1 /*
2 * minimum free energy
3 * RNA secondary structure prediction
4 *
5 * c Ivo Hofacker, Chrisoph Flamm
6 * original implementation by
7 * Walter Fontana
8 *
9 * Vienna RNA package
10 */
11
12 #ifdef HAVE_CONFIG_H
13 #include "config.h"
14 #endif
15
16 #include <stdio.h>
17 #include <stdlib.h>
18 #include <math.h>
19 #include <ctype.h>
20 #include <string.h>
21 #include "ViennaRNA/utils/basic.h"
22 #include "ViennaRNA/utils/structures.h"
23 #include "ViennaRNA/params/default.h"
24 #include "ViennaRNA/fold_vars.h"
25 #include "ViennaRNA/pair_mat.h"
26 #include "ViennaRNA/params/basic.h"
27 #include "ViennaRNA/snofold.h"
28 #include "ViennaRNA/loops/all.h"
29
30 #ifdef __GNUC__
31 #define INLINE inline
32 #else
33 #define INLINE
34 #endif
35
36 #define PAREN
37
38 #define PUBLIC
39 #define PRIVATE static
40
41 #define STACK_BULGE1 1 /* stacking energies for bulges of size 1 */
42 #define NEW_NINIO 1 /* new asymetry penalty */
43
44 /*@unused@*/
45 PRIVATE void
46 get_arrays(unsigned int size);
47
48
49 /* PRIVATE int stack_energy(int i, const char *string); */
50 PRIVATE void
51 make_ptypes(const short *S,
52 const char *structure);
53
54
55 PRIVATE void
56 encode_seq(const char *sequence);
57
58
59 PRIVATE void
60 backtrack(const char *sequence,
61 int s);
62
63
64 PRIVATE int
65 fill_arrays(const char *sequence,
66 const int max_asymm,
67 const int threshloop,
68 const int min_s2,
69 const int max_s2,
70 const int half_stem,
71 const int max_half_stem);
72
73
74 /*@unused@*/
75
76
77 /* alifold */
78 PRIVATE void
79 alisnoinitialize_fold(const int length);
80
81
82 PRIVATE void
83 make_pscores(int length,
84 const short *const *S,
85 const char *const *AS,
86 int n_seq,
87 const char *structure);
88
89
90 PRIVATE int
91 alifill_arrays(const char **string,
92 const int max_asymm,
93 const int threshloop,
94 const int min_s2,
95 const int max_s2,
96 const int half_stem,
97 const int max_half_stem);
98
99
100 PRIVATE void
101 aliget_arrays(unsigned int size);
102
103
104 PRIVATE short *
105 aliencode_seq(const char *sequence);
106
107
108 PRIVATE int
109 alibacktrack(const char **strings,
110 int s);
111
112
113 #define UNIT 100
114 #define MINPSCORE -2 * UNIT
115 /* end alifold */
116
117 #define MAXSECTORS 500 /* dimension for a backtrack array */
118 #define LOCALITY 0. /* locality parameter for base-pairs */
119
120 #define MIN2(A, B) ((A) < (B) ? (A) : (B))
121 #define MAX2(A, B) ((A) > (B) ? (A) : (B))
122 #define SAME_STRAND(I, J) (((I) >= cut_point) || ((J) < cut_point))
123
124 PRIVATE int *pscore; /* precomputed array of pair types */
125 PRIVATE short **Sali;
126 PRIVATE vrna_param_t *P = NULL;
127
128 PRIVATE int *indx = NULL; /* index for moving in the triangle matrices c[] and fMl[]*/
129
130 PRIVATE int *c = NULL; /* energy array, given that i-j pair */
131 PRIVATE int *cc = NULL; /* linear array for calculating canonical structures */
132 PRIVATE int *cc1 = NULL; /* " " */
133 PRIVATE int *Fmi = NULL; /* holds row i of fML (avoids jumps in memory) */
134 PRIVATE int *DMLi = NULL; /* DMLi[j] holds MIN(fML[i,k]+fML[k+1,j]) */
135 PRIVATE int *DMLi1 = NULL; /* MIN(fML[i+1,k]+fML[k+1,j]) */
136 PRIVATE int *DMLi2 = NULL; /* MIN(fML[i+2,k]+fML[k+1,j]) */
137 PRIVATE char *ptype = NULL; /* precomputed array of pair types */
138 PRIVATE short *S = NULL, *S1 = NULL;
139 PRIVATE int init_length = -1;
140 PRIVATE int *mLoop = NULL; /*contains the minimum of c for a xy range*/
141 PRIVATE folden **foldlist = NULL;
142 PRIVATE folden **foldlist_XS = NULL;
143
144 PRIVATE int *BP = NULL; /* contains the structure constrainsts: BP[i]
145 * -1: | = base must be paired
146 * -2: < = base must be paired with j<i
147 * -3: > = base must be paired with j>i
148 * -4: x = base must not pair
149 * positive int: base is paired with int */
150
151
152 static sect sector[MAXSECTORS]; /* stack of partial structures for backtracking */
153
154 PRIVATE char alpha[] = "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz";
155 /* needed by cofold/eval */
156 /* PRIVATE int cut_in_loop(int i); */
157 /* PRIVATE int min_hairpin = TURN; */
158
159 /* some definitions to take circfold into account... */
160 /* PRIVATE int *fM2 = NULL;*/ /* fM2 = multiloop region with exactly two stems, extending to 3' end */
161 PUBLIC int Fc, FcH, FcI, FcM; /* parts of the exterior loop energies */
162 /*--------------------------------------------------------------------------*/
163
164 void
snoinitialize_fold(const int length)165 snoinitialize_fold(const int length)
166 {
167 unsigned int n;
168
169 if (length < 1)
170 vrna_message_error("snoinitialize_fold: argument must be greater 0");
171
172 if (init_length > 0)
173 snofree_arrays(length);
174
175 get_arrays((unsigned)length);
176 init_length = length;
177 for (n = 1; n <= (unsigned)length; n++)
178 indx[n] = (n * (n - 1)) >> 1; /* n(n-1)/2 */
179
180 snoupdate_fold_params();
181 }
182
183
184 PRIVATE void
alisnoinitialize_fold(const int length)185 alisnoinitialize_fold(const int length)
186 {
187 unsigned int n;
188
189 if (length < 1)
190 vrna_message_error("snoinitialize_fold: argument must be greater 0");
191
192 if (init_length > 0)
193 snofree_arrays(length);
194
195 aliget_arrays((unsigned)length);
196 make_pair_matrix();
197 init_length = length;
198 for (n = 1; n <= (unsigned)length; n++)
199 indx[n] = (n * (n - 1)) >> 1; /* n(n-1)/2 */
200 snoupdate_fold_params();
201 }
202
203
204 /*--------------------------------------------------------------------------*/
205
206 PRIVATE void
get_arrays(unsigned int size)207 get_arrays(unsigned int size)
208 {
209 indx = (int *)vrna_alloc(sizeof(int) * (size + 1));
210 c = (int *)vrna_alloc(sizeof(int) * ((size * (size + 1)) / 2 + 2));
211 mLoop = (int *)vrna_alloc(sizeof(int) * ((size * (size + 1)) / 2 + 2));
212
213 ptype = (char *)vrna_alloc(sizeof(char) * ((size * (size + 1)) / 2 + 2));
214 cc = (int *)vrna_alloc(sizeof(int) * (size + 2));
215 cc1 = (int *)vrna_alloc(sizeof(int) * (size + 2));
216 Fmi = (int *)vrna_alloc(sizeof(int) * (size + 1));
217 DMLi = (int *)vrna_alloc(sizeof(int) * (size + 1));
218 DMLi1 = (int *)vrna_alloc(sizeof(int) * (size + 1));
219 DMLi2 = (int *)vrna_alloc(sizeof(int) * (size + 1));
220 if (base_pair)
221 free(base_pair);
222
223 base_pair = (vrna_bp_stack_t *)vrna_alloc(sizeof(vrna_bp_stack_t) * (1 + size / 2));
224 /* extra array(s) for circfold() */
225 }
226
227
228 PRIVATE void
aliget_arrays(unsigned int size)229 aliget_arrays(unsigned int size)
230 {
231 indx = (int *)vrna_alloc(sizeof(int) * (size + 1));
232 c = (int *)vrna_alloc(sizeof(int) * ((size * (size + 1)) / 2 + 2));
233 mLoop = (int *)vrna_alloc(sizeof(int) * ((size * (size + 1)) / 2 + 2));
234 pscore = (int *)vrna_alloc(sizeof(int) * ((size * (size + 1)) / 2 + 2));
235 ptype = (char *)vrna_alloc(sizeof(char) * ((size * (size + 1)) / 2 + 2));
236 cc = (int *)vrna_alloc(sizeof(int) * (size + 2));
237 cc1 = (int *)vrna_alloc(sizeof(int) * (size + 2));
238 Fmi = (int *)vrna_alloc(sizeof(int) * (size + 1));
239 DMLi = (int *)vrna_alloc(sizeof(int) * (size + 1));
240 DMLi1 = (int *)vrna_alloc(sizeof(int) * (size + 1));
241 DMLi2 = (int *)vrna_alloc(sizeof(int) * (size + 1));
242 if (base_pair)
243 free(base_pair);
244
245 base_pair = (vrna_bp_stack_t *)vrna_alloc(sizeof(vrna_bp_stack_t) * (1 + size / 2));
246 /* extra array(s) for circfold() */
247 }
248
249
250 /*--------------------------------------------------------------------------*/
251
252
253 void
snofree_arrays(const int length)254 snofree_arrays(const int length)
255 {
256 free(indx);
257 free(c);
258 free(cc);
259 free(cc1);
260 free(ptype);
261 free(mLoop);
262 int i;
263 for (i = length; i > -1; i--) {
264 while (foldlist[i] != NULL) {
265 folden *n = foldlist[i];
266 foldlist[i] = foldlist[i]->next;
267 free(n);
268 }
269 free(foldlist[i]);
270 }
271 free(foldlist);
272 for (i = length; i > -1; i--) {
273 while (foldlist_XS[i] != NULL) {
274 folden *n = foldlist_XS[i];
275 foldlist_XS[i] = foldlist_XS[i]->next;
276 free(n);
277 }
278 free(foldlist_XS[i]);
279 }
280 free(foldlist_XS);
281 free(base_pair);
282 base_pair = NULL;
283 free(Fmi);
284 free(DMLi);
285 free(DMLi1);
286 free(DMLi2);
287 free(BP);
288 init_length = 0;
289 }
290
291
292 void
alisnofree_arrays(const int length)293 alisnofree_arrays(const int length)
294 {
295 free(indx);
296 free(c);
297 free(cc);
298 free(cc1);
299 free(ptype);
300 free(mLoop);
301 free(pscore);
302 int i;
303 for (i = length - 1; i > -1; i--) {
304 while (foldlist[i] != NULL) {
305 folden *n = foldlist[i];
306 foldlist[i] = foldlist[i]->next;
307 free(n);
308 }
309 free(foldlist[i]);
310 }
311 free(foldlist);
312 free(base_pair);
313 base_pair = NULL;
314 free(Fmi);
315 free(DMLi);
316 free(DMLi1);
317 free(DMLi2);
318 free(BP);
319 init_length = 0;
320 }
321
322
323 /*--------------------------------------------------------------------------*/
324
325 void
snoexport_fold_arrays(int ** indx_p,int ** mLoop_p,int ** cLoop,folden *** fold_p,folden *** fold_p_XS)326 snoexport_fold_arrays(int **indx_p,
327 int **mLoop_p,
328 int **cLoop,
329 folden ***fold_p,
330 folden ***fold_p_XS)
331 {
332 /* make the DP arrays available to routines such as subopt() */
333 *indx_p = indx;
334 *mLoop_p = mLoop;
335 *cLoop = c;
336 *fold_p = foldlist;
337 *fold_p_XS = foldlist_XS;
338 }
339
340
341 /* void alisnoexport_fold_arrays(int **indx_p, int **mLoop_p, int **cLoop, folden ***fold_p, int **pscores) { */
342 /* /\* make the DP arrays available to routines such as subopt() *\/ */
343 /* *indx_p = indx; *mLoop_p = mLoop; */
344 /* *cLoop = c; *fold_p = foldlist; */
345 /* *pscores=pscore; */
346 /* } */
347
348 /*--------------------------------------------------------------------------*/
349
350
351 int
snofold(const char * string,char * structure,const int max_assym,const int threshloop,const int min_s2,const int max_s2,const int half_stem,const int max_half_stem)352 snofold(const char *string,
353 char *structure,
354 const int max_assym,
355 const int threshloop,
356 const int min_s2,
357 const int max_s2,
358 const int half_stem,
359 const int max_half_stem)
360 {
361 int length, energy, bonus, bonus_cnt, s;
362
363 /* Variable initialization */
364 bonus = 0;
365 bonus_cnt = 0;
366 s = 0;
367 length = (int)strlen(string);
368
369 S = encode_sequence(string, 0);
370 S1 = encode_sequence(string, 1);
371
372
373 /* structure = (char *) vrna_alloc((unsigned) length+1); */
374
375 if (length > init_length)
376 snoinitialize_fold(length);
377 else if (fabs(P->temperature - temperature) > 1e-6)
378 snoupdate_fold_params();
379
380 /* encode_seq(string); */
381 BP = (int *)vrna_alloc(sizeof(int) * (length + 2));
382 make_ptypes(S, structure);
383 energy = fill_arrays(string, max_assym, threshloop, min_s2, max_s2, half_stem, max_half_stem);
384 backtrack(string, s);
385
386 free(structure);
387 free(S);
388 free(S1); /* free(BP); */
389 return energy;
390 }
391
392
393 PRIVATE void
make_pscores(int n,const short * const * S,const char * const * AS,int n_seq,const char * structure)394 make_pscores(int n,
395 const short *const *S,
396 const char *const *AS,
397 int n_seq,
398 const char *structure)
399 {
400 /* calculate co-variance bonus for each pair depending on */
401 /* compensatory/consistent mutations and incompatible seqs */
402 /* should be 0 for conserved pairs, >0 for good pairs */
403 #define NONE -10000 /* score for forbidden pairs */
404 int i, j, k, l, s, score;
405 int dm[7][7] = { { 0, 0, 0, 0, 0, 0, 0 }, /* hamming distance between pairs */
406 { 0, 0, 2, 2, 1, 2, 2 } /* CG */,
407 { 0, 2, 0, 1, 2, 2, 2 } /* GC */,
408 { 0, 2, 1, 0, 2, 1, 2 } /* GU */,
409 { 0, 1, 2, 2, 0, 2, 1 } /* UG */,
410 { 0, 2, 2, 1, 2, 0, 2 } /* AU */,
411 { 0, 2, 2, 2, 1, 2, 0 } /* UA */ };
412
413 for (i = 1; i < n; i++) {
414 for (j = i + 1; (j < i + TURN + 1) && (j <= n); j++)
415 pscore[indx[j] + i] = NONE;
416 for (j = i + TURN + 1; j <= n; j++) {
417 int pfreq[8] = {
418 0, 0, 0, 0, 0, 0, 0, 0
419 };
420 for (s = 0; s < n_seq; s++) {
421 int type;
422 if (Sali[s][i] == 0 && Sali[s][j] == 0) {
423 type = 7; /* gap-gap */
424 } else {
425 if ((AS[s][i] == '~') || (AS[s][j] == '~'))
426 type = 7;
427 else
428 type = pair[Sali[s][i]][Sali[s][j]];
429 }
430
431 pfreq[type]++;
432 }
433 if (pfreq[0] * 2 > n_seq) {
434 pscore[indx[j] + i] = NONE;
435 continue;
436 }
437
438 for (k = 1, score = 0; k <= 6; k++) /* ignore pairtype 7 (gap-gap) */
439 for (l = k + 1; l <= 6; l++)
440 /* scores for replacements between pairtypes */
441 /* consistent or compensatory mutations score 1 or 2 */
442 score += pfreq[k] * pfreq[l] * dm[k][l];
443 /* counter examples score -1, gap-gap scores -0.25 */
444 pscore[indx[j] + i] = cv_fact *
445 ((UNIT * score) / n_seq - nc_fact * UNIT *
446 (pfreq[0] + pfreq[7] * 0.25));
447 }
448 }
449
450 if (noLonelyPairs) {
451 /* remove unwanted pairs */
452 for (k = 1; k < n - TURN - 1; k++)
453 for (l = 1; l <= 2; l++) {
454 int type, ntype = 0, otype = 0;
455 i = k;
456 j = i + TURN + l;
457 type = pscore[indx[j] + i];
458 while ((i >= 1) && (j <= n)) {
459 if ((i > 1) && (j < n))
460 ntype = pscore[indx[j + 1] + i - 1];
461
462 if ((otype < -4 * UNIT) && (ntype < -4 * UNIT)) /* worse than 2 counterex */
463 pscore[indx[j] + i] = NONE; /* i.j can only form isolated pairs */
464
465 otype = type;
466 type = ntype;
467 i--;
468 j++;
469 }
470 }
471 }
472
473 if (fold_constrained && (structure != NULL)) {
474 int psij, hx, hx2, *stack, *stack2;
475 stack = (int *)vrna_alloc(sizeof(int) * (n + 1));
476 stack2 = (int *)vrna_alloc(sizeof(int) * (n + 1));
477
478 for (hx = hx2 = 0, j = 1; j <= n; j++) {
479 switch (structure[j - 1]) {
480 case 'x': /* can't pair */
481 for (l = 1; l < j - TURN; l++)
482 pscore[indx[j] + l] = NONE;
483 for (l = j + TURN + 1; l <= n; l++)
484 pscore[indx[l] + j] = NONE;
485 break;
486 case '(':
487 stack[hx++] = j;
488 /* fallthrough */
489 case '[':
490 stack2[hx2++] = j;
491 /* fallthrough */
492 case '<': /* pairs upstream */
493 for (l = 1; l < j - TURN; l++)
494 pscore[indx[j] + l] = NONE;
495 break;
496 case ']':
497 if (hx2 <= 0)
498 vrna_message_error("unbalanced brackets in constraints\n%s", structure);
499
500 i = stack2[--hx2];
501 pscore[indx[j] + i] = NONE;
502 break;
503 case ')':
504 if (hx <= 0)
505 vrna_message_error("unbalanced brackets in constraints\n%s", structure);
506
507 i = stack[--hx];
508 psij = pscore[indx[j] + i]; /* store for later */
509 for (k = j; k <= n; k++)
510 for (l = i; l <= j; l++)
511 pscore[indx[k] + l] = NONE;
512 for (l = i; l <= j; l++)
513 for (k = 1; k <= i; k++)
514 pscore[indx[l] + k] = NONE;
515 for (k = i + 1; k < j; k++)
516 pscore[indx[k] + i] = pscore[indx[j] + k] = NONE;
517 pscore[indx[j] + i] = (psij > 0) ? psij : 0;
518 /* fallthrough */
519 case '>': /* pairs downstream */
520 for (l = j + TURN + 1; l <= n; l++)
521 pscore[indx[l] + j] = NONE;
522 break;
523 }
524 }
525 if (hx != 0)
526 vrna_message_error("unbalanced brackets in constraint string\n%s", structure);
527
528 free(stack);
529 free(stack2);
530 }
531 }
532
533
534 float
alisnofold(const char ** strings,const int max_assym,const int threshloop,const int min_s2,const int max_s2,const int half_stem,const int max_half_stem)535 alisnofold(const char **strings,
536 const int max_assym,
537 const int threshloop,
538 const int min_s2,
539 const int max_s2,
540 const int half_stem,
541 const int max_half_stem)
542 {
543 int s, n_seq, length, energy;
544 char *structure;
545
546 length = (int)strlen(strings[0]);
547 /* structure = (char *) vrna_alloc((unsigned) length+1); */
548 structure = NULL;
549 if (length > init_length)
550 alisnoinitialize_fold(length);
551
552 if (fabs(P->temperature - temperature) > 1e-6)
553 snoupdate_fold_params();
554
555 for (s = 0; strings[s] != NULL; s++) ;
556 n_seq = s;
557 Sali = (short **)vrna_alloc(n_seq * sizeof(short *));
558 for (s = 0; s < n_seq; s++) {
559 if (strlen(strings[s]) != length)
560 vrna_message_error("uneqal seqence lengths");
561
562 Sali[s] = aliencode_seq(strings[s]);
563 }
564 make_pscores(length, (const short **)Sali, (const char *const *)strings, n_seq, structure);
565 energy = alifill_arrays(strings, max_assym, threshloop, min_s2, max_s2, half_stem, max_half_stem);
566 alibacktrack((const char **)strings, 0);
567 for (s = 0; s < n_seq; s++)
568 free(Sali[s]);
569 free(Sali);
570 /* free(structure); */
571 /* free(S)*/;
572 free(S1); /* free(BP); */
573 return (float)energy / 100.;
574 }
575
576
577 PRIVATE int
alifill_arrays(const char ** strings,const int max_asymm,const int threshloop,const int min_s2,const int max_s2,const int half_stem,const int max_half_stem)578 alifill_arrays(const char **strings,
579 const int max_asymm,
580 const int threshloop,
581 const int min_s2,
582 const int max_s2,
583 const int half_stem,
584 const int max_half_stem)
585 {
586 int i, j, length, energy;
587 /* int decomp, new_fML; */
588 int *type, type_2;
589 int bonus, n_seq, s;
590
591
592 for (n_seq = 0; strings[n_seq] != NULL; n_seq++) ;
593 type = (int *)vrna_alloc(n_seq * sizeof(int));
594 length = strlen(strings[0]);
595 bonus = 0;
596 /* max_separation = (int) ((1.-LOCALITY)*(double)(length-2));*/ /* not in use */
597
598 /* for (i=(j>TURN?(j-TURN):1); i<j; i++) { */
599 /* } */
600 for (i = (length) - TURN - 1; i >= 1; i--) {
601 /* i,j in [1..length] */
602 for (j = i + TURN + 1; j <= length; j++) {
603 int p, q, ij, psc;
604 ij = indx[j] + i;
605 for (s = 0; s < n_seq; s++) {
606 type[s] = pair[Sali[s][i]][Sali[s][j]];
607 if (type[s] == 0)
608 type[s] = 7;
609 }
610 psc = pscore[indx[j] + i];
611 if (psc >= MINPSCORE) {
612 /* we have a pair */
613 int new_c = 0, stackEnergy = INF; /* seems that new_c immer den minimum von cij enthaelt */
614 /* hairpin ----------------------------------------------*/
615
616 for (new_c = s = 0; s < n_seq; s++)
617 new_c += E_Hairpin(j - i - 1,
618 type[s],
619 Sali[s][i + 1],
620 Sali[s][j - 1],
621 strings[s] + i - 1,
622 P);
623 /*--------------------------------------------------------
624 * check for elementary structures involving more than one
625 * closing pair (interior loop).
626 * --------------------------------------------------------*/
627
628 for (p = i + 1; p <= MIN2(j - 2 - TURN, i + MAXLOOP + 1); p++) {
629 int minq = j - i + p - MAXLOOP - 2;
630 if (minq < p + 1 + TURN)
631 minq = p + 1 + TURN;
632
633 for (q = minq; q < j; q++) {
634 if (pscore[indx[q] + p] < MINPSCORE)
635 continue;
636
637 if (abs((p - i) - (j - q)) > max_asymm)
638 continue;
639
640 for (energy = s = 0; s < n_seq; s++) {
641 type_2 = pair[Sali[s][q]][Sali[s][p]]; /* q,p not p,q! */
642 if (type_2 == 0)
643 type_2 = 7;
644
645 energy += E_IntLoop(p - i - 1, j - q - 1, type[s], type_2,
646 Sali[s][i + 1], Sali[s][j - 1],
647 Sali[s][p - 1], Sali[s][q + 1], P);
648 }
649 new_c = MIN2(energy + c[indx[q] + p], new_c);
650 if ((p == i + 1) && (j == q + 1))
651 stackEnergy = energy; /* remember stack energy */
652 } /* end q-loop */
653 } /* end p-loop */
654
655 /* coaxial stacking of (i.j) with (i+1.k) or (k+1.j-1) */
656
657 new_c = MIN2(new_c, cc1[j - 1] + stackEnergy);
658 cc[j] = new_c - psc; /* add covariance bonnus/penalty */
659 c[ij] = cc[j];
660 } /* end >> if (pair) << */
661 else {
662 c[ij] = INF;
663 }
664
665 /* done with c[i,j], now compute fML[i,j] */
666 /* free ends ? -----------------------------------------*/
667 }
668
669 {
670 int *FF; /* rotate the auxilliary arrays */
671 FF = DMLi2;
672 DMLi2 = DMLi1;
673 DMLi1 = DMLi;
674 DMLi = FF;
675 FF = cc1;
676 cc1 = cc;
677 cc = FF;
678 for (j = 1; j <= length; j++)
679 cc[j] = Fmi[j] = DMLi[j] = INF;
680 }
681 }
682 foldlist = (folden **)vrna_alloc((length) * sizeof(folden *));
683
684 for (i = 0; i < length; i++) {
685 foldlist[i] = (folden *)vrna_alloc(sizeof(folden));
686 foldlist[i]->next = NULL;
687 foldlist[i]->k = INF + 1;
688 foldlist[i]->energy = INF;
689 }
690 folden *head; /* we save the stem loop information in a list like structure */
691
692 for (i = length - TURN - 1; i >= 1; i--) {
693 /* i,j in [1..length] */
694 int max_k, min_k;
695 max_k = MIN2(length - min_s2, i + max_half_stem + 1);
696 min_k = MAX2(i + half_stem + 1, length - max_s2);
697 for (j = i + TURN + 1; j <= length; j++) {
698 int ij, a, b;
699 ij = indx[j] + i;
700 for (a = 0; a < MISMATCH; a++)
701 for (b = 0; b < MISMATCH; b++)
702 mLoop[ij] = MIN2(mLoop[ij], c[indx[j - a] + i + b]);
703 if (mLoop[ij] >= n_seq * threshloop) {
704 mLoop[ij] = INF;
705 } else {
706 if (j >= min_k - 1 && j < max_k) {
707 /* comment if out to recover the known behaviour */
708 head = (folden *)vrna_alloc(sizeof(folden));
709 head->k = j;
710 head->energy = mLoop[ij];
711 head->next = foldlist[i];
712 foldlist[i] = head;
713 }
714 }
715 }
716 }
717 free(type);
718 return mLoop[indx[length] + 1];/* mLoop; */
719 }
720
721
722 PRIVATE int
alibacktrack(const char ** strings,int s)723 alibacktrack(const char **strings,
724 int s)
725 {
726 /*------------------------------------------------------------------
727 * trace back through the "c", "f5" and "fML" arrays to get the
728 * base pairing list. No search for equivalent structures is done.
729 * This is fast, since only few structure elements are recalculated.
730 * ------------------------------------------------------------------*/
731
732 /* normally s=0.
733 * If s>0 then s items have been already pushed onto the sector stack */
734 int i, j, length, energy; /* , new; */
735 int type_2;
736 int bonus, n_seq, *type;
737 int b = 0, cov_en = 0;
738
739 length = strlen(strings[0]);
740 for (n_seq = 0; strings[n_seq] != NULL; n_seq++) ;
741 type = (int *)vrna_alloc(n_seq * sizeof(int));
742 if (s == 0) {
743 sector[++s].i = 1;
744 sector[s].j = length;
745 sector[s].ml = 2;
746 }
747
748 while (s > 0) {
749 int ml, ss, cij, traced, i1, j1, p, q;
750 int canonical = 1; /* (i,j) closes a canonical structure */
751 i = sector[s].i;
752 j = sector[s].j;
753 ml = sector[s--].ml; /* ml is a flag indicating if backtracking is to
754 * occur in the fML- (1) or in the f-array (0) */
755 if (ml == 2) {
756 base_pair[++b].i = i;
757 base_pair[b].j = j;
758 goto repeat1;
759 }
760
761 if (j < i + TURN + 1)
762 continue; /* no more pairs in this interval */
763
764 repeat1:
765
766 /*----- begin of "repeat:" -----*/
767 if (canonical)
768 cij = c[indx[j] + i];
769
770 for (ss = 0; ss < n_seq; ss++) {
771 type[ss] = pair[Sali[ss][i]][Sali[ss][j]];
772 if (type[ss] == 0)
773 type[ss] = 7;
774 }
775 bonus = 0;
776
777 if (noLonelyPairs) {
778 if (cij == c[indx[j] + i]) {
779 /* (i.j) closes canonical structures, thus
780 * (i+1.j-1) must be a pair */
781 for (ss = 0; ss < n_seq; ss++) {
782 type_2 = pair[Sali[ss][j - 1]][Sali[ss][i + 1]]; /* j,i not i,j */
783 if (type_2 == 0)
784 type_2 = 7;
785
786 cij -= P->stack[type[ss]][type_2];
787 }
788 cij += pscore[indx[j] + i];
789 base_pair[++b].i = i + 1;
790 base_pair[b].j = j - 1;
791 cov_en += pscore[indx[j - 1] + i + 1];
792 i++;
793 j--;
794 canonical = 0;
795 goto repeat1;
796 }
797 }
798
799 canonical = 1;
800 cij += pscore[indx[j] + i];
801 {
802 int cc = 0;
803 for (ss = 0; ss < n_seq; ss++)
804 cc += E_Hairpin(j - i - 1,
805 type[ss],
806 Sali[ss][i + 1],
807 Sali[ss][j - 1],
808 strings[ss] + i - 1,
809 P);
810 if (cij == cc) /* found hairpin */
811 continue;
812 }
813 for (p = i + 1; p <= MIN2(j - 2 - TURN, i + MAXLOOP + 1); p++) {
814 int minq;
815 minq = j - i + p - MAXLOOP - 2;
816 if (minq < p + 1 + TURN)
817 minq = p + 1 + TURN;
818
819 for (q = j - 1; q >= minq; q--) {
820 for (ss = energy = 0; ss < n_seq; ss++) {
821 type_2 = pair[Sali[ss][q]][Sali[ss][p]]; /* q,p not p,q */
822 if (type_2 == 0)
823 type_2 = 7;
824
825 energy += E_IntLoop(p - i - 1, j - q - 1, type[ss], type_2,
826 Sali[ss][i + 1], Sali[ss][j - 1],
827 Sali[ss][p - 1], Sali[ss][q + 1], P);
828 }
829 traced = (cij == energy + c[indx[q] + p]);
830 if (traced) {
831 base_pair[++b].i = p;
832 base_pair[b].j = q;
833 cov_en += pscore[indx[q] + p];
834 i = p, j = q;
835 goto repeat1;
836 }
837 }
838 }
839
840 /* end of repeat: --------------------------------------------------*/
841
842 /* (i.j) must close a multi-loop */
843 /* tt = rtype[type]; */
844 /* mm = bonus+P->MLclosing+P->MLintern[tt]; */
845 /* d5 = P->dangle5[tt][S1[j-1]]; */
846 /* d3 = P->dangle3[tt][S1[i+1]]; */
847 i1 = i + 1;
848 j1 = j - 1;
849 sector[s + 1].ml = sector[s + 2].ml = 1;
850
851 /* if (k<=j-3-TURN) { /\\* found the decomposition *\\/ *\/ */
852 /* sector[++s].i = i1; */
853 /* sector[s].j = k; */
854 /* sector[++s].i = k+1; */
855 /* sector[s].j = j1; */
856 /* } /\* else { *\/ */
857 /* vrna_message_error("backtracking failed in repeat"); */
858 /* } */
859 }
860 base_pair[0].i = b; /* save the total number of base pairs */
861 free(type);
862 return cov_en;
863 }
864
865
866 PRIVATE int
fill_arrays(const char * string,const int max_asymm,const int threshloop,const int min_s2,const int max_s2,const int half_stem,const int max_half_stem)867 fill_arrays(const char *string,
868 const int max_asymm,
869 const int threshloop,
870 const int min_s2,
871 const int max_s2,
872 const int half_stem,
873 const int max_half_stem)
874 {
875 int i, j, length, energy;
876 /* int decomp;*/ /*, new_fML; */
877 int no_close, type, type_2;
878 int bonus;
879 int min_c;
880
881 min_c = INF;
882 length = (int)strlen(string);
883 bonus = 0;
884 /* max_separation = (int) ((1.-LOCALITY)*(double)(length-2)); */ /* not in use */
885
886
887 for (i = length - TURN - 1; i >= 1; i--) {
888 /* i,j in [1..length] */
889 /* printf("i=%d\t",i); */
890 for (j = i + TURN + 1; j <= length; j++) {
891 /* printf("j=%d,",j); */
892 int p, q, ij;
893 ij = indx[j] + i;
894 type = ptype[ij];
895 bonus = 0;
896 energy = INF;
897
898 if ((BP[i] == j) || (BP[i] == -1) || (BP[i] == -2))
899 bonus -= BONUS;
900
901 if ((BP[j] == -1) || (BP[j] == -3))
902 bonus -= BONUS;
903
904 if ((BP[i] == -4) || (BP[j] == -4))
905 type = 0;
906
907 no_close = (((type == 3) || (type == 4)) && no_closingGU);
908
909 /* if (j-i-1 > max_separation) type = 0; */ /* forces locality degree */
910
911 if (type) {
912 /* we have a pair */
913 int new_c = 0, stackEnergy = INF; /* seems that new_c immer den minimum von cij enthaelt */
914 /* hairpin ----------------------------------------------*/
915
916 if (no_close)
917 new_c = FORBIDDEN;
918 else
919 new_c = E_Hairpin(j - i - 1, type, S1[i + 1], S1[j - 1], string + i - 1, P); /* computes hair pin structure for subsequence i...j */
920
921 /*--------------------------------------------------------
922 * check for elementary structures involving more than one
923 * closing pair (interior loop).
924 * --------------------------------------------------------*/
925
926 for (p = i + 1; p <= MIN2(j - 2 - TURN, i + MAXLOOP + 1); p++) {
927 int minq = j - i + p - MAXLOOP - 2;
928 if (minq < p + 1 + TURN)
929 minq = p + 1 + TURN;
930
931 for (q = minq; q < j; q++) {
932 if (abs((p - i) - (j - q)) > max_asymm)
933 continue;
934
935 type_2 = ptype[indx[q] + p];
936
937 if (type_2 == 0)
938 continue;
939
940 type_2 = rtype[type_2];
941
942 if (no_closingGU)
943 if (no_close || (type_2 == 3) || (type_2 == 4))
944 if ((p > i + 1) || (q < j - 1))
945 continue;
946
947 /* continue unless stack */
948
949 energy = E_IntLoop(p - i - 1, j - q - 1, type, type_2,
950 S1[i + 1], S1[j - 1], S1[p - 1], S1[q + 1], P);
951 new_c = MIN2(energy + c[indx[q] + p], new_c);
952 if ((p == i + 1) && (j == q + 1))
953 stackEnergy = energy; /* remember stack energy */
954 } /* end q-loop */
955 } /* end p-loop */
956
957
958 /* coaxial stacking of (i.j) with (i+1.k) or (k+1.j-1) */
959
960
961 new_c = MIN2(new_c, cc1[j - 1] + stackEnergy);
962 cc[j] = new_c;
963 c[ij] = new_c;
964 /* min_c=MIN2(min_c, c[ij]); */
965 } /* end >> if (pair) << */
966 else {
967 c[ij] = INF;
968 }
969
970 /* done with c[i,j], now compute fML[i,j] */
971 /* free ends ? -----------------------------------------*/
972 }
973
974 {
975 int *FF; /* rotate the auxilliary arrays */
976 FF = DMLi2;
977 DMLi2 = DMLi1;
978 DMLi1 = DMLi;
979 DMLi = FF;
980 FF = cc1;
981 cc1 = cc;
982 cc = FF;
983 for (j = 1; j <= length; j++)
984 cc[j] = Fmi[j] = DMLi[j] = INF;
985 }
986 }
987 foldlist = (folden **)vrna_alloc((length + 1) * sizeof(folden *));
988 foldlist_XS = (folden **)vrna_alloc((length + 1) * sizeof(folden *));
989 /* linked list initialization*/
990 for (i = 0; i <= length; i++) {
991 foldlist[i] = (folden *)vrna_alloc(sizeof(folden));
992 foldlist[i]->next = NULL;
993 foldlist[i]->k = INF + 1;
994 foldlist[i]->energy = INF;
995 foldlist_XS[i] = (folden *)vrna_alloc(sizeof(folden));
996 foldlist_XS[i]->next = NULL;
997 foldlist_XS[i]->k = INF + 1;
998 foldlist_XS[i]->energy = INF;
999 }
1000 folden *head; /* we save the stem loop information in a list like structure */
1001 folden *head_XS;
1002 for (i = length - TURN - 1; i >= 1; i--) {
1003 /* i,j in [1..length] */
1004 int max_k, min_k;
1005 max_k = MIN2(length - min_s2, i + max_half_stem + 1);
1006 min_k = MAX2(i + half_stem + 1, length - max_s2);
1007
1008
1009 for (j = i + TURN + 1; j <= length; j++) {
1010 int ij, a, b;
1011 ij = indx[j] + i;
1012 for (a = 0; a < MISMATCH; a++) {
1013 for (b = 0; b < MISMATCH; b++) {
1014 mLoop[ij] = MIN2(mLoop[ij], c[indx[j - a] + i + b]);
1015 /* #mLoop[ij]=MIN2(mLoop[ij], c[indx[j-2]+i]); */
1016 /* #mLoop[ij]=MIN2(mLoop[ij], c[indx[j]+i+1]); */
1017 /* #mLoop[ij]=MIN2(mLoop[ij], c[indx[j-1]+i+1]); */
1018 /* #mLoop[ij]=MIN2(mLoop[ij], c[indx[j-2]+i+1]); */
1019 /* #mLoop[ij]=MIN2(mLoop[ij], c[indx[j]+i+2]); */
1020 /* #mLoop[ij]=MIN2(mLoop[ij], c[indx[j-1]+i+2]); */
1021 /* #mLoop[ij]=MIN2(mLoop[ij], c[indx[j-2]+i+2]); */
1022 }
1023 }
1024 min_c = MIN2(mLoop[ij], min_c);
1025
1026 if (mLoop[ij] >= threshloop) {
1027 mLoop[ij] = INF;
1028 } else {
1029 if (j >= min_k - 1 && j <= max_k) {
1030 /* comment if out to recover the known behaviour */
1031 head = (folden *)vrna_alloc(sizeof(folden));
1032 head->k = j;
1033 head->energy = mLoop[ij];
1034 head->next = foldlist[i];
1035 foldlist[i] = head;
1036 head_XS = (folden *)vrna_alloc(sizeof(folden));
1037 head_XS->k = i;
1038 head_XS->energy = mLoop[ij];
1039 head_XS->next = foldlist_XS[j];
1040 foldlist_XS[j] = head_XS;
1041 }
1042 }
1043 }
1044 }
1045 /* int count=0; */
1046 /* for(i=0; i< length; i++){ */
1047 /* folden *temp; */
1048 /* temp = foldlist[i]; */
1049 /* while(temp->next){ */
1050 /* count++; */
1051 /* printf("count %d: i%d j%d energy %d \n", count, i, temp->k, temp->energy); */
1052 /* temp=temp->next; */
1053 /* } */
1054 /* } */
1055 /* printf("Count %d \n", count); */
1056 /* count=0; */
1057 /* for(i=length-1; i>=0; i--){ */
1058 /* folden *temp; */
1059 /* temp = foldlist_XS[i]; */
1060 /* while(temp->next){ */
1061 /* count++; */
1062 /* printf("count %d: i%d j%d energy %d \n", count, temp->k,i, temp->energy); */
1063 /* temp=temp->next; */
1064 /* } */
1065 /* } */
1066 /* printf("Count %d \n", count); */
1067 /* return mLoop[indx[length]+1]; */ /* mLoop; */
1068 return min_c;
1069 /* printf("\nmin_array = %d\n", min_c); */
1070 /* return f5[length]; */
1071 }
1072
1073
1074 PRIVATE void
backtrack(const char * string,int s)1075 backtrack(const char *string,
1076 int s)
1077 {
1078 /*------------------------------------------------------------------
1079 * trace back through the "c", "f5" and "fML" arrays to get the
1080 * base pairing list. No search for equivalent structures is done.
1081 * This is fast, since only few structure elements are recalculated.
1082 * ------------------------------------------------------------------*/
1083
1084 /* normally s=0.
1085 * If s>0 then s items have been already pushed onto the sector stack */
1086 int i, j, /*k,*/ length, energy, new;
1087 int no_close, type, type_2; /* , tt; */
1088 int bonus;
1089 int b = 0;
1090
1091 length = strlen(string);
1092 if (s == 0) {
1093 sector[++s].i = 1;
1094 sector[s].j = length;
1095 sector[s].ml = 2;
1096 }
1097
1098 while (s > 0) {
1099 int ml, cij, traced, i1, j1, /*d3, d5, mm,*/ p, q;
1100 int canonical = 1; /* (i,j) closes a canonical structure */
1101 i = sector[s].i;
1102 j = sector[s].j;
1103 ml = sector[s--].ml; /* ml is a flag indicating if backtracking is to
1104 * occur in the fML- (1) or in the f-array (0) */
1105 if (ml == 2) {
1106 base_pair[++b].i = i;
1107 base_pair[b].j = j;
1108 goto repeat1;
1109 }
1110
1111 if (j < i + TURN + 1)
1112 continue; /* no more pairs in this interval */
1113
1114 repeat1:
1115
1116 /*----- begin of "repeat:" -----*/
1117 if (canonical)
1118 cij = c[indx[j] + i];
1119
1120 type = ptype[indx[j] + i];
1121 bonus = 0;
1122 if (fold_constrained) {
1123 if ((BP[i] == j) || (BP[i] == -1) || (BP[i] == -2))
1124 bonus -= BONUS;
1125
1126 if ((BP[j] == -1) || (BP[j] == -3))
1127 bonus -= BONUS;
1128 }
1129
1130 if (noLonelyPairs) {
1131 if (cij == c[indx[j] + i]) {
1132 /* (i.j) closes canonical structures, thus
1133 * (i+1.j-1) must be a pair */
1134 type_2 = ptype[indx[j - 1] + i + 1];
1135 type_2 = rtype[type_2];
1136 cij -= P->stack[type][type_2] + bonus;
1137 base_pair[++b].i = i + 1;
1138 base_pair[b].j = j - 1;
1139 i++;
1140 j--;
1141 canonical = 0;
1142 goto repeat1;
1143 }
1144 }
1145
1146 canonical = 1;
1147 no_close = (((type == 3) || (type == 4)) && no_closingGU && (bonus == 0));
1148 if (no_close) {
1149 if (cij == FORBIDDEN)
1150 continue;
1151 } else
1152 if (cij == E_Hairpin(j - i - 1, type, S1[i + 1], S1[j - 1], string + i - 1, P) + bonus) {
1153 continue;
1154 }
1155
1156 for (p = i + 1; p <= MIN2(j - 2 - TURN, i + MAXLOOP + 1); p++) {
1157 int minq;
1158 minq = j - i + p - MAXLOOP - 2;
1159 if (minq < p + 1 + TURN)
1160 minq = p + 1 + TURN;
1161
1162 for (q = j - 1; q >= minq; q--) {
1163 type_2 = ptype[indx[q] + p];
1164 if (type_2 == 0)
1165 continue;
1166
1167 type_2 = rtype[type_2];
1168 if (no_closingGU)
1169 if (no_close || (type_2 == 3) || (type_2 == 4))
1170 if ((p > i + 1) || (q < j - 1))
1171 continue;
1172
1173 /* continue unless stack */
1174
1175 energy = E_IntLoop(p - i - 1, j - q - 1, type, type_2,
1176 S1[i + 1], S1[j - 1], S1[p - 1], S1[q + 1], P);
1177 new = energy + c[indx[q] + p] + bonus;
1178 traced = (cij == new);
1179 if (traced) {
1180 base_pair[++b].i = p;
1181 base_pair[b].j = q;
1182 i = p, j = q;
1183 goto repeat1;
1184 }
1185 }
1186 }
1187
1188 /* end of repeat: --------------------------------------------------*/
1189
1190 /* (i.j) must close a multi-loop */
1191 /* tt = rtype[type]; */
1192 /* mm = bonus+P->MLclosing+P->MLintern[tt]; */
1193 /* d5 = P->dangle5[tt][S1[j-1]]; */
1194 /* d3 = P->dangle3[tt][S1[i+1]]; */
1195 i1 = i + 1;
1196 j1 = j - 1;
1197 sector[s + 1].ml = sector[s + 2].ml = 1;
1198
1199 /* if (k<=j-3-TURN) { */ /* found the decomposition */
1200 /* sector[++s].i = i1; */
1201 /* sector[s].j = k; */
1202 /* sector[++s].i = k+1; */
1203 /* sector[s].j = j1; */
1204 /* } else { */
1205 /* vrna_message_error("backtracking failed in repeat"); */
1206 /* } */
1207 /* */
1208 }
1209
1210 base_pair[0].i = b; /* save the total number of base pairs */
1211 }
1212
1213
1214 char *
snobacktrack_fold_from_pair(const char * sequence,int i,int j)1215 snobacktrack_fold_from_pair(const char *sequence,
1216 int i,
1217 int j)
1218 {
1219 char *structure;
1220
1221 sector[1].i = i;
1222 sector[1].j = j;
1223 sector[1].ml = 2;
1224 base_pair[0].i = 0;
1225 encode_seq(sequence);
1226 backtrack(sequence, 1);
1227 structure = vrna_db_from_bp_stack(base_pair, strlen(sequence));
1228 free(S);
1229 free(S1);
1230 return structure;
1231 }
1232
1233
1234 char *
alisnobacktrack_fold_from_pair(const char ** strings,int i,int j,int * cov)1235 alisnobacktrack_fold_from_pair(const char **strings,
1236 int i,
1237 int j,
1238 int *cov)
1239 {
1240 char *structure;
1241 int n_seq, s, length;
1242
1243 length = (int)strlen(strings[0]);
1244 for (s = 0; strings[s] != NULL; s++) ;
1245 n_seq = s;
1246 sector[1].i = i;
1247 sector[1].j = j;
1248 sector[1].ml = 2;
1249 base_pair[0].i = 0;
1250 /* encode_seq(sequence); */
1251 Sali = (short **)vrna_alloc(n_seq * sizeof(short *));
1252 for (s = 0; s < n_seq; s++) {
1253 if (strlen(strings[s]) != length)
1254 vrna_message_error("uneqal seqence lengths");
1255
1256 Sali[s] = aliencode_seq(strings[s]);
1257 }
1258 *cov = alibacktrack(strings, 1);
1259 structure = vrna_db_from_bp_stack(base_pair, length);
1260 free(S);
1261 free(S1);
1262 for (s = 0; s < n_seq; s++)
1263 free(Sali[s]);
1264 free(Sali);
1265 return structure;
1266 }
1267
1268
1269 /*---------------------------------------------------------------------------*/
1270
1271
1272 /*---------------------------------------------------------------------------*/
1273
1274
1275 /*--------------------------------------------------------------------------*/
1276
1277
1278 /*---------------------------------------------------------------------------*/
1279
1280
1281 PRIVATE void
encode_seq(const char * sequence)1282 encode_seq(const char *sequence)
1283 {
1284 unsigned int i, l;
1285
1286 l = strlen(sequence);
1287 S = (short *)vrna_alloc(sizeof(short) * (l + 2));
1288 S1 = (short *)vrna_alloc(sizeof(short) * (l + 2));
1289 /* S1 exists only for the special X K and I bases and energy_set!=0 */
1290 S[0] = (short)l;
1291
1292 for (i = 1; i <= l; i++) {
1293 /* make numerical encoding of sequence */
1294 S[i] = (short)encode_char(toupper(sequence[i - 1]));
1295 S1[i] = alias[S[i]]; /* for mismatches of nostandard bases */
1296 }
1297 /* for circular folding add first base at position n+1 and last base at
1298 * position 0 in S1 */
1299 S[l + 1] = S[1];
1300 S1[l + 1] = S1[1];
1301 S1[0] = S1[l];
1302 }
1303
1304
1305 PRIVATE short *
aliencode_seq(const char * sequence)1306 aliencode_seq(const char *sequence)
1307 {
1308 unsigned int i, l;
1309 short *Stemp;
1310
1311 l = strlen(sequence);
1312 Stemp = (short *)vrna_alloc(sizeof(short) * (l + 2));
1313 Stemp[0] = (short)l;
1314
1315 /* make numerical encoding of sequence */
1316 for (i = 1; i <= l; i++)
1317 Stemp[i] = (short)encode_char(toupper(sequence[i - 1]));
1318
1319 /* for circular folding add first base at position n+1 */
1320 /* Stemp[l+1] = Stemp[1]; */
1321
1322 return Stemp;
1323 }
1324
1325
1326 /*---------------------------------------------------------------------------*/
1327
1328 PUBLIC void
snoupdate_fold_params(void)1329 snoupdate_fold_params(void)
1330 {
1331 vrna_md_t md;
1332
1333 if (P)
1334 free(P);
1335
1336 set_model_details(&md);
1337 P = vrna_params(&md);
1338 make_pair_matrix();
1339 if (init_length < 0)
1340 init_length = 0;
1341 }
1342
1343
1344 /*---------------------------------------------------------------------------*/
1345
1346 PRIVATE void
make_ptypes(const short * S,const char * structure)1347 make_ptypes(const short *S,
1348 const char *structure)
1349 {
1350 int n, i, j, k, l;
1351
1352 n = S[0];
1353 for (k = 1; k < n - TURN; k++)
1354 for (l = 1; l <= 2; l++) {
1355 int type, ntype = 0, otype = 0;
1356 i = k;
1357 j = i + TURN + l;
1358 if (j > n)
1359 continue;
1360
1361 type = pair[S[i]][S[j]];
1362 while ((i >= 1) && (j <= n)) {
1363 if ((i > 1) && (j < n))
1364 ntype = pair[S[i - 1]][S[j + 1]];
1365
1366 if (noLonelyPairs && (!otype) && (!ntype))
1367 type = 0; /* i.j can only form isolated pairs */
1368
1369 ptype[indx[j] + i] = (char)type;
1370 otype = type;
1371 type = ntype;
1372 i--;
1373 j++;
1374 }
1375 }
1376
1377 if (fold_constrained && (structure != NULL))
1378 constrain_ptypes(structure, (unsigned int)n, ptype, BP, TURN, 0);
1379 }
1380