diff --git a/Assembly.cpp b/Assembly.cpp index d4f0921..6aa1cbd 100644 --- a/Assembly.cpp +++ b/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); diff --git a/hist.cpp b/hist.cpp index 5021c4b..9a1e31e 100644 --- a/hist.cpp +++ b/hist.cpp @@ -1,7 +1,7 @@ #include #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 diff --git a/htab.cpp b/htab.cpp index 421fbda..6b696b9 100644 --- a/htab.cpp +++ b/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<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<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<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; diff --git a/htab.h b/htab.h index 4f68d96..d7232cf 100644 --- a/htab.h +++ b/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_gen_mzidx(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_from_store, All_reads *rs); -void trio_partition(void); - +void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs); 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 { diff --git a/main.cpp b/main.cpp index e257b2b..842c0f9 100644 --- a/main.cpp +++ b/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); - idx = ha_gen_mzidx(&asm_opt, flt_tab, 0, &R_INF); - ha_idx_destroy(idx); + flt_tab = ha_ft_gen(&asm_opt, &R_INF); + idx = ha_pt_gen(&asm_opt, flt_tab, 0, &R_INF); + ha_pt_destroy(idx); ha_ft_destroy(flt_tab); if (0) { Correct_Reads(asm_opt.number_of_round);