xxay/hifiasm
1
0
Fork 0
mirror of https://github.com/chhylp123/hifiasm.git synced 2026-10-01 13:28:12 +08:00
Code Issues Packages Projects Releases Wiki Activity

fall back to avx2

Browse Source
This commit is contained in:
chhylp123
2026-04-01 22:28:54 -04:00
parent a137a6b06f
commit 815d709fb9
16 changed files with 1122 additions and 390 deletions
Download Patch File Download Diff File Expand all files Collapse all files
Assembly.cpp
+4 -2
View File
@@ -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;
}
CommandLines.cpp
+13 -3
View File
@@ -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 <INT minimizers [%ld]\n", asm_opt->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 <INT in size in contig graphs [%lld]\n", asm_opt->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);
}
CommandLines.h
+6 -2
View File
@@ -5,7 +5,7 @@
#include <pthread.h>
#include <stdint.h>
#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;
Correct.cpp
+264 -46
View File
@@ -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);
}
}
}
Levenshtein_avx2.cpp
+162
View File
@@ -0,0 +1,162 @@
#include "Levenshtein_distance.h"
#include <immintrin.h>
#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)<<AVX_GS2)-1; int32_t bd, ibd[AVX_GS2], i, last_high = (thre<<1), tn0 = tn - 1, cut = thre+last_high;
lone = _mm256_set1_epi64x(1);
VP = _mm256_setzero_si256();
VN_mm = Peq_mm[0];
VN_mm[0] = (((Word)1)<<(abs_diag_a[0]))-1; VN_mm[1] = (((Word)1)<<(abs_diag_a[1]))-1; VN_mm[2] = (((Word)1)<<(abs_diag_a[2]))-1; VN_mm[3] = (((Word)1)<<(abs_diag_a[3]))-1;
VN = _mm256_loadu_si256((const __m256i *)VN_mm);
VN_mm[0] = abs_diag_a[0]; VN_mm[1] = abs_diag_a[1]; VN_mm[2] = abs_diag_a[2]; VN_mm[3] = abs_diag_a[3];
E = _mm256_loadu_si256((const __m256i *)VN_mm);
memset(VN_mm, 0, (sizeof((*VN_mm))*AVX_GS2)); VN_mm = NULL;///reset
init_simd_ed4(psa, pna, thre, abs_diag_a, r_err, r_pe, 0, tn, cut, bd, i, mm, Peq_mm, lz, ibd);
init_simd_ed4(psa, pna, thre, abs_diag_a, r_err, r_pe, 1, tn, cut, bd, i, mm, Peq_mm, lz, ibd);
init_simd_ed4(psa, pna, thre, abs_diag_a, r_err, r_pe, 2, tn, cut, bd, i, mm, Peq_mm, lz, ibd);
init_simd_ed4(psa, pna, thre, abs_diag_a, r_err, r_pe, 3, tn, cut, bd, i, mm, Peq_mm, lz, ibd);
ht &= lz;
if(ht == 0) return;
Peq[0] = _mm256_loadu_si256((const __m256i *)Peq_mm[0]);
Peq[1] = _mm256_loadu_si256((const __m256i *)Peq_mm[1]);
Peq[2] = _mm256_loadu_si256((const __m256i *)Peq_mm[2]);
Peq[3] = _mm256_loadu_si256((const __m256i *)Peq_mm[3]);
Peq[4] = _mm256_setzero_si256();
C = _mm256_set1_epi64x(cut + 1);///_mm512_set1_epi64x(cut);
mm = ((Word)1 << (thre<<1));///for the incoming char/last char**
mmk[0] = _mm256_set_epi64x(0, 0, 0, mm);
mmk[1] = _mm256_set_epi64x(0, 0, mm, 0);
mmk[2] = _mm256_set_epi64x(0, mm, 0, 0);
mmk[3] = _mm256_set_epi64x(mm, 0, 0, 0);
i = 0;
while (i < tn0) {
ed_core_64x4(Peq[seq_nt4_table[(uint8_t)tstr[i]]], VP, VN, X, D0, HN, HP);
E = _mm256_add_epi64(_mm256_xor_si256(lone, _mm256_and_si256(D0, lone)), E);
// ht = _mm512_cmple_epi64_mask(E, C);
ht =_mm256_movemask_pd(_mm256_castsi256_pd(_mm256_cmpgt_epi64(C, E)));
ht &= lz;
if(ht == 0) return;
Peq[0] = _mm256_srli_epi64(Peq[0], 1);
Peq[1] = _mm256_srli_epi64(Peq[1], 1);
Peq[2] = _mm256_srli_epi64(Peq[2], 1);
Peq[3] = _mm256_srli_epi64(Peq[3], 1);
i++; ///c = 4;
ed_core_upx4(Peq, psa, pna, ibd, ht, c, mmk, 0);
ed_core_upx4(Peq, psa, pna, ibd, ht, c, mmk, 1);
ed_core_upx4(Peq, psa, pna, ibd, ht, c, mmk, 2);
ed_core_upx4(Peq, psa, pna, ibd, ht, c, mmk, 3);
}
ed_core_64x4(Peq[seq_nt4_table[(uint8_t)tstr[i]]], VP, VN, X, D0, HN, HP);
E = _mm256_add_epi64(_mm256_xor_si256(lone, _mm256_and_si256(D0, lone)), E);
// ht = _mm512_cmple_epi64_mask(E, C);
ht =_mm256_movemask_pd(_mm256_castsi256_pd(_mm256_cmpgt_epi64(C, E)));
ht &= lz;
if(ht == 0) return;
int32_t st, ai; int64_t err_mm[AVX_GS2], uge;
VN_mm = Peq_mm[0]; VP_mm = Peq_mm[1];
_mm256_storeu_si256((__m256i *)VN_mm, VN); _mm256_storeu_si256((__m256i *)VP_mm, VP); _mm256_storeu_si256((__m256i *)err_mm, E);
ED_TAIL_LANE(0);
ED_TAIL_LANE(1);
ED_TAIL_LANE(2);
ED_TAIL_LANE(3);
}
Levenshtein_avx512.cpp
+276
View File
@@ -0,0 +1,276 @@
#include "Levenshtein_distance.h"
#include <immintrin.h>
#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)<<AVX_GS)-1, mtf, mt, mbest; int32_t bd, ibd[AVX_GS], i, last_high = (thre<<1), tn0 = tn - 1, cut = thre+last_high;
lone = _mm512_set1_epi64(1);
VP = _mm512_setzero_si512();
VN_mm = Peq_mm[0];
VN_mm[0] = (((Word)1)<<(abs_diag_a[0]))-1; VN_mm[1] = (((Word)1)<<(abs_diag_a[1]))-1; VN_mm[2] = (((Word)1)<<(abs_diag_a[2]))-1; VN_mm[3] = (((Word)1)<<(abs_diag_a[3]))-1;
VN_mm[4] = (((Word)1)<<(abs_diag_a[4]))-1; VN_mm[5] = (((Word)1)<<(abs_diag_a[5]))-1; VN_mm[6] = (((Word)1)<<(abs_diag_a[6]))-1; VN_mm[7] = (((Word)1)<<(abs_diag_a[7]))-1;
VN = _mm512_loadu_si512(VN_mm);
VN_mm[0] = abs_diag_a[0]; VN_mm[1] = abs_diag_a[1]; VN_mm[2] = abs_diag_a[2]; VN_mm[3] = abs_diag_a[3];
VN_mm[4] = abs_diag_a[4]; VN_mm[5] = abs_diag_a[5]; VN_mm[6] = abs_diag_a[6]; VN_mm[7] = abs_diag_a[7];
E = _mm512_loadu_si512(VN_mm);
memset(VN_mm, 0, (sizeof((*VN_mm))*AVX_GS)); VN_mm = NULL;///reset
init_simd_ed(psa, pna, thre, abs_diag_a, r_err, r_pe, 0, tn, cut, bd, i, mm, Peq_mm, lz, ibd);
init_simd_ed(psa, pna, thre, abs_diag_a, r_err, r_pe, 1, tn, cut, bd, i, mm, Peq_mm, lz, ibd);
init_simd_ed(psa, pna, thre, abs_diag_a, r_err, r_pe, 2, tn, cut, bd, i, mm, Peq_mm, lz, ibd);
init_simd_ed(psa, pna, thre, abs_diag_a, r_err, r_pe, 3, tn, cut, bd, i, mm, Peq_mm, lz, ibd);
init_simd_ed(psa, pna, thre, abs_diag_a, r_err, r_pe, 4, tn, cut, bd, i, mm, Peq_mm, lz, ibd);
init_simd_ed(psa, pna, thre, abs_diag_a, r_err, r_pe, 5, tn, cut, bd, i, mm, Peq_mm, lz, ibd);
init_simd_ed(psa, pna, thre, abs_diag_a, r_err, r_pe, 6, tn, cut, bd, i, mm, Peq_mm, lz, ibd);
init_simd_ed(psa, pna, thre, abs_diag_a, r_err, r_pe, 7, tn, cut, bd, i, mm, Peq_mm, lz, ibd);
ht &= lz;
if(ht == 0) return;
Peq[0] = _mm512_loadu_si512(Peq_mm[0]);
Peq[1] = _mm512_loadu_si512(Peq_mm[1]);
Peq[2] = _mm512_loadu_si512(Peq_mm[2]);
Peq[3] = _mm512_loadu_si512(Peq_mm[3]);
Peq[4] = _mm512_setzero_si512();
C = _mm512_set1_epi64(cut);
mm = ((Word)1 << (thre<<1));///for the incoming char/last char**
mmk[0] = _mm512_mask_set1_epi64(VP, 1, mm);
mmk[1] = _mm512_mask_set1_epi64(VP, 2, mm);
mmk[2] = _mm512_mask_set1_epi64(VP, 4, mm);
mmk[3] = _mm512_mask_set1_epi64(VP, 8, mm);
mmk[4] = _mm512_mask_set1_epi64(VP, 16, mm);
mmk[5] = _mm512_mask_set1_epi64(VP, 32, mm);
mmk[6] = _mm512_mask_set1_epi64(VP, 64, mm);
mmk[7] = _mm512_mask_set1_epi64(VP, 128, mm);
i = 0;
while (i < tn0) {
ed_core_64x8(Peq[seq_nt4_table[(uint8_t)tstr[i]]], VP, VN, X, D0, HN, HP);
E = _mm512_add_epi64(_mm512_xor_si512(lone, _mm512_and_si512(D0, lone)), E);
ht = _mm512_cmple_epi64_mask(E, C);
ht &= lz;
if(ht == 0) return;
Peq[0] = _mm512_srli_epi64(Peq[0], 1);
Peq[1] = _mm512_srli_epi64(Peq[1], 1);
Peq[2] = _mm512_srli_epi64(Peq[2], 1);
Peq[3] = _mm512_srli_epi64(Peq[3], 1);
i++; ///c = 4;
ed_core_upx8(Peq, psa, pna, ibd, ht, c, mmk, 0);
ed_core_upx8(Peq, psa, pna, ibd, ht, c, mmk, 1);
ed_core_upx8(Peq, psa, pna, ibd, ht, c, mmk, 2);
ed_core_upx8(Peq, psa, pna, ibd, ht, c, mmk, 3);
ed_core_upx8(Peq, psa, pna, ibd, ht, c, mmk, 4);
ed_core_upx8(Peq, psa, pna, ibd, ht, c, mmk, 5);
ed_core_upx8(Peq, psa, pna, ibd, ht, c, mmk, 6);
ed_core_upx8(Peq, psa, pna, ibd, ht, c, mmk, 7);
}
ed_core_64x8(Peq[seq_nt4_table[(uint8_t)tstr[i]]], VP, VN, X, D0, HN, HP);
E = _mm512_add_epi64(_mm512_xor_si512(lone, _mm512_and_si512(D0, lone)), E);
ht = _mm512_cmple_epi64_mask(E, C);
ht &= lz;
if(ht == 0) return;
// site = tn - 1 - abs_diag;/**up bound**/
// ai = pn - tn + abs_diag; /**in most cases, ai = (thre<<1)**/
int32_t st[AVX_GS] = {tn-1, tn-1, tn-1, tn-1, tn-1, tn-1, tn-1, tn-1};
int32_t ai[AVX_GS] = {-tn, -tn, -tn, -tn, -tn, -tn, -tn, -tn};
int32_t k[AVX_GS] = {0}; i = -1; bd = 0;
int64_t err_mm[AVX_GS], uge_mm[AVX_GS] = {INT32_MAX, INT32_MAX, INT32_MAX, INT32_MAX, INT32_MAX, INT32_MAX, INT32_MAX, INT32_MAX};
VN_mm = Peq_mm[0]; VP_mm = Peq_mm[1];
_mm512_storeu_si512(VN_mm, VN); _mm512_storeu_si512(VP_mm, VP); _mm512_storeu_si512(err_mm, E);
ed_tail_upx8(ht, 0, st, ai, pna, abs_diag_a, k, err_mm, VP_mm, VN_mm, thre, r_err, r_pe, bd, i);
ed_tail_upx8(ht, 1, st, ai, pna, abs_diag_a, k, err_mm, VP_mm, VN_mm, thre, r_err, r_pe, bd, i);
ed_tail_upx8(ht, 2, st, ai, pna, abs_diag_a, k, err_mm, VP_mm, VN_mm, thre, r_err, r_pe, bd, i);
ed_tail_upx8(ht, 3, st, ai, pna, abs_diag_a, k, err_mm, VP_mm, VN_mm, thre, r_err, r_pe, bd, i);
ed_tail_upx8(ht, 4, st, ai, pna, abs_diag_a, k, err_mm, VP_mm, VN_mm, thre, r_err, r_pe, bd, i);
ed_tail_upx8(ht, 5, st, ai, pna, abs_diag_a, k, err_mm, VP_mm, VN_mm, thre, r_err, r_pe, bd, i);
ed_tail_upx8(ht, 6, st, ai, pna, abs_diag_a, k, err_mm, VP_mm, VN_mm, thre, r_err, r_pe, bd, i);
ed_tail_upx8(ht, 7, st, ai, pna, abs_diag_a, k, err_mm, VP_mm, VN_mm, thre, r_err, r_pe, bd, i);
if(bd <= 0) return;
if(bd > 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;
}
}
Levenshtein_distance.h
+3 -268
View File
@@ -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)<<AVX_GS)-1, mtf, mt, mbest; int32_t bd, ibd[AVX_GS], i, last_high = (thre<<1), tn0 = tn - 1, cut = thre+last_high;
lone = _mm512_set1_epi64(1);
VP = _mm512_setzero_si512();
VN_mm = Peq_mm[0];
VN_mm[0] = (((Word)1)<<(abs_diag_a[0]))-1; VN_mm[1] = (((Word)1)<<(abs_diag_a[1]))-1; VN_mm[2] = (((Word)1)<<(abs_diag_a[2]))-1; VN_mm[3] = (((Word)1)<<(abs_diag_a[3]))-1;
VN_mm[4] = (((Word)1)<<(abs_diag_a[4]))-1; VN_mm[5] = (((Word)1)<<(abs_diag_a[5]))-1; VN_mm[6] = (((Word)1)<<(abs_diag_a[6]))-1; VN_mm[7] = (((Word)1)<<(abs_diag_a[7]))-1;
VN = _mm512_loadu_si512(VN_mm);
VN_mm[0] = abs_diag_a[0]; VN_mm[1] = abs_diag_a[1]; VN_mm[2] = abs_diag_a[2]; VN_mm[3] = abs_diag_a[3];
VN_mm[4] = abs_diag_a[4]; VN_mm[5] = abs_diag_a[5]; VN_mm[6] = abs_diag_a[6]; VN_mm[7] = abs_diag_a[7];
E = _mm512_loadu_si512(VN_mm);
memset(VN_mm, 0, (sizeof((*VN_mm))*AVX_GS)); VN_mm = NULL;///reset
init_simd_ed(psa, pna, thre, abs_diag_a, r_err, r_pe, 0, tn, cut, bd, i, mm, Peq_mm, lz, ibd);
init_simd_ed(psa, pna, thre, abs_diag_a, r_err, r_pe, 1, tn, cut, bd, i, mm, Peq_mm, lz, ibd);
init_simd_ed(psa, pna, thre, abs_diag_a, r_err, r_pe, 2, tn, cut, bd, i, mm, Peq_mm, lz, ibd);
init_simd_ed(psa, pna, thre, abs_diag_a, r_err, r_pe, 3, tn, cut, bd, i, mm, Peq_mm, lz, ibd);
init_simd_ed(psa, pna, thre, abs_diag_a, r_err, r_pe, 4, tn, cut, bd, i, mm, Peq_mm, lz, ibd);
init_simd_ed(psa, pna, thre, abs_diag_a, r_err, r_pe, 5, tn, cut, bd, i, mm, Peq_mm, lz, ibd);
init_simd_ed(psa, pna, thre, abs_diag_a, r_err, r_pe, 6, tn, cut, bd, i, mm, Peq_mm, lz, ibd);
init_simd_ed(psa, pna, thre, abs_diag_a, r_err, r_pe, 7, tn, cut, bd, i, mm, Peq_mm, lz, ibd);
ht &= lz;
if(ht == 0) return;
Peq[0] = _mm512_loadu_si512(Peq_mm[0]);
Peq[1] = _mm512_loadu_si512(Peq_mm[1]);
Peq[2] = _mm512_loadu_si512(Peq_mm[2]);
Peq[3] = _mm512_loadu_si512(Peq_mm[3]);
Peq[4] = _mm512_setzero_si512();
C = _mm512_set1_epi64(cut);
mm = ((Word)1 << (thre<<1));///for the incoming char/last char**
mmk[0] = _mm512_mask_set1_epi64(VP, 1, mm);
mmk[1] = _mm512_mask_set1_epi64(VP, 2, mm);
mmk[2] = _mm512_mask_set1_epi64(VP, 4, mm);
mmk[3] = _mm512_mask_set1_epi64(VP, 8, mm);
mmk[4] = _mm512_mask_set1_epi64(VP, 16, mm);
mmk[5] = _mm512_mask_set1_epi64(VP, 32, mm);
mmk[6] = _mm512_mask_set1_epi64(VP, 64, mm);
mmk[7] = _mm512_mask_set1_epi64(VP, 128, mm);
i = 0;
while (i < tn0) {
ed_core_64x8(Peq[seq_nt4_table[(uint8_t)tstr[i]]], VP, VN, X, D0, HN, HP);
E = _mm512_add_epi64(_mm512_xor_si512(lone, _mm512_and_si512(D0, lone)), E);
ht = _mm512_cmple_epi64_mask(E, C);
ht &= lz;
if(ht == 0) return;
Peq[0] = _mm512_srli_epi64(Peq[0], 1);
Peq[1] = _mm512_srli_epi64(Peq[1], 1);
Peq[2] = _mm512_srli_epi64(Peq[2], 1);
Peq[3] = _mm512_srli_epi64(Peq[3], 1);
i++; ///c = 4;
ed_core_upx8(Peq, psa, pna, ibd, ht, c, mmk, 0);
ed_core_upx8(Peq, psa, pna, ibd, ht, c, mmk, 1);
ed_core_upx8(Peq, psa, pna, ibd, ht, c, mmk, 2);
ed_core_upx8(Peq, psa, pna, ibd, ht, c, mmk, 3);
ed_core_upx8(Peq, psa, pna, ibd, ht, c, mmk, 4);
ed_core_upx8(Peq, psa, pna, ibd, ht, c, mmk, 5);
ed_core_upx8(Peq, psa, pna, ibd, ht, c, mmk, 6);
ed_core_upx8(Peq, psa, pna, ibd, ht, c, mmk, 7);
}
ed_core_64x8(Peq[seq_nt4_table[(uint8_t)tstr[i]]], VP, VN, X, D0, HN, HP);
E = _mm512_add_epi64(_mm512_xor_si512(lone, _mm512_and_si512(D0, lone)), E);
ht = _mm512_cmple_epi64_mask(E, C);
ht &= lz;
if(ht == 0) return;
// site = tn - 1 - abs_diag;/**up bound**/
// ai = pn - tn + abs_diag; /**in most cases, ai = (thre<<1)**/
int32_t st[AVX_GS] = {tn-1, tn-1, tn-1, tn-1, tn-1, tn-1, tn-1, tn-1};
int32_t ai[AVX_GS] = {-tn, -tn, -tn, -tn, -tn, -tn, -tn, -tn};
int32_t k[AVX_GS] = {0}; i = -1; bd = 0;
int64_t err_mm[AVX_GS], uge_mm[AVX_GS] = {INT32_MAX, INT32_MAX, INT32_MAX, INT32_MAX, INT32_MAX, INT32_MAX, INT32_MAX, INT32_MAX};
VN_mm = Peq_mm[0]; VP_mm = Peq_mm[1];
_mm512_storeu_si512(VN_mm, VN); _mm512_storeu_si512(VP_mm, VP); _mm512_storeu_si512(err_mm, E);
ed_tail_upx8(ht, 0, st, ai, pna, abs_diag_a, k, err_mm, VP_mm, VN_mm, thre, r_err, r_pe, bd, i);
ed_tail_upx8(ht, 1, st, ai, pna, abs_diag_a, k, err_mm, VP_mm, VN_mm, thre, r_err, r_pe, bd, i);
ed_tail_upx8(ht, 2, st, ai, pna, abs_diag_a, k, err_mm, VP_mm, VN_mm, thre, r_err, r_pe, bd, i);
ed_tail_upx8(ht, 3, st, ai, pna, abs_diag_a, k, err_mm, VP_mm, VN_mm, thre, r_err, r_pe, bd, i);
ed_tail_upx8(ht, 4, st, ai, pna, abs_diag_a, k, err_mm, VP_mm, VN_mm, thre, r_err, r_pe, bd, i);
ed_tail_upx8(ht, 5, st, ai, pna, abs_diag_a, k, err_mm, VP_mm, VN_mm, thre, r_err, r_pe, bd, i);
ed_tail_upx8(ht, 6, st, ai, pna, abs_diag_a, k, err_mm, VP_mm, VN_mm, thre, r_err, r_pe, bd, i);
ed_tail_upx8(ht, 7, st, ai, pna, abs_diag_a, k, err_mm, VP_mm, VN_mm, thre, r_err, r_pe, bd, i);
if(bd <= 0) return;
if(bd > 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)
{
Makefile
+12 -3
View File
@@ -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)
Overlaps.cpp
+15 -2
View File
@@ -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) {
Process_Read.cpp
+49
View File
@@ -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;
Process_Read.h
+6
View File
@@ -8,6 +8,7 @@
#include <zlib.h>
#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);
ecovlp.cpp
+255 -59
View File
@@ -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;
ecovlp.h
+3
View File
@@ -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
htab.cpp
+26 -3
View File
@@ -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);
htab.h
+1 -1
View File
@@ -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);
main.cpp
+27 -1
View File
@@ -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;
}