diff --git a/CommandLines.cpp b/CommandLines.cpp
index 75048c6..b773a54 100644
--- a/CommandLines.cpp
+++ b/CommandLines.cpp
@@ -34,6 +34,10 @@ static ko_longopt_t long_options[] = {
     { "m-rate",      ko_required_argument, 319 },
     { "b-partition",      ko_no_argument, 320 },
     { "t-occ",      ko_required_argument, 321 },
+    { "seed",      ko_required_argument, 322 },
+    { "n-perturb",      ko_required_argument, 323 },
+    { "f-perturb",      ko_required_argument, 324 },
+    { "n-hap",      ko_required_argument, 325 },
 	{ 0, 0, 0 }
 };
 
@@ -49,64 +53,71 @@ void Print_H(hifiasm_opt_t* asm_opt)
     fprintf(stderr, "Usage: hifiasm [options] <in_1.fq> <in_2.fq> <...>\n");
     fprintf(stderr, "Options:\n");
 	fprintf(stderr, "  Input/Output:\n");
-    fprintf(stderr, "    -o STR      prefix of output files [%s]\n", asm_opt->output_file_name);
-    fprintf(stderr, "    -i          ignore saved read correction and overlaps\n");
-    fprintf(stderr, "    -t INT      number of threads [%d]\n", asm_opt->thread_num);
-    fprintf(stderr, "    -z INT      length of adapters that should be removed [%d]\n", asm_opt->adapterLen);
-    fprintf(stderr, "    --version   show version number\n");
+    fprintf(stderr, "    -o STR       prefix of output files [%s]\n", asm_opt->output_file_name);
+    fprintf(stderr, "    -i           ignore saved read correction and overlaps\n");
+    fprintf(stderr, "    -t INT       number of threads [%d]\n", asm_opt->thread_num);
+    fprintf(stderr, "    -z INT       length of adapters that should be removed [%d]\n", asm_opt->adapterLen);
+    fprintf(stderr, "    --version    show version number\n");
 	fprintf(stderr, "  Overlap/Error correction:\n");
-    fprintf(stderr, "    -k INT      k-mer length (must be <64) [%d]\n", asm_opt->k_mer_length);
-	fprintf(stderr, "    -w INT      minimizer window size [%d]\n", asm_opt->mz_win);
-	fprintf(stderr, "    -f INT      number of bits for bloom filter; 0 to disable [%d]\n", asm_opt->bf_shift);
-	fprintf(stderr, "    -D FLOAT    drop k-mers occurring >FLOAT*coverage times [%.1f]\n", asm_opt->high_factor);
-	fprintf(stderr, "    -N INT      consider up to max(-D*coverage,-N) overlaps for each oriented read [%d]\n", asm_opt->max_n_chain);
-    fprintf(stderr, "    -r INT      round of correction [%d]\n", asm_opt->number_of_round);
+    fprintf(stderr, "    -k INT       k-mer length (must be <64) [%d]\n", asm_opt->k_mer_length);
+	fprintf(stderr, "    -w INT       minimizer window size [%d]\n", asm_opt->mz_win);
+	fprintf(stderr, "    -f INT       number of bits for bloom filter; 0 to disable [%d]\n", asm_opt->bf_shift);
+	fprintf(stderr, "    -D FLOAT     drop k-mers occurring >FLOAT*coverage times [%.1f]\n", asm_opt->high_factor);
+	fprintf(stderr, "    -N INT       consider up to max(-D*coverage,-N) overlaps for each oriented read [%d]\n", asm_opt->max_n_chain);
+    fprintf(stderr, "    -r INT       round of correction [%d]\n", asm_opt->number_of_round);
     fprintf(stderr, "  Assembly:\n");
-    fprintf(stderr, "    -a INT      round of assembly cleaning [%d]\n", asm_opt->clean_round);
-    fprintf(stderr, "    -m INT      pop bubbles of <INT in size in contig graphs [%lld]\n", asm_opt->large_pop_bubble_size);
-    fprintf(stderr, "    -p INT      pop bubbles of <INT in size in unitig graphs [%lld]\n", asm_opt->small_pop_bubble_size);
-    fprintf(stderr, "    -n INT      remove tip unitigs composed of <=INT reads [%d]\n", asm_opt->max_short_tip);
-    fprintf(stderr, "    -x FLOAT    max overlap drop ratio [%.2g]\n", asm_opt->max_drop_rate);
-    fprintf(stderr, "    -y FLOAT    min overlap drop ratio [%.2g]\n", asm_opt->min_drop_rate);
-    fprintf(stderr, "    -u          disable post join contigs step which may improve N50\n");
-    fprintf(stderr, "    --lowQ      INT\n");
-    fprintf(stderr, "                output contig regions with >=INT%% inconsistency in BED format; 0 to disable [%d]\n", asm_opt->bed_inconsist_rate);
-    fprintf(stderr, "    --b-cov     INT\n");
-    fprintf(stderr, "                break contigs at positions with <INT-fold coverage; work with '--m-rate'; 0 to disable [%d]\n", asm_opt->b_low_cov);
-    fprintf(stderr, "    --h-cov     INT\n");
-    fprintf(stderr, "                break contigs at positions with >INT-fold coverage; work with '--m-rate'; -1 to disable [%d]\n", asm_opt->b_high_cov);
-    fprintf(stderr, "    --m-rate    FLOAT\n");
-    fprintf(stderr, "                break contigs at positions with <=FLOAT*coverage exact overlaps;\n");
-    fprintf(stderr, "                only work with '--b-cov' or '--h-cov'[%.2f]\n", asm_opt->m_rate);
+    fprintf(stderr, "    -a INT       round of assembly cleaning [%d]\n", asm_opt->clean_round);
+    fprintf(stderr, "    -m INT       pop bubbles of <INT in size in contig graphs [%lld]\n", asm_opt->large_pop_bubble_size);
+    fprintf(stderr, "    -p INT       pop bubbles of <INT in size in unitig graphs [%lld]\n", asm_opt->small_pop_bubble_size);
+    fprintf(stderr, "    -n INT       remove tip unitigs composed of <=INT reads [%d]\n", asm_opt->max_short_tip);
+    fprintf(stderr, "    -x FLOAT     max overlap drop ratio [%.2g]\n", asm_opt->max_drop_rate);
+    fprintf(stderr, "    -y FLOAT     min overlap drop ratio [%.2g]\n", asm_opt->min_drop_rate);
+    fprintf(stderr, "    -u           disable post join contigs step which may improve N50\n");
+    fprintf(stderr, "    --lowQ       INT\n");
+    fprintf(stderr, "                 output contig regions with >=INT%% inconsistency in BED format; 0 to disable [%d]\n", asm_opt->bed_inconsist_rate);
+    fprintf(stderr, "    --b-cov      INT\n");
+    fprintf(stderr, "                 break contigs at positions with <INT-fold coverage; work with '--m-rate'; 0 to disable [%d]\n", asm_opt->b_low_cov);
+    fprintf(stderr, "    --h-cov      INT\n");
+    fprintf(stderr, "                 break contigs at positions with >INT-fold coverage; work with '--m-rate'; -1 to disable [%d]\n", asm_opt->b_high_cov);
+    fprintf(stderr, "    --m-rate     FLOAT\n");
+    fprintf(stderr, "                 break contigs at positions with <=FLOAT*coverage exact overlaps;\n");
+    fprintf(stderr, "                 only work with '--b-cov' or '--h-cov'[%.2f]\n", asm_opt->m_rate);
     
 //	fprintf(stderr, "    --pri-range INT1[,INT2]\n");
 //	fprintf(stderr, "                keep contigs with coverage in this range in p_ctg.gfa; -1 to disable [auto,inf]\n");
 
     fprintf(stderr, "  Trio-partition:\n");
-    fprintf(stderr, "    -1 FILE     hap1/paternal k-mer dump generated by \"yak count\" []\n");
-    fprintf(stderr, "    -2 FILE     hap2/maternal k-mer dump generated by \"yak count\" []\n");
-    fprintf(stderr, "    -c INT      lower bound of the binned k-mer's frequency [%d]\n", asm_opt->min_cnt);
-    fprintf(stderr, "    -d INT      upper bound of the binned k-mer's frequency [%d]\n", asm_opt->mid_cnt);
-	fprintf(stderr, "    -3 FILE     list of hap1/paternal read names []\n");
-	fprintf(stderr, "    -4 FILE     list of hap2/maternal read names []\n");
-    fprintf(stderr, "    --t-occ     INT\n");
-    fprintf(stderr, "                force remove unitigs with >INT unexpected haplotype-specific reads;\n");
-    fprintf(stderr, "                ignore graph topology; [%d]\n", asm_opt->trio_flag_occ_thres);
+    fprintf(stderr, "    -1 FILE      hap1/paternal k-mer dump generated by \"yak count\" []\n");
+    fprintf(stderr, "    -2 FILE      hap2/maternal k-mer dump generated by \"yak count\" []\n");
+    fprintf(stderr, "    -c INT       lower bound of the binned k-mer's frequency [%d]\n", asm_opt->min_cnt);
+    fprintf(stderr, "    -d INT       upper bound of the binned k-mer's frequency [%d]\n", asm_opt->mid_cnt);
+	fprintf(stderr, "    -3 FILE      list of hap1/paternal read names []\n");
+	fprintf(stderr, "    -4 FILE      list of hap2/maternal read names []\n");
+    fprintf(stderr, "    --t-occ      INT\n");
+    fprintf(stderr, "                 force remove unitigs with >INT unexpected haplotype-specific reads;\n");
+    fprintf(stderr, "                 ignore graph topology; [%d]\n", asm_opt->trio_flag_occ_thres);
 
     fprintf(stderr, "  Purge-dups:\n");
-    fprintf(stderr, "    -l INT      purge level. 0: no purging; 1: light; 2/3: aggressive [0 for trio; 2 for unzip]\n");
-    fprintf(stderr, "    -s FLOAT    similarity threshold for duplicate haplotigs [%g for -l1/-l2, %g for -l3]\n", 
-                                     asm_opt->purge_simi_rate_l2, asm_opt->purge_simi_rate_l3);
-    fprintf(stderr, "    -O INT      min number of overlapped reads for duplicate haplotigs [%d]\n", 
-                                     asm_opt->purge_overlap_len);
-    fprintf(stderr, "    --purge-cov INT\n");
-    fprintf(stderr, "                coverage upper bound of Purge-dups [auto]\n");
-
-    ///fprintf(stderr, "    --high-het  enable this mode for high heterozygosity sample [experimental, not stable]\n");
+    fprintf(stderr, "    -l INT       purge level. 0: no purging; 1: light; 2/3: aggressive [0 for trio; 3 for unzip]\n");
+    fprintf(stderr, "    -s FLOAT     similarity threshold for duplicate haplotigs [%g for -l1/-l2, %g for -l3]\n", 
+                                      asm_opt->purge_simi_rate_l2, asm_opt->purge_simi_rate_l3);
+    fprintf(stderr, "    -O INT       min number of overlapped reads for duplicate haplotigs [%d]\n", 
+                                      asm_opt->purge_overlap_len);
+    fprintf(stderr, "    --purge-cov  INT\n");
+    fprintf(stderr, "                 coverage upper bound of Purge-dups [auto]\n");
+    fprintf(stderr, "    --n-hap      INT\n");
+    fprintf(stderr, "                 number of haplotypes [%d]\n", asm_opt->polyploidy);
     
-    fprintf(stderr, "  Hi-C-partition [experimental, not stable]:\n");
+    // fprintf(stderr, "  Hi-C-partition [experimental, not stable]:\n");
+    fprintf(stderr, "  Hi-C-partition:\n");
     fprintf(stderr, "    --h1 FILEs   file names of Hi-C R1  [r1_1.fq,r1_2.fq,...]\n");
     fprintf(stderr, "    --h2 FILEs   file names of Hi-C R2  [r2_1.fq,r2_2.fq,...]\n");
+    fprintf(stderr, "    --seed INT   RNG seed [%lu]\n", asm_opt->seed);
+	fprintf(stderr, "    --n-perturb  INT\n");
+    fprintf(stderr, "                 rounds of perturbation [%d]\n", asm_opt->n_perturb);
+    fprintf(stderr, "    --f-perturb  FLOAT\n");
+	fprintf(stderr, "                 fraction to flip for perturbation [%.3g]\n", asm_opt->f_perturb);
+
 
     fprintf(stderr, "Example: ./hifiasm -o NA12878.asm -t 32 NA12878.fq.gz\n");
     fprintf(stderr, "See `man ./hifiasm.1' for detailed description of these command-line options.\n");
@@ -154,7 +165,7 @@ void init_opt(hifiasm_opt_t* asm_opt)
     asm_opt->max_short_tip = 3;
     asm_opt->min_cnt = 2;
     asm_opt->mid_cnt = 5;
-    asm_opt->purge_level_primary = 2;
+    asm_opt->purge_level_primary = 3;
     asm_opt->purge_level_trio = 0;
     asm_opt->purge_simi_rate_l2 = 0.75;
     asm_opt->purge_simi_rate_l3 = 0.55;
@@ -175,6 +186,9 @@ void init_opt(hifiasm_opt_t* asm_opt)
     asm_opt->hap_occ = 1;
     asm_opt->polyploidy = 2;
     asm_opt->trio_flag_occ_thres = 60;
+    asm_opt->seed = 11;
+    asm_opt->n_perturb = 50000;
+    asm_opt->f_perturb = 0.1;
 }
 
 void destory_enzyme(enzyme* f)
@@ -469,27 +483,7 @@ int check_option(hifiasm_opt_t* asm_opt)
         fprintf(stderr, "[ERROR] [-s] must >= 0\n");
         return 0;
     }
-
-    // fprintf(stderr, "input file num: %d\n", asm_opt->num_reads);
-    // fprintf(stderr, "output file: %s\n", asm_opt->output_file_name);
-    // fprintf(stderr, "number of threads: %d\n", asm_opt->thread_num);
-    // fprintf(stderr, "number of rounds for correction: %d\n", asm_opt->number_of_round);
-    // fprintf(stderr, "number of rounds for assembly cleaning: %d\n", asm_opt->clean_round);
-    // fprintf(stderr, "length of removed adapters: %d\n", asm_opt->adapterLen);
-    // fprintf(stderr, "length of k_mer: %d\n", asm_opt->k_mer_length);
-    // fprintf(stderr, "min overlap drop ratio: %.2g\n", asm_opt->min_drop_rate);
-    // fprintf(stderr, "max overlap drop ratio: %.2g\n", asm_opt->max_drop_rate);
-    // fprintf(stderr, "size of popped small bubbles: %lld\n", asm_opt->small_pop_bubble_size);
-    // fprintf(stderr, "size of popped large bubbles: %lld\n", asm_opt->large_pop_bubble_size);
-    // fprintf(stderr, "small removed unitig threshold: %d\n", asm_opt->max_short_tip);
-    // fprintf(stderr, "small removed unitig threshold: %d\n", asm_opt->max_short_tip);
-    // fprintf(stderr, "min_cnt: %d\n", asm_opt->min_cnt);
-    // fprintf(stderr, "mid_cnt: %d\n", asm_opt->mid_cnt);
-    // fprintf(stderr, "purge_level_primary: %d\n", asm_opt->purge_level_primary);
-    // fprintf(stderr, "purge_level_trio: %d\n", asm_opt->purge_level_trio);
-    // fprintf(stderr, "purge_simi_thres: %f\n", asm_opt->purge_simi_thres);
-    // fprintf(stderr, "purge_overlap_len: %d\n", asm_opt->purge_overlap_len);
-
+    
     return 1;
 }
 
@@ -652,6 +646,10 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
         else if (c == 319) asm_opt->m_rate = atof(opt.arg);
         else if (c == 320) asm_opt->flag |= HA_F_PARTITION;
         else if (c == 321) asm_opt->trio_flag_occ_thres = atoi(opt.arg);
+        else if (c == 322) asm_opt->seed = atol(opt.arg);
+        else if (c == 323) asm_opt->n_perturb = atoi(opt.arg);
+        else if (c == 324) asm_opt->f_perturb = atof(opt.arg);
+        else if (c == 325) asm_opt->polyploidy = atoi(opt.arg);
         else if (c == 'l')
         {   ///0: disable purge_dup; 1: purge containment; 2: purge overlap
             asm_opt->purge_level_primary = asm_opt->purge_level_trio = atoi(opt.arg);
diff --git a/CommandLines.h b/CommandLines.h
index 44c2bdc..aac7329 100644
--- a/CommandLines.h
+++ b/CommandLines.h
@@ -2,8 +2,9 @@
 #define __COMMAND_LINE_PARSER__
 
 #include <pthread.h>
+#include <stdint.h>
 
-#define HA_VERSION "0.14.2-r326"
+#define HA_VERSION "0.15-r327"
 
 #define VERBOSE 0
 
@@ -97,7 +98,9 @@ typedef struct {
     int hap_occ;
     int polyploidy;
     int trio_flag_occ_thres;
-    
+    uint64_t seed;
+    int32_t n_perturb;
+    double f_perturb;
 } hifiasm_opt_t;
 
 extern hifiasm_opt_t asm_opt;
diff --git a/Overlaps.cpp b/Overlaps.cpp
index 8e3c5f0..fe372f0 100644
--- a/Overlaps.cpp
+++ b/Overlaps.cpp
@@ -12578,7 +12578,7 @@ trans_chain* load_hc_hits(const char *fn)
 
     fclose(fp);
     free(buf);
-    fprintf(stderr, "[M::%s::] ==> Hi-C cov have been loaded\n", __func__);
+    // fprintf(stderr, "[M::%s::] ==> Hi-C cov have been loaded\n", __func__);
     return t_ch;
 }
 
@@ -12630,30 +12630,20 @@ bub_label_t* b_mask_t)
         write_trans_chain(cov->t_ch, output_file_name);
     }
 
-    
     hic_analysis(ug, sg, cov?cov->t_ch:t_ch);
 
-    fprintf(stderr, "sa-0-sa\n");
-
-
+    /**
     char* gfa_name = (char*)malloc(strlen(output_file_name)+25);
     sprintf(gfa_name, "%s.d_utg.noseq.gfa", output_file_name);
     FILE* output_file = fopen(gfa_name, "w");
     ma_ug_print_simple(ug, &R_INF, sg, coverage_cut, sources, ruIndex, "utg", output_file);
     fclose(output_file);
     free(gfa_name);
-
-
-    fprintf(stderr, "sa-1-sa\n");
-
-
-
-
+    **/
 
     if(cov) destory_hap_cov_t(&cov);
     if(t_ch) destory_trans_chain(&t_ch);
 
-    fprintf(stderr, "sa-2-sa\n");
     ma_ug_destroy(ug);
     kv_destroy(new_rtg_edges.a); 
     
diff --git a/Overlaps.h b/Overlaps.h
index 3bdb77f..6caf5c2 100644
--- a/Overlaps.h
+++ b/Overlaps.h
@@ -115,6 +115,7 @@ typedef struct {
 typedef struct {
 	uint32_t qs, qe, qn;
 	uint32_t ts, te, tn;
+	uint32_t occ;
 	double nw;
 	uint8_t f:6, rev:1, del:1;
 } u_trans_t;
diff --git a/Purge_Dups.cpp b/Purge_Dups.cpp
index 7835452..4fa7f18 100644
--- a/Purge_Dups.cpp
+++ b/Purge_Dups.cpp
@@ -5237,7 +5237,7 @@ uint32_t just_coverage, hap_cov_t *cov, uint32_t collect_p_trans)
 
     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);
+        mc_solve(&all_ovlp, cov->t_ch, NULL, ug, read_g, 0.8, R_INF.trio_flag, 1, NULL, 1, NULL);
         ///pt_solve(&all_ovlp, cov->t_ch, ug, read_g, 0.8, R_INF.trio_flag);
     } 
     
diff --git a/hic.cpp b/hic.cpp
index d98f2e5..3b72cf4 100644
--- a/hic.cpp
+++ b/hic.cpp
@@ -190,6 +190,7 @@ typedef struct {
 typedef struct {
     kvec_t(pe_hit) a;
     kvec_t(uint64_t) idx;
+    kvec_t(uint64_t) occ;
 } kvec_pe_hit;
 
 typedef struct {
@@ -523,7 +524,7 @@ int write_hc_pt_index(ha_ug_index* idx, char* file_name)
 int load_hc_pt_index(ha_ug_index** r_idx, char* file_name)
 {
     uint64_t flag = 0;
-    double index_time = yak_realtime();
+    // double index_time = yak_realtime();
     char* gfa_name = (char*)malloc(strlen(file_name)+25);
     sprintf(gfa_name, "%s.hic.tlb.bin", file_name);
     FILE* fp = fopen(gfa_name, "r");
@@ -558,7 +559,7 @@ int load_hc_pt_index(ha_ug_index** r_idx, char* file_name)
 
     free(gfa_name);
     fclose(fp);
-    fprintf(stderr, "[M::%s::%.3f] ==> HiC index has been loaded\n", __func__, yak_realtime()-index_time);
+    // fprintf(stderr, "[M::%s::%.3f] ==> HiC index has been loaded\n", __func__, yak_realtime()-index_time);
     return 1;
 }
 
@@ -3226,7 +3227,6 @@ void get_shortest_path(uint32_t src, pdq* pq, asg_t *sg, uint32_t* pre)
 
 void all_pair_shortest_path(asg_t *sg, hc_links* link, MT* M)
 {
-    double index_time = yak_realtime();
     hc_linkeage* t = NULL;
     pdq pq;
     init_pdq(&pq, sg->n_seq<<1);
@@ -3249,7 +3249,6 @@ void all_pair_shortest_path(asg_t *sg, hc_links* link, MT* M)
     }
     
     destory_pdq(&pq);
-    fprintf(stderr, "[M::%s::%.3f]\n", __func__, yak_realtime()-index_time);
 }
 
 uint64_t LCA_distance(long long d_x, long long d_y, long long xLen, long long yLen, uint8_t* rev)
@@ -3512,7 +3511,7 @@ static void worker_for_dis(void *data, long i, int tid)
 
 void fill_utg_distance_multi(asg_t *sg, hc_links* link, MT* M, bubble_type* bub)
 {
-    double index_time = yak_realtime();
+    // double index_time = yak_realtime();
     uint32_t i;
     utg_d_t s;
     s.sg = sg; s.link = link; s.M = M; s.bub = bub;
@@ -3530,7 +3529,7 @@ void fill_utg_distance_multi(asg_t *sg, hc_links* link, MT* M, bubble_type* bub)
         free(s.dis_buf[i]);
     }
     free(s.dis_buf);
-    fprintf(stderr, "[M::%s::%.3f]\n", __func__, yak_realtime()-index_time);
+    // fprintf(stderr, "[M::%s::%.3f]\n", __func__, yak_realtime()-index_time);
 }
 
 void init_MT(MT* M, uint32_t n_vtx)
@@ -4168,7 +4167,7 @@ void collect_hc_links(const ha_ug_index* idx, kvec_pe_hit* hits, hc_links* link,
 
 void measure_distance(const ma_ug_t* ug, kvec_pe_hit* hits, hc_links* link, bubble_type* bub, kv_u_trans_t *ta)
 {
-    double index_time = yak_realtime();
+    // double index_time = yak_realtime();
     MT M;
     init_MT(&M, ug->g->n_seq<<1);
     uint64_t uID_bits;
@@ -4217,7 +4216,7 @@ void measure_distance(const ma_ug_t* ug, kvec_pe_hit* hits, hc_links* link, bubb
 
 
     destory_MT(&M);
-    fprintf(stderr, "[M::%s::%.3f] ==> Hi-C linkages have been counted\n", __func__, yak_realtime()-index_time);
+    // fprintf(stderr, "[M::%s::%.3f] ==> Hi-C linkages have been counted\n", __func__, yak_realtime()-index_time);
     return;
 }
 
@@ -4669,7 +4668,7 @@ int load_hc_hits(kvec_pe_hit* hits, const char *fn)
 
     fclose(fp);
     free(buf);
-    fprintf(stderr, "[M::%s::] ==> Hi-C linkages have been loaded\n", __func__);
+    // fprintf(stderr, "[M::%s::] ==> Hi-C linkages have been loaded\n", __func__);
     return 1;
 }
 
@@ -8783,12 +8782,12 @@ void clean_bubble_chain_by_HiC(ma_ug_t* ug, hc_links* link, bubble_type* bub)
     kvec_asg_arc_t_warp edges; kv_init(edges.a);
     ma_ug_t *back_bs_ug = copy_untig_graph(bs_ug);
 
-    for (i = 0; i < bs_ug->g->n_arc; i++)
+    for (i = 0; i < bs_ug->g->n_arc; i++)///weight of bs_ug's edges
     {
         e_w[i] = -1;
     }
 
-    for (i = 0; i < bs_ug->g->n_seq; i++)
+    for (i = 0; i < bs_ug->g->n_seq; i++)///init all chain with flag_aux 
     {
         set_b_utg_weight_flag(bub, &b, i<<1, vis, flag_aux, NULL);
     }
@@ -9563,12 +9562,12 @@ H_partition* hap, int8_t *s, trans_idx* dis)
         dis->a[i+1].beg = dis->a[i].end;
     }
     
-    for (i = 0; i < dis->n; i++)
-    {
-        if(i > 0 && dis->a[i].beg != dis->a[i-1].end) fprintf(stderr, "ERROR: dis->a[i].beg: %lu, dis->a[i-1].end: %lu\n", dis->a[i].beg, dis->a[i-1].end);
-        fprintf(stderr, "beg: %lu, end: %lu, cnt_0: %lu, cnt_1: %lu, error_rate: %f\n", 
-        dis->a[i].beg, dis->a[i].end, dis->a[i].cnt_0, dis->a[i].cnt_1, (double)(dis->a[i].cnt_1)/(double)(dis->a[i].cnt_1 + dis->a[i].cnt_0));
-    }
+    // for (i = 0; i < dis->n; i++)
+    // {
+    //     if(i > 0 && dis->a[i].beg != dis->a[i-1].end) fprintf(stderr, "ERROR: dis->a[i].beg: %lu, dis->a[i-1].end: %lu\n", dis->a[i].beg, dis->a[i-1].end);
+    //     fprintf(stderr, "beg: %lu, end: %lu, cnt_0: %lu, cnt_1: %lu, error_rate: %f\n", 
+    //     dis->a[i].beg, dis->a[i].end, dis->a[i].cnt_0, dis->a[i].cnt_1, (double)(dis->a[i].cnt_1)/(double)(dis->a[i].cnt_1 + dis->a[i].cnt_0));
+    // }
     
     LeastSquare_advance(dis, idx, med);
     // fprintf(stderr, "idx->a: %f, idx->b: %f, idx->frac: %f, med: %lu\n", 
@@ -11405,7 +11404,7 @@ typedef struct{
     block_phase_type* x;
     uint32_t n_thread;
     bubble_type* bub;
-    uint64_t* chain_idx;
+    uint64_t* chain_idx;///index of each chain
     uint64_t chain_idx_n;
     uint64_t chain_ele_occ;
     block_res_type* res;
@@ -11469,7 +11468,7 @@ void init_mul_block_phase_type(mul_block_phase_type* x, G_partition* g_p, bubble
         n += shift_block_phase_type(u, g_p, bub, &b, (uint32_t)-1);
         x->chain_idx_n++;
     }
-    x->chain_idx[i] = n;
+    x->chain_idx[i] = n;///index of chain
     x->chain_ele_occ = n;
 }
 
@@ -12401,7 +12400,7 @@ uint32_t init_contig_partition(H_partition* hap, ha_ug_index* idx, bubble_type*
 
     link_phase_group(hap, bub);
 
-    assign_per_unitig_G_partition(&(hap->g_p), hap->n, link, bub, 0);
+    assign_per_unitig_G_partition(&(hap->g_p), hap->n, link, bub, 0);///warp unitigs
 
     adjust_contig_partition(hap, link);
 
@@ -12419,7 +12418,1083 @@ uint32_t init_contig_partition(H_partition* hap, ha_ug_index* idx, bubble_type*
     return 1;
 }
 
+double get_path_phasing_weight(uint32_t query, uint32_t v0, uint32_t root, bub_p_t_warp *b, kv_u_trans_t *ta)
+{
+    if(v0 == root) return 0;
+    uint32_t v, u;
+    u_trans_t *p = NULL;
+    double nw = 0;
+    v = v0;
+    do {
+        u = b->a[v].p; // u->v
+        get_u_trans_spec(ta, query>>1, v>>1, &p, NULL);
+        if(p && (!p->del)) nw += p->nw;
+        v = u;
+    } while (v != root);
 
+    return nw;
+}
+
+void get_related_phasing_weight(uint32_t x, kv_u_trans_t *ta, double* w0, double* w1, int8_t *s)
+{
+    (*w0) = (*w1) = 0;
+    if(x >= ta->idx.n) return;
+    uint32_t e_n, k;
+    u_trans_t* e = u_trans_a(*ta, x);
+    e_n = u_trans_n(*ta, x);
+    for (k = 0; k < e_n; k++)
+    {
+        if(e[k].del) continue;
+        if(s[e[k].tn] > 0) (*w0) += e[k].nw;
+        else if(s[e[k].tn] < 0) (*w1) += e[k].nw;
+    }
+}
+
+void set_phase_path(bub_p_t_warp *b, uint32_t root, kv_u_trans_t *ta, ps_t *s)
+{
+    int8_t f;
+    uint32_t v, u;
+    double z[2], cur_w[2];
+    z[0] = z[1] = 0;
+
+    v = b->S.a[0];
+    do {
+        u = b->a[v].p; // u->v
+        if(v != b->S.a[0]) 
+        {
+            get_related_phasing_weight(v>>1, ta, &cur_w[0], &cur_w[1], s->s);
+            z[0] += cur_w[0]; z[1] += cur_w[1];
+        }
+        v = u;
+    } while (v != root);
+
+    if(z[0] - z[1] < 0)
+    {
+        f = 1;
+    }
+    else if(z[0] - z[1] > 0)
+    {
+        f = -1;
+    }
+    else
+    {
+        s->xs = kr_splitmix64(s->xs);
+        f = s->xs&1? 1 : -1;
+    }
+
+    v = b->S.a[0];
+    do {
+        u = b->a[v].p; // u->v
+        if(v != b->S.a[0]) s->s[v>>1] = f;
+        v = u;
+    } while (v != root);
+}
+
+uint32_t bub_phase(ma_ug_t *ug, uint32_t beg, uint32_t end, bub_p_t_warp *b, kv_u_trans_t *ta, ps_t *s)
+{
+    asg_t *g = ug->g;
+    if(g->seq[beg>>1].del) return 0; // already deleted
+    if(get_real_length(g, beg, NULL)<2) return 0;
+    uint32_t i, is_end, n_pending, to_replace, cur_nc, cur_uc, cur_ac, n_tips, tip_end, n_pop;
+    double cur_nh, cur_w0, cur_w1, cur_rate, max_rate, cur_weight, min_weight;
+    ///S saves nodes with all incoming edges visited
+    b->S.n = b->T.n = b->b.n = b->e.n = 0;
+    ///for each node, b->a saves all related information
+    b->a[beg].d = b->a[beg].nc = b->a[beg].ac = b->a[beg].uc = 0; 
+    b->a[beg].nh = b->a[beg].w[0] = b->a[beg].w[1] = 0;
+    b->a[beg].p = (uint32_t)-1;
+    ///b->S is the nodes with all incoming edges visited
+    kv_push(uint32_t, b->S, beg);
+    n_pop = n_tips = n_pending = 0;
+    tip_end = (uint32_t)-1;
+
+    do {
+        ///v is a node that all incoming edges have been visited
+        ///d is the distance from v0 to v
+        uint32_t v = kv_pop(b->S);
+        uint32_t d = b->a[v].d, nc = b->a[v].nc, uc = b->a[v].uc, ac = b->a[v].ac; 
+        double nh = b->a[v].nh;///path weight
+        double nw_0 = b->a[v].w[0], nw_1 = b->a[v].w[1];///weight to haplotype 1/2
+
+        uint32_t nv = asg_arc_n(g, v);
+        asg_arc_t *av = asg_arc_a(g, v);
+        for (i = 0; i < nv; ++i) {
+            uint32_t w = av[i].v, l = (uint32_t)av[i].ul; // v->w with length l, not overlap length
+            bub_p_t *t = &b->a[w];
+            is_end = 0;
+            if((w>>1) == (end>>1)) is_end = 1;
+            //got a circle
+            if ((w>>1) == (beg>>1)) goto pop_reset;
+            if (av[i].del) continue;
+            ///push the edge
+            kv_push(uint32_t, b->e, (g->idx[v]>>32) + i);
+
+            if (t->s == 0) 
+            { // this vertex has never been visited
+                kv_push(uint32_t, b->b, w); // save it for revert
+                ///t->p is the parent node of 
+                ///t->s = 1 means w has been visited
+                ///d is len(v0->v), l is len(v->w), so t->d is len(v0->w)
+                t->p = v, t->s = 1, t->d = d + l; 
+                t->r = get_real_length(g, w^1, NULL);
+
+                if(is_end == 0)
+                {
+                    t->nc = nc + ug->u.a[(w>>1)].n;
+                    t->nh = nh + get_path_phasing_weight(w, v, beg, b, ta);
+                    get_related_phasing_weight(w>>1, ta, &(t->w[0]), &(t->w[1]), s->s);
+                    t->w[0] += nw_0; t->w[1] += nw_1; 
+                    t->ac = ac + ((s->s[w>>1] == 0)? ug->u.a[(w>>1)].n : 0);
+                    t->uc = uc + ((s->s[w>>1] != 0)? ug->u.a[(w>>1)].n : 0);
+                }
+
+                ++n_pending;
+            }
+            else {
+                to_replace = 0;
+
+                if(is_end)
+                {
+                    cur_nc = nc; cur_nh = nh; 
+                    cur_w0 = nw_0; cur_w1= nw_1;
+                    cur_ac = ac; cur_uc = uc;
+                }
+                else
+                {
+                    cur_nc = nc + ug->u.a[(w>>1)].n;
+                    cur_nh = nh + get_path_phasing_weight(w, v, beg, b, ta);
+                    get_related_phasing_weight(w>>1, ta, &cur_w0, &cur_w1, s->s);
+                    cur_w0 += nw_0; cur_w1 += nw_1;
+                    cur_ac = ac + ((s->s[w>>1] == 0)? ug->u.a[(w>>1)].n : 0);
+                    cur_uc = uc + ((s->s[w>>1] != 0)? ug->u.a[(w>>1)].n : 0);
+                }
+                
+                cur_weight = cur_nh + MIN(cur_w0, cur_w1) - MAX(cur_w0, cur_w1);
+                min_weight = t->nh + MIN(t->w[0], t->w[1]) - MAX(t->w[0], t->w[1]);
+                cur_rate = ((cur_ac+cur_uc == 0)? -1 : ((double)(cur_ac)/(double)(cur_ac+cur_uc)));
+                max_rate = ((t->ac+t->uc == 0)? -1 : ((double)(t->ac)/(double)(t->ac+t->uc)));
+
+                if(cur_rate > max_rate)
+                {
+                    to_replace = 1;
+                }
+                else if(cur_rate == max_rate)
+                {
+                    if(cur_weight < min_weight)
+                    {
+                        to_replace = 1;
+                    }
+                    else if(cur_weight == min_weight)
+                    {
+                        if(cur_nc > t->nc)
+                        {
+                            to_replace = 1;
+                        }
+                        else if(cur_nc == t->nc)
+                        {
+                            if(d + l > t->d)
+                            {
+                                to_replace = 1;
+                            }
+                        }
+                    }
+                }
+
+
+                if(to_replace)
+                {
+                    t->p = v;
+                    t->nc = cur_nc;
+                    t->nh = cur_nh;
+                    t->ac = cur_ac;
+                    t->uc = cur_uc;
+                    t->w[0] = cur_w0;
+                    t->w[1] = cur_w1;
+                }
+
+
+                if (d + l < t->d) t->d = d + l; // update dist
+            }
+
+            if (--(t->r) == 0) {
+                uint32_t x = get_real_length(g, w, NULL);
+                if(x > 0)
+                {
+                    kv_push(uint32_t, b->S, w);
+                }
+                else
+                {
+                    ///at most one tip
+                    if(n_tips != 0) goto pop_reset;
+                    n_tips++;
+                    tip_end = w;
+                }
+                --n_pending;
+            }
+        }
+
+        if(n_tips == 1)
+        {
+            if(tip_end != (uint32_t)-1 && n_pending == 0 && b->S.n == 0)
+            {
+                ///sink is b.S.a[0]
+                kv_push(uint32_t, b->S, tip_end);
+                break;
+            }
+            else
+            {
+                goto pop_reset;
+            }
+        }
+
+        if (i < nv || b->S.n == 0) goto pop_reset;
+    }while (b->S.n > 1 || n_pending);
+
+    n_pop = 1;
+    /**need fix**/
+    set_phase_path(b, beg, ta, s);
+
+    pop_reset:
+
+    for (i = 0; i < b->b.n; ++i) { // clear the states of visited vertices
+        bub_p_t *t = &b->a[b->b.a[i]];
+        t->p = t->d = t->nc = t->ac = t->uc = t->r = t->s = 0;
+        t->nh = t->w[0] = t->w[1] = 0;
+    }
+    return n_pop;
+}
+
+uint32_t get_weightest_node(kv_u_trans_t *ta, bubble_type* bub, ma_ug_t* ug, int8_t *s, uint8_t *vis)
+{
+    double w_a, w_n, max_w_a, max_w_n;
+    uint32_t i, occ, k, m, v, *a = NULL, a_n, e_n, max_w_a_i, max_w_n_i;
+    u_trans_t *e = NULL;
+    max_w_a = max_w_n = -1; max_w_a_i = max_w_n_i = (uint32_t)-1;
+    for (i = 0; i < bub->f_bub; i++)
+    {
+        get_bubbles(bub, i, NULL, NULL, &a, &a_n, NULL);
+        if(vis[i]) continue;
+        for (k = 0, w_a = w_n = 0; k < a_n; k++)
+        {
+            v = a[k]>>1;
+            if(s[v] != 0) break; 
+            e = u_trans_a(*ta, v);
+            e_n = u_trans_n(*ta, v);
+            for (m = 0; m < e_n; m++)
+            {
+                if(s[e[m].tn] != 0)
+                {
+                    w_a += (e[m].nw>=0?e[m].nw:-e[m].nw);
+                }
+                else
+                {
+                    w_n += (e[m].nw>=0?e[m].nw:-e[m].nw);
+                }
+            }
+        }
+
+        if(k >= a_n && a_n > 0)//unset whole bubble
+        {
+            if(w_a > max_w_a)
+            {
+                max_w_a_i = i<<1;
+                max_w_a = w_a;
+            }
+
+            if(w_n > max_w_n)
+            {
+                max_w_n_i = i<<1;
+                max_w_n = w_n;
+            }
+        }
+        else
+        {
+            for (k = occ = 0; k < a_n; k++)
+            {
+                v = a[k]>>1;
+                if(s[v] != 0)
+                {
+                    occ++;
+                    continue;
+                } 
+                
+                e = u_trans_a(*ta, v);
+                e_n = u_trans_n(*ta, v);
+                for (m = 0, w_a = w_n = 0; m < e_n; m++)
+                {
+                    if(s[e[m].tn] != 0)
+                    {
+                        w_a += (e[m].nw>=0?e[m].nw:-e[m].nw);
+                    }
+                    else
+                    {
+                        w_n += (e[m].nw>=0?e[m].nw:-e[m].nw);
+                    }
+                }
+
+                if(w_a > max_w_a)
+                {
+                    max_w_a_i = (v<<1)+1;
+                    max_w_a = w_a;
+                }
+
+                if(w_n > max_w_n)
+                {
+                    max_w_n_i = (v<<1)+1;
+                    max_w_n = w_n;
+                }
+            }
+
+            if(occ == a_n) vis[i] = 1;
+        }
+    }
+
+    for (i = 0; i < ug->g->n_seq; i++)
+    {
+        if(s[i] != 0 || (IF_BUB(i, *bub))) continue;
+        v = i;
+        e = u_trans_a(*ta, v);
+        e_n = u_trans_n(*ta, v);
+        for (m = 0, w_a = w_n = 0; m < e_n; m++)
+        {
+            if(s[e[m].tn] != 0)
+            {
+                w_a += (e[m].nw>=0?e[m].nw:-e[m].nw);
+            }
+            else
+            {
+                w_n += (e[m].nw>=0?e[m].nw:-e[m].nw);
+            }
+        }
+
+        if(w_a > max_w_a)
+        {
+            max_w_a_i = (i<<1)+1;
+            max_w_a = w_a;
+        }
+
+        if(w_n > max_w_n)
+        {
+            max_w_n_i = (i<<1)+1;
+            max_w_n = w_n;
+        }
+    }
+    
+    if(max_w_a_i != (uint32_t)-1) return max_w_a_i;
+    return max_w_n_i;
+}
+
+void init_phase(ha_ug_index* idx, kv_u_trans_t *ta, bubble_type* bub, ps_t *st)
+{
+    double index_time = yak_realtime();
+    uint8_t *vis = NULL; CALLOC(vis, idx->ug->g->n_seq);
+    bub_p_t_warp b; memset(&b, 0, sizeof(bub_p_t_warp));
+    CALLOC(b.a, idx->ug->g->n_seq*2);
+    uint32_t i, k, *a = NULL, n, beg, end;
+    memset(st->s, 0, sizeof(int8_t)*idx->ug->g->n_seq);
+    
+    double z[2];
+    u_trans_t *e = NULL;
+    uint32_t e_n;
+    while(1)
+    {
+        i = get_weightest_node(ta, bub, idx->ug, st->s, vis);
+        if(i == (uint32_t)-1) break;
+        if(i&1)
+        {
+            i>>=1;
+            z[0] = z[1] = 0;
+            e = u_trans_a(*ta, i);
+            e_n = u_trans_n(*ta, i);
+            for (k = 0; k < e_n; k++)
+            {
+                if(e[k].del) continue;
+                if(st->s[e[k].tn] > 0) z[0] += e[k].nw;
+                else if(st->s[e[k].tn] < 0) z[1] += e[k].nw;
+            }
+            if(z[0] - z[1] < 0)
+            {
+                st->s[i] = 1;
+            } 
+            else if(z[0] - z[1] > 0)
+            {
+                st->s[i] = -1;
+            }
+            else
+            {
+                st->xs = kr_splitmix64(st->xs);
+                st->s[i] = st->xs&1? 1 : -1;
+            }
+        }
+        else
+        {
+            i>>=1;
+            get_bubbles(bub, i, &beg, &end, &a, &n, NULL);
+            bub_phase(idx->ug, beg, end, &b, ta, st);
+            bub_phase(idx->ug, beg, end, &b, ta, st);
+        }
+    }
+    free(b.a); free(b.S.a); free(b.T.a); free(b.b.a); free(b.e.a);
+    free(vis);
+
+    for (i = 0; i < idx->ug->g->n_seq; i++)
+    {
+        if(st->s[i] == 0) fprintf(stderr, "ERROR\n");
+        if(u_trans_n(*ta, i) == 0) st->s[i] = 0;
+    }
+
+    fprintf(stderr, "[M::%s::%.3f]\n", __func__, yak_realtime()-index_time);
+}
+
+
+double dfs_weight_hic(uint32_t v, uint8_t* vis_flag, uint8_t* is_vis, kv_u_trans_t *ta, 
+kvec_t_u32_warp* stack, kvec_t_u32_warp* result, uint32_t e_flag, uint32_t ava_flag, uint32_t* link_occ)
+{
+    u_trans_t *e = NULL;
+    uint32_t cur, i, next = (uint32_t)-1, e_n;
+    stack->a.n = 0;
+    kv_push(uint32_t, stack->a, v);
+    double w = 0;
+    if(link_occ) (*link_occ) = 0;
+    while (stack->a.n > 0)
+    {
+        stack->a.n--;
+        cur = stack->a.a[stack->a.n];
+        if(is_vis[cur]) continue;
+        is_vis[cur] = 1;
+        if(cur!=v && vis_flag[cur] != ava_flag) continue;
+        e = u_trans_a(*ta, cur); e_n = u_trans_n(*ta, cur);
+        for (i = 0; i < e_n; i++)
+        {
+            if(e[i].del) continue;
+            next = e[i].tn;
+            if(vis_flag[next]&e_flag)
+            {
+                w += (e[i].nw >= 0? e[i].nw : -e[i].nw);
+                if(link_occ) (*link_occ) += e[i].occ;
+                continue;
+            }
+            if(is_vis[next]) continue;
+            if(vis_flag[next] != ava_flag) continue;
+            kv_push(uint32_t, stack->a, next);
+        }
+    }
+    return w;
+}
+
+double get_chain_weight_hic(bubble_type* bub, ma_ug_t *bub_ug, buf_t* b, uint32_t v, uint32_t convex_source, 
+kv_u_trans_t *ta, uint8_t* vis_flag, uint8_t* is_vis, ma_ug_t* ug, kvec_t_u32_warp* stack, 
+kvec_t_u32_warp* result, uint32_t e_flag, uint32_t ava_flag, kvec_t_u32_warp* res_utg, uint32_t* link_occ)
+{
+    long long nodeLen, baseLen, max_stop_nodeLen, max_stop_baseLen;
+    ma_utg_t *u = NULL;
+    uint32_t convex, k, k_i, k_j, *a, n, beg, sink, uID, root, cur, ncur, n_vx = ug->g->n_seq<<1, occ;
+    asg_arc_t *acur = NULL;
+    double w = 0;
+    b->b.n = 0;
+    get_unitig(bub_ug->g, NULL, v, &convex, &nodeLen, &baseLen, &max_stop_nodeLen, 
+                                                                    &max_stop_baseLen, 1, b);
+    memset(is_vis, 0, n_vx);
+
+    for (k = 0; k < b->b.n; k++)
+    {
+        u = &(bub_ug->u.a[b->b.a[k]>>1]);
+        if(u->n == 0) continue;
+
+        for (k_i = 0; k_i < u->n; k_i++)
+        {
+            get_bubbles(bub, u->a[k_i]>>33, &beg, &sink, &a, &n, NULL);
+            for(k_j = 0; k_j < n; k_j++) is_vis[a[k_j]] = is_vis[a[k_j]^1] = 1;
+            if(beg != (uint32_t)-1) is_vis[beg] = is_vis[beg^1] = 1; 
+            if(sink != (uint32_t)-1) is_vis[sink] = is_vis[sink^1] = 1; 
+        }
+    }
+
+    u = &(bub_ug->u.a[v>>1]);
+    if((v&1)==0)
+    {
+        get_bubbles(bub, (u->a[0]>>32)>>1, (((u->a[0]>>32)&1)^1)==1?&root:NULL, 
+                                    (((u->a[0]>>32)&1)^1) == 0?&root:NULL, NULL, NULL, NULL);
+    } 
+    else
+    {
+        get_bubbles(bub, (u->a[u->n-1]>>32)>>1, ((u->a[u->n-1]>>32)&1)==1?&root:NULL, 
+                                    ((u->a[u->n-1]>>32)&1) == 0?&root:NULL, NULL, NULL, NULL);
+    }
+    
+    root ^= 1;
+    ///fprintf(stderr, "root=utg%.6dl\n", (root>>1)+1);
+    is_vis[root] = 0; 
+    stack->a.n = 0;
+    kv_push(uint32_t, stack->a, root);
+    while (stack->a.n > 0)
+    {
+        stack->a.n--;
+        cur = stack->a.a[stack->a.n];
+        if(is_vis[cur]) continue;
+        is_vis[cur] = 1;
+        if(vis_flag[cur>>1] == 0) vis_flag[cur>>1] = ava_flag;///unitig not in any chain
+       
+        ncur = asg_arc_n(ug->g, cur);
+        acur = asg_arc_a(ug->g, cur);
+        for (k = 0; k < ncur; k++)
+        {
+            if(acur[k].del) continue;
+            if(is_vis[acur[k].v]) continue;
+            if(vis_flag[acur[k].v>>1] != 0 && vis_flag[acur[k].v>>1] != ava_flag) continue;
+            kv_push(uint32_t, stack->a, acur[k].v);
+        }
+    }
+
+
+    ///vis_flag keeps isloated nodes
+    uint32_t aim_0, aim_1, root_source;
+    aim_0 = root>>1;
+    u = &(bub_ug->u.a[convex_source>>1]);
+    if((convex_source&1)==1)
+    {
+        get_bubbles(bub, (u->a[0]>>32)>>1, (((u->a[0]>>32)&1)^1)==1?&root_source:NULL, 
+                                    (((u->a[0]>>32)&1)^1) == 0?&root_source:NULL, NULL, NULL, NULL);
+    } 
+    else
+    {
+        get_bubbles(bub, (u->a[u->n-1]>>32)>>1, ((u->a[u->n-1]>>32)&1)==1?&root_source:NULL, 
+                                    ((u->a[u->n-1]>>32)&1) == 0?&root_source:NULL, NULL, NULL, NULL);
+    }
+    root_source ^= 1;
+    aim_1 = root_source>>1;
+
+    ///fprintf(stderr, "aim_0=utg%.6ul, aim_1=utg%.6ul\n", aim_0+1, aim_1+1);
+
+    cur = root_source;///scan nodes that cannot be reached from root but can be reached from root_source
+    ncur = asg_arc_n(ug->g, cur);
+    acur = asg_arc_a(ug->g, cur);
+    for (k_i = 0; k_i < ncur; k_i++)
+    {
+        if(acur[k_i].del) continue;
+        if(vis_flag[acur[k_i].v>>1] != 0) continue;///skip nodes that are already reachable
+        if_conflict_utg(acur[k_i].v, &aim_0, &aim_1, ug, vis_flag, is_vis, ava_flag, stack);
+    }
+
+
+    for (k = 0; k < ug->g->n_seq; k++)
+    {
+        if(vis_flag[k] == ava_flag)
+        {
+            cur = k<<1;
+            ncur = asg_arc_n(ug->g, cur);
+            acur = asg_arc_a(ug->g, cur);
+            for (k_i = 0; k_i < ncur; k_i++)
+            {
+                if(acur[k_i].del) continue;
+                if(vis_flag[acur[k_i].v>>1] != 0) continue;
+                if_conflict_utg(acur[k_i].v, &aim_0, &aim_1, ug, vis_flag, is_vis, ava_flag, stack);
+            }
+
+
+
+            cur = (k<<1)+1;
+            ncur = asg_arc_n(ug->g, cur);
+            acur = asg_arc_a(ug->g, cur);
+            for (k_i = 0; k_i < ncur; k_i++)
+            {
+                if(acur[k_i].del) continue;
+                if(vis_flag[acur[k_i].v>>1] != 0) continue;
+                if_conflict_utg(acur[k_i].v, &aim_0, &aim_1, ug, vis_flag, is_vis, ava_flag, stack);
+            }
+        }
+    }
+    
+
+
+    memset(is_vis, 0, n_vx);
+    if(link_occ) (*link_occ) = 0;
+    for (k = result->a.n = 0, w = 0; k < b->b.n; k++)
+    {
+        u = &(bub_ug->u.a[b->b.a[k]>>1]);
+        if(u->n == 0) continue;
+        for (k_i = 0; k_i < u->n; k_i++)
+        {
+            get_bubbles(bub, u->a[k_i]>>33, NULL, NULL, &a, &n, NULL);
+
+            for (k_j = 0; k_j < n; k_j++)
+            {
+                uID = a[k_j]>>1;
+                w += dfs_weight_hic(uID, vis_flag, is_vis, ta, stack, result, e_flag, ava_flag, &occ);
+                if(link_occ) (*link_occ) += occ;
+            }
+        }
+    }
+
+
+    for (k = 0; k < ug->g->n_seq; k++)
+    {
+        if(vis_flag[k] == ava_flag)
+        {
+            vis_flag[k] = 0;
+            if(res_utg && (!IF_HOM(k, *bub)))
+            {
+                kv_push(uint32_t, res_utg->a, k<<1);
+            }
+        }
+    }
+
+    return w;
+}
+
+
+void clean_bubble_chain_by_hic(ma_ug_t* ug, kv_u_trans_t *ta, bubble_type* bub)
+{   
+    double index_time = yak_realtime();
+    ma_ug_t *bs_ug = bub->b_ug;
+    uint32_t v, u, i, m, max_i, nv, rv, n_vx, root, flag_pri = 1, flag_aux = 2, flag_ava = 4, occ;
+    double w, cutoff = 2;
+    uint32_t max_w_occ = 4;
+    asg_arc_t *av = NULL;
+    n_vx = bs_ug->g->n_seq << 1;
+    uint8_t *vis = NULL; CALLOC(vis, ug->g->n_seq<<1);
+    uint8_t *is_vis = NULL; CALLOC(is_vis, ug->g->n_seq<<1);
+    uint8_t *is_used = NULL; CALLOC(is_used, n_vx);
+    uint8_t *dedup = NULL; CALLOC(dedup, ug->g->n_seq<<1);
+    buf_t b; memset(&b, 0, sizeof(buf_t));
+    kvec_t_u32_warp stack, result, res_utg;
+    kv_init(stack.a); kv_init(result.a); kv_init(res_utg.a);
+    double *e_w = NULL; MALLOC(e_w, bs_ug->g->n_arc);
+    uint32_t *e_occ = NULL, *a_occ = NULL; CALLOC(e_occ, bs_ug->g->n_arc);
+    double *aw = NULL, max_w = 0;
+    kvec_asg_arc_t_warp edges; kv_init(edges.a);
+    ma_ug_t *back_bs_ug = copy_untig_graph(bs_ug);
+
+    for (i = 0; i < bs_ug->g->n_arc; i++)///weight of bs_ug's edges
+    {
+        e_w[i] = -1;
+    }
+
+    for (i = 0; i < bs_ug->g->n_seq; i++)///init all chain with flag_aux 
+    {
+        set_b_utg_weight_flag(bub, &b, i<<1, vis, flag_aux, NULL);
+    }
+
+
+    for (v = 0; v < n_vx; v++)
+    {
+        av = asg_arc_a(bs_ug->g, v);
+        nv = asg_arc_n(bs_ug->g, v);
+        aw = (&e_w[bs_ug->g->idx[v]>>32]);
+        a_occ = (&e_occ[bs_ug->g->idx[v]>>32]);
+        if(nv <= 1 || get_real_length(bs_ug->g, v, NULL) <= 1) continue;
+        set_b_utg_weight_flag_xor(bub, bs_ug, &b, v^1, vis, flag_pri, NULL);
+
+        for (i = 0; i < nv; i++)
+        {
+            if(av[i].del) continue;
+            w = get_chain_weight_hic(bub, bs_ug, &b, av[i].v, v, ta, vis, is_vis, ug, &stack, &result, flag_pri, flag_ava, NULL, &occ);
+            aw[i] = w;
+            a_occ[i] = occ;
+        }
+
+        set_b_utg_weight_flag_xor(bub, bs_ug, &b, v^1, vis, flag_pri, NULL);
+    }
+
+
+    for (v = 0; v < n_vx; v++)
+    {
+        av = asg_arc_a(bs_ug->g, v);
+        nv = asg_arc_n(bs_ug->g, v);
+        aw = (&e_w[bs_ug->g->idx[v]>>32]);
+        a_occ = (&e_occ[bs_ug->g->idx[v]>>32]);
+        if(nv <= 1 || get_real_length(bs_ug->g, v, NULL) <= 1) continue;
+
+        for (i = rv = 0, max_i = (uint32_t)-1; i < nv; i++)
+        {
+            if(av[i].del) continue;
+            if(max_i == (uint32_t)-1)
+            {
+                max_i = i;
+                max_w = aw[i];
+            }
+            else if(max_w < aw[i])
+            {
+                max_i = i;
+                max_w = aw[i];
+            }
+            rv++;
+        }
+
+        if(max_i == (uint32_t)-1) continue;
+        ///if(max_w <= max_w_cutoff) continue; //must be <=
+        if(a_occ[max_i] <= max_w_occ) continue; //must be <=
+        if(rv < 2) continue;
+
+        for (i = 0; i < nv; i++)
+        {
+            if(av[i].del) continue;
+            if(i == max_i) continue;
+            ///if((av[i].v>>1) == (v>>1) && aw[i] <= max_w_cutoff) continue; ///might be not reasonable
+            if((av[i].v>>1) == (v>>1) && a_occ[i] <= max_w_occ) continue; ///might be not reasonable
+            if(aw[i]*cutoff < max_w && double_check_bub_branch(&av[i], bs_ug, e_w, e_occ, cutoff, max_w_occ))
+            {
+                av[i].del = 1; asg_arc_del(bs_ug->g, (av[i].v)^1, (av[i].ul>>32)^1, 1);
+            }
+        }
+    }
+
+    uint32_t rId_0, ori_0, rId_1, ori_1, root_0, root_1, new_bub;
+    if(bub->num.n > 0) bub->num.n--;
+    new_bub = bub->b_g->n_seq;
+
+    for (v = 0; v < n_vx; v++)
+    {
+        av = asg_arc_a(bs_ug->g, v);
+        nv = asg_arc_n(bs_ug->g, v);
+        rv = get_real_length(bs_ug->g, v, NULL);
+        if(nv == rv) continue;
+        if(rv != 1 || nv <= 1) continue;
+        get_real_length(bs_ug->g, v, &u);
+        u ^= 1;
+        if(get_real_length(bs_ug->g, u, NULL) != 1) continue;
+        drop_g_edges_by_utg(bub, bub->b_g, bs_ug, NULL, v, u);
+        if(is_used[v] || is_used[u]) continue;
+
+        is_used[v] = is_used[u] = 1;
+        root = get_utg_end_from_btg(bub, bs_ug, v);
+        rId_0 = root>>1;
+        ori_0 = root&1;
+        get_bubbles(bub, rId_0, ori_0 == 1?&root_0:NULL, ori_0 == 0?&root_0:NULL, NULL, NULL, NULL);
+
+        root = get_utg_end_from_btg(bub, bs_ug, u);
+        rId_1 = root>>1;
+        ori_1 = root&1;
+        get_bubbles(bub, rId_1, ori_1 == 1?&root_1:NULL, ori_1 == 0?&root_1:NULL, NULL, NULL, NULL);
+
+        res_utg.a.n = 0;
+        set_b_utg_weight_flag_xor(bub, back_bs_ug, &b, v^1, vis, flag_pri, NULL);
+        get_chain_weight_hic(bub, back_bs_ug, &b, u^1, v, ta, vis, is_vis, ug, &stack, &result, flag_pri, flag_ava, &res_utg, NULL);
+        set_b_utg_weight_flag_xor(bub, back_bs_ug, &b, v^1, vis, flag_pri, NULL);
+        for (i = 0; i < res_utg.a.n; i++) dedup[res_utg.a.a[i]>>1] |= 1;
+
+        set_b_utg_weight_flag_xor(bub, back_bs_ug, &b, u^1, vis, flag_pri, NULL);
+        get_chain_weight_hic(bub, back_bs_ug, &b, v^1, u, ta, vis, is_vis, ug, &stack, &result, flag_pri, flag_ava, &res_utg, NULL);
+        set_b_utg_weight_flag_xor(bub, back_bs_ug, &b, u^1, vis, flag_pri, NULL);
+        for (; i < res_utg.a.n; i++) dedup[res_utg.a.a[i]>>1] |= 2;
+
+        
+        for (i = m = 0; i < res_utg.a.n; i++)
+        {
+            if(dedup[res_utg.a.a[i]>>1] == 3)
+            {
+                res_utg.a.a[m] = res_utg.a.a[i];
+                m++;
+            }
+            dedup[res_utg.a.a[i]>>1] = 0;
+        }
+        res_utg.a.n = m;
+
+        // fprintf(stderr, "res_utg.a.n: %u, m: %u, beg-utg%.6ul, sink-utg%.6ul\n", 
+        //             res_utg.a.n, m, (root_0>>1)+1, (root_1>>1)+1);
+        
+        if(!IF_HOM(root_0>>1, *bub)) kv_push(uint32_t, res_utg.a, root_0);
+        if(!IF_HOM(root_1>>1, *bub)) kv_push(uint32_t, res_utg.a, root_1);
+
+        update_bubble_graph(&res_utg, root_0^1, rId_0, root_1^1, rId_1, bub, &edges, bub->b_g, NULL, NULL, ug, NULL, 0);
+
+        ///fprintf(stderr, "\n******src-btg%.6ul------>dest-btg%.6ul\n", (v>>1)+1, (u>>1)+1);
+    }
+
+    kv_push(uint32_t, bub->num, bub->list.n);
+    new_bub = bub->b_g->n_seq - new_bub;
+    bub->cross_bub += new_bub;
+    ///actually not useful, and may have bug when one bubble at multipe chains
+    if(new_bub) update_bub_b_s_idx(bub);
+    
+    ///debug_tangle_bubble(bub, bub->b_g->n_seq - bub->cross_bub, bub->b_g->n_seq - 1, "Cross-tangle");
+    update_bsg(bub->b_g, &edges);
+
+    ma_ug_destroy(bs_ug);
+    bs_ug = ma_ug_gen(bub->b_g);
+    bub->b_ug = bs_ug;
+    kv_destroy(bub->chain_weight);
+    ma_utg_t *u_x = NULL;
+    bs_ug = bub->b_ug;
+    kv_malloc(bub->chain_weight, bs_ug->u.n); bub->chain_weight.n = bs_ug->u.n;
+    for (i = 0; i < bs_ug->u.n; i++)
+    {
+        u_x = &(bs_ug->u.a[i]);
+        bub->chain_weight.a[i].id = i;
+        // if(u->n <= 1) ///not a chain
+        // {
+        //     bub->chain_weight.a[i].b_occ = bub->chain_weight.a[i].g_occ = 0;
+        //     bub->chain_weight.a[i].del = 1;
+        // } 
+        // else
+        {
+            bub->chain_weight.a[i].del = 0;
+            calculate_chain_weight(u_x, bub, ug, &(bub->chain_weight.a[i]));
+        }
+    }
+    qsort(bub->chain_weight.a, bub->chain_weight.n, sizeof(chain_w_type), cmp_chain_weight);
+
+    
+    
+    free(vis); free(is_vis); free(is_used); free(dedup); free(b.b.a); free(e_w); free(e_occ);
+    kv_destroy(stack.a); kv_destroy(result.a); kv_destroy(res_utg.a); kv_destroy(edges.a);
+    ma_ug_destroy(back_bs_ug);
+    fprintf(stderr, "[M::%s::%.3f]\n", __func__, yak_realtime()-index_time);
+}
+
+void append_boundary_chain_hic(ma_ug_t* ug, kv_u_trans_t *ta, bubble_type* bub)
+{
+    double index_time = yak_realtime();
+    ma_ug_t *bs_ug = bub->b_ug;
+    uint32_t v, u, i, k, beg_idx, m, nv, n_vx, flag_pri = 1, flag_aux = 2, flag_ava = 4;
+    uint32_t root, rId_0, ori_0, root_0, new_bub;
+    asg_arc_t *av = NULL;
+    n_vx = bs_ug->g->n_seq << 1;
+    uint8_t *vis = NULL; CALLOC(vis, ug->g->n_seq<<1);
+    uint8_t *is_vis = NULL; CALLOC(is_vis, ug->g->n_seq<<1);
+    uint8_t *is_used = NULL; CALLOC(is_used, n_vx);
+    uint8_t *dedup = NULL; CALLOC(dedup, ug->g->n_seq<<1);
+    buf_t b; memset(&b, 0, sizeof(buf_t));
+    kvec_t_u32_warp stack, result, res_utg;
+    kv_init(stack.a); kv_init(result.a); kv_init(res_utg.a);
+    kvec_asg_arc_t_warp edges; kv_init(edges.a);
+
+    for (i = 0; i < bs_ug->g->n_seq; i++)
+    {
+        set_b_utg_weight_flag(bub, &b, i<<1, vis, flag_aux, NULL);
+    }
+
+    if(bub->num.n > 0) bub->num.n--;
+    new_bub = bub->b_g->n_seq;
+    for (v = 0; v < n_vx; v++)
+    {
+        av = asg_arc_a(bs_ug->g, v);
+        nv = asg_arc_n(bs_ug->g, v);
+        if(nv == 0 || get_real_length(bs_ug->g, v, NULL) == 0) continue;
+        res_utg.a.n = 0;
+        for (i = 0; i < nv; i++)
+        {
+            if(av[i].del) continue;
+            u = av[i].v^1;
+            beg_idx = res_utg.a.n;
+            set_b_utg_weight_flag_xor(bub, bs_ug, &b, u^1, vis, flag_pri, NULL);
+            get_chain_weight_hic(bub, bs_ug, &b, v^1, u, ta, vis, is_vis, ug, &stack, &result, flag_pri, flag_ava, &res_utg, NULL);
+            set_b_utg_weight_flag_xor(bub, bs_ug, &b, u^1, vis, flag_pri, NULL);
+            for (k = m = beg_idx; k < res_utg.a.n; k++)
+            {
+                if(dedup[res_utg.a.a[k]>>1] != 0) continue;
+                dedup[res_utg.a.a[k]>>1] = 1;
+                res_utg.a.a[m] = res_utg.a.a[k];
+                m++;
+            }
+            res_utg.a.n = m;
+        }
+
+        for (k = 0; k < res_utg.a.n; k++) dedup[res_utg.a.a[k]>>1] = 0;
+
+        root = get_utg_end_from_btg(bub, bs_ug, v);
+        rId_0 = root>>1;
+        ori_0 = root&1;
+        get_bubbles(bub, rId_0, ori_0 == 1?&root_0:NULL, ori_0 == 0?&root_0:NULL, NULL, NULL, NULL);
+
+        if(root_0 != (uint32_t)-1 && (!IF_HOM(root_0>>1, *bub))) kv_push(uint32_t, res_utg.a, root_0);
+
+        if(v&1)
+        {
+            update_bubble_graph(&res_utg, root_0^1, rId_0, (uint32_t)-1, (uint32_t)-1, bub, &edges, bub->b_g, NULL, NULL, ug, NULL, 0);
+        }
+        else
+        {
+            update_bubble_graph(&res_utg, (uint32_t)-1, (uint32_t)-1, root_0^1, rId_0, bub, &edges, bub->b_g, NULL, NULL, ug, NULL, 0);
+        }
+    }
+    kv_push(uint32_t, bub->num, bub->list.n);
+    new_bub = bub->b_g->n_seq - new_bub;
+    bub->mess_bub += new_bub;
+    ///actually not useful, and may have bug when one bubble at multipe chains
+    if(new_bub) update_bub_b_s_idx(bub);
+
+
+    for (v = 0; v < n_vx; v++)
+    {
+        av = asg_arc_a(bs_ug->g, v);
+        nv = asg_arc_n(bs_ug->g, v);
+        if(nv == 0 || get_real_length(bs_ug->g, v, NULL) == 0) continue;
+        drop_g_edges_by_utg(bub, bub->b_g, bs_ug, NULL, v, (uint32_t)-1);
+    }
+
+    update_bsg(bub->b_g, &edges);
+
+    ma_ug_destroy(bs_ug);
+    bs_ug = ma_ug_gen(bub->b_g);
+    bub->b_ug = bs_ug;
+    kv_destroy(bub->chain_weight);
+    ma_utg_t *u_x = NULL;
+    bs_ug = bub->b_ug;
+    kv_malloc(bub->chain_weight, bs_ug->u.n); bub->chain_weight.n = bs_ug->u.n;
+    for (i = 0; i < bs_ug->u.n; i++)
+    {
+        u_x = &(bs_ug->u.a[i]);
+        bub->chain_weight.a[i].id = i;
+        // if(u->n <= 1) ///not a chain
+        // {
+        //     bub->chain_weight.a[i].b_occ = bub->chain_weight.a[i].g_occ = 0;
+        //     bub->chain_weight.a[i].del = 1;
+        // } 
+        // else
+        {
+            bub->chain_weight.a[i].del = 0;
+            calculate_chain_weight(u_x, bub, ug, &(bub->chain_weight.a[i]));
+        }
+    }
+    qsort(bub->chain_weight.a, bub->chain_weight.n, sizeof(chain_w_type), cmp_chain_weight);
+
+
+
+    free(vis); free(is_vis); free(is_used); free(dedup); free(b.b.a);
+    kv_destroy(stack.a); kv_destroy(result.a); kv_destroy(res_utg.a);
+    kv_destroy(edges.a);
+    fprintf(stderr, "[M::%s::%.3f]\n", __func__, yak_realtime()-index_time);
+}
+
+double get_specific_hic_weight_by_chain(uint32_t uid, kv_u_trans_t *ta, uint8_t* vis, uint8_t flag)
+{
+    u_trans_t *e = u_trans_a(*ta, uid);
+    uint64_t k, e_n = u_trans_n(*ta, uid);
+    double w = 0;
+    
+    for (k = 0; k < e_n; k++)
+    {
+        if(e[k].del) continue;
+        if(vis[e[k].tn] != flag) continue;
+        w += (e[k].nw>=0?e[k].nw:-e[k].nw);
+    }
+    return w;
+}
+
+void update_bubble_weight(bub_sort_vec* w_stack, uint32_t idx, kv_u_trans_t *ta, uint8_t* vis, uint32_t flag_cur)
+{
+    w_stack->a[idx].used = 1;
+    uint32_t k, m, uid = w_stack->a[idx].p_id>>1;
+    u_trans_t *e = u_trans_a(*ta, uid);
+    uint32_t e_n = u_trans_n(*ta, uid);
+    for (k = 0; k < e_n; k++)
+    {
+        if(e[k].del) continue;
+        if(vis[e[k].tn] != flag_cur) continue;
+        for (m = 0; m < w_stack->n; m++)
+        {
+            if(e[k].tn == (w_stack->a[m].p_id>>1)) break;
+        }
+        if(m >= w_stack->n) continue;
+        if(w_stack->a[m].used) continue;
+        w_stack->a[m].weight += (e[k].nw>=0?e[k].nw:-e[k].nw);
+    }
+}
+
+void reorder_bubbble_chain(kv_u_trans_t *ta, bubble_type* bub, bub_sort_vec* w_stack, 
+uint8_t* vis, uint32_t n_utg, uint32_t chain_id)
+{
+    uint32_t i, k, m, max_idx, flag_cur = 3, flag_right = 2, flag_left = 1, flag_unset = 0, *a, n;
+    uint64_t bid, uid;
+    ma_utg_t *u = &(bub->b_ug->u.a[chain_id]);
+    memset(vis, flag_unset, n_utg);
+    for (i = 0; i < u->n; i++)
+    {
+        bid = u->a[i]>>33;
+        get_bubbles(bub, bid, NULL, NULL, &a, &n, NULL);
+        for (k = 0; k < n; k++)
+        {
+            uid = a[k]>>1;
+            vis[uid] = flag_right;
+        }
+    }
+
+    for (i = 0; i < u->n; i++)
+    {
+        w_stack->n = 0;
+        bid = u->a[i]>>33;
+        get_bubbles(bub, bid, NULL, NULL, &a, &n, NULL);
+        kv_resize(bub_sort_type, *w_stack, n); 
+        w_stack->n = n;
+        for (k = 0; k < n; k++)
+        {
+            uid = a[k]>>1;
+            vis[uid] = flag_cur;
+            w_stack->a[k].p_id = a[k];
+            w_stack->a[k].weight = 0;
+            w_stack->a[k].used = 0;
+        }
+        
+        for (k = 0; k < w_stack->n; k++)
+        {
+            w_stack->a[k].weight += get_specific_hic_weight_by_chain(w_stack->a[k].p_id>>1, 
+            ta, vis, flag_left);
+            w_stack->a[k].weight -= get_specific_hic_weight_by_chain(w_stack->a[k].p_id>>1, 
+            ta, vis, flag_right);
+        }
+        m = 0;
+        while ((max_idx = get_max_hap_g(w_stack, NULL)) != (uint32_t)-1)
+        {
+            a[m] = w_stack->a[max_idx].p_id;
+            m++;
+            update_bubble_weight(w_stack, max_idx, ta, vis, flag_cur);
+        }
+
+        while ((max_idx = get_max_hap_g(w_stack, &max_idx)) != (uint32_t)-1)
+        {
+            a[m] = w_stack->a[max_idx].p_id;
+            m++;
+            update_bubble_weight(w_stack, max_idx, ta, vis, flag_cur);
+        }
+        if(m != n) fprintf(stderr, "ERROR\n");
+
+
+        for (k = 0; k < n; k++)
+        {
+            uid = a[k]>>1;
+            vis[uid] = flag_left;
+        }
+    }
+}
+
+void reorder_bubbles(bubble_type* bub, kv_u_trans_t *ta, uint32_t n_utg)
+{
+    double index_time = yak_realtime();
+    uint8_t* vis = NULL; MALLOC(vis, n_utg);
+    bub_sort_vec w_stack; kv_init(w_stack);
+    uint32_t i;
+
+    for (i = 0; i < bub->chain_weight.n; i++)
+    {
+        if(bub->chain_weight.a[i].del) continue;
+        reorder_bubbble_chain(ta, bub, &w_stack, vis, n_utg, bub->chain_weight.a[i].id);
+    }
+
+    free(vis); kv_destroy(w_stack);
+    fprintf(stderr, "[M::%s::%.3f]\n", __func__, yak_realtime()-index_time);
+}
+
+void update_trans_g(ha_ug_index* idx, kv_u_trans_t *ta, bubble_type* bub)
+{
+    double index_time = yak_realtime();
+    update_bubble_chain(idx->ug, bub, 0, 1);
+
+    // print_debug_bubble_graph(bub, idx->ug, "bub-0");
+
+    resolve_bubble_chain_tangle(idx->ug, bub);
+
+    // print_debug_bubble_graph(bub, idx->ug, "bub-1");
+
+    clean_bubble_chain_by_hic(idx->ug, ta, bub);
+
+    // print_debug_bubble_graph(bub, idx->ug, "bub-2");
+
+    append_boundary_chain_hic(idx->ug, ta, bub);
+
+    reorder_bubbles(bub, ta, idx->ug->g->n_seq);
+    fprintf(stderr, "[M::%s::%.3f]\n", __func__, yak_realtime()-index_time);
+}
 
 uint32_t get_max_unitig(H_partition* h, G_partition* g_p, hc_links* link, bubble_type* bub)
 {
@@ -13409,13 +14484,19 @@ void debug_gfa_space(ma_ug_t* ug, trans_chain* t_ch)
     destory_hc_links(&link);
 }
 
-void idx_hc_links(kvec_pe_hit* hits, ha_ug_index* idx)
+void idx_hc_links(kvec_pe_hit* hits, ha_ug_index* idx, bubble_type* bub)
 {
     uint64_t k, l;
+    uint32_t qn, tn;
+
     kv_resize(uint64_t, hits->idx, idx->ug->g->n_seq);
     hits->idx.n = idx->ug->g->n_seq;
     memset(hits->idx.a, 0, hits->idx.n*sizeof(uint64_t));
 
+    kv_resize(uint64_t, hits->occ, idx->ug->g->n_seq);
+    hits->occ.n = idx->ug->g->n_seq;
+    memset(hits->occ.a, 0, hits->occ.n*sizeof(uint64_t));
+
     radix_sort_pe_hit_idx_an1(hits->a.a, hits->a.a + hits->a.n);
     for (k = 1, l = 0; k <= hits->a.n; ++k) 
     {   
@@ -13426,6 +14507,17 @@ void idx_hc_links(kvec_pe_hit* hits, ha_ug_index* idx)
             
             hits->idx.a[((hits->a.a[l].s<<1)>>(64 - idx->uID_bits))] 
                                                 = (uint64_t)l << 32 | (k - l);
+
+            for (; l < k; l++)
+            {
+                qn = ((hits->a.a[l].s<<1)>>(64 - idx->uID_bits));
+                tn = ((hits->a.a[l].e<<1)>>(64 - idx->uID_bits));
+                if(IF_HOM(qn, *bub)) continue;
+                if(IF_HOM(tn, *bub)) continue;
+                hits->occ.a[qn]++;
+                hits->occ.a[tn]++;
+            }
+            
             l = k;
         }
     }
@@ -13451,7 +14543,7 @@ kv_u_trans_t *ta, trans_idx* dis)
             kv_pushp(u_trans_t, *ta, &p);
             memset(p, 0, sizeof(u_trans_t));
             p->qn = i; p->tn = link->a.a[i].e.a[k].uID;
-            p->nw = 0; 
+            p->nw = 0; p->occ = 0;
         }
     }
     kt_u_trans_t_idx(ta, idx->ug->g->n_seq);
@@ -13479,8 +14571,8 @@ kv_u_trans_t *ta, trans_idx* dis)
             weight = get_trans_weight_advance(idx, t_d, dis);
         }
         
-        e1->nw -= weight;
-        e2->nw -= weight;
+        e1->nw -= weight; e1->occ++;
+        e2->nw -= weight; e2->occ++;
     }    
 }
 
@@ -13567,27 +14659,9 @@ pe_hit *hits, uint32_t occ, uint32_t qid, uint32_t qs, uint32_t qe, uint32_t tid
 
         w += weight;
     }
-
-    // for (k = 0, w = 0; k < occ; k++)///all hits of qid
-    // {   
-    //     interpr_hit(idx, hits[k].s, hits[k].len>>32, &s_uid, &s_beg, &s_end);
-    //     if(s_uid != qid) continue;
-    //     if(!(qs <= s_beg && qe >= s_end)) continue;
-
-    //     interpr_hit(idx, hits[k].e, (uint32_t)hits[k].len, &e_uid, &e_beg, &e_end);
-    //     if(e_uid != tid) continue;
-    //     if(!(ts <= e_beg && te >= e_end)) continue;
-
-    //     t_d = get_hic_distance(&hits[k], link, idx);
-    //     if(t_d == (uint64_t)-1) continue;
-
-    //     weight = 1;
-    //     if(dis) weight = get_trans_weight_advance(idx, t_d, dis);
-
-    //     w += weight;
-    // }
     return w;
 }
+
 double get_hits_weight(ha_ug_index* idx, kvec_pe_hit* hits, hc_links* link, trans_idx* dis, 
 u_trans_t *t_a, uint32_t t_n, uint32_t qid, kv_u_trans_t *ta_idx)
 {
@@ -13597,7 +14671,8 @@ u_trans_t *t_a, uint32_t t_n, uint32_t qid, kv_u_trans_t *ta_idx)
     /****************************may have bugs********************************/
     uint32_t k, i, x_n, y;
     u_trans_t *x_a = NULL;
-    double w;
+    double occ_q, occ_t;
+    double w, i_w;
     for (k = 0, w = 0; k < t_n; k++)
     {
         /****************************may have bugs********************************/
@@ -13617,16 +14692,24 @@ u_trans_t *t_a, uint32_t t_n, uint32_t qid, kv_u_trans_t *ta_idx)
             if(i >= x_n) continue; ///no hit bridging tn and qid
         }
         /****************************may have bugs********************************/
+        occ_q = hits->occ.a[qid];
+        occ_t = hits->occ.a[t_a[k].tn] * ((double)(t_a[k].te - t_a[k].ts) / (double)(idx->ug->g->seq[t_a[k].tn].len));
+        if(occ_t < 1) occ_t = 1;
+
         ///q--->t hits
-        w += get_interval_weight(idx, link, dis, hits->a.a + (hits->idx.a[qid]>>32),
+        i_w = get_interval_weight(idx, link, dis, hits->a.a + (hits->idx.a[qid]>>32),
         (uint32_t)(hits->idx.a[qid]), qid, 0, idx->ug->g->seq[qid].len, t_a[k].tn, 
         t_a[k].ts, t_a[k].te);
+        if(i_w != 0) i_w /= (MIN(occ_q, occ_t));
+        w += i_w;
         if(t_a[k].tn == qid) continue;
 
         ///t--->q hits
-        w += get_interval_weight(idx, link, dis, hits->a.a + (hits->idx.a[t_a[k].tn]>>32),
+        i_w = get_interval_weight(idx, link, dis, hits->a.a + (hits->idx.a[t_a[k].tn]>>32),
         (uint32_t)(hits->idx.a[t_a[k].tn]), t_a[k].tn, t_a[k].ts, t_a[k].te, 
-        qid, 0, idx->ug->g->seq[qid].len);        
+        qid, 0, idx->ug->g->seq[qid].len);    
+        if(i_w != 0) i_w /= (MIN(occ_q, occ_t));
+        w += i_w;
     }
     return w;
 }
@@ -13638,27 +14721,23 @@ kv_u_trans_t *ta, kv_u_trans_t *ref, trans_idx* dis)
     uint32_t n, k, i, m, qn, tn;
     uint8_t *vis = NULL; CALLOC(vis, idx->ug->g->n_seq);
     double w;
-    // for (k = 0; k < ta->idx.n; k++)
-    // {
-    //     a = u_trans_a(*ta, k);
-    //     n = u_trans_n(*ta, k);
-    //     ///for each pair qn, tn
-    //     ///count hic pairs between (qn, tn^) and (qn^1, tn)
-    //     for (i = 0; i < n; i++)
-    //     {
-    //         if(a[i].qn == a[i].tn) continue;
-    //         if(IF_HOM(a[i].qn, *bub)) continue;
-    //         if(IF_HOM(a[i].tn, *bub)) continue;
-    //         ///(qn, tn^)
-    //         a[i].nw += get_hits_weight(idx, hits, link, dis, 
-    //                                 u_trans_a(*ref, a[i].tn), u_trans_n(*ref, a[i].tn), a[i].qn);
-    //         ///(qn^1, tn)
-    //         a[i].nw += get_hits_weight(idx, hits, link, dis, 
-    //                                 u_trans_a(*ref, a[i].qn), u_trans_n(*ref, a[i].qn), a[i].tn);
-    //     }
-    // }
+
+    for (i = m = 0; i < ta->n; i++)
+    {
+        if(ta->a[i].del) continue;
+        if(IF_HOM(ta->a[i].qn, *bub)) continue;
+        if(IF_HOM(ta->a[i].tn, *bub)) continue;
+        ta->a[m] = ta->a[i];
+        if(ta->a[m].nw != 0)
+        {
+            ta->a[m].nw /= (double)(MIN(hits->occ.a[ta->a[m].qn], hits->occ.a[ta->a[m].tn]));
+        }
+        
+        m++;
+    }
+    ta->n = m;
     
-    fprintf(stderr, "+++++ta->n=%u\n", (uint32_t)ta->n);
+    // fprintf(stderr, "+++++ta->n=%u\n", (uint32_t)ta->n);
     for (k = 0; k < ta->idx.n; k++)///all nodes
     {
         if(IF_HOM(k, *bub)) continue;
@@ -13680,15 +14759,14 @@ kv_u_trans_t *ta, kv_u_trans_t *ref, trans_idx* dis)
                                     u_trans_a(*ref, a[i].qn), u_trans_n(*ref, a[i].qn), a[i].tn, ta);
         }
         
+        
         for (i = 0; i < ta->idx.n; i++)///all edges
         {
-            // fprintf(stderr, "+a+k=%u, i=%u, ta->idx.n=%u\n", k, i, (uint32_t)ta->idx.n);
             qn = k; tn = i; w = 0;
             if(IF_HOM(qn, *bub)) continue;
             if(IF_HOM(tn, *bub)) continue;
             if(vis[tn]) continue;
             if(qn == tn) continue;
-            // fprintf(stderr, "+b+k=%u, i=%u, ta->idx.n=%u\n", k, i, (uint32_t)ta->idx.n);
             ///(qn, tn^)
             w += get_hits_weight(idx, hits, link, dis, 
                                 u_trans_a(*ref, tn), u_trans_n(*ref, tn), qn, ta);
@@ -13699,13 +14777,31 @@ kv_u_trans_t *ta, kv_u_trans_t *ref, trans_idx* dis)
             kv_pushp(u_trans_t, *ta, &p);
             memset(p, 0, sizeof(u_trans_t));///extra edges
             p->nw = w; p->qn = qn; p->tn = tn;
-            // fprintf(stderr, "+c+k=%u, i=%u, ta->idx.n=%u\n", k, i, (uint32_t)ta->idx.n);
         }
 
+        // for (i = 0; i < u_trans_n(*ref, k); i++)///only trans edges
+        // {
+        //     qn = k; tn = (u_trans_a(*ref, k))[i].tn; w = 0;
+        //     if(IF_HOM(qn, *bub)) continue;
+        //     if(IF_HOM(tn, *bub)) continue;
+        //     if(vis[tn]) continue;
+        //     if(qn == tn) continue;
+        //     ///(qn, tn^)
+        //     w += get_hits_weight(idx, hits, link, dis, 
+        //                         u_trans_a(*ref, tn), u_trans_n(*ref, tn), qn, ta);
+        //     ///(qn^1, tn)
+        //     w += get_hits_weight(idx, hits, link, dis, 
+        //                         u_trans_a(*ref, qn), u_trans_n(*ref, qn), tn, ta);
+        //     if(w == 0) continue;
+        //     kv_pushp(u_trans_t, *ta, &p);
+        //     memset(p, 0, sizeof(u_trans_t));///extra edges
+        //     p->nw = w; p->qn = qn; p->tn = tn;
+        // }
+
         for (i = 0, vis[k] = 0; i < n; i++) vis[a[i].tn] = 0;
     }
 
-    fprintf(stderr, "------ta->n=%u\n", (uint32_t)ta->n);
+    // fprintf(stderr, "------ta->n=%u\n", (uint32_t)ta->n);
 
     for (i = m = 0; i < ta->n; i++)
     {
@@ -13758,11 +14854,10 @@ ha_ug_index* idx, bubble_type* bub, int8_t *s, uint32_t ignore_dis)
         lk->a.a[i].e.n = m;
     }
 
-    if(hits->idx.n == 0) idx_hc_links(hits, idx);
+    if(hits->idx.n == 0) idx_hc_links(hits, idx, bub);
 
     weight_kv_u_trans(idx, hits, lk, bub, ta, is_comples_weight == 1? &dis : NULL);
     adjust_weight_kv_u_trans(idx, hits, lk, bub, ta, ref, is_comples_weight == 1? &dis : NULL);
-    // memset(s, 0, sizeof(int8_t)*idx->ug->g->n_seq);
     kv_destroy(dis);
 }
 
@@ -13783,6 +14878,20 @@ void print_kv_u_trans(kv_u_trans_t *ta, hc_links* lk, int8_t *s)
     
 }
 
+ps_t* init_ps_t(uint64_t seed, uint64_t n)
+{
+    ps_t *s = NULL; CALLOC(s, 1);
+    s->xs = seed;
+    CALLOC(s->s, n);
+    return s;
+}
+
+void destory_ps_t(ps_t **s)
+{
+    free((*s)->s);
+    free((*s));
+}
+
 int hic_short_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx)
 {
     double index_time = yak_realtime();
@@ -13797,7 +14906,7 @@ int hic_short_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx)
     idx->hap_cnt = asm_opt.hap_occ;
     ///int_kvec_pe_hit_hap(&sl.hits);
     ///int_kvec_pe_hit(&sl.hits);
-    kv_init(sl.hits.a); kv_init(sl.hits.idx);
+    kv_init(sl.hits.a); kv_init(sl.hits.idx); kv_init(sl.hits.occ);
 
     
     if(!load_hc_hits(&sl.hits, asm_opt.output_file_name))
@@ -13831,7 +14940,7 @@ int hic_short_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx)
     bubble_type bub; 
     kv_u_trans_t k_trans; 
     kv_init(k_trans); kv_init(k_trans.idx);
-    int8_t *s = NULL; CALLOC(s, idx->ug->g->n_seq);
+    ps_t *s = init_ps_t(11, idx->ug->g->n_seq);
     memset(&bub, 0, sizeof(bubble_type));
     bub.round_id = 0; bub.n_round = 2;
     for (bub.round_id = 0; bub.round_id < bub.n_round; bub.round_id++)
@@ -13842,11 +14951,15 @@ int hic_short_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx)
             measure_distance(idx->ug, &sl.hits, &link, &bub, &(idx->t_ch->k_trans));
             collect_hc_reverse_links(&link, idx->ug, &bub);
         }
-        renew_kv_u_trans(&k_trans, &link, &sl.hits, &(idx->t_ch->k_trans), idx, &bub, s, 0);
-        mc_solve(NULL, NULL, &k_trans, idx->ug, idx->read_g, 0.8, R_INF.trio_flag, (bub.round_id == 0? 1 : 0), s, 0);
-        fprintf(stderr, "sb-0-sb\n");
-        label_unitigs_sm(s, idx->ug);
-        fprintf(stderr, "sb-1-sb\n");
+
+        renew_kv_u_trans(&k_trans, &link, &sl.hits, &(idx->t_ch->k_trans), idx, &bub, s->s, 0);
+        // if(bub.round_id == 0) init_phase(idx, &k_trans, &bub, s); 
+        update_trans_g(idx, &k_trans, &bub);
+        /*******************************for debug************************************/
+        mc_solve(NULL, NULL, &k_trans, idx->ug, idx->read_g, 0.8, R_INF.trio_flag, 
+        (bub.round_id == 0? 1 : 0), s->s, 0, /**&bub**/NULL);
+        /*******************************for debug************************************/
+        label_unitigs_sm(s->s, idx->ug);
         /**
         init_hic_advance((ha_ug_index*)sl.idx, &sl.hits, &link, &bub, &hap, 0);
         reset_H_partition(&hap, (bub.round_id == 0? 1 : 0));
@@ -13857,7 +14970,7 @@ int hic_short_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx)
     }
 
     ///print_hc_links(&link, 0, &hap);
-    print_kv_u_trans(&k_trans, &link, s);
+    // print_kv_u_trans(&k_trans, &link, s->s);
 
     // cluster_contigs(&bub, idx, &sl.hits, &M, &hap, &link);
     // destory_MT(&M);
@@ -13880,16 +14993,16 @@ int hic_short_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx)
     
 
 
-    fprintf(stderr, "sb-2-sb\n");
     // destory_contig_partition(&hap);
     kv_destroy(back_hc_edge.a);
     ///destory_kvec_pe_hit_hap(&sl.hits);
     kv_destroy(sl.hits.a);
+    kv_destroy(sl.hits.idx);
+    kv_destroy(sl.hits.occ);
     destory_hc_links(&link);
     kv_destroy(k_trans); 
     kv_destroy(k_trans.idx);
-    free(s);
-    fprintf(stderr, "sb-3-sb\n");
+    destory_ps_t(&s);
     return 1;
     
     /*******************************for debug************************************/
diff --git a/hic.h b/hic.h
index 9eec15c..93e6a72 100644
--- a/hic.h
+++ b/hic.h
@@ -49,6 +49,11 @@ typedef struct {
     kvec_t(chain_w_type) chain_weight;
     chain_hic_warp c_w;
 } bubble_type;
+
+typedef struct {
+    int8_t *s;
+    uint64_t xs;
+} ps_t;
 #define P_het(B) ((B).num.n)
 #define M_het(B) ((B).num.n + 1)
 // #define IF_BUB(ID, B) ((B).index[(ID)] < (B).num.n)
diff --git a/rcut.cpp b/rcut.cpp
index 51c0e70..4efb56c 100644
--- a/rcut.cpp
+++ b/rcut.cpp
@@ -6,6 +6,8 @@
 #include "Purge_Dups.h"
 #include "Correct.h"
 #include "ksort.h"
+#include "kthread.h"
+#include "hic.h"
 
 #define mc_edge_key(e) ((e).x)
 KRADIX_SORT_INIT(mce, mc_edge_t, mc_edge_key, member_size(mc_edge_t, x))
@@ -17,6 +19,8 @@ KRADIX_SORT_INIT(mc64, uint64_t, mc_generic_key, 8)
 #define ma_x(z) (((z).x>>32))
 #define ma_y(z) (((uint32_t)((z).x)))
 
+uint8_t bit_filed[8] = {1, 2, 4, 8, 16, 32, 64, 128};
+#define is_bit_set(id, a) ((a)[(id)>>3]&bit_filed[(id)&7])
 
 typedef struct {
 	int32_t max_iter;
@@ -25,6 +29,7 @@ typedef struct {
     uint64_t seed;
 } mc_opt_t;
 
+
 typedef struct {
 	t_w_t z[2];
 } mc_pairsc_t;
@@ -40,14 +45,66 @@ typedef struct {
 	uint8_t *f;
 } mc_svaux_t;
 
-void mc_opt_init(mc_opt_t *opt)
+typedef struct{
+	uint64_t chain_id, bid, uid;
+}mc_bp_iter;
+
+typedef struct{
+    uint32_t chain_id;
+    uint32_t f_bid;
+    uint32_t f_uid;
+    uint32_t l_bid; 
+    uint32_t l_uid;
+    uint32_t id;
+    t_w_t w;
+}mc_bp_res;
+
+typedef struct{
+	size_t n, m;
+	uint8_t *a;
+}bits_p;
+
+typedef struct{
+    bubble_type* b_b;
+	uint64_t *idx, idx_n, occ, n_thread;
+	mc_bp_res *res;
+	uint8_t *lock;
+	mc_svaux_t *b_aux;
+	mc_match_t *ma;
+	bits_p *vis;
+}mc_bp_t;
+
+typedef struct {
+	uint64_t x; // RNG
+	uint32_t cc_off, cc_size;
+	kvec_t(uint32_t) cc_node;
+	kvec_t(uint32_t) bfs;
+	///share
+	uint32_t *bfs_mark;
+	mc_pairsc_t *z, *z_opt;///keep scores to nodes(1) and nodes(-1)
+	int8_t *s, *s_opt;
+} mc_svaux_t_s;
+
+typedef struct {
+	kvec_t(mc_svaux_t_s) bs;
+	kvec_t(uint64_t) cc_edge;
+	uint32_t *bfs_mark;
+	mc_pairsc_t *z, *z_opt;///keep scores to nodes(1) and nodes(-1)
+	int8_t *s, *s_opt;
+	uint32_t n_t, n_g;
+} mc_svaux_t_mul;
+
+
+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 = 5000;
-	opt->n_perturb = 100000;
-	opt->f_perturb = 0.1;
+	// 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 = 11;
+	opt->seed = seed;
 }
 
 mc_g_t *init_mc_g_t(ma_ug_t *ug, asg_t *read_g, int8_t *s, uint32_t renew_s)
@@ -87,22 +144,6 @@ void destory_mc_g_t(mc_g_t **p)
 }
 
 
-static inline uint64_t kr_splitmix64(uint64_t x)
-{
-	uint64_t z = (x += 0x9E3779B97F4A7C15ULL);
-	z = (z ^ (z >> 30)) * 0xBF58476D1CE4E5B9ULL;
-	z = (z ^ (z >> 27)) * 0x94D049BB133111EBULL;
-	return z ^ (z >> 31);
-}
-
-static inline double kr_drand_r(uint64_t *x)
-{
-    union { uint64_t i; double d; } u;
-	*x = kr_splitmix64(*x);
-    u.i = 0x3FFULL << 52 | (*x) >> 12;
-    return u.d - 1.0;
-}
-
 static void ks_shuffle_uint32_t(size_t n, uint32_t a[], uint64_t *x)
 {
 	size_t i, j;
@@ -567,6 +608,57 @@ void mc_g_cc(mc_match_t *ma)
 {
 	ma->cc = mc_g_cc_core(ma);
 }
+mc_bp_t *mc_bp_t_init(mc_match_t *ma, mc_svaux_t *b_aux, bubble_type* bub, uint64_t n_thread)
+{
+	uint32_t i, k, n;
+	mc_bp_t *bp = NULL; 
+	ma_utg_t *u = NULL;
+	CALLOC(bp, 1);
+	bp->b_aux = b_aux;
+	bp->ma = ma;
+	bp->b_b = bub;
+	bp->n_thread = n_thread;
+	CALLOC(bp->lock, bub->ug->g->n_seq);
+	CALLOC(bp->res, n_thread);
+	CALLOC(bp->vis, n_thread);
+	for (i = 0; i < n_thread; i++)
+	{
+		bp->vis[i].m = bp->vis[i].n = bub->ug->g->n_seq;
+		CALLOC(bp->vis[i].a, bp->vis[i].n);
+	}
+	
+
+	MALLOC(bp->idx, bub->chain_weight.n+1);
+	bp->idx_n = bp->occ = 0;
+	for (i = 0; i < bub->chain_weight.n; i++)
+	{
+		bp->idx[i] = bp->occ;
+		bp->idx_n++;
+		if(bub->chain_weight.a[i].del) continue;
+		u = &(bub->b_ug->u.a[bub->chain_weight.a[i].id]);///list of bubbles
+		for (k = 0; k < u->n; k++)
+		{
+			get_bubbles(bub, u->a[k]>>33, NULL, NULL, NULL, &n, NULL);
+			bp->occ += n;
+		}
+	}
+	bp->idx[i] = bp->occ;
+	fprintf(stderr, "# nodes in chains: %lu, # chains: %lu\n", bp->occ, bp->idx_n);
+	return bp;
+}
+
+void destroy_mc_bp_t(mc_bp_t **bp)
+{
+	uint32_t i;
+	for (i = 0; i < (*bp)->n_thread; i++)
+	{
+		free((*bp)->vis[i].a);
+	}
+	free((*bp)->idx);
+	free((*bp)->res);
+	free((*bp)->lock);
+	free((*bp));
+}
 
 mc_svaux_t *mc_svaux_init(const mc_g_t *mg, uint64_t x)
 {
@@ -584,7 +676,10 @@ mc_svaux_t *mc_svaux_init(const mc_g_t *mg, uint64_t x)
 	b->s = mg->s.a;
 	CALLOC(b->s_opt, ma->n_seq);
 	MALLOC(b->bfs, ma->n_seq);
-	MALLOC(b->bfs_mark, ma->n_seq);
+
+	MALLOC(b->bfs_mark, ma->n_seq); 
+	memset(b->bfs_mark, -1, ma->n_seq*sizeof(uint32_t));
+	
 	CALLOC(b->z, ma->n_seq);
 	CALLOC(b->z_opt, ma->n_seq);
 	CALLOC(b->f, ma->n_seq);
@@ -601,6 +696,64 @@ void mc_svaux_destroy(mc_svaux_t *b)
 	free(b);
 }
 
+
+mc_svaux_t_mul *init_mc_svaux_t_mul(const mc_g_t *mg, uint64_t n_threads)
+{
+	uint32_t st, i;
+	mc_match_t *ma = mg->e;
+	mc_svaux_t_mul *b; CALLOC(b, 1);
+	for (st = 0, i = 1, b->n_t = 0; i <= ma->n_seq; ++i)
+	{
+		if (i == ma->n_seq || ma->cc[st]>>32 != ma->cc[i]>>32)
+		{
+			b->n_t++;
+		}
+	}
+	b->n_g = b->n_t;
+	if(n_threads < b->n_t) b->n_t = n_threads;
+	
+	kv_init(b->cc_edge);
+	MALLOC(b->bfs_mark, ma->n_seq);
+	memset(b->bfs_mark, -1, ma->n_seq*sizeof(uint32_t));
+	CALLOC(b->s_opt, ma->n_seq);
+	CALLOC(b->z, ma->n_seq);
+	CALLOC(b->z_opt, ma->n_seq);
+	b->s = mg->s.a;
+
+	kv_init(b->bs); CALLOC(b->bs.a, b->n_t); 
+	b->bs.n = b->bs.m = b->n_t;
+	for (i = 0; i < b->bs.n; i++)
+	{
+		kv_init(b->bs.a[i].cc_node);
+		kv_init(b->bs.a[i].bfs);
+		b->bs.a[i].bfs_mark = b->bfs_mark;
+		b->bs.a[i].z = b->z;
+		b->bs.a[i].z_opt = b->z_opt;
+		b->bs.a[i].s = b->s;
+		b->bs.a[i].s_opt = b->s_opt;
+	}
+	return b;
+}
+
+
+void destroy_mc_svaux_t_mul(mc_svaux_t_mul *b)
+{
+	uint32_t i;
+	kv_destroy(b->cc_edge);
+	free(b->bfs_mark);
+	free(b->s_opt);
+	free(b->z);
+	free(b->z_opt);
+	b->s = NULL;
+	for (i = 0; i < b->bs.n; i++)
+	{
+		kv_destroy(b->bs.a[i].cc_node);
+		kv_destroy(b->bs.a[i].bfs);
+	}
+	kv_destroy(b->bs);
+	free(b);
+}
+
 uint32_t mc_best(const mc_match_t *ma, mc_svaux_t *b)
 {
 	uint32_t i, max_i = (uint32_t)-1;
@@ -634,6 +787,17 @@ t_w_t mc_score(const mc_match_t *ma, mc_svaux_t *b)
 	return z;
 }
 
+t_w_t mc_score_all(const mc_match_t *ma, mc_svaux_t *b)
+{
+	uint32_t k;
+	t_w_t z = 0;
+	for (k = 0; k < ma->n_seq; ++k) 
+	{
+		z += -((t_w_t)(b->s[k])) * (b->z[k].z[0] - b->z[k].z[1]);
+	}
+	return z;
+}
+
 void mc_reset_z(const mc_match_t *ma, mc_svaux_t *b)
 {
 	uint32_t i;
@@ -651,10 +815,42 @@ void mc_reset_z(const mc_match_t *ma, mc_svaux_t *b)
 	}
 }
 
+void mc_reset_z_debug(const mc_match_t *ma, mc_svaux_t *b)
+{
+	uint32_t i;
+	t_w_t z[2];
+	for (i = 0; i < b->cc_size; ++i) {
+		uint32_t k = (uint32_t)ma->cc[b->cc_off + i];///uid
+		uint32_t o = ma->idx.a[k] >> 32;
+		uint32_t j, n = (uint32_t)ma->idx.a[k];
+		z[0] = b->z[k].z[0]; z[1] = b->z[k].z[1];
+		b->z[k].z[0] = b->z[k].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 (b->s[t] > 0) b->z[k].z[0] += e->w;
+			else if (b->s[t] < 0) b->z[k].z[1] += e->w;
+		}
+		if(z[0] != b->z[k].z[0]) fprintf(stderr, "ERROR1\n");
+		if(z[1] != b->z[k].z[1]) fprintf(stderr, "ERROR2\n");
+	}
+}
+
 t_w_t mc_init_spin(const mc_opt_t *opt, const mc_match_t *ma, mc_svaux_t *b)
 {
 	uint32_t i;
 	b->cc_edge.n = 0;
+	for (i = 0; i < b->cc_size; ++i) {///how many nodes
+		uint32_t k = (uint32_t)ma->cc[b->cc_off + i];///node id
+		if(b->s[k] == 0) break;
+	}
+	if(i >= b->cc_size)
+	{
+		// fprintf(stderr, "------Set\n");
+		mc_reset_z(ma, b);
+		return mc_score(ma, b);
+	}
+	// fprintf(stderr, "++++++UnSet\n");
 	for (i = 0; i < b->cc_size; ++i) {///how many nodes
 		uint32_t k = (uint32_t)ma->cc[b->cc_off + i];///node id
 		uint32_t o = ma->idx.a[k] >> 32;///cc group id
@@ -664,7 +860,6 @@ t_w_t mc_init_spin(const mc_opt_t *opt, const mc_match_t *ma, mc_svaux_t *b)
 			w_t w = ma->ma.a[o + j].w;
 			w = w > 0? w : -w;
 			kv_push(uint64_t, b->cc_edge, (uint64_t)((uint32_t)-1 - ((uint32_t)w)) << 32 | (o + j));
-			///b->cc_edge[b->n_cc_edge++] = (uint64_t)((uint32_t)-1 - w) << 32 | (o + j);
 		}
 	}
 	radix_sort_mc64(b->cc_edge.a, b->cc_edge.a + b->cc_edge.n);
@@ -706,33 +901,23 @@ static void mc_set_spin(const mc_match_t *ma, mc_svaux_t *b, uint32_t k, int8_t
 	b->s[k] = s;
 }
 
-t_w_t mc_best_flip(const mc_match_t *ma, mc_svaux_t *b, t_w_t *sc_max)
+void mc_best_flip(const mc_match_t *ma, mc_svaux_t *b)
 {
-	uint32_t idx, k;
-	t_w_t z = 0, w = 0;
+	uint32_t idx;
 	for (idx = 0; idx < b->cc_size; ++idx) {
-		k = (uint32_t)ma->cc[b->cc_off + idx];
-		b->f[k] = 0;
-		z += -((t_w_t)(b->s[k])) * (b->z[k].z[0] - b->z[k].z[1]);
-		///b->f[(uint32_t)ma->cc[b->cc_off + idx]] = 0;
+		b->f[(uint32_t)ma->cc[b->cc_off + idx]] = 0;
         ///uint32_t k = (uint32_t)ma->cc[b->cc_off + idx];///uid
 	}
 	while (1)
 	{
 		idx = mc_best(ma, b);
 		if(idx == (uint32_t)-1) break;
-
-		w = ((t_w_t)(b->s[k])) * (b->z[k].z[0] - b->z[k].z[1]) * 4;
-		if(sc_max && (*sc_max) >= (z + w)) break;
-		z += w;
-
 		mc_set_spin(ma, b, idx, -b->s[idx]);
 		b->f[idx] = 1;
 	}
-	return z;
 }
 
-static t_w_t mc_optimize_local(const mc_opt_t *opt, const mc_match_t *ma, mc_svaux_t *b, uint32_t *n_iter, t_w_t *sc_max)
+static t_w_t mc_optimize_local(const mc_opt_t *opt, const mc_match_t *ma, mc_svaux_t *b, uint32_t *n_iter)
 {
 	uint32_t i, n_flip = 0;
 	int32_t n_iter_local = 0;
@@ -753,13 +938,7 @@ static t_w_t mc_optimize_local(const mc_opt_t *opt, const mc_match_t *ma, mc_sva
 		if (n_flip == 0) break;
 	}
 
-	if(n_flip != 0) 
-	{
-		t_w_t z_debug =  mc_best_flip(ma, b, sc_max);
-		if(z_debug != mc_score(ma, b)) fprintf(stderr, "ERROR\n");
-		return z_debug;
-	}
-
+	// if(n_flip != 0) mc_best_flip(ma, b);
 	return mc_score(ma, b);
 }
 
@@ -808,6 +987,249 @@ static void mc_perturb_node(const mc_opt_t *opt, const mc_match_t *ma, mc_svaux_
 		mc_set_spin(ma, b, b->bfs[i], -b->s[b->bfs[i]]);
 }
 
+void clean_mc_bp_res(mc_bp_res *res)
+{
+	res->chain_id = (uint32_t)-1;
+	res->f_bid = res->f_uid = res->l_bid = res->l_uid = (uint32_t)-1;
+	res->id = (uint32_t)-1; res->w = -1;
+}
+
+void reset_mc_bp_iter(mc_bp_t* bp, mc_bp_iter *x, uint64_t id)
+{
+	ma_utg_t *u = NULL;
+	uint32_t i, n, occ;
+	for (i = 0; i < bp->idx_n; i++)
+	{
+		if(id >= bp->idx[i]) break;
+	}
+
+	id -= bp->idx[i];
+	x->chain_id = bp->b_b->chain_weight.a[i].id;
+	u = &(bp->b_b->b_ug->u.a[x->chain_id]);
+	for (i = occ = 0; i < u->n; i++)
+	{
+		get_bubbles(bp->b_b, u->a[i]>>33, NULL, NULL, NULL, &n, NULL);
+		occ += n;
+		if(id < occ)
+		{
+			x->bid = i;
+			x->uid = id - (occ -n);
+			break;
+		}
+	}
+}
+inline uint32_t next_uid(mc_bp_iter *iter, bubble_type* bub, uint32_t *c_bid, uint32_t *c_uid)
+{
+	ma_utg_t *u = &(bub->b_ug->u.a[iter->chain_id]);
+	uint32_t *a, n, uid;
+	while (1)
+	{
+		if(iter->bid >= u->n) break;
+		get_bubbles(bub, u->a[iter->bid]>>33, NULL, NULL, &a, &n, NULL);
+		while (1)
+		{
+			if(iter->uid >= n) break;
+			uid = a[iter->uid]>>1;
+			if(c_bid) (*c_bid) = iter->bid;
+			if(c_uid) (*c_uid) = iter->uid;
+			iter->uid++;
+			return uid;
+		}
+		iter->bid++, iter->uid = 0;
+	}
+	return (uint32_t)-1;
+}
+
+t_w_t incre_weight(mc_svaux_t *b_aux, mc_match_t *ma, uint8_t* vis, uint32_t uid)
+{
+	mc_edge_t *o = NULL;
+	uint32_t n, i, t;
+	t_w_t w = ((t_w_t)(b_aux->s[uid])) * (b_aux->z[uid].z[0] - b_aux->z[uid].z[1]) * 2;
+	t_w_t w_off = 0;
+	o = pt_a(*ma, uid);
+	n = pt_n(*ma, uid);
+	for (i = 0; i < n; ++i) 
+	{
+		t = ma_y(o[i]);
+		if(vis[t] == 0) continue;
+        if(t == uid) continue;
+		w_off += (b_aux->s[uid]*b_aux->s[t]*o[i].w);
+	}
+	return w - (w_off*4);//2 for self; 4 for both directions
+}
+
+void select_min_bp(bubble_type* bub, mc_match_t *ma, mc_svaux_t *b_aux, mc_bp_t* bp,
+uint8_t *lock, bits_p *vis, uint64_t id, mc_bp_res* r)
+{
+	uint32_t uid, val = 0, max_bid, max_uid, c_bid, c_uid, f_bid, f_uid;
+	t_w_t w = 0, max_w = -1;
+	mc_bp_iter i;
+	reset_mc_bp_iter(bp, &i, id);
+	memset(vis->a, 0, vis->n);
+	max_bid = max_uid = (uint32_t)-1;
+	f_bid = i.bid; f_uid = i.uid;
+
+	while (1)
+	{
+		uid = next_uid(&i, bub, &c_bid, &c_uid);
+		if(uid == (uint32_t)-1) break;
+		if(vis->a[uid]) continue; ///already flip uid
+		///update w
+		w += incre_weight(b_aux, ma, vis->a, uid);
+		vis->a[uid] = 1;
+		if(lock[uid] == 0) val = 1;
+        if(val == 0) continue;
+		///update max_w
+		if(max_w < w)
+		{
+			max_w = w;
+			max_bid = c_bid;
+			max_uid = c_uid;
+		}
+	}
+
+	if(max_w <= 0 || max_bid == (uint32_t)-1 || max_uid == (uint32_t)-1) return;
+	if((max_w > r->w) || (max_w == r->w && id < r->id))
+	{
+		r->w = max_w;
+		r->id = id;
+		r->chain_id = i.chain_id;
+		r->l_bid = max_bid;
+		r->l_uid = max_uid;
+		r->f_bid = f_bid;
+		r->f_uid = f_uid;
+	}
+	///if((res->min_w > i_b->weight) || (res->min_w == i_b->weight && id < res->min_idx))
+}
+
+static void worker_for_min_bp(void *data, long i, int tid) // callback for kt_for()
+{
+    mc_bp_t* bp = (mc_bp_t *)data;
+	select_min_bp(bp->b_b, bp->ma, bp->b_aux, bp, bp->lock, &(bp->vis[tid]), i, &(bp->res[tid]));
+}
+
+uint32_t best_bp(mc_bp_t *bp, mc_bp_res *res)
+{
+	uint32_t i;
+	clean_mc_bp_res(res);
+	for (i = 0; i < bp->n_thread; i++)
+	{
+		clean_mc_bp_res(&(bp->res[i]));
+	} 
+	kt_for(bp->n_thread, worker_for_min_bp, bp, bp->occ);
+	
+	for (i = 0; i < bp->n_thread; i++)
+	{
+		if(bp->res[i].chain_id == (uint32_t)-1) continue;
+		if(bp->res[i].w <= 0) continue;
+		if((bp->res[i].w > res->w) || (bp->res[i].w == res->w && bp->res[i].id < res->id))
+		{
+			(*res) = bp->res[i];
+		}
+	}
+	if(res->chain_id != (uint32_t)-1) return 1;
+	return 0;
+}
+
+void mc_set_bp_spin(bubble_type* bub, mc_match_t *ma, mc_svaux_t *b_aux, mc_bp_t* bp,
+uint8_t *lock, bits_p *vis, mc_bp_res *res)
+{
+	uint32_t uid, val = 0, c_bid, c_uid;
+	mc_bp_iter i;
+	i.chain_id = res->chain_id;
+	i.bid = res->f_bid;
+	i.uid = res->f_uid;
+	memset(vis->a, 0, vis->n);
+	while (1)
+    {
+        uid = next_uid(&i, bub, &c_bid, &c_uid);
+        if(uid == (uint32_t)-1) break;
+		if(lock[uid] == 0)
+		{
+			val = 1;
+			break;
+		}
+        if(c_bid == res->l_bid && c_uid == res->l_uid) break;
+    }
+
+	if(val == 0)
+	{
+		fprintf(stderr, "ERROR-1\n");
+		return;
+	}
+
+	i.bid = res->f_bid;
+	i.uid = res->f_uid;
+	while (1)
+	{
+		uid = next_uid(&i, bub, &c_bid, &c_uid);
+        if(uid == (uint32_t)-1) break;
+        if(vis->a[uid] == 1) continue;
+        lock[uid] = 1;
+        vis->a[uid] = 1;
+		mc_set_spin(ma, b_aux, uid, -b_aux->s[uid]);
+        if(c_bid == res->l_bid && c_uid == res->l_uid) break;
+	}
+}
+
+double mc_solve_bp_cc(mc_bp_t *bp)
+{
+	mc_bp_res res;
+	memset(bp->lock, 0, bp->b_b->ug->g->n_seq);
+	while (best_bp(bp, &res))
+	{
+		mc_set_bp_spin(bp->b_b, bp->ma, bp->b_aux, bp, bp->lock, &(bp->vis[0]), &res);
+	}
+
+	return mc_score_all(bp->ma, bp->b_aux);
+}
+
+void mc_reset_z_all(const mc_match_t *ma, mc_svaux_t *b)
+{
+	t_w_t z[2];
+	uint32_t k;
+	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] = b->z[k].z[0]; z[1] = b->z[k].z[1];
+        b->z[k].z[0] = b->z[k].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 (b->s[t] > 0) b->z[k].z[0] += e->w;
+            else if (b->s[t] < 0) b->z[k].z[1] += e->w;
+        }
+		if(z[0] != b->z[k].z[0]) fprintf(stderr, "ERROR1-all\n");
+        if(z[1] != b->z[k].z[1]) fprintf(stderr, "ERROR2-all\n");
+	}
+}
+
+void mc_solve_bp(mc_bp_t *bp)
+{
+	double index_time = yak_realtime();
+	uint32_t r = 1;
+	double sc_opt, sc;
+	mc_reset_z_all(bp->ma, bp->b_aux);
+	sc_opt = mc_score_all(bp->ma, bp->b_aux);
+
+	while (1)
+	{
+		sc = mc_solve_bp_cc(bp);
+		fprintf(stderr, "[M::%s::# round: %u] sc_opt: %f, sc: %f\n", __func__, r, sc_opt, sc);
+		if(sc <= sc_opt) break;
+		sc_opt = sc;
+		r++;
+	}
+	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)
+{
+	fprintf(stderr, "# iter: %u, sc_opt: %f, sc-local: %f, sc-global: %f\n", 
+					n_iter, sc_opt, mc_score(ma, b), mc_score_all(ma, b));
+}
+
 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)
 {
 	uint32_t j, k, n_iter = 0;
@@ -816,14 +1238,15 @@ uint32_t mc_solve_cc(const mc_opt_t *opt, const mc_g_t *mg, mc_svaux_t *b, uint3
 	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(opt, mg->e, b);
 	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
 		b->z_opt[b->cc_node[j]] = b->z[b->cc_node[j]]; ///z[0]: positive weight; z[1]: positive weight
 	}
-
-	sc = mc_optimize_local(opt, mg->e, b, &n_iter, &sc_opt);
+	// print_sc(opt, mg->e, b, sc_opt, n_iter);
+	sc = mc_optimize_local(opt, mg->e, b, &n_iter);
 	if (sc > sc_opt) {
 		for (j = 0; j < b->cc_size; ++j) {
 			b->s_opt[b->cc_node[j]] = b->s[b->cc_node[j]];
@@ -836,11 +1259,13 @@ 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]];
 		}
 	}
+	// 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);
 	for (k = 0; k < (uint32_t)opt->n_perturb; ++k) {
 		if (k&1) mc_perturb(opt, mg->e, b);
 		else mc_perturb_node(opt, mg->e, b, 3);
-		sc = mc_optimize_local(opt, mg->e, b, &n_iter, &sc_opt);
+		sc = mc_optimize_local(opt, mg->e, b, &n_iter);
 		// fprintf(stderr, "(%u) sc_opt: %f, sc: %f\n", k, sc_opt, sc);
 		if (sc > sc_opt) {
 			for (j = 0; j < b->cc_size; ++j) {
@@ -854,27 +1279,60 @@ 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->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]];
 	return n_iter;
 }
 
-void mc_solve_core(const mc_opt_t *opt, mc_g_t *mg)
+void mc_init_spin_all(const mc_opt_t *opt, mc_g_t *mg, mc_svaux_t *b)
+{
+	uint32_t st, i;
+	for (st = 0, i = 1; i <= mg->e->n_seq; ++i) {
+		if (i == mg->e->n_seq || mg->e->cc[st]>>32 != mg->e->cc[i]>>32) {
+			b->cc_off = st, b->cc_size = i - st;
+			if (b->cc_size >= 2)
+			{
+				mc_init_spin(opt, mg->e, b);
+			}
+			st = i;
+		}
+	}
+}
+
+
+
+void mc_solve_core(const mc_opt_t *opt, mc_g_t *mg, bubble_type* bub)
 {
 	double index_time = yak_realtime();
 	uint32_t st, i;
 	mc_svaux_t *b;
+	mc_bp_t *bp = NULL;
 	mc_g_cc(mg->e);
+
 	b = mc_svaux_init(mg, opt->seed);
+	if(bub) bp = mc_bp_t_init(mg->e, b, bub, asm_opt.thread_num);
+	/*******************************for debug************************************/
+	if(bp)
+	{
+		mc_init_spin_all(opt, mg, b);
+		mc_solve_bp(bp);
+	} 
+	/*******************************for debug************************************/
+	// fprintf(stderr, "\n\n\n\n\n*************beg-[M::%s::score->%f] ==> Partition\n", __func__, mc_score_all(mg->e, b));
 	for (st = 0, i = 1; i <= mg->e->n_seq; ++i) {
 		if (i == mg->e->n_seq || mg->e->cc[st]>>32 != mg->e->cc[i]>>32) {
 			mc_solve_cc(opt, mg, b, st, i - st);
 			st = i;
 		}
 	}
+	// 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);
+	if(bp) mc_svaux_destroy(b);
+	destroy_mc_bp_t(&bp);
 	fprintf(stderr, "[M::%s::%.3f] ==> Partition\n", __func__, yak_realtime()-index_time);
 }
 
@@ -981,14 +1439,14 @@ void p_nodes(mc_g_t *mg, trans_chain* t_ch, uint8_t* trio_flag)
 	}
 }
 
-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)
+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)
 {
 	mc_opt_t opt;
-	mc_opt_init(&opt);
+	mc_opt_init(&opt, asm_opt.n_perturb, asm_opt.f_perturb, asm_opt.seed);
 	mc_g_t *mg = init_mc_g_t(ug, read_g, s, renew_s);
 	update_mc_edges(mg, ovlp, ta, t_ch, f_rate, is_sys);
 	///debug_mc_g_t(mg);
-	mc_solve_core(&opt, mg);
+	mc_solve_core(&opt, mg, bub);
 
 	if((asm_opt.flag & HA_F_PARTITION) && t_ch)
 	{
diff --git a/rcut.h b/rcut.h
index 2571119..b0dd888 100644
--- a/rcut.h
+++ b/rcut.h
@@ -5,6 +5,7 @@
 #include "kvec.h"
 #include "Overlaps.h"
 #include "Purge_Dups.h"
+#include "hic.h"
 
 typedef struct {
     uint32_t bS, bE;
@@ -42,5 +43,20 @@ typedef struct {
     mc_match_t* e;
 }mc_g_t;
 
-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);
+static inline uint64_t kr_splitmix64(uint64_t x)
+{
+	uint64_t z = (x += 0x9E3779B97F4A7C15ULL);
+	z = (z ^ (z >> 30)) * 0xBF58476D1CE4E5B9ULL;
+	z = (z ^ (z >> 27)) * 0x94D049BB133111EBULL;
+	return z ^ (z >> 31);
+}
+
+static inline double kr_drand_r(uint64_t *x)
+{
+    union { uint64_t i; double d; } u;
+	*x = kr_splitmix64(*x);
+    u.i = 0x3FFULL << 52 | (*x) >> 12;
+    return u.d - 1.0;
+}
+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);
 #endif
\ No newline at end of file