r878

CommandLines.h

@@ -5,7 +5,7 @@
 #include <pthread.h>
 #include <stdint.h>
 
-#define HA_VERSION "0.25.0-r877"
+#define HA_VERSION "0.25.1-r878"
 
 #define VERBOSE 0
 

Correct.cpp

@@ -35989,6 +35989,192 @@ 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_adv_adp_smp_0(gen_hc_aln_t *ez, uint8_t set_match)
+{
+    uint64_t i, bs, k, w, ql = ez->qu->length, *wsrt = NULL, wsrt_n = 0, spn0 = 0, tot_b = 0; double err, e_max, align_gap_max, align_gap_rate, chem_drop; overlap_region t, *z; asg64_v *sp = ez->srt;
+    ez->ol->mapped_overlaps_length = 0; uint32_t *ocn = ez->v32->a, *osc = ez->v32->a + ez->ol->length; 
+    if(ez->ol->length <= 0) return;
+
+    ///base alignment
+    bs = (MAX(ez->wl[0], ez->wl[1]))+(THRESHOLD_MAX_SIZE<<1)+1;
+    resize_UC_Read(ez->tu, bs<<1); spn0 = sp->n; 
+
+    wsrt = mmp_chn_select(ez->ol, ez->cl, sp, ez->ocw, ocn, osc, ql, &wsrt_n, set_match, ez->max_n_chain, ez->max_n_chain_f, ez->chain_cutoff, ez->ave_cov_min, 0.333333, 16, 16);
+
+    for (i = 0; i < wsrt_n; i++) {
+        z = &(ez->ol->list[(uint32_t)wsrt[i]]); 
+
+        if(z->y_id < ez->t_cut) {
+            err = ez->e_rate[1]; w = ez->wl[1]; align_gap_max = ez->align_gap_max[1]; align_gap_rate = ez->align_gap_rate[1]; chem_drop = ez->chem_drop[1];
+        } else {
+            err = ez->e_rate[0]; w = ez->wl[0]; align_gap_max = ez->align_gap_max[0]; align_gap_rate = ez->align_gap_rate[0]; chem_drop = ez->chem_drop[0];
+        }
+        e_max = err * 1.5;
+
+
+        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, NULL/**ez->hpz->a**/, ez->buf, 0, &tot_b))) {
+            continue;
+        }
+
+        z->is_match = 1; z->strong = z->without_large_indel = 0;
+    }
+
+    for (i = k = 0; i < ez->ol->length; i++) {///primary chain
+        z = &(ez->ol->list[i]); 
+        if(z->is_match == 0) continue;
+        if(k != i) {
+            t = ez->ol->list[k];
+            ez->ol->list[k] = ez->ol->list[i];
+            ez->ol->list[i] = t;
+        }
+        k++;
+    }
+
+    ez->ol->length = k; sp->n = spn0; 
+}
+
+void batch_simd_ck(gen_hc_aln_t *ez, uint64_t *wsrt, uint64_t wsrt_n, uint64_t ql, double err, uint64_t w, double align_gap_max, double align_gap_rate, double chem_drop)
+{
+    if(!wsrt_n) return;
+    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];
+    uint64_t i, k, zk, in0, s, e, q[2], os, oe; overlap_region *z;
+
+    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)(ez->ol->list[(uint32_t)wsrt[zk]].x_pos_e + 1))<<32)|((uint32_t)wsrt[zk]);
+                }
+                radix_sort_bc64(wsrt + i, wsrt + k);
+            }
+            i = k;
+        }
+    }
+
+    ez->v32->n = (ez->ol->length<<1); in0 = ez->v32->n;
+    i = 0; s = 0; e = w; e = ((e<=ql)?e:ql); rr = 0;
+    for (; s < ql; ) {///[s, e)
+        if(rr) {
+            for (k = zk = in0; k < ez->v32->n; k++) {
+                if((wsrt[ez->v32->a[k]]>>32) == UINT32_MAX) continue;///passed
+                z = &(ez->ol->list[(uint32_t)wsrt[ez->v32->a[k]]]);
+                q[0] = z->x_pos_s;
+                q[1] = z->x_pos_e + 1;
+                ///[s, e) && [q[0], q[1])
+                os = MAX(q[0], s); oe = MIN(q[1], e);
+                if(oe > os) ez->v32->a[zk++] = ez->v32->a[k];
+            }
+            ez->v32->n = zk;
+        }
+
+        for (; i < wsrt_n; ++i) {
+            z = &(ez->ol->list[(uint32_t)wsrt[i]]);
+            q[0] = z->x_pos_s;
+            q[1] = z->x_pos_e + 1;
+            if(q[0] >= e) break;
+            if(z->is_match == 0) {
+                os = MAX(q[0], s); oe = MIN(q[1], e);
+                if(oe > os) kv_push(uint32_t, *(ez->v32), i);
+            }
+        }
+
+        rr = hc_aln_simd(ez->ol->list, wsrt, ez->v32->a + in0, ez->v32->n - in0, s, e, w, err, ez->rref, ez->qu->seq, ez->tu, fi, aux_beg, aux_end, t_s, t_pri_l, ez->exz, OVERLAP_THRESHOLD_HIFI_FILTER, 0, NULL/**&tot_b**/);            
+        s += w; e += w; e = ((e<=ql)?e:ql);
+    }
+}
+
+void gen_hc_r_alin_adv_adp_smp_1(gen_hc_aln_t *ez, uint8_t set_match)
+{
+    uint64_t i, bs, k, w, ql = ez->qu->length, nol_1, *wsrt = NULL, wsrt_n = 0, spn0 = 0; double err, e_max, align_gap_max, align_gap_rate, chem_drop; overlap_region t, *z; asg64_v *sp = ez->srt;
+    ez->ol->mapped_overlaps_length = 0; uint32_t *ocn = ez->v32->a, *osc = ez->v32->a + ez->ol->length; 
+    if(ez->ol->length <= 0) return;
+
+    ///base alignment
+    bs = (MAX(ez->wl[0], ez->wl[1]))+(THRESHOLD_MAX_SIZE<<1)+1;
+    resize_UC_Read(ez->tu, bs<<1); spn0 = sp->n; 
+
+    wsrt = mmp_chn_select(ez->ol, ez->cl, sp, ez->ocw, ocn, osc, ql, &wsrt_n, set_match, ez->max_n_chain, ez->max_n_chain_f, ez->chain_cutoff, ez->ave_cov_min, 0.333333, 16, 16);
+
+
+    for (i = nol_1 = 0; i < wsrt_n; i++) {
+        z = &(ez->ol->list[(uint32_t)wsrt[i]]);
+        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]);
+
+        if(z->y_id < ez->t_cut) {
+            if(nol_1 != i) {
+                bs = wsrt[i]; wsrt[i] = wsrt[nol_1]; wsrt[nol_1] = bs;
+            }
+            nol_1++;
+        }  
+    }
+
+    batch_simd_ck(ez, wsrt, nol_1, ql, ez->e_rate[1], ez->wl[1], ez->align_gap_max[1], ez->align_gap_rate[1], ez->chem_drop[1]);
+    batch_simd_ck(ez, wsrt + nol_1, wsrt_n - nol_1, ql, ez->e_rate[0], ez->wl[0], ez->align_gap_max[0], ez->align_gap_rate[0], ez->chem_drop[0]);
+    ocn = ez->v32->a; osc = ez->v32->a + ez->ol->length;
+
+    for (i = 0; i < wsrt_n; i++) {
+        if((wsrt[i]>>32) == UINT32_MAX) {
+            continue;
+        }
+        z = &(ez->ol->list[(uint32_t)wsrt[i]]); 
+
+        if(z->y_id < ez->t_cut) {
+            err = ez->e_rate[1]; w = ez->wl[1]; align_gap_max = ez->align_gap_max[1]; align_gap_rate = ez->align_gap_rate[1]; chem_drop = ez->chem_drop[1];
+        } else {
+            err = ez->e_rate[0]; w = ez->wl[0]; align_gap_max = ez->align_gap_max[0]; align_gap_rate = ez->align_gap_rate[0]; chem_drop = ez->chem_drop[0];
+        }
+        e_max = err * 1.5;
+
+        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, 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, NULL/**ez->hpz->a**/, ez->buf, 1, NULL/**&tot_b**/))) {
+                continue;
+            }
+        }
+        // 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, NULL/**ez->hpz->a**/, ez->buf, 0, &tot_b))) {
+        //     continue;
+        // }
+
+        z->is_match = 1; z->strong = z->without_large_indel = 0;
+    }
+
+    /**
+    // uint64_t kqs, kqe, kts, kte, kerr;
+    for (i = 0; i < wsrt_n; i++) {
+        if((wsrt[i]>>32) == UINT32_MAX) {
+            continue;
+        }
+        z = &(ol->list[(uint32_t)wsrt[i]]); 
+
+        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))) {
+                continue;
+            }
+        }
+
+        z->is_match = 1; z->strong = z->without_large_indel = 0;
+    }
+    **/
+
+    for (i = k = 0; i < ez->ol->length; i++) {///primary chain
+        z = &(ez->ol->list[i]); 
+        if(z->is_match == 0) continue;
+        if(k != i) {
+            t = ez->ol->list[k];
+            ez->ol->list[k] = ez->ol->list[i];
+            ez->ol->list[i] = t;
+        }
+        k++;
+    }
+
+    ez->ol->length = k; sp->n = spn0; 
+}
+
 
 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)
 {

Correct.h

@@ -1442,6 +1442,8 @@ uint64_t gen_hc_r_alin_adp_mmp_1(overlap_region_alloc* ol, Candidates_list *cl,
 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 gen_hc_r_alin_adv_adp_smp_0(gen_hc_aln_t *ez, uint8_t set_match);
+void gen_hc_r_alin_adv_adp_smp_1(gen_hc_aln_t *ez, uint8_t set_match);
 void pp_chn_a(overlap_region *z, Candidates_list *cl, uint8_t is_raw);
 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, 

ecovlp.cpp

@@ -3686,6 +3686,73 @@ void gen_hc_r_alin_ea_adv_flt(gen_hc_aln_t *ez)
 }
 
 
+void gen_hc_r_alin_ea_adv_flt_mmp(gen_hc_aln_t *ez)
+{
+    if(ez->ol->length <= 0) return;
+
+    uint64_t i, k, m, *ei, en, *oi, on, tid, trev, nec; overlap_region *z; ma_hit_t *p;
+    ez->v32->n = ez->ol->length<<1; kv_resize(uint32_t, *(ez->v32), ez->v32->n); 
+    for (i = 0; i < ez->ol->length; i++) {
+        ez->v32->a[i] = ez->ol->list[i].align_length;
+        ez->v32->a[i+ez->ol->length] = ez->ol->list[i].shared_seed;
+        ez->ol->list[i].align_length = 0;
+    }
+
+    // uint64_t k, i, m, *ei, en, *oi, on, tid, trev, nec; overlap_region *z; ma_hit_t *p;
+    for (k = ez->srt->n = 0; k < ez->in->length; k++) {
+        if(ez->in->buffer[k].el) {
+            m = ez->in->buffer[k].tn; m <<= 1; m |= ez->in->buffer[k].rev; 
+            m <<= 32; m |= k; kv_push(uint64_t, (*(ez->srt)), m);
+        }
+    }
+
+    if(!(ez->srt->n)) {
+        // gen_hc_r_alin_adv_adp_smp(ez, a_cu, a_ci, ocn, osc, idx_cu, n_cu, 1);
+        gen_hc_r_alin_adv_adp_smp_1(ez, 1);
+    } 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);
+        kv_resize(uint64_t, *(ez->srt), (ez->srt->n + ez->ol->length));
+        ei = ez->srt->a; en = ez->srt->n; oi = ez->srt->a + ez->srt->n; on = ez->ol->length;
+        for (k = 0; k < on; k++) {
+            z = &(ez->ol->list[k]); z->is_match = z->strong = z->without_large_indel = 0;
+            oi[k] = z->y_id; oi[k] <<= 1; oi[k] |= z->y_pos_strand;
+            oi[k] <<= 32; oi[k] |= k;
+        }
+
+        radix_sort_ec64(ei, ei + en); radix_sort_ec64(oi, oi + on);
+        for (k = i = nec = 0; k < on; k++) {
+            z = &(ez->ol->list[(uint32_t)oi[k]]); tid = z->y_id; trev = z->y_pos_strand;
+            for (; (i < en) && ((ei[i]>>32) < ((tid<<1)|trev)); i++);
+            if((i < en) && ((ei[i]>>32) == ((tid<<1)|trev))) {
+                p = &(ez->in->buffer[(uint32_t)ei[i]]);
+                if((z->x_pos_s == ((uint32_t)p->qns)) && (z->x_pos_e + 1 == p->qe) && 
+                            (z->y_pos_s == p->ts) && (z->y_pos_e + 1 == p->te)) {
+                    resize_UC_Read(ez->tu, p->te - p->ts); recover_UC_Read_sub_region(ez->tu->seq, p->ts, p->te - p->ts, trev, ez->rref, tid);
+                    if(exact_ec_check(ez->qu->seq, ez->qu->length, ez->tu->seq, p->te - p->ts, ((uint32_t)p->qns), p->qe, 0, p->te - p->ts)) {
+                        z->is_match = 1; z->shared_seed = z->non_homopolymer_errors;///for index
+                        z->non_homopolymer_errors = 0; z->strong = z->without_large_indel = 0;
+                        set_exact_exz(ez->exz, z->x_pos_s, z->x_pos_e + 1, z->y_pos_s, z->y_pos_e + 1); push_alnw(z, ez->exz);
+                        nec++;
+                    }
+                }
+            }
+        }
+        ///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);
+
+        if(on > nec) {
+            // gen_hc_r_alin_adv_adp_smp(ez, a_cu, a_ci, ocn, osc, idx_cu, n_cu, 0);
+            gen_hc_r_alin_adv_adp_smp_1(ez, 0);
+        }
+        // fprintf(stderr, "[M::%s] srt->n::%u, nec::%lu, on::%lu\n", __func__, (uint32_t)srt->n, nec, on);
+        ///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);
+    }
+}
+
+
+
 void prt_ovlp_sam_0(char *cm, FILE *fp, char *ref_id, int32_t ref_id_n, char *qry_id, int32_t qry_id_n, char *qry_seq, uint64_t qry_seq_n, uint64_t rs, uint64_t re, uint64_t qs, uint64_t qe, uint64_t flag, uint64_t err0, bit_extz_t *ez)
 {
     uint64_t ci = 0, err1 = 0; uint16_t c; uint32_t cl, cl0 = 0; char c0 = (char)-1;
@@ -5129,14 +5196,15 @@ static void worker_hap_ec_hybrid(void *data, long i, int tid)
     aux_o = fetch_aux_ovlp(&b->olist, NULL/**&rse_o**/);///must be here
 
     copy_asg_arr(buf0, b->sp); 
-    init_gen_hc_aln_t(&ez, &b->olist, &b->clist, &R_INF, &b->self_read, &b->ovlp_read, &b->exz, aux_o, rse_o, &b->v8q,
+    init_gen_hc_aln_t(&ez, &b->olist, &b->clist, &R_INF, &b->self_read, &b->ovlp_read, &b->exz, aux_o, rse_o, /**&b->v8q**/NULL,
         e_l, e_h, (qid < R_INF.tqn)?(WINDOW_OHC):(WINDOW_HC), WINDOW_OHC, i, E_KHIT, 1, &b->v16, &b->v64, &(R_INF.paf[i]),
         ((qid < R_INF.tqn)?(1):(0)), 1, (qid < R_INF.tqn)?(0.006):(-1), 0.006, (qid < R_INF.tqn)?(64):(-1), 64, 
         (qid < R_INF.tqn)?(512):(0), (qid < R_INF.tqn)?(6):(0), (qid < R_INF.tqn)?(1.5):(-1), (qid < R_INF.tqn)?(0.1):(-1), NULL, 
         &(b->v32), &buf0, b->ab, (asm_opt.max_n_chain>0)?(asm_opt.max_n_chain):(1), ((asm_opt.max_n_chain*HC_MF_R)>0)?(asm_opt.max_n_chain*HC_MF_R):1, 
         asm_opt.chn_occ, ((asm_opt.hom_cov*HC_AV_MIN)>0)?(asm_opt.hom_cov*HC_AV_MIN):(1), qw, R_INF.tqn, asm_opt.hom_cov); 
         // gen_hc_r_alin_ea_adv(&ez);
-    gen_hc_r_alin_ea_adv_flt(&ez);
+    // gen_hc_r_alin_ea_adv_flt(&ez);
+    gen_hc_r_alin_ea_adv_flt_mmp(&ez);
     copy_asg_arr(b->sp, buf0);
 
     // fprintf(stderr, "-c-[M::%s] rid::%ld\n", __func__, i);