implemented ha_pt_get()

Assembly.cpp

@@ -12,7 +12,7 @@
 #include "htab.h"
 
 void *ha_flt_tab;
-void *ha_idx;
+ha_pt_t *ha_idx;
 
 All_reads R_INF;
 pthread_mutex_t statistics;
@@ -817,7 +817,7 @@ void Overlap_calculate_multipe_thr()
     }
 
 
-	ha_idx_destroy(ha_idx);
+	ha_pt_destroy(ha_idx);
 	ha_idx = 0;
 
     fprintf(stderr, "All overlaps have been calculated.\n");
@@ -1391,7 +1391,7 @@ void generate_overlaps(int last_round)
     asm_opt.roundID = asm_opt.number_of_round - last_round;
     fprintf(stderr, "Begin calculting final overlaps ...\n");
 
-	ha_idx = ha_gen_mzidx(&asm_opt, ha_flt_tab, asm_opt.roundID == 0? 0 : 1, &R_INF);
+	ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, asm_opt.roundID == 0? 0 : 1, &R_INF);
 
     pthread_t *_r_threads;
 
@@ -1412,7 +1412,7 @@ 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);
-	ha_idx_destroy(ha_idx);
+	ha_pt_destroy(ha_idx);
 	ha_idx = 0;
 
     fprintf(stderr, "Final overlaps have been calculated.\n");
@@ -1453,9 +1453,9 @@ void Correct_Reads(int last_round)
 
     fprintf(stderr, "Error correction: Start the %d-th round ...\n", asm_opt.roundID);
 
-	ha_idx = ha_gen_mzidx(&asm_opt, ha_flt_tab, asm_opt.roundID == 0? 0 : 1, &R_INF);
+	ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, asm_opt.roundID == 0? 0 : 1, &R_INF);
     Overlap_calculate_multipe_thr();
-	ha_idx_destroy(ha_idx);
+	ha_pt_destroy(ha_idx);
 
     fprintf(stderr, "Error correction: The %d-th round has been completed.\n", asm_opt.roundID);
 

hist.cpp

@@ -1,7 +1,7 @@
 #include <stdio.h>
 #include "htab.h"
 
-static void yak_hist_line(int c, int x, int exceed, int64_t cnt)
+static void ha_hist_line(int c, int x, int exceed, int64_t cnt)
 {
 	int j;
 	if (c >= 0) fprintf(stderr, "[M::%s] %5d: ", __func__, c);
@@ -11,7 +11,7 @@ static void yak_hist_line(int c, int x, int exceed, int64_t cnt)
 	fprintf(stderr, " %lld\n", (long long)cnt);
 }
 
-int yak_analyze_count(int n_cnt, const int64_t *cnt, int *peak_het)
+int ha_analyze_count(int n_cnt, const int64_t *cnt, int *peak_het)
 {
 	const int hist_max = 100;
 	int i, start, low_i, max_i, max2_i, max3_i;
@@ -40,7 +40,7 @@ int yak_analyze_count(int n_cnt, const int64_t *cnt, int *peak_het)
 		x = (int)((double)hist_max * cnt[i] / cnt[max_i] + .499);
 		if (x > hist_max) exceed = 1, x = hist_max; // may happen if cnt[2] is higher
 		if (i > max_i && x == 0) break;
-		yak_hist_line(i, x, exceed, cnt[i]);
+		ha_hist_line(i, x, exceed, cnt[i]);
 	}
 	{
 		int x, exceed = 0;
@@ -48,7 +48,7 @@ int yak_analyze_count(int n_cnt, const int64_t *cnt, int *peak_het)
 		for (; i < n_cnt; ++i) rest += cnt[i];
 		x = (int)((double)hist_max * rest / cnt[max_i] + .499);
 		if (x > hist_max) exceed = 1, x = hist_max;
-		yak_hist_line(-1, x, exceed, rest);
+		ha_hist_line(-1, x, exceed, rest);
 	}
 
 	// look for smaller peak on the low end

htab.cpp

@@ -281,11 +281,11 @@ typedef struct {
 	uint64_t *a;
 } ha_pt1_t;
 
-typedef struct {
+struct ha_pt_s {
 	int k, pre;
 	uint64_t tot, tot_pos;
 	ha_pt1_t *h;
-} ha_pt_t;
+};
 
 typedef struct {
 	const ha_ct_t *ct;
@@ -302,7 +302,7 @@ static void worker_pt_gen(void *data, long i, int tid) // callback for kt_for()
 		if (kh_exist(g, k)) {
 			int absent;
 			khint_t l;
-			l = yak_pt_put(b->h, kh_key(g, k) >> YAK_COUNTER_BITS << YAK_COUNTER_BITS, &absent);
+			l = yak_pt_put(b->h, kh_key(g, k) >> a->ct->pre << YAK_COUNTER_BITS, &absent);
 			kh_val(b->h, l) = b->n;
 			b->n += kh_key(g, k) & YAK_MAX_COUNT;
 		}
@@ -331,7 +331,7 @@ ha_pt_t *ha_pt_gen(ha_ct_t *ct, int n_thread)
 	return pt;
 }
 
-static int ha_pt_insert_list(ha_pt_t *h, int n, const ha_mz1_t *a)
+int ha_pt_insert_list(ha_pt_t *h, int n, const ha_mz1_t *a)
 {
 	int j, mask = (1<<h->pre) - 1, n_ins = 0;
 	ha_pt1_t *g;
@@ -367,12 +367,12 @@ static void worker_pt_sort(void *data, long i, int tid)
 	}
 }
 
-static void ha_pt_sort(ha_pt_t *h, int n_thread)
+void ha_pt_sort(ha_pt_t *h, int n_thread)
 {
 	kt_for(n_thread, worker_pt_sort, h, 1<<h->pre);
 }
 
-static void ha_pt_destroy(ha_pt_t *h)
+void ha_pt_destroy(ha_pt_t *h)
 {
 	int i;
 	if (h == 0) return;
@@ -383,9 +383,15 @@ static void ha_pt_destroy(ha_pt_t *h)
 	free(h->h); free(h);
 }
 
-void ha_idx_destroy(void *h)
+const uint64_t *ha_pt_get(const ha_pt_t *h, uint64_t hash, int *n)
 {
-	ha_pt_destroy((ha_pt_t*)h);
+	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->a[kh_val(g->h, k)];
 }
 
 /**********************************
@@ -729,7 +735,7 @@ static yak_ft_t *gen_hh(const ha_ct_t *h)
 		khint_t k;
 		for (k = 0; k < kh_end(ht); ++k) {
 			if (kh_exist(ht, k)) {
-				uint64_t y = kh_key(ht, k) >> YAK_COUNTER_BITS << h->pre | i;
+				uint64_t y = kh_key(ht, k) >> h->pre << YAK_COUNTER_BITS | i;
 				int absent;
 				yak_ft_put(hh, y, &absent);
 			}
@@ -755,7 +761,7 @@ void ha_ft_destroy(void *h)
  * High-level interfaces *
  *************************/
 
-void *ha_gen_flt_tab(const hifiasm_opt_t *asm_opt, All_reads *rs)
+void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs)
 {
 	yak_ft_t *flt_tab;
 	int64_t cnt[YAK_N_COUNTS];
@@ -763,7 +769,7 @@ void *ha_gen_flt_tab(const hifiasm_opt_t *asm_opt, All_reads *rs)
 	ha_ct_t *h;
 	h = ha_count(asm_opt, HAF_COUNT_ALL|HAF_RS_WRITE_LEN, NULL, NULL, rs);
 	ha_ct_hist(h, cnt, asm_opt->thread_num);
-	peak_hom = yak_analyze_count(YAK_N_COUNTS, cnt, &peak_het);
+	peak_hom = ha_analyze_count(YAK_N_COUNTS, cnt, &peak_het);
 	if (peak_hom > 0) fprintf(stderr, "[M::%s] peak_hom: %d; peak_het: %d\n", __func__, peak_hom, peak_het);
 	cutoff = (int)(peak_hom * asm_opt->high_factor);
 	if (cutoff > YAK_MAX_COUNT - 1) cutoff = YAK_MAX_COUNT - 1;
@@ -775,7 +781,7 @@ void *ha_gen_flt_tab(const hifiasm_opt_t *asm_opt, All_reads *rs)
 	return (void*)flt_tab;
 }
 
-void *ha_gen_mzidx(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_from_store, All_reads *rs)
+ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_from_store, All_reads *rs)
 {
 	int64_t cnt[YAK_N_COUNTS], tot_cnt;
 	int peak_hom, peak_het, i, extra_flag = read_from_store? HAF_RS_READ : HAF_RS_WRITE_SEQ;
@@ -786,7 +792,7 @@ void *ha_gen_mzidx(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_f
 			yak_realtime(), yak_cputime() / yak_realtime(), (long)ct->tot);
 	ha_ct_hist(ct, cnt, asm_opt->thread_num);
 	fprintf(stderr, "[M::%s] count[%d] = %ld (for sanity check)\n", __func__, YAK_MAX_COUNT, (long)cnt[YAK_MAX_COUNT]);
-	peak_hom = yak_analyze_count(YAK_N_COUNTS, cnt, &peak_het);
+	peak_hom = ha_analyze_count(YAK_N_COUNTS, cnt, &peak_het);
 	if (peak_hom > 0) fprintf(stderr, "[M::%s] peak_hom: %d; peak_het: %d\n", __func__, peak_hom, peak_het);
 	ha_ct_shrink(ct, 2, YAK_MAX_COUNT - 1, asm_opt->thread_num);
 	for (i = 2, tot_cnt = 0; i <= YAK_MAX_COUNT - 1; ++i) tot_cnt += cnt[i] * i;

htab.h

@@ -12,23 +12,28 @@ typedef struct {
 
 typedef struct { uint32_t n, m; ha_mz1_t *a; } ha_mz1_v;
 
+struct ha_pt_s;
+typedef struct ha_pt_s ha_pt_t;
+
 extern const unsigned char seq_nt4_table[256];
 
-void *ha_gen_flt_tab(const hifiasm_opt_t *asm_opt, All_reads *rs);
+void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs);
-void *ha_gen_mzidx(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_from_store, All_reads *rs);
-void trio_partition(void);
-
 int ha_ft_isflt(const void *hh, uint64_t y);
 void ha_ft_destroy(void *h);
-void ha_idx_destroy(void *h);
+
+ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_from_store, All_reads *rs);
+void ha_pt_destroy(ha_pt_t *h);
+const uint64_t *ha_pt_get(const ha_pt_t *h, uint64_t hash, int *n);
 
 double yak_cputime(void);
 void yak_reset_realtime(void);
 double yak_realtime(void);
 long yak_peakrss(void);
 
+void trio_partition(void);
+
 void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mz1_v *p, const void *hf);
-int yak_analyze_count(int n_cnt, const int64_t *cnt, int *peak_het);
+int ha_analyze_count(int n_cnt, const int64_t *cnt, int *peak_het);
 
 static inline uint64_t yak_hash64(uint64_t key, uint64_t mask) // invertible integer hash function
 {

main.cpp

@@ -10,16 +10,16 @@ int main(int argc, char *argv[])
 {
 	int i;
 	void *flt_tab;
-	void *idx;
+	ha_pt_t *idx;
 
     init_opt(&asm_opt);
 
     if (!CommandLine_process(argc, argv, &asm_opt)) return 1;
 
 	yak_reset_realtime();
-	flt_tab = ha_gen_flt_tab(&asm_opt, &R_INF);
+	flt_tab = ha_ft_gen(&asm_opt, &R_INF);
-	idx = ha_gen_mzidx(&asm_opt, flt_tab, 0, &R_INF);
+	idx = ha_pt_gen(&asm_opt, flt_tab, 0, &R_INF);
-	ha_idx_destroy(idx);
+	ha_pt_destroy(idx);
 	ha_ft_destroy(flt_tab);
 	if (0) {
 	    Correct_Reads(asm_opt.number_of_round);