29 using namespace shogun;
51 -1, 0, 0, 2, -5, 3, 7,
52 0, -3, -1, -2, -4, -3, -3, 8,
53 -2, 0, 1, -2, -4, 1, 0, -3, 11,
54 -2, -5, -5, -5, -2, -4, -5, -6, -5, 6,
55 -2, -3, -5, -5, -2, -3, -4, -5, -4, 2, 6,
56 -1, 3, 0, -1, -5, 2, 1, -2, -1, -4, -4, 7,
57 -1, -2, -3, -5, -2, -1, -3, -4, -2, 2, 3, -2, 8,
58 -3, -4, -5, -5, -4, -5, -5, -5, -2, 0, 1, -5, 0, 9,
59 -1, -3, -3, -2, -4, -2, -2, -3, -3, -4, -4, -2, -4, -5, 11,
60 2, -1, 1, 0, -1, 0, 0, 0, -1, -4, -4, 0, -2, -4, -1, 6,
61 0, -2, 0, -2, -1, -1, -1, -2, -3, -1, -2, -1, -1, -3, -2, 2, 7,
62 -4, -4, -6, -6, -3, -3, -4, -4, -4, -4, -2, -4, -2, 1, -6, -4, -4, 16,
63 -3, -3, -3, -5, -4, -2, -3, -5, 3, -2, -2, -3, -1, 4, -4, -3, -2, 3, 10,
64 0, -4, -4, -5, -1, -3, -4, -5, -5, 4, 1, -3, 1, -1, -4, -2, 0, -4, -2, 6};
67 #define BINDEX(i,j) (((i)>(j))?(j)+(((i)*(i+1))/2):(i)+(((j)*(j+1))/2))
81 #define LOGP(x,y) LogSum(x,y)
90 #define INTSCALE 1000.0
98 init_static_variables();
109 init_static_variables();
141 void CLocalAlignmentStringKernel::init_logsum(){
148 int32_t CLocalAlignmentStringKernel::LogSum(int32_t p1, int32_t p2)
151 static int32_t firsttime=1;
160 if (diff>=LOGSUM_TBL)
return p1;
161 else if (diff<=-LOGSUM_TBL)
return p2;
163 else return p2+logsum_lookup[-diff];
170 return (p1-p2 > 50.) ? p1 : p1 + log(1. + exp(p2-p1));
172 return (p2-p1 > 50.) ? p2 : p2 + log(1. + exp(p1-p2));
176 void CLocalAlignmentStringKernel::init_static_variables()
182 if ((
aaIndex=(int32_t *)calloc(
NLET,
sizeof(int32_t))) == NULL)
188 if ((
isAA=(int32_t *)calloc(256,
sizeof(int32_t))) == NULL)
194 for (i=0 ; i<NAA*(NAA+1)/2; i++)
208 float64_t CLocalAlignmentStringKernel::LAkernelcompute(
209 int32_t* aaX, int32_t* aaY,
210 int32_t nX, int32_t nY )
270 for (i=1;i<=nX;i++) {
277 logX2[curpos] =
LOG0;
278 logY2[curpos] =
LOG0;
281 for (j=1;j<=nY;j++) {
288 frompos = old*cl + j;
295 logX2[curpos] =
LOGP( logM[frompos] , logX2[frompos] );
300 frompos = cur*cl + j-1;
307 aux =
LOGP( logM[frompos] , logY2[frompos] );
308 logY2[curpos] =
LOGP( aux , logX2[frompos] );
313 frompos = old*cl + j-1;
315 aux =
LOGP( logX[frompos] , logY[frompos] );
316 aux2 =
LOGP( 0 , logM[frompos] );
335 curpos = old*cl + nY;
336 aux =
LOGP( logX2[curpos] , logY2[curpos] );
337 aux2 =
LOGP( 0 , logM[curpos] );
369 if ((lx<1) || (ly<1))
374 if ((aax=(int32_t *)calloc(lx,
sizeof(int32_t))) == NULL)
376 if ((aay=(int32_t *)calloc(ly,
sizeof(int32_t))) == NULL)
382 for (i=0 ; i<lx ; i++)
383 if (
isAA[toupper(x[i])])
384 aax[j++] =
aaIndex[toupper(x[i])-
'A'];
387 for (i=0 ; i<ly ; i++)
388 if (
isAA[toupper(y[i])])
389 aay[j++] =
aaIndex[toupper(y[i])-
'A'];
394 float64_t result=LAkernelcompute(aax,aay,lx,ly);
406 void CLocalAlignmentStringKernel::init()