diff --git a/CommandLines.cpp b/CommandLines.cpp index aa63f9b..9c6833e 100644 --- a/CommandLines.cpp +++ b/CommandLines.cpp @@ -88,6 +88,8 @@ static ko_longopt_t long_options[] = { { "resume", ko_no_argument, 370}, { "flt-kocc", ko_required_argument, 371}, { "chn-occ", ko_required_argument, 372}, + { "dbg-in1", ko_required_argument, 373}, + { "dbg-in2", ko_required_argument, 374}, // { "path-round", ko_required_argument, 348}, { 0, 0, 0 } }; @@ -415,6 +417,9 @@ void init_opt(hifiasm_opt_t* asm_opt) asm_opt->hmo_cov_ss = -1; asm_opt->het_cov_ss = -1; asm_opt->chn_occ = 2; + + asm_opt->dbg_run_1 = NULL; + asm_opt->dbg_run_2 = NULL; } void destory_enzyme(enzyme* f) @@ -1075,8 +1080,11 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt) fprintf(stderr, "[ERROR] chain cutoff should be > 0 (--chn-occ)"); return 1; } - } - else if (c == 'l') { ///0: disable purge_dup; 1: purge containment; 2: purge overlap + } else if (c == 373) { + asm_opt->dbg_run_1 = opt.arg; + } else if (c == 374) { + asm_opt->dbg_run_2 = 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); } else if (c == 's') asm_opt->purge_simi_rate_l2 = asm_opt->purge_simi_rate_l3 = atof(opt.arg); diff --git a/CommandLines.h b/CommandLines.h index 196acde..84fe77f 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -5,7 +5,7 @@ #include #include -#define HA_VERSION "0.25.0-r856" +#define HA_VERSION "0.25.0-r857" #define VERBOSE 0 @@ -192,6 +192,9 @@ typedef struct { int64_t hmo_cov_ss; int64_t het_cov_ss; int64_t chn_occ; + + char *dbg_run_1, *dbg_run_2; + } hifiasm_opt_t; extern hifiasm_opt_t asm_opt; diff --git a/Correct.cpp b/Correct.cpp index 66f4768..59c210e 100644 --- a/Correct.cpp +++ b/Correct.cpp @@ -15378,7 +15378,7 @@ uint32_t align_ul_ed_post_extz(overlap_region *z, const ul_idx_t *uref, hpc_t *h uint32_t align_hc_ed_post_extz(overlap_region *z, All_reads *rref, char* qstr, char *tstr, bit_extz_t *exz, double e_rate, int64_t w_l, double ovlp_cut, int64_t force_aln, uint64_t *tot_b) { - int64_t q_s, q_e, nw, k, q_l, t_l, t_tot_l, aux_beg, aux_end, t_s, thre, aln_l, t_pri_l; + int64_t q_s, q_e, nw, k, q_l, t_l, t_tot_l, aux_beg, aux_end, t_s, thre, aln_l, t_pri_l; ///uint64_t tot_b0 = 0; char *q_string, *t_string; z->w_list.n = 0; z->is_match = 0; z->align_length = 0; nw = get_num_wins(z->x_pos_s, z->x_pos_e+1, w_l); @@ -15416,6 +15416,13 @@ uint32_t align_hc_ed_post_extz(overlap_region *z, All_reads *rref, char* qstr, c if(tot_b) (*tot_b) += q_l; + // if(z->y_id == 5544) { + // fprintf(stderr, "-z-[M::%s]\ttid::%u\t%.*s(%c)\tq::[%ld,%ld)\terr::%d\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], + // q_s, q_e + 1, exz->err, (*tot_b) - tot_b0); + // tot_b0 = (*tot_b); + // } + + // 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); @@ -31255,12 +31262,12 @@ uint64_t qick_realign(Candidates_list *cl, overlap_region *a, uint64_t n_a, uint return k; } -void gen_hc_r_alin(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 *kp, uint8_t *hpf) +uint64_t gen_hc_r_alin(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 *kp, uint8_t *hpf) { - uint64_t i, bs, k, ql = qu->length, ks = 0, ke = sec_aln_win, kcnt = sec_aln_cov; Window_Pool w; double err, e_max, rr; int64_t re, kk; + uint64_t i, bs, k, ql = qu->length, ks = 0, ke = sec_aln_win, kcnt = sec_aln_cov, tot_b = 0/**, tot_b0 = 0**/; Window_Pool w; double err, e_max, rr; int64_t re, kk; overlap_region t; overlap_region *z; //asg64_v iidx, buf, buf1; ol->mapped_overlaps_length = 0; - if(ol->length <= 0) return; + if(ol->length <= 0) return tot_b; // if(ol->length && ol->list[0].x_id == 19350) e_rate = 0.1; ///base alignment @@ -31271,12 +31278,22 @@ void gen_hc_r_alin(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rre // fprintf(stderr, "+[M::%s] window_length::%lld, err::%f\n", __func__, w.window_length, err); // if(kp) kp->n = 0; + // fprintf(stderr, "-0-[M::%s] tot_b::%lu\n", __func__, tot_b); + for (i = k = 0; i < ol->length; i++) { z = &(ol->list[i]); z->shared_seed = z->non_homopolymer_errors; z->non_homopolymer_errors = (uint32_t)-1;///for index // if(z->y_id == 234) fprintf(stderr, "-z-[M::%s] tid::%u\t%.*s\n", __func__, z->y_id, (int)Get_NAME_LENGTH(R_INF, z->y_id), Get_NAME(R_INF, z->y_id)); - - if(!align_hc_ed_post_extz(z, rref, qu->seq, tu->seq, exz, err, w.window_length, OVERLAP_THRESHOLD_HIFI_FILTER, 0, NULL)) continue; + //tot_b0 = tot_b; + if(!align_hc_ed_post_extz(z, rref, qu->seq, tu->seq, exz, err, w.window_length, OVERLAP_THRESHOLD_HIFI_FILTER, 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; + } + // else { + // fprintf(stderr, "-im-[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); + // } // if(z->y_id == 234) fprintf(stderr, "-m-[M::%s] tid::%u\t%.*s\n", __func__, z->y_id, (int)Get_NAME_LENGTH(R_INF, z->y_id), Get_NAME(R_INF, z->y_id)); @@ -31328,6 +31345,8 @@ void gen_hc_r_alin(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rre } } + // fprintf(stderr, "-1-[M::%s] tot_b::%lu\n", __func__, tot_b); + // fprintf(stderr, "-[M::%s] window_length::%lld, err::%f\n", __func__, w.window_length, err); if((kp) && (kcnt < sec_aln_cov)) { @@ -31338,7 +31357,8 @@ void gen_hc_r_alin(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rre // fprintf(stderr, "--[M::%s] window_length::%lld, err::%f\n", __func__, w.window_length, err); ol->length = k; - if(ol->length <= 0) return; + // if(ol->length <= 0) return; + return tot_b; } @@ -32445,7 +32465,7 @@ 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\n", __func__, z->y_id, (int)Get_NAME_LENGTH(R_INF, z->y_id), Get_NAME(R_INF, z->y_id)); // } - f = align_hc_ed_post_extz(z, rref, qu->seq, tu->seq, exz, err, wsl, OVERLAP_THRESHOLD_NOSI_FILTER, 0, tot_b); + 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), @@ -33292,15 +33312,15 @@ void rescue_cu_aln_adv(gen_hc_aln_t *ez, int64_t ql, uint64_t *wcut, uint64_t wc ///need to consider coverage, this information is missing right now (currently only use numbers) -void 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) +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 i, k, bs, ql = qu->length, *wcut = NULL, wcut_n = 0, *wsrt = NULL, wsrt_n = 0, lch = 0, wsrt_n0, spn0 = 0, tz, tot_b = 0, t_cov0, z_cov, focv = 0, fcov_0 = 0, rescue_b_cut = 0; Window_Pool w; double err, e_max; + uint64_t i, k, bs, ql = qu->length, *wcut = NULL, wcut_n = 0, *wsrt = NULL, wsrt_n = 0, lch = 0, wsrt_n0, spn0 = 0, tz, tot_b = 0, t_cov0, z_cov, focv = 0, fcov_0 = 0, rescue_b_cut = 0/**, tot_b0 = 0**/; Window_Pool w; double err, e_max; overlap_region *z, t; uint8_t fc = (((ol->length > max_n_chain) && (n_cu > 0))?(1):(0)), wf, ff; int32_t s[4]; s[0] = s[1] = s[2] = s[3] = 0; uint32_t scn[4]; scn[0] = scn[1] = scn[2] = scn[3] = 0; ol->mapped_overlaps_length = 0; - if(ol->length <= 0 || ql <= 0) return; + 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); @@ -33329,7 +33349,7 @@ void gen_hc_r_alin_adp_smp(overlap_region_alloc* ol, Candidates_list *cl, All_re rescue_b_cut = (fcov_0*0.75)*ql; } - + // fprintf(stderr, "-0-[M::%s] tot_b::%lu\n", __func__, tot_b); // fprintf(stderr, "-0-[M::%s] focv::%lu, ave_cov_min::%lu, fcov_0::%lu, max_n_chain::%lu, wcut_n::%lu, ql::%lu, chain_cutoff::%lu\n", __func__, focv, ave_cov_min, fcov_0, max_n_chain, wcut_n, ql, chain_cutoff); ////this is very important!!! main issues coming from uneven read length, some short reads may have too many chains need to be calculated? @@ -33412,10 +33432,16 @@ void gen_hc_r_alin_adp_smp(overlap_region_alloc* ol, Candidates_list *cl, All_re // osc[i], ocn[i], z->is_match, z->non_homopolymer_errors); // } - + // 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))) { + // 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; - } + } + // else { + // fprintf(stderr, "-im-[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); + // } // if(z->x_id == 142 && z->y_id == 207) { // fprintf(stderr, "-5-[M::%s::]\tqid::%u\ttid::%u\t%.*s\tq::[%u,%u)\tt::[%u,%u)\tsc::%u\tcn::%u\tis_match::%u\terr::%u\n", __func__, z->x_id, @@ -33436,6 +33462,8 @@ void gen_hc_r_alin_adp_smp(overlap_region_alloc* ol, Candidates_list *cl, All_re // (uint32_t)wsrt[i], wsrt[i]>>32); // } + // fprintf(stderr, "-1-[M::%s] tot_b::%lu\n", __func__, tot_b); + if(fc && wsrt_n) { t_cov0 = cal_mm_wins_cov(wcut, wcut_n, ocw, ql, fcov_0); @@ -33511,16 +33539,27 @@ void gen_hc_r_alin_adp_smp(overlap_region_alloc* ol, Candidates_list *cl, All_re tot_b += z_cov; } + // fprintf(stderr, "-2-[M::%s] tot_b::%lu\n", __func__, tot_b); + if(lch) { gen_hc_aln_small_chn_smp(ol, 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, wcut, wcut_n, ocw, wsrt, wsrt_n, ocn, osc, chain_cutoff, /**3, 3**/16, 16); } + // fprintf(stderr, "-3-[M::%s] tot_b::%lu\n", __func__, tot_b); + if(rescue_b_cut > 0 && n_cu > 0) { tot_b += rescue_cu_aln(ol, 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, wcut, wcut_n, ocw, wsrt, a_cu, a_ci, idx_cu, n_cu, rescue_b_cut); } + // fprintf(stderr, "-4-[M::%s] tot_b::%lu\n", __func__, tot_b); + + if(!is_dedup) { + sp->n = spn0; + return tot_b; + } + for (i = k = 0; i < ol->length; i++) {///primary chain z = &(ol->list[i]); @@ -33550,7 +33589,8 @@ void gen_hc_r_alin_adp_smp(overlap_region_alloc* ol, Candidates_list *cl, All_re // 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; + // if(ol->length <= 0) return tot_b; + return tot_b; } @@ -33908,12 +33948,12 @@ void gen_hc_r_alin_adv_adp_smp(gen_hc_aln_t *ez, uint32_t *a_cu, uint32_t *a_ci, -void gen_hc_r_alin_nec(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 *kp, uint8_t *hpf) +uint64_t gen_hc_r_alin_nec(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 *kp, uint8_t *hpf) { - uint64_t i, bs, k, ql = qu->length, ks = 0, ke = sec_aln_win, kcnt = sec_aln_cov; Window_Pool w; double err, e_max, rr; int64_t re, kk; + uint64_t i, bs, k, ql = qu->length, ks = 0, ke = sec_aln_win, kcnt = sec_aln_cov, tot_b = 0;; Window_Pool w; double err, e_max, rr; int64_t re, kk; overlap_region t; overlap_region *z; //asg64_v iidx, buf, buf1; ol->mapped_overlaps_length = 0; - if(ol->length <= 0) return; + if(ol->length <= 0) return tot_b; // if(ol->length && ol->list[0].x_id == 19350) e_rate = 0.1; ///base alignment @@ -33933,7 +33973,7 @@ void gen_hc_r_alin_nec(overlap_region_alloc* ol, Candidates_list *cl, All_reads // if(z->x_id == 19350 && z->y_id == 19324) fprintf(stderr, "-z-[M::%s] tid::%u\t%.*s\n", __func__, z->y_id, (int)Get_NAME_LENGTH(R_INF, z->y_id), Get_NAME(R_INF, z->y_id)); - if(!align_hc_ed_post_extz(z, rref, qu->seq, tu->seq, exz, err, w.window_length, OVERLAP_THRESHOLD_HIFI_FILTER, 0, NULL)) continue; + if(!align_hc_ed_post_extz(z, rref, qu->seq, tu->seq, exz, err, w.window_length, OVERLAP_THRESHOLD_HIFI_FILTER, 0, &tot_b)) continue; // if(z->x_id == 19350 && z->y_id == 19324) fprintf(stderr, "-m-[M::%s] tid::%u\t%.*s\n", __func__, z->y_id, (int)Get_NAME_LENGTH(R_INF, z->y_id), Get_NAME(R_INF, z->y_id)); @@ -33993,7 +34033,8 @@ void gen_hc_r_alin_nec(overlap_region_alloc* ol, Candidates_list *cl, All_reads ol->length = k; - if(ol->length <= 0) return; + // if(ol->length <= 0) return; + return tot_b; } void gen_hc_r_alin_nec_adv(gen_hc_aln_t *ez) diff --git a/Correct.h b/Correct.h index 933c59c..fa1fb5f 100644 --- a/Correct.h +++ b/Correct.h @@ -1433,17 +1433,17 @@ const ul_idx_t *uref, char* qstr, UC_Read *tu, overlap_region_alloc *ol, overlap bit_extz_t *exz, double e_rate, int64_t qs); -void 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); +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); 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); -void gen_hc_r_alin(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 *kp, uint8_t *hpf); +uint64_t gen_hc_r_alin(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 *kp, uint8_t *hpf); void gen_hc_r_alin_flt(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, overlap_region *aux_b, 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 *kp, 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); void gen_hc_r_alin_adp(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, overlap_region *aux_b, 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 *kp, 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); void gen_hc_r_alin_adv(gen_hc_aln_t *ez); -void gen_hc_r_alin_nec(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 *kp, uint8_t *hpf); +uint64_t gen_hc_r_alin_nec(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 *kp, uint8_t *hpf); void gen_hc_r_alin_nec_adv(gen_hc_aln_t *ez); uint64_t gen_hc_r_alin_re(overlap_region* z, Candidates_list *cl, char* qstr, uint64_t ql, char* tstr, uint64_t tl, 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); void rphase_hc(overlap_region_alloc* ol, All_reads *rref, haplotype_evdience_alloc* hp, UC_Read* qu, UC_Read* tu, kv_ul_ov_t *c_idx, asg64_v* idx, asg64_v* buf, int64_t bd, int64_t wl, int64_t ql, uint8_t occ_thres/**, uint8_t is_dbg**/, uint64_t rid, uint64_t hpc_len, uint64_t std_bs, Chain_Data *dp, asg8_v *q8, asg8_v *t8, uint8_t lindel, uint64_t tcut, uint64_t site_sc, int64_t h0_w, asg32_v *b32); diff --git a/Overlaps.h b/Overlaps.h index 7417a78..25816b4 100644 --- a/Overlaps.h +++ b/Overlaps.h @@ -263,6 +263,7 @@ typedef struct { size_t n, m; uint64_t *a; } asg64_v; typedef struct { size_t n, m; uint32_t *a; } asg32_v; typedef struct { size_t n, m; ma_utg_t *a;} ma_utg_v; typedef struct { asg64_v idx; kv_ul_ov_t srt;} mask_ul_ov_t; +typedef struct { uint64_t n, m; char *a; } asgchr_v; typedef struct { ma_utg_v u; diff --git a/ecovlp.cpp b/ecovlp.cpp index e16b031..204cf31 100644 --- a/ecovlp.cpp +++ b/ecovlp.cpp @@ -27,6 +27,17 @@ typedef struct { } dbg_cnt_ss; dbg_cnt_ss *dbg_a = NULL; +typedef struct { + dbg_cnt_ss *ssa, *ssb; + uint64_t nrid; + asg64_v srt_aa; ///uint64_t srt_aa_cut; + asg64_v srt_ba; ///uint64_t srt_ba_cut; + asgchr_v *spt_mul; uint64_t mul_n; + uint64_t *fa, fn; asg32_v fthr; +} dbg_cmp_ss; +dbg_cmp_ss *dbgss = NULL; + + KDQ_INIT(uint32_t) typedef struct { @@ -3019,12 +3030,12 @@ void ggen_chain_clus_0(overlap_region_alloc* ol, asg32_v *v32, asg64_v *v64, uin } } -void gen_hc_r_alin_ea_flt(ha_abuf_t *ab, overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, uint64_t max_n_chain, uint64_t max_n_chain_f, uint64_t chain_cutoff, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v *buf, asg64_v *srt, ma_hit_t_alloc *in, 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 *bp, uint64_t ocw, asg8_v *hpz, asg32_v *v32) +uint64_t gen_hc_r_alin_ea_flt(ha_abuf_t *ab, overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, uint64_t max_n_chain, uint64_t max_n_chain_f, uint64_t chain_cutoff, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v *buf, asg64_v *srt, ma_hit_t_alloc *in, 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 *bp, uint64_t ocw, asg8_v *hpz, asg32_v *v32, uint8_t is_dedup) { - if(ol->length <= 0) return; + if(ol->length <= 0) return 0; - uint32_t *a_cu = NULL, *a_ci = NULL, *ocn = NULL, *osc = NULL; uint64_t k, *idx_cu = NULL, n_cu = 0; v32->n = 0; + uint32_t *a_cu = NULL, *a_ci = NULL, *ocn = NULL, *osc = NULL; uint64_t k, *idx_cu = NULL, n_cu = 0, tot_b = 0; v32->n = 0; if(ol->length > max_n_chain) { gen_chain_clus(ab, ol, cl, v32);///ol->align_length has not been set to 0 } @@ -3042,8 +3053,8 @@ void gen_hc_r_alin_ea_flt(ha_abuf_t *ab, overlap_region_alloc* ol, Candidates_li if(!(srt->n)) { gen_ff_hpc(hpz, qu->seq, qu->length, HPC_RR_Q, HPC_CC_Q); - gen_hc_r_alin_adp_smp(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, a_cu, a_ci, ocn, osc, - idx_cu, n_cu, 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); + tot_b = gen_hc_r_alin_adp_smp(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, a_cu, a_ci, ocn, osc, + idx_cu, n_cu, 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, is_dedup); } else { kv_resize(uint64_t, *srt, (srt->n + ol->length)); ei = srt->a; en = srt->n; oi = srt->a + srt->n; on = ol->length; @@ -3074,22 +3085,24 @@ void gen_hc_r_alin_ea_flt(ha_abuf_t *ab, overlap_region_alloc* ol, Candidates_li if(on > nec) { gen_ff_hpc(hpz, qu->seq, qu->length, HPC_RR_Q, HPC_CC_Q); - gen_hc_r_alin_adp_smp(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, a_cu, a_ci, ocn, osc, - idx_cu, n_cu, 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); + tot_b = gen_hc_r_alin_adp_smp(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, a_cu, a_ci, ocn, osc, + idx_cu, n_cu, 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, is_dedup); } } if(ol->length) srt_olst(ol); + + return tot_b; } -void gen_hc_r_alin_ea(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, asg64_v *srt, ma_hit_t_alloc *in, uint8_t chem_drop, double align_gap_rate, int64_t align_gap_max, +uint64_t gen_hc_r_alin_ea(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, asg64_v *srt, ma_hit_t_alloc *in, 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 *kp, asg8_v *hpz) { - if(ol->length <= 0) return; + if(ol->length <= 0) return 0; // uint64_t k, l, i, s, m, mm_k, *ei, en, *oi, on, tid, trev, nec; int64_t sc, mm_sc, plus, minus; overlap_region *z, t; ma_hit_t *p; - uint64_t k, i, m, *ei, en, *oi, on, tid, trev, nec; overlap_region *z; ma_hit_t *p; + uint64_t k, i, m, *ei, en, *oi, on, tid, trev, nec, tot_b = 0; overlap_region *z; ma_hit_t *p; srt->n = 0; for (k = 0; k < in->length; k++) { if(in->buffer[k].el) { @@ -3100,7 +3113,7 @@ void gen_hc_r_alin_ea(overlap_region_alloc* ol, Candidates_list *cl, All_reads * if(!(srt->n)) { gen_ff_hpc(hpz, qu->seq, qu->length, HPC_RR_Q, HPC_CC_Q); - gen_hc_r_alin(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, kp, hpz->a); + tot_b = gen_hc_r_alin(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, kp, hpz->a); } 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); @@ -3136,7 +3149,7 @@ void gen_hc_r_alin_ea(overlap_region_alloc* ol, Candidates_list *cl, All_reads * if(on > nec) { gen_ff_hpc(hpz, qu->seq, qu->length, HPC_RR_Q, HPC_CC_Q); - gen_hc_r_alin_nec(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, kp, hpz->a); + tot_b = gen_hc_r_alin_nec(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, kp, hpz->a); } // fprintf(stderr, "[M::%s] srt->n::%u, nec::%lu, on::%lu\n", __func__, (uint32_t)srt->n, nec, on); @@ -3144,6 +3157,8 @@ void gen_hc_r_alin_ea(overlap_region_alloc* ol, Candidates_list *cl, All_reads * // snprintf(NULL, 0, "dwn::%u\tdcn::%u", (uint32_t)aux_o->w_list.n, (uint32_t)aux_o->w_list.c.n); } + return tot_b; + /** if(ol->length > 1) {///for duplicated chains overlap_region_sort_y_id(ol->list, ol->length); @@ -3435,6 +3450,72 @@ void stderr_phase_ovlp(overlap_region_alloc* ol) } } + +void stderr_phase_ovlp_buf(uint32_t srid, overlap_region_alloc* ol, Candidates_list *cl, uint32_t *ocn, uint32_t *osc, asgchr_v *ou, dbg_cnt_ss *mz, const char *cmd) +{ + int64_t on = ol->length, k, t, prt_len; overlap_region *z; + if(!on) return; + uint64_t qry_n = 0, rid, ref_n, qid, sc, chn; int32_t err, mm; + rid = ol->list[0].x_id; ref_n = Get_READ_LENGTH(R_INF, rid); + + + for (k = qry_n = 0; k < on; k++) { + z = &(ol->list[k]); qid = ol->list[k].y_id; + + if(ocn && osc) { + sc = osc[k]; chn = ocn[k]; err = ((z->non_homopolymer_errors >= (UINT32_MAX-1))?(-1):(z->non_homopolymer_errors)); + mm = z->is_match; if((mm == 0) && (z->non_homopolymer_errors == UINT32_MAX)) mm = -1;///skipped + } else { + sc = z->shared_seed; + for (t = z->non_homopolymer_errors; (t < cl->length) && (cl->list[t].readID == cl->list[z->non_homopolymer_errors].readID); t++); + chn = t - z->non_homopolymer_errors; assert(chn); + err = -1; + mm = -2;///not base alignment performed + } + if(mm == -1) continue; + qry_n += z->x_pos_e + 1 - z->x_pos_s; + } + + ///label + // fprintf(fn, "rid::%lu\ttot_bs::%lu\taln_bs::%lu\tchn_tm::%f\taln_tm::%f\tphs_tm::%f\tcns_tm::%f\n", k, p[k].fbs, p[k].faln, p[k].chn_tm, p[k].aln_tm, p[k].phs_tm, p[k].cns_tm); + prt_len = snprintf(NULL, 0, "(%s)srid::%u\trid::%lu\tcov_bs::%lu\ttot_bs::%lu\taln_bs::%lu\tchn_tm::%f\taln_tm::%f\tphs_tm::%f\tcns_tm::%f\n", cmd, srid, rid, + qry_n, mz->fbs, mz->faln, mz->chn_tm, mz->aln_tm, mz->phs_tm, mz->cns_tm);///snprintf exclude \0 + ou->n += prt_len; if((ou->n + 1) > ou->m) kv_resize(char, *ou, (ou->n + 1)); + snprintf(ou->a + ou->n - prt_len, prt_len + 1, "(%s)srid::%u\trid::%lu\tcov_bs::%lu\ttot_bs::%lu\taln_bs::%lu\tchn_tm::%f\taln_tm::%f\tphs_tm::%f\tcns_tm::%f\n", cmd, srid, rid, + qry_n, mz->fbs, mz->faln, mz->chn_tm, mz->aln_tm, mz->phs_tm, mz->cns_tm); + + + + + + for (k = 0; k < on; k++) { + z = &(ol->list[k]); qid = ol->list[k].y_id; + qry_n = Get_READ_LENGTH(R_INF, qid); + + if(ocn && osc) { + sc = osc[k]; chn = ocn[k]; err = ((z->non_homopolymer_errors >= (UINT32_MAX-1))?(-1):(z->non_homopolymer_errors)); + mm = z->is_match; if((mm == 0) && (z->non_homopolymer_errors == UINT32_MAX)) mm = -1;///skipped + } else { + sc = z->shared_seed; + for (t = z->non_homopolymer_errors; (t < cl->length) && (cl->list[t].readID == cl->list[z->non_homopolymer_errors].readID); t++); + chn = t - z->non_homopolymer_errors; assert(chn); + err = -1; + mm = -2;///not base alignment performed + } + if(mm == -1) continue; + + ///fprintf(fn, "rid::%lu\ttot_bs::%lu\taln_bs::%lu\tchn_tm::%f\taln_tm::%f\tphs_tm::%f\tcns_tm::%f\n", k, p[k].fbs, p[k].faln, p[k].chn_tm, p[k].aln_tm, p[k].phs_tm, p[k].cns_tm); + prt_len = snprintf(NULL, 0, "%.*s(qid::%lu)\tql::%lu\tq::[%u,%u)\t%c\t%.*s(tid::%lu)\ttl::%lu\tt::[%u,%u)\tmatch::%d\terr::%d\tchn::%lu\tosc::%lu\n", + (int32_t)Get_NAME_LENGTH(R_INF, rid), Get_NAME(R_INF, rid), rid, ref_n, z->x_pos_s, z->x_pos_e + 1, "+-"[z->y_pos_strand], + (int32_t)Get_NAME_LENGTH(R_INF, qid), Get_NAME(R_INF, qid), qid, qry_n, z->y_pos_s, z->y_pos_e + 1, mm, err, chn, sc); + ou->n += prt_len; if((ou->n + 1) > ou->m) kv_resize(char, *ou, (ou->n + 1)); + snprintf(ou->a + ou->n - prt_len, prt_len + 1, "%.*s(qid::%lu)\tql::%lu\tq::[%u,%u)\t%c\t%.*s(tid::%lu)\ttl::%lu\tt::[%u,%u)\tmatch::%d\terr::%d\tchn::%lu\tosc::%lu\n", + (int32_t)Get_NAME_LENGTH(R_INF, rid), Get_NAME(R_INF, rid), rid, ref_n, z->x_pos_s, z->x_pos_e + 1, "+-"[z->y_pos_strand], + (int32_t)Get_NAME_LENGTH(R_INF, qid), Get_NAME(R_INF, qid), qid, qry_n, z->y_pos_s, z->y_pos_e + 1, mm, err, chn, sc); + } + +} + void dedup_chains(overlap_region_alloc* ol) { uint64_t k, l, s, m, mm_k, mm_m, sf; int64_t sc, mm_sc, plus, minus; overlap_region *z, t; @@ -3735,13 +3816,21 @@ void init_gen_hc_aln_t(gen_hc_aln_t *ez, overlap_region_alloc *ol, Candidates_li ez->hpz = hpz; } +uint64_t cal_aln_bs(overlap_region_alloc *ol) +{ + uint64_t tot = 0, k; + for (k = 0; k < ol->length; k++) { + tot += ol->list[k].x_pos_e + 1 - ol->list[k].x_pos_s; + } + return tot; +} static void worker_hap_ec(void *data, long i, int tid) { ec_ovec_buf_t0 *b = &(((ec_ovec_buf_t*)data)->a[tid]); uint32_t high_occ = asm_opt.hom_cov * (2.0 - HA_KMER_GOOD_RATIO); uint32_t low_occ = asm_opt.hom_cov * HA_KMER_GOOD_RATIO; ///gen_hc_aln_t ez; - overlap_region *aux_o = NULL/**, *rse_o = NULL**/; asg64_v buf0; uint32_t qlen = 0, qw = 0; + overlap_region *aux_o = NULL/**, *rse_o = NULL**/; asg64_v buf0; uint32_t qlen = 0, qw = 0; uint64_t tot_b = 0; double tt0 = 0, tt1 = 0; b->v8q.n = b->v8t.n = 0; // if((i != 733166) && (i != 858708) && (i != 858732) && (i != 859819) && (i != 859899) && (i != 863486) && (i != 872165) && (i != 899887) && (i != 902298) && @@ -3804,21 +3893,33 @@ static void worker_hap_ec(void *data, long i, int tid) // return; // } + if(DBG_TIME && dbg_a) { + dbg_a[i].chn_tm = dbg_a[i].aln_tm = dbg_a[i].phs_tm = dbg_a[i].cns_tm = 0; + } + // debug_retrive_bqual(D, &b->v8t, i, 256); return; - + if(DBG_TIME && dbg_a) { + tt0 = yak_realtime_0(); + } recover_UC_Read(&b->self_read, &R_INF, i); qlen = b->self_read.length; qw = ((qlen < (COV_W_AC<<1))?(qlen>>1):(COV_W_AC)); if(!qw) qw = 1; // if(qlen <= 0) return; - + h_ec_lchain(b->ab, i, b->self_read.seq, b->self_read.length, asm_opt.mz_win, asm_opt.k_mer_length, &R_INF, &b->olist, &b->clist, ((asm_opt.is_ont)?(0.05):(0.02)), asm_opt.max_n_chain, 1, NULL, NULL, &(b->sp), &high_occ, &low_occ, 1, 1, 3, 0.7, 2, 32, COV_W, 1);///ONT high error - + /** + h_ec_lchain(b->ab, i, b->self_read.seq, b->self_read.length, asm_opt.mz_win, asm_opt.k_mer_length, &R_INF, &b->olist, &b->clist, ((asm_opt.is_ont)?(0.05):(0.02)), asm_opt.max_n_chain, 1, NULL, NULL, &(b->sp), &high_occ, &low_occ, 1, 1, 3, 0.7, 2, 32, COV_W, 0);///ONT high error + **/ // b->num_read_base += b->olist.length; b->cnt[0] += b->self_read.length; aux_o = fetch_aux_ovlp(&b->olist, NULL/**&rse_o**/);///must be here - + if(DBG_TIME && dbg_a) { + tt1 = yak_realtime_0(); + dbg_a[i].chn_tm = tt1 - tt0; + tt0 = tt1; + } // stderr_phase_ovlp(&b->olist); @@ -3829,33 +3930,43 @@ static void worker_hap_ec(void *data, long i, int tid) // fprintf(stderr, "\n+[M::%s]\trid::%ld\t%.*s\n", __func__, i, (int)Get_NAME_LENGTH(R_INF, i), Get_NAME(R_INF, i)); - + ///r769: kp (gen_hc_r_alin_ea) -> NULL; site_sc (rphase_hc) -> 0 ///r770: kp (gen_hc_r_alin_ea) -> buf0; site_sc (rphase_hc) -> 0 ///r789: kp (gen_hc_r_alin_ea) -> NULL; site_sc (rphase_hc) -> 0 ///r791: kp (gen_hc_r_alin_ea) -> buf0; site_sc (rphase_hc) -> 0 copy_asg_arr(buf0, b->sp); - gen_hc_r_alin_ea_flt(b->ab, &b->olist, &b->clist, &R_INF, &b->self_read, &b->ovlp_read, &b->exz, aux_o, asm_opt.max_n_chain, asm_opt.max_n_chain*HC_MF_R, asm_opt.chn_occ, asm_opt.max_ov_diff_ec, (asm_opt.is_ont)?(WINDOW_OHC):(WINDOW_HC), i, E_KHIT, + tot_b = gen_hc_r_alin_ea_flt(b->ab, &b->olist, &b->clist, &R_INF, &b->self_read, &b->ovlp_read, &b->exz, aux_o, asm_opt.max_n_chain, asm_opt.max_n_chain*HC_MF_R, asm_opt.chn_occ, asm_opt.max_ov_diff_ec, (asm_opt.is_ont)?(WINDOW_OHC):(WINDOW_HC), i, E_KHIT, 1, &b->v16, &b->v64, &(R_INF.paf[i]), asm_opt.is_ont, (asm_opt.is_ont)?(0.006):(-1), (asm_opt.is_ont)?(64):(-1), (asm_opt.is_ont)?(512):(0), (asm_opt.is_ont)?(6):(0), - (asm_opt.is_ont)?(1.5):(-1), (asm_opt.is_ont)?(0.1):(-1), &buf0, qw, &b->v8q, &b->v32); + (asm_opt.is_ont)?(1.5):(-1), (asm_opt.is_ont)?(0.1):(-1), &buf0, qw, &b->v8q, &b->v32, 1); copy_asg_arr(b->sp, buf0); /** copy_asg_arr(buf0, b->sp); //kp: r763 -> r765: buf0 -> NULL //kp: r766 -> r767: NULL -> buf0 - gen_hc_r_alin_ea(&b->olist, &b->clist, &R_INF, &b->self_read, &b->ovlp_read, &b->exz, aux_o, asm_opt.max_ov_diff_ec, (asm_opt.is_ont)?(WINDOW_OHC):(WINDOW_HC), i, E_KHIT, + tot_b = gen_hc_r_alin_ea(&b->olist, &b->clist, &R_INF, &b->self_read, &b->ovlp_read, &b->exz, aux_o, asm_opt.max_ov_diff_ec, (asm_opt.is_ont)?(WINDOW_OHC):(WINDOW_HC), i, E_KHIT, 1, &b->v16, &b->v64, &(R_INF.paf[i]), asm_opt.is_ont, (asm_opt.is_ont)?(0.006):(-1), (asm_opt.is_ont)?(64):(-1), (asm_opt.is_ont)?(512):(0), (asm_opt.is_ont)?(6):(0), (asm_opt.is_ont)?(1.5):(-1), (asm_opt.is_ont)?(0.1):(-1), (asm_opt.is_ont)?(&buf0):(NULL), &b->v8q); - + copy_asg_arr(b->sp, buf0); + **/ // init_gen_hc_aln_t(&ez, &b->olist, &b->clist, &R_INF, &b->self_read, &b->ovlp_read, &b->exz, aux_o, // asm_opt.max_ov_diff_ec, asm_opt.max_ov_diff_ec, (asm_opt.is_ont)?(WINDOW_OHC):(WINDOW_HC), (asm_opt.is_ont)?(WINDOW_OHC):(WINDOW_HC), i, E_KHIT, 1, &b->v16, &b->v64, &(R_INF.paf[i]), // asm_opt.is_ont, asm_opt.is_ont, (asm_opt.is_ont)?(0.006):(-1), (asm_opt.is_ont)?(0.006):(-1), (asm_opt.is_ont)?(64):(-1), (asm_opt.is_ont)?(64):(-1), // (asm_opt.is_ont)?(512):(0), (asm_opt.is_ont)?(6):(0), (asm_opt.is_ont)?(1.5):(-1), (asm_opt.is_ont)?(0.1):(-1), (asm_opt.is_ont)?(&buf0):(NULL), (uint64_t)-1); // gen_hc_r_alin_ea_adv(&ez); - copy_asg_arr(b->sp, buf0); - **/ + + + + if(DBG_TIME && dbg_a) { + tt1 = yak_realtime_0(); + dbg_a[i].aln_tm = tt1 - tt0; + tt0 = tt1; + + dbg_a[i].faln = cal_aln_bs(&b->olist); + dbg_a[i].fbs = tot_b; + } // fprintf(stderr, "-[M::%s] rid::%ld\n", __func__, i); //for debug indel @@ -3883,6 +3994,12 @@ static void worker_hap_ec(void *data, long i, int tid) ///for debug indel // stderr_phase_ovlp(&b->olist); + if(DBG_TIME && dbg_a) { + tt1 = yak_realtime_0(); + dbg_a[i].phs_tm = tt1 - tt0; + tt0 = tt1; + } + dedup_chains(&b->olist); @@ -3892,6 +4009,12 @@ static void worker_hap_ec(void *data, long i, int tid) b->cnt[1] += wcns_gen(&b->olist, &R_INF, i, &b->self_read, &b->ovlp_read, &b->exz, &b->pidx, &b->v64, &buf0, 0, 512, b->self_read.length, 3, 0.500001, aux_o, &b->v32, &b->cns, 256, i, ((uint64_t)-1), NULL); copy_asg_arr(b->sp, buf0); + if(DBG_TIME && dbg_a) { + tt1 = yak_realtime_0(); + dbg_a[i].cns_tm = tt1 - tt0; + tt0 = tt1; + } + push_nec_re(aux_o, &(scc.a[i])); push_nec_re(aux_o, &(scb.a[i])); @@ -3981,6 +4104,132 @@ static void worker_hap_ec(void *data, long i, int tid) //fprintf(stderr, "-[M::%s]\trid::%ld\t%.*s\n", __func__, i, (int)Get_NAME_LENGTH(R_INF, i), Get_NAME(R_INF, i)); } + + +static void worker_hap_ec_ss(void *data, long i, int tid) +{ + ec_ovec_buf_t0 *b = &(((ec_ovec_buf_t*)data)->a[tid]); + uint32_t high_occ = asm_opt.hom_cov * (2.0 - HA_KMER_GOOD_RATIO); + uint32_t low_occ = asm_opt.hom_cov * HA_KMER_GOOD_RATIO; ///gen_hc_aln_t ez; + overlap_region *aux_o = NULL/**, *rse_o = NULL**/; asg64_v buf0; uint32_t qlen = 0, qw = 0; uint64_t /**tot_b = 0,**/ i0 = i, prt_n0; + b->v8q.n = b->v8t.n = 0; + + // if(((uint64_t)i) < dbgss->fn) return; + i = (uint32_t)dbgss->fa[i0]; prt_n0 = dbgss->spt_mul[tid].n; + + // if (memcmp("4e144e93-4653-4ebf-8920-7943e378cf9a", Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0) { + // fprintf(stderr, "-a-[M::%s-beg] rid->%ld\n", __func__, i); + // } else { + // return; + // } + + // if(i != 6814) return; + + // fprintf(stderr, "\n[M::%s] rid::%ld\t%.*s\tlen::%lld\n", __func__, i, (int)Get_NAME_LENGTH(R_INF, i), Get_NAME(R_INF, i), b->self_read.length); + + + recover_UC_Read(&b->self_read, &R_INF, i); qlen = b->self_read.length; + qw = ((qlen < (COV_W_AC<<1))?(qlen>>1):(COV_W_AC)); if(!qw) qw = 1; + // if(qlen <= 0) return; + + + ////new version + aux_o = NULL; ///tot_b = 0; + h_ec_lchain(b->ab, i, b->self_read.seq, b->self_read.length, asm_opt.mz_win, asm_opt.k_mer_length, &R_INF, &b->olist, &b->clist, ((asm_opt.is_ont)?(0.05):(0.02)), asm_opt.max_n_chain, 1, NULL, NULL, &(b->sp), &high_occ, &low_occ, 1, 1, 3, 0.7, 2, 32, COV_W, 1);///ONT high error + aux_o = fetch_aux_ovlp(&b->olist, NULL/**&rse_o**/);///must be here + + copy_asg_arr(buf0, b->sp); + /**tot_b =**/ gen_hc_r_alin_ea_flt(b->ab, &b->olist, &b->clist, &R_INF, &b->self_read, &b->ovlp_read, &b->exz, aux_o, asm_opt.max_n_chain, asm_opt.max_n_chain*HC_MF_R, asm_opt.chn_occ, asm_opt.max_ov_diff_ec, (asm_opt.is_ont)?(WINDOW_OHC):(WINDOW_HC), i, E_KHIT, + 1, &b->v16, &b->v64, &(R_INF.paf[i]), asm_opt.is_ont, (asm_opt.is_ont)?(0.006):(-1), (asm_opt.is_ont)?(64):(-1), (asm_opt.is_ont)?(512):(0), (asm_opt.is_ont)?(6):(0), + (asm_opt.is_ont)?(1.5):(-1), (asm_opt.is_ont)?(0.1):(-1), &buf0, qw, &b->v8q, &b->v32, 0); + copy_asg_arr(b->sp, buf0); + + stderr_phase_ovlp_buf(i0, &b->olist, &b->clist, b->v32.a + b->v32.n - b->olist.length - b->olist.length, b->v32.a + b->v32.n - b->olist.length, &(dbgss->spt_mul[tid]), &(dbgss->ssa[i]), "++"); + + + + ////old version + aux_o = NULL; ///tot_b = 0; + h_ec_lchain(b->ab, i, b->self_read.seq, b->self_read.length, asm_opt.mz_win, asm_opt.k_mer_length, &R_INF, &b->olist, &b->clist, ((asm_opt.is_ont)?(0.05):(0.02)), asm_opt.max_n_chain, 1, NULL, NULL, &(b->sp), &high_occ, &low_occ, 1, 1, 3, 0.7, 2, 32, COV_W, 0);///ONT high error + + aux_o = fetch_aux_ovlp(&b->olist, NULL);///must be here + copy_asg_arr(buf0, b->sp); + /**tot_b =**/ gen_hc_r_alin_ea(&b->olist, &b->clist, &R_INF, &b->self_read, &b->ovlp_read, &b->exz, aux_o, asm_opt.max_ov_diff_ec, (asm_opt.is_ont)?(WINDOW_OHC):(WINDOW_HC), i, E_KHIT, + 1, &b->v16, &b->v64, &(R_INF.paf[i]), asm_opt.is_ont, (asm_opt.is_ont)?(0.006):(-1), (asm_opt.is_ont)?(64):(-1), (asm_opt.is_ont)?(512):(0), (asm_opt.is_ont)?(6):(0), + (asm_opt.is_ont)?(1.5):(-1), (asm_opt.is_ont)?(0.1):(-1), (asm_opt.is_ont)?(&buf0):(NULL), &b->v8q); + copy_asg_arr(b->sp, buf0); + + + stderr_phase_ovlp_buf(i0, &b->olist, &b->clist, NULL, NULL, &(dbgss->spt_mul[tid]), &(dbgss->ssb[i]), "--"); + + if(dbgss->spt_mul[tid].n > prt_n0) { + dbgss->fa[i0] = prt_n0; + dbgss->fthr.a[i0] = tid; + kv_push(char, dbgss->spt_mul[tid], '\0'); + } else { + dbgss->fa[i0] = (uint64_t)-1; + dbgss->fthr.a[i0] = tid; + } + + return; + + + + b->cnt[0] += b->self_read.length; + + + // fprintf(stderr, "-[M::%s] rid::%ld\n", __func__, i); + //for debug indel + // prt_ovlp_sam(&b->olist, &b->ovlp_read, b->self_read.seq, b->self_read.length); + // stderr_phase_ovlp(&b->olist); + + + // fprintf(stderr, "\n[M::%s] rid::%ld\t%.*s\tlen::%lld\tocc::%lu\n", __func__, i, (int)Get_NAME_LENGTH(R_INF, i), Get_NAME(R_INF, i), b->self_read.length, b->olist.length); + + // fprintf(stderr, "[M::%s] rid::%ld\n", __func__, i); + // debug_mm_exact_cigar(&b->olist, i, &b->self_read, &b->ovlp_read); + + // b->num_correct_base += b->olist.length; + /** + * ///r779: enable this + copy_asg_arr(buf0, b->sp); + gen_reseed_re(&b->olist, &b->clist, aux_o, rse_o, &R_INF, &b->self_read, &b->ovlp_read, &b->exz, &b->pidx, &b->v64, &buf0, 0, asm_opt.mz_win, 19, i, asm_opt.max_ov_diff_ec, asm_opt.max_ov_diff_ec, &b->v16, R_INF.tqn, b->v8q.a); + copy_asg_arr(b->sp, buf0); + **/ + + 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); + copy_asg_arr(b->sp, buf0); + ///for debug indel + // stderr_phase_ovlp(&b->olist); + + + dedup_chains(&b->olist); + + + + copy_asg_arr(buf0, b->sp); + b->cnt[1] += wcns_gen(&b->olist, &R_INF, i, &b->self_read, &b->ovlp_read, &b->exz, &b->pidx, &b->v64, &buf0, 0, 512, b->self_read.length, 3, 0.500001, aux_o, &b->v32, &b->cns, 256, i, ((uint64_t)-1), NULL); + copy_asg_arr(b->sp, buf0); + + + push_nec_re(aux_o, &(scc.a[i])); + push_nec_re(aux_o, &(scb.a[i])); + + + push_ne_ovlp(&(R_INF.paf[i]), &b->olist, 1, &R_INF, &(scc.a[i])/**, i, &b->self_read, &b->ovlp_read**/); + push_ne_ovlp(&(R_INF.reverse_paf[i]), &b->olist, 2, &R_INF, NULL/**, i, NULL, NULL**/); + + + check_well_cal(&(scc.a[i]), &b->v64, &(R_INF.paf[i].is_fully_corrected), &(R_INF.paf[i].is_abnormal), qlen, (MIN_COVERAGE_THRESHOLD*2), &(R_INF.paf[i])); + R_INF.trio_flag[i] = AMBIGU; + + refresh_ec_ovec_buf_t0(b, REFRESH_N); + +} + + static void worker_hap_ec_hybrid(void *data, long i, int tid) { ec_ovec_buf_t0 *b = &(((ec_ovec_buf_t*)data)->a[tid]); @@ -6775,6 +7024,165 @@ static void worker_sl_ec(void *data, long i, int tid) } } +dbg_cnt_ss* fetch_dbg_cnt_ss(char *in, int64_t n_a, dbg_cnt_ss *in_p) +{ + FILE *fp = fopen(in, "r"); + if (!fp) { + fprintf(stderr, "No %s.", in); + return NULL; + } + + char buffer[8192], *pch = NULL; int64_t k = 0; + dbg_cnt_ss *p = in_p; + if(!p) MALLOC(p, n_a); + + while (fgets(buffer, sizeof(buffer), fp)) { + pch = strtok(buffer, "\t"); + while (pch != NULL) { + if((strlen(pch) >= 6) && (!memcmp(pch, "rid::", 5))) { + assert(k == atoll(pch+5)); + } else if(strlen(pch) >= 9) { + if(!memcmp(pch, "tot_bs::", 8)) { + p[k].fbs = atoll(pch+8); + } else if(!memcmp(pch, "aln_bs::", 8)) { + p[k].faln = atoll(pch+8); + } else if(!memcmp(pch, "chn_tm::", 8)) { + p[k].chn_tm = atof(pch+8); + } else if(!memcmp(pch, "aln_tm::", 8)) { + p[k].aln_tm = atof(pch+8); + } else if(!memcmp(pch, "phs_tm::", 8)) { + p[k].phs_tm = atof(pch+8); + } else if(!memcmp(pch, "cns_tm::", 8)) { + p[k].cns_tm = atof(pch+8); + } + } + pch = strtok (NULL, "\t"); + } + k++; + } + assert(n_a == k); + + fclose(fp); + return p; +} + +void prt_dbg_stats(dbg_cnt_ss *p, uint64_t n_a, char *fn_n, uint64_t rr) +{ + uint64_t k; + char* ga_n = (char*)malloc(strlen(fn_n)+64); + sprintf(ga_n, "%s.%lu.run.stat.log", fn_n, rr); + FILE *fn = fopen(ga_n, "w"); + free(ga_n); + + for (k = 0; k < n_a; k++) { + fprintf(fn, "rid::%lu\ttot_bs::%lu\taln_bs::%lu\tchn_tm::%f\taln_tm::%f\tphs_tm::%f\tcns_tm::%f\n", k, p[k].fbs, p[k].faln, p[k].chn_tm, p[k].aln_tm, p[k].phs_tm, p[k].cns_tm); + } + + fclose(fn); +} + + +dbg_cmp_ss* init_dbgss(uint64_t n_thre, uint64_t n_a, dbg_cmp_ss *in) +{ + dbg_cmp_ss *p = NULL; uint64_t k; + + if(!in) { + CALLOC(p, 1); + p->nrid = n_a; CALLOC(p->ssa, n_a); CALLOC(p->ssb, n_a); + p->mul_n = n_thre; CALLOC(p->spt_mul, n_thre); + // p->srt_aa_cut = p->srt_ba_cut = 0; + } else { + p = in; + free(p->ssa); free(p->ssb); + free(p->srt_aa.a); free(p->srt_ba.a); + for (k = 0; k < p->mul_n; k++) { + free(p->spt_mul[k].a); + } + free(p->spt_mul); free(p->fthr.a); free(p); + p = NULL; + } + + return p; +} + +void prt_dbgss_cmd(asgchr_v *buf_a, uint64_t *a, uint64_t an, uint32_t *tid_a, const char *fn_n) +{ + uint64_t k; + char* ga_n = (char*)malloc(strlen(fn_n)+64); + sprintf(ga_n, "%s.cmp.diff.stat.log", fn_n); + FILE *fn = fopen(ga_n, "w"); + free(ga_n); + + for (k = 0; k < an; k++) { + if(a[k] == ((uint64_t)-1)) continue; + // fprintf(stderr, "k::%lu, tid::%u, t_aid::%lu\n", k, tid_a[k], a[k]); + // fprintf(stderr, "t_aid_n::%lu\n", buf_a[tid_a[k]].n); + fputs(buf_a[tid_a[k]].a + a[k], fn); + } + + fclose(fn); +} + +void cal_ec_multiple_stat_cmp(ec_ovec_buf_t *b, uint64_t n_thre, uint64_t n_a, char *stat_1, char *stat_2, double ssa_cut, double ssb_cut) +{ + uint64_t k = 0, sp, tot_a, tot_b, tt, tcut, *sa = NULL, sn; uint8_t sf; + dbgss = init_dbgss(n_thre, n_a, NULL); + + fetch_dbg_cnt_ss(stat_1, n_a, dbgss->ssa); + fetch_dbg_cnt_ss(stat_2, n_a, dbgss->ssb); + + // prt_dbg_stats(dbg_a, n_a, (char *)"dbg_a", 768); + // prt_dbg_stats(dbg_b, n_a, (char *)"dbg_b", 768); + for (k = dbgss->srt_aa.n = dbgss->srt_ba.n = tot_a = tot_b = 0; k < n_a; k++) { + if(dbgss->ssa[k].fbs >= dbgss->ssb[k].fbs) { + sp = dbgss->ssa[k].fbs - dbgss->ssb[k].fbs; sf = 1; + tot_a += sp; + } else { + sp = dbgss->ssb[k].fbs - dbgss->ssa[k].fbs; sf = 0; + tot_b += sp; + } + if(sp > UINT32_MAX) sp = UINT32_MAX; + sp = UINT32_MAX - sp; sp <<= 32; sp |= k; + if(sf) { + kv_push(uint64_t, dbgss->srt_aa, sp); + } else { + kv_push(uint64_t, dbgss->srt_ba, sp); + } + } + + radix_sort_ec64(dbgss->srt_aa.a, dbgss->srt_aa.a + dbgss->srt_aa.n); + radix_sort_ec64(dbgss->srt_ba.a, dbgss->srt_ba.a + dbgss->srt_ba.n); + + + + for (k = 0; k < dbgss->mul_n; k++) dbgss->spt_mul[k].n = 0; + sa = dbgss->srt_aa.a; sn = dbgss->srt_aa.n; tcut = tot_a*ssa_cut; + for (k = tt = 0; k < sn && tt <= tcut; k++) { + tt += (UINT32_MAX - (sa[k]>>32)); + } + dbgss->fa = sa; dbgss->fn = k; kv_resize(uint32_t, dbgss->fthr, dbgss->fn); + fprintf(stderr, "[M::a>=b] # top->%f diff bases::%lu(# read::%lu); # tot diff bases::%lu; # reads::%lu\n", ssa_cut, tt, k, tot_a, n_a); + kt_for(n_thre, worker_hap_ec_ss, b, dbgss->fn);///debug_for_fix + prt_dbgss_cmd(dbgss->spt_mul, sa, k, dbgss->fthr.a, "a_b"); + + + + for (k = 0; k < dbgss->mul_n; k++) dbgss->spt_mul[k].n = 0; + sa = dbgss->srt_ba.a; sn = dbgss->srt_ba.n; tcut = tot_b*ssb_cut; + for (k = tt = 0; k < sn && tt <= tcut; k++) { + tt += (UINT32_MAX - (sa[k]>>32)); + } + dbgss->fa = sa; dbgss->fn = k; kv_resize(uint32_t, dbgss->fthr, dbgss->fn); + fprintf(stderr, "[M:::a%f diff bases::%lu(# read::%lu); # tot diff bases::%lu; # reads::%lu\n", ssb_cut, tt, k, tot_b, n_a); + kt_for(n_thre, worker_hap_ec_ss, b, dbgss->fn);///debug_for_fix + prt_dbgss_cmd(dbgss->spt_mul, sa, k, dbgss->fthr.a, "b_a"); + + + + init_dbgss(n_thre, n_a, dbgss); + exit(1); +} + uint64_t cal_ec_multiple(ec_ovec_buf_t *b, uint64_t n_thre, uint64_t n_a, uint64_t *r_base) { double tt0 = yak_realtime_0(); @@ -6790,6 +7198,9 @@ uint64_t cal_ec_multiple(ec_ovec_buf_t *b, uint64_t n_thre, uint64_t n_a, uint64 for (k = 0; k < n_thre; ++k) b->a[k].cnt[0] = b->a[k].cnt[1] = 0; + // fprintf(stderr, "[M::%s] n_thre->%lu\n", __func__, n_thre); + if(asm_opt.dbg_run_1 && asm_opt.dbg_run_2) cal_ec_multiple_stat_cmp(b, n_thre, n_a, asm_opt.dbg_run_1, asm_opt.dbg_run_2, 0.1, 0.1); + if(!(asm_opt.hf)) kt_for(n_thre, worker_hap_ec, b, n_a);///debug_for_fix else kt_for(n_thre, worker_hap_ec_hybrid, b, n_a);///debug_for_fix @@ -6980,11 +7391,11 @@ void write_ec_reads(const char *suffix_ou) fclose(ou); destory_UC_Read(&qstr); destory_UC_Read(&tstr); } -// dbg_cnt_ss* gen_dbg_cnt_ss(uint64_t n_a) -// { -// dbg_cnt_ss *p = NULL; CALLOC(p, n_a); -// return p; -// } +dbg_cnt_ss* gen_dbg_cnt_ss(uint64_t n_a) +{ + dbg_cnt_ss *p = NULL; CALLOC(p, n_a); + return p; +} 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) { @@ -6996,9 +7407,15 @@ void cal_ec_r(uint64_t n_thre, uint64_t round, uint64_t n_round, uint64_t n_a, u ec_ovec_buf_t *b = NULL; uint64_t k, is_cr = (round&1); (*tot_b) = (*tot_e) = 0; - // if(DBG_TIME) dbg_a = gen_dbg_cnt_ss(n_a); + if(DBG_TIME && ((!asm_opt.dbg_run_1) || (!asm_opt.dbg_run_2))) { + dbg_a = gen_dbg_cnt_ss(n_a); + } b = gen_ec_ovec_buf_t(n_thre); (*tot_e) += cal_ec_multiple(b, n_thre, n_a, tot_b); ///exit(1); + if(DBG_TIME && dbg_a) { + prt_dbg_stats(dbg_a, n_a, asm_opt.output_file_name, round); + free(dbg_a); dbg_a = NULL; + } sl_ec_r(n_thre, n_a); for (k = 0; k < n_round; k++) {