1 /* truncyk.c
2 * DLK
3 *
4 * Fully local alignment of target (sub)sequence to a CM
5 * using truncated-CYK algorithm
6 */
7
8 /************************************************************
9 *
10 * truncyk external API:
11 *
12 * TrCYK_DnC() - Divide and conquer
13 * TrCYK_Inside() - Inside with or without traceback
14 *
15 ************************************************************/
16
17 #include "esl_config.h"
18 #include "p7_config.h"
19 #include "config.h"
20
21 #include <stdio.h>
22 #include <stdlib.h>
23 #include <math.h>
24
25 #include "easel.h"
26 #include "esl_alphabet.h"
27 #include "esl_stack.h"
28 #include "esl_vectorops.h"
29
30 #include "hmmer.h"
31
32 #include "infernal.h"
33
34 /*
35 struct deckpool_s {
36 float ***pool;
37 int n;
38 int nalloc;
39 int block;
40 };
41 */
42
43 /* Structure: AlphaMats_t */
44 typedef struct alphamats_s {
45 float ***J;
46 float ***L;
47 float ***R;
48 float ***T;
49 } AlphaMats_t;
50
51 /* structure: BetaMats_t */
52 typedef struct betamats_s {
53 float ***J;
54 float **L;
55 float **R;
56 /* no T because T only applies at bifurcations, and beta/outside is only calculated on unbifurcated subgraphs */
57 } BetaMats_t;
58
59 /* Structure: ShadowMats_t */
60 typedef struct shadowmats_s {
61 void ***J;
62 void ***L;
63 void ***Lmode;
64 void ***R;
65 void ***Rmode;
66 void ***T;
67 } ShadowMats_t;
68
69 /* Divide and conquer */
70 float tr_generic_splitter(CM_t *cm, ESL_DSQ *dsq, int L, Parsetree_t *tr,
71 int r, int vend, int i0, int j0,
72 int r_allow_J, int r_allow_L, int r_allow_R);
73 float tr_wedge_splitter(CM_t *cm, ESL_DSQ *dsq, int L, Parsetree_t *tr,
74 int r, int z, int i0, int j0,
75 int r_allow_J, int r_allow_L, int r_allow_R);
76 void tr_v_splitter(CM_t *cm, ESL_DSQ *dsq, int L, Parsetree_t *tr,
77 int r, int z, int i0, int i1, int j1, int j0,
78 int useEL, int r_allow_J, int r_allow_L, int r_allow_R,
79 int z_allow_J, int z_allow_L, int z_allow_R);
80
81 /* Alignment engines */
82 /* trinside is legacy, aviod use! */
83 float trinside (CM_t *cm, ESL_DSQ *dsq, int L, int vroot, int vend, int i0, int j0, int do_full,
84 void ****ret_shadow, void ****ret_L_shadow, void ****ret_R_shadow,
85 void ****ret_T_shadow, void ****ret_Lmode_shadow, void ****ret_Rmode_shadow,
86 int *ret_mode, int *ret_v, int *ret_i, int *ret_j);
87 float tr_inside(CM_t *cm, ESL_DSQ *dsq, int L, int vroot, int vend, int i0, int j0, int do_full,
88 int allow_begin, int r_allow_J, int r_allow_L, int r_allow_R, int lenCORREX,
89 AlphaMats_t *arg_alpha, AlphaMats_t *ret_alpha,
90 struct deckpool_s *dpool, struct deckpool_s **ret_dpool,
91 ShadowMats_t *ret_shadow, int *ret_mode, int *ret_v, int *ret_i, int *ret_j);
92 float tr_outside(CM_t *cm, ESL_DSQ *dsq, int L, int vroot, int vend, int i0, int j0, int do_full,
93 int r_allow_J, int r_allow_L, int r_allow_R,
94 BetaMats_t *arg_beta, BetaMats_t *ret_beta,
95 struct deckpool_s *dpool, struct deckpool_s **ret_dpool,
96 int *ret_mode, int *ret_v, int *ret_j);
97 float tr_vinside(CM_t *cm, ESL_DSQ *dsq, int L, int r, int z, int i0, int i1, int j1, int j0,
98 int useEL, int do_full, int allow_begin,
99 int r_allow_J, int r_allow_L, int r_allow_R,
100 int z_allow_J, int z_allow_L, int z_allow_R,
101 AlphaMats_t *arg_alpha, AlphaMats_t *ret_alpha,
102 struct deckpool_s *dpool, struct deckpool_s **ret_dpool,
103 ShadowMats_t *ret_shadow, int *ret_mode, int *ret_v, int *ret_i, int *ret_j);
104 void tr_voutside(CM_t *cm, ESL_DSQ *dsq, int L, int r, int z, int i0, int i1, int j1, int j0,
105 int useEL, int do_full, int r_allow_J, int r_allow_L, int r_allow_R,
106 int z_allow_J, int z_allow_L, int z_allow_R, BetaMats_t *arg_beta,
107 BetaMats_t *ret_beta, struct deckpool_s *dpool, struct deckpool_s **ret_dpool);
108
109 /* Traceback routine */
110 float tr_insideT(CM_t *cm, ESL_DSQ *dsq, int L, Parsetree_t *tr, int r, int z, int i0, int j0,
111 int r_allow_J, int r_allow_L, int r_allow_R, int lenCORREX);
112 float tr_vinsideT(CM_t *cm, ESL_DSQ *dsq, int L, Parsetree_t *tr, int r, int z,
113 int i0, int i1, int j1, int j0, int useEL,
114 int r_allow_J, int r_allow_L, int r_allow_R,
115 int z_allow_J, int z_allow_L, int z_allow_R);
116
117 /* Function: SetMarginalScores_reproduce_i27()
118 * Author: DLK
119 *
120 * Purpose: Given an otherwise initialized CM,
121 * set marginalized emission score vectors.
122 * Requires cm->abc and cm->esc.
123 *
124 * EPN, Thu Aug 25 09:21:13 2011
125 * This function contains an unfixed version of bug i27 from
126 * BUGTRAX in that it calls
127 * LeftMarginalScore_reproduce_i27 and
128 * RightMarginalScore_reproduce_i27. See those functions
129 * below for more information.
130 *
131 * Args: cm
132 *
133 * Returns: none (cm is modified)
134 */
135 void
SetMarginalScores_reproduce_i27(CM_t * cm)136 SetMarginalScores_reproduce_i27(CM_t *cm)
137 {
138 int i,v;
139
140 cm->lmesc = malloc(sizeof(float *) * (cm->M));
141 cm->rmesc = malloc(sizeof(float *) * (cm->M));
142 cm->ilmesc = malloc(sizeof(float *) * (cm->M));
143 cm->irmesc = malloc(sizeof(float *) * (cm->M));
144
145 cm->lmesc[0] = malloc(sizeof(float) * (cm->M*cm->abc->Kp));
146 cm->rmesc[0] = malloc(sizeof(float) * (cm->M*cm->abc->Kp));
147
148 for (v = 0; v < cm->M; v++)
149 {
150 cm->lmesc[v] = cm->lmesc[0] + v*cm->abc->Kp;
151 cm->rmesc[v] = cm->rmesc[0] + v*cm->abc->Kp;
152
153 if (cm->sttype[v] == MP_st)
154 for (i = 0; i < cm->abc->Kp; i++)
155 {
156 cm->lmesc[v][i] = LeftMarginalScore_reproduce_i27(cm->abc, cm->esc[v], i);
157 cm->rmesc[v][i] = RightMarginalScore_reproduce_i27(cm->abc, cm->esc[v], i);
158 }
159 else if (cm->sttype[v] == ML_st || cm->sttype[v] == IL_st)
160 for (i = 0; i < cm->abc->Kp; i++)
161 {
162 cm->lmesc[v][i] = cm->esc[v][i];
163 cm->rmesc[v][i] = 0.0;
164 }
165 else if (cm->sttype[v] == MR_st || cm->sttype[v] == IR_st)
166 for (i = 0; i < cm->abc->Kp; i++)
167 {
168 cm->lmesc[v][i] = 0.0;
169 cm->rmesc[v][i] = cm->esc[v][i];
170 }
171 else
172 for (i = 0; i < cm->abc->Kp; i++)
173 {
174 cm->lmesc[v][i] = 0.0;
175 cm->rmesc[v][i] = 0.0;
176 }
177 }
178
179 return;
180 }
181
182
183 /* Function: LeftMarginalScore_reproduce_i27()
184 * Author: DLK
185 *
186 * Purpose: Calculate marginal probability for left half
187 * of an emission pair. Implicitly assumes
188 * a uniform background distribution
189 *
190 * EPN, Thu Aug 25 09:20:32 2011 This function contains an
191 * unfixed version of bug i27 from BUGTRAX. The esc vector
192 * is log_2 odds scores but esl_vec_FLogSum() works in log
193 * base e, therefore the scores returned from this function
194 * were wrong and corresponded to emission probabilities
195 * that do not sum to 1.0 across the canonical residues of
196 * the alphabet (when called successively for
197 * dres=0..abc->K-1). I left this in to allow reproduction
198 * of the benchmark from Kolbe, Eddy 2009 paper on truncated
199 * CYK, the code for which included this bug. The fixed
200 * version is now incorporated into CMLogoddsify, which now
201 * calculates marginal scores and thus removes the need for
202 * this function.
203 */
204 float
LeftMarginalScore_reproduce_i27(const ESL_ALPHABET * abc,float * esc,ESL_DSQ dres)205 LeftMarginalScore_reproduce_i27(const ESL_ALPHABET *abc, float *esc, ESL_DSQ dres)
206 {
207 float *left = NULL;
208 int status;
209 ESL_ALLOC(left, (sizeof(float) * (abc->K+1)));
210 int i;
211 float sc;
212
213 if (dres < abc->K)
214 {
215 sc = esl_vec_FLogSum(&(esc[dres*abc->K]),abc->K);
216 sc -= sreLOG2(abc->K);
217 }
218 else /* degenerate */
219 {
220 esl_vec_FSet(left, abc->K, 0.);
221 esl_abc_FCount(abc, left, dres, 1.);
222
223 sc = 0.;
224 for (i = 0; i < abc->K; i++)
225 {
226 sc += esl_vec_FLogSum(&(esc[i*abc->K]),abc->K)*left[i];
227 sc -= sreLOG2(abc->K)*left[i];
228 }
229 }
230
231 free(left);
232
233 return sc;
234
235 ERROR:
236 cm_Fail("Memory allocation error.");
237 return 0.0; /* never reached */
238 }
239
240 /* Function: RightMarginalScore_reproduce_i27()
241 * Author: DLK
242 *
243 * Purpose: Calculate marginal probability for right half
244 * of an emission pair. Implicitly assumes
245 * a uniform background distribution.
246 *
247 * EPN, Thu Aug 25 09:18:43 2011
248 * This contains bug i27 from BUGTRAX. The esc vector is
249 * log_2 odds scores but esl_vec_FLogSum() works in log base
250 * e, therefore the scores returned from this function were
251 * wrong and corresponded to emission probabilities that do
252 * not sum to 1.0 across the canonical residues of the
253 * alphabet (when called successively for dres=0..abc->K-1).
254 * I left this in to allow reproduction of the benchmark
255 * from Kolbe, Eddy 2009 paper on truncated CYK, the code
256 * for which included this bug. The fixed version is now
257 * incorporated into CMLogoddsify, which now calculates
258 * marginal scores and thus removes the need for this
259 * function.
260 */
261 float
RightMarginalScore_reproduce_i27(const ESL_ALPHABET * abc,float * esc,ESL_DSQ dres)262 RightMarginalScore_reproduce_i27(const ESL_ALPHABET *abc, float *esc, ESL_DSQ dres)
263 {
264 float *right = NULL;
265 int status;
266 int i,j;
267 float sc;
268 float row[abc->K];
269 ESL_ALLOC(right, (sizeof(float) * (abc->K+1)));
270
271 if (dres < abc->K)
272 {
273 for (i=0; i<abc->K; i++)
274 row[i] = esc[i*abc->K+dres];
275 sc = esl_vec_FLogSum(row,abc->K);
276 sc -= sreLOG2(abc->K);
277 }
278 else /* degenerate */
279 {
280 esl_vec_FSet(right, abc->K, 0.);
281 esl_abc_FCount(abc, right, dres, 1.);
282
283 sc = 0.;
284 for (i=0; i < abc->K; i++)
285 {
286 for (j=0; j<abc->K; j++)
287 row[j] = esc[j*abc->K+dres];
288 sc += esl_vec_FLogSum(row,abc->K)*right[i];
289 sc -= sreLOG2(abc->K)*right[i];
290 }
291 }
292
293 free(right);
294
295 return sc;
296
297 ERROR:
298 cm_Fail("Memory allocation error.");
299 return 0.0; /* never reached */
300 }
301
302 /* Function: TrCYK_DnC()
303 * Author: DLK
304 *
305 * Purpose: Divide-and-conquer CYK alignment
306 * for truncated sequences with traceback
307 *
308 * Args: cm - the covariance model
309 * dsq - the sequence, 1..L
310 * L - length of the sequence
311 * r - root of subgraph to align to target subseq (usually 0, the model's root)
312 * i0 - start of target subsequence (usually 1, beginning of dsq)
313 * j0 - end of target subsequence (usually L, end of dsq)
314 * pass_idx - pipeline pass index, tr->pass_idx set as this
315 * do_1p0 - TRUE: behave like the 1.0 version;
316 * FALSE: be consistent with cm_dpalign_trunc.c functions
317 * ret_tr - RETURN: traceback (pass NULL if trace isn't wanted)
318 *
319 * Returns: score of the alignment in bits
320 */
321 float
TrCYK_DnC(CM_t * cm,ESL_DSQ * dsq,int L,int r,int i0,int j0,int pass_idx,int do_1p0,Parsetree_t ** ret_tr)322 TrCYK_DnC(CM_t *cm, ESL_DSQ *dsq, int L, int r, int i0, int j0, int pass_idx, int do_1p0, Parsetree_t **ret_tr)
323 {
324 Parsetree_t *tr;
325 int z;
326 float sc, bsc;
327 int pty_idx; /* index for truncation penalty, determined by pass_idx */
328
329 /* Check input parameters */
330 if ( cm->stid[r] != ROOT_S )
331 {
332 if (! (cm->flags & CMH_LOCAL_BEGIN)) cm_Die("internal error: we're not in local mode, but r is not root");
333 if ( (cm->stid[r] != MATP_MP) &&
334 (cm->stid[r] != MATL_ML) &&
335 (cm->stid[r] != MATR_MR) &&
336 (cm->stid[r] != BIF_B ) ) cm_Die("internal error: trying to do a local begin at a non-mainline start");
337 }
338
339 /* Create parse tree and initialize */
340 tr = CreateParsetree(100);
341 tr->is_std = FALSE; /* lower is_std flag, now we'll know this parsetree was created by a truncated (non-standard) alignment function */
342 tr->pass_idx = pass_idx;
343
344 if(! do_1p0) InsertTraceNode(tr, -1, TRACE_LEFT_CHILD, i0, j0, 0); /* trcyk 1.0 did not add state 0 */
345 z = cm->M-1;
346
347 /* If local begin is known */
348 if ( r != 0 )
349 {
350 InsertTraceNode(tr, -1, TRACE_LEFT_CHILD, i0, j0, r);
351 z = CMSubtreeFindEnd(cm, r);
352 }
353
354 /* Solve by calling tr_generic_splitter() */
355 sc = tr_generic_splitter(cm, dsq, L, tr, r, z, i0, j0, TRUE, TRUE, TRUE);
356
357 if(! do_1p0) {
358 /* modify mode of initial node(s) now that we know the mode of the full alignment */
359 if(r == 0) {
360 tr->mode[0] = tr->mode[1];
361 }
362 else {
363 tr->mode[0] = tr->mode[2];
364 tr->mode[1] = tr->mode[2];
365 }
366 }
367
368 if(do_1p0) {
369 /* 2.0 instead of 2 to force floating point division, not integer division */
370 /* This is the fragment penalty */
371 bsc = sreLOG2(2.0/(cm->clen*(cm->clen+1)));
372 /* printf("Best truncated score: %.4f (%.4f)\n", sc, (sc+bsc)); */
373 }
374 else {
375 if((pty_idx = cm_tr_penalties_IdxForPass(pass_idx)) == -1) cm_Die("TrCYK_DnC, unexpected pass idx: %d", pass_idx);
376 bsc = (cm->flags & CMH_LOCAL_BEGIN) ? cm->trp->l_ptyAA[pty_idx][tr->state[1]] : cm->trp->g_ptyAA[pty_idx][tr->state[1]];
377 }
378 tr->trpenalty = bsc;
379 sc += bsc;
380
381 if ( ret_tr != NULL ) { *ret_tr = tr; }
382 else { FreeParsetree(tr); }
383
384 return sc;
385 }
386
387 /* Function: TrCYK_Inside()
388 * Author: DLK
389 *
390 * Purpose: Full CYK alignment for truncated sequences
391 * with traceback
392 *
393 * Based on CYKInside()
394 *
395 * Args: cm - the covariance model
396 * dsq - the sequence, 1..L
397 * L - length of the sequence
398 * r - root of subgraph to align to target subseq (usually 0, the model's root)
399 * i0 - start of target subsequence (usually 1, beginning of dsq)
400 * j0 - end of target subsequence (usually L, end of dsq)
401 * pass_idx - pipeline pass index, tr->pass_idx set as this
402 * ret_tr - RETURN: traceback (pass NULL if trace isn't wanted)
403 *
404 * Returns; score of the alignment in bits
405 */
406 float
TrCYK_Inside(CM_t * cm,ESL_DSQ * dsq,int L,int r,int i0,int j0,int pass_idx,int do_1p0,int lenCORREX,Parsetree_t ** ret_tr)407 TrCYK_Inside(CM_t *cm, ESL_DSQ *dsq, int L, int r, int i0, int j0, int pass_idx, int do_1p0, int lenCORREX, Parsetree_t **ret_tr)
408 {
409 Parsetree_t *tr = NULL;
410 int z;
411 float sc, bsc, rsc, nullsc;
412 int pty_idx; /* index for truncation penalty, determined by pass_idx */
413
414 /* Check input parameters */
415 if ( cm->stid[r] != ROOT_S )
416 {
417 if (! (cm->flags & CMH_LOCAL_BEGIN)) cm_Die("internal error: we're not in local mode, but r is not root");
418 if ( (cm->stid[r] != MATP_MP) &&
419 (cm->stid[r] != MATL_ML) &&
420 (cm->stid[r] != MATR_MR) &&
421 (cm->stid[r] != BIF_B ) ) cm_Die("internal error: trying to do a local begin at a non-mainline start");
422 }
423
424 if ( ret_tr != NULL)
425 {
426 /* Create parse tree and initialize */
427 tr = CreateParsetree(100);
428 tr->is_std = FALSE; /* lower is_std flag, now we'll know this parsetree was created by a truncated (non-standard) alignment function */
429 tr->pass_idx = pass_idx;
430
431 if(! do_1p0) InsertTraceNode(tr, -1, TRACE_LEFT_CHILD, i0, j0, 0); /* trcyk 1.0 did not add state 0 */
432 z = cm->M-1;
433
434 /* If local begin is known */
435 if ( r != 0 )
436 {
437 InsertTraceNode(tr, -1, TRACE_LEFT_CHILD, i0, j0, r);
438 z = CMSubtreeFindEnd(cm, r);
439 }
440
441 /* Solve by calling tr_insideT() */
442 sc = tr_insideT(cm, dsq, L, tr, r, z, i0, j0, TRUE, TRUE, TRUE, lenCORREX);
443
444 if(! do_1p0) {
445 /* modify mode of initial node(s) now that we know the mode of the full alignment */
446 if(r == 0) {
447 tr->mode[0] = tr->mode[1]; /* modify mode of first trace node now that we know the mode of the full alignment */
448 }
449 else {
450 tr->mode[0] = tr->mode[2];
451 tr->mode[1] = tr->mode[2];
452 }
453 }
454 }
455 else
456 {
457 z = (r == 0) ? cm->M-1 : CMSubtreeFindEnd(cm, r);
458 sc = tr_inside(cm, dsq, L, r, z, i0, j0, BE_EFFICIENT,
459 TRUE, TRUE, TRUE, TRUE, lenCORREX,
460 NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL);
461 }
462
463 if(do_1p0) {
464 /* 2.0 instead of 2 to force floating point division, not integer division */
465 /* This is the fragment penalty */
466 bsc = sreLOG2(2.0/(cm->clen*(cm->clen+1)));
467 /* printf("Best truncated score: %.4f (%.4f)\n", sc, (sc+bsc)); */
468 }
469 else {
470 if((pty_idx = cm_tr_penalties_IdxForPass(pass_idx)) == -1) cm_Die("TrCYK_DnC, unexpected pass idx: %d", pass_idx);
471 bsc = (cm->flags & CMH_LOCAL_BEGIN) ? cm->trp->l_ptyAA[pty_idx][tr->state[1]] : cm->trp->g_ptyAA[pty_idx][tr->state[1]];
472 }
473 if(tr != NULL) tr->trpenalty = bsc;
474
475 /* null model length correction */
476 if (lenCORREX)
477 {
478 rsc = (float) L/(float) (L+1);
479 nullsc = L*sreLOG2(rsc) + sreLOG2(1 - rsc);
480 sc -= nullsc;
481 }
482 sc += bsc;
483
484 if ( ret_tr != NULL ) *ret_tr = tr;
485 else if(tr != NULL) { FreeParsetree(tr); }
486
487 return sc;
488 }
489
490 /* Function: tr_generic_splitter()
491 * Author: DLK
492 *
493 * Purpose: Generic problem for divide-and-conquer
494 * Based closely on generic_splitter()
495 *
496 * Args:
497 *
498 * Returns:
499 */
500 float
tr_generic_splitter(CM_t * cm,ESL_DSQ * dsq,int L,Parsetree_t * tr,int r,int z,int i0,int j0,int r_allow_J,int r_allow_L,int r_allow_R)501 tr_generic_splitter(CM_t *cm, ESL_DSQ *dsq, int L, Parsetree_t *tr,
502 int r, int z, int i0, int j0,
503 int r_allow_J, int r_allow_L, int r_allow_R)
504 {
505 AlphaMats_t *alpha;
506 BetaMats_t *beta;
507 struct deckpool_s *pool;
508 int v,w,y;
509 int wend, yend;
510 int tv;
511 int jp;
512 int W;
513 float sc;
514 int j,d,k;
515 float best_sc;
516 int best_j, best_d, best_k;
517 int v_mode, w_mode, y_mode;
518 int b1_mode, b2_mode, b3_mode;
519 int b1_v, b1_i, b1_j;
520 int b2_v, b2_i, b2_j;
521 int b3_v, b3_j;
522 float b1_sc, b2_sc, b3_sc;
523 int useEL;
524
525 int v_allow_T = FALSE;
526
527 if (r == 0) v_allow_T = TRUE;
528 if (!r_allow_J && !r_allow_L && !r_allow_R) v_allow_T = TRUE;
529
530 /* Case 1: problem size is small; solve with tr_insideT()
531 * size calculation is heuristic based on size of insideT() */
532 if (5*insideT_size(cm, L, r, z, i0, j0) < RAMLIMIT)
533 {
534 sc = tr_insideT(cm, dsq, L, tr, r, z, i0, j0, r_allow_J, r_allow_L, r_allow_R, FALSE);
535 return sc;
536 }
537
538 /* Case 2: find a bifurcation */
539 for (v = r; v <= z-5; v++)
540 { if (cm->sttype[v] == B_st) break; }
541
542 /* Case 3: no bifurcations -> wedge problem */
543 if (cm->sttype[v] != B_st)
544 {
545 if (cm->sttype[z] != E_st) cm_Die("z in tr_generic_splitter not E_st - that ain't right");
546 sc = tr_wedge_splitter(cm, dsq, L, tr, r, z, i0, j0, r_allow_J, r_allow_L, r_allow_R);
547 return sc;
548 }
549
550 alpha = malloc(sizeof(AlphaMats_t));
551 beta = malloc(sizeof(BetaMats_t));
552
553 /* Unusual cases dispatched, back to case 2 (bifurcation) */
554 w = cm->cfirst[v];
555 y = cm->cnum[v];
556 if (w < y) { wend = y-1; yend = z; }
557 else { yend = w-1; wend = z; }
558
559 /* Calculate alphas for w and y
560 * also pick up best local begins in each subtree */
561 b1_sc = tr_inside(cm, dsq, L, w, wend, i0, j0, BE_EFFICIENT,
562 (r == 0), TRUE, r_allow_L, r_allow_R, FALSE,
563 NULL, alpha, NULL, &pool, NULL, &b1_mode, &b1_v, &b1_i, &b1_j);
564 if (r != 0) b1_sc = IMPOSSIBLE;
565 b2_sc = tr_inside(cm, dsq, L, y, yend, i0, j0, BE_EFFICIENT,
566 (r == 0), TRUE, r_allow_L, r_allow_R, FALSE,
567 alpha, alpha, pool, NULL, NULL, &b2_mode, &b2_v, &b2_i, &b2_j);
568 if (r != 0) b2_sc = IMPOSSIBLE;
569
570 /* Calculate beta; release pool */
571 b3_sc = tr_outside(cm, dsq, L, r, v, i0, j0, BE_EFFICIENT,
572 r_allow_J, r_allow_L, r_allow_R,
573 NULL, beta, NULL, NULL, &b3_mode, &b3_v, &b3_j);
574
575 /* OK, to the point of actually finding the best split
576 * We have a lot more types of splits than the non-truncated
577 * version, so we need a better way to keep track of them */
578 W = j0 - i0 + 1;
579 best_sc = IMPOSSIBLE;
580 for (jp = 0; jp <= W; jp ++)
581 {
582 j = i0 - 1 + jp;
583 for (d = 0; d <= jp; d++)
584 {
585 for (k = 0; k <= d; k++)
586 {
587 /* Attempted bug fix cases - these cases have priority */
588 if ( v_allow_T && k > 0 && k < d)
589 if ( (sc = alpha->J[w][j-k][d-k] + alpha->L[y][j][k]) > best_sc )
590 {
591 best_sc = sc;
592 best_k = k;
593 best_j = j;
594 best_d = d;
595 v_mode = 0; w_mode = 3; y_mode = 2;
596 }
597 if ( v_allow_T && k > 0 && k < d)
598 if ( (sc = alpha->R[w][j-k][d-k] + alpha->J[y][j][k]) > best_sc )
599 {
600 best_sc = sc;
601 best_k = k;
602 best_j = j;
603 best_d = d;
604 v_mode = 0; w_mode = 1; y_mode = 3;
605 }
606 if ( (sc = alpha->J[w][j-k][d-k] + alpha->J[y][j][k] + beta->J[v][j][d]) > best_sc )
607 {
608 best_sc = sc;
609 best_k = k;
610 best_j = j;
611 best_d = d;
612 v_mode = 3; w_mode = 3; y_mode = 3;
613 }
614 if ( r_allow_L && k > 0 /* && j-d+1 > i0 */ )
615 if ( (sc = alpha->J[w][j-k][d-k] + alpha->L[y][j][k] + beta->L[v][j-d+1]) > best_sc )
616 {
617 best_sc = sc;
618 best_k = k;
619 best_j = j;
620 best_d = d;
621 v_mode = 2; w_mode = 3; y_mode = 2;
622 }
623 /* Attempted bug fix case */
624 if ( r_allow_L && k > 0 )
625 if ( (sc = alpha->J[w][j-k][d-k] + alpha->J[y][j][k] + beta->L[v][j-d+1]) > best_sc )
626 {
627 best_sc = sc;
628 best_k = k;
629 best_j = j;
630 best_d = d;
631 v_mode = 2; w_mode = 3; y_mode = 3;
632 }
633 /* j < j0 test causes problems if there are no R emitters between r and v */
634 if ( r_allow_R && k < d /* && j < j0*/ )
635 if ( (sc = alpha->R[w][j-k][d-k] + alpha->J[y][j][k] + beta->R[v][j]) > best_sc )
636 {
637 best_sc = sc;
638 best_k = k;
639 best_j = j;
640 best_d = d;
641 v_mode = 1; w_mode = 1; y_mode = 3;
642 }
643 /* Attempted bug fix case */
644 if ( r_allow_R && k < d )
645 if ( (sc = alpha->J[w][j-k][d-k] + alpha->J[y][j][k] + beta->R[v][j]) > best_sc )
646 {
647 best_sc = sc;
648 best_k = k;
649 best_j = j;
650 best_d = d;
651 v_mode = 1; w_mode = 3; y_mode = 3;
652 }
653 if ( v_allow_T && k > 0 && k < d )
654 if ( (sc = alpha->R[w][j-k][d-k] + alpha->L[y][j][k]) > best_sc )
655 {
656 best_sc = sc;
657 best_k = k;
658 best_j = j;
659 best_d = d;
660 v_mode = 0; w_mode = 1; y_mode = 2;
661 }
662 }
663
664 if ( r_allow_L )
665 if ( (sc = alpha->L[w][j][d] + beta->L[v][j-d+1]) > best_sc )
666 {
667 best_sc = sc;
668 best_k = 0;
669 best_j = j;
670 best_d = d;
671 v_mode = 2; w_mode = 2; y_mode = 0;
672 }
673 if ( r_allow_L )
674 if ( (sc = alpha->J[w][j][d] + beta->L[v][j-d+1]) > best_sc )
675 {
676 best_sc = sc;
677 best_k = 0;
678 best_j = j;
679 best_d = d;
680 v_mode = 2; w_mode = 3; y_mode = 0;
681 }
682 if ( r_allow_R )
683 if ( (sc = alpha->R[y][j][d] + beta->R[v][j]) > best_sc )
684 {
685 best_sc = sc;
686 best_k = d;
687 best_j = j;
688 best_d = d;
689 v_mode = 1; w_mode = 0; y_mode = 1;
690 }
691 if ( r_allow_R )
692 if ( (sc = alpha->J[y][j][d] + beta->R[v][j]) > best_sc )
693 {
694 best_sc = sc;
695 best_k = d;
696 best_j = j;
697 best_d = d;
698 v_mode = 1; w_mode = 0; y_mode = 3;
699 }
700 if ( (sc = beta->J[cm->M][j][d]) > best_sc) /* Joint parent to EL */
701 {
702 best_sc = sc;
703 best_k = -1;
704 best_j = j;
705 best_d = d;
706 v_mode = 3; w_mode = 0; y_mode = 0;
707 useEL = TRUE;
708 }
709 }
710 }
711
712 /* Check for local entry in one of the child sub-trees */
713 if (r == 0)
714 {
715 if (b1_sc > best_sc)
716 {
717 best_sc = b1_sc;
718 best_k = b1_v;
719 best_j = b1_j;
720 best_d = b1_j - b1_i + 1;
721 v_mode = 0; w_mode = b1_mode; y_mode = 0;
722 }
723 if (b2_sc > best_sc)
724 {
725 best_sc = b2_sc;
726 best_k = b2_v;
727 best_j = b2_j;
728 best_d = b2_j - b2_i + 1;
729 v_mode = 0; w_mode = 0; y_mode = b2_mode;
730 }
731 }
732
733 /* local hit in parent (must be marginal) */
734 if (b3_sc > best_sc)
735 {
736 best_sc = b3_sc;
737 best_k = b3_v;
738 best_j = b3_j;
739 v_mode = b3_mode; w_mode = 0; y_mode = 0;
740 useEL = FALSE;
741 }
742
743 /* Free alphas */
744 free_vjd_matrix(alpha->J, cm->M, i0, j0);
745 free_vjd_matrix(alpha->L, cm->M, i0, j0);
746 free_vjd_matrix(alpha->R, cm->M, i0, j0);
747 free_vjd_matrix(alpha->T, cm->M, i0, j0);
748 free_vjd_matrix( beta->J, cm->M, i0, j0);
749 free(beta->L[0]); free(beta->L);
750 free(beta->R[0]); free(beta->R);
751 free(alpha);
752 free(beta);
753
754 /* Found the best path, now to interpret and sub-divide */
755 if ( v_mode ) /* parent graph is non-empty */
756 {
757 if ( w_mode == TRMODE_T && y_mode == TRMODE_T ) /* local hit in parent (marginal) */
758 {
759 if (!useEL && b3_v == -1)
760 cm_Die("1Superbad: passing z = -1!\n");
761 tr_v_splitter(cm, dsq, L, tr, r, (useEL ? v : b3_v), i0, best_j, best_j, j0,
762 useEL, r_allow_J, r_allow_L, r_allow_R, (v_mode == TRMODE_J), (v_mode == TRMODE_L), (v_mode == TRMODE_R));
763 return best_sc;
764 }
765 else
766 {
767 if (v == -1) cm_Die("2Superbad: passing z = -1!\n");
768 tr_v_splitter(cm, dsq, L, tr, r, v, i0, best_j-best_d+1, best_j, j0,
769 FALSE, r_allow_J, r_allow_L, r_allow_R, (v_mode == TRMODE_J), (v_mode == TRMODE_L), (v_mode == TRMODE_R));
770 }
771 }
772 else if ( w_mode == TRMODE_T || y_mode == TRMODE_T ) /* local entry to one of the children */
773 {
774 if ( b1_sc > b2_sc )
775 {
776 InsertTraceNodewithMode(tr, tr->n-1, TRACE_LEFT_CHILD, b1_i, b1_j, b1_v, b1_mode);
777 z = CMSubtreeFindEnd(cm, b1_v);
778 tr_generic_splitter(cm, dsq, L, tr, b1_v, z, b1_i, b1_j, (b1_mode == TRMODE_J), (b1_mode == TRMODE_L), (b1_mode == TRMODE_R));
779 return best_sc;
780 }
781 else
782 {
783 InsertTraceNodewithMode(tr, tr->n-1, TRACE_LEFT_CHILD, b2_i, b2_j, b2_v, b2_mode);
784 z = CMSubtreeFindEnd(cm, b2_v);
785 tr_generic_splitter(cm, dsq, L, tr, b2_v, z, b2_i, b2_j, (b2_mode == TRMODE_J), (b2_mode == TRMODE_L), (b2_mode == TRMODE_R));
786 return best_sc;
787 }
788 }
789 else /* case T: parent is empty, but both children are non-empty */
790 {
791 InsertTraceNodewithMode(tr, tr->n-1, TRACE_LEFT_CHILD, best_j - best_d + 1, best_j, v, 0);
792 }
793
794 tv = tr->n - 1;
795 if ( w_mode )
796 {
797 InsertTraceNodewithMode(tr, tv, TRACE_LEFT_CHILD, best_j - best_d + 1, best_j - best_k, w, w_mode);
798 tr_generic_splitter(cm, dsq, L, tr, w, wend, best_j - best_d + 1, best_j - best_k, (w_mode == TRMODE_J), (w_mode == TRMODE_L), (w_mode == TRMODE_R)); }
799 else
800 {
801 InsertTraceNodewithMode(tr, tr->n-1, TRACE_LEFT_CHILD, best_j - best_d + 1, best_j - best_d, w, w_mode);
802 /*InsertTraceNodewithMode(tr, tr->n-1, TRACE_LEFT_CHILD, best_j - best_d + 1, best_j - best_d, cm->M, 3);*/
803 }
804
805 if ( y_mode )
806 {
807 InsertTraceNodewithMode(tr, tv, TRACE_RIGHT_CHILD, best_j - best_k + 1, best_j, y, y_mode);
808 tr_generic_splitter(cm, dsq, L, tr, y, yend, best_j - best_k + 1, best_j, (y_mode == TRMODE_J), (y_mode == TRMODE_L), (y_mode == TRMODE_R));
809 }
810 else
811 {
812 InsertTraceNodewithMode(tr, tv, TRACE_RIGHT_CHILD, best_j + 1, best_j, y, y_mode);
813 /*InsertTraceNodewithMode(tr, tv, TRACE_RIGHT_CHILD, best_j + 1, best_j, cm->M, 3);*/
814 }
815
816 return best_sc;
817 }
818
819 /* Function: tr_wedge_splitter()
820 * Author: DLK
821 *
822 * Purpose: Wedge problem for divide-and-conquer
823 * Based closely on wedge_splitter()
824 *
825 * Args:
826 *
827 * Returns:
828 */
829 float
tr_wedge_splitter(CM_t * cm,ESL_DSQ * dsq,int L,Parsetree_t * tr,int r,int z,int i0,int j0,int r_allow_J,int r_allow_L,int r_allow_R)830 tr_wedge_splitter(CM_t *cm, ESL_DSQ *dsq, int L, Parsetree_t *tr,
831 int r, int z, int i0, int j0,
832 int r_allow_J, int r_allow_L, int r_allow_R)
833 {
834 AlphaMats_t *alpha;
835 BetaMats_t *beta;
836 float sc;
837 float best_sc;
838 int v,w,y;
839 int W;
840 int jp, j, d;
841 int best_v, best_j, best_d;
842 int p_mode, c_mode;
843 int midnode;
844 float b1_sc, b2_sc;
845 int b1_mode, b1_v, b1_i, b1_j;
846 int b2_mode, b2_v, b2_j;
847
848 /* Special case: problem is small enough to be solved with traceback */
849 if ( (cm->ndidx[z] == cm->ndidx[r] + 1) ||
850 (5 * insideT_size(cm, L, r, z, i0, j0) < RAMLIMIT) )
851 {
852 sc = tr_insideT(cm, dsq, L, tr, r, z, i0, j0, r_allow_J, r_allow_L, r_allow_R, FALSE);
853 return sc;
854 }
855
856 alpha = malloc(sizeof(AlphaMats_t));
857 beta = malloc(sizeof(BetaMats_t));
858
859 /* Calculate a midpoint to split at */
860 midnode = cm->ndidx[r] + ((cm->ndidx[z] - cm->ndidx[r])/2);
861 w = cm->nodemap[midnode];
862 y = cm->cfirst[w] - 1;
863
864 /* Get alphas and betas */
865 b1_sc = tr_inside(cm, dsq, L, w, z, i0, j0, BE_EFFICIENT,
866 (r == 0), TRUE, r_allow_L, r_allow_R, FALSE,
867 NULL, alpha, NULL, NULL, NULL, &b1_mode, &b1_v, &b1_i, &b1_j);
868 if (r != 0) b1_sc = IMPOSSIBLE;
869 b2_sc = tr_outside(cm, dsq, L, r, y, i0, j0, BE_EFFICIENT,
870 r_allow_J, r_allow_L, r_allow_R,
871 NULL, beta, NULL, NULL, &b2_mode, &b2_v, &b2_j);
872 if ( b2_mode == TRMODE_L && !r_allow_L ) b2_sc = IMPOSSIBLE;
873 if ( b2_mode == TRMODE_R && !r_allow_R ) b2_sc = IMPOSSIBLE;
874
875 /* Find the split */
876 W = j0 - i0 + 1;
877 best_sc = IMPOSSIBLE;
878
879 /* Special case: parent empty, child has local hit */
880 if (b1_sc > best_sc)
881 {
882 best_sc = b1_sc;
883 best_v = b1_v;
884 best_j = b1_j;
885 best_d = b1_j - b1_i + 1;
886 p_mode = 0; c_mode = b1_mode;
887 }
888
889 /* Special case: child empty, parent has local hit */
890 /* 1 and 2 are the only appropriate values for b2_mode */
891 if (b2_sc > best_sc)
892 {
893 best_sc = b2_sc;
894 best_v = b2_v;
895 best_j = b2_j;
896 best_d = 1;
897 p_mode = b2_mode; c_mode = 0;
898 }
899
900 /* Standard cases */
901 for (v = w; v <= y; v++)
902 {
903 for (jp = 0; jp <= W; jp++)
904 {
905 j = i0 - 1 + jp;
906 for (d = 0; d <= jp; d++)
907 {
908 if ( (sc = alpha->J[v][j][d] + beta->J[v][j][d]) > best_sc )
909 {
910 best_sc = sc;
911 best_v = v;
912 best_d = d;
913 best_j = j;
914 p_mode = 3; c_mode = 3;
915 }
916 if ( r_allow_L )
917 if ( (sc = alpha->J[v][j][d] + beta->L[v][j-d+1]) > best_sc )
918 {
919 best_sc = sc;
920 best_v = v;
921 best_d = d;
922 best_j = j;
923 p_mode = 2; c_mode = 3;
924 }
925 if ( r_allow_L )
926 if ( (sc = alpha->L[v][j][d] + beta->L[v][j-d+1]) > best_sc )
927 {
928 best_sc = sc;
929 best_v = v;
930 best_d = d;
931 best_j = j;
932 p_mode = 2; c_mode = 2;
933 }
934 if ( r_allow_R )
935 if ( (sc = alpha->J[v][j][d] + beta->R[v][j]) > best_sc )
936 {
937 best_sc = sc;
938 best_v = v;
939 best_d = d;
940 best_j = j;
941 p_mode = 1; c_mode = 3;
942 }
943 if ( r_allow_R )
944 if ( (sc = alpha->R[v][j][d] + beta->R[v][j]) > best_sc )
945 {
946 best_sc = sc;
947 best_v = v;
948 best_d = d;
949 best_j = j;
950 p_mode = 1; c_mode = 1;
951 }
952 }
953 }
954 }
955
956 /* Special case: joint parent to EL */
957 for (jp = 0; jp <= W; jp++)
958 {
959 j = i0 - 1 + jp;
960 for (d = 0; d <= jp; d++)
961 {
962 if ( (sc = beta->J[cm->M][j][d]) > best_sc )
963 {
964 best_sc = sc;
965 best_v = -1;
966 best_j = j;
967 best_d = d;
968 p_mode = 3; c_mode = 0;
969 }
970 }
971 }
972
973 /* Free alpha and beta */
974 free_vjd_matrix(alpha->J, cm->M, i0, j0);
975 free_vjd_matrix(alpha->L, cm->M, i0, j0);
976 free_vjd_matrix(alpha->R, cm->M, i0, j0);
977 free_vjd_matrix(alpha->T, cm->M, i0, j0);
978 free_vjd_matrix( beta->J, cm->M, i0, j0);
979 free(beta->L[0]); free(beta->L);
980 free(beta->R[0]); free(beta->R);
981 free(alpha);
982 free(beta);
983
984 if ( p_mode )
985 {
986 if ( c_mode == TRMODE_T ) /* child empty */
987 {
988 if ((p_mode == TRMODE_J ? w : b2_v) == -1) cm_Die("3Superbad: passing z = -1!\n");
989 tr_v_splitter(cm, dsq, L, tr, r, (p_mode == TRMODE_J ? w : b2_v), i0, best_j - best_d + 1, best_j, j0,
990 (p_mode == TRMODE_J), r_allow_J, r_allow_L, r_allow_R, (p_mode == TRMODE_J), (p_mode == TRMODE_L), (p_mode == TRMODE_R));
991 return best_sc;
992 }
993 else
994 {
995 if (best_v == -1) cm_Die("4Superbad: passing z = -1!\n");
996 tr_v_splitter(cm, dsq, L, tr, r, best_v, i0, best_j - best_d + 1, best_j, j0,
997 FALSE, r_allow_J, r_allow_L, r_allow_R, (c_mode == TRMODE_J), (c_mode == TRMODE_L), (c_mode == TRMODE_R));
998 }
999 }
1000
1001 if ( c_mode )
1002 {
1003 if ( p_mode == TRMODE_T )
1004 {
1005 InsertTraceNodewithMode(tr, tr->n-1, TRACE_LEFT_CHILD, best_j - best_d + 1, best_j, best_v, c_mode);
1006 }
1007 tr_wedge_splitter(cm, dsq, L, tr, best_v, z, best_j - best_d + 1, best_j, (c_mode == TRMODE_J), (c_mode == TRMODE_L), (c_mode == TRMODE_R));
1008 }
1009 else /* parent and child both empty */
1010 cm_Die("Danger, danger! p_mode = %d c_mode = %d\n",p_mode,c_mode);
1011
1012 return best_sc;
1013 }
1014
1015 /* Function: tr_v_splitter()
1016 * Author: DLK
1017 *
1018 * Purpose: 'V problem' - closely based on v_splitter()
1019 *
1020 * Args:
1021 *
1022 * Returns; none
1023 */
1024 void
tr_v_splitter(CM_t * cm,ESL_DSQ * dsq,int L,Parsetree_t * tr,int r,int z,int i0,int i1,int j1,int j0,int useEL,int r_allow_J,int r_allow_L,int r_allow_R,int z_allow_J,int z_allow_L,int z_allow_R)1025 tr_v_splitter(CM_t *cm, ESL_DSQ *dsq, int L, Parsetree_t *tr, int r, int z, int i0, int i1,
1026 int j1, int j0, int useEL, int r_allow_J, int r_allow_L, int r_allow_R,
1027 int z_allow_J, int z_allow_L, int z_allow_R)
1028 {
1029 AlphaMats_t *alpha;
1030 BetaMats_t *beta;
1031 float sc, best_sc;
1032 int v, w, y;
1033 int best_v, best_i, best_j;
1034 int midnode;
1035 int p_mode, c_mode;
1036 int jp, ip;
1037 float b_sc;
1038 int b_mode, b_v, b_i, b_j;
1039
1040 /* Recommend a special handler for the fully marginal cases (linear alg.)*/
1041 /*
1042 if ( force_LM)
1043 {
1044
1045 }
1046 else if ( force_RM )
1047 {
1048
1049 }
1050 */
1051
1052 /* Special case: solve without splitting for small problems and boundary conditions */
1053 if (cm->ndidx[z] == cm->ndidx[r] + 1 || r == z ||
1054 5*vinsideT_size(cm, r, z, i0, i1, j1, j0) < RAMLIMIT)
1055 {
1056 tr_vinsideT(cm, dsq, L, tr, r, z, i0, i1, j1, j0, useEL,
1057 r_allow_J, r_allow_L, r_allow_R, z_allow_J, z_allow_L, z_allow_R);
1058 return;
1059 }
1060
1061 alpha = malloc(sizeof(AlphaMats_t));
1062 beta = malloc(sizeof(BetaMats_t));
1063
1064 /* Find split set */
1065 midnode = cm->ndidx[r] + ((cm->ndidx[z] - cm->ndidx[r])/2);
1066 w = cm->nodemap[midnode];
1067 y = cm->cfirst[w] - 1;
1068
1069 /* Calculate alphas and betas */
1070 b_sc = tr_vinside(cm, dsq, L, w, z, i0, i1, j1, j0, useEL, BE_EFFICIENT, (r == 0),
1071 z_allow_J, r_allow_L, r_allow_R, z_allow_J, z_allow_L, z_allow_R,
1072 NULL, alpha, NULL, NULL, NULL, &b_mode, &b_v, &b_i, &b_j);
1073 if (r != 0) b_sc = IMPOSSIBLE;
1074 tr_voutside(cm, dsq, L, r, y, i0, i1, j1, j0, useEL, BE_EFFICIENT,
1075 r_allow_J, r_allow_L, r_allow_R, z_allow_J, z_allow_L, z_allow_R,
1076 NULL, beta, NULL, NULL);
1077
1078 best_sc = IMPOSSIBLE;
1079
1080 /* check local begin */
1081 if (b_sc > best_sc)
1082 {
1083 best_sc = b_sc;
1084 best_v = b_v;
1085 best_i = b_i;
1086 best_j = b_j;
1087 p_mode = 0; c_mode = b_mode;
1088 }
1089
1090 /* Find our best split */
1091 for (v = w; v <= y; v++)
1092 {
1093 for (ip = 0; ip <= i1-i0; ip++)
1094 {
1095 for (jp = 0; jp <= j0-j1; jp++)
1096 {
1097 if ( z_allow_J )
1098 if ( (sc = alpha->J[v][jp][ip] + beta->J[v][jp][ip]) > best_sc )
1099 {
1100 best_sc = sc;
1101 best_v = v;
1102 best_i = ip + i0;
1103 best_j = jp + j1;
1104 p_mode = 3; c_mode = 3;
1105 }
1106 if ( z_allow_J && r_allow_L )
1107 if ( (sc = alpha->J[v][jp][ip] + beta->L[v][ip]) > best_sc )
1108 {
1109 best_sc = sc;
1110 best_v = v;
1111 best_i = ip + i0;
1112 best_j = jp + j1;
1113 p_mode = 2; c_mode = 3;
1114 }
1115 if ( r_allow_L )
1116 if ( (sc = alpha->L[v][jp][ip] + beta->L[v][ip]) > best_sc )
1117 {
1118 best_sc = sc;
1119 best_v = v;
1120 best_i = ip + i0;
1121 best_j = jp + j1;
1122 p_mode = 2; c_mode = 2;
1123 }
1124 if ( z_allow_J && r_allow_R )
1125 if ( (sc = alpha->J[v][jp][ip] + beta->R[v][jp]) > best_sc )
1126 {
1127 best_sc = sc;
1128 best_v = v;
1129 best_i = ip + i0;
1130 best_j = jp + j1;
1131 p_mode = 1; c_mode = 3;
1132 }
1133 if ( r_allow_R )
1134 if ( (sc = alpha->R[v][jp][ip] + beta->R[v][jp]) > best_sc )
1135 {
1136 best_sc = sc;
1137 best_v = v;
1138 best_i = ip + i0;
1139 best_j = jp + j1;
1140 p_mode = 1; c_mode = 1;
1141 }
1142 }
1143 }
1144 }
1145
1146 /* check EL */
1147 if ( useEL )
1148 {
1149 for (ip = 0; ip <= i1-i0; ip++)
1150 {
1151 for (jp = 0; jp <= j0-j1; jp++)
1152 {
1153 if ( (sc = beta->J[cm->M][jp][ip]) > best_sc )
1154 {
1155 best_sc = sc;
1156 best_v = cm->M;
1157 best_i = ip + i0;
1158 best_j = jp + j1;
1159 p_mode = 3; c_mode = 0;
1160 }
1161 }
1162 }
1163 }
1164
1165 /* Free memory */
1166 free_vji_matrix(alpha->J, cm->M, j1, j0);
1167 free_vji_matrix(alpha->L, cm->M, j1, j0);
1168 free_vji_matrix(alpha->R, cm->M, j1, j0);
1169 free_vji_matrix( beta->J, cm->M, j1, j0);
1170 free(beta->L[0]); free(beta->L);
1171 free(beta->R[0]); free(beta->R);
1172 free(alpha);
1173 free(beta);
1174
1175 /* Interpret and subdivide */
1176 if ( p_mode )
1177 {
1178 if ( c_mode )
1179 {
1180 if (best_v == -1) cm_Die("5Superbad: passing z = -1!\n");
1181 tr_v_splitter(cm, dsq, L, tr, r, best_v, i0, best_i, best_j, j0,
1182 FALSE, r_allow_J, r_allow_L, r_allow_R, (c_mode == TRMODE_J), (c_mode == TRMODE_L), (c_mode == TRMODE_R));
1183 if (z == -1) cm_Die("6Superbad: passing z = -1!\n");
1184 tr_v_splitter(cm, dsq, L, tr, best_v, z, best_i, i1, j1, best_j,
1185 useEL, (c_mode == TRMODE_J), (c_mode == TRMODE_L), (c_mode == TRMODE_R), z_allow_J, z_allow_L, z_allow_R);
1186 }
1187 else
1188 {
1189 tr_v_splitter(cm, dsq, L, tr, r, w, i0, best_i, best_j, j0,
1190 TRUE, r_allow_J, r_allow_L, r_allow_R, TRUE, FALSE, FALSE);
1191 }
1192 }
1193 else
1194 {
1195 if (best_v != z)
1196 {
1197 InsertTraceNodewithMode(tr, tr->n-1, TRACE_LEFT_CHILD, best_i, best_j, best_v, c_mode);
1198 }
1199 if (z == -1) cm_Die("7Superbad: passing z = -1!\n");
1200 tr_v_splitter(cm, dsq, L, tr, best_v, z, best_i, i1, j1, best_j,
1201 useEL, (c_mode == TRMODE_J), (c_mode == TRMODE_L), (c_mode == TRMODE_R), z_allow_J, z_allow_L, z_allow_R);
1202 }
1203
1204 return;
1205 }
1206
1207 /* Legacy wrapper */
1208 float
trinside(CM_t * cm,ESL_DSQ * dsq,int L,int vroot,int vend,int i0,int j0,int do_full,void **** ret_shadow,void **** ret_L_shadow,void **** ret_R_shadow,void **** ret_T_shadow,void **** ret_Lmode_shadow,void **** ret_Rmode_shadow,int * ret_mode,int * ret_v,int * ret_i,int * ret_j)1209 trinside (CM_t *cm, ESL_DSQ *dsq, int L, int vroot, int vend, int i0, int j0, int do_full,
1210 void ****ret_shadow, void ****ret_L_shadow, void ****ret_R_shadow,
1211 void ****ret_T_shadow, void ****ret_Lmode_shadow, void ****ret_Rmode_shadow,
1212 int *ret_mode, int *ret_v, int *ret_i, int *ret_j)
1213 {
1214 float sc;
1215 ShadowMats_t shadow;
1216
1217 shadow.J = *ret_shadow;
1218 shadow.L = *ret_L_shadow;
1219 shadow.R = *ret_R_shadow;
1220 shadow.T = *ret_T_shadow;
1221 shadow.Lmode = *ret_Lmode_shadow;
1222 shadow.Rmode = *ret_Rmode_shadow;
1223
1224 sc = tr_inside(cm, dsq, L, vroot, vend, i0, j0, do_full,
1225 TRUE, TRUE, TRUE, TRUE, FALSE,
1226 NULL, NULL, NULL, NULL, &shadow,
1227 ret_mode, ret_v, ret_i, ret_j);
1228 *ret_shadow = shadow.J;
1229 *ret_L_shadow = shadow.L;
1230 *ret_R_shadow = shadow.R;
1231 *ret_T_shadow = shadow.T;
1232 *ret_Lmode_shadow = shadow.Lmode;
1233 *ret_Rmode_shadow = shadow.Rmode;
1234
1235 return sc;
1236 }
1237
1238 /* Function: tr_inside()
1239 * Author: DLK
1240 *
1241 * Purpose: inside phase of CYK on truncated sequence
1242 * based on inside()
1243 *
1244 * Score matrices and deckpool are only managed
1245 * within this func, not available to caller
1246 *
1247 * Args: cm - the covariance model
1248 * dsq - the sequence, 1..L
1249 * L - length of the sequence
1250 * r - root of subgraph to align to target subseq (usually 0, the model's root)
1251 * z - last state of the subgraph
1252 * i0 - start of target subsequence (usually 1, beginning of dsq)
1253 * j0 - end of target subsequence (usually L, end of dsq)
1254 * do_full - if TRUE, save all decks rather than re-using
1255 *
1256 * Returns: Score of the optimal alignment
1257 */
1258 float
tr_inside(CM_t * cm,ESL_DSQ * dsq,int L,int vroot,int vend,int i0,int j0,int do_full,int allow_begin,int r_allow_J,int r_allow_L,int r_allow_R,int lenCORREX,AlphaMats_t * arg_alpha,AlphaMats_t * ret_alpha,struct deckpool_s * dpool,struct deckpool_s ** ret_dpool,ShadowMats_t * ret_shadow,int * ret_mode,int * ret_v,int * ret_i,int * ret_j)1259 tr_inside(CM_t *cm, ESL_DSQ *dsq, int L, int vroot, int vend, int i0, int j0, int do_full,
1260 int allow_begin, int r_allow_J, int r_allow_L, int r_allow_R, int lenCORREX,
1261 AlphaMats_t *arg_alpha, AlphaMats_t *ret_alpha,
1262 struct deckpool_s *dpool, struct deckpool_s **ret_dpool,
1263 ShadowMats_t *ret_shadow, int *ret_mode, int *ret_v, int *ret_i, int *ret_j)
1264 {
1265 float **end;
1266 int nends;
1267 int *touch;
1268 int v,y,z;
1269 int j,d,i,k;
1270 float sc;
1271 int yoffset;
1272 int W;
1273 int jp;
1274 void ***shadow;
1275 void ***L_shadow;
1276 void ***R_shadow;
1277 void ***T_shadow;
1278 int ***Lmode_shadow;
1279 int ***Rmode_shadow;
1280 int **kshad;
1281 char **yshad;
1282 int r_v, r_i, r_j, r_mode;
1283 float r_sc;
1284 float p1, p2, psc;
1285
1286 float ***alpha;
1287 float ***L_alpha;
1288 float ***R_alpha;
1289 float ***T_alpha;
1290
1291 if ( arg_alpha == NULL )
1292 {
1293 alpha = NULL;
1294 L_alpha = NULL;
1295 R_alpha = NULL;
1296 T_alpha = NULL;
1297 }
1298 else
1299 {
1300 alpha = arg_alpha->J;
1301 L_alpha = arg_alpha->L;
1302 R_alpha = arg_alpha->R;
1303 T_alpha = arg_alpha->T;
1304 }
1305
1306 /*Initialization */
1307 r_v = -1;
1308 r_i = i0;
1309 r_j = j0;
1310 r_mode = 3;
1311 r_sc = IMPOSSIBLE;
1312 W = j0-i0+1;
1313 p1 = (float) L / ((float) L + 2.);
1314 p2 = sreLOG2(p1);
1315 p1 = 2*sreLOG2(1.-p1);
1316
1317 /* Make a deckpool */
1318 if ( dpool == NULL ) dpool = deckpool_create();
1319 if (! deckpool_pop(dpool, &end) )
1320 { end = alloc_vjd_deck(L, i0, j0); }
1321 nends = CMSubtreeCountStatetype(cm, vroot, E_st);
1322 for ( jp=0; jp<=W; jp++ )
1323 {
1324 j = i0+jp-1;
1325 end[j][0] = 0.0;
1326 for ( d=1; d<=jp; d++ ) { end[j][d] = IMPOSSIBLE; }
1327 }
1328
1329 /* Create score matrices */
1330 if ( alpha == NULL )
1331 {
1332 alpha = malloc(sizeof(float **) * (cm->M+1));
1333 for ( v=0; v<=cm->M; v++ ) { alpha[v] = NULL; }
1334 }
1335 if ( L_alpha == NULL )
1336 {
1337 L_alpha = malloc(sizeof(float **) * (cm->M+1));
1338 for ( v=0; v<=cm->M; v++ ) { L_alpha[v] = NULL; }
1339 }
1340 if ( R_alpha == NULL )
1341 {
1342 R_alpha = malloc(sizeof(float **) * (cm->M+1));
1343 for ( v=0; v<=cm->M; v++ ) { R_alpha[v] = NULL; }
1344 }
1345 if ( T_alpha == NULL )
1346 {
1347 T_alpha = malloc(sizeof(float **) * (cm->M+1));
1348 for ( v=0; v<=cm->M; v++ ) { T_alpha[v] = NULL; }
1349 }
1350
1351 touch = malloc(sizeof(int) *cm->M);
1352 for ( v=0; v<vroot; v++ ) { touch[v] = 0; }
1353 for ( v=vroot; v<=vend; v++ ) { touch[v] = cm->pnum[v]; }
1354 for ( v=vend+1; v<cm->M; v++ ) {touch[v] = 0; }
1355
1356 /* Create shadow matrices */
1357 if ( ret_shadow != NULL )
1358 {
1359 shadow = (void ***) malloc(sizeof(void **) * cm->M);
1360 for ( v=0; v<cm->M; v++ ) { shadow[v] = NULL; }
1361
1362 L_shadow = (void ***) malloc(sizeof(void **) * cm->M);
1363 for ( v=0; v<cm->M; v++ ) { L_shadow[v] = NULL; }
1364
1365 R_shadow = (void ***) malloc(sizeof(void **) * cm->M);
1366 for ( v=0; v<cm->M; v++ ) { R_shadow[v] = NULL; }
1367
1368 T_shadow = (void ***) malloc(sizeof(void **) * cm->M);
1369 for ( v=0; v<cm->M; v++ ) { T_shadow[v] = NULL; }
1370
1371 Lmode_shadow = (int ***) malloc(sizeof(int **) * cm->M);
1372 for ( v=0; v<cm->M; v++ ) { Lmode_shadow[v] = NULL; }
1373
1374 Rmode_shadow = (int ***) malloc(sizeof(int **) * cm->M);
1375 for ( v=0; v<cm->M; v++ ) { Rmode_shadow[v] = NULL; }
1376 }
1377
1378 /* Main recursion */
1379 for ( v = vend; v >= vroot; v-- )
1380 {
1381 if ( cm->sttype[v] == E_st )
1382 {
1383 alpha[v] = end;
1384 L_alpha[v] = end;
1385 R_alpha[v] = end;
1386 continue;
1387 }
1388 /* Assign alpha decks */
1389 if (! deckpool_pop(dpool, &(alpha[v])) )
1390 { alpha[v] = alloc_vjd_deck(L, i0, j0); }
1391 if (! deckpool_pop(dpool, &(L_alpha[v])) )
1392 { L_alpha[v] = alloc_vjd_deck(L, i0, j0); }
1393 if (! deckpool_pop(dpool, &(R_alpha[v])) )
1394 { R_alpha[v] = alloc_vjd_deck(L, i0, j0); }
1395 if ( (cm->sttype[v] == B_st) && (! deckpool_pop(dpool, &(T_alpha[v])) ) )
1396 { T_alpha[v] = alloc_vjd_deck(L, i0, j0); }
1397
1398 /* Assign shadow decks */
1399 if ( ret_shadow != NULL )
1400 {
1401 if ( cm->sttype[v] == B_st )
1402 {
1403 kshad = alloc_vjd_kshadow_deck(L, i0, j0);
1404 shadow[v] = (void **) kshad;
1405 }
1406 else
1407 {
1408 yshad = alloc_vjd_yshadow_deck(L, i0, j0);
1409 shadow[v] = (void **) yshad;
1410 }
1411
1412 if ( cm->sttype[v] == B_st )
1413 {
1414 kshad = alloc_vjd_kshadow_deck(L, i0, j0);
1415 L_shadow[v] = (void **) kshad;
1416 }
1417 else
1418 {
1419 yshad = alloc_vjd_yshadow_deck(L, i0, j0);
1420 L_shadow[v] = (void **) yshad;
1421 }
1422
1423 if ( cm->sttype[v] == B_st )
1424 {
1425 kshad = alloc_vjd_kshadow_deck(L, i0, j0);
1426 R_shadow[v] = (void **) kshad;
1427 }
1428 else
1429 {
1430 yshad = alloc_vjd_yshadow_deck(L, i0, j0);
1431 R_shadow[v] = (void **) yshad;
1432 }
1433
1434 if ( cm->sttype[v] == B_st )
1435 {
1436 kshad = alloc_vjd_kshadow_deck(L, i0, j0);
1437 T_shadow[v] = (void **) kshad;
1438 }
1439
1440 kshad = alloc_vjd_kshadow_deck(L, i0, j0);
1441 Lmode_shadow[v] = (int **) kshad;
1442
1443 kshad = alloc_vjd_kshadow_deck(L, i0, j0);
1444 Rmode_shadow[v] = (int **) kshad;
1445 }
1446
1447 if ( cm->sttype[v] == D_st || cm->sttype[v] == S_st )
1448 {
1449 for ( jp=0; jp<=W; jp++ )
1450 {
1451 j = i0-1+jp;
1452 for ( d=0; d<=jp; d++ )
1453 {
1454 y = cm->cfirst[v];
1455 alpha[v][j][d] = cm->endsc[v] + (cm->el_selfsc * (d));
1456 /*
1457 if (cm->sttype[v] == S_st) alpha[v][j][d] = cm->endsc[v] + (cm->el_selfsc * (d));
1458 else alpha[v][j][d] = IMPOSSIBLE;
1459 */
1460 L_alpha[v][j][d] = IMPOSSIBLE;
1461 R_alpha[v][j][d] = IMPOSSIBLE;
1462 if ( ret_shadow != NULL )
1463 {
1464 ((char **) shadow[v])[j][d] = USED_EL;
1465 /* Set USED_EL to prevent a traceback bug */
1466 ((char **)L_shadow[v])[j][d] = USED_EL;
1467 ((char **)R_shadow[v])[j][d] = USED_EL;
1468 }
1469
1470 for ( yoffset=0; yoffset<cm->cnum[v]; yoffset++ )
1471 {
1472 if ( (sc = alpha[y+yoffset][j][d] + cm->tsc[v][yoffset]) > alpha[v][j][d] )
1473 {
1474 alpha[v][j][d] = sc;
1475 if ( ret_shadow != NULL ) ((char **)shadow[v])[j][d] = yoffset;
1476 }
1477 if ( (sc = L_alpha[y+yoffset][j][d] + cm->tsc[v][yoffset]) > L_alpha[v][j][d] )
1478 {
1479 L_alpha[v][j][d] = sc;
1480 if ( ret_shadow != NULL ) ((char **)L_shadow[v])[j][d] = yoffset;
1481 if ( ret_shadow != NULL ) Lmode_shadow[v][j][d] = 2;
1482 }
1483 if ( (sc = R_alpha[y+yoffset][j][d] + cm->tsc[v][yoffset]) > R_alpha[v][j][d] )
1484 {
1485 R_alpha[v][j][d] = sc;
1486 if ( ret_shadow != NULL ) ((char **)R_shadow[v])[j][d] = yoffset;
1487 if ( ret_shadow != NULL ) Rmode_shadow[v][j][d] = 1;
1488 }
1489 }
1490
1491 if ( d == 0 )
1492 {
1493 L_alpha[v][j][d] = IMPOSSIBLE;
1494 R_alpha[v][j][d] = IMPOSSIBLE;
1495 }
1496
1497 if ( alpha[v][j][d] < IMPOSSIBLE ) { alpha[v][j][d] = IMPOSSIBLE; }
1498 if ( L_alpha[v][j][d] < IMPOSSIBLE ) { L_alpha[v][j][d] = IMPOSSIBLE; }
1499 if ( R_alpha[v][j][d] < IMPOSSIBLE ) { R_alpha[v][j][d] = IMPOSSIBLE; }
1500 }
1501 }
1502 }
1503 else if ( cm->sttype[v] == B_st )
1504 {
1505 for ( jp=0; jp<=W; jp++ )
1506 {
1507 j = i0-1+jp;
1508 for ( d=0; d<=jp; d++ )
1509 {
1510 int allow_L_exit = 0;
1511 int allow_R_exit = 0;
1512 int allow_J_exit = 0;
1513 y = cm->cfirst[v];
1514 z = cm->cnum[v];
1515
1516 alpha[v][j][d] = alpha[y][j][d] + alpha[z][j][0];
1517 L_alpha[v][j][d] = L_alpha[y][j][d] ;
1518 R_alpha[v][j][d] = + R_alpha[z][j][d];
1519 if ( ret_shadow != NULL )
1520 {
1521 ((int **) shadow[v])[j][d] = 0;
1522 ((int **)L_shadow[v])[j][d] = 0;
1523 ((int **)R_shadow[v])[j][d] = d;
1524 Lmode_shadow[v][j][d] = 2;
1525 Rmode_shadow[v][j][d] = 1;
1526 }
1527
1528 if ( (sc = alpha[y][j][d] ) > L_alpha[v][j][d] )
1529 {
1530 L_alpha[v][j][d] = sc;
1531 if ( ret_shadow != NULL ) { ((int **)L_shadow[v])[j][d] = 0; }
1532 if ( ret_shadow != NULL ) { Lmode_shadow[v][j][d] = 3; }
1533 }
1534
1535 for ( k=1; k<=d; k++ )
1536 {
1537 if ( (sc = alpha[y][j-k][d-k] + alpha[z][j][k]) > alpha[v][j][d] )
1538 {
1539 alpha[v][j][d] = sc;
1540 if (ret_shadow != NULL ) { ((int **)shadow[v])[j][d] = k; }
1541 if (k == d) allow_J_exit = 0;
1542 else allow_J_exit = 1;
1543 }
1544 if ( (sc = alpha[y][j-k][d-k] + L_alpha[z][j][k]) > L_alpha[v][j][d] )
1545 {
1546 L_alpha[v][j][d] = sc;
1547 if ( ret_shadow != NULL ) { ((int **)L_shadow[v])[j][d] = k; }
1548 if ( ret_shadow != NULL ) { Lmode_shadow[v][j][d] = 3; }
1549 allow_L_exit = 1;
1550 }
1551 }
1552 if ( (sc = alpha[z][j][d]) > R_alpha[v][j][d] )
1553 {
1554 R_alpha[v][j][d] = sc;
1555 if ( ret_shadow != NULL ) { ((int **)R_shadow[v])[j][d] = d; }
1556 if ( ret_shadow != NULL ) { Rmode_shadow[v][j][d] = 3; }
1557 }
1558 for ( k=0; k< d; k++ )
1559 {
1560 if ( (sc = R_alpha[y][j-k][d-k] + alpha[z][j][k]) > R_alpha[v][j][d] )
1561 {
1562 R_alpha[v][j][d] = sc;
1563 if ( ret_shadow != NULL ) { ((int **)R_shadow[v])[j][d] = k; }
1564 if ( ret_shadow != NULL ) { Rmode_shadow[v][j][d] = 3; }
1565 allow_R_exit = 1;
1566 }
1567 }
1568
1569 if (d == 0) {
1570 L_alpha[v][j][d] = IMPOSSIBLE;
1571 R_alpha[v][j][d] = IMPOSSIBLE;
1572 }
1573
1574 if (d >= 2) {
1575 T_alpha[v][j][d] = R_alpha[y][j-1][d-1] + L_alpha[z][j][1];
1576 if ( ret_shadow != NULL ) { ((int **)T_shadow[v])[j][d] = 1; }
1577 for ( k=2; k<d; k++ )
1578 {
1579 if ( (sc = R_alpha[y][j-k][d-k] + L_alpha[z][j][k]) > T_alpha[v][j][d] )
1580 {
1581 T_alpha[v][j][d] = sc;
1582 if ( ret_shadow != NULL) { ((int **)T_shadow[v])[j][d] = k; }
1583 }
1584 }
1585 }
1586 else {
1587 T_alpha[v][j][d] = IMPOSSIBLE;
1588 }
1589
1590 if ( alpha[v][j][d] < IMPOSSIBLE ) { alpha[v][j][d] = IMPOSSIBLE; }
1591 if ( L_alpha[v][j][d] < IMPOSSIBLE ) { L_alpha[v][j][d] = IMPOSSIBLE; }
1592 if ( R_alpha[v][j][d] < IMPOSSIBLE ) { R_alpha[v][j][d] = IMPOSSIBLE; }
1593
1594 if ( allow_begin )
1595 {
1596 /* Shouldn't allow exit from marginal B if one of the children is NULL, sinee that is covered by the */
1597 /* root of the other child, and we haven't added anything above the bifurcation */
1598 if (!lenCORREX)
1599 {
1600 if (( alpha[v][j][d] > r_sc) && (allow_J_exit) )
1601 { r_mode = 3; r_v = v; r_j = j; r_i = j-d+1; r_sc = alpha[v][j][d]; }
1602 if ((L_alpha[v][j][d] > r_sc) && (allow_L_exit) )
1603 { r_mode = 2; r_v = v; r_j = j; r_i = j-d+1; r_sc = L_alpha[v][j][d]; }
1604 if ((R_alpha[v][j][d] > r_sc) && (allow_R_exit) )
1605 { r_mode = 1; r_v = v; r_j = j; r_i = j-d+1; r_sc = R_alpha[v][j][d]; }
1606 if ( T_alpha[v][j][d] > r_sc )
1607 { r_mode = 0; r_v = v; r_j = j; r_i = j-d+1; r_sc = T_alpha[v][j][d]; }
1608 }
1609 else
1610 {
1611 psc = p1 + (L-d)*p2;
1612 if (( alpha[v][j][d] + psc > r_sc) && (allow_J_exit) )
1613 { r_mode = 3; r_v = v; r_j = j; r_i = j-d+1; r_sc = alpha[v][j][d] + psc; }
1614 if ((L_alpha[v][j][d] + psc > r_sc) && (allow_L_exit) )
1615 { r_mode = 2; r_v = v; r_j = j; r_i = j-d+1; r_sc = L_alpha[v][j][d] + psc; }
1616 if ((R_alpha[v][j][d] + psc > r_sc) && (allow_R_exit) )
1617 { r_mode = 1; r_v = v; r_j = j; r_i = j-d+1; r_sc = R_alpha[v][j][d] + psc; }
1618 if ( T_alpha[v][j][d] + psc > r_sc )
1619 { r_mode = 0; r_v = v; r_j = j; r_i = j-d+1; r_sc = T_alpha[v][j][d] + psc; }
1620 }
1621 }
1622 }
1623 }
1624 }
1625 else if ( cm->sttype[v] == MP_st )
1626 {
1627 for ( jp=0; jp<=W; jp++ )
1628 {
1629 j = i0-1+jp;
1630 y = cm->cfirst[v];
1631 alpha[v][j][0] = IMPOSSIBLE;
1632 L_alpha[v][j][0] = IMPOSSIBLE;
1633 R_alpha[v][j][0] = IMPOSSIBLE;
1634 if ( jp > 0 ) {
1635 alpha[v][j][1] = IMPOSSIBLE;
1636 L_alpha[v][j][1] = cm->lmesc[v][(int) dsq[j]];
1637 R_alpha[v][j][1] = cm->rmesc[v][(int) dsq[j]];
1638 if ( ret_shadow != NULL ) { ((char **)L_shadow[v])[j][1] = USED_EL; }
1639 if ( ret_shadow != NULL ) { ((char **)R_shadow[v])[j][1] = USED_EL; }
1640 }
1641 for ( d=2; d<=jp; d++)
1642 {
1643 alpha[v][j][d] = cm->endsc[v] + (cm->el_selfsc * (d-2));
1644 L_alpha[v][j][d] = IMPOSSIBLE;
1645 R_alpha[v][j][d] = IMPOSSIBLE;
1646 if ( ret_shadow != NULL ) { ((char **)shadow[v])[j][d] = USED_EL; }
1647
1648 for ( yoffset=0; yoffset<cm->cnum[v]; yoffset++ )
1649 {
1650 if ( (sc = alpha[y+yoffset][j-1][d-2] + cm->tsc[v][yoffset]) > alpha[v][j][d] )
1651 {
1652 alpha[v][j][d] = sc;
1653 if ( ret_shadow != NULL ) { ((char **)shadow[v])[j][d] = yoffset; }
1654 }
1655
1656 if ( (sc = alpha[y+yoffset][j][d-1] + cm->tsc[v][yoffset]) > L_alpha[v][j][d])
1657 {
1658 L_alpha[v][j][d] = sc;
1659 if ( ret_shadow != NULL ) { ((char **)L_shadow[v])[j][d] = yoffset; }
1660 if ( ret_shadow != NULL ) { Lmode_shadow[v][j][d] = 3; }
1661 }
1662 if ( (sc = L_alpha[y+yoffset][j][d-1] + cm->tsc[v][yoffset]) > L_alpha[v][j][d])
1663 {
1664 L_alpha[v][j][d] = sc;
1665 if ( ret_shadow != NULL ) { ((char **)L_shadow[v])[j][d] = yoffset; }
1666 if ( ret_shadow != NULL ) { Lmode_shadow[v][j][d] = 2; }
1667 }
1668
1669 if ( (sc = alpha[y+yoffset][j-1][d-1] + cm->tsc[v][yoffset]) > R_alpha[v][j][d])
1670 {
1671 R_alpha[v][j][d] = sc;
1672 if ( ret_shadow != NULL ) { ((char **)R_shadow[v])[j][d] = yoffset; }
1673 if ( ret_shadow != NULL ) { Rmode_shadow[v][j][d] = 3; }
1674 }
1675 if ( (sc = R_alpha[y+yoffset][j-1][d-1] + cm->tsc[v][yoffset]) > R_alpha[v][j][d])
1676 {
1677 R_alpha[v][j][d] = sc;
1678 if ( ret_shadow != NULL ) { ((char **)R_shadow[v])[j][d] = yoffset; }
1679 if ( ret_shadow != NULL ) { Rmode_shadow[v][j][d] = 1; }
1680 }
1681 }
1682
1683 i = j-d+1;
1684 if ( dsq[i] < cm->abc->K && dsq[j] < cm->abc->K )
1685 { alpha[v][j][d] += cm->esc[v][(int) (dsq[i]*cm->abc->K+dsq[j])]; }
1686 else
1687 { alpha[v][j][d] += DegeneratePairScore(cm->abc, cm->esc[v], dsq[i], dsq[j]); }
1688 { L_alpha[v][j][d] += cm->lmesc[v][(int) dsq[i]]; }
1689 { R_alpha[v][j][d] += cm->rmesc[v][(int) dsq[j]]; }
1690
1691 if ( alpha[v][j][d] < IMPOSSIBLE ) { alpha[v][j][d] = IMPOSSIBLE; }
1692 if ( L_alpha[v][j][d] < IMPOSSIBLE ) { L_alpha[v][j][d] = IMPOSSIBLE; }
1693 if ( R_alpha[v][j][d] < IMPOSSIBLE ) { R_alpha[v][j][d] = IMPOSSIBLE; }
1694 }
1695
1696 for ( d = 1; d <= jp; d++ )
1697 {
1698 if ( allow_begin )
1699 {
1700 if (!lenCORREX)
1701 {
1702 if ( alpha[v][j][d] > r_sc ) { r_mode = 3; r_v = v; r_j = j; r_i = j-d+1; r_sc = alpha[v][j][d]; }
1703 if ( L_alpha[v][j][d] > r_sc ) { r_mode = 2; r_v = v; r_j = j; r_i = j-d+1; r_sc = L_alpha[v][j][d]; }
1704 if ( R_alpha[v][j][d] > r_sc ) { r_mode = 1; r_v = v; r_j = j; r_i = j-d+1; r_sc = R_alpha[v][j][d]; }
1705 }
1706 else
1707 {
1708 psc = p1 + (L-d)*p2;
1709 if ( alpha[v][j][d] + psc > r_sc ) { r_mode = 3; r_v = v; r_j = j; r_i = j-d+1; r_sc = alpha[v][j][d] + psc; }
1710 if ( L_alpha[v][j][d] + psc > r_sc ) { r_mode = 2; r_v = v; r_j = j; r_i = j-d+1; r_sc = L_alpha[v][j][d] + psc; }
1711 if ( R_alpha[v][j][d] + psc > r_sc ) { r_mode = 1; r_v = v; r_j = j; r_i = j-d+1; r_sc = R_alpha[v][j][d] + psc; }
1712 }
1713 }
1714 }
1715 }
1716 }
1717 else if ( cm->sttype[v] == IL_st || cm->sttype[v] == ML_st )
1718 {
1719 for ( jp = 0; jp <= W; jp++ )
1720 {
1721 j = i0-1+jp;
1722 y = cm->cfirst[v];
1723
1724 alpha[v][j][0] = IMPOSSIBLE;
1725 L_alpha[v][j][0] = IMPOSSIBLE;
1726 R_alpha[v][j][0] = IMPOSSIBLE;
1727
1728 for ( d = 1; d <= jp; d++ )
1729 {
1730 alpha[v][j][d] = cm->endsc[v] + (cm->el_selfsc * (d-1));
1731 if (d == 1) L_alpha[v][j][d] = 0.0;
1732 else L_alpha[v][j][d] = IMPOSSIBLE;
1733 R_alpha[v][j][d] = IMPOSSIBLE;
1734 if ( ret_shadow != NULL ) { ((char **) shadow[v])[j][d] = USED_EL; }
1735 if ( ret_shadow != NULL ) { ((char **)L_shadow[v])[j][d] = USED_EL; }
1736 if ( ret_shadow != NULL ) { Lmode_shadow[v][j][d] = 3; }
1737
1738 for ( yoffset = 0; yoffset < cm->cnum[v]; yoffset++ )
1739 {
1740 if ( (sc = alpha[y+yoffset][j][d-1] + cm->tsc[v][yoffset]) > alpha[v][j][d] )
1741 {
1742 alpha[v][j][d] = sc;
1743 if ( ret_shadow != NULL ) { ((char **)shadow[v])[j][d] = yoffset; }
1744 }
1745
1746 if ( d > 1 )
1747 if ( (sc = L_alpha[y+yoffset][j][d-1] + cm->tsc[v][yoffset]) > L_alpha[v][j][d] )
1748 {
1749 L_alpha[v][j][d] = sc;
1750 if ( ret_shadow != NULL ) { ((char **)L_shadow[v])[j][d] = yoffset; }
1751 if ( ret_shadow != NULL ) { Lmode_shadow[v][j][d] = 2; }
1752 }
1753 }
1754
1755 /* we need a separate 'for(yoffset...' loop for the R matrix
1756 * because it depends on a fully calculated J matrix cell
1757 * alpha[v][j][d] in some cases (when v is an IL, and yoffset is 0)
1758 */
1759 for ( yoffset = 0; yoffset < cm->cnum[v]; yoffset++ )
1760 {
1761 if ( (sc = alpha[y+yoffset][j][d] + cm->tsc[v][yoffset]) > R_alpha[v][j][d] )
1762 {
1763 R_alpha[v][j][d] = sc;
1764 if ( ret_shadow != NULL ) { ((char **)R_shadow[v])[j][d] = yoffset; }
1765 if ( ret_shadow != NULL ) { Rmode_shadow[v][j][d] = 3; }
1766 }
1767
1768 if ( (sc = R_alpha[y+yoffset][j][d] + cm->tsc[v][yoffset]) > R_alpha[v][j][d] )
1769 {
1770 R_alpha[v][j][d] = sc;
1771 if ( ret_shadow != NULL ) { ((char **)R_shadow[v])[j][d] = yoffset; }
1772 if ( ret_shadow != NULL ) { Rmode_shadow[v][j][d] = 1; }
1773 }
1774 }
1775 #if 0
1776 for ( yoffset = 0; yoffset < cm->cnum[v]; yoffset++ )
1777 {
1778 if ( (sc = alpha[y+yoffset][j][d-1] + cm->tsc[v][yoffset]) > alpha[v][j][d] )
1779 {
1780 alpha[v][j][d] = sc;
1781 if ( ret_shadow != NULL ) { ((char **)shadow[v])[j][d] = yoffset; }
1782 }
1783
1784 if ( d > 1 )
1785 if ( (sc = L_alpha[y+yoffset][j][d-1] + cm->tsc[v][yoffset]) > L_alpha[v][j][d] )
1786 {
1787 L_alpha[v][j][d] = sc;
1788 if ( ret_shadow != NULL ) { ((char **)L_shadow[v])[j][d] = yoffset; }
1789 if ( ret_shadow != NULL ) { Lmode_shadow[v][j][d] = 2; }
1790 }
1791
1792 /* EPN, Wed Aug 17 09:34:13 2011
1793 * I think this is a minor bug:
1794 * If y == v and yoffset is 0, then
1795 * 1. alpha[y+yoffset][j][d] == alpha[v][j][d] This
1796 * is a problem because alpha[v][j][d] hasn't
1797 * been fully calculated yet, see line about 12
1798 * lines up: 'alpha[v][j][d] = sc;', it could be
1799 * changed for next yoffset value. This means
1800 * that the if statement below: 'if ( (sc =
1801 * alpha[y+yoffset][j][d] + cm->tsc[v][yoffset])
1802 * > R_alpha[v][j][d] )' is using a
1803 * not-yet-finished alpha[y+yoffset][j][d]
1804 * (remember this equals alpha[v][j][d]) value,
1805 * which is bad. The way to fix this is to
1806 * use separate 'for(yoffset...' loops for each
1807 * alpha and R_alpha.
1808 *
1809 * 2. R_alpha[y+yoffset][j][d] == R_alpha[v][j][d]
1810 * I think this is okay because R_alpha[v][j][d] is
1811 * init'ed to IMPOSSIBLE and transitions are always
1812 * <= 0, so R_alpha[v][j][d] will stay as IMPOSSIBLE
1813 * below.
1814 */
1815 if ( (sc = alpha[y+yoffset][j][d] + cm->tsc[v][yoffset]) > R_alpha[v][j][d] )
1816 {
1817 R_alpha[v][j][d] = sc;
1818 if ( ret_shadow != NULL ) { ((char **)R_shadow[v])[j][d] = yoffset; }
1819 if ( ret_shadow != NULL ) { Rmode_shadow[v][j][d] = 3; }
1820 }
1821
1822 if ( (sc = R_alpha[y+yoffset][j][d] + cm->tsc[v][yoffset]) > R_alpha[v][j][d] )
1823 {
1824 R_alpha[v][j][d] = sc;
1825 if ( ret_shadow != NULL ) { ((char **)R_shadow[v])[j][d] = yoffset; }
1826 if ( ret_shadow != NULL ) { Rmode_shadow[v][j][d] = 1; }
1827 }
1828 }
1829 #endif
1830 i = j-d+1;
1831 if ( dsq[i] < cm->abc->K )
1832 {
1833 alpha[v][j][d] += cm->esc[v][(int) dsq[i]];
1834 L_alpha[v][j][d] += cm->esc[v][(int) dsq[i]];
1835 }
1836 else
1837 {
1838 alpha[v][j][d] += esl_abc_FAvgScore(cm->abc, dsq[i], cm->esc[v]);
1839 L_alpha[v][j][d] += esl_abc_FAvgScore(cm->abc, dsq[i], cm->esc[v]);
1840 }
1841
1842 if ( alpha[v][j][d] < IMPOSSIBLE ) { alpha[v][j][d] = IMPOSSIBLE; }
1843 if ( L_alpha[v][j][d] < IMPOSSIBLE ) { L_alpha[v][j][d] = IMPOSSIBLE; }
1844 if ( R_alpha[v][j][d] < IMPOSSIBLE ) { R_alpha[v][j][d] = IMPOSSIBLE; }
1845 }
1846
1847 for ( d = 1; d <= jp; d++ )
1848 {
1849 if ( cm->sttype[v] == ML_st && allow_begin )
1850 {
1851 if (!lenCORREX)
1852 {
1853 if ( alpha[v][j][d] > r_sc ) { r_mode = 3; r_v = v; r_j = j; r_i = j-d+1; r_sc = alpha[v][j][d]; }
1854 if ( L_alpha[v][j][d] > r_sc ) { r_mode = 2; r_v = v; r_j = j; r_i = j-d+1; r_sc = L_alpha[v][j][d]; }
1855 }
1856 else
1857 {
1858 psc = p1 + (L-d)*p2;
1859 if ( alpha[v][j][d] + psc > r_sc ) { r_mode = 3; r_v = v; r_j = j; r_i = j-d+1; r_sc = alpha[v][j][d] + psc; }
1860 if ( L_alpha[v][j][d] + psc > r_sc ) { r_mode = 2; r_v = v; r_j = j; r_i = j-d+1; r_sc = L_alpha[v][j][d] + psc; }
1861 }
1862 }
1863 }
1864 }
1865 }
1866 else if ( cm->sttype[v] == IR_st || cm->sttype[v] == MR_st )
1867 {
1868 for ( jp = 0; jp <= W; jp++ )
1869 {
1870 j = i0-1+jp;
1871 y = cm->cfirst[v];
1872
1873 alpha[v][j][0] = IMPOSSIBLE;
1874 L_alpha[v][j][0] = IMPOSSIBLE;
1875 R_alpha[v][j][0] = IMPOSSIBLE;
1876
1877 for ( d = 1; d <= jp; d++ )
1878 {
1879 alpha[v][j][d] = cm->endsc[v] + (cm->el_selfsc * (d-1));
1880 L_alpha[v][j][d] = IMPOSSIBLE;
1881 if (d == 1) R_alpha[v][j][d] = 0.0;
1882 else R_alpha[v][j][d] = IMPOSSIBLE;
1883 if ( ret_shadow != NULL ) { ((char **) shadow[v])[j][d] = USED_EL; }
1884 if ( ret_shadow != NULL ) { ((char **)R_shadow[v])[j][d] = USED_EL; }
1885 if ( ret_shadow != NULL ) { Rmode_shadow[v][j][d] = 3; }
1886
1887 for ( yoffset = 0; yoffset < cm->cnum[v]; yoffset++ )
1888 {
1889 if ( (sc = alpha[y+yoffset][j-1][d-1] + cm->tsc[v][yoffset]) > alpha[v][j][d] )
1890 {
1891 alpha[v][j][d] = sc;
1892 if ( ret_shadow != NULL ) { ((char **)shadow[v])[j][d] = yoffset; }
1893 }
1894 if ( d > 1 )
1895 if ( (sc = R_alpha[y+yoffset][j-1][d-1] + cm->tsc[v][yoffset]) > R_alpha[v][j][d] )
1896 {
1897 R_alpha[v][j][d] = sc;
1898 if ( ret_shadow != NULL ) { ((char **)R_shadow[v])[j][d] = yoffset; }
1899 if ( ret_shadow != NULL ) { Rmode_shadow[v][j][d] = 1; }
1900 }
1901 }
1902
1903 /* we need a separate 'for(yoffset...' loop for the L matrix
1904 * because it depends on a fully calculated J matrix cell
1905 * alpha[v][j][d] in some cases (when v is an IR, and yoffset is 0)
1906 */
1907 for ( yoffset = 0; yoffset < cm->cnum[v]; yoffset++ )
1908 {
1909 if ( (sc = alpha[y+yoffset][j][d] + cm->tsc[v][yoffset]) > L_alpha[v][j][d] )
1910 {
1911 L_alpha[v][j][d] = sc;
1912 if ( ret_shadow != NULL ) { ((char **)L_shadow[v])[j][d] = yoffset; }
1913 if ( ret_shadow != NULL ) { Lmode_shadow[v][j][d] = 3; }
1914 }
1915
1916 if ( (sc = L_alpha[y+yoffset][j][d] + cm->tsc[v][yoffset]) > L_alpha[v][j][d] )
1917 {
1918 L_alpha[v][j][d] = sc;
1919 if ( ret_shadow != NULL ) { ((char **)L_shadow[v])[j][d] = yoffset; }
1920 if ( ret_shadow != NULL ) { Lmode_shadow[v][j][d] = 2; }
1921 }
1922 }
1923
1924 #if 0
1925 /* EPN, Wed Aug 17 09:34:13 2011 I think this is a
1926 * minor bug, analogous to the one in the same spot
1927 * above for IL, ML states, see that comment for more
1928 * info. (The bug is actually only for self-looping
1929 * IL/IR states.)
1930 */
1931
1932 for ( yoffset = 0; yoffset < cm->cnum[v]; yoffset++ )
1933 {
1934 if ( (sc = alpha[y+yoffset][j-1][d-1] + cm->tsc[v][yoffset]) > alpha[v][j][d] )
1935 {
1936 alpha[v][j][d] = sc;
1937 if ( ret_shadow != NULL ) { ((char **)shadow[v])[j][d] = yoffset; }
1938 }
1939
1940 if ( (sc = alpha[y+yoffset][j][d] + cm->tsc[v][yoffset]) > L_alpha[v][j][d] )
1941 {
1942 L_alpha[v][j][d] = sc;
1943 if ( ret_shadow != NULL ) { ((char **)L_shadow[v])[j][d] = yoffset; }
1944 if ( ret_shadow != NULL ) { Lmode_shadow[v][j][d] = 3; }
1945 }
1946
1947 if ( (sc = L_alpha[y+yoffset][j][d] + cm->tsc[v][yoffset]) > L_alpha[v][j][d] )
1948 {
1949 L_alpha[v][j][d] = sc;
1950 if ( ret_shadow != NULL ) { ((char **)L_shadow[v])[j][d] = yoffset; }
1951 if ( ret_shadow != NULL ) { Lmode_shadow[v][j][d] = 2; }
1952 }
1953
1954 if ( d > 1 )
1955 if ( (sc = R_alpha[y+yoffset][j-1][d-1] + cm->tsc[v][yoffset]) > R_alpha[v][j][d] )
1956 {
1957 R_alpha[v][j][d] = sc;
1958 if ( ret_shadow != NULL ) { ((char **)R_shadow[v])[j][d] = yoffset; }
1959 if ( ret_shadow != NULL ) { Rmode_shadow[v][j][d] = 1; }
1960 }
1961 }
1962 #endif
1963
1964 if ( dsq[j] < cm->abc->K )
1965 {
1966 alpha[v][j][d] += cm->esc[v][(int) dsq[j]];
1967 R_alpha[v][j][d] += cm->esc[v][(int) dsq[j]];
1968 }
1969 else
1970 {
1971 alpha[v][j][d] += esl_abc_FAvgScore(cm->abc, dsq[j], cm->esc[v]);
1972 R_alpha[v][j][d] += esl_abc_FAvgScore(cm->abc, dsq[j], cm->esc[v]);
1973 }
1974
1975 if ( alpha[v][j][d] < IMPOSSIBLE ) { alpha[v][j][d] = IMPOSSIBLE; }
1976 if ( L_alpha[v][j][d] < IMPOSSIBLE ) { L_alpha[v][j][d] = IMPOSSIBLE; }
1977 if ( R_alpha[v][j][d] < IMPOSSIBLE ) { R_alpha[v][j][d] = IMPOSSIBLE; }
1978 }
1979
1980 for ( d = 1; d <= jp; d++ )
1981 {
1982 if ( cm->sttype[v] == MR_st && allow_begin )
1983 {
1984 if (!lenCORREX)
1985 {
1986 if ( alpha[v][j][d] > r_sc ) { r_mode = 3; r_v = v; r_j = j; r_i = j-d+1; r_sc = alpha[v][j][d]; }
1987 if ( R_alpha[v][j][d] > r_sc ) { r_mode = 1; r_v = v; r_j = j; r_i = j-d+1; r_sc = R_alpha[v][j][d]; }
1988 }
1989 else
1990 {
1991 psc = p1 + (L-d)*p2;
1992 if ( alpha[v][j][d] + psc > r_sc ) { r_mode = 3; r_v = v; r_j = j; r_i = j-d+1; r_sc = alpha[v][j][d] + psc; }
1993 if ( R_alpha[v][j][d] + psc > r_sc ) { r_mode = 1; r_v = v; r_j = j; r_i = j-d+1; r_sc = R_alpha[v][j][d] + psc; }
1994 }
1995 }
1996 }
1997 }
1998 }
1999 else
2000 {
2001 cm_Die("'Inconceivable!'\n'You keep using that word...'");
2002 }
2003
2004 #if 0
2005 /* TEMP PRINTING */
2006 if(cm->sttype[v] == B_st) {
2007 for (jp = 0; jp <= W; jp++) {
2008 j = i0-1+jp;
2009 if(j > 0) {
2010 for ( d = 0; d <= jp; d++ ) {
2011 printf("D j: %3d v: %3d d: %3d J: %10.4f L: %10.4f R: %10.4f T: %10.4f\n",
2012 j, v, d,
2013 NOT_IMPOSSIBLE( alpha[v][j][d]) ? alpha[v][j][d] : -9999.9,
2014 NOT_IMPOSSIBLE(L_alpha[v][j][d]) ? L_alpha[v][j][d] : -9999.9,
2015 NOT_IMPOSSIBLE(R_alpha[v][j][d]) ? R_alpha[v][j][d] : -9999.9,
2016 NOT_IMPOSSIBLE(T_alpha[v][j][d]) ? T_alpha[v][j][d] : -9999.9);
2017 }
2018 }
2019 }
2020 }
2021 else {
2022 for (jp = 0; jp <= W; jp++) {
2023 j = i0-1+jp;
2024 if(j > 0) {
2025 for ( d = 0; d <= jp; d++ ) {
2026 printf("D j: %3d v: %3d d: %3d J: %10.4f L: %10.4f R: %10.4f T: %10.4f\n",
2027 j, v, d,
2028 NOT_IMPOSSIBLE( alpha[v][j][d]) ? alpha[v][j][d] : -9999.9,
2029 NOT_IMPOSSIBLE(L_alpha[v][j][d]) ? L_alpha[v][j][d] : -9999.9,
2030 NOT_IMPOSSIBLE(R_alpha[v][j][d]) ? R_alpha[v][j][d] : -9999.9,
2031 -9999.9);
2032 }
2033 }
2034 }
2035 }
2036
2037 printf("\n");
2038 /* TEMP PRINTING */
2039 #endif
2040
2041 if ( v == vroot )
2042 {
2043 if (!lenCORREX)
2044 {
2045 if ( alpha[v][j0][W] > r_sc ) { r_mode = 3; r_v = v; r_j = j0; r_i = j0-W+1; r_sc = alpha[v][j0][W]; }
2046 if ( L_alpha[v][j0][W] > r_sc ) { r_mode = 2; r_v = v; r_j = j0; r_i = j0-W+1; r_sc = L_alpha[v][j0][W]; }
2047 if ( R_alpha[v][j0][W] > r_sc ) { r_mode = 1; r_v = v; r_j = j0; r_i = j0-W+1; r_sc = R_alpha[v][j0][W]; }
2048 }
2049 else
2050 {
2051 psc = p1 + (L-W)*p2;
2052 if ( alpha[v][j0][W] + psc > r_sc ) { r_mode = 3; r_v = v; r_j = j0; r_i = j0-W+1; r_sc = alpha[v][j0][W] + psc; }
2053 if ( L_alpha[v][j0][W] + psc > r_sc ) { r_mode = 2; r_v = v; r_j = j0; r_i = j0-W+1; r_sc = L_alpha[v][j0][W] + psc; }
2054 if ( R_alpha[v][j0][W] + psc > r_sc ) { r_mode = 1; r_v = v; r_j = j0; r_i = j0-W+1; r_sc = R_alpha[v][j0][W] + psc; }
2055 }
2056 }
2057
2058 if ( v==0 )
2059 {
2060 alpha[0][j0][W] = r_sc;
2061 L_alpha[0][j0][W] = r_sc;
2062 R_alpha[0][j0][W] = r_sc;
2063 if ( ret_shadow != NULL ) { ((char **) shadow[0])[j0][W] = USED_LOCAL_BEGIN; }
2064 if ( ret_shadow != NULL ) { ((char **)L_shadow[0])[j0][W] = USED_LOCAL_BEGIN; }
2065 if ( ret_shadow != NULL ) { ((char **)R_shadow[0])[j0][W] = USED_LOCAL_BEGIN; }
2066 if ( ret_shadow != NULL ) { Lmode_shadow[0][j0][W] = r_mode; }
2067 if ( ret_shadow != NULL ) { Rmode_shadow[0][j0][W] = r_mode; }
2068 }
2069
2070 if (! do_full)
2071 {
2072 if ( cm->sttype[v] == B_st )
2073 {
2074 y = cm->cfirst[v];
2075 deckpool_push(dpool, alpha[y]); alpha[y] = NULL;
2076 deckpool_push(dpool, L_alpha[y]); L_alpha[y] = NULL;
2077 deckpool_push(dpool, R_alpha[y]); R_alpha[y] = NULL;
2078 z = cm->cnum[v];
2079 deckpool_push(dpool, alpha[z]); alpha[z] = NULL;
2080 deckpool_push(dpool, L_alpha[z]); L_alpha[z] = NULL;
2081 deckpool_push(dpool, R_alpha[z]); R_alpha[z] = NULL;
2082 }
2083 else
2084 {
2085 for ( y = cm->cfirst[v]; y < cm->cfirst[v]+cm->cnum[v]; y++ )
2086 {
2087 touch[y]--;
2088 if ( touch[y] == 0 )
2089 {
2090 if ( cm->sttype[y] == E_st )
2091 {
2092 nends--;
2093 if ( nends == 0 ) { deckpool_push(dpool, end); end = NULL; }
2094 }
2095 else
2096 {
2097 deckpool_push(dpool, alpha[y]);
2098 deckpool_push(dpool, L_alpha[y]);
2099 deckpool_push(dpool, R_alpha[y]);
2100 if ( cm->sttype[y] == B_st ) { deckpool_push(dpool, T_alpha[y]); }
2101 }
2102
2103 alpha[y] = NULL;
2104 L_alpha[y] = NULL;
2105 R_alpha[y] = NULL;
2106 T_alpha[y] = NULL;
2107 }
2108 }
2109 }
2110 }
2111 } /* end loop over all v */
2112
2113 sc = r_sc;
2114 if ( ret_v != NULL ) { *ret_v = r_v; }
2115 if ( ret_i != NULL ) { *ret_i = r_i; }
2116 if ( ret_j != NULL ) { *ret_j = r_j; }
2117 if ( ret_mode != NULL ) { *ret_mode = r_mode; }
2118
2119 /* Free or return score matrices */
2120 if ( ret_alpha == NULL )
2121 {
2122 for ( v = vroot; v <= vend; v++ )
2123 {
2124 if ( alpha[v] != NULL )
2125 {
2126 if ( cm->sttype[v] != E_st )
2127 {
2128 deckpool_push(dpool, alpha[v]); alpha[v] = NULL;
2129 deckpool_push(dpool, L_alpha[v]); L_alpha[v] = NULL;
2130 deckpool_push(dpool, R_alpha[v]); R_alpha[v] = NULL;
2131 if ( T_alpha[v] != NULL )
2132 { deckpool_push(dpool, T_alpha[v]); T_alpha[v] = NULL; }
2133 }
2134 else
2135 { end = alpha[v]; }
2136 }
2137 }
2138 if ( end != NULL) {deckpool_push(dpool, end); end = NULL; }
2139 free( alpha);
2140 free(L_alpha);
2141 free(R_alpha);
2142 free(T_alpha);
2143 }
2144 else
2145 {
2146 ret_alpha->J = alpha;
2147 ret_alpha->L = L_alpha;
2148 ret_alpha->R = R_alpha;
2149 ret_alpha->T = T_alpha;
2150 }
2151
2152 /* Free or return deckpool */
2153 if ( ret_dpool == NULL )
2154 {
2155 while ( deckpool_pop(dpool, &end)) free_vjd_deck(end, i0, j0);
2156 deckpool_free(dpool);
2157 }
2158 else
2159 {
2160 *ret_dpool = dpool;
2161 }
2162
2163 free(touch);
2164 if ( ret_shadow != NULL ) ret_shadow->J = shadow;
2165 if ( ret_shadow != NULL ) ret_shadow->L = L_shadow;
2166 if ( ret_shadow != NULL ) ret_shadow->R = R_shadow;
2167 if ( ret_shadow != NULL ) ret_shadow->T = T_shadow;
2168 if ( ret_shadow != NULL ) ret_shadow->Lmode = (void ***)Lmode_shadow;
2169 if ( ret_shadow != NULL ) ret_shadow->Rmode = (void ***)Rmode_shadow;
2170 return sc;
2171 }
2172
2173 /* Function: tr_outside()
2174 * Author: DLK
2175 *
2176 * Purpose: outside version of truncated CYK run on
2177 * an unbifurcated model segment vroot..vend
2178 * Closely based on outside()
2179 * Args;
2180 *
2181 * Returns: Score of best local hit (not extending to vend)
2182 */
2183 float
tr_outside(CM_t * cm,ESL_DSQ * dsq,int L,int vroot,int vend,int i0,int j0,int do_full,int r_allow_J,int r_allow_L,int r_allow_R,BetaMats_t * arg_beta,BetaMats_t * ret_beta,struct deckpool_s * dpool,struct deckpool_s ** ret_dpool,int * ret_mode,int * ret_v,int * ret_j)2184 tr_outside(CM_t *cm, ESL_DSQ *dsq, int L, int vroot, int vend, int i0, int j0, int do_full,
2185 int r_allow_J, int r_allow_L, int r_allow_R,
2186 BetaMats_t *arg_beta, BetaMats_t *ret_beta,
2187 struct deckpool_s *dpool, struct deckpool_s **ret_dpool,
2188 int *ret_mode, int *ret_v, int *ret_j)
2189 {
2190 int v,y;
2191 int j,d,i;
2192 float sc;
2193 int *touch;
2194 float esc;
2195 int W;
2196 int jp;
2197 int voffset;
2198 int w1, w2;
2199 int allow_begin;
2200
2201 float b_sc;
2202 int b_mode, b_v, b_j;
2203
2204 BetaMats_t *beta;
2205
2206 W = j0 - i0 + 1;
2207 if ( dpool == NULL ) dpool = deckpool_create();
2208
2209 beta = malloc(sizeof(BetaMats_t));
2210 if ( arg_beta == NULL )
2211 {
2212 beta->J = malloc(sizeof(float **) * (cm->M+1));
2213 beta->L = malloc(sizeof(float *) * (cm->M+1));
2214 beta->R = malloc(sizeof(float *) * (cm->M+1));
2215 beta->L[0] = malloc(sizeof(float) * (cm->M+1)*(L+2));
2216 beta->R[0] = malloc(sizeof(float) * (cm->M+1)*(L+2));
2217 for ( v = 0; v < cm->M+1; v++ )
2218 {
2219 beta->J[v] = NULL;
2220 beta->L[v] = beta->L[0] + v*(L+2);
2221 beta->R[v] = beta->R[0] + v*(L+2);
2222 }
2223 }
2224 else
2225 {
2226 beta->J = arg_beta->J;
2227 beta->L = arg_beta->L;
2228 beta->R = arg_beta->R;
2229 for ( v = 0; v < cm->M+1; v++ )
2230 {
2231 beta->J[v] = arg_beta->J[v];
2232 beta->L[v] = arg_beta->L[v];
2233 beta->R[v] = arg_beta->R[v];
2234 }
2235 }
2236
2237 /* Initialize the root deck, and its split set if applicable */
2238 w1 = cm->nodemap[cm->ndidx[vroot]];
2239 if (cm->sttype[vroot] == B_st)
2240 {
2241 w2 = w1;
2242 if (vend != vroot) cm_Die("vroot B but vroot != vend!\n");
2243 }
2244 else
2245 w2 = cm->cfirst[w1] - 1;
2246
2247 for (v = w1; v <= w2; v++)
2248 {
2249 allow_begin = TRUE;
2250 if ( vroot != 0 ) allow_begin = FALSE;
2251 if ( cm->sttype[v] == IL_st ||
2252 cm->sttype[v] == IR_st ||
2253 cm->sttype[v] == S_st ||
2254 cm->sttype[v] == D_st ||
2255 cm->sttype[v] == E_st ) allow_begin = FALSE;
2256 if (! deckpool_pop(dpool, &(beta->J[v])) )
2257 beta->J[v] = alloc_vjd_deck(L, i0, j0);
2258 for (jp = 0; jp <= W; jp++)
2259 {
2260 j = i0 + jp - 1;
2261 for (d = 0; d <= jp; d++)
2262 if ( allow_begin )
2263 beta->J[v][j][d] = 0.0;
2264 else
2265 beta->J[v][j][d] = IMPOSSIBLE;
2266 if ( allow_begin )
2267 {
2268 beta->L[v][j] = 0.0;
2269 beta->R[v][j] = 0.0;
2270 }
2271 else
2272 {
2273 beta->L[v][j] = IMPOSSIBLE;
2274 beta->R[v][j] = IMPOSSIBLE;
2275 }
2276 }
2277 if ( allow_begin )
2278 {
2279 beta->L[v][i0+W] = 0.0;
2280 beta->R[v][i0+W] = 0.0;
2281 }
2282 else
2283 {
2284 beta->L[v][i0+W] = IMPOSSIBLE;
2285 beta->R[v][i0+W] = IMPOSSIBLE;
2286 }
2287 }
2288 beta->J[vroot][j0][W] = 0.0;
2289 if (r_allow_L) beta->L[vroot][i0] = 0.0; else beta->L[vroot][i0] = IMPOSSIBLE;
2290 if (r_allow_R) beta->R[vroot][j0] = 0.0; else beta->R[vroot][j0] = IMPOSSIBLE;
2291
2292 /* Initialize EL */
2293 if (! deckpool_pop(dpool, &(beta->J[cm->M])) )
2294 beta->J[cm->M] = alloc_vjd_deck(L, i0, j0);
2295 for (jp = 0; jp <= W; jp++)
2296 {
2297 j = i0 + jp - 1;
2298 for (d = 0; d <= jp; d++)
2299 beta->J[cm->M][j][d] = IMPOSSIBLE;
2300 beta->L[cm->M][j] = IMPOSSIBLE;
2301 beta->R[cm->M][j] = IMPOSSIBLE;
2302 }
2303
2304 /* deal with vroot->EL */
2305 /* Marginal modes don't transition to EL,
2306 * so beta->L and beta->R remain at their
2307 * initialization values of IMPOSSIBLE */
2308 if (NOT_IMPOSSIBLE(cm->endsc[vroot]))
2309 {
2310 switch (cm->sttype[vroot])
2311 {
2312 case MP_st:
2313 if (W < 2) break;
2314 if (dsq[i0] < cm->abc->K && dsq[j0] < cm->abc->K)
2315 esc = cm->esc[vroot][(int) (dsq[i0]*cm->abc->K+dsq[j0])];
2316 else
2317 esc = DegeneratePairScore(cm->abc, cm->esc[vroot], dsq[i0], dsq[j0]);
2318 beta->J[cm->M][j0-1][W-2] = cm->endsc[vroot] + (cm->el_selfsc * (W-2)) + esc;
2319 if (beta->J[cm->M][j0-1][W-2] < IMPOSSIBLE) beta->J[cm->M][j0-1][W-2] = IMPOSSIBLE;
2320 break;
2321 case ML_st:
2322 case IL_st:
2323 if (W < 1) break;
2324 if (dsq[i0] < cm->abc->K)
2325 esc = cm->esc[vroot][(int) dsq[i0]];
2326 else
2327 esc = esl_abc_FAvgScore(cm->abc, dsq[i0], cm->esc[vroot]);
2328 beta->J[cm->M][j0][W-1] = cm->endsc[vroot] + (cm->el_selfsc * (W-1)) + esc;
2329 if (beta->J[cm->M][j0][W-1] < IMPOSSIBLE) beta->J[cm->M][j0][W-1] = IMPOSSIBLE;
2330 break;
2331 case MR_st:
2332 case IR_st:
2333 if (W < 1) break;
2334 if (dsq[j0] < cm->abc->K)
2335 esc = cm->esc[vroot][(int) dsq[j0]];
2336 else
2337 esc = esl_abc_FAvgScore(cm->abc, dsq[j0], cm->esc[vroot]);
2338 beta->J[cm->M][j0-1][W-1] = cm->endsc[vroot] + (cm->el_selfsc * (W-1)) + esc;
2339 if (beta->J[cm->M][j0-1][W-1] < IMPOSSIBLE) beta->J[cm->M][j0][W-1] = IMPOSSIBLE;
2340 break;
2341 case S_st:
2342 case D_st:
2343 beta->J[cm->M][j0][W] = cm->endsc[vroot] + (cm->el_selfsc * W);
2344 if (beta->J[cm->M][j0][W] < IMPOSSIBLE) beta->J[cm->M][j0][W] = IMPOSSIBLE;
2345 break;
2346 case B_st: /* B_st can't go to EL? */
2347 default:
2348 cm_Die("bogus parent state %d\n",cm->sttype[vroot]);
2349 }
2350 }
2351
2352 /* Initialize touch vector for controlling deck de-allocation */
2353 touch = malloc(sizeof(int) * cm->M);
2354 for (v = 0; v < w1; v++) touch[v] = 0;
2355 for (v = vend+1; v < cm->M; v++) touch[v] = 0;
2356 for (v = w1; v <= vend; v++)
2357 {
2358 if (cm->sttype[v] == B_st) touch[v] = 2;
2359 else touch[v] = cm->cnum[v];
2360 }
2361
2362 b_sc = IMPOSSIBLE;
2363 b_v = -1;
2364 b_j = -1;
2365 b_mode = -1;
2366
2367 /* Main loop through decks */
2368 for (v = w2+1; v <= vend; v++)
2369 {
2370 allow_begin = TRUE;
2371 if ( vroot != 0 ) allow_begin = FALSE;
2372 if ( cm->sttype[v] == IL_st ||
2373 cm->sttype[v] == IR_st ||
2374 cm->sttype[v] == S_st ||
2375 cm->sttype[v] == D_st ||
2376 cm->sttype[v] == E_st ) allow_begin = FALSE;
2377
2378 /* Get a deck */
2379 if (! deckpool_pop(dpool, &(beta->J[v])) )
2380 beta->J[v] = alloc_vjd_deck(L, i0, j0);
2381 for (jp = W; jp >= 0; jp--)
2382 {
2383 j = i0 + jp - 1;
2384 for (d = jp; d >= 0; d--)
2385 {
2386 if ( allow_begin )
2387 beta->J[v][j][d] = 0.0;
2388 else
2389 beta->J[v][j][d] = IMPOSSIBLE;
2390 }
2391 if ( allow_begin )
2392 {
2393 beta->L[v][j] = 0.0;
2394 beta->R[v][j] = 0.0;
2395 }
2396 else
2397 {
2398 beta->L[v][j] = IMPOSSIBLE;
2399 beta->R[v][j] = IMPOSSIBLE;
2400 }
2401 }
2402 beta->L[v][i0+W] = IMPOSSIBLE;
2403
2404 /* mini-recursion for beta->L */
2405 if ( r_allow_L )
2406 for (j = i0; j <= j0+1; j++)
2407 {
2408 for (y = cm->plast[v]; y > cm->plast[v] - cm->pnum[v]; y--)
2409 {
2410 if (y < vroot) continue;
2411 voffset = v - cm->cfirst[y];
2412
2413 switch (cm->sttype[y])
2414 {
2415 case MP_st:
2416 if (j > i0)
2417 {
2418 esc = cm->lmesc[y][(int) dsq[j-1]];
2419 if ( (sc = beta->L[y][j-1] + cm->tsc[y][voffset] + esc) > beta->L[v][j] )
2420 beta->L[v][j] = sc;
2421 }
2422 break;
2423 case ML_st:
2424 case IL_st:
2425 if (j > i0)
2426 {
2427 if (dsq[j-1] < cm->abc->K)
2428 esc = cm->esc[y][(int) dsq[j-1]];
2429 else
2430 esc = esl_abc_FAvgScore(cm->abc, dsq[j-1], cm->esc[v]);
2431 if ( (sc = beta->L[y][j-1] + cm->tsc[y][voffset] + esc) > beta->L[v][j] )
2432 beta->L[v][j] = sc;
2433 }
2434 break;
2435 case MR_st:
2436 case IR_st:
2437 case S_st:
2438 case E_st:
2439 case D_st:
2440 if ( (sc = beta->L[y][j] + cm->tsc[y][voffset]) > beta->L[v][j] )
2441 beta->L[v][j] = sc;
2442 break;
2443 default:
2444 cm_Die("Bogus parent type %d for y = %d, v = %d\n",cm->sttype[y],y,v);
2445 }
2446 }
2447 esc = 0.0;
2448 if ( j <= j0 )
2449 {
2450 if (cm->sttype[v] == MP_st)
2451 {
2452 esc = cm->lmesc[v][(int) dsq[j]];
2453 }
2454 if (cm->sttype[v] == ML_st || cm->sttype[v] == IL_st)
2455 {
2456 if (dsq[j] < cm->abc->K) esc = cm->esc[v][(int) dsq[j]];
2457 else esc = esl_abc_FAvgScore(cm->abc, dsq[j], cm->esc[v]);
2458 }
2459 }
2460 if (beta->L[v][j] + esc > b_sc)
2461 {
2462 b_sc = beta->L[v][j] + esc;
2463 b_v = v;
2464 b_j = j;
2465 b_mode = 2;
2466 }
2467 }
2468
2469 /* mini-recursion for beta->R */
2470 if ( r_allow_R )
2471 for (j = j0; j >= i0-1; j--)
2472 {
2473 for (y = cm->plast[v]; y > cm->plast[v] - cm->pnum[v]; y--)
2474 {
2475 if (y < vroot) continue;
2476 voffset = v - cm->cfirst[y];
2477
2478 switch (cm->sttype[y])
2479 {
2480 case MP_st:
2481 if (j < j0)
2482 {
2483 esc = cm->rmesc[y][(int) dsq[j+1]];
2484 if ( (sc = beta->R[y][j+1] + cm->tsc[y][voffset] + esc) > beta->R[v][j] )
2485 beta->R[v][j] = sc;
2486 }
2487 break;
2488 case MR_st:
2489 case IR_st:
2490 if (j < j0)
2491 {
2492 if (dsq[j+1] < cm->abc->K)
2493 esc = cm->esc[y][(int) dsq[j+1]];
2494 else
2495 esc = esl_abc_FAvgScore(cm->abc, dsq[j+1], cm->esc[y]);
2496 if ( (sc = beta->R[y][j+1] + cm->tsc[y][voffset] + esc) > beta->R[v][j] )
2497 beta->R[v][j] = sc;
2498 }
2499 break;
2500 case ML_st:
2501 case IL_st:
2502 case S_st:
2503 case E_st:
2504 case D_st:
2505 if ( (sc = beta->R[y][j] + cm->tsc[y][voffset]) > beta->R[v][j] )
2506 beta->R[v][j] = sc;
2507 break;
2508 default:
2509 cm_Die("Bogus parent type %d for y = %d, v = %d\n",cm->sttype[y],y,v);
2510 }
2511 }
2512 esc = 0.0;
2513 if ( j >= i0 )
2514 {
2515 if (cm->sttype[v] == MP_st)
2516 {
2517 esc = cm->rmesc[v][(int) dsq[j]];
2518 }
2519 if (cm->sttype[v] == MR_st || cm->sttype[v] == IR_st)
2520 {
2521 if (dsq[j] < cm->abc->K) esc = cm->esc[v][(int) dsq[j]];
2522 else esc = esl_abc_FAvgScore(cm->abc, dsq[j], cm->esc[v]);
2523 }
2524 }
2525 if (beta->R[v][j] + esc > b_sc)
2526 {
2527 b_sc = beta->R[v][j] + esc;
2528 b_v = v;
2529 b_j = j;
2530 b_mode = 1;
2531 }
2532 }
2533
2534 /* main recursion */
2535 for (jp = W; jp >= 0; jp--)
2536 {
2537 j = i0 + jp - 1;
2538 for (d = jp; d >= 0; d--)
2539 {
2540 i = j - d + 1;
2541 for (y = cm->plast[v]; y > cm->plast[v] - cm->pnum[v]; y--)
2542 {
2543 if (y < vroot) continue;
2544 voffset = v - cm->cfirst[y];
2545
2546 switch (cm->sttype[y])
2547 {
2548 case MP_st:
2549 if (j != j0 && d != jp)
2550 {
2551 if (dsq[i-1] < cm->abc->K && dsq[j+1] < cm->abc->K)
2552 esc = cm->esc[y][(int) (dsq[i-1]*cm->abc->K + dsq[j+1])];
2553 else
2554 esc = DegeneratePairScore(cm->abc, cm->esc[y], dsq[i-1], dsq[j+1]);
2555 if ( (sc = beta->J[y][j+1][d+2] + cm->tsc[y][voffset] + esc) > beta->J[v][j][d] )
2556 beta->J[v][j][d] = sc;
2557 if ( (sc = beta->L[y][i-1] + cm->tsc[y][voffset] + esc) > beta->J[v][j][d] )
2558 beta->J[v][j][d] = sc;
2559 if ( (sc = beta->R[y][j+1] + cm->tsc[y][voffset] + esc) > beta->J[v][j][d] )
2560 beta->J[v][j][d] = sc;
2561 }
2562 break;
2563 case ML_st:
2564 case IL_st:
2565 if (d != jp)
2566 {
2567 if (dsq[i-1] < cm->abc->K)
2568 esc = cm->esc[y][(int) dsq[i-1]];
2569 else
2570 esc = esl_abc_FAvgScore(cm->abc, dsq[i-1], cm->esc[y]);
2571 if ( (sc = beta->J[y][j][d+1] + cm->tsc[y][voffset] + esc) > beta->J[v][j][d] )
2572 beta->J[v][j][d] = sc;
2573 if ( (sc = beta->L[y][i-1] + cm->tsc[y][voffset] + esc) > beta->J[v][j][d] )
2574 beta->J[v][j][d] = sc;
2575 if ( (sc = beta->R[y][j] + cm->tsc[y][voffset] + esc) > beta->J[v][j][d] )
2576 beta->J[v][j][d] = sc;
2577 }
2578 break;
2579 case MR_st:
2580 case IR_st:
2581 if (j != j0)
2582 {
2583 if (dsq[j+1] < cm->abc->K)
2584 esc = cm->esc[y][(int) dsq[j+1]];
2585 else
2586 esc = esl_abc_FAvgScore(cm->abc, dsq[j+1], cm->esc[y]);
2587 if ( (sc = beta->J[y][j+1][d+1] + cm->tsc[y][voffset] + esc) > beta->J[v][j][d] )
2588 beta->J[v][j][d] = sc;
2589 if ( (sc = beta->L[y][i] + cm->tsc[y][voffset] + esc) > beta->J[v][j][d] )
2590 beta->J[v][j][d] = sc;
2591 if ( (sc = beta->R[y][j+1] + cm->tsc[y][voffset] + esc) > beta->J[v][j][d] )
2592 beta->J[v][j][d] = sc;
2593 }
2594 break;
2595 case S_st:
2596 case E_st:
2597 case D_st:
2598 if ( (sc = beta->J[y][j][d] + cm->tsc[y][voffset]) > beta->J[v][j][d] )
2599 beta->J[v][j][d] = sc;
2600 break;
2601 default:
2602 cm_Die("Bogus parent type %d for y = %d, v = %d\n",cm->sttype[y],y,v);
2603 } /* End switch over state types */
2604 } /* End loop over parent states - cell done */
2605 } /* End loop over d - row done */
2606 } /* End loop over jp - deck done */
2607
2608 /* v->EL transitions (beta->J only */
2609 if (NOT_IMPOSSIBLE(cm->endsc[v]))
2610 {
2611 for (jp = 0; jp <= W; jp++)
2612 {
2613 j = i0-1+jp;
2614 for (d = 0; d <= jp; d++)
2615 {
2616 i = j-d+1;
2617 switch (cm->sttype[v])
2618 {
2619 case MP_st:
2620 if (j == j0 || d == jp) continue; /* boundary condition */
2621 if (dsq[i-1] < cm->abc->K && dsq[j+1] < cm->abc->K)
2622 esc = cm->esc[v][(int) (dsq[i-1]*cm->abc->K+dsq[j+1])];
2623 else
2624 esc = DegeneratePairScore(cm->abc, cm->esc[v], dsq[i-1], dsq[j+1]);
2625 if ((sc = beta->J[v][j+1][d+2] + cm->endsc[v] + (cm->el_selfsc * d) + esc) > beta->J[cm->M][j][d])
2626 beta->J[cm->M][j][d] = sc;
2627 break;
2628 case ML_st:
2629 case IL_st:
2630 if (d == jp) continue;
2631 if (dsq[i-1] < cm->abc->K)
2632 esc = cm->esc[v][(int) dsq[i-1]];
2633 else
2634 esc = esl_abc_FAvgScore(cm->abc, dsq[i-1], cm->esc[v]);
2635 if ((sc = beta->J[v][j][d+1] + cm->endsc[v] + (cm->el_selfsc * d) + esc) > beta->J[cm->M][j][d])
2636 /*(cm->el_selfsc * (d+1)) + esc) > beta[cm->M][j][d])*/
2637 beta->J[cm->M][j][d] = sc;
2638 break;
2639 case MR_st:
2640 case IR_st:
2641 if (j == j0) continue;
2642 if (dsq[j+1] < cm->abc->K)
2643 esc = cm->esc[v][(int) dsq[j+1]];
2644 else
2645 esc = esl_abc_FAvgScore(cm->abc, dsq[j+1], cm->esc[v]);
2646 if ((sc = beta->J[v][j+1][d+1] + cm->endsc[v] + (cm->el_selfsc * d) + esc) > beta->J[cm->M][j][d])
2647 /*(cm->el_selfsc * (d+1)) + esc) > beta[cm->M][j][d])*/
2648 beta->J[cm->M][j][d] = sc;
2649 break;
2650 case S_st:
2651 case D_st:
2652 case E_st:
2653 if ((sc = beta->J[v][j][d] + cm->endsc[v] + (cm->el_selfsc * d)) > beta->J[cm->M][j][d])
2654 beta->J[cm->M][j][d] = sc;
2655 break;
2656 case B_st:
2657 default: cm_Die("bogus parent state %d\n", cm->sttype[v]);
2658 /* note that although B is a valid vend for a segment we'd do
2659 outside on, B->EL is set to be impossible, by the local alignment
2660 config. There's no point in having a B->EL because B is a nonemitter
2661 (indeed, it would introduce an alignment ambiguity). The same
2662 alignment case is handled by the X->EL transition where X is the
2663 parent consensus state (S, MP, ML, or MR) above the B. Thus,
2664 this code is relying on the NOT_IMPOSSIBLE() test, above,
2665 to make sure the sttype[vend]=B case gets into this switch.
2666 */
2667 } /* end switch over parent state type v */
2668 } /* end inner loop over d */
2669 } /* end outer loop over jp */
2670 } /* end conditional section for dealing w/ v->EL local end transitions */
2671
2672 /* Recycle memory */
2673 if (! do_full)
2674 {
2675 for (y = cm->plast[v]; y > cm->plast[v] - cm->pnum[v]; y--)
2676 {
2677 touch[y]--;
2678 if (touch[y] == 0) { deckpool_push(dpool, beta->J[y]); beta->J[y] = NULL; }
2679 }
2680 }
2681 } /* end loop over decks v */
2682
2683 /* Clean-up */
2684 if (ret_beta == NULL)
2685 {
2686 for (v = w1; v <= vend; v++)
2687 if (beta->J[v] != NULL) { deckpool_push(dpool, beta->J[v]); beta->J[v] = NULL; }
2688 deckpool_push(dpool, beta->J[cm->M]); beta->J[cm->M] = NULL;
2689 free(beta->L[0]); free(beta->L);
2690 free(beta->R[0]); free(beta->R);
2691 }
2692 else
2693 {
2694 ret_beta->J = beta->J;
2695 ret_beta->L = beta->L;
2696 ret_beta->R = beta->R;
2697 }
2698 free(beta);
2699
2700 if (ret_dpool == NULL)
2701 {
2702 float **a;
2703 while (deckpool_pop(dpool,&a)) free_vjd_deck(a, i0, j0);
2704 deckpool_free(dpool);
2705 }
2706 else
2707 {
2708 *ret_dpool = dpool;
2709 }
2710 free(touch);
2711
2712 if (ret_mode != NULL) *ret_mode = b_mode;
2713 if (ret_v != NULL) *ret_v = b_v;
2714 if (ret_j != NULL) *ret_j = b_j;
2715
2716 return b_sc;
2717 }
2718
2719 /* Function: tr_vinside()
2720 * Author: DLK
2721 *
2722 * Purpose: Inside-type tr-CYK for a v-problem
2723 * Closely modeled on vinside() and tr_inside()
2724 * Note use of vji coordinates rather than vjd
2725 * Args:
2726 *
2727 * Returns:
2728 */
2729 float
tr_vinside(CM_t * cm,ESL_DSQ * dsq,int L,int r,int z,int i0,int i1,int j1,int j0,int useEL,int do_full,int allow_begin,int r_allow_J,int r_allow_L,int r_allow_R,int z_allow_J,int z_allow_L,int z_allow_R,AlphaMats_t * arg_alpha,AlphaMats_t * ret_alpha,struct deckpool_s * dpool,struct deckpool_s ** ret_dpool,ShadowMats_t * ret_shadow,int * ret_mode,int * ret_v,int * ret_i,int * ret_j)2730 tr_vinside(CM_t *cm, ESL_DSQ *dsq, int L, int r, int z, int i0, int i1, int j1, int j0,
2731 int useEL, int do_full, int allow_begin,
2732 int r_allow_J, int r_allow_L, int r_allow_R,
2733 int z_allow_J, int z_allow_L, int z_allow_R,
2734 AlphaMats_t *arg_alpha, AlphaMats_t *ret_alpha,
2735 struct deckpool_s *dpool, struct deckpool_s **ret_dpool,
2736 ShadowMats_t *ret_shadow, int *ret_mode, int *ret_v, int *ret_i, int *ret_j)
2737 {
2738 int v,i,j;
2739 int w1, w2;
2740 int jp, ip;
2741 int *touch;
2742 int y, yoffset;
2743 float sc;
2744
2745 float b_sc;
2746 int b_v, b_i, b_j, b_mode;
2747
2748 AlphaMats_t *alpha;
2749 ShadowMats_t *shadow;
2750
2751 /* Initialization */
2752 b_v = -1;
2753 b_i = i0;
2754 b_j = j0;
2755 b_mode = 3;
2756 b_sc = IMPOSSIBLE;
2757
2758 if ( dpool == NULL ) dpool = deckpool_create();
2759
2760 alpha = malloc(sizeof(AlphaMats_t));
2761 shadow = malloc(sizeof(ShadowMats_t));
2762
2763 /* Create and initialize score matrices */
2764 if ( arg_alpha == NULL )
2765 {
2766 alpha->J = NULL; alpha->L = NULL; alpha->R = NULL;
2767 }
2768 else
2769 {
2770 alpha->J = arg_alpha->J; alpha->L = arg_alpha->L; alpha->R = arg_alpha->R;
2771 }
2772
2773 if ( alpha->J == NULL )
2774 {
2775 alpha->J = malloc(sizeof(float **) * (cm->M+1));
2776 for (v = 0; v <= cm->M; v++) { alpha->J[v] = NULL; }
2777 }
2778 if ( alpha->L == NULL )
2779 {
2780 alpha->L = malloc(sizeof(float **) * (cm->M+1));
2781 for (v = 0; v <= cm->M; v++) { alpha->L[v] = NULL; }
2782 }
2783 if ( alpha->R == NULL )
2784 {
2785 alpha->R = malloc(sizeof(float **) * (cm->M+1));
2786 for (v = 0; v <= cm->M; v++) { alpha->R[v] = NULL; }
2787 }
2788
2789 w1 = cm->nodemap[cm->ndidx[z]];
2790 w2 = cm->cfirst[w1]-1;
2791 for (v = w1; v <= w2; v++)
2792 {
2793 if (! deckpool_pop(dpool, &(alpha->J[v])) )
2794 alpha->J[v] = alloc_vji_deck(i0, i1, j1, j0);
2795 if (! deckpool_pop(dpool, &(alpha->L[v])) )
2796 alpha->L[v] = alloc_vji_deck(i0, i1, j1, j0);
2797 if (! deckpool_pop(dpool, &(alpha->R[v])) )
2798 alpha->R[v] = alloc_vji_deck(i0, i1, j1, j0);
2799 for (jp = 0; jp <= j0-j1; jp++)
2800 {
2801 for (ip = 0; ip <= i1-i0; ip++)
2802 {
2803 alpha->J[v][jp][ip] = IMPOSSIBLE;
2804 alpha->L[v][jp][ip] = IMPOSSIBLE;
2805 alpha->R[v][jp][ip] = IMPOSSIBLE;
2806 }
2807 }
2808 }
2809
2810 touch = malloc(sizeof(int) * cm->M);
2811 for (v = 0; v < r; v++) { touch[v] = 0; }
2812 for (v = r; v <= w2; v++) { touch[v] = cm->pnum[v]; }
2813 for (v = w2+1; v < cm->M; v++) { touch[v] = 0; }
2814
2815 /* Create shadow matrices */
2816 if (ret_shadow != NULL)
2817 {
2818 shadow->J = malloc(sizeof(char **) * cm->M);
2819 shadow->L = malloc(sizeof(char **) * cm->M);
2820 shadow->R = malloc(sizeof(char **) * cm->M);
2821 shadow->Lmode = malloc(sizeof(char **) * cm->M);
2822 shadow->Rmode = malloc(sizeof(char **) * cm->M);
2823 for (v = 0; v < cm->M; v++)
2824 {
2825 shadow->J[v] = NULL;
2826 shadow->L[v] = NULL;
2827 shadow->R[v] = NULL;
2828 shadow->Lmode[v] = NULL;
2829 shadow->Rmode[v] = NULL;
2830 }
2831 }
2832
2833 /* Initialize our non-IMPOSSIBLE boundary condition */
2834 /* (Includes an unroll of the main recursion to handle EL) */
2835 ip = i1 - i0;
2836 jp = 0;
2837 if (! useEL)
2838 {
2839 if ( z_allow_J )
2840 alpha->J[z][jp][ip] = 0.0;
2841 if ( z_allow_L )
2842 alpha->L[z][jp][ip] = 0.0;
2843 if ( z_allow_R )
2844 alpha->R[z][jp][ip] = 0.0;
2845 }
2846 else if ( z_allow_J )
2847 {
2848 if (ret_shadow != NULL)
2849 {
2850 shadow->J[z] = (void **) alloc_vji_shadow_deck(i0,i1,j1,j0);
2851 shadow->L[z] = (void **) alloc_vji_shadow_deck(i0,i1,j1,j0);
2852 shadow->R[z] = (void **) alloc_vji_shadow_deck(i0,i1,j1,j0);
2853 shadow->Lmode[z] = (void **) alloc_vji_shadow_deck(i0,i1,j1,j0);
2854 shadow->Rmode[z] = (void **) alloc_vji_shadow_deck(i0,i1,j1,j0);
2855 }
2856
2857 switch (cm->sttype[z])
2858 {
2859 case D_st:
2860 case S_st:
2861 alpha->J[z][jp][ip] = cm->endsc[z] + (cm->el_selfsc * ((jp+j1)-(ip+i0)+1));
2862 if (ret_shadow != NULL) ((char **)shadow->J[z])[jp][ip] = USED_EL;
2863 break;
2864 case MP_st:
2865 if (i0 == i1 || j1 == j0) break;
2866 alpha->J[z][jp+1][ip-1] = cm->endsc[z] + (cm->el_selfsc * ((jp+j1)-(ip+i0)+1));
2867 if (dsq[i1-1] < cm->abc->K && dsq[j1+1] < cm->abc->K)
2868 alpha->J[z][jp+1][ip-1] += cm->esc[z][(int) (dsq[i1-1]*cm->abc->K+dsq[j1+1])];
2869 else
2870 alpha->J[z][jp+1][ip-1] += DegeneratePairScore(cm->abc, cm->esc[z], dsq[i1-1], dsq[j1+1]);
2871 if (ret_shadow != NULL) ((char **)shadow->J[z])[jp+1][ip-1] = USED_EL;
2872 if (alpha->J[z][jp+1][ip-1] < IMPOSSIBLE) alpha->J[z][jp+1][ip-1] = IMPOSSIBLE;
2873 break;
2874 case ML_st:
2875 case IL_st:
2876 if (i0 == i1 ) break;
2877 alpha->J[z][jp][ip-1] = cm->endsc[z] + (cm->el_selfsc * ((jp+j1)-(ip+i0)+1));
2878 if (dsq[i1-1] < cm->abc->K)
2879 alpha->J[z][jp][ip-1] += cm->esc[z][(int) dsq[i1-1]];
2880 else
2881 alpha->J[z][jp][ip-1] += esl_abc_FAvgScore(cm->abc, dsq[i1-1], cm->esc[z]);
2882 if (ret_shadow != NULL) ((char **)shadow->J[z])[jp][ip-1] = USED_EL;
2883 if (alpha->J[z][jp][ip-1] < IMPOSSIBLE) alpha->J[z][jp][ip-1] = IMPOSSIBLE;
2884 break;
2885 case MR_st:
2886 case IR_st:
2887 if (j1 == j0) break;
2888 alpha->J[z][jp+1][ip] = cm->endsc[z] + (cm->el_selfsc * ((jp+j1)-(ip+i0)+1));
2889 if (dsq[j1+1] < cm->abc->K)
2890 alpha->J[z][jp+1][ip] += cm->esc[z][(int) dsq[j1+1]];
2891 else
2892 alpha->J[z][jp+1][ip] += esl_abc_FAvgScore(cm->abc, dsq[j1+1], cm->esc[z]);
2893 if (ret_shadow != NULL) ((char **)shadow->J[z])[jp+1][ip] = USED_EL;
2894 if (alpha->J[z][jp+1][ip] < IMPOSSIBLE) alpha->J[z][jp+1][ip] = IMPOSSIBLE;
2895 break;
2896 /*
2897 default:
2898 cm_Die("Bad input combination in tr_vinside: useEL TRUE, but cm->sttype[z] = %d\n",cm->sttype[z]);
2899 */
2900 }
2901
2902 alpha->L[z][jp][ip] = IMPOSSIBLE;
2903 alpha->R[z][jp][ip] = IMPOSSIBLE;
2904 if (ret_shadow != NULL)
2905 {
2906 /* didn't actually use EL, but this prevents a traceback bug */
2907 ((char **)shadow->L[z])[jp][ip] = USED_EL;
2908 ((char **)shadow->R[z])[jp][ip] = USED_EL;
2909 }
2910 }
2911 else
2912 cm_Die("Bad input combination in tr_vinside: useEL %d z_allow_J %d \n",useEL,z_allow_J);
2913
2914 /* Special case: empty sequence */
2915 if (r == 0)
2916 {
2917 b_v = z; b_i = i1; b_j = j1;
2918 b_sc = IMPOSSIBLE; b_mode = 0;
2919 if (z_allow_J && alpha->J[z][0][i1-i0] > b_sc)
2920 {
2921 b_sc = alpha->J[z][0][i1-i0];
2922 b_mode = 3;
2923 }
2924 if (z_allow_L && alpha->L[z][0][i1-i0] > b_sc)
2925 {
2926 b_sc = alpha->L[z][0][i1-i0];
2927 b_mode = 2;
2928 }
2929 if (z_allow_R && alpha->R[z][0][i1-i0] > b_sc)
2930 {
2931 b_sc = alpha->R[z][0][i1-i0];
2932 b_mode = 1;
2933 }
2934
2935 if (z == 0)
2936 {
2937 /* FIXME */
2938 /* I don't understand what exactly Sean's doing in this block */
2939 cm_Die("Potentially unhandled case!\n");
2940 }
2941 }
2942
2943 /* Main recursion */
2944 for (v = w1-1; v >= r; v--)
2945 {
2946 /* Get decks */
2947 if (! deckpool_pop(dpool, &(alpha->J[v])) )
2948 alpha->J[v] = alloc_vji_deck(i0,i1,j1,j0);
2949 if (! deckpool_pop(dpool, &(alpha->L[v])) )
2950 alpha->L[v] = alloc_vji_deck(i0,i1,j1,j0);
2951 if (! deckpool_pop(dpool, &(alpha->R[v])) )
2952 alpha->R[v] = alloc_vji_deck(i0,i1,j1,j0);
2953
2954 if (ret_shadow != NULL)
2955 {
2956 shadow->J[v] = (void **) alloc_vji_shadow_deck(i0,i1,j1,j0);
2957 shadow->L[v] = (void **) alloc_vji_shadow_deck(i0,i1,j1,j0);
2958 shadow->R[v] = (void **) alloc_vji_shadow_deck(i0,i1,j1,j0);
2959 shadow->Lmode[v] = (void **) alloc_vji_shadow_deck(i0,i1,j1,j0);
2960 shadow->Rmode[v] = (void **) alloc_vji_shadow_deck(i0,i1,j1,j0);
2961 }
2962
2963 /* Full initialization of the deck */
2964 /* This may be wasteful, since it could be folded into the rest
2965 * of the DP */
2966 for (jp = 0; jp <= j0-j1; jp++)
2967 for (ip = i1-i0; ip >= 0; ip--)
2968 {
2969 alpha->J[v][jp][ip] = IMPOSSIBLE;
2970 alpha->L[v][jp][ip] = IMPOSSIBLE;
2971 alpha->R[v][jp][ip] = IMPOSSIBLE;
2972 if (ret_shadow != NULL)
2973 {
2974 /* Didn't really use EL, but trying to eliminate uninitialized values */
2975 ((char **)shadow->J[v])[jp][ip] = USED_EL;
2976 ((char **)shadow->L[v])[jp][ip] = USED_EL;
2977 ((char **)shadow->R[v])[jp][ip] = USED_EL;
2978 ((char **)shadow->Lmode[v])[jp][ip] = 0;
2979 ((char **)shadow->Rmode[v])[jp][ip] = 0;
2980 }
2981 }
2982
2983 /* Double-check problem type */
2984 if (cm->sttype[v] == E_st || cm->sttype[v] == B_st || (cm->sttype[v] == S_st && v > r))
2985 cm_Die("Non-V problem in tr_vinside(); cm->sttype[%d] = %d\n",v,cm->sttype[v]);
2986
2987 if (cm->sttype[v] == D_st || cm->sttype[v] == S_st)
2988 {
2989 for (jp = 0; jp <= j0-j1; jp++)
2990 {
2991 for (ip = i1-i0; ip >= 0; ip--)
2992 {
2993 y = cm->cfirst[v];
2994 if (useEL && NOT_IMPOSSIBLE(cm->endsc[v]) && z_allow_J)
2995 if ( (sc = cm->endsc[v] + (cm->el_selfsc * ((jp+j1)-(ip+i0)+1))) > alpha->J[v][jp][ip])
2996 {
2997 alpha->J[v][jp][ip] = sc;
2998 if (ret_shadow != NULL) ((char **)shadow->J[v])[jp][ip] = USED_EL;
2999 }
3000
3001 for (yoffset = 0; yoffset < cm->cnum[v]; yoffset++)
3002 {
3003 if ( z_allow_J )
3004 if ( (sc = alpha->J[y+yoffset][jp][ip] + cm->tsc[v][yoffset]) > alpha->J[v][jp][ip])
3005 {
3006 alpha->J[v][jp][ip] = sc;
3007 if (ret_shadow != NULL) ((char **)shadow->J[v])[jp][ip] = (char) yoffset;
3008 }
3009 if ( r_allow_L )
3010 if ( (sc = alpha->L[y+yoffset][jp][ip] + cm->tsc[v][yoffset]) > alpha->L[v][jp][ip])
3011 {
3012 alpha->L[v][jp][ip] = sc;
3013 if (ret_shadow != NULL) { ((char **)shadow->L[v])[jp][ip] = (char) yoffset; ((char **)shadow->Lmode[v])[jp][ip] = 2; }
3014 }
3015 if ( r_allow_R )
3016 if ( (sc = alpha->R[y+yoffset][jp][ip] + cm->tsc[v][yoffset]) > alpha->R[v][jp][ip])
3017 {
3018 alpha->R[v][jp][ip] = sc;
3019 if (ret_shadow != NULL) { ((char **)shadow->R[v])[jp][ip] = (char) yoffset; ((char **)shadow->Rmode[v])[jp][ip] = 1; }
3020 }
3021 }
3022
3023 if ( alpha->J[v][jp][ip] < IMPOSSIBLE ) alpha->J[v][jp][ip] = IMPOSSIBLE;
3024 if ( alpha->L[v][jp][ip] < IMPOSSIBLE ) alpha->L[v][jp][ip] = IMPOSSIBLE;
3025 if ( alpha->R[v][jp][ip] < IMPOSSIBLE ) alpha->R[v][jp][ip] = IMPOSSIBLE;
3026 }
3027 }
3028 }
3029 else if (cm->sttype[v] == MP_st)
3030 {
3031 for (jp = 0; jp <= j0-j1; jp++)
3032 {
3033 j = jp+j1;
3034 for (ip = i1-i0; ip >= 0; ip--)
3035 {
3036 i = ip+i0;
3037 y = cm->cfirst[v];
3038
3039 if (useEL && NOT_IMPOSSIBLE(cm->endsc[v]) && z_allow_J && jp > 0 && ip < i1-i0)
3040 if ( (sc = cm->endsc[v] + (cm->el_selfsc * ((jp+j1)-(ip+i0)+1 - 2))) > alpha->J[v][jp][ip] )
3041 {
3042 alpha->J[v][jp][ip] = sc;
3043 if (ret_shadow != NULL) ((char **)shadow->J[v])[jp][ip] = USED_EL;
3044 }
3045
3046 for (yoffset = 0; yoffset < cm->cnum[v]; yoffset++)
3047 {
3048 if (z_allow_J && jp > 0 && ip < i1-i0)
3049 if ( (sc = alpha->J[y+yoffset][jp-1][ip+1] + cm->tsc[v][yoffset]) > alpha->J[v][jp][ip])
3050 {
3051 alpha->J[v][jp][ip] = sc;
3052 if (ret_shadow != NULL) ((char **)shadow->J[v])[jp][ip] = (char) yoffset;
3053 }
3054 if (r_allow_L && ip < i1-i0)
3055 if ( (sc = alpha->J[y+yoffset][jp][ip+1] + cm->tsc[v][yoffset]) > alpha->L[v][jp][ip])
3056 {
3057 alpha->L[v][jp][ip] = sc;
3058 if (ret_shadow != NULL) { ((char **)shadow->L[v])[jp][ip] = (char) yoffset; ((char **)shadow->Lmode[v])[jp][ip] = 3; }
3059 }
3060 if (r_allow_L && ip < i1-i0)
3061 if ( (sc = alpha->L[y+yoffset][jp][ip+1] + cm->tsc[v][yoffset]) > alpha->L[v][jp][ip])
3062 {
3063 alpha->L[v][jp][ip] = sc;
3064 if (ret_shadow != NULL) { ((char **)shadow->L[v])[jp][ip] = (char) yoffset; ((char **)shadow->Lmode[v])[jp][ip] = 2; }
3065 }
3066 if (r_allow_R && jp > 0)
3067 if ( (sc = alpha->J[y+yoffset][jp-1][ip] + cm->tsc[v][yoffset]) > alpha->R[v][jp][ip])
3068 {
3069 alpha->R[v][jp][ip] = sc;
3070 if (ret_shadow != NULL) { ((char **)shadow->R[v])[jp][ip] = (char) yoffset; ((char **)shadow->Rmode[v])[jp][ip] = 3; }
3071 }
3072 if (r_allow_R && jp > 0)
3073 if ( (sc = alpha->R[y+yoffset][jp-1][ip] + cm->tsc[v][yoffset]) > alpha->R[v][jp][ip])
3074 {
3075 alpha->R[v][jp][ip] = sc;
3076 if (ret_shadow != NULL) { ((char **)shadow->R[v])[jp][ip] = (char) yoffset; ((char **)shadow->Rmode[v])[jp][ip] = 1; }
3077 }
3078 }
3079
3080 if (jp > 0 && ip < i1-i0)
3081 {
3082 if (dsq[i] < cm->abc->K && dsq[j] < cm->abc->K)
3083 alpha->J[v][jp][ip] += cm->esc[v][(int) (dsq[i]*cm->abc->K+dsq[j])];
3084 else
3085 alpha->J[v][jp][ip] += DegeneratePairScore(cm->abc, cm->esc[v], dsq[i], dsq[j]);
3086 }
3087 if (ip < i1-i0)
3088 {
3089 alpha->L[v][jp][ip] += cm->lmesc[v][(int) dsq[i]];
3090 }
3091 if (jp > 0)
3092 {
3093 alpha->R[v][jp][ip] += cm->rmesc[v][(int) dsq[j]];
3094 }
3095
3096 if ( alpha->J[v][jp][ip] < IMPOSSIBLE) alpha->J[v][jp][ip] = IMPOSSIBLE;
3097 if ( alpha->L[v][jp][ip] < IMPOSSIBLE) alpha->L[v][jp][ip] = IMPOSSIBLE;
3098 if ( alpha->R[v][jp][ip] < IMPOSSIBLE) alpha->R[v][jp][ip] = IMPOSSIBLE;
3099
3100 if ( allow_begin )
3101 {
3102 if ( r_allow_J )
3103 if ( alpha->J[v][jp][ip] > b_sc ) { b_mode = 3; b_v = v; b_j = j1+jp; b_i = i0+ip; b_sc = alpha->J[v][jp][ip]; }
3104 if ( r_allow_L )
3105 if ( alpha->L[v][jp][ip] > b_sc ) { b_mode = 2; b_v = v; b_j = j1+jp; b_i = i0+ip; b_sc = alpha->L[v][jp][ip]; }
3106 if ( r_allow_R )
3107 if ( alpha->R[v][jp][ip] > b_sc ) { b_mode = 1; b_v = v; b_j = j1+jp; b_i = i0+ip; b_sc = alpha->R[v][jp][ip]; }
3108 }
3109 }
3110 }
3111 }
3112 else if (cm->sttype[v] == ML_st || cm->sttype[v] == IL_st)
3113 {
3114 for (jp = 0; jp <= j0-j1; jp++)
3115 {
3116 for (ip = i1-i0; ip >= 0; ip--)
3117 {
3118 i = i0+ip;
3119 y = cm->cfirst[v];
3120
3121 if (useEL && NOT_IMPOSSIBLE(cm->endsc[v]) && z_allow_J && ip < i1-i0)
3122 if ( (sc = cm->endsc[v] + (cm->el_selfsc * ((jp+j1)-(ip+i0)+1 - 1))) > alpha->J[v][jp][ip] )
3123 {
3124 alpha->J[v][jp][ip] = sc;
3125 if (ret_shadow != NULL) ((char **)shadow->J[v])[jp][ip] = USED_EL;
3126 }
3127
3128 for (yoffset = 0; yoffset < cm->cnum[v]; yoffset++)
3129 {
3130 if (z_allow_J && ip < i1-i0)
3131 if ( (sc = alpha->J[y+yoffset][jp][ip+1] + cm->tsc[v][yoffset]) > alpha->J[v][jp][ip])
3132 {
3133 alpha->J[v][jp][ip] = sc;
3134 if (ret_shadow != NULL) ((char **)shadow->J[v])[jp][ip] = (char) yoffset;
3135 }
3136 if (r_allow_L && ip < i1-i0)
3137 if ( (sc = alpha->L[y+yoffset][jp][ip+1] + cm->tsc[v][yoffset]) > alpha->L[v][jp][ip])
3138 {
3139 alpha->L[v][jp][ip] = sc;
3140 if (ret_shadow != NULL) { ((char **)shadow->L[v])[jp][ip] = (char) yoffset; ((char **)shadow->Lmode[v])[jp][ip] = 2; }
3141 }
3142 if ( r_allow_R )
3143 if ( (sc = alpha->J[y+yoffset][jp][ip] + cm->tsc[v][yoffset]) > alpha->R[v][jp][ip])
3144 {
3145 alpha->R[v][jp][ip] = sc;
3146 if (ret_shadow != NULL) { ((char **)shadow->R[v])[jp][ip] = (char) yoffset; ((char **)shadow->Rmode[v])[jp][ip] = 3; }
3147 }
3148 if ( r_allow_R )
3149 if ( (sc = alpha->R[y+yoffset][jp][ip] + cm->tsc[v][yoffset]) > alpha->R[v][jp][ip])
3150 {
3151 alpha->R[v][jp][ip] = sc;
3152 if (ret_shadow != NULL) { ((char **)shadow->R[v])[jp][ip] = (char) yoffset; ((char **)shadow->Rmode[v])[jp][ip] = 1; }
3153 }
3154 }
3155
3156 if (ip < i1-i0)
3157 {
3158 if (dsq[i] < cm->abc->K)
3159 {
3160 alpha->J[v][jp][ip] += cm->esc[v][(int) dsq[i]];
3161 alpha->L[v][jp][ip] += cm->esc[v][(int) dsq[i]];
3162 }
3163 else
3164 {
3165 alpha->J[v][jp][ip] += esl_abc_FAvgScore(cm->abc, dsq[i], cm->esc[v]);
3166 alpha->L[v][jp][ip] += esl_abc_FAvgScore(cm->abc, dsq[i], cm->esc[v]);
3167 }
3168 }
3169
3170 if ( alpha->J[v][jp][ip] < IMPOSSIBLE) alpha->J[v][jp][ip] = IMPOSSIBLE;
3171 if ( alpha->L[v][jp][ip] < IMPOSSIBLE) alpha->L[v][jp][ip] = IMPOSSIBLE;
3172 if ( alpha->R[v][jp][ip] < IMPOSSIBLE) alpha->R[v][jp][ip] = IMPOSSIBLE;
3173
3174 if ( cm->sttype[v] == ML_st && allow_begin )
3175 {
3176 if ( r_allow_J )
3177 if ( alpha->J[v][jp][ip] > b_sc ) { b_mode = 3; b_v = v; b_j = j1+jp; b_i = i0+ip; b_sc = alpha->J[v][jp][ip]; }
3178 if ( r_allow_L )
3179 if ( alpha->L[v][jp][ip] > b_sc ) { b_mode = 2; b_v = v; b_j = j1+jp; b_i = i0+ip; b_sc = alpha->L[v][jp][ip]; }
3180 }
3181 }
3182 }
3183 }
3184 else if (cm->sttype[v] == MR_st || cm->sttype[v] == IR_st)
3185 {
3186 for (jp = 0; jp <= j0-j1; jp++)
3187 {
3188 j = j1+jp;
3189 for (ip = i1-i0; ip >= 0; ip--)
3190 {
3191 y = cm->cfirst[v];
3192
3193 if (useEL && NOT_IMPOSSIBLE(cm->endsc[v]) && z_allow_J && jp > 0)
3194 if ( (sc = cm->endsc[v] + (cm->el_selfsc * ((j1+jp)-(i0+ip)+1 - 1))) > alpha->J[v][jp][ip] )
3195 {
3196 alpha->J[v][jp][ip] = sc;
3197 if (ret_shadow != NULL) ((char **)shadow->J[v])[jp][ip] = USED_EL;
3198 }
3199
3200 for (yoffset = 0; yoffset < cm->cnum[v]; yoffset++)
3201 {
3202 if (z_allow_J && jp > 0)
3203 if ( (sc = alpha->J[y+yoffset][jp-1][ip] + cm->tsc[v][yoffset]) > alpha->J[v][jp][ip])
3204 {
3205 alpha->J[v][jp][ip] = sc;
3206 if (ret_shadow != NULL) ((char **)shadow->J[v])[jp][ip] = (char) yoffset;
3207 }
3208 if ( r_allow_L )
3209 if ( (sc = alpha->J[y+yoffset][jp][ip] + cm->tsc[v][yoffset]) > alpha->L[v][jp][ip])
3210 {
3211 alpha->L[v][jp][ip] = sc;
3212 if (ret_shadow != NULL) { ((char **)shadow->L[v])[jp][ip] = (char) yoffset; ((char **)shadow->Lmode[v])[jp][ip] = 3; }
3213 }
3214 if ( r_allow_L )
3215 if ( (sc = alpha->L[y+yoffset][jp][ip] + cm->tsc[v][yoffset]) > alpha->L[v][jp][ip])
3216 {
3217 alpha->L[v][jp][ip] = sc;
3218 if (ret_shadow != NULL) { ((char **)shadow->L[v])[jp][ip] = (char) yoffset; ((char **)shadow->Lmode[v])[jp][ip] = 2; }
3219 }
3220 if (r_allow_R && jp > 0)
3221 if ( (sc = alpha->R[y+yoffset][jp-1][ip] + cm->tsc[v][yoffset]) > alpha->R[v][jp][ip])
3222 {
3223 alpha->R[v][jp][ip] = sc;
3224 if (ret_shadow != NULL) { ((char **)shadow->R[v])[jp][ip] = (char) yoffset; ((char **)shadow->Rmode[v])[jp][ip] = 1; }
3225 }
3226 }
3227
3228 if (jp > 0)
3229 {
3230 if (dsq[j] < cm->abc->K)
3231 {
3232 alpha->J[v][jp][ip] += cm->esc[v][(int) dsq[j]];
3233 alpha->R[v][jp][ip] += cm->esc[v][(int) dsq[j]];
3234 }
3235 else
3236 {
3237 alpha->J[v][jp][ip] += esl_abc_FAvgScore(cm->abc, dsq[j], cm->esc[v]);
3238 alpha->R[v][jp][ip] += esl_abc_FAvgScore(cm->abc, dsq[j], cm->esc[v]);
3239 }
3240 }
3241
3242 if ( alpha->J[v][jp][ip] < IMPOSSIBLE) alpha->J[v][jp][ip] = IMPOSSIBLE;
3243 if ( alpha->L[v][jp][ip] < IMPOSSIBLE) alpha->L[v][jp][ip] = IMPOSSIBLE;
3244 if ( alpha->R[v][jp][ip] < IMPOSSIBLE) alpha->R[v][jp][ip] = IMPOSSIBLE;
3245
3246 if ( cm->sttype[v] == MR_st && allow_begin )
3247 {
3248 if ( r_allow_J )
3249 if ( alpha->J[v][jp][ip] > b_sc ) { b_mode = 3; b_v = v; b_j = j1+jp; b_i = i0+ip; b_sc = alpha->J[v][jp][ip]; }
3250 if ( r_allow_R )
3251 if ( alpha->R[v][jp][ip] > b_sc ) { b_mode = 1; b_v = v; b_j = j1+jp; b_i = i0+ip; b_sc = alpha->R[v][jp][ip]; }
3252 }
3253 }
3254 }
3255 }
3256 else
3257 {
3258 cm_Die("There's no way we could have gotten here - should have died before now\n");
3259 }
3260
3261 if (v == r)
3262 {
3263 if ( r_allow_J )
3264 if ( alpha->J[v][j0-j1][0] > b_sc) { b_mode = 3; b_v = v; b_j = j0; b_i = i0; b_sc = alpha->J[v][j0-j1][0]; }
3265 if ( r_allow_L )
3266 if ( alpha->L[v][j0-j1][0] > b_sc) { b_mode = 2; b_v = v; b_j = j0; b_i = i0; b_sc = alpha->L[v][j0-j1][0]; }
3267 if ( r_allow_R )
3268 if ( alpha->R[v][j0-j1][0] > b_sc) { b_mode = 1; b_v = v; b_j = j0; b_i = i0; b_sc = alpha->R[v][j0-j1][0]; }
3269 }
3270
3271 /* If we're at root, give it the best (local) score */
3272 if (v == 0)
3273 {
3274 alpha->J[v][j0-j1][0] = b_sc;
3275 alpha->L[v][j0-j1][0] = b_sc;
3276 alpha->R[v][j0-j1][0] = b_sc;
3277 if (ret_shadow != NULL)
3278 {
3279 ((char **)shadow->J[v])[j0-j1][0] = USED_LOCAL_BEGIN;
3280 ((char **)shadow->L[v])[j0-j1][0] = USED_LOCAL_BEGIN;
3281 ((char **)shadow->R[v])[j0-j1][0] = USED_LOCAL_BEGIN;
3282 ((char **)shadow->Lmode[v])[j0-j1][0] = b_mode;
3283 ((char **)shadow->Rmode[v])[j0-j1][0] = b_mode;
3284 }
3285 }
3286
3287 /* Recycle memory */
3288 if (! do_full)
3289 {
3290 for (y = cm->cfirst[v]; y < cm->cfirst[v]+cm->cnum[v]; y++)
3291 {
3292 touch[y]--;
3293 if (touch[y] == 0)
3294 {
3295 deckpool_push(dpool, alpha->J[y]);
3296 deckpool_push(dpool, alpha->L[y]);
3297 deckpool_push(dpool, alpha->R[y]);
3298 alpha->J[y] = NULL;
3299 alpha->L[y] = NULL;
3300 alpha->R[y] = NULL;
3301 }
3302 }
3303 }
3304 } /* end loop over v */
3305
3306 sc = b_sc;
3307 if (ret_v != NULL ) *ret_v = b_v;
3308 if (ret_i != NULL ) *ret_i = b_i;
3309 if (ret_j != NULL ) *ret_j = b_j;
3310 if (ret_mode != NULL ) *ret_mode = b_mode;
3311
3312 /* Free or return score matrices */
3313 if (ret_alpha == NULL)
3314 {
3315 for (v = r; v <= w2; v++)
3316 {
3317 if (alpha->J[v] != NULL)
3318 {
3319 deckpool_push(dpool, alpha->J[v]);
3320 alpha->J[v] = NULL;
3321 }
3322 if (alpha->L[v] != NULL)
3323 {
3324 deckpool_push(dpool, alpha->L[v]);
3325 alpha->L[v] = NULL;
3326 }
3327 if (alpha->R[v] != NULL)
3328 {
3329 deckpool_push(dpool, alpha->R[v]);
3330 alpha->R[v] = NULL;
3331 }
3332 }
3333 free(alpha->J);
3334 free(alpha->L);
3335 free(alpha->R);
3336 }
3337 else
3338 {
3339 ret_alpha->J = alpha->J;
3340 ret_alpha->L = alpha->L;
3341 ret_alpha->R = alpha->R;
3342 }
3343 free(alpha);
3344
3345 /* Free or return deck pool */
3346 if (ret_dpool == NULL)
3347 {
3348 float **foo;
3349 while (deckpool_pop(dpool, &foo))
3350 free_vji_deck(foo, j1, j0);
3351 deckpool_free(dpool);
3352 }
3353 else
3354 {
3355 *ret_dpool = dpool;
3356 }
3357
3358 free(touch);
3359 if (ret_shadow != NULL)
3360 {
3361 ret_shadow->J = shadow->J;
3362 ret_shadow->L = shadow->L;
3363 ret_shadow->R = shadow->R;
3364 ret_shadow->Lmode = shadow->Lmode;
3365 ret_shadow->Rmode = shadow->Rmode;
3366 }
3367 free(shadow);
3368
3369 return sc;
3370 }
3371
3372 /* Function: tr_voutside()
3373 * Author: DLK
3374 *
3375 * Purpose: Outside direction TrCYK for a v-problem
3376 * Based closely on voutside() and tr_outside()
3377 * Note use of vji instead of vjd coordinates
3378 *
3379 * Args:
3380 *
3381 * Returns:
3382 */
3383 void
tr_voutside(CM_t * cm,ESL_DSQ * dsq,int L,int r,int z,int i0,int i1,int j1,int j0,int useEL,int do_full,int r_allow_J,int r_allow_L,int r_allow_R,int z_allow_J,int z_allow_L,int z_allow_R,BetaMats_t * arg_beta,BetaMats_t * ret_beta,struct deckpool_s * dpool,struct deckpool_s ** ret_dpool)3384 tr_voutside(CM_t *cm, ESL_DSQ *dsq, int L, int r, int z, int i0, int i1, int j1, int j0,
3385 int useEL, int do_full, int r_allow_J, int r_allow_L, int r_allow_R,
3386 int z_allow_J, int z_allow_L, int z_allow_R, BetaMats_t *arg_beta,
3387 BetaMats_t *ret_beta, struct deckpool_s *dpool, struct deckpool_s **ret_dpool)
3388 {
3389 int v,y;
3390 int i,j;
3391 int ip, jp;
3392 float sc, esc;
3393 int voffset;
3394 int *touch;
3395 int allow_begin;
3396
3397 BetaMats_t *beta;
3398
3399 /* Initialization */
3400 if (dpool == NULL) dpool = deckpool_create();
3401
3402 beta = malloc(sizeof(BetaMats_t));
3403 if (arg_beta == NULL)
3404 {
3405 beta->J = malloc(sizeof(float **) * (cm->M+1));
3406 beta->L = malloc(sizeof(float *) * (cm->M+1));
3407 beta->R = malloc(sizeof(float *) * (cm->M+1));
3408 beta->L[0] = malloc(sizeof(float) * (cm->M+1)*(i1-i0+1));
3409 beta->R[0] = malloc(sizeof(float) * (cm->M+1)*(j0-j1+1));
3410 for (v = 0; v < cm->M+1; v++)
3411 {
3412 beta->J[v] = NULL;
3413 beta->L[v] = beta->L[0] + (v * (i1-i0+1));
3414 beta->R[v] = beta->R[0] + (v * (j0-j1+1));
3415 }
3416 }
3417 else
3418 {
3419 beta->J = arg_beta->J;
3420 beta->L = arg_beta->L;
3421 beta->R = arg_beta->R;
3422 }
3423
3424 /* Initialize root deck */
3425 /* outside()/tr_outside() also initialize the root's
3426 split set, if it has one, while voutside() doesn't
3427 I think that this is because in calls to voutside()
3428 (and analagously, tr_voutside() ) we've already
3429 determined that the root state is actually used in
3430 the solution, whereas for the more generic outside
3431 that's not necessarily the case. Not sure, though.
3432 outside()/tr_outside() might not even need to worry
3433 about the split set, but do it anyway (legacy code?) */
3434 if (! deckpool_pop(dpool, &(beta->J[r])) )
3435 beta->J[r] = alloc_vji_deck(i0,i1,j1,j0);
3436 for (jp = 0; jp <= j0-j1; jp++)
3437 {
3438 for (ip = 0; ip <= i1-i0; ip++)
3439 if (r == 0 && r_allow_J )
3440 beta->J[r][jp][ip] = 0.0;
3441 else
3442 beta->J[r][jp][ip] = IMPOSSIBLE;
3443 if ( r == 0 && r_allow_R )
3444 beta->R[r][jp] = 0.0;
3445 else
3446 beta->R[r][jp] = IMPOSSIBLE;
3447 }
3448 for (ip = 0; ip <= i1-i0; ip++)
3449 {
3450 if (r == 0 && r_allow_L )
3451 beta->L[r][ip] = 0.0;
3452 else
3453 beta->L[r][ip] = IMPOSSIBLE;
3454 }
3455 if ( r_allow_J )
3456 beta->J[r][j0-j1][0] = 0.0;
3457 if ( r_allow_L )
3458 beta->L[r][0] = 0.0;
3459 if ( r_allow_R )
3460 beta->R[r][j0-j1] = 0.0;
3461
3462 /* Deal with vroot->EL; marginal modes don't use EL */
3463 if (useEL)
3464 {
3465 if (! deckpool_pop(dpool, &(beta->J[cm->M])) )
3466 beta->J[cm->M] = alloc_vji_deck(i0,i1,j1,j0);
3467 for (jp = 0; jp <= j0-j1; jp++)
3468 for (ip = 0; ip <= i1-i0; ip++)
3469 beta->J[cm->M][jp][ip] = IMPOSSIBLE;
3470 }
3471 if (useEL && NOT_IMPOSSIBLE(cm->endsc[r]))
3472 {
3473 switch(cm->sttype[r])
3474 {
3475 case MP_st:
3476 if (i0 == i1 || j1 == j0) break;
3477 if (dsq[i0] < cm->abc->K && dsq[j0] < cm->abc->K)
3478 esc = cm->esc[r][(int) (dsq[i0]*cm->abc->K+dsq[j0])];
3479 else
3480 esc = DegeneratePairScore(cm->abc, cm->esc[r], dsq[i0], dsq[j0]);
3481 beta->J[cm->M][j0-j1-1][1] = cm->endsc[r] + (cm->el_selfsc * ((j0-1)-(i0+1)+1)) + esc;
3482 if (beta->J[cm->M][j0-j1-1][1] < IMPOSSIBLE) beta->J[cm->M][j0-j1-1][1] = IMPOSSIBLE;
3483 break;
3484 case ML_st:
3485 case IL_st:
3486 if (i0 == i1) break;
3487 if (dsq[i0] < cm->abc->K)
3488 esc = cm->esc[r][(int) dsq[i0]];
3489 else
3490 esc = esl_abc_FAvgScore(cm->abc, dsq[i0], cm->esc[r]);
3491 beta->J[cm->M][j0-j1][1] = cm->endsc[r] + (cm->el_selfsc * ((j0)-(i0+1)+1)) + esc;
3492 if (beta->J[cm->M][j0-j1][1] < IMPOSSIBLE) beta->J[cm->M][j0-j1][1] = IMPOSSIBLE;
3493 break;
3494 case MR_st:
3495 case IR_st:
3496 if (j1 == j0) break;
3497 if (dsq[j0] < cm->abc->K)
3498 esc = cm->esc[r][(int) dsq[j0]];
3499 else
3500 esc = esl_abc_FAvgScore(cm->abc, dsq[j0], cm->esc[r]);
3501 beta->J[cm->M][j0-j1-1][0] = cm->endsc[r] + (cm->el_selfsc * ((j0-1)-(i0)+1)) + esc;
3502 if (beta->J[cm->M][j0-j1-1][0] < IMPOSSIBLE) beta->J[cm->M][j0-j1-1][0] = IMPOSSIBLE;
3503 break;
3504 case S_st:
3505 case D_st:
3506 beta->J[cm->M][j0-j1][0] = cm->endsc[r] + (cm->el_selfsc * ((j0)-(i0)+1));
3507 break;
3508 default:
3509 cm_Die("bogus parent state %d\n",cm->sttype[r]);
3510 }
3511 }
3512
3513 /* Initialize touch vector for controlling deck recycling */
3514 touch = malloc(sizeof(int) * cm->M);
3515 for (v = 0; v < r; v++) touch[v] = 0;
3516 for (v = z+1; v < cm->M; v++) touch[v] = 0;
3517 for (v = r; v <= z; v++)
3518 {
3519 if (cm->sttype[v] == B_st) touch[v] = 2;
3520 else touch[v] = cm->cnum[v];
3521 }
3522
3523 /* Main loop through decks */
3524 for (v = r+1; v <= z; v++)
3525 {
3526 allow_begin = TRUE;
3527 if (r != 0) allow_begin = FALSE;
3528 if ( cm->sttype[v] == IL_st ||
3529 cm->sttype[v] == IR_st ||
3530 cm->sttype[v] == S_st ||
3531 cm->sttype[v] == D_st ||
3532 cm->sttype[v] == E_st ) allow_begin = FALSE;
3533 /* Get a deck */
3534 if (! deckpool_pop(dpool, &(beta->J[v])) )
3535 beta->J[v] = alloc_vji_deck(i0,i1,j1,j0);
3536 for (jp = j0-j1; jp >= 0; jp--)
3537 {
3538 for (ip = 0; ip <= i1-i0; ip++)
3539 {
3540 if (allow_begin && r_allow_J )
3541 beta->J[v][jp][ip] = 0.0;
3542 else
3543 beta->J[v][jp][ip] = IMPOSSIBLE;
3544 }
3545 if (allow_begin && r_allow_R )
3546 beta->R[v][jp] = 0.0;
3547 else
3548 beta->R[v][jp] = IMPOSSIBLE;
3549 }
3550 for (ip = 0; ip <= i1-i0; ip++)
3551 {
3552 if (allow_begin && r_allow_L )
3553 beta->L[v][ip] = 0.0;
3554 else
3555 beta->L[v][ip] = IMPOSSIBLE;
3556 }
3557
3558 /* mini-recursion for beta->L */
3559 if ( r_allow_L )
3560 for (ip = 0; ip <= i1-i0; ip++)
3561 {
3562 i = i0+ip;
3563 for (y = cm->plast[v]; y > cm->plast[v] - cm->pnum[v]; y--)
3564 {
3565 if (y < r) continue;
3566 voffset = v - cm->cfirst[y];
3567
3568 switch (cm->sttype[y])
3569 {
3570 case MP_st:
3571 if (ip > 0)
3572 {
3573 esc = cm->lmesc[y][(int) dsq[i-1]];
3574 if ( (sc = beta->L[y][ip-1] + cm->tsc[y][voffset] + esc) > beta->L[v][ip] )
3575 beta->L[v][ip] = sc;
3576 }
3577 break;
3578 case ML_st:
3579 case IL_st:
3580 if (ip > 0)
3581 {
3582 if (dsq[i-1] < cm->abc->K)
3583 esc = cm->esc[y][(int) dsq[i-1]];
3584 else
3585 esc = esl_abc_FAvgScore(cm->abc, dsq[i-1], cm->esc[y]);
3586 if ( (sc = beta->L[y][ip-1] + cm->tsc[y][voffset] + esc) > beta->L[v][ip] )
3587 beta->L[v][ip] = sc;
3588 }
3589 break;
3590 case MR_st:
3591 case IR_st:
3592 case S_st:
3593 case E_st:
3594 case D_st:
3595 if ( (sc = beta->L[y][ip] + cm->tsc[y][voffset]) > beta->L[v][ip] )
3596 beta->L[v][ip] = sc;
3597 break;
3598 default:
3599 cm_Die("Bogus parent type %d for y = %d, v = %d\n",cm->sttype[y],y,v);
3600 }
3601 }
3602
3603 if (beta->L[v][ip] < IMPOSSIBLE) beta->L[v][ip] = IMPOSSIBLE;
3604 }
3605
3606 /* mini-recursion for beta->R */
3607 if ( r_allow_R )
3608 for (jp = j0-j1; jp >= 0; jp--)
3609 {
3610 j = j1+jp;
3611 for (y = cm->plast[v]; y > cm->plast[v] - cm->pnum[v]; y--)
3612 {
3613 if (y < r) continue;
3614 voffset = v - cm->cfirst[y];
3615
3616 switch (cm->sttype[y])
3617 {
3618 case MP_st:
3619 if (jp < j0-j1)
3620 {
3621 esc = cm->rmesc[y][(int) dsq[j+1]];
3622 if ( (sc = beta->R[y][jp+1] + cm->tsc[y][voffset] + esc) > beta->R[v][jp] )
3623 beta->R[v][jp] = sc;
3624 }
3625 break;
3626 case MR_st:
3627 case IR_st:
3628 if (jp < j0-j1)
3629 {
3630 if (dsq[j+1] < cm->abc->K)
3631 esc = cm->esc[y][(int) dsq[j+1]];
3632 else
3633 esc = esl_abc_FAvgScore(cm->abc, dsq[j+1], cm->esc[y]);
3634 if ( (sc = beta->R[y][jp+1] + cm->tsc[y][voffset] + esc) > beta->R[v][jp] )
3635 beta->R[v][jp] = sc;
3636 }
3637 break;
3638 case ML_st:
3639 case IL_st:
3640 case S_st:
3641 case E_st:
3642 case D_st:
3643 if ( (sc = beta->R[y][jp] + cm->tsc[y][voffset]) > beta->R[v][jp] )
3644 beta->R[v][jp] = sc;
3645 break;
3646 default:
3647 cm_Die("Bogus parent type %d for y = %d, v = %d\n",cm->sttype[y],y,v);
3648 }
3649 }
3650
3651 if (beta->R[v][jp] < IMPOSSIBLE) beta->R[v][jp] = IMPOSSIBLE;
3652 }
3653
3654 /* Main recursion */
3655 if ( z_allow_J )
3656 for (jp = j0-j1; jp >= 0; jp--)
3657 {
3658 j = j1+jp;
3659 for (ip = 0; ip <= i1-i0; ip++)
3660 {
3661 i = i0+ip;
3662 for (y = cm->plast[v]; y > cm->plast[v] - cm->pnum[v]; y--)
3663 {
3664 if (y < r) continue;
3665 voffset = v - cm->cfirst[y];
3666
3667 switch(cm->sttype[y])
3668 {
3669 case MP_st:
3670 if (j != j0 && i != i0)
3671 {
3672 if (dsq[i-1] < cm->abc->K && dsq[j+1] < cm->abc->K)
3673 esc = cm->esc[y][(int) (dsq[i-1]*cm->abc->K+dsq[j+1])];
3674 else
3675 esc = DegeneratePairScore(cm->abc, cm->esc[y], dsq[i-1], dsq[j+1]);
3676 if ( (sc = beta->J[y][jp+1][ip-1] + cm->tsc[y][voffset] + esc) > beta->J[v][jp][ip] )
3677 beta->J[v][jp][ip] = sc;
3678 if ( (sc = beta->L[y][ip-1] + cm->tsc[y][voffset] + esc) > beta->J[v][jp][ip] )
3679 beta->J[v][jp][ip] = sc;
3680 if ( (sc = beta->R[y][jp+1] + cm->tsc[y][voffset] + esc) > beta->J[v][jp][ip] )
3681 beta->J[v][jp][ip] = sc;
3682 }
3683 break;
3684 case ML_st:
3685 case IL_st:
3686 if (i != i0)
3687 {
3688 if (dsq[i-1] < cm->abc->K)
3689 esc = cm->esc[y][(int) dsq[i-1]];
3690 else
3691 esc = esl_abc_FAvgScore(cm->abc, dsq[i-1], cm->esc[y]);
3692 if ( (sc = beta->J[y][jp][ip-1] + cm->tsc[y][voffset] + esc) > beta->J[v][jp][ip] )
3693 beta->J[v][jp][ip] = sc;
3694 if ( (sc = beta->L[y][ip-1] + cm->tsc[y][voffset] + esc) > beta->J[v][jp][ip] )
3695 beta->J[v][jp][ip] = sc;
3696 if ( (sc = beta->R[y][jp] + cm->tsc[y][voffset] + esc) > beta->J[v][jp][ip] )
3697 beta->J[v][jp][ip] = sc;
3698 }
3699 break;
3700 case MR_st:
3701 case IR_st:
3702 if (j != j0)
3703 {
3704 if (dsq[j+1] < cm->abc->K)
3705 esc = cm->esc[y][(int) dsq[j+1]];
3706 else
3707 esc = esl_abc_FAvgScore(cm->abc, dsq[j+1], cm->esc[y]);
3708 if ( (sc = beta->J[y][jp+1][ip] + cm->tsc[y][voffset] + esc) > beta->J[v][jp][ip] )
3709 beta->J[v][jp][ip] = sc;
3710 if ( (sc = beta->L[y][ip] + cm->tsc[y][voffset] + esc) > beta->J[v][jp][ip] )
3711 beta->J[v][jp][ip] = sc;
3712 if ( (sc = beta->R[y][jp+1] + cm->tsc[y][voffset] + esc) > beta->J[v][jp][ip] )
3713 beta->J[v][jp][ip] = sc;
3714 }
3715 break;
3716 case S_st:
3717 case E_st:
3718 case D_st:
3719 if ( (sc = beta->J[y][jp][ip] + cm->tsc[y][voffset]) > beta->J[v][jp][ip] )
3720 beta->J[v][jp][ip] = sc;
3721 break;
3722 default:
3723 cm_Die("Bogus parent type %d for y = %d, v = %d\n",cm->sttype[y],y,v);
3724 }
3725 }
3726
3727 if (beta->J[v][jp][ip] < IMPOSSIBLE) beta->J[v][jp][ip] = IMPOSSIBLE;
3728 }
3729 }
3730
3731 /* v->EL transitions (beta->J only) */
3732 if (useEL && NOT_IMPOSSIBLE(cm->endsc[v]))
3733 {
3734 for (jp = j0-j1; jp >= 0; jp--)
3735 {
3736 j = j1+jp;
3737 for (ip = 0; ip <= i1-i0; ip++)
3738 {
3739 i = i0+ip;
3740
3741 switch (cm->sttype[v])
3742 {
3743 case MP_st:
3744 if (j != j0 && i != i0)
3745 {
3746 if (dsq[i-1] < cm->abc->K && dsq[j+1] < cm->abc->K)
3747 esc = cm->esc[v][(int) (dsq[i-1]*cm->abc->K+dsq[j+1])];
3748 else
3749 esc = DegeneratePairScore(cm->abc, cm->esc[v], dsq[i-1], dsq[j+1]);
3750 if ( (sc = beta->J[v][jp+1][ip-1] + cm->endsc[v] + (cm->el_selfsc* (j-i+1)) + esc) > beta->J[cm->M][jp][ip] )
3751 beta->J[cm->M][jp][ip] = sc;
3752 }
3753 break;
3754 case ML_st:
3755 case IL_st:
3756 if (i != i0)
3757 {
3758 if (dsq[i-1] < cm->abc->K)
3759 esc = cm->esc[v][(int) dsq[i-1]];
3760 else
3761 esc = esl_abc_FAvgScore(cm->abc, dsq[i-1], cm->esc[v]);
3762 if ( (sc = beta->J[v][jp][ip-1] + cm->endsc[v] + (cm->el_selfsc* (j-i+1)) + esc) > beta->J[cm->M][jp][ip] )
3763 beta->J[cm->M][jp][ip] = sc;
3764 }
3765 break;
3766 case MR_st:
3767 case IR_st:
3768 if (j != j0)
3769 {
3770 if (dsq[j+1] < cm->abc->K)
3771 esc = cm->esc[v][(int) dsq[j+1]];
3772 else
3773 esc = esl_abc_FAvgScore(cm->abc, dsq[j+1], cm->esc[v]);
3774 if ( (sc = beta->J[v][jp+1][ip] + cm->endsc[v] + (cm->el_selfsc* (j-i+1)) + esc) > beta->J[cm->M][jp][ip] )
3775 beta->J[cm->M][jp][ip] = sc;
3776 }
3777 break;
3778 case S_st:
3779 case E_st:
3780 case D_st:
3781 if ( (sc = beta->J[v][jp][ip] + cm->endsc[v] + (cm->el_selfsc * (j-1+1)) + esc) > beta->J[cm->M][jp][ip] )
3782 beta->J[cm->M][jp][ip] = sc;
3783 default:
3784 cm_Die("Bogus parent type %d for y = %d, v = %d\n",cm->sttype[y],y,v);
3785 }
3786
3787 if (beta->J[cm->M][jp][ip] < IMPOSSIBLE) beta->J[cm->M][jp][ip] = IMPOSSIBLE;
3788 }
3789 }
3790 }
3791
3792 /* Recycle memory */
3793 if (! do_full)
3794 {
3795 for (y = cm->plast[v]; y > cm->plast[v] - cm->pnum[v]; y--)
3796 {
3797 touch[y]--;
3798 if (touch[y] == 0) { deckpool_push(dpool, beta->J[y]); beta->J[y] = NULL; }
3799 }
3800 }
3801 } /* end loop over decks v */
3802
3803 /* Clean-up */
3804 if (ret_beta == NULL)
3805 {
3806 for (v = r; v <= z; v++)
3807 if (beta->J[v] != NULL) { deckpool_push(dpool, beta->J[v]); beta->J[v] = NULL; }
3808 deckpool_push(dpool, beta->J[cm->M]); beta->J[cm->M] = NULL;
3809 free(beta->L[0]); free(beta->L);
3810 free(beta->R[0]); free(beta->R);
3811 }
3812 else
3813 {
3814 ret_beta->J = beta->J;
3815 ret_beta->L = beta->L;
3816 ret_beta->R = beta->R;
3817 }
3818 free(beta);
3819
3820 if (ret_dpool == NULL)
3821 {
3822 float **a;
3823 while (deckpool_pop(dpool, &a))
3824 {
3825 if (a == NULL) { fprintf(stderr,"WARNING: We've got issues: popped from deckpool but it's NULL!\n"); continue; }
3826 free_vji_deck(a,j1,j0);
3827 }
3828 deckpool_free(dpool);
3829 }
3830 else
3831 {
3832 *ret_dpool = dpool;
3833 }
3834
3835 free(touch);
3836
3837 return;
3838 }
3839
3840 /* Function: tr_insideT()
3841 * Author: DLK
3842 *
3843 * Purpose: inside with traceback on truncated sequence
3844 * based on insideT()
3845 *
3846 * Args: cm - the covariance model
3847 * dsq - the sequence, 1..L
3848 * L - length of the sequence
3849 * tr - Parsetree for traceback
3850 * r - root of subgraph to align to target subseq (usually 0, the model's root)
3851 * z - last state of the subgraph
3852 * i0 - start of target subsequence (usually 1, beginning of dsq)
3853 * j0 - end of target subsequence (usually L, end of dsq)
3854 *
3855 * Returns: score of optimal alignment (float)
3856 */
3857 float
tr_insideT(CM_t * cm,ESL_DSQ * dsq,int L,Parsetree_t * tr,int r,int z,int i0,int j0,int r_allow_J,int r_allow_L,int r_allow_R,int lenCORREX)3858 tr_insideT(CM_t *cm, ESL_DSQ *dsq, int L, Parsetree_t *tr, int r, int z,
3859 int i0, int j0, int r_allow_J, int r_allow_L, int r_allow_R, int lenCORREX)
3860 {
3861 int status; /* easel status code */
3862 void ***shadow; /* standard shadow matrix with state information */
3863 void ***L_shadow; /* left marginal shadow matrix with state information */
3864 void ***R_shadow; /* right marginal shadow matrix with state information */
3865 void ***T_shadow; /* terminal shadow matrix with state information */
3866 void ***Lmode_shadow; /* left marginal shadow matrix with mode information */
3867 void ***Rmode_shadow; /* right marginal shadow matrix with mode information */
3868 float sc; /* score of the CYK alignment */
3869 ESL_STACK *pda; /* stack for storing info of 2nd child at B_st */
3870 int v,i,j,d; /* indices for state, position, & distance */
3871 int mode,nxtmode;
3872 int k;
3873 int y, yoffset;
3874 int bifparent;
3875
3876 ShadowMats_t *all_shadow;
3877 all_shadow = malloc(sizeof(ShadowMats_t));
3878
3879 /*
3880 sc = trinside(cm, dsq, L, r, z, i0, j0,
3881 BE_EFFICIENT,
3882 &shadow,
3883 &L_shadow, &R_shadow, &T_shadow,
3884 &Lmode_shadow, &Rmode_shadow,
3885 &mode, &v, &i, &j );
3886 */
3887
3888 sc = tr_inside(cm, dsq, L, r, z, i0, j0, BE_EFFICIENT,
3889 (r == 0), r_allow_J, r_allow_L, r_allow_R, lenCORREX,
3890 NULL, NULL, NULL, NULL, all_shadow,
3891 &mode, &v, &i, &j);
3892 shadow = all_shadow->J;
3893 L_shadow = all_shadow->L;
3894 R_shadow = all_shadow->R;
3895 T_shadow = all_shadow->T;
3896 Lmode_shadow = all_shadow->Lmode;
3897 Rmode_shadow = all_shadow->Rmode;
3898
3899 if((pda = esl_stack_ICreate()) == NULL) goto ERROR;
3900 d = j-i+1;
3901
3902 if (r == 0)
3903 {
3904 InsertTraceNodewithMode(tr, tr->n-1, TRACE_LEFT_CHILD, i, j, v, mode);
3905 }
3906
3907 while (1)
3908 {
3909 if ( cm->sttype[v] == B_st )
3910 {
3911 if ( mode == TRMODE_J )
3912 {
3913 k = ((int **) shadow[v])[j][d];
3914
3915 if((status = esl_stack_IPush(pda, j)) != eslOK) goto ERROR;
3916 if((status = esl_stack_IPush(pda, k)) != eslOK) goto ERROR;
3917 if((status = esl_stack_IPush(pda, mode)) != eslOK) goto ERROR;
3918 if((status = esl_stack_IPush(pda, tr->n-1)) != eslOK) goto ERROR;
3919 }
3920 else if ( mode == TRMODE_L )
3921 {
3922 k = ((int **) L_shadow[v])[j][d];
3923
3924 /* In left marginal mode, right child is always left marginal */
3925 if((status = esl_stack_IPush(pda, j)) != eslOK) goto ERROR;
3926 if((status = esl_stack_IPush(pda, k)) != eslOK) goto ERROR;
3927 if((status = esl_stack_IPush(pda, mode)) != eslOK) goto ERROR;
3928 if((status = esl_stack_IPush(pda, tr->n-1)) != eslOK) goto ERROR;
3929
3930 /* Retrieve mode of left child (should be 3 or 2) */
3931 mode = ((int **)Lmode_shadow[v])[j][d];
3932 }
3933 else if ( mode == TRMODE_R )
3934 {
3935 k = ((int **) R_shadow[v])[j][d];
3936
3937 /* Retrieve mode of right child (should be 3 or 1) */
3938 mode = ((int **)Rmode_shadow[v])[j][d];
3939 if((status = esl_stack_IPush(pda, j)) != eslOK) goto ERROR;
3940 if((status = esl_stack_IPush(pda, k)) != eslOK) goto ERROR;
3941 if((status = esl_stack_IPush(pda, mode)) != eslOK) goto ERROR;
3942 if((status = esl_stack_IPush(pda, tr->n-1)) != eslOK) goto ERROR;
3943
3944 /* In right marginal mode, left child is always right marginal */
3945 mode = 1;
3946 }
3947 else if ( mode == TRMODE_T )
3948 {
3949 k = ((int **) T_shadow[v])[j][d];
3950
3951 mode = 2;
3952 if((status = esl_stack_IPush(pda, j)) != eslOK) goto ERROR;
3953 if((status = esl_stack_IPush(pda, k)) != eslOK) goto ERROR;
3954 if((status = esl_stack_IPush(pda, mode)) != eslOK) goto ERROR;
3955 if((status = esl_stack_IPush(pda, tr->n-1)) != eslOK) goto ERROR;
3956
3957 mode = 1;
3958 }
3959 else { cm_Die("Unknown mode in traceback!"); }
3960
3961 j = j-k;
3962 d = d-k;
3963 i = j-d+1;
3964 v = cm->cfirst[v];
3965 InsertTraceNodewithMode(tr, tr->n-1, TRACE_LEFT_CHILD, i, j, v, mode);
3966 }
3967 else if ( (cm->sttype[v] == E_st) || (cm->sttype[v] == EL_st) )
3968 {
3969 if (esl_stack_IPop(pda, &bifparent) == eslEOD) break;
3970 esl_stack_IPop(pda, &mode);
3971 esl_stack_IPop(pda, &d);
3972 esl_stack_IPop(pda, &j);
3973 v = tr->state[bifparent];
3974 y = cm->cnum[v];
3975 i = j-d+1;
3976
3977 v = y;
3978 InsertTraceNodewithMode(tr, bifparent, TRACE_RIGHT_CHILD, i, j, v, mode);
3979 }
3980 else
3981 {
3982 if ( mode == TRMODE_J )
3983 {
3984 yoffset = ((char **) shadow[v])[j][d];
3985 nxtmode = 3;
3986 }
3987 else if ( mode == TRMODE_L )
3988 {
3989 yoffset = ((char **) L_shadow[v])[j][d];
3990 nxtmode = ((int **)Lmode_shadow[v])[j][d];
3991 }
3992 else if ( mode == TRMODE_R )
3993 {
3994 yoffset = ((char **) R_shadow[v])[j][d];
3995 nxtmode = ((int **)Rmode_shadow[v])[j][d];
3996 }
3997 else { cm_Die("Unknown mode in traceback!"); }
3998
3999 switch (cm->sttype[v])
4000 {
4001 case D_st:
4002 break;
4003 case MP_st:
4004 if ( mode == TRMODE_J ) i++;
4005 if ( mode == TRMODE_L && d > 0 ) i++;
4006 if ( mode == TRMODE_J ) j--;
4007 if ( mode == TRMODE_R && d > 0 ) j--;
4008 break;
4009 case ML_st:
4010 if ( mode == TRMODE_J ) i++;
4011 if ( mode == TRMODE_L && d > 0 ) i++;
4012 break;
4013 case MR_st:
4014 if ( mode == TRMODE_J ) j--;
4015 if ( mode == TRMODE_R && d > 0 ) j--;
4016 break;
4017 case IL_st:
4018 if ( mode == TRMODE_J ) i++;
4019 if ( mode == TRMODE_L && d > 0 ) i++;
4020 break;
4021 case IR_st:
4022 if ( mode == TRMODE_J ) j--;
4023 if ( mode == TRMODE_R && d > 0 ) j--;
4024 break;
4025 case S_st:
4026 break;
4027 default:
4028 cm_Die("'Inconceivable!'\n'You keep using that word...'");
4029 }
4030 d = j-i+1;
4031
4032 if ( yoffset == USED_EL )
4033 {
4034 v = cm->M;
4035 if (mode == TRMODE_J)
4036 {
4037 InsertTraceNodewithMode(tr, tr->n-1, TRACE_LEFT_CHILD, i, j, v, mode);
4038 }
4039 }
4040 else if ( yoffset == USED_LOCAL_BEGIN )
4041 { /* local begin, can only happen once, from root */
4042 /* However, all hits from truncyk() are local hits, and this should have
4043 been dealt with immediately after return from the DP function.
4044 If we've reached this point, there's a major problem */
4045 cm_Die("Impossible local begin in traceback\n");
4046 }
4047 else
4048 {
4049 mode = nxtmode;
4050 v = cm->cfirst[v] + yoffset;
4051 InsertTraceNodewithMode(tr, tr->n-1, TRACE_LEFT_CHILD, i, j, v, mode);
4052 }
4053 }
4054 }
4055
4056 esl_stack_Destroy(pda);
4057 free_vjd_shadow_matrix(shadow, cm, i0, j0);
4058 free_vjd_shadow_matrix(L_shadow, cm, i0, j0);
4059 free_vjd_shadow_matrix(R_shadow, cm, i0, j0);
4060 free_vjd_shadow_matrix(T_shadow, cm, i0, j0);
4061 free_vjd_shadow_matrix(Lmode_shadow, cm, i0, j0);
4062 free_vjd_shadow_matrix(Rmode_shadow, cm, i0, j0);
4063
4064 return sc;
4065
4066 ERROR:
4067 cm_Fail("Memory error.");
4068 return 0.; /* NEVERREACHED */
4069 }
4070
4071 /* Function: tr_vinsideT()
4072 * Author: DLK
4073 *
4074 * Purpose: Traceback wrapper for tr_vinside()
4075 * Appends trace to a traceback which
4076 * already has state r at t->n-1
4077 * Args:
4078 *
4079 * Returns:
4080 */
4081 float
tr_vinsideT(CM_t * cm,ESL_DSQ * dsq,int L,Parsetree_t * tr,int r,int z,int i0,int i1,int j1,int j0,int useEL,int r_allow_J,int r_allow_L,int r_allow_R,int z_allow_J,int z_allow_L,int z_allow_R)4082 tr_vinsideT(CM_t *cm, ESL_DSQ *dsq, int L, Parsetree_t *tr, int r, int z,
4083 int i0, int i1, int j1, int j0, int useEL,
4084 int r_allow_J, int r_allow_L, int r_allow_R,
4085 int z_allow_J, int z_allow_L, int z_allow_R)
4086 {
4087 float sc;
4088 int v, i, j;
4089 int ip, jp;
4090 int mode, nxtmode;
4091 int yoffset;
4092
4093 AlphaMats_t *alpha;
4094 ShadowMats_t *shadow;
4095 alpha = malloc(sizeof(AlphaMats_t));
4096 shadow = malloc(sizeof(ShadowMats_t));
4097
4098 if (r == z)
4099 {
4100 if ( r_allow_J ) mode = 3;
4101 else if ( r_allow_L ) mode = 2;
4102 else mode = 1;
4103
4104 InsertTraceNodewithMode(tr, tr->n-1, TRACE_LEFT_CHILD, i0, j0, r, mode);
4105 return 0.0;
4106 }
4107
4108 sc = tr_vinside(cm, dsq, L, r, z, i0, i1, j1, j0, useEL, BE_EFFICIENT, (r == 0),
4109 r_allow_J, r_allow_L, r_allow_R, z_allow_J, z_allow_L, z_allow_R,
4110 NULL, alpha, NULL, NULL, shadow, &mode, &v, &i, &j);
4111
4112 if (r == 0)
4113 {
4114 InsertTraceNodewithMode(tr, tr->n-1, TRACE_LEFT_CHILD, i, j, v, mode);
4115 }
4116
4117 if (r != 0 && r != v)
4118 {
4119 v = r;
4120 i = i0;
4121 j = j0;
4122 ip = 0; jp = j0-j1;
4123 mode = 3;
4124 if (alpha->L[v][jp][ip] > alpha->J[v][jp][ip])
4125 mode = 2;
4126 if (alpha->R[v][jp][ip] > alpha->J[v][jp][ip] && alpha->R[v][jp][ip] > alpha->L[v][jp][ip])
4127 mode = 1;
4128 }
4129
4130 free_vji_matrix(alpha->J, cm->M, j1, j0);
4131 free_vji_matrix(alpha->L, cm->M, j1, j0);
4132 free_vji_matrix(alpha->R, cm->M, j1, j0);
4133 free(alpha);
4134
4135 /* start traceback */
4136 while (v != z)
4137 {
4138 jp = j-j1;
4139 ip = i-i0;
4140
4141 if ( mode == TRMODE_J )
4142 {
4143 yoffset = ((char **) shadow->J[v])[jp][ip];
4144 nxtmode = 3;
4145 }
4146 else if ( mode == TRMODE_L )
4147 {
4148 yoffset = ((char **) shadow->L[v])[jp][ip];
4149 nxtmode = ((char **) shadow->Lmode[v])[jp][ip];
4150 }
4151 else if ( mode == TRMODE_R )
4152 {
4153 yoffset = ((char **) shadow->R[v])[jp][ip];
4154 nxtmode = ((char **) shadow->Rmode[v])[jp][ip];
4155 }
4156 else
4157 cm_Die("Unknown mode in traceback!\n");
4158
4159 switch (cm->sttype[v])
4160 {
4161 case S_st:
4162 case D_st:
4163 break;
4164 case MP_st:
4165 if ( mode == TRMODE_J || mode == TRMODE_L ) i++;
4166 if ( mode == TRMODE_J || mode == TRMODE_R ) j--;
4167 break;
4168 case ML_st:
4169 case IL_st:
4170 if ( mode == TRMODE_J || mode == TRMODE_L ) i++;
4171 break;
4172 case MR_st:
4173 case IR_st:
4174 if ( mode == TRMODE_J || mode == TRMODE_R ) j--;
4175 break;
4176 default:
4177 cm_Die("'Inconceivable!'\n'Youu keep using that word...'");
4178 }
4179 mode = nxtmode;
4180
4181 if (yoffset == USED_EL)
4182 {
4183 v = cm->M;
4184 InsertTraceNodewithMode(tr, tr->n-1, TRACE_LEFT_CHILD, i, j, v, mode);
4185 break;
4186 }
4187 else if (yoffset == USED_LOCAL_BEGIN)
4188 { /* local begin, can only happen once, from root */
4189 if (v != 0)
4190 cm_Die("Impossible local begin in traceback!\n");
4191 else
4192 cm_Die("DEV: you actually need to deal with this local begin case\n");
4193 }
4194 else
4195 {
4196 v = cm->cfirst[v] + yoffset;
4197 InsertTraceNodewithMode(tr, tr->n-1, TRACE_LEFT_CHILD, i, j, v, mode);
4198 }
4199 }
4200
4201 if (useEL)
4202 {
4203 switch (cm->sttype[z])
4204 {
4205 case MP_st: i++; j--; break;
4206 case ML_st:
4207 case IL_st: i++; break;
4208 case MR_st:
4209 case IR_st: j--; break;
4210 }
4211
4212 InsertTraceNodewithMode(tr, tr->n-1, TRACE_LEFT_CHILD, i, j, cm->M, 3);
4213 }
4214
4215 free_vji_shadow_matrix((char ***) shadow->J, cm->M, j1, j0);
4216 free_vji_shadow_matrix((char ***) shadow->L, cm->M, j1, j0);
4217 free_vji_shadow_matrix((char ***) shadow->R, cm->M, j1, j0);
4218 free_vji_shadow_matrix((char ***) shadow->Lmode, cm->M, j1, j0);
4219 free_vji_shadow_matrix((char ***) shadow->Rmode, cm->M, j1, j0);
4220 free(shadow);
4221
4222 return sc;
4223 }
4224