Trie.cpp

Go to the documentation of this file.
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     //printf( "W_helper( %d )\n", i );
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   // --- init
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       // --- go thru direct children
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   // --- collect precalced values from left neighbor tree
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   // --- add w and return results
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   // --- init
00135   node->S = 0;
00136   node->L = 0;
00137   node->R = 0;
00138 
00139   // --- recurse thru direct children
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   // --- collect precalced values from left neighbor tree
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               // - second order correction for S
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   // --- add w and return results
00216   const DREAL w0 = node->weight;
00217   //printf( "  d=%d, node=%d, dS=%.3f, w=%.3f\n", depth, nodeIdx, node->S, w0 );
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     //printf( "i=%d, d=%d, y=%d:  w=%.3f \n", i, k, y, valW );
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     // --- add superstring score; subtract "w", as it was counted twice
00288     if( debug==0 || debug==2 )
00289     {
00290         poims[ k-1 ][ i*nk + y ] += valS - valW;
00291     }
00292 
00293     // --- left partial overlaps
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     // --- right partial overlaps
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 }

SHOGUN Machine Learning Toolbox - Documentation