19 const int32_t nodeIdx,
const int32_t depth,
const int32_t offset,
20 const int32_t y0,
float64_t*
const*
const W,
const int32_t K )
24 float64_t*
const W_kiy = & W[ depth ][ offset + y0 ];
25 POIMTrie*
const node = &TreeMem[ nodeIdx ];
28 if( depth < degree-1 )
30 const int32_t offset1 = offset * NUM_SYMS;
31 for( sym = 0; sym < NUM_SYMS; ++sym )
34 const int32_t childIdx = node->children[ sym ];
35 if( childIdx != NO_CHILD )
37 W_kiy[ sym ] = TreeMem[ childIdx ].weight;
41 const int32_t y1 = ( y0 + sym ) * NUM_SYMS;
42 POIMs_extract_W_helper( childIdx, depth+1, offset1, y1, W, K );
50 for( sym = 0; sym < NUM_SYMS; ++sym )
53 W_kiy[ sym ] = node->child_weights[ sym ];
60 float64_t*
const*
const W,
const int32_t K)
64 const int32_t N = length;
66 for( i = 0; i < N; ++i ) {
68 POIMs_extract_W_helper( trees[i], 0, i*NUM_SYMS, 0*NUM_SYMS, W, K );
74 const float64_t*
const distrib,
const int32_t i,
const int32_t nodeIdx,
75 int32_t left_tries_idx[4],
const int32_t depth, int32_t
const lastSym,
81 const float64_t*
const distribLeft = & distrib[ (i-1) * NUM_SYMS ];
82 POIMTrie*
const node = &TreeMem[ nodeIdx ];
93 const float64_t*
const distribRight = & distrib[ (i+depth) * NUM_SYMS ];
96 for( symRight = 0; symRight < NUM_SYMS; ++symRight ) {
97 const float64_t w1 = node->child_weights[ symRight ];
98 const float64_t pRight = distribRight[ symRight ];
106 for( symLeft = 0; symLeft < NUM_SYMS; ++symLeft )
108 if (left_tries_idx[symLeft] != NO_CHILD)
110 POIMTrie* nodeLeft = &TreeMem[left_tries_idx[symLeft]];
113 const float64_t w2 = nodeLeft->child_weights[ lastSym ];
114 const float64_t pLeft = distribLeft[ symLeft ];
134 const float64_t*
const distrib,
const int32_t i,
const int32_t nodeIdx,
138 ASSERT(0<=depth && depth<=degree-2)
141 const float64_t*
const distribLeft = & distrib[ (i-1) * NUM_SYMS ];
142 POIMTrie*
const node = &TreeMem[ nodeIdx ];
155 for( symRight = 0; symRight < NUM_SYMS; ++symRight )
157 const int32_t childIdx = node->
children[ symRight ];
158 if( childIdx != NO_CHILD )
160 if( depth < degree-2 )
162 int32_t new_left_tries_idx[4];
164 for( symLeft = 0; symLeft < NUM_SYMS; ++symLeft )
166 new_left_tries_idx[symLeft]=NO_CHILD;
168 if (left_tries_idx[symLeft] != NO_CHILD)
170 POIMTrie* nodeLeft = &TreeMem[left_tries_idx[symLeft]];
172 new_left_tries_idx[ symLeft ]=nodeLeft->children[ symRight ];
175 POIMs_calc_SLR_helper2( distrib, i, childIdx, new_left_tries_idx, depth+1, &auxS, &dummy, &auxR );
178 POIMs_calc_SLR_helper1( distrib, i, childIdx, left_tries_idx, depth+1, symRight, &auxS, &dummy, &auxR );
182 const float64_t*
const distribRight = & distrib[ (i+depth) * NUM_SYMS ];
183 const float64_t pRight = distribRight[ symRight ];
184 node->S += pRight * auxS;
185 node->R += pRight * auxR;
191 for( symLeft = 0; symLeft < NUM_SYMS; ++symLeft )
193 if (left_tries_idx[symLeft] != NO_CHILD)
195 const POIMTrie* nodeLeft = &TreeMem[left_tries_idx[symLeft]];
197 const float64_t pLeft = distribLeft[ symLeft ];
199 node->S += pLeft * nodeLeft->S;
200 node->L += pLeft * nodeLeft->L;
204 const float64_t*
const distribRight = & distrib[ (i+depth) * NUM_SYMS ];
207 if( depth < degree-2 )
209 for( symRight = 0; symRight < NUM_SYMS; ++symRight )
211 const int32_t childIdxLeft = nodeLeft->children[ symRight ];
212 if( childIdxLeft != NO_CHILD )
214 const POIMTrie*
const childLeft = &TreeMem[ childIdxLeft ];
215 auxS += distribRight[symRight] * childLeft->S;
221 for( symRight = 0; symRight < NUM_SYMS; ++symRight ) {
222 auxS += distribRight[symRight] * nodeLeft->child_weights[ symRight ];
225 node->S -= pLeft* auxS;
251 const int32_t N = length;
254 int32_t leftSubtrees[4];
255 for( symLeft = 0; symLeft < NUM_SYMS; ++symLeft )
256 leftSubtrees[ symLeft ] = NO_CHILD;
258 for(int32_t i = 0; i < N; ++i )
260 POIMs_calc_SLR_helper2( distrib, i, trees[i], leftSubtrees, 0, &dummy, &dummy, &dummy );
262 const POIMTrie*
const node = &TreeMem[ trees[i] ];
263 ASSERT(trees[i]!=NO_CHILD)
265 for(symLeft = 0; symLeft < NUM_SYMS; ++symLeft )
266 leftSubtrees[ symLeft ] = node->children[ symLeft ];
272 const int32_t parentIdx,
const int32_t sym,
const int32_t depth,
275 ASSERT(parentIdx!=NO_CHILD)
276 const POIMTrie*
const parent = &TreeMem[ parentIdx ];
277 if( depth < degree ) {
278 const int32_t nodeIdx = parent->children[ sym ];
279 const POIMTrie*
const node = &TreeMem[ nodeIdx ];
285 const float64_t w = parent->child_weights[ sym ];
294 float64_t*
const*
const poims,
const int32_t K,
const int32_t k,
295 const int32_t i,
const int32_t y,
const float64_t valW,
300 const int32_t nk = nofsKmers[ k ];
307 if( debug==0 || debug==2 )
309 poims[ k-1 ][ i*nk + y ] += valS - valW;
313 if( debug==0 || debug==3 )
316 for( r = 1; k+r <= K; ++r )
318 const int32_t nr = nofsKmers[ r ];
319 const int32_t nz = nofsKmers[ k+r ];
320 float64_t*
const poim = & poims[ k+r-1 ][ i*nz ];
322 for( j = 0; j < nr; ++j )
324 if( !( 0 <= z && z < nz ) ) {
325 SG_PRINT(
"k=%d, nk=%d, r=%d, nr=%d, nz=%d \n", k, nk, r, nr, nz )
326 SG_PRINT(
" j=%d, y=%d, z=%d \n", j, y, z )
329 poim[ z ] += valL - valW;
335 if( debug==0 || debug==4 )
338 for( l = 1; k+l <= K && l <= i; ++l )
340 const int32_t nl = nofsKmers[ l ];
341 const int32_t nz = nofsKmers[ k+l ];
342 float64_t*
const poim = & poims[ k+l-1 ][ (i-l)*nz ];
344 for( j = 0; j < nl; ++j ) {
346 poim[ z ] += valR - valW;
355 const int32_t nodeIdx,
const int32_t depth,
const int32_t i,
356 const int32_t y0,
float64_t*
const*
const poims,
const int32_t K,
361 POIMTrie*
const node = &TreeMem[ nodeIdx ];
363 if( depth < degree-1 )
367 for( sym = 0; sym < NUM_SYMS; ++sym )
369 const int32_t childIdx = node->
children[ sym ];
370 if( childIdx != NO_CHILD )
372 POIMTrie*
const child = &TreeMem[ childIdx ];
373 const int32_t y = y0 + sym;
374 const int32_t y1 = y * NUM_SYMS;
376 POIMs_add_SLR_helper2( poims, K, depth+1, i, y, w, child->S, child->L, child->R, debug );
377 POIMs_add_SLR_helper1( childIdx, depth+1, i, y1, poims, K, debug );
384 for( sym = 0; sym < NUM_SYMS; ++sym )
386 const int32_t childIdx = node->children[ sym ];
387 if( childIdx != NO_CHILD ) {
388 POIMTrie*
const child = &TreeMem[ childIdx ];
389 const int32_t y = y0 + sym;
391 POIMs_add_SLR_helper2( poims, K, depth+1, i, y, w, child->S, child->L, child->R, debug );
399 for( sym = 0; sym < NUM_SYMS; ++sym )
401 const float64_t w = node->child_weights[ sym ];
402 const int32_t y = y0 + sym;
403 POIMs_add_SLR_helper2( poims, K, depth+1, i, y, w, w, w, w, debug );
412 float64_t*
const*
const poims,
const int32_t K,
const int32_t debug)
417 const int32_t m = ( degree > K ) ? degree : K;
418 nofsKmers = SG_MALLOC(int32_t, m+1 );
422 for( k = 0; k < m+1; ++k ) {
427 const int32_t N = length;
429 for( i = 0; i < N; ++i ) {
430 POIMs_add_SLR_helper1( trees[i], 0, i, 0*NUM_SYMS, poims, K, debug );
Template class Trie implements a suffix trie, i.e. a tree in which all suffixes up to a certain lengt...
void POIMs_extract_W_helper(const int32_t nodeIdx, const int32_t depth, const int32_t offset, const int32_t y0, float64_t *const *const W, const int32_t K)