disable debug

Merge branch 'master' of https://github.com/chhylp123/Long_read_assembly
Update README
2026-09-16 05:07:55 +08:00 · 2020-03-23 00:25:20 -04:00 · 2020-03-22 18:32:47 -04:00 · 2020-03-22 18:31:59 -04:00 · 2020-03-22 18:01:04 -04:00 · 2020-03-22 17:06:36 -04:00
16 changed files with 12362 additions and 803 deletions
--- a/Assembly.cpp
+++ b/Assembly.cpp
@@ -9,6 +9,7 @@
 #include "POA.h"
 #include "Correct.h"
 #include "Output.h"
+#include "Trio.h"

 Total_Count_Table TCB;
 Total_Pos_Table PCB;
@@ -1257,7 +1258,7 @@ void Overlap_calculate_multipe_thr()
    double start_time = Get_T();
    
    fprintf(stderr, "Begin calculating overlaps... \n");
-
+    
    pthread_t *_r_threads;

 	_r_threads = (pthread_t *)malloc(sizeof(pthread_t)*asm_opt.thread_num);
@@ -1327,19 +1328,6 @@ void Overlap_calculate_multipe_thr()
 }


-// int load_pre_cauculated_index()
-// {
-//     if(load_Total_Pos_Table(&PCB, asm_opt.read_file_name) && load_All_reads(&R_INF, asm_opt.read_file_name))
-//     {
-//         return 1;
-//     }
-//     else
-//     {
-//         return 0;
-//     }
-    
-// }
-

 void update_overlaps(overlap_region_alloc* overlap_list, ma_hit_t_alloc* paf, 
 UC_Read* g_read, UC_Read* overlap_read, int is_match, int is_exact)
@@ -1758,7 +1746,7 @@ void Output_PAF()
    free(paf_name);
    fclose(output_file);

-    fprintf(stderr, "PAF has beem written.\n");
+    fprintf(stderr, "PAF has been written.\n");
 }


@@ -1838,9 +1826,6 @@ long long readNum, long long rescue_threshold, float cluster_threshold)
        if((long long)edge_vector.n >= rescue_threshold && 
        check_cluster(edge_vector.a, edge_vector.n, paf, cluster_threshold) == 1)
        {
-            // fprintf(stderr,"\nremove following %u edges...\n", edge_vector.n);
-            // print_revise_edges(&(paf[i]), edge_vector_index.a, edge_vector_index.n);
-
            add_overlaps(&(paf[i]), &(rev_paf[i]), edge_vector_index.a, edge_vector_index.n);
            remove_overlaps(&(paf[i]), edge_vector_index.a, edge_vector_index.n);
            revises = revises + edge_vector.n;
@@ -1868,10 +1853,6 @@ long long readNum, long long rescue_threshold, float cluster_threshold)
        if((long long)edge_vector.n >= rescue_threshold && 
        check_cluster(edge_vector.a, edge_vector.n, paf, cluster_threshold) == 1)
        {
-            // fprintf(stderr,"\nadd following %u edges...\n", edge_vector.n);
-            // print_revise_edges(&(rev_paf[i]), edge_vector_index.a, edge_vector_index.n);
-
-
            remove_overlaps(&(rev_paf[i]), edge_vector_index.a, edge_vector_index.n);
            add_overlaps_from_different_sources(paf, &(paf[i]), edge_vector.a, edge_vector.n);
            revises = revises + edge_vector.n;
@@ -1886,6 +1867,7 @@ long long readNum, long long rescue_threshold, float cluster_threshold)
 }


+
 void generate_overlaps(int last_round)
 {
    double start_time = Get_T();
@@ -1914,11 +1896,14 @@ void generate_overlaps(int last_round)
    free(_r_threads);

    ///rescue_edges(R_INF.paf, R_INF.reverse_paf, R_INF.total_reads, 4, 0.985);
+    destory_Total_Pos_Table(&PCB);

    fprintf(stderr, "Final overlaps have been calculated.\n");
    fprintf(stderr, "%-30s%18.2f\n\n", "Final overlaps calculation time:", Get_T() - start_time); 

    Output_PAF();
+
+    trio_partition();
    
    build_string_graph_without_clean(asm_opt.min_overlap_coverage, R_INF.paf, R_INF.reverse_paf, 
    R_INF.total_reads, R_INF.read_length, asm_opt.min_overlap_Len, asm_opt.max_hang_Len, asm_opt.clean_round, 
--- a/Assembly.h
+++ b/Assembly.h
@@ -10,7 +10,6 @@

 void Counting_multiple_thr();
 void Build_hash_table_multiple_thr();
-///int load_pre_cauculated_index();
 void Overlap_calculate_multipe_thr();
 void Correct_Reads(int last_round);
 #endif
--- a/CommandLines.cpp
+++ b/CommandLines.cpp
@@ -5,7 +5,7 @@
 #include "ketopt.h"
 #include <sys/time.h>

-#define VERSION "0.2.0"
+#define VERSION "0.3.0"
 #define DEFAULT_OUTPUT "hifiasm.asm"

 hifiasm_opt_t asm_opt;
@@ -21,7 +21,9 @@ void Print_H(hifiasm_opt_t* asm_opt)
 {
    fprintf(stderr, "Usage: hifiasm [options] <in_1.fq> <in_2.fq> <...>\n");
    fprintf(stderr, "Options:\n");
+    fprintf(stderr, "  Assembly:\n");
    fprintf(stderr, "    -o FILE       prefix of output files [%s]\n", asm_opt->output_file_name);
+    ///fprintf(stderr, "    -c FILE       file including trio information\n");
    fprintf(stderr, "    -t INT        number of threads [%d]\n", asm_opt->thread_num);
    fprintf(stderr, "    -r INT        round of correction [%d]\n", asm_opt->number_of_round);
    fprintf(stderr, "    -a INT        round of assembly cleaning [%d]\n", asm_opt->clean_round);
@@ -40,6 +42,13 @@ void Print_H(hifiasm_opt_t* asm_opt)
    fprintf(stderr, "    -y FLOAT      min overlap drop ratio [%.2g]\n", asm_opt->min_drop_rate);
    fprintf(stderr, "    -v            show version number\n");
    fprintf(stderr, "    -h            show help information\n");
+
+    fprintf(stderr, "  Trio-partition:\n");
+    fprintf(stderr, "    -P FILE       paternal trio index generated by \"yak count\" [NULL]\n");
+    fprintf(stderr, "    -M FILE       maternal trio index generated by \"yak count\" [NULL]\n");
+    fprintf(stderr, "    -c INT        lower bound of the binned k-mer's frequency [%d]\n", asm_opt->min_cnt);
+    fprintf(stderr, "    -d INT        upper bound of the binned k-mer's frequency [%d]\n", asm_opt->mid_cnt);
+
    fprintf(stderr, "Example: ./hifiasm -o NA12878.asm -t 32 NA12878.fq.gz\n");
    fprintf(stderr, "See `man ./hifiasm.1' for detailed description of these command-line options.\n");
 }
@@ -51,6 +60,8 @@ void init_opt(hifiasm_opt_t* asm_opt)
    asm_opt->read_file_names = NULL;
    asm_opt->output_file_name = (char*)(DEFAULT_OUTPUT);
    asm_opt->required_read_name = NULL;
+    asm_opt->pat_index = NULL;
+    asm_opt->mat_index = NULL;
    asm_opt->thread_num = 1;
    asm_opt->k_mer_length = 40;
    asm_opt->k_mer_min_freq = 3;
@@ -71,6 +82,8 @@ void init_opt(hifiasm_opt_t* asm_opt)
    asm_opt->min_overlap_Len = 50;
    asm_opt->min_overlap_coverage = 0;
    asm_opt->max_short_tip = 3;
+    asm_opt->min_cnt = 2;
+    asm_opt->mid_cnt = 5;
 }

 void destory_opt(hifiasm_opt_t* asm_opt)
@@ -90,6 +103,24 @@ void clear_opt(hifiasm_opt_t* asm_opt, int last_round)
    asm_opt->roundID = asm_opt->number_of_round - last_round;
 }

+int check_file(char* name, const char* opt)
+{
+    if(!name)
+    {
+        fprintf(stderr, "[ERROR] file does not exist (-%s)\n", opt);
+        return 0;
+    } 
+    FILE* is_exist = NULL;
+    is_exist = fopen(name,"r");
+    if(!is_exist)
+    {
+        fprintf(stderr, "[ERROR] %s does not exist (-%s)\n", name, opt);
+        return 0;
+    } 
+
+    fclose(is_exist);
+    return 1;
+}

 int check_option(hifiasm_opt_t* asm_opt)
 {
@@ -206,6 +237,10 @@ int check_option(hifiasm_opt_t* asm_opt)
        return 0;
    }

+
+    if(asm_opt->pat_index != NULL && check_file(asm_opt->pat_index, "P") == 0) return 0;
+    if(asm_opt->mat_index != NULL && check_file(asm_opt->mat_index, "M") == 0) return 0;
+
    // fprintf(stderr, "input file num: %d\n", asm_opt->num_reads);
    // fprintf(stderr, "output file: %s\n", asm_opt->output_file_name);
    // fprintf(stderr, "number of threads: %d\n", asm_opt->thread_num);
@@ -218,6 +253,11 @@ int check_option(hifiasm_opt_t* asm_opt)
    // fprintf(stderr, "size of popped small bubbles: %lld\n", asm_opt->small_pop_bubble_size);
    // fprintf(stderr, "size of popped large bubbles: %lld\n", asm_opt->large_pop_bubble_size);
    // fprintf(stderr, "small removed unitig threshold: %d\n", asm_opt->max_short_tip);
+    // fprintf(stderr, "small removed unitig threshold: %d\n", asm_opt->max_short_tip);
+    // fprintf(stderr, "pat_index: %s\n", asm_opt->pat_index);
+    // fprintf(stderr, "mat_index: %s\n", asm_opt->mat_index);
+    // fprintf(stderr, "min_cnt: %d\n", asm_opt->min_cnt);
+    // fprintf(stderr, "mid_cnt: %d\n", asm_opt->mid_cnt);

    return 1;
 }
@@ -255,7 +295,7 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)

    int c;

-    while ((c = ketopt(&opt, argc, argv, 1, "hvt:o:k:lwm:n:r:a:b:z:x:y:p:i", 0)) >= 0) {
+    while ((c = ketopt(&opt, argc, argv, 1, "hvt:o:k:lwm:n:r:a:b:z:x:y:p:c:d:M:P:i", 0)) >= 0) {
        if (c == 'h')
        {
            Print_H(asm_opt);
@@ -276,6 +316,10 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
        else if (c == 'a') asm_opt->clean_round = atoi(opt.arg); 
        else if (c == 'z') asm_opt->adapterLen = atoi(opt.arg);
        else if (c == 'b') asm_opt->required_read_name = opt.arg;
+        else if (c == 'c') asm_opt->min_cnt = atoi(opt.arg);
+        else if (c == 'd') asm_opt->mid_cnt = atoi(opt.arg);
+        else if (c == 'P') asm_opt->pat_index = opt.arg;
+        else if (c == 'M') asm_opt->mat_index = opt.arg;
        else if (c == 'x') asm_opt->max_drop_rate = atof(opt.arg);
        else if (c == 'y') asm_opt->min_drop_rate = atof(opt.arg);
        else if (c == 'p') asm_opt->small_pop_bubble_size = atoll(opt.arg);
--- a/CommandLines.h
+++ b/CommandLines.h
@@ -4,12 +4,15 @@
 #include <pthread.h>

 #define VERBOSE 0
+#define VERBOSE_GFA 0

 typedef struct {
    int num_reads;
    char** read_file_names;
    char* output_file_name;
    char* required_read_name;
+    char* pat_index;
+    char* mat_index;
    int thread_num;
    int k_mer_length;
    int k_mer_min_freq;
@@ -26,6 +29,8 @@ typedef struct {
    int min_overlap_Len;
    int min_overlap_coverage;
    int max_short_tip;
+    int min_cnt;
+    int mid_cnt;

    float max_hang_rate;
    float min_drop_rate;
--- a/6
+++ b/6
@@ -3,7 +3,7 @@ CXXFLAGS=	-g -O3 -msse4.2 -mpopcnt -fomit-frame-pointer -Wall #-fsanitize=addres
 CPPFLAGS=
 INCLUDES=
 OBJS=		Output.o CommandLines.o Process_Read.o Assembly.o kmer.o Hash_Table.o \
-			POA.o Correct.o Levenshtein_distance.o Overlaps.o #ksw2_extz2_sse.o
+			POA.o Correct.o Levenshtein_distance.o Overlaps.o Trio.o kthread.o #ksw2_extz2_sse.o
 EXE=		hifiasm
 LIBS=		-lz -lpthread -lm #-fsanitize=address -fno-omit-frame-pointer

@@ -33,7 +33,7 @@ depend:

 Assembly.o: Assembly.h Process_Read.h kseq.h Overlaps.h kvec.h kdq.h
 Assembly.o: CommandLines.h kmer.h Hash_Table.h khash.h POA.h Correct.h
-Assembly.o: Levenshtein_distance.h Output.h
+Assembly.o: Levenshtein_distance.h Output.h Trio.h
 CommandLines.o: CommandLines.h ketopt.h
 Correct.o: Correct.h Hash_Table.h khash.h kmer.h Process_Read.h kseq.h
 Correct.o: Overlaps.h kvec.h kdq.h CommandLines.h Levenshtein_distance.h
@@ -51,4 +51,6 @@ Process_Read.o: Process_Read.h kseq.h Overlaps.h kvec.h kdq.h CommandLines.h
 kmer.o: kmer.h Process_Read.h kseq.h Overlaps.h kvec.h kdq.h CommandLines.h
 main.o: CommandLines.h Process_Read.h kseq.h Overlaps.h kvec.h kdq.h
 main.o: Assembly.h Levenshtein_distance.h
+Trio.o: Trio.h khashl.h kthread.h Process_Read.h CommandLines.h
+kthread.o: kthread.h
 #ksw2_extz2_sse.o: ksw2.h
--- a/Overlaps.cpp
+++ b/Overlaps.cpp
--- a/Overlaps.h
+++ b/Overlaps.h
@@ -1,5 +1,6 @@
 #ifndef __OVERLAPS__
 #define __OVERLAPS__
+#include <stdio.h>
 #include <stdint.h>
 #include "kvec.h"
 #include "kdq.h"
@@ -20,6 +21,7 @@
 #define PRIMARY_LABLE 0
 #define ALTER_LABLE 1
 #define HAP_LABLE 2
+#define TRIO_THRES 0.9
 // #define PRIMARY_LABLE 1
 // #define ALTER_LABLE 2
 // #define HAP_LABLE 4
@@ -40,6 +42,11 @@
 #define LONG_TIPS_UNDER_MAX_EXT 6
 #define LOOP 7

+#define TRIM 10
+#define CUT 11
+#define CUT_DIF_HAP 12
+
+



@@ -65,6 +72,7 @@ typedef struct {
 void init_ma_hit_t_alloc(ma_hit_t_alloc* x);
 void clear_ma_hit_t_alloc(ma_hit_t_alloc* x);
 void resize_ma_hit_t_alloc(ma_hit_t_alloc* x, uint64_t size);
+void destory_all_ma_hit_t_alloc(ma_hit_t_alloc* x, uint64_t n_read);
 void destory_ma_hit_t_alloc(ma_hit_t_alloc* x);
 void add_ma_hit_t_alloc(ma_hit_t_alloc* x, ma_hit_t* element);
 void ma_hit_sort_tn(ma_hit_t *a, long long n);
@@ -73,7 +81,6 @@ void ma_hit_sort_qns(ma_hit_t *a, long long n);
 int load_all_data_from_disk(ma_hit_t_alloc **sources, ma_hit_t_alloc **reverse_sources, 
 char* output_file_name);

-void normalize_ma_hit_t(ma_hit_t_alloc* sources, long long num_sources);


 typedef struct {
@@ -99,6 +106,17 @@ typedef struct {
 	uint8_t no_l_indel;
 } asg_arc_t;

+
+typedef struct {
+	uint32_t len:31, circ:1; // len: length of the unitig; circ: circular if non-zero
+	uint32_t start, end; // start: starting vertex in the string graph; end: ending vertex
+	uint32_t m, n; // number of reads
+	uint64_t *a; // list of reads
+	char *s; // unitig sequence is not null
+} ma_utg_t;
+
+
+
 typedef struct {
 	uint32_t len:31, del:1;
 	uint8_t c;
@@ -114,11 +132,26 @@ typedef struct {
 	uint64_t *idx;

 	uint8_t* seq_vis;
+
+	uint32_t n_F_seq;
+	ma_utg_t* F_seq;
 } asg_t;

 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 {
+	ma_utg_v u;
+	asg_t *g;
+} ma_ug_t;
+
+typedef struct {
+	uint32_t utg:31, ori:1, start, len;
+} utg_intv_t;
+
+
 #define MA_HT_INT        (-1)
 #define MA_HT_QCONT      (-2)
 #define MA_HT_TCONT      (-3)
@@ -209,6 +242,23 @@ static inline int ma_hit2arc(const ma_hit_t *h, int ql, int tl, int max_hang, fl
 #define asg_arc_n(g, v) ((uint32_t)(g)->idx[(v)])
 #define asg_arc_a(g, v) (&(g)->arc[(g)->idx[(v)]>>32])

+static inline uint32_t asg_get_arc(asg_t *g, uint32_t v, uint32_t w, asg_arc_t* t)
+{
+	uint32_t i, nv = asg_arc_n(g, v);
+	asg_arc_t *av = asg_arc_a(g, v);
+	for (i = 0; i < nv; ++i)
+	{
+		if(av[i].del) continue;
+		if(av[i].v == w)
+		{
+			(*t) = av[i];
+			return 1;
+		}
+	}
+
+	return 0;
+}
+
 // append an arc
 static inline asg_arc_t *asg_arc_pushp(asg_t *g)
 {
@@ -262,7 +312,8 @@ static inline void asg_seq_drop(asg_t *g, uint32_t s)
 				///if output node is at primary
 				/****************************may have hap bugs********************************/
 				///if(g->seq[(av[i].v>>1)].c == PRIMARY_LABLE)
-				if(g->seq[(av[i].v>>1)].c == PRIMARY_LABLE || g->seq[(av[i].v>>1)].c == HAP_LABLE)
+				///if(g->seq[(av[i].v>>1)].c == PRIMARY_LABLE || g->seq[(av[i].v>>1)].c == HAP_LABLE)
+				if(g->seq[(av[i].v>>1)].c != ALTER_LABLE)
 				{/****************************may have hap bugs********************************/
 					av[i].del = 1;
 					asg_arc_del(g, av[i].v^1, v^1, 1);
@@ -273,24 +324,7 @@ static inline void asg_seq_drop(asg_t *g, uint32_t s)
 }


-typedef struct {
-	uint32_t len:31, circ:1; // len: length of the unitig; circ: circular if non-zero
-	uint32_t start, end; // start: starting vertex in the string graph; end: ending vertex
-	uint32_t m, n; // number of reads
-	uint64_t *a; // list of reads
-	char *s; // unitig sequence is not null
-} ma_utg_t;

-typedef struct { size_t n, m; ma_utg_t *a; } ma_utg_v;
-
-typedef struct {
-	ma_utg_v u;
-	asg_t *g;
-} ma_ug_t;
-
-typedef struct {
-	uint32_t utg:31, ori:1, start, len;
-} utg_intv_t;


 /******************
@@ -300,7 +334,10 @@ typedef struct {
 typedef struct {
 	uint32_t p; // the optimal parent vertex
 	uint32_t d; // the shortest distance from the initial vertex
-	uint32_t c; // max count of reads
+	uint32_t c; // max count of positive reads
+	uint32_t m; // max count of negative reads
+	uint32_t np; // max count of non-positive reads
+	uint32_t nc; // max count of reads, no matter positive or negative
 	uint32_t r:31, s:1; // r: the number of remaining incoming arc; s: state
 	//s: state, s=0, this edge has not been visited, otherwise, s=1
 } binfo_t;
@@ -331,6 +368,12 @@ typedef struct {
 	uint64_t i;
 } kvec_t_u64_warp;

+typedef struct {
+	kvec_t(asg_arc_t) a;
+	uint64_t i;
+}kvec_asg_arc_t_warp;
+
+

 typedef struct {
 	uint32_t q_pos;
@@ -387,7 +430,6 @@ void add_overlaps(ma_hit_t_alloc* source_paf, ma_hit_t_alloc* dest_paf, uint64_t
 void remove_overlaps(ma_hit_t_alloc* source_paf, uint64_t* source_index, long long listLen);
 void add_overlaps_from_different_sources(ma_hit_t_alloc* source_paf_list, ma_hit_t_alloc* dest_paf, 
 uint64_t* source_index, long long listLen);
-void print_revise_edges(ma_hit_t_alloc* source_paf, uint64_t* source_index, long long listLen);

 #define EvaluateLen(U, id) ((U).a[(id)].start)
 #define IsMerge(U, id) ((U).a[(id)].end)
@@ -421,5 +463,574 @@ void destory_R_to_U(R_to_U* x);
 void set_R_to_U(R_to_U* x, uint32_t rID, uint32_t uID, uint32_t is_Unitig);
 void get_R_to_U(R_to_U* x, uint32_t rID, uint32_t* uID, uint32_t* is_Unitig);
 void transfor_R_to_U(R_to_U* x);
+void debug_utg_graph(ma_ug_t *ug, asg_t* read_g, int require_equal_nv, int test_tangle);
+void clean_untig_graph(ma_ug_t *ug, asg_t *read_g, ma_hit_t_alloc* reverse_sources,
+long long bubble_dist, long long tipsLen, float tip_drop_ratio, long long stops_threshold, 
+R_to_U* ruIndex, buf_t* b_0, uint8_t* visit, float density, uint32_t miniHapLen, 
+uint32_t miniBiGraph, float chimeric_rate, int is_final_clean);
+int asg_pop_bubble_primary(asg_t *g, int max_dist);
+long long asg_arc_del_simple_circle_untig(ma_hit_t_alloc* sources, ma_sub_t* coverage_cut, asg_t *g, long long circleLen, int is_drop);
+
+typedef struct {
+	asg_t* g;
+
+	asg_arc_t *av;
+	uint32_t nv;
+	uint32_t av_i;
+
+	asg_arc_t* new_edges;
+	uint32_t new_edges_n;
+	uint32_t new_edges_i;
+} Edge_iter;
+
+void init_Edge_iter(asg_t* g, uint32_t v, asg_arc_t* new_edges, uint32_t new_edges_n, Edge_iter* x);
+int get_arc_t(Edge_iter* x, asg_arc_t* get);
+int asg_pop_bubble_primary_trio(ma_ug_t *ug, int max_dist, uint32_t positive_flag, uint32_t negative_flag);
+
+
+inline int get_real_length(asg_t *g, uint32_t v, uint32_t* v_s)
+{
+    uint32_t i, kv = 0;
+    for (i = 0, kv = 0; i < asg_arc_n(g, v); i++)
+    {
+        if(!asg_arc_a(g, v)[i].del)
+        {
+            if(v_s) v_s[kv] = asg_arc_a(g, v)[i].v;
+            kv++;
+        }
+    }
+
+    return kv;
+}
+
+inline uint32_t check_tip(asg_t *sg, uint32_t begNode, uint32_t* endNode, buf_t* b, uint32_t max_ext)
+{
+    ///cut tip of length <= max_ext
+    uint32_t v = begNode, w;
+    uint32_t kv;
+    uint32_t eLen = 0;
+    (*endNode) = (uint32_t)-1;
+    b->b.n = 0;
+    while (1)
+    {
+        kv = get_real_length(sg, v, NULL);
+        (*endNode) = v;
+        eLen++;
+        if(b) kv_push(uint32_t, b->b, v);
+        if(kv == 0) return END_TIPS;
+        if(kv > 1) return MUL_OUTPUT;
+        ///if(eLen > max_ext) return LONG_TIPS;
+        ///kv must be 1 here
+        kv = get_real_length(sg, v, &w);
+        ///here this value must be >= 1
+        if(get_real_length(sg, w^1, NULL)!=1) return MUL_INPUT;
+        v = w;
+        if(v == begNode) return LOOP;
+        if(eLen >= max_ext) return LONG_TIPS;
+    }
+}
+
+inline uint32_t get_unitig_back(asg_t *sg, ma_ug_t *ug, uint32_t begNode, uint32_t* endNode, 
+long long* nodeLen, long long* baseLen, buf_t* b)
+{
+	ma_utg_v* u = NULL; 
+    uint32_t v = begNode, w, k;
+    uint32_t kv;
+	(*nodeLen) = (*baseLen) = 0;
+    (*endNode) = (uint32_t)-1;
+	if(ug!=NULL) u = &(ug->u);
+
+    while (1)
+    {
+        kv = get_real_length(sg, v, NULL);
+        (*endNode) = v;
+		if(u == NULL) 
+		{
+			(*nodeLen)++;
+		}
+		else
+		{
+			(*nodeLen) += EvaluateLen((*u), v>>1);
+		}
+        if(b) kv_push(uint32_t, b->b, v);
+		///means reach the end of a unitig
+		if(kv!=1) (*baseLen) += sg->seq[v>>1].len;
+		if(kv==0) return END_TIPS;
+        if(kv>1) return MUL_OUTPUT;
+        ///kv must be 1 here
+        kv = get_real_length(sg, v, &w);
+		///means reach the end of a unitig
+        if(get_real_length(sg, w^1, NULL)!=1)
+		{
+			(*baseLen) += sg->seq[v>>1].len;
+			return MUL_INPUT;
+		}
+
+		for (k = 0; k < asg_arc_n(sg, v); k++)
+        {
+			if(asg_arc_a(sg, v)[k].del) continue;
+			///here is just one undeleted edge
+			(*baseLen) += asg_arc_len(asg_arc_a(sg, v)[k]);
+			break;
+		}
+
+        v = w;
+        if(v == begNode) return LOOP;
+    }
+}
+
+inline uint32_t get_unitig(asg_t *sg, ma_ug_t *ug, uint32_t begNode, uint32_t* endNode, 
+long long* nodeLen, long long* baseLen, long long* max_stop_nodeLen, long long* max_stop_baseLen, 
+uint32_t stops_threshold, buf_t* b)
+{
+	ma_utg_v* u = NULL; 
+    uint32_t v = begNode, w, k;
+    uint32_t kv, return_flag, n_stops = 0;
+	long long pre_baseLen = 0, pre_nodeLen = 0;
+	long long cur_baseLen = 0, cur_nodeLen = 0;
+	(*max_stop_nodeLen) = (*max_stop_baseLen) = (*nodeLen) = (*baseLen) = 0;
+    (*endNode) = (uint32_t)-1;
+	if(ug!=NULL) u = &(ug->u);
+
+    while (1)
+    {
+        kv = get_real_length(sg, v, NULL);
+        (*endNode) = v;
+		if(u == NULL) 
+		{
+			(*nodeLen)++;
+		}
+		else
+		{
+			(*nodeLen) += EvaluateLen((*u), v>>1);
+		}
+        if(b) kv_push(uint32_t, b->b, v);
+		///means reach the end of a unitig
+		if(kv!=1) (*baseLen) += sg->seq[v>>1].len;
+		if(kv==0)
+		{
+			return_flag = END_TIPS;
+			break;
+			///return END_TIPS;
+		} 
+		if(kv>1)
+		{
+			return_flag = MUL_OUTPUT;
+			break;
+			///return MUL_OUTPUT;
+		}
+        ///kv must be 1 here
+        kv = get_real_length(sg, v, &w);
+		///means reach the end of a unitig
+        if(get_real_length(sg, w^1, NULL)!=1)
+		{
+
+			n_stops++;
+			if(n_stops >= stops_threshold)
+			{
+				(*baseLen) += sg->seq[v>>1].len;
+				return_flag = MUL_INPUT;
+				break;
+				///return MUL_INPUT;
+			}
+			else
+			{
+				for (k = 0; k < asg_arc_n(sg, v); k++)
+				{
+					if(asg_arc_a(sg, v)[k].del) continue;
+					///here is just one undeleted edge
+					(*baseLen) += asg_arc_len(asg_arc_a(sg, v)[k]);
+					break;
+				}
+			}
+
+			cur_baseLen = (*baseLen) - pre_baseLen;
+			pre_baseLen = (*baseLen);
+			if(cur_baseLen > (*max_stop_baseLen))
+			{
+				(*max_stop_baseLen) = cur_baseLen;
+			}
+
+
+			cur_nodeLen = (*nodeLen) - pre_nodeLen;
+			pre_nodeLen = (*nodeLen);
+			if(cur_nodeLen > (*max_stop_nodeLen))
+			{
+				(*max_stop_nodeLen) = cur_nodeLen;
+			}
+		}
+		else
+		{
+			for (k = 0; k < asg_arc_n(sg, v); k++)
+			{
+				if(asg_arc_a(sg, v)[k].del) continue;
+				///here is just one undeleted edge
+				(*baseLen) += asg_arc_len(asg_arc_a(sg, v)[k]);
+				break;
+			}
+		}
+		
+
+        v = w;
+        if(v == begNode)
+		{
+			return_flag = LOOP;
+			break;
+			///return LOOP;
+		} 
+    }
+
+
+
+
+	cur_baseLen = (*baseLen) - pre_baseLen;
+	pre_baseLen = (*baseLen);
+	if(cur_baseLen > (*max_stop_baseLen))
+	{
+		(*max_stop_baseLen) = cur_baseLen;
+	}
+
+
+	cur_nodeLen = (*nodeLen) - pre_nodeLen;
+	pre_nodeLen = (*nodeLen);
+	if(cur_nodeLen > (*max_stop_nodeLen))
+	{
+		(*max_stop_nodeLen) = cur_nodeLen;
+	}
+
+	return return_flag;
+}
+
+#define UNAVAILABLE (uint32_t)-1
+#define PLOID 0
+#define NON_PLOID 1
+#define DIFF_HAP_RATE 0.75
+#define TRIO_DROP_THRES 0.9
+#define TRIO_DROP_LENGTH_THRES 0.8
+#define MAX_STOP_RATE 0.6
+#define TANGLE_MISSED_THRES 0.6
+///if ug == NULL, nsg should be equal to read_sg
+inline uint32_t check_different_haps(asg_t *nsg, ma_ug_t *ug, asg_t *read_sg, 
+uint32_t v_0, uint32_t v_1, ma_hit_t_alloc* reverse_sources, buf_t* b_0, buf_t* b_1, 
+R_to_U* ruIndex, uint32_t min_edge_length, uint32_t stops_threshold)
+{
+	uint32_t vEnd, qn, tn, j, is_Unitig, uId;
+	long long ELen_0, ELen_1, tmp, max_stop_nodeLen, max_stop_baseLen;
+
+	b_0->b.n = b_1->b.n = 0;
+	if(get_unitig(nsg, ug, v_0, &vEnd, &ELen_0, &tmp, &max_stop_nodeLen, &max_stop_baseLen, 
+	stops_threshold, b_0) == LOOP)
+	{
+		return UNAVAILABLE;
+	} 
+	if(get_unitig(nsg, ug, v_1, &vEnd, &ELen_1, &tmp, &max_stop_nodeLen, &max_stop_baseLen, 
+	stops_threshold, b_1) == LOOP)
+	{
+		return UNAVAILABLE;
+	}
+
+	if(ELen_0<=min_edge_length || ELen_1<=min_edge_length) return UNAVAILABLE;
+
+	rIdContig b_max, b_min;
+	b_max.b_0 = b_min.b_0 = NULL;
+	b_max.offset = b_max.readI = b_max.untigI = 0;
+	b_min.offset = b_min.readI = b_min.untigI = 0;
+
+	if(ELen_0<=ELen_1)
+	{
+        b_min.b_0 = b_0;
+        b_max.b_0 = b_1;
+    }
+    else
+    {
+        b_min.b_0 = b_1;
+        b_max.b_0 = b_0;
+    }
+
+	uint32_t max_count = 0, min_count = 0;
+	ma_utg_t *node_min = NULL, *node_max = NULL;
+
+	if(ug != NULL)
+	{
+		/*****************************label all unitigs****************************************/
+		for (b_max.untigI = 0; b_max.untigI < b_max.b_0->b.n; b_max.untigI++)
+		{
+			node_max = &(ug->u.a[b_max.b_0->b.a[b_max.untigI]>>1]);
+			///each read
+			for (b_max.readI = 0; b_max.readI < node_max->n; b_max.readI++)
+			{
+				qn = (node_max->a[b_max.readI]>>33);
+				set_R_to_U(ruIndex, qn, (b_max.b_0->b.a[b_max.untigI]>>1), 1);
+			}
+		}
+		/*****************************label all unitigs****************************************/
+
+
+		///each unitig
+		for (b_min.untigI = 0; b_min.untigI < b_min.b_0->b.n; b_min.untigI++)
+		{
+			
+			node_min = &(ug->u.a[(b_min.b_0->b.a[b_min.untigI]>>1)]);
+
+			///each read
+			for (b_min.readI = 0; b_min.readI < node_min->n; b_min.readI++)
+        	{
+				qn = node_min->a[b_min.readI]>>33;
+
+				if(reverse_sources[qn].length>=0) min_count++;
+				for (j = 0; j < (long long)reverse_sources[qn].length; j++)
+				{
+					tn = Get_tn(reverse_sources[qn].buffer[j]);
+					if(read_sg->seq[tn].del == 1)
+					{
+						get_R_to_U(ruIndex, tn, &tn, &is_Unitig);
+						if(tn == (uint32_t)-1 || is_Unitig == 1 || read_sg->seq[tn].del == 1) continue;
+					}
+
+					get_R_to_U(ruIndex, tn, &uId, &is_Unitig);
+					if(uId!=(uint32_t)-1 && is_Unitig == 1)
+					{
+						max_count++;
+						break;
+					}
+				}
+			}
+		}
+
+		/*****************************label all unitigs****************************************/
+		for (b_max.untigI = 0; b_max.untigI < b_max.b_0->b.n; b_max.untigI++)
+		{
+			node_max = &(ug->u.a[b_max.b_0->b.a[b_max.untigI]>>1]);
+			///each read
+			for (b_max.readI = 0; b_max.readI < node_max->n; b_max.readI++)
+			{
+				qn = (node_max->a[b_max.readI]>>33);
+				ruIndex->index[qn] = (uint32_t)-1;
+			}
+		}
+		/*****************************label all unitigs****************************************/
+		
+	}
+	else
+	{
+		/*****************************label all reads****************************************/
+		for (b_max.untigI = 0; b_max.untigI < b_max.b_0->b.n; b_max.untigI++)
+		{
+			qn = (b_max.b_0->b.a[b_max.untigI]>>1);
+			set_R_to_U(ruIndex, qn, 1, 1);
+		}
+		/*****************************label all reads****************************************/
+
+		///each read
+		for (b_min.untigI = 0; b_min.untigI < b_min.b_0->b.n; b_min.untigI++)
+		{
+			qn = (b_min.b_0->b.a[b_min.untigI]>>1);
+
+			if(reverse_sources[qn].length>=0) min_count++;
+
+			for (j = 0; j < (long long)reverse_sources[qn].length; j++)
+        	{
+            	tn = Get_tn(reverse_sources[qn].buffer[j]);
+				if(nsg->seq[tn].del == 1)
+				{
+					get_R_to_U(ruIndex, tn, &tn, &is_Unitig);
+					if(tn == (uint32_t)-1 || is_Unitig == 1 || nsg->seq[tn].del == 1) continue;
+				}
+
+
+				get_R_to_U(ruIndex, tn, &uId, &is_Unitig);
+				if(uId!=(uint32_t)-1 && is_Unitig == 1)
+				{
+					max_count++;
+					break;
+				}
+			}
+		}
+
+		/*****************************label all reads****************************************/
+		for (b_max.untigI = 0; b_max.untigI < b_max.b_0->b.n; b_max.untigI++)
+		{
+			qn = (b_max.b_0->b.a[b_max.untigI]>>1);
+			ruIndex->index[qn] = (uint32_t)-1;
+		}
+		/*****************************label all reads****************************************/
+	}
+
+	// if(((v_0==7707) && (v_1==26867))||((v_1==7707) && (v_0==26867)))
+	// {
+	// 	fprintf(stderr, "******\nv_0>>1: %u, v_0&1: %u, ELen_0: %u\n", v_0>>1, v_0&1, (uint32_t)ELen_0);
+	// 	fprintf(stderr, "v_1>>1: %u, v_1&1: %u, ELen_1: %u\n", v_1>>1, v_1&1, (uint32_t)ELen_1);
+	// 	fprintf(stderr, "min_count: %u, max_count: %u, DIFF_HAP_RATE: %f\n\n", 
+	// 	min_count, max_count, DIFF_HAP_RATE);
+	// }
+	
+	if(min_count == 0) return UNAVAILABLE;
+	if(max_count > min_count*DIFF_HAP_RATE) return PLOID;
+	return NON_PLOID;
+}
+
+
+inline uint32_t check_different_haps_naive(asg_t *nsg, ma_ug_t *ug, asg_t *read_sg, 
+uint32_t v_0, uint32_t v_1, ma_hit_t_alloc* reverse_sources, buf_t* b_0, buf_t* b_1, 
+R_to_U* ruIndex, uint32_t min_edge_length, uint32_t stops_threshold)
+{
+	uint32_t vEnd, qn, tn, j, is_Unitig;
+	long long ELen_0, ELen_1, tmp, max_stop_nodeLen, max_stop_baseLen;
+
+	b_0->b.n = b_1->b.n = 0;
+	if(get_unitig(nsg, ug, v_0, &vEnd, &ELen_0, &tmp, &max_stop_nodeLen, &max_stop_baseLen, 
+	stops_threshold, b_0) == LOOP)
+	{
+		return UNAVAILABLE;
+	} 
+	if(get_unitig(nsg, ug, v_1, &vEnd, &ELen_1, &tmp, &max_stop_nodeLen, &max_stop_baseLen, 
+	stops_threshold, b_1) == LOOP)
+	{
+		return UNAVAILABLE;
+	}
+
+	if(ELen_0<=min_edge_length || ELen_1<=min_edge_length) return UNAVAILABLE;
+
+	rIdContig b_max, b_min;
+	b_max.b_0 = b_min.b_0 = NULL;
+	b_max.offset = b_max.readI = b_max.untigI = 0;
+	b_min.offset = b_min.readI = b_min.untigI = 0;
+
+	if(ELen_0<=ELen_1)
+	{
+        b_min.b_0 = b_0;
+        b_max.b_0 = b_1;
+    }
+    else
+    {
+        b_min.b_0 = b_1;
+        b_max.b_0 = b_0;
+    }
+
+	uint32_t max_count = 0, min_count = 0;
+	ma_utg_t *node_min = NULL, *node_max = NULL;
+
+	if(ug != NULL)
+	{
+		///each unitig
+		for (b_min.untigI = 0; b_min.untigI < b_min.b_0->b.n; b_min.untigI++)
+		{
+			
+			node_min = &(ug->u.a[(b_min.b_0->b.a[b_min.untigI]>>1)]);
+
+			///each read
+			for (b_min.readI = 0; b_min.readI < node_min->n; b_min.readI++)
+        	{
+				qn = node_min->a[b_min.readI]>>33;
+
+				if(reverse_sources[qn].length>=0) min_count++;
+				for (j = 0; j < (long long)reverse_sources[qn].length; j++)
+				{
+					tn = Get_tn(reverse_sources[qn].buffer[j]);
+					if(read_sg->seq[tn].del == 1)
+					{
+						get_R_to_U(ruIndex, tn, &tn, &is_Unitig);
+						if(tn == (uint32_t)-1 || is_Unitig == 1 || read_sg->seq[tn].del == 1) continue;
+					}
+
+					///each unitig
+					for (b_max.untigI = 0; b_max.untigI < b_max.b_0->b.n; b_max.untigI++)
+					{
+						node_max = &(ug->u.a[b_max.b_0->b.a[b_max.untigI]>>1]);
+						///each read
+						for (b_max.readI = 0; b_max.readI < node_max->n; b_max.readI++)
+						{
+							if(tn == (node_max->a[b_max.readI]>>33))
+							{
+								max_count++;
+								goto end_check_different_haps_ug;
+							}
+						}
+					}
+				}
+
+				end_check_different_haps_ug:;
+			}
+		}
+	}
+	else
+	{
+		///each read
+		for (b_min.untigI = 0; b_min.untigI < b_min.b_0->b.n; b_min.untigI++)
+		{
+			qn = (b_min.b_0->b.a[b_min.untigI]>>1);
+
+			if(reverse_sources[qn].length>=0) min_count++;
+
+			for (j = 0; j < (long long)reverse_sources[qn].length; j++)
+        	{
+            	tn = Get_tn(reverse_sources[qn].buffer[j]);
+				if(nsg->seq[tn].del == 1)
+				{
+					get_R_to_U(ruIndex, tn, &tn, &is_Unitig);
+					if(tn == (uint32_t)-1 || is_Unitig == 1 || nsg->seq[tn].del == 1) continue;
+				}
+
+				///each read
+				for (b_max.untigI = 0; b_max.untigI < b_max.b_0->b.n; b_max.untigI++)
+                {
+					if((b_max.b_0->b.a[b_max.untigI]>>1) == tn)
+					{
+						max_count++;
+						goto end_check_different_haps_non_ug;
+					}
+				}
+			}
+
+			end_check_different_haps_non_ug:;
+		}
+	}
+
+	// if(((v_0==7707) && (v_1==26867))||((v_1==7707) && (v_0==26867)))
+	// {
+	// 	fprintf(stderr, "******\nv_0>>1: %u, v_0&1: %u, ELen_0: %u\n", v_0>>1, v_0&1, (uint32_t)ELen_0);
+	// 	fprintf(stderr, "v_1>>1: %u, v_1&1: %u, ELen_1: %u\n", v_1>>1, v_1&1, (uint32_t)ELen_1);
+	// 	fprintf(stderr, "min_count: %u, max_count: %u, DIFF_HAP_RATE: %f\n\n", 
+	// 	min_count, max_count, DIFF_HAP_RATE);
+	// }
+	
+	if(min_count == 0) return UNAVAILABLE;
+	if(max_count > min_count*DIFF_HAP_RATE) return PLOID;
+	return NON_PLOID;
+}
+
+
+
+typedef struct {
+	uint32_t father_occ;
+	uint32_t mother_occ;
+	uint32_t ambig_occ;
+	uint32_t drop_occ;
+	uint32_t total;
+} Trio_counter;
+
+void resolve_tangles(ma_ug_t *src, asg_t *read_g, ma_hit_t_alloc* reverse_sources, long long minLongUntig, 
+long long maxShortUntig, float l_untig_rate, float max_node_threshold, R_to_U* ruIndex, uint32_t trio_flag, 
+float drop_ratio);
+void adjust_utg_advance(asg_t *sg, ma_ug_t *ug, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex);
+void rescue_contained_reads_aggressive(ma_ug_t *i_ug, asg_t *r_g, ma_hit_t_alloc* sources, ma_sub_t *coverage_cut,
+R_to_U* ruIndex, int max_hang, int min_ovlp, long long bubble_dist, uint32_t chainLenThres, uint32_t is_bubble_check, 
+uint32_t is_primary_check, kvec_asg_arc_t_warp* new_rtg_edges, kvec_t_u32_warp* new_rtg_nodes);
+void rescue_missing_overlaps_aggressive(ma_ug_t *i_ug, asg_t *r_g, ma_hit_t_alloc* sources, ma_sub_t *coverage_cut,
+R_to_U* ruIndex, int max_hang, int min_ovlp, long long bubble_dist, uint32_t is_bubble_check, 
+uint32_t is_primary_check, kvec_asg_arc_t_warp* new_rtg_edges);
+void deduplicate(ma_ug_t *src, asg_t *read_g, ma_hit_t_alloc* reverse_sources, long long minLongUntig, 
+long long maxShortUntig, float l_untig_rate, float max_node_threshold, R_to_U* ruIndex);
+void all_to_all_deduplicate(ma_ug_t* ug, uint8_t postive_flag, float drop_rate, 
+ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex);
+void drop_semi_circle(ma_ug_t *ug, asg_t* nsg, asg_t* read_g, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex);
+void rescue_wrong_overlaps_to_unitigs(ma_ug_t *i_ug, asg_t *r_g,  ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, 
+ma_sub_t *coverage_cut, R_to_U* ruIndex, int max_hang, int min_ovlp, long long bubble_dist, kvec_asg_arc_t_warp* keep_edges);
+void get_unitig_trio_flag(ma_utg_t* nsu, uint32_t flag, uint32_t* require, uint32_t* non_require, uint32_t* ambigious);
+void rescue_missing_overlaps_backward(ma_ug_t *i_ug, asg_t *r_g, ma_hit_t_alloc* sources, ma_sub_t *coverage_cut,
+R_to_U* ruIndex, int max_hang, int min_ovlp, long long bubble_dist, uint32_t backward_steps, 
+uint32_t is_bubble_check, uint32_t is_primary_check);
+uint32_t get_edge_from_source(ma_hit_t_alloc* sources, ma_sub_t *coverage_cut, 
+R_to_U* ruIndex, int max_hang, int min_ovlp, uint32_t query, uint32_t target, asg_arc_t* t);

 #endif
--- a/Process_Read.cpp
+++ b/Process_Read.cpp
@@ -59,6 +59,7 @@ void init_All_reads(All_reads* r)
 	r->name = NULL;
 	r->total_name_length = 0;
 	r->total_reads = 0;
+	r->trio_flag = NULL;
 }

 void destory_All_reads(All_reads* r)
@@ -77,6 +78,7 @@ void destory_All_reads(All_reads* r)
 	free(r->name);
 	free(r->name_index);
 	free(r->read_length);
+	free(r->trio_flag);
 }


@@ -123,6 +125,7 @@ void write_All_reads(All_reads* r, char* read_file_name)
 	
 	fwrite(r->name, sizeof(char), r->total_name_length, fp);
 	fwrite(r->name_index, sizeof(uint64_t), r->name_index_size, fp);
+	fwrite(r->trio_flag, sizeof(uint8_t), r->total_reads, fp);
    free(index_name);    
 	fflush(fp);
    fclose(fp);
@@ -201,6 +204,10 @@ int load_All_reads(All_reads* r, char* read_file_name)
 	r->name_index = (uint64_t*)malloc(sizeof(uint64_t)*r->name_index_size);
 	f_flag += fread(r->name_index, sizeof(uint64_t), r->name_index_size, fp);

+	/****************************may have bugs********************************/
+	r->trio_flag = (uint8_t*)malloc(sizeof(uint8_t)*r->total_reads);
+	f_flag += fread(r->trio_flag, sizeof(uint8_t), r->total_reads, fp);
+	/****************************may have bugs********************************/

 	r->cigars = (Compressed_Cigar_record*)malloc(sizeof(Compressed_Cigar_record)*r->total_reads);
 	r->second_round_cigar = (Compressed_Cigar_record*)malloc(sizeof(Compressed_Cigar_record)*r->total_reads);
@@ -279,7 +286,8 @@ void malloc_All_reads(All_reads* r)

 	r->name = (char*)malloc(sizeof(char)*r->total_name_length);
 	r->N_site = (uint64_t**)calloc(r->total_reads, sizeof(uint64_t*));
-
+	r->trio_flag = (uint8_t*)malloc(r->total_reads*sizeof(uint8_t));
+	memset(r->trio_flag, AMBIGU, r->total_reads*sizeof(uint8_t));
 }

 void destory_UC_Read(UC_Read* r)
@@ -1065,3 +1073,59 @@ void reverse_complement(char* pattern, uint64_t length)
 }


+
+typedef struct {
+  char* tmp;
+  long long tmpSize;
+  char* dest;
+  long long destSize;
+  FILE* fp;
+} LineReader;
+
+int get_single_line(LineReader* line)
+{
+    long long currentLen = 0, getLen = 0;
+    line->tmp[0] = '\0';
+    while (fgets(line->tmp, line->tmpSize, line->fp) != NULL)
+    {
+        getLen = strlen(line->tmp);
+        if(getLen + currentLen >= line->destSize)
+        {
+            line->destSize = getLen + currentLen + 1;
+            line->dest = (char*)realloc(line->dest, line->destSize);
+        }
+        memcpy(line->dest + currentLen, line->tmp, getLen+1);
+        currentLen = currentLen + getLen;
+        if(currentLen > 0 && line->dest[currentLen - 1] == '\n')
+        {
+            return 1;
+        }
+    }
+        
+    if(currentLen > 0)
+    {
+        return 1;
+    }
+    else
+    {
+        return 0;
+    }
+}
+
+void get_trio_info(char* input, uint8_t* pm)
+{
+    uint32_t i;
+    (*pm) = AMBIGU;
+    for (i = 0; input[i] != '\0'; i++)
+    {
+        if(input[i] == '\t')
+        {
+            break;
+        }
+    }
+
+    input[i] = '\0';
+    i++;
+    if(input[i] == 'p') (*pm) = FATHER;
+    if(input[i] == 'm') (*pm) = MOTHER;
+}
--- a/Process_Read.h
+++ b/Process_Read.h
@@ -106,7 +106,12 @@ typedef struct
 	uint32_t new_length;
 }Compressed_Cigar_record;

-
+#define AMBIGU 0
+#define FATHER 1
+#define MOTHER 2
+#define MIX_TRIO 3
+#define NON_TRIO 4
+#define DROP 5
 typedef struct
 {
 	uint64_t** N_site;
@@ -117,6 +122,7 @@ typedef struct
 	uint8_t** read_sperate;
 	uint64_t* read_length;
 	uint64_t* read_size;
+	uint8_t* trio_flag;

 	///seq start pos in uint8_t* read
 	///do not need it
--- a/README.md
+++ b/README.md
@@ -13,7 +13,8 @@ cd hifiasm && make
 Hifiasm is a fast haplotype-resolved de novo assembler for PacBio
 Hifi reads. Unlike most existing assemblers, hifiasm starts from uncollapsed
 genome. Thus, it is able to keep the haplotype information as much as possible.
-The input of hifiasm is the PacBio Hifi reads in fasta/fastq format, and its
+
+For non-trio assembly, the input of hifiasm is the PacBio Hifi reads in fasta/fastq format, and its
 outputs consist of: 

 1. Haplotype-resolved raw [unitig][unitig] graph in [GFA][gfa] format
@@ -30,9 +31,18 @@ outputs consist of:
 5. Haplotype-aware error corrected reads in fasta format (*prefix*.ec.fa).
 6. All-to-all overlaps in the [PAF][paf] format (*prefix*.ovlp.paf).

-So far hifiasm is still in early development stage, it will output phased
-chromosome-level high-quality assembly in the near future. In addition, hifiasm
-also outputs three binary files that save all overlap inforamtion
+For trio assembly, the input of hifiasm is the PacBio Hifi reads in fasta/fastq format, and the paternal/maternal trio indexes generated by `yak count` (see https://github.com/lh3/yak). The outputs consist of:
+1. Haplotype-resolved raw [unitig][unitig] graph in [GFA][gfa] format
+   (*prefix*.r\_utg.gfa). This graph keeps all haplotype information. 
+
+2. Phased maternal unitig graph (*prefix*.m.r\_utg.gfa).
+This graph keeps the phased maternal assembly.
+
+3. Phased paternal unitig graph (*prefix*.p.r\_utg.gfa).
+This graph keeps the phased paternal assembly.
+
+
+In addition, hifiasm also outputs three binary files that save all overlap inforamtion
 (hifiasm.asm.ovlp, hifiasm.asm.ovlp.source, hifiasm.asm.ovlp.reverse in default). With these files, hifiasm can avoid the time-consuming all-to-all overlap calculation step, and do the assembly
 directly and quickly. This might be helpful when you want to get an optimized
 assembly by multiple rounds of experiments with different parameters.
@@ -44,9 +54,10 @@ assembly in a few hours. Hifiasm has been tested on the following datasets:
 |<sub>Dataset<sub>|<sub>GSize<sub>|<sub>Cov<sub>|<sub>Asm options<sub>|<sub>CPU time<sub>|<sub>Wall time<sub>|<sub>RAM<sub>|<sub>[unitig][unitig]/[contig][unitig] N50<sup>[1]</sup><sub>|
 |:---------------|-----:|-----:|:---------------------|-------:|--------:|----:|----------------:|
 |<sub>[Human NA12878]<sub>|<sub>3Gb<sub>|<sub>x28<sub>|<sub>-k 40 -t 42 -r 2<sub>|<sub>200h<sub>|    <sub>5h32m<sub>|<sub>114G<sub>|<sub>93.5Kb/28.2Mb<sub>|
-|<sub>[Human HG002]<sub>|<sub>3Gb<sub>|<sub>x43<sub>|<sub>-k 40 -t 42 -r 2<sub>|<sub>405h10m<sub>|<sub>12h7m<sub>|<sub>146G<sub>|<sub>320kb/35.3Mb<sub>|
+|<sub>[Human HG002]<sub>|<sub>3Gb<sub>|<sub>x43<sub>|<sub>-k 40 -t 42 -r 2<sub>|<sub>405h10m<sub>|<sub>12h7m<sub>|<sub>146G<sub>|<sub>320kb/46.0Mb<sub>|
 |<sub>[Human CHM13]<sub>|<sub>3Gb<sub>|<sub>x27<sub>|<sub>-k 40 -t 42 -r 2<sub>|<sub>157h28m<sub>|<sub>5h10m<sub>|<sub>85.8G<sub>|<sub>NA<sup>[2]</sup>/41.4Mb<sub>|
 |<sub>[Butterfly]<sub>|<sub>358Mb<sub>|<sub>x35<sub>|<sub>-k 40 -t 42 -r 2 -z 20<sub>|<sub>17h6m<sub>|<sub>36m<sub>|<sub>16G<sub>|<sub>7.5Mb/NA<sup>[3]</sup><sub>|
+|<sub>[\[Redwood\]](https://downloads.pacbcloud.com/public/dataset/redwood2020/)<sub>|<sub>26.5Gb<sub>|<sub>x23<sub>|<sub>-k 40 -t 64 -r 2<sub>|<sub>7274h30m<sub>|<sub>141h30m<sub>|<sub>512G<sub>|<sub>1.7Mb/1.9Mb<sub>|

 <sub>[1] unitig N50 is the N50 of assembly graph with haplotype information (i.e., bubbles), while the contig N50 is the N50 of haplotype collapsed assembly (i.e., without bubbles).
 [2] CHM13 is a homozygous sample, so that unitig N50 makes no sense.
@@ -85,6 +96,22 @@ the assembly quality, adapters should be removed by `-z` as follow:

 In this example, hifiasm will remove 20 bases from both ends of each read.

+For trio assembly, first the trio indexes of paternal/maternal should be generated by 
+`yak count` (see https://github.com/lh3/yak):
+
+```sh
+./yak count -k31 -b37 -t16 -o mat.yak mat.fq.gz
+```
+```sh
+./yak count -k31 -b37 -t16 -o pat.yak pat.fq.gz
+```
+
+and then run hifiasm as follow:
+
+```sh
+./hifiasm -o NA12878.asm -t 32 -P pat.yak -M mat.yak NA12878_1.fq.gz NA12878_2.fq.gz
+```
+
 [unitig]: http://wgs-assembler.sourceforge.net/wiki/index.php/Celera_Assembler_Terminology
 [gfa]: https://github.com/pmelsted/GFA-spec/blob/master/GFA-spec.md
 [paf]: https://github.com/lh3/miniasm/blob/master/PAF.md
--- a/Trio.cpp
+++ b/Trio.cpp
@@ -0,0 +1,364 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <stdarg.h>
+#include <string.h>
+#include <assert.h>
+#include "khashl.h" // hash table
+#include "kthread.h"
+#include "Process_Read.h"
+#include "Trio.h"
+#include "CommandLines.h"
+
+#define CALLOC(ptr, len) ((ptr) = (__typeof__(ptr))calloc((len), sizeof(*(ptr))))
+#define MALLOC(ptr, len) ((ptr) = (__typeof__(ptr))malloc((len) * sizeof(*(ptr))))
+#define REALLOC(ptr, len) ((ptr) = (__typeof__(ptr))realloc((ptr), (len) * sizeof(*(ptr))))
+#define yak_ch_eq(a, b) ((a)>>YAK_COUNTER_BITS == (b)>>YAK_COUNTER_BITS) // lower 8 bits for counts; higher bits for k-mer
+#define yak_ch_hash(a) ((a)>>YAK_COUNTER_BITS)
+KHASHL_SET_INIT(, yak_ht_t, yak_ht, uint64_t, yak_ch_hash, yak_ch_eq)
+
+///#define CHUNK_SIZE 200000
+
+unsigned char seq_nt4_table[256] = { // translate ACGT to 0123
+	0, 1, 2, 3,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
+	4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
+	4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
+	4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
+	4, 0, 4, 1,  4, 4, 4, 2,  4, 4, 4, 4,  4, 4, 4, 4,
+	4, 4, 4, 4,  3, 3, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
+	4, 0, 4, 1,  4, 4, 4, 2,  4, 4, 4, 4,  4, 4, 4, 4,
+	4, 4, 4, 4,  3, 3, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
+	4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
+	4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
+	4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
+	4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
+	4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
+	4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
+	4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
+	4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4
+};
+
+static inline uint64_t yak_hash64(uint64_t key, uint64_t mask) // invertible integer hash function
+{
+	key = (~key + (key << 21)) & mask; // key = (key << 21) - key - 1;
+	key = key ^ key >> 24;
+	key = ((key + (key << 3)) + (key << 8)) & mask; // key * 265
+	key = key ^ key >> 14;
+	key = ((key + (key << 2)) + (key << 4)) & mask; // key * 21
+	key = key ^ key >> 28;
+	key = (key + (key << 31)) & mask;
+	return key;
+}
+
+static inline uint64_t yak_hash64_64(uint64_t key)
+{
+	key = ~key + (key << 21);
+	key = key ^ key >> 24;
+	key = (key + (key << 3)) + (key << 8);
+	key = key ^ key >> 14;
+	key = (key + (key << 2)) + (key << 4);
+	key = key ^ key >> 28;
+	key = key + (key << 31);
+	return key;
+}
+
+static inline uint64_t yak_hash_long(uint64_t x[4])
+{
+	int j = x[1] < x[3]? 0 : 1;
+	return yak_hash64_64(x[j<<1|0]) + yak_hash64_64(x[j<<1|1]);
+}
+
+int yak_ch_get(const yak_ch_t *h, uint64_t x)
+{
+	int mask = (1<<h->pre) - 1;
+	yak_ht_t *g = h->h[x&mask].h;
+	khint_t k;
+	k = yak_ht_get(g, x >> h->pre << YAK_COUNTER_BITS);
+	return k == kh_end(g)? -1 : kh_key(g, k)&YAK_MAX_COUNT;
+}
+
+yak_bf_t *yak_bf_init(int n_shift, int n_hashes)
+{
+	yak_bf_t *b;
+	void *ptr = 0;
+	if (n_shift + YAK_BLK_SHIFT > 64 || n_shift < YAK_BLK_SHIFT) return 0;
+	b = (yak_bf_t*)calloc(1, sizeof(yak_bf_t));
+	b->n_shift = n_shift;
+	b->n_hashes = n_hashes;
+	posix_memalign(&ptr, 1<<(YAK_BLK_SHIFT-3), 1ULL<<(n_shift-3));
+	b->b = (uint8_t*)ptr;
+	bzero(b->b, 1ULL<<(n_shift-3));
+	return b;
+}
+
+yak_ch_t *yak_ch_init(int k, int pre, int n_hash, int n_shift)
+{
+	yak_ch_t *h;
+	int i;
+	if (pre < YAK_COUNTER_BITS) return 0;
+	CALLOC(h, 1);
+	h->k = k, h->pre = pre;
+	CALLOC(h->h, 1<<h->pre);
+	for (i = 0; i < 1<<h->pre; ++i)
+		h->h[i].h = yak_ht_init();
+	if (n_hash > 0 && n_shift > h->pre) {
+		h->n_hash = n_hash, h->n_shift = n_shift;
+		for (i = 0; i < 1<<h->pre; ++i)
+			h->h[i].b = yak_bf_init(h->n_shift - h->pre, h->n_hash);
+	}
+	return h;
+}
+
+void yak_bf_destroy(yak_bf_t *b)
+{
+	if (b == 0) return;
+	free(b->b); free(b);
+}
+
+void yak_ch_destroy_bf(yak_ch_t *h)
+{
+	int i;
+	for (i = 0; i < 1<<h->pre; ++i) {
+		if (h->h[i].b)
+			yak_bf_destroy(h->h[i].b);
+		h->h[i].b = 0;
+	}
+}
+
+
+yak_ch_t *yak_ch_restore_core(yak_ch_t *ch0, const char *fn, int mode, ...)
+{
+	va_list ap;
+	FILE *fp;
+	uint32_t t[3];
+	char magic[4];
+	int i, j, absent, min_cnt = 0, mid_cnt = 0, mode_err = 0;
+	uint64_t mask = (1ULL<<YAK_COUNTER_BITS) - 1, n_ins = 0, n_new = 0;
+	yak_ch_t *ch;
+
+	va_start(ap, mode);
+	if (mode == YAK_LOAD_ALL) { // do nothing
+	} else if (mode == YAK_LOAD_TRIOBIN1 || mode == YAK_LOAD_TRIOBIN2) {
+		assert(YAK_COUNTER_BITS >= 4);
+		min_cnt = va_arg(ap, int);
+		mid_cnt = va_arg(ap, int);
+		if (ch0 == 0 && mode == YAK_LOAD_TRIOBIN2)
+			mode_err = 1;
+	} else mode_err = 1;
+	va_end(ap);
+	if (mode_err) return 0;
+
+	if ((fp = fopen(fn, "rb")) == 0) return 0;
+	if (fread(magic, 1, 4, fp) != 4) return 0;
+	if (strncmp(magic, YAK_MAGIC, 4) != 0) {
+		fprintf(stderr, "ERROR: wrong file magic.\n");
+		fclose(fp);
+		return 0;
+	}
+	fread(t, 4, 3, fp);
+	if (t[2] != YAK_COUNTER_BITS) {
+		fprintf(stderr, "ERROR: saved counter bits: %d; compile-time counter bits: %d\n", t[2], YAK_COUNTER_BITS);
+		fclose(fp);
+		return 0;
+	}
+
+	ch = ch0 == 0? yak_ch_init(t[0], t[1], 0, 0) : ch0;
+	assert((int)t[0] == ch->k && (int)t[1] == ch->pre);
+	for (i = 0; i < 1<<ch->pre; ++i) {
+		yak_ht_t *h = ch->h[i].h;
+		fread(t, 4, 2, fp);
+		if (ch0 == 0) yak_ht_resize(h, t[0]);
+		for (j = 0; j < t[1]; ++j) {
+			uint64_t key;
+			fread(&key, 8, 1, fp);
+			if (mode == YAK_LOAD_ALL) {
+				++n_ins;
+				yak_ht_put(h, key, &absent);
+				if (absent) ++n_new;
+			} else if (mode == YAK_LOAD_TRIOBIN1 || mode == YAK_LOAD_TRIOBIN2) {
+				int cnt = key & mask, x, shift = mode == YAK_LOAD_TRIOBIN1? 0 : 2;
+				if (cnt >= mid_cnt) x = 2<<shift;
+				else if (cnt >= min_cnt) x = 1<<shift;
+				else x = -1;
+				if (x >= 0) {
+					khint_t k;
+					key = (key & ~mask) | x;
+					++n_ins;
+					k = yak_ht_put(h, key, &absent);
+					if (absent) ++n_new;
+					else kh_key(h, k) = kh_key(h, k) | x;
+				}
+			}
+		}
+	}
+	fclose(fp);
+	///fprintf(stderr, "[M::%s] inserted %ld k-mers, of which %ld are new\n", __func__, (long)n_ins, (long)n_new);
+	return ch;
+}
+
+
+void yak_ch_destroy(yak_ch_t *h)
+{
+	int i;
+	if (h == 0) return;
+	yak_ch_destroy_bf(h);
+	for (i = 0; i < 1<<h->pre; ++i)
+		yak_ht_destroy(h->h[i].h);
+	free(h->h); free(h);
+}
+
+
+static char tb_classify(const int sc[2], const int *c, int k, double ratio_thres)
+{
+	char type;
+	if (sc[0] == 0 && sc[1] == 0) {
+		if (c[0<<2|2] == c[2<<2|0]) type = '0';
+		else if (c[0<<2|2] >= k - 4 + c[2<<2|0] && (c[2<<2|0] <= 1 || c[0<<2|2] * 0.05 > c[2<<2|0])) type = 'p';
+		else if (c[2<<2|0] >= k - 4 + c[0<<2|2] && (c[0<<2|2] <= 1 || c[2<<2|0] * 0.05 > c[0<<2|2])) type = 'm';
+		else type = '0';
+	} else if (sc[0] > k && sc[1] > k) {
+		type = 'a';
+	} else if (sc[0] >= k - 4 + sc[1] && sc[0] * 0.05 >= sc[1] && c[0<<2|2] * ratio_thres > c[2<<2|0]) {
+		type = 'p';
+	} else if (sc[1] >= k - 4 + sc[0] && sc[1] * 0.05 >= sc[0] && c[2<<2|0] * ratio_thres > c[0<<2|2]) {
+		type = 'm';
+	} else {
+		type = 'a';
+	}
+	return type;
+}
+
+static void tb_worker(void *_data, long k, int tid)
+{
+	///tb_step_t *t = (tb_step_t*)_data;
+	///tb_shared_t *aux = t->aux;
+	tb_shared_t *aux = (tb_shared_t*)_data;
+	UC_Read *s = &aux->bseq[tid]; 
+	recover_UC_Read(s, aux->seq, k);
+	tb_buf_t *b = &aux->buf[tid];
+	tb_cnt_t cnt; memset(&cnt, 0, sizeof(tb_cnt_t));
+	uint64_t x[4], mask;
+	int i, l, shift;
+	if (aux->ch->k < 32) {
+		mask = (1ULL<<2*aux->ch->k) - 1;
+		shift = 2 * (aux->ch->k - 1);
+	} else {
+		mask = (1ULL<<aux->ch->k) - 1;
+		shift = aux->ch->k - 1;
+	}
+	if (s->length > b->max) {
+		b->max = s->length;
+		kroundup32(b->max);
+		b->s = (uint32_t*)realloc(b->s, b->max * sizeof(uint32_t));
+	}
+	memset(b->s, 0, s->length * sizeof(uint32_t));
+	for (i = l = 0, x[0] = x[1] = x[2] = x[3] = 0; i < s->length; ++i) {
+		int flag, c = seq_nt4_table[(uint8_t)s->seq[i]];
+		if (c < 4) {
+			if (aux->ch->k < 32) {
+				x[0] = (x[0] << 2 | c) & mask;
+				x[1] = x[1] >> 2 | (uint64_t)(3 - c) << shift;
+			} else {
+				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 >= aux->k) {
+				int type = 0, c1, c2;
+				uint64_t y;
+
+
+				//++t->cnt[k].nk;
+				++cnt.nk;
+
+
+				if (aux->ch->k < 32)
+					y = yak_hash64(x[0] < x[1]? x[0] : x[1], mask);
+				else
+					y = yak_hash_long(x);
+				flag = yak_ch_get(aux->ch, y);
+				if (flag < 0) flag = 0;
+				c1 = flag&3, c2 = flag>>2&3;
+				if (c1 == 2 && c2 == 0) type = 1;
+				else if (c2 == 2 && c1 == 0) type = 2;
+				b->s[i] = type;
+
+
+				///++t->cnt[k].c[flag];
+				++cnt.c[flag];
+
+				// if (aux->print_diff && (flag>>2&3) != (flag&3))
+				// 	printf("D\t%s\t%d\t%d\t%d\n", s->name, i, flag&3, flag>>2&3);
+			}
+		} else l = 0, x[0] = x[1] = x[2] = x[3] = 0;
+	}
+	for (l = 0, i = 1; i <= s->length; ++i) {
+		if (i == s->length || b->s[i] != b->s[l]) {
+			if (b->s[l] > 0 && i - l >= aux->k - 4)
+			{
+				///t->cnt[k].sc[b->s[l] - 1] += i - l;
+				cnt.sc[b->s[l] - 1] += i - l;
+			}
+			l = i;
+		}
+	}
+
+
+	int *c = cnt.c;
+	char type;
+	type = tb_classify(cnt.sc, c, aux->k, aux->ratio_thres);
+	aux->seq->trio_flag[k] = AMBIGU;
+	if(type == 'p') aux->seq->trio_flag[k] = FATHER;
+	if(type == 'm') aux->seq->trio_flag[k] = MOTHER;
+}
+
+
+
+void trio_partition()
+{
+    if(asm_opt.pat_index == NULL || asm_opt.mat_index == NULL)
+	{
+		memset(R_INF.trio_flag, AMBIGU, R_INF.total_reads*sizeof(uint8_t));
+		return;
+	}
+
+	double start_time = Get_T();
+	fprintf(stderr, "Start trio binning ...... \n");
+
+    yak_ch_t *ch;
+    int i/**, min_cnt = 2, mid_cnt = 5**/;
+    tb_shared_t aux;
+    memset(&aux, 0, sizeof(tb_shared_t));
+	aux.n_threads = asm_opt.thread_num, aux.print_diff = 0;
+	aux.ratio_thres = 0.33;
+	aux.seq = &R_INF;
+
+
+    ch = yak_ch_restore_core(0,  asm_opt.pat_index, YAK_LOAD_TRIOBIN1, asm_opt.min_cnt, asm_opt.mid_cnt);
+	ch = yak_ch_restore_core(ch, asm_opt.mat_index, YAK_LOAD_TRIOBIN2, asm_opt.min_cnt, asm_opt.mid_cnt);
+
+
+
+    aux.k = ch->k;
+    aux.ch = ch;
+	aux.buf = (tb_buf_t*)calloc(aux.n_threads, sizeof(tb_buf_t));
+	aux.bseq = (UC_Read*)calloc(aux.n_threads, sizeof(UC_Read));
+	for (i = 0; i < aux.n_threads; ++i)
+	{
+		init_UC_Read(&aux.bseq[i]);
+	} 
+	
+	kt_for(aux.n_threads, tb_worker, &aux, aux.seq->total_reads);
+	
+    for (i = 0; i < aux.n_threads; ++i)
+	{
+		free(aux.buf[i].s);
+		destory_UC_Read(&aux.bseq[i]);
+	} 
+	free(aux.buf);
+	free(aux.bseq);
+    yak_ch_destroy(ch);
+
+	fprintf(stderr, "Trio binning has been done.\n");
+	fprintf(stderr, "%-30s%18.2f\n\n", "Trio binning time:", Get_T() - start_time); 
+}
--- a/Trio.h
+++ b/Trio.h
@@ -0,0 +1,71 @@
+#ifndef __TRIO__
+#define __TRIO__
+#include <stdint.h>
+
+#define YAK_MAX_KMER     31
+#define YAK_COUNTER_BITS 10
+#define YAK_N_COUNTS     (1<<YAK_COUNTER_BITS)
+#define YAK_MAX_COUNT    ((1<<YAK_COUNTER_BITS)-1)
+
+#define YAK_BLK_SHIFT  9 // 64 bytes, the size of a cache line
+#define YAK_BLK_MASK   ((1<<(YAK_BLK_SHIFT)) - 1)
+
+#define YAK_LOAD_ALL       1
+#define YAK_LOAD_TRIOBIN1  2
+#define YAK_LOAD_TRIOBIN2  3
+
+#define YAK_MAGIC "YAK\2"
+
+typedef struct {
+	int n_shift, n_hashes;
+	uint8_t *b;
+} yak_bf_t;
+
+struct yak_ht_t;
+
+typedef struct {
+	struct yak_ht_t *h;
+	yak_bf_t *b;
+} yak_ch1_t;
+
+typedef struct {
+	int k, pre, n_hash, n_shift;
+	uint64_t tot;
+	yak_ch1_t *h;
+} yak_ch_t;
+
+
+typedef struct {
+	int max;
+	uint32_t *s;
+} tb_buf_t;
+
+
+typedef struct {
+	int k, n_threads, print_diff;
+	double ratio_thres;
+	///bseq_file_t *fp;
+	const yak_ch_t *ch;
+	tb_buf_t *buf;
+	UC_Read *bseq;
+	All_reads* seq; 
+} tb_shared_t;
+
+typedef struct {
+	int c[16];
+	int sc[2];
+	int nk;
+} tb_cnt_t;
+
+
+typedef struct {
+	int n_seq;
+	tb_shared_t *aux;
+	///bseq1_t *seq;
+	///All_reads* seq;
+	///tb_cnt_t *cnt;
+} tb_step_t;
+
+
+void trio_partition();
+#endif
--- a/hifiasm.1
+++ b/hifiasm.1
@@ -1,4 +1,4 @@
-.TH hifiasm 1 "3 Jan 2020" "hifiasm-0.1.0" "Bioinformatics tools"
+.TH hifiasm 1 "22 Mar 2020" "hifiasm-0.3.0" "Bioinformatics tools"

 .SH NAME
 .PP
@@ -164,6 +164,47 @@ This might be helpful when users want to get an optimized assembly by multiple r
 with different parameters.


+.SS Trio-partition options
+
+.TP 10
+.BI -P \ FILE
+Paternal trio index. This index should be generated by
+.I [yak count]
+with the paternal short reads. For details of yak, please see
+.I [https://github.com/lh3/yak]
+
+
+.TP 10
+.BI -M \ FILE
+Maternal trio index. This index should be generated by
+.I [yak count]
+with the maternal short reads. For details of yak, please see
+.I [https://github.com/lh3/yak]
+
+.TP 10
+.BI -c \ INT
+Lower bound of the binned k-mer's frequency [2]. When doing trio binning, 
+a k-mer is said to be differentiating if it occurs >=
+.I [-d upper_bound] 
+times in one sample 
+but occurs <
+.I [-c lower_bound]  
+times in the other sample.
+
+.TP 10
+.BI -d \ INT
+Upper bound of the binned k-mer's frequency [5]. When doing trio binning, 
+a k-mer is said to be differentiating if it occurs >=
+.I [-d upper_bound] 
+times in one sample 
+but occurs <
+.I [-c lower_bound]  
+times in the other sample.
+
+
+
+
+
 .SH EXAMPLES

 .TP
@@ -180,6 +221,22 @@ With
 .I [-z 20],
 hifiasm will remove 20 bases from both ends of each read.

+.SH EXAMPLES FRO TRIO
+.TP
+.BR ./yak " " count " " \-k31 " " \-b37 " " \-t16 " " \-o " " mat.yak " " mat.fq.gz
+Build maternal trio index from mat.fq.gz.
+
+.TP
+.BR ./yak " " count " " \-k31 " " \-b37 " " \-t16 " " \-o " " pat.yak " " pat.fq.gz
+Build paternal trio index from pat.fq.gz.
+
+.TP
+.BR ./hifiasm " " \-o " " NA12878.asm " " \-t " " 32 " " \-P " " pat.yak " " \-M " " mat.yak " " NA12878_1.fq.gz " " NA12878_2.fq.gz
+In this example, hifiasm will do trio assembly with 32 CPU threads. The paternal assembly can be found at [NA12878.asm.p.r_utg.gfa],
+and the maternal assembly can be found at [NA12878.asm.m.r_utg.gfa].
+
+
+
 .SH OUTPUTS


@@ -194,7 +251,7 @@ During the error correction step, hifiasm outputs the following two files:
 2. All-to-all overlaps in paf format [outPrefix.ovlp.paf].

 .PP
-During the assembly step, hifiasm outputs the following four assembly graphs in GFA format:
+During the non-trio assembly step, hifiasm outputs the following four assembly graphs in GFA format:


 .IP
@@ -214,6 +271,21 @@ This graph collapses different haplotypes.
 4. Alternate assembly contig graph [outPrefix.a_ctg.gfa].
 This graph consists of all assemblies that are discarded in primary assembly contig graph.

+
+.PP
+If you have trio information, hifiasm outputs the following three assembly graphs in GFA format:
+
+.IP
+1. Phased maternal unitig graph [outPrefix.m.r_utg.gfa].
+This graph keeps the phased maternal assembly.
+
+2. Phased paternal unitig graph [outPrefix.p.r_utg.gfa].
+This graph keeps the phased paternal assembly.
+
+3. Haplotype-resolved raw unitig graph [outPrefix.r_utg.gfa].
+This graph keeps all haplotype information.
+
+
 .PP
 For each graph, hifiasm also outputs a simplified version without sequences. These simplified
 graphs can be easily visualized.
--- a/khashl.h
+++ b/khashl.h
@@ -0,0 +1,351 @@
+/* The MIT License
+
+   Copyright (c) 2019 by Attractive Chaos <attractor@live.co.uk>
+
+   Permission is hereby granted, free of charge, to any person obtaining
+   a copy of this software and associated documentation files (the
+   "Software"), to deal in the Software without restriction, including
+   without limitation the rights to use, copy, modify, merge, publish,
+   distribute, sublicense, and/or sell copies of the Software, and to
+   permit persons to whom the Software is furnished to do so, subject to
+   the following conditions:
+
+   The above copyright notice and this permission notice shall be
+   included in all copies or substantial portions of the Software.
+
+   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+   EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+   MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+   NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
+   BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
+   ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
+   CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+   SOFTWARE.
+*/
+
+#ifndef __AC_KHASHL_H
+#define __AC_KHASHL_H
+
+#define AC_VERSION_KHASHL_H "0.1"
+
+#include <stdlib.h>
+#include <string.h>
+#include <limits.h>
+
+/************************************
+ * Compiler specific configurations *
+ ************************************/
+
+#if UINT_MAX == 0xffffffffu
+typedef unsigned int khint32_t;
+#elif ULONG_MAX == 0xffffffffu
+typedef unsigned long khint32_t;
+#endif
+
+#if ULONG_MAX == ULLONG_MAX
+typedef unsigned long khint64_t;
+#else
+typedef unsigned long long khint64_t;
+#endif
+
+#ifndef kh_inline
+#ifdef _MSC_VER
+#define kh_inline __inline
+#else
+#define kh_inline inline
+#endif
+#endif /* kh_inline */
+
+#ifndef klib_unused
+#if (defined __clang__ && __clang_major__ >= 3) || (defined __GNUC__ && __GNUC__ >= 3)
+#define klib_unused __attribute__ ((__unused__))
+#else
+#define klib_unused
+#endif
+#endif /* klib_unused */
+
+#define KH_LOCAL static kh_inline klib_unused
+
+typedef khint32_t khint_t;
+
+/******************
+ * malloc aliases *
+ ******************/
+
+#ifndef kcalloc
+#define kcalloc(N,Z) calloc(N,Z)
+#endif
+#ifndef kmalloc
+#define kmalloc(Z) malloc(Z)
+#endif
+#ifndef krealloc
+#define krealloc(P,Z) realloc(P,Z)
+#endif
+#ifndef kfree
+#define kfree(P) free(P)
+#endif
+
+/****************************
+ * Simple private functions *
+ ****************************/
+
+#define __kh_used(flag, i)       (flag[i>>5] >> (i&0x1fU) & 1U)
+#define __kh_set_used(flag, i)   (flag[i>>5] |= 1U<<(i&0x1fU))
+#define __kh_set_unused(flag, i) (flag[i>>5] &= ~(1U<<(i&0x1fU)))
+
+#define __kh_fsize(m) ((m) < 32? 1 : (m)>>5)
+
+static kh_inline khint_t __kh_h2b(khint_t hash, khint_t bits) { return hash * 2654435769U >> (32 - bits); }
+
+/*******************
+ * Hash table base *
+ *******************/
+
+#define __KHASHL_TYPE(HType, khkey_t) \
+	typedef struct HType { \
+		khint_t bits, count; \
+		khint32_t *used; \
+		khkey_t *keys; \
+	} HType;
+
+#define __KHASHL_PROTOTYPES(HType, prefix, khkey_t) \
+	extern HType *prefix##_init(void); \
+	extern void prefix##_destroy(HType *h); \
+	extern void prefix##_clear(HType *h); \
+	extern khint_t prefix##_getp(const HType *h, const khkey_t *key); \
+	extern int prefix##_resize(HType *h, khint_t new_n_buckets); \
+	extern khint_t prefix##_putp(HType *h, const khkey_t *key, int *absent); \
+	extern void prefix##_del(HType *h, khint_t k);
+
+#define __KHASHL_IMPL_BASIC(SCOPE, HType, prefix) \
+	SCOPE HType *prefix##_init(void) { \
+		return (HType*)kcalloc(1, sizeof(HType)); \
+	} \
+	SCOPE void prefix##_destroy(HType *h) { \
+		if (!h) return; \
+		kfree((void *)h->keys); kfree(h->used); \
+		kfree(h); \
+	} \
+	SCOPE void prefix##_clear(HType *h) { \
+		if (h && h->used) { \
+			uint32_t n_buckets = 1U << h->bits; \
+			memset(h->used, 0, __kh_fsize(n_buckets) * sizeof(khint32_t)); \
+			h->count = 0; \
+		} \
+	}
+
+#define __KHASHL_IMPL_GET(SCOPE, HType, prefix, khkey_t, __hash_fn, __hash_eq) \
+	SCOPE khint_t prefix##_getp(const HType *h, const khkey_t *key) { \
+		khint_t i, last, n_buckets, mask; \
+		if (h->keys == 0) return 0; \
+		n_buckets = 1U << h->bits; \
+		mask = n_buckets - 1U; \
+		i = last = __kh_h2b(__hash_fn(*key), h->bits); \
+		while (__kh_used(h->used, i) && !__hash_eq(h->keys[i], *key)) { \
+			i = (i + 1U) & mask; \
+			if (i == last) return n_buckets; \
+		} \
+		return !__kh_used(h->used, i)? n_buckets : i; \
+	} \
+	SCOPE khint_t prefix##_get(const HType *h, khkey_t key) { return prefix##_getp(h, &key); }
+
+#define __KHASHL_IMPL_RESIZE(SCOPE, HType, prefix, khkey_t, __hash_fn, __hash_eq) \
+	SCOPE int prefix##_resize(HType *h, khint_t new_n_buckets) { \
+		khint32_t *new_used = 0; \
+		khint_t j = 0, x = new_n_buckets, n_buckets, new_bits, new_mask; \
+		while ((x >>= 1) != 0) ++j; \
+		if (new_n_buckets & (new_n_buckets - 1)) ++j; \
+		new_bits = j > 2? j : 2; \
+		new_n_buckets = 1U << new_bits; \
+		if (h->count > (new_n_buckets>>1) + (new_n_buckets>>2)) return 0; /* requested size is too small */ \
+		new_used = (khint32_t*)kmalloc(__kh_fsize(new_n_buckets) * sizeof(khint32_t)); \
+		memset(new_used, 0, __kh_fsize(new_n_buckets) * sizeof(khint32_t)); \
+		if (!new_used) return -1; /* not enough memory */ \
+		n_buckets = h->keys? 1U<<h->bits : 0U; \
+		if (n_buckets < new_n_buckets) { /* expand */ \
+			khkey_t *new_keys = (khkey_t*)krealloc((void*)h->keys, new_n_buckets * sizeof(khkey_t)); \
+			if (!new_keys) { kfree(new_used); return -1; } \
+			h->keys = new_keys; \
+		} /* otherwise shrink */ \
+		new_mask = new_n_buckets - 1; \
+		for (j = 0; j != n_buckets; ++j) { \
+			khkey_t key; \
+			if (!__kh_used(h->used, j)) continue; \
+			key = h->keys[j]; \
+			__kh_set_unused(h->used, j); \
+			while (1) { /* kick-out process; sort of like in Cuckoo hashing */ \
+				khint_t i; \
+				i = __kh_h2b(__hash_fn(key), new_bits); \
+				while (__kh_used(new_used, i)) i = (i + 1) & new_mask; \
+				__kh_set_used(new_used, i); \
+				if (i < n_buckets && __kh_used(h->used, i)) { /* kick out the existing element */ \
+					{ khkey_t tmp = h->keys[i]; h->keys[i] = key; key = tmp; } \
+					__kh_set_unused(h->used, i); /* mark it as deleted in the old hash table */ \
+				} else { /* write the element and jump out of the loop */ \
+					h->keys[i] = key; \
+					break; \
+				} \
+			} \
+		} \
+		if (n_buckets > new_n_buckets) /* shrink the hash table */ \
+			h->keys = (khkey_t*)krealloc((void *)h->keys, new_n_buckets * sizeof(khkey_t)); \
+		kfree(h->used); /* free the working space */ \
+		h->used = new_used, h->bits = new_bits; \
+		return 0; \
+	}
+
+#define __KHASHL_IMPL_PUT(SCOPE, HType, prefix, khkey_t, __hash_fn, __hash_eq) \
+	SCOPE khint_t prefix##_putp(HType *h, const khkey_t *key, int *absent) { \
+		khint_t n_buckets, i, last, mask; \
+		n_buckets = h->keys? 1U<<h->bits : 0U; \
+		*absent = -1; \
+		if (h->count >= (n_buckets>>1) + (n_buckets>>2)) { /* rehashing */ \
+			if (prefix##_resize(h, n_buckets + 1U) < 0) \
+				return n_buckets; \
+			n_buckets = 1U<<h->bits; \
+		} /* TODO: to implement automatically shrinking; resize() already support shrinking */ \
+		mask = n_buckets - 1; \
+		i = last = __kh_h2b(__hash_fn(*key), h->bits); \
+		while (__kh_used(h->used, i) && !__hash_eq(h->keys[i], *key)) { \
+			i = (i + 1U) & mask; \
+			if (i == last) break; \
+		} \
+		if (!__kh_used(h->used, i)) { /* not present at all */ \
+			h->keys[i] = *key; \
+			__kh_set_used(h->used, i); \
+			++h->count; \
+			*absent = 1; \
+		} else *absent = 0; /* Don't touch h->keys[i] if present */ \
+		return i; \
+	} \
+	SCOPE khint_t prefix##_put(HType *h, khkey_t key, int *absent) { return prefix##_putp(h, &key, absent); }
+
+#define __KHASHL_IMPL_DEL(SCOPE, HType, prefix, khkey_t, __hash_fn) \
+	SCOPE int prefix##_del(HType *h, khint_t i) { \
+		khint_t j = i, k, mask, n_buckets; \
+		if (h->keys == 0) return 0; \
+		n_buckets = 1U<<h->bits; \
+		mask = n_buckets - 1U; \
+		while (1) { \
+			j = (j + 1U) & mask; \
+			if (j == i || !__kh_used(h->used, j)) break; /* j==i only when the table is completely full */ \
+			k = __kh_h2b(__hash_fn(h->keys[j]), h->bits); \
+			if ((j > i && (k <= i || k > j)) || (j < i && (k <= i && k > j))) \
+				h->keys[i] = h->keys[j], i = j; \
+		} \
+		__kh_set_unused(h->used, i); \
+		--h->count; \
+		return 1; \
+	}
+
+#define KHASHL_DECLARE(HType, prefix, khkey_t) \
+	__KHASHL_TYPE(HType, khkey_t) \
+	__KHASHL_PROTOTYPES(HType, prefix, khkey_t)
+
+#define KHASHL_INIT(SCOPE, HType, prefix, khkey_t, __hash_fn, __hash_eq) \
+	__KHASHL_TYPE(HType, khkey_t) \
+	__KHASHL_IMPL_BASIC(SCOPE, HType, prefix) \
+	__KHASHL_IMPL_GET(SCOPE, HType, prefix, khkey_t, __hash_fn, __hash_eq) \
+	__KHASHL_IMPL_RESIZE(SCOPE, HType, prefix, khkey_t, __hash_fn, __hash_eq) \
+	__KHASHL_IMPL_PUT(SCOPE, HType, prefix, khkey_t, __hash_fn, __hash_eq) \
+	__KHASHL_IMPL_DEL(SCOPE, HType, prefix, khkey_t, __hash_fn)
+
+/*****************************
+ * More convenient interface *
+ *****************************/
+
+#define __kh_packed __attribute__ ((__packed__))
+#define __kh_cached_hash(x) ((x).hash)
+
+#define KHASHL_SET_INIT(SCOPE, HType, prefix, khkey_t, __hash_fn, __hash_eq) \
+	typedef struct { khkey_t key; } __kh_packed HType##_s_bucket_t; \
+	static kh_inline khint_t prefix##_s_hash(HType##_s_bucket_t x) { return __hash_fn(x.key); } \
+	static kh_inline int prefix##_s_eq(HType##_s_bucket_t x, HType##_s_bucket_t y) { return __hash_eq(x.key, y.key); } \
+	KHASHL_INIT(KH_LOCAL, HType, prefix##_s, HType##_s_bucket_t, prefix##_s_hash, prefix##_s_eq) \
+	SCOPE HType *prefix##_init(void) { return prefix##_s_init(); } \
+	SCOPE void prefix##_destroy(HType *h) { prefix##_s_destroy(h); } \
+	SCOPE void prefix##_resize(HType *h, khint_t new_n_buckets) { prefix##_s_resize(h, new_n_buckets); } \
+	SCOPE khint_t prefix##_get(const HType *h, khkey_t key) { HType##_s_bucket_t t; t.key = key; return prefix##_s_getp(h, &t); } \
+	SCOPE int prefix##_del(HType *h, khint_t k) { return prefix##_s_del(h, k); } \
+	SCOPE khint_t prefix##_put(HType *h, khkey_t key, int *absent) { HType##_s_bucket_t t; t.key = key; return prefix##_s_putp(h, &t, absent); }
+
+#define KHASHL_MAP_INIT(SCOPE, HType, prefix, khkey_t, kh_val_t, __hash_fn, __hash_eq) \
+	typedef struct { khkey_t key; kh_val_t val; } __kh_packed HType##_m_bucket_t; \
+	static kh_inline khint_t prefix##_m_hash(HType##_m_bucket_t x) { return __hash_fn(x.key); } \
+	static kh_inline int prefix##_m_eq(HType##_m_bucket_t x, HType##_m_bucket_t y) { return __hash_eq(x.key, y.key); } \
+	KHASHL_INIT(KH_LOCAL, HType, prefix##_m, HType##_m_bucket_t, prefix##_m_hash, prefix##_m_eq) \
+	SCOPE HType *prefix##_init(void) { return prefix##_m_init(); } \
+	SCOPE void prefix##_destroy(HType *h) { prefix##_m_destroy(h); } \
+	SCOPE khint_t prefix##_get(const HType *h, khkey_t key) { HType##_m_bucket_t t; t.key = key; return prefix##_m_getp(h, &t); } \
+	SCOPE int prefix##_del(HType *h, khint_t k) { return prefix##_m_del(h, k); } \
+	SCOPE khint_t prefix##_put(HType *h, khkey_t key, int *absent) { HType##_m_bucket_t t; t.key = key; return prefix##_m_putp(h, &t, absent); }
+
+#define KHASHL_CSET_INIT(SCOPE, HType, prefix, khkey_t, __hash_fn, __hash_eq) \
+	typedef struct { khkey_t key; khint_t hash; } __kh_packed HType##_cs_bucket_t; \
+	static kh_inline int prefix##_cs_eq(HType##_cs_bucket_t x, HType##_cs_bucket_t y) { return x.hash == y.hash && __hash_eq(x.key, y.key); } \
+	KHASHL_INIT(KH_LOCAL, HType, prefix##_cs, HType##_cs_bucket_t, __kh_cached_hash, prefix##_cs_eq) \
+	SCOPE HType *prefix##_init(void) { return prefix##_cs_init(); } \
+	SCOPE void prefix##_destroy(HType *h) { prefix##_cs_destroy(h); } \
+	SCOPE khint_t prefix##_get(const HType *h, khkey_t key) { HType##_cs_bucket_t t; t.key = key; t.hash = __hash_fn(key); return prefix##_cs_getp(h, &t); } \
+	SCOPE int prefix##_del(HType *h, khint_t k) { return prefix##_cs_del(h, k); } \
+	SCOPE khint_t prefix##_put(HType *h, khkey_t key, int *absent) { HType##_cs_bucket_t t; t.key = key, t.hash = __hash_fn(key); return prefix##_cs_putp(h, &t, absent); }
+
+#define KHASHL_CMAP_INIT(SCOPE, HType, prefix, khkey_t, kh_val_t, __hash_fn, __hash_eq) \
+	typedef struct { khkey_t key; kh_val_t val; khint_t hash; } __kh_packed HType##_cm_bucket_t; \
+	static kh_inline int prefix##_cm_eq(HType##_cm_bucket_t x, HType##_cm_bucket_t y) { return x.hash == y.hash && __hash_eq(x.key, y.key); } \
+	KHASHL_INIT(KH_LOCAL, HType, prefix##_cm, HType##_cm_bucket_t, __kh_cached_hash, prefix##_cm_eq) \
+	SCOPE HType *prefix##_init(void) { return prefix##_cm_init(); } \
+	SCOPE void prefix##_destroy(HType *h) { prefix##_cm_destroy(h); } \
+	SCOPE khint_t prefix##_get(const HType *h, khkey_t key) { HType##_cm_bucket_t t; t.key = key; t.hash = __hash_fn(key); return prefix##_cm_getp(h, &t); } \
+	SCOPE int prefix##_del(HType *h, khint_t k) { return prefix##_cm_del(h, k); } \
+	SCOPE khint_t prefix##_put(HType *h, khkey_t key, int *absent) { HType##_cm_bucket_t t; t.key = key, t.hash = __hash_fn(key); return prefix##_cm_putp(h, &t, absent); }
+
+/**************************
+ * Public macro functions *
+ **************************/
+
+#define kh_bucket(h, x) ((h)->keys[x])
+#define kh_size(h) ((h)->count)
+#define kh_capacity(h) ((h)->keys? 1U<<(h)->bits : 0U)
+#define kh_end(h) kh_capacity(h)
+
+#define kh_key(h, x) ((h)->keys[x].key)
+#define kh_val(h, x) ((h)->keys[x].val)
+#define kh_exist(h, x) __kh_used((h)->used, (x))
+
+/**************************************
+ * Common hash and equality functions *
+ **************************************/
+
+#define kh_eq_generic(a, b) ((a) == (b))
+#define kh_eq_str(a, b) (strcmp((a), (b)) == 0)
+#define kh_hash_dummy(x) ((khint_t)(x))
+
+static kh_inline khint_t kh_hash_uint32(khint_t key) {
+	key += ~(key << 15);
+	key ^=  (key >> 10);
+	key +=  (key << 3);
+	key ^=  (key >> 6);
+	key += ~(key << 11);
+	key ^=  (key >> 16);
+	return key;
+}
+
+static kh_inline khint_t kh_hash_uint64(khint64_t key) {
+	key = ~key + (key << 21);
+	key = key ^ key >> 24;
+	key = (key + (key << 3)) + (key << 8);
+	key = key ^ key >> 14;
+	key = (key + (key << 2)) + (key << 4);
+	key = key ^ key >> 28;
+	key = key + (key << 31);
+	return (khint_t)key;
+}
+
+static kh_inline khint_t kh_hash_str(const char *s) {
+	khint_t h = (khint_t)*s;
+	if (h) for (++s ; *s; ++s) h = (h << 5) - h + (khint_t)*s;
+	return h;
+}
+
+#endif /* __AC_KHASHL_H */
--- a/kthread.cpp
+++ b/kthread.cpp
@@ -0,0 +1,159 @@
+#include <pthread.h>
+#include <stdlib.h>
+#include <limits.h>
+#include <stdint.h>
+#include "kthread.h"
+
+#if (defined(WIN32) || defined(_WIN32)) && defined(_MSC_VER)
+#define __sync_fetch_and_add(ptr, addend)     _InterlockedExchangeAdd((void*)ptr, addend)
+#endif
+
+/************
+ * kt_for() *
+ ************/
+
+struct kt_for_t;
+
+typedef struct {
+	struct kt_for_t *t;
+	long i;
+} ktf_worker_t;
+
+typedef struct kt_for_t {
+	int n_threads;
+	long n;
+	ktf_worker_t *w;
+	void (*func)(void*,long,int);
+	void *data;
+} kt_for_t;
+
+static inline long steal_work(kt_for_t *t)
+{
+	int i, min_i = -1;
+	long k, min = LONG_MAX;
+	for (i = 0; i < t->n_threads; ++i)
+		if (min > t->w[i].i) min = t->w[i].i, min_i = i;
+	k = __sync_fetch_and_add(&t->w[min_i].i, t->n_threads);
+	return k >= t->n? -1 : k;
+}
+
+static void *ktf_worker(void *data)
+{
+	ktf_worker_t *w = (ktf_worker_t*)data;
+	long i;
+	for (;;) {
+		i = __sync_fetch_and_add(&w->i, w->t->n_threads);
+		if (i >= w->t->n) break;
+		w->t->func(w->t->data, i, w - w->t->w);
+	}
+	while ((i = steal_work(w->t)) >= 0)
+		w->t->func(w->t->data, i, w - w->t->w);
+	pthread_exit(0);
+}
+
+void kt_for(int n_threads, void (*func)(void*,long,int), void *data, long n)
+{
+	if (n_threads > 1) {
+		int i;
+		kt_for_t t;
+		pthread_t *tid;
+		t.func = func, t.data = data, t.n_threads = n_threads, t.n = n;
+		t.w = (ktf_worker_t*)calloc(n_threads, sizeof(ktf_worker_t));
+		tid = (pthread_t*)calloc(n_threads, sizeof(pthread_t));
+		for (i = 0; i < n_threads; ++i)
+			t.w[i].t = &t, t.w[i].i = i;
+		for (i = 0; i < n_threads; ++i) pthread_create(&tid[i], 0, ktf_worker, &t.w[i]);
+		for (i = 0; i < n_threads; ++i) pthread_join(tid[i], 0);
+		free(tid); free(t.w);
+	} else {
+		long j;
+		for (j = 0; j < n; ++j) func(data, j, 0);
+	}
+}
+
+/*****************
+ * kt_pipeline() *
+ *****************/
+
+struct ktp_t;
+
+typedef struct {
+	struct ktp_t *pl;
+	int64_t index;
+	int step;
+	void *data;
+} ktp_worker_t;
+
+typedef struct ktp_t {
+	void *shared;
+	void *(*func)(void*, int, void*);
+	int64_t index;
+	int n_workers, n_steps;
+	ktp_worker_t *workers;
+	pthread_mutex_t mutex;
+	pthread_cond_t cv;
+} ktp_t;
+
+static void *ktp_worker(void *data)
+{
+	ktp_worker_t *w = (ktp_worker_t*)data;
+	ktp_t *p = w->pl;
+	while (w->step < p->n_steps) {
+		// test whether we can kick off the job with this worker
+		pthread_mutex_lock(&p->mutex);
+		for (;;) {
+			int i;
+			// test whether another worker is doing the same step
+			for (i = 0; i < p->n_workers; ++i) {
+				if (w == &p->workers[i]) continue; // ignore itself
+				if (p->workers[i].step <= w->step && p->workers[i].index < w->index)
+					break;
+			}
+			if (i == p->n_workers) break; // no workers with smaller indices are doing w->step or the previous steps
+			pthread_cond_wait(&p->cv, &p->mutex);
+		}
+		pthread_mutex_unlock(&p->mutex);
+
+		// working on w->step
+		w->data = p->func(p->shared, w->step, w->step? w->data : 0); // for the first step, input is NULL
+
+		// update step and let other workers know
+		pthread_mutex_lock(&p->mutex);
+		w->step = w->step == p->n_steps - 1 || w->data? (w->step + 1) % p->n_steps : p->n_steps;
+		if (w->step == 0) w->index = p->index++;
+		pthread_cond_broadcast(&p->cv);
+		pthread_mutex_unlock(&p->mutex);
+	}
+	pthread_exit(0);
+}
+
+void kt_pipeline(int n_threads, void *(*func)(void*, int, void*), void *shared_data, int n_steps)
+{
+	ktp_t aux;
+	pthread_t *tid;
+	int i;
+
+	if (n_threads < 1) n_threads = 1;
+	aux.n_workers = n_threads;
+	aux.n_steps = n_steps;
+	aux.func = func;
+	aux.shared = shared_data;
+	aux.index = 0;
+	pthread_mutex_init(&aux.mutex, 0);
+	pthread_cond_init(&aux.cv, 0);
+
+	aux.workers = (ktp_worker_t*)calloc(n_threads, sizeof(ktp_worker_t));
+	for (i = 0; i < n_threads; ++i) {
+		ktp_worker_t *w = &aux.workers[i];
+		w->step = 0; w->pl = &aux; w->data = 0;
+		w->index = aux.index++;
+	}
+
+	tid = (pthread_t*)calloc(n_threads, sizeof(pthread_t));
+	for (i = 0; i < n_threads; ++i) pthread_create(&tid[i], 0, ktp_worker, &aux.workers[i]);
+	for (i = 0; i < n_threads; ++i) pthread_join(tid[i], 0);
+	free(tid); free(aux.workers);
+
+	pthread_mutex_destroy(&aux.mutex);
+	pthread_cond_destroy(&aux.cv);
+}
--- a/kthread.h
+++ b/kthread.h
@@ -0,0 +1,15 @@
+#ifndef KTHREAD_H
+#define KTHREAD_H
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+void kt_for(int n_threads, void (*func)(void*,long,int), void *data, long n);
+void kt_pipeline(int n_threads, void *(*func)(void*, int, void*), void *shared_data, int n_steps);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif
Author	SHA1	Message	Date
chhylp123	3c1d3cf6a1	disable debug	2020-03-23 00:25:20 -04:00
chhylp123	ee2573a05c	Merge branch 'master' of https://github.com/chhylp123/Long_read_assembly	2020-03-22 18:32:47 -04:00
chhylp123	d3cb016aba	Update README	2020-03-22 18:31:59 -04:00
chhylp123	670bd10093	update man	2020-03-22 18:01:04 -04:00
chhylp123	90e290636a	code clean	2020-03-22 17:06:36 -04:00
chhylp123	c72d419711	for hifiasm 3.0	2020-03-22 16:54:32 -04:00
chhylp123	4536141b9a	update README for Redwood update README for Redwood	2020-03-04 14:58:04 -05:00