diff --git a/Correct.cpp b/Correct.cpp index 6107f88..1e36199 100644 --- a/Correct.cpp +++ b/Correct.cpp @@ -3657,16 +3657,19 @@ char *qstr, char *tstr, char *tstr1, Correct_dumy* dumy, uint32_t rev, uint32_t // fprintf(stderr, "[qstr] %.*s\n", (int32_t)ql, q_string); // } // assert(dbg_e <= (int32_t)error); - // bit_extz_t exz; ///ed_band_cal_global_128bit(t_string+r_ts, t_end+1-r_ts, q_string, ql, thres, &exz); + // bit_extz_t exz; // ed_band_cal_extension_128bit(t_string+r_ts, t_end+1-r_ts, q_string, ql, thres, &exz); + // assert(exz.err <= (int32_t)error); + + // ed_band_cal_global_128bit(t_string+r_ts, t_end+1-r_ts, q_string, ql, thres, &exz); + // assert(exz.err <= (int32_t)error); // if(exz.err > (int32_t)error && ql == 1) { // fprintf(stderr, "[M::%s::] error::%u, ed_extension::%d, ql::%ld, thres::%ld\n", // __func__, error, exz.err, ql, thres); // 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(exz.err <= (int32_t)error); - + // 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)); diff --git a/Levenshtein_distance.h b/Levenshtein_distance.h index 4968266..ba5a0fd 100644 --- a/Levenshtein_distance.h +++ b/Levenshtein_distance.h @@ -612,6 +612,11 @@ typedef struct { // asg16_v cigar; w64_trace_t path; } bit_extz_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;\ +} + #define w128_bit(x, b) ((x).a[((b)>>bitw)]|=(((w_sig)1)<<((b)&bitz))) #define w128_get_bit(x, b) (((x).a[((b)>>bitw)]>>((b)&bitz))&((w_sig)1)) @@ -650,7 +655,9 @@ typedef struct { else (x).a[1] = (((w_sig)1)<<((l)-bitwbit))-1;\ } while (0) \ -#define ed_core_w128(Peq, VP, VN, X, D0, HN, HP) { \ + + +#define ed_core_w128b(Peq, VP, VN, X, D0, HN, HP) { \ /**X = Peq[seq_nt4_table[(uint8_t)tstr[i]]] | VN;**/\ c = seq_nt4_table[(uint8_t)tstr[i]]; w128_or(X, Peq[c], VN);\ /**D0 = ((VP + (X&VP)) ^ VP) | X;**/\ @@ -672,228 +679,227 @@ typedef struct { w128_or(VP, X, HP);\ w128_self_not(VP);\ w128_self_or(VP, HN);\ - } +} -#define ed_core_w128_reshift(Peq) {\ +#define ed_core(sf, Peq, VP, VN, X, D0, HN, HP) { \ + /**X = Peq[seq_nt4_table[(uint8_t)tstr[i]]] | VN;**/\ + c = seq_nt4_table[(uint8_t)tstr[i]]; ##sf##or(X, Peq[c], VN);\ + /**D0 = ((VP + (X&VP)) ^ VP) | X;**/\ + ##sf##and(D0, X, VP);\ + ##sf##self_add(D0, VP);\ + ##sf##self_xor(D0, VP);\ + ##sf##self_or(D0, X);\ + /**HN = VP&D0;**/\ + ##sf##and(HN, VP, D0);\ + /**HP = VN | ~(VP | D0);**/\ + ##sf##or(HP, VP, D0);\ + ##sf##self_not(HP);\ + ##sf##self_or(HP, VN);\ + /**X = D0 >> 1;**/\ + X = D0; ##sf##self_rsft_1(X);\ + /**VN = X&HP;**/\ + ##sf##and(VN, X, HP);\ + /**VP = HN | ~(X | HP);**/\ + ##sf##or(VP, X, HP);\ + ##sf##self_not(VP);\ + ##sf##self_or(VP, HN);\ +} + + +#define HA_ED_INIT(sf)\ +inline void ed_band_cal_global_##sf##bit(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;\ + 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 = thre+1, i_bd = thre; uint8_t c;\ + 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->mm); w##sf##_bit(ez->mm, thre); /**mm = (((Word)1)<Peq[seq_nt4_table[(uint8_t)pstr[i]]], ez->mm); w##sf##_self_lsft_1(ez->mm);\ + /** Peq[seq_nt4_table[(uint8_t)pstr[i]]] |= mm; mm <<= 1;**/\ + }\ + 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; \ + /**for the incoming char/last char**/\ + w##sf##_clear(ez->mm); w##sf##_bit(ez->mm, (thre<<1)); /**mm = ((Word)1 << (thre<<1));**/\ + while (i < tn0) {\ + ed_core_w##sf##b(ez->Peq, ez->VP, ez->VN, ez->X, ez->D0, ez->HN, ez->HP);\ + if (!(ez->D0.a[0]&(1ULL))) {\ + ++err; if (err>cut) return;\ + }\ /** Peq[0] >>= 1; Peq[1] >>= 1; Peq[2] >>= 1; Peq[3] >>= 1;**/\ - w128_self_rsft_1(Peq[0]); w128_self_rsft_1(Peq[1]);\ - w128_self_rsft_1(Peq[2]); w128_self_rsft_1(Peq[3]);\ -} - -#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;\ -} - -#define tst_band_err(p_off, t_off, p_end, dif, ez, k) {\ - if((p_off) + (((ez).thre)<<1) >= (p_end)) { \ - for ((k) = 0; (p_off) < (p_end); (p_off)++) {\ - (dif) += w128_get_bit((ez).VP, (k));\ - (dif) -= w128_get_bit((ez).VN, (k));\ - (k)++;}\ - if((dif) <= (ez).thre && (dif) < (ez).err) {\ - (ez).err = dif; (ez).pe = p_off; (ez).te = t_off;}\ + 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);\ + }\ + }\ + ed_core_w##sf##b(ez->Peq, ez->VP, ez->VN, ez->X, ez->D0, ez->HN, ez->HP);\ + if (!(ez->D0.a[0]&(1ULL))) {\ + ++err; if (err>cut) return;\ }\ + \ + 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); \ + }\ + \ + if (site == cut && err <= thre) {\ + ez->err = err; \ + ez->pe = pn-1; ez->te = tn-1;\ + }\ + return;\ +}\ +inline void ed_band_cal_semi_##sf##bit(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->mm); w##sf##_bit(ez->mm, 0);\ + 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 = (thre<<1)+1, i, err = 0, i_bd = (thre<<1), last_high = (thre<<1), tn0 = tn - 1;\ + int32_t cut = thre+last_high; uint8_t c;\ + \ + for (i = 0; i < bd; i++) {\ + w##sf##_self_or(ez->Peq[seq_nt4_table[(uint8_t)pstr[i]]], ez->mm); w##sf##_self_lsft_1(ez->mm);\ + /** Peq[seq_nt4_table[(uint8_t)pstr[i]]] |= mm; mm <<= 1;**/\ + }\ + /**should make Peq[4] = 0 if N is always an error**/\ + /** 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##sf##_bit(ez->mm, (thre<<1));\ + \ + i = 0;\ + while (i < tn0) {\ + ed_core_w##sf##b(ez->Peq, ez->VP, ez->VN, ez->X, ez->D0, ez->HN, ez->HP);\ + 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; \ + /** 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);\ + }\ + ed_core_w##sf##b(ez->Peq, ez->VP, ez->VN, ez->X, ez->D0, ez->HN, ez->HP);\ + if (!(ez->D0.a[0]&(1ULL))) {\ + ++err; if (err>cut) return;\ + }\ + \ + int32_t site = tn - 1;/**up bound**/\ + /**in most cases, ai = (thre<<1)**/\ + int32_t ai = pn - tn, uge = INT32_MAX;\ + if ((err <= thre) && (err <= ez->err)) {\ + ez->err = err; ez->pe = site;\ + }\ + 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); \ + ++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;\ +}\ +/**require:: (pn >= tn - thre && pn <= tn + thre)**/\ +inline void ed_band_cal_extension_##sf##bit(char *pstr, int32_t pn, char *tstr, int32_t tn, int32_t thre, bit_extz_t *ez)\ +{\ + 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 = thre+1, i_bd = thre; uint8_t 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]);\ + \ + w##sf##_clear(ez->mm); w##sf##_bit(ez->mm, thre); /**mm = (((Word)1)<Peq[seq_nt4_table[(uint8_t)pstr[i]]], ez->mm); w##sf##_self_lsft_1(ez->mm);\ + /**Peq[seq_nt4_table[(uint8_t)pstr[i]]] |= mm; mm <<= 1;**/\ + }\ + 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**/\ + w##sf##_clear(ez->mm); w##sf##_bit(ez->mm, (thre<<1)); /**mm = ((Word)1 << (thre<<1));**/\ + while (i < tn0) {\ + ed_core_w##sf##b(ez->Peq, ez->VP, ez->VN, ez->X, ez->D0, ez->HN, ez->HP);\ + 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##sf##_get_bit((*ez).VP, (k));\ + (tmp_e) -= w##sf##_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;}\ + }\ + }\ + \ + /** 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);\ + }\ + }\ + ed_core_w##sf##b(ez->Peq, ez->VP, ez->VN, ez->X, ez->D0, ez->HN, ez->HP);\ + 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); \ + 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;\ + }\ + return;\ } - -inline void ed_band_cal_global_128bit(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; - 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 = thre+1, i_bd = thre; uint8_t c; - w128_clear(ez->Peq[0]); w128_clear(ez->Peq[1]); w128_clear(ez->Peq[2]); w128_clear(ez->Peq[3]); w128_clear(ez->Peq[4]); - - w128_clear(ez->mm); w128_bit(ez->mm, thre); ///mm = (((Word)1)<Peq[seq_nt4_table[(uint8_t)pstr[i]]], ez->mm); w128_self_lsft_1(ez->mm); - // Peq[seq_nt4_table[(uint8_t)pstr[i]]] |= mm; mm <<= 1; - } - w128_clear(ez->Peq[4]); - err = thre; - w128_set_bit_lsub(ez->VN, thre); ///VN = (((Word)1)<<(thre))-1; - w128_set_bit_lsub(ez->VP, (thre<<1)+1); ///VP = (((Word)1)<<((thre<<1)+1))-1; - w128_self_xor(ez->VP, ez->VN); ///VP ^= VN; - - // print_bits(Peq[0].a, (thre<<1)+1, "-Peq[A]"); - // print_bits(Peq[1].a, (thre<<1)+1, "-Peq[C]"); - // print_bits(Peq[2].a, (thre<<1)+1, "-Peq[G]"); - // print_bits(Peq[3].a, (thre<<1)+1, "-Peq[T]"); - // print_bits(VN.a, (thre<<1)+1, "-VN"); - // print_bits(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 - w128_clear(ez->mm); w128_bit(ez->mm, (thre<<1)); ///mm = ((Word)1 << (thre<<1)); - //VP + ((Peq|VN)&VP) - while (i < tn0) { - ed_core_w128(ez->Peq, ez->VP, ez->VN, ez->X, ez->D0, ez->HN, ez->HP); - if (!(ez->D0.a[0]&(1ULL))) { - ++err; if (err>cut) return; - } - ed_core_w128_reshift(ez->Peq); - ++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) w128_self_or(ez->Peq[c], ez->mm); - } - } - ed_core_w128(ez->Peq, ez->VP, ez->VN, ez->X, ez->D0, ez->HN, ez->HP); - if (!(ez->D0.a[0]&(1ULL))) { - ++err; if (err>cut) return; - } - - 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); w128_self_rsft_1(ez->VP); - // err -= ((VN >> i)&(1ULL)); - err -= ez->VN.a[0]&(1ULL); w128_self_rsft_1(ez->VN); - } - - if (site == cut && err <= thre) { - ez->err = err; - ez->pe = pn-1; ez->te = tn-1; - } - return; -} - -inline void ed_band_cal_semi_128bit(char *pstr, int32_t pn, char *tstr, int32_t tn, int32_t thre, bit_extz_t *ez) -{ - // (*re_err) = INT32_MAX; - init_base_ed(*ez, thre, pn, tn); ez->ps = ez->pe = -1; ez->ts = 0; ez->te = tn-1; - w128_clear(ez->VP); w128_clear(ez->VN); w128_clear(ez->mm); w128_bit(ez->mm, 0); - w128_clear(ez->Peq[0]); w128_clear(ez->Peq[1]); w128_clear(ez->Peq[2]); w128_clear(ez->Peq[3]); w128_clear(ez->Peq[4]); - int32_t bd = (thre<<1)+1, i, err = 0, i_bd = (thre<<1), last_high = (thre<<1), tn0 = tn - 1; - int32_t cut = thre+last_high; uint8_t c; - - for (i = 0; i < bd; i++) { - w128_self_or(ez->Peq[seq_nt4_table[(uint8_t)pstr[i]]], ez->mm); w128_self_lsft_1(ez->mm); - // Peq[seq_nt4_table[(uint8_t)pstr[i]]] |= mm; mm <<= 1; - } - ///should make Peq[4] = 0 if N is always an error - // Peq[4] = 0; - w128_clear(ez->Peq[4]); - //mm = ((Word)1 << (thre<<1));///for the incoming char/last char - w128_clear(ez->mm); w128_bit(ez->mm, (thre<<1)); - - i = 0; - while (i < tn0) { - ed_core_w128(ez->Peq, ez->VP, ez->VN, ez->X, ez->D0, ez->HN, ez->HP); - if (!(ez->D0.a[0]&(1ULL))) { - ++err; if (err>cut) return; - } - ed_core_w128_reshift(ez->Peq); - - ++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) w128_self_or(ez->Peq[c], ez->mm); - } - ed_core_w128(ez->Peq, ez->VP, ez->VN, ez->X, ez->D0, ez->HN, ez->HP); - if (!(ez->D0.a[0]&(1ULL))) { - ++err; if (err>cut) return; - } - - int32_t site = tn - 1;///up bound - ///in most cases, ai = (thre<<1) - int32_t ai = pn - tn, uge = INT32_MAX; - if ((err <= thre) && (err <= ez->err)) { - ez->err = err; ez->pe = site; - } - i = 0; - - while (i < ai) { - // err += ((VP >> i)&(1ULL)); - err += ez->VP.a[0]&(1ULL); w128_self_rsft_1(ez->VP); - // err -= ((VN >> i)&(1ULL)); - err -= ez->VN.a[0]&(1ULL); w128_self_rsft_1(ez->VN); - ++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; -} - -///require:: (pn >= tn - thre && pn <= tn + thre) -inline void ed_band_cal_extension_128bit(char *pstr, int32_t pn, char *tstr, int32_t tn, int32_t thre, bit_extz_t *ez) -{ - // fprintf(stderr, "\n[M::%s::] pn::%d, tn::%d, thre::%d\n", __func__, pn, tn, thre); - if(pn > tn + thre) pn = tn + thre; - else if(tn > pn + thre) tn = pn + thre; - // fprintf(stderr, "[M::%s::] pn::%d, tn::%d, thre::%d\n", __func__, pn, tn, thre); - init_base_ed(*ez, thre, pn, tn); ez->ps = ez->ts = 0; ez->pe = ez->te = -1; - // 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 = thre+1, i_bd = thre; uint8_t c; - int32_t poff, pe = pn-1, tmp_e, k; - w128_clear(ez->Peq[0]); w128_clear(ez->Peq[1]); w128_clear(ez->Peq[2]); w128_clear(ez->Peq[3]); w128_clear(ez->Peq[4]); - - w128_clear(ez->mm); w128_bit(ez->mm, thre); ///mm = (((Word)1)<Peq[seq_nt4_table[(uint8_t)pstr[i]]], ez->mm); w128_self_lsft_1(ez->mm); - // Peq[seq_nt4_table[(uint8_t)pstr[i]]] |= mm; mm <<= 1; - } - w128_clear(ez->Peq[4]); - err = thre; - w128_set_bit_lsub(ez->VN, thre); ///VN = (((Word)1)<<(thre))-1; - w128_set_bit_lsub(ez->VP, (thre<<1)+1); ///VP = (((Word)1)<<((thre<<1)+1))-1; - w128_self_xor(ez->VP, ez->VN); ///VP ^= VN; - - // print_bits(ez->Peq[0].a, (thre<<1)+1, "-Peq[A]"); - // print_bits(ez->Peq[1].a, (thre<<1)+1, "-Peq[C]"); - // print_bits(ez->Peq[2].a, (thre<<1)+1, "-Peq[G]"); - // print_bits(ez->Peq[3].a, (thre<<1)+1, "-Peq[T]"); - // print_bits(ez->VN.a, (thre<<1)+1, "-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 - w128_clear(ez->mm); w128_bit(ez->mm, (thre<<1)); ///mm = ((Word)1 << (thre<<1)); - //VP + ((Peq|VN)&VP) - while (i < tn0) { - ed_core_w128(ez->Peq, ez->VP, ez->VN, ez->X, ez->D0, ez->HN, ez->HP); - if (!(ez->D0.a[0]&(1ULL))) { - ++err; if (err>cut) return; - } - poff = i-thre; tmp_e = err; //poff:[i-thre, i+thre] - tst_band_err(poff, i, pe, tmp_e, *ez, k); - - ed_core_w128_reshift(ez->Peq); - ++i; ++i_bd; - if(i_bd < pn) { - c = seq_nt4_table[(uint8_t)pstr[i_bd]]; - if(c < 4) w128_self_or(ez->Peq[c], ez->mm); - } - } - ed_core_w128(ez->Peq, ez->VP, ez->VN, ez->X, ez->D0, ez->HN, ez->HP); - 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); w128_self_rsft_1(ez->VP); - // err -= ((VN >> i)&(1ULL)); - err -= ez->VN.a[0]&(1ULL); w128_self_rsft_1(ez->VN); - 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; - } - return; -} +HA_ED_INIT(128) /** diff --git a/main.cpp b/main.cpp index 621d4dc..4030b4b 100644 --- a/main.cpp +++ b/main.cpp @@ -12,11 +12,11 @@ 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; ///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); + //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); + // fprintf(stderr, "ed_extension::%d, pe::%d, te::%d\n", exz.err, exz.pe, exz.te); + // 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),