r373

Assembly.cpp

@@ -1636,8 +1636,8 @@ void ug_idx_build(ma_ug_t *ug, int hap_n)
 {
     int flag = asm_opt.flag&HA_F_NO_HPC, i;
     asm_opt.flag -= flag;
-    ha_flt_tab = ha_ft_ug_gen(&asm_opt, &(ug->u), hap_n);
+    // ha_flt_tab = ha_ft_ug_gen(&asm_opt, &(ug->u), hap_n, hap_n);
-    ha_idx = ha_pt_ug_gen(&asm_opt, ha_flt_tab, &(ug->u), hap_n);
+    // ha_idx = ha_pt_ug_gen(&asm_opt, ha_flt_tab, &(ug->u), hap_n);
 
     ha_ovec_buf_t **b = NULL;
     // overlap and correct reads

CommandLines.cpp

@@ -44,6 +44,8 @@ static ko_longopt_t long_options[] = {
     { "dp-er",     ko_required_argument, 330},
     { "max-kocc",     ko_required_argument, 331},
     { "hg-size",     ko_required_argument, 332},
+    { "ul",     ko_required_argument, 333},
+    { "unskew",     ko_no_argument, 334},
 	{ 0, 0, 0 }
 };
 
@@ -155,6 +157,7 @@ void init_opt(hifiasm_opt_t* asm_opt)
     asm_opt->hic_enzymes = NULL;
     asm_opt->hic_reads[0] = NULL;
     asm_opt->hic_reads[1] = NULL;
+    asm_opt->ar = NULL;
     asm_opt->thread_num = 1;
     asm_opt->k_mer_length = 51;
     asm_opt->hic_mer_length = 31;
@@ -245,6 +248,7 @@ void destory_opt(hifiasm_opt_t* asm_opt)
     if(asm_opt->hic_enzymes != NULL) destory_enzyme(asm_opt->hic_enzymes);
     if(asm_opt->hic_reads[0] != NULL) destory_enzyme(asm_opt->hic_reads[0]);
     if(asm_opt->hic_reads[1] != NULL) destory_enzyme(asm_opt->hic_reads[1]);
+    if(asm_opt->ar != NULL) destory_enzyme(asm_opt->ar);
 }
 
 void ha_opt_reset_to_round(hifiasm_opt_t* asm_opt, int round)
@@ -493,6 +497,13 @@ int check_option(hifiasm_opt_t* asm_opt)
         return 0;
     }
 
+    if(asm_opt->ar != NULL && check_hic_reads(asm_opt->ar, "UL") == 0) return 0;
+    if(asm_opt->ar != NULL && asm_opt->ar->n == 0)
+    {
+        fprintf(stderr, "[ERROR] wrong UL reads (--ul)\n");
+        return 0;
+    }
+
     if(asm_opt->b_low_cov < 0)
     {
         fprintf(stderr, "[ERROR] must >= 0 (--b-cov)\n");
@@ -725,6 +736,8 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
         else if (c == 330) asm_opt->dp_e = atof(opt.arg);
         else if (c == 331) asm_opt->max_kmer_cnt = atol(opt.arg);
         else if (c == 332) asm_opt->hg_size = inter_gsize(opt.arg);      
+        else if (c == 333) get_hic_enzymes(opt.arg, &(asm_opt->ar), 0);
+        else if (c == 334) asm_opt->flag |= HA_F_USKEW;
         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);

CommandLines.h

@@ -4,7 +4,7 @@
 #include <pthread.h>
 #include <stdint.h>
 
-#define HA_VERSION "0.16.0-r369"
+#define HA_VERSION "0.16.1-r373"
 
 #define VERBOSE 0
 
@@ -21,6 +21,7 @@
 #define HA_F_HIGH_HET        0x400
 #define HA_F_PARTITION       0x800
 #define HA_F_FAST            0x1000
+#define HA_F_USKEW            0x2000
 
 #define HA_MIN_OV_DIFF       0.02 // min sequence divergence in an overlap
 
@@ -40,6 +41,7 @@ typedef struct {
 	char *extract_list;
     enzyme *hic_reads[2];
     enzyme *hic_enzymes;
+    enzyme *ar;
 	int extract_iter;
     int thread_num;
     int k_mer_length;

Overlaps.cpp

@@ -1475,14 +1475,6 @@ R_to_U* ruIndex, int max_hang, int min_ovlp)
         
                 delete_all_edges(sources, coverage_cut, Get_qn(*h));
                 set_R_to_U(ruIndex, Get_qn(*h), Get_tn(*h), 0, NULL);
-
-                // if(delete_all_edges_carefully(sources, coverage_cut, max_hang, min_ovlp, 
-                // Get_qn(*h))==0)
-                // {
-                //     set_R_to_U(ruIndex, Get_qn(*h), Get_tn(*h), 0);
-                // }
-                // sq->del = 1;
-                // set_R_to_U(ruIndex, Get_qn(*h), Get_tn(*h), 0);
             }
 		    else if (r == MA_HT_TCONT) 
             {
@@ -1491,15 +1483,6 @@ R_to_U* ruIndex, int max_hang, int min_ovlp)
 
                 delete_all_edges(sources, coverage_cut, Get_tn(*h));
                 set_R_to_U(ruIndex, Get_tn(*h), Get_qn(*h), 0, NULL);
-
-                // if(delete_all_edges_carefully(sources, coverage_cut, max_hang, 
-                // min_ovlp, Get_tn(*h)) == 0)
-                // {
-                //     set_R_to_U(ruIndex, Get_tn(*h), Get_qn(*h), 0);
-                //     no_fully_tn_num++;
-                // }
-                // st->del = 1;
-                // set_R_to_U(ruIndex, Get_tn(*h), Get_qn(*h), 0);
             }
         }
     }
@@ -9673,7 +9656,7 @@ ma_hit_t_alloc* sources, R_to_U* ruIndex, uint8_t* r_flag)
         r_flag[rId] = 0;
     }
 
-    return C_bases/R_bases;
+    return R_bases == 0? 0:C_bases/R_bases;
 }
 
 uint32_t get_ug_coverage_aggressive(ma_ug_t *ug, uint32_t uID, asg_t* read_g, 
@@ -13158,7 +13141,7 @@ long long gap_fuzz, bub_label_t* b_mask_t)
         ///asm_opt.purge_simi_thres = asm_opt.purge_simi_rate_hic;
         adjust_utg_by_primary(&copy_ug, copy_sg, TRIO_THRES, sources, reverse_sources, coverage_cut, 
         tipsLen, tip_drop_ratio, stops_threshold, ruIndex, chimeric_rate, drop_ratio, 
-        max_hang, min_ovlp, &new_rtg_edges, &cov, b_mask_t, 1, 1/**0**/);
+        max_hang, min_ovlp, &new_rtg_edges, &cov, b_mask_t, 1, 0);
 
         print_utg(copy_ug, copy_sg, coverage_cut, output_file_name, sources, ruIndex, max_hang, 
         min_ovlp, &new_rtg_edges);
@@ -13182,7 +13165,6 @@ long long gap_fuzz, bub_label_t* b_mask_t)
     }
 
     // char* gfa_name = (char*)malloc(strlen(output_file_name)+50);
-    // FILE* output_file = NULL;
     // sprintf(gfa_name, "%s.pre.clean_d_utg.noseq.gfa", output_file_name);
     // FILE* output_file = fopen(gfa_name, "w");
     // ma_ug_print_simple(ug, sg, coverage_cut, sources, ruIndex, "utg", output_file);
@@ -13744,7 +13726,7 @@ void kt_u_trans_t_symm(kv_u_trans_t *ta, ma_ug_t *ug)
             if (i == n || a[i].tn != a[st].tn)
             {
                 get_u_trans_spec(ta, a[st].tn, a[st].qn, &r_a, &r_n);
-                if(i - st == 1 && r_n == 0)
+                if(i - st == 1 && r_n == 0)///should be always here
                 {
                     st = i;
                     continue;
@@ -14811,9 +14793,37 @@ void label_r_set(buf_t* b, R_to_U* ruIndex, ma_ug_t *ug, uint32_t flag)
     }
 }
 
+inline int trio_check(ma_ug_t *ug, uint32_t *a, uint32_t a_n, uint32_t flag)
+{
+    if(flag != FATHER && flag != MOTHER) return 0;
+    uint32_t flag_occ = 0, non_flag_occ = 0, ambigious = 0, u_n = 0, f, nf, ab, k;
+    for (k = 0; k < a_n; k++) {
+        get_unitig_trio_flag(&(ug->u.a[a[k]>>1]), flag, &f, &nf, &ab);
+        flag_occ += f;
+        non_flag_occ += nf;
+        ambigious += ab;
+        u_n += ug->u.a[a[k]>>1].n;
+    }    
+    if((flag_occ+non_flag_occ) == 0) return 0;
+    if(flag_occ <= ((non_flag_occ+flag_occ)*0.75)) return 0;
+    if(non_flag_occ == 0 && flag_occ >= 20) return 1;
+    if(u_n >= 100)
+    {
+        if(flag_occ < u_n*DOUBLE_CHECK_THRES) return 0;
+    }
+    else if(u_n >= 50)
+    {
+        if(flag_occ < u_n*DOUBLE_CHECK_THRES*0.5) return 0;
+    }
+    else
+    {
+        if(flag_occ < u_n*DOUBLE_CHECK_THRES*0.25) return 0;
+    }
+    return 1;
+}
 
 int asg_arc_cut_trio_long_tip_primary(asg_t *g, ma_ug_t *ug, asg_t *read_sg, ma_hit_t_alloc* reverse_sources,
-R_to_U* ruIndex, uint32_t min_edge_length, float drop_ratio, hap_cov_t *cov, utg_trans_t *o)
+R_to_U* ruIndex, uint32_t min_edge_length, float drop_ratio, uint32_t trio_flag, hap_cov_t *cov, utg_trans_t *o)
 {
     double startTime = Get_T();
     ///the reason is that each read has two direction (query->target, target->query)
@@ -14874,6 +14884,7 @@ R_to_U* ruIndex, uint32_t min_edge_length, float drop_ratio, hap_cov_t *cov, utg
                 } 
 
                 if(ll >= (v_maxLen*drop_ratio)) continue;
+                if(trio_check(ug, b.b.a, b.b.n, trio_flag)) continue;
                 n_reduced++;
 
                 operation = TRIM;
@@ -14934,7 +14945,7 @@ R_to_U* ruIndex, uint32_t min_edge_length, float drop_ratio, hap_cov_t *cov, utg
 }
 
 int asg_arc_cut_trio_long_tip_primary_complex(asg_t *g, ma_ug_t *ug, asg_t *read_sg, ma_hit_t_alloc* reverse_sources,
-R_to_U* ruIndex, uint32_t min_edge_length, float drop_ratio, uint32_t stops_threshold, hap_cov_t *cov, utg_trans_t *o)
+R_to_U* ruIndex, uint32_t min_edge_length, float drop_ratio, uint32_t stops_threshold, hap_cov_t *cov, utg_trans_t *o, uint32_t trio_flag)
 {
     double startTime = Get_T();
     uint32_t v, n_vtx = g->n_seq * 2, n_reduced = 0, convex, in, flag, operation;
@@ -14958,6 +14969,7 @@ R_to_U* ruIndex, uint32_t min_edge_length, float drop_ratio, uint32_t stops_thre
         &max_stop_baseLen, 1, &b);
 
         if(return_flag != MUL_INPUT) continue;
+        if(trio_check(ug, b.b.a, b.b.n, trio_flag)) continue;
         in = convex^1;
         get_real_length(g, convex, &convex);
         convex = convex^1;
@@ -15177,17 +15189,11 @@ hap_cov_t *cov, utg_trans_t *o)
             &max_stop_baseLen, 1, &b);
 
             if(return_flag != END_TIPS) continue;
+            if(trio_check(ug, b.b.a, b.b.n, trio_flag)) continue;
 
             flag = check_different_haps(g, ug, read_sg, av[base_maxLen_i].v, av[i].v,
             reverse_sources, &b_0, &b_1, ruIndex, cov->is_r_het, miniedgeLen, 1);
 
-            // if((av[i].v>>1) == 255 && (av[base_maxLen_i].v>>1) == 33)
-            // if((av[i].v>>1) == 1852 && (av[base_maxLen_i].v>>1) == 2441)
-            // {
-            //     fprintf(stderr, "max-utg%.6ul (%u), p-utg%.6ul (%u)\n", 
-            //     (av[base_maxLen_i].v>>1)+1, av[base_maxLen_i].v, (av[i].v>>1)+1, av[i].v);
-            // }
-
             // #define UNAVAILABLE (uint32_t)-1
             // #define PLOID 0
             // #define NON_PLOID 1
@@ -15445,7 +15451,7 @@ R_to_U* ruIndex, uint32_t positive_flag, float drop_rate)
 
 int asg_arc_cut_trio_long_equal_tips_assembly_complex(asg_t *g, ma_ug_t *ug, asg_t *read_sg, 
 ma_hit_t_alloc* reverse_sources, long long miniedgeLen, R_to_U* ruIndex, uint32_t stops_threshold, 
-hap_cov_t *cov, utg_trans_t *o)
+hap_cov_t *cov, utg_trans_t *o, uint32_t trio_flag)
 {
     double startTime = Get_T();
     uint32_t v, n_vtx = g->n_seq * 2, n_reduced = 0, convex, in, flag;
@@ -15471,6 +15477,7 @@ hap_cov_t *cov, utg_trans_t *o)
         &max_stop_baseLen, 1, &b);
 
         if(return_flag != MUL_INPUT) continue;
+        if(trio_check(ug, b.b.a, b.b.n, trio_flag)) continue;
         in = convex^1;
         get_real_length(g, convex, &convex);
         convex = convex^1;
@@ -15955,10 +15962,11 @@ float drop_ratio, uint32_t trio_flag, float trio_drop_rate, hap_cov_t *cov)
         if(just_bubble_pop == 0)
         {
             ///need consider tangles
-            asg_arc_cut_trio_long_tip_primary(g, ug, read_g, reverse_sources, ruIndex, 2, tip_drop_ratio, cov, NULL);            
+            asg_arc_cut_trio_long_tip_primary(g, ug, read_g, reverse_sources, ruIndex, 2, tip_drop_ratio, trio_flag, cov, NULL);            
+            // if(trio_flag == MOTHER) print_debug_gfa(read_g, ug, coverage_cut, "debug_dups", sources, ruIndex, asm_opt.max_hang_Len, asm_opt.min_overlap_Len);
             asg_arc_cut_trio_long_equal_tips_assembly(g, ug, read_g, reverse_sources, 2, ruIndex, trio_flag, cov, NULL);            
-            asg_arc_cut_trio_long_tip_primary_complex(g, ug, read_g, reverse_sources, ruIndex, 2, tip_drop_ratio, stops_threshold, cov, NULL);
+            asg_arc_cut_trio_long_tip_primary_complex(g, ug, read_g, reverse_sources, ruIndex, 2, tip_drop_ratio, stops_threshold, cov, NULL, trio_flag);
-            asg_arc_cut_trio_long_equal_tips_assembly_complex(g, ug, read_g, reverse_sources, 2, ruIndex, stops_threshold, cov, NULL);
+            asg_arc_cut_trio_long_equal_tips_assembly_complex(g, ug, read_g, reverse_sources, 2, ruIndex, stops_threshold, cov, NULL, trio_flag);
             detect_chimeric_by_topo(g, ug, read_g, reverse_sources, 2, stops_threshold, chimeric_rate, ruIndex, NULL, cov->is_r_het);
             ///need consider tangles
             ///note we need both the read graph and the untig graph
@@ -15972,10 +15980,11 @@ float drop_ratio, uint32_t trio_flag, float trio_drop_rate, hap_cov_t *cov)
     }
 
     ///print_debug_gfa(read_g, ug, coverage_cut, "debug_dups", sources, ruIndex, asm_opt.max_hang_Len, asm_opt.min_overlap_Len);
-
+    magic_trio_phasing(g, ug, read_g, coverage_cut, sources, reverse_sources, 2, ruIndex, trio_flag, trio_drop_rate); 
     resolve_tangles(ug, read_g, reverse_sources, 20, 100, 0.05, 0.2, ruIndex, cov->is_r_het, trio_flag, drop_ratio);    
     drop_semi_circle(ug, g, read_g, reverse_sources, ruIndex, cov->is_r_het);
     all_to_all_deduplicate(ug, read_g, coverage_cut, sources, trio_flag, trio_drop_rate, reverse_sources, ruIndex, cov->is_r_het, DOUBLE_CHECK_THRES, asm_opt.trio_flag_occ_thres);
+    // if(trio_flag == MOTHER) print_untig_by_read(ug, "m54329U_190827_173812/30214441/ccs", (uint32_t)-1, NULL, NULL, "bf-16");
     if(is_first)
     {
         is_first = 0;
@@ -16031,10 +16040,10 @@ int just_bubble_pop, float drop_ratio, hap_cov_t *cov)
         if(just_bubble_pop == 0)
         {
             ///need consider tangles  
-            asg_arc_cut_trio_long_tip_primary(g, ug, read_g, reverse_sources, ruIndex, 2, tip_drop_ratio, cov, NULL); 
+            asg_arc_cut_trio_long_tip_primary(g, ug, read_g, reverse_sources, ruIndex, 2, tip_drop_ratio, (uint32_t)-1, cov, NULL); 
             asg_arc_cut_trio_long_equal_tips_assembly(g, ug, read_g, reverse_sources, 2, ruIndex, (uint32_t)-1, cov, NULL);   
-            asg_arc_cut_trio_long_tip_primary_complex(g, ug, read_g, reverse_sources, ruIndex, 2, tip_drop_ratio, stops_threshold, cov, NULL);   
+            asg_arc_cut_trio_long_tip_primary_complex(g, ug, read_g, reverse_sources, ruIndex, 2, tip_drop_ratio, stops_threshold, cov, NULL, (uint32_t)-1);   
-            asg_arc_cut_trio_long_equal_tips_assembly_complex(g, ug, read_g, reverse_sources, 2, ruIndex, stops_threshold, cov, NULL);
+            asg_arc_cut_trio_long_equal_tips_assembly_complex(g, ug, read_g, reverse_sources, 2, ruIndex, stops_threshold, cov, NULL, (uint32_t)-1);
             detect_chimeric_by_topo(g, ug, read_g, reverse_sources, 2, stops_threshold, chimeric_rate, ruIndex, NULL, cov->is_r_het);    
             if(round != T_ROUND)
             {
@@ -16094,10 +16103,10 @@ int min_ovlp, hap_cov_t *cov)
                 asg_pop_bubble_primary_trio(ug, NULL, (uint32_t)-1, DROP, cov, o, 1);
                 
                 ///need consider tangles  
-                asg_arc_cut_trio_long_tip_primary(g, ug, read_g, reverse_sources, ruIndex, 2, tip_drop_ratio, cov, o);              
+                asg_arc_cut_trio_long_tip_primary(g, ug, read_g, reverse_sources, ruIndex, 2, tip_drop_ratio, (uint32_t)-1, cov, o);              
                 asg_arc_cut_trio_long_equal_tips_assembly(g, ug, read_g, reverse_sources, 2, ruIndex, (uint32_t)-1, cov, o);   
-                asg_arc_cut_trio_long_tip_primary_complex(g, ug, read_g, reverse_sources, ruIndex, 2, tip_drop_ratio, stops_threshold, cov, o);   
+                asg_arc_cut_trio_long_tip_primary_complex(g, ug, read_g, reverse_sources, ruIndex, 2, tip_drop_ratio, stops_threshold, cov, o, (uint32_t)-1);   
-                asg_arc_cut_trio_long_equal_tips_assembly_complex(g, ug, read_g, reverse_sources, 2, ruIndex, stops_threshold, cov, o);
+                asg_arc_cut_trio_long_equal_tips_assembly_complex(g, ug, read_g, reverse_sources, 2, ruIndex, stops_threshold, cov, o, (uint32_t)-1);
                 detect_chimeric_by_topo(g, ug, read_g, reverse_sources, 2, stops_threshold, chimeric_rate, ruIndex, o, cov->is_r_het);    
                 
                 cur_cons = get_graph_statistic(g);
@@ -23920,6 +23929,8 @@ int write_ruIndex(R_to_U* ruIndex, char* read_file_name)
     fwrite(ruIndex->index, sizeof(ruIndex->index[0]), ruIndex->len, fp);
     fwrite(R_INF.trio_flag, sizeof(R_INF.trio_flag[0]), ruIndex->len, fp);
     // fwrite(ruIndex->is_het, 1, ruIndex->len, fp);
+    fwrite(&(asm_opt.hom_global_coverage_set), sizeof(asm_opt.hom_global_coverage_set), 1, fp);
+    fwrite(&(asm_opt.hom_global_coverage), sizeof(asm_opt.hom_global_coverage), 1, fp);
     free(index_name);    
     fflush(fp);
     fclose(fp);
@@ -23947,6 +23958,9 @@ int load_ruIndex(R_to_U* ruIndex, char* read_file_name)
     // CALLOC(ruIndex->is_het, ruIndex->len);
     // f_flag += fread(ruIndex->is_het, 1, ruIndex->len, fp);
 
+    f_flag += fread(&(asm_opt.hom_global_coverage_set), sizeof(asm_opt.hom_global_coverage_set), 1, fp);
+    f_flag += fread(&(asm_opt.hom_global_coverage), sizeof(asm_opt.hom_global_coverage), 1, fp);
+
     free(index_name);    
     fflush(fp);
     fclose(fp);
@@ -30869,12 +30883,10 @@ ma_sub_t **coverage_cut_ptr, int debug_g)
 
     output_contig_graph_primary_pre(sg, coverage_cut, o_file, sources, reverse_sources, 
         asm_opt.small_pop_bubble_size, asm_opt.max_short_tip, ruIndex, max_hang_length, mini_overlap_length);
-
     if (asm_opt.flag & HA_F_VERBOSE_GFA)
     {
         write_debug_graph(sg, sources, coverage_cut, output_file_name, n_read, reverse_sources, ruIndex);
         debug_gfa:;
-        set_hom_global_coverage(&asm_opt, sg, coverage_cut, sources, reverse_sources, ruIndex, max_hang_length, mini_overlap_length);
     }
     
     if (ha_opt_triobin(&asm_opt) && ha_opt_hic(&asm_opt))

Overlaps.h

@@ -190,7 +190,7 @@ void print_gfa(asg_t *g);
 typedef struct { size_t n, m; uint64_t *a; } asg64_v;
 
 
-typedef struct { size_t n, m; ma_utg_t *a; } ma_utg_v;
+typedef struct { size_t n, m; ma_utg_t *a; int h;} ma_utg_v;
 
 typedef struct {
 	ma_utg_v u;
@@ -891,6 +891,12 @@ typedef struct {
 	buf_t b0, b1;
 } utg_trans_t;
 
+typedef struct {
+	ma_ug_t *ug;
+	kvec_t(uint64_t) idx;
+	kvec_t(uint32_t) dst;
+} spg_t;
+
 void init_hc_links(hc_links* link, uint64_t ug_num, trans_chain* t_ch);
 void destory_hc_links(hc_links* link);
 uint64_t get_bub_pop_max_dist(asg_t *g, buf_t *b);

Purge_Dups.cpp

@@ -51,7 +51,7 @@ typedef struct {
 typedef struct {
     hap_alignment_struct* buf;
     uint32_t num_threads;
-
+    uint8_t *hh;
     ma_ug_t *ug;
     asg_t *read_g;
     ma_hit_t_alloc* sources;
@@ -471,6 +471,30 @@ void destory_hap_alignment_struct(hap_alignment_struct* x)
     kv_destroy(x->u_can.a);
 }
 
+uint8_t *init_pip_hh(asg_t *rg, ma_hit_t_alloc* reverse_sources, ma_sub_t *coverage_cut, long long sc)
+{
+    uint8_t *c = NULL; CALLOC(c, rg->n_seq);
+    ma_hit_t *h = NULL;
+    uint32_t i, k;
+    long long R_Base, C_Base;
+    for (i = 0; i < rg->n_seq; i++) {
+        if(sc <= 0){
+            c[i] = 1;
+            continue;
+        }
+        R_Base = (coverage_cut[i].e - coverage_cut[i].s);
+        for (k = 0, C_Base = 0; k < reverse_sources[i].length; k++){
+            h = &(reverse_sources[i].buffer[k]);
+            C_Base += (Get_qe((*h)) - Get_qs((*h)));
+        }
+        C_Base = (R_Base!=0?(C_Base/R_Base):0);
+        C_Base /= sc;
+        c[i] = 1;
+        if(C_Base < REV_W) c[i] = REV_W -  C_Base;
+    }
+    return c;
+}
+
 void init_hap_alignment_struct_pip(hap_alignment_struct_pip* x, uint32_t num_threads, uint32_t n_seq,
 ma_ug_t *ug, asg_t *read_g,  ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, ma_sub_t *coverage_cut, 
 uint64_t* position_index, float Hap_rate, int max_hang, int min_ovlp, float chain_rate, hap_overlaps_list* all_ovlp, hap_cov_t *cov)
@@ -496,6 +520,17 @@ uint64_t* position_index, float Hap_rate, int max_hang, int min_ovlp, float chai
     x->chain_rate = chain_rate;
     x->all_ovlp = all_ovlp;
     x->cov = cov;
+    long long sc = -1;
+    if(asm_opt.hom_global_coverage_set) {
+        sc = asm_opt.hom_global_coverage*1.75;
+    }
+    else {
+        if(asm_opt.hom_global_coverage > 0){
+            sc = ((int)(((double)asm_opt.hom_global_coverage)/((double)HOM_PEAK_RATE)))*1.75;
+        }
+    }
+    if(sc <= 0) sc = -1;
+    x->hh = init_pip_hh(read_g, reverse_sources, coverage_cut, sc);
 }
 
 
@@ -508,6 +543,7 @@ void destory_hap_alignment_struct_pip(hap_alignment_struct_pip* x)
     }
 
     free(x->buf);
+    free(x->hh);
 }
 
 void init_hap_overlaps_list(hap_overlaps_list* x, uint32_t num)
@@ -677,16 +713,10 @@ void deduplicate_edge(kvec_asg_arc_t_offset* u_buffer)
     long long i = u_buffer->a.n - 1, k, i_off;
     uint32_t v = u_buffer->a.a[i].x.ul>>33, m;
 
-    for (; i >= 0; i--)
+    for (; i >= 0; i--) {
-    {
+        if((u_buffer->a.a[i].x.ul>>33) != v) break;
-        if((u_buffer->a.a[i].x.ul>>33) != v)
-        {
-            break;
-        } 
     }
 
-    ///fprintf(stderr, "u_buffer->a.n: %u, i: %lld\n", u_buffer->a.n, i);
-
     i = i + 1;
     for (m = i; i < (long long)u_buffer->a.n; i++)
     {
@@ -960,7 +990,7 @@ ma_hit_t_alloc* reverse_sources, asg_t *read_g, R_to_U* ruIndex, double* Match,
     uint32_t i, j, qn, tn, is_Unitig, uId, min_count = 0, max_count = 0, cutoff = 0;;
     for (i = 0; i < Len; i++)
     {
-        if(cutoff > CUTOFF_THRES)
+        if(cutoff > CUTOFF_THRES && cutoff > (Len>>1))
         {
             max_count = 0;
             min_count = Len;
@@ -2254,6 +2284,7 @@ uint32_t* xBeg, uint32_t* xEnd, uint32_t* yBeg, uint32_t* yEnd)
 
 #define generic_key(x) (x)
 KRADIX_SORT_INIT(i32, int32_t, generic_key, sizeof(int32_t))
+KRADIX_SORT_INIT(ru32, uint32_t, generic_key, sizeof(uint32_t))
 
 long long get_chain_score(ma_utg_t *xReads, asg_t *read_g, kvec_asg_arc_t_offset* u_buffer, kvec_t_i32_warp* tailIndex, kvec_t_i32_warp* idx, 
 ma_hit_t_alloc* reverse_sources, long long xBegPos, long long xEndPos)
@@ -2820,7 +2851,7 @@ long long y_readLen)
 {
     u_can->a.n = 0;
     if(u_buffer->a.n == 0) return;
-    uint32_t i = 0, anchor_i = 0, m = 1, break_point = (uint32_t)-1, is_merge;
+    uint32_t i = 0, /**anchor_i = 0, **/m = 1, break_point = (uint32_t)-1, is_merge;
 
     qsort(u_buffer->a.a, u_buffer->a.n, sizeof(asg_arc_t_offset), cmp_hap_alignment_chaining);
 
@@ -2832,14 +2863,14 @@ long long y_readLen)
         if(u_buffer->a.a[m-1].x.el == u_buffer->a.a[i].x.el)
         {
             if(u_buffer->a.a[m-1].Off == u_buffer->a.a[i].Off) is_merge = 1;
-            if(is_merge == 0 && (Get_xOff(u_buffer->a.a[m-1].Off)==Get_xOff(u_buffer->a.a[i].Off)))
+            ///I think we don't need the following merging
-            {
+            // if(is_merge == 0 && (Get_xOff(u_buffer->a.a[m-1].Off)==Get_xOff(u_buffer->a.a[i].Off)))
-                if((Get_yOff(u_buffer->a.a[i].Off)-(Get_yOff(u_buffer->a.a[m-1].Off))) ==
+            // {
-                  (i-anchor_i))///not sure why, does it use for tolerate indels in overlaps? 
+            //     if((Get_yOff(u_buffer->a.a[i].Off)-(Get_yOff(u_buffer->a.a[m-1].Off))) == (i-anchor_i))///not sure why, does it use for tolerate indels in overlaps? 
-                {
+            //     {
-                    is_merge = 1;
+            //         is_merge = 1;
-                }
+            //     }
-            }
+            // }
 
             if(is_merge)
             {
@@ -2848,7 +2879,7 @@ long long y_readLen)
             }
         } 
         u_buffer->a.a[m] = u_buffer->a.a[i];
-        anchor_i = i;
+        // anchor_i = i;
         if(u_buffer->a.a[m].x.el != u_buffer->a.a[m-1].x.el) break_point = m;
         m++;
     }
@@ -2957,6 +2988,7 @@ static void hap_alignment_advance_worker(void *_data, long eid, int tid)
     kvec_t_u8_warp* flag_vec = &(hap_buf->buf[tid].u_buffer_flag);
     uint64_t cov_threshold = hap_buf->cov_threshold;
     hap_cov_t *cov = hap_buf->cov;
+    uint8_t *hh = hap_buf->hh, hhc;
     if(hap_buf->cov_threshold < 0) cov_threshold = (uint64_t)-1;
     ma_utg_t *xReads = NULL, *yReads = NULL;
     ma_hit_t_alloc *xR = NULL;
@@ -3036,7 +3068,7 @@ static void hap_alignment_advance_worker(void *_data, long eid, int tid)
         for (k = 0; k < xReads->n; k++)
         {
             xR = &(reverse_sources[xReads->a[k]>>33]);
-
+            hhc = hh[xReads->a[k]>>33];
             for (j = 0; j < xR->length; j++)
             {
                 h = &(xR->buffer[j]);
@@ -3078,7 +3110,7 @@ static void hap_alignment_advance_worker(void *_data, long eid, int tid)
                 if(((t_offset.Off>>32) == (uint32_t)-1) || (((uint32_t)t_offset.Off) == (uint32_t)-1)) continue;
 
                 t_offset.x = t;
-                t_offset.weight = 1;
+                t_offset.weight = hhc;
 
                 kv_push(asg_arc_t_offset, u_buffer->a, t_offset);
             }
@@ -5255,6 +5287,243 @@ void collect_purge_trans_cov(ma_ug_t *ug, hap_overlaps_list* ha, hap_cov_t *cov,
     }
 }
 
+/**
+typedef struct {
+    uint32_t qn, qs, qe;
+    uint32_t tn, ts, te;
+    uint32_t oid;
+    uint8_t rev;
+}scg_hits;
+
+typedef struct {
+    scg_hits *a;
+    size_t n,m;
+    kvec_t(uint64_t) idx;
+}scg_hits_v;
+
+#define scg_key_qtn(a) ((((uint64_t)(a).qn)<<32)|((uint64_t)(a).tn))
+KRADIX_SORT_INIT(scg_qtn, scg_hits, scg_key_qtn, 8)
+#define scg_key_qts(a) ((((uint64_t)(a).qs)<<32)|((uint64_t)(a).ts))
+KRADIX_SORT_INIT(scg_qts, scg_hits, scg_key_qts, 8)
+#define scg_key_qte(a) ((((uint64_t)(a).qe)<<32)|((uint64_t)(a).te))
+KRADIX_SORT_INIT(scg_qte, scg_hits, scg_key_qte, 8)
+#define scg_key_rev(a) ((a).rev)
+KRADIX_SORT_INIT(scg_rev, scg_hits, scg_key_rev, member_size(scg_hits, rev))
+
+inline void rev_scg_hits(scg_hits *p, spg_t *scg)
+{
+    if(p->rev){
+        uint32_t t;
+        p->ts = scg->ug->u.a[p->tn].len - p->ts - 1;
+        p->te = scg->ug->u.a[p->tn].len - (p->te - 1) - 1;
+        t = p->ts; p->ts = p->te; p->te = t; p->te++;
+    }
+}
+
+#define arc_first(g, v) ((g)->arc[(g)->idx[(v)]>>32])
+scg_hits_v *get_scg_hits_v(scg_hits_v *vp, spg_t *scg)
+{   
+    scg_hits_v *hh = NULL; CALLOC(hh, 1);
+    ma_utg_v *u = &(scg->ug->u);
+    uint32_t i, mn, *ma = NULL;
+    uint64_t offset, *idx = NULL;
+
+    for (i = 0; i < scg->idx.n; i++) {
+        mn = (uint32_t)scg->idx.a[i];
+        ma = scg->dst.a + (scg->idx.a[i]>>32);
+    }
+    
+
+
+
+    return hh;
+}
+
+void refine_scg(spg_t *scg, ma_ug_t *lug, hap_overlaps_list *ha, hap_cov_t *cov, uint64_t* position_index)
+{
+    scg_hits_v vp; kv_init(vp);
+    uint32_t v, i, k, st, c[2];
+    hap_overlaps *x = NULL;
+    u_trans_t *z = NULL;
+    scg_hits *p = NULL;
+
+    for (v = 0; v < cov->t_ch->k_trans.n; v++){
+        z = &(cov->t_ch->k_trans.a[v]);
+        if(z->del) continue;
+        kv_pushp(scg_hits, vp, &p);
+        p->rev = z->rev; p->oid = (uint32_t)-1;
+        p->qn = z->qn; p->qs = z->qs; p->qe = z->qe;
+        p->tn = z->tn; p->ts = z->ts; p->te = z->te;
+        // rev_scg_hits(p, scg);
+        kv_pushp(scg_hits, vp, &p);
+        p->rev = z->rev; p->oid = (uint32_t)-1;
+        p->qn = z->tn; p->qs = z->ts; p->qe = z->te;
+        p->tn = z->qn; p->ts = z->qs; p->te = z->qe;
+        // rev_scg_hits(p, scg);
+    }
+    
+    for (v = 0; v < ha->num; v++){
+        for (i = 0; i < ha->x[v].a.n; i++){
+            x = &(ha->x[v].a.a[i]);
+            st = cov->t_ch->k_trans.n;
+            chain_origin_trans_uid_by_purge(x, lug, cov, position_index);
+            for (k = st; k < cov->t_ch->k_trans.n; k++){
+                z = &(cov->t_ch->k_trans.a[k]);
+                if(z->del) continue;
+                kv_pushp(scg_hits, vp, &p);
+                p->rev = z->rev; p->oid = v;
+                p->qn = z->qn; p->qs = z->qs; p->qe = z->qe;
+                p->tn = z->tn; p->ts = z->ts; p->te = z->te;
+                // rev_scg_hits(p, scg);
+                kv_pushp(scg_hits, vp, &p);
+                p->rev = z->rev; p->oid = v;
+                p->qn = z->tn; p->qs = z->ts; p->qe = z->te;
+                p->tn = z->qn; p->ts = z->qs; p->te = z->qe;
+                // rev_scg_hits(p, scg);
+            }
+            cov->t_ch->k_trans.n = st;
+        }
+    }
+
+    ///two scg_hits might be totally equal; must remove first
+    radix_sort_scg_qtn(vp.a, vp.a + vp.n);
+    for (st = 0, i = 1; i <= vp.n; ++i){
+        if (i == vp.n || vp.a[i].qn != vp.a[st].qn || vp.a[i].tn != vp.a[st].tn){
+           if(i - st > 1) radix_sort_scg_rev(vp.a+st, vp.a+i);
+           for (v = st, c[0] = c[1] = 0; v < i; v++) c[vp.a[v].rev]++;
+           if(c[0]>1) radix_sort_scg_qts(vp.a+st, vp.a+st+c[0]);
+           if(c[1]>1) radix_sort_scg_qts(vp.a+st+c[0], vp.a+st+c[0]+c[1]);
+           st = i;
+        }
+    }
+    for (st = 0, i = 1; i <= vp.n; ++i){
+        if (i == vp.n || vp.a[i].rev != vp.a[st].rev || 
+                vp.a[i].qn != vp.a[st].qn || vp.a[i].tn != vp.a[st].tn || 
+                vp.a[i].qs != vp.a[st].qs || vp.a[i].ts != vp.a[st].ts)
+        {
+           if(i - st > 1) radix_sort_scg_qte(vp.a+st, vp.a+i);
+           st = i;
+        }
+    }
+    for (st = 0, i = 1, k = 0; i <= vp.n; ++i){
+        if (i == vp.n || vp.a[i].rev != vp.a[st].rev ||
+                vp.a[i].qn != vp.a[st].qn || vp.a[i].tn != vp.a[st].tn || 
+                vp.a[i].qs != vp.a[st].qs || vp.a[i].ts != vp.a[st].ts ||
+                vp.a[i].qe != vp.a[st].qe || vp.a[i].te != vp.a[st].te)
+        {
+           vp.a[k] = vp.a[st];
+           k++;
+           st = i;
+        }
+    }
+    ///build idx
+    vp.n = k; kv_resize(uint64_t, vp.idx, scg->ug->u.n); vp.idx.n = scg->ug->u.n; 
+    memset(vp.idx.a, 0, vp.idx.n*sizeof(uint64_t));
+    for (st = 0, i = 1; i <= vp.n; ++i)
+    {
+        if (i == vp.n || vp.a[i].qn != vp.a[st].qn)
+        {
+            vp.idx.a[vp.a[st].qn] = (uint64_t)st << 32 | (i - st);
+            st = i;
+        }
+    }
+
+
+    kv_destroy(vp); kv_destroy(vp.idx);
+}
+**/
+uint32_t seed_uid(ma_utg_t *vu, uint64_t* ps_idx, R_to_U* ruIndex, ma_ug_t *rug)
+{
+    int64_t v_i, v, w, w_i, wb, we, vb, ve, k;
+    uint32_t uid, is_u;
+    ma_utg_t *wu = NULL;
+    for (v_i = 0; v_i < vu->n; v_i++) {
+        v = vu->a[v_i]>>32;
+        get_R_to_U(ruIndex, v>>1, &uid, &is_u);
+        if(is_u == 0 || uid == (uint32_t)-1 || ps_idx[v>>1] == (uint64_t)-1) continue;
+        w_i = (uint32_t)ps_idx[v>>1];
+        wu = &(rug->u.a[uid]);
+        w = wu->a[w_i]>>32;
+        if((v>>1)!=(w>>1)) continue;
+        vb = 0; ve = vu->n; ///[vb, ve)
+        if(v == w){ ///[wb, we)
+            wb = w_i - v_i;
+            we = wb + vu->n;
+            if(wb < 0 || we > wu->n) continue;
+            for (k = 0; k < vu->n; k++){
+                if((vu->a[k+vb]>>32) != (wu->a[k+wb]>>32)) break;
+            }
+            if(k >= vu->n) return uid;
+        } else {
+            wb = w_i + 1 - (ve - v_i);
+            we = wb + vu->n;
+            if(wb < 0 || we > wu->n) continue;
+            for (k = 0; k < vu->n; k++){
+                if((vu->a[k+vb]>>32) != ((wu->a[we-k-1]>>32)^1)) break;
+            }
+            if(k >= vu->n) return uid;
+        }
+    }
+    return (uint32_t)-1;
+}
+
+void filter_ovlp_vecs(hap_overlaps_list* ha, uint32_t *a, uint32_t a_n)
+{
+    uint32_t st, k, i, m, v, w;
+    int idx;
+    radix_sort_ru32(a, a + a_n);
+    for (st = 0, m = 0, k = 1; k <= a_n; k++){
+        if(k == a_n || a[k] != a[st]){
+            a[m++] = a[st];
+            st = k;
+        }
+    }
+    a_n = m;
+    if(a_n < 2) return;
+    for (k = 0; k < a_n; k++){
+        v = a[k];
+        for (i = k+1; i < a_n; i++) {
+            w = a[i];
+            idx = get_specific_hap_overlap(&(ha->x[v]), v, w);
+            if(idx != -1) ha->x[v].a.a[idx].status = DELETE;
+            idx = get_specific_hap_overlap(&(ha->x[w]), w, v);
+            if(idx != -1) ha->x[w].a.a[idx].status = DELETE;
+        }
+    }
+}
+
+void filter_ovlp_scg(hap_overlaps_list* ha, uint64_t* ps_idx, R_to_U* ruIndex, ma_ug_t *rug, spg_t *scg)
+{
+    uint32_t i, k, v, *ma = NULL, mn, luid;
+    ma_utg_v *pp = &(scg->ug->u);
+    kvec_t(uint32_t) vv; kv_init(vv);
+    for (i = 0; i < scg->idx.n; i++){
+        ma = scg->dst.a + (scg->idx.a[i]>>32);
+        mn = (uint64_t)scg->idx.a[i];
+        if(mn < 2) continue;
+        for (k = 0, vv.n = 0; k < mn; k++){
+            luid = seed_uid(&(pp->a[ma[k]>>1]), ps_idx, ruIndex, rug);
+            if(luid == (uint32_t)-1) {
+                fprintf(stderr, "ERROR-scg\n");
+                continue;
+            }
+            kv_push(uint32_t, vv, luid);
+        }
+        if(vv.n < 2) continue;
+        filter_ovlp_vecs(ha, vv.a, vv.n);
+    }
+
+    for (v = 0; v < ha->num; v++){
+        for (i = 0, k = 0; i < ha->x[v].a.n; i++){
+            if(ha->x[v].a.a[i].status == DELETE) continue;
+            ha->x[v].a.a[k++] = ha->x[v].a.a[i];
+        }
+        ha->x[v].a.n = k;
+    }
+
+    kv_destroy(vv);
+}
+
 void purge_dups(ma_ug_t *ug, asg_t *read_g, ma_sub_t* coverage_cut, ma_hit_t_alloc* sources, 
 ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, kvec_asg_arc_t_warp* edge, float density, 
 uint32_t purege_minLen, int max_hang, int min_ovlp, float drop_ratio, uint32_t just_contain, 
@@ -5343,7 +5612,9 @@ uint32_t just_coverage, hap_cov_t *cov, uint32_t collect_p_trans, uint32_t colle
     ///if(debug_enable) print_all_purge_ovlp(ug, &all_ovlp);
     filter_hap_overlaps_by_length(&all_ovlp, purege_minLen);
 
-    normalize_hap_overlaps_advance(&all_ovlp, &back_all_ovlp, ug, read_g, reverse_sources, ruIndex);
+    // normalize_hap_overlaps_advance(&all_ovlp, &back_all_ovlp, ug, read_g, reverse_sources, ruIndex);
+    pg = init_p_g_t(ug, cov, read_g);
+    normalize_hap_overlaps_advance_by_p_g_t(&all_ovlp, &back_all_ovlp, ug, read_g, reverse_sources, ruIndex, pg, cov, 0.8);
 
     if(collect_p_trans && collect_p_trans_f == 0)
     {
@@ -5352,7 +5623,7 @@ uint32_t just_coverage, hap_cov_t *cov, uint32_t collect_p_trans, uint32_t colle
     
     if(asm_opt.polyploidy <= 2)
     {
-        mc_solve(&all_ovlp, cov->t_ch, NULL, ug, read_g, 0.8, R_INF.trio_flag, 1, NULL, 1, NULL, NULL);
+        mc_solve(&all_ovlp, cov->t_ch, NULL, ug, read_g, 0.8, R_INF.trio_flag, 1, NULL, 1, NULL, NULL, 1);
     } 
     
     if(collect_p_trans && collect_p_trans_f == 1)
@@ -5360,11 +5631,6 @@ uint32_t just_coverage, hap_cov_t *cov, uint32_t collect_p_trans, uint32_t colle
         collect_purge_trans_cov(ug, &all_ovlp, cov, position_index);
     } 
 
-    pg = init_p_g_t(ug, cov, read_g);
-    
-    normalize_hap_overlaps_advance_by_p_g_t(&all_ovlp, &back_all_ovlp, ug, read_g, reverse_sources, 
-    ruIndex, pg, cov, 0.8);
-
     ///normalize_hap_overlaps_advance(&all_ovlp, &back_all_ovlp, ug, read_g, reverse_sources, ruIndex);
     ///debug_hap_overlaps(&all_ovlp, &back_all_ovlp);
 

Purge_Dups.h

@@ -12,6 +12,7 @@
 #define ALTER_COV_THRES 0.9
 #define REAL_ALTER_THRES 0.25
 #define CHAIN_FILTER_RATE 0.7
+#define REV_W 8
 
 #define SELF_EXIST 0
 #define REVE_EXIST 1

anchor.cpp

@@ -76,7 +76,7 @@ void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_reg
 	rlen = Get_READ_LENGTH(R_INF, rid); // read length
 
 	// get the list of anchors
-	ha_sketch(ucr->seq, ucr->length, asm_opt.mz_win, asm_opt.k_mer_length, 0, !(asm_opt.flag & HA_F_NO_HPC), &ab->mz, ha_flt_tab, asm_opt.mz_sample_dist, k_flag, dbg_ct, NULL, -1, asm_opt.dp_min_len, -1, sp, asm_opt.mz_rewin, 0);
+	mz1_ha_sketch(ucr->seq, ucr->length, asm_opt.mz_win, asm_opt.k_mer_length, 0, !(asm_opt.flag & HA_F_NO_HPC), &ab->mz, ha_flt_tab, asm_opt.mz_sample_dist, k_flag, dbg_ct, NULL, -1, asm_opt.dp_min_len, -1, sp, asm_opt.mz_rewin, 0);
 	// minimizer of queried read
 	if (ab->mz.m > ab->old_mz_m) {
 		ab->old_mz_m = ab->mz.m;

hic.cpp

@@ -2378,7 +2378,7 @@ void dfs_bubble(asg_t *g, kvec_t_u32_warp* stack, kvec_t_u32_warp* result, uint3
 uint32_t get_unitig_het_arb(ma_ug_t* ug, uint32_t uid, uint8_t *r_het_flag, kv_u_trans_t *ref, uint32_t m_het_occ,
 uint32_t m_het_label, uint32_t p_het_label, uint32_t n_het_label)
 {
-    if(ref && u_trans_n(*ref, uid) > 0) return m_het_label;
+    if(ref && ref->idx.n > 0 && u_trans_n(*ref, uid) > 0) return m_het_label;
     ma_utg_t *u = &(ug->u.a[uid]);
     uint32_t k, rId;
     uint32_t het_occ, hom_occ;
@@ -2411,7 +2411,6 @@ void identify_bubbles(ma_ug_t* ug, bubble_type* bub, uint8_t *r_het_flag, kv_u_t
     uint64_t pathLen;
     bub->ug = ug; 
     bub->b_bub = bub->b_end_bub = bub->tangle_bub = bub->cross_bub = bub->mess_bub = 0;
-
     if(bub->round_id == 0)
     {
         buf_t b; memset(&b, 0, sizeof(buf_t)); b.a = (binfo_t*)calloc(n_vtx, sizeof(binfo_t));
@@ -2421,10 +2420,8 @@ void identify_bubbles(ma_ug_t* ug, bubble_type* bub, uint8_t *r_het_flag, kv_u_t
         bub->b_ug = NULL; kv_init(bub->chain_weight);
         bub->b_s_idx.n = ug->g->n_seq;
         memset(bub->b_s_idx.a, -1, bub->b_s_idx.n * sizeof(uint64_t));
-
         CALLOC(bub->index, n_vtx);
         for (i = 0; i < ug->g->n_seq; i++) ug->g->seq[i].c = 0;
-
         for (v = 0; v < n_vtx; ++v) 
         {
             if(ug->g->seq[v>>1].del) continue;
@@ -2445,7 +2442,6 @@ void identify_bubbles(ma_ug_t* ug, bubble_type* bub, uint8_t *r_het_flag, kv_u_t
             }
         }
 
-
         kvec_t_u32_warp stack, result;
         kv_init(stack.a); kv_init(result.a);
         for (v = 0; v < n_vtx; ++v) 
@@ -2498,16 +2494,14 @@ void identify_bubbles(ma_ug_t* ug, bubble_type* bub, uint8_t *r_het_flag, kv_u_t
                 }
             }
         }
-        
         kv_push(uint32_t, bub->num, bub->list.n);
         free(b.a); free(b.S.a); free(b.T.a); free(b.b.a); free(b.e.a);
         bub->f_bub = bub->num.n - 1; ///bub->s_bub = bub->num.n - 1;
         
-       for (i = 0; i < ug->g->n_seq; i++)
+        for (i = 0; i < ug->g->n_seq; i++)
-       {
+        {
-           bub->index[i] = get_unitig_het_arb(ug, i, r_het_flag, ref, 20, M_het(*bub), P_het(*bub), (uint32_t)-1);
+            bub->index[i] = get_unitig_het_arb(ug, i, r_het_flag, ref, 20, M_het(*bub), P_het(*bub), (uint32_t)-1);
-       }
+        }
-
         for (i = 0; i < bub->f_bub; i++)
         {
             get_bubbles(bub, i, &beg, &sink, &a, &n, &pathLen);
@@ -2555,7 +2549,6 @@ void identify_bubbles(ma_ug_t* ug, bubble_type* bub, uint8_t *r_het_flag, kv_u_t
                 bub->b_s_idx.a[v] |= i;
             }        
         }
-        
         for (i = 0; i < ug->g->n_seq; i++)
         {
             if(bub->index[i] == M_het(*bub)) bub->index[i] = P_het(*bub);
@@ -4923,7 +4916,7 @@ int load_hc_hits(kvec_pe_hit* hits, ma_ug_t* ug, const char *fn)
     }
 
     fclose(fp);
-    // 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;
 }
 
@@ -14883,14 +14876,56 @@ inline uint32_t trans_checking_pass(bubble_type* bub, kv_u_trans_t *ref, uint32_
     return 1;
 }
 
+uint32_t get_u_trans_spec_idx(kv_u_trans_t *ta, uint32_t qn, uint32_t tn, u_trans_t **r_a, uint32_t *occ, uint32_t *idx)
+{
+    if(r_a) (*r_a) = NULL; 
+    if(occ) (*occ) = 0;
+    if(idx) (*idx) = (uint32_t)-1;
+    u_trans_t *a = NULL;
+    uint32_t n, st, i;
+    a = u_trans_a(*ta, qn);
+    n = u_trans_n(*ta, qn);
+    for (st = 0, i = 1; i <= n; ++i)
+    {
+        if (i == n || a[i].tn != a[st].tn)
+        {
+            if(a[st].tn == tn)
+            {
+                if(r_a) (*r_a) = a + st;
+                if(occ) (*occ) = i - st;
+                if(idx) (*idx) = st + ((*ta).idx.a[(qn)]>>32);
+                return 1;
+            }
+            st = i;
+        }
+    }
+    return 0;
+}
+
+void interpr_hit(ha_ug_index* idx, uint64_t x, uint32_t rLen, uint32_t *uid, uint32_t *beg, uint32_t *end);
+int hic_sc_type(ha_ug_index* idx, kvec_pe_hit* hits, uint64_t k)
+{
+    uint32_t s_uid, s_beg, s_end, e_uid, e_beg, e_end, slen, elen, x = 0;
+
+    interpr_hit(idx, hits->a.a[k].s, hits->a.a[k].len>>32, &s_uid, &s_beg, &s_end); 
+    s_beg = (s_beg+s_end)>>1; slen = idx->ug->u.a[s_uid].len;
+    if(s_beg >= (slen>>1)) x+=1;
+
+    interpr_hit(idx, hits->a.a[k].e, (uint32_t)hits->a.a[k].len, &e_uid, &e_beg, &e_end); 
+    e_beg = (e_beg+e_end)>>1; elen = idx->ug->u.a[e_uid].len;
+    if(e_beg >= (elen>>1)) x+=2;
+    return x;
+}
+
 void weight_kv_u_trans(ha_ug_index* idx, kvec_pe_hit* hits, hc_links* link, bubble_type* bub, 
-kv_u_trans_t *ta, trans_idx* dis)
+kv_u_trans_t *ta, trans_idx* dis, int sc_weight)
 {
     uint64_t k, i, shif = 64 - idx->uID_bits, beg, end, t_d;
     u_trans_t *e1 = NULL, *e2 = NULL;
     long double weight;
+    double *sw = NULL;
     u_trans_t *p = NULL;
-    uint32_t is_cc;
+    uint32_t is_cc, ii1, ii2;
 
     for (i = 0, ta->idx.n = ta->n = 0; i < link->a.n; i++)
     {
@@ -14908,6 +14943,11 @@ kv_u_trans_t *ta, trans_idx* dis)
         }
     }
     kt_u_trans_t_idx(ta, idx->ug->g->n_seq);
+    if(sc_weight) {
+        k = ta->n*3;
+        MALLOC(sw, k);
+        for (i = 0; i < k; i++) sw[i] = 0;
+    }
 
 
     for (k = 0; k < hits->a.n; ++k) 
@@ -14925,15 +14965,37 @@ kv_u_trans_t *ta, trans_idx* dis)
         // if(t_d == (uint64_t)-1 && is_cc == 0) continue;
         // if(t_d == (uint64_t)-1 && !dis) continue;
 
-        get_u_trans_spec(ta, beg, end, &e1, NULL);
+        get_u_trans_spec_idx(ta, beg, end, &e1, NULL, &ii1);
-        get_u_trans_spec(ta, end, beg, &e2, NULL);
+        get_u_trans_spec_idx(ta, end, beg, &e2, NULL, &ii2);
 
         if(e1 == NULL || e2 == NULL) continue;
         weight = 1;
         if(dis) weight = get_trans_weight_advance(idx, t_d, dis);
         
-        e1->nw -= weight; e1->occ++;
+        if(sc_weight){
-        e2->nw -= weight; e2->occ++;
+            i = hic_sc_type(idx, hits, k);
+            if(i == 0){
+                e1->nw -= weight; e2->nw -= weight;
+            }
+            else{
+                i--;
+                sw[(ii1*3)+i] -= weight;
+                sw[(ii2*3)+i] -= weight;
+            }
+        } else{
+            e1->nw -= weight; e2->nw -= weight; 
+        }
+        e1->occ++; e2->occ++;
+    }    
+
+    if(sc_weight) {
+        for (i = 0; i < ta->n; ++i){
+            if(ta->a[i].nw > sw[(i*3)]) ta->a[i].nw = sw[(i*3)];
+            if(ta->a[i].nw > sw[(i*3)+1]) ta->a[i].nw = sw[(i*3)+1];
+            if(ta->a[i].nw > sw[(i*3)+2]) ta->a[i].nw = sw[(i*3)+2];
+            ta->a[i].nw *= 4/**2**/;
+        }
+        free(sw);
     }
 }
 
@@ -15450,7 +15512,7 @@ ha_ug_index* idx, bubble_type* bub, int8_t *s, mc_gg_status *sa, uint32_t ignore
 
     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);
+    weight_kv_u_trans(idx, hits, lk, bub, ta, is_comples_weight == 1? &dis : NULL, asm_opt.flag&HA_F_USKEW?0:1);
     // if(bub->round_id == bub->n_round-1)
     // {
     //     print_debug_hc_links(idx, bub, lk, ta, hits);
@@ -15686,12 +15748,9 @@ void resolve_tangles_hic(ha_ug_index *idx, bubble_type *bub, kvec_pe_hit *hits,
     uint64_t shif = 64 - idx->uID_bits, qn, tn;
     pe_hit *h_a = NULL;
     u_trans_t *p = NULL;
-
     identify_bubbles(idx->ug, bub, idx->t_ch->ir_het, &(idx->t_ch->k_trans));
-
     ta->idx.n = ta->n = 0;
     if(hits->idx.n == 0) idx_hc_links(hits, idx, NULL);
-
     for (qn = 0; qn < hits->idx.n; qn++)
     {
         h_a = hits->a.a + (hits->idx.a[qn]>>32);
@@ -15713,9 +15772,7 @@ void resolve_tangles_hic(ha_ug_index *idx, bubble_type *bub, kvec_pe_hit *hits,
             }
         }
     }
-
     radix_sort_u_trans_m(ta->a, ta->a + ta->n);
-
     for (k = 1, l = 0, m = 0; k <= ta->n; ++k)
     {
         if (k == ta->n || ta->a[l].qn != ta->a[k].qn || ta->a[l].tn != ta->a[k].tn) //same qn and tn
@@ -15737,22 +15794,20 @@ void resolve_tangles_hic(ha_ug_index *idx, bubble_type *bub, kvec_pe_hit *hits,
     ta->n = m;
     kt_u_trans_t_idx(ta, idx->ug->g->n_seq);
     
-    
     resolve_bubble_chain_by_hic(idx, ta, bub);
-
     ta->idx.n = ta->n = 0;
 }
 
 
-void print_kv_u_trans_t(kv_u_trans_t *ta)
+void print_kv_u_trans_t(kv_u_trans_t *ta, ma_ug_t* ug)
 {
     uint32_t i;
     u_trans_t *p = NULL;
     for (i = 0; i < ta->n; i++)
     {
         p = &(ta->a[i]);
-        fprintf(stderr, "q-utg%.6ul\tqs(%u)\tqe(%u)\tt-utg%.6ul\tts(%u)\tte(%u)\trev(%u)\tw(%f)\tf(%u)\n", 
+        fprintf(stderr, "q-utg%.6ul\tql(%u)\tqs(%u)\tqe(%u)\tt-utg%.6ul\ttl(%u)\tts(%u)\tte(%u)\trev(%u)\tw(%f)\tf(%u)\n", 
-        p->qn+1, p->qs, p->qe, p->tn+1, p->ts, p->te, p->rev, p->nw, p->f);
+        p->qn+1, ug->u.a[p->qn].len, p->qs, p->qe, p->tn+1, ug->u.a[p->tn].len, p->ts, p->te, p->rev, p->nw, p->f);
     }
     fprintf(stderr, "[M::%s::] \n", __func__);
 }
@@ -16191,6 +16246,61 @@ void renew_idx_para(ha_ug_index* idx, ma_ug_t* ug)
     idx->rev_mode = ((uint64_t)1) << 63;
 }
 
+uint32_t get_oe_occ(uint32_t qn, uint32_t tn, kvec_pe_hit* hits, ha_ug_index* idx)
+{
+    uint64_t shif = 64 - idx->uID_bits, occ = 0;
+    pe_hit *h_a = hits->a.a + (hits->idx.a[qn]>>32);
+    uint32_t h_occ = (uint32_t)(hits->idx.a[qn]), i;
+    for (i = 0; i < h_occ; i++) {
+        if(((h_a[i].s<<1)>>shif)==qn && ((h_a[i].e<<1)>>shif)==tn) occ++;
+    }
+    return occ;
+}
+
+void optimize_u_trans(kv_u_trans_t *ovlp, kvec_pe_hit* hits, ha_ug_index* idx)
+{
+    if(hits->idx.n == 0) idx_hc_links(hits, idx, NULL);
+    uint64_t i, m, occ;
+    u_trans_t *x = NULL, *p = NULL;
+    kv_u_trans_t k_trans; 
+    kv_init(k_trans); kv_init(k_trans.idx);
+    for (i = 0; i < ovlp->n; i++){
+        x = &(ovlp->a[i]);
+        if(x->qn > x->tn) continue;
+        if(x->f != RC_2 || x->del) continue;
+        occ = get_oe_occ(x->qn, x->tn, hits, idx) + get_oe_occ(x->tn, x->qn, hits, idx);
+        kv_pushp(u_trans_t, k_trans, &p);
+        (*p) = (*x); p->nw = (x->nw*(1-(((double)(occ<<1))/((double)(hits->occ.a[x->qn]+hits->occ.a[x->tn])))));
+        if(p->nw < 0) fprintf(stderr, "ERROR-nw\n");
+        if(p->nw == 0) p->nw = x->nw*0.005;
+        if(p->nw == 0) {
+            k_trans.n--;
+        }
+        else {
+            kv_pushp(u_trans_t, k_trans, &p);
+            (*p) = k_trans.a[k_trans.n-2];
+            p->qn = k_trans.a[k_trans.n-2].tn; p->qs = k_trans.a[k_trans.n-2].ts; p->qe = k_trans.a[k_trans.n-2].te;
+            p->tn = k_trans.a[k_trans.n-2].qn; p->ts = k_trans.a[k_trans.n-2].qs; p->te = k_trans.a[k_trans.n-2].qe;
+        }
+    }
+    kt_u_trans_t_idx(&k_trans, idx->ug->g->n_seq);
+    mc_solve(NULL, NULL, &k_trans, idx->ug, idx->read_g, 0.8, R_INF.trio_flag, 1, NULL, 1, NULL, NULL, 1);
+    for (i = m = 0; i < ovlp->n; i++){
+        x = &(ovlp->a[i]);
+        if(x->del) continue;
+        if(x->f == RC_2){
+            get_u_trans_spec(&k_trans, x->qn, x->tn, &p, NULL);
+            if(!p) continue;
+        }
+        ovlp->a[m++] = ovlp->a[i];
+    }
+    ovlp->n = m;
+    kt_u_trans_t_idx(ovlp, idx->ug->g->n_seq);
+    free(hits->idx.a); hits->idx.a = NULL; hits->idx.n = hits->idx.m = 0;
+    free(hits->occ.a); hits->occ.a = NULL; hits->occ.n = hits->occ.m = 0;
+    kv_destroy(k_trans); kv_destroy(k_trans.idx);
+}
+
 int hic_short_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx, ug_opt_t *opt, kvec_pe_hit **rhits)
 {
     double index_time = yak_realtime();
@@ -16209,15 +16319,16 @@ int hic_short_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx, ug_o
         alignment_worker_pipeline(&sl, fn1, fn2);
         write_hc_hits(&sl.hits, idx->ug, asm_opt.output_file_name);
     }
+    sl.hits.uID_bits = idx->uID_bits; sl.hits.pos_mode = idx->pos_mode;
+    optimize_u_trans(&(idx->t_ch->k_trans), &sl.hits, idx);
     filter_kv_u_trans_t(&(idx->t_ch->k_trans), idx->ug, 0.5);
+    // print_kv_u_trans_t(&(idx->t_ch->k_trans), idx->ug);
     // flter_by_cov(idx, &sl.hits, 2);
     // update_hits(idx, &sl.hits, idx->t_ch->is_r_het);
     ///debug_hc_hits_v14(&sl.hits, asm_opt.output_file_name, sl.idx);
     ////dedup_hits(&(sl.hits), sl.idx);   
     ///write_hc_hits_v14(&sl.hits, asm_opt.output_file_name);
     
-
-    sl.hits.uID_bits = idx->uID_bits; sl.hits.pos_mode = idx->pos_mode;
     if(asm_opt.misjoin_len > 0)
     {
         update_switch_unitig(idx->ug, idx->read_g, &(sl.hits), &(idx->t_ch->k_trans), 10, 20, asm_opt.misjoin_len, 0.15);
@@ -16255,7 +16366,7 @@ int hic_short_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx, ug_o
         // 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, 1, /**&bub**/NULL, &(idx->t_ch->k_trans));
+        (bub.round_id == 0? 1 : 0), s->s, 1, /**&bub**/NULL, &(idx->t_ch->k_trans), 0);
         /*******************************for debug************************************/
         label_unitigs_sm(s->s, NULL, idx->ug);
 
@@ -16321,6 +16432,49 @@ int hic_short_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx, ug_o
     return 1;
 }
 
+spg_t *hic_short_pre_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx, ug_opt_t *opt, kvec_pe_hit **rhits)
+{
+    // double index_time = yak_realtime();
+    sldat_t sl;
+    sl.idx = idx;
+    sl.t_ch = idx->t_ch;
+    sl.chunk_size = 20000000;
+    sl.n_thread = asm_opt.thread_num;
+    sl.total_base = sl.total_pair = 0;
+    idx->hap_cnt = asm_opt.hap_occ;
+    kv_init(sl.hits.a); kv_init(sl.hits.idx); kv_init(sl.hits.occ);
+
+    if(!load_hc_hits(&sl.hits, idx->ug, asm_opt.output_file_name))
+    {
+        alignment_worker_pipeline(&sl, fn1, fn2);
+        fprintf(stderr, "sb0sb\n");
+        write_hc_hits(&sl.hits, idx->ug, asm_opt.output_file_name);
+        fprintf(stderr, "sb1sb\n");
+    }
+    sl.hits.uID_bits = idx->uID_bits; sl.hits.pos_mode = idx->pos_mode;
+    kv_u_trans_t k_trans; 
+    kv_init(k_trans); kv_init(k_trans.idx);
+    bubble_type bub; 
+    memset(&bub, 0, sizeof(bubble_type));
+    bub.round_id = 0; bub.n_round = asm_opt.n_weight;
+    fprintf(stderr, "sb2sb\n");
+    resolve_tangles_hic(idx, &bub, &sl.hits, &k_trans);
+    fprintf(stderr, "sb3sb\n");
+    spg_t *scg = horder_utg(&sl.hits, idx->uID_bits, idx->pos_mode, idx->read_g, idx->ug, &bub, opt);
+    fprintf(stderr, "sb4sb\n");
+    if(rhits){
+        CALLOC(*rhits, 1);
+        (**rhits) = sl.hits;
+        sl.hits.a.a = NULL;
+        sl.hits.idx.a = NULL;
+        sl.hits.occ.a = NULL;
+    }
+    kv_destroy(sl.hits.a); kv_destroy(sl.hits.idx); kv_destroy(sl.hits.occ);
+    kv_destroy(k_trans); kv_destroy(k_trans.idx);
+    destory_bubbles(&bub);
+    return scg;
+}
+
 int load_psg_t(psg_t **sg, const char *fn)
 {
     uint64_t flag = 0;
@@ -16433,7 +16587,7 @@ int hic_short_align_poy(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx,
         // }
         /*******************************for debug************************************/
     }
-    write_psg_t(s, asm_opt.output_file_name);
+    if((asm_opt.flag & HA_F_VERBOSE_GFA)) write_psg_t(s, asm_opt.output_file_name);
 
     skip_flipping:
     verbose_het_stat(&bub);
@@ -16488,6 +16642,19 @@ void hic_analysis(ma_ug_t *ug, asg_t* read_g, trans_chain* t_ch, ug_opt_t *opt,
     destory_hc_pt_index(ug_index);free(ug_index);
 }
 
+spg_t *hic_pre_analysis(ma_ug_t *ug, asg_t* read_g, trans_chain* t_ch, ug_opt_t *opt, kvec_pe_hit **rhits)
+{
+    ug_index = NULL;
+    int exist = (asm_opt.load_index_from_disk? 
+                    load_hc_pt_index(&ug_index, ug, asm_opt.output_file_name) : 0);
+    if(exist == 0) ug_index = build_unitig_index(ug, asm_opt.hic_mer_length, asm_opt.hap_occ, 0, asm_opt.thread_num);
+    if(exist == 0) write_hc_pt_index(ug_index, asm_opt.output_file_name);
+    ug_index->ug = ug;
+    ug_index->read_g = read_g;
+    ug_index->t_ch = t_ch;
+    return hic_short_pre_align(asm_opt.hic_reads[0], asm_opt.hic_reads[1], ug_index, opt, rhits);    
+}
+
 
 void init_ug_idx(ma_ug_t *ug, uint64_t k, uint64_t up_bound, uint64_t low_bound, uint64_t build_idx)
 {

hic.h

@@ -75,6 +75,7 @@ typedef struct{
 }pdq;
 
 
+
 #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)
@@ -106,5 +107,5 @@ long long *dis);
 void set_utg_by_dis(uint32_t v, pdq* pq, asg_t *g, kvec_t_u32_warp *res, uint32_t dis);
 void dedup_hits(kvec_pe_hit* hits, uint64_t is_dup);
 void hic_analysis(ma_ug_t *ug, asg_t* read_g, trans_chain* t_ch, ug_opt_t *opt, uint32_t is_poy, kvec_pe_hit **rhits);
-
+spg_t *hic_pre_analysis(ma_ug_t *ug, asg_t* read_g, trans_chain* t_ch, ug_opt_t *opt, kvec_pe_hit **rhits);
 #endif

horder.cpp

@@ -2178,9 +2178,9 @@ uint64_t rs, uint64_t re, uint64_t limit_s, uint64_t limit_e, int unique_only)
     return cnt;
 }
 
-void detect_lowNs(kvec_pe_hit *hit, uint64_t sHit, uint64_t eHit, kvec_t_u64_warp *b, 
+int detect_lowNs(kvec_pe_hit *hit, uint64_t sHit, uint64_t eHit, kvec_t_u64_warp *b, 
-h_cov_t *Np, uint64_t len, uint64_t cutoff_s, uint64_t cutoff_e, h_covs *res,
+h_cov_t *Np, uint64_t len, uint64_t cutoff_s, uint64_t cutoff_e, uint64_t force_cutoff, uint64_t force_cutoff_cov,
-h_covs *cov_buf, h_covs *b_points, uint64_t local_bound, int unique_only)
+h_covs *res, h_covs *cov_buf, h_covs *b_points, uint64_t local_bound, int unique_only)
 {
     uint64_t cov_hic, cov_utg, cov_ava, i, p0s, p0e, p1s, p1e, span_s, span_e, cutoff, bs, be, occ = 0;
     uint64_t sPos, ePos, min_cutoff;
@@ -2219,6 +2219,18 @@ h_covs *cov_buf, h_covs *b_points, uint64_t local_bound, int unique_only)
     radix_sort_ho64(b->a.a, b->a.a+b->a.n);
     cov_utg = get_hic_cov_interval(b->a.a, b->a.n, 1, NULL, NULL, NULL);
     cov_ava = (cov_utg? cov_hic/cov_utg:0);
+    if(force_cutoff != (uint64_t)-1 || force_cutoff_cov != (uint64_t)-1)
+    {
+        min_cutoff = get_sub_cov(hit, sHit, eHit, len, Np->s, Np->e, sPos, ePos, unique_only);
+        if((force_cutoff != (uint64_t)-1 && min_cutoff <= (cov_ava/force_cutoff)) || 
+                (force_cutoff_cov != (uint64_t)-1 && min_cutoff <= force_cutoff_cov))
+        {
+            kv_pushp(h_cov_t, *b_points, &p);
+            p->s = get_hit_suid(*hit, sHit); p->e = Np->dp; p->dp = 0;
+            return 1;
+        }
+    }
+    
     ///if cov_ava == 0, do nothing or break?
     /*******************************for debug************************************/
     // fprintf(stderr, "\n[M::%s::] utg%.6lul, ulen: %lu, # hic hits: %lu, map cov: %lu, utg cov: %lu, average: %lu\n", 
@@ -2269,24 +2281,29 @@ h_covs *cov_buf, h_covs *b_points, uint64_t local_bound, int unique_only)
         {
             kv_pushp(h_cov_t, *b_points, &p);
             p->s = get_hit_suid(*hit, sHit); p->e = Np->dp; p->dp = 0;
+            return 1;
             /*******************************for debug************************************/
             // fprintf(stderr, "consensus_break-rid: %lu\n", p->e);
             /*******************************for debug************************************/
         }
     }
+    return 0;
 }
 
-uint64_t break_scaffold(horder_t *h, uint64_t cutoff_s, uint64_t cutoff_e, uint64_t local_bound, int unique_only)
+uint64_t break_scaffold(horder_t *h, uint64_t cutoff_s, uint64_t cutoff_e, uint64_t force_cutoff, uint64_t force_cutoff_cov,
+uint64_t local_bound, int unique_only, h_covs *r_b_points)
 {
     uint64_t k, l, i, ulen;
     kvec_t_u64_warp b; kv_init(b.a);
     h_covs cov_buf; kv_init(cov_buf);
     h_covs res; kv_init(res);
-    h_covs b_points; kv_init(b_points);
+    h_covs *b_points = NULL; 
+    if(r_b_points) b_points = r_b_points;
+    else CALLOC(b_points, 1);
+    b_points->n = 0;
     h_covs Ns; kv_init(Ns);
     ma_ug_t *ug = h->ug;
     kvec_pe_hit *hits = &(h->u_hits);
-    b_points.n = 0; 
     for (k = 1, l = 0; k <= hits->a.n; ++k) 
     {   
         if (k == hits->a.n || (get_hit_suid(*hits, k) != get_hit_suid(*hits, l))) 
@@ -2302,8 +2319,10 @@ uint64_t break_scaffold(horder_t *h, uint64_t cutoff_s, uint64_t cutoff_e, uint6
                 {
                     for (i = 0; i < Ns.n; i++)
                     {
-                        detect_lowNs(hits, l, k, &b, &(Ns.a[i]), ulen, cutoff_s, cutoff_e, 
+                        if(detect_lowNs(hits, l, k, &b, &(Ns.a[i]), ulen, cutoff_s, cutoff_e, force_cutoff,
-                        &res, &cov_buf, &b_points, local_bound, unique_only);
+                        force_cutoff_cov, &res, &cov_buf, b_points, local_bound, unique_only) && r_b_points){
+                            b_points->a[b_points->n-1].dp = i;
+                        }
                     }
                 }
             }  
@@ -2311,15 +2330,68 @@ uint64_t break_scaffold(horder_t *h, uint64_t cutoff_s, uint64_t cutoff_e, uint6
         }
     }
 
-    break_utg_horder(h, &b_points);
+    if(!r_b_points){
+        break_utg_horder(h, b_points);
+        kv_destroy(*b_points);
+    }
     
     kv_destroy(b.a);
     kv_destroy(cov_buf);
     kv_destroy(res);
-    kv_destroy(b_points);
     kv_destroy(Ns);
 
-    return b_points.n;
+    return b_points->n;
+}
+
+uint64_t break_scaffold_mean(horder_t *h, uint64_t cutoff_s, uint64_t cutoff_e, uint64_t force_cutoff, uint64_t force_cutoff_cov,
+uint64_t local_bound, int unique_only, h_covs *r_b_points)
+{
+    uint64_t k, l, i, ulen;
+    kvec_t_u64_warp b; kv_init(b.a);
+    h_covs cov_buf; kv_init(cov_buf);
+    h_covs res; kv_init(res);
+    h_covs *b_points = NULL; 
+    if(r_b_points) b_points = r_b_points;
+    else CALLOC(b_points, 1);
+    b_points->n = 0;
+    h_covs Ns; kv_init(Ns);
+    ma_ug_t *ug = h->ug;
+    kvec_pe_hit *hits = &(h->u_hits);
+    for (k = 1, l = 0; k <= hits->a.n; ++k) 
+    {   
+        if (k == hits->a.n || (get_hit_suid(*hits, k) != get_hit_suid(*hits, l))) 
+        {
+            ulen = ug->u.a[get_hit_suid(*hits, l)].len;
+            Ns.n = 0;
+            // if(ulen >= BREAK_THRES)
+            {
+                get_Ns(&(ug->u.a[get_hit_suid(*hits, l)]), &Ns);
+                if(Ns.n)
+                {
+                    for (i = 0; i < Ns.n; i++)
+                    {
+                        if(detect_lowNs(hits, l, k, &b, &(Ns.a[i]), ulen, cutoff_s, cutoff_e, force_cutoff,
+                        force_cutoff_cov, &res, &cov_buf, b_points, local_bound, unique_only) && r_b_points){
+                            b_points->a[b_points->n-1].dp = i;
+                        }
+                    }
+                }
+            }  
+            l = k;
+        }
+    }
+
+    if(!r_b_points){
+        break_utg_horder(h, b_points);
+        kv_destroy(*b_points);
+    }
+    
+    kv_destroy(b.a);
+    kv_destroy(cov_buf);
+    kv_destroy(res);
+    kv_destroy(Ns);
+
+    return b_points->n;
 }
 
 
@@ -2729,7 +2801,7 @@ void update_scg(horder_t *h, trans_col_t *t_idx)
         h->sg.g->seq[i].ez[0] = ug->u.a[i].len>>1;
         h->sg.g->seq[i].ez[1] = ug->u.a[i].len - (ug->u.a[i].len>>1);
     }
-    idx = build_interval_idx(hits, ug);
+    idx = build_interval_idx(hits, ug);///idx is used to get density
     
     for (i = 0, e.n = 0; i < hits->a.n; i++)
     {
@@ -3397,10 +3469,12 @@ uint32_t get_sl_occ(sc_lay_t *sl)
 
 void refine_layout(horder_t *h, sc_lay_t *sl, uint8_t *vis)
 {
-    uint32_t k, m/**, max_utg, max_sc**/;
+    uint32_t k/**m, max_utg, max_sc**/;
     lay_t *p = NULL;
     uint8_t *sgv = NULL; MALLOC(sgv, sl->n);
     double *w = NULL; MALLOC(w, sl->n);
+    
+    /**
     uint32_t *idx = NULL; MALLOC(idx, h->sg.g->n_seq);
     memset(idx, -1, sizeof(uint32_t)*h->sg.g->n_seq);
 
@@ -3412,8 +3486,6 @@ void refine_layout(horder_t *h, sc_lay_t *sl, uint8_t *vis)
             idx[p->a[m]>>1] = k;
         }
     }
-
-    /**
     while (get_max_anchor(h, sl, vis, w, sgv, idx, &max_utg, &max_sc))
     {
         
@@ -3422,6 +3494,7 @@ void refine_layout(horder_t *h, sc_lay_t *sl, uint8_t *vis)
         vis[max_utg<<1] = vis[(max_utg<<1)+1] = 1;
         idx[max_utg] = max_sc;
     }
+    free(idx); 
     **/
 
     for (k = 0; k < h->sg.g->n_seq; k++)
@@ -3434,8 +3507,7 @@ void refine_layout(horder_t *h, sc_lay_t *sl, uint8_t *vis)
         vis[(k<<1)] = vis[(k<<1)+1] = 1;
     }
 
-    
+    free(w); free(sgv);
-    free(w); free(idx); free(sgv);
 }
 
 
@@ -3591,7 +3663,7 @@ void update_avoids(horder_t *h, sc_lay_t *sl)
     h->avoid.n = m;
 }
 
-void update_ug_by_layout(horder_t *h, sc_lay_t *sl)
+void update_ug_by_layout(horder_t *h, sc_lay_t *sl, ma_ug_t* i_ug)
 {
     uint32_t i;
     lay_t *p = NULL;
@@ -3603,7 +3675,7 @@ void update_ug_by_layout(horder_t *h, sc_lay_t *sl)
     {
         p = &(sl->a[i]);
         kv_pushp(ma_utg_t, sug->u, &pu);
-        generate_scaffold(pu, p, h->ug, h->r_g);
+        generate_scaffold(pu, p, i_ug?i_ug:h->ug, h->r_g);
     }
     ma_ug_destroy(h->ug);
     h->ug = sug;
@@ -3677,7 +3749,14 @@ void get_long_switch_scaffolds(horder_t *h, sc_lay_t *sl, osg_t *lg)
     kv_destroy(idx);
 }
 
-void layout_scg(horder_t *h, double nw_thres, uint32_t occ_thres)
+void destory_sc_lay_t(sc_lay_t *sl)
+{
+    uint32_t i;
+    for (i = 0; i < sl->n; i++) kv_destroy(sl->a[i]);
+    kv_destroy(*sl);
+}
+
+void layout_scg(horder_t *h, double nw_thres, uint32_t occ_thres, sc_lay_t *r_sl)
 {
     uint32_t k;
     osg_arc_t *p = NULL, *lp = NULL;
@@ -3708,17 +3787,23 @@ void layout_scg(horder_t *h, double nw_thres, uint32_t occ_thres)
 
     get_backbone_layout(h, &sl, lg, vis);
 
-    get_long_switch_scaffolds(h, &sl, lg);
+    // get_long_switch_scaffolds(h, &sl, lg);
     
     refine_layout(h, &sl, vis);
     
     // print_N50_layout(h->ug, &sl);
 
-    update_ug_by_layout(h, &sl);
+    update_ug_by_layout(h, &sl, NULL);
 
     print_N50(h->ug);
 
-    kv_destroy(sl);
+    if(r_sl){
+        r_sl->a = sl.a; sl.a = NULL;
+        r_sl->m = sl.m; sl.m = 0;
+        r_sl->n = sl.n; sl.n = 0;
+
+    }
+    destory_sc_lay_t(&sl);
     osg_destroy(lg);
     free(vis);
 }
@@ -3729,7 +3814,7 @@ void renew_scaffold(horder_t *h)
     while (1)
     {
         update_u_hits(&(h->u_hits), &(h->r_hits), h->ug, h->r_g);
-        if(!break_scaffold(h, 5, 15, 2500000, 1)) break;
+        if(!break_scaffold(h, 5, 15, (uint64_t)-1, (uint64_t)-1, 2500000, 1, NULL)) break;
         print_N50(h->ug);
     }
     fprintf(stderr, "[M::%s::%.3f] \n", __func__, yak_realtime()-index_time);
@@ -3789,7 +3874,7 @@ void scaffold_hap(horder_t *h, ug_opt_t *opt, trans_col_t *t_idx, uint32_t round
     {
         update_u_hits(&(h->u_hits), &(h->r_hits), h->ug, h->r_g);
         update_scg(h, t_idx);
-        layout_scg(h, 1.001, 19);
+        layout_scg(h, 1.001, 19, NULL);
         renew_scaffold(h);
     }
     
@@ -3832,7 +3917,7 @@ void scaffold_ug(horder_t *h, ma_ug_t *ug, ug_opt_t *opt, uint32_t round, char *
         fprintf(stderr, "[M::%s::]**i->%u**\n", __func__, i);
         update_scg(h, NULL);
         fprintf(stderr, "[M::%s::]***i->%u***\n", __func__, i);
-        layout_scg(h, 1.001, 19);
+        layout_scg(h, 1.001, 19, NULL);
         fprintf(stderr, "[M::%s::]****i->%u****\n", __func__, i);
         renew_scaffold(h);
         fprintf(stderr, "[M::%s::]*****i->%u*****\n", __func__, i);
@@ -3887,7 +3972,7 @@ asg_t *i_rg, ma_ug_t* i_ug, bubble_type* bub, kv_u_trans_t *ref, ug_opt_t *opt,
     {
         update_u_hits(&(h->u_hits), &(h->r_hits), h->ug, h->r_g);
         update_scg(h, t_idx);
-        layout_scg(h, 1.001, 19);
+        layout_scg(h, 1.001, 19, NULL);
         renew_scaffold(h);
     }
     
@@ -3902,6 +3987,113 @@ asg_t *i_rg, ma_ug_t* i_ug, bubble_type* bub, kv_u_trans_t *ref, ug_opt_t *opt,
     return h;
 }
 
+void cpy_u_hits(kvec_pe_hit *u_hits, kvec_pe_hit *i_hits, uint32_t u_n)
+{
+    uint64_t i;
+    memset(u_hits, 0, sizeof(kvec_pe_hit));
+    u_hits->pos_mode = i_hits->pos_mode;
+    u_hits->uID_bits = i_hits->uID_bits;
+    u_hits->a.n = u_hits->a.m = i_hits->a.n;
+    MALLOC(u_hits->a.a, u_hits->a.n);
+    memcpy(u_hits->a.a, i_hits->a.a, u_hits->a.n*sizeof(pe_hit));
+    for (i = 0; i < u_hits->a.n; i++) u_hits->a.a[i].id = 1;
+    idx_hits(u_hits, u_n);
+}
+
+void update_sc_lay(sc_lay_t *sl, h_covs *b)
+{
+    uint64_t k, l, i, pidx, cidx;
+    lay_t *s = NULL;
+    lay_t *p = NULL;
+    for (k = 1, l = 0; k <= b->n; ++k) 
+    {   
+        if (k == b->n || (b->a[k].s != b->a[l].s)) 
+        {
+            s = &(sl->a[b->a[l].s]);
+            for (i = l, pidx = 0; i < k; i++){
+                cidx = b->a[i].dp;
+                kv_pushp(lay_t, *sl, &p);
+                kv_init(*p);
+                p->n = p->m = (cidx - pidx + 1)<<1;
+                MALLOC(p->a, p->n);
+                mempcpy(p->a, s->a + (pidx<<1), sizeof(*(p->a))*p->n);
+                pidx = cidx + 1;
+            }
+
+            if(pidx >= (s->n>>1)) fprintf(stderr, "ERROR-update\n");
+            cidx = (s->n>>1)-1;
+            kv_pushp(lay_t, *sl, &p);
+            kv_init(*p);
+            p->n = p->m = (cidx - pidx + 1)<<1;
+            MALLOC(p->a, p->n);
+            mempcpy(p->a, s->a + (pidx<<1), sizeof(*(p->a))*p->n);
+            free(s->a); s->n = s->m = 0;
+            l = k;
+        }
+    }
+
+    for (i = k = 0; i < sl->n; i++){
+        if(!sl->a[i].a) continue;
+        sl->a[k] = sl->a[i];
+        sl->a[i].a = NULL; sl->a[i].n = sl->a[i].m = 0;
+        if(sl->a[k].n == 2){
+            sl->a[k].a[0] >>= 1; sl->a[k].a[0] <<= 1;
+            sl->a[k].a[1] >>= 1; sl->a[k].a[1] <<= 1; sl->a[k].a[1]++;
+        }
+        k++;
+    }
+    sl->n = k;
+}
+
+void renew_scaffold_utg(horder_t *h, sc_lay_t *sl, ma_ug_t* i_ug)
+{
+    double index_time = yak_realtime();
+    h_covs b; kv_init(b);
+    while (1)
+    {
+        update_u_hits(&(h->u_hits), &(h->r_hits), h->ug, h->r_g);
+        if(!break_scaffold(h, 5, 15, 15, 25, 2500000, 1, &b)) break;
+        update_sc_lay(sl, &b);
+        update_ug_by_layout(h, sl, i_ug);
+        print_N50(h->ug);
+    }
+    fprintf(stderr, "[M::%s::%.3f] \n", __func__, yak_realtime()-index_time);
+    kv_destroy(b);
+}
+
+spg_t *scf_g(sc_lay_t *sl, ma_ug_t* ug)
+{
+    spg_t *scg = NULL; CALLOC(scg, 1); scg->ug = ug;    
+    uint64_t i, k;
+    lay_t *p = NULL;
+    for (i = 0; i < sl->n; i++){
+        p = &(sl->a[i]);
+        kv_push(uint64_t, scg->idx, (uint64_t)scg->dst.n << 32 | (p->n>>1));
+        for (k = 0; k < p->n; k+=2) kv_push(uint32_t, scg->dst, p->a[k]);        
+    }
+    return scg;
+}
+
+spg_t *horder_utg(kvec_pe_hit *i_hits, uint64_t i_hits_uid_bits, uint64_t i_hits_pos_mode, 
+asg_t *i_rg, ma_ug_t* i_ug, bubble_type* bub, ug_opt_t *opt)
+{
+    horder_t *h = NULL; CALLOC(h, 1);
+    sc_lay_t sl; kv_init(sl);
+    get_r_hits(i_hits, &(h->r_hits), i_rg, i_ug, bub, i_hits_uid_bits, i_hits_pos_mode);
+    h->r_g = copy_read_graph(i_rg);
+    horder_clean_sg_by_utg(h->r_g, i_ug);
+    h->ug = copy_untig_graph(i_ug); asg_destroy(h->ug->g); h->ug->g = NULL;
+    cpy_u_hits(&(h->u_hits), i_hits, h->ug->u.n);
+
+    update_scg(h, NULL);
+    layout_scg(h, ((double)1)/((double)0.75), 19, &sl);
+    renew_scaffold_utg(h, &sl, i_ug);
+    
+    spg_t *scg = scf_g(&sl, i_ug);
+    destory_sc_lay_t(&sl);
+    destory_horder_t(&h);
+    return scg;
+}
 
 
 void ha_aware_order(kvec_pe_hit *r_hits, asg_t *rg, ma_ug_t *ug_fa, ma_ug_t *ug_mo, kv_u_trans_t *ref, 

horder.h

@@ -67,4 +67,6 @@ kvec_pe_hit *get_r_hits_order(kvec_pe_hit *uhits, uint64_t hits_uid_bits, uint64
 asg_t *rg, ma_ug_t* ug, bubble_type* bub);
 void ha_aware_order(kvec_pe_hit *r_hits, asg_t *rg, ma_ug_t *ug_fa, ma_ug_t *ug_mo, kv_u_trans_t *ref, 
 ug_opt_t *opt, uint32_t round);
+spg_t *horder_utg(kvec_pe_hit *i_hits, uint64_t i_hits_uid_bits, uint64_t i_hits_pos_mode, 
+asg_t *i_rg, ma_ug_t* i_ug, bubble_type* bub, ug_opt_t *opt);
 #endif

htab.cpp

@@ -38,75 +38,6 @@ void *ha_flt_tab_hp;
 ha_pt_t *ha_idx_hp;
 void *ha_ct_table;
 
-#define MZ_FUNC_INIT(sf, HType) \
-static inline void sf##_init_kuf(pl_data_t *p, st_data_t *s){\
-	int i, n_pre = 1<<p->opt->pre, m;\
-	/**allocate the k-mer buffer**/\
-	CALLOC(s->buf, n_pre);\
-	m = (int)(s->nk * 1.2 / n_pre) + 1;\
-	/**pre-allocate memory for each of 4096 buffer**/\
-	for (i = 0; i < n_pre; ++i) {\
-		s->buf[i].m = m;\
-		/**for 0-th counting, p->pt = NULL**/\
-		if (p->pt && !(p->flag&HAF_COUNT_REFINE)) MALLOC(s->buf[i].b_##sf, m);\
-		else MALLOC(s->buf[i].a, m);\
-	}\
-}\
-static inline void sf##_destory_kuf(pl_data_t *p, st_data_t *s, int n){\
-	int i;\
-	uint64_t n_ins = 0;\
-	/**n_ins is number of distinct k-mers**/\
-	for (i = 0; i < n; ++i) {\
-		n_ins += s->buf[i].n_ins;\
-		if (p->pt && !(p->flag&HAF_COUNT_REFINE)) free(s->buf[i].b_##sf);\
-		else free(s->buf[i].a);\
-	}\
-	if (p->ct) p->ct->tot += n_ins, p->ct->bs += s->sum_len;\
-	if (p->pt) p->pt->tot_pos += n_ins;\
-	free(s->buf);\
-	/**#if 0\
-	fprintf(stderr, "[M::%s::%.3f*%.2f] processed %ld sequences; %ld %s in the hash table\n", __func__,\
-			yak_realtime(), yak_cpu_usage(), (long)s->n_seq0 + s->n_seq,\
-			(long)(p->pt? p->pt->tot_pos : p->ct->tot), p->pt? "positions" : "distinct k-mers");\
-	#endif**/\
-	free(s);\
-}\
-static inline void sf##_pt_insert_buf(ch_buf_t *buf, int p, const HType *y){\
-	/**assign minimizer to one of 4096 bins by low 12 bits**/\
-	int pre = y->x & ((1<<p) - 1);\
-	ch_buf_t *b = &buf[pre];\
-	if (b->n == b->m) {\
-		b->m = b->m < 8? 8 : b->m + (b->m>>1);\
-		REALLOC(b->b_##sf, b->m);\
-	}\
-	b->b_##sf[b->n++] = *y;\
-}\
-static inline void sf##_mselect(pl_data_t *p, st_data_t *s){\
-	int i; uint32_t j;\
-	/**s->n_seq is how many reads at this buffer**/\
-	/**s->mz && s->mz_buf are lists of minimzer vectors**/\
-	CALLOC(s->sf, s->n_seq), CALLOC(s->sf##_buf, p->opt->n_thread), CALLOC(s->mt, p->opt->n_thread);\
-	/**calculate minimzers for each read, each read corresponds to one thread**/\
-	kt_for(p->opt->n_thread, worker_for_mz, s, s->n_seq);\
-	for (i = 0; i < p->opt->n_thread; ++i) free(s->mt[i].a), free(s->sf##_buf[i].a);\
-	free(s->mt), free(s->sf##_buf);\
-	/**insert minimizers**/\
-	if (p->pt && !(p->flag&HAF_COUNT_REFINE)) {/**insert whole minimizer**/\
-		for (i = 0; i < s->n_seq; ++i)\
-			for (j = 0; j < s->sf[i].n; ++j)\
-				sf##_pt_insert_buf(s->buf, p->opt->pre, &s->sf[i].a[j]);\
-	} else {/**just insert the hash key of minimizer**/\
-		for (i = 0; i < s->n_seq; ++i)\
-			for (j = 0; j < s->sf[i].n; ++j)\
-				ct_insert_buf(s->buf, p->opt->pre, s->sf[i].a[j].x);\
-	}\
-	for (i = 0; i < s->n_seq; ++i) {\
-		p->n_mz += s->sf[i].n;\
-		free(s->sf[i].a);\
-		if (!p->is_store) free(s->seq[i]);\
-	}\
-	free(s->sf);}
-
 /***************************
  * Yak specific parameters *
  ***************************/
@@ -601,77 +532,7 @@ const int ha_pt_cnt(const ha_pt_t *h, uint64_t hash)
 /**********************************
  * Buffer for counting all k-mers *
  **********************************/
-
-typedef struct {
-	int n, m;
-	uint64_t n_ins;
-	uint64_t *a;
-	ha_mz1_t *b_mz;
-	ha_mzl_t *b_mzl;
-} ch_buf_t;
-
-///p = 12
-static inline void ct_insert_buf(ch_buf_t *buf, int p, uint64_t y) // insert a k-mer $y to a linear buffer
-{
-	///assign k-mer to one of the 4096 bins
-	///using low 12 bits for assigning
-	///so all elements at b have the same low 12 bits
-	int pre = y & ((1<<p) - 1);
-	ch_buf_t *b = &buf[pre];
-	if (b->n == b->m) {
-		b->m = b->m < 8? 8 : b->m + (b->m>>1);
-		REALLOC(b->a, b->m);
-	}
-	b->a[b->n++] = y;
-}
-
-///buf is the read block, k is the k-mer length, p = 12, len is the read length, seq is the read
-static void count_seq_buf(ch_buf_t *buf, int k, int p, int len, const char *seq) // insert k-mers in $seq to linear buffer $buf
-{
-	int i, l;
-	uint64_t x[4], mask = (1ULL<<k) - 1, shift = k - 1;
-	for (i = l = 0, x[0] = x[1] = x[2] = x[3] = 0; i < len; ++i) {
-		int c = seq_nt4_table[(uint8_t)seq[i]];
-		///c = 00, 01, 10, 11
-		if (c < 4) { // not an "N" base
-			///x[0] & x[1] are the forward k-mer
-			///x[2] & x[3] are the reverse complementary k-mer
-			x[0] = (x[0] << 1 | (c&1))  & mask;
-			x[1] = (x[1] << 1 | (c>>1)) & mask;
-			x[2] = x[2] >> 1 | (uint64_t)(1 - (c&1))  << shift;
-			x[3] = x[3] >> 1 | (uint64_t)(1 - (c>>1)) << shift;
-			if (++l >= k)
-				ct_insert_buf(buf, p, yak_hash_long(x));
-		} else l = 0, x[0] = x[1] = x[2] = x[3] = 0; // if there is an "N", restart
-	}
-}
-
-static void count_seq_buf_HPC(ch_buf_t *buf, int k, int p, int len, const char *seq) // insert k-mers in $seq to linear buffer $buf
-{
-	int i, l, last = -1;
-	uint64_t x[4], mask = (1ULL<<k) - 1, shift = k - 1;
-	for (i = l = 0, x[0] = x[1] = x[2] = x[3] = 0; i < len; ++i) {
-		int c = seq_nt4_table[(uint8_t)seq[i]];
-		if (c < 4) { // not an "N" base
-			if (c != last) {
-				x[0] = (x[0] << 1 | (c&1))  & mask;
-				x[1] = (x[1] << 1 | (c>>1)) & mask;
-				x[2] = x[2] >> 1 | (uint64_t)(1 - (c&1))  << shift;
-				x[3] = x[3] >> 1 | (uint64_t)(1 - (c>>1)) << shift;
-				if (++l >= k)
-					ct_insert_buf(buf, p, yak_hash_long(x));
-				last = c;
-			}
-		} else l = 0, last = -1, x[0] = x[1] = x[2] = x[3] = 0; // if there is an "N", restart
-	}
-}
-
-/******************
- * K-mer counting *
- ******************/
-
 KSEQ_INIT(gzFile, gzread)
-
 #define HAF_COUNT_EXACT  0x1
 #define HAF_COUNT_ALL    0x2
 #define HAF_RS_WRITE_LEN 0x4
@@ -696,178 +557,285 @@ typedef struct { // global data structure for kt_pipeline()
 	const ma_utg_v *us_in;
 } pl_data_t;
 
-typedef struct { // data structure for each step in kt_pipeline()
+#define MZ_TEST_INIT(sf, HType, VType, IType, Ia) \
-	pl_data_t *p;
+typedef struct {int n, m; uint64_t n_ins; uint64_t *a; HType *b;} sf##_ch_buf_t;\
-	uint64_t n_seq0;  ///the start index of current buffer block at R_INF
+static inline void sf##_ct_insert_buf(sf##_ch_buf_t *buf, int p, uint64_t y) /** insert a k-mer $y to a linear buffer**/\
-	///sum_len = total bases, nk = number of k-mers
+{\
-	int n_seq, m_seq, sum_len, nk, uq;
+	/**assign k-mer to one of the 4096 bins**/\
-	int *len;
+	/**using low 12 bits for assigning**/\
-	char **seq;
+	/**so all elements at b have the same low 12 bits**/\
-	ha_mz1_v *mz_buf;
+	int pre = y & ((1<<p) - 1);\
-	ha_mz1_v *mz;
+	sf##_ch_buf_t *b = &buf[pre];\
-	ha_mzl_v *mzl_buf;
+	if (b->n == b->m) {\
-	ha_mzl_v *mzl;
+		b->m = b->m < 8? 8 : b->m + (b->m>>1);\
-	ch_buf_t *buf;
+		REALLOC(b->a, b->m);\
-	st_mt_t *mt;
+	}\
-} st_data_t;
+	b->a[b->n++] = y;\
-
+}\
-static void worker_for_insert(void *data, long i, int tid) // callback for kt_for()
+/**buf is the read block, k is the k-mer length, p = 12, len is the read length, seq is the read**/\
-{
+static void sf##_count_seq_buf(sf##_ch_buf_t *buf, int k, int p, int len, const char *seq) /**insert k-mers in $seq to linear buffer $buf**/\
-	st_data_t *s = (st_data_t*)data;
+{\
-	ch_buf_t *b = &s->buf[i];
+	int i, l;\
-	if (s->p->pt)
+	uint64_t x[4], mask = (1ULL<<k) - 1, shift = k - 1;\
-	{
+	for (i = l = 0, x[0] = x[1] = x[2] = x[3] = 0; i < len; ++i) {\
-		if(s->p->flag&HAF_COUNT_REFINE) b->n_ins += ha_pt_cnt_insert_list(s->p->pt, b->n, b->a);
+		int c = seq_nt4_table[(uint8_t)seq[i]];\
-		else b->n_ins += ha_pt_insert_list(s->p->pt, b->n, b->b_mz);
+		/**c = 00, 01, 10, 11**/\
-	}
+		if (c < 4) { /** not an "N" base**/\
-	else///for 0-th count, go into here
+			/**x[0] & x[1] are the forward k-mer**/\
-	{
+			/**x[2] & x[3] are the reverse complementary k-mer**/\
-		b->n_ins += ha_ct_insert_list(s->p->ct, s->p->create_new, b->n, b->a);
+			x[0] = (x[0] << 1 | (c&1))  & mask;\
-	}		
+			x[1] = (x[1] << 1 | (c>>1)) & mask;\
+			x[2] = x[2] >> 1 | (uint64_t)(1 - (c&1))  << shift;\
+			x[3] = x[3] >> 1 | (uint64_t)(1 - (c>>1)) << shift;\
+			if (++l >= k)\
+				sf##_ct_insert_buf(buf, p, yak_hash_long(x));\
+		} else l = 0, x[0] = x[1] = x[2] = x[3] = 0; /** if there is an "N", restart**/\
+	}\
+}\
+static void sf##_count_seq_buf_HPC(sf##_ch_buf_t *buf, int k, int p, int len, const char *seq) /**insert k-mers in $seq to linear buffer $buf**/\
+{\
+	int i, l, last = -1;\
+	uint64_t x[4], mask = (1ULL<<k) - 1, shift = k - 1;\
+	for (i = l = 0, x[0] = x[1] = x[2] = x[3] = 0; i < len; ++i) {\
+		int c = seq_nt4_table[(uint8_t)seq[i]];\
+		if (c < 4) { /** not an "N" base**/\
+			if (c != last) {\
+				x[0] = (x[0] << 1 | (c&1))  & mask;\
+				x[1] = (x[1] << 1 | (c>>1)) & mask;\
+				x[2] = x[2] >> 1 | (uint64_t)(1 - (c&1))  << shift;\
+				x[3] = x[3] >> 1 | (uint64_t)(1 - (c>>1)) << shift;\
+				if (++l >= k)\
+					sf##_ct_insert_buf(buf, p, yak_hash_long(x));\
+				last = c;\
+			}\
+		} else l = 0, last = -1, x[0] = x[1] = x[2] = x[3] = 0; /**if there is an "N", restart**/\
+	}\
+}\
+int sf##_ha_pt_insert_list(ha_pt_t *h, int n, const HType *a)\
+{\
+	int j, mask = (1<<h->pre) - 1, n_ins = 0;\
+	ha_pt1_t *g;\
+	if (n == 0) return 0;\
+	g = &h->h[a[0].x&mask];\
+	for (j = 0; j < n; ++j) {\
+		uint64_t x = a[j].x >> h->pre;\
+		khint_t k;\
+		int n;\
+		IType *p;\
+		assert((a[j].x&mask) == (a[0].x&mask));\
+		k = yak_pt_get(g->h, x<<YAK_COUNTER_BITS);\
+		if (k == kh_end(g->h)) continue; \
+		n = kh_key(g->h, k) & YAK_MAX_COUNT;\
+		assert(n < YAK_MAX_COUNT);\
+		p = &g->Ia[kh_val(g->h, k) + n];\
+		p->rid = a[j].rid, p->rev = a[j].rev, p->pos = a[j].pos, p->span = a[j].span;\
+		/**(uint64_t)a[j].rid<<36 | (uint64_t)a[j].rev<<35 | (uint64_t)a[j].pos<<8 | (uint64_t)a[j].span;**/\
+		++kh_key(g->h, k);\
+		++n_ins;\
+	}\
+	return n_ins;\
+}\
+/** data structure for each step in kt_pipeline()**/\
+typedef struct {pl_data_t *p;uint64_t n_seq0; int n_seq, m_seq, sum_len, nk, uq, *len; char **seq; VType *mz_buf; VType *mz;sf##_ch_buf_t *buf;st_mt_t *mt;} sf##_st_data_t;\
+static void sf##_worker_for_insert(void *data, long i, int tid) /** callback for kt_for()**/\
+{\
+	sf##_st_data_t *s = (sf##_st_data_t*)data;\
+	sf##_ch_buf_t *b = &s->buf[i];\
+	if (s->p->pt){\
+		if(s->p->flag&HAF_COUNT_REFINE) b->n_ins += ha_pt_cnt_insert_list(s->p->pt, b->n, b->a);\
+		else b->n_ins += sf##_ha_pt_insert_list(s->p->pt, b->n, b->b);\
+	}else{\
+		b->n_ins += ha_ct_insert_list(s->p->ct, s->p->create_new, b->n, b->a);\
+	}\
+}\
+static void sf##_worker_for_mz(void *data, long i, int tid)\
+{\
+	sf##_st_data_t *s = (sf##_st_data_t*)data;\
+	/**get the corresponding minimzer vector of this read**/\
+	VType *b = &s->mz_buf[tid];\
+	s->mz_buf[tid].n = 0;\
+	sf##_ha_sketch(s->seq[i], s->len[i], s->p->opt->w, s->p->opt->k, s->n_seq0 + i, s->p->opt->is_HPC, b, s->p->flt_tab, asm_opt.mz_sample_dist, 0, 0, \
+	(s->p->pt&&(s->p->flag&HAF_COUNT_REFINE))?s->p->pt:NULL, s->p->opt->min_rcnt, asm_opt.dp_min_len, asm_opt.dp_e, &(s->mt[tid]), asm_opt.mz_rewin, s->uq);\
+	s->mz[i].n = s->mz[i].m = b->n;\
+	MALLOC(s->mz[i].a, b->n);\
+	memcpy(s->mz[i].a, b->a, b->n * sizeof(VType));\
+}\
+static inline void sf##_pt_insert_buf(sf##_ch_buf_t *buf, int p, const HType *y){\
+	/**assign minimizer to one of 4096 bins by low 12 bits**/\
+	int pre = y->x & ((1<<p) - 1);\
+	sf##_ch_buf_t *b = &buf[pre];\
+	if (b->n == b->m) {\
+		b->m = b->m < 8? 8 : b->m + (b->m>>1);\
+		REALLOC(b->b, b->m);\
+	}\
+	b->b[b->n++] = *y;\
+}\
+static void *sf##_worker_count(void *data, int step, void *in) /** callback for kt_pipeline()**/\
+{\
+	pl_data_t *p = (pl_data_t*)data;\
+	if (step == 0) { /** step 1: read a block of sequences**/\
+		int ret;\
+		sf##_st_data_t *s;\
+		CALLOC(s, 1);\
+		s->p = p;\
+		s->n_seq0 = p->n_seq;\
+		if (p->rs_in && (p->flag & HAF_RS_READ)) {\
+			while (p->n_seq < p->rs_in->total_reads) {\
+				if ((p->flag & HAF_SKIP_READ) && p->rs_in->trio_flag[p->n_seq] != AMBIGU) {\
+					++p->n_seq;\
+					continue;\
+				}\
+				int l;\
+				recover_UC_Read(&p->ucr, p->rs_in, p->n_seq);\
+				l = p->ucr.length;\
+				if (s->n_seq == s->m_seq) {\
+					s->m_seq = s->m_seq < 16? 16 : s->m_seq + (s->m_seq>>1);\
+					REALLOC(s->len, s->m_seq);\
+					REALLOC(s->seq, s->m_seq);\
+				}\
+				MALLOC(s->seq[s->n_seq], l);\
+				memcpy(s->seq[s->n_seq], p->ucr.seq, l);\
+				s->len[s->n_seq++] = l;\
+				++p->n_seq;\
+				s->sum_len += l;\
+				s->nk += l >= p->opt->k? l - p->opt->k + 1 : 0;\
+				if (s->sum_len >= p->opt->chunk_size)\
+					break;\
+			}\
+		} else if(p->us_in) {\
+			ma_utg_t *u; s->uq = p->us_in->h;\
+			while (p->n_seq < p->us_in->n) {\
+				u = &(p->us_in->a[p->n_seq]);\
+				if (s->n_seq == s->m_seq) {\
+					s->m_seq = s->m_seq < 16? 16 : s->m_seq + (s->m_seq>>1);\
+					REALLOC(s->len, s->m_seq);\
+					REALLOC(s->seq, s->m_seq);\
+				}\
+				MALLOC(s->seq[s->n_seq], u->len);\
+				memcpy(s->seq[s->n_seq], u->s, u->len);\
+				s->len[s->n_seq++] = u->len;\
+				++p->n_seq;\
+				s->sum_len += u->len;\
+				s->nk += u->len >= p->opt->k? u->len - p->opt->k + 1 : 0;\
+				if (s->sum_len >= p->opt->chunk_size)\
+					break;\
+			}\
+		} else {\
+			while ((ret = kseq_read(p->ks)) >= 0) {\
+				int l = (int)(p->ks->seq.l) - (int)(p->opt->adaLen) - (int)(p->opt->adaLen);\
+				if(l <= 0) continue;\
+				if (p->n_seq >= 1<<28) {\
+					fprintf(stderr, "ERROR: this implementation supports no more than %d reads\n", 1<<28);\
+					exit(1);\
+				}\
+				if (p->rs_out) {\
+					/**for 0-th count, just insert read length to R_INF, instead of read**/\
+					if (p->flag & HAF_RS_WRITE_LEN) {\
+						assert(p->n_seq == p->rs_out->total_reads);\
+						ha_insert_read_len(p->rs_out, l, p->ks->name.l);\
+					} else if (p->flag & HAF_RS_WRITE_SEQ) {\
+						int i, n_N;\
+						assert(l == (int)p->rs_out->read_length[p->n_seq]);\
+						for (i = n_N = 0; i < l; ++i) /** count number of ambiguous bases**/\
+							if (seq_nt4_table[(uint8_t)p->ks->seq.s[i+p->opt->adaLen]] >= 4)\
+								++n_N;\
+						ha_compress_base(Get_READ(*p->rs_out, p->n_seq), p->ks->seq.s+p->opt->adaLen, l, &p->rs_out->N_site[p->n_seq], n_N);\
+						memcpy(&p->rs_out->name[p->rs_out->name_index[p->n_seq]], p->ks->name.s, p->ks->name.l);\
+					}\
+				}\
+				if (s->n_seq == s->m_seq) {\
+					s->m_seq = s->m_seq < 16? 16 : s->m_seq + (s->m_seq>>1);\
+					REALLOC(s->len, s->m_seq);\
+					REALLOC(s->seq, s->m_seq);\
+				}\
+				MALLOC(s->seq[s->n_seq], l);\
+				memcpy(s->seq[s->n_seq], p->ks->seq.s+p->opt->adaLen, l);\
+				s->len[s->n_seq++] = l;\
+				++p->n_seq;\
+				s->sum_len += l;\
+				s->nk += l >= p->opt->k? l - p->opt->k + 1 : 0;\
+				/**p->opt->chunk_size is the block max size**/\
+				if (s->sum_len >= p->opt->chunk_size)\
+					break;\
+			}\
+		}\
+		if (s->sum_len == 0) free(s);\
+		else return s;\
+	} else if (step == 1) { /** step 2: extract k-mers**/\
+		/**s is the block of reads**/\
+		sf##_st_data_t *s = (sf##_st_data_t*)in;\
+		int i, n_pre = 1<<p->opt->pre, m;\
+		/**allocate the k-mer buffer**/\
+		CALLOC(s->buf, n_pre);\
+		m = (int)(s->nk * 1.2 / n_pre) + 1;\
+		/**pre-allocate memory for each of 4096 buffer**/\
+		for (i = 0; i < n_pre; ++i) {\
+			s->buf[i].m = m;\
+			/**for 0-th counting, p->pt = NULL**/\
+			if (p->pt && !(p->flag&HAF_COUNT_REFINE)) MALLOC(s->buf[i].b, m);\
+			else MALLOC(s->buf[i].a, m);\
+		}\
+		if (p->opt->w == 1) { /** enumerate all k-mers**/\
+			int i;\
+			for (i = 0; i < s->n_seq; ++i) {\
+				if (p->opt->is_HPC)\
+					sf##_count_seq_buf_HPC(s->buf, p->opt->k, p->opt->pre, s->len[i], s->seq[i]);\
+				else\
+					sf##_count_seq_buf(s->buf, p->opt->k, p->opt->pre, s->len[i], s->seq[i]);\
+				if (!p->is_store) free(s->seq[i]);\
+			}\
+		} else { /** minimizers only**/\
+			uint32_t j;\
+			/**s->n_seq is how many reads at this buffer**/\
+			/**s->mz && s->mz_buf are lists of minimzer vectors**/\
+			CALLOC(s->mz, s->n_seq), CALLOC(s->mz_buf, p->opt->n_thread), CALLOC(s->mt, p->opt->n_thread);\
+			/**calculate minimzers for each read, each read corresponds to one thread**/\
+			kt_for(p->opt->n_thread, sf##_worker_for_mz, s, s->n_seq);\
+			for (i = 0; i < p->opt->n_thread; ++i) free(s->mt[i].a), free(s->mz_buf[i].a);\
+			free(s->mt), free(s->mz_buf);\
+			/**insert minimizers**/\
+			if (p->pt && !(p->flag&HAF_COUNT_REFINE)) {/**insert whole minimizer**/\
+				for (i = 0; i < s->n_seq; ++i)\
+					for (j = 0; j < s->mz[i].n; ++j)\
+						sf##_pt_insert_buf(s->buf, p->opt->pre, &s->mz[i].a[j]);\
+			} else {/**just insert the hash key of minimizer**/\
+				for (i = 0; i < s->n_seq; ++i)\
+					for (j = 0; j < s->mz[i].n; ++j)\
+						sf##_ct_insert_buf(s->buf, p->opt->pre, s->mz[i].a[j].x);\
+			}\
+			for (i = 0; i < s->n_seq; ++i) {\
+				p->n_mz += s->mz[i].n;\
+				free(s->mz[i].a);\
+				if (!p->is_store) free(s->seq[i]);\
+			}\
+			free(s->mz);\
+		}\
+		/**just clean seq**/\
+		free(s->seq); free(s->len);\
+		s->seq = 0, s->len = 0;\
+		return s;\
+	} else if (step == 2) { /** step 3: insert k-mers to hash table**/\
+		sf##_st_data_t *s = (sf##_st_data_t*)in;\
+		int i, n = 1<<p->opt->pre;uint64_t n_ins = 0;\
+		/**for 0-th counting, p->pt = NULL**/\
+		kt_for(p->opt->n_thread, sf##_worker_for_insert, s, n);\
+		/**n_ins is number of distinct k-mers**/\
+		for (i = 0; i < n; ++i) {\
+			n_ins += s->buf[i].n_ins;\
+			if (p->pt && !(p->flag&HAF_COUNT_REFINE)) free(s->buf[i].b);\
+			else free(s->buf[i].a);\
+		}\
+		if (p->ct) p->ct->tot += n_ins, p->ct->bs += s->sum_len;\
+		if (p->pt) p->pt->tot_pos += n_ins;\
+		free(s->buf);\
+		free(s);\
+	}\
+	return 0;\
 }
 
-static void worker_for_mz(void *data, long i, int tid)
+MZ_TEST_INIT(mz1, ha_mz1_t, ha_mz1_v, ha_idxpos_t, a)
-{
+MZ_TEST_INIT(mz2, ha_mzl_t, ha_mzl_v, ha_idxposl_t, al)
-	st_data_t *s = (st_data_t*)data;
-	///get the corresponding minimzer vector of this read
-	ha_mz1_v *b = &s->mz_buf[tid];
-	s->mz_buf[tid].n = 0;
-	ha_sketch(s->seq[i], s->len[i], s->p->opt->w, s->p->opt->k, s->n_seq0 + i, s->p->opt->is_HPC, b, s->p->flt_tab, asm_opt.mz_sample_dist, 0, 0, 
-	(s->p->pt&&(s->p->flag&HAF_COUNT_REFINE))?s->p->pt:NULL, s->p->opt->min_rcnt, asm_opt.dp_min_len, asm_opt.dp_e, &(s->mt[tid]), asm_opt.mz_rewin, s->uq);
-	s->mz[i].n = s->mz[i].m = b->n;
-	MALLOC(s->mz[i].a, b->n);
-	memcpy(s->mz[i].a, b->a, b->n * sizeof(ha_mz1_t));
-}
 
-MZ_FUNC_INIT(mz, ha_mz1_t)
-MZ_FUNC_INIT(mzl, ha_mzl_t)
-static void *worker_count(void *data, int step, void *in) // callback for kt_pipeline()
-{
-	pl_data_t *p = (pl_data_t*)data;
-	if (step == 0) { // step 1: read a block of sequences
-		int ret;
-		st_data_t *s;
-		CALLOC(s, 1);
-		s->p = p;
-		s->n_seq0 = p->n_seq;
-		if (p->rs_in && (p->flag & HAF_RS_READ)) {
-			while (p->n_seq < p->rs_in->total_reads) {
-				if ((p->flag & HAF_SKIP_READ) && p->rs_in->trio_flag[p->n_seq] != AMBIGU) {
-					++p->n_seq;
-					continue;
-				}
-				int l;
-				recover_UC_Read(&p->ucr, p->rs_in, p->n_seq);
-				l = p->ucr.length;
-				if (s->n_seq == s->m_seq) {
-					s->m_seq = s->m_seq < 16? 16 : s->m_seq + (s->m_seq>>1);
-					REALLOC(s->len, s->m_seq);
-					REALLOC(s->seq, s->m_seq);
-				}
-				MALLOC(s->seq[s->n_seq], l);
-				memcpy(s->seq[s->n_seq], p->ucr.seq, l);
-				s->len[s->n_seq++] = l;
-				++p->n_seq;
-				s->sum_len += l;
-				s->nk += l >= p->opt->k? l - p->opt->k + 1 : 0;
-				if (s->sum_len >= p->opt->chunk_size)
-					break;
-			}
-		} else if(p->us_in) {
-			ma_utg_t *u; s->uq = 1;
-			while (p->n_seq < p->us_in->n) {
-				u = &(p->us_in->a[p->n_seq]);
-				if (s->n_seq == s->m_seq) {
-					s->m_seq = s->m_seq < 16? 16 : s->m_seq + (s->m_seq>>1);
-					REALLOC(s->len, s->m_seq);
-					REALLOC(s->seq, s->m_seq);
-				}
-				MALLOC(s->seq[s->n_seq], u->len);
-				memcpy(s->seq[s->n_seq], u->s, u->len);
-				s->len[s->n_seq++] = u->len;
-				++p->n_seq;
-				s->sum_len += u->len;
-				s->nk += u->len >= p->opt->k? u->len - p->opt->k + 1 : 0;
-				if (s->sum_len >= p->opt->chunk_size)
-					break;
-			}
-		} else {
-			while ((ret = kseq_read(p->ks)) >= 0) {
-				int l = (int)(p->ks->seq.l) - (int)(p->opt->adaLen) - (int)(p->opt->adaLen);
-				if(l <= 0) continue;
-
-				if (p->n_seq >= 1<<28) {
-					fprintf(stderr, "ERROR: this implementation supports no more than %d reads\n", 1<<28);
-					exit(1);
-				}
-				if (p->rs_out) {
-					///for 0-th count, just insert read length to R_INF, instead of read
-					if (p->flag & HAF_RS_WRITE_LEN) {
-						assert(p->n_seq == p->rs_out->total_reads);
-						ha_insert_read_len(p->rs_out, l, p->ks->name.l);
-					} else if (p->flag & HAF_RS_WRITE_SEQ) {
-						int i, n_N;
-						assert(l == (int)p->rs_out->read_length[p->n_seq]);
-						for (i = n_N = 0; i < l; ++i) // count number of ambiguous bases
-							if (seq_nt4_table[(uint8_t)p->ks->seq.s[i+p->opt->adaLen]] >= 4)
-								++n_N;
-						ha_compress_base(Get_READ(*p->rs_out, p->n_seq), p->ks->seq.s+p->opt->adaLen, l, &p->rs_out->N_site[p->n_seq], n_N);
-						memcpy(&p->rs_out->name[p->rs_out->name_index[p->n_seq]], p->ks->name.s, p->ks->name.l);
-					}
-				}
-				///for 0-th count, insert both seq and length to local block
-				if (s->n_seq == s->m_seq) {
-					s->m_seq = s->m_seq < 16? 16 : s->m_seq + (s->m_seq>>1);
-					REALLOC(s->len, s->m_seq);
-					REALLOC(s->seq, s->m_seq);
-				}
-				MALLOC(s->seq[s->n_seq], l);
-				memcpy(s->seq[s->n_seq], p->ks->seq.s+p->opt->adaLen, l);
-				s->len[s->n_seq++] = l;
-				++p->n_seq;
-				s->sum_len += l;
-				s->nk += l >= p->opt->k? l - p->opt->k + 1 : 0;
-				///p->opt->chunk_size is the block max size
-				if (s->sum_len >= p->opt->chunk_size)
-					break;
-			}
-		}
-		if (s->sum_len == 0) free(s);
-		else return s;
-	} else if (step == 1) { // step 2: extract k-mers
-		///s is the block of reads
-		st_data_t *s = (st_data_t*)in;
-		if(p->us_in) mzl_init_kuf(p, s);
-		else mz_init_kuf(p, s);
-		// fill the buffer
-		///for 0-th counting, p->opt->w == 1
-		if (p->opt->w == 1) { // enumerate all k-mers
-			///scan all reads
-			int i;
-			for (i = 0; i < s->n_seq; ++i) {
-				if (p->opt->is_HPC)
-					count_seq_buf_HPC(s->buf, p->opt->k, p->opt->pre, s->len[i], s->seq[i]);
-				else
-					count_seq_buf(s->buf, p->opt->k, p->opt->pre, s->len[i], s->seq[i]);
-				if (!p->is_store) free(s->seq[i]);
-			}
-		} else { // minimizers only
-			if(p->us_in) mzl_mselect(p, s);
-			else mz_mselect(p, s);
-		}
-		///just clean seq
-		free(s->seq); free(s->len);
-		s->seq = 0, s->len = 0;
-		return s;
-	} else if (step == 2) { // step 3: insert k-mers to hash table
-		st_data_t *s = (st_data_t*)in;
-		///for 0-th counting, p->pt = NULL
-		kt_for(p->opt->n_thread, worker_for_insert, s, 1<<p->opt->pre);
-		if(p->us_in) mzl_destory_kuf(p, s, 1<<p->opt->pre);
-		else mz_destory_kuf(p, s, 1<<p->opt->pre);
-	}
-	return 0;
-}
 
 void debug_adapter(const hifiasm_opt_t *asm_opt, All_reads *rs)
 {
@@ -951,7 +919,8 @@ static ha_ct_t *yak_count(const yak_copt_t *opt, const char *fn, int flag, ha_pt
 		pl.ct = ha_ct_init(opt->k, opt->pre, opt->bf_n_hash, opt->bf_shift);
 	}
 	if(pl.ct) pl.ct->bs = 0;
-	kt_pipeline(3, worker_count, &pl, 3);
+	if(ug_rs) kt_pipeline(3, mz2_worker_count, &pl, 3);
+	else kt_pipeline(3, mz1_worker_count, &pl, 3);
 	if (read_rs) {
 		destory_UC_Read(&pl.ucr);
 	} else if(!read_rs && !ug_rs) {
@@ -963,7 +932,7 @@ static ha_ct_t *yak_count(const yak_copt_t *opt, const char *fn, int flag, ha_pt
 	return pl.ct;
 }
 
-ha_ct_t *ha_count(const hifiasm_opt_t *asm_opt, int flag, ha_pt_t *p0, const void *flt_tab, All_reads *rs, ma_utg_v *us, int keep_adapter, int *low_freq)
+ha_ct_t *ha_count(const hifiasm_opt_t *asm_o, int flag, int HPC, int k, int w, ha_pt_t *p0, const void *flt_tab, All_reads *rs, ma_utg_v *us, int keep_adapter, int *low_freq)
 {
 	int i;
 	int64_t n_seq = 0;
@@ -979,20 +948,19 @@ ha_ct_t *ha_count(const hifiasm_opt_t *asm_opt, int flag, ha_pt_t *p0, const voi
 			malloc_All_reads(rs);
 	}
 	yak_copt_init(&opt);
-	opt.k = us? asm_opt->ul_mer_length:asm_opt->k_mer_length;
+	opt.k = k;
-	///always 0
+	opt.is_HPC = HPC;
-	opt.is_HPC = !(asm_opt->flag&HA_F_NO_HPC);
 	///for ft-counting, shoud be 1
-	opt.w = flag & HAF_COUNT_ALL? 1 : (us? asm_opt->ul_mz_win:asm_opt->mz_win);
+	opt.w = flag & HAF_COUNT_ALL? 1 : w;
 	///for ft-counting, shoud be 37
 	///for ha_pt_gen, shoud be 0
-	opt.bf_shift = flag & HAF_COUNT_EXACT? 0 : asm_opt->bf_shift;
+	opt.bf_shift = flag & HAF_COUNT_EXACT? 0 : asm_o->bf_shift;
-	opt.n_thread = asm_opt->thread_num;
+	opt.n_thread = asm_o->thread_num;
-	opt.adaLen = (keep_adapter? asm_opt->adapterLen : 0);
+	opt.adaLen = (keep_adapter? asm_o->adapterLen : 0);
 	opt.min_rcnt = (low_freq?*low_freq:-1);
 	///asm_opt->num_reads is the number of fastq files
-	for (i = n_bs = 0; i < (us?1:asm_opt->num_reads); ++i){
+	for (i = n_bs = 0; i < (us?1:asm_o->num_reads); ++i){
-		h = yak_count(&opt, asm_opt->read_file_names[i], flag|HAF_CREATE_NEW, p0, h, flt_tab, rs, us, &n_seq);
+		h = yak_count(&opt, asm_o->read_file_names[i], flag|HAF_CREATE_NEW, p0, h, flt_tab, rs, us, &n_seq);
 		if(h) n_bs += h->bs;
 	}
 	if(h) h->bs = n_bs;	
@@ -1077,32 +1045,32 @@ void debug_ct_index(void* q_ct_idx, void* r_ct_idx)
  * High-level interfaces *
  *************************/
 
-void *ha_ft_ug_gen(const hifiasm_opt_t *asm_opt, ma_utg_v *us, int hap_n)
+void *ha_ft_ug_gen(const hifiasm_opt_t *asm_opt, ma_utg_v *us, int is_HPC, int k, int w, int min_freq, int max_freq)
 {
 	yak_ft_t *flt_tab;
 	ha_ct_t *h;
 	///HAF_COUNT_EXACT ---> no bf; HAF_COUNT_ALL ---> no minimizer
-	h = ha_count(asm_opt, HAF_COUNT_ALL|HAF_UG_READ|HAF_COUNT_EXACT, NULL, NULL, NULL, us, 0, NULL);
+	h = ha_count(asm_opt, HAF_COUNT_ALL|HAF_UG_READ|HAF_COUNT_EXACT, is_HPC, k, w, NULL, NULL, NULL, us, 0, NULL);
-	ha_ct_shrink(h, 1, YAK_MAX_COUNT-1, asm_opt->thread_num);
+	ha_ct_shrink(h, min_freq, max_freq>YAK_MAX_COUNT-1?YAK_MAX_COUNT-1:max_freq, asm_opt->thread_num);
 	flt_tab = gen_hh(h, asm_opt->max_kmer_cnt);
 	ha_ct_destroy(h);
 	return (void*)flt_tab;
 }
 
 
-ha_pt_t *ha_pt_ug_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, ma_utg_v *us, int hap_n)
+ha_pt_t *ha_pt_ug_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, ma_utg_v *us, int is_HPC, int k, int w, int min_freq)
 {
 	ha_ct_t *ct;
 	ha_pt_t *pt;
 	///HAF_COUNT_EXACT: no bf
-	ct = ha_count(asm_opt, HAF_COUNT_EXACT|HAF_UG_READ, NULL, flt_tab, NULL, us, 0, NULL);
+	ct = ha_count(asm_opt, HAF_COUNT_EXACT|HAF_UG_READ, is_HPC, k, w, NULL, flt_tab, NULL, us, 0, NULL);
 	fprintf(stderr, "[M::%s::%.3f*%.2f] ==> counted %ld distinct minimizer k-mers\n", __func__,
 			yak_realtime(), yak_cpu_usage(), (long)ct->tot);
 	///minimizer with YAK_MAX_COUNT occ may apper > YAK_MAX_COUNT times, so it may lead to overflow at ha_pt_gen
-	ha_ct_shrink(ct, 1, YAK_MAX_COUNT - 1, asm_opt->thread_num);
+	ha_ct_shrink(ct, min_freq, YAK_MAX_COUNT - 1, asm_opt->thread_num);
 
 	pt = ha_pt_gen(ct, asm_opt->thread_num, 1);
-	ha_count(asm_opt, HAF_COUNT_EXACT|HAF_UG_READ, pt, flt_tab, NULL, us, 0, NULL);
+	ha_count(asm_opt, HAF_COUNT_EXACT|HAF_UG_READ, is_HPC, k, w, pt, flt_tab, NULL, us, 0, NULL);
 	//ha_pt_sort(pt, asm_opt->thread_num);
 	fprintf(stderr, "[M::%s::%.3f*%.2f] ==> indexed %ld positions\n", __func__,
 			yak_realtime(), yak_cpu_usage(), (long)pt->tot_pos);
@@ -1116,7 +1084,7 @@ void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov, int i
 	int peak_hom, peak_het, cutoff = YAK_MAX_COUNT - 1, ex_flag = 0;
 	if(is_hp_mode) ex_flag = HAF_RS_READ|HAF_SKIP_READ;
 	ha_ct_t *h;
-	h = ha_count(asm_opt, HAF_COUNT_ALL|HAF_RS_WRITE_LEN|ex_flag, NULL, NULL, rs, NULL, 1, NULL);
+	h = ha_count(asm_opt, HAF_COUNT_ALL|HAF_RS_WRITE_LEN|ex_flag, !(asm_opt->flag&HA_F_NO_HPC), asm_opt->k_mer_length, asm_opt->mz_win, NULL, NULL, rs, NULL, 1, NULL);
 	if((asm_opt->flag & HA_F_VERBOSE_GFA))
 	{
 		write_ct_index((void*)h, asm_opt->output_file_name);
@@ -1148,12 +1116,12 @@ ha_pt_t *ha_pt_gen_dp(const hifiasm_opt_t *asm_opt, ha_ct_t *ct, int flag, int n
 {
 	int low_freq = mz_low_b(peak_hom, peak_het);
 	ha_pt_t *pt = ha_pt_gen_count(ct, n_thread); ///key = cnt, val = 0
-	ha_count(asm_opt, HAF_COUNT_EXACT|HAF_COUNT_REFINE|flag, pt, flt_tab, rs, NULL, 1, &low_freq);
+	ha_count(asm_opt, HAF_COUNT_EXACT|HAF_COUNT_REFINE|flag, !(asm_opt->flag&HA_F_NO_HPC), asm_opt->k_mer_length, asm_opt->mz_win, pt, flt_tab, rs, NULL, 1, &low_freq);
 	uint64_t occ = ha_pt_shrink(pt, n_thread);
 	if(flag&HAF_RS_WRITE_LEN) flag -= HAF_RS_WRITE_LEN;
 	if(flag&HAF_RS_WRITE_SEQ) flag -= HAF_RS_WRITE_SEQ;
 	flag |= HAF_RS_READ; pt->tot_pos = 0;
-	ha_count(asm_opt, HAF_COUNT_EXACT|flag, pt, flt_tab, rs, NULL, 1, NULL);
+	ha_count(asm_opt, HAF_COUNT_EXACT|flag, !(asm_opt->flag&HA_F_NO_HPC), asm_opt->k_mer_length, asm_opt->mz_win, pt, flt_tab, rs, NULL, 1, NULL);
 	// fprintf(stderr, "[M::%s::] counted %lu distinct minimizer k-mers\n", __func__, pt->tot);
 	// fprintf(stderr, "[M::%s::] collected %lu minimizers\n\n\n", __func__, pt->tot_pos);
 	assert(occ == pt->tot_pos);
@@ -1177,7 +1145,7 @@ ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_f
 	}
 	if(is_hp_mode) extra_flag1 |= HAF_SKIP_READ, extra_flag2 |= HAF_SKIP_READ;
 
-	ct = ha_count(asm_opt, HAF_COUNT_EXACT|extra_flag1, NULL, flt_tab, rs, NULL, 1, NULL);
+	ct = ha_count(asm_opt, HAF_COUNT_EXACT|extra_flag1,  !(asm_opt->flag&HA_F_NO_HPC), asm_opt->k_mer_length, asm_opt->mz_win, NULL, flt_tab, rs, NULL, 1, NULL);
 	fprintf(stderr, "[M::%s::%.3f*%.2f] ==> counted %ld distinct minimizer k-mers\n", __func__,
 			yak_realtime(), yak_cpu_usage(), (long)ct->tot);
 	ha_ct_hist(ct, cnt, asm_opt->thread_num);
@@ -1203,7 +1171,7 @@ ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_f
 	{
 		fprintf(stderr, "[M::%s::] counting in normal mode\n", __func__);
 		pt = ha_pt_gen(ct, asm_opt->thread_num, 0);
-		ha_count(asm_opt, HAF_COUNT_EXACT|extra_flag2, pt, flt_tab, rs, NULL, 1, NULL);
+		ha_count(asm_opt, HAF_COUNT_EXACT|extra_flag2,  !(asm_opt->flag&HA_F_NO_HPC), asm_opt->k_mer_length, asm_opt->mz_win, pt, flt_tab, rs, NULL, 1, NULL);
 		assert((uint64_t)tot_cnt == pt->tot_pos);
 	}
 	else

htab.h

@@ -25,14 +25,12 @@ typedef struct { uint32_t n, m; ha_mz1_t *a; } ha_mz1_v;
 
 typedef struct {
 	uint64_t x; ///x is the hash key
-	uint64_t rid:31, rev:1;
+	uint64_t rid:31, rev:1, pos:32;
-	uint32_t pos;
 	uint8_t span;
 } ha_mzl_t;
 
 typedef struct {
-	uint64_t rid:31, rev:1;
+	uint64_t rid:31, rev:1, pos:32;
-	uint32_t pos;
 	uint8_t span;
 } ha_idxposl_t;
 
@@ -71,12 +69,12 @@ extern void *ha_flt_tab_hp;
 extern ha_pt_t *ha_idx_hp;
 extern void *ha_ct_table;
 
-void *ha_ft_ug_gen(const hifiasm_opt_t *asm_opt, ma_utg_v *us, int hap_n);
+void *ha_ft_ug_gen(const hifiasm_opt_t *asm_opt, ma_utg_v *us, int is_HPC, int k, int w, int min_freq, int max_freq);
 void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov, int is_hp_mode);
 int32_t ha_ft_cnt(const void *hh, uint64_t y);
 void ha_ft_destroy(void *h);
 
-ha_pt_t *ha_pt_ug_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, ma_utg_v *us, int hap_n);
+ha_pt_t *ha_pt_ug_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, ma_utg_v *us, int is_HPC, int k, int w, int min_freq);
 ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_from_store, int is_hp_mode, All_reads *rs, int *hom_cov, int *het_cov);
 void ha_pt_destroy(ha_pt_t *h);
 const ha_idxpos_t *ha_pt_get(const ha_pt_t *h, uint64_t hash, int *n);
@@ -101,7 +99,8 @@ double yak_cpu_usage(void);
 
 void ha_triobin(const hifiasm_opt_t *opt);
 
-void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mz1_v *p, const void *hf, int sample_dist, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct, ha_pt_t *pt, int min_freq, int32_t dp_min_len, float dp_e, st_mt_t *mt, int32_t ws, int32_t is_unique);
+void mz1_ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mz1_v *p, const void *hf, int sample_dist, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct, ha_pt_t *pt, int min_freq, int32_t dp_min_len, float dp_e, st_mt_t *mt, int32_t ws, int32_t is_unique);
+void mz2_ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mzl_v *p, const void *hf, int sample_dist, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct, ha_pt_t *pt, int min_freq, int32_t dp_min_len, float dp_e, st_mt_t *mt, int32_t ws, int32_t is_unique);
 int ha_analyze_count(int n_cnt, int start_cnt, int m_peak_hom, const int64_t *cnt, int *peak_het);
 int adj_m_peak_hom(int m_peak_hom, int max_i, int max2_i, int max3_i, int *peak_het);
 void print_hist_lines(int n_cnt, int start_cnt, const int64_t *cnt);

inter.cpp

@@ -7,19 +7,26 @@
 #include "CommandLines.h"
 #include "htab.h"
 
-void uidx_build(ma_ug_t *ug, int hap_n)
+void uidx_build(ma_ug_t *ug, int is_HPC, int k, int w, int hap_n)
 {
     int flag = asm_opt.flag;
     asm_opt.flag |= HA_F_NO_HPC;
-    ha_flt_tab = ha_ft_ug_gen(&asm_opt, &(ug->u), hap_n);
+    ug->u.h = hap_n;
-    ha_idx = ha_pt_ug_gen(&asm_opt, ha_flt_tab, &(ug->u), hap_n);
+    ha_flt_tab = ha_ft_ug_gen(&asm_opt, &(ug->u), is_HPC, k, w, hap_n, hap_n*10);
-
+    ha_idx = ha_pt_ug_gen(&asm_opt, ha_flt_tab, &(ug->u), is_HPC, k, w, hap_n);
-
-
-
-
-
     asm_opt.flag = flag;
+}
+
+void uidx_destory()
+{
     ha_ft_destroy(ha_flt_tab); 
     ha_pt_destroy(ha_idx);
 }
+
+
+void ul_resolve(ma_ug_t *ug, int hap_n)
+{
+    uidx_build(ug, 1, 63, 63, hap_n);
+    
+    uidx_destory();
+}

rcut.cpp

@@ -2767,6 +2767,48 @@ void clean_ovlp_by_mc(mc_g_t *mg, hap_overlaps_list* ha)
 }
 
 
+void filter_ta_by_mc(ma_ug_t *ug, asg_t *read_g, uint32_t uID, kv_u_trans_t* ta, mc_match_t *ma, int8_t *s)
+{
+	mc_edge_t *o = pt_a(*ma, uID);
+	uint32_t n = pt_n(*ma, uID), k, qn, tn;
+	u_trans_t *p = NULL;
+	for (k = 0; k < n; ++k)
+	{
+		qn = ma_x(o[k]); tn = ma_y(o[k]); p = NULL;
+		if((s[qn]*s[tn])!=-1) continue;
+		get_u_trans_spec(ta, qn, tn, &p, NULL);
+		if(p && p->nw == o[k].w) {
+			p->del = 0; 
+		}
+		else {
+			get_u_trans_spec(ta, tn, qn, &p, NULL);
+			if(p && p->nw == o[k].w) p->del = 0;
+		}
+		if(!p) fprintf(stderr, "ERROR-ta-p\n");
+	}
+}
+
+
+void clean_ta_by_mc(mc_g_t *mg, kv_u_trans_t *ta)
+{
+	uint32_t v, i;
+	u_trans_t *p = NULL;
+	for (i = 0; i < ta->n; i++) ta->a[i].del = 1;
+	for (i = 0; i < mg->e->n_seq; ++i) filter_ta_by_mc(mg->ug, mg->rg, i, ta, mg->e, mg->s.a);
+	for (i = v = 0; i < ta->n; i++) {
+		if(ta->a[i].del) continue;
+		ta->a[v++] = ta->a[i];
+	}
+	ta->n = v;
+	for (i = 0; i < v; i++) {
+		kv_pushp(u_trans_t, *ta, &p);
+		(*p) = ta->a[i];
+		p->qn = ta->a[i].tn; p->qs = ta->a[i].ts; p->qe = ta->a[i].te;
+        p->tn = ta->a[i].qn; p->ts = ta->a[i].qs; p->te = ta->a[i].qe;
+	}
+	kt_u_trans_t_idx(ta, mg->ug->g->n_seq);
+}
+
 void p_nodes(mc_g_t *mg, trans_chain* t_ch, uint8_t* trio_flag)
 {
 	uint32_t i;
@@ -2844,7 +2886,7 @@ void print_hap_s(int8_t *s, uint32_t sn)
 	}
 }
 
-void mc_solve(hap_overlaps_list* ovlp, trans_chain* t_ch, kv_u_trans_t *ta, ma_ug_t *ug, asg_t *read_g, double f_rate, uint8_t* trio_flag, uint32_t renew_s, int8_t *s, uint32_t is_sys, bubble_type* bub, kv_u_trans_t *ref)
+void mc_solve(hap_overlaps_list* ovlp, trans_chain* t_ch, kv_u_trans_t *ta, ma_ug_t *ug, asg_t *read_g, double f_rate, uint8_t* trio_flag, uint32_t renew_s, int8_t *s, uint32_t is_sys, bubble_type* bub, kv_u_trans_t *ref, int clean_ov)
 {
 	mc_opt_t opt;
 	mc_opt_init(&opt, asm_opt.n_perturb, asm_opt.f_perturb, asm_opt.seed);
@@ -2863,8 +2905,10 @@ void mc_solve(hap_overlaps_list* ovlp, trans_chain* t_ch, kv_u_trans_t *ta, ma_u
 		p_nodes(mg, t_ch, trio_flag);
 	}
 
-	if(ovlp) clean_ovlp_by_mc(mg, ovlp);
+	if(clean_ov){
-
+		if(ovlp) clean_ovlp_by_mc(mg, ovlp);
+		if(ta) clean_ta_by_mc(mg, ta);
+	}
 	// print_hap_s(s, ug->u.n);
 
 	destory_mc_g_t(&mg);	

rcut.h

@@ -118,7 +118,7 @@ static inline double kr_drand_r(uint64_t *x)
     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, kv_u_trans_t *ref);
+void mc_solve(hap_overlaps_list* ovlp, trans_chain* t_ch, kv_u_trans_t *ta, ma_ug_t *ug, asg_t *read_g, double f_rate, uint8_t* trio_flag, uint32_t renew_s, int8_t *s, uint32_t is_sys, bubble_type* bub, kv_u_trans_t *ref, int clean_ov);
 void debug_mc_g_t(const char* name);
 void mc_solve_general(kv_u_trans_t *ta, uint32_t un, kv_gg_status *s, uint16_t hapN, uint16_t update_ta, uint16_t write_dump);
 kv_gg_status *init_mc_gg_status(ma_ug_t *ug, asg_t *read_g, ma_sub_t* coverage_cut, 

sketch.cpp

@@ -10,14 +10,8 @@
 #define MAX_HIGH_OCC     8   // TODO: don't hard code if we need to tune this parameter
 #define MAX_MAX_HIGH_OCC 16
 #define GMC(a, x,y,xn) ((a)[(x)*(xn)+(y)])
-
+#define GL(x, i) ((int64_t)((uint32_t)((x).a[(i)])))
-static inline int mzcmp(const ha_mz1_t *a, const ha_mz1_t *b)
+#define A_M(p, i) ((i) >= 0 && (p).a[(i)].rid > 0)
-{
-	return a->rid < b->rid? -1 : a->rid > b->rid? 1 : ((a->x > b->x) - (a->x < b->x));
-}
-
-#define mz_lt(a, b) (mzcmp(&(a), &(b)) < 0)
-KSORT_INIT(mz, ha_mz1_t, mz_lt)
 
 void debug_refine(ha_mz1_t *ma, uint64_t *mmt, int32_t sn, int32_t n, int32_t m, int32_t end)
 {
@@ -42,287 +36,6 @@ void debug_refine(ha_mz1_t *ma, uint64_t *mmt, int32_t sn, int32_t n, int32_t m,
 	if(nt != tot) fprintf(stderr, "ERROR-TOT, nt: %ld, tot: %ld\n", nt, tot);
 }
 
-void refine_select(ha_mz1_v *mz, int32_t sidx, int32_t eidx, int32_t sn, int32_t min_freq, st_mt_t *mm, 
-int32_t *rsi, int32_t *rei)
-{
-	int32_t n = sn, m = eidx + 1 - sidx, i, k, t, mk=-1;
-	uint64_t ix, kx, ks;
-	kv_resize(uint64_t, *mm, mm->n+n*m);
-	ha_mz1_t *ma = mz->a + sidx;
-	uint64_t *mmt = mm->a + mm->n;
-	// fprintf(stderr, "[M::%s::] ==> +n: %d, m: %d, sn: %d, sidx: %d, eidx: %d\n", __func__, n, m, sn, sidx, eidx);
-
-	for (i = 0; i < n; i++) ///how many selected minimizers 
-	{
-		for (k = 0, mk = -1; k < m; k++) ///how many minimizers in total
-		{
-			if((int32_t)(ma[k].rid)<min_freq) continue;
-			
-			ks = ma[k].pos + 1 - ma[k].span; t = -1;
-			if(i > 0)
-			{
-				for (t = k-1; t >= 0 && (ma[t].pos >= ks||(int32_t)(ma[t].rid)<min_freq); t--);
-			} 
-			
-			ix = (i <= 0?0:(t<0?0xffffffff:(GMC(mmt, i-1,t,m)&0xffffffff)));
-			if(ix < 0xffffffff) ix += (ma[k].rid);
-			kx = (mk < 0?0xffffffff:(GMC(mmt, i, mk,m)&0xffffffff));
-
-			
-
-			ks = MIN(ix, kx);
-			// fprintf(stderr, "ks: %lu, i: %d (n-%d), k: %d (m-%d), ix: %lu, kx: %lu, t: %d, mk: %d\n", 
-			// ks, i, n, k, m, ix, kx, t, mk);
-			if((ks&0xffffffff) == 0xffffffff) ks |= ((uint64_t)0xffffffff)<<32;
-			else if(ks == ix) ks |= (uint64_t)(i>0?(i-1)*m+t:0xffffffff)<<32;
-			else if(ks == kx) ks |= (uint64_t)(mk>=0?i*m+mk:0xffffffff)<<32;
-			
-			GMC(mmt, i,k,m) = ks;	
-			mk = k;
-		}
-	}
-	// fprintf(stderr, "[M::%s::] ==> ++n: %d, m: %d, sn: %d, sidx: %d, eidx: %d\n", __func__, n, m, sn, sidx, eidx);
-
-	ks = (n-1)*m + mk; ix = (uint64_t)-1; kx = 0;
-	while (ks != 0xffffffff)
-	{
-		i = ks/m; k = ks%m;
-		ks = mmt[ks]>>32;
-		// fprintf(stderr, "i: %d, k: %d, ks: %lu\n", i, k, ks);
-		if(ks == 0xffffffff || (int32_t)(ks/m) == (i-1))
-		{
-			mm->a[sidx+k] = 1;
-			ix = MIN((uint64_t)k, ix); kx = MAX((uint64_t)k, kx);
-		}
-	}
-	///debug
-	// debug_refine(ma, mmt, sn, n, m, (n-1)*m + mk);
-
-	if(rsi) (*rsi) = ix + sidx;
-	if(rei) (*rei) = kx + sidx;	
-}
-
-void refine_sketch(ha_mz1_v *p, ha_pt_t *pt, int32_t rlen, int32_t dp_min_len, float er, int32_t min_freq, st_mt_t *mt)
-{
-	// fprintf(stderr, "[M::%s::] ==> #########10#########, rlen: %d\n", __func__, rlen);
-
-	int32_t i, n = p->n, bd, len = MIN(rlen, dp_min_len), sublen, cnt, ei, li, ri;
-	int32_t sn =  len*er + 1;
-	kv_resize(uint64_t, *mt, (int64_t)p->n); 
-	mt->n = p->n; memset(mt->a, 0, sizeof(uint64_t)*p->n);
-	for (i = 0; i < n; i++) p->a[i].rid = ha_pt_cnt(pt, p->a[i].x);
-
-	for (i = cnt = 0, bd = -1, ei = -1; i < n; i++)
-	{
-		if((int32_t)(p->a[i].rid)<min_freq) continue;
-		sublen = p->a[i].pos + 1;
-		if(sublen > len) break;
-		else ei = i;
-		
-		if((int32_t)(p->a[i].pos + 1 - p->a[i].span) > bd)
-		{
-			bd = p->a[i].pos;
-			cnt++;
-		}
-	}
-
-	// fprintf(stderr, "[M::%s::] ==> +cnt: %d, sn: %d, ei: %d, n: %d\n", __func__, cnt, sn, ei, n);
-
-	
-	if(cnt >= sn) refine_select(p, 0, ei, sn, min_freq, mt, NULL, &li);
-	else
-	{
-		li = i-1; 
-		for (i = 0; i <= li; i++) mt->a[i] = 1;
-	} 
-	
-	
-	if(len < rlen)
-	{
-		for (i = n-1, cnt = 0, bd = rlen+1, ei = -1; i >= 0; i--)
-		{
-			if((int32_t)(p->a[i].rid)<min_freq) continue;
-			sublen = rlen - (p->a[i].pos + 1 - p->a[i].span);
-			if(sublen > len) break;
-			else ei = i;
-
-			if((int32_t)(p->a[i].pos) < bd)
-			{
-				bd = p->a[i].pos + 1 - p->a[i].span;
-				cnt++;
-			}
-		}
-
-		// fprintf(stderr, "[M::%s::] ==> -cnt: %d, sn: %d, ei: %d, n: %d\n", __func__, cnt, sn, ei, n);
-
-		
-		if(cnt >= sn) refine_select(p, ei, n-1, sn, min_freq, mt, &ri, NULL);
-		else 
-		{
-			ri = i+1;
-			for (i = ri; i <= n-1; i++) mt->a[i] = 1;
-		} 
-		
-		// fprintf(stderr, "[M::%s::] ==> --cnt: %d, sn: %d, ei: %d, n: %d\n", __func__, cnt, sn, ei, n);
-
-		if(ri - li >= 2)
-		{
-			li++; ri--;
-			sn =  (p->a[ri].pos - p->a[li].pos + p->a[li].span)*er + 1;
-			for (i = li, cnt = 0, bd = -1; i <= ri; i++)
-			{
-				if((int32_t)(p->a[i].rid)<min_freq) continue;
-				if((int32_t)(p->a[i].pos + 1 - p->a[i].span) > bd)
-				{
-					bd = p->a[i].pos;
-					cnt++;
-					if(cnt >= sn) break;
-				}
-			}
-
-			if(cnt >= sn) refine_select(p, li, ri, sn, min_freq, mt, NULL, NULL);
-			else for (i = li; i <= ri; i++) mt->a[i] = 1;
-		}
-	}
-
-	// fprintf(stderr, "[M::%s::] ==> #########20#########, p->n: %u, n: %d\n", __func__, p->n, n);
-	for (i = sn = 0; i < n; i++)
-	{
-		if(mt->a[i])
-		{
-			p->a[sn] = p->a[i];
-			sn++;
-		}
-	}
-	// if(p->n != sn) fprintf(stderr, "[M::%s::] ==> #########21#########, p->n: %u, sn: %d\n", __func__, p->n, sn);
-	p->n = sn;
-	
-}
-
-inline int hf_dp(ha_mz1_v *mz, int32_t sidx, int32_t eidx, int32_t sn, int32_t min_freq, st_mt_t *mm, 
-int32_t *rsi, int32_t *rei)
-{
-	return 0;
-}
-
-inline void hf_select(ha_mz1_v *p, int32_t si, int32_t ei, int32_t n, int32_t len, int32_t sample_dist, ha_mz1_t *b, int32_t force)
-{
-	if(ei - si <= 1) return;
-	int32_t ps = si < 0? 0 : p->a[si].pos;
-	int32_t pe = ei == n? len : p->a[ei].pos;
-	int32_t j, k, st = si + 1, en = ei;
-	int32_t max_high_occ = (int32_t)((double)(pe - ps) / sample_dist + .499);
-	if (max_high_occ > MAX_MAX_HIGH_OCC)
-		max_high_occ = MAX_MAX_HIGH_OCC;
-	for (j = st, k = 0; j < en && k < max_high_occ; ++j, ++k)
-		b[k] = p->a[j], b[k].pos = j; // b[].pos keeps the index in p->a[]
-	ks_heapmake_mz(k, b); // initialize the binomial heap
-	for (; j < en; ++j) { // if there are more, choose top max_high_occ
-		if (mz_lt(p->a[j], b[0])) { // then update the heap
-			b[0] = p->a[j], b[0].pos = j;
-			ks_heapdown_mz(0, k, b);
-		}
-	}
-	//ks_heapsort_mz(k, b); // sorting is not needed for now
-	for (j = 0; j < k; ++j)
-		if (b[j].rid < pe - ps || force)
-			p->a[b[j].pos].rid = 0;
-}
-
-void select_mz(ha_mz1_v *p, int len, int sample_dist, int32_t dp_min_len)
-{ // for high-occ minimizers, choose up to max_high_occ in each high-occ streak
-	int32_t i, last0 = -1, n = (int32_t)p->n, m = 0, nw[2], min_len;
-	ha_mz1_t b[MAX_MAX_HIGH_OCC]; // this is to avoid a heap allocation
-
-	if (n == 0 || n == 1) return;
-	assert(n < 1<<27); // 27 is the number of bits for ha_mz1_t::pos; this should be safe as there are more bases than minimizers
-	for (i = 0; i < n; ++i)
-		if (p->a[i].rid != 0) ++m;
-	if (m == 0) return; // no high-frequency k-mers; do nothing
-	for (i = 0; i <= n; ++i) {
-		if (i == n || p->a[i].rid == 0) {
-			if (i - last0 > 1) {
-				hf_select(p, last0, i, n, len, sample_dist, b, 0);
-				// int32_t ps = last0 < 0? 0 : p->a[last0].pos;
-				// int32_t pe = i == n? len : p->a[i].pos;
-				// int32_t j, k, st = last0 + 1, en = i;
-				// int32_t max_high_occ = (int32_t)((double)(pe - ps) / sample_dist + .499);
-				// if (max_high_occ > MAX_MAX_HIGH_OCC)
-				// 	max_high_occ = MAX_MAX_HIGH_OCC;
-				// for (j = st, k = 0; j < en && k < max_high_occ; ++j, ++k)
-				// 	b[k] = p->a[j], b[k].pos = j; // b[].pos keeps the index in p->a[]
-				// ks_heapmake_mz(k, b); // initialize the binomial heap
-				// for (; j < en; ++j) { // if there are more, choose top max_high_occ
-				// 	if (mz_lt(p->a[j], b[0])) { // then update the heap
-				// 		b[0] = p->a[j], b[0].pos = j;
-				// 		ks_heapdown_mz(0, k, b);
-				// 	}
-				// }
-				// //ks_heapsort_mz(k, b); // sorting is not needed for now
-				// for (j = 0; j < k; ++j)
-				// 	if (b[j].rid < pe - ps)
-				// 		p->a[b[j].pos].rid = 0;
-			}
-			last0 = i;
-		}
-	}
-
-	min_len = MAX(dp_min_len, (p->a[0].pos+1)+sample_dist);
-	for (i = 0, nw[0] = nw[1] = 0; i < n; i++)
-    {
-		nw[(p->a[i].rid!=0)]++;
-        if((p->a[i].pos + 1) > min_len) break;
-    }
-	if(nw[0]==0 && nw[1]>0) hf_select(p, -1, i, n, len, sample_dist, b, 1);
-
-	min_len = MAX(dp_min_len, (len - (p->a[n-1].pos + 1 - p->a[n-1].span))+sample_dist);
-	for (i = n-1, nw[0] = nw[1] = 0; i >= 0; i--)
-	{
-		nw[(p->a[i].rid!=0)]++;
-		if((len - (p->a[i].pos + 1 - p->a[i].span)) > min_len) break;
-	}
-	if(nw[0]==0 && nw[1]>0) hf_select(p, i, n, n, len, sample_dist, b, 1);
-
-	for (i = n = 0; i < (int32_t)p->n; ++i) // squeeze out filtered minimizers
-		if (p->a[i].rid == 0)
-			p->a[n++] = p->a[i];
-	p->n = n;
-}
-
-static inline int mzcmp_l(const ha_mz1_v *p, int32_t ai, int32_t bi)
-{
-	if(ai >= 0 && bi >= 0){
-		ha_mz1_t *a = &(p->a[ai]), *b = &(p->a[bi]);
-		if(a->rid > 0 && b->rid > 0) return mzcmp(a, b);
-		return (a->rid == 0) - (b->rid == 0);
-	}
-	return (ai < 0) - (bi < 0);
-}
-
-#define GL(x, i) ((int64_t)((uint32_t)((x).a[(i)])))
-#define A_M(p, i) ((i) >= 0 && (p).a[(i)].rid > 0)
-int32_t qfw(ha_mz1_v *p, st_mt_t *mt, int32_t n, int32_t tot_l, int32_t ws, int32_t i, int32_t *mi)
-{
-	int32_t m, si;
-	for (si = i, (*mi) = -1; i < n; i++){
-		if(GL(*mt, i) >= ws || (i+1 < n && GL(*mt, i) < ws && GL(*mt, i+1) > ws) || 
-				(i+1 == n && tot_l >= ws && GL(*mt, i) < ws)){
-			for (m = si; m <= i; m++){
-				if(!A_M(*p, m)) continue;
-				if(mzcmp_l(p, *mi, m) >= 0) (*mi) = m;
-			}
-			if((*mi) >= 0 && A_M(*p, *mi)){
-				for (m = si; m <= i; m++){
-					if(!A_M(*p, m)) continue;
-					if(mzcmp_l(p, *mi, m) == 0) mt->a[m] |= 0x100000000;
-				}
-			}
-			break;
-		}
-	}
-	return i;
-}
-
 void dbg_boundary(ha_mz1_v *p, st_mt_t *mt, int32_t w, int32_t k, int32_t tot_l)
 {
 	if(tot_l < w + k -1) return;
@@ -401,110 +114,6 @@ void dbg_boundary(ha_mz1_v *p, st_mt_t *mt, int32_t w, int32_t k, int32_t tot_l)
 	}
 }
 
-static void select_mz_h(ha_mz1_v *p, st_mt_t *mt, int len, int sample_dist, int32_t w, int32_t k, int32_t tot_l)
-{ // for high-occ minimizers, choose up to max_high_occ in each high-occ streak
-	int32_t i, mi = -1, si, last0 = -1, n = (int32_t)p->n, m = 0, ws = w + k - 1;
-
-	if (n == 0) return;
-	assert(n < 1<<27); // 27 is the number of bits for ha_mz1_t::pos; this should be safe as there are more bases than minimizers
-
-	for (i = m = 0, last0 = -1; i <= n; ++i) {
-		if (i == n || p->a[i].rid == 0) {
-			if (i - last0 > 1) {
-				int32_t ps = last0 < 0? 0 : p->a[last0].pos;
-				int32_t pe = i == n? len : p->a[i].pos;
-				if(((int32_t)((double)(pe - ps) / sample_dist + .499)) > 0){
-					last0 = -2;
-					m++;
-					break;
-				}
-			}
-			last0 = i;
-		}
-	}
-	if (m == 0) return; // no high-frequency k-mers; do nothing
-	if(last0 >= -1) goto ff;
-	i = 0;
-	i = qfw(p, mt, n, tot_l, ws, i, &mi);
-	if(i == n) goto ff;
-	
-	for (si = 0, i++; i < n; i++){
-		for (; si < i; si++){
-			if(GL(*mt, si) + w > GL(*mt, i)) break;
-		}
-
-		// a new minimum; then write the old min
-		if(mzcmp_l(p, i, mi) <= 0) {
-			if(A_M(*p, mi)) mt->a[mi] |= 0x100000000;
-			mi = i;
-		}// old min has moved outside the window
-		else if(si > mi){
-			if(A_M(*p, mi)) mt->a[mi] |= 0x100000000;
-			for (m = si, mi = -1; m <= i; m++){
-				if(mzcmp_l(p, mi, m) >= 0) mi = m;
-            }
-			if(A_M(*p, mi)){
-				for (m = si; m <= i; m++){
-					if(!A_M(*p, m)) continue;
-					if(mzcmp_l(p, mi, m) == 0) mt->a[m] |= 0x100000000;
-				}
-			}
-		}
-	}	
-
-	if(A_M(*p, mi)) mt->a[mi] |= 0x100000000;
-
-	for (i = n - 1; si < n && GL(*mt, si) + w <= tot_l + 1; si++){
-		if(si > mi){
-			if(A_M(*p, mi)) mt->a[mi] |= 0x100000000;
-			for (m = si, mi = -1; m <= i; m++){
-				if(mzcmp_l(p, mi, m) >= 0) mi = m;
-            }
-			if(A_M(*p, mi)){
-				for (m = si; m <= i; m++){
-					if(!A_M(*p, m)) continue;
-					if(mzcmp_l(p, mi, m) == 0) mt->a[m] |= 0x100000000;
-				}
-			}
-		}
-	}
-	
-	/**
-	dbg_boundary(p, mt, w, k, tot_l);
-	fprintf(stderr, "\n");
-	for (i = 0; i < (int32_t)p->n; ++i){
-		if(p->a[i].rid == 0) continue;
-		fprintf(stderr, "%cl: %u, pos: %lu, cnt: %lu, key: %lu, i: %d\n", "+-"[!!(mt->a[i]&0x100000000)], 
-													(uint32_t)mt->a[i], p->a[i].pos, p->a[i].rid, p->a[i].x, i);
-		// if (mt->a[i]&0x100000000){
-		// 	fprintf(stderr, "+l: %u, pos: %lu, cnt: %lu\n", (uint32_t)mt->a[i], p->a[i].pos, p->a[i].rid);
-		// }
-	}
-	**/
-	ha_mz1_t b[MAX_MAX_HIGH_OCC];
-	for (i = 0, last0 = -1; i <= n; ++i) {
-		if (i == n || p->a[i].rid == 0) {
-			if (i - last0 > 1) {
-				int32_t ps = last0 < 0? 0 : p->a[last0].pos;
-				int32_t pe = i == n? len : p->a[i].pos;
-				if(((int32_t)((double)(pe - ps) / sample_dist + .499)) > 0){
-					for (m = last0 + 1, mi = 0; m < i; ++m){
-						if(mt->a[m]&0x100000000) p->a[m].rid = 0, mi++;
-					}
-					if(mi == 0) hf_select(p, last0, i, n, len, sample_dist, b, 0);
-				} 
-			}
-			last0 = i;
-		}
-	}
-
-	ff:
-	for (i = n = 0; i < (int32_t)p->n; ++i) // squeeze out filtered minimizers
-		if (p->a[i].rid == 0)
-			p->a[n++] = p->a[i];
-	p->n = n;
-}
-
 void debug_pl(const char *str, int len, int w, int k, int is_hpc, ha_mz1_v *p, const void *hf, st_mt_t *mt)
 {
 	int i, l, dbi, dbcnt = 0, kmer_span = 0;
@@ -571,153 +180,402 @@ void debug_pl(const char *str, int len, int w, int k, int is_hpc, ha_mz1_v *p, c
 	}
 }
 
-/**
+static inline int mz1_mzcmp(const ha_mz1_t *a, const ha_mz1_t *b){return a->rid < b->rid? -1 : a->rid > b->rid? 1 : ((a->x > b->x) - (a->x < b->x));}
- * Find symmetric (w,k)-minimizers on a DNA sequence
+#define mz1_mz_lt(a, b) (mz1_mzcmp(&(a), &(b)) < 0)
- *
+KSORT_INIT(mz1_mz, ha_mz1_t, mz1_mz_lt)
- * @param str    DNA sequence
- * @param len    length of $str
- * @param w      find a minimizer for every $w consecutive k-mers
- * @param k      k-mer size
- * @param rid    reference ID; will be copied to the output $p array
- * @param is_hpc homopolymer-compressed or not
- * @param p      minimizers
- */
-void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mz1_v *p, const void *hf, int sample_dist, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct, ha_pt_t *pt, int min_freq, int32_t dp_min_len, float dp_e, st_mt_t *mt, int32_t ws, int32_t is_unique)
-{   ///in default, w = 51, k = 51, is_hpc = 1
-    /**
-     uint64_t x;
-     uint64_t rid:28, pos:27, rev:1, span:8;
-     **/
-    extern void *ha_ct_table;
-    static const ha_mz1_t dummy = { UINT64_MAX, (1<<28) - 1, 0, 0, 0};
-    uint64_t shift1 = k - 1, mask = (1ULL<<k) - 1, kmer[4] = {0,0,0,0};
-    int i, j, l, tl = 0, buf_pos, min_pos, kmer_span = 0;
-    ha_mz1_t buf[256], min = dummy;
-    uint32_t buf_p[256], min_s = (uint32_t)-1;
-    tiny_queue_t tq;
 
-    assert(len > 0 && len < 1<<27 && rid < 1<<28 && (w > 0 && w < 256) && (k > 0 && k <= 63));
+static inline int mz2_mzcmp(const ha_mzl_t *a, const ha_mzl_t *b){return a->rid < b->rid? -1 : a->rid > b->rid? 1 : ((a->x > b->x) - (a->x < b->x));}
-    if (dbg_ct != NULL) dbg_ct->a.n = 0;
+#define mz2_mz_lt(a, b) (mz2_mzcmp(&(a), &(b)) < 0)
-    if (k_flag != NULL) {
+KSORT_INIT(mz2_mz, ha_mzl_t, mz2_mz_lt)
-        kv_resize(uint8_t, k_flag->a, (uint64_t)len);
-        k_flag->a.n = len; 
-        memset(k_flag->a.a, 0, k_flag->a.n);
-    } 
 
-    memset(buf, 0xff, w * sizeof(ha_mz1_t));
-    memset(&tq, 0, sizeof(tiny_queue_t));
-    ///len/w is the evaluated minimizer numbers
-    kv_resize(ha_mz1_t, *p, p->n + len/w);
-    kv_resize(uint64_t, *mt, (int64_t)p->m); mt->n = p->n;
 
-    for (i = l = tl = buf_pos = min_pos = 0; i < len; ++i) {
+#define HA_SC_INIT(sf, HType, VType, RidBits, PosBits)\
-        int c = seq_nt4_table[(uint8_t)str[i]];
+inline void sf##_hf_select(VType *p, int32_t si, int32_t ei, int32_t n, int32_t len, int32_t sample_dist, HType *b, int32_t force)\
-        ha_mz1_t info = dummy;
+{\
-        if (c < 4) { // not an ambiguous base
+	if(ei - si <= 1) return;\
-            int z;
+	int32_t ps = si < 0? 0 : p->a[si].pos;\
-            if (is_hpc) {
+	int32_t pe = ei == n? len : p->a[ei].pos;\
-                int skip_len = 1;
+	int32_t j, k, st = si + 1, en = ei;\
-                if (i + 1 < len && seq_nt4_table[(uint8_t)str[i + 1]] == c) {
+	int32_t max_high_occ = (int32_t)((double)(pe - ps) / sample_dist + .499);\
-                    for (skip_len = 2; i + skip_len < len; ++skip_len)
+	if (max_high_occ > MAX_MAX_HIGH_OCC)\
-                        if (seq_nt4_table[(uint8_t)str[i + skip_len]] != c)
+		max_high_occ = MAX_MAX_HIGH_OCC;\
-                            break;
+	for (j = st, k = 0; j < en && k < max_high_occ; ++j, ++k)\
-                    i += skip_len - 1; // put $i at the end of the current homopolymer run
+		b[k] = p->a[j], b[k].pos = j; /** b[].pos keeps the index in p->a[]**/\
-                }
+	ks_heapmake_##sf##_mz(k, b); /** initialize the binomial heap**/\
-                tq_push(&tq, skip_len);
+	for (; j < en; ++j) { /** if there are more, choose top max_high_occ**/\
-                kmer_span += skip_len;
+		if (sf##_mz_lt(p->a[j], b[0])) { /** then update the heap**/\
-                ///how many bases that are covered by this HPC k-mer
+			b[0] = p->a[j], b[0].pos = j;\
-                ///kmer_span includes at most k HPC elements
+			ks_heapdown_##sf##_mz(0, k, b);\
-                if (tq.count > k) kmer_span -= tq_shift(&tq);
+		}\
-            } else kmer_span = l + 1 < k? l + 1 : k; 
+	}\
-            ///kmer_span should be used for HPC k-mer
+	/**ks_heapsort_mz(k, b); // sorting is not needed for now**/\
-            ///non-HPC k-mer, kmer_span should be k
+	for (j = 0; j < k; ++j)\
-            ///kmer_span is used to calculate anchor pos on reverse complementary strand
+		if (b[j].rid < pe - ps || force)\
-
+			p->a[b[j].pos].rid = 0;\
-            if (k_flag != NULL) k_flag->a.a[i] = 1;///lable all useful base, which are not ignored by HPC
+}\
-
+static inline int sf##_mzcmp_l(const VType *p, int32_t ai, int32_t bi)\
-            kmer[0] = (kmer[0] << 1 | (c&1))  & mask;                  // forward k-mer
+{\
-            kmer[1] = (kmer[1] << 1 | (c>>1)) & mask;
+	if(ai >= 0 && bi >= 0){\
-            kmer[2] = kmer[2] >> 1 | (uint64_t)(1 - (c&1))  << shift1; // reverse k-mer
+		HType *a = &(p->a[ai]), *b = &(p->a[bi]);\
-            kmer[3] = kmer[3] >> 1 | (uint64_t)(1 - (c>>1)) << shift1;
+		if(a->rid > 0 && b->rid > 0) return sf##_mzcmp(a, b);\
-            if (kmer[1] == kmer[3]) continue; // skip "symmetric k-mers" as we don't know it strand
+		return (a->rid == 0) - (b->rid == 0);\
-            z = kmer[1] < kmer[3]? 0 : 1; // strand
+	}\
-            ++l; tl++;
+	return (ai < 0) - (bi < 0);\
-            if (l >= k && kmer_span < 256) {
+}\
-                uint64_t y;
+int32_t sf##_qfw(VType *p, st_mt_t *mt, int32_t n, int32_t tot_l, int32_t ws, int32_t i, int32_t *mi)\
-                int32_t cnt, filtered;
+{\
-                y = yak_hash64_64(kmer[z<<1|0]) + yak_hash64_64(kmer[z<<1|1]);
+	int32_t m, si;\
-                cnt = hf? ha_ft_cnt(hf, y) : 0;
+	for (si = i, (*mi) = -1; i < n; i++){\
-                filtered = (cnt >= 1<<28);
+		if(GL(*mt, i) >= ws || (i+1 < n && GL(*mt, i) < ws && GL(*mt, i+1) > ws) || \
-				if(is_unique){
+				(i+1 == n && tot_l >= ws && GL(*mt, i) < ws)){\
-					filtered = (cnt == 0);
+			for (m = si; m <= i; m++){\
-					cnt = (cnt == 1? 0:cnt);
+				if(!A_M(*p, m)) continue;\
-				}
+				if(sf##_mzcmp_l(p, *mi, m) >= 0) (*mi) = m;\
-                if (dbg_ct != NULL) kv_push(uint64_t, dbg_ct->a, ((((uint64_t)(query_ct_index(ha_ct_table, y))<<1)|filtered)<<32)|(uint64_t)(i));
+			}\
-                if (!filtered) info.x = y, info.rid = cnt, info.pos = i, info.rev = z, info.span = kmer_span; // initially ha_mz1_t::rid keeps the k-mer count
+			if((*mi) >= 0 && A_M(*p, *mi)){\
-                if (k_flag != NULL) k_flag->a.a[i]++;
+				for (m = si; m <= i; m++){\
-                if (k_flag != NULL && filtered > 0) k_flag->a.a[i]++;
+					if(!A_M(*p, m)) continue;\
-            }
+					if(sf##_mzcmp_l(p, *mi, m) == 0) mt->a[m] |= 0x100000000;\
-        } else l = 0, tq.count = tq.front = 0, kmer_span = 0;
+				}\
-
+			}\
-        buf[buf_pos] = info; // need to do this here as appropriate buf_pos and buf[buf_pos] are needed below
+			break;\
-        buf_p[buf_pos] = l;
+		}\
-        if (l == w + k - 1 && min.x != UINT64_MAX) { // special case for the first window - because identical k-mers are not stored yet
+	}\
-			for (j = buf_pos + 1; j < w; ++j){
+	return i;\
-                if (mzcmp(&min, &buf[j]) == 0 && buf[j].pos != min.pos){
+}\
-                    kv_push(ha_mz1_t, *p, buf[j]); kv_push(uint64_t, *mt, buf_p[j]);
+static void sf##_select_mz_h(VType *p, st_mt_t *mt, int len, int sample_dist, int32_t w, int32_t k, int32_t tot_l)\
-                }
+{ /**for high-occ minimizers, choose up to max_high_occ in each high-occ streak**/\
-            }
+	int32_t i, mi = -1, si, last0 = -1, n = (int32_t)p->n, m = 0, ws = w + k - 1;\
-            for (j = 0; j < buf_pos; ++j){
+	if (n == 0) return;\
-                if (mzcmp(&min, &buf[j]) == 0 && buf[j].pos != min.pos){
+	assert((int64_t)(n) < (int64_t)((((uint64_t)1)<<PosBits)));\
-                    kv_push(ha_mz1_t, *p, buf[j]); kv_push(uint64_t, *mt, buf_p[j]);
+	for (i = m = 0, last0 = -1; i <= n; ++i) {\
-                } 
+		if (i == n || p->a[i].rid == 0) {\
-            }
+			if (i - last0 > 1) {\
-        }
+				int32_t ps = last0 < 0? 0 : p->a[last0].pos;\
-
+				int32_t pe = i == n? len : p->a[i].pos;\
-        /**
+				if(((int32_t)((double)(pe - ps) / sample_dist + .499)) > 0){\
-         * There are three cases:
+					last0 = -2;\
-         * 1. info.x <= min.x, means info is a new minimizer
+					m++;\
-         * 2. info.x > min.x, info is not a new minimizer
+					break;\
-         *    (1) buf_pos != min_pos, do nothing
+				}\
-         *    (2) buf_pos == min_pos, means current minimizer has moved outside the window
+			}\
-         * **/
+			last0 = i;\
-        ///three cases: 1. 
+		}\
-        if (mzcmp(&min, &info) >= 0) { // a new minimum; then write the old min
+	}\
-            if (l >= w + k && min.x != UINT64_MAX){
+	if (m == 0) return; /**no high-frequency k-mers; do nothing**/\
-                kv_push(ha_mz1_t, *p, min); kv_push(uint64_t, *mt, min_s);
+	if(last0 >= -1) goto sf##_ff;\
-            } 
+	i = 0;\
-            min = info, min_pos = buf_pos, min_s = buf_p[buf_pos];
+	i = sf##_qfw(p, mt, n, tot_l, ws, i, &mi);\
-        } else if (buf_pos == min_pos) { // old min has moved outside the window
+	if(i == n) goto sf##_ff;\
-            if (l >= w + k - 1 && min.x != UINT64_MAX){
+	for (si = 0, i++; i < n; i++){\
-                kv_push(ha_mz1_t, *p, min); kv_push(uint64_t, *mt, min_s);
+		for (; si < i; si++){\
-            } 
+			if(GL(*mt, si) + w > GL(*mt, i)) break;\
-            ///buf_pos == min_pos, means current minimizer has moved outside the window
+		}\
-            ///so for now we need to find a new minimizer at the current window (w k-mers)
+		/**a new minimum; then write the old min**/\
-            for (j = buf_pos + 1, min = dummy; j < w; ++j) // the two loops are necessary when there are identical k-mers
+		if(sf##_mzcmp_l(p, i, mi) <= 0) {\
-                if (mzcmp(&min, &buf[j]) >= 0) min = buf[j], min_pos = j, min_s = buf_p[j]; // >= is important s.t. min is always the closest k-mer
+			if(A_M(*p, mi)) mt->a[mi] |= 0x100000000;\
-            for (j = 0; j <= buf_pos; ++j)
+			mi = i;\
-                if (mzcmp(&min, &buf[j]) >= 0) min = buf[j], min_pos = j, min_s = buf_p[j];
+		}/**old min has moved outside the window**/\
-
+		else if(si > mi){\
-            if (l >= w + k - 1 && min.x != UINT64_MAX) { // write identical k-mers
+			if(A_M(*p, mi)) mt->a[mi] |= 0x100000000;\
-                for (j = buf_pos + 1; j < w; ++j) // these two loops make sure the output is sorted
+			for (m = si, mi = -1; m <= i; m++){\
-                    if (mzcmp(&min, &buf[j]) == 0 && min.pos != buf[j].pos){
+				if(sf##_mzcmp_l(p, mi, m) >= 0) mi = m;\
-                        kv_push(ha_mz1_t, *p, buf[j]); kv_push(uint64_t, *mt, buf_p[j]);
+            }\
-                    }
+			if(A_M(*p, mi)){\
-                for (j = 0; j <= buf_pos; ++j)
+				for (m = si; m <= i; m++){\
-                    if (mzcmp(&min, &buf[j]) == 0 && min.pos != buf[j].pos){
+					if(!A_M(*p, m)) continue;\
-                        kv_push(ha_mz1_t, *p, buf[j]); kv_push(uint64_t, *mt, buf_p[j]);
+					if(sf##_mzcmp_l(p, mi, m) == 0) mt->a[m] |= 0x100000000;\
-                    } 
+				}\
-            }
+			}\
-        }
+		}\
-        if (++buf_pos == w) buf_pos = 0;
+	}\
-    }
+	if(A_M(*p, mi)) mt->a[mi] |= 0x100000000;\
-    if (min.x != UINT64_MAX){
+	for (i = n - 1; si < n && GL(*mt, si) + w <= tot_l + 1; si++){\
-        kv_push(ha_mz1_t, *p, min); kv_push(uint64_t, *mt, min_s); 
+		if(si > mi){\
-    }
+			if(A_M(*p, mi)) mt->a[mi] |= 0x100000000;\
-	// debug_pl(str, len, w, k, is_hpc, p, hf, mt);
+			for (m = si, mi = -1; m <= i; m++){\
-    if (sample_dist > w) select_mz_h(p, mt, len, sample_dist, ws, k, tl);
+				if(sf##_mzcmp_l(p, mi, m) >= 0) mi = m;\
-    if (dp_min_len > 0 && pt && mt) refine_sketch(p, pt, len, dp_min_len, dp_e, min_freq, mt);
+            }\
-	for (i = 0; i < (int)p->n; ++i) // populate .rid as this was keeping counts
+			if(A_M(*p, mi)){\
-		p->a[i].rid = rid;
+				for (m = si; m <= i; m++){\
+					if(!A_M(*p, m)) continue;\
+					if(sf##_mzcmp_l(p, mi, m) == 0) mt->a[m] |= 0x100000000;\
+				}\
+			}\
+		}\
+	}\
+	/**dbg_boundary(p, mt, w, k, tot_l);**/\
+	HType b[MAX_MAX_HIGH_OCC];\
+	for (i = 0, last0 = -1; i <= n; ++i) {\
+		if (i == n || p->a[i].rid == 0) {\
+			if (i - last0 > 1) {\
+				int32_t ps = last0 < 0? 0 : p->a[last0].pos;\
+				int32_t pe = i == n? len : p->a[i].pos;\
+				if(((int32_t)((double)(pe - ps) / sample_dist + .499)) > 0){\
+					for (m = last0 + 1, mi = 0; m < i; ++m){\
+						if(mt->a[m]&0x100000000) p->a[m].rid = 0, mi++;\
+					}\
+					if(mi == 0) sf##_hf_select(p, last0, i, n, len, sample_dist, b, 0);\
+				}\
+			}\
+			last0 = i;\
+		}\
+	}\
+	sf##_ff:\
+	for (i = n = 0; i < (int32_t)p->n; ++i) /**squeeze out filtered minimizers**/\
+		if (p->a[i].rid == 0)\
+			p->a[n++] = p->a[i];\
+	p->n = n;\
+}\
+void sf##_refine_select(VType *mz, int32_t sidx, int32_t eidx, int32_t sn, int32_t min_freq, st_mt_t *mm, int32_t *rsi, int32_t *rei)\
+{\
+	int32_t n = sn, m = eidx + 1 - sidx, i, k, t, mk=-1;\
+	uint64_t ix, kx, ks;\
+	kv_resize(uint64_t, *mm, mm->n+n*m);\
+	HType *ma = mz->a + sidx;\
+	uint64_t *mmt = mm->a + mm->n;\
+	/**fprintf(stderr, "[M::%s::] ==> +n: %d, m: %d, sn: %d, sidx: %d, eidx: %d\n", __func__, n, m, sn, sidx, eidx);**/\
+	for (i = 0; i < n; i++) /**how many selected minimizers**/\
+	{\
+		for (k = 0, mk = -1; k < m; k++) /**how many minimizers in total**/\
+		{\
+			if((int32_t)(ma[k].rid)<min_freq) continue;\
+			ks = ma[k].pos + 1 - ma[k].span; t = -1;\
+			if(i > 0){\
+				for (t = k-1; t >= 0 && (ma[t].pos >= ks||(int32_t)(ma[t].rid)<min_freq); t--);\
+			}\
+			ix = (i <= 0?0:(t<0?0xffffffff:(GMC(mmt, i-1,t,m)&0xffffffff)));\
+			if(ix < 0xffffffff) ix += (ma[k].rid);\
+			kx = (mk < 0?0xffffffff:(GMC(mmt, i, mk,m)&0xffffffff));\
+			ks = MIN(ix, kx);\
+			/**fprintf(stderr, "ks: %lu, i: %d (n-%d), k: %d (m-%d), ix: %lu, kx: %lu, t: %d, mk: %d\n", ks, i, n, k, m, ix, kx, t, mk);**/\
+			if((ks&0xffffffff) == 0xffffffff) ks |= ((uint64_t)0xffffffff)<<32;\
+			else if(ks == ix) ks |= (uint64_t)(i>0?(i-1)*m+t:0xffffffff)<<32;\
+			else if(ks == kx) ks |= (uint64_t)(mk>=0?i*m+mk:0xffffffff)<<32;\
+			GMC(mmt, i,k,m) = ks;\
+			mk = k;\
+		}\
+	}\
+	/**fprintf(stderr, "[M::%s::] ==> ++n: %d, m: %d, sn: %d, sidx: %d, eidx: %d\n", __func__, n, m, sn, sidx, eidx);**/\
+	ks = (n-1)*m + mk; ix = (uint64_t)-1; kx = 0;\
+	while (ks != 0xffffffff)\
+	{\
+		i = ks/m; k = ks%m;\
+		ks = mmt[ks]>>32;\
+		/**fprintf(stderr, "i: %d, k: %d, ks: %lu\n", i, k, ks);**/\
+		if(ks == 0xffffffff || (int32_t)(ks/m) == (i-1)){\
+			mm->a[sidx+k] = 1;\
+			ix = MIN((uint64_t)k, ix); kx = MAX((uint64_t)k, kx);\
+		}\
+	}\
+	/**debug_refine(ma, mmt, sn, n, m, (n-1)*m + mk);**/\
+	if(rsi) (*rsi) = ix + sidx;\
+	if(rei) (*rei) = kx + sidx;\
+}\
+void sf##_refine_sketch(VType *p, ha_pt_t *pt, int32_t rlen, int32_t dp_min_len, float er, int32_t min_freq, st_mt_t *mt)\
+{\
+	/**fprintf(stderr, "[M::%s::] ==> #########10#########, rlen: %d\n", __func__, rlen);**/\
+	int32_t i, n = p->n, bd, len = MIN(rlen, dp_min_len), sublen, cnt, ei, li, ri;\
+	int32_t sn =  len*er + 1;\
+	kv_resize(uint64_t, *mt, (int64_t)p->n);\
+	mt->n = p->n; memset(mt->a, 0, sizeof(uint64_t)*p->n);\
+	for (i = 0; i < n; i++) p->a[i].rid = ha_pt_cnt(pt, p->a[i].x);\
+	for (i = cnt = 0, bd = -1, ei = -1; i < n; i++){\
+		if((int32_t)(p->a[i].rid)<min_freq) continue;\
+		sublen = p->a[i].pos + 1;\
+		if(sublen > len) break;\
+		else ei = i;\
+		if((int32_t)(p->a[i].pos + 1 - p->a[i].span) > bd){\
+			bd = p->a[i].pos;\
+			cnt++;\
+		}\
+	}\
+	/**fprintf(stderr, "[M::%s::] ==> +cnt: %d, sn: %d, ei: %d, n: %d\n", __func__, cnt, sn, ei, n);**/\
+	if(cnt >= sn) sf##_refine_select(p, 0, ei, sn, min_freq, mt, NULL, &li);\
+	else{\
+		li = i-1;\
+		for (i = 0; i <= li; i++) mt->a[i] = 1;\
+	}\
+	if(len < rlen){\
+		for (i = n-1, cnt = 0, bd = rlen+1, ei = -1; i >= 0; i--){\
+			if((int32_t)(p->a[i].rid)<min_freq) continue;\
+			sublen = rlen - (p->a[i].pos + 1 - p->a[i].span);\
+			if(sublen > len) break;\
+			else ei = i;\
+			if((int32_t)(p->a[i].pos) < bd){\
+				bd = p->a[i].pos + 1 - p->a[i].span;\
+				cnt++;\
+			}\
+		}\
+		/**fprintf(stderr, "[M::%s::] ==> -cnt: %d, sn: %d, ei: %d, n: %d\n", __func__, cnt, sn, ei, n);**/\
+		if(cnt >= sn) sf##_refine_select(p, ei, n-1, sn, min_freq, mt, &ri, NULL);\
+		else {\
+			ri = i+1;\
+			for (i = ri; i <= n-1; i++) mt->a[i] = 1;\
+		}\
+		/**fprintf(stderr, "[M::%s::] ==> --cnt: %d, sn: %d, ei: %d, n: %d\n", __func__, cnt, sn, ei, n);**/\
+		if(ri - li >= 2){\
+			li++; ri--;\
+			sn =  (p->a[ri].pos - p->a[li].pos + p->a[li].span)*er + 1;\
+			for (i = li, cnt = 0, bd = -1; i <= ri; i++){\
+				if((int32_t)(p->a[i].rid)<min_freq) continue;\
+				if((int32_t)(p->a[i].pos + 1 - p->a[i].span) > bd){\
+					bd = p->a[i].pos;\
+					cnt++;\
+					if(cnt >= sn) break;\
+				}\
+			}\
+			if(cnt >= sn) sf##_refine_select(p, li, ri, sn, min_freq, mt, NULL, NULL);\
+			else for (i = li; i <= ri; i++) mt->a[i] = 1;\
+		}\
+	}\
+	/**fprintf(stderr, "[M::%s::] ==> #########20#########, p->n: %u, n: %d\n", __func__, p->n, n);**/\
+	for (i = sn = 0; i < n; i++){\
+		if(mt->a[i]){\
+			p->a[sn] = p->a[i];\
+			sn++;\
+		}\
+	}\
+	/**if(p->n != sn) fprintf(stderr, "[M::%s::] ==> #########21#########, p->n: %u, sn: %d\n", __func__, p->n, sn);**/\
+	p->n = sn;\
+}\
+/**\
+ * Find symmetric (w,k)-minimizers on a DNA sequence\
+ *\
+ * @param str    DNA sequence\
+ * @param len    length of $str\
+ * @param w      find a minimizer for every $w consecutive k-mers\
+ * @param k      k-mer size\
+ * @param rid    reference ID; will be copied to the output $p array\
+ * @param is_hpc homopolymer-compressed or not\
+ * @param p      minimizers\
+ */\
+void sf##_ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, VType *p, const void *hf, int sample_dist, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct, ha_pt_t *pt, int min_freq, int32_t dp_min_len, float dp_e, st_mt_t *mt, int32_t ws, int32_t is_unique)\
+{   /**in default, w = 51, k = 51, is_hpc = 1**/\
+    extern void *ha_ct_table;\
+    static const HType dummy = { UINT64_MAX, (((uint64_t)1)<<RidBits) - 1, 0, 0, 0};\
+    uint64_t shift1 = k - 1, mask = (1ULL<<k) - 1, kmer[4] = {0,0,0,0};\
+    int i, j, l, tl = 0, buf_pos, min_pos, kmer_span = 0;\
+    HType buf[256], min = dummy;\
+    uint32_t buf_p[256], min_s = (uint32_t)-1;\
+    tiny_queue_t tq;\
+    assert(len > 0 && (int64_t)(len) < (int64_t)((((uint64_t)1)<<PosBits)) && (int64_t)(rid) < (int64_t)((((uint64_t)1)<<RidBits)) && (w > 0 && w < 256) && (k > 0 && k <= 63));\
+    if (dbg_ct != NULL) dbg_ct->a.n = 0;\
+    if (k_flag != NULL) {\
+        kv_resize(uint8_t, k_flag->a, (uint64_t)len);\
+        k_flag->a.n = len;\
+        memset(k_flag->a.a, 0, k_flag->a.n);\
+    }\
+    memset(buf, 0xff, w * sizeof(HType));\
+    memset(&tq, 0, sizeof(tiny_queue_t));\
+    /**len/w is the evaluated minimizer numbers**/\
+    kv_resize(HType, *p, p->n + len/w);\
+    kv_resize(uint64_t, *mt, (int64_t)p->m); mt->n = p->n;\
+    for (i = l = tl = buf_pos = min_pos = 0; i < len; ++i) {\
+        int c = seq_nt4_table[(uint8_t)str[i]];\
+        HType info = dummy;\
+        if (c < 4) { /**not an ambiguous base**/\
+            int z;\
+            if (is_hpc) {\
+                int skip_len = 1;\
+                if (i + 1 < len && seq_nt4_table[(uint8_t)str[i + 1]] == c) {\
+                    for (skip_len = 2; i + skip_len < len; ++skip_len)\
+                        if (seq_nt4_table[(uint8_t)str[i + skip_len]] != c)\
+                            break;\
+                    i += skip_len - 1; /**put $i at the end of the current homopolymer run**/\
+                }\
+                tq_push(&tq, skip_len);\
+                kmer_span += skip_len;\
+                /**how many bases that are covered by this HPC k-mer\
+                kmer_span includes at most k HPC elements**/\
+                if (tq.count > k) kmer_span -= tq_shift(&tq);\
+            } else kmer_span = l + 1 < k? l + 1 : k;\
+            /**kmer_span should be used for HPC k-mer\
+            non-HPC k-mer, kmer_span should be k\
+            kmer_span is used to calculate anchor pos on reverse complementary strand**/\
+            if (k_flag != NULL) k_flag->a.a[i] = 1;/**lable all useful base, which are not ignored by HPC**/\
+            kmer[0] = (kmer[0] << 1 | (c&1))  & mask;/**forward k-mer**/\
+            kmer[1] = (kmer[1] << 1 | (c>>1)) & mask;\
+            kmer[2] = kmer[2] >> 1 | (uint64_t)(1 - (c&1))  << shift1; /**reverse k-mer**/\
+            kmer[3] = kmer[3] >> 1 | (uint64_t)(1 - (c>>1)) << shift1;\
+            if (kmer[1] == kmer[3]) continue; /** skip "symmetric k-mers" as we don't know it strand**/\
+            z = kmer[1] < kmer[3]? 0 : 1; /** strand**/\
+            ++l; tl++;\
+            if (l >= k && kmer_span < 256) {\
+                uint64_t y;\
+                int32_t cnt, filtered;\
+                y = yak_hash64_64(kmer[z<<1|0]) + yak_hash64_64(kmer[z<<1|1]);\
+                cnt = hf? ha_ft_cnt(hf, y) : 0;\
+                filtered = (cnt >= 1<<28);\
+				if(is_unique){\
+					filtered = (cnt < is_unique);\
+					cnt = (cnt == is_unique? 0:cnt);\
+				}\
+                if (dbg_ct != NULL) kv_push(uint64_t, dbg_ct->a, ((((uint64_t)(query_ct_index(ha_ct_table, y))<<1)|filtered)<<32)|(uint64_t)(i));\
+                if (!filtered) info.x = y, info.rid = cnt, info.pos = i, info.rev = z, info.span = kmer_span; /** initially ha_mz1_t::rid keeps the k-mer count**/\
+                if (k_flag != NULL) k_flag->a.a[i]++;\
+                if (k_flag != NULL && filtered > 0) k_flag->a.a[i]++;\
+            }\
+        } else l = 0, tq.count = tq.front = 0, kmer_span = 0;\
+        buf[buf_pos] = info; /**need to do this here as appropriate buf_pos and buf[buf_pos] are needed below**/\
+        buf_p[buf_pos] = l;\
+        if (l == w + k - 1 && min.x != UINT64_MAX) { /**special case for the first window - because identical k-mers are not stored yet**/\
+			for (j = buf_pos + 1; j < w; ++j){\
+                if (sf##_mzcmp(&min, &buf[j]) == 0 && buf[j].pos != min.pos){\
+                    kv_push(HType, *p, buf[j]); kv_push(uint64_t, *mt, buf_p[j]);\
+                }\
+            }\
+            for (j = 0; j < buf_pos; ++j){\
+                if (sf##_mzcmp(&min, &buf[j]) == 0 && buf[j].pos != min.pos){\
+                    kv_push(HType, *p, buf[j]); kv_push(uint64_t, *mt, buf_p[j]);\
+                }\
+            }\
+        }\
+        /**\
+         * There are three cases:\
+         * 1. info.x <= min.x, means info is a new minimizer\
+         * 2. info.x > min.x, info is not a new minimizer\
+         *    (1) buf_pos != min_pos, do nothing\
+         *    (2) buf_pos == min_pos, means current minimizer has moved outside the window\
+         * **/\
+        /**three cases: 1.**/\
+        if (sf##_mzcmp(&min, &info) >= 0) { /**a new minimum; then write the old min**/\
+            if (l >= w + k && min.x != UINT64_MAX){\
+                kv_push(HType, *p, min); kv_push(uint64_t, *mt, min_s);\
+            }\
+            min = info, min_pos = buf_pos, min_s = buf_p[buf_pos];\
+        } else if (buf_pos == min_pos) { /**old min has moved outside the window**/\
+            if (l >= w + k - 1 && min.x != UINT64_MAX){\
+                kv_push(HType, *p, min); kv_push(uint64_t, *mt, min_s);\
+            }\
+            /**buf_pos == min_pos, means current minimizer has moved outside the window\
+            so for now we need to find a new minimizer at the current window (w k-mers)**/\
+            for (j = buf_pos + 1, min = dummy; j < w; ++j) /**the two loops are necessary when there are identical k-mers**/\
+                if (sf##_mzcmp(&min, &buf[j]) >= 0) min = buf[j], min_pos = j, min_s = buf_p[j]; /** >= is important s.t. min is always the closest k-mer**/\
+            for (j = 0; j <= buf_pos; ++j)\
+                if (sf##_mzcmp(&min, &buf[j]) >= 0) min = buf[j], min_pos = j, min_s = buf_p[j];\
+            if (l >= w + k - 1 && min.x != UINT64_MAX) { /**write identical k-mers**/\
+                for (j = buf_pos + 1; j < w; ++j) /**these two loops make sure the output is sorted**/\
+                    if (sf##_mzcmp(&min, &buf[j]) == 0 && min.pos != buf[j].pos){\
+                        kv_push(HType, *p, buf[j]); kv_push(uint64_t, *mt, buf_p[j]);\
+                    }\
+                for (j = 0; j <= buf_pos; ++j)\
+                    if (sf##_mzcmp(&min, &buf[j]) == 0 && min.pos != buf[j].pos){\
+                        kv_push(HType, *p, buf[j]); kv_push(uint64_t, *mt, buf_p[j]);\
+                    }\
+            }\
+        }\
+        if (++buf_pos == w) buf_pos = 0;\
+    }\
+    if (min.x != UINT64_MAX){\
+        kv_push(HType, *p, min); kv_push(uint64_t, *mt, min_s);\
+    }\
+	/**debug_pl(str, len, w, k, is_hpc, p, hf, mt);**/\
+    if (sample_dist > w) sf##_select_mz_h(p, mt, len, sample_dist, ws, k, tl);\
+    if (dp_min_len > 0 && pt && mt) sf##_refine_sketch(p, pt, len, dp_min_len, dp_e, min_freq, mt);\
+	for (i = 0; i < (int)p->n; ++i) /**populate .rid as this was keeping counts**/\
+		p->a[i].rid = rid;\
 }
 
+HA_SC_INIT(mz1, ha_mz1_t, ha_mz1_v, 28, 27)
+HA_SC_INIT(mz2, ha_mzl_t, ha_mzl_v, 31, 32)