Permalink
Browse files

This branch keeps a backup of xmatch-like strategy

  • Loading branch information...
1 parent 8c13339 commit f60f222606a092d2a8dcb7fb64a82bf74c397c9e Heng Li committed Sep 11, 2011
Showing with 20 additions and 7 deletions.
  1. +20 −7 ksw.c
View
27 ksw.c
@@ -92,9 +92,9 @@ ksw_query_t *ksw_qinit(int size, int qlen, const uint8_t *query, int m, const in
int ksw_sse2_16(ksw_query_t *q, int tlen, const uint8_t *target, ksw_aux_t *a) // the first gap costs -(_o+_e)
{
- int slen, i, m_b, n_b, te = -1, gmax = 0;
+ int sum = 0, slen, i, m_b, n_b, te = -1, gmax = 0, _thres = 254 - q->mdiff;
uint64_t *b;
- __m128i zero, gapoe, gape, shift, *H0, *H1, *E, *Hmax;
+ __m128i zero, gapoe, gape, shift, *H0, *H1, *E, *Hmax, thres, reduce;
#define __max_16(ret, xx) do { \
(xx) = _mm_max_epu8((xx), _mm_srli_si128((xx), 8)); \
@@ -110,6 +110,8 @@ int ksw_sse2_16(ksw_query_t *q, int tlen, const uint8_t *target, ksw_aux_t *a) /
gapoe = _mm_set1_epi8(a->gapo + a->gape);
gape = _mm_set1_epi8(a->gape);
shift = _mm_set1_epi8(q->shift);
+ thres = _mm_set1_epi8(_thres);
+ reduce = _mm_set1_epi8(127);
H0 = q->H0; H1 = q->H1; E = q->E; Hmax = q->Hmax;
slen = q->slen;
for (i = 0; i < slen; ++i) {
@@ -164,24 +166,35 @@ int ksw_sse2_16(ksw_query_t *q, int tlen, const uint8_t *target, ksw_aux_t *a) /
end_loop16:
//int k;for (k=0;k<16;++k)printf("%d ", ((uint8_t*)&max)[k]);printf("\n");
__max_16(imax, max); // imax is the maximum number in max
- if (imax >= a->T) { // write the b array; this condition adds branching unfornately
+ if (imax + sum >= a->T) { // write the b array; this condition adds branching unfornately
if (n_b == 0 || (int32_t)b[n_b-1] + 1 != i) { // then append
if (n_b == m_b) {
m_b = m_b? m_b<<1 : 8;
b = realloc(b, 8 * m_b);
}
- b[n_b++] = (uint64_t)imax<<32 | i;
- } else if ((int)(b[n_b-1]>>32) < imax) b[n_b-1] = (uint64_t)imax<<32 | i; // modify the last
+ b[n_b++] = (uint64_t)(imax + sum)<<32 | i;
+ } else if ((int)(b[n_b-1]>>32) < imax + sum) b[n_b-1] = (uint64_t)(imax + sum)<<32 | i; // modify the last
}
if (imax > gmax) {
gmax = imax; te = i; // te is the end position on the target
for (j = 0; LIKELY(j < slen); ++j) // keep the H1 vector
_mm_store_si128(Hmax + j, _mm_load_si128(H1 + j));
if (gmax + q->shift >= 255) break;
}
+ if (gmax >= _thres) { // reduce
+ // When overflow is going to happen, subtract 127 from all scores. This heuristic
+ // may miss the best alignment, but in practice this should happen very rarely.
+ sum += 127; gmax -= 127;
+ for (j = 0; LIKELY(j < slen); ++j) {
+ h = _mm_subs_epu8(_mm_load_si128(H1 + j), reduce);
+ _mm_store_si128(H1 + j, h);
+ e = _mm_subs_epu8(_mm_load_si128(E + j), reduce);
+ _mm_store_si128(E + j, e);
+ }
+ }
S = H1; H1 = H0; H0 = S; // swap H0 and H1
}
- a->score = gmax; a->te = te;
+ a->score = gmax + sum; a->te = te;
{ // get a->qe, the end of query match; find the 2nd best score
int max = -1, low, high, qlen = slen * 16;
uint8_t *t = (uint8_t*)Hmax;
@@ -197,7 +210,7 @@ int ksw_sse2_16(ksw_query_t *q, int tlen, const uint8_t *target, ksw_aux_t *a) /
}
}
free(b);
- return a->score + q->shift >= 255? 255 : a->score;
+ return a->score;
}
int ksw_sse2_8(ksw_query_t *q, int tlen, const uint8_t *target, ksw_aux_t *a) // the first gap costs -(_o+_e)

0 comments on commit f60f222

Please sign in to comment.