From 815d709fb9868ac5f5f526cbf947d7bb30422745 Mon Sep 17 00:00:00 2001 From: chhylp123 Date: Wed, 1 Apr 2026 22:28:54 -0400 Subject: [PATCH] fall back to avx2 --- Assembly.cpp | 6 +- CommandLines.cpp | 16 ++- CommandLines.h | 8 +- Correct.cpp | 310 ++++++++++++++++++++++++++++++++++------ Levenshtein_avx2.cpp | 162 +++++++++++++++++++++ Levenshtein_avx512.cpp | 276 ++++++++++++++++++++++++++++++++++++ Levenshtein_distance.h | 271 +---------------------------------- Makefile | 15 +- Overlaps.cpp | 17 ++- Process_Read.cpp | 49 +++++++ Process_Read.h | 6 + ecovlp.cpp | 314 +++++++++++++++++++++++++++++++++-------- ecovlp.h | 3 + htab.cpp | 29 +++- htab.h | 2 +- main.cpp | 28 +++- 16 files changed, 1122 insertions(+), 390 deletions(-) create mode 100644 Levenshtein_avx2.cpp create mode 100644 Levenshtein_avx512.cpp diff --git a/Assembly.cpp b/Assembly.cpp index 99bac20..02165dd 100644 --- a/Assembly.cpp +++ b/Assembly.cpp @@ -2135,7 +2135,9 @@ int ha_assemble(void) fprintf(stderr, "[M::%s::%.3f*%.2f@%.3fGB] ==> found overlaps for the final round\n", __func__, yak_realtime(), yak_cpu_usage(), yak_peakrss_in_gb()); // fprintf(stderr, "\n[M::%s::%.3f*%.2f@%.3fGB] ==> found overlaps for the final round\n", __func__, yak_realtime(), yak_cpu_usage(), yak_peakrss_in_gb()); // ha_print_ovlp_stat(R_INF.paf, R_INF.reverse_paf, R_INF.total_reads); - if(!(asm_opt.write_pos_idx)) ha_ft_destroy(ha_flt_tab); + if(!(asm_opt.write_pos_idx)) { + ha_ft_destroy(ha_flt_tab); ha_flt_tab = NULL; + } if (asm_opt.flag & HA_F_WRITE_PAF) Output_PAF(); ha_triobin(&asm_opt); @@ -2147,7 +2149,7 @@ int ha_assemble(void) 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, asm_opt.gap_fuzz, asm_opt.min_drop_rate, asm_opt.max_drop_rate, asm_opt.output_file_name, asm_opt.large_pop_bubble_size, 0, !ovlp_loaded); - destory_All_reads(&R_INF); + destory_All_reads(&R_INF); if(asm_opt.dbg_bam) destroy_cc_v(&scb); return 0; } diff --git a/CommandLines.cpp b/CommandLines.cpp index f015622..c50258d 100644 --- a/CommandLines.cpp +++ b/CommandLines.cpp @@ -92,6 +92,7 @@ static ko_longopt_t long_options[] = { { "dbg-in2", ko_required_argument, 374}, { "ec-only", ko_no_argument, 375}, { "hyb-syn", ko_required_argument, 376}, + { "simd-m", ko_required_argument, 377}, // { "path-round", ko_required_argument, 348}, { 0, 0, 0 } }; @@ -136,6 +137,8 @@ void Print_H(hifiasm_opt_t* asm_opt) fprintf(stderr, " --chn-occ INT\n"); fprintf(stderr, " discard overlaps supported by chn_occ); fprintf(stderr, " --ec-only error correction only; disable overlapping and assembly\n"); + fprintf(stderr, " --simd-m use SIMD acceleration when supported: AVX-512 (2), AVX2 (1), or non-SIMD (0)\n"); + fprintf(stderr, " Assembly:\n"); fprintf(stderr, " -a INT round of assembly cleaning [%d]\n", asm_opt->clean_round); fprintf(stderr, " -m INT pop bubbles of large_pop_bubble_size); @@ -247,7 +250,8 @@ void Print_H(hifiasm_opt_t* asm_opt) fprintf(stderr, " --ont assemble ONT Simplex reads in fastq format\n"); // fprintf(stderr, " --sc-n consider base qual value for assembly\n"); fprintf(stderr, " --chem-c INT\n"); - fprintf(stderr, " detect chimeric reads with <=INT other reads support [%lu]\n", asm_opt->chemical_cov); + // fprintf(stderr, " detect chimeric reads with <=INT other reads support [%lu]\n", asm_opt->chemical_cov); + fprintf(stderr, " detect chimeric reads with <=INT other reads support [auto]\n"); fprintf(stderr, " --chem-f INT\n"); fprintf(stderr, " length of flanking regions for chimeric read detection [%lu]\n", asm_opt->chemical_flank); fprintf(stderr, " --rl-cut INT\n"); @@ -396,7 +400,7 @@ void init_opt(hifiasm_opt_t* asm_opt) asm_opt->is_ont = 0; asm_opt->is_sc = 0; - asm_opt->chemical_cov = 1; + asm_opt->chemical_cov = -1/**1**/; asm_opt->chemical_flank = 256; asm_opt->ul_mod = 0; @@ -416,7 +420,7 @@ void init_opt(hifiasm_opt_t* asm_opt) asm_opt->hf_cutoff = -1; - asm_opt->write_pos_idx = 1; + asm_opt->write_pos_idx = 0/**1**/; asm_opt->hom_cov_0 = -1; asm_opt->het_cov_0 = -1; @@ -433,6 +437,10 @@ void init_opt(hifiasm_opt_t* asm_opt) asm_opt->ec_only = 0; asm_opt->hyb_syn = 1; asm_opt->step_rd = -1/**128**/; + + asm_opt->dbg_bam = 0/**1**/; + + asm_opt->simd_mm = -1; } void destory_enzyme(enzyme* f) @@ -1106,6 +1114,8 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt) asm_opt->ec_only = 1; } else if (c == 376) { asm_opt->hyb_syn = atoi(opt.arg); + } else if (c == 377) { + asm_opt->simd_mm = atoi(opt.arg); } else if (c == 'l') { ///0: disable purge_dup; 1: purge containment; 2: purge overlap asm_opt->purge_level_primary = asm_opt->purge_level_trio = atoi(opt.arg); } diff --git a/CommandLines.h b/CommandLines.h index bb7eb72..7d2e16f 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -5,7 +5,7 @@ #include #include -#define HA_VERSION "0.25.1-r878" +#define HA_VERSION "0.25.1-r891" #define VERBOSE 0 @@ -168,7 +168,7 @@ typedef struct { uint64_t is_ont; uint64_t is_sc; - uint64_t chemical_cov; + int64_t chemical_cov; uint64_t chemical_flank; int64_t rl_cut; @@ -199,6 +199,10 @@ typedef struct { int32_t hyb_syn; int64_t step_rd; + + uint8_t dbg_bam; + + int8_t simd_mm; } hifiasm_opt_t; extern hifiasm_opt_t asm_opt; diff --git a/Correct.cpp b/Correct.cpp index e76fe68..cde2e0d 100644 --- a/Correct.cpp +++ b/Correct.cpp @@ -35114,7 +35114,7 @@ uint64_t gen_hc_r_alin_adp_mmp_0(overlap_region_alloc* ol, Candidates_list *cl, ///[ws, we) -uint8_t hc_aln_simd(overlap_region* ol, uint64_t *ffa, uint32_t *ia, uint64_t in, int64_t ws, int64_t we, int64_t wl, double e_rate, All_reads *rref, char *qu, UC_Read* tu, +uint8_t hc_aln_simd_avx512(overlap_region* ol, uint64_t *ffa, uint32_t *ia, uint64_t in, int64_t ws, int64_t we, int64_t wl, double e_rate, All_reads *rref, char *qu, UC_Read* tu, uint64_t *fi, int32_t *baux_beg, int32_t *baux_end, int32_t *bt_s, int32_t *bt_pri_l, bit_extz_t *exz, double ovlp_cut, int64_t force_aln, uint64_t *tot_b) { uint64_t k, mm_k = 0, *msk = NULL; uint8_t rr = 0, fr = 0; overlap_region *z; int64_t bl = we - ws, bthre, thre, q[2], wqs, wqe, wql, wts, wte, aux_beg, aux_end, aln_l, t_tot_l, t_pri_l, mbl = 0; @@ -35311,6 +35311,156 @@ uint8_t hc_aln_simd(overlap_region* ol, uint64_t *ffa, uint32_t *ia, uint64_t in return ((rr || fr)?1:0); } +///[ws, we) +uint8_t hc_aln_simd_avx2(overlap_region* ol, uint64_t *ffa, uint32_t *ia, uint64_t in, int64_t ws, int64_t we, int64_t wl, double e_rate, All_reads *rref, char *qu, UC_Read* tu, + uint64_t *fi, int32_t *baux_beg, int32_t *baux_end, int32_t *bt_s, int32_t *bt_pri_l, bit_extz_t *exz, double ovlp_cut, int64_t force_aln, uint64_t *tot_b) +{ + uint64_t k, mm_k = 0, *msk = NULL; uint8_t rr = 0, fr = 0; overlap_region *z; int64_t bl = we - ws, bthre, thre, q[2], wqs, wqe, wql, wts, wte, aux_beg, aux_end, aln_l, t_tot_l, t_pri_l, mbl = 0; + char *qstr = NULL, *tstr[AVX_GS2]; int64_t ez_er[AVX_GS2], ez_pe[AVX_GS2]; + + mbl = wl+(THRESHOLD_MAX_SIZE<<1)+1; resize_UC_Read(tu, mbl<<1); + bthre = bl*e_rate; bthre = Adjust_Threshold(bthre, bl); + if(bthre > THRESHOLD_MAX_SIZE) bthre = THRESHOLD_MAX_SIZE; + + for (k = mbl = mm_k = 0; k < in; k++) { + z = &(ol[(uint32_t)ffa[ia[k]]]); + q[0] = z->x_pos_s; q[1] = z->x_pos_e + 1; + if(q[1] <= we) rr = 1; + wqs = MAX(q[0], ws); wqe = MIN(q[1], we); + if(wqe <= wqs) continue; + + wql = wqe - wqs; aux_beg = aux_end = 0; + thre = wql*e_rate; thre = Adjust_Threshold(thre, wql); + if(thre > THRESHOLD_MAX_SIZE) thre = THRESHOLD_MAX_SIZE; + + wts = (wqs - z->x_pos_s) + z->y_pos_s; + wts += y_start_offset(wqs, &(z->f_cigar)); + + aln_l = wql + (thre<<1); t_tot_l = Get_READ_LENGTH((*rref), z->y_id); + if(init_waln(thre, wts, t_tot_l, aln_l, &aux_beg, &aux_end, &wts, &t_pri_l)) { + if(wql == bl) { + fi[mm_k] = ia[k]; baux_beg[mm_k] = aux_beg; baux_end[mm_k] = aux_end; + bt_s[mm_k] = wts; bt_pri_l[mm_k] = t_pri_l; mm_k++; mbl += t_pri_l; + + if(mm_k == AVX_GS2) { + qstr = qu + wqs; + if(mbl > tu->size) resize_UC_Read(tu, mbl); + + // fprintf(stderr, "[M::%s]\tw::[%ld,%ld)\twq::[%ld,%ld)\tmbl::%ld\tbt_pri_l[0]::%ld\tbt_pri_l[1]::%ld\tbt_pri_l[2]::%ld\tbt_pri_l[3]::%ld\tbt_pri_l[4]::%ld\tbt_pri_l[5]::%ld\tbt_pri_l[6]::%ld\tbt_pri_l::%ld\n", __func__, ws, we, wqs, wqe, mbl, + // bt_pri_l[0], bt_pri_l[1], bt_pri_l[2], bt_pri_l[3], bt_pri_l[4], bt_pri_l[5], bt_pri_l[6], bt_pri_l[7]); + + + tstr[0] = tu->seq; recover_UC_Read_sub_region(tstr[0], bt_s[0], bt_pri_l[0], ol[(uint32_t)ffa[fi[0]]].y_pos_strand, rref, ol[(uint32_t)ffa[fi[0]]].y_id); + tstr[1] = tstr[0] + bt_pri_l[0]; recover_UC_Read_sub_region(tstr[1], bt_s[1], bt_pri_l[1], ol[(uint32_t)ffa[fi[1]]].y_pos_strand, rref, ol[(uint32_t)ffa[fi[1]]].y_id); + tstr[2] = tstr[1] + bt_pri_l[1]; recover_UC_Read_sub_region(tstr[2], bt_s[2], bt_pri_l[2], ol[(uint32_t)ffa[fi[2]]].y_pos_strand, rref, ol[(uint32_t)ffa[fi[2]]].y_id); + tstr[3] = tstr[2] + bt_pri_l[2]; recover_UC_Read_sub_region(tstr[3], bt_s[3], bt_pri_l[3], ol[(uint32_t)ffa[fi[3]]].y_pos_strand, rref, ol[(uint32_t)ffa[fi[3]]].y_id); + + ///same thre due to the same wql + ed_band_cal_semi_64_w_absent_diag_avx4(tstr, bt_pri_l, qstr, wql, thre, baux_beg, ez_er, ez_pe); + + // ed_band_cal_semi_64_w_absent_diag(tstr[0], bt_pri_l[0], qstr, wql, thre, baux_beg[0], exz); ez_er[0] = exz->err; ez_pe[0] = exz->pe; + // ed_band_cal_semi_64_w_absent_diag(tstr[1], bt_pri_l[1], qstr, wql, thre, baux_beg[1], exz); ez_er[1] = exz->err; ez_pe[1] = exz->pe; + // ed_band_cal_semi_64_w_absent_diag(tstr[2], bt_pri_l[2], qstr, wql, thre, baux_beg[2], exz); ez_er[2] = exz->err; ez_pe[2] = exz->pe; + // ed_band_cal_semi_64_w_absent_diag(tstr[3], bt_pri_l[3], qstr, wql, thre, baux_beg[3], exz); ez_er[3] = exz->err; ez_pe[3] = exz->pe; + + exz->ps = exz->pe = -1; exz->ts = 0; exz->te = wql-1; + + init_base_ed((*exz), thre, bt_pri_l[0], wql); exz->err = ez_er[0]; exz->pe = ez_pe[0]; z = &(ol[(uint32_t)ffa[fi[0]]]); + wts = bt_s[0]; wte = bt_s[0] + exz->pe; t_tot_l = Get_READ_LENGTH((*rref), z->y_id); + aux_beg = baux_beg[0]; aux_end = baux_end[0]; msk = &(ffa[fi[0]]); + // if(z->x_id == 1 && z->y_id == 24) { + // fprintf(stderr, "-0-[M::%s]\tw::[%ld,%ld)\twq::[%ld,%ld)\tzql::%u\talign_length::%u\tthre::%ld\taln::%u\n", __func__, ws, we, wqs, wqe, z->x_pos_e+1-z->x_pos_s, z->align_length, thre, (is_align(*exz))); + // } + if ((is_align(*exz)) && (!push_hc_wlst_exz(NULL, NULL, rref, z, qu, tu->seq, exz, THRESHOLD_MAX_SIZE, wqs, wqe-1, wts, wte, t_tot_l, aux_beg, aux_end, e_rate, wl, ovlp_cut, force_aln, tot_b, 0))) { + (*msk) |= (((uint64_t)UINT32_MAX)<<32); fr = 1; ///disable this chain + } + + init_base_ed((*exz), thre, bt_pri_l[1], wql); exz->err = ez_er[1]; exz->pe = ez_pe[1]; z = &(ol[(uint32_t)ffa[fi[1]]]); + wts = bt_s[1]; wte = bt_s[1] + exz->pe; t_tot_l = Get_READ_LENGTH((*rref), z->y_id); + aux_beg = baux_beg[1]; aux_end = baux_end[1]; msk = &(ffa[fi[1]]); + // if(z->x_id == 1 && z->y_id == 24) { + // fprintf(stderr, "-1-[M::%s]\tw::[%ld,%ld)\twq::[%ld,%ld)\tzql::%u\talign_length::%u\tthre::%ld\taln::%u\n", __func__, ws, we, wqs, wqe, z->x_pos_e+1-z->x_pos_s, z->align_length, thre, (is_align(*exz))); + // } + if ((is_align(*exz)) && (!push_hc_wlst_exz(NULL, NULL, rref, z, qu, tu->seq, exz, THRESHOLD_MAX_SIZE, wqs, wqe-1, wts, wte, t_tot_l, aux_beg, aux_end, e_rate, wl, ovlp_cut, force_aln, tot_b, 0))) { + (*msk) |= (((uint64_t)UINT32_MAX)<<32); fr = 1; ///disable this chain + } + + init_base_ed((*exz), thre, bt_pri_l[2], wql); exz->err = ez_er[2]; exz->pe = ez_pe[2]; z = &(ol[(uint32_t)ffa[fi[2]]]); + wts = bt_s[2]; wte = bt_s[2] + exz->pe; t_tot_l = Get_READ_LENGTH((*rref), z->y_id); + aux_beg = baux_beg[2]; aux_end = baux_end[2]; msk = &(ffa[fi[2]]); + // if(z->x_id == 1 && z->y_id == 24) { + // fprintf(stderr, "-2-[M::%s]\tw::[%ld,%ld)\twq::[%ld,%ld)\tzql::%u\talign_length::%u\tthre::%ld\taln::%u\n", __func__, ws, we, wqs, wqe, z->x_pos_e+1-z->x_pos_s, z->align_length, thre, (is_align(*exz))); + // } + if ((is_align(*exz)) && (!push_hc_wlst_exz(NULL, NULL, rref, z, qu, tu->seq, exz, THRESHOLD_MAX_SIZE, wqs, wqe-1, wts, wte, t_tot_l, aux_beg, aux_end, e_rate, wl, ovlp_cut, force_aln, tot_b, 0))) { + (*msk) |= (((uint64_t)UINT32_MAX)<<32); fr = 1; ///disable this chain + } + + init_base_ed((*exz), thre, bt_pri_l[3], wql); exz->err = ez_er[3]; exz->pe = ez_pe[3]; z = &(ol[(uint32_t)ffa[fi[3]]]); + wts = bt_s[3]; wte = bt_s[3] + exz->pe; t_tot_l = Get_READ_LENGTH((*rref), z->y_id); + aux_beg = baux_beg[3]; aux_end = baux_end[3]; msk = &(ffa[fi[3]]); + // if(z->x_id == 1 && z->y_id == 24) { + // fprintf(stderr, "-3-[M::%s]\tw::[%ld,%ld)\twq::[%ld,%ld)\tzql::%u\talign_length::%u\tthre::%ld\taln::%u\n", __func__, ws, we, wqs, wqe, z->x_pos_e+1-z->x_pos_s, z->align_length, thre, (is_align(*exz))); + // } + if ((is_align(*exz)) && (!push_hc_wlst_exz(NULL, NULL, rref, z, qu, tu->seq, exz, THRESHOLD_MAX_SIZE, wqs, wqe-1, wts, wte, t_tot_l, aux_beg, aux_end, e_rate, wl, ovlp_cut, force_aln, tot_b, 0))) { + (*msk) |= (((uint64_t)UINT32_MAX)<<32); fr = 1; ///disable this chain + } + + mm_k = mbl = 0; + } + + } else { + qstr = qu + wqs; + recover_UC_Read_sub_region(tu->seq, wts, t_pri_l, z->y_pos_strand, rref, z->y_id); + ed_band_cal_semi_64_w_absent_diag(tu->seq, t_pri_l, qstr, wql, thre, aux_beg, exz); + wte = wts + exz->pe; msk = &(ffa[ia[k]]); + // if(z->x_id == 1 && z->y_id == 24) { + // fprintf(stderr, "-s-[M::%s]\tw::[%ld,%ld)\twq::[%ld,%ld)\tzql::%u\talign_length::%u\tthre::%ld\taln::%u\twql::%ld\te_rate::%f\n", __func__, ws, we, wqs, wqe, z->x_pos_e+1-z->x_pos_s, z->align_length, + // thre, (is_align(*exz)), wql, e_rate); + // } + if ((is_align(*exz)) && (!push_hc_wlst_exz(NULL, NULL, rref, z, qu, tu->seq, exz, THRESHOLD_MAX_SIZE, wqs, wqe-1, wts, wte, t_tot_l, aux_beg, aux_end, e_rate, wl, ovlp_cut, force_aln, tot_b, 0))) { + (*msk) |= (((uint64_t)UINT32_MAX)<<32); fr = 1; ///disable this chain + } + } + } + } + + if(mm_k > 0) { + wqs = ws; wqe = we; wql = bl; thre = bthre; + qstr = qu + wqs; + if(mbl > tu->size) resize_UC_Read(tu, mbl); + + for (k = mbl = 0; k < mm_k; k++) { + tstr[k] = tu->seq + mbl; + recover_UC_Read_sub_region(tstr[k], bt_s[k], bt_pri_l[k], ol[(uint32_t)ffa[fi[k]]].y_pos_strand, rref, ol[(uint32_t)ffa[fi[k]]].y_id); + mbl += bt_pri_l[k]; + } + for (; k < AVX_GS2; k++) { + tstr[k] = NULL; bt_pri_l[k] = 0; baux_beg[k] = 0; + } + + + if(mm_k > 1) { + ed_band_cal_semi_64_w_absent_diag_avx4(tstr, bt_pri_l, qstr, wql, thre, baux_beg, ez_er, ez_pe); + } else { + ed_band_cal_semi_64_w_absent_diag(tstr[0], bt_pri_l[0], qstr, wql, thre, baux_beg[0], exz); ez_er[0] = exz->err; ez_pe[0] = exz->pe; + } + + for (k = 0; k < mm_k; k++) { + init_base_ed((*exz), thre, bt_pri_l[k], wql); exz->err = ez_er[k]; exz->pe = ez_pe[k]; z = &(ol[(uint32_t)ffa[fi[k]]]); + wts = bt_s[k]; wte = bt_s[k] + exz->pe; t_tot_l = Get_READ_LENGTH((*rref), z->y_id); + aux_beg = baux_beg[k]; aux_end = baux_end[k]; msk = &(ffa[fi[k]]); + // if(z->x_id == 1 && z->y_id == 24) { + // fprintf(stderr, "-mm_k::%lu-[M::%s]\tw::[%ld,%ld)\twq::[%ld,%ld)\tzql::%u\talign_length::%u\tthre::%ld\taln::%u\n", mm_k, __func__, ws, we, wqs, wqe, z->x_pos_e+1-z->x_pos_s, z->align_length, thre, (is_align(*exz))); + // } + if ((is_align(*exz)) && (!push_hc_wlst_exz(NULL, NULL, rref, z, qu, tu->seq, exz, THRESHOLD_MAX_SIZE, wqs, wqe-1, wts, wte, t_tot_l, aux_beg, aux_end, e_rate, wl, ovlp_cut, force_aln, tot_b, 0))) { + (*msk) |= (((uint64_t)UINT32_MAX)<<32); fr = 1; ///disable this chain + } + } + } + + return ((rr || fr)?1:0); +} + uint64_t gen_hc_r_alin_adp_mmp_1(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t wl0, int64_t rid, int64_t khit, int64_t move_gap, asg16_v* buf, uint8_t chem_drop, double align_gap_rate, int64_t align_gap_max, uint64_t sec_aln_win, uint64_t sec_aln_cov, double sec_aln_err_rate, double sec_aln_max, asg64_v *sp, uint64_t ocw, uint8_t *hpf, asg32_v *v32, asg64_v *bp, uint64_t max_n_chain, uint64_t max_n_chain_f, uint64_t chain_cutoff, uint64_t ave_cov_min, uint8_t set_match) { @@ -35353,37 +35503,73 @@ uint64_t gen_hc_r_alin_adp_mmp_1(overlap_region_alloc* ol, Candidates_list *cl, v32->n = (ol->length<<1); in0 = v32->n; i = 0; s = 0; e = nwl; e = ((e<=ql)?e:ql); rr = 0; - for (; s < ql; ) {///[s, e) - if(rr) { - for (k = zk = in0; k < v32->n; k++) { - if((wsrt[v32->a[k]]>>32) == UINT32_MAX) continue;///passed - z = &(ol->list[(uint32_t)wsrt[v32->a[k]]]); + + if(asm_opt.simd_mm == 2) { + for (; s < ql; ) {///[s, e) + if(rr) { + for (k = zk = in0; k < v32->n; k++) { + if((wsrt[v32->a[k]]>>32) == UINT32_MAX) continue;///passed + z = &(ol->list[(uint32_t)wsrt[v32->a[k]]]); + q[0] = z->x_pos_s; + q[1] = z->x_pos_e + 1; + ///[s, e) && [q[0], q[1]) + os = MAX(q[0], s); oe = MIN(q[1], e); + if(oe > os) v32->a[zk++] = v32->a[k]; + } + v32->n = zk; + } + + for (; i < wsrt_n; ++i) { + z = &(ol->list[(uint32_t)wsrt[i]]); q[0] = z->x_pos_s; q[1] = z->x_pos_e + 1; - ///[s, e) && [q[0], q[1]) - os = MAX(q[0], s); oe = MIN(q[1], e); - if(oe > os) v32->a[zk++] = v32->a[k]; + if(q[0] >= e) break; + if(z->is_match == 0) { + os = MAX(q[0], s); oe = MIN(q[1], e); + if(oe > os) kv_push(uint32_t, *v32, i); + } } - v32->n = zk; - } - for (; i < wsrt_n; ++i) { - z = &(ol->list[(uint32_t)wsrt[i]]); - q[0] = z->x_pos_s; - q[1] = z->x_pos_e + 1; - if(q[0] >= e) break; - if(z->is_match == 0) { - os = MAX(q[0], s); oe = MIN(q[1], e); - if(oe > os) kv_push(uint32_t, *v32, i); + rr = hc_aln_simd_avx512(ol->list, wsrt, v32->a + in0, v32->n - in0, s, e, nwl, err, rref, qu->seq, tu, fi, aux_beg, aux_end, t_s, t_pri_l, exz, OVERLAP_THRESHOLD_HIFI_FILTER, 0, &tot_b); + // hc_est_robust_rr(ol->list, ql, ix->a + srt_n, ix->n - srt_n, ix->a + ix->n, ix->a + ix->n + ix->n - srt_n, s, e, 2.0, c_idx->a, &est_bd, &est_e, min_dp, 0); + // fprintf(stderr, "[M::%s-0-]\tq::[%ld,%ld)\test_bd::%ld\test_e::%ld\n", __func__, s, e, est_bd, est_e); + + s += nwl; e += nwl; e = ((e<=ql)?e:ql); + } + } else if (asm_opt.simd_mm == 1) { + for (; s < ql; ) {///[s, e) + if(rr) { + for (k = zk = in0; k < v32->n; k++) { + if((wsrt[v32->a[k]]>>32) == UINT32_MAX) continue;///passed + z = &(ol->list[(uint32_t)wsrt[v32->a[k]]]); + q[0] = z->x_pos_s; + q[1] = z->x_pos_e + 1; + ///[s, e) && [q[0], q[1]) + os = MAX(q[0], s); oe = MIN(q[1], e); + if(oe > os) v32->a[zk++] = v32->a[k]; + } + v32->n = zk; } + + for (; i < wsrt_n; ++i) { + z = &(ol->list[(uint32_t)wsrt[i]]); + q[0] = z->x_pos_s; + q[1] = z->x_pos_e + 1; + if(q[0] >= e) break; + if(z->is_match == 0) { + os = MAX(q[0], s); oe = MIN(q[1], e); + if(oe > os) kv_push(uint32_t, *v32, i); + } + } + + rr = hc_aln_simd_avx2(ol->list, wsrt, v32->a + in0, v32->n - in0, s, e, nwl, err, rref, qu->seq, tu, fi, aux_beg, aux_end, t_s, t_pri_l, exz, OVERLAP_THRESHOLD_HIFI_FILTER, 0, &tot_b); + // hc_est_robust_rr(ol->list, ql, ix->a + srt_n, ix->n - srt_n, ix->a + ix->n, ix->a + ix->n + ix->n - srt_n, s, e, 2.0, c_idx->a, &est_bd, &est_e, min_dp, 0); + // fprintf(stderr, "[M::%s-0-]\tq::[%ld,%ld)\test_bd::%ld\test_e::%ld\n", __func__, s, e, est_bd, est_e); + + s += nwl; e += nwl; e = ((e<=ql)?e:ql); } - - rr = hc_aln_simd(ol->list, wsrt, v32->a + in0, v32->n - in0, s, e, nwl, err, rref, qu->seq, tu, fi, aux_beg, aux_end, t_s, t_pri_l, exz, OVERLAP_THRESHOLD_HIFI_FILTER, 0, &tot_b); - // hc_est_robust_rr(ol->list, ql, ix->a + srt_n, ix->n - srt_n, ix->a + ix->n, ix->a + ix->n + ix->n - srt_n, s, e, 2.0, c_idx->a, &est_bd, &est_e, min_dp, 0); - // fprintf(stderr, "[M::%s-0-]\tq::[%ld,%ld)\test_bd::%ld\test_e::%ld\n", __func__, s, e, est_bd, est_e); - - s += nwl; e += nwl; e = ((e<=ql)?e:ql); } + ocn = v32->a; osc = v32->a + ol->length; // fprintf(stderr, "\n\n\n"); @@ -36054,33 +36240,65 @@ void batch_simd_ck(gen_hc_aln_t *ez, uint64_t *wsrt, uint64_t wsrt_n, uint64_t q ez->v32->n = (ez->ol->length<<1); in0 = ez->v32->n; i = 0; s = 0; e = w; e = ((e<=ql)?e:ql); rr = 0; - for (; s < ql; ) {///[s, e) - if(rr) { - for (k = zk = in0; k < ez->v32->n; k++) { - if((wsrt[ez->v32->a[k]]>>32) == UINT32_MAX) continue;///passed - z = &(ez->ol->list[(uint32_t)wsrt[ez->v32->a[k]]]); + + if(asm_opt.simd_mm == 2) { + for (; s < ql; ) {///[s, e) + if(rr) { + for (k = zk = in0; k < ez->v32->n; k++) { + if((wsrt[ez->v32->a[k]]>>32) == UINT32_MAX) continue;///passed + z = &(ez->ol->list[(uint32_t)wsrt[ez->v32->a[k]]]); + q[0] = z->x_pos_s; + q[1] = z->x_pos_e + 1; + ///[s, e) && [q[0], q[1]) + os = MAX(q[0], s); oe = MIN(q[1], e); + if(oe > os) ez->v32->a[zk++] = ez->v32->a[k]; + } + ez->v32->n = zk; + } + + for (; i < wsrt_n; ++i) { + z = &(ez->ol->list[(uint32_t)wsrt[i]]); q[0] = z->x_pos_s; q[1] = z->x_pos_e + 1; - ///[s, e) && [q[0], q[1]) - os = MAX(q[0], s); oe = MIN(q[1], e); - if(oe > os) ez->v32->a[zk++] = ez->v32->a[k]; + if(q[0] >= e) break; + if(z->is_match == 0) { + os = MAX(q[0], s); oe = MIN(q[1], e); + if(oe > os) kv_push(uint32_t, *(ez->v32), i); + } } - ez->v32->n = zk; - } - for (; i < wsrt_n; ++i) { - z = &(ez->ol->list[(uint32_t)wsrt[i]]); - q[0] = z->x_pos_s; - q[1] = z->x_pos_e + 1; - if(q[0] >= e) break; - if(z->is_match == 0) { - os = MAX(q[0], s); oe = MIN(q[1], e); - if(oe > os) kv_push(uint32_t, *(ez->v32), i); + rr = hc_aln_simd_avx512(ez->ol->list, wsrt, ez->v32->a + in0, ez->v32->n - in0, s, e, w, err, ez->rref, ez->qu->seq, ez->tu, fi, aux_beg, aux_end, t_s, t_pri_l, ez->exz, OVERLAP_THRESHOLD_HIFI_FILTER, 0, NULL/**&tot_b**/); + s += w; e += w; e = ((e<=ql)?e:ql); + } + } else if (asm_opt.simd_mm == 1) { + for (; s < ql; ) {///[s, e) + if(rr) { + for (k = zk = in0; k < ez->v32->n; k++) { + if((wsrt[ez->v32->a[k]]>>32) == UINT32_MAX) continue;///passed + z = &(ez->ol->list[(uint32_t)wsrt[ez->v32->a[k]]]); + q[0] = z->x_pos_s; + q[1] = z->x_pos_e + 1; + ///[s, e) && [q[0], q[1]) + os = MAX(q[0], s); oe = MIN(q[1], e); + if(oe > os) ez->v32->a[zk++] = ez->v32->a[k]; + } + ez->v32->n = zk; } - } - rr = hc_aln_simd(ez->ol->list, wsrt, ez->v32->a + in0, ez->v32->n - in0, s, e, w, err, ez->rref, ez->qu->seq, ez->tu, fi, aux_beg, aux_end, t_s, t_pri_l, ez->exz, OVERLAP_THRESHOLD_HIFI_FILTER, 0, NULL/**&tot_b**/); - s += w; e += w; e = ((e<=ql)?e:ql); + for (; i < wsrt_n; ++i) { + z = &(ez->ol->list[(uint32_t)wsrt[i]]); + q[0] = z->x_pos_s; + q[1] = z->x_pos_e + 1; + if(q[0] >= e) break; + if(z->is_match == 0) { + os = MAX(q[0], s); oe = MIN(q[1], e); + if(oe > os) kv_push(uint32_t, *(ez->v32), i); + } + } + + rr = hc_aln_simd_avx2(ez->ol->list, wsrt, ez->v32->a + in0, ez->v32->n - in0, s, e, w, err, ez->rref, ez->qu->seq, ez->tu, fi, aux_beg, aux_end, t_s, t_pri_l, ez->exz, OVERLAP_THRESHOLD_HIFI_FILTER, 0, NULL/**&tot_b**/); + s += w; e += w; e = ((e<=ql)?e:ql); + } } } diff --git a/Levenshtein_avx2.cpp b/Levenshtein_avx2.cpp new file mode 100644 index 0000000..879813d --- /dev/null +++ b/Levenshtein_avx2.cpp @@ -0,0 +1,162 @@ +#include "Levenshtein_distance.h" +#include + + +#define init_simd_ed4(PSA, PNA, THRE, ABS_DIAG, R_ERR, R_PE, SI, TN, CUT, BD, I, MM, PEQ_MM, LZ, IBD) {\ + (R_ERR)[(SI)] = INT32_MAX; (R_PE)[(SI)] = -1; (IBD)[(SI)] = ((THRE)<<1) - (ABS_DIAG)[(SI)];\ + if(((PNA)[(SI)] <= (TN) + (CUT)) && ((TN) <= (PNA)[(SI)] + (CUT))) {\ + (BD) = (((THRE)<<1)+1)-(ABS_DIAG)[(SI)]; (BD) = (((BD)<=(PNA)[(SI)])?(BD):(PNA)[(SI)]); (LZ) |= (((int32_t)1u) << (SI));\ + for ((I) = 0, (MM) = (((Word)1)<<((ABS_DIAG)[(SI)])); (I) < (BD); (I)++) {\ + (PEQ_MM)[seq_nt4_table[(uint8_t)(PSA)[(SI)][(I)]]][(SI)] |= (MM); (MM) <<= 1;\ + }\ + }\ +} + +#define ed_core_64x4(PEQz, VPz, VNz, Xz, D0z, HNz, HPz) { \ + /**(X) = (Peq)|(VN);**/\ + (Xz) = _mm256_or_si256((PEQz), (VNz)); \ + /**(D0) = (((VP) + ((X)&(VP))) ^ (VP)) | (X);**/\ + (D0z) = _mm256_or_si256(_mm256_xor_si256(_mm256_add_epi64((VPz), _mm256_and_si256((Xz), (VPz))), (VPz)), (Xz)); \ + /**(HN) = (VP)&(D0);**/\ + (HNz) = _mm256_and_si256((VPz), (D0z)); \ + /**(HP) = (VN) | ~((VP) | (D0));**/\ + (HPz) = _mm256_or_si256((VNz), _mm256_andnot_si256(_mm256_or_si256((VPz), (D0z)), _mm256_set1_epi64x(-1))); \ + /**(X) = (D0) >> 1;**/\ + (Xz) = _mm256_srli_epi64((D0z), 1); \ + /**(VN) = (X)&(HP);**/\ + (VNz) = _mm256_and_si256((Xz), (HPz)); \ + /**(VP) = (HN) | ~((X) | (HP));**/\ + (VPz) = _mm256_or_si256((HNz), _mm256_andnot_si256(_mm256_or_si256((Xz), (HPz)), _mm256_set1_epi64x(-1))); \ +} + +#define ed_core_upx4(PEQz, PSA, PNA, IBD, HT, CC, MMK, SI) { \ + if((HT) & (((int32_t)1u) << (SI))) {\ + (IBD)[(SI)]++;\ + if((IBD)[(SI)] < (PNA)[(SI)]) {\ + (CC) = seq_nt4_table[(uint8_t)(PSA)[(SI)][(IBD)[(SI)]]];\ + if((CC) < 4) (PEQz)[(CC)] = _mm256_or_si256((PEQz)[(CC)], (MMK)[(SI)]);\ + }\ + }\ +} + +#define ed_tail_upx4(HT, SI, ST, AI, PNA, ABS_DIAG, K, ERR_MM, VP_MM, VN_MM, THRE, R_ERR, R_PE, BD, I) {\ + if((HT) & (((int32_t)1u) << (SI))) {\ + (ST)[(SI)] -= (ABS_DIAG)[(SI)]; (AI)[(SI)] += (PNA)[(SI)] + (ABS_DIAG)[(SI)];\ + for ((K)[(SI)] = 0; (ST)[(SI)] < 0 && (K)[(SI)] < (AI)[(SI)]; (K)[(SI)]++, (ST)[(SI)]++) {\ + (ERR_MM)[(SI)] += ((VP_MM)[(SI)]&(1ULL)); (VP_MM)[(SI)]>>=1;\ + (ERR_MM)[(SI)] -= ((VN_MM)[(SI)]&(1ULL)); (VN_MM)[(SI)]>>=1;\ + }\ + if (((ERR_MM)[(SI)] <= (THRE)) && ((ERR_MM)[(SI)] <= (R_ERR)[(SI)])) {\ + (R_ERR)[(SI)] = (ERR_MM)[(SI)]; (R_PE)[(SI)] = (ST)[(SI)];\ + }\ + (ST)[(SI)] -= (K)[(SI)]; (BD)++; (I) = (SI);\ + }\ +} + +#define ED_TAIL_LANE(K) do { \ + if ((ht) & (((int32_t)1u) << (K))) { \ + st = tn - 1 - abs_diag_a[(K)]; \ + ai = pna[(K)] - tn + abs_diag_a[(K)]; \ + for (i = 0, uge = INT64_MAX; st < 0 && i < ai; i++, st++) { \ + err_mm[(K)] += ((VP_mm[(K)] >> i) & 1ULL); \ + err_mm[(K)] -= ((VN_mm[(K)] >> i) & 1ULL); \ + } \ + if ((err_mm[(K)] <= thre) && (err_mm[(K)] <= r_err[(K)])) { \ + r_err[(K)] = err_mm[(K)]; \ + r_pe[(K)] = st; \ + } \ + st -= i; \ + while (i < ai) { \ + err_mm[(K)] += ((VP_mm[(K)] >> i) & 1ULL); \ + err_mm[(K)] -= ((VN_mm[(K)] >> i) & 1ULL); \ + ++i; \ + if ((err_mm[(K)] <= thre) && (err_mm[(K)] <= r_err[(K)])) { \ + r_err[(K)] = err_mm[(K)]; \ + r_pe[(K)] = st + i; \ + } \ + if (i == thre) uge = err_mm[(K)]; \ + } \ + if ((uge <= thre) && (uge == r_err[(K)])) r_pe[(K)] = st + thre; \ + } \ +} while (0) + +void ed_band_cal_semi_64_w_absent_diag_avx4(char **psa, int32_t *pna, char *tstr, int32_t tn, int32_t thre, int32_t *abs_diag_a, int64_t *r_err, int64_t *r_pe) +{ + // r_err[0] = r_err[1] = r_err[2] = r_err[3] = r_err[4] = r_err[5] = r_err[6] = r_err[7] = thre+1; + // r_pe[0] = r_pe[1] = r_pe[2] = r_pe[3] = r_pe[4] = r_pe[5] = r_pe[6] = r_pe[7] = -1; + + Word mm, Peq_mm[5][AVX_GS2] = {{0}}, *VN_mm = NULL, *VP_mm = NULL, c = 0; __m256i Peq[5], VP, VN, X, D0, HN, HP, lone, E, C, mmk[AVX_GS2]; + int32_t lz = 0, ht = (((int32_t)1u)< + +#define init_simd_ed(PSA, PNA, THRE, ABS_DIAG, R_ERR, R_PE, SI, TN, CUT, BD, I, MM, PEQ_MM, LZ, IBD) {\ + (R_ERR)[(SI)] = INT32_MAX; (R_PE)[(SI)] = -1; (IBD)[(SI)] = ((THRE)<<1) - (ABS_DIAG)[(SI)];\ + if(((PNA)[(SI)] <= (TN) + (CUT)) && ((TN) <= (PNA)[(SI)] + (CUT))) {\ + (BD) = (((THRE)<<1)+1)-(ABS_DIAG)[(SI)]; (BD) = (((BD)<=(PNA)[(SI)])?(BD):(PNA)[(SI)]); (LZ) |= (((__mmask8)1u) << (SI));\ + for ((I) = 0, (MM) = (((Word)1)<<((ABS_DIAG)[(SI)])); (I) < (BD); (I)++) {\ + (PEQ_MM)[seq_nt4_table[(uint8_t)(PSA)[(SI)][(I)]]][(SI)] |= (MM); (MM) <<= 1;\ + }\ + }\ +} + +#define ed_core_64x8(PEQz, VPz, VNz, Xz, D0z, HNz, HPz) { \ + /**(X) = (Peq)|(VN);**/\ + (Xz) = _mm512_or_si512((PEQz), (VNz));\ + /**(D0) = (((VP) + ((X)&(VP))) ^ (VP)) | (X);**/\ + (D0z) = _mm512_or_si512(_mm512_xor_si512(_mm512_add_epi64((VPz), _mm512_and_si512((Xz), (VPz))), (VPz)), (Xz));\ + /**(HN) = (VP)&(D0);**/\ + (HNz) = _mm512_and_si512((VPz), (D0z));\ + /**(HP) = (VN) | ~((VP) | (D0));**/\ + (HPz) = _mm512_or_si512((VNz), _mm512_andnot_si512(_mm512_or_si512((VPz), (D0z)), _mm512_set1_epi64(-1)));\ + /**(X) = (D0) >> 1;**/\ + (Xz) = _mm512_srli_epi64((D0z), 1);\ + /**(VN) = (X)&(HP);**/\ + (VNz) = _mm512_and_si512((Xz), (HPz));\ + /**(VP) = (HN) | ~((X) | (HP));**/\ + (VPz) = _mm512_or_si512((HNz), _mm512_andnot_si512(_mm512_or_si512((Xz), (HPz)), _mm512_set1_epi64(-1)));\ +} + +#define ed_core_upx8(PEQz, PSA, PNA, IBD, HT, CC, MMK, SI) { \ + if((HT) & (((__mmask8)1u) << (SI))) {\ + (IBD)[(SI)]++;\ + if((IBD)[(SI)] < (PNA)[(SI)]) {\ + (CC) = seq_nt4_table[(uint8_t)(PSA)[(SI)][(IBD)[(SI)]]];\ + if((CC) < 4) (PEQz)[(CC)] = _mm512_or_si512((PEQz)[(CC)], (MMK)[(SI)]);\ + }\ + }\ +} + +#define ed_tail_upx8(HT, SI, ST, AI, PNA, ABS_DIAG, K, ERR_MM, VP_MM, VN_MM, THRE, R_ERR, R_PE, BD, I) {\ + if((HT) & (((__mmask8)1u) << (SI))) {\ + (ST)[(SI)] -= (ABS_DIAG)[(SI)]; (AI)[(SI)] += (PNA)[(SI)] + (ABS_DIAG)[(SI)];\ + for ((K)[(SI)] = 0; (ST)[(SI)] < 0 && (K)[(SI)] < (AI)[(SI)]; (K)[(SI)]++, (ST)[(SI)]++) {\ + (ERR_MM)[(SI)] += ((VP_MM)[(SI)]&(1ULL)); (VP_MM)[(SI)]>>=1;\ + (ERR_MM)[(SI)] -= ((VN_MM)[(SI)]&(1ULL)); (VN_MM)[(SI)]>>=1;\ + }\ + if (((ERR_MM)[(SI)] <= (THRE)) && ((ERR_MM)[(SI)] <= (R_ERR)[(SI)])) {\ + (R_ERR)[(SI)] = (ERR_MM)[(SI)]; (R_PE)[(SI)] = (ST)[(SI)];\ + }\ + (ST)[(SI)] -= (K)[(SI)]; (BD)++; (I) = (SI);\ + }\ +} + +#define ed_tail_ck8(MBEST, SI, R_PE, ST, K, THRE, UGE_MM, ERR_MM, AI, HT) {\ + (K)[(SI)]++;\ + if((MBEST) & (((__mmask8)1u) << (SI))) {\ + (R_PE)[(SI)] = (ST)[(SI)] + (K)[(SI)];\ + }\ + if((K)[(SI)] >= (AI)[(SI)]) (HT) &= ~(((__mmask8)1u) << (SI));\ + if((K)[(SI)] == (THRE)) (UGE_MM)[(SI)] = (ERR_MM)[(SI)];\ +} + + +void ed_band_cal_semi_64_w_absent_diag_avx8(char **psa, int32_t *pna, char *tstr, int32_t tn, int32_t thre, int32_t *abs_diag_a, int64_t *r_err, int64_t *r_pe) +{ + // r_err[0] = r_err[1] = r_err[2] = r_err[3] = r_err[4] = r_err[5] = r_err[6] = r_err[7] = thre+1; + // r_pe[0] = r_pe[1] = r_pe[2] = r_pe[3] = r_pe[4] = r_pe[5] = r_pe[6] = r_pe[7] = -1; + /** + ed_band_cal_semi_64_w_absent_diag_avx4(psa, pna, tstr, tn, thre, abs_diag_a, r_err, r_pe); + ed_band_cal_semi_64_w_absent_diag_avx4(psa + 4, pna + 4, tstr, tn, thre, abs_diag_a + 4, r_err + 4, r_pe + 4); + return; + **/ + + + Word mm, Peq_mm[5][AVX_GS] = {{0}}, *VN_mm = NULL, *VP_mm = NULL, c = 0; __m512i Peq[5], VP, VN, X, D0, HN, HP, lone, E, C, bestE, bestPE, curPE, cutPE, threPE, ugE, mmk[AVX_GS]; + __mmask8 lz = ((__mmask8)0u), ht = (((__mmask8)1u)< 1) { + VN = _mm512_loadu_si512(VN_mm); VP = _mm512_loadu_si512(VP_mm); E = _mm512_loadu_si512(err_mm); i = 0; + + bestE = _mm512_loadu_si512(r_err); ///threE = _mm512_set1_epi64(thre); + if(k[0] >= ai[0]) ht &= ((__mmask8)(255-1)); + if(k[1] >= ai[1]) ht &= ((__mmask8)(255-2)); + if(k[2] >= ai[2]) ht &= ((__mmask8)(255-4)); + if(k[3] >= ai[3]) ht &= ((__mmask8)(255-8)); + if(k[4] >= ai[4]) ht &= ((__mmask8)(255-16)); + if(k[5] >= ai[5]) ht &= ((__mmask8)(255-32)); + if(k[6] >= ai[6]) ht &= ((__mmask8)(255-64)); + if(k[7] >= ai[7]) ht &= ((__mmask8)(255-128)); + + err_mm[0] = r_pe[0]; err_mm[1] = r_pe[1]; err_mm[2] = r_pe[2]; err_mm[3] = r_pe[3]; + err_mm[4] = r_pe[4]; err_mm[5] = r_pe[5]; err_mm[6] = r_pe[6]; err_mm[7] = r_pe[7]; + bestPE = _mm512_loadu_si512(err_mm); + err_mm[0] = st[0] + k[0]; err_mm[1] = st[1] + k[1]; err_mm[2] = st[2] + k[2]; err_mm[3] = st[3] + k[3]; + err_mm[4] = st[4] + k[4]; err_mm[5] = st[5] + k[5]; err_mm[6] = st[6] + k[6]; err_mm[7] = st[7] + k[7]; + curPE = _mm512_loadu_si512(err_mm); + err_mm[0] = st[0] + thre; err_mm[1] = st[1] + thre; err_mm[2] = st[2] + thre; err_mm[3] = st[3] + thre; + err_mm[4] = st[4] + thre; err_mm[5] = st[5] + thre; err_mm[6] = st[6] + thre; err_mm[7] = st[7] + thre; + threPE = _mm512_loadu_si512(err_mm); + err_mm[0] = st[0] + ai[0]; err_mm[1] = st[1] + ai[1]; err_mm[2] = st[2] + ai[2]; err_mm[3] = st[3] + ai[3]; + err_mm[4] = st[4] + ai[4]; err_mm[5] = st[5] + ai[5]; err_mm[6] = st[6] + ai[6]; err_mm[7] = st[7] + ai[7]; + cutPE = _mm512_loadu_si512(err_mm); + + ugE = _mm512_loadu_si512(uge_mm); + + // mtf = _mm512_cmpge_epi64_mask(curPE, threPE) | ((__mmask8)(~ht)); + mtf = _mm512_cmpge_epi64_mask(curPE, threPE); + + while ((ht != 0) && ((mtf|((__mmask8)(~ht))) != (__mmask8)255)) { + E = _mm512_add_epi64(E, _mm512_and_si512(VP, lone)); VP = _mm512_srli_epi64(VP, 1); + E = _mm512_sub_epi64(E, _mm512_and_si512(VN, lone)); VN = _mm512_srli_epi64(VN, 1); + // i++; + + curPE = _mm512_add_epi64(curPE, lone); + ht &= _mm512_cmple_epi64_mask(curPE, cutPE); + if (ht == 0) break; + + mbest = _mm512_cmple_epi64_mask(E, bestE) & ht; + + bestE = _mm512_mask_mov_epi64(bestE, mbest, E); + bestPE = _mm512_mask_mov_epi64(bestPE, mbest, curPE); + + mt = _mm512_cmpeq_epi64_mask(curPE, threPE); + ugE = _mm512_mask_mov_epi64(ugE, mt&ht, E); + + mtf |= mt; + + // if(mbest && i < thre) _mm512_storeu_si512(err_mm, E); + + // ed_tail_ck8(mbest, 0, r_pe, st, k, thre, uge_mm, err_mm, ai, ht); + // ed_tail_ck8(mbest, 1, r_pe, st, k, thre, uge_mm, err_mm, ai, ht); + // ed_tail_ck8(mbest, 2, r_pe, st, k, thre, uge_mm, err_mm, ai, ht); + // ed_tail_ck8(mbest, 3, r_pe, st, k, thre, uge_mm, err_mm, ai, ht); + // ed_tail_ck8(mbest, 4, r_pe, st, k, thre, uge_mm, err_mm, ai, ht); + // ed_tail_ck8(mbest, 5, r_pe, st, k, thre, uge_mm, err_mm, ai, ht); + // ed_tail_ck8(mbest, 6, r_pe, st, k, thre, uge_mm, err_mm, ai, ht); + // ed_tail_ck8(mbest, 7, r_pe, st, k, thre, uge_mm, err_mm, ai, ht); + } + + + while (ht != 0) { + E = _mm512_add_epi64(E, _mm512_and_si512(VP, lone)); VP = _mm512_srli_epi64(VP, 1); + E = _mm512_sub_epi64(E, _mm512_and_si512(VN, lone)); VN = _mm512_srli_epi64(VN, 1); + // i++; + + curPE = _mm512_add_epi64(curPE, lone); + ht &= _mm512_cmple_epi64_mask(curPE, cutPE); + if (ht == 0) break; + + mbest = _mm512_cmple_epi64_mask(E, bestE) & ht; + + bestE = _mm512_mask_mov_epi64(bestE, mbest, E); + bestPE = _mm512_mask_mov_epi64(bestPE, mbest, curPE); + } + + cutPE = _mm512_set1_epi64(thre); + ht = _mm512_cmpgt_epi64_mask(bestE, cutPE); + bestE = _mm512_mask_set1_epi64(bestE, ht, INT32_MAX); + bestPE = _mm512_mask_set1_epi64(bestPE, ht, -1); + + ht = _mm512_cmple_epi64_mask(ugE, cutPE) & _mm512_cmpeq_epi64_mask(ugE, bestE); + bestPE = _mm512_mask_mov_epi64(bestPE, ht, threPE); + + _mm512_storeu_si512(r_err, bestE); + _mm512_storeu_si512(r_pe, bestPE); + _mm512_storeu_si512(uge_mm, ugE); + } else {///bd == 1 + while (k[i] < ai[i]) { + err_mm[i] += (VP_mm[i]&(1ULL)); VP_mm[i]>>=1; + err_mm[i] -= (VN_mm[i]&(1ULL)); VN_mm[i]>>=1; + ++k[i]; + if ((err_mm[i] <= thre) && (err_mm[i] <= r_err[i])) { + r_err[i] = err_mm[i]; r_pe[i] = st[i] + k[i]; + } + if(k[i] == thre) uge_mm[i] = err_mm[i]; + } + if((uge_mm[i] <= thre) && (uge_mm[i] == r_err[i])) r_pe[i] = st[i] + thre; + } +} diff --git a/Levenshtein_distance.h b/Levenshtein_distance.h index 039700c..d53b2a9 100644 --- a/Levenshtein_distance.h +++ b/Levenshtein_distance.h @@ -13,6 +13,7 @@ #include "kvec.h" #define AVX_GS 8 +#define AVX_GS2 4 extern const unsigned char seq_nt4_table[256]; typedef uint64_t Word; @@ -3736,275 +3737,9 @@ inline void ed_band_cal_extension_64_1_w_trace(char *pstr, int32_t pn, char *tst return; } +void ed_band_cal_semi_64_w_absent_diag_avx4(char **psa, int32_t *pna, char *tstr, int32_t tn, int32_t thre, int32_t *abs_diag_a, int64_t *r_err, int64_t *r_pe); -#define init_simd_ed(PSA, PNA, THRE, ABS_DIAG, R_ERR, R_PE, SI, TN, CUT, BD, I, MM, PEQ_MM, LZ, IBD) {\ - (R_ERR)[(SI)] = INT32_MAX; (R_PE)[(SI)] = -1; (IBD)[(SI)] = ((THRE)<<1) - (ABS_DIAG)[(SI)];\ - if(((PNA)[(SI)] <= (TN) + (CUT)) && ((TN) <= (PNA)[(SI)] + (CUT))) {\ - (BD) = (((THRE)<<1)+1)-(ABS_DIAG)[(SI)]; (BD) = (((BD)<=(PNA)[(SI)])?(BD):(PNA)[(SI)]); (LZ) |= (((__mmask8)1u) << (SI));\ - for ((I) = 0, (MM) = (((Word)1)<<((ABS_DIAG)[(SI)])); (I) < (BD); (I)++) {\ - (PEQ_MM)[seq_nt4_table[(uint8_t)(PSA)[(SI)][(I)]]][(SI)] |= (MM); (MM) <<= 1;\ - }\ - }\ -} - -#define ed_core_64x8(PEQz, VPz, VNz, Xz, D0z, HNz, HPz) { \ - /**(X) = (Peq)|(VN);**/\ - (Xz) = _mm512_or_si512((PEQz), (VNz));\ - /**(D0) = (((VP) + ((X)&(VP))) ^ (VP)) | (X);**/\ - (D0z) = _mm512_or_si512(_mm512_xor_si512(_mm512_add_epi64((VPz), _mm512_and_si512((Xz), (VPz))), (VPz)), (Xz));\ - /**(HN) = (VP)&(D0);**/\ - (HNz) = _mm512_and_si512((VPz), (D0z));\ - /**(HP) = (VN) | ~((VP) | (D0));**/\ - (HPz) = _mm512_or_si512((VNz), _mm512_andnot_si512(_mm512_or_si512((VPz), (D0z)), _mm512_set1_epi64(-1)));\ - /**(X) = (D0) >> 1;**/\ - (Xz) = _mm512_srli_epi64((D0z), 1);\ - /**(VN) = (X)&(HP);**/\ - (VNz) = _mm512_and_si512((Xz), (HPz));\ - /**(VP) = (HN) | ~((X) | (HP));**/\ - (VPz) = _mm512_or_si512((HNz), _mm512_andnot_si512(_mm512_or_si512((Xz), (HPz)), _mm512_set1_epi64(-1)));\ -} - -#define ed_core_upx8(PEQz, PSA, PNA, IBD, HT, CC, MMK, SI) { \ - if((HT) & (((__mmask8)1u) << (SI))) {\ - (IBD)[(SI)]++;\ - if((IBD)[(SI)] < (PNA)[(SI)]) {\ - (CC) = seq_nt4_table[(uint8_t)(PSA)[(SI)][(IBD)[(SI)]]];\ - if((CC) < 4) (PEQz)[(CC)] = _mm512_or_si512((PEQz)[(CC)], (MMK)[(SI)]);\ - }\ - }\ -} - -#define ed_tail_upx8(HT, SI, ST, AI, PNA, ABS_DIAG, K, ERR_MM, VP_MM, VN_MM, THRE, R_ERR, R_PE, BD, I) {\ - if((HT) & (((__mmask8)1u) << (SI))) {\ - (ST)[(SI)] -= (ABS_DIAG)[(SI)]; (AI)[(SI)] += (PNA)[(SI)] + (ABS_DIAG)[(SI)];\ - for ((K)[(SI)] = 0; (ST)[(SI)] < 0 && (K)[(SI)] < (AI)[(SI)]; (K)[(SI)]++, (ST)[(SI)]++) {\ - (ERR_MM)[(SI)] += ((VP_MM)[(SI)]&(1ULL)); (VP_MM)[(SI)]>>=1;\ - (ERR_MM)[(SI)] -= ((VN_MM)[(SI)]&(1ULL)); (VN_MM)[(SI)]>>=1;\ - }\ - if (((ERR_MM)[(SI)] <= (THRE)) && ((ERR_MM)[(SI)] <= (R_ERR)[(SI)])) {\ - (R_ERR)[(SI)] = (ERR_MM)[(SI)]; (R_PE)[(SI)] = (ST)[(SI)];\ - }\ - (ST)[(SI)] -= (K)[(SI)]; (BD)++; (I) = (SI);\ - }\ -} - -#define ed_tail_ck8(MBEST, SI, R_PE, ST, K, THRE, UGE_MM, ERR_MM, AI, HT) {\ - (K)[(SI)]++;\ - if((MBEST) & (((__mmask8)1u) << (SI))) {\ - (R_PE)[(SI)] = (ST)[(SI)] + (K)[(SI)];\ - }\ - if((K)[(SI)] >= (AI)[(SI)]) (HT) &= ~(((__mmask8)1u) << (SI));\ - if((K)[(SI)] == (THRE)) (UGE_MM)[(SI)] = (ERR_MM)[(SI)];\ -} - - -inline void ed_band_cal_semi_64_w_absent_diag_avx8(char **psa, int32_t *pna, char *tstr, int32_t tn, int32_t thre, int32_t *abs_diag_a, int64_t *r_err, int64_t *r_pe) -{ - // r_err[0] = r_err[1] = r_err[2] = r_err[3] = r_err[4] = r_err[5] = r_err[6] = r_err[7] = thre+1; - // r_pe[0] = r_pe[1] = r_pe[2] = r_pe[3] = r_pe[4] = r_pe[5] = r_pe[6] = r_pe[7] = -1; - - Word mm, Peq_mm[5][AVX_GS] = {{0}}, *VN_mm = NULL, *VP_mm = NULL, c = 0; __m512i Peq[5], VP, VN, X, D0, HN, HP, lone, E, C, bestE, bestPE, curPE, cutPE, threPE, ugE, mmk[AVX_GS]; - __mmask8 lz = ((__mmask8)0u), ht = (((__mmask8)1u)< 1) { - VN = _mm512_loadu_si512(VN_mm); VP = _mm512_loadu_si512(VP_mm); E = _mm512_loadu_si512(err_mm); i = 0; - - bestE = _mm512_loadu_si512(r_err); ///threE = _mm512_set1_epi64(thre); - if(k[0] >= ai[0]) ht &= ((__mmask8)(255-1)); - if(k[1] >= ai[1]) ht &= ((__mmask8)(255-2)); - if(k[2] >= ai[2]) ht &= ((__mmask8)(255-4)); - if(k[3] >= ai[3]) ht &= ((__mmask8)(255-8)); - if(k[4] >= ai[4]) ht &= ((__mmask8)(255-16)); - if(k[5] >= ai[5]) ht &= ((__mmask8)(255-32)); - if(k[6] >= ai[6]) ht &= ((__mmask8)(255-64)); - if(k[7] >= ai[7]) ht &= ((__mmask8)(255-128)); - - err_mm[0] = r_pe[0]; err_mm[1] = r_pe[1]; err_mm[2] = r_pe[2]; err_mm[3] = r_pe[3]; - err_mm[4] = r_pe[4]; err_mm[5] = r_pe[5]; err_mm[6] = r_pe[6]; err_mm[7] = r_pe[7]; - bestPE = _mm512_loadu_si512(err_mm); - err_mm[0] = st[0] + k[0]; err_mm[1] = st[1] + k[1]; err_mm[2] = st[2] + k[2]; err_mm[3] = st[3] + k[3]; - err_mm[4] = st[4] + k[4]; err_mm[5] = st[5] + k[5]; err_mm[6] = st[6] + k[6]; err_mm[7] = st[7] + k[7]; - curPE = _mm512_loadu_si512(err_mm); - err_mm[0] = st[0] + thre; err_mm[1] = st[1] + thre; err_mm[2] = st[2] + thre; err_mm[3] = st[3] + thre; - err_mm[4] = st[4] + thre; err_mm[5] = st[5] + thre; err_mm[6] = st[6] + thre; err_mm[7] = st[7] + thre; - threPE = _mm512_loadu_si512(err_mm); - err_mm[0] = st[0] + ai[0]; err_mm[1] = st[1] + ai[1]; err_mm[2] = st[2] + ai[2]; err_mm[3] = st[3] + ai[3]; - err_mm[4] = st[4] + ai[4]; err_mm[5] = st[5] + ai[5]; err_mm[6] = st[6] + ai[6]; err_mm[7] = st[7] + ai[7]; - cutPE = _mm512_loadu_si512(err_mm); - - ugE = _mm512_loadu_si512(uge_mm); - - // mtf = _mm512_cmpge_epi64_mask(curPE, threPE) | ((__mmask8)(~ht)); - mtf = _mm512_cmpge_epi64_mask(curPE, threPE); - - while ((ht != 0) && ((mtf|((__mmask8)(~ht))) != (__mmask8)255)) { - E = _mm512_add_epi64(E, _mm512_and_si512(VP, lone)); VP = _mm512_srli_epi64(VP, 1); - E = _mm512_sub_epi64(E, _mm512_and_si512(VN, lone)); VN = _mm512_srli_epi64(VN, 1); - // i++; - - curPE = _mm512_add_epi64(curPE, lone); - ht &= _mm512_cmple_epi64_mask(curPE, cutPE); - if (ht == 0) break; - - mbest = _mm512_cmple_epi64_mask(E, bestE) & ht; - - bestE = _mm512_mask_mov_epi64(bestE, mbest, E); - bestPE = _mm512_mask_mov_epi64(bestPE, mbest, curPE); - - mt = _mm512_cmpeq_epi64_mask(curPE, threPE); - ugE = _mm512_mask_mov_epi64(ugE, mt&ht, E); - - mtf |= mt; - - // if(mbest && i < thre) _mm512_storeu_si512(err_mm, E); - - // ed_tail_ck8(mbest, 0, r_pe, st, k, thre, uge_mm, err_mm, ai, ht); - // ed_tail_ck8(mbest, 1, r_pe, st, k, thre, uge_mm, err_mm, ai, ht); - // ed_tail_ck8(mbest, 2, r_pe, st, k, thre, uge_mm, err_mm, ai, ht); - // ed_tail_ck8(mbest, 3, r_pe, st, k, thre, uge_mm, err_mm, ai, ht); - // ed_tail_ck8(mbest, 4, r_pe, st, k, thre, uge_mm, err_mm, ai, ht); - // ed_tail_ck8(mbest, 5, r_pe, st, k, thre, uge_mm, err_mm, ai, ht); - // ed_tail_ck8(mbest, 6, r_pe, st, k, thre, uge_mm, err_mm, ai, ht); - // ed_tail_ck8(mbest, 7, r_pe, st, k, thre, uge_mm, err_mm, ai, ht); - } - - - while (ht != 0) { - E = _mm512_add_epi64(E, _mm512_and_si512(VP, lone)); VP = _mm512_srli_epi64(VP, 1); - E = _mm512_sub_epi64(E, _mm512_and_si512(VN, lone)); VN = _mm512_srli_epi64(VN, 1); - // i++; - - curPE = _mm512_add_epi64(curPE, lone); - ht &= _mm512_cmple_epi64_mask(curPE, cutPE); - if (ht == 0) break; - - mbest = _mm512_cmple_epi64_mask(E, bestE) & ht; - - bestE = _mm512_mask_mov_epi64(bestE, mbest, E); - bestPE = _mm512_mask_mov_epi64(bestPE, mbest, curPE); - } - - cutPE = _mm512_set1_epi64(thre); - ht = _mm512_cmpgt_epi64_mask(bestE, cutPE); - bestE = _mm512_mask_set1_epi64(bestE, ht, INT32_MAX); - bestPE = _mm512_mask_set1_epi64(bestPE, ht, -1); - - ht = _mm512_cmple_epi64_mask(ugE, cutPE) & _mm512_cmpeq_epi64_mask(ugE, bestE); - bestPE = _mm512_mask_mov_epi64(bestPE, ht, threPE); - - _mm512_storeu_si512(r_err, bestE); - _mm512_storeu_si512(r_pe, bestPE); - _mm512_storeu_si512(uge_mm, ugE); - } else {///bd == 1 - while (k[i] < ai[i]) { - err_mm[i] += (VP_mm[i]&(1ULL)); VP_mm[i]>>=1; - err_mm[i] -= (VN_mm[i]&(1ULL)); VN_mm[i]>>=1; - ++k[i]; - if ((err_mm[i] <= thre) && (err_mm[i] <= r_err[i])) { - r_err[i] = err_mm[i]; r_pe[i] = st[i] + k[i]; - } - if(k[i] == thre) uge_mm[i] = err_mm[i]; - } - if((uge_mm[i] <= thre) && (uge_mm[i] == r_err[i])) r_pe[i] = st[i] + thre; - } -} - +void ed_band_cal_semi_64_w_absent_diag_avx8(char **psa, int32_t *pna, char *tstr, int32_t tn, int32_t thre, int32_t *abs_diag_a, int64_t *r_err, int64_t *r_pe); inline void ed_band_cal_semi_64_w_absent_diag(char *pstr, int32_t pn, char *tstr, int32_t tn, int32_t thre, int32_t abs_diag, bit_extz_t *ez) { diff --git a/Makefile b/Makefile index c33c5ce..1d450a2 100644 --- a/Makefile +++ b/Makefile @@ -1,11 +1,11 @@ CXX= g++ CC= gcc -CXXFLAGS= -g -O3 -mavx512f -msse4.2 -mpopcnt -fomit-frame-pointer -Wall +CXXFLAGS= -g -O3 -msse4.2 -mpopcnt -fomit-frame-pointer -Wall CFLAGS= $(CXXFLAGS) CPPFLAGS= INCLUDES= OBJS= CommandLines.o Process_Read.o Assembly.o Hash_Table.o \ - POA.o Correct.o Levenshtein_distance.o Overlaps.o Trio.o kthread.o Purge_Dups.o \ + POA.o Correct.o Levenshtein_distance.o Levenshtein_avx2.o Levenshtein_avx512.o Overlaps.o Trio.o kthread.o Purge_Dups.o \ htab.o hist.o sketch.o anchor.o extract.o sys.o hic.o rcut.o horder.o ecovlp.o\ tovlp.o inter.o kalloc.o gfa_ut.o gchain_map.o EXE= hifiasm @@ -19,13 +19,22 @@ endif .SUFFIXES:.cpp .c .o .PHONY:all clean depend +all:$(EXE) + .cpp.o: $(CXX) -c $(CXXFLAGS) $(CPPFLAGS) $(INCLUDES) $< -o $@ .c.o: $(CC) -c $(CFLAGS) $(CPPFLAGS) $(INCLUDES) $< -o $@ -all:$(EXE) +# compiled only with AVX2 +Levenshtein_avx2.o: Levenshtein_avx2.cpp Levenshtein_distance.h + $(CXX) -c $(CXXFLAGS) -mavx2 $(CPPFLAGS) $(INCLUDES) $< -o $@ + +# compiled only with AVX512 +Levenshtein_avx512.o: Levenshtein_avx512.cpp Levenshtein_distance.h + $(CXX) -c $(CXXFLAGS) -mavx512f $(CPPFLAGS) $(INCLUDES) $< -o $@ + $(EXE):$(OBJS) main.o $(CXX) $(CXXFLAGS) $^ -o $@ $(LIBS) diff --git a/Overlaps.cpp b/Overlaps.cpp index 3d5716a..c0a3596 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -13918,6 +13918,10 @@ ma_hit_t_alloc* sources, R_to_U* ruIndex, int max_hang, int min_ovlp) ug = ma_ug_gen(sg); ma_ug_seq(ug, sg, coverage_cut, sources, &new_rtg_edges, max_hang, min_ovlp, 0, 1); + // if(asm_opt.dbg_bam) { + // gen_gfa_bam(ug); + // } + fprintf(stderr, "Writing raw unitig GFA to disk... \n"); char* gfa_name = (char*)malloc(strlen(output_file_name)+25); sprintf(gfa_name, "%s.r_utg.gfa", output_file_name); @@ -23743,6 +23747,11 @@ void write_all_data_to_disk(ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sou sprintf(gfa_name, "%s.ovlp.reverse", output_file_name); write_ma_hit_ts(reverse_sources, RNF->total_reads, gfa_name); + if(asm_opt.dbg_bam) { + sprintf(gfa_name, "%s.rec", output_file_name); + write_cc_v(&scb, gfa_name); + } + free(gfa_name); fprintf(stderr, "bin files have been written.\n"); if(asm_opt.bin_only) exit(0); @@ -23767,8 +23776,7 @@ int load_all_data_from_disk(ma_hit_t_alloc **sources, ma_hit_t_alloc **reverse_s } } - if((asm_opt.flag & HA_F_VERBOSE_GFA) && load_debug_graph(NULL, NULL, NULL, output_file_name, NULL, NULL, NULL)) - { + if((asm_opt.flag & HA_F_VERBOSE_GFA) && load_debug_graph(NULL, NULL, NULL, output_file_name, NULL, NULL, NULL)) { (*sources) = NULL; (*reverse_sources) = NULL; free(gfa_name); @@ -23785,6 +23793,11 @@ int load_all_data_from_disk(ma_hit_t_alloc **sources, ma_hit_t_alloc **reverse_s free(gfa_name); return 0; } + + if(asm_opt.dbg_bam) { + sprintf(gfa_name, "%s.rec", output_file_name); + load_cc_v(&scb, gfa_name); ///write_ec_reads("lec.raw.fa", &scb, 0); + } free(gfa_name); if(asm_opt.write_pos_idx) { diff --git a/Process_Read.cpp b/Process_Read.cpp index 673e657..b88d50d 100644 --- a/Process_Read.cpp +++ b/Process_Read.cpp @@ -247,6 +247,55 @@ int load_All_reads(All_reads* r, char* read_file_name) return 1; } + +void write_cc_v(cc_v* r, char* read_file_name) +{ + fprintf(stderr, "Writing raw reads to disk... \n"); + char* index_name = (char*)malloc(strlen(read_file_name)+32); + sprintf(index_name, "%s.bin", read_file_name); + FILE* fp = fopen(index_name, "w"); free(index_name); + + ///typedef struct {size_t n, m; asg16_v *a; uint8_t *f; uint16_t *er; uint64_t bid;} cc_v; + asg16_v *z; + uint64_t k, rn = r->n; uint32_t zn; fwrite(&rn, sizeof(rn), 1, fp); + for (k = 0; k < r->n; k++) { + z = &(r->a[k]); + zn = z->n; + fwrite(&zn, sizeof(zn), 1, fp); + fwrite(z->a, sizeof((*(z->a))), zn, fp); + } + + fflush(fp); fclose(fp); + fprintf(stderr, "Raw reads has been written.\n"); +} + +uint8_t load_cc_v(cc_v* r, char* read_file_name) +{ + fprintf(stderr, "Loading raw reads... \n"); + char* index_name = (char*)malloc(strlen(read_file_name)+32); + sprintf(index_name, "%s.bin", read_file_name); + FILE* fp = fopen(index_name, "r"); free(index_name); + if (!fp) { + fprintf(stderr, "No raw read bin.\n"); + return 0; + } + + ///typedef struct {size_t n, m; asg16_v *a; uint8_t *f; uint16_t *er; uint64_t bid;} cc_v; + asg16_v *z; int f_flag = 0; + uint64_t k, rn; uint32_t zn; f_flag += fread(&rn, sizeof(rn), 1, fp); + r->n = r->m = rn; MALLOC(r->a, r->n); + for (k = 0; k < r->n; k++) { + z = &(r->a[k]); + f_flag += fread(&zn, sizeof(zn), 1, fp); + z->n = z->m = zn; MALLOC(z->a, z->n); + f_flag += fread(z->a, sizeof((*(z->a))), zn, fp); + } + + fflush(fp); fclose(fp); + fprintf(stderr, "Raw reads has been loaded.\n"); + return 1; +} + void read_ma(ma_hit_t* x, FILE* fp) { int f_flag; diff --git a/Process_Read.h b/Process_Read.h index 2205242..b8486f5 100644 --- a/Process_Read.h +++ b/Process_Read.h @@ -8,6 +8,7 @@ #include #include "Overlaps.h" #include "CommandLines.h" +#include "Levenshtein_distance.h" ///#include "Hash_Table.h" #define READ_INIT_NUMBER 1000 @@ -151,6 +152,9 @@ typedef struct extern All_reads R_INF; +typedef struct {size_t n, m; asg16_v *a; uint8_t *f; uint16_t *er; uint64_t bid;} cc_v; +extern cc_v scb; + typedef struct { char* seq; @@ -240,6 +244,8 @@ void destory_UC_Read(UC_Read* r); void reverse_complement(char* pattern, uint64_t length); void write_All_reads(All_reads* r, char* read_file_name); int load_All_reads(All_reads* r, char* read_file_name); +uint8_t load_cc_v(cc_v* r, char* read_file_name); +void write_cc_v(cc_v* r, char* read_file_name); int append_All_reads(All_reads* r, char *idx, uint32_t id); void destory_All_reads(All_reads* r); int destory_read_bin(All_reads* r); diff --git a/ecovlp.cpp b/ecovlp.cpp index 50213fe..579988e 100644 --- a/ecovlp.cpp +++ b/ecovlp.cpp @@ -120,7 +120,6 @@ KRADIX_SORT_INIT(ec64, uint64_t, generic_key, 8) #define kdq_clear(q) ((q)->count = (q)->front = 0) -typedef struct {size_t n, m; asg16_v *a; uint8_t *f; uint16_t *er; uint64_t bid;} cc_v; cc_v scc = {0, 0, NULL, NULL, NULL, 0}; cc_v scb = {0, 0, NULL, NULL, NULL, 0}; cc_v sca = {0, 0, NULL, NULL, NULL, 0}; @@ -2842,7 +2841,7 @@ void push_ne_ovlp_flt(ma_hit_t_alloc* paf, overlap_region_alloc* ov, uint32_t fl } } -void push_ne_ovlp(ma_hit_t_alloc* paf, overlap_region_alloc* ov, uint32_t flag, All_reads* R_INF, asg16_v *ec/**, uint64_t qid, UC_Read *qu, UC_Read *tu**/) +void push_ne_ovlp_back(ma_hit_t_alloc* paf, overlap_region_alloc* ov, uint32_t flag, All_reads* R_INF, asg16_v *ec/**, uint64_t qid, UC_Read *qu, UC_Read *tu**/) { // if(qu && tu) { // debug_extract_max_exact_sub(qid, qu, tu); @@ -2899,6 +2898,67 @@ void push_ne_ovlp(ma_hit_t_alloc* paf, overlap_region_alloc* ov, uint32_t flag, } } +void push_ne_ovlp(ma_hit_t_alloc* paf, overlap_region_alloc* ov, uint32_t flag, All_reads* R_INF, asg16_v *ec/**, uint64_t qid, UC_Read *qu, UC_Read *tu**/) +{ + // if(qu && tu) { + // debug_extract_max_exact_sub(qid, qu, tu); + // } + uint64_t k, n; ma_hit_t *z; uint32_t rxs, rxe, rys, rye; + for (k = n = 0; k < ov->length; k++) { + if(ov->list[k].is_match == flag) n++; + } + + if(n > paf->size) { + paf->size = n; + REALLOC(paf->buffer, paf->size); + } + + for (k = paf->length = 0; k < ov->length; k++) { + if(ov->list[k].is_match == flag) { + // fprintf(stderr, "@%s\tSN:%.*s(id::%u)\terr::%u\n", flag==1?"SQ":"RQ", (int32_t)Get_NAME_LENGTH((*R_INF), ov->list[k].y_id), Get_NAME((*R_INF), ov->list[k].y_id), ov->list[k].y_id, ov->list[k].non_homopolymer_errors); + + z = &(paf->buffer[paf->length++]); + + z->qns = ov->list[k].x_id; + z->qns = z->qns << 32; + z->tn = ov->list[k].y_id; + + z->qns = z->qns | (uint64_t)(ov->list[k].x_pos_s); + z->qe = ov->list[k].x_pos_e + 1; + z->ts = ov->list[k].y_pos_s; + z->te = ov->list[k].y_pos_e + 1; + + ///for overlap_list, the x_strand of all overlaps are 0, so the tmp.rev is the same as the y_strand + z->rev = ov->list[k].y_pos_strand; + + z->bl = Get_READ_LENGTH((*R_INF), ov->list[k].y_id); + z->ml = ov->list[k].strong; + z->no_l_indel = ov->list[k].without_large_indel; + + z->el = 0; + if(ec) { + extract_max_exact(&ov->list[k], ec, /**qu, tu,**/ &rxs, &rxe, &rys, &rye); + // z->el = 0; + // fprintf(stderr, "[M::%s]\tq::[%u,\t%u)\tt::[%u,\t%u)\teq::[%u,\t%u)\tet::[%u,\t%u)\n", __func__, ov->list[k].x_pos_s, ov->list[k].x_pos_e + 1, ov->list[k].y_pos_s, ov->list[k].y_pos_e + 1, rxs, rxe, rys, rye); + if(rxe > rxs) { + // assert(rxe - rxs == rye - rys); + // z->qns = ov->list[k].x_id; + // z->qns = z->qns << 32; + // z->qns = z->qns | (uint64_t)(rxs); + // z->qe = rxe; + // z->ts = rys; + // z->te = rye; + + z->qns = (((uint64_t)(rxs))<<32) | ((uint64_t)((uint32_t)z->qns)); + z->bl = rys; z->cc = rxe - rxs; + + z->el = 1; + } + } + } + } +} + void pull_ovlp_syn(ma_hit_t_alloc* paf, uint64_t tqn) { @@ -3369,8 +3429,14 @@ uint64_t gen_hc_r_alin_ea_flt_mmp(ha_abuf_t *ab, overlap_region_alloc* ol, Candi if(!(srt->n)) { // gen_ff_hpc(hpz, qu->seq, qu->length, HPC_RR_Q, HPC_CC_Q); - tot_b = gen_hc_r_alin_adp_mmp_1(ol, cl, rref, qu, tu, exz, aux_o, e_rate, wl, rid, khit, move_gap, buf, chem_drop, align_gap_rate, align_gap_max, sec_aln_win, sec_aln_cov, sec_aln_err_rate, sec_aln_max, srt, ocw, NULL/**hpz->a**/, - v32, bp, max_n_chain>0?max_n_chain:1, max_n_chain_f>0?max_n_chain_f:1, chain_cutoff, ((asm_opt.hom_cov*HC_AV_MIN)>0)?(asm_opt.hom_cov*HC_AV_MIN):(1), 1); + if(asm_opt.simd_mm > 0) { + tot_b = gen_hc_r_alin_adp_mmp_1(ol, cl, rref, qu, tu, exz, aux_o, e_rate, wl, rid, khit, move_gap, buf, chem_drop, align_gap_rate, align_gap_max, sec_aln_win, sec_aln_cov, sec_aln_err_rate, sec_aln_max, srt, ocw, NULL/**hpz->a**/, + v32, bp, max_n_chain>0?max_n_chain:1, max_n_chain_f>0?max_n_chain_f:1, chain_cutoff, ((asm_opt.hom_cov*HC_AV_MIN)>0)?(asm_opt.hom_cov*HC_AV_MIN):(1), 1); + } else { + tot_b = gen_hc_r_alin_adp_mmp_0(ol, cl, rref, qu, tu, exz, aux_o, e_rate, wl, rid, khit, move_gap, buf, chem_drop, align_gap_rate, align_gap_max, sec_aln_win, sec_aln_cov, sec_aln_err_rate, sec_aln_max, srt, ocw, NULL/**hpz->a**/, + v32, bp, max_n_chain>0?max_n_chain:1, max_n_chain_f>0?max_n_chain_f:1, chain_cutoff, ((asm_opt.hom_cov*HC_AV_MIN)>0)?(asm_opt.hom_cov*HC_AV_MIN):(1), 1); + } + } else { kv_resize(uint64_t, *srt, (srt->n + ol->length)); ei = srt->a; en = srt->n; oi = srt->a + srt->n; on = ol->length; @@ -3401,8 +3467,13 @@ uint64_t gen_hc_r_alin_ea_flt_mmp(ha_abuf_t *ab, overlap_region_alloc* ol, Candi if(on > nec) { // gen_ff_hpc(hpz, qu->seq, qu->length, HPC_RR_Q, HPC_CC_Q); - tot_b = gen_hc_r_alin_adp_mmp_1(ol, cl, rref, qu, tu, exz, aux_o, e_rate, wl, rid, khit, move_gap, buf, chem_drop, align_gap_rate, align_gap_max, sec_aln_win, sec_aln_cov, sec_aln_err_rate, sec_aln_max, srt, ocw, NULL/**hpz->a**/, - v32, bp, max_n_chain>0?max_n_chain:1, max_n_chain_f>0?max_n_chain_f:1, chain_cutoff, ((asm_opt.hom_cov*HC_AV_MIN)>0)?(asm_opt.hom_cov*HC_AV_MIN):(1), 0); + if(asm_opt.simd_mm > 0) { + tot_b = gen_hc_r_alin_adp_mmp_1(ol, cl, rref, qu, tu, exz, aux_o, e_rate, wl, rid, khit, move_gap, buf, chem_drop, align_gap_rate, align_gap_max, sec_aln_win, sec_aln_cov, sec_aln_err_rate, sec_aln_max, srt, ocw, NULL/**hpz->a**/, + v32, bp, max_n_chain>0?max_n_chain:1, max_n_chain_f>0?max_n_chain_f:1, chain_cutoff, ((asm_opt.hom_cov*HC_AV_MIN)>0)?(asm_opt.hom_cov*HC_AV_MIN):(1), 0); + } else { + tot_b = gen_hc_r_alin_adp_mmp_0(ol, cl, rref, qu, tu, exz, aux_o, e_rate, wl, rid, khit, move_gap, buf, chem_drop, align_gap_rate, align_gap_max, sec_aln_win, sec_aln_cov, sec_aln_err_rate, sec_aln_max, srt, ocw, NULL/**hpz->a**/, + v32, bp, max_n_chain>0?max_n_chain:1, max_n_chain_f>0?max_n_chain_f:1, chain_cutoff, ((asm_opt.hom_cov*HC_AV_MIN)>0)?(asm_opt.hom_cov*HC_AV_MIN):(1), 0); + } } } @@ -3708,7 +3779,8 @@ void gen_hc_r_alin_ea_adv_flt_mmp(gen_hc_aln_t *ez) if(!(ez->srt->n)) { // gen_hc_r_alin_adv_adp_smp(ez, a_cu, a_ci, ocn, osc, idx_cu, n_cu, 1); - gen_hc_r_alin_adv_adp_smp_1(ez, 1); + if(asm_opt.simd_mm > 0) gen_hc_r_alin_adv_adp_smp_1(ez, 1); + else gen_hc_r_alin_adv_adp_smp_0(ez, 1); } else { ///debug for memory // snprintf(NULL, 0, "dwn::%u\tdcn::%u", (uint32_t)aux_o->w_list.n, (uint32_t)aux_o->w_list.c.n); @@ -3743,7 +3815,8 @@ void gen_hc_r_alin_ea_adv_flt_mmp(gen_hc_aln_t *ez) if(on > nec) { // gen_hc_r_alin_adv_adp_smp(ez, a_cu, a_ci, ocn, osc, idx_cu, n_cu, 0); - gen_hc_r_alin_adv_adp_smp_1(ez, 0); + if(asm_opt.simd_mm > 0) gen_hc_r_alin_adv_adp_smp_1(ez, 0); + else gen_hc_r_alin_adv_adp_smp_0(ez, 0); } // fprintf(stderr, "[M::%s] srt->n::%u, nec::%lu, on::%lu\n", __func__, (uint32_t)srt->n, nec, on); ///debug for memory @@ -4221,6 +4294,7 @@ static void worker_init_ec_step(void *data, long i, int tid) if(((het_r) < 0) || ((het_r) > ((hom_r)/(n_hap)))) {(het_r) = (hom_r)/(n_hap);}\ } while (0) +void update_scb(All_reads *R_INF, asg16_v *scc, asg16_v *scb, asg16_v *scb_res, UC_Read *qu, UC_Read *tu, asg64_v *srt, bit_extz_t *exz, uint64_t rid); static void worker_hap_ec(void *data, long i, int tid) @@ -4386,7 +4460,7 @@ static void worker_hap_ec(void *data, long i, int tid) copy_asg_arr(buf0, b->sp); //site_sc: r765 -> r766: 1 -> 0 - rphase_hc(&b->olist, &R_INF, &b->hap, &b->self_read, &b->ovlp_read, &b->pidx, &b->v64, &buf0, 0, WINDOW_MAX_SIZE, b->self_read.length, 1/**, 0**/, i, (asm_opt.is_ont)?HPC_PL:0, asm_opt.is_ont, ((asm_opt.is_ont)?&(b->clist.chainDP):NULL), ((asm_opt.is_sc)?&(b->v8q):NULL), /**((asm_opt.is_sc)?&(b->v8t):NULL)**/&(b->v8t), (asm_opt.is_ont)?1:0, ((uint64_t)-1), 0, HC0_W, &b->v32, + rphase_hc(&b->olist, &R_INF, &b->hap, &b->self_read, &b->ovlp_read, &b->pidx, &b->v64, &buf0, 0, WINDOW_MAX_SIZE, b->self_read.length, 1/**, 0**/, i, (asm_opt.is_ont)?HPC_PL:0, asm_opt.is_ont, ((asm_opt.is_ont)?&(b->clist.chainDP):NULL), ((asm_opt.is_sc)?&(b->v8q):NULL), ((asm_opt.is_sc)?&(b->v8t):NULL)/**&(b->v8t)**/, (asm_opt.is_ont)?1:0, ((uint64_t)-1), 0, HC0_W, &b->v32, asm_opt.s_hap_cov, asm_opt.infor_cov, het_a, hom_a, asm_opt.polyploidy, -1.0); copy_asg_arr(b->sp, buf0); ///for debug indel @@ -4418,7 +4492,11 @@ static void worker_hap_ec(void *data, long i, int tid) push_nec_re(aux_o, &(scc.a[i])); - push_nec_re(aux_o, &(scb.a[i])); + // push_nec_re(aux_o, &(scb.a[i])); + if(asm_opt.dbg_bam) { + update_scb(&R_INF, &(scc.a[i]), &(scb.a[i]), &(b->v16), &b->self_read, &b->ovlp_read, &b->v64, &b->exz, i); + kv_resize(uint16_t, scb.a[i], b->v16.n); scb.a[i].n = b->v16.n; memcpy(scb.a[i].a, b->v16.a, b->v16.n*sizeof(*(scb.a[i].a))); + } // if((asm_opt.is_ont) && is_chemical_r_qual(&b->olist, &b->v64, qlen, 1, 16, &(b->v8q), i)/**(is_uncorrected_read(&b->olist, &b->v64, qlen, 1600))**/) { // // b->olist.length = 0; @@ -5258,7 +5336,11 @@ static void worker_hap_ec_hybrid(void *data, long i, int tid) copy_asg_arr(b->sp, buf0); copy_asg_arr(b->hap.snp_srt, buf1); push_nec_re(aux_o, &(scc.a[i])); - push_nec_re(aux_o, &(scb.a[i])); + // push_nec_re(aux_o, &(scb.a[i])); + if(asm_opt.dbg_bam) { + update_scb(&R_INF, &(scc.a[i]), &(scb.a[i]), &(b->v16), &b->self_read, &b->ovlp_read, &b->v64, &b->exz, i); + kv_resize(uint16_t, scb.a[i], b->v16.n); memcpy(scb.a[i].a, b->v16.a, b->v16.n); + } // if((asm_opt.is_ont) && is_chemical_r_qual(&b->olist, &b->v64, qlen, 1, 16, &(b->v8q), i)/**(is_uncorrected_read(&b->olist, &b->v64, qlen, 1600))**/) { // // b->olist.length = 0; @@ -5571,10 +5653,11 @@ uint32_t adjust_exact_match(asg16_v *in, int64_t xs0, int64_t xe0, int64_t ys0, return (*rxe) - (*rxs); } -uint32_t quick_exact_match(ma_hit_t *z, All_reads *rref, UC_Read* qu, UC_Read* tu, cc_v *sc) +uint32_t quick_exact_match(ma_hit_t *z, All_reads *rref, UC_Read* qu, UC_Read* tu, cc_v *sc, uint64_t qid) { uint64_t rts, rte, rqs, rqe, f = 0; int64_t ql, tl, qr, tr, qs, qe, ts, te; - + qs = z->qns>>32; qe = qs + z->cc; ts = z->bl; te = ts + z->cc; + z->qe = qe; z->ts = ts; z->te = te; z->qns = qid; z->qns <<= 32; z->qns |= qs; if((R_INF.is_syn == 1) && ((z->qns>>32) >= R_INF.tqn) && (z->tn < R_INF.tqn)) {///HiFi-to-ONT f = 1; } else if(adjust_exact_match(&(sc->a[z->tn]), z->ts, z->te, ((uint32_t)(z->qns)), z->qe, &rts, &rte, &rqs, &rqe, z->rev)) { @@ -5671,7 +5754,7 @@ uint64_t gen_hap_dc_cov(asg64_v *be, asg64_v *ba, ma_hit_t_alloc *paf, All_reads // (uint32_t)z->qns, z->qe, z->ts, z->te, z->el); // } - if((z->el) && (quick_exact_match(z, rref, qu, tu, sc))) { + if((z->el) && (quick_exact_match(z, rref, qu, tu, sc, rid))) { s = ((uint32_t)(z->qns)); e = z->qe; kv_push(uint64_t, (*be), (s<<1)); kv_push(uint64_t, (*be), (e<<1)|1); @@ -5892,7 +5975,7 @@ static void worker_update_dc_ec(void *data, long i, int tid) // (int32_t)Get_NAME_LENGTH(R_INF, z->tn), Get_NAME(R_INF, z->tn), z->tn, Get_READ_LENGTH(R_INF, (z->tn)), // z->ts, z->te, z->el?1:0); // } - if((z->el) && (quick_exact_match(z, &R_INF, &b->self_read, &b->ovlp_read, &scc))) { + if((z->el) && (quick_exact_match(z, &R_INF, &b->self_read, &b->ovlp_read, &scc, i))) { z->el = 1; b->cnt[0]++; } else { z->el = 0; b->cnt[1]++; @@ -5967,6 +6050,45 @@ static void worker_hap_post_rev(void *data, long i, int tid) for (k = 0; k < l; k++) a[l - k - 1] = b->v8q.a[k]; ha_compress_qual_bit(Get_QUAL(R_INF, i), a, l, sc_bn); } + + if(scb.a[i].n) { + uint16_t *z = NULL, c; + kl = scb.a[i].n>>1; + for (k = 0; k < kl; k++) { + nn = scb.a[i].a[k]; + scb.a[i].a[k] = scb.a[i].a[scb.a[i].n - k - 1]; + scb.a[i].a[scb.a[i].n - k - 1] = nn; + + z = &(scb.a[i].a[k]); + c = (*z)>>14; + if(c > 0) { + *z ^= (0x3u << 12); + if(c == 1) { + *z ^= (0x3u << 10); + } + } + + z = &(scb.a[i].a[scb.a[i].n - k - 1]); + c = (*z)>>14; + if(c > 0) { + *z ^= (0x3u << 12); + if(c == 1) { + *z ^= (0x3u << 10); + } + } + } + + if(((uint32_t)scb.a[i].n)&1) { + z = &(scb.a[i].a[k]); + c = (*z)>>14; + if(c > 0) { + *z ^= (0x3u << 12); + if(c == 1) { + *z ^= (0x3u << 10); + } + } + } + } } static void worker_hap_dc_ec_gen(void *data, long i, int tid) @@ -7849,9 +7971,21 @@ void gen_ori_seq0(char *tstr, uint64_t tl, UC_Read *qu, asg16_v *sc, uint64_t ri if(c != 3) tk += len; wq[1] = qk; wt[1] = tk; } - if(!(tk == tl)) { - fprintf(stderr, "[M::%s] rid::%lu, tk::%lu, tl::%lu\n", __func__, rid, tk, tl); - } + // if(!(tk == tl)) { + // if(rid == 8) { + // fprintf(stderr, "[M::%s] rid::%lu, tk::%lu, tl::%lu\n", __func__, rid, tk, tl); + // ck = 0; + // while (ck < sc->n) { + // wq[0] = qk; wt[0] = tk; + // ck = pop_trace_bp_f(sc, ck, &c, &bq, &bt, &len); + // if(c != 2) qk += len; + // if(c != 3) tk += len; + // wq[1] = qk; wt[1] = tk; + // fprintf(stderr, "[M::%s] c::%u, len::%u\n", __func__, c, len); + // } + // } + + // } assert(tk == tl); resize_UC_Read(qu, qk); qstr = qu->seq; qu->length = qk; @@ -8169,44 +8303,47 @@ void update_scb(All_reads *R_INF, asg16_v *scc, asg16_v *scb, asg16_v *scb_res, { char *qstr = NULL, *tstr = NULL; uint64_t ql = 0, tl = 0; uint64_t ck, qk, tk, k, wq[2], wt[2]; uint32_t len; uint16_t c, bq, bt; - gen_ori_seq0(qu->seq, qu->length, tu, scb, rid); ///tstr = tu->seq; tl = tu->length; + // if(!scb->n) push_trace_bp_f(scb, 0, (uint16_t)-1, (uint16_t)-1, qu->length, 0);///init + if(scb->n) { + gen_ori_seq0(qu->seq, qu->length, tu, scb, rid); ///tstr = tu->seq; tl = tu->length; - ck = qk = tk = 0; ql = qu->length; - while (ck < scc->n) { - wq[0] = qk; wt[0] = tk; - ck = pop_trace_bp_f(scc, ck, &c, &bq, &bt, &len); - if(c != 2) qk += len; - if(c != 3) tk += len; - wq[1] = qk; wt[1] = tk; - } - assert(qk == ql); - tl = tk; resize_UC_Read(qu, ql + tl); - qstr = qu->seq; tstr = qu->seq + ql; - - ck = 0; qk = tk = 0; - while (ck < scc->n) { - wq[0] = qk; wt[0] = tk; - ck = pop_trace_bp_f(scc, ck, &c, &bq, &bt, &len); - if(c != 2) qk += len; - if(c != 3) tk += len; - wq[1] = qk; wt[1] = tk; - // if(xk > (uint32_t)p->z.length) fprintf(stderr, "[M::%s] xk::%u, len::%u, c::%u, rid::%ld\n", __func__, xk, (uint32_t)p->z.length, c, i); - if(c == 0) { - memcpy(tstr + wt[0], qstr + wq[0], (wq[1]-wq[0])*sizeof((*qstr))); - } else if(c == 1 || c == 2) { - for (k = wt[0]; k < wt[1]; k++) tstr[k] = s_H[bt]; + ck = qk = tk = 0; ql = qu->length; + while (ck < scc->n) { + wq[0] = qk; wt[0] = tk; + ck = pop_trace_bp_f(scc, ck, &c, &bq, &bt, &len); + if(c != 2) qk += len; + if(c != 3) tk += len; + wq[1] = qk; wt[1] = tk; } - // if(i == 700) fprintf(stderr, "|%u%c(%c)(x::%u)(y::%u)", len, cm[c], ((c==1)||(c==2))?(cc[b]):('*'), wx[1], wy[1]); // s_H - } + assert(qk == ql); + tl = tk; resize_UC_Read(qu, ql + tl); + qstr = qu->seq; tstr = qu->seq + ql; - qstr = tstr; ql = tl; - tstr = tu->seq; tl = tu->length; + ck = 0; qk = tk = 0; + while (ck < scc->n) { + wq[0] = qk; wt[0] = tk; + ck = pop_trace_bp_f(scc, ck, &c, &bq, &bt, &len); + if(c != 2) qk += len; + if(c != 3) tk += len; + wq[1] = qk; wt[1] = tk; + // if(xk > (uint32_t)p->z.length) fprintf(stderr, "[M::%s] xk::%u, len::%u, c::%u, rid::%ld\n", __func__, xk, (uint32_t)p->z.length, c, i); + if(c == 0) { + memcpy(tstr + wt[0], qstr + wq[0], (wq[1]-wq[0])*sizeof((*qstr))); + } else if(c == 1 || c == 2) { + for (k = wt[0]; k < wt[1]; k++) tstr[k] = s_H[bt]; + } + // if(i == 700) fprintf(stderr, "|%u%c(%c)(x::%u)(y::%u)", len, cm[c], ((c==1)||(c==2))?(cc[b]):('*'), wx[1], wy[1]); // s_H + } - // fprintf(stderr, "\n[M::%s] ql::%lu, tl::%lu, rid::%lu\n", __func__, ql, tl, rid); + qstr = tstr; ql = tl; + tstr = tu->seq; tl = tu->length; - gen_updated_trace(scc, scb, scb_res, qstr, ql, tstr, tl, srt, exz, rid); + // fprintf(stderr, "\n[M::%s] ql::%lu, tl::%lu, rid::%lu\n", __func__, ql, tl, rid); - + gen_updated_trace(scc, scb, scb_res, qstr, ql, tstr, tl, srt, exz, rid); + } else { + kv_resize(uint16_t, *scb_res, scc->n); scb_res->n = scc->n; memcpy(scb_res->a, scc->a, scc->n*sizeof(*(scb_res->a))); + } ///debug // resize_UC_Read(tu, ql + tl); @@ -9181,23 +9318,35 @@ uint64_t cal_sec_ec_multiple(ec_ovec_buf_t *b, uint64_t n_thre, uint64_t n_a, in } -void write_ec_reads(const char *suffix_ou) +void write_ec_reads(const char *suffix_ou, cc_v *cvt, uint8_t is_rev) { - uint64_t k, strl; UC_Read qstr, tstr; char *nn = NULL, *str = NULL; + uint64_t k, strl, zk, zn; UC_Read qstr, tstr; char *nn = NULL, *str = NULL, zc; init_UC_Read(&qstr); init_UC_Read(&tstr); MALLOC(nn, strlen(suffix_ou) + strlen(asm_opt.output_file_name) + 36); sprintf(nn, "%s.%s", asm_opt.output_file_name, suffix_ou); + // fprintf(stderr, "[M::%s]\tnn::%s\n", __func__, nn); FILE *ou = fopen(nn, "w"); free(nn); for (k = 0; k < R_INF.total_reads; k++) { recover_UC_Read(&qstr, &R_INF, k); - if(scb.a) { - gen_ori_seq0(qstr.seq, qstr.length, &tstr, &(scb.a[k]), k); str = tstr.seq; strl = tstr.length; + if(cvt) { + gen_ori_seq0(qstr.seq, qstr.length, &tstr, &(cvt->a[k]), k); str = tstr.seq; strl = tstr.length; } else { str = qstr.seq; strl = qstr.length; } + if(is_rev) { + zn = strl>>1; + for (zk = 0; zk < zn; zk++) { + zc = str[strl-zk-1]; str[strl-zk-1] = RC_CHAR(str[zk]); str[zk] = RC_CHAR(zc); + } + if(strl&1) { + zc = str[strl-zk-1]; str[strl-zk-1] = RC_CHAR(str[zk]); str[zk] = RC_CHAR(zc); + } + } + + fwrite(">", 1, 1, ou); fwrite(Get_NAME(R_INF, k), 1, Get_NAME_LENGTH(R_INF, k), ou); fwrite("\n", 1, 1, ou); @@ -9228,9 +9377,19 @@ void prt_nel_ovlp(ma_hit_t_alloc *pa, uint64_t p_n) exit(1); } +void dbg_write_ec_reads(const char* i_cmd, uint64_t round, cc_v *cvt, uint8_t is_rev) +{ + char *nn = NULL; MALLOC(nn, strlen(i_cmd) + 64); + sprintf(nn, "r%lu.%s", round, i_cmd); + write_ec_reads((const char*)nn, cvt, is_rev); + free(nn); +} + void cal_ec_r(uint64_t n_thre, uint64_t round, uint64_t n_round, uint64_t n_a, uint64_t is_sv, uint64_t *tot_b, uint64_t *tot_e) { - // write_ec_reads("ec0.fa"); + // if(round == 0) { + // write_ec_reads("raw.ec.fa", NULL, 0); + // } // fprintf(stderr, "[M::%s]\tn_thre::%lu, round::%lu, n_round::%lu, n_a::%lu, is_sv::%lu\n", __func__, n_thre, round, n_round, n_a, is_sv); fprintf(stderr, "-0-[M::%s]\t# tqn::%lu, Ont base::%lu, # HiFi bases::%lu\n", __func__, R_INF.tqn, R_INF.tr[0], R_INF.tr[1]); @@ -9252,6 +9411,13 @@ void cal_ec_r(uint64_t n_thre, uint64_t round, uint64_t n_round, uint64_t n_a, u prt_dbg_stats(dbg_a, n_a, asm_opt.output_file_name, round, 0); free(dbg_a); dbg_a = NULL; } + + // uint32_t ck, len; uint16_t c; uint16_t bq, bt; + // ck = 0; fprintf(stderr, "\n-a-[M::%s] scb.a[8].n::%u\n", __func__, (uint32_t)scb.a[8].n); + // while (ck < scb.a[8].n) { + // ck = pop_trace_bp_f(&scb.a[8], ck, &c, &bq, &bt, &len); + // fprintf(stderr, "-a-[M::%s] c::%u, len::%u\n", __func__, c, len); + // } sl_ec_r(n_thre, n_a); for (k = 0; k < n_round; k++) { @@ -9269,6 +9435,8 @@ void cal_ec_r(uint64_t n_thre, uint64_t round, uint64_t n_round, uint64_t n_a, u // prt_nel_ovlp(R_INF.paf, n_a); // exit(1); + // dbg_write_ec_reads("ec12.fa", round, &scb, is_cr); + if((!is_sv) || (is_sv && is_cr)) { kt_for(n_thre, worker_hap_post_rev, b, n_a); } @@ -9282,7 +9450,8 @@ void cal_ec_r(uint64_t n_thre, uint64_t round, uint64_t n_round, uint64_t n_a, u fprintf(stderr, "-4-[M::%s]\t# tqn::%lu, Ont base::%lu, # HiFi bases::%lu\n", __func__, R_INF.tqn, R_INF.tr[0], R_INF.tr[1]); - // write_ec_reads("ec16.fa"); + // dbg_write_ec_reads("ec16.fa", round, &scb, !is_cr); + // exit(1); // uint64_t z; // for (z = 0; z < scc.n; z++) { @@ -9291,6 +9460,19 @@ void cal_ec_r(uint64_t n_thre, uint64_t round, uint64_t n_round, uint64_t n_a, u // } } +void gen_gfa_bam(ma_ug_t *ug) +{ + if(ha_flt_tab) { + ha_ft_destroy(ha_flt_tab); ha_flt_tab = NULL; + } + if(ha_idx) { + ha_pt_destroy(ha_idx); ha_idx = NULL; + } + + ha_flt_tab = ha_ft_ug_gen(&asm_opt, &(ug->u), 0, asm_opt.k_mer_length, asm_opt.mz_win, -1, -1); + // ha_idx = ha_pt_ug_gen(&asm_opt, ha_flt_tab, &(ug->u), hap_n); +} + void print_ov_dbg_paf(FILE *fp, char *ref_str, char *ref_id, int32_t ref_id_n, char *qry_str, char *qry_id, int32_t qry_id_n, uint64_t rs, uint64_t re, uint64_t rl, uint64_t qs, uint64_t qe, uint64_t ql, uint64_t rev, bit_extz_t *ez, char *ezh) { uint64_t ci = 0; uint16_t c; uint32_t cl; @@ -9371,12 +9553,12 @@ void cal_ov_r(uint64_t n_thre, uint64_t n_a, uint64_t new_idx) b = gen_ec_ovec_buf_t(n_thre); if(new_idx) { // kt_for(n_thre, worker_hap_dc_ec, b, n_a);///update overlaps - destroy_cc_v(&scc); destroy_cc_v(&scb); destroy_cc_v(&sca); + destroy_cc_v(&scc); if(!asm_opt.dbg_bam) destroy_cc_v(&scb); destroy_cc_v(&sca); ha_print_ovlp_stat_0(b, n_thre, n_a); } else { ha_print_ovlp_stat_1(b, n_thre, n_a); - destroy_cc_v(&scc); destroy_cc_v(&scb); destroy_cc_v(&sca); + destroy_cc_v(&scc); if(!asm_opt.dbg_bam) destroy_cc_v(&scb); destroy_cc_v(&sca); } destroy_ec_ovec_buf_t(b); @@ -9451,7 +9633,21 @@ void handle_chemical_arc(uint64_t n_thre, uint64_t n_a) uint8_t* gen_chemical_arc_rf(uint64_t n_thre, uint64_t n_a) { - ec_ovec_buf_t *b = NULL; uint64_t k, chem_n = 0; uint8_t *ra = NULL; + ec_ovec_buf_t *b = NULL; uint64_t k, chem_n = 0; uint8_t *ra = NULL; int64_t auto_chem_c = 0, hom_a; + + hom_a = asm_opt.hom_cov; + if(asm_opt.hom_global_coverage_set) hom_a = asm_opt.hom_global_coverage; + + auto_chem_c = 1 + (((hom_a>50)?(hom_a-50):(0))+5)/25; + fprintf(stderr, "[M::%s::auto] inferred chimeric threshold: %ld\n", __func__, auto_chem_c); + if(asm_opt.chemical_cov < 0) { + asm_opt.chemical_cov = auto_chem_c; + } else { + fprintf(stderr, "[M::%s::user] override requested chimeric threshold: %ld\n", __func__, asm_opt.chemical_cov); + } + if(asm_opt.chemical_cov < 0) asm_opt.chemical_cov = 0; + fprintf(stderr, "[M::%s::final] using chimeric threshold: %ld\n", __func__, asm_opt.chemical_cov); + b = gen_ec_ovec_buf_t(n_thre); for (k = 0; k < n_thre; ++k) { b->a[k].cnt[0] = 0; b->a[k].cnt[1] = 0; diff --git a/ecovlp.h b/ecovlp.h index 3e1fad8..7ddf18a 100644 --- a/ecovlp.h +++ b/ecovlp.h @@ -16,5 +16,8 @@ void handle_chemical_r(uint64_t n_thre, uint64_t n_a); void handle_chemical_arc(uint64_t n_thre, uint64_t n_a); uint8_t* gen_chemical_arc_rf(uint64_t n_thre, uint64_t n_a); void cal_ec_r_dbg(uint64_t n_thre, uint64_t n_a); +void write_ec_reads(const char *suffix_ou, cc_v *cvt, uint8_t is_rev); +void destroy_cc_v(cc_v *z); +void gen_gfa_bam(ma_ug_t *ug); #endif \ No newline at end of file diff --git a/htab.cpp b/htab.cpp index af9802b..e44fbf6 100644 --- a/htab.cpp +++ b/htab.cpp @@ -1138,12 +1138,35 @@ void *ha_ft_ul_gen(const hifiasm_opt_t *asm_opt, ma_utg_v *us, int k, int w, int return (void*)flt_tab; } -void *ha_ft_ug_gen(const hifiasm_opt_t *asm_opt, ma_utg_v *us, int is_HPC, int k, int w, int min_freq, int max_freq) +// int64_t ha_ct_ug_cutoff(ha_ct_t *h, int64_t num_thre, double cut_rate, double cut_an_rate) +// { +// int64_t cnt[YAK_N_COUNTS], k, tot_n = 0, tot_cutn = 0, tot_an = 0, tot_cut_an = 0, tot_k = 0, tot_ak = 0; +// ha_ct_hist(h, cnt, num_thre); +// for (k = tot_n = tot_an = 0; k < YAK_N_COUNTS; k++) { +// tot_n += cnt[k]; tot_an += cnt[k]*k; +// } + +// tot_cutn = tot_n - (tot_n*cut_rate); +// tot_cut_an = tot_an - (tot_an*cut_an_rate); + + + +// } + +void *ha_ft_ug_gen(hifiasm_opt_t *asm_opt, ma_utg_v *us, int is_HPC, int k, int w, int min_freq, int max_freq) { - yak_ft_t *flt_tab; - ha_ct_t *h; + yak_ft_t *flt_tab; + ha_ct_t *h; ///int32_t b0; + // if(min_freq < 0 || max_freq < 0) { + // b0 = asm_opt->bf_shift; asm_opt->bf_shift = 0; + // } ///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, 0); + // if(min_freq < 0 || max_freq < 0) { + // asm_opt->bf_shift = b0; + + // } + 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, YAK_MAX_COUNT); ha_ct_destroy(h); diff --git a/htab.h b/htab.h index 425dd07..d7ca3d6 100644 --- a/htab.h +++ b/htab.h @@ -73,7 +73,7 @@ extern ha_pt_t *ha_idx_hp; extern void *ha_ct_table; void *ha_ft_ul_gen(const hifiasm_opt_t *asm_opt, ma_utg_v *us, int k, int w, int cutoff); -void *ha_ft_ug_gen(const hifiasm_opt_t *asm_opt, ma_utg_v *us, int is_HPC, int k, int w, int min_freq, int max_freq); +void *ha_ft_ug_gen(hifiasm_opt_t *asm_opt, ma_utg_v *us, int is_HPC, int k, int w, int min_freq, int max_freq); void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov, int is_hp_mode, int read_from_store); int32_t ha_ft_cnt(const void *hh, uint64_t y); void ha_ft_destroy(void *h); diff --git a/main.cpp b/main.cpp index a6270e2..3ce4a0c 100644 --- a/main.cpp +++ b/main.cpp @@ -62,6 +62,32 @@ int main(int argc, char *argv[]) // ed_band_cal_global((char*)"ACTTTTTT", 8, (char*)"AATTTT", 6, 3), // ed_band_cal_global_128bit((char*)"ACTTTTTT", 8, (char*)"AATTTT", 6, 3)); // exit(1); + + int8_t simd_auto = 0; + #if defined(__x86_64__) || defined(__i386__) + __builtin_cpu_init(); + if (__builtin_cpu_supports("avx512f")) simd_auto = 2; + else if (__builtin_cpu_supports("avx2")) simd_auto = 1; + else simd_auto = 0; + #endif + + if (simd_auto) { + fprintf(stderr, "[M::%s::auto] detected CPU support for %s\n", __func__, (simd_auto == 2) ? "AVX-512" : "AVX2"); + } else { + fprintf(stderr, "[M::%s::auto] no supported SIMD extension detected; falling back to non-SIMD\n", __func__); + } + + if (asm_opt.simd_mm == 0 || asm_opt.simd_mm == 1 || asm_opt.simd_mm == 2) { + fprintf(stderr, "[M::%s::user] user requested %s mode\n", __func__, (asm_opt.simd_mm == 2) ? "AVX-512" : ((asm_opt.simd_mm == 1) ? "AVX2" : "non-SIMD")); + if(asm_opt.simd_mm > simd_auto) asm_opt.simd_mm = simd_auto; + } else { + asm_opt.simd_mm = simd_auto; + if(asm_opt.simd_mm >= 1) asm_opt.simd_mm = 1; ///use avx2 rather than avx512, looks like avx512 still has issues right now + } + + fprintf(stderr, "[M::%s::final] using %s mode\n", __func__, (asm_opt.simd_mm == 2) ? "AVX-512" : ((asm_opt.simd_mm == 1) ? "AVX2" : "non-SIMD")); + + if(asm_opt.sec_in) ret = ha_assemble_pair(); else if(asm_opt.dbg_ovec_cal) ret = ha_ec_dbg(); else ret = ha_assemble(); @@ -71,6 +97,6 @@ int main(int argc, char *argv[]) fprintf(stderr, "[M::%s] CMD:", __func__); for (i = 0; i < argc; ++i) fprintf(stderr, " %s", argv[i]); - fprintf(stderr, "\n[M::%s] Real time: %.3f sec; CPU: %.3f sec; Peak RSS: %.3f GB\n", __func__, yak_realtime(), yak_cputime(), yak_peakrss_in_gb()); + fprintf(stderr, "\n[M::%s] Real time: %.3f sec; CPU: %.3f sec; Peak RSS: %.3f GB; SIMD: %u\n", __func__, yak_realtime(), yak_cputime(), yak_peakrss_in_gb(), asm_opt.simd_mm); return ret; }