mempcpy -> memcpy

Assembly.cpp

@@ -425,6 +425,7 @@ typedef struct {
 	overlap_region_alloc olist;
     overlap_region_alloc olist_hp;
 	ha_abuf_t *ab;
+    ha_abufl_t *abl;
 	// error correction related buffers
 	int64_t num_read_base, num_correct_base, num_recorrect_base;
 	Cigar_record cigar1;
@@ -441,7 +442,7 @@ typedef struct {
     st_mt_t sp;
 } ha_ovec_buf_t;
 
-ha_ovec_buf_t *ha_ovec_init(int is_final, int save_ov)
+ha_ovec_buf_t *ha_ovec_init(int is_final, int save_ov, int is_ug)
 {
 	ha_ovec_buf_t *b;
 	CALLOC(b, 1);
@@ -456,7 +457,8 @@ ha_ovec_buf_t *ha_ovec_init(int is_final, int save_ov)
     kv_init(b->r_buf.a);
     kv_init(b->k_flag.a);
     kv_init(b->sp);
-	b->ab = ha_abuf_init();
+	if(!is_ug) b->ab = ha_abuf_init();
+    else b->abl = ha_abufl_init();
 	if (!b->is_final) {
 		init_Cigar_record(&b->cigar1);
 		init_Graph(&b->POA_Graph);
@@ -476,6 +478,7 @@ void ha_ovec_destroy(ha_ovec_buf_t *b)
 	destory_overlap_region_alloc(&b->olist);
     destory_overlap_region_alloc(&b->olist_hp);
 	ha_abuf_destroy(b->ab);
+    ha_abufl_destroy(b->abl);
     destory_fake_cigar(&(b->tmp_region.f_cigar));
     kv_destroy(b->b_buf.a);
     kv_destroy(b->r_buf.a);
@@ -526,7 +529,8 @@ int64_t ha_ovec_mem(const ha_ovec_buf_t *b)
 		mem_olist += r->boundary_cigars.size * sizeof(window_list);
 	}
 
-	mem = ha_abuf_mem(b->ab) + mem_clist + mem_olist;
+    if(b->ab) mem = ha_abuf_mem(b->ab) + mem_clist + mem_olist;
+    if(b->abl) mem = ha_abufl_mem(b->abl) + mem_clist + mem_olist;
 	if (!b->is_final) {
 		mem += sizeof(Cigar_record) + b->cigar1.lost_base_size + b->cigar1.size * 4;
 		mem += sizeof(Correct_dumy) + b->correct.size * 8;
@@ -852,7 +856,7 @@ void ha_overlap_and_correct(int round)
 	// overlap and correct reads
 	CALLOC(b, asm_opt.thread_num);
 	for (i = 0; i < asm_opt.thread_num; ++i)
-		b[i] = ha_ovec_init(0, (round == asm_opt.number_of_round - 1));
+		b[i] = ha_ovec_init(0, (round == asm_opt.number_of_round - 1),0);
     if(ha_idx) hom_cov = asm_opt.hom_cov;
 	if(ha_idx == NULL) ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, round == 0? 0 : 1, 0, &R_INF, &hom_cov, &het_cov); // build the index
 	///debug_adapter(&asm_opt, &R_INF);
@@ -1588,7 +1592,7 @@ void ha_overlap_final(void)
 
 	CALLOC(b, asm_opt.thread_num);
 	for (i = 0; i < asm_opt.thread_num; ++i)
-		b[i] = ha_ovec_init(asm_opt.flag & HA_F_HIGH_HET, 1);///b[i] = ha_ovec_init(1, 1);
+		b[i] = ha_ovec_init(asm_opt.flag & HA_F_HIGH_HET, 1,0);///b[i] = ha_ovec_init(1, 1);
 	ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, 1, 0, &R_INF, &hom_cov, &het_cov); // build the index
     if(asm_opt.flag & HA_F_HIGH_HET)
     {
@@ -1644,7 +1648,7 @@ void ug_idx_build(ma_ug_t *ug, int hap_n)
     CALLOC(b, asm_opt.thread_num);
     for (i = 0; i < asm_opt.thread_num; ++i)
     {
-        b[i] = ha_ovec_init(1, 1);
+        b[i] = ha_ovec_init(1, 1, 0);
         b[i]->ua = &(ug->u);
     }
             

CommandLines.h

@@ -4,7 +4,7 @@
 #include <pthread.h>
 #include <stdint.h>
 
-#define HA_VERSION "0.16.1-r374"
+#define HA_VERSION "0.16.1-r375"
 
 #define VERBOSE 0
 

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; int h;} ma_utg_v;
+typedef struct { size_t n, m; ma_utg_t *a;} ma_utg_v;
 
 typedef struct {
 	ma_utg_v u;

anchor.cpp

@@ -29,6 +29,11 @@ typedef struct {
 	const ha_idxpos_t *a;
 } seed1_t;
 
+typedef struct {
+	int n, cnt;
+	const ha_idxposl_t *a;
+} seedl_t;
+
 struct ha_abuf_s {
 	uint64_t n_a, m_a;///number of anchors (seed positions)
 	uint32_t old_mz_m;///number of seeds
@@ -37,20 +42,23 @@ struct ha_abuf_s {
 	anchor1_t *a;
 };
 
-ha_abuf_t *ha_abuf_init(void)
+struct ha_abufl_s {
-{
+	uint64_t n_a, m_a;///number of anchors (seed positions)
-	return (ha_abuf_t*)calloc(1, sizeof(ha_abuf_t));
+	uint32_t old_mz_m;///number of seeds
+	ha_mzl_v mz;
+	seedl_t *seed;
+	anchor1_t *a;
+};
+
+#define HA_ABUF_INIT(HType, MZType, SDType, sf) \
+HType *sf##_init(void){return (HType*)calloc(1, sizeof(HType));}\
+void sf##_destroy(HType *ab){if(ab){free(ab->seed); free(ab->a); free(ab->mz.a); free(ab);}}\
+uint64_t sf##_mem(const HType *ab){\
+	return ab->m_a * sizeof(anchor1_t) + ab->mz.m * (sizeof(MZType) + sizeof(SDType)) + sizeof(HType);\
 }
 
-void ha_abuf_destroy(ha_abuf_t *ab)
+HA_ABUF_INIT(ha_abuf_s, ha_mz1_t, seed1_t, ha_abuf)
-{
+HA_ABUF_INIT(ha_abufl_s, ha_mzl_t, seedl_t, ha_abufl)
-	free(ab->seed); free(ab->a); free(ab->mz.a); free(ab);
-}
-
-uint64_t ha_abuf_mem(const ha_abuf_t *ab)
-{
-	return ab->m_a * sizeof(anchor1_t) + ab->mz.m * (sizeof(ha_mz1_t) + sizeof(seed1_t)) + sizeof(ha_abuf_t);
-}
 
 int ha_ov_type(const overlap_region *r, uint32_t len)
 {
@@ -257,122 +265,248 @@ void calculate_ug_chaining(Candidates_list* candidates, overlap_region_alloc* ov
     }
 }
 
+/**
+void ha_get_inter_candidates(ha_abufl_t *ab, int64_t uid, ma_utg_v *ua, overlap_region_alloc *ovlp, 
+Candidates_list *cl, double bw_thres, int max_n_chain, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* chain_idx, 
+void *ha_flt_tab, ha_pt_t *ha_idx, overlap_region* f_cigar, kvec_t_u64_warp* dbg_ct, st_mt_t *sp, 
+double chain_match_rate, uint32_t uk, uint32_t uw, uint32_t is_hpc, uint32_t h_occ)
+{
+
+	uint32_t i, rlen;
+    uint64_t k, l;
+	ma_utg_t *u = &(ua->a[uid]);
+	// prepare
+    clear_Candidates_list(cl);
+    clear_overlap_region_alloc(ovlp);
+	ab->mz.n = 0, ab->n_a = 0;
+
+	mz2_ha_sketch(u->s, u->len, uw, uk, uid, is_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, 1);
+	// minimizer of queried read
+	if (ab->mz.m > ab->old_mz_m) {
+        ab->old_mz_m = ab->mz.m;
+        REALLOC(ab->seed, ab->old_mz_m);
+    }
+    for (i = 0, ab->n_a = 0; i < ab->mz.n; ++i) {
+        int n;
+        ab->seed[i].a = ha_ptl_get(ha_idx, ab->mz.a[i].x, &n);
+        ab->seed[i].n = n;
+        ab->seed[i].cnt = ha_ft_cnt(ha_flt_tab, ab->mz.a[i].x);
+        ab->n_a += n;
+    }
+    if (ab->n_a > ab->m_a) {
+        ab->m_a = ab->n_a;
+        kroundup64(ab->m_a);
+        REALLOC(ab->a, ab->m_a);
+    }
+    for (i = 0, k = 0; i < ab->mz.n; ++i) {
+        int j;
+        ///z is one of the minimizer
+        ha_mzl_t *z = &ab->mz.a[i];
+        seedl_t *s = &ab->seed[i];
+        for (j = 0; j < s->n; ++j) {
+            const ha_idxposl_t *y = &s->a[j];
+            anchor1_t *an = &ab->a[k++];
+            uint8_t rev = z->rev == y->rev? 0 : 1;
+            an->other_off = y->pos;
+            an->self_off = rev? u->len - 1 - (z->pos + 1 - z->span) : z->pos;
+            an->cnt = s->cnt;
+            an->srt = (uint64_t)y->rid<<33 | (uint64_t)rev<<32 | an->other_off;
+        }
+    }
+
+    // sort anchors
+    radix_sort_ha_an1(ab->a, ab->a + ab->n_a);
+    for (k = 1, l = 0; k <= ab->n_a; ++k) {
+        if (k == ab->n_a || ab->a[k].srt != ab->a[l].srt) {
+            if (k - l > 1)
+                radix_sort_ha_an2(ab->a + l, ab->a + k);
+            l = k;
+        }
+    }
+
+
+    // copy over to _cl_
+    if (ab->m_a >= (uint64_t)cl->size) {
+        cl->size = ab->m_a;
+        REALLOC(cl->list, cl->size);
+    }
+    for (k = 0; k < ab->n_a; ++k) {
+        k_mer_hit *p = &cl->list[k];
+        p->readID = ab->a[k].srt >> 33;
+        p->strand = ab->a[k].srt >> 32 & 1;
+        p->offset = ab->a[k].other_off;
+        p->self_offset = ab->a[k].self_off;
+        p->cnt = (ab->a[k].cnt == 1? 1 : 16);
+    }
+    cl->length = ab->n_a;
+    // calculate_overlap_region_by_chaining(cl, ovlp, chain_idx, rid, ucr->length, &R_INF, bw_thres, keep_whole_chain, f_cigar);
+	// calculate_ug_chaining(cl, overlap_list, chain_idx, rid, ua, bw_thres, keep_whole_chain, f_cigar, ab->mz.n, chain_match_rate);
+
+
+    #if 0
+    if (overlap_list->length > 0) {
+        fprintf(stderr, "B\t%ld\t%ld\t%d\n", (long)rid, (long)overlap_list->length, rlen);
+        for (int i = 0; i < (int)overlap_list->length; ++i) {
+            overlap_region *r = &overlap_list->list[i];
+            fprintf(stderr, "C\t%d\t%d\t%d\t%c\t%d\t%ld\t%d\t%d\t%c\t%d\t%d\n", (int)r->x_id, (int)r->x_pos_s, (int)r->x_pos_e, "+-"[r->x_pos_strand],
+                    (int)r->y_id, (long)Get_READ_LENGTH(R_INF, r->y_id), (int)r->y_pos_s, (int)r->y_pos_e, "+-"[r->y_pos_strand], (int)r->shared_seed, ha_ov_type(r, rlen));
+        }
+    }
+    #endif
+    
+    if ((int)ovlp->length > max_n_chain) {
+        int32_t w, n[4], s[4];
+        n[0] = n[1] = n[2] = n[3] = 0, s[0] = s[1] = s[2] = s[3] = 0;
+        ks_introsort_or_ss(ovlp->length, ovlp->list);
+        for (i = 0; i < (uint32_t)ovlp->length; ++i) {
+            const overlap_region *r = &ovlp->list[i];
+            w = ha_ov_type(r, rlen);
+            ++n[w];
+            if ((int)n[w] == max_n_chain) s[w] = r->shared_seed;
+        }
+        if (s[0] > 0 || s[1] > 0 || s[2] > 0 || s[3] > 0) {
+            // n[0] = n[1] = n[2] = n[3] = 0;
+            for (i = 0, k = 0; i < (uint32_t)ovlp->length; ++i) {
+                overlap_region *r = &ovlp->list[i];
+                w = ha_ov_type(r, rlen);
+                // ++n[w];
+                // if (((int)n[w] <= max_n_chain) || (r->shared_seed >= s[w] && s[w] >= (asm_opt.k_mer_length<<1))) {
+                if (r->shared_seed >= s[w]) {
+                    if ((uint32_t)k != i) {
+                        overlap_region t;
+                        t = ovlp->list[k];
+                        ovlp->list[k] = ovlp->list[i];
+                        ovlp->list[i] = t;
+                    }
+                    ++k;
+                }
+            }
+            ovlp->length = k;
+        }
+    }
+    ///ks_introsort_or_xs(overlap_list->length, overlap_list->list);
+}
+**/
+
+
 void ha_get_ug_candidates(ha_abuf_t *ab, int64_t rid, ma_utg_t *u, ma_utg_v *ua, overlap_region_alloc *overlap_list, Candidates_list *cl, double bw_thres, int max_n_chain, int keep_whole_chain, kvec_t_u8_warp* k_flag,
 kvec_t_u64_warp* chain_idx, void *ha_flt_tab, ha_pt_t *ha_idx, overlap_region* f_cigar, kvec_t_u64_warp* dbg_ct, double chain_match_rate)
 {
-	uint32_t i;
+    uint32_t i;
-	uint64_t k, l;
+    uint64_t k, l;
 
-	// prepare
+    // prepare
     clear_Candidates_list(cl);
     clear_overlap_region_alloc(overlap_list);
-	ab->mz.n = 0, ab->n_a = 0;
+    ab->mz.n = 0, ab->n_a = 0;
 
-	// get the list of anchors
+    // get the list of anchors
-	//should use the new version...
+    //should use the new version...
-	///ha_sketch_query(u->s, u->len, asm_opt.mz_win, asm_opt.k_mer_length, 0, !(asm_opt.flag & HA_F_NO_HPC), &ab->mz, ha_flt_tab, k_flag, dbg_ct);
+    ///ha_sketch_query(u->s, u->len, asm_opt.mz_win, asm_opt.k_mer_length, 0, !(asm_opt.flag & HA_F_NO_HPC), &ab->mz, ha_flt_tab, k_flag, dbg_ct);
-	// minimizer of queried read
+    // minimizer of queried read
-	if (ab->mz.m > ab->old_mz_m) {
+    if (ab->mz.m > ab->old_mz_m) {
-		ab->old_mz_m = ab->mz.m;
+        ab->old_mz_m = ab->mz.m;
-		REALLOC(ab->seed, ab->old_mz_m);
+        REALLOC(ab->seed, ab->old_mz_m);
-	}
+    }
-	for (i = 0, ab->n_a = 0; i < ab->mz.n; ++i) {
+    for (i = 0, ab->n_a = 0; i < ab->mz.n; ++i) {
-		int n;
+        int n;
-		ab->seed[i].a = ha_pt_get(ha_idx, ab->mz.a[i].x, &n);
+        ab->seed[i].a = ha_pt_get(ha_idx, ab->mz.a[i].x, &n);
-		ab->seed[i].n = n;
+        ab->seed[i].n = n;
-		ab->n_a += n;
+        ab->n_a += n;
-	}
+    }
-	if (ab->n_a > ab->m_a) {
+    if (ab->n_a > ab->m_a) {
-		ab->m_a = ab->n_a;
+        ab->m_a = ab->n_a;
-		kroundup64(ab->m_a);
+        kroundup64(ab->m_a);
-		REALLOC(ab->a, ab->m_a);
+        REALLOC(ab->a, ab->m_a);
-	}
+    }
-	for (i = 0, k = 0; i < ab->mz.n; ++i) {
+    for (i = 0, k = 0; i < ab->mz.n; ++i) {
-		int j;
+        int j;
-		///z is one of the minimizer
+        ///z is one of the minimizer
-		ha_mz1_t *z = &ab->mz.a[i];
+        ha_mz1_t *z = &ab->mz.a[i];
-		seed1_t *s = &ab->seed[i];
+        seed1_t *s = &ab->seed[i];
-		for (j = 0; j < s->n; ++j) {
+        for (j = 0; j < s->n; ++j) {
-			const ha_idxpos_t *y = &s->a[j];
+            const ha_idxpos_t *y = &s->a[j];
-			anchor1_t *an = &ab->a[k++];
+            anchor1_t *an = &ab->a[k++];
-			uint8_t rev = z->rev == y->rev? 0 : 1;
+            uint8_t rev = z->rev == y->rev? 0 : 1;
-			an->other_off = y->pos;
+            an->other_off = y->pos;
-			an->self_off = rev? u->len - 1 - (z->pos + 1 - z->span) : z->pos;
+            an->self_off = rev? u->len - 1 - (z->pos + 1 - z->span) : z->pos;
-			an->cnt = 1;
+            an->cnt = 1;
-			an->srt = (uint64_t)y->rid<<33 | (uint64_t)rev<<32 | an->other_off;
+            an->srt = (uint64_t)y->rid<<33 | (uint64_t)rev<<32 | an->other_off;
-		}
+        }
-	}
+    }
 
-	// sort anchors
+    // sort anchors
-	radix_sort_ha_an1(ab->a, ab->a + ab->n_a);
+    radix_sort_ha_an1(ab->a, ab->a + ab->n_a);
-	for (k = 1, l = 0; k <= ab->n_a; ++k) {
+    for (k = 1, l = 0; k <= ab->n_a; ++k) {
-		if (k == ab->n_a || ab->a[k].srt != ab->a[l].srt) {
+        if (k == ab->n_a || ab->a[k].srt != ab->a[l].srt) {
-			if (k - l > 1)
+            if (k - l > 1)
-				radix_sort_ha_an2(ab->a + l, ab->a + k);
+                radix_sort_ha_an2(ab->a + l, ab->a + k);
-			l = k;
+            l = k;
-		}
+        }
-	}
+    }
 
 
-	// copy over to _cl_
+    // copy over to _cl_
-	if (ab->m_a >= (uint64_t)cl->size) {
+    if (ab->m_a >= (uint64_t)cl->size) {
-		cl->size = ab->m_a;
+        cl->size = ab->m_a;
-		REALLOC(cl->list, cl->size);
+        REALLOC(cl->list, cl->size);
-	}
+    }
-	for (k = 0; k < ab->n_a; ++k) {
+    for (k = 0; k < ab->n_a; ++k) {
-		k_mer_hit *p = &cl->list[k];
+        k_mer_hit *p = &cl->list[k];
-		p->readID = ab->a[k].srt >> 33;
+        p->readID = ab->a[k].srt >> 33;
-		p->strand = ab->a[k].srt >> 32 & 1;
+        p->strand = ab->a[k].srt >> 32 & 1;
-		p->offset = ab->a[k].other_off;
+        p->offset = ab->a[k].other_off;
-		p->self_offset = ab->a[k].self_off;
+        p->self_offset = ab->a[k].self_off;
-		p->cnt = 1;
+        p->cnt = 1;
-	}
+    }
-	cl->length = ab->n_a;
+    cl->length = ab->n_a;
 
-	calculate_ug_chaining(cl, overlap_list, chain_idx, rid, ua, bw_thres, keep_whole_chain, f_cigar, ab->mz.n, chain_match_rate);
+    calculate_ug_chaining(cl, overlap_list, chain_idx, rid, ua, bw_thres, keep_whole_chain, f_cigar, ab->mz.n, chain_match_rate);
 
-	#if 0
+    #if 0
-	if (overlap_list->length > 0) {
+    if (overlap_list->length > 0) {
-		fprintf(stderr, "B\t%ld\t%ld\t%d\n", (long)rid, (long)overlap_list->length, rlen);
+        fprintf(stderr, "B\t%ld\t%ld\t%d\n", (long)rid, (long)overlap_list->length, rlen);
-		for (int i = 0; i < (int)overlap_list->length; ++i) {
+        for (int i = 0; i < (int)overlap_list->length; ++i) {
-			overlap_region *r = &overlap_list->list[i];
+            overlap_region *r = &overlap_list->list[i];
-			fprintf(stderr, "C\t%d\t%d\t%d\t%c\t%d\t%ld\t%d\t%d\t%c\t%d\t%d\n", (int)r->x_id, (int)r->x_pos_s, (int)r->x_pos_e, "+-"[r->x_pos_strand],
+            fprintf(stderr, "C\t%d\t%d\t%d\t%c\t%d\t%ld\t%d\t%d\t%c\t%d\t%d\n", (int)r->x_id, (int)r->x_pos_s, (int)r->x_pos_e, "+-"[r->x_pos_strand],
-					(int)r->y_id, (long)Get_READ_LENGTH(R_INF, r->y_id), (int)r->y_pos_s, (int)r->y_pos_e, "+-"[r->y_pos_strand], (int)r->shared_seed, ha_ov_type(r, rlen));
+                    (int)r->y_id, (long)Get_READ_LENGTH(R_INF, r->y_id), (int)r->y_pos_s, (int)r->y_pos_e, "+-"[r->y_pos_strand], (int)r->shared_seed, ha_ov_type(r, rlen));
-		}
+        }
-	}
+    }
-	#endif
+    #endif
 
-	if ((int)overlap_list->length > max_n_chain) {
+    if ((int)overlap_list->length > max_n_chain) {
-		int32_t w, n[4], s[4];
+        int32_t w, n[4], s[4];
-		n[0] = n[1] = n[2] = n[3] = 0, s[0] = s[1] = s[2] = s[3] = 0;
+        n[0] = n[1] = n[2] = n[3] = 0, s[0] = s[1] = s[2] = s[3] = 0;
-		ks_introsort_or_ss(overlap_list->length, overlap_list->list);
+        ks_introsort_or_ss(overlap_list->length, overlap_list->list);
-		for (i = 0; i < (uint32_t)overlap_list->length; ++i) {
+        for (i = 0; i < (uint32_t)overlap_list->length; ++i) {
-			const overlap_region *r = &overlap_list->list[i];
+            const overlap_region *r = &overlap_list->list[i];
-			w = ha_ov_type(r, u->len);
+            w = ha_ov_type(r, u->len);
-			++n[w];
+            ++n[w];
-			if ((int)n[w] == max_n_chain) s[w] = r->shared_seed;
+            if ((int)n[w] == max_n_chain) s[w] = r->shared_seed;
-		}
+        }
-		if (s[0] > 0 || s[1] > 0 || s[2] > 0 || s[3] > 0) {
+        if (s[0] > 0 || s[1] > 0 || s[2] > 0 || s[3] > 0) {
-			for (i = 0, k = 0; i < (uint32_t)overlap_list->length; ++i) {
+            for (i = 0, k = 0; i < (uint32_t)overlap_list->length; ++i) {
-				overlap_region *r = &overlap_list->list[i];
+                overlap_region *r = &overlap_list->list[i];
-				w = ha_ov_type(r, u->len);
+                w = ha_ov_type(r, u->len);
-				if (r->shared_seed >= s[w]) {
+                if (r->shared_seed >= s[w]) {
-					if ((uint32_t)k != i) {
+                    if ((uint32_t)k != i) {
-						overlap_region t;
+                        overlap_region t;
-						t = overlap_list->list[k];
+                        t = overlap_list->list[k];
-						overlap_list->list[k] = overlap_list->list[i];
+                        overlap_list->list[k] = overlap_list->list[i];
-						overlap_list->list[i] = t;
+                        overlap_list->list[i] = t;
-					}
+                    }
-					++k;
+                    ++k;
-				}
+                }
-			}
+            }
-			overlap_list->length = k;
+            overlap_list->length = k;
-		}
+        }
-	}
+    }
 
-	///ks_introsort_or_xs(overlap_list->length, overlap_list->list);
+    ///ks_introsort_or_xs(overlap_list->length, overlap_list->list);
 }
 
+
 void lable_matched_ovlp(overlap_region_alloc* overlap_list, ma_hit_t_alloc* paf)
 {
 	uint64_t j = 0, inner_j = 0;

horder.cpp

@@ -4016,7 +4016,7 @@ void update_sc_lay(sc_lay_t *sl, h_covs *b)
                 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);
+                memcpy(p->a, s->a + (pidx<<1), sizeof(*(p->a))*p->n);
                 pidx = cidx + 1;
             }
 
@@ -4026,7 +4026,7 @@ void update_sc_lay(sc_lay_t *sl, h_covs *b)
             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);
+            memcpy(p->a, s->a + (pidx<<1), sizeof(*(p->a))*p->n);
             free(s->a); s->n = s->m = 0;
             l = k;
         }

htab.cpp

@@ -504,7 +504,13 @@ void ha_pt_destroy(ha_pt_t *h)
 	if (h == 0) return;
 	for (i = 0; i < 1<<h->pre; ++i) {
 		yak_pt_destroy(h->h[i].h);
-		free(h->h[i].a);
+		if(h->h[i].a){
+			free(h->h[i].a); h->h[i].a = NULL;
+		}
+		if(h->h[i].al){
+			free(h->h[i].al); h->h[i].al = NULL;
+		}
+		
 	}
 	free(h->h); free(h);
 }
@@ -520,6 +526,17 @@ const ha_idxpos_t *ha_pt_get(const ha_pt_t *h, uint64_t hash, int *n)
 	return &g->a[kh_val(g->h, k)];
 }
 
+const ha_idxposl_t *ha_ptl_get(const ha_pt_t *h, uint64_t hash, int *n)
+{
+	khint_t k;
+	const ha_pt1_t *g = &h->h[hash & ((1ULL<<h->pre) - 1)];
+	*n = 0;
+	k = yak_pt_get(g->h, hash >> h->pre << YAK_COUNTER_BITS);
+	if (k == kh_end(g->h)) return 0;
+	*n = kh_key(g->h, k) & YAK_MAX_COUNT;
+	return &g->al[kh_val(g->h, k)];
+}
+
 const int ha_pt_cnt(const ha_pt_t *h, uint64_t hash)
 {
 	khint_t k;
@@ -703,7 +720,7 @@ static void *sf##_worker_count(void *data, int step, void *in) /** callback for
 					break;\
 			}\
 		} else if(p->us_in) {\
-			ma_utg_t *u; s->uq = p->us_in->h;\
+			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) {\
@@ -1052,7 +1069,7 @@ void *ha_ft_ug_gen(const hifiasm_opt_t *asm_opt, ma_utg_v *us, int is_HPC, int k
 	///HAF_COUNT_EXACT ---> no bf; HAF_COUNT_ALL ---> no minimizer
 	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, 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);
+	flt_tab = gen_hh(h, YAK_MAX_COUNT);
 	ha_ct_destroy(h);
 	return (void*)flt_tab;
 }

htab.h

@@ -62,6 +62,9 @@ typedef struct ha_pt_s ha_pt_t;
 struct ha_abuf_s;
 typedef struct ha_abuf_s ha_abuf_t;
 
+struct ha_abufl_s;
+typedef struct ha_abufl_s ha_abufl_t;
+
 extern const unsigned char seq_nt4_table[256];
 extern void *ha_flt_tab;
 extern ha_pt_t *ha_idx;
@@ -78,6 +81,7 @@ ha_pt_t *ha_pt_ug_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, ma_utg_
 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);
+const ha_idxposl_t *ha_ptl_get(const ha_pt_t *h, uint64_t hash, int *n);
 const int ha_pt_cnt(const ha_pt_t *h, uint64_t hash);
 
 int write_pt_index(void *flt_tab, ha_pt_t *ha_idx, All_reads* r, hifiasm_opt_t* opt, char* file_name);
@@ -89,6 +93,9 @@ int query_ct_index(void* ct_idx, uint64_t hash);
 ha_abuf_t *ha_abuf_init(void);
 void ha_abuf_destroy(ha_abuf_t *ab);
 uint64_t ha_abuf_mem(const ha_abuf_t *ab);
+ha_abufl_t *ha_abufl_init(void);
+void ha_abufl_destroy(ha_abufl_t *ab);
+uint64_t ha_abufl_mem(const ha_abufl_t *ab);
 
 double yak_cputime(void);
 void yak_reset_realtime(void);

inter.cpp

@@ -2,18 +2,45 @@
 #include <stdlib.h>
 #include <assert.h>
 #include <zlib.h>
+#include "kseq.h" // FASTA/Q parser
+#include "kthread.h"
 #include "inter.h"
 #include "Overlaps.h"
 #include "CommandLines.h"
 #include "htab.h"
+KSEQ_INIT(gzFile, gzread)
 
-void uidx_build(ma_ug_t *ug, int is_HPC, int k, int w, int hap_n)
+typedef struct {
+	int w, k, bw, max_gap, is_HPC, hap_n;
+} mg_idxopt_t;
+
+typedef struct { // global data structure for kt_pipeline()
+    const void *ha_flt_tab;
+    const ha_pt_t *ha_idx;
+    const mg_idxopt_t *opt;
+	kseq_t *ks;
+    int64_t chunk_size;
+    uint64_t n_thread;
+    uint64_t total_base;
+    uint64_t total_pair;
+} uldat_t;
+
+void init_mg_opt(mg_idxopt_t *opt, int is_HPC, int k, int w, int hap_n)
+{
+    opt->k = k; 
+    opt->w = w;
+    opt->hap_n = hap_n;
+    opt->is_HPC = is_HPC;
+    opt->bw = 2000;
+    opt->max_gap = 5000;
+}
+
+void uidx_build(ma_ug_t *ug, mg_idxopt_t *opt)
 {
     int flag = asm_opt.flag;
     asm_opt.flag |= HA_F_NO_HPC;
-    ug->u.h = hap_n;
+    ha_flt_tab = ha_ft_ug_gen(&asm_opt, &(ug->u), opt->is_HPC, opt->k, opt->w, 1, opt->hap_n*10);
-    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), opt->is_HPC, opt->k, opt->w, 1);
-    ha_idx = ha_pt_ug_gen(&asm_opt, ha_flt_tab, &(ug->u), is_HPC, k, w, hap_n);
     asm_opt.flag = flag;
 }
 
@@ -23,10 +50,110 @@ void uidx_destory()
     ha_pt_destroy(ha_idx);
 }
 
+static void *worker_ul_pipeline(void *data, int step, void *in) // callback for kt_pipeline()
+{
+    /**
+    uldat_t *p = (uldat_t*)data;
+    ///uint64_t total_base = 0, total_pair = 0;
+    if (step == 0) { // step 1: read a block of sequences
+        int ret1, ret2;
+        uint64_t l1, l2;
+		stepdat_t *s;
+		CALLOC(s, 1);
+        s->idx = p->idx; s->id = p->total_pair; s->t_ch = p->t_ch;
+        while (((ret1 = kseq_read(p->ks1)) >= 0)&&((ret2 = kseq_read(p->ks2)) >= 0)) 
+        {
+            if (p->ks1->seq.l < p->idx->k || p->ks2->seq.l < p->idx->k) continue;
+            if (s->n == s->m) {
+                s->m = s->m < 16? 16 : s->m + (s->n>>1);
+                REALLOC(s->len, s->m);
+                REALLOC(s->seq, s->m);
+            }
+
+            l1 = p->ks1->seq.l; l2 = p->ks2->seq.l;
+            MALLOC(s->seq[s->n], l1+l2);
+            s->sum_len += l1+l2;
+            memcpy(s->seq[s->n], p->ks1->seq.s, l1);
+            memcpy(s->seq[s->n]+l1, p->ks2->seq.s, l2);
+            s->len[s->n++] = (uint64_t)(l1<<32)|(uint64_t)l2;
+
+            if (s->sum_len >= p->chunk_size) break;            
+        }
+        p->total_pair += s->n;
+        if (s->sum_len == 0) free(s);
+		else return s;
+    }
+    else if (step == 1) { // step 2: alignment
+        stepdat_t *s = (stepdat_t*)in;
+        CALLOC(s->pos_buf, p->n_thread);
+        CALLOC(s->pos, s->n);
+        int i;
+        kt_for(p->n_thread, worker_for_alignment, s, s->n);
+        for (i = 0; i < s->n; ++i) {
+            free(s->seq[i]);
+            p->total_base += (s->len[i]>>32) + (uint32_t)s->len[i];
+        }
+        
+        free(s->seq); free(s->len);
+        for (i = 0; i < (int)p->n_thread; ++i) {
+            free(s->pos_buf[i].a.a);
+        }
+        free(s->pos_buf);
+		return s;
+    }
+    else if (step == 2) { // step 3: dump
+        stepdat_t *s = (stepdat_t*)in;
+        int i;
+        for (i = 0; i < s->n; ++i) {
+            // if(s->pos[i].a == NULL) continue;
+            // kv_push(pe_hit_hap, p->hits, s->pos[i]);
+            if(s->pos[i].s == (uint64_t)-1) continue;
+            kv_push(pe_hit, p->hits.a, s->pos[i]);
+        }
+        free(s->pos);
+        free(s);
+    }
+    **/
+    return 0;
+}
+
+int alignment_ul_pipeline(uldat_t* sl, const enzyme *fn)
+{
+    double index_time = yak_realtime();
+    int i;
+    for (i = 0; i < fn->n; i++){
+        gzFile fp;
+        if ((fp = gzopen(fn->a[i], "r")) == 0) return 0;
+        sl->ks = kseq_init(fp);
+        kt_pipeline(3, worker_ul_pipeline, sl, 3);
+        kseq_destroy(sl->ks);
+        gzclose(fp);
+    }
+    fprintf(stderr, "[M::%s::%.3f] ==> Qualification\n", __func__, yak_realtime()-index_time);
+    return 1;
+}
+
+int ul_align(mg_idxopt_t *opt, const enzyme *fn, void *ha_flt_tab, ha_pt_t *ha_idx)
+{
+    uldat_t sl; memset(&sl, 0, sizeof(sl));
+    sl.ha_flt_tab = ha_flt_tab;
+    sl.ha_idx = ha_idx;
+    sl.opt = opt;   
+    sl.chunk_size = 20000000;
+    sl.n_thread = asm_opt.thread_num;
+    alignment_ul_pipeline(&sl, fn);
+    return 1;
+}
 
 void ul_resolve(ma_ug_t *ug, int hap_n)
 {
-    uidx_build(ug, 1, 63, 63, hap_n);
+    mg_idxopt_t opt;
-    
+    init_mg_opt(&opt, 0, 19, 10, hap_n);
+    uidx_build(ug, &opt);
+    ul_align(&opt, asm_opt.ar, ha_flt_tab, ha_idx);
+
+
+
+
     uidx_destory();
 }

sketch.cpp

@@ -506,10 +506,10 @@ void sf##_ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is
                 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);\
+				filtered = (cnt >= 1<<28);\
-				if(is_unique){\
+                if(is_unique && (!filtered)) {\
-					filtered = (cnt < is_unique);\
+					filtered = (cnt == 0);\
-					cnt = (cnt == is_unique? 0:cnt);\
+					cnt = (cnt == 1? 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**/\