00001 #include "lib/common.h"
00002 #include "lib/io.h"
00003 #include "lib/Trie.h"
00004 #include "lib/Mathematics.h"
00005
00006 template <>
00007 void CTrie<POIMTrie>::POIMs_extract_W_helper( const INT nodeIdx, const int depth, const INT offset, const INT y0,
00008 DREAL* const* const W, const INT K )
00009 {
00010 ASSERT(nodeIdx!=NO_CHILD);
00011 ASSERT(depth<K);
00012 DREAL* const W_kiy = & W[ depth ][ offset + y0 ];
00013 POIMTrie* const node = &TreeMem[ nodeIdx ];
00014 INT sym;
00015
00016 if( depth < degree-1 )
00017 {
00018 const INT offset1 = offset * NUM_SYMS;
00019 for( sym = 0; sym < NUM_SYMS; ++sym )
00020 {
00021 ASSERT(W_kiy[sym]==0);
00022 const INT childIdx = node->children[ sym ];
00023 if( childIdx != NO_CHILD )
00024 {
00025 W_kiy[ sym ] = TreeMem[ childIdx ].weight;
00026
00027 if (depth < K-1)
00028 {
00029 const INT y1 = ( y0 + sym ) * NUM_SYMS;
00030 POIMs_extract_W_helper( childIdx, depth+1, offset1, y1, W, K );
00031 }
00032 }
00033 }
00034 }
00035 else
00036 {
00037 ASSERT(depth==degree-1);
00038 for( sym = 0; sym < NUM_SYMS; ++sym )
00039 {
00040 ASSERT(W_kiy[sym]==0);
00041 W_kiy[ sym ] = node->child_weights[ sym ];
00042 }
00043 }
00044 }
00045
00046 template <>
00047 void CTrie<POIMTrie>::POIMs_extract_W( DREAL* const* const W, const INT K )
00048 {
00049 ASSERT(degree>=1);
00050 ASSERT(K>=1);
00051 const INT N = length;
00052 INT i;
00053 for( i = 0; i < N; ++i ) {
00054
00055 POIMs_extract_W_helper( trees[i], 0, i*NUM_SYMS, 0*NUM_SYMS, W, K );
00056 }
00057 }
00058
00059 template <>
00060 void CTrie<POIMTrie>::POIMs_calc_SLR_helper1( const DREAL* const distrib, const INT i,
00061 const INT nodeIdx, INT left_tries_idx[4], const int depth, INT const lastSym,
00062 DREAL* S, DREAL* L, DREAL* R )
00063 {
00064 ASSERT(depth==degree-1);
00065 ASSERT(nodeIdx!=NO_CHILD);
00066
00067 const DREAL* const distribLeft = & distrib[ (i-1) * NUM_SYMS ];
00068 POIMTrie* const node = &TreeMem[ nodeIdx ];
00069 INT symRight;
00070 INT symLeft;
00071
00072
00073 node->S = 0;
00074 node->L = 0;
00075 node->R = 0;
00076
00077 if (i+depth<length)
00078 {
00079 const DREAL* const distribRight = & distrib[ (i+depth) * NUM_SYMS ];
00080
00081
00082 for( symRight = 0; symRight < NUM_SYMS; ++symRight ) {
00083 const DREAL w1 = node->child_weights[ symRight ];
00084 const DREAL pRight = distribRight[ symRight ];
00085 const DREAL incr1 = pRight * w1;
00086 node->S += incr1;
00087 node->R += incr1;
00088 }
00089 }
00090
00091
00092 for( symLeft = 0; symLeft < NUM_SYMS; ++symLeft )
00093 {
00094 if (left_tries_idx[symLeft] != NO_CHILD)
00095 {
00096 POIMTrie* nodeLeft = &TreeMem[left_tries_idx[symLeft]];
00097
00098 ASSERT(nodeLeft);
00099 const DREAL w2 = nodeLeft->child_weights[ lastSym ];
00100 const DREAL pLeft = distribLeft[ symLeft ];
00101 const DREAL incr2 = pLeft * w2;
00102 node->S += incr2;
00103 node->L += incr2;
00104 }
00105 }
00106
00107
00108 const DREAL w0 = node->weight;
00109 node->S += w0;
00110 node->L += w0;
00111 node->R += w0;
00112 *S = node->S;
00113 *L = node->L;
00114 *R = node->R;
00115 }
00116
00117
00118 template <>
00119 void CTrie<POIMTrie>::POIMs_calc_SLR_helper2( const DREAL* const distrib, const INT i,
00120 const INT nodeIdx, INT left_tries_idx[4], const int depth,
00121 DREAL* S, DREAL* L, DREAL* R )
00122 {
00123 ASSERT(0<=depth && depth<=degree-2);
00124 ASSERT(nodeIdx!=NO_CHILD);
00125
00126 const DREAL* const distribLeft = & distrib[ (i-1) * NUM_SYMS ];
00127 POIMTrie* const node = &TreeMem[ nodeIdx ];
00128 DREAL dummy;
00129 DREAL auxS;
00130 DREAL auxR;
00131 INT symRight;
00132 INT symLeft;
00133
00134
00135 node->S = 0;
00136 node->L = 0;
00137 node->R = 0;
00138
00139
00140 for( symRight = 0; symRight < NUM_SYMS; ++symRight )
00141 {
00142 const INT childIdx = node->children[ symRight ];
00143 if( childIdx != NO_CHILD )
00144 {
00145 if( depth < degree-2 )
00146 {
00147 INT new_left_tries_idx[4];
00148
00149 for( symLeft = 0; symLeft < NUM_SYMS; ++symLeft )
00150 {
00151 new_left_tries_idx[symLeft]=NO_CHILD;
00152
00153 if (left_tries_idx[symLeft] != NO_CHILD)
00154 {
00155 POIMTrie* nodeLeft = &TreeMem[left_tries_idx[symLeft]];
00156 ASSERT(nodeLeft);
00157 new_left_tries_idx[ symLeft ]=nodeLeft->children[ symRight ];
00158 }
00159 }
00160 POIMs_calc_SLR_helper2( distrib, i, childIdx, new_left_tries_idx, depth+1, &auxS, &dummy, &auxR );
00161 }
00162 else
00163 POIMs_calc_SLR_helper1( distrib, i, childIdx, left_tries_idx, depth+1, symRight, &auxS, &dummy, &auxR );
00164
00165 if (i+depth<length)
00166 {
00167 const DREAL* const distribRight = & distrib[ (i+depth) * NUM_SYMS ];
00168 const DREAL pRight = distribRight[ symRight ];
00169 node->S += pRight * auxS;
00170 node->R += pRight * auxR;
00171 }
00172 }
00173 }
00174
00175
00176 for( symLeft = 0; symLeft < NUM_SYMS; ++symLeft )
00177 {
00178 if (left_tries_idx[symLeft] != NO_CHILD)
00179 {
00180 const POIMTrie* nodeLeft = &TreeMem[left_tries_idx[symLeft]];
00181 ASSERT(nodeLeft);
00182 const DREAL pLeft = distribLeft[ symLeft ];
00183
00184 node->S += pLeft * nodeLeft->S;
00185 node->L += pLeft * nodeLeft->L;
00186
00187 if (i+depth<length)
00188 {
00189 const DREAL* const distribRight = & distrib[ (i+depth) * NUM_SYMS ];
00190
00191 auxS = 0;
00192 if( depth < degree-2 )
00193 {
00194 for( symRight = 0; symRight < NUM_SYMS; ++symRight )
00195 {
00196 const INT childIdxLeft = nodeLeft->children[ symRight ];
00197 if( childIdxLeft != NO_CHILD )
00198 {
00199 const POIMTrie* const childLeft = &TreeMem[ childIdxLeft ];
00200 auxS += distribRight[symRight] * childLeft->S;
00201 }
00202 }
00203 }
00204 else
00205 {
00206 for( symRight = 0; symRight < NUM_SYMS; ++symRight ) {
00207 auxS += distribRight[symRight] * nodeLeft->child_weights[ symRight ];
00208 }
00209 }
00210 node->S -= pLeft* auxS;
00211 }
00212 }
00213 }
00214
00215
00216 const DREAL w0 = node->weight;
00217
00218 node->S += w0;
00219 node->L += w0;
00220 node->R += w0;
00221 *S = node->S;
00222 *L = node->L;
00223 *R = node->R;
00224 }
00225
00226
00227
00228 template <>
00229 void CTrie<POIMTrie>::POIMs_precalc_SLR( const DREAL* const distrib )
00230 {
00231 if( degree == 1 ) {
00232 return;
00233 }
00234
00235 ASSERT(degree>=2);
00236 const INT N = length;
00237 DREAL dummy;
00238 INT symLeft;
00239 INT leftSubtrees[4];
00240 for( symLeft = 0; symLeft < NUM_SYMS; ++symLeft )
00241 leftSubtrees[ symLeft ] = NO_CHILD;
00242
00243 for(INT i = 0; i < N; ++i )
00244 {
00245 POIMs_calc_SLR_helper2( distrib, i, trees[i], leftSubtrees, 0, &dummy, &dummy, &dummy );
00246
00247 const POIMTrie* const node = &TreeMem[ trees[i] ];
00248 ASSERT(trees[i]!=NO_CHILD);
00249
00250 for(symLeft = 0; symLeft < NUM_SYMS; ++symLeft )
00251 leftSubtrees[ symLeft ] = node->children[ symLeft ];
00252 }
00253 }
00254
00255 template <>
00256 void CTrie<POIMTrie>::POIMs_get_SLR( const INT parentIdx, const INT sym, const int depth,
00257 DREAL* S, DREAL* L, DREAL* R )
00258 {
00259 ASSERT(parentIdx!=NO_CHILD);
00260 const POIMTrie* const parent = &TreeMem[ parentIdx ];
00261 if( depth < degree ) {
00262 const INT nodeIdx = parent->children[ sym ];
00263 const POIMTrie* const node = &TreeMem[ nodeIdx ];
00264 *S = node->S;
00265 *L = node->L;
00266 *R = node->R;
00267 } else {
00268 ASSERT(depth==degree);
00269 const DREAL w = parent->child_weights[ sym ];
00270 *S = w;
00271 *L = w;
00272 *R = w;
00273 }
00274 }
00275
00276 template <>
00277 void CTrie<POIMTrie>::POIMs_add_SLR_helper2( DREAL* const* const poims, const int K, const int k, const INT i, const INT y,
00278 const DREAL valW, const DREAL valS, const DREAL valL, const DREAL valR, const INT debug )
00279 {
00280
00281 const INT nk = nofsKmers[ k ];
00282 ASSERT(1<=k && k<=K);
00283 ASSERT(0<=y && y<nk);
00284 INT z;
00285 INT j;
00286
00287
00288 if( debug==0 || debug==2 )
00289 {
00290 poims[ k-1 ][ i*nk + y ] += valS - valW;
00291 }
00292
00293
00294 if( debug==0 || debug==3 )
00295 {
00296 INT r;
00297 for( r = 1; k+r <= K; ++r )
00298 {
00299 const INT nr = nofsKmers[ r ];
00300 const INT nz = nofsKmers[ k+r ];
00301 DREAL* const poim = & poims[ k+r-1 ][ i*nz ];
00302 z = y * nr;
00303 for( j = 0; j < nr; ++j )
00304 {
00305 if( !( 0 <= z && z < nz ) ) {
00306 printf( "k=%d, nk=%d, r=%d, nr=%d, nz=%d \n", k, nk, r, nr, nz );
00307 printf( " j=%d, y=%d, z=%d \n", j, y, z );
00308 }
00309 ASSERT(0<=z && z<nz);
00310 poim[ z ] += valL - valW;
00311 ++z;
00312 }
00313 }
00314 }
00315
00316 if( debug==0 || debug==4 )
00317 {
00318 INT l;
00319 for( l = 1; k+l <= K && l <= i; ++l )
00320 {
00321 const INT nl = nofsKmers[ l ];
00322 const INT nz = nofsKmers[ k+l ];
00323 DREAL* const poim = & poims[ k+l-1 ][ (i-l)*nz ];
00324 z = y;
00325 for( j = 0; j < nl; ++j ) {
00326 ASSERT(0<=z && z<nz);
00327 poim[ z ] += valR - valW;
00328 z += nk;
00329 }
00330 }
00331 }
00332 }
00333
00334 template <>
00335 void CTrie<POIMTrie>::POIMs_add_SLR_helper1( const INT nodeIdx, const int depth, const INT i, const INT y0,
00336 DREAL* const* const poims, const INT K, const INT debug )
00337 {
00338 ASSERT(nodeIdx!=NO_CHILD);
00339 ASSERT(depth<K);
00340 POIMTrie* const node = &TreeMem[ nodeIdx ];
00341 INT sym;
00342 if( depth < degree-1 )
00343 {
00344 if( depth < K-1 )
00345 {
00346 for( sym = 0; sym < NUM_SYMS; ++sym )
00347 {
00348 const INT childIdx = node->children[ sym ];
00349 if( childIdx != NO_CHILD )
00350 {
00351 POIMTrie* const child = &TreeMem[ childIdx ];
00352 const INT y = y0 + sym;
00353 const INT y1 = y * NUM_SYMS;
00354 const DREAL w = child->weight;
00355 POIMs_add_SLR_helper2( poims, K, depth+1, i, y, w, child->S, child->L, child->R, debug );
00356 POIMs_add_SLR_helper1( childIdx, depth+1, i, y1, poims, K, debug );
00357 }
00358 }
00359 }
00360 else
00361 {
00362 ASSERT(depth==K-1);
00363 for( sym = 0; sym < NUM_SYMS; ++sym )
00364 {
00365 const INT childIdx = node->children[ sym ];
00366 if( childIdx != NO_CHILD ) {
00367 POIMTrie* const child = &TreeMem[ childIdx ];
00368 const INT y = y0 + sym;
00369 const DREAL w = child->weight;
00370 POIMs_add_SLR_helper2( poims, K, depth+1, i, y, w, child->S, child->L, child->R, debug );
00371 }
00372 }
00373 }
00374 }
00375 else
00376 {
00377 ASSERT(depth==degree-1);
00378 for( sym = 0; sym < NUM_SYMS; ++sym )
00379 {
00380 const DREAL w = node->child_weights[ sym ];
00381 const INT y = y0 + sym;
00382 POIMs_add_SLR_helper2( poims, K, depth+1, i, y, w, w, w, w, debug );
00383 }
00384 }
00385 }
00386
00387
00388
00389 template <>
00390 void CTrie<POIMTrie>::POIMs_add_SLR( DREAL* const* const poims, const INT K, const INT debug )
00391 {
00392 ASSERT(degree>=1);
00393 ASSERT(K>=1);
00394 {
00395 const int m = ( degree > K ) ? degree : K;
00396 nofsKmers = new INT[ m+1 ];
00397 INT n;
00398 INT k;
00399 n = 1;
00400 for( k = 0; k < m+1; ++k ) {
00401 nofsKmers[ k ] = n;
00402 n *= NUM_SYMS;
00403 }
00404 }
00405 const INT N = length;
00406 INT i;
00407 for( i = 0; i < N; ++i ) {
00408 POIMs_add_SLR_helper1( trees[i], 0, i, 0*NUM_SYMS, poims, K, debug );
00409 }
00410 delete[] nofsKmers;
00411 }