1 #ifdef HAVE_CONFIG_H
2 #include "config.h"
3 #endif
4 
5 #include <stdio.h>
6 #include <stdlib.h>
7 #include <math.h>
8 #include <ctype.h>
9 #include <string.h>
10 #include <limits.h>
11 
12 #include "ViennaRNA/fold_vars.h"
13 #include "ViennaRNA/params/default.h"
14 #include "ViennaRNA/utils/basic.h"
15 #include "ViennaRNA/alphabet.h"
16 #include "ViennaRNA/constraints/hard.h"
17 #include "ViennaRNA/constraints/soft.h"
18 #include "ViennaRNA/gquad.h"
19 #include "ViennaRNA/structured_domains.h"
20 #include "ViennaRNA/unstructured_domains.h"
21 #include "ViennaRNA/loops/external.h"
22 
23 #ifdef __GNUC__
24 # define INLINE inline
25 #else
26 # define INLINE
27 #endif
28 
29 #ifdef VRNA_WITH_SVM
30 #include "ViennaRNA/zscore_dat.inc"
31 #endif
32 
33 #include "external_hc.inc"
34 #include "external_sc.inc"
35 
36 /*
37  #################################
38  # PRIVATE FUNCTION DECLARATIONS #
39  #################################
40  */
41 
42 PRIVATE int
43 BT_ext_loop_f5(vrna_fold_compound_t *fc,
44                int                  *k,
45                int                  *i,
46                int                  *j,
47                vrna_bp_stack_t      *bp_stack,
48                int                  *stack_count);
49 
50 
51 PRIVATE int
52 BT_ext_loop_f5_comparative(vrna_fold_compound_t *fc,
53                            int                  *k,
54                            int                  *i,
55                            int                  *j,
56                            vrna_bp_stack_t      *bp_stack,
57                            int                  *stack_count);
58 
59 
60 PRIVATE int
61 BT_ext_loop_f3(vrna_fold_compound_t *fc,
62                int                  *k,
63                int                  maxdist,
64                int                  *i,
65                int                  *j,
66                vrna_bp_stack_t      *bp_stack,
67                int                  *stack_count);
68 
69 
70 PRIVATE int
71 BT_ext_loop_f3_comparative(vrna_fold_compound_t *fc,
72                            int                  *k,
73                            int                  maxdist,
74                            int                  *i,
75                            int                  *j,
76                            vrna_bp_stack_t      *bp_stack,
77                            int                  *stack_count);
78 
79 
80 PRIVATE int
81 BT_ext_loop_f3_pp(vrna_fold_compound_t  *fc,
82                   int                   *i,
83                   int                   maxj);
84 
85 
86 PRIVATE int
87 BT_ext_loop_f3_pp_comparative(vrna_fold_compound_t  *fc,
88                               int                   *i,
89                               int                   maxj);
90 
91 
92 /*
93  #################################
94  # BEGIN OF FUNCTION DEFINITIONS #
95  #################################
96  */
97 PUBLIC int
vrna_BT_ext_loop_f5(vrna_fold_compound_t * fc,int * k,int * i,int * j,vrna_bp_stack_t * bp_stack,int * stack_count)98 vrna_BT_ext_loop_f5(vrna_fold_compound_t  *fc,
99                     int                   *k,
100                     int                   *i,
101                     int                   *j,
102                     vrna_bp_stack_t       *bp_stack,
103                     int                   *stack_count)
104 {
105   if (fc) {
106     switch (fc->type) {
107       case VRNA_FC_TYPE_SINGLE:
108         return BT_ext_loop_f5(fc, k, i, j, bp_stack, stack_count);
109         break;
110 
111       case VRNA_FC_TYPE_COMPARATIVE:
112         return BT_ext_loop_f5_comparative(fc, k, i, j, bp_stack, stack_count);
113         break;
114     }
115   }
116 
117   return -1;
118 }
119 
120 
121 PUBLIC int
vrna_BT_ext_loop_f3(vrna_fold_compound_t * fc,int * k,int maxdist,int * i,int * j,vrna_bp_stack_t * bp_stack,int * stack_count)122 vrna_BT_ext_loop_f3(vrna_fold_compound_t  *fc,
123                     int                   *k,
124                     int                   maxdist,
125                     int                   *i,
126                     int                   *j,
127                     vrna_bp_stack_t       *bp_stack,
128                     int                   *stack_count)
129 {
130   if (fc) {
131     switch (fc->type) {
132       case VRNA_FC_TYPE_SINGLE:
133         return BT_ext_loop_f3(fc, k, maxdist, i, j, bp_stack, stack_count);
134         break;
135 
136       case VRNA_FC_TYPE_COMPARATIVE:
137         return BT_ext_loop_f3_comparative(fc, k, maxdist, i, j, bp_stack, stack_count);
138         break;
139     }
140   }
141 
142   return -1;
143 }
144 
145 
146 PUBLIC int
vrna_BT_ext_loop_f3_pp(vrna_fold_compound_t * fc,int * i,int maxj)147 vrna_BT_ext_loop_f3_pp(vrna_fold_compound_t *fc,
148                        int                  *i,
149                        int                  maxj)
150 {
151   if (fc) {
152     switch (fc->type) {
153       case VRNA_FC_TYPE_SINGLE:
154         return BT_ext_loop_f3_pp(fc, i, maxj);
155         break;
156 
157       case VRNA_FC_TYPE_COMPARATIVE:
158         return BT_ext_loop_f3_pp_comparative(fc, i, maxj);
159         break;
160     }
161   }
162 
163   return -1;
164 }
165 
166 
167 PRIVATE int
BT_ext_loop_f5(vrna_fold_compound_t * fc,int * k,int * i,int * j,vrna_bp_stack_t * bp_stack,int * stack_count)168 BT_ext_loop_f5(vrna_fold_compound_t *fc,
169                int                  *k,
170                int                  *i,
171                int                  *j,
172                vrna_bp_stack_t      *bp_stack,
173                int                  *stack_count)
174 {
175   char                      *ptype;
176   short                     mm5, mm3, *S1;
177   unsigned int              *sn, type;
178   int                       length, fij, fi, jj, u, en, e, *my_f5, *my_c, *my_ggg, *idx,
179                             dangle_model, turn, with_gquad, cnt, ii, with_ud;
180   vrna_param_t              *P;
181   vrna_md_t                 *md;
182   vrna_sc_t                 *sc;
183   vrna_ud_t                 *domains_up;
184   vrna_callback_hc_evaluate *evaluate;
185   struct hc_ext_def_dat     hc_dat_local;
186 
187   length        = fc->length;
188   P             = fc->params;
189   md            = &(P->model_details);
190   sn            = fc->strand_number;
191   sc            = fc->sc;
192   my_f5         = fc->matrices->f5;
193   my_c          = fc->matrices->c;
194   my_ggg        = fc->matrices->ggg;
195   domains_up    = fc->domains_up;
196   idx           = fc->jindx;
197   ptype         = fc->ptype;
198   S1            = fc->sequence_encoding;
199   dangle_model  = md->dangles;
200   turn          = md->min_loop_size;
201   with_gquad    = md->gquad;
202   with_ud       = (domains_up && domains_up->energy_cb) ? 1 : 0;
203   evaluate      = prepare_hc_ext_def(fc, &hc_dat_local);
204 
205   jj = *k;
206 
207   /* nibble off unpaired 3' stretches harboring bound ligands (interspersed with unpaired nucleotides) */
208   if (with_ud) {
209     do {
210       fij = my_f5[jj];
211       fi  = INF;
212 
213       /* try nibble off one unpaired nucleotide first */
214       if (evaluate(1, jj, 1, jj - 1, VRNA_DECOMP_EXT_EXT, &hc_dat_local)) {
215         fi = my_f5[jj - 1];
216 
217         if (sc) {
218           if (sc->energy_up)
219             fi += sc->energy_up[jj][1];
220 
221           if (sc->f)
222             fi += sc->f(1, jj, 1, jj - 1, VRNA_DECOMP_EXT_EXT, sc->data);
223         }
224 
225         if (jj == 1) {
226           /* no more pairs */
227           *i  = *j = -1;
228           *k  = 0;
229           return 1;
230         }
231 
232         if (fij == fi) {
233           jj--;
234           continue;
235         }
236       }
237 
238       /* next, try nibble off a ligand */
239       for (cnt = 0; cnt < domains_up->uniq_motif_count; cnt++) {
240         u   = domains_up->uniq_motif_size[cnt];
241         ii  = jj - u + 1;
242         if ((ii > 0) && evaluate(1, jj, 1, jj - u, VRNA_DECOMP_EXT_EXT, &hc_dat_local)) {
243           en = domains_up->energy_cb(fc,
244                                      ii,
245                                      jj,
246                                      VRNA_UNSTRUCTURED_DOMAIN_EXT_LOOP | VRNA_UNSTRUCTURED_DOMAIN_MOTIF,
247                                      domains_up->data);
248           if (sc) {
249             if (sc->energy_up)
250               en += sc->energy_up[ii][u];
251 
252             if (sc->f)
253               en += sc->f(1, jj, 1, jj - u, VRNA_DECOMP_EXT_EXT, sc->data);
254           }
255 
256           fi  = my_f5[ii - 1];
257           fi  += en;
258 
259           if (fij == fi) {
260             /* skip remaining motifs after first hit */
261             jj = ii - 1;
262             break;
263           }
264         }
265       }
266 
267       if (jj == 0) {
268         /* no more pairs */
269         *i  = *j = -1;
270         *k  = 0;
271         return 1;
272       }
273     } while (fij == fi);
274   } else {
275     /* nibble off unpaired 3' bases */
276     do {
277       fij = my_f5[jj];
278       fi  = INF;
279 
280       if (evaluate(1, jj, 1, jj - 1, VRNA_DECOMP_EXT_EXT, &hc_dat_local)) {
281         fi = my_f5[jj - 1];
282 
283         if (sc) {
284           if (sc->energy_up)
285             fi += sc->energy_up[jj][1];
286 
287           if (sc->f)
288             fi += sc->f(1, jj, 1, jj - 1, VRNA_DECOMP_EXT_EXT, sc->data);
289         }
290       }
291 
292       if (--jj == 0)
293         break;
294     } while (fij == fi);
295     jj++;
296   }
297 
298   if (jj < turn + 2) {
299     /* no more pairs */
300     *i  = *j = -1;
301     *k  = 0;
302     return 1;
303   }
304 
305   /* must have found a decomposition */
306   switch (dangle_model) {
307     case 0:   /* j is paired. Find pairing partner */
308       for (u = jj - turn - 1; u >= 1; u--) {
309         if (with_gquad) {
310           if (fij == my_f5[u - 1] + my_ggg[idx[jj] + u]) {
311             *i  = *j = -1;
312             *k  = u - 1;
313             return vrna_BT_gquad_mfe(fc, u, jj, bp_stack, stack_count);
314           }
315         }
316 
317         if (evaluate(1, jj, u - 1, u, VRNA_DECOMP_EXT_EXT_STEM, &hc_dat_local)) {
318           type = vrna_get_ptype(idx[jj] + u, ptype);
319 
320           en = my_c[idx[jj] + u];
321           if (sc)
322             if (sc->f)
323               en += sc->f(1, jj, u - 1, u, VRNA_DECOMP_EXT_EXT_STEM, sc->data);
324 
325           if (sn[jj] != sn[u])
326             en += P->DuplexInit;
327 
328           if (fij == vrna_E_ext_stem(type, -1, -1, P) + en + my_f5[u - 1]) {
329             *i                            = u;
330             *j                            = jj;
331             *k                            = u - 1;
332             bp_stack[++(*stack_count)].i  = u;
333             bp_stack[(*stack_count)].j    = jj;
334             return 1;
335           }
336         }
337       }
338       break;
339 
340     case 2:
341       mm3 = ((jj < length) && (sn[jj + 1] == sn[jj])) ? S1[jj + 1] : -1;
342       for (u = jj - turn - 1; u >= 1; u--) {
343         if (with_gquad) {
344           if (fij == my_f5[u - 1] + my_ggg[idx[jj] + u]) {
345             *i  = *j = -1;
346             *k  = u - 1;
347             return vrna_BT_gquad_mfe(fc, u, jj, bp_stack, stack_count);
348           }
349         }
350 
351         if (evaluate(1, jj, u - 1, u, VRNA_DECOMP_EXT_EXT_STEM, &hc_dat_local)) {
352           mm5   = ((u > 1) && (sn[u] == sn[u - 1])) ? S1[u - 1] : -1;
353           type  = vrna_get_ptype(idx[jj] + u, ptype);
354 
355           en = my_c[idx[jj] + u];
356           if (sc)
357             if (sc->f)
358               en += sc->f(1, jj, u - 1, u, VRNA_DECOMP_EXT_EXT_STEM, sc->data);
359 
360           if (sn[jj] != sn[u])
361             en += P->DuplexInit;
362 
363           if (fij == vrna_E_ext_stem(type, mm5, mm3, P) + en + my_f5[u - 1]) {
364             *i                            = u;
365             *j                            = jj;
366             *k                            = u - 1;
367             bp_stack[++(*stack_count)].i  = u;
368             bp_stack[(*stack_count)].j    = jj;
369             return 1;
370           }
371         }
372       }
373       break;
374 
375     default:
376       if (with_gquad) {
377         if (fij == my_ggg[idx[jj] + 1]) {
378           *i  = *j = -1;
379           *k  = 0;
380           return vrna_BT_gquad_mfe(fc, 1, jj, bp_stack, stack_count);
381         }
382       }
383 
384       if (evaluate(1, jj, 1, jj, VRNA_DECOMP_EXT_STEM, &hc_dat_local)) {
385         type = vrna_get_ptype(idx[jj] + 1, ptype);
386 
387         en = my_c[idx[jj] + 1];
388         if (sc)
389           if (sc->f)
390             en += sc->f(1, jj, 1, jj, VRNA_DECOMP_EXT_STEM, sc->data);
391 
392         if (sn[jj] != sn[1])
393           en += P->DuplexInit;
394 
395         if (fij == en + vrna_E_ext_stem(type, -1, -1, P)) {
396           *i                            = 1;
397           *j                            = jj;
398           *k                            = 0;
399           bp_stack[++(*stack_count)].i  = 1;
400           bp_stack[(*stack_count)].j    = jj;
401           return 1;
402         }
403       }
404 
405       if (evaluate(1, jj, 1, jj - 1, VRNA_DECOMP_EXT_STEM, &hc_dat_local)) {
406         if (sn[jj] == sn[jj - 1]) {
407           mm3   = S1[jj];
408           type  = vrna_get_ptype(idx[jj - 1] + 1, ptype);
409 
410           en = my_c[idx[jj - 1] + 1];
411           if (sc) {
412             if (sc->energy_up)
413               en += sc->energy_up[jj][1];
414 
415             if (sc->f)
416               en += sc->f(1, jj, 1, jj - 1, VRNA_DECOMP_EXT_STEM, sc->data);
417           }
418 
419           if (sn[jj - 1] != sn[1])
420             en += P->DuplexInit;
421 
422           if (fij == en + vrna_E_ext_stem(type, -1, mm3, P)) {
423             *i                            = 1;
424             *j                            = jj - 1;
425             *k                            = 0;
426             bp_stack[++(*stack_count)].i  = 1;
427             bp_stack[(*stack_count)].j    = jj - 1;
428             return 1;
429           }
430         }
431       }
432 
433       for (u = jj - turn - 1; u > 1; u--) {
434         if (with_gquad) {
435           if (fij == my_f5[u - 1] + my_ggg[idx[jj] + u]) {
436             *i  = *j = -1;
437             *k  = u - 1;
438             return vrna_BT_gquad_mfe(fc, u, jj, bp_stack, stack_count);
439           }
440         }
441 
442         type = vrna_get_ptype(idx[jj] + u, ptype);
443 
444         en = my_c[idx[jj] + u];
445         if (sn[jj] != sn[u])
446           en += P->DuplexInit;
447 
448         if (evaluate(1, jj, u - 1, u, VRNA_DECOMP_EXT_EXT_STEM, &hc_dat_local)) {
449           e = my_f5[u - 1] +
450               en +
451               vrna_E_ext_stem(type, -1, -1, P);
452 
453           if (sc)
454             if (sc->f)
455               e += sc->f(1, jj, u - 1, u, VRNA_DECOMP_EXT_EXT_STEM, sc->data);
456 
457           if (fij == e) {
458             *i                            = u;
459             *j                            = jj;
460             *k                            = u - 1;
461             bp_stack[++(*stack_count)].i  = u;
462             bp_stack[(*stack_count)].j    = jj;
463             return 1;
464           }
465         }
466 
467         if (evaluate(1, jj, u - 2, u, VRNA_DECOMP_EXT_EXT_STEM, &hc_dat_local)) {
468           if (sn[u] == sn[u - 1]) {
469             mm5 = S1[u - 1];
470             e   = my_f5[u - 2] +
471                   en +
472                   vrna_E_ext_stem(type, mm5, -1, P);
473 
474             if (sc) {
475               if (sc->energy_up)
476                 e += sc->energy_up[u - 1][1];
477 
478               if (sc->f)
479                 e += sc->f(1, jj, u - 2, u, VRNA_DECOMP_EXT_EXT_STEM, sc->data);
480             }
481 
482             if (fij == e) {
483               *i                            = u;
484               *j                            = jj;
485               *k                            = u - 2;
486               bp_stack[++(*stack_count)].i  = u;
487               bp_stack[(*stack_count)].j    = jj;
488               return 1;
489             }
490           }
491         }
492 
493         type = vrna_get_ptype(idx[jj - 1] + u, ptype);
494 
495         en = my_c[idx[jj - 1] + u];
496         if (sn[jj - 1] != sn[u])
497           en += P->DuplexInit;
498 
499         mm5 = (sn[u] == sn[u - 1]) ? S1[u - 1] : -1;
500         mm3 = (sn[jj] == sn[jj - 1]) ? S1[jj] : -1;
501 
502         if (evaluate(1, jj, u - 1, u, VRNA_DECOMP_EXT_EXT_STEM1, &hc_dat_local)) {
503           e = my_f5[u - 1] +
504               en +
505               vrna_E_ext_stem(type, -1, mm3, P);
506 
507           if (sc) {
508             if (sc->energy_up)
509               e += sc->energy_up[jj][1];
510 
511             if (sc->f)
512               e += sc->f(1, jj, u - 1, u, VRNA_DECOMP_EXT_EXT_STEM1, sc->data);
513           }
514 
515           if (fij == e) {
516             *i                            = u;
517             *j                            = jj - 1;
518             *k                            = u - 1;
519             bp_stack[++(*stack_count)].i  = u;
520             bp_stack[(*stack_count)].j    = jj - 1;
521             return 1;
522           }
523         }
524 
525         if (evaluate(1, jj, u - 2, u, VRNA_DECOMP_EXT_EXT_STEM1, &hc_dat_local)) {
526           e = my_f5[u - 2] + en + vrna_E_ext_stem(type, mm5, mm3, P);
527           if (sc) {
528             if (sc->energy_up)
529               e += sc->energy_up[jj][1] +
530                    sc->energy_up[u - 1][1];
531 
532             if (sc->f)
533               e += sc->f(1, jj, u - 2, u, VRNA_DECOMP_EXT_EXT_STEM1, sc->data);
534           }
535 
536           if (fij == e) {
537             *i                            = u;
538             *j                            = jj - 1;
539             *k                            = u - 2;
540             bp_stack[++(*stack_count)].i  = u;
541             bp_stack[(*stack_count)].j    = jj - 1;
542             return 1;
543           }
544         }
545       }
546       break;
547   }
548 
549   return 0;
550 }
551 
552 
553 PRIVATE int
BT_ext_loop_f5_comparative(vrna_fold_compound_t * fc,int * k,int * i,int * j,vrna_bp_stack_t * bp_stack,int * stack_count)554 BT_ext_loop_f5_comparative(vrna_fold_compound_t *fc,
555                            int                  *k,
556                            int                  *i,
557                            int                  *j,
558                            vrna_bp_stack_t      *bp_stack,
559                            int                  *stack_count)
560 {
561   unsigned int              **a2s, n;
562   short                     **S, **S5, **S3;
563   unsigned int              tt;
564   int                       fij, fi, jj, u, en, *my_f5, *my_c, *my_ggg, *idx,
565                             dangle_model, turn, with_gquad, n_seq, ss, mm5, mm3;
566   vrna_param_t              *P;
567   vrna_md_t                 *md;
568   vrna_sc_t                 **scs;
569   vrna_callback_hc_evaluate *evaluate;
570   struct hc_ext_def_dat     hc_dat_local;
571 
572   n_seq         = fc->n_seq;
573   n             = fc->length;
574   S             = fc->S;
575   S5            = fc->S5;
576   S3            = fc->S3;
577   a2s           = fc->a2s;
578   P             = fc->params;
579   md            = &(P->model_details);
580   scs           = fc->scs;
581   my_f5         = fc->matrices->f5;
582   my_c          = fc->matrices->c;
583   my_ggg        = fc->matrices->ggg;
584   idx           = fc->jindx;
585   dangle_model  = md->dangles;
586   turn          = md->min_loop_size;
587   with_gquad    = md->gquad;
588   evaluate      = prepare_hc_ext_def(fc, &hc_dat_local);
589 
590   jj = *k;
591 
592   /* nibble off unpaired 3' bases */
593   do {
594     fij = my_f5[jj];
595     fi  = INF;
596 
597     if (evaluate(1, jj, 1, jj - 1, VRNA_DECOMP_EXT_EXT, &hc_dat_local)) {
598       fi = my_f5[jj - 1];
599 
600       if (scs) {
601         for (ss = 0; ss < n_seq; ss++)
602           if (scs[ss]) {
603             if (scs[ss]->energy_up)
604               fi += scs[ss]->energy_up[a2s[ss][jj]][1];
605 
606             if (scs[ss]->f)
607               fi += scs[ss]->f(1, jj, 1, jj - 1, VRNA_DECOMP_EXT_EXT, scs[ss]->data);
608           }
609       }
610     }
611 
612     if (--jj == 0)
613       break;
614   } while (fij == fi);
615   jj++;
616 
617   if (jj < turn + 2) {
618     /* no more pairs */
619     *i  = *j = -1;
620     *k  = 0;
621     return 1;
622   }
623 
624   /* must have found a decomposition */
625   switch (dangle_model) {
626     case 0:   /* j is paired. Find pairing partner */
627       for (u = jj - turn - 1; u >= 1; u--) {
628         if (with_gquad) {
629           if (fij == my_f5[u - 1] + my_ggg[idx[jj] + u]) {
630             *i  = *j = -1;
631             *k  = u - 1;
632             return vrna_BT_gquad_mfe(fc, u, jj, bp_stack, stack_count);
633           }
634         }
635 
636         if (evaluate(1, jj, u - 1, u, VRNA_DECOMP_EXT_EXT_STEM, &hc_dat_local)) {
637           en = my_c[idx[jj] + u] +
638                my_f5[u - 1];
639 
640           for (ss = 0; ss < n_seq; ss++) {
641             tt  = vrna_get_ptype_md(S[ss][u], S[ss][jj], md);
642             en  += vrna_E_ext_stem(tt, -1, -1, P);
643           }
644 
645           if (scs) {
646             for (ss = 0; ss < n_seq; ss++)
647               if (scs[ss] && scs[ss]->f)
648                 en += scs[ss]->f(1, jj, u - 1, u, VRNA_DECOMP_EXT_EXT_STEM, scs[ss]->data);
649           }
650 
651           if (fij == en) {
652             *i                            = u;
653             *j                            = jj;
654             *k                            = u - 1;
655             bp_stack[++(*stack_count)].i  = u;
656             bp_stack[(*stack_count)].j    = jj;
657             return 1;
658           }
659         }
660       }
661       break;
662 
663     case 2:
664       for (u = jj - turn - 1; u >= 1; u--) {
665         if (with_gquad) {
666           if (fij == my_f5[u - 1] + my_ggg[idx[jj] + u]) {
667             *i  = *j = -1;
668             *k  = u - 1;
669             return vrna_BT_gquad_mfe(fc, u, jj, bp_stack, stack_count);
670           }
671         }
672 
673         if (evaluate(1, jj, u - 1, u, VRNA_DECOMP_EXT_EXT_STEM, &hc_dat_local)) {
674           en = my_c[idx[jj] + u] +
675                my_f5[u - 1];
676 
677           for (ss = 0; ss < n_seq; ss++) {
678             tt  = vrna_get_ptype_md(S[ss][u], S[ss][jj], md);
679             mm5 = (a2s[ss][u] > 1) ? S5[ss][u] : -1;
680             mm3 = (a2s[ss][jj] < a2s[ss][n]) ? S3[ss][jj] : -1;
681             en  += vrna_E_ext_stem(tt, mm5, mm3, P);
682           }
683 
684           if (scs) {
685             for (ss = 0; ss < n_seq; ss++)
686               if (scs[ss] && scs[ss]->f)
687                 en += scs[ss]->f(1, jj, u - 1, u, VRNA_DECOMP_EXT_EXT_STEM, scs[ss]->data);
688           }
689 
690           if (fij == en) {
691             *i                            = u;
692             *j                            = jj;
693             *k                            = u - 1;
694             bp_stack[++(*stack_count)].i  = u;
695             bp_stack[(*stack_count)].j    = jj;
696             return 1;
697           }
698         }
699       }
700       break;
701   }
702 
703   return 0;
704 }
705 
706 
707 PRIVATE int
BT_ext_loop_f3(vrna_fold_compound_t * fc,int * k,int maxdist,int * i,int * j,vrna_bp_stack_t * bp_stack,int * stack_count)708 BT_ext_loop_f3(vrna_fold_compound_t *fc,
709                int                  *k,
710                int                  maxdist,
711                int                  *i,
712                int                  *j,
713                vrna_bp_stack_t      *bp_stack,
714                int                  *stack_count)
715 {
716   char                      **ptype;
717   short                     mm5, mm3, *S1;
718   unsigned int              type;
719   int                       length, fij, fj, ii, u, *f3, **c, **ggg,
720                             dangle_model, turn, with_gquad, en;
721   vrna_param_t              *P;
722   vrna_md_t                 *md;
723   vrna_sc_t                 *sc;
724   vrna_callback_hc_evaluate *evaluate;
725   struct hc_ext_def_dat     hc_dat_local;
726 
727   length        = fc->length;
728   P             = fc->params;
729   md            = &(P->model_details);
730   sc            = fc->sc;
731   f3            = fc->matrices->f3_local;
732   c             = fc->matrices->c_local;
733   ggg           = fc->matrices->ggg_local;
734   ptype         = fc->ptype_local;
735   S1            = fc->sequence_encoding;
736   dangle_model  = md->dangles;
737   turn          = md->min_loop_size;
738   with_gquad    = md->gquad;
739   evaluate      = prepare_hc_ext_def_window(fc, &hc_dat_local);
740 
741   ii = *k;
742 
743   /* nibble off unpaired 5' bases */
744   do {
745     fij = f3[ii];
746     fj  = INF;
747 
748     if (evaluate(ii, length, ii + 1, length, VRNA_DECOMP_EXT_EXT, &hc_dat_local)) {
749       fj = f3[ii + 1];
750       if (sc) {
751         if (sc->energy_up)
752           fj += sc->energy_up[ii][1];
753 
754         if (sc->f)
755           fj += sc->f(ii, length, ii + 1, length, VRNA_DECOMP_EXT_EXT, sc->data);
756       }
757     }
758 
759     if (++ii > maxdist)
760       break;
761   } while (fij == fj);
762   ii--;
763 
764   if (ii > maxdist - turn + 1) {
765     /* no more pairs */
766     *i  = *j = -1;
767     *k  = 0;
768     return 1;
769   }
770 
771   /*
772    *  must have found a decomposition
773    *  i is paired. Find pairing partner
774    */
775   switch (dangle_model) {
776     /* no dangles */
777     case 0:
778       for (u = maxdist; u > ii + turn; u--) {
779         if (with_gquad) {
780           if (fij == ggg[ii][u - ii] + f3[u + 1]) {
781             *i  = *j = -1;
782             *k  = u + 1;
783             return vrna_BT_gquad_mfe(fc, ii, u, bp_stack, stack_count);
784           }
785         }
786 
787         if (evaluate(ii, length, u, u + 1, VRNA_DECOMP_EXT_STEM_EXT, &hc_dat_local)) {
788           type = vrna_get_ptype_window(ii, u, ptype);
789 
790           en = c[ii][u - ii] +
791                vrna_E_ext_stem(type, -1, -1, P) +
792                f3[u + 1];
793 
794           if ((sc) && (sc->f))
795             en += sc->f(ii, length, u, u + 1, VRNA_DECOMP_EXT_STEM_EXT, sc->data);
796 
797           if (fij == en) {
798             *i                            = ii;
799             *j                            = u;
800             *k                            = u + 1;
801             bp_stack[++(*stack_count)].i  = ii;
802             bp_stack[(*stack_count)].j    = u;
803             return 1;
804           }
805         }
806       }
807       break;
808 
809     /* dangles on both sides */
810     case 2:
811       mm5 = (ii > 1) ? S1[ii - 1] : -1;
812       for (u = maxdist; u > ii + turn; u--) {
813         if (with_gquad) {
814           if (fij == ggg[ii][u - ii] + f3[u + 1]) {
815             *i  = *j = -1;
816             *k  = u + 1;
817             return vrna_BT_gquad_mfe(fc, ii, u, bp_stack, stack_count);
818           }
819         }
820 
821         if (evaluate(ii, length, u, u + 1, VRNA_DECOMP_EXT_STEM_EXT, &hc_dat_local)) {
822           mm3   = (u < length) ? S1[u + 1] : -1;
823           type  = vrna_get_ptype_window(ii, u, ptype);
824 
825           en = c[ii][u - ii] +
826                vrna_E_ext_stem(type, mm5, mm3, P) +
827                f3[u + 1];
828 
829           if ((sc) && (sc->f))
830             en += sc->f(ii, length, u, u + 1, VRNA_DECOMP_EXT_STEM_EXT, sc->data);
831 
832           if (fij == en) {
833             *i                            = ii;
834             *j                            = u;
835             *k                            = u + 1;
836             bp_stack[++(*stack_count)].i  = ii;
837             bp_stack[(*stack_count)].j    = u;
838             return 1;
839           }
840         }
841       }
842       break;
843 
844     default:
845       mm5 = S1[ii];
846       for (u = maxdist; u > ii + turn; u--) {
847         if (with_gquad) {
848           if (fij == ggg[ii][u - ii] + f3[u + 1]) {
849             *i  = *j = -1;
850             *k  = u + 1;
851             return vrna_BT_gquad_mfe(fc, ii, u, bp_stack, stack_count);
852           }
853         }
854 
855         if (u + 2 <= length) {
856           if (evaluate(ii, length, u, u + 2, VRNA_DECOMP_EXT_STEM_EXT1, &hc_dat_local)) {
857             mm3   = S1[u + 1];
858             type  = vrna_get_ptype_window(ii + 1, u, ptype);
859 
860             en = c[ii + 1][u - ii - 1] + vrna_E_ext_stem(type, mm5, mm3, P) + f3[u + 2];
861 
862             if (sc) {
863               if (sc->energy_up)
864                 en += sc->energy_up[u + 1][1] +
865                       sc->energy_up[ii][1];
866 
867               if (sc->f)
868                 en += sc->f(ii, length, u, u + 2, VRNA_DECOMP_EXT_STEM_EXT1, sc->data);
869             }
870 
871             if (fij == en) {
872               *i                            = ii + 1;
873               *j                            = u;
874               *k                            = u + 2;
875               bp_stack[++(*stack_count)].i  = ii + 1;
876               bp_stack[(*stack_count)].j    = u;
877               return 1;
878             }
879           }
880         } else {
881           if (evaluate(ii, length, ii + 1, u, VRNA_DECOMP_EXT_STEM, &hc_dat_local)) {
882             mm3   = (u < length) ? S1[u + 1] : -1;
883             type  = vrna_get_ptype_window(ii + 1, u, ptype);
884 
885             en = c[ii + 1][u - ii - 1] +
886                  vrna_E_ext_stem(type, mm5, mm3, P);
887 
888             if (sc) {
889               if (sc->energy_up) {
890                 en += sc->energy_up[ii][1];
891                 if (u < length)
892                   en += sc->energy_up[u + 1][1];
893               }
894 
895               if (sc->f)
896                 en += sc->f(ii, length, ii + 1, u, VRNA_DECOMP_EXT_STEM, sc->data);
897             }
898 
899             if (fij == en) {
900               *i                            = ii + 1;
901               *j                            = u;
902               *k                            = u + 2;
903               bp_stack[++(*stack_count)].i  = ii + 1;
904               bp_stack[(*stack_count)].j    = u;
905               return 1;
906             }
907           }
908         }
909 
910         if (evaluate(ii, length, u, u + 1, VRNA_DECOMP_EXT_STEM_EXT1, &hc_dat_local)) {
911           type = vrna_get_ptype_window(ii + 1, u, ptype);
912 
913           en = c[ii + 1][u - ii - 1] +
914                vrna_E_ext_stem(type, mm5, -1, P) +
915                f3[u + 1];
916 
917           if (sc) {
918             if (sc->energy_up)
919               en += sc->energy_up[ii][1];
920 
921             if (sc->f)
922               en += sc->f(ii, length, u, u + 1, VRNA_DECOMP_EXT_STEM_EXT1, sc->data);
923           }
924 
925           if (fij == en) {
926             *i                            = ii + 1;
927             *j                            = u;
928             *k                            = u + 1;
929             bp_stack[++(*stack_count)].i  = ii + 1;
930             bp_stack[(*stack_count)].j    = u;
931             return 1;
932           }
933         }
934 
935         if (u + 2 <= length) {
936           if (evaluate(ii, length, u, u + 2, VRNA_DECOMP_EXT_STEM_EXT, &hc_dat_local)) {
937             mm3   = S1[u + 1];
938             type  = vrna_get_ptype_window(ii, u, ptype);
939 
940             en = c[ii][u - ii] +
941                  vrna_E_ext_stem(type, -1, mm3, P) +
942                  f3[u + 2];
943 
944             if (sc) {
945               if (sc->energy_up)
946                 en += sc->energy_up[u + 1][1];
947 
948               if (sc->f)
949                 en += sc->f(ii, length, u, u + 2, VRNA_DECOMP_EXT_STEM_EXT, sc->data);
950             }
951 
952             if (fij == en) {
953               *i                            = ii;
954               *j                            = u;
955               *k                            = u + 2;
956               bp_stack[++(*stack_count)].i  = ii;
957               bp_stack[(*stack_count)].j    = u;
958               return 1;
959             }
960           }
961         } else {
962           if (evaluate(ii, length, ii, u, VRNA_DECOMP_EXT_STEM, &hc_dat_local)) {
963             mm3   = (u < length) ? S1[u + 1] : -1;
964             type  = vrna_get_ptype_window(ii, u, ptype);
965 
966             en = c[ii][u - ii] +
967                  vrna_E_ext_stem(type, -1, mm3, P);
968 
969             if (sc) {
970               if ((sc->energy_up) && (u < length))
971                 en += sc->energy_up[u + 1][1];
972 
973               if (sc->f)
974                 en += sc->f(ii, length, ii, u, VRNA_DECOMP_EXT_STEM, sc->data);
975             }
976 
977             if (fij == en) {
978               *i                            = ii;
979               *j                            = u;
980               *k                            = u + 2;
981               bp_stack[++(*stack_count)].i  = ii;
982               bp_stack[(*stack_count)].j    = u;
983               return 1;
984             }
985           }
986         }
987 
988         if (evaluate(ii, length, u, u + 1, VRNA_DECOMP_EXT_STEM_EXT, &hc_dat_local)) {
989           type = vrna_get_ptype_window(ii, u, ptype);
990 
991           en = c[ii][u - ii] +
992                vrna_E_ext_stem(type, -1, -1, P) + f3[u + 1];
993 
994           if (sc)
995             if (sc->f)
996               en += sc->f(ii, length, u, u + 1, VRNA_DECOMP_EXT_STEM_EXT, sc->data);
997 
998           if (fij == en) {
999             *i                            = ii;
1000             *j                            = u;
1001             *k                            = u + 1;
1002             bp_stack[++(*stack_count)].i  = ii;
1003             bp_stack[(*stack_count)].j    = u;
1004             return 1;
1005           }
1006         }
1007       }
1008 
1009       break;
1010   }
1011 
1012   return 0;
1013 }
1014 
1015 
1016 PRIVATE int
BT_ext_loop_f3_comparative(vrna_fold_compound_t * fc,int * k,int maxdist,int * i,int * j,vrna_bp_stack_t * bp_stack,int * stack_count)1017 BT_ext_loop_f3_comparative(vrna_fold_compound_t *fc,
1018                            int                  *k,
1019                            int                  maxdist,
1020                            int                  *i,
1021                            int                  *j,
1022                            vrna_bp_stack_t      *bp_stack,
1023                            int                  *stack_count)
1024 {
1025   short                     **S, **S5, **S3;
1026   unsigned int              type, ss, n_seq, **a2s;
1027   int                       n, fij, cc, fj, ii, u, *f3, **c, **ggg,
1028                             dangle_model, turn, with_gquad;
1029   vrna_param_t              *P;
1030   vrna_md_t                 *md;
1031   vrna_sc_t                 **scs;
1032   vrna_callback_hc_evaluate *evaluate;
1033   struct hc_ext_def_dat     hc_dat_local;
1034 
1035   n             = fc->length;
1036   n_seq         = fc->n_seq;
1037   S             = fc->S;
1038   S5            = fc->S5;   /* S5[s][i] holds next base 5' of i in sequence s */
1039   S3            = fc->S3;   /* Sl[s][i] holds next base 3' of i in sequence s */
1040   a2s           = fc->a2s;
1041   P             = fc->params;
1042   md            = &(P->model_details);
1043   scs           = fc->scs;
1044   f3            = fc->matrices->f3_local;
1045   c             = fc->matrices->c_local;
1046   ggg           = fc->matrices->ggg_local;
1047   dangle_model  = md->dangles;
1048   turn          = md->min_loop_size;
1049   with_gquad    = md->gquad;
1050   evaluate      = prepare_hc_ext_def_window(fc, &hc_dat_local);
1051 
1052   ii = *k;
1053 
1054   /* nibble off unpaired 5' bases */
1055   do {
1056     fij = f3[ii];
1057     fj  = INF;
1058 
1059     if (evaluate(ii, n, ii + 1, n, VRNA_DECOMP_EXT_EXT, &hc_dat_local)) {
1060       fj = f3[ii + 1];
1061       if (scs) {
1062         for (ss = 0; ss < n_seq; ss++)
1063           if (scs[ss]) {
1064             if (scs[ss]->energy_up)
1065               fj += scs[ss]->energy_up[ii][1];
1066 
1067             if (scs[ss]->f)
1068               fj += scs[ss]->f(ii, n, ii + 1, n, VRNA_DECOMP_EXT_EXT, scs[ss]->data);
1069           }
1070       }
1071     }
1072 
1073     if (++ii > maxdist)
1074       break;
1075   } while (fij == fj);
1076   ii--;
1077 
1078   if (ii > maxdist - turn + 1) {
1079     /* no more pairs */
1080     *i  = *j = -1;
1081     *k  = 0;
1082     return 1;
1083   }
1084 
1085   /*
1086    *  must have found a decomposition
1087    *  i is paired. Find pairing partner
1088    */
1089   switch (dangle_model) {
1090     /* no dangles */
1091     case 0:
1092       for (u = maxdist; u > ii + turn; u--) {
1093         if (with_gquad) {
1094           if (fij == ggg[ii][u - ii] + f3[u + 1]) {
1095             *i  = *j = -1;
1096             *k  = u + 1;
1097             return vrna_BT_gquad_mfe(fc, ii, u, bp_stack, stack_count);
1098           }
1099         }
1100 
1101         if (evaluate(ii, n, u, u + 1, VRNA_DECOMP_EXT_STEM_EXT, &hc_dat_local)) {
1102           cc = c[ii][u - ii];
1103           for (ss = 0; ss < n_seq; ss++) {
1104             type  = vrna_get_ptype_md(S[ss][ii], S[ss][u], md);
1105             cc    += vrna_E_ext_stem(type, -1, -1, P);
1106           }
1107 
1108           if (fij == cc + f3[u + 1]) {
1109             *i                            = ii;
1110             *j                            = u;
1111             *k                            = u + 1;
1112             bp_stack[++(*stack_count)].i  = ii;
1113             bp_stack[(*stack_count)].j    = u;
1114             return 1;
1115           }
1116         }
1117       }
1118       break;
1119 
1120     /* dangles on both sides */
1121     case 2:
1122       for (u = maxdist; u > ii + turn; u--) {
1123         if (with_gquad) {
1124           if (fij == ggg[ii][u - ii] + f3[u + 1]) {
1125             *i  = *j = -1;
1126             *k  = u + 1;
1127             return vrna_BT_gquad_mfe(fc, ii, u, bp_stack, stack_count);
1128           }
1129         }
1130 
1131         if (evaluate(ii, n, u, u + 1, VRNA_DECOMP_EXT_STEM_EXT, &hc_dat_local)) {
1132           cc = c[ii][u - ii];
1133           for (ss = 0; ss < n_seq; ss++) {
1134             type  = vrna_get_ptype_md(S[ss][ii], S[ss][u], md);
1135             cc    +=
1136               vrna_E_ext_stem(type, (a2s[ss][ii] > 1) ? S5[ss][ii] : -1,
1137                               (a2s[ss][u] < a2s[ss][n]) ? S3[ss][u] : -1, P);
1138           }
1139 
1140           if (fij == cc + f3[u + 1]) {
1141             *i                            = ii;
1142             *j                            = u;
1143             *k                            = u + 1;
1144             bp_stack[++(*stack_count)].i  = ii;
1145             bp_stack[(*stack_count)].j    = u;
1146             return 1;
1147           }
1148         }
1149       }
1150       break;
1151   }
1152 
1153   return 0;
1154 }
1155 
1156 
1157 PRIVATE int
BT_ext_loop_f3_pp(vrna_fold_compound_t * fc,int * i,int maxj)1158 BT_ext_loop_f3_pp(vrna_fold_compound_t  *fc,
1159                   int                   *i,
1160                   int                   maxj)
1161 {
1162   int j, start;
1163 
1164   j     = -1;
1165   start = *i;
1166 
1167   if (fc) {
1168     char                      **ptype;
1169     short                     *S1;
1170     unsigned int              type;
1171     int                       traced2, length, turn, dangle_model, with_gquad, maxdist, cc,
1172                               **c, **ggg, *f3, fij, ii;
1173     vrna_param_t              *P;
1174     vrna_md_t                 *md;
1175     vrna_sc_t                 *sc;
1176     vrna_callback_hc_evaluate *evaluate;
1177     struct hc_ext_def_dat     hc_dat_local;
1178 #ifdef VRNA_WITH_SVM
1179     int                       zsc_pre_filter;
1180     vrna_zsc_dat_t            zsc_data;
1181 #endif
1182 
1183     length        = fc->length;
1184     S1            = fc->sequence_encoding;
1185     ptype         = fc->ptype_local;
1186     f3            = fc->matrices->f3_local;
1187     c             = fc->matrices->c_local;
1188     ggg           = fc->matrices->ggg_local;
1189     sc            = fc->sc;
1190     P             = fc->params;
1191     md            = &(P->model_details);
1192     turn          = md->min_loop_size;
1193     dangle_model  = md->dangles;
1194     with_gquad    = md->gquad;
1195     maxdist       = MIN2(fc->window_size, maxj);
1196     traced2       = 0;
1197     ii            = start;
1198     evaluate      = prepare_hc_ext_def_window(fc, &hc_dat_local);
1199 #ifdef VRNA_WITH_SVM
1200     zsc_data        = fc->zscore_data;
1201     zsc_pre_filter  = ((zsc_data) &&
1202                        (zsc_data->filter_on) &&
1203                        (zsc_data->pre_filter) &&
1204                        (zsc_data->current_i == ii)) ? 1 : 0;
1205 #endif
1206 
1207     fij = f3[start];
1208 
1209     /* try to nibble off unpaired 5' bases */
1210     if ((sc) && (evaluate(start, length, start + 1, length, VRNA_DECOMP_EXT_EXT, &hc_dat_local))) {
1211       cc = f3[start + 1];
1212 
1213       if (sc->energy_up)
1214         cc += sc->energy_up[start][1];
1215 
1216       if (sc->f)
1217         cc += sc->f(start, length, start + 1, length, VRNA_DECOMP_EXT_EXT, sc->data);
1218 
1219       if (fij == cc)
1220         /* simple 5' unpaired extensions, so we skip this hit */
1221         return 0;
1222     }
1223 
1224     /* get pairing partner j */
1225     switch (dangle_model) {
1226       case 0:
1227         for (j = start + turn; j <= ii + maxdist; j++) {
1228           if (evaluate(start, length, j, j + 1, VRNA_DECOMP_EXT_STEM_EXT, &hc_dat_local)) {
1229 #ifdef VRNA_WITH_SVM
1230             if ((zsc_pre_filter) &&
1231                 (zsc_data->current_z[j] > zsc_data->min_z))
1232               continue;
1233 
1234 #endif
1235             type = vrna_get_ptype_window(start, j, ptype);
1236 
1237             cc = c[start][j - start] +
1238                  vrna_E_ext_stem(type, -1, -1, P) +
1239                  f3[j + 1];
1240 
1241             if ((sc) && (sc->f))
1242               cc += sc->f(start, length, j, j + 1, VRNA_DECOMP_EXT_STEM_EXT, sc->data);
1243 
1244             if (fij == cc) {
1245               traced2 = 1;
1246               break;
1247             }
1248           }
1249 
1250           if (with_gquad) {
1251             cc = ggg[start][j - start] + f3[j + 1];
1252             if (fij == cc) {
1253               traced2 = 1;
1254               break;
1255             }
1256           }
1257         }
1258         break;
1259 
1260       case 2:
1261         for (j = start + turn; j <= ii + maxdist; j++) {
1262           if (evaluate(start, length, j, j + 1, VRNA_DECOMP_EXT_STEM_EXT, &hc_dat_local)) {
1263 #ifdef VRNA_WITH_SVM
1264             if ((zsc_pre_filter) &&
1265                 (zsc_data->current_z[j] > zsc_data->min_z))
1266               continue;
1267 
1268 #endif
1269             type = vrna_get_ptype_window(start, j, ptype);
1270 
1271             cc = c[start][j - start] +
1272                  vrna_E_ext_stem(type, (start > 1) ? S1[start - 1] : -1,
1273                                  (j < length) ? S1[j + 1] : -1, P) +
1274                  f3[j + 1];
1275 
1276             if ((sc) && (sc->f))
1277               cc += sc->f(start, length, j, j + 1, VRNA_DECOMP_EXT_STEM_EXT, sc->data);
1278 
1279             if (fij == cc) {
1280               traced2 = 1;
1281               break;
1282             }
1283           }
1284 
1285           if (with_gquad) {
1286             cc = ggg[start][j - start] +
1287                  f3[j + 1];
1288 
1289             if (fij == cc) {
1290               traced2 = 1;
1291               break;
1292             }
1293           }
1294         }
1295         break;
1296 
1297       default:
1298         for (j = start + turn; j <= ii + maxdist; j++) {
1299           if (evaluate(start, length, j, j + 1, VRNA_DECOMP_EXT_STEM_EXT, &hc_dat_local)) {
1300             type = vrna_get_ptype_window(start, j, ptype);
1301 
1302             cc = c[start][j - start] +
1303                  vrna_E_ext_stem(type, -1, -1, P) +
1304                  f3[j + 1];
1305 
1306             if ((sc) && (sc->f))
1307               cc += sc->f(start, length, j, j + 1, VRNA_DECOMP_EXT_STEM_EXT, sc->data);
1308 
1309             if (fij == cc) {
1310               traced2 = 1;
1311               break;
1312             }
1313           }
1314 
1315           if (j < length) {
1316             if (evaluate(start, length, j, j + 2, VRNA_DECOMP_EXT_STEM_EXT, &hc_dat_local)) {
1317               type  = vrna_get_ptype_window(start, j, ptype);
1318               cc    = c[start][j - start] +
1319                       vrna_E_ext_stem(type, -1, S1[j + 1], P) +
1320                       f3[j + 2];
1321 
1322               if (sc) {
1323                 if (sc->energy_up)
1324                   cc += sc->energy_up[j + 1][1];
1325 
1326                 if (sc->f)
1327                   cc += sc->f(start, length, j, j + 2, VRNA_DECOMP_EXT_STEM_EXT, sc->data);
1328               }
1329 
1330               if (fij == cc) {
1331                 traced2 = 1;
1332                 break;
1333               }
1334             }
1335           }
1336 
1337           if (with_gquad) {
1338             cc = ggg[start][j - start] +
1339                  f3[j + 1];
1340 
1341             if (fij == cc) {
1342               traced2 = 1;
1343               break;
1344             }
1345           }
1346 
1347           if (evaluate(start, length, j, j + 1, VRNA_DECOMP_EXT_STEM_EXT1, &hc_dat_local)) {
1348             type = vrna_get_ptype_window(start + 1, j, ptype);
1349 
1350             cc = c[start + 1][j - (start + 1)] +
1351                  vrna_E_ext_stem(type, S1[start], -1, P) +
1352                  f3[j + 1];
1353 
1354             if (sc) {
1355               if (sc->energy_up)
1356                 cc += sc->energy_up[start][1];
1357 
1358               if (sc->f)
1359                 cc += sc->f(start, length, j, j + 1, VRNA_DECOMP_EXT_STEM_EXT1, sc->data);
1360             }
1361 
1362             if (fij == cc) {
1363               traced2 = 1;
1364               break;
1365             }
1366           }
1367 
1368           if (j < length) {
1369             if (evaluate(start, length, j, j + 2, VRNA_DECOMP_EXT_STEM_EXT1, &hc_dat_local)) {
1370               type = vrna_get_ptype_window(start + 1, j, ptype);
1371 
1372               cc = c[start + 1][j - (start + 1)] +
1373                    vrna_E_ext_stem(type, S1[start], S1[j + 1], P) +
1374                    f3[j + 2];
1375 
1376               if (sc) {
1377                 if (sc->energy_up)
1378                   cc += sc->energy_up[start][1] +
1379                         sc->energy_up[j + 1][1];
1380 
1381                 if (sc->f)
1382                   cc += sc->f(start, length, j, j + 2, VRNA_DECOMP_EXT_STEM_EXT1, sc->data);
1383               }
1384 
1385               if (fij == cc) {
1386                 traced2 = 1;
1387                 break;
1388               }
1389             }
1390           }
1391         }
1392         break;
1393     }
1394 
1395     if (!traced2)
1396       j = -1;
1397   }
1398 
1399   *i = start;
1400   return j;
1401 }
1402 
1403 
1404 PRIVATE int
BT_ext_loop_f3_pp_comparative(vrna_fold_compound_t * fc,int * i,int maxj)1405 BT_ext_loop_f3_pp_comparative(vrna_fold_compound_t  *fc,
1406                               int                   *i,
1407                               int                   maxj)
1408 {
1409   int j, start;
1410 
1411   j = -1;
1412 
1413   if (fc) {
1414     short                     **S, **S5, **S3;
1415     unsigned int              tt, s, n_seq, **a2s;
1416     int                       traced2, length, turn, dangle_model, with_gquad, maxdist, cc, **c,
1417                               **ggg, *f3, fij;
1418     vrna_param_t              *P;
1419     vrna_md_t                 *md;
1420     vrna_sc_t                 **scs;
1421     vrna_callback_hc_evaluate *evaluate;
1422     struct hc_ext_def_dat     hc_dat_local;
1423 
1424     length        = fc->length;
1425     n_seq         = fc->n_seq;
1426     S             = fc->S;
1427     S5            = fc->S5;   /* S5[s][start] holds next base 5' of start in sequence s */
1428     S3            = fc->S3;   /* Sl[s][start] holds next base 3' of start in sequence s */
1429     a2s           = fc->a2s;
1430     f3            = fc->matrices->f3_local;
1431     c             = fc->matrices->c_local;
1432     ggg           = fc->matrices->ggg_local;
1433     scs           = fc->scs;
1434     P             = fc->params;
1435     md            = &(P->model_details);
1436     turn          = md->min_loop_size;
1437     dangle_model  = md->dangles;
1438     with_gquad    = md->gquad;
1439     maxdist       = MIN2(fc->window_size, maxj);
1440     traced2       = 0;
1441     start         = *i;
1442     evaluate      = prepare_hc_ext_def_window(fc, &hc_dat_local);
1443 
1444     fij = f3[start];
1445 
1446     /* try to nibble off unpaired 5' bases */
1447     if ((scs) && (evaluate(start, length, start + 1, length, VRNA_DECOMP_EXT_EXT, &hc_dat_local))) {
1448       cc = f3[start + 1];
1449 
1450       for (s = 0; s < n_seq; s++)
1451         if (scs[s]) {
1452           if (scs[s]->energy_up)
1453             cc += scs[s]->energy_up[start][1];
1454 
1455           if (scs[s]->f)
1456             cc +=scs[s]->f(start,
1457                            length,
1458                            start + 1,
1459                            length,
1460                            VRNA_DECOMP_EXT_EXT,
1461                            scs[s]->data);
1462         }
1463 
1464       if (fij == cc)
1465         /* simple 5' unpaired extensions, so we skip this hit */
1466         return 0;
1467     }
1468 
1469     /* get pairing partner j */
1470     switch (dangle_model) {
1471       case 0:
1472         for (j = start + turn; j <= MIN2(start + maxdist, length - 1); j++) {
1473           if (evaluate(start, length, j, j + 1, VRNA_DECOMP_EXT_STEM_EXT, &hc_dat_local)) {
1474             cc = c[start][j - start] +
1475                  f3[j + 1];
1476 
1477             for (s = 0; s < n_seq; s++) {
1478               tt  = vrna_get_ptype_md(S[s][start], S[s][j], md);
1479               cc  += vrna_E_ext_stem(tt, -1, -1, P);
1480             }
1481 
1482             if (scs) {
1483               for (s = 0; s < n_seq; s++) {
1484                 if ((scs[s]) && (scs[s]->f))
1485                   cc += scs[s]->f(start, length, j, j + 1, VRNA_DECOMP_EXT_STEM_EXT, scs[s]->data);
1486               }
1487             }
1488 
1489             if (fij == cc) {
1490               traced2 = 1;
1491               break;
1492             }
1493           }
1494 
1495           if (with_gquad) {
1496             cc = ggg[start][j - start] +
1497                  f3[j + 1];
1498 
1499             if (fij == cc) {
1500               traced2 = 1;
1501               break;
1502             }
1503           }
1504         }
1505 
1506         if ((!traced2) && (length <= start + maxdist)) {
1507           j = length;
1508           if (evaluate(start, length, start, j, VRNA_DECOMP_EXT_STEM, &hc_dat_local)) {
1509             cc = c[start][j - start];
1510 
1511             for (s = 0; s < n_seq; s++) {
1512               tt  = vrna_get_ptype_md(S[s][start], S[s][j], md);
1513               cc  += vrna_E_ext_stem(tt, -1, -1, P);
1514             }
1515 
1516             if (scs) {
1517               for (s = 0; s < n_seq; s++) {
1518                 if ((scs[s]) && (scs[s]->f))
1519                   cc += scs[s]->f(start, length, start, j, VRNA_DECOMP_EXT_STEM, scs[s]->data);
1520               }
1521             }
1522 
1523             if (fij == cc) {
1524               traced2 = 1;
1525               break;
1526             }
1527           }
1528 
1529           if (with_gquad) {
1530             cc = ggg[start][j - start];
1531 
1532             if (fij == cc) {
1533               traced2 = 1;
1534               break;
1535             }
1536           }
1537         }
1538 
1539         break;
1540 
1541       case 2:
1542         for (j = start + turn; j <= MIN2(start + maxdist, length - 1); j++) {
1543           if (evaluate(start, length, j, j + 1, VRNA_DECOMP_EXT_STEM_EXT, &hc_dat_local)) {
1544             cc = c[start][j - start] +
1545                  f3[j + 1];
1546 
1547             for (s = 0; s < n_seq; s++) {
1548               tt  = vrna_get_ptype_md(S[s][start], S[s][j], md);
1549               cc  +=
1550                 vrna_E_ext_stem(tt,
1551                                 (a2s[s][start] > 1) ? S5[s][start] : -1,
1552                                 (a2s[s][j] < a2s[s][length]) ? S3[s][j] : -1,
1553                                 P);
1554             }
1555 
1556             if (scs) {
1557               for (s = 0; s < n_seq; s++) {
1558                 if ((scs[s]) && (scs[s]->f))
1559                   cc += scs[s]->f(start, length, j, j + 1, VRNA_DECOMP_EXT_STEM_EXT, scs[s]->data);
1560               }
1561             }
1562 
1563             if (fij == cc) {
1564               traced2 = 1;
1565               break;
1566             }
1567           }
1568 
1569           if (with_gquad) {
1570             cc = ggg[start][j - start] +
1571                  f3[j + 1];
1572 
1573             if (fij == cc) {
1574               traced2 = 1;
1575               break;
1576             }
1577           }
1578         }
1579 
1580         if ((!traced2) && (length <= start + maxdist)) {
1581           j = length;
1582           if (evaluate(start, length, start, j, VRNA_DECOMP_EXT_STEM, &hc_dat_local)) {
1583             cc = c[start][j - start];
1584             for (s = 0; s < n_seq; s++) {
1585               tt  = vrna_get_ptype_md(S[s][start], S[s][j], md);
1586               cc  += vrna_E_ext_stem(tt, (a2s[s][start] > 1) ? S5[s][start] :  -1, -1, P);
1587             }
1588 
1589             if (scs) {
1590               for (s = 0; s < n_seq; s++) {
1591                 if ((scs[s]) && (scs[s]->f))
1592                   cc += scs[s]->f(start, length, start, j, VRNA_DECOMP_EXT_STEM, scs[s]->data);
1593               }
1594             }
1595 
1596             if (fij == cc) {
1597               traced2 = 1;
1598               break;
1599             }
1600           }
1601 
1602           if (with_gquad) {
1603             cc = ggg[start][j - start];
1604 
1605             if (fij == cc) {
1606               traced2 = 1;
1607               break;
1608             }
1609           }
1610         }
1611 
1612         break;
1613 
1614 
1615     }
1616 
1617     if (!traced2)
1618       j = -1;
1619   }
1620 
1621   return j;
1622 }
1623