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

simd version 0

Browse Source
This commit is contained in:
chhylp123
2026-03-24 17:54:36 -04:00
parent 9a2884f518
commit aa4d12b1c6
6 changed files with 671 additions and 27 deletions
Download Patch File Download Diff File Expand all files Collapse all files
CommandLines.h
+1 -1
View File
@@ -5,7 +5,7 @@
#include <pthread.h>
#include <stdint.h>
#define HA_VERSION "0.25.0-r869"
#define HA_VERSION "0.25.0-r873"
#define VERBOSE 0
Correct.cpp
+390 -20
View File
@@ -15123,6 +15123,9 @@ uint32_t push_hc_wlst_exz(const ul_idx_t *uref, hpc_t *hpc_g, All_reads *rref, o
// if((!force_aln) && (!simi_pass(ovl, aln, 0, ovlp_cut, &e_rate)) && (!simi_pass(ovl, aln, sec_check, ovlp_cut, NULL))) {
if((!force_aln) && (!pass_qovlp(ovl, aln, ovlp_cut))) {
rf = 0;
// if(ol->x_id == 1 && ol->y_id == 24) {
// fprintf(stderr, "-a-[M::%s]\tqid::%u\ttid::%uovl::%ld\taln::%ld\tovlp_cut::%f\n", __func__, ol->x_id, ol->y_id, ovl, aln, ovlp_cut);
// }
if(!is_srt) {
kv_push(window_list, ol->w_list, p);
return 0;
@@ -15307,18 +15310,21 @@ uint32_t align_hc_ed_post_extz(overlap_region *z, All_reads *rref, char* qstr, c
// tot_b0 = (*tot_b);
// }
// if(z->x_id == 1 && z->y_id == 24) {
// fprintf(stderr, "-z-[M::%s]\twq::[%ld,%ld)\tzql::%u\talign_length::%u\tthre::%ld\taln::%u\n", __func__, q_s, q_e+1, z->x_pos_e+1-z->x_pos_s, z->align_length, thre, (is_align(*exz)));
// }
// if(z->x_id == 5569 && z->y_id == 5557 && q_s == 10075 && q_e == 10849) {
// fprintf(stderr, "\n[M::%s::semi::t_s->%ld::t_pri_l->%ld::aux_beg->%ld::aux_end->%ld::thre->%ld] exz->ps::%d, exz->pe::%d, exz->ts::%d, exz->te::%d, exz->err::%d, exz->cigar.n::%d, thre::%ld\n",
// __func__, t_s, t_pri_l, aux_beg, aux_end, thre, exz->ps, exz->pe, exz->ts, exz->te, exz->err, (int32_t)exz->cigar.n, thre);
// fprintf(stderr, "[tstr::len->%ld] %.*s\n", t_l, (int32_t)t_l, t_string);
// fprintf(stderr, "[qstr::len->%ld] %.*s\n", q_l, (int32_t)q_l, q_string);
// if(z->x_id == 1 && z->y_id == 24) {
// fprintf(stderr, "\n[M::%s::semi::t_s->%ld::t_pri_l->%ld::aux_beg->%ld::aux_end->%ld::thre->%ld] exz->ps::%d, exz->pe::%d, exz->ts::%d, exz->te::%d, exz->err::%d, exz->cigar.n::%d, thre::%ld, e_rate::%f, q::[%ld,%ld), is_aln::%u\n",
// __func__, t_s, t_pri_l, aux_beg, aux_end, thre, exz->ps, exz->pe, exz->ts, exz->te, exz->err, (int32_t)exz->cigar.n, thre, e_rate, q_s, q_e + 1, is_align(*exz));
// fprintf(stderr, "[tstr::len->%ld] %.*s\n", t_l, (int32_t)t_l, t_string);
// fprintf(stderr, "[qstr::len->%ld] %.*s\n", q_l, (int32_t)q_l, q_string);
// }
if (is_align(*exz)) {
// ed_band_cal_semi_64_w(t_string, aln_l, q_string, q_l, thre, exz);
// assert(exz->err <= exz->thre);
// if(z->x_id == 19350 && z->y_id == 19324) {
// fprintf(stderr, "+[M::%s]\tq::[%ld,%ld)\tt::[%ld,%ld)\texz->err::%d\n", __func__, q_s, q_e + 1, t_s, t_s + exz->pe + 1, exz->err);
// if(z->x_id == 63 && z->y_id == 9) {
// fprintf(stderr, "\n+[M::%s]\tq::[%ld,%ld)\tt::[%ld,%ld)\texz->err::%d\tzaln::%u\n", __func__, q_s, q_e + 1, t_s, t_s + exz->pe + 1, exz->err, z->align_length);
// }
///t_s do not have aux_beg, while t_s + t_end (aka, te) has
if(!push_hc_wlst_exz(NULL, NULL, rref, z, qstr, tstr, exz, THRESHOLD_MAX_SIZE, q_s, q_e, t_s, t_s + exz->pe,
@@ -15326,6 +15332,9 @@ uint32_t align_hc_ed_post_extz(overlap_region *z, All_reads *rref, char* qstr, c
// if(z->y_id == 234) fprintf(stderr, "-b-[M::%s] tid::%u(%c)\tq::[%u,%u)\tt::[%u,%u)\n", __func__, z->y_id, "+-"[z->y_pos_strand], z->x_pos_s, z->x_pos_e+1, z->y_pos_s, z->y_pos_e+1);
return 0;
}
// if(z->x_id == 63 && z->y_id == 9) {
// fprintf(stderr, "-[M::%s]\tq::[%ld,%ld)\tt::[%ld,%ld)\texz->err::%d\tzaln::%u\n", __func__, q_s, q_e + 1, t_s, t_s + exz->pe + 1, exz->err, z->align_length);
// }
// append_window_list(z, q_s, q_e, t_s, t_s + t_end, error, aux_beg, aux_end, thre, w_l, km);
}
// else {
@@ -16022,6 +16031,7 @@ char *tstr, bit_extz_t *exz, uint64_t *v_idx, int64_t block_s, double ovlp_cut,
}
}
// fprintf(stderr, "-[M::%s]\tz->x_id::%u\tz->y_id::%u\tzq::[%u,%u)\ttot_l::%ld\tovl::%ld\n", __func__, z->x_id, z->y_id, z->x_pos_s, z->x_pos_e + 1, tot_l, ovl);
assert(tot_l == ovl); if(r_e) (*r_e) = tot_e;
return (double)(tot_e)/(double)(tot_l);
}
@@ -33561,7 +33571,7 @@ uint8_t inline gen_hc_r_alin_flt_1(overlap_region *oa, uint64_t zi, Candidates_l
uint8_t inline gen_hc_r_alin_flt_1_smp(overlap_region *z, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o,
double err, double e_max, double e_rate, int64_t wsl, int64_t ql, int64_t rid, int64_t khit, uint8_t chem_drop, double align_gap_rate, int64_t align_gap_max, uint8_t *hpf, asg16_v* buf, uint64_t *tot_b)
double err, double e_max, double e_rate, int64_t wsl, int64_t ql, int64_t rid, int64_t khit, uint8_t chem_drop, double align_gap_rate, int64_t align_gap_max, uint8_t *hpf, asg16_v* buf, uint8_t pre_win, uint64_t *tot_b)
{
uint8_t f = 1; double rr; int64_t re;
@@ -33573,7 +33583,9 @@ uint8_t inline gen_hc_r_alin_flt_1_smp(overlap_region *z, Candidates_list *cl, A
// fprintf(stderr, "-0-[M::%s]\ttid::%u\t%.*s\ttlen::%lu\n", __func__, z->y_id, (int)Get_NAME_LENGTH(R_INF, z->y_id), Get_NAME(R_INF, z->y_id), Get_READ_LENGTH(R_INF, z->y_id));
// }
f = align_hc_ed_post_extz(z, rref, qu->seq, tu->seq, exz, err, wsl, OVERLAP_THRESHOLD_HIFI_FILTER/**OVERLAP_THRESHOLD_NOSI_FILTER**/, 0, tot_b);
if(pre_win == 0) {
f = align_hc_ed_post_extz(z, rref, qu->seq, tu->seq, exz, err, wsl, OVERLAP_THRESHOLD_HIFI_FILTER/**OVERLAP_THRESHOLD_NOSI_FILTER**/, 0, tot_b);
}
// if(z->x_id == 142 && z->y_id == 207) {
// fprintf(stderr, "-a-[M::%s::]\tqid::%u\ttid::%u\t%.*s\tq::[%u,%u)\tt::[%u,%u)\tis_match::%u\terr::%u\n", __func__, z->x_id,
// z->y_id, (int)Get_NAME_LENGTH(R_INF, z->y_id), Get_NAME(R_INF, z->y_id),
@@ -34182,7 +34194,7 @@ void gen_hc_aln_small_chn_smp(overlap_region_alloc* ol, Candidates_list *cl, All
// osc[(uint32_t)srt_a[m]], ocn[(uint32_t)srt_a[m]], zm->is_match, zm->non_homopolymer_errors, f);
if(f) continue;
if(!gen_hc_r_alin_flt_1_smp(zm, cl, rref, qu, tu, exz, aux_o, err, e_max, e_rate, wsl, ql, rid, khit, chem_drop, align_gap_rate, align_gap_max, hpf, buf, tot_b)) {
if(!gen_hc_r_alin_flt_1_smp(zm, cl, rref, qu, tu, exz, aux_o, err, e_max, e_rate, wsl, ql, rid, khit, chem_drop, align_gap_rate, align_gap_max, hpf, buf, 0, tot_b)) {
continue;
}
@@ -34412,7 +34424,7 @@ void gen_hc_aln_small_chn_smp_adv(gen_hc_aln_t *ez, uint64_t ql, uint32_t *a_cu,
}
e_max = err * 1.5;
if(!gen_hc_r_alin_flt_1_smp(zm, ez->cl, ez->rref, ez->qu, ez->tu, ez->exz, ez->aux_o, err, e_max, err, w, ql, ez->rid, ez->khit, chem_drop, align_gap_rate, align_gap_max, ez->hpz->a, ez->buf, tot_b)) {
if(!gen_hc_r_alin_flt_1_smp(zm, ez->cl, ez->rref, ez->qu, ez->tu, ez->exz, ez->aux_o, err, e_max, err, w, ql, ez->rid, ez->khit, chem_drop, align_gap_rate, align_gap_max, ez->hpz->a, ez->buf, 0, tot_b)) {
continue;
}
@@ -34464,7 +34476,7 @@ uint64_t rescue_cu_aln(overlap_region_alloc* ol, Candidates_list *cl, All_reads
// z->x_pos_s, z->x_pos_e + 1, z->y_pos_s, z->y_pos_e + 1,
// ((uint32_t)-1) - (srt_a[k]>>32), z->is_match, z->non_homopolymer_errors, zcov, tot_b_cut);
if(!gen_hc_r_alin_flt_1_smp(z, cl, rref, qu, tu, exz, aux_o, err, e_max, e_rate, wsl, ql, rid, khit, chem_drop, align_gap_rate, align_gap_max, hpf, buf, &zcov)) {
if(!gen_hc_r_alin_flt_1_smp(z, cl, rref, qu, tu, exz, aux_o, err, e_max, e_rate, wsl, ql, rid, khit, chem_drop, align_gap_rate, align_gap_max, hpf, buf, 0, &zcov)) {
continue;
}
@@ -34522,7 +34534,7 @@ void rescue_cu_aln_adv(gen_hc_aln_t *ez, int64_t ql, uint64_t *wcut, uint64_t wc
}
e_max = err * 1.5;
if(!gen_hc_r_alin_flt_1_smp(z, ez->cl, ez->rref, ez->qu, ez->tu, ez->exz, ez->aux_o, err, e_max, err, w, ql, ez->rid, ez->khit, chem_drop, align_gap_rate, align_gap_max, ez->hpz->a, ez->buf, &zcov)) {
if(!gen_hc_r_alin_flt_1_smp(z, ez->cl, ez->rref, ez->qu, ez->tu, ez->exz, ez->aux_o, err, e_max, err, w, ql, ez->rid, ez->khit, chem_drop, align_gap_rate, align_gap_max, ez->hpz->a, ez->buf, 0, &zcov)) {
continue;
}
@@ -34659,7 +34671,7 @@ uint64_t gen_hc_r_alin_adp_smp(overlap_region_alloc* ol, Candidates_list *cl, Al
// }
// tot_b0 = tot_b;
if((z->is_match == 0) && (!gen_hc_r_alin_flt_1_smp(&(ol->list[i]), cl, rref, qu, tu, exz, aux_o, err, e_max, e_rate, w.window_length, ql, rid, khit, chem_drop, align_gap_rate, align_gap_max, hpf, buf, &tot_b))) {
if((z->is_match == 0) && (!gen_hc_r_alin_flt_1_smp(&(ol->list[i]), cl, rref, qu, tu, exz, aux_o, err, e_max, e_rate, w.window_length, ql, rid, khit, chem_drop, align_gap_rate, align_gap_max, hpf, buf, 0, &tot_b))) {
// fprintf(stderr, "-m-[M::%s]\ttid::%u\t%.*s(%c)\tq::[%u,%u)\tt::[%u,%u)\ttot_b::%lu\n", __func__, z->y_id, (int)Get_NAME_LENGTH(R_INF, z->y_id), Get_NAME(R_INF, z->y_id), "+-"[z->y_pos_strand],
// z->x_pos_s, z->x_pos_e + 1, z->y_pos_s, z->y_pos_e + 1, tot_b - tot_b0);
continue;
@@ -34734,7 +34746,7 @@ uint64_t gen_hc_r_alin_adp_smp(overlap_region_alloc* ol, Candidates_list *cl, Al
// z->x_pos_s, z->x_pos_e + 1, z->y_pos_s, z->y_pos_e + 1,
// osc[i], ocn[i], z->is_match, z->non_homopolymer_errors);
if(!gen_hc_r_alin_flt_1_smp(&(ol->list[i]), cl, rref, qu, tu, exz, aux_o, err, e_max, e_rate, w.window_length, ql, rid, khit, chem_drop, align_gap_rate, align_gap_max, hpf, buf, &z_cov)) {
if(!gen_hc_r_alin_flt_1_smp(&(ol->list[i]), cl, rref, qu, tu, exz, aux_o, err, e_max, e_rate, w.window_length, ql, rid, khit, chem_drop, align_gap_rate, align_gap_max, hpf, buf, 0, &z_cov)) {
continue;
}
@@ -35045,11 +35057,11 @@ double fcov_rat, uint64_t ch_occ, uint64_t ch_sc)
}
uint64_t gen_hc_r_alin_adp_mmp(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 wl, 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,
uint64_t gen_hc_r_alin_adp_mmp_0(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 wl, 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)
{
uint64_t i, k, bs, ql = qu->length, *wsrt = NULL, wsrt_n = 0, spn0 = 0, tot_b = 0; Window_Pool w; double err, e_max;
overlap_region *z, t; uint32_t *ocn = v32->a, *osc = v32->a + ol->length;
overlap_region *z, t; uint32_t *ocn = v32->a, *osc = v32->a + ol->length;
ol->mapped_overlaps_length = 0;
if(ol->length <= 0 || ql <= 0) return tot_b;
@@ -35069,7 +35081,7 @@ uint64_t gen_hc_r_alin_adp_mmp(overlap_region_alloc* ol, Candidates_list *cl, Al
// fprintf(stderr, "[M::%s]\trid::%u(%c)\t%.*s\tq::[%u,%u)\tt::[%u,%u)\n", __func__,
// z->y_id, "+-"[z->y_pos_strand], (int)Get_NAME_LENGTH(R_INF, z->y_id), Get_NAME(R_INF, z->y_id),
// z->x_pos_s, z->x_pos_e + 1, z->y_pos_s, z->y_pos_e + 1);
if((z->is_match == 0) && (!gen_hc_r_alin_flt_1_smp(z, cl, rref, qu, tu, exz, aux_o, err, e_max, e_rate, w.window_length, ql, rid, khit, chem_drop, align_gap_rate, align_gap_max, hpf, buf, &tot_b))) {
if((z->is_match == 0) && (!gen_hc_r_alin_flt_1_smp(z, cl, rref, qu, tu, exz, aux_o, err, e_max, e_rate, w.window_length, ql, rid, khit, chem_drop, align_gap_rate, align_gap_max, hpf, buf, 0, &tot_b))) {
continue;
}
z->is_match = 1; z->strong = z->without_large_indel = 0;
@@ -35099,6 +35111,364 @@ uint64_t gen_hc_r_alin_adp_mmp(overlap_region_alloc* ol, Candidates_list *cl, Al
}
///[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,
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_GS]; int64_t ez_er[AVX_GS], ez_pe[AVX_GS];
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_GS) {
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);
tstr[4] = tstr[3] + bt_pri_l[3]; recover_UC_Read_sub_region(tstr[4], bt_s[4], bt_pri_l[4], ol[(uint32_t)ffa[fi[4]]].y_pos_strand, rref, ol[(uint32_t)ffa[fi[4]]].y_id);
tstr[5] = tstr[4] + bt_pri_l[4]; recover_UC_Read_sub_region(tstr[5], bt_s[5], bt_pri_l[5], ol[(uint32_t)ffa[fi[5]]].y_pos_strand, rref, ol[(uint32_t)ffa[fi[5]]].y_id);
tstr[6] = tstr[5] + bt_pri_l[5]; recover_UC_Read_sub_region(tstr[6], bt_s[6], bt_pri_l[6], ol[(uint32_t)ffa[fi[6]]].y_pos_strand, rref, ol[(uint32_t)ffa[fi[6]]].y_id);
tstr[7] = tstr[6] + bt_pri_l[6]; recover_UC_Read_sub_region(tstr[7], bt_s[7], bt_pri_l[7], ol[(uint32_t)ffa[fi[7]]].y_pos_strand, rref, ol[(uint32_t)ffa[fi[7]]].y_id);
///same thre due to the same wql
ed_band_cal_semi_64_w_absent_diag_avx8(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;
// ed_band_cal_semi_64_w_absent_diag(tstr[4], bt_pri_l[4], qstr, wql, thre, baux_beg[4], exz); ez_er[4] = exz->err; ez_pe[4] = exz->pe;
// ed_band_cal_semi_64_w_absent_diag(tstr[5], bt_pri_l[5], qstr, wql, thre, baux_beg[5], exz); ez_er[5] = exz->err; ez_pe[5] = exz->pe;
// ed_band_cal_semi_64_w_absent_diag(tstr[6], bt_pri_l[6], qstr, wql, thre, baux_beg[6], exz); ez_er[6] = exz->err; ez_pe[6] = exz->pe;
// ed_band_cal_semi_64_w_absent_diag(tstr[7], bt_pri_l[7], qstr, wql, thre, baux_beg[7], exz); ez_er[7] = exz->err; ez_pe[7] = 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
}
init_base_ed((*exz), thre, bt_pri_l[4], wql); exz->err = ez_er[4]; exz->pe = ez_pe[4]; z = &(ol[(uint32_t)ffa[fi[4]]]);
wts = bt_s[4]; wte = bt_s[4] + exz->pe; t_tot_l = Get_READ_LENGTH((*rref), z->y_id);
aux_beg = baux_beg[4]; aux_end = baux_end[4]; msk = &(ffa[fi[4]]);
// if(z->x_id == 1 && z->y_id == 24) {
// fprintf(stderr, "-4-[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[5], wql); exz->err = ez_er[5]; exz->pe = ez_pe[5]; z = &(ol[(uint32_t)ffa[fi[5]]]);
wts = bt_s[5]; wte = bt_s[5] + exz->pe; t_tot_l = Get_READ_LENGTH((*rref), z->y_id);
aux_beg = baux_beg[5]; aux_end = baux_end[5]; msk = &(ffa[fi[5]]);
// if(z->x_id == 1 && z->y_id == 24) {
// fprintf(stderr, "-5-[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[6], wql); exz->err = ez_er[6]; exz->pe = ez_pe[6]; z = &(ol[(uint32_t)ffa[fi[6]]]);
wts = bt_s[6]; wte = bt_s[6] + exz->pe; t_tot_l = Get_READ_LENGTH((*rref), z->y_id);
aux_beg = baux_beg[6]; aux_end = baux_end[6]; msk = &(ffa[fi[6]]);
// if(z->x_id == 1 && z->y_id == 24) {
// fprintf(stderr, "-6-[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[7], wql); exz->err = ez_er[7]; exz->pe = ez_pe[7]; z = &(ol[(uint32_t)ffa[fi[7]]]);
wts = bt_s[7]; wte = bt_s[7] + exz->pe; t_tot_l = Get_READ_LENGTH((*rref), z->y_id);
aux_beg = baux_beg[7]; aux_end = baux_end[7]; msk = &(ffa[fi[7]]);
// if(z->x_id == 1 && z->y_id == 24) {
// fprintf(stderr, "-7-[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_GS; k++) {
tstr[k] = NULL; bt_pri_l[k] = 0; baux_beg[k] = 0;
}
// for (k = 0; k < mm_k; k++) {
// ed_band_cal_semi_64_w_absent_diag(tstr[k], bt_pri_l[k], qstr, wql, thre, baux_beg[k], exz); ez_er[k] = exz->err; ez_pe[k] = exz->pe;
// }
if(mm_k > 1) {
ed_band_cal_semi_64_w_absent_diag_avx8(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)
{
uint64_t i, k, zk, bs, ql = qu->length, *wsrt = NULL, wsrt_n = 0, spn0 = 0, tot_b = 0, s, e, nwl, in0, q[2], os, oe; Window_Pool w; double err, e_max;
overlap_region *z, t; uint32_t *ocn = v32->a, *osc = v32->a + ol->length; uint8_t rr; uint64_t fi[AVX_GS]; int32_t aux_beg[AVX_GS], aux_end[AVX_GS], t_s[AVX_GS], t_pri_l[AVX_GS];
ol->mapped_overlaps_length = 0;
if(ol->length <= 0 || ql <= 0) return tot_b;
// prt_chain_cluster(ol, cl, a_cu, a_ci, ocn, osc, idx_cu, n_cu, 0, NULL);
///base alignment
err = e_rate; e_max = err * 1.5;
init_Window_Pool(&w, ql, wl0, (int)(1.0/err));
bs = (w.window_length)+(THRESHOLD_MAX_SIZE<<1)+1;
resize_UC_Read(tu, bs<<1); spn0 = sp->n; nwl = w.window_length;
wsrt = mmp_chn_select(ol, cl, sp, ocw, ocn, osc, ql, &wsrt_n, set_match, max_n_chain, max_n_chain_f, chain_cutoff, ave_cov_min, 0.333333, 16, 16);
for (i = 0; i < wsrt_n; i++) {
z = &(ol->list[(uint32_t)wsrt[i]]);
// if(z->x_id == 63 && z->y_id == 9) {
// fprintf(stderr, "-a-[M::%s]\tqid::%u\ttid::%u\tqa::[%u,%u)\tta::[%u,%u)\n", __func__, z->x_id, z->y_id, z->x_pos_s, z->x_pos_e + 1, z->y_pos_s, z->y_pos_e + 1);
// }
if(z->is_match == 0) {
z->w_list.n = 0; z->align_length = 0;
}
wsrt[i] = ((((uint64_t)z->x_pos_s))<<32)|((uint32_t)wsrt[i]);
}
radix_sort_bc64(wsrt, wsrt + wsrt_n);
for (k = 1, i = 0; k <= wsrt_n; k++) {
if (k == wsrt_n || (wsrt[k]>>32) != (wsrt[i]>>32)) {
if(k - i > 1) {
for (zk = i; zk < k; zk++) {
wsrt[zk] = (((uint64_t)(ol->list[(uint32_t)wsrt[zk]].x_pos_e + 1))<<32)|((uint32_t)wsrt[zk]);
}
radix_sort_bc64(wsrt + i, wsrt + k);
}
i = k;
}
}
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]]]);
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(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");
// uint64_t kqs, kqe, kts, kte, kerr;
for (i = 0; i < wsrt_n; i++) {
if((wsrt[i]>>32) == UINT32_MAX) {
// z = &(ol->list[(uint32_t)wsrt[i]]);
// kqs = z->x_pos_s; kqe = z->x_pos_e; kts = z->y_pos_s; kte = z->y_pos_e; kerr = z->non_homopolymer_errors;
// if((z->is_match == 1) || (gen_hc_r_alin_flt_1_smp(z, cl, rref, qu, tu, exz, aux_o, err, e_max, e_rate, w.window_length, ql, rid, khit, chem_drop, align_gap_rate, align_gap_max, hpf, buf, 0, &tot_b))) {
// fprintf(stderr, "-m-[M::%s]\tqid::%u\tq::[%lu,%lu)\ttid::%u\tt::[%lu,%lu)\tzlan::%u\terr::%u\n", __func__, z->x_id, kqs, kqe + 1, z->y_id, kts, kte + 1, z->align_length, z->non_homopolymer_errors);
// exit(1);
// }
continue;
}
z = &(ol->list[(uint32_t)wsrt[i]]);
// kqs = z->x_pos_s; kqe = z->x_pos_e; kts = z->y_pos_s; kte = z->y_pos_e; kerr = z->non_homopolymer_errors;
// if(z->x_id == 1 && z->y_id == 24) {
// fprintf(stderr, "[M::%s]\tzql::%u\talign_length::%u\n", __func__, z->x_pos_e+1-z->x_pos_s, z->align_length);
// }
if(z->is_match == 0) {
if((!pass_qovlp(z->x_pos_e+1-z->x_pos_s, z->align_length, OVERLAP_THRESHOLD_HIFI_FILTER)) || (!gen_hc_r_alin_flt_1_smp(z, cl, rref, qu, tu, exz, aux_o, err, e_max,
e_rate, w.window_length, ql, rid, khit, chem_drop, align_gap_rate, align_gap_max, hpf, buf, 1, &tot_b))) {
// z->x_pos_s = kqs; z->x_pos_e = kqe; z->y_pos_s = kts; z->y_pos_e = kte; z->non_homopolymer_errors = kerr;
// if(gen_hc_r_alin_flt_1_smp(z, cl, rref, qu, tu, exz, aux_o, err, e_max, e_rate, w.window_length, ql, rid, khit, chem_drop, align_gap_rate, align_gap_max, hpf, buf, 0, &tot_b)) {
// fprintf(stderr, "-s-[M::%s]\tqid::%u\ttid::%u\t\n", __func__, z->x_id, z->y_id);
// exit(1);
// }
continue;
}
}
// z->x_pos_s = kqs; z->x_pos_e = kqe; z->y_pos_s = kts; z->y_pos_e = kte; z->non_homopolymer_errors = kerr;
// if(!gen_hc_r_alin_flt_1_smp(z, cl, rref, qu, tu, exz, aux_o, err, e_max, e_rate, w.window_length, ql, rid, khit, chem_drop, align_gap_rate, align_gap_max, hpf, buf, 0, &tot_b)) {
// fprintf(stderr, "-e-[M::%s]\tqid::%u\ttid::%u\tqa::[%u,%u)\tta::[%u,%u)\n", __func__, z->x_id, z->y_id, z->x_pos_s, z->x_pos_e + 1, z->y_pos_s, z->y_pos_e + 1);
// continue;
// }
// if(z->x_pos_s != kqs || z->x_pos_e != kqe || z->y_pos_s != kts || z->y_pos_e != kte || z->non_homopolymer_errors != kerr) {
// fprintf(stderr, "-dbg-[M::%s]\tqid::%u\ttid::%u\tzqs::%u(%lu)\tzqe::%u(%lu)\tzts::%u(%lu)\tzte::%u(%lu)\terr::%u(%lu)\n", __func__, z->x_id, z->y_id,
// z->x_pos_s, kqs, z->x_pos_e + 1, kqe + 1, z->y_pos_s, kts, z->y_pos_e + 1, kte + 1, z->non_homopolymer_errors, kerr);
// exit(1);
// }
z->is_match = 1; z->strong = z->without_large_indel = 0;
}
// for (i = 0; i < wsrt_n; i++) {
// z = &(ol->list[(uint32_t)wsrt[i]]);
// if((z->is_match == 0) && (!gen_hc_r_alin_flt_1_smp(z, cl, rref, qu, tu, exz, aux_o, err, e_max, e_rate, w.window_length, ql, rid, khit, chem_drop, align_gap_rate, align_gap_max, hpf, buf, &tot_b))) {
// continue;
// }
// z->is_match = 1; z->strong = z->without_large_indel = 0;
// }
for (i = k = 0; i < ol->length; i++) {///primary chain
z = &(ol->list[i]);
if(z->is_match == 0) continue;
if(k != i) {
t = ol->list[k];
ol->list[k] = ol->list[i];
ol->list[i] = t;
}
k++;
}
// print_mm_wins_all(wcut, wcut_n, ocw, ql, 0);
// print_mm_wins_all(wcut, wcut_n, ocw, ql, 1);
ol->length = k;
// prt_chain_cluster(ol, cl, a_cu, a_ci, ocn, osc, idx_cu, n_cu, 1, NULL);
// fprintf(stderr, "-[M::%s]\trid::%ld\ttot_b::%lu\tql::%lu\tmax_n_chain::%ld\ttot_b_cov:::%ld\n", __func__, rid, tot_b, ql, max_n_chain, tot_b/ql);
// exit(1);
sp->n = spn0;
// if(ol->length <= 0) return tot_b;
return tot_b;
}
///need to consider coverage, this information is missing right now (currently only use numbers)
uint64_t gen_hc_r_alin_adp_smp_ff_ec(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, int64_t rid, asg64_v *sp, uint64_t ocw, uint32_t *ocn, uint32_t *osc, uint64_t max_n_chain, uint64_t max_n_chain_f, uint64_t chain_cutoff, uint64_t ave_cov_min)
{
@@ -35479,7 +35849,7 @@ void gen_hc_r_alin_adv_adp_smp(gen_hc_aln_t *ez, uint32_t *a_cu, uint32_t *a_ci,
// osc[i], ocn[i], z->is_match, z->non_homopolymer_errors);
if((z->is_match == 0) && (!gen_hc_r_alin_flt_1_smp(z, ez->cl, ez->rref, ez->qu, ez->tu, ez->exz, ez->aux_o, err, e_max, err, w, ql, ez->rid, ez->khit, chem_drop, align_gap_rate, align_gap_max, ez->hpz->a, ez->buf, &tot_b))) {
if((z->is_match == 0) && (!gen_hc_r_alin_flt_1_smp(z, ez->cl, ez->rref, ez->qu, ez->tu, ez->exz, ez->aux_o, err, e_max, err, w, ql, ez->rid, ez->khit, chem_drop, align_gap_rate, align_gap_max, ez->hpz->a, ez->buf, 0, &tot_b))) {
continue;
}
@@ -35565,7 +35935,7 @@ void gen_hc_r_alin_adv_adp_smp(gen_hc_aln_t *ez, uint32_t *a_cu, uint32_t *a_ci,
// z->x_pos_s, z->x_pos_e + 1, z->y_pos_s, z->y_pos_e + 1,
// osc[i], ocn[i], z->is_match, z->non_homopolymer_errors);
if(!gen_hc_r_alin_flt_1_smp(z, ez->cl, ez->rref, ez->qu, ez->tu, ez->exz, ez->aux_o, err, e_max, err, w, ql, ez->rid, ez->khit, chem_drop, align_gap_rate, align_gap_max, ez->hpz->a, ez->buf, &tot_b)) {
if(!gen_hc_r_alin_flt_1_smp(z, ez->cl, ez->rref, ez->qu, ez->tu, ez->exz, ez->aux_o, err, e_max, err, w, ql, ez->rid, ez->khit, chem_drop, align_gap_rate, align_gap_max, ez->hpz->a, ez->buf, 0, &tot_b)) {
continue;
}
Correct.h
+3 -1
View File
@@ -1437,7 +1437,9 @@ bit_extz_t *exz, double e_rate, int64_t qs);
uint64_t gen_hc_r_alin_adp_smp_ff_ec(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, int64_t rid, asg64_v *sp, uint64_t ocw, uint32_t *ocn, uint32_t *osc, uint64_t max_n_chain, uint64_t max_n_chain_f, uint64_t chain_cutoff, uint64_t ave_cov_min);
uint64_t gen_hc_r_alin_adp_smp(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 wl, 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, uint32_t *a_cu, uint32_t *a_ci, uint32_t *ocn, uint32_t *osc, uint64_t *idx_cu, uint64_t n_cu, 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, uint8_t is_dedup);
uint64_t gen_hc_r_alin_adp_mmp(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 wl, 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,
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 wl, 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);
uint64_t gen_hc_r_alin_adp_mmp_0(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 wl, 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);
void gen_hc_r_alin_adv_adp_smp(gen_hc_aln_t *ez, uint32_t *a_cu, uint32_t *a_ci, uint32_t *ocn, uint32_t *osc, uint64_t *idx_cu, uint64_t n_cu, uint8_t set_match);
void pp_chn_a(overlap_region *z, Candidates_list *cl, uint8_t is_raw);
Levenshtein_distance.h
+272
View File
@@ -12,6 +12,8 @@
#include <stdio.h>
#include "kvec.h"
#define AVX_GS 8
extern const unsigned char seq_nt4_table[256];
typedef uint64_t Word;
typedef uint32_t Word_32;
@@ -3734,6 +3736,276 @@ inline void ed_band_cal_extension_64_1_w_trace(char *pstr, int32_t pn, char *tst
return;
}
#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;
}
}
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)
{
init_base_ed(*ez, thre, pn, tn); ez->ps = ez->pe = -1; ez->ts = 0; ez->te = tn-1;
Makefile
+1 -1
View File
@@ -1,6 +1,6 @@
CXX= g++
CC= gcc
CXXFLAGS= -g -O3 -msse4.2 -mpopcnt -fomit-frame-pointer -Wall
CXXFLAGS= -g -O3 -mavx512f -msse4.2 -mpopcnt -fomit-frame-pointer -Wall
CFLAGS= $(CXXFLAGS)
CPPFLAGS=
INCLUDES=
ecovlp.cpp
+4 -4
View File
@@ -3369,7 +3369,7 @@ 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(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, hpz->a,
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, 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));
@@ -3401,7 +3401,7 @@ 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(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, hpz->a,
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, 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);
}
}
@@ -4211,7 +4211,7 @@ static void worker_hap_ec(void *data, long i, int tid)
// if(i != 339646) return;
// if(i!=854835) return;
// if(i != 533) return;
// if(i != 1) return;
// if(i % 100000 == 0) fprintf(stderr, "-a-[M::%s-beg] rid->%ld\n", __func__, i);
@@ -8045,7 +8045,7 @@ static void update_scb0(void *data, long i, int tid)
if(!scc.f[i]) return;
ma_hit_t_alloc *ov, *os; uint64_t k, kr, qn, tn, ql, tl, qs, qe, ts, te; ma_hit_t *z, *r;
ma_hit_t_alloc *ov, *os; uint64_t k, kr, qn, tn, ql, tl, qs, qe, ts, te; ma_hit_t *z, *r = NULL;
uint64_t ck; uint16_t op, bq, bt; uint32_t cl;
ov = &(R_INF.paf[i]);