binaries/src/ViennaRNA/lib/snofold.c

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