diff --git a/Correct.cpp b/Correct.cpp index 9c3776d..2fb0c9a 100644 --- a/Correct.cpp +++ b/Correct.cpp @@ -1787,7 +1787,7 @@ inline int move_gap_greedy(char* path, int path_i, int path_length, char* x, int char oper = path[path_i]; - if(oper == 3) + if(oper == 3)///there are more x { path_i++; y_i--; @@ -1817,7 +1817,7 @@ inline int move_gap_greedy(char* path, int path_i, int path_length, char* x, int } } - else if(oper == 2) + else if(oper == 2)///there are more y { path_i++; x_i--; @@ -1925,11 +1925,11 @@ inline void generate_cigar(char* path, int path_length, window_list *idx, window else if(path[i] == 1) { x_i++; y_i++; } - else if(path[i] == 2) { + else if(path[i] == 2) {///there are more y move_gap_greedy(path, i, path_length, x, x_i, y, y_i, old_error); y_i++; } - else if(path[i] == 3) { + else if(path[i] == 3) {///there are more x move_gap_greedy(path, i, path_length, x, x_i, y, y_i, old_error); x_i++; } @@ -3623,6 +3623,12 @@ int32_t y_strand, int32_t y_id) return 0; } +void gen_rev_str(char *in, char **out, uint32_t len) +{ + char *r; uint32_t k; MALLOC(r, len); (*out) = r; + for (k = 0; k < len; k++) r[k] = in[len-k-1]; +} + inline uint32_t gen_backtrace_adv(window_list *p, overlap_region *z, All_reads *rref, hpc_t *hpc_g, const ul_idx_t *uref, char *qstr, char *tstr, char *tstr1, Correct_dumy* dumy, uint32_t rev, uint32_t id) { @@ -3650,58 +3656,146 @@ char *qstr, char *tstr, char *tstr1, Correct_dumy* dumy, uint32_t rev, uint32_t // assert(error != (unsigned int)-1); if(error != (unsigned int)-1) { - // int32_t dbg_e = ed_band_cal_global(t_string+r_ts, t_end+1-r_ts, q_string, ql, thres); - // if(dbg_e != (int32_t)error && dbg_e <= 4) { - // fprintf(stderr, "\n[M::%s::] error::%u, ed_global::%d\n", __func__, error, dbg_e); - // fprintf(stderr, "[tstr] %.*s\n", t_end+1-r_ts, t_string+r_ts); - // fprintf(stderr, "[qstr] %.*s\n", (int32_t)ql, q_string); - // } - // assert(dbg_e <= (int32_t)error); + /** bit_extz_t exz, exz64; init_bit_extz_t(&exz, thres); init_bit_extz_t(&exz64, thres); - - // if(!(exz.err <= (int32_t)error)) { - // fprintf(stderr, "[M::%s::] error::%u, ed_extension::%d, ql::%ld, thres::%ld\n", - // __func__, error, exz.err, ql, thres); - - // } - // exit(1); - // - // ed_band_cal_extension_256_w(t_string+r_ts, t_end+1-r_ts, q_string, ql, thres, &exz); - ed_band_cal_extension_infi0_w(t_string+r_ts, t_end+1-r_ts, q_string, ql, thres, NULL, &exz); - ed_band_cal_extension_64_w(t_string+r_ts, t_end+1-r_ts, q_string, ql, thres, &exz64); - - if(!(exz.err == exz64.err && exz.ps == exz64.ps && exz.pe == exz64.pe && exz.ts == exz64.ts && exz.te == exz64.te)) { - fprintf(stderr, "\n[M::%s::] ql::%ld, thres::%ld, exz.err::%d, exz64.err::%d, exz.ps::%d, exz64.ps::%d, exz.pe::%d, exz64.pe::%d, exz.ts::%d, exz64.ts::%d, exz.te::%d, exz64.te::%d\n", __func__, - ql, thres, exz.err, exz64.err, exz.ps, exz64.ps, exz.pe, exz64.pe, exz.ts, exz64.ts, exz.te, exz64.te); - } + + ed_band_cal_extension_64_0_w(t_string+r_ts, t_end+1-r_ts, q_string, ql, thres, &exz64); + ed_band_cal_extension_64_0_w_trace(t_string+r_ts, t_end+1-r_ts, q_string, ql, thres, &exz64); + cigar_check(t_string+r_ts, q_string, &(exz64)); + exz64.err = INT32_MAX; + ed_band_cal_extension_64_0_w_trace(t_string+r_ts, t_end+1-r_ts, q_string, ql, thres, &exz64); + cigar_check(t_string+r_ts, q_string, &(exz64)); + + ed_band_cal_extension_infi_0_w(t_string+r_ts, t_end+1-r_ts, q_string, ql, thres, NULL, &exz); + ed_band_cal_extension_infi_0_w_trace(t_string+r_ts, t_end+1-r_ts, q_string, ql, thres, NULL, &exz); + cigar_check(t_string+r_ts, q_string, &exz); + exz.err = INT32_MAX; + ed_band_cal_extension_infi_0_w_trace(t_string+r_ts, t_end+1-r_ts, q_string, ql, thres, NULL, &exz); + cigar_check(t_string+r_ts, q_string, &exz); assert(exz.err <= (int32_t)error && exz.err >= 0); assert(exz.err == exz64.err && exz.ps == exz64.ps && exz.pe == exz64.pe && exz.ts == exz64.ts && exz.te == exz64.te); - ed_band_cal_extension_infi1_w(t_string+r_ts, t_end+1-r_ts, q_string, ql, thres, NULL, &exz); - assert(exz.err <= (int32_t)error); + ed_band_cal_extension_256_0_w(t_string+r_ts, t_end+1-r_ts, q_string, ql, thres, &exz); + ed_band_cal_extension_256_0_w_trace(t_string+r_ts, t_end+1-r_ts, q_string, ql, thres, &exz); + cigar_check(t_string+r_ts, q_string, &exz); + exz.err = INT32_MAX; + ed_band_cal_extension_256_0_w_trace(t_string+r_ts, t_end+1-r_ts, q_string, ql, thres, &exz); + cigar_check(t_string+r_ts, q_string, &exz); + assert(exz.err <= (int32_t)error && exz.err >= 0); + assert(exz.err == exz64.err && exz.ps == exz64.ps && exz.pe == exz64.pe && exz.ts == exz64.ts && exz.te == exz64.te); - // ed_band_cal_global_256_w(t_string+r_ts, t_end+1-r_ts, q_string, ql, thres, &exz); - ed_band_cal_global_infi_w(t_string+r_ts, t_end+1-r_ts, q_string, ql, thres, NULL, &exz); - ed_band_cal_global_64_w(t_string+r_ts, t_end+1-r_ts, q_string, ql, thres, &exz64); // if(!(exz.err == exz64.err && exz.ps == exz64.ps && exz.pe == exz64.pe && exz.ts == exz64.ts && exz.te == exz64.te)) { // fprintf(stderr, "\n[M::%s::] ql::%ld, thres::%ld, exz.err::%d, exz64.err::%d, exz.ps::%d, exz64.ps::%d, exz.pe::%d, exz64.pe::%d, exz.ts::%d, exz64.ts::%d, exz.te::%d, exz64.te::%d\n", __func__, // ql, thres, exz.err, exz64.err, exz.ps, exz64.ps, exz.pe, exz64.pe, exz.ts, exz64.ts, exz.te, exz64.te); // } + + + + + char *qr, *tr; + gen_rev_str(q_string, &qr, ql); gen_rev_str(t_string+r_ts, &tr, t_end+1-r_ts); + + ed_band_cal_extension_64_1_w(tr, t_end+1-r_ts, qr, ql, thres, &exz64); + ed_band_cal_extension_64_1_w_trace(tr, t_end+1-r_ts, qr, ql, thres, &exz64); + cigar_check(tr, qr, &exz64); + exz64.err = INT32_MAX; + ed_band_cal_extension_64_1_w_trace(tr, t_end+1-r_ts, qr, ql, thres, &exz64); + cigar_check(tr, qr, &exz64); + assert(exz.err == exz64.err && exz.ps == (exz64.pl-exz64.pe-1) && exz.pe == (exz64.pl-exz64.ps-1) && exz.ts == (exz64.tl-exz64.te-1) && exz.te == (exz64.tl-exz64.ts-1)); + + ed_band_cal_extension_infi_1_w(tr, t_end+1-r_ts, qr, ql, thres, NULL, &exz); + ed_band_cal_extension_infi_1_w_trace(tr, t_end+1-r_ts, qr, ql, thres, NULL, &exz); + cigar_check(tr, qr, &exz); + exz.err = INT32_MAX; + ed_band_cal_extension_infi_1_w_trace(tr, t_end+1-r_ts, qr, ql, thres, NULL, &exz); + cigar_check(tr, qr, &exz); + assert(exz.err == exz64.err && exz.ps == exz64.ps && exz.pe == exz64.pe && exz.ts == exz64.ts && exz.te == exz64.te); + + ed_band_cal_extension_256_1_w(tr, t_end+1-r_ts, qr, ql, thres, &exz); + ed_band_cal_extension_256_1_w_trace(tr, t_end+1-r_ts, qr, ql, thres, &exz); + cigar_check(tr, qr, &exz); + exz.err = INT32_MAX; + ed_band_cal_extension_256_1_w_trace(tr, t_end+1-r_ts, qr, ql, thres, &exz); + cigar_check(tr, qr, &exz); + assert(exz.err == exz64.err && exz.ps == exz64.ps && exz.pe == exz64.pe && exz.ts == exz64.ts && exz.te == exz64.te); + + free(qr); free(tr); + + + ed_band_cal_global_64_w(t_string+r_ts, t_end+1-r_ts, q_string, ql, thres, &exz64); + ed_band_cal_global_64_w_trace(t_string+r_ts, t_end+1-r_ts, q_string, ql, thres, &exz64); + cigar_check(t_string+r_ts, q_string, &(exz64)); + exz64.err = INT32_MAX; + ed_band_cal_global_64_w_trace(t_string+r_ts, t_end+1-r_ts, q_string, ql, thres, &exz64); + cigar_check(t_string+r_ts, q_string, &(exz64)); + + // if(!cigar_check(t_string+r_ts, t_end+1-r_ts, q_string, ql, &(exz64.cigar), exz64.err) && ql <= 12) { + // fprintf(stderr, "\n[M::%s::] ql::%ld, thres::%ld, exz.err::%d, exz64.err::%d, exz.ps::%d, exz64.ps::%d, exz.pe::%d, exz64.pe::%d, exz.ts::%d, exz64.ts::%d, exz.te::%d, exz64.te::%d\n", __func__, + // ql, thres, exz.err, exz64.err, exz.ps, exz64.ps, exz.pe, exz64.pe, exz.ts, exz64.ts, exz.te, exz64.te); + // fprintf(stderr, "[tstr] %.*s\n", (int32_t)(t_end+1-r_ts), t_string+r_ts); + // fprintf(stderr, "[qstr] %.*s\n", (int32_t)ql, q_string); + // } + ed_band_cal_global_infi_w(t_string+r_ts, t_end+1-r_ts, q_string, ql, thres, NULL, &exz); + ed_band_cal_global_infi_w_trace(t_string+r_ts, t_end+1-r_ts, q_string, ql, thres, NULL, &exz); + cigar_check(t_string+r_ts, q_string, &exz); + exz.err = INT32_MAX; + ed_band_cal_global_infi_w_trace(t_string+r_ts, t_end+1-r_ts, q_string, ql, thres, NULL, &exz); + cigar_check(t_string+r_ts, q_string, &exz); assert(exz.err <= (int32_t)error && exz.err >= 0); assert(exz.err == exz64.err && exz.ps == exz64.ps && exz.pe == exz64.pe && exz.ts == exz64.ts && exz.te == exz64.te); + ed_band_cal_global_256_w(t_string+r_ts, t_end+1-r_ts, q_string, ql, thres, &exz); + ed_band_cal_global_256_w_trace(t_string+r_ts, t_end+1-r_ts, q_string, ql, thres, &exz); + cigar_check(t_string+r_ts, q_string, &exz); + exz.err = INT32_MAX; + ed_band_cal_global_256_w_trace(t_string+r_ts, t_end+1-r_ts, q_string, ql, thres, &exz); + cigar_check(t_string+r_ts, q_string, &exz); + assert(exz.err <= (int32_t)error && exz.err >= 0); + assert(exz.err == exz64.err && exz.ps == exz64.ps && exz.pe == exz64.pe && exz.ts == exz64.ts && exz.te == exz64.te); - // ed_band_cal_semi_256_w(t_string, aln_l, q_string, ql, thres, &exz); - ed_band_cal_semi_infi_w(t_string, aln_l, q_string, ql, thres, NULL, &exz); + // if(!(exz.err == exz64.err && exz.ps == exz64.ps && exz.pe == exz64.pe && exz.ts == exz64.ts && exz.te == exz64.te)) { + // fprintf(stderr, "\n[M::%s::] ql::%ld, thres::%ld, exz.err::%d, exz64.err::%d, exz.ps::%d, exz64.ps::%d, exz.pe::%d, exz64.pe::%d, exz.ts::%d, exz64.ts::%d, exz.te::%d, exz64.te::%d\n", __func__, + // ql, thres, exz.err, exz64.err, exz.ps, exz64.ps, exz.pe, exz64.pe, exz.ts, exz64.ts, exz.te, exz64.te); + // } + ed_band_cal_semi_64_w(t_string, aln_l, q_string, ql, thres, &exz64); + ed_band_cal_semi_64_w_trace(t_string, aln_l, q_string, ql, thres, &exz64); + cigar_check(t_string, q_string, &(exz64)); + exz64.err = INT32_MAX; + ed_band_cal_semi_64_w_trace(t_string, aln_l, q_string, ql, thres, &exz64); + cigar_check(t_string, q_string, &(exz64)); + + // if(!cigar_check(t_string, q_string, &(exz64)) && ql <= 4) { + // fprintf(stderr, "\n[M::%s::] ql::%ld, thres::%ld, exz64.err::%d, exz64.ps::%d, exz64.pe::%d, exz64.ts::%d, exz64.te::%d\n", __func__, + // ql, thres, exz64.err, exz64.ps, exz64.pe, exz64.ts, exz64.te); + // fprintf(stderr, "[tstr] %.*s\n", (int32_t)aln_l, t_string); + // fprintf(stderr, "[qstr] %.*s\n", (int32_t)ql, q_string); + // } + + ed_band_cal_semi_infi_w(t_string, aln_l, q_string, ql, thres, NULL, &exz); + ed_band_cal_semi_infi_w_trace(t_string, aln_l, q_string, ql, thres, NULL, &exz); + cigar_check(t_string, q_string, &exz); + exz.err = INT32_MAX; + ed_band_cal_semi_infi_w_trace(t_string, aln_l, q_string, ql, thres, NULL, &exz); + cigar_check(t_string, q_string, &exz); + assert(exz.err <= (int32_t)error && exz.err >= 0); + assert(exz.err == exz64.err && exz.pe == exz64.pe && exz.ts == exz64.ts && exz.te == exz64.te); + + ed_band_cal_semi_256_w(t_string, aln_l, q_string, ql, thres, &exz); + ed_band_cal_semi_256_w_trace(t_string, aln_l, q_string, ql, thres, &exz); + cigar_check(t_string, q_string, &exz); + exz.err = INT32_MAX; + ed_band_cal_semi_256_w_trace(t_string, aln_l, q_string, ql, thres, &exz); + cigar_check(t_string, q_string, &exz); + assert(exz.err <= (int32_t)error && exz.err >= 0); + assert(exz.err == exz64.err && exz.pe == exz64.pe && exz.ts == exz64.ts && exz.te == exz64.te); + // if((!(exz.err == exz64.err && exz.ps == exz64.ps && exz.pe == exz64.pe && exz.ts == exz64.ts && exz.te == exz64.te)) || (exz.err != (int32_t)error)) { // fprintf(stderr, "\n[M::%s::semi] error::%u, ql::%ld, thres::%ld, exz.err::%d, exz64.err::%d, exz.ps::%d, exz64.ps::%d, exz.pe::%d, exz64.pe::%d, exz.ts::%d, exz64.ts::%d, exz.te::%d, exz64.te::%d\n", __func__, // error, ql, thres, exz.err, exz64.err, exz.ps, exz64.ps, exz.pe, exz64.pe, exz.ts, exz64.ts, exz.te, exz64.te); // fprintf(stderr, "[tstr] %.*s\n", (int32_t)aln_l, t_string); // fprintf(stderr, "[qstr] %.*s\n", (int32_t)ql, q_string); // } - assert(exz.err <= (int32_t)error && exz.err >= 0); - assert(exz.err == exz64.err && exz.ps == exz64.ps && exz.pe == exz64.pe && exz.ts == exz64.ts && exz.te == exz64.te); destroy_bit_extz_t(&exz); destroy_bit_extz_t(&exz64); // if(exz.err > (int32_t)error && ql == 1) { @@ -3713,6 +3807,7 @@ char *qstr, char *tstr, char *tstr1, Correct_dumy* dumy, uint32_t rev, uint32_t // assert(ed_band_cal_global(t_string+r_ts, t_end+1-r_ts, q_string, ql, thres) == // ed_band_cal_global_128bit(t_string+r_ts, t_end+1-r_ts, q_string, ql, thres)); + **/ ///this condition is always wrong ///in best case, r_ts = threshold, t_end = aln_l - thres - 1 @@ -11908,12 +12003,22 @@ void ul_lalign_hpc(overlap_region_alloc* ol, Candidates_list *cl, const ul_idx_t **/ } +// void ul_phase(overlap_region *oa, int64_t on, uint64_t *idx, uint64_t *buf) +// { +// int64_t i, oi; overlap_region *z; +// for (i = 0; i < on; i++) { +// oi = (uint32_t)idx[i]; z = &(oa[oi]); +// if(z->is_match != 1) continue; +// } +// } + void ul_lalign(overlap_region_alloc* ol, Candidates_list *cl, const ul_idx_t *uref, char *qstr, uint64_t ql, UC_Read* qu, UC_Read* tu, Correct_dumy* dumy, haplotype_evdience_alloc* hap, kvec_t_u64_warp* v_idx, double e_rate, int64_t wl, uint64_t is_base, void *km) { - uint64_t i, bs, k, ovl; Window_Pool w; double err; overlap_region t; overlap_region *z; + uint64_t i, bs, k, ovl/**, on**/; Window_Pool w; double err; + /**int64_t sc;**/ overlap_region t; overlap_region *z; ol->mapped_overlaps_length = 0; if(ol->length <= 0) return; @@ -11962,6 +12067,24 @@ void ul_lalign(overlap_region_alloc* ol, Candidates_list *cl, const ul_idx_t *ur } else { // fprintf(stderr, "-[M::%s] on::%lu\n", __func__, ol->length); if(ol->length <= 1) return; + // for (i = 0; (i < ol->length) && (ol->list[i].is_match == 1); i++); on = i; + // if(on <= 1) return; + // kv_resize(uint64_t, v_idx->a, (on<<1)); v_idx->a.n = on; + // for (i = 0; i < on; i++) { + // sc = ol->list[i].x_pos_e+1-ol->list[i].x_pos_s; + // sc -= ((int64_t)(ol->list[i].non_homopolymer_errors*ERROR_RATE)); + // if(sc < 0) sc = 0; + // v_idx->a.a[i] = sc; v_idx->a.a[i] <<= 32; v_idx->a.a[i] += i; + // } + // radix_sort_bc64(v_idx->a.a, v_idx->a.a+on); + // for (i = 0; i < on; i++) { + // ol->list[i].non_homopolymer_errors = 0; + // k = ((uint32_t)v_idx->a.a[i]); + // v_idx->a.a[k]<<=32; v_idx->a.a[k]>>=32; v_idx->a.a[k]|=i; + // } + // ul_phase(ol->list, on, v_idx->a.a, v_idx->a.a+on); + + for (i = 0; i < ol->length; i++) { z = &(ol->list[i]); ovl = z->x_pos_e+1-z->x_pos_s; z->is_match = 1; for (k = 0; k < z->w_list.n; k++) { diff --git a/Correct.h b/Correct.h index 76cd241..0f73994 100644 --- a/Correct.h +++ b/Correct.h @@ -17,6 +17,7 @@ #define MISMATCH 1 #define INSERTION 2 #define DELETION 3 +#define ERROR_RATE 1.25 #define WINDOW_MAX_SIZE (WINDOW + (int)(1.0 / HA_MIN_OV_DIFF) + 3) // TODO: why 1/max_ov_diff? diff --git a/Levenshtein_distance.h b/Levenshtein_distance.h index 635cd1a..4f9ad24 100644 --- a/Levenshtein_distance.h +++ b/Levenshtein_distance.h @@ -512,27 +512,49 @@ inline int32_t ed_band_cal_global(char *pstr, int32_t pn, char *tstr, int32_t tn return INT32_MAX; } +#define EAC_M 0 +#define MIS_M 1 +#define MOR_YP 2 +#define MOR_XT 3 -// typedef uint64_t w_sig; -// typedef struct {w_sig *a;} w128_t; -// #define bitw (6) -// #define bitwbit (64) -// #define bitz (63) +inline void push_trace(asg16_v *res, uint16_t c, uint32_t len) +{ + uint16_t p; c <<= 14; + while (len >= (0x3fff)) { + p = (c + (0x3fff)); kv_push(uint16_t, *res, p); len -= (0x3fff); + } + if(len) { + p = (c + len); kv_push(uint16_t, *res, p); + } +} + +inline uint32_t pop_trace(asg16_v *res, uint32_t i, uint16_t *c, uint32_t *len) +{ + (*c) = (res->a[i]>>14); (*len) = (res->a[i]&(0x3fff)); + for (i++; (i < res->n) && ((*c) == (res->a[i]>>14)); i++) { + (*len) += (res->a[i]&(0x3fff)); + } + return i; +} + + +typedef uint64_t w_sig; +#define bitw (6) +#define bitwbit (64) +#define bitz (63) // typedef uint32_t w_sig; -// typedef struct {w_sig a[3];} w128_t; // #define bitw (5) // #define bitwbit (32) // #define bitz (31) // typedef uint16_t w_sig; -// typedef struct {w_sig *a;} w128_t; // #define bitw (4) // #define bitwbit (16) // #define bitz (15) -typedef uint8_t w_sig; +// typedef uint8_t w_sig; +// #define bitw (3) +// #define bitwbit (8) +// #define bitz (7) typedef struct {w_sig *a;} w128_t; -#define bitw (3) -#define bitwbit (8) -#define bitz (7) typedef struct {size_t n, m; w_sig *a;} w64_trace_t; typedef struct { int32_t done_cigar, cigar_n, done_path, path_n; @@ -540,9 +562,59 @@ typedef struct { int32_t thre, err, nword, mword; uint32_t m, mm_thres; w_sig *a; w128_t Peq[5], mm, VP, VN, X, D0, HN, HP; - // asg16_v cigar; w64_trace_t path; + asg16_v cigar; w64_trace_t path; } bit_extz_t; +inline uint32_t cigar_check(char *pstr, char *tstr, bit_extz_t *ez) +{ + int32_t pi = ez->ps, ti = ez->ts, err = 0; uint32_t ci = 0, cl, k; uint16_t c; + while (ci < ez->cigar.n) { + ci = pop_trace(&(ez->cigar), ci, &c, &cl); + // fprintf(stderr, "# %u = %u\n", c, cl); + if(c == 0) { + for (k=0;(kerr) { + fprintf(stderr, "ERROR-err\n"); + return 0; + } + return 1; +} + +inline int32_t dbg_ext_err(int32_t i, int32_t thre, int32_t err, int32_t pe, w128_t *VP, w128_t *VN) +{ + int32_t poff = i-thre, tmp_e = INT32_MAX, k, bd, k_bd; + if((poff) + (thre<<1) >= (pe)) { + tmp_e = err; + for ((k) = 0; (poff) < (pe); (poff)++) { + bd = (k>>bitw); k_bd = (k&bitz); + (tmp_e) += ((VP->a[bd]>>k_bd)&((w_sig)1)); + (tmp_e) -= ((VN->a[bd]>>k_bd)&((w_sig)1)); + (k)++; + } + } + return tmp_e; +} + + inline void print_bits(w_sig *az, int64_t w, const char *cmd) { int64_t k, m, s = (sizeof(*az)<<3), sw = (w/s) + (!!(w%s)), ks; @@ -602,10 +674,95 @@ inline void init_bit_extz_t(bit_extz_t *ex, uint64_t thres) { } inline void destroy_bit_extz_t(bit_extz_t *ex) { - free((*ex).a); + free((*ex).a); free((*ex).path.a); free((*ex).cigar.a); } +inline void gen_trace(bit_extz_t *ez, int32_t ptrim, int32_t reverse)///ptrim = thre for global and extension; = 0 for semi +{ + ez->cigar.n = 0; + int32_t V, H, D, min, cur, tn = (ez->te+1-ez->ts), pn = tn + (ez->thre<<1), bd = (ez->thre<<1)+1; + int32_t bs = ez->path.n/tn, bbs = bs/5, poff = ez->pe, sft = bd - (pn - ez->pe - ptrim); + int32_t i = tn, low = bd-1, wi, wm, d = 0, pd = -1, pdn = 0; cur = ez->err; + w_sig *D0, *VP, *VN, *HP, *HN; + + // fprintf(stderr, "\n[M::%s::] bs::%d, bbs::%d, cur::%d, poff::%d, sft::%d\n", __func__, bs, bbs, cur, poff, sft); + + while (i > 0 && cur > 0) { + D0 = ez->path.a + ((i-1)*bs); VP = D0 + bbs; VN = VP + bbs; HP = VN + bbs; HN = HP + bbs; + + wi = (sft>>bitw); wm = sft&bitz; + D = cur - ((~(D0[wi]>>wm))&(1ULL)); d = 0; min = D; + H = V = INT32_MAX; + if(sft!=low) { + H = cur + ((HN[wi]>> wm)&(1ULL)) - ((HP[wi]>> wm)&(1ULL)); + if ((H+1) == cur && H <= min) {//prefer indels + min = H; d = 3; + } + } + if(sft!=0) { + wi = ((sft-1)>>bitw); wm = (sft-1)&bitz; + V = cur + ((VN[wi]>> wm)&(1ULL)) - ((VP[wi]>> wm)&(1ULL)); + if ((V+1) == cur && V <= min) {//prefer indels + min = V; d = 2; + } + } + // fprintf(stderr, "[M::%s::] cur::%d, D::%d, V::%d, H::%d, poff::%d, toff::%d, d::%d\n", __func__, cur, D, V, H, poff, i, d); + if(d == 0) { + if(D != cur) d = 1; + i--; poff--; + } else if(d == 2) {//more pstr + sft--; poff--; + } else if(d == 3) {///more tstr + i--; sft++; + } + + if(d == pd) { + pdn++; + } else { + if(pdn > 0) push_trace(&(ez->cigar), pd, pdn); + pd = d; pdn = 1; + } + cur = min; + } + + if (i > 0) { + d = 0; poff -= i; + if(d == pd) { + pdn += i; + } else { + if(pdn > 0) push_trace(&(ez->cigar), pd, pdn); + pd = d; pdn = i; + } + } + + poff++; + // fprintf(stderr, "[M::%s::] poff::%d, ez->ps::%d\n", __func__, poff, ez->ps); + if(ez->ps < 0 || ez->ps >= ez->pl) {//ps is unavailable + ez->ps = poff; + } else if(poff > ez->ps){ + d = 2; i = poff - ez->ps; + if(d == pd) { + pdn += i; + } else { + if(pdn > 0) push_trace(&(ez->cigar), pd, pdn); + pd = d; pdn = i; + } + } + if(pdn > 0) push_trace(&(ez->cigar), pd, pdn); + + if(reverse) { + uint16_t t; pdn = ez->cigar.n>>1; + for (i = 0; i < pdn; i++) { + t = ez->cigar.a[i]; + ez->cigar.a[i] = ez->cigar.a[ez->cigar.n-i-1]; + ez->cigar.a[ez->cigar.n-i-1] = t; + } + } +} + + + #define init_base_ed(ez, thre, pn, tn) {\ (ez).thre = (thre), (ez).err = INT32_MAX, (ez).pl = pn, (ez).tl = tn;\ (ez).done_cigar = (ez).done_path = (ez).cigar_n = (ez).path_n = 0;\ @@ -618,6 +775,8 @@ inline void destroy_bit_extz_t(bit_extz_t *ex) { /***********************2 words***********************/ #define w_128_clear(x) ((x).a[0]=(x).a[1]=0) +#define w_128_word (2) + #define w_128_self_not(x) ((x).a[0]=~(x).a[0], \ (x).a[1]=~(x).a[1]) @@ -655,6 +814,8 @@ inline void destroy_bit_extz_t(bit_extz_t *ex) { /***********************3 words***********************/ #define w_192_clear(x) ((x).a[0]=(x).a[1]=(x).a[2]=0) +#define w_192_word (3) + #define w_192_self_not(x) ((x).a[0]=~(x).a[0],\ (x).a[1]=~(x).a[1],\ (x).a[2]=~(x).a[2]) @@ -702,6 +863,8 @@ inline void destroy_bit_extz_t(bit_extz_t *ex) { /***********************4 words***********************/ #define w_256_clear(x) ((x).a[0]=(x).a[1]=(x).a[2]=(x).a[3]=0) +#define w_256_word (4) + #define w_256_self_not(x) ((x).a[0]=~(x).a[0],\ (x).a[1]=~(x).a[1],\ (x).a[2]=~(x).a[2],\ @@ -904,8 +1067,7 @@ inline void ed_band_cal_global_##sf##_w(char *pstr, int32_t pn, char *tstr, int3 {\ init_base_ed(*ez, thre, pn, tn); ez->ps = ez->ts = 0;\ if((pn > tn + thre) || (tn > pn + thre)) return;\ - /**if((pn < thre + 1) || (tn < thre + 1)) return;**/\ - int32_t i, err, tn0 = tn - 1, cut = thre+(thre<<1), bd, i_bd; w_sig c;\ + int32_t i, err, tn0 = tn - 1, cut = thre+(thre<<1), bd, i_bd, Peq_i; w_sig c, Peq_m;\ w_##sf##_clear(ez->Peq[0]); w_##sf##_clear(ez->Peq[1]); w_##sf##_clear(ez->Peq[2]); w_##sf##_clear(ez->Peq[3]); w_##sf##_clear(ez->Peq[4]);\ \ bd = thre+1; bd = ((bd<=pn)?bd:pn); i_bd = thre;\ @@ -922,7 +1084,7 @@ inline void ed_band_cal_global_##sf##_w(char *pstr, int32_t pn, char *tstr, int3 /**should make Peq[4] = 0 if N is always an error**/\ i = 0; \ /**for the incoming char/last char**/\ - w_##sf##_clear(ez->mm); w_bit(ez->mm, (thre<<1)); /**mm = ((Word)1 << (thre<<1));**/\ + Peq_i = (((thre<<1))>>bitw); Peq_m = (((w_sig)1)<<(((thre<<1))&bitz));\ while (i < tn0) {\ ed_core(_##sf##_, ez->Peq, ez->VP, ez->VN, ez->X, ez->D0, ez->HN, ez->HP, (uint8_t)tstr[i], c);\ if (!(ez->D0.a[0]&(1ULL))) {\ @@ -931,12 +1093,10 @@ inline void ed_band_cal_global_##sf##_w(char *pstr, int32_t pn, char *tstr, int3 /** Peq[0] >>= 1; Peq[1] >>= 1; Peq[2] >>= 1; Peq[3] >>= 1;**/\ w_##sf##_self_rsft_1(ez->Peq[0]); w_##sf##_self_rsft_1(ez->Peq[1]);\ w_##sf##_self_rsft_1(ez->Peq[2]); w_##sf##_self_rsft_1(ez->Peq[3]);\ - ++i; ++i_bd;\ - if(i_bd < pn) {\ - c = seq_nt4_table[(uint8_t)pstr[i_bd]];\ - /**if(c < 4) Peq[c] |= mm;**/\ - if(c < 4) w_##sf##_self_or(ez->Peq[c], ez->mm);\ - }\ + \ + ++i; ++i_bd; c = 4;\ + if(i_bd < pn) c = seq_nt4_table[(uint8_t)pstr[i_bd]];\ + if(c < 4) ez->Peq[c].a[Peq_i]|=Peq_m;\ }\ ed_core(_##sf##_, ez->Peq, ez->VP, ez->VN, ez->X, ez->D0, ez->HN, ez->HP, (uint8_t)tstr[i], c);\ if (!(ez->D0.a[0]&(1ULL))) {\ @@ -944,11 +1104,10 @@ inline void ed_band_cal_global_##sf##_w(char *pstr, int32_t pn, char *tstr, int3 }\ \ int32_t site = tn - 1 - thre;/**up bound**/\ - for (cut = pn - 1; site < cut; site++) {\ - /**err += ((VP >> i)&(1ULL));**/\ - err += ez->VP.a[0]&(1ULL); w_##sf##_self_rsft_1(ez->VP); \ - /**err -= ((VN >> i)&(1ULL));**/\ - err -= ez->VN.a[0]&(1ULL); w_##sf##_self_rsft_1(ez->VN); \ + for (cut = pn - 1, i = 0; site < cut; site++, i++) {\ + bd = (i>>bitw); i_bd = (i&bitz);\ + err += (((*ez).VP.a[bd]>>i_bd)&((w_sig)1));\ + err -= (((*ez).VN.a[bd]>>i_bd)&((w_sig)1));\ }\ \ if (site == cut && err <= thre) {\ @@ -960,20 +1119,24 @@ inline void ed_band_cal_global_##sf##_w(char *pstr, int32_t pn, char *tstr, int3 inline void ed_band_cal_semi_##sf##_w(char *pstr, int32_t pn, char *tstr, int32_t tn, int32_t thre, bit_extz_t *ez)\ {\ init_base_ed(*ez, thre, pn, tn); ez->ps = ez->pe = -1; ez->ts = 0; ez->te = tn-1;\ - w_##sf##_clear(ez->VP); w_##sf##_clear(ez->VN);\ - w_##sf##_clear(ez->Peq[0]); w_##sf##_clear(ez->Peq[1]); w_##sf##_clear(ez->Peq[2]); w_##sf##_clear(ez->Peq[3]); w_##sf##_clear(ez->Peq[4]);\ - int32_t bd, i, err = 0, i_bd, last_high = (thre<<1), tn0 = tn - 1;\ - int32_t cut = thre+last_high; w_sig c;\ + int32_t bd, i, err = 0, i_bd, cut = thre+(thre<<1), tn0 = tn - 1, Peq_i; w_sig c, Peq_m;\ + \ + w_##sf##_clear(ez->Peq[0]);\ + w_##sf##_clear(ez->Peq[1]);\ + w_##sf##_clear(ez->Peq[2]);\ + w_##sf##_clear(ez->Peq[3]);\ + w_##sf##_clear(ez->Peq[4]);\ + w_##sf##_clear(ez->VP);\ + w_##sf##_clear(ez->VN);\ \ bd = (thre<<1)+1; bd = ((bd<=pn)?bd:pn); i_bd = 0;\ ed_init_core(i, bd, i_bd, pstr, ez->Peq);\ bd = (thre<<1)+1, i_bd = (thre<<1);\ - /** Peq[4] = 0;**/\ - w_##sf##_clear(ez->Peq[4]);\ - /**mm = ((Word)1 << (thre<<1));///for the incoming char/last char**/\ - w_##sf##_clear(ez->mm); w_bit(ez->mm, (thre<<1));\ \ + w_##sf##_clear(ez->Peq[4]);\ i = 0;\ + /**for the incoming char/last char; mm = ((Word)1 << (thre<<1))**/\ + Peq_i = (((thre<<1))>>bitw); Peq_m = (((w_sig)1)<<(((thre<<1))&bitz));\ while (i < tn0) {\ ed_core(_##sf##_, ez->Peq, ez->VP, ez->VN, ez->X, ez->D0, ez->HN, ez->HP, (uint8_t)tstr[i], c);\ if (!(ez->D0.a[0]&(1ULL))) {\ @@ -982,10 +1145,9 @@ inline void ed_band_cal_semi_##sf##_w(char *pstr, int32_t pn, char *tstr, int32_ /** Peq[0] >>= 1; Peq[1] >>= 1; Peq[2] >>= 1; Peq[3] >>= 1;**/\ w_##sf##_self_rsft_1(ez->Peq[0]); w_##sf##_self_rsft_1(ez->Peq[1]);\ w_##sf##_self_rsft_1(ez->Peq[2]); w_##sf##_self_rsft_1(ez->Peq[3]);\ - ++i; ++i_bd; \ - /** Peq[seq_nt4_table[(uint8_t)pstr[i_bd]]] |= mm; Peq[4] = 0;**/\ - c = seq_nt4_table[(uint8_t)pstr[i_bd]]; \ - if(c < 4) w_##sf##_self_or(ez->Peq[c], ez->mm);\ + ++i; ++i_bd; c = 4;\ + if(i_bd < pn) c = seq_nt4_table[(uint8_t)pstr[i_bd]];\ + if(c < 4) ez->Peq[c].a[Peq_i]|=Peq_m;\ }\ ed_core(_##sf##_, ez->Peq, ez->VP, ez->VN, ez->X, ez->D0, ez->HN, ez->HP, (uint8_t)tstr[i], c);\ if (!(ez->D0.a[0]&(1ULL))) {\ @@ -998,13 +1160,12 @@ inline void ed_band_cal_semi_##sf##_w(char *pstr, int32_t pn, char *tstr, int32_ if ((err <= thre) && (err <= ez->err)) {\ ez->err = err; ez->pe = site;\ }\ - i = 0;\ \ + i = 0;\ while (i < ai) {\ - /** err += ((VP >> i)&(1ULL)); **/\ - err += ez->VP.a[0]&(1ULL); w_##sf##_self_rsft_1(ez->VP); \ - /** err -= ((VN >> i)&(1ULL)); **/\ - err -= ez->VN.a[0]&(1ULL); w_##sf##_self_rsft_1(ez->VN); \ + bd = (i>>bitw); i_bd = (i&bitz);\ + err += (((*ez).VP.a[bd]>>i_bd)&((w_sig)1));\ + err -= (((*ez).VN.a[bd]>>i_bd)&((w_sig)1));\ ++i;\ if ((err <= thre) && (err <= ez->err)) {\ ez->err = err; ez->pe = site + i;\ @@ -1015,18 +1176,24 @@ inline void ed_band_cal_semi_##sf##_w(char *pstr, int32_t pn, char *tstr, int32_ if((uge <= thre) && (uge == ez->err)) ez->pe = site + thre;\ }\ /**require:: (pn >= tn - thre && pn <= tn + thre)**/\ -inline void ed_band_cal_extension_##sf##_w(char *pstr, int32_t pn, char *tstr, int32_t tn, int32_t thre, bit_extz_t *ez)\ +inline void ed_band_cal_extension_##sf##_0_w(char *pstr, int32_t pn, char *tstr, int32_t tn, int32_t thre, bit_extz_t *ez)\ {\ + init_base_ed(*ez, thre, pn, tn); ez->ps = ez->ts = 0; ez->pe = ez->te = -1;\ if(pn > tn + thre) pn = tn + thre;\ else if(tn > pn + thre) tn = pn + thre;\ - init_base_ed(*ez, thre, pn, tn); ez->ps = ez->ts = 0; ez->pe = ez->te = -1;\ - int32_t i, err, tn0 = tn - 1, cut = thre+(thre<<1), bd, i_bd; w_sig c;\ - int32_t poff, pe = pn-1, tmp_e, k;\ - w_##sf##_clear(ez->Peq[0]); w_##sf##_clear(ez->Peq[1]); w_##sf##_clear(ez->Peq[2]); w_##sf##_clear(ez->Peq[3]); w_##sf##_clear(ez->Peq[4]);\ + int32_t i, err, tn0 = tn - 1, cut = thre+(thre<<1), bd, i_bd, k_bd, Peq_i; w_sig c, Peq_m;\ + int32_t poff, pe = pn-1, tmp_e = INT32_MAX, k;\ + \ + w_##sf##_clear(ez->Peq[0]);\ + w_##sf##_clear(ez->Peq[1]);\ + w_##sf##_clear(ez->Peq[2]);\ + w_##sf##_clear(ez->Peq[3]);\ + w_##sf##_clear(ez->Peq[4]);\ \ bd = thre+1; bd = ((bd<=pn)?bd:pn); i_bd = thre;\ ed_init_core(i, bd, i_bd, pstr, ez->Peq);\ bd = thre+1, i_bd = thre;\ + \ w_##sf##_clear(ez->Peq[4]);\ err = thre;\ w_##sf##_set_bit_lsub(ez->VN, thre); /**VN = (((Word)1)<<(thre))-1; **/\ @@ -1036,46 +1203,56 @@ inline void ed_band_cal_extension_##sf##_w(char *pstr, int32_t pn, char *tstr, i /** print_bits(ez->VP.a, (thre<<1)+1, "-VP");**/\ /**should make Peq[4] = 0 if N is always an error**/\ i = 0; \ - /**for the incoming char/last char**/\ - w_##sf##_clear(ez->mm); w_bit(ez->mm, (thre<<1)); /**mm = ((Word)1 << (thre<<1));**/\ + /**for the incoming char/last char; mm = ((Word)1 << (thre<<1))**/\ + Peq_i = (((thre<<1))>>bitw); Peq_m = (((w_sig)1)<<(((thre<<1))&bitz));\ while (i < tn0) {\ ed_core(_##sf##_, ez->Peq, ez->VP, ez->VN, ez->X, ez->D0, ez->HN, ez->HP, (uint8_t)tstr[i], c);\ if (!(ez->D0.a[0]&(1ULL))) {\ ++err; if (err>cut) return;\ }\ \ - {\ - poff = i-thre; tmp_e = err; /**poff:[i-thre, i+thre]**/\ - if((poff) + (((*ez).thre)<<1) >= (pe)) { \ - for ((k) = 0; (poff) < (pe); (poff)++) {\ - (tmp_e) += w_get_bit((*ez).VP, (k));\ - (tmp_e) -= w_get_bit((*ez).VN, (k));\ - (k)++;}\ - if((tmp_e) <= (*ez).thre && (tmp_e) < (*ez).err) {\ - (*ez).err = tmp_e; (*ez).pe = poff; (*ez).te = i;}\ - }\ - }\ + {\ + poff = i-thre; k = i+thre-pe;/**poff:[i-thre, i+thre]**/\ + if(k >= 0) {\ + if(k == 0) {\ + tmp_e = err;\ + for ((k) = 0; (poff) < (pe); (poff)++, (k)++) {\ + bd = (k>>bitw); k_bd = (k&bitz);\ + (tmp_e) += (((*ez).VP.a[bd]>>k_bd)&((w_sig)1));\ + (tmp_e) -= (((*ez).VN.a[bd]>>k_bd)&((w_sig)1));\ + }\ + } else {\ + k = (thre<<1) - k;\ + if(k >= 0) {\ + bd = (k>>bitw); k_bd = (k&bitz);\ + (tmp_e) += (((*ez).HP.a[bd]>>k_bd)&((w_sig)1));\ + (tmp_e) -= (((*ez).HN.a[bd]>>k_bd)&((w_sig)1));\ + }\ + }\ + if((tmp_e) <= (*ez).thre && (tmp_e) < (*ez).err) {\ + (*ez).err = tmp_e; (*ez).pe = pe; (*ez).te = i;\ + }\ + }\ + /**if(dbg_ext_err(i, thre, err, pe, &((*ez).VP), &((*ez).VN))!=tmp_e)**/\ + }\ \ /** Peq[0] >>= 1; Peq[1] >>= 1; Peq[2] >>= 1; Peq[3] >>= 1;**/\ w_##sf##_self_rsft_1(ez->Peq[0]); w_##sf##_self_rsft_1(ez->Peq[1]);\ w_##sf##_self_rsft_1(ez->Peq[2]); w_##sf##_self_rsft_1(ez->Peq[3]);\ - ++i; ++i_bd;\ - if(i_bd < pn) {\ - c = seq_nt4_table[(uint8_t)pstr[i_bd]]; \ - if(c < 4) w_##sf##_self_or(ez->Peq[c], ez->mm);\ - }\ + ++i; ++i_bd; c = 4;\ + if(i_bd < pn) c = seq_nt4_table[(uint8_t)pstr[i_bd]];\ + if(c < 4) ez->Peq[c].a[Peq_i]|=Peq_m;\ }\ ed_core(_##sf##_, ez->Peq, ez->VP, ez->VN, ez->X, ez->D0, ez->HN, ez->HP, (uint8_t)tstr[i], c);\ if (!(ez->D0.a[0]&(1ULL))) {\ ++err; if (err>cut) return;\ }\ - /**i = tn - 1**/\ + \ int32_t site = tn - 1 - thre;/**up bound; site:[tn - 1 - thre, tn - 1 + thre]**/\ - for (cut = pn - 1; site < cut; ) {\ - /** err += ((VP >> i)&(1ULL)); **/\ - err += ez->VP.a[0]&(1ULL); w_##sf##_self_rsft_1(ez->VP); \ - /** err -= ((VN >> i)&(1ULL));**/\ - err -= ez->VN.a[0]&(1ULL); w_##sf##_self_rsft_1(ez->VN); \ + for (cut = pn - 1, i = 0; site < cut; i++) {\ + bd = (i>>bitw); i_bd = (i&bitz);\ + err += (((*ez).VP.a[bd]>>i_bd)&((w_sig)1));\ + err -= (((*ez).VN.a[bd]>>i_bd)&((w_sig)1));\ site++;\ if(err <= thre && err < ez->err) {\ ez->err = err; ez->pe = site; ez->te = tn-1;\ @@ -1085,7 +1262,491 @@ inline void ed_band_cal_extension_##sf##_w(char *pstr, int32_t pn, char *tstr, i ez->err = err; ez->pe = site; ez->te = tn-1;\ }\ return;\ -} +}\ +/**require:: (pn >= tn - thre && pn <= tn + thre)**/\ +inline void ed_band_cal_extension_##sf##_1_w(char *pstr, int32_t pn, char *tstr, int32_t tn, int32_t thre, bit_extz_t *ez)\ +{\ + init_base_ed(*ez, thre, pn, tn); ez->ps = ez->ts = INT32_MAX; ez->pe = pn-1; ez->te = tn-1;\ + if(pn > tn + thre) pn = tn + thre;\ + else if(tn > pn + thre) tn = pn + thre;\ + int32_t i, err, tn0 = tn - 1, cut = thre+(thre<<1), bd, i_bd, k_bd, Peq_i; w_sig c, Peq_m;\ + int32_t poff, pe = pn-1, tmp_e = INT32_MAX, k, pidx = ez->pe, tidx = ez->te;\ + \ + w_##sf##_clear(ez->Peq[0]);\ + w_##sf##_clear(ez->Peq[1]);\ + w_##sf##_clear(ez->Peq[2]);\ + w_##sf##_clear(ez->Peq[3]);\ + w_##sf##_clear(ez->Peq[4]);\ + \ + bd = thre+1; bd = ((bd<=pn)?bd:pn); i_bd = thre;\ + for (i = 0; i < bd; i++, i_bd++) {\ + w_bit((ez->Peq)[seq_nt4_table[(uint8_t)(pstr)[pidx-i]]], (i_bd));\ + }\ + bd = thre+1, i_bd = thre;\ + \ + w_##sf##_clear(ez->Peq[4]);\ + err = thre;\ + w_##sf##_set_bit_lsub(ez->VN, thre); /**VN = (((Word)1)<<(thre))-1; **/\ + w_##sf##_set_bit_lsub(ez->VP, (thre<<1)+1); /**VP = (((Word)1)<<((thre<<1)+1))-1;**/\ + w_##sf##_self_xor(ez->VP, ez->VN); /**VP ^= VN;**/\ + \ + /** print_bits(ez->VP.a, (thre<<1)+1, "-VP");**/\ + /**should make Peq[4] = 0 if N is always an error**/\ + i = 0; \ + /**for the incoming char/last char; mm = ((Word)1 << (thre<<1))**/\ + Peq_i = (((thre<<1))>>bitw); Peq_m = (((w_sig)1)<<(((thre<<1))&bitz));\ + while (i < tn0) {\ + ed_core(_##sf##_, ez->Peq, ez->VP, ez->VN, ez->X, ez->D0, ez->HN, ez->HP, (uint8_t)tstr[tidx-i], c);\ + if (!(ez->D0.a[0]&(1ULL))) {\ + ++err; if (err>cut) return;\ + }\ + \ + {\ + poff = i-thre; k = i+thre-pe;/**poff:[i-thre, i+thre]**/\ + if(k >= 0) {\ + if(k == 0) {\ + tmp_e = err;\ + for ((k) = 0; (poff) < (pe); (poff)++, (k)++) {\ + bd = (k>>bitw); k_bd = (k&bitz);\ + (tmp_e) += (((*ez).VP.a[bd]>>k_bd)&((w_sig)1));\ + (tmp_e) -= (((*ez).VN.a[bd]>>k_bd)&((w_sig)1));\ + }\ + } else {\ + k = (thre<<1) - k;\ + if(k >= 0) {\ + bd = (k>>bitw); k_bd = (k&bitz);\ + (tmp_e) += (((*ez).HP.a[bd]>>k_bd)&((w_sig)1));\ + (tmp_e) -= (((*ez).HN.a[bd]>>k_bd)&((w_sig)1));\ + }\ + }\ + if((tmp_e) <= (*ez).thre && (tmp_e) < (*ez).err) {\ + (*ez).err = tmp_e; (*ez).ps = pidx - pe; (*ez).ts = tidx-i;\ + }\ + }\ + /**if(dbg_ext_err(i, thre, err, pe, &((*ez).VP), &((*ez).VN))!=tmp_e)**/\ + }\ + \ + /** Peq[0] >>= 1; Peq[1] >>= 1; Peq[2] >>= 1; Peq[3] >>= 1;**/\ + w_##sf##_self_rsft_1(ez->Peq[0]); w_##sf##_self_rsft_1(ez->Peq[1]);\ + w_##sf##_self_rsft_1(ez->Peq[2]); w_##sf##_self_rsft_1(ez->Peq[3]);\ + ++i; ++i_bd; c = 4;\ + if(i_bd < pn) c = seq_nt4_table[(uint8_t)pstr[pidx-i_bd]];\ + if(c < 4) ez->Peq[c].a[Peq_i]|=Peq_m;\ + }\ + ed_core(_##sf##_, ez->Peq, ez->VP, ez->VN, ez->X, ez->D0, ez->HN, ez->HP, (uint8_t)tstr[tidx-i], c);\ + if (!(ez->D0.a[0]&(1ULL))) {\ + ++err; if (err>cut) return;\ + }\ + \ + int32_t site = tn - 1 - thre;/**up bound; site:[tn - 1 - thre, tn - 1 + thre]**/\ + for (cut = pn - 1, i = 0; site < cut; i++) {\ + bd = (i>>bitw); i_bd = (i&bitz);\ + err += (((*ez).VP.a[bd]>>i_bd)&((w_sig)1));\ + err -= (((*ez).VN.a[bd]>>i_bd)&((w_sig)1));\ + site++;\ + if(err <= thre && err < ez->err) {\ + ez->err = err; ez->ps = pidx-site; ez->ts = tidx+1-tn;\ + }\ + }\ + if(err <= thre && err < ez->err) {\ + ez->err = err; ez->ps = pidx-site; ez->ts = tidx+1-tn;\ + }\ + return;\ +}\ +inline void ed_band_cal_global_##sf##_w_trace(char *pstr, int32_t pn, char *tstr, int32_t tn, int32_t thre, bit_extz_t *ez)\ +{\ + ez->cigar.n = 0; ez->nword = w_##sf##_word;\ + if(ez->err > thre) {\ + init_base_ed(*ez, thre, pn, tn); ez->ps = ez->ts = 0;\ + } else if(ez->err == 0) {\ + push_trace(&(ez->cigar), 0, ez->te+1-ez->ts); return;\ + }\ + if((pn > tn + thre) || (tn > pn + thre)) return;\ + int32_t i, err, tn0 = tn - 1, cut = thre+(thre<<1), bd, i_bd, ws, Peq_i; w_sig c, Peq_m;\ + w_##sf##_clear(ez->Peq[0]); w_##sf##_clear(ez->Peq[1]); w_##sf##_clear(ez->Peq[2]); w_##sf##_clear(ez->Peq[3]); w_##sf##_clear(ez->Peq[4]);\ + \ + bd = thre+1; bd = ((bd<=pn)?bd:pn); i_bd = thre;\ + ed_init_core(i, bd, i_bd, pstr, ez->Peq);\ + bd = thre+1, i_bd = thre;\ + \ + w_##sf##_clear(ez->Peq[4]);\ + err = thre;\ + w_##sf##_set_bit_lsub(ez->VN, thre); /**VN = (((Word)1)<<(thre))-1;**/\ + w_##sf##_set_bit_lsub(ez->VP, (thre<<1)+1); /**VP = (((Word)1)<<((thre<<1)+1))-1;**/\ + w_##sf##_self_xor(ez->VP, ez->VN); /**VP ^= VN;**/\ + \ + /** print_bits(VP.a, (thre<<1)+1, "-VP");**/\ + /**should make Peq[4] = 0 if N is always an error**/\ + i = 0; ws = sizeof(*(ez->a))*(ez->nword);\ + ez->path.n=(ez->nword*tn*5);\ + kv_resize(w_sig, ez->path, ez->path.n); ez->path.n=0;\ + \ + /**for the incoming char/last char**/\ + Peq_i = (((thre<<1))>>bitw); Peq_m = (((w_sig)1)<<(((thre<<1))&bitz));\ + while (i < tn0) {\ + ed_core(_##sf##_, ez->Peq, ez->VP, ez->VN, ez->X, ez->D0, ez->HN, ez->HP, (uint8_t)tstr[i], c);\ + if (!(ez->D0.a[0]&(1ULL))) {\ + ++err; if (err>cut) return;\ + }\ + /** Peq[0] >>= 1; Peq[1] >>= 1; Peq[2] >>= 1; Peq[3] >>= 1;**/\ + w_##sf##_self_rsft_1(ez->Peq[0]); w_##sf##_self_rsft_1(ez->Peq[1]);\ + w_##sf##_self_rsft_1(ez->Peq[2]); w_##sf##_self_rsft_1(ez->Peq[3]);\ + \ + ++i; ++i_bd; c = 4;\ + if(i_bd < pn) c = seq_nt4_table[(uint8_t)pstr[i_bd]];\ + if(c < 4) ez->Peq[c].a[Peq_i]|=Peq_m;\ + \ + memcpy(ez->path.a+ez->path.n, ez->D0.a, ws); ez->path.n += ez->nword;\ + memcpy(ez->path.a+ez->path.n, ez->VP.a, ws); ez->path.n += ez->nword;\ + memcpy(ez->path.a+ez->path.n, ez->VN.a, ws); ez->path.n += ez->nword;\ + memcpy(ez->path.a+ez->path.n, ez->HP.a, ws); ez->path.n += ez->nword;\ + memcpy(ez->path.a+ez->path.n, ez->HN.a, ws); ez->path.n += ez->nword;\ + }\ + ed_core(_##sf##_, ez->Peq, ez->VP, ez->VN, ez->X, ez->D0, ez->HN, ez->HP, (uint8_t)tstr[i], c);\ + if (!(ez->D0.a[0]&(1ULL))) {\ + ++err; if (err>cut) return;\ + }\ + memcpy(ez->path.a+ez->path.n, ez->D0.a, ws); ez->path.n += ez->nword;\ + memcpy(ez->path.a+ez->path.n, ez->VP.a, ws); ez->path.n += ez->nword;\ + memcpy(ez->path.a+ez->path.n, ez->VN.a, ws); ez->path.n += ez->nword;\ + memcpy(ez->path.a+ez->path.n, ez->HP.a, ws); ez->path.n += ez->nword;\ + memcpy(ez->path.a+ez->path.n, ez->HN.a, ws); ez->path.n += ez->nword;\ + \ + int32_t site = tn - 1 - thre;/**up bound**/\ + if(ez->err > thre) {\ + for (cut = pn - 1, i = 0; site < cut; site++, i++) {\ + bd = (i>>bitw); i_bd = (i&bitz);\ + err += (((*ez).VP.a[bd]>>i_bd)&((w_sig)1));\ + err -= (((*ez).VN.a[bd]>>i_bd)&((w_sig)1));\ + }\ + \ + if (site == cut && err <= thre) {\ + ez->err = err; \ + ez->pe = pn-1; ez->te = tn-1;\ + }\ + }\ + gen_trace(ez, thre, 1);\ + return;\ +}\ +inline void ed_band_cal_semi_##sf##_w_trace(char *pstr, int32_t pn, char *tstr, int32_t tn, int32_t thre, bit_extz_t *ez)\ +{\ + ez->cigar.n = 0; ez->nword = w_##sf##_word;\ + if(ez->err > thre) {\ + init_base_ed(*ez, thre, pn, tn); ez->ps = ez->pe = -1; ez->ts = 0; ez->te = tn-1;\ + } else if(ez->err == 0) {\ + push_trace(&(ez->cigar), 0, ez->te+1-ez->ts);\ + ez->ps = ez->pe - (ez->te-ez->ts);\ + return;\ + }\ + int32_t bd, i, err = 0, i_bd, cut = thre+(thre<<1), tn0 = tn - 1, Peq_i, ws; w_sig c, Peq_m;\ + \ + w_##sf##_clear(ez->Peq[0]);\ + w_##sf##_clear(ez->Peq[1]);\ + w_##sf##_clear(ez->Peq[2]);\ + w_##sf##_clear(ez->Peq[3]);\ + w_##sf##_clear(ez->Peq[4]);\ + w_##sf##_clear(ez->VP);\ + w_##sf##_clear(ez->VN);\ + \ + bd = (thre<<1)+1; bd = ((bd<=pn)?bd:pn); i_bd = 0;\ + ed_init_core(i, bd, i_bd, pstr, ez->Peq);\ + bd = (thre<<1)+1, i_bd = (thre<<1);\ + \ + w_##sf##_clear(ez->Peq[4]);\ + i = 0; ws = sizeof(*(ez->a))*(ez->nword);\ + ez->path.n=(ez->nword*tn*5);\ + kv_resize(w_sig, ez->path, ez->path.n); ez->path.n=0;\ + /**for the incoming char/last char; mm = ((Word)1 << (thre<<1))**/\ + Peq_i = (((thre<<1))>>bitw); Peq_m = (((w_sig)1)<<(((thre<<1))&bitz));\ + while (i < tn0) {\ + ed_core(_##sf##_, ez->Peq, ez->VP, ez->VN, ez->X, ez->D0, ez->HN, ez->HP, (uint8_t)tstr[i], c);\ + if (!(ez->D0.a[0]&(1ULL))) {\ + ++err; if (err>cut) return;\ + }\ + /** Peq[0] >>= 1; Peq[1] >>= 1; Peq[2] >>= 1; Peq[3] >>= 1;**/\ + w_##sf##_self_rsft_1(ez->Peq[0]); w_##sf##_self_rsft_1(ez->Peq[1]);\ + w_##sf##_self_rsft_1(ez->Peq[2]); w_##sf##_self_rsft_1(ez->Peq[3]);\ + ++i; ++i_bd; c = 4;\ + if(i_bd < pn) c = seq_nt4_table[(uint8_t)pstr[i_bd]];\ + if(c < 4) ez->Peq[c].a[Peq_i]|=Peq_m;\ + \ + memcpy(ez->path.a+ez->path.n, ez->D0.a, ws); ez->path.n += ez->nword;\ + memcpy(ez->path.a+ez->path.n, ez->VP.a, ws); ez->path.n += ez->nword;\ + memcpy(ez->path.a+ez->path.n, ez->VN.a, ws); ez->path.n += ez->nword;\ + memcpy(ez->path.a+ez->path.n, ez->HP.a, ws); ez->path.n += ez->nword;\ + memcpy(ez->path.a+ez->path.n, ez->HN.a, ws); ez->path.n += ez->nword;\ + }\ + ed_core(_##sf##_, ez->Peq, ez->VP, ez->VN, ez->X, ez->D0, ez->HN, ez->HP, (uint8_t)tstr[i], c);\ + if (!(ez->D0.a[0]&(1ULL))) {\ + ++err; if (err>cut) return;\ + }\ + memcpy(ez->path.a+ez->path.n, ez->D0.a, ws); ez->path.n += ez->nword;\ + memcpy(ez->path.a+ez->path.n, ez->VP.a, ws); ez->path.n += ez->nword;\ + memcpy(ez->path.a+ez->path.n, ez->VN.a, ws); ez->path.n += ez->nword;\ + memcpy(ez->path.a+ez->path.n, ez->HP.a, ws); ez->path.n += ez->nword;\ + memcpy(ez->path.a+ez->path.n, ez->HN.a, ws); ez->path.n += ez->nword;\ + \ + int32_t site = tn - 1;/**up bound**/\ + /**in most cases, ai = (thre<<1)**/\ + int32_t ai = pn - tn, uge = INT32_MAX;\ + if(ez->err > thre) {\ + if ((err <= thre) && (err <= ez->err)) {\ + ez->err = err; ez->pe = site;\ + }\ + \ + i = 0;\ + while (i < ai) {\ + bd = (i>>bitw); i_bd = (i&bitz);\ + err += (((*ez).VP.a[bd]>>i_bd)&((w_sig)1));\ + err -= (((*ez).VN.a[bd]>>i_bd)&((w_sig)1));\ + ++i;\ + if ((err <= thre) && (err <= ez->err)) {\ + ez->err = err; ez->pe = site + i;\ + }\ + if(i == thre) uge = err;\ + }\ + \ + if((uge <= thre) && (uge == ez->err)) ez->pe = site + thre;\ + }\ + gen_trace(ez, 0, 1);\ +}\ +inline void ed_band_cal_extension_##sf##_0_w_trace(char *pstr, int32_t pn, char *tstr, int32_t tn, int32_t thre, bit_extz_t *ez)\ +{\ + int32_t done = 1;\ + ez->cigar.n = 0; ez->nword = w_##sf##_word;\ + if(ez->err > thre) {\ + init_base_ed(*ez, thre, pn, tn); ez->ps = ez->ts = 0; ez->pe = ez->te = -1;\ + if(pn > tn + thre) pn = tn + thre;\ + else if(tn > pn + thre) tn = pn + thre;\ + done = 0;\ + } else if(ez->err == 0) {\ + push_trace(&(ez->cigar), 0, ez->te+1-ez->ts); return;\ + } else {\ + pn = ez->pe + 1 - ez->ps; tn = ez->te + 1 - ez->ts;\ + }\ + int32_t i, err, tn0 = tn - 1, cut = thre+(thre<<1), bd, i_bd, k_bd, Peq_i, ws; w_sig c, Peq_m;\ + int32_t poff, pe = pn-1, tmp_e = INT32_MAX, k;\ + \ + w_##sf##_clear(ez->Peq[0]);\ + w_##sf##_clear(ez->Peq[1]);\ + w_##sf##_clear(ez->Peq[2]);\ + w_##sf##_clear(ez->Peq[3]);\ + w_##sf##_clear(ez->Peq[4]);\ + \ + bd = thre+1; bd = ((bd<=pn)?bd:pn); i_bd = thre;\ + ed_init_core(i, bd, i_bd, pstr, ez->Peq);\ + bd = thre+1, i_bd = thre;\ + \ + w_##sf##_clear(ez->Peq[4]);\ + err = thre;\ + w_##sf##_set_bit_lsub(ez->VN, thre); /**VN = (((Word)1)<<(thre))-1; **/\ + w_##sf##_set_bit_lsub(ez->VP, (thre<<1)+1); /**VP = (((Word)1)<<((thre<<1)+1))-1;**/\ + w_##sf##_self_xor(ez->VP, ez->VN); /**VP ^= VN;**/\ + \ + /** print_bits(ez->VP.a, (thre<<1)+1, "-VP");**/\ + /**should make Peq[4] = 0 if N is always an error**/\ + i = 0; ws = sizeof(*(ez->a))*(ez->nword);\ + ez->path.n=(ez->nword*tn*5);\ + kv_resize(w_sig, ez->path, ez->path.n); ez->path.n=0;\ + /**for the incoming char/last char; mm = ((Word)1 << (thre<<1))**/\ + Peq_i = (((thre<<1))>>bitw); Peq_m = (((w_sig)1)<<(((thre<<1))&bitz));\ + while (i < tn0) {\ + ed_core(_##sf##_, ez->Peq, ez->VP, ez->VN, ez->X, ez->D0, ez->HN, ez->HP, (uint8_t)tstr[i], c);\ + if (!(ez->D0.a[0]&(1ULL))) {\ + ++err; if (err>cut) return;\ + }\ + \ + if(!done) {\ + poff = i-thre; k = i+thre-pe;/**poff:[i-thre, i+thre]**/\ + if(k >= 0) {\ + if(k == 0) {\ + tmp_e = err;\ + for ((k) = 0; (poff) < (pe); (poff)++, (k)++) {\ + bd = (k>>bitw); k_bd = (k&bitz);\ + (tmp_e) += (((*ez).VP.a[bd]>>k_bd)&((w_sig)1));\ + (tmp_e) -= (((*ez).VN.a[bd]>>k_bd)&((w_sig)1));\ + }\ + } else {\ + k = (thre<<1) - k;\ + if(k >= 0) {\ + bd = (k>>bitw); k_bd = (k&bitz);\ + (tmp_e) += (((*ez).HP.a[bd]>>k_bd)&((w_sig)1));\ + (tmp_e) -= (((*ez).HN.a[bd]>>k_bd)&((w_sig)1));\ + }\ + }\ + if((tmp_e) <= (*ez).thre && (tmp_e) < (*ez).err) {\ + (*ez).err = tmp_e; (*ez).pe = pe; (*ez).te = i;\ + }\ + }\ + /**if(dbg_ext_err(i, thre, err, pe, &((*ez).VP), &((*ez).VN))!=tmp_e)**/\ + }\ + \ + /** Peq[0] >>= 1; Peq[1] >>= 1; Peq[2] >>= 1; Peq[3] >>= 1;**/\ + w_##sf##_self_rsft_1(ez->Peq[0]); w_##sf##_self_rsft_1(ez->Peq[1]);\ + w_##sf##_self_rsft_1(ez->Peq[2]); w_##sf##_self_rsft_1(ez->Peq[3]);\ + ++i; ++i_bd; c = 4;\ + if(i_bd < pn) c = seq_nt4_table[(uint8_t)pstr[i_bd]];\ + if(c < 4) ez->Peq[c].a[Peq_i]|=Peq_m;\ + \ + memcpy(ez->path.a+ez->path.n, ez->D0.a, ws); ez->path.n += ez->nword;\ + memcpy(ez->path.a+ez->path.n, ez->VP.a, ws); ez->path.n += ez->nword;\ + memcpy(ez->path.a+ez->path.n, ez->VN.a, ws); ez->path.n += ez->nword;\ + memcpy(ez->path.a+ez->path.n, ez->HP.a, ws); ez->path.n += ez->nword;\ + memcpy(ez->path.a+ez->path.n, ez->HN.a, ws); ez->path.n += ez->nword;\ + }\ + ed_core(_##sf##_, ez->Peq, ez->VP, ez->VN, ez->X, ez->D0, ez->HN, ez->HP, (uint8_t)tstr[i], c);\ + if (!(ez->D0.a[0]&(1ULL))) {\ + ++err; if (err>cut) return;\ + }\ + memcpy(ez->path.a+ez->path.n, ez->D0.a, ws); ez->path.n += ez->nword;\ + memcpy(ez->path.a+ez->path.n, ez->VP.a, ws); ez->path.n += ez->nword;\ + memcpy(ez->path.a+ez->path.n, ez->VN.a, ws); ez->path.n += ez->nword;\ + memcpy(ez->path.a+ez->path.n, ez->HP.a, ws); ez->path.n += ez->nword;\ + memcpy(ez->path.a+ez->path.n, ez->HN.a, ws); ez->path.n += ez->nword;\ + \ + int32_t site = tn - 1 - thre;/**up bound; site:[tn - 1 - thre, tn - 1 + thre]**/\ + if(!done) {\ + for (cut = pn - 1, i = 0; site < cut; i++) {\ + bd = (i>>bitw); i_bd = (i&bitz);\ + err += (((*ez).VP.a[bd]>>i_bd)&((w_sig)1));\ + err -= (((*ez).VN.a[bd]>>i_bd)&((w_sig)1));\ + site++;\ + if(err <= thre && err < ez->err) {\ + ez->err = err; ez->pe = site; ez->te = tn-1;\ + }\ + }\ + if(err <= thre && err < ez->err) {\ + ez->err = err; ez->pe = site; ez->te = tn-1;\ + }\ + }\ + \ + if((ez->te-ez->ts+1) != tn) {\ + ez->path.n /= tn; ez->path.n *= (ez->te+1-ez->ts);\ + }\ + gen_trace(ez, thre, 1);\ + return;\ +}\ +inline void ed_band_cal_extension_##sf##_1_w_trace(char *pstr, int32_t pn, char *tstr, int32_t tn, int32_t thre, bit_extz_t *ez)\ +{\ + int32_t done = 1;\ + ez->cigar.n = 0; ez->nword = w_##sf##_word;\ + if(ez->err > thre) {\ + init_base_ed(*ez, thre, pn, tn); ez->ps = ez->ts = INT32_MAX; ez->pe = pn-1; ez->te = tn-1;\ + if(pn > tn + thre) pn = tn + thre;\ + else if(tn > pn + thre) tn = pn + thre;\ + done = 0;\ + } else if(ez->err == 0) {\ + push_trace(&(ez->cigar), 0, ez->te+1-ez->ts); return;\ + } else {\ + pn = ez->pe + 1 - ez->ps; tn = ez->te + 1 - ez->ts;\ + }\ + int32_t i, err, tn0 = tn - 1, cut = thre+(thre<<1), bd, i_bd, k_bd, Peq_i, ws; w_sig c, Peq_m;\ + int32_t poff, pe = pn-1, tmp_e = INT32_MAX, k, pidx = ez->pe, tidx = ez->te;\ + \ + w_##sf##_clear(ez->Peq[0]);\ + w_##sf##_clear(ez->Peq[1]);\ + w_##sf##_clear(ez->Peq[2]);\ + w_##sf##_clear(ez->Peq[3]);\ + w_##sf##_clear(ez->Peq[4]);\ + \ + bd = thre+1; bd = ((bd<=pn)?bd:pn); i_bd = thre;\ + for (i = 0; i < bd; i++, i_bd++) {\ + w_bit((ez->Peq)[seq_nt4_table[(uint8_t)(pstr)[pidx-i]]], (i_bd));\ + }\ + bd = thre+1, i_bd = thre;\ + \ + w_##sf##_clear(ez->Peq[4]);\ + err = thre;\ + w_##sf##_set_bit_lsub(ez->VN, thre); /**VN = (((Word)1)<<(thre))-1; **/\ + w_##sf##_set_bit_lsub(ez->VP, (thre<<1)+1); /**VP = (((Word)1)<<((thre<<1)+1))-1;**/\ + w_##sf##_self_xor(ez->VP, ez->VN); /**VP ^= VN;**/\ + \ + /** print_bits(ez->VP.a, (thre<<1)+1, "-VP");**/\ + /**should make Peq[4] = 0 if N is always an error**/\ + i = 0; ws = sizeof(*(ez->a))*(ez->nword);\ + ez->path.n=(ez->nword*tn*5);\ + kv_resize(w_sig, ez->path, ez->path.n); ez->path.n=0;\ + /**for the incoming char/last char; mm = ((Word)1 << (thre<<1))**/\ + Peq_i = (((thre<<1))>>bitw); Peq_m = (((w_sig)1)<<(((thre<<1))&bitz));\ + while (i < tn0) {\ + ed_core(_##sf##_, ez->Peq, ez->VP, ez->VN, ez->X, ez->D0, ez->HN, ez->HP, (uint8_t)tstr[tidx-i], c);\ + if (!(ez->D0.a[0]&(1ULL))) {\ + ++err; if (err>cut) return;\ + }\ + \ + if(!done) {\ + poff = i-thre; k = i+thre-pe;/**poff:[i-thre, i+thre]**/\ + if(k >= 0) {\ + if(k == 0) {\ + tmp_e = err;\ + for ((k) = 0; (poff) < (pe); (poff)++, (k)++) {\ + bd = (k>>bitw); k_bd = (k&bitz);\ + (tmp_e) += (((*ez).VP.a[bd]>>k_bd)&((w_sig)1));\ + (tmp_e) -= (((*ez).VN.a[bd]>>k_bd)&((w_sig)1));\ + }\ + } else {\ + k = (thre<<1) - k;\ + if(k >= 0) {\ + bd = (k>>bitw); k_bd = (k&bitz);\ + (tmp_e) += (((*ez).HP.a[bd]>>k_bd)&((w_sig)1));\ + (tmp_e) -= (((*ez).HN.a[bd]>>k_bd)&((w_sig)1));\ + }\ + }\ + if((tmp_e) <= (*ez).thre && (tmp_e) < (*ez).err) {\ + (*ez).err = tmp_e; (*ez).ps = pidx - pe; (*ez).ts = tidx-i;\ + }\ + }\ + /**if(dbg_ext_err(i, thre, err, pe, &((*ez).VP), &((*ez).VN))!=tmp_e)**/\ + }\ + \ + /** Peq[0] >>= 1; Peq[1] >>= 1; Peq[2] >>= 1; Peq[3] >>= 1;**/\ + w_##sf##_self_rsft_1(ez->Peq[0]); w_##sf##_self_rsft_1(ez->Peq[1]);\ + w_##sf##_self_rsft_1(ez->Peq[2]); w_##sf##_self_rsft_1(ez->Peq[3]);\ + ++i; ++i_bd; c = 4;\ + if(i_bd < pn) c = seq_nt4_table[(uint8_t)pstr[pidx-i_bd]];\ + if(c < 4) ez->Peq[c].a[Peq_i]|=Peq_m;\ + \ + memcpy(ez->path.a+ez->path.n, ez->D0.a, ws); ez->path.n += ez->nword;\ + memcpy(ez->path.a+ez->path.n, ez->VP.a, ws); ez->path.n += ez->nword;\ + memcpy(ez->path.a+ez->path.n, ez->VN.a, ws); ez->path.n += ez->nword;\ + memcpy(ez->path.a+ez->path.n, ez->HP.a, ws); ez->path.n += ez->nword;\ + memcpy(ez->path.a+ez->path.n, ez->HN.a, ws); ez->path.n += ez->nword;\ + }\ + ed_core(_##sf##_, ez->Peq, ez->VP, ez->VN, ez->X, ez->D0, ez->HN, ez->HP, (uint8_t)tstr[tidx-i], c);\ + if (!(ez->D0.a[0]&(1ULL))) {\ + ++err; if (err>cut) return;\ + }\ + memcpy(ez->path.a+ez->path.n, ez->D0.a, ws); ez->path.n += ez->nword;\ + memcpy(ez->path.a+ez->path.n, ez->VP.a, ws); ez->path.n += ez->nword;\ + memcpy(ez->path.a+ez->path.n, ez->VN.a, ws); ez->path.n += ez->nword;\ + memcpy(ez->path.a+ez->path.n, ez->HP.a, ws); ez->path.n += ez->nword;\ + memcpy(ez->path.a+ez->path.n, ez->HN.a, ws); ez->path.n += ez->nword;\ + \ + int32_t site = tn - 1 - thre;/**up bound; site:[tn - 1 - thre, tn - 1 + thre]**/\ + if(!done) {\ + for (cut = pn - 1, i = 0; site < cut; i++) {\ + bd = (i>>bitw); i_bd = (i&bitz);\ + err += (((*ez).VP.a[bd]>>i_bd)&((w_sig)1));\ + err -= (((*ez).VN.a[bd]>>i_bd)&((w_sig)1));\ + site++;\ + if(err <= thre && err < ez->err) {\ + ez->err = err; ez->ps = pidx-site; ez->ts = tidx+1-tn;\ + }\ + }\ + if(err <= thre && err < ez->err) {\ + ez->err = err; ez->ps = pidx-site; ez->ts = tidx+1-tn;\ + }\ + }\ + \ + if((ez->te-ez->ts+1) != tn) {\ + ez->path.n /= tn; ez->path.n *= (ez->te+1-ez->ts);\ + }\ + \ + poff = ez->ps; ez->ps = pidx - ez->pe; ez->pe = pidx - poff;\ + gen_trace(ez, thre, 0);\ + poff = ez->ps; ez->ps = pidx - ez->pe; ez->pe = pidx - poff;\ + return;\ +}\ HA_ED_INIT(128) HA_ED_INIT(192) @@ -1279,13 +1940,13 @@ inline void ed_band_cal_semi_infi_w(char *pstr, int32_t pn, char *tstr, int32_t if((uge <= thre) && (uge == ez->err)) ez->pe = site + thre; } -inline void ed_band_cal_extension_infi0_w(char *pstr, int32_t pn, char *tstr, int32_t tn, int32_t thre, int32_t *nword, bit_extz_t *ez) +inline void ed_band_cal_extension_infi_0_w(char *pstr, int32_t pn, char *tstr, int32_t tn, int32_t thre, int32_t *nword, bit_extz_t *ez) { init_base_ed(*ez, thre, pn, tn); ez->ps = ez->ts = 0; ez->pe = ez->te = -1; if(pn > tn + thre) pn = tn + thre; else if(tn > pn + thre) tn = pn + thre; int32_t i, err, tn0 = tn - 1, cut = thre+(thre<<1), bd, i_bd, k_bd, wz, wz1, Peq_i; w_sig c, ad, Peq_m; - int32_t poff, pe = pn-1, tmp_e, k; + int32_t poff, pe = pn-1, tmp_e = INT32_MAX, k; if(nword) { ez->nword = (*nword); } else { @@ -1326,16 +1987,30 @@ inline void ed_band_cal_extension_infi0_w(char *pstr, int32_t pn, char *tstr, in } { - poff = i-thre; tmp_e = err; /**poff:[i-thre, i+thre]**/ - if((poff) + (((*ez).thre)<<1) >= (pe)) { - for ((k) = 0; (poff) < (pe); (poff)++) { - bd = (k>>bitw); k_bd = (k&bitz); - (tmp_e) += (((*ez).VP.a[bd]>>k_bd)&((w_sig)1)); - (tmp_e) -= (((*ez).VN.a[bd]>>k_bd)&((w_sig)1)); - (k)++;} - if((tmp_e) <= (*ez).thre && (tmp_e) < (*ez).err) { - (*ez).err = tmp_e; (*ez).pe = poff; (*ez).te = i;} - } + poff = i-thre; k = i+thre-pe;/**poff:[i-thre, i+thre]**/ + if(k >= 0) { + if(k == 0) { + tmp_e = err; + for ((k) = 0; (poff) < (pe); (poff)++, (k)++) { + bd = (k>>bitw); k_bd = (k&bitz); + (tmp_e) += (((*ez).VP.a[bd]>>k_bd)&((w_sig)1)); + (tmp_e) -= (((*ez).VN.a[bd]>>k_bd)&((w_sig)1)); + } + } else { + k = (thre<<1) - k; + if(k >= 0) { + bd = (k>>bitw); k_bd = (k&bitz); + (tmp_e) += (((*ez).HP.a[bd]>>k_bd)&((w_sig)1)); + (tmp_e) -= (((*ez).HN.a[bd]>>k_bd)&((w_sig)1)); + } + } + if((tmp_e) <= (*ez).thre && (tmp_e) < (*ez).err) { + (*ez).err = tmp_e; (*ez).pe = pe; (*ez).te = i; + } + } + // if(dbg_ext_err(i, thre, err, pe, &((*ez).VP), &((*ez).VN))!=tmp_e) { + // fprintf(stderr, "i::%d, tmp_e::%d, dbg_ext_err::%d\n", i, tmp_e, dbg_ext_err(i, thre, err, pe, &((*ez).VP), &((*ez).VN))); + // } } ed_infi_post_Peq(ez->Peq, wz, wz1, ez->nword); @@ -1364,13 +2039,13 @@ inline void ed_band_cal_extension_infi0_w(char *pstr, int32_t pn, char *tstr, in return; } -inline void ed_band_cal_extension_infi1_w(char *pstr, int32_t pn, char *tstr, int32_t tn, int32_t thre, int32_t *nword, bit_extz_t *ez) +inline void ed_band_cal_extension_infi_1_w(char *pstr, int32_t pn, char *tstr, int32_t tn, int32_t thre, int32_t *nword, bit_extz_t *ez) { init_base_ed(*ez, thre, pn, tn); ez->ps = ez->ts = INT32_MAX; ez->pe = pn-1; ez->te = tn-1; if(pn > tn + thre) pn = tn + thre; else if(tn > pn + thre) tn = pn + thre; int32_t i, err, tn0 = tn - 1, cut = thre+(thre<<1), bd, i_bd, k_bd, wz, wz1, Peq_i; - int32_t poff, pe = pn-1, tmp_e, k, pidx = ez->pe, tidx = ez->te; w_sig c, ad, Peq_m; + int32_t poff, pe = pn-1, tmp_e = INT32_MAX, k, pidx = ez->pe, tidx = ez->te; w_sig c, ad, Peq_m; if(nword) { ez->nword = (*nword); } else { @@ -1406,17 +2081,31 @@ inline void ed_band_cal_extension_infi1_w(char *pstr, int32_t pn, char *tstr, in } { - poff = i-thre; tmp_e = err; /**poff:[i-thre, i+thre]**/ - if((poff) + (((*ez).thre)<<1) >= (pe)) { - for ((k) = 0; (poff) < (pe); (poff)++) { - bd = (k>>bitw); k_bd = (k&bitz); - (tmp_e) += (((*ez).VP.a[bd]>>k_bd)&((w_sig)1)); - (tmp_e) -= (((*ez).VN.a[bd]>>k_bd)&((w_sig)1)); - (k)++;} - if((tmp_e) <= (*ez).thre && (tmp_e) < (*ez).err) { - (*ez).err = tmp_e; (*ez).pe = pidx-poff; (*ez).te = tidx-i;} - } + poff = i-thre; k = i+thre-pe;/**poff:[i-thre, i+thre]**/ + if(k >= 0) { + if(k == 0) { + tmp_e = err; + for ((k) = 0; (poff) < (pe); (poff)++, (k)++) { + bd = (k>>bitw); k_bd = (k&bitz); + (tmp_e) += (((*ez).VP.a[bd]>>k_bd)&((w_sig)1)); + (tmp_e) -= (((*ez).VN.a[bd]>>k_bd)&((w_sig)1)); + } + } else { + k = (thre<<1) - k; + if(k >= 0) { + bd = (k>>bitw); k_bd = (k&bitz); + (tmp_e) += (((*ez).HP.a[bd]>>k_bd)&((w_sig)1)); + (tmp_e) -= (((*ez).HN.a[bd]>>k_bd)&((w_sig)1)); + } + } + if((tmp_e) <= (*ez).thre && (tmp_e) < (*ez).err) { + (*ez).err = tmp_e; (*ez).ps = pidx - pe; (*ez).ts = tidx-i; + } + } } + // if(dbg_ext_err(i, thre, err, pe, &((*ez).VP), &((*ez).VN))!=tmp_e) { + // fprintf(stderr, "i::%d, tmp_e::%d, dbg_ext_err::%d\n", i, tmp_e, dbg_ext_err(i, thre, err, pe, &((*ez).VP), &((*ez).VN))); + // } ed_infi_post_Peq(ez->Peq, wz, wz1, ez->nword); ++i; ++i_bd; c = 4; @@ -1435,15 +2124,448 @@ inline void ed_band_cal_extension_infi1_w(char *pstr, int32_t pn, char *tstr, in err -= (((*ez).VN.a[bd]>>i_bd)&((w_sig)1)); site++; if(err <= thre && err < ez->err) { - ez->err = err; ez->pe = pidx-site; ez->te = tidx+1-tn; + ez->err = err; ez->ps = pidx-site; ez->ts = tidx+1-tn; } } if(err <= thre && err < ez->err) { - ez->err = err; ez->pe = pidx-site; ez->te = tidx+1-tn; + ez->err = err; ez->ps = pidx-site; ez->ts = tidx+1-tn; } return; } +inline void ed_band_cal_global_infi_w_trace(char *pstr, int32_t pn, char *tstr, int32_t tn, int32_t thre, int32_t *nword, bit_extz_t *ez) +{ + ez->cigar.n = 0;//diff + if(ez->err > thre) {//diff + init_base_ed(*ez, thre, pn, tn); ez->ps = ez->ts = 0; + } else if(ez->err == 0) {//diff + push_trace(&(ez->cigar), 0, ez->te+1-ez->ts); return;//diff + } + if((pn > tn + thre) || (tn > pn + thre)) return; + int32_t i, err, tn0 = tn - 1, cut = thre+(thre<<1), bd, i_bd, wz, wz1, Peq_i, ws; w_sig c, ad, Peq_m; + if(nword) { + ez->nword = (*nword); + } else { + bd = (((thre)<<1)+1); ez->nword = ((bd>>bitw)+(!!(bd&bitz))); + } + resize_bit_extz_t((*ez), (thre)); + wz = ws = sizeof(*(ez->a))*(ez->nword);//diff + memset((ez->Peq[0]).a, 0, wz); + memset((ez->Peq[1]).a, 0, wz); + memset((ez->Peq[2]).a, 0, wz); + memset((ez->Peq[3]).a, 0, wz); + memset((ez->Peq[4]).a, 0, wz); + + bd = thre+1; bd = ((bd<=pn)?bd:pn); i_bd = thre; + ed_init_core(i, bd, i_bd, pstr, ez->Peq); + bd = thre+1, i_bd = thre; + + memset((ez->Peq[4]).a, 0, wz); + err = thre; + w_infi_set_bit_lsub(ez->VN, thre, ez->nword); /**VN = (((Word)1)<<(thre))-1;**/ + w_infi_set_bit_lsub(ez->VP, (thre<<1)+1, ez->nword); /**VP = (((Word)1)<<((thre<<1)+1))-1;**/ + w_infi_self_xor(ez->VP, ez->VN, ez->nword, wz) /**VP ^= VN;**/ + + ez->path.n=(ez->nword*tn*5);//diff + kv_resize(w_sig, ez->path, ez->path.n); ez->path.n=0;//diff + + i = 0; + /**for the incoming char/last char; mm = ((Word)1 << (thre<<1))**/ + Peq_i = (((thre<<1))>>bitw); Peq_m = (((w_sig)1)<<(((thre<<1))&bitz)); + while (i < tn0) { + // fprintf(stderr, "\ni::%d\n", i); + // prt_bit_extz_t((*ez), (((thre<<1))+1)); + ed_infi_core(ez->Peq, ez->VP, ez->VN, ez->X, ez->D0, ez->HN, ez->HP, (uint8_t)tstr[i], c, ad, wz, ez->nword); + if (!(ez->D0.a[0]&(1ULL))) { + ++err; if (err>cut) return; + } + + ed_infi_post_Peq(ez->Peq, wz, wz1, ez->nword); + ++i; ++i_bd; c = 4; + if(i_bd < pn) c = seq_nt4_table[(uint8_t)pstr[i_bd]]; + if(c < 4) ez->Peq[c].a[Peq_i]|=Peq_m; + // fprintf(stderr, "c::%u\n", c); + // print_bits(ez->Peq[c].a, (((thre<<1))+1), "Peq[c]"); + + memcpy(ez->path.a+ez->path.n, ez->D0.a, ws); ez->path.n += ez->nword;//diff + memcpy(ez->path.a+ez->path.n, ez->VP.a, ws); ez->path.n += ez->nword;//diff + memcpy(ez->path.a+ez->path.n, ez->VN.a, ws); ez->path.n += ez->nword;//diff + memcpy(ez->path.a+ez->path.n, ez->HP.a, ws); ez->path.n += ez->nword;//diff + memcpy(ez->path.a+ez->path.n, ez->HN.a, ws); ez->path.n += ez->nword;//diff + } + ed_infi_core(ez->Peq, ez->VP, ez->VN, ez->X, ez->D0, ez->HN, ez->HP, (uint8_t)tstr[i], c, ad, wz, ez->nword); + if (!(ez->D0.a[0]&(1ULL))) { + ++err; if (err>cut) return; + } + memcpy(ez->path.a+ez->path.n, ez->D0.a, ws); ez->path.n += ez->nword;//diff + memcpy(ez->path.a+ez->path.n, ez->VP.a, ws); ez->path.n += ez->nword;//diff + memcpy(ez->path.a+ez->path.n, ez->VN.a, ws); ez->path.n += ez->nword;//diff + memcpy(ez->path.a+ez->path.n, ez->HP.a, ws); ez->path.n += ez->nword;//diff + memcpy(ez->path.a+ez->path.n, ez->HN.a, ws); ez->path.n += ez->nword;//diff + + int32_t site = tn - 1 - thre;/**up bound**/ + if(ez->err > thre) {//diff + for (cut = pn - 1, i = 0; site < cut; site++, i++) { + bd = (i>>bitw); i_bd = (i&bitz); + err += (((*ez).VP.a[bd]>>i_bd)&((w_sig)1)); + err -= (((*ez).VN.a[bd]>>i_bd)&((w_sig)1)); + } + + if (site == cut && err <= thre) { + ez->err = err; + ez->pe = pn-1; ez->te = tn-1; + } + } + gen_trace(ez, thre, 1);//diff + return; +} + +inline void ed_band_cal_semi_infi_w_trace(char *pstr, int32_t pn, char *tstr, int32_t tn, int32_t thre, int32_t *nword, bit_extz_t *ez) +{ + ez->cigar.n = 0;//diff + if(ez->err > thre) {//diff + init_base_ed(*ez, thre, pn, tn); ez->ps = ez->pe = -1; ez->ts = 0; ez->te = tn-1; + } else if(ez->err == 0) {//diff + push_trace(&(ez->cigar), 0, ez->te+1-ez->ts); //diff + ez->ps = ez->pe - (ez->te-ez->ts);//diff + return;//diff + } + int32_t bd, i, err = 0, i_bd, cut = thre+(thre<<1), tn0 = tn - 1, wz, wz1, ws, Peq_i; w_sig c, ad, Peq_m; + if(nword) { + ez->nword = (*nword); + } else { + bd = (((thre)<<1)+1); ez->nword = ((bd>>bitw)+(!!(bd&bitz))); + } + resize_bit_extz_t((*ez), (thre)); + wz = ws = sizeof(*(ez->a))*(ez->nword);//diff + memset((ez->Peq[0]).a, 0, wz); + memset((ez->Peq[1]).a, 0, wz); + memset((ez->Peq[2]).a, 0, wz); + memset((ez->Peq[3]).a, 0, wz); + memset((ez->Peq[4]).a, 0, wz); + memset((ez->VP).a, 0, wz); + memset((ez->VN).a, 0, wz); + + bd = (thre<<1)+1; bd = ((bd<=pn)?bd:pn); i_bd = 0; + ed_init_core(i, bd, i_bd, pstr, ez->Peq); + bd = (thre<<1)+1, i_bd = (thre<<1); + + memset((ez->Peq[4]).a, 0, wz); + + ez->path.n=(ez->nword*tn*5);//diff + kv_resize(w_sig, ez->path, ez->path.n); ez->path.n=0;//diff + + i = 0; + /**for the incoming char/last char; mm = ((Word)1 << (thre<<1))**/ + Peq_i = (((thre<<1))>>bitw); Peq_m = (((w_sig)1)<<(((thre<<1))&bitz)); + while (i < tn0) { + ed_infi_core(ez->Peq, ez->VP, ez->VN, ez->X, ez->D0, ez->HN, ez->HP, (uint8_t)tstr[i], c, ad, wz, ez->nword); + if (!(ez->D0.a[0]&(1ULL))) { + ++err; if (err>cut) return; + } + + ed_infi_post_Peq(ez->Peq, wz, wz1, ez->nword); + ++i; ++i_bd; c = 4; + if(i_bd < pn) c = seq_nt4_table[(uint8_t)pstr[i_bd]]; + if(c < 4) ez->Peq[c].a[Peq_i]|=Peq_m; + + memcpy(ez->path.a+ez->path.n, ez->D0.a, ws); ez->path.n += ez->nword;//diff + memcpy(ez->path.a+ez->path.n, ez->VP.a, ws); ez->path.n += ez->nword;//diff + memcpy(ez->path.a+ez->path.n, ez->VN.a, ws); ez->path.n += ez->nword;//diff + memcpy(ez->path.a+ez->path.n, ez->HP.a, ws); ez->path.n += ez->nword;//diff + memcpy(ez->path.a+ez->path.n, ez->HN.a, ws); ez->path.n += ez->nword;//diff + } + ed_infi_core(ez->Peq, ez->VP, ez->VN, ez->X, ez->D0, ez->HN, ez->HP, (uint8_t)tstr[i], c, ad, wz, ez->nword); + if (!(ez->D0.a[0]&(1ULL))) { + ++err; if (err>cut) return; + } + memcpy(ez->path.a+ez->path.n, ez->D0.a, ws); ez->path.n += ez->nword;//diff + memcpy(ez->path.a+ez->path.n, ez->VP.a, ws); ez->path.n += ez->nword;//diff + memcpy(ez->path.a+ez->path.n, ez->VN.a, ws); ez->path.n += ez->nword;//diff + memcpy(ez->path.a+ez->path.n, ez->HP.a, ws); ez->path.n += ez->nword;//diff + memcpy(ez->path.a+ez->path.n, ez->HN.a, ws); ez->path.n += ez->nword;//diff + + int32_t site = tn - 1;/**up bound**/ + /**in most cases, ai = (thre<<1)**/ + int32_t ai = pn - tn, uge = INT32_MAX; + if(ez->err > thre) {//diff + if ((err <= thre) && (err <= ez->err)) { + ez->err = err; ez->pe = site; + } + + i = 0; + while (i < ai) { + bd = (i>>bitw); i_bd = (i&bitz); + err += (((*ez).VP.a[bd]>>i_bd)&((w_sig)1)); + err -= (((*ez).VN.a[bd]>>i_bd)&((w_sig)1)); + ++i; + if ((err <= thre) && (err <= ez->err)) { + ez->err = err; ez->pe = site + i; + } + if(i == thre) uge = err; + } + if((uge <= thre) && (uge == ez->err)) ez->pe = site + thre; + } + gen_trace(ez, 0, 1);//diff +} + +inline void ed_band_cal_extension_infi_0_w_trace(char *pstr, int32_t pn, char *tstr, int32_t tn, int32_t thre, int32_t *nword, bit_extz_t *ez) +{ + int32_t done = 1;//diff + ez->cigar.n = 0; //diff + if(ez->err > thre) {//diff + init_base_ed(*ez, thre, pn, tn); ez->ps = ez->ts = 0; ez->pe = ez->te = -1; + if(pn > tn + thre) pn = tn + thre; + else if(tn > pn + thre) tn = pn + thre; + done = 0;//diff + } else if(ez->err == 0) {//diff + push_trace(&(ez->cigar), 0, ez->te+1-ez->ts); return;//diff + } else { + pn = ez->pe + 1 - ez->ps; tn = ez->te + 1 - ez->ts; //diff + } + int32_t i, err, tn0 = tn - 1, cut = thre+(thre<<1), bd, i_bd, k_bd, wz, wz1, ws, Peq_i; w_sig c, ad, Peq_m; + int32_t poff, pe = pn-1, tmp_e = INT32_MAX, k; + if(nword) { + ez->nword = (*nword); + } else { + bd = (((thre)<<1)+1); ez->nword = ((bd>>bitw)+(!!(bd&bitz))); + } + resize_bit_extz_t((*ez), (thre)); + wz = ws = sizeof(*(ez->a))*(ez->nword);//diff + memset((ez->Peq[0]).a, 0, wz); + memset((ez->Peq[1]).a, 0, wz); + memset((ez->Peq[2]).a, 0, wz); + memset((ez->Peq[3]).a, 0, wz); + memset((ez->Peq[4]).a, 0, wz); + + bd = thre+1; bd = ((bd<=pn)?bd:pn); i_bd = thre; + ed_init_core(i, bd, i_bd, pstr, ez->Peq); + bd = thre+1, i_bd = thre; + + memset((ez->Peq[4]).a, 0, wz); + err = thre; + w_infi_set_bit_lsub(ez->VN, thre, ez->nword); /**VN = (((Word)1)<<(thre))-1;**/ + w_infi_set_bit_lsub(ez->VP, (thre<<1)+1, ez->nword); /**VP = (((Word)1)<<((thre<<1)+1))-1;**/ + w_infi_self_xor(ez->VP, ez->VN, ez->nword, wz) /**VP ^= VN;**/ + + ez->path.n=(ez->nword*tn*5);//diff + kv_resize(w_sig, ez->path, ez->path.n); ez->path.n=0;//diff + + i = 0; + /**for the incoming char/last char; mm = ((Word)1 << (thre<<1))**/ + Peq_i = (((thre<<1))>>bitw); Peq_m = (((w_sig)1)<<(((thre<<1))&bitz)); + //for debug + // memset((ez->X).a, 0, sizeof(*(ez->a))*(ez->nword)); + // memset((ez->D0).a, 0, sizeof(*(ez->a))*(ez->nword)); + // memset((ez->HN).a, 0, sizeof(*(ez->a))*(ez->nword)); + // memset((ez->HP).a, 0, sizeof(*(ez->a))*(ez->nword)); + while (i < tn0) { + // fprintf(stderr, "i::%d\n", i); + // prt_bit_extz_t((*ez), ((thre<<1)+1)); + ed_infi_core(ez->Peq, ez->VP, ez->VN, ez->X, ez->D0, ez->HN, ez->HP, (uint8_t)tstr[i], c, ad, wz, ez->nword); + if (!(ez->D0.a[0]&(1ULL))) { + ++err; if (err>cut) return; + } + + if(!done) {//diff + poff = i-thre; k = i+thre-pe;/**poff:[i-thre, i+thre]**/ + if(k >= 0) { + if(k == 0) { + tmp_e = err; + for ((k) = 0; (poff) < (pe); (poff)++, (k)++) { + bd = (k>>bitw); k_bd = (k&bitz); + (tmp_e) += (((*ez).VP.a[bd]>>k_bd)&((w_sig)1)); + (tmp_e) -= (((*ez).VN.a[bd]>>k_bd)&((w_sig)1)); + } + } else { + k = (thre<<1) - k; + if(k >= 0) { + bd = (k>>bitw); k_bd = (k&bitz); + (tmp_e) += (((*ez).HP.a[bd]>>k_bd)&((w_sig)1)); + (tmp_e) -= (((*ez).HN.a[bd]>>k_bd)&((w_sig)1)); + } + } + if((tmp_e) <= (*ez).thre && (tmp_e) < (*ez).err) { + (*ez).err = tmp_e; (*ez).pe = pe; (*ez).te = i; + } + } + // if(dbg_ext_err(i, thre, err, pe, &((*ez).VP), &((*ez).VN))!=tmp_e) { + // fprintf(stderr, "i::%d, tmp_e::%d, dbg_ext_err::%d\n", i, tmp_e, dbg_ext_err(i, thre, err, pe, &((*ez).VP), &((*ez).VN))); + // } + } + + ed_infi_post_Peq(ez->Peq, wz, wz1, ez->nword); + ++i; ++i_bd; c = 4; + if(i_bd < pn) c = seq_nt4_table[(uint8_t)pstr[i_bd]]; + if(c < 4) ez->Peq[c].a[Peq_i]|=Peq_m; + + memcpy(ez->path.a+ez->path.n, ez->D0.a, ws); ez->path.n += ez->nword;//diff + memcpy(ez->path.a+ez->path.n, ez->VP.a, ws); ez->path.n += ez->nword;//diff + memcpy(ez->path.a+ez->path.n, ez->VN.a, ws); ez->path.n += ez->nword;//diff + memcpy(ez->path.a+ez->path.n, ez->HP.a, ws); ez->path.n += ez->nword;//diff + memcpy(ez->path.a+ez->path.n, ez->HN.a, ws); ez->path.n += ez->nword;//diff + } + ed_infi_core(ez->Peq, ez->VP, ez->VN, ez->X, ez->D0, ez->HN, ez->HP, (uint8_t)tstr[i], c, ad, wz, ez->nword); + if (!(ez->D0.a[0]&(1ULL))) { + ++err; if (err>cut) return; + } + memcpy(ez->path.a+ez->path.n, ez->D0.a, ws); ez->path.n += ez->nword;//diff + memcpy(ez->path.a+ez->path.n, ez->VP.a, ws); ez->path.n += ez->nword;//diff + memcpy(ez->path.a+ez->path.n, ez->VN.a, ws); ez->path.n += ez->nword;//diff + memcpy(ez->path.a+ez->path.n, ez->HP.a, ws); ez->path.n += ez->nword;//diff + memcpy(ez->path.a+ez->path.n, ez->HN.a, ws); ez->path.n += ez->nword;//diff + + int32_t site = tn - 1 - thre;/**up bound; site:[tn - 1 - thre, tn - 1 + thre]**/ + if(!done) {//diff + for (cut = pn - 1, i = 0; site < cut; i++) { + bd = (i>>bitw); i_bd = (i&bitz); + err += (((*ez).VP.a[bd]>>i_bd)&((w_sig)1)); + err -= (((*ez).VN.a[bd]>>i_bd)&((w_sig)1)); + site++; + if(err <= thre && err < ez->err) { + ez->err = err; ez->pe = site; ez->te = tn-1; + } + } + if(err <= thre && err < ez->err) { + ez->err = err; ez->pe = site; ez->te = tn-1; + } + } + + if((ez->te-ez->ts+1) != tn) {//diff + ez->path.n /= tn; ez->path.n *= (ez->te+1-ez->ts);//diff + }//diff + gen_trace(ez, thre, 1);//diff + return; +} + +inline void ed_band_cal_extension_infi_1_w_trace(char *pstr, int32_t pn, char *tstr, int32_t tn, int32_t thre, int32_t *nword, bit_extz_t *ez) +{ + int32_t done = 1;//diff + ez->cigar.n = 0;//diff + if(ez->err > thre) {//diff + init_base_ed(*ez, thre, pn, tn); ez->ps = ez->ts = INT32_MAX; ez->pe = pn-1; ez->te = tn-1; + if(pn > tn + thre) pn = tn + thre; + else if(tn > pn + thre) tn = pn + thre; + done = 0;//diff + } else if(ez->err == 0) {//diff + push_trace(&(ez->cigar), 0, ez->te+1-ez->ts); return;//diff + } else { + pn = ez->pe + 1 - ez->ps; tn = ez->te + 1 - ez->ts; //diff + } + int32_t i, err, tn0 = tn - 1, cut = thre+(thre<<1), bd, i_bd, k_bd, wz, wz1, ws, Peq_i; + int32_t poff, pe = pn-1, tmp_e = INT32_MAX, k, pidx = ez->pe, tidx = ez->te; w_sig c, ad, Peq_m; + if(nword) { + ez->nword = (*nword); + } else { + bd = (((thre)<<1)+1); ez->nword = ((bd>>bitw)+(!!(bd&bitz))); + } + resize_bit_extz_t((*ez), (thre)); + wz = ws = sizeof(*(ez->a))*(ez->nword);//diff + memset((ez->Peq[0]).a, 0, wz); + memset((ez->Peq[1]).a, 0, wz); + memset((ez->Peq[2]).a, 0, wz); + memset((ez->Peq[3]).a, 0, wz); + memset((ez->Peq[4]).a, 0, wz); + + bd = thre+1; bd = ((bd<=pn)?bd:pn); i_bd = thre; + for (i = 0; i < bd; i++, i_bd++) { + w_bit((ez->Peq)[seq_nt4_table[(uint8_t)(pstr)[pidx-i]]], (i_bd)); + } + bd = thre+1, i_bd = thre; + + memset((ez->Peq[4]).a, 0, wz); + err = thre; + w_infi_set_bit_lsub(ez->VN, thre, ez->nword); /**VN = (((Word)1)<<(thre))-1;**/ + w_infi_set_bit_lsub(ez->VP, (thre<<1)+1, ez->nword); /**VP = (((Word)1)<<((thre<<1)+1))-1;**/ + w_infi_self_xor(ez->VP, ez->VN, ez->nword, wz) /**VP ^= VN;**/ + + ez->path.n=(ez->nword*tn*5);//diff + kv_resize(w_sig, ez->path, ez->path.n); ez->path.n=0;//diff + + i = 0; + /**for the incoming char/last char; mm = ((Word)1 << (thre<<1))**/ + Peq_i = (((thre<<1))>>bitw); Peq_m = (((w_sig)1)<<(((thre<<1))&bitz)); + while (i < tn0) { + ed_infi_core(ez->Peq, ez->VP, ez->VN, ez->X, ez->D0, ez->HN, ez->HP, (uint8_t)tstr[tidx-i], c, ad, wz, ez->nword); + if (!(ez->D0.a[0]&(1ULL))) { + ++err; if (err>cut) return; + } + + if(!done) {//diff + poff = i-thre; k = i+thre-pe;/**poff:[i-thre, i+thre]**/ + if(k >= 0) { + if(k == 0) { + tmp_e = err; + for ((k) = 0; (poff) < (pe); (poff)++, (k)++) { + bd = (k>>bitw); k_bd = (k&bitz); + (tmp_e) += (((*ez).VP.a[bd]>>k_bd)&((w_sig)1)); + (tmp_e) -= (((*ez).VN.a[bd]>>k_bd)&((w_sig)1)); + } + } else { + k = (thre<<1) - k; + if(k >= 0) { + bd = (k>>bitw); k_bd = (k&bitz); + (tmp_e) += (((*ez).HP.a[bd]>>k_bd)&((w_sig)1)); + (tmp_e) -= (((*ez).HN.a[bd]>>k_bd)&((w_sig)1)); + } + } + if((tmp_e) <= (*ez).thre && (tmp_e) < (*ez).err) { + (*ez).err = tmp_e; (*ez).ps = pidx - pe; (*ez).ts = tidx-i; + } + } + } + // if(dbg_ext_err(i, thre, err, pe, &((*ez).VP), &((*ez).VN))!=tmp_e) { + // fprintf(stderr, "i::%d, tmp_e::%d, dbg_ext_err::%d\n", i, tmp_e, dbg_ext_err(i, thre, err, pe, &((*ez).VP), &((*ez).VN))); + // } + + ed_infi_post_Peq(ez->Peq, wz, wz1, ez->nword); + ++i; ++i_bd; c = 4; + if(i_bd < pn) c = seq_nt4_table[(uint8_t)pstr[pidx-i_bd]]; + if(c < 4) ez->Peq[c].a[Peq_i]|=Peq_m; + + memcpy(ez->path.a+ez->path.n, ez->D0.a, ws); ez->path.n += ez->nword;//diff + memcpy(ez->path.a+ez->path.n, ez->VP.a, ws); ez->path.n += ez->nword;//diff + memcpy(ez->path.a+ez->path.n, ez->VN.a, ws); ez->path.n += ez->nword;//diff + memcpy(ez->path.a+ez->path.n, ez->HP.a, ws); ez->path.n += ez->nword;//diff + memcpy(ez->path.a+ez->path.n, ez->HN.a, ws); ez->path.n += ez->nword;//diff + } + ed_infi_core(ez->Peq, ez->VP, ez->VN, ez->X, ez->D0, ez->HN, ez->HP, (uint8_t)tstr[tidx-i], c, ad, wz, ez->nword); + if (!(ez->D0.a[0]&(1ULL))) { + ++err; if (err>cut) return; + } + memcpy(ez->path.a+ez->path.n, ez->D0.a, ws); ez->path.n += ez->nword;//diff + memcpy(ez->path.a+ez->path.n, ez->VP.a, ws); ez->path.n += ez->nword;//diff + memcpy(ez->path.a+ez->path.n, ez->VN.a, ws); ez->path.n += ez->nword;//diff + memcpy(ez->path.a+ez->path.n, ez->HP.a, ws); ez->path.n += ez->nword;//diff + memcpy(ez->path.a+ez->path.n, ez->HN.a, ws); ez->path.n += ez->nword;//diff + + int32_t site = tn - 1 - thre;/**up bound; site:[tn - 1 - thre, tn - 1 + thre]**/ + if(!done) {//diff + for (cut = pn - 1, i = 0; site < cut; i++) { + bd = (i>>bitw); i_bd = (i&bitz); + err += (((*ez).VP.a[bd]>>i_bd)&((w_sig)1)); + err -= (((*ez).VN.a[bd]>>i_bd)&((w_sig)1)); + site++; + if(err <= thre && err < ez->err) { + ez->err = err; ez->ps = pidx-site; ez->ts = tidx+1-tn; + } + } + if(err <= thre && err < ez->err) { + ez->err = err; ez->ps = pidx-site; ez->ts = tidx+1-tn; + } + } + + if((ez->te-ez->ts+1) != tn) {//diff + ez->path.n /= tn; ez->path.n *= (ez->te+1-ez->ts);//diff + }//diff + + poff = ez->ps; ez->ps = pidx - ez->pe; ez->pe = pidx - poff;//diff + gen_trace(ez, thre, 0);//diff + poff = ez->ps; ez->ps = pidx - ez->pe; ez->pe = pidx - poff;//diff + return; +} + + #define ed_core_64(Peq, VP, VN, X, D0, HN, HP, z) { \ /**X = Peq[seq_nt4_table[(uint8_t)tstr[i]]] | VN;**/\ (X) = (Peq)[(seq_nt4_table[(z)])]|(VN);\ @@ -1551,12 +2673,12 @@ inline void ed_band_cal_semi_64_w(char *pstr, int32_t pn, char *tstr, int32_t tn if((uge <= thre) && (uge == ez->err)) ez->pe = site + thre; } -inline void ed_band_cal_extension_64_w(char *pstr, int32_t pn, char *tstr, int32_t tn, int32_t thre, bit_extz_t *ez) +inline void ed_band_cal_extension_64_0_w(char *pstr, int32_t pn, char *tstr, int32_t tn, int32_t thre, bit_extz_t *ez) { + init_base_ed(*ez, thre, pn, tn); ez->ps = ez->ts = 0; ez->pe = ez->te = -1; if(pn > tn + thre) pn = tn + thre; else if(tn > pn + thre) tn = pn + thre; - init_base_ed(*ez, thre, pn, tn); ez->ps = ez->ts = 0; ez->pe = ez->te = -1; - int32_t i, err, tn0 = tn - 1, cut = thre+(thre<<1), bd, i_bd, poff, pe = pn-1, tmp_e, k; + int32_t i, err, tn0 = tn - 1, cut = thre+(thre<<1), bd, i_bd, poff, pe = pn-1, tmp_e = INT32_MAX, k; Word c, Peq[5] = {0}, VP, VN, X, D0, HN, HP, mm; bd = thre+1; bd = ((bd<=pn)?bd:pn); i_bd = thre; @@ -1577,16 +2699,25 @@ inline void ed_band_cal_extension_64_w(char *pstr, int32_t pn, char *tstr, int32 } { - poff = i-thre; tmp_e = err; /**poff:[i-thre, i+thre]**/ - if(poff + (((*ez).thre)<<1) >= pe) { - for (k = 0; poff < pe; poff++) { - tmp_e += ((VP>>k)&(1ULL)); tmp_e -= ((VN>>k)&(1ULL)); (k)++; - } + poff = i-thre; k = i+thre-pe;/**poff:[i-thre, i+thre]**/ + if(k >= 0) { + if(k == 0) { + tmp_e = err; + for (k = 0; poff < pe; poff++, (k)++) { + tmp_e += ((VP>>k)&(1ULL)); tmp_e -= ((VN>>k)&(1ULL)); + } + } else { + k = (thre<<1) - k; + if(k >= 0) { + tmp_e += ((HP>>k)&(1ULL)); tmp_e -= ((HN>>k)&(1ULL)); + } + } if(tmp_e <= (*ez).thre && tmp_e < (*ez).err) { - (*ez).err = tmp_e; (*ez).pe = poff; (*ez).te = i; - } + (*ez).err = tmp_e; (*ez).pe = pe; (*ez).te = i; + } } } + // fprintf(stderr, "i::%d, tmp_e::%d, err::%d\n", i, tmp_e, err); Peq[0] >>= 1; Peq[1] >>= 1; Peq[2] >>= 1; Peq[3] >>= 1; ++i; ++i_bd; @@ -1614,6 +2745,429 @@ inline void ed_band_cal_extension_64_w(char *pstr, int32_t pn, char *tstr, int32 return; } +inline void ed_band_cal_extension_64_1_w(char *pstr, int32_t pn, char *tstr, int32_t tn, int32_t thre, bit_extz_t *ez) +{ + init_base_ed(*ez, thre, pn, tn); ez->ps = ez->ts = INT32_MAX; ez->pe = pn-1; ez->te = tn-1; + if(pn > tn + thre) pn = tn + thre; + else if(tn > pn + thre) tn = pn + thre; + int32_t i, err, tn0 = tn - 1, cut = thre+(thre<<1), bd, i_bd; + int32_t poff, pe = pn-1, tmp_e = INT32_MAX, k, pidx = ez->pe, tidx = ez->te; + Word c, Peq[5] = {0}, VP, VN, X, D0, HN, HP, mm; + + bd = thre+1; bd = ((bd<=pn)?bd:pn); i_bd = thre; + for (i = 0, mm = (((Word)1)<cut) return; + } + + { + poff = i-thre; k = i+thre-pe;/**poff:[i-thre, i+thre]**/ + if(k >= 0) { + if(k == 0) { + tmp_e = err; + for (k = 0; poff < pe; poff++, (k)++) { + tmp_e += ((VP>>k)&(1ULL)); tmp_e -= ((VN>>k)&(1ULL)); + } + } else { + k = (thre<<1) - k; + if(k >= 0) { + tmp_e += ((HP>>k)&(1ULL)); tmp_e -= ((HN>>k)&(1ULL)); + } + } + if(tmp_e <= (*ez).thre && tmp_e < (*ez).err) { + (*ez).err = tmp_e; (*ez).ps = pidx - pe; (*ez).ts = tidx-i; + } + } + } + // fprintf(stderr, "i::%d, tmp_e::%d, err::%d\n", i, tmp_e, err); + + Peq[0] >>= 1; Peq[1] >>= 1; Peq[2] >>= 1; Peq[3] >>= 1; + ++i; ++i_bd; + if(i_bd < pn) { + c = seq_nt4_table[(uint8_t)pstr[pidx-i_bd]]; + if(c < 4) Peq[c] |= mm; + } + } + ed_core_64(Peq, VP, VN, X, D0, HN, HP, (uint8_t)tstr[tidx-i]); + if (!(D0&(1ULL))) { + ++err; if (err>cut) return; + } + + int32_t site = tn - 1 - thre;/**up bound; site:[tn - 1 - thre, tn - 1 + thre]**/\ + for (cut = pn - 1; site < cut; ) { + err += VP&(1ULL); VP >>= 1; err -= VN&(1ULL); VN >>= 1; + site++; + if(err <= thre && err < ez->err) { + ez->err = err; ez->ps = pidx-site; ez->ts = tidx+1-tn; + } + } + if(err <= thre && err < ez->err) { + ez->err = err; ez->ps = pidx-site; ez->ts = tidx+1-tn; + } + return; +} + +inline void ed_band_cal_global_64_w_trace(char *pstr, int32_t pn, char *tstr, int32_t tn, int32_t thre, bit_extz_t *ez) +{ + ez->cigar.n = 0; ez->nword = 1;//diff + if(ez->err > thre) {//diff + init_base_ed(*ez, thre, pn, tn); ez->ps = ez->ts = 0;//diff + } else if(ez->err == 0) {//diff + push_trace(&(ez->cigar), 0, ez->te+1-ez->ts); return;//diff + } + if((pn > tn + thre) || (tn > pn + thre)) return; + int32_t i, err, tn0 = tn - 1, cut = thre+(thre<<1), bd, i_bd; + w_sig c, Peq[5] = {0}, VP, VN, X, D0, HN, HP, mm; + + bd = thre+1; bd = ((bd<=pn)?bd:pn); i_bd = thre; + for (i = 0, mm = (((w_sig)1)<path.n=(ez->nword*tn*5);//diff + kv_resize(w_sig, ez->path, ez->path.n); ez->path.n=0;//diff + + i = 0; mm = ((w_sig)1 << (thre<<1));///for the incoming char/last char + while (i < tn0) { + // fprintf(stderr, "\ni::%d\n", i); + // prt_vector(Peq, VP, VN, X, D0, HN, HP, (((thre<<1))+1)); + ed_core_64(Peq, VP, VN, X, D0, HN, HP, (uint8_t)tstr[i]); + if (!(D0&(1ULL))) { + ++err; if (err>cut) return; + } + + Peq[0] >>= 1; Peq[1] >>= 1; Peq[2] >>= 1; Peq[3] >>= 1; + ++i; ++i_bd; + if(i_bd < pn) { + c = seq_nt4_table[(uint8_t)pstr[i_bd]]; + if(c < 4) Peq[c] |= mm; + } + ez->path.a[ez->path.n++] = D0;//diff + ez->path.a[ez->path.n++] = VP;//diff + ez->path.a[ez->path.n++] = VN;//diff + ez->path.a[ez->path.n++] = HP;//diff + ez->path.a[ez->path.n++] = HN;//diff + } + ed_core_64(Peq, VP, VN, X, D0, HN, HP, (uint8_t)tstr[i]); + if (!(D0&(1ULL))) { + ++err; if (err>cut) return; + } + ez->path.a[ez->path.n++] = D0;//diff + ez->path.a[ez->path.n++] = VP;//diff + ez->path.a[ez->path.n++] = VN;//diff + ez->path.a[ez->path.n++] = HP;//diff + ez->path.a[ez->path.n++] = HN;//diff + + int32_t site = tn - 1 - thre;/**up bound**/\ + if(ez->err > thre) {//diff + for (cut = pn - 1; site < cut; site++) { + err += VP&(1ULL); VP >>= 1; err -= VN&(1ULL); VN >>= 1; + } + + if (site == cut && err <= thre) { + ez->err = err; + ez->pe = pn-1; ez->te = tn-1; + } + } + ///should update ez->path.n for extension + gen_trace(ez, thre, 1);//diff + return; +} + +inline void ed_band_cal_semi_64_w_trace(char *pstr, int32_t pn, char *tstr, int32_t tn, int32_t thre, bit_extz_t *ez) +{ + ez->cigar.n = 0; ez->nword = 1;//diff + if(ez->err > thre) {//diff + init_base_ed(*ez, thre, pn, tn); ez->ps = ez->pe = -1; ez->ts = 0; ez->te = tn-1;//diff + } else if(ez->err == 0) {//diff + push_trace(&(ez->cigar), 0, ez->te+1-ez->ts); //diff + ez->ps = ez->pe - (ez->te-ez->ts);//diff + return;//diff + } + Word c, Peq[5] = {0}, VP = 0, VN = 0, X, D0, HN, HP, mm; + int32_t bd, i, err = 0, i_bd, last_high = (thre<<1), tn0 = tn - 1, cut = thre+last_high; + + bd = (thre<<1)+1; bd = ((bd<=pn)?bd:pn); + for (i = 0, mm = 1; i < bd; i++) { + Peq[seq_nt4_table[(uint8_t)pstr[i]]] |= mm; mm <<= 1; + } + bd = (thre<<1)+1, i_bd = (thre<<1); + + ez->path.n=(ez->nword*tn*5);//diff + kv_resize(w_sig, ez->path, ez->path.n); ez->path.n=0;//diff + + i = 0; Peq[4] = 0; mm = ((Word)1 << (thre<<1));///for the incoming char/last char** + while (i < tn0) { + ed_core_64(Peq, VP, VN, X, D0, HN, HP, (uint8_t)tstr[i]); + if (!(D0&(1ULL))) { + ++err; if (err>cut) return; + } + + Peq[0] >>= 1; Peq[1] >>= 1; Peq[2] >>= 1; Peq[3] >>= 1; + ++i; ++i_bd; c = seq_nt4_table[(uint8_t)pstr[i_bd]]; + if(c < 4) Peq[c] |= mm; + + ez->path.a[ez->path.n++] = D0;//diff + ez->path.a[ez->path.n++] = VP;//diff + ez->path.a[ez->path.n++] = VN;//diff + ez->path.a[ez->path.n++] = HP;//diff + ez->path.a[ez->path.n++] = HN;//diff + } + ed_core_64(Peq, VP, VN, X, D0, HN, HP, (uint8_t)tstr[i]); + if (!(D0&(1ULL))) { + ++err; if (err>cut) return; + } + ez->path.a[ez->path.n++] = D0;//diff + ez->path.a[ez->path.n++] = VP;//diff + ez->path.a[ez->path.n++] = VN;//diff + ez->path.a[ez->path.n++] = HP;//diff + ez->path.a[ez->path.n++] = HN;//diff + + int32_t site = tn - 1;/**up bound**/ + /**in most cases, ai = (thre<<1)**/ + int32_t ai = pn - tn, uge = INT32_MAX; + if(ez->err > thre) {//diff + if ((err <= thre) && (err <= ez->err)) { + ez->err = err; ez->pe = site; + } + + i = 0; + while (i < ai) { + err += ((VP >> i)&(1ULL)); err -= ((VN >> i)&(1ULL)); ++i; + if ((err <= thre) && (err <= ez->err)) { + ez->err = err; ez->pe = site + i; + } + if(i == thre) uge = err; + } + if((uge <= thre) && (uge == ez->err)) ez->pe = site + thre; + } + ///should update ez->path.n for extension + gen_trace(ez, 0, 1);//diff +} + +inline void ed_band_cal_extension_64_0_w_trace(char *pstr, int32_t pn, char *tstr, int32_t tn, int32_t thre, bit_extz_t *ez) +{ + int32_t done = 1;//diff + ez->cigar.n = 0; ez->nword = 1;//diff + if(ez->err > thre) {//diff + init_base_ed(*ez, thre, pn, tn); ez->ps = ez->ts = 0; ez->pe = ez->te = -1; + if(pn > tn + thre) pn = tn + thre; + else if(tn > pn + thre) tn = pn + thre; + done = 0;//diff + } else if(ez->err == 0) {//diff + push_trace(&(ez->cigar), 0, ez->te+1-ez->ts); return;//diff + } else { + pn = ez->pe + 1 - ez->ps; tn = ez->te + 1 - ez->ts; //diff + } + int32_t i, err, tn0 = tn - 1, cut = thre+(thre<<1), bd, i_bd, poff, pe = pn-1, tmp_e = INT32_MAX, k; + Word c, Peq[5] = {0}, VP, VN, X, D0, HN, HP, mm; + + bd = thre+1; bd = ((bd<=pn)?bd:pn); i_bd = thre; + for (i = 0, mm = (((Word)1)<path.n=(ez->nword*tn*5);//diff + kv_resize(w_sig, ez->path, ez->path.n); ez->path.n=0;//diff + + i = 0; mm = ((Word)1 << (thre<<1)); + // X = D0 = HN = HP = mm = 0;//for debug + while (i < tn0) { + // fprintf(stderr, "i::%d\n", i); + // prt_vector(Peq, VP, VN, X, D0, HN, HP, ((thre<<1)+1)); + ed_core_64(Peq, VP, VN, X, D0, HN, HP, (uint8_t)tstr[i]); + if (!(D0&(1ULL))) { + ++err; if (err>cut) return; + } + + if(!done) {//diff + poff = i-thre; k = i+thre-pe;/**poff:[i-thre, i+thre]**/ + if(k >= 0) { + if(k == 0) { + tmp_e = err; + for (k = 0; poff < pe; poff++, (k)++) { + tmp_e += ((VP>>k)&(1ULL)); tmp_e -= ((VN>>k)&(1ULL)); + } + } else { + k = (thre<<1) - k; + if(k >= 0) { + tmp_e += ((HP>>k)&(1ULL)); tmp_e -= ((HN>>k)&(1ULL)); + } + } + if(tmp_e <= (*ez).thre && tmp_e < (*ez).err) { + (*ez).err = tmp_e; (*ez).pe = pe; (*ez).te = i; + } + } + } + // fprintf(stderr, "i::%d, tmp_e::%d, err::%d\n", i, tmp_e, err); + + Peq[0] >>= 1; Peq[1] >>= 1; Peq[2] >>= 1; Peq[3] >>= 1; + ++i; ++i_bd; + if(i_bd < pn) { + c = seq_nt4_table[(uint8_t)pstr[i_bd]]; + if(c < 4) Peq[c] |= mm; + } + + ez->path.a[ez->path.n++] = D0;//diff + ez->path.a[ez->path.n++] = VP;//diff + ez->path.a[ez->path.n++] = VN;//diff + ez->path.a[ez->path.n++] = HP;//diff + ez->path.a[ez->path.n++] = HN;//diff + } + ed_core_64(Peq, VP, VN, X, D0, HN, HP, (uint8_t)tstr[i]); + if (!(D0&(1ULL))) { + ++err; if (err>cut) return; + } + ez->path.a[ez->path.n++] = D0;//diff + ez->path.a[ez->path.n++] = VP;//diff + ez->path.a[ez->path.n++] = VN;//diff + ez->path.a[ez->path.n++] = HP;//diff + ez->path.a[ez->path.n++] = HN;//diff + + int32_t site = tn - 1 - thre;/**up bound; site:[tn - 1 - thre, tn - 1 + thre]**/ + if(!done) {//diff + for (cut = pn - 1; site < cut; ) { + err += VP&(1ULL); VP >>= 1; err -= VN&(1ULL); VN >>= 1; + site++; + if(err <= thre && err < ez->err) { + ez->err = err; ez->pe = site; ez->te = tn-1; + } + } + if(err <= thre && err < ez->err) { + ez->err = err; ez->pe = site; ez->te = tn-1; + } + } + + if((ez->te-ez->ts+1) != tn) {//diff + ez->path.n /= tn; ez->path.n *= (ez->te+1-ez->ts);//diff + }//diff + gen_trace(ez, thre, 1);//diff + return; +} + +inline void ed_band_cal_extension_64_1_w_trace(char *pstr, int32_t pn, char *tstr, int32_t tn, int32_t thre, bit_extz_t *ez) +{ + int32_t done = 1;//diff + ez->cigar.n = 0; ez->nword = 1;//diff + if(ez->err > thre) {//diff + init_base_ed(*ez, thre, pn, tn); ez->ps = ez->ts = INT32_MAX; ez->pe = pn-1; ez->te = tn-1; + if(pn > tn + thre) pn = tn + thre; + else if(tn > pn + thre) tn = pn + thre; + done = 0;//diff + } else if(ez->err == 0) {//diff + push_trace(&(ez->cigar), 0, ez->te+1-ez->ts); return;//diff + } else { + pn = ez->pe + 1 - ez->ps; tn = ez->te + 1 - ez->ts; //diff + } + + int32_t i, err, tn0 = tn - 1, cut = thre+(thre<<1), bd, i_bd; + int32_t poff, pe = pn-1, tmp_e = INT32_MAX, k, pidx = ez->pe, tidx = ez->te; + Word c, Peq[5] = {0}, VP, VN, X, D0, HN, HP, mm; + + bd = thre+1; bd = ((bd<=pn)?bd:pn); i_bd = thre; + for (i = 0, mm = (((Word)1)<path.n=(ez->nword*tn*5);//diff + kv_resize(w_sig, ez->path, ez->path.n); ez->path.n=0;//diff + + i = 0; mm = ((Word)1 << (thre<<1)); + // X = D0 = HN = HP = mm = 0;//for debug + while (i < tn0) { + // fprintf(stderr, "i::%d\n", i); + // prt_vector(Peq, VP, VN, X, D0, HN, HP, ((thre<<1)+1)); + ed_core_64(Peq, VP, VN, X, D0, HN, HP, (uint8_t)tstr[tidx-i]); + if (!(D0&(1ULL))) { + ++err; if (err>cut) return; + } + + if(!done) {//diff + poff = i-thre; k = i+thre-pe;/**poff:[i-thre, i+thre]**/ + if(k >= 0) { + if(k == 0) { + tmp_e = err; + for (k = 0; poff < pe; poff++, (k)++) { + tmp_e += ((VP>>k)&(1ULL)); tmp_e -= ((VN>>k)&(1ULL)); + } + } else { + k = (thre<<1) - k; + if(k >= 0) { + tmp_e += ((HP>>k)&(1ULL)); tmp_e -= ((HN>>k)&(1ULL)); + } + } + if(tmp_e <= (*ez).thre && tmp_e < (*ez).err) { + (*ez).err = tmp_e; (*ez).ps = pidx - pe; (*ez).ts = tidx-i; + } + } + } + // fprintf(stderr, "i::%d, tmp_e::%d, err::%d\n", i, tmp_e, err); + + Peq[0] >>= 1; Peq[1] >>= 1; Peq[2] >>= 1; Peq[3] >>= 1; + ++i; ++i_bd; + if(i_bd < pn) { + c = seq_nt4_table[(uint8_t)pstr[pidx-i_bd]]; + if(c < 4) Peq[c] |= mm; + } + + ez->path.a[ez->path.n++] = D0;//diff + ez->path.a[ez->path.n++] = VP;//diff + ez->path.a[ez->path.n++] = VN;//diff + ez->path.a[ez->path.n++] = HP;//diff + ez->path.a[ez->path.n++] = HN;//diff + } + ed_core_64(Peq, VP, VN, X, D0, HN, HP, (uint8_t)tstr[tidx-i]); + if (!(D0&(1ULL))) { + ++err; if (err>cut) return; + } + ez->path.a[ez->path.n++] = D0;//diff + ez->path.a[ez->path.n++] = VP;//diff + ez->path.a[ez->path.n++] = VN;//diff + ez->path.a[ez->path.n++] = HP;//diff + ez->path.a[ez->path.n++] = HN;//diff + + int32_t site = tn - 1 - thre;/**up bound; site:[tn - 1 - thre, tn - 1 + thre]**/\ + if(!done) {//diff + for (cut = pn - 1; site < cut; ) { + err += VP&(1ULL); VP >>= 1; err -= VN&(1ULL); VN >>= 1; + site++; + if(err <= thre && err < ez->err) { + ez->err = err; ez->ps = pidx-site; ez->ts = tidx+1-tn; + } + } + if(err <= thre && err < ez->err) { + ez->err = err; ez->ps = pidx-site; ez->ts = tidx+1-tn; + } + } + + if((ez->te-ez->ts+1) != tn) {//diff + ez->path.n /= tn; ez->path.n *= (ez->te+1-ez->ts);//diff + }//diff + + poff = ez->ps; ez->ps = pidx - ez->pe; ez->pe = pidx - poff;//diff + gen_trace(ez, thre, 0);//diff + poff = ez->ps; ez->ps = pidx - ez->pe; ez->pe = pidx - poff;//diff + return; +} + + /** pattern is the longer one, while text is the shorter one @@ -2546,22 +4100,6 @@ inline int Reserve_Banded_BPM_4_SSE_only(char *pattern1, char *pattern2, char *p } -#define EAC_M 0 -#define MIS_M 1 -#define MOR_YP 2 -#define MOR_XT 3 - -inline void push_trace(asg16_v *res, uint16_t c, uint32_t len) -{ - uint16_t p; c <<= 14; - while (len >= (0x3fff)) { - p = (c + (0x3fff)); kv_push(uint16_t, *res, p); len -= (0x3fff); - } - if(len) { - p = (c + len); kv_push(uint16_t, *res, p); - } -} - // void move_trace_gap(uint16_t *trace, int32_t trace_n, int32_t trace_i, // char *pstr, int32_t pi, char *tstr, int32_t ti, int32_t *err) // { diff --git a/inter.cpp b/inter.cpp index 0db5bee..80a8e42 100644 --- a/inter.cpp +++ b/inter.cpp @@ -5376,7 +5376,7 @@ int64_t filter_sec(overlap_region_alloc *ol, ul_ov_t *idx, int64_t idx_n, ul_ov_ set_w_e(&(ol->list[a[z].qn]), w_idx, wl, ql); // fprintf(stderr, "[M::%s::utg%.6dl]\n", __func__, (int32_t)a[z].tn+1); } - on_contain += (((idx[k].te-idx[k].ts)==1)?1:0); + on_contain += (((idx[k].te-idx[k].ts)==1)?1:0); } if(on_contain == idx_n) {///each primary chain only has one alignment on_contain = 0; @@ -5413,6 +5413,8 @@ int64_t filter_sec(overlap_region_alloc *ol, ul_ov_t *idx, int64_t idx_n, ul_ov_ } // fprintf(stderr, "-[M::%s] oln::%ld\n", __func__, ol->length); if(alt_occ == 0 || ol->length == 1) return 1;//if all alignments are primary or there is only one alignment + // for (k = ol->length; k < on; k++) ol->list[k].is_match = 2;//reover trans alignments + // ol->length = on; return 0; } @@ -5441,9 +5443,12 @@ int64_t gl_chain_flter(overlap_region_alloc* olist, Correct_dumy* dumy, st_mt_t kv_resize(uint64_t, *sps, idx->n); kv_resize(ul_ov_t, ll->tk, idx->n); - occ = ed_dp_c(olist, idx, ll->tk.a, uref, uopt, G_CHAIN_BW, N_GCHAIN_RATE, ql, 75, dumy->overlapID, ll->srt.a.a, sps->a, 1.25, 0, NULL, uref->ug); + occ = ed_dp_c(olist, idx, ll->tk.a, uref, uopt, G_CHAIN_BW, N_GCHAIN_RATE, ql, 75, dumy->overlapID, ll->srt.a.a, sps->a, ERROR_RATE, 0, NULL, uref->ug); if((!occ) || (!idx->n)) return 0; idx_n = idx->n; p = &(idx->a[idx_n-1]); + if(idx_n <= 1) {//one chain; nothing to do + (*need_phase) = 0; return 1; + } // fprintf(stderr, "[M::%s] qs::%u, qe::%u, ql::%ld, occ::%u\n", __func__, p->qs, p->qe, ql, p->te - p->ts); if(p->qe-p->qs <= (ql*0.25)) return 0;///primary chain is too short i = idx_n-1; occ = p->te - p->ts; @@ -5462,8 +5467,9 @@ int64_t gl_chain_flter(overlap_region_alloc* olist, Correct_dumy* dumy, st_mt_t i++; } - // fprintf(stderr, "[M::%s] i::%ld, idx_n::%ld\n", __func__, i, ((int64_t)idx->n)); - if(occ == (int64_t)olist->length) return 1; + if(occ == (int64_t)olist->length) {//all alignments are primary chains; nothing to do + (*need_phase) = 0; return 1; + } // if(i >= ((int64_t)idx->n)) return 0; nw = get_num_wins(0, ql, wl); kv_resize(uint64_t, ll->srt.a, (uint64_t)nw); if(filter_sec(olist, idx->a+i, idx->n-i, ll->tk.a, ll->srt.a.a, nw, wl, ql)) { diff --git a/main.cpp b/main.cpp index 4030b4b..69e10ff 100644 --- a/main.cpp +++ b/main.cpp @@ -12,11 +12,38 @@ int main(int argc, char *argv[]) yak_reset_realtime(); init_opt(&asm_opt); if (!CommandLine_process(argc, argv, &asm_opt)) return 0; + /** + bit_extz_t exz, exz64; init_bit_extz_t(&exz, 2); init_bit_extz_t(&exz64, 2); + + char *pstr = "GTGTGG", *tsrt = "CTGT"; int32_t thre = 1; + ed_band_cal_semi_infi_w((char*)pstr, strlen(pstr), (char*)tsrt, strlen(tsrt), thre, NULL, &exz); + ed_band_cal_semi_64_w((char*)pstr, strlen(pstr), (char*)tsrt, strlen(tsrt), thre, &exz64); + fprintf(stderr, "\n[M::%s::] exz.err::%d, exz64.err::%d, exz.ps::%d, exz64.ps::%d, exz.pe::%d, exz64.pe::%d, exz.ts::%d, exz64.ts::%d, exz.te::%d, exz64.te::%d\n", __func__, + exz.err, exz64.err, exz.ps, exz64.ps, exz.pe, exz64.pe, exz.ts, exz64.ts, exz.te, exz64.te); + ed_band_cal_semi_64_w_trace((char*)pstr, strlen(pstr), (char*)tsrt, strlen(tsrt), thre, &exz64); + cigar_check((char*)pstr, (char*)tsrt, &(exz64)); + + + // char *pstr = "TGT", *tsrt = "CTGT"; int32_t thre = 1; + // ed_band_cal_global_infi_w((char*)pstr, strlen(pstr), (char*)tsrt, strlen(tsrt), thre, NULL, &exz); + // ed_band_cal_global_64_w((char*)pstr, strlen(pstr), (char*)tsrt, strlen(tsrt), thre, &exz64); + // ed_band_cal_global_64_w_trace((char*)pstr, strlen(pstr), (char*)tsrt, strlen(tsrt), thre, &exz64); + // fprintf(stderr, "\n[M::%s::] exz.err::%d, exz64.err::%d, exz.ps::%d, exz64.ps::%d, exz.pe::%d, exz64.pe::%d, exz.ts::%d, exz64.ts::%d, exz.te::%d, exz64.te::%d\n", __func__, + // exz.err, exz64.err, exz.ps, exz64.ps, exz.pe, exz64.pe, exz.ts, exz64.ts, exz.te, exz64.te); + // cigar_check((char*)pstr, (char*)tsrt, &(exz64)); + + // ed_band_cal_extension_infi0_w((char *)"AAT", 3, (char *)"ACTTTTTT", 8, 2, NULL, &exz); + // ed_band_cal_extension_64_w((char *)"AAT", 3, (char *)"ACTTTTTT", 8, 2, &exz64); + // fprintf(stderr, "\n[M::%s::] exz.err::%d, exz64.err::%d, exz.ps::%d, exz64.ps::%d, exz.pe::%d, exz64.pe::%d, exz.ts::%d, exz64.ts::%d, exz.te::%d, exz64.te::%d\n", __func__, + // exz.err, exz64.err, exz.ps, exz64.ps, exz.pe, exz64.pe, exz.ts, exz64.ts, exz.te, exz64.te); + //bit_extz_t exz; ///ed_band_cal_global_128bit(t_string+r_ts, t_end+1-r_ts, q_string, ql, thres, &exz); // ed_band_cal_extension_128bit((char *)"AAGTTTA", 7, (char *)"CCTTTTTT", 8, 4, &exz); // ed_band_cal_extension_128bit((char *)"AA", 2, (char *)"ACTTTTTT", 8, 1, &exz); // fprintf(stderr, "ed_extension::%d, pe::%d, te::%d\n", exz.err, exz.pe, exz.te); - // exit(1); + exit(1); + **/ + // fprintf(stderr, "[M::%s::] ed_global::%d, ed_global_128bit::%d\n", __func__, // ed_band_cal_global((char *)"ACT", 3, (char *)"AAT", 3, 1),