binaries/src/ViennaRNA/lib/plex.c

   1
   2 /* Last changed Time-stamp: <2007-10-30 14:06:22 htafer> */
   3 /*
   4            compute the duplex structure of two RNA strands,
   5                 allowing only inter-strand base pairs.
   6          see cofold() for computing hybrid structures without
   7                              restriction.
   8                              Ivo Hofacker
   9                           Vienna RNA package
  10
  11 */
  12
  13
  14 /*
  15   library containing the function used in rnaplex
  16   the program rnaplex uses the following function
  17   Lduplexfold: finds high scoring segments
  18   it stores the end-position of these segments in an array
  19   and call then for each of these positions the duplexfold function
  20   which allows to make backtracking for each of the high scoring position
  21   It allows to find suboptimal partially overlapping (depends on a a parameter)
  22   duplexes between a long RNA and a shorter one.
  23   Contrarly to RNAduplex, the energy model is not in E~log(N),
  24   where N is the length of an interial loop but used an affine model,
  25   where the extension and begin parameter are fitted to the energy
  26   parameter used by RNAduplex. This allows to check for duplex between a short RNA(20nt)
  27   and a long one at the speed of 1Mnt/s. At this speed the whole genome (3Gnt) can be analyzed for one siRNA
  28   in about 50 minutes.
  29   The algorithm is based on an idea by Durbin and Eddy:when the alginment reach a value larger than a
  30   given threshold this value is stored in an array. When the alignment score goes
  31   then under this threshold, the alignemnent begin from this value, in that way the backtracking allow us
  32   to find all non-overlapping high-scoring segments.
  33   For more information check "durbin, biological sequence analysis"
  34 */
  35
  36 #include <config.h>
  37 #include <stdio.h>
  38 #include <stdlib.h>
  39 #include <math.h>
  40 #include <ctype.h>
  41 #include <string.h>
  42 #include "utils.h"
  43 #include "energy_par.h"
  44 #include "fold_vars.h"
  45 #include "fold.h"
  46 #include "pair_mat.h"
  47 #include "params.h"
  48 #include "plex.h"
  49 #include "ali_plex.h"
  50 #include "loop_energies.h"
  51 /* #################SIMD############### */
  52
  53 /*@unused@*/
  54 static char rcsid[] UNUSED = "$Id: plex.c,v 1.14 2007/06/12 12:50:16 htafer Exp $";
  55 /* int subopt_sorted=0; */
  56
  57 #define PUBLIC
  58 #define PRIVATE static
  59
  60 #define STACK_BULGE1  1   /* stacking energies for bulges of size 1 */
  61 #define NEW_NINIO     1   /* new asymetry penalty */
  62 #define ARRAY 32          /*array size*/
  63 #define UNIT 100
  64 #define MINPSCORE -2 * UNIT
  65
  66 /**
  67 *** Macro that define indices for the Single Array approach defined in FLduplexfold_XS->gain of 20% in runtime
  68 *** so that everything is done in a 1D array.
  69 *** input is idx for i, j for j and the length of the query RNA
  70 *** 1D is divided in 6 subarrays, one for each number of allowed state
  71 *** The length of each subarray is 5*L. 5 the maximal stored distance on the target sequence,
  72 *** L is the length of the query sequence
  73 **/
  74 #define LCI(i,j,l)      ((i     )*l + j)
  75 #define LINI(i,j,l)     ((i +  5)*l + j)
  76 #define LBXI(i,j,l)     ((i + 10)*l + j)
  77 #define LBYI(i,j,l)     ((i + 15)*l + j)
  78 #define LINIX(i,j,l)    ((i + 20)*l + j)
  79 #define LINIY(i,j,l)    ((i + 25)*l + j)
  80
  81 PRIVATE void  encode_seqs(const char *s1, const char *s2);
  82 PRIVATE short *encode_seq(const char *seq);
  83 PRIVATE void  update_dfold_params(void);
  84 /**
  85 *** duplexfold(_XS)/backtrack(_XS) computes duplex interaction with standard energy and considers extension_cost
  86 *** find_max(_XS)/plot_max(_XS) find suboptimals and MFE
  87 *** fduplexfold(_XS) computes duplex in a plex way
  88 **/
  89 PRIVATE duplexT duplexfold(const char *s1, const char *s2, const int extension_cost);
  90 PRIVATE char * backtrack(int i, int j, const int extension_cost);
  91 PRIVATE void   find_max(const int *position, const int *position_j, const int delta, const int threshold, const int length, const char *s1, const char *s2, const int extension_cost, const int fast,const int il_a, const int il_b, const int b_a, const int b_b);
  92 PRIVATE void   plot_max(const int max, const int max_pos, const int max_pos_j, const int alignment_length, const char *s1, const char *s2, const int extension_cost, const int fast,const int il_a, const int il_b, const int b_a, const int b_b);
  93
  94 /* PRIVATE duplexT duplexfold_XS(const char *s1, const char *s2,const int **access_s1, const int **access_s2, const int i_pos, const int j_pos, const int threshold); */
  95 PRIVATE duplexT duplexfold_XS(const char *s1, const char *s2,const int **access_s1, const int **access_s2, const int i_pos, const int j_pos, const int threshold, const int i_flag, const int j_flag);
  96 /* PRIVATE char *   backtrack_XS(int i, int j, const int** access_s1, const int** access_s2); */
  97 PRIVATE char *backtrack_XS(int i, int j, const int **access_s1,const int **access_s2,const int i_flag, const int j_flag );
  98 PRIVATE void find_max_XS(const int *position, const int *position_j,const int delta, const int threshold, const int alignment_length,
  99                          const char *s1, const char *s2, const int **access_s1, const int **access_s2, const int fast,const int il_a, const int il_b, const int b_a, const int b_b);
 100 PRIVATE void plot_max_XS(const int max, const int max_pos, const int max_pos_j, const int alignment_length, const char *s1, const char *s2, const int **access_s1, const int **access_s2, const int fast,const int il_a, const int il_b, const int b_a, const int b_b);
 101 PRIVATE duplexT fduplexfold(const char *s1, const char *s2, const int extension_cost, const int il_a, const int il_b, const int b_a, const int b_b);
 102 PRIVATE char *fbacktrack(int i, int j, const int extension_cost,const int il_a, const int il_b, const int b_a, const int b_b, int *dG);
 103 PRIVATE duplexT fduplexfold_XS(const char *s1, const char *s2, const int **access_s1, const int **access_s2, const int i_pos, const int j_pos, const int threshold,const int il_a, const int il_b, const int b_a, const int b_b);
 104 PRIVATE char * fbacktrack_XS(int i, int j, const int **access_s1, const int **access_s2, const int i_pos, const int j_pos, const int il_a, const int il_b, const int b_a, const int b_b, int *dGe, int *dGeplex, int *dGx, int *dGy);
 105
 106
 107
 108 /*@unused@*/
 109
 110 #define MAXSECTORS      500     /* dimension for a backtrack array */
 111 #define LOCALITY        0.      /* locality parameter for base-pairs */
 112
 113 #define MIN2(A, B)      ((A) < (B) ? (A) : (B))
 114 #define MAX2(A, B)      ((A) > (B) ? (A) : (B))
 115
 116 PRIVATE paramT *P = NULL;
 117
 118 /**
 119 *** energy array used in fduplexfold and fduplexfold_XS
 120 *** We do not use the 1D array here as it is not time critical
 121 *** It also makes the code more readable
 122 *** c -> stack;in -> interior loop;bx/by->bulge;inx/iny->1xn loops
 123 **/
 124
 125 PRIVATE int   **c=NULL, **in=NULL, **bx=NULL, **by=NULL, **inx=NULL, **iny=NULL;
 126
 127 /**
 128 *** S1, SS1, ... contains the encoded sequence for target and query
 129 *** n1, n2, n3, n4 contains target and query length
 130 **/
 131
 132 PRIVATE short  *S1=NULL, *SS1=NULL, *S2=NULL, *SS2=NULL;/*contains the sequences*/
 133 PRIVATE int   n1,n2;    /* sequence lengths */
 134 PRIVATE int n3, n4; /*sequence length for the duplex*/;
 135
 136
 137 /*-----------------------------------------------------------------------duplexfold_XS---------------------------------------------------------------------------*/
 138
 139 /**
 140 *** duplexfold_XS is the pendant to the duplex function as defined in duplex.c
 141 *** but takes the accessibility into account. It is similar to the MFE version of RNAup
 142 *** The only approximation made is that target 3' end - query 5' end base pair is known
 143 *** s1,s2 are the query and target sequence; access_s1, access_s2 are the accessibility
 144 *** profiles, i_pos, j_pos are the coordinates of the closing pair.
 145 **/
 146
 147
 148
 149 PRIVATE duplexT duplexfold_XS(const char *s1, const char *s2, const int **access_s1, const int **access_s2, const int i_pos, const int j_pos, const int threshold, const int i_flag, const int j_flag) {
 150   int i, j,p,q, Emin=INF, l_min=0, k_min=0;
 151   char *struc;
 152   struc=NULL;
 153   duplexT mfe;
 154   n3 = (int) strlen(s1);
 155   n4 = (int) strlen(s2);
 156   if ((!P) || (fabs(P->temperature - temperature)>1e-6)) {
 157     update_fold_params();  if(P) free(P); P = scale_parameters();
 158     make_pair_matrix();
 159   }
 160
 161   c = (int **) space(sizeof(int *) * (n3+1));
 162   for (i=0; i<=n3; i++) c[i] = (int *) space(sizeof(int) * (n4+1));
 163   for (i=0; i<=n3; i++){
 164     for(j=0;j<=n4;j++){
 165       c[i][j]=INF;
 166     }
 167   }
 168   encode_seqs(s1, s2);
 169   int type, type2, type3, E, k,l;
 170   i=n3-i_flag; j=1+j_flag;
 171   type = pair[S1[i]][S2[j]];
 172   if(!type){
 173     printf("Error during initialization of the duplex in duplexfold_XS\n");
 174     mfe.structure=NULL;
 175     mfe.energy = INF;
 176     return mfe;
 177   }
 178   c[i][j] = P->DuplexInit;
 179   /**  if (type>2) c[i][j] += P->TerminalAU;
 180   ***  c[i][j]+=P->dangle3[rtype[type]][SS1[i+1]];
 181   ***  c[i][j]+=P->dangle5[rtype[type]][SS2[j-1]];
 182   *** The three above lines are replaced by the line below
 183   **/
 184
 185
 186   c[i][j] += E_ExtLoop(rtype[type], (j_flag ? SS2[j-1] : -1) , (i_flag ? SS1[i+1] : -1),  P);
 187
 188 /*   if(j_flag ==0 && i_flag==0){ */
 189 /*     c[i][j] += E_ExtLoop(rtype[type], -1 , -1 , P); */
 190 /*   }else if(j_flag ==0 && i_flag==1){ */
 191 /*     c[i][j] += E_ExtLoop(rtype[type], -1 , SS1[i+1], P); */
 192 /*   }else if(j_flag ==1 && i_flag==0){ */
 193 /*     c[i][j] += E_ExtLoop(rtype[type], SS2[j-1] , -1, P); */
 194 /*   }else { */
 195 /*     c[i][j] += E_ExtLoop(rtype[type], SS2[j-1] , SS1[i+1], P); */
 196 /*   } */
 197   /*  Just in case we have only one bp, we initialize ... */
 198   /*  k_min, l_min and Emin */
 199   k_min=i; l_min=j;Emin=c[i][j];
 200   for (k=i; k>1 ; k--) {
 201     if(k<i) c[k+1][0]=INF;
 202     for (l=j; l<=n4-1; l++) {
 203       if(!(k==i && l==j)){
 204         c[k][l]=INF;
 205       }
 206       type2 = pair[S1[k]][S2[l]];
 207       if (!type2) continue;
 208       for (p=k+1; p<= n3 - i_flag && p<k+MAXLOOP-1; p++) {
 209         for (q = l-1; q >= 1+j_flag; q--) {
 210           if (p-k+l-q-2>MAXLOOP) break;
 211           type3=pair[S1[p]][S2[q]];
 212           if(!type3) continue;
 213           E = E_IntLoop(p-k-1, l-q-1, type2, rtype[type3],SS1[k+1], SS2[l-1], SS1[p-1], SS2[q+1],P);
 214           c[k][l] = MIN2(c[k][l], c[p][q]+E);
 215         }
 216       }
 217       E = c[k][l];
 218       E+=access_s1[i-k+1][i_pos]+access_s2[l-1][j_pos+(l-1)-1];
 219       /**if (type2>2) E += P->TerminalAU;
 220       ***if (k>1) E += P->dangle5[type2][SS1[k-1]];
 221       ***if (l<n4) E += P->dangle3[type2][SS2[l+1]];
 222       *** Replaced by the line below
 223       **/
 224       E+=E_ExtLoop(type2, (k>1) ? SS1[k-1] : -1, (l<n4) ? SS2[l+1] : -1, P);
 225
 226       if (E<Emin) {
 227         Emin=E; k_min=k; l_min=l;
 228       }
 229     }
 230   }
 231
 232   if(Emin  > threshold){
 233     mfe.energy=INF;
 234     mfe.ddG=INF;
 235     mfe.structure=NULL;
 236     for (i=0; i<=n3; i++) free(c[i]);
 237     free(c);
 238     free(S1); free(S2); free(SS1); free(SS2);
 239     return mfe;
 240   } else{
 241     struc = backtrack_XS(k_min, l_min, access_s1, access_s2, i_flag, j_flag);
 242   }
 243
 244
 245   /**
 246   *** find best dangles combination
 247   **/
 248   int dx_5, dx_3, dy_5, dy_3,dGx,dGy,bonus_x;
 249   dx_5=0; dx_3=0; dy_5=0; dy_3=0;dGx=0;dGy=0;bonus_x=0;
 250   /* x--------x */
 251   /*  |||||||| */
 252   /* x--------x */
 253    dGx = access_s1[i-k_min+1][i_pos];dx_3=0; dx_5=0;bonus_x=0;
 254    dGy = access_s2[l_min-j+1][j_pos + (l_min-1)];
 255
 256   mfe.tb=i_pos -9 - i + k_min -1 -dx_5;
 257   mfe.te=i_pos -9 -1 + dx_3;
 258   mfe.qb=j_pos -9 -1 - dy_5;
 259   mfe.qe=j_pos + l_min -3 -9 + dy_3;
 260   mfe.ddG=(double) Emin * 0.01;
 261   mfe.dG1=(double) dGx*0.01 ;
 262   mfe.dG2=(double) dGy*0.01 ;
 263
 264   mfe.energy= mfe.ddG - mfe.dG1 - mfe.dG2;
 265
 266   mfe.structure = struc;
 267   for (i=0; i<=n3; i++) free(c[i]);
 268   free(c);
 269   free(S1); free(S2); free(SS1); free(SS2);
 270   return mfe;
 271 }
 272
 273
 274
 275
 276
 277
 278 PRIVATE char *backtrack_XS(int i, int j, const int **access_s1,const int **access_s2, const int i_flag, const int j_flag) {
 279   /* backtrack structure going backwards from i, and forwards from j
 280      return structure in bracket notation with & as separator */
 281   int k, l, type, type2, E, traced, i0, j0;
 282   char *st1, *st2, *struc;
 283   st1 = (char *) space(sizeof(char)*(n3+1));
 284   st2 = (char *) space(sizeof(char)*(n4+1));
 285   i0=i;/*MAX2(i-1,1);*/j0=j;/*MIN2(j+1,n4);*/
 286   while (i<=n3-i_flag && j>=1+j_flag) {
 287     E = c[i][j]; traced=0;
 288     st1[i-1] = '(';
 289     st2[j-1] = ')';
 290     type = pair[S1[i]][S2[j]];
 291     if (!type) nrerror("backtrack failed in fold duplex bli");
 292     for (k=i+1; k<=n3 && k>i-MAXLOOP-2; k++) {
 293       for (l=j-1; l>=1; l--) {
 294         int LE;
 295         if (i-k+l-j-2>MAXLOOP) break;
 296         type2 = pair[S1[k]][S2[l]];
 297         if (!type2) continue;
 298         LE = E_IntLoop(k-i-1, j-l-1, type, rtype[type2], SS1[i+1], SS2[j-1], SS1[k-1], SS2[l+1],P);
 299         if (E == c[k][l]+LE) {
 300           traced=1;
 301           i=k; j=l;
 302           break;
 303         }
 304       }
 305       if (traced) break;
 306     }
 307     if (!traced) {
 308 #if 0
 309       /**
 310       ***if (i<n3) E -= P->dangle3[rtype[type]][SS1[i+1]];//+access_s1[1][i+1];
 311       ***if (j>1)  E -= P->dangle5[rtype[type]][SS2[j-1]];//+access_s2[1][j+1];
 312       ***if (type>2) E -= P->TerminalAU;
 313       *** Changes with line below
 314       **/
 315 #endif
 316       E -= E_ExtLoop(rtype[type], SS2[j-1] , SS1[i+1], P);
 317       break;
 318       if (E != P->DuplexInit) {
 319         nrerror("backtrack failed in fold duplex bal");
 320       } else break;
 321     }
 322   }
 323   /* if (i<n3)  i++; */
 324   /* if (j>1)   j--; */
 325   struc = (char *) space(i-i0+1+j0-j+1+2);
 326   for (k=MAX2(i0,1); k<=i; k++) if (!st1[k-1]) st1[k-1] = '.';
 327   for (k=j; k<=j0; k++) if (!st2[k-1]) st2[k-1] = '.';
 328   strcpy(struc, st1+MAX2(i0-1,0)); strcat(struc, "&");
 329   strcat(struc, st2+j-1);
 330   free(st1); free(st2);
 331   return struc;
 332 }
 333
 334
 335 /**
 336 *** fduplexfold(_XS) computes the interaction based on the plex energy model.
 337 *** Faster than duplex approach, but not standard model compliant
 338 *** We use the standard matrix (c, in, etc..., because we backtrack)
 339 **/
 340
 341 PRIVATE duplexT fduplexfold_XS(const char *s1, const char *s2, const int **access_s1, const int **access_s2, const int i_pos, const int j_pos, const int threshold,const int il_a, const int il_b, const int b_a, const int b_b) {
 342   /**
 343   *** i,j recursion index
 344   *** Emin, i_min, j_min MFE position and energy
 345   *** mfe struc duplex structure
 346   **/
 347   int i, j, Emin, i_min, j_min,l1;
 348   duplexT mfe;
 349   char *struc;
 350   /**
 351   *** bext=b_a bulge extension parameter for linear model
 352   *** iopen=il_b interior opening for linear model
 353   *** iext_s=2*il_a asymmetric extension for interior loop
 354   *** iext_ass=60+il_a symmetric extension for interior loop
 355   *** min_colonne=INF; max score of a row
 356   *** i_length;
 357   *** max_pos; position of best hit during recursion on target
 358   *** max_pos_j; position of best hit during recursion on query
 359   *** temp; temp variable for min_colonne
 360   *** min_j_colonne; position of the minimum on query in row j
 361   *** max=INF; absolute MFE
 362   *** n3,n4 length of target and query
 363   *** DJ contains the accessibility penalty for the query sequence
 364   *** maxPenalty contains the maximum penalty
 365   **/
 366   int bopen=b_b;
 367   int bext=b_a;
 368   int iopen=il_b;
 369   int iext_s=2*il_a;/* iext_s 2 nt nucleotide extension of interior loop, on i and j side */
 370   int iext_ass=50+il_a;/* iext_ass assymetric extension of interior loop, either on i or on j side. */
 371   int min_colonne=INF; /* enthaelt das maximum einer kolonne */
 372   int i_length;
 373   int max_pos;/* get position of the best hit */
 374   int max_pos_j;
 375   int temp=INF;
 376   int min_j_colonne;
 377   int max=INF;
 378   int **DJ;
 379   int maxPenalty[4];
 380   /**
 381   *** variable initialization
 382   **/
 383   n3 = (int) strlen(s1);
 384   n4 = (int) strlen(s2);
 385   if ((!P) || (fabs(P->temperature - temperature)>1e-6)) {
 386     update_fold_params();  if(P) free(P); P = scale_parameters();
 387     make_pair_matrix();
 388   }
 389   /**
 390   *** array initialization
 391   **/
 392   c  = (int**) space(sizeof(int *)  * (n3+1));
 393   in = (int**) space (sizeof(int *) * (n3+1));
 394   bx = (int**) space (sizeof(int *) * (n3+1));
 395   by = (int**) space (sizeof(int *) * (n3+1));
 396   inx= (int**) space (sizeof(int *) * (n3+1));
 397   iny= (int**) space (sizeof(int *) * (n3+1));
 398   for (i=0; i<=n3; i++){
 399     c[i]  = (int *) space(sizeof(int) * (n4+1));
 400     in[i] = (int *) space(sizeof(int) * (n4+1));
 401     bx[i] = (int *) space(sizeof(int) * (n4+1));
 402     by[i] = (int *) space(sizeof(int) * (n4+1));
 403     inx[i]= (int *) space(sizeof(int) * (n4+1));
 404     iny[i]= (int *) space(sizeof(int) * (n4+1));
 405   }
 406   for(i=0; i<n3; i++){
 407     for(j=0; j<n4; j++){
 408       in[i][j]=INF;/* no in before  1 */
 409       c[i][j] =INF; /* no bulge and no in before n2 */
 410       bx[i][j]=INF;/* no bulge before 1 */
 411       by[i][j]=INF;
 412       inx[i][j]=INF;/* no bulge before 1 */
 413       iny[i][j]=INF;
 414     }
 415   }
 416   /**
 417   *** sequence encoding
 418   **/
 419   encode_seqs(s1,s2);
 420   /**
 421   *** Compute accessibility penalty for the query only once
 422   **/
 423   maxPenalty[0]=(int) -1*P->stack[2][2]/2;
 424   maxPenalty[1]=(int) -1*P->stack[2][2];
 425   maxPenalty[2]=(int) -3*P->stack[2][2]/2;
 426   maxPenalty[3]=(int) -2*P->stack[2][2];
 427
 428
 429   DJ=(int **)   space(4*sizeof(int*));
 430   DJ[0]=(int *) space((1+n4)*sizeof(int));
 431   DJ[1]=(int *) space((1+n4)*sizeof(int));
 432   DJ[2]=(int *) space((1+n4)*sizeof(int));
 433   DJ[3]=(int *) space((1+n4)*sizeof(int));
 434
 435   j=n4-9;
 436   while(--j>9){
 437     int jdiff=j_pos+j-11;
 438     DJ[0][j] = access_s2[5][jdiff+4] - access_s2[4][jdiff+4]           ;
 439     DJ[1][j] = access_s2[5][jdiff+5] - access_s2[4][jdiff+5] + DJ[0][j];
 440     DJ[2][j] = access_s2[5][jdiff+6] - access_s2[4][jdiff+6] + DJ[1][j];
 441     DJ[3][j] = access_s2[5][jdiff+7] - access_s2[4][jdiff+7] + DJ[2][j];
 442     DJ[0][j] = MIN2(DJ[0][j],maxPenalty[0]);
 443     DJ[1][j] = MIN2(DJ[1][j],maxPenalty[1]);
 444     DJ[2][j] = MIN2(DJ[2][j],maxPenalty[2]);
 445     DJ[3][j] = MIN2(DJ[3][j],maxPenalty[3]);
 446   }
 447
 448   /**
 449   *** Start of the recursion
 450   *** first and last 10 nucleotides on target and query are dummy nucleotides
 451   *** allow to reduce number of if test
 452   **/
 453   i=11;
 454   i_length=n3-9;
 455   while(i < i_length) {
 456     int di1,di2,di3,di4;
 457     int idiff=i_pos-(n3-10-i);
 458     di1 = access_s1[5][idiff]   - access_s1[4][idiff-1];
 459     di2 = access_s1[5][idiff-1] - access_s1[4][idiff-2] + di1;
 460     di3 = access_s1[5][idiff-2] - access_s1[4][idiff-3] + di2;
 461     di4 = access_s1[5][idiff-3] - access_s1[4][idiff-4] + di3;
 462     di1=MIN2(di1,maxPenalty[0]);
 463     di2=MIN2(di2,maxPenalty[1]);
 464     di3=MIN2(di3,maxPenalty[2]);
 465     di4=MIN2(di4,maxPenalty[3]);
 466     j=n4-9;
 467     min_colonne=INF;
 468     while (10 < --j) {
 469       int dj1,dj2,dj3,dj4;
 470       int jdiff=j_pos+j-11;
 471       dj1 = DJ[0][j];
 472       dj2 = DJ[1][j];
 473       dj3 = DJ[2][j];
 474       dj4 = DJ[3][j];
 475       int type, type2;
 476       type = pair[S1[i]][S2[j]];
 477       /**
 478       *** Start duplex
 479       **/
 480       c[i][j]=type ? P->DuplexInit + access_s1[1][idiff]+access_s2[1][jdiff] : INF;
 481       /**
 482       *** update lin bx by linx liny matrix
 483       **/
 484       type2=pair[S2[j+1]][S1[i-1]];
 485       /**
 486       *** start/extend interior loop
 487       **/
 488       in[i][j]=MIN2(c[i - 1][j+1]+P->mismatchI[type2][SS2[j]][SS1[i]]+iopen+iext_s+di1+dj1,
 489                     in[i - 1][j]+iext_ass + di1);
 490
 491       /**
 492       *** start/extend nx1 target
 493       *** use same type2 as for in
 494       **/
 495       inx[i][j]=MIN2(c[i-1][j+1]+P->mismatch1nI[type2][SS2[j]][SS1[i]]+iopen+iext_s+di1+dj1,
 496                      inx[i-1][j]+iext_ass+di1);
 497       /**
 498       *** start/extend 1xn target
 499       *** use same type2 as for in
 500       **/
 501       iny[i][j]=MIN2(c[i-1][j+1]+P->mismatch1nI[type2][SS2[j]][SS1[i]]+iopen+iext_s+di1+dj1,
 502                      iny[i][j+1]+iext_ass+dj1);
 503       /**
 504       *** extend interior loop
 505       **/
 506       in[i][j]=MIN2(in[i][j],in[i][j+1]+iext_ass + dj1);
 507       in[i][j]=MIN2(in[i][j],in[i - 1][j+1]+iext_s + di1 + dj1);
 508       /**
 509       *** start/extend bulge target
 510       **/
 511       type2=pair[S2[j]][S1[i-1]];
 512       bx[i][j]=MIN2(bx[i - 1][j]+bext + di1, c[i - 1][j]+bopen+bext+(type2>2?P->TerminalAU:0) + di1);
 513       /**
 514       *** start/extend bulge query
 515       **/
 516       type2=pair[S2[j+1]][S1[i]];
 517       by[i][j]=MIN2(by[i][j+1]+bext+dj1, c[i][j+1]+bopen+bext+(type2>2?P->TerminalAU:0)+dj1);
 518       /**
 519       ***end update recursion
 520       ***######################## Start stack extension##############################
 521       **/
 522       if(!type){continue;}
 523       c[i][j]+=E_ExtLoop(type, SS1[i-1], SS2[j+1],P);
 524       /**
 525       *** stack extension
 526       **/
 527       if((type2=pair[S1[i-1]][S2[j+1]]))
 528         c[i][j]=MIN2(c[i - 1][j+1]+P->stack[rtype[type]][type2]+di1+dj1, c[i][j]);
 529       /**
 530       *** 1x0 / 0x1 stack extension
 531       **/
 532       if((type2=pair[S1[i-1]][S2[j+2]]))
 533         c[i][j]=MIN2(c[i - 1][j+2]+P->bulge[1]+P->stack[rtype[type]][type2]+di1+dj2,c[i][j]);
 534       if((type2=pair[S1[i-2]][S2[j+1]]))
 535         c[i][j]=MIN2(c[i - 2][j+1]+P->bulge[1]+P->stack[type2][rtype[type]]+di2+dj1,c[i][j]);
 536       /**
 537       *** 1x1 / 2x2 stack extension
 538       **/
 539       if((type2=pair[S1[i-2]][S2[j+2]]))
 540         c[i][j]=MIN2(c[i - 2][j+2]+P->int11[type2][rtype[type]][SS1[i-1]][SS2[j+1]]+di2+dj2, c[i][j]);
 541       if((type2 = pair[S1[i-3]][S2[j+3]]))
 542         c[i][j]=MIN2(c[i - 3][j+3]+P->int22[type2][rtype[type]][SS1[i-2]][SS1[i-1]][SS2[j+1]][SS2[j+2]]+di3+dj3,c[i][j]);
 543       /**
 544       *** 1x2 / 2x1 stack extension
 545       *** E_IntLoop(1,2,type2, rtype[type],SS1[i-1], SS2[j+2], SS1[i-1], SS2[j+1], P) corresponds to
 546       *** P->int21[rtype[type]][type2][SS2[j+2]][SS1[i-1]][SS1[i-1]]
 547       **/
 548       if((type2 = pair[S1[i-3]][S2[j+2]]))
 549         c[i][j]=MIN2(c[i - 3][j+2]+P->int21[rtype[type]][type2][SS2[j+1]][SS1[i-2]][SS1[i-1]]+di3+dj2, c[i][j]);
 550       if((type2 = pair[S1[i-2]][S2[j+3]]))
 551         c[i][j]=MIN2(c[i - 2][j+3]+P->int21[type2][rtype[type]][SS1[i-1]][SS2[j+1]][SS2[j+2]]+di2+dj3, c[i][j]);
 552
 553       /**
 554       *** 2x3 / 3x2 stack extension
 555       **/
 556       if((type2 = pair[S1[i-4]][S2[j+3]]))
 557         c[i][j]=MIN2(c[i - 4][j+3]+P->internal_loop[5]+P->ninio[2]+
 558                      P->mismatch23I[type2][SS1[i-3]][SS2[j+2]]+P->mismatch23I[rtype[type]][SS2[j+1]][SS1[i-1]]+di4+dj3, c[i][j]);
 559       if((type2 = pair[S1[i-3]][S2[j+4]]))
 560         c[i][j]=MIN2(c[i - 3][j+4]+P->internal_loop[5]+P->ninio[2]+
 561                      P->mismatch23I[type2][SS1[i-2]][SS2[j+3]]+P->mismatch23I[rtype[type]][SS2[j+1]][SS1[i-1]]+di3+dj4, c[i][j]);
 562
 563       /**
 564       *** So now we have to handle 1x3, 3x1, 3x3, and mxn m,n > 3
 565       **/
 566       /**
 567       *** 3x3 or more
 568       **/
 569       c[i][j]=MIN2(in[i - 3][j+3]+P->mismatchI[rtype[type]][SS1[i-1]][SS2[j+1]]+2*iext_s+di3+dj3, c[i][j]);
 570       /**
 571       *** 2xn or more
 572       **/
 573       c[i][j]=MIN2(in[i - 4][j+2]+P->mismatchI[rtype[type]][SS1[i-1]][SS2[j+1]]+iext_s+2*iext_ass+di4+dj2, c[i][j]);
 574       /**
 575       *** nx2 or more
 576       **/
 577       c[i][j]=MIN2(in[i - 2][j+4]+P->mismatchI[rtype[type]][SS1[i-1]][SS2[j+1]]+iext_s+2*iext_ass+di2+dj4, c[i][j]);
 578       /**
 579       *** nx1 n>2
 580       **/
 581       c[i][j]=MIN2(inx[i - 3][j+1]+P->mismatch1nI[rtype[type]][SS1[i-1]][SS2[j+1]]+iext_ass+iext_ass+di3+dj1, c[i][j]);
 582       /**
 583       *** 1xn n>2
 584       **/
 585       c[i][j]=MIN2(iny[i - 1][j+3]+P->mismatch1nI[rtype[type]][SS1[i-1]][SS2[j+1]]+iext_ass+iext_ass+dj3+di1, c[i][j]);
 586       /**
 587       *** nx0 n>1
 588       **/
 589       int bAU;
 590       bAU=(type>2?P->TerminalAU:0);
 591       c[i][j]=MIN2(bx[i - 2][j+1]+di2+dj1+bext+bAU, c[i][j]);
 592       /**
 593       *** 0xn n>1
 594       **/
 595       c[i][j]=MIN2(by[i - 1][j+2]+di1+dj2+bext+bAU, c[i][j]);
 596       temp=min_colonne;
 597       min_colonne=MIN2(c[i][j]+ E_ExtLoop(rtype[type], SS2[j-1], SS1[i+1], P),min_colonne);
 598       if(temp>min_colonne){
 599         min_j_colonne=j;
 600       }
 601       /* ---------------------------------------------------------------------end update */
 602     }
 603     if(max>=min_colonne){
 604       max=min_colonne;
 605       max_pos=i;
 606       max_pos_j=min_j_colonne;
 607     }
 608     i++;
 609   }
 610   Emin=max;
 611   if(Emin>threshold){
 612     free(S1); free(S2); free(SS1); free(SS2);
 613     for (i=0; i<=n3; i++) {
 614       free(c[i]);
 615       free(in[i]);
 616       free(bx[i]);
 617       free(by[i]);
 618       free(inx[i]);
 619       free(iny[i]);
 620     }
 621     for (i=0; i<=3; i++) {
 622       free(DJ[i]);
 623     }
 624     free(c);free(in);free(bx);free(by);free(inx);free(iny);free(DJ);
 625     mfe.energy=0;
 626     mfe.structure=NULL;
 627     return mfe;
 628   }
 629   i_min=max_pos;
 630   j_min=max_pos_j;
 631   int dGe, dGeplex, dGx, dGy;
 632   dGe=dGeplex=dGx=dGy=0;
 633   /* printf("MAX fduplexfold_XS %d\n",Emin); */
 634   struc = fbacktrack_XS(i_min, j_min, access_s1, access_s2, i_pos, j_pos,il_a, il_b,b_a,b_b,&dGe, &dGeplex, &dGx, &dGy);
 635
 636   l1 = strchr(struc, '&')-struc;
 637   int size;
 638   size=strlen(struc)-1;
 639   int lengthx; int endx; int lengthy; int endy;
 640   lengthx=l1;
 641   lengthx-=(struc[0]=='.'?1:0);
 642   lengthx-=(struc[l1-1]=='.'?1:0);
 643   endx=(i_pos-(n3-i_min));
 644   lengthy=size-l1;
 645   lengthy-=(struc[size]=='.'?1:0);
 646   lengthy-=(struc[l1+1]=='.'?1:0);
 647   endy=j_pos+j_min+lengthy -22;
 648   if (i_min<n3-10) i_min++;
 649   if (j_min>11 ) j_min--;
 650   mfe.i = i_min;
 651   mfe.j = j_min;
 652   mfe.ddG = (double) Emin*0.01;
 653   mfe.structure = struc;
 654   mfe.energy_backtrack = (double) dGe * 0.01;
 655   mfe.energy = (double) dGeplex * 0.01;
 656   mfe.opening_backtrack_x = (double) dGx * 0.01;
 657   mfe.opening_backtrack_y = (double) dGy * 0.01;
 658   mfe.dG1=(double) access_s1[lengthx][endx+10] * 0.01;
 659   mfe.dG2=(double) access_s2[lengthy][endy+10] * 0.01;
 660   free(S1); free(S2); free(SS1); free(SS2);
 661   for (i=0; i<=n3; i++) {
 662     free(c[i]);
 663     free(in[i]);
 664     free(bx[i]);
 665     free(by[i]);
 666     free(inx[i]);
 667     free(iny[i]);
 668   }
 669   for (i=0; i<=3; i++) {
 670     free(DJ[i]);
 671   }
 672   free(DJ);
 673   free(c);free(in);free(bx);free(by);free(iny);free(inx);
 674   return mfe;
 675 }
 676
 677 PRIVATE char *fbacktrack_XS(int i, int j, const int** access_s1, const int** access_s2, const int i_pos, const int j_pos,const int il_a, const int il_b, const int b_a, const int b_b, int *dG, int *dGplex, int *dGx, int *dGy) {
 678   /* backtrack structure going backwards from i, and forwards from j
 679      return structure in bracket notation with & as separator */
 680   int k, l, type, type2, E, traced, i0, j0;
 681   char *st1, *st2, *struc;
 682   int bopen=b_b;
 683   int bext=b_a;
 684   int iopen=il_b;
 685   int iext_s=2*il_a;/* iext_s 2 nt nucleotide extension of interior loop, on i and j side */
 686   int iext_ass=50+il_a;/* iext_ass assymetric extension of interior loop, either on i or on j side. */
 687   st1 = (char *) space(sizeof(char)*(n3+1));
 688   st2 = (char *) space(sizeof(char)*(n4+1));
 689   i0=MIN2(i+1,n3-10); j0=MAX2(j-1,11);
 690   int state;
 691   state=1; /* we start backtracking from a a pair , i.e. c-matrix */
 692   /* state 1 -> base pair, c
 693      state 2 -> interior loop, in
 694      state 3 -> bx loop, bx
 695      state 4 -> by loop, by
 696   */
 697   traced=1;
 698   k=i; l=j; /* stores the i,j information for subsequence usage see * */
 699   int idiff,jdiff;
 700   /**
 701   *** (type>2?P->TerminalAU:0)+P->dangle3[rtype[type]][SS1[i+1]]+P->dangle5[rtype[type]][SS2[j-1]];
 702   **/
 703
 704   int maxPenalty[4];
 705   if ((!P) || (fabs(P->temperature - temperature)>1e-6)){
 706     update_dfold_params();
 707     if(P) free(P); P = scale_parameters();
 708     make_pair_matrix();
 709   }
 710   maxPenalty[0]=(int) -1*P->stack[2][2]/2;
 711   maxPenalty[1]=(int) -1*P->stack[2][2];
 712   maxPenalty[2]=(int) -3*P->stack[2][2]/2;
 713   maxPenalty[3]=(int) -2*P->stack[2][2];
 714
 715   type = pair[S1[i]][S2[j]];
 716   *dG+= E_ExtLoop(rtype[type], SS2[j-1] , SS1[i+1] , P);
 717   *dGplex=*dG;
 718
 719   while (i>10 && j<=n4-9 && traced) {
 720     int di1,di2,di3,di4;
 721     idiff=i_pos-(n3-10-i);
 722     di1 = access_s1[5][idiff]   - access_s1[4][idiff-1];
 723     di2 = access_s1[5][idiff-1] - access_s1[4][idiff-2] + di1;
 724     di3 = access_s1[5][idiff-2] - access_s1[4][idiff-3] + di2;
 725     di4 = access_s1[5][idiff-3] - access_s1[4][idiff-4] + di3;
 726     di1=MIN2(di1,maxPenalty[0]);
 727     di2=MIN2(di2,maxPenalty[1]);
 728     di3=MIN2(di3,maxPenalty[2]);
 729     di4=MIN2(di4,maxPenalty[3]);
 730     int dj1,dj2,dj3,dj4;
 731     jdiff=j_pos+j-11;
 732     dj1 = access_s2[5][jdiff+4] - access_s2[4][jdiff+4];
 733     dj2 = access_s2[5][jdiff+5] - access_s2[4][jdiff+5] + dj1;
 734     dj3 = access_s2[5][jdiff+6] - access_s2[4][jdiff+6] + dj2;
 735     dj4 = access_s2[5][jdiff+7] - access_s2[4][jdiff+7] + dj3;
 736     dj1=MIN2(dj1,maxPenalty[0]);
 737     dj2=MIN2(dj2,maxPenalty[1]);
 738     dj3=MIN2(dj3,maxPenalty[2]);
 739     dj4=MIN2(dj4,maxPenalty[3]);
 740     traced=0;
 741     switch(state){
 742     case 1:
 743       type = pair[S1[i]][S2[j]];
 744       int bAU;
 745       bAU=(type>2?P->TerminalAU:0);
 746       if(!type) nrerror("backtrack failed in fold duplex");
 747       type2=pair[S1[i-1]][S2[j+1]];
 748       if(type2 && c[i][j]== (c[i - 1][j+1]+P->stack[rtype[type]][type2]+di1+dj1)){
 749         k=i-1;
 750         l=j+1;
 751         (*dG)+=E_IntLoop(i-k-1, l-j-1, type2, rtype[type],SS1[k+1], SS2[l-1], SS1[i-1], SS2[j+1],P);
 752         *dGplex+=E_IntLoop(i-k-1, l-j-1, type2, rtype[type],SS1[k+1], SS2[l-1], SS1[i-1], SS2[j+1],P);
 753         *dGx+=di1;
 754         *dGy+=dj1;
 755         st1[i-1] = '(';
 756         st2[j-1] = ')';
 757         i=k;
 758         j=l;
 759         state=1;
 760         traced=1;
 761         break;
 762       }
 763       type2=pair[S1[i-1]][S2[j+2]];
 764       if(type2 && c[i][j]==(c[i - 1][j+2]+P->bulge[1]+P->stack[rtype[type]][type2]+di1+dj2)){
 765         k=i-1;
 766         l=j+2;
 767         *dG+=E_IntLoop(i-k-1, l-j-1, type2, rtype[type],SS1[k+1], SS2[l-1], SS1[i-1], SS2[j+1],P);
 768         *dGplex+=E_IntLoop(i-k-1, l-j-1, type2, rtype[type],SS1[k+1], SS2[l-1], SS1[i-1], SS2[j+1],P);
 769         *dGx+=di1;
 770         *dGy+=dj2;
 771         st1[i-1] = '(';
 772         st2[j-1] = ')';
 773         i=k;
 774         j=l;
 775         state=1;
 776         traced=1;
 777         break;
 778       }
 779       type2=pair[S1[i-2]][S2[j+1]];
 780       if(type2 && c[i][j]==(c[i - 2][j+1]+P->bulge[1]+P->stack[type2][rtype[type]]+di2+dj1)){
 781         k=i-2;
 782         l=j+1;
 783         *dG+=E_IntLoop(i-k-1, l-j-1, type2, rtype[type],SS1[k+1], SS2[l-1], SS1[i-1], SS2[j+1],P);
 784         *dGplex+=E_IntLoop(i-k-1, l-j-1, type2, rtype[type],SS1[k+1], SS2[l-1], SS1[i-1], SS2[j+1],P);
 785         *dGx+=di2;
 786         *dGy+=dj1;
 787         st1[i-1] = '(';
 788         st2[j-1] = ')';
 789         i=k;
 790         j=l;
 791         state=1;
 792         traced=1;
 793         break;
 794       }
 795       type2=pair[S1[i-2]][S2[j+2]];
 796       if(type2 && c[i][j]==(c[i - 2][j+2]+P->int11[type2][rtype[type]][SS1[i-1]][SS2[j+1]]+di2+dj2)){
 797         k=i-2;
 798         l=j+2;
 799         *dG+=E_IntLoop(i-k-1, l-j-1, type2, rtype[type],SS1[k+1], SS2[l-1], SS1[i-1], SS2[j+1],P);
 800         *dGplex+=E_IntLoop(i-k-1, l-j-1, type2, rtype[type],SS1[k+1], SS2[l-1], SS1[i-1], SS2[j+1],P);
 801         *dGx+=di2;
 802         *dGy+=dj2;
 803         st1[i-1] = '(';
 804         st2[j-1] = ')';
 805         i=k;
 806         j=l;
 807         state=1;
 808         traced=1;
 809         break;
 810       }
 811       type2 = pair[S1[i-3]][S2[j+3]];
 812       if(type2 && c[i][j]==(c[i - 3][j+3]+P->int22[type2][rtype[type]][SS1[i-2]][SS1[i-1]][SS2[j+1]][SS2[j+2]]+di3+dj3)){
 813         k=i-3;
 814         l=j+3;
 815         *dG+=E_IntLoop(i-k-1, l-j-1, type2, rtype[type],SS1[k+1], SS2[l-1], SS1[i-1], SS2[j+1],P);
 816         *dGplex+=E_IntLoop(i-k-1, l-j-1, type2, rtype[type],SS1[k+1], SS2[l-1], SS1[i-1], SS2[j+1],P);
 817         *dGx+=di3;
 818         *dGy+=dj3;
 819         st1[i-1] = '(';
 820         st2[j-1] = ')';
 821         i=k;
 822         j=l;
 823         state=1;
 824         traced=1;
 825         break;
 826       }
 827       type2 = pair[S1[i-3]][S2[j+2]];
 828       if(type2 && c[i][j]==(c[i - 3][j+2]+P->int21[rtype[type]][type2][SS2[j+1]][SS1[i-2]][SS1[i-1]]+di3+dj2)){
 829         k=i-3;
 830         l=j+2;
 831         *dG+=E_IntLoop(i-k-1, l-j-1, type2, rtype[type],SS1[k+1], SS2[l-1], SS1[i-1], SS2[j+1],P);
 832         *dGplex+=E_IntLoop(i-k-1, l-j-1, type2, rtype[type],SS1[k+1], SS2[l-1], SS1[i-1], SS2[j+1],P);
 833         *dGx+=di3;
 834         *dGy+=dj2;
 835         st1[i-1] = '(';
 836         st2[j-1] = ')';
 837         i=k;
 838         j=l;
 839         state=1;
 840         traced=1;
 841         break;
 842       }
 843       type2 = pair[S1[i-2]][S2[j+3]];
 844       if(type2 && c[i][j]==(c[i - 2][j+3]+P->int21[type2][rtype[type]][SS1[i-1]][SS2[j+1]][SS2[j+2]]+di2+dj3)){
 845         k=i-2;
 846         l=j+3;
 847         *dG+=E_IntLoop(i-k-1, l-j-1, type2, rtype[type],SS1[k+1], SS2[l-1], SS1[i-1], SS2[j+1],P);
 848         *dGplex+=E_IntLoop(i-k-1, l-j-1, type2, rtype[type],SS1[k+1], SS2[l-1], SS1[i-1], SS2[j+1],P);
 849         *dGx+=di2;
 850         *dGy+=dj3;
 851         st1[i-1] = '(';
 852         st2[j-1] = ')';
 853         i=k;
 854         j=l;
 855         state=1;
 856         traced=1;
 857         break;
 858       }
 859       type2 = pair[S1[i-4]][S2[j+3]];
 860       if(type2 && c[i][j]==(c[i - 4][j+3]+P->internal_loop[5]+P->ninio[2]+
 861                             P->mismatch23I[type2][SS1[i-3]][SS2[j+2]]+P->mismatch23I[rtype[type]][SS2[j+1]][SS1[i-1]]+di4+dj3)){
 862         k=i-4;
 863         l=j+3;
 864         *dG+=E_IntLoop(i-k-1, l-j-1, type2, rtype[type],SS1[k+1], SS2[l-1], SS1[i-1], SS2[j+1],P);
 865         *dGplex+=E_IntLoop(i-k-1, l-j-1, type2, rtype[type],SS1[k+1], SS2[l-1], SS1[i-1], SS2[j+1],P);
 866         *dGx+=di2;
 867         *dGy+=dj3;
 868         st1[i-1] = '(';
 869         st2[j-1] = ')';
 870         i=k;
 871         j=l;
 872         state=1;
 873         traced=1;
 874         break;
 875       }
 876       type2 = pair[S1[i-3]][S2[j+4]];
 877       if(type2 && c[i][j]==(c[i - 3][j+4]+P->internal_loop[5]+P->ninio[2]+
 878                             P->mismatch23I[type2][SS1[i-2]][SS2[j+3]]+P->mismatch23I[rtype[type]][SS2[j+1]][SS1[i-1]]+di3+dj4)){
 879         k=i-3;
 880         l=j+4;
 881         *dG+=E_IntLoop(i-k-1, l-j-1, type2, rtype[type],SS1[k+1], SS2[l-1], SS1[i-1], SS2[j+1],P);
 882         *dGplex+=E_IntLoop(i-k-1, l-j-1, type2, rtype[type],SS1[k+1], SS2[l-1], SS1[i-1], SS2[j+1],P);
 883         *dGx+=di2;
 884         *dGy+=dj3;
 885         st1[i-1] = '(';
 886         st2[j-1] = ')';
 887         i=k;
 888         j=l;
 889         state=1;
 890         traced=1;
 891         break;
 892       }
 893       if(c[i][j]==(in[i - 3][j+3]+P->mismatchI[rtype[type]][SS1[i-1]][SS2[j+1]]+di3+dj3+2*iext_s)){
 894         k=i;
 895         l=j;
 896         *dGplex+=P->mismatchI[rtype[type]][SS1[i-1]][SS2[j+1]]+2*iext_s;
 897         *dGx+=di3;
 898         *dGy+=dj3;
 899         st1[i-1] = '(';
 900         st2[j-1] = ')';
 901         i=i-3;
 902         j=j+3;
 903         state=2;
 904         traced=1;
 905         break;
 906       }
 907       if(c[i][j]==(in[i - 4][j+2]+P->mismatchI[rtype[type]][SS1[i-1]][SS2[j+1]]+di4+dj2+iext_s+2*iext_ass)){
 908         k=i;
 909         l=j;
 910         *dGplex+=P->mismatchI[rtype[type]][SS1[i-1]][SS2[j+1]]+iext_s+2*iext_ass;
 911         *dGx+=di4;
 912         *dGy+=dj2;
 913         st1[i-1] = '(';
 914         st2[j-1] = ')';
 915         i=i-4;
 916         j=j+2;
 917         state=2;
 918         traced=1;
 919         break;
 920       }
 921       if(c[i][j]==(in[i - 2][j+4]+P->mismatchI[rtype[type]][SS1[i-1]][SS2[j+1]]+di2+dj4+iext_s+2*iext_ass)){
 922         k=i;
 923         l=j;
 924         *dGplex+=P->mismatchI[rtype[type]][SS1[i-1]][SS2[j+1]]+iext_s+2*iext_ass;
 925         *dGx+=di2;
 926         *dGy+=dj4;
 927         st1[i-1] = '(';
 928         st2[j-1] = ')';
 929         i=i-2;
 930         j=j+4;
 931         state=2;
 932         traced=1;
 933         break;
 934       }
 935       if(c[i][j]==(inx[i - 3][j+1]+P->mismatch1nI[rtype[type]][SS1[i-1]][SS2[j+1]]+iext_ass+iext_ass+di3+dj1)){
 936         k=i;
 937         l=j;
 938         *dGplex+=P->mismatch1nI[rtype[type]][SS1[i-1]][SS2[j+1]]+iext_ass+iext_ass+di3+dj1;
 939         *dGx+=di3;
 940         *dGy+=dj1;
 941         st1[i-1] = '(';
 942         st2[j-1] = ')';
 943         i=i-3;
 944         j=j+1;
 945         state=5;
 946         traced=1;
 947         break;
 948       }
 949       if(c[i][j]==(iny[i - 1][j+3]+P->mismatch1nI[rtype[type]][SS1[i-1]][SS2[j+1]]+iext_ass+iext_ass+di1+dj3)){
 950         k=i;
 951         l=j;
 952         *dGplex+=P->mismatch1nI[rtype[type]][SS1[i-1]][SS2[j+1]]+iext_ass+iext_ass+di1+dj3;
 953         *dGx+=di1;
 954         *dGy+=dj3;
 955         st1[i-1] = '(';
 956         st2[j-1] = ')';
 957         i=i-1;
 958         j=j+3;
 959         state=6;
 960         traced=1;
 961         break;
 962       }
 963       if(c[i][j]==(bx[i - 2][j+1]+di2+dj1+bext+bAU)){
 964         k=i;
 965         l=j;
 966         st1[i-1] = '(';
 967         st2[j-1] = ')';
 968         *dGplex+=bext+bAU;
 969         *dGx+=di2;
 970         *dGy+=dj1;
 971         i=i-2;
 972         j=j+1;
 973         state=3;
 974         traced=1;
 975         break;
 976       }
 977       if(c[i][j]==(by[i - 1][j+2]+di1+dj2+bext+bAU)){
 978         k=i;
 979         l=j;
 980         *dGplex+=bext+bAU;
 981         *dGx+=di1;
 982         *dGy+=dj2;
 983         st1[i-1] = '(';
 984         st2[j-1] = ')';
 985         i=i-1;
 986         j=j+2;
 987         state=4;
 988         traced=1;
 989         break;
 990       }
 991       break;
 992     case 2:
 993       if(in[i][j]==(in[i - 1][j+1]+iext_s + di1 + dj1)){
 994         i--;
 995         j++;
 996         *dGplex+=iext_s;
 997         *dGx+=di1;
 998         *dGy+=dj1;
 999         state=2;
1000         traced=1;
1001         break;
1002       }
1003       if(in[i][j]==(in[i - 1][j]+iext_ass + di1)){
1004         i=i-1;
1005         *dGplex+=iext_ass;
1006         *dGx+=di1;
1007         state=2;
1008         traced=1;
1009         break;
1010       }
1011       if(in[i][j]==(in[i][j+1]+iext_ass + dj1)){
1012         j++;
1013         state=2;
1014         *dGy+=dj1;
1015         *dGplex+=iext_ass;
1016         traced=1;
1017         break;
1018       }
1019       type2=pair[SS2[j+1]][SS1[i-1]];
1020       if(type2 && in[i][j]==(c[i - 1][j+1]+P->mismatchI[type2][SS2[j]][SS1[i]]+iopen+iext_s + di1 +dj1)){
1021         *dGplex+=P->mismatchI[type2][SS2[j]][SS1[i]]+iopen+iext_s;
1022         int temp; temp=k; k=i-1; i=temp;
1023         temp=l; l=j+1; j=temp;
1024         type=pair[S1[i]][S2[j]];
1025         *dG+=E_IntLoop(i-k-1, l-j-1, type2, rtype[type],SS1[k+1], SS2[l-1], SS1[i-1], SS2[j+1],P);
1026         *dGx+=di1;
1027         *dGy+=dj1;
1028         i=k;
1029         j=l;
1030         state=1;
1031         traced=1;
1032         break;
1033       }
1034     case 3:
1035       if(bx[i][j]==(bx[i - 1][j]+bext+di1)){
1036         i--;
1037         *dGplex+=bext;
1038         *dGx+=di1;
1039         state=3;
1040         traced=1;
1041         break;
1042       }
1043       type2=pair[S2[j]][S1[i-1]];
1044       if(type2 && bx[i][j]==(c[i - 1][j]+bopen+bext+(type2>2?P->TerminalAU:0)+di1)){
1045         int temp; temp=k; k=i-1; i=temp;
1046         temp=l; l=j; j=temp;
1047         type=pair[S1[i]][S2[j]];
1048         *dG+=E_IntLoop(i-k-1, l-j-1, type2, rtype[type],SS1[k+1], SS2[l-1], SS1[i-1], SS2[j+1],P);
1049         *dGplex+=bopen+bext+(type2>2?P->TerminalAU:0);
1050         *dGx+=di1;
1051         i=k;
1052         j=l;
1053         state=1;
1054         traced=1;
1055         break;
1056       }
1057     case 4:
1058       if(by[i][j]==(by[i][j+1] + bext +dj1)){
1059         j++;
1060         *dGplex+=bext;
1061         state=4;
1062         traced=1;
1063         break;
1064       }
1065       type2=pair[S2[j+1]][S1[i]];
1066       if(type2 && by[i][j]==(c[i][j+1]+bopen+bext+(type2>2?P->TerminalAU:0) + dj1)){
1067         int temp; temp=k; k=i; i=temp;
1068         temp=l; l=j+1; j=temp;
1069         type=pair[S1[i]][S2[j]];
1070         *dG+=E_IntLoop(i-k-1, l-j-1, type2, rtype[type],SS1[k+1], SS2[l-1], SS1[i-1], SS2[j+1],P);
1071         *dGplex+=bopen+bext+(type2>2?P->TerminalAU:0);
1072         *dGy+=dj1;
1073         i=k;
1074         j=l;
1075         state=1;
1076         traced=1;
1077         break;
1078       }
1079     case 5:
1080       if(inx[i][j]==(inx[i-1][j]+iext_ass+di1)) {
1081         i--;
1082         *dGplex+=iext_ass;
1083         *dGx+=di1;
1084         state=5;
1085         traced=1;
1086         break;
1087       }
1088       type2=pair[S2[j+1]][S1[i-1]];
1089       if(type2 && inx[i][j]==(c[i-1][j+1]+P->mismatch1nI[type2][SS2[j]][SS1[i]]+iopen+iext_s+di1+dj1)){
1090         *dGplex+=P->mismatch1nI[type2][SS2[j]][SS1[i]]+iopen+iext_s;
1091         int temp; temp=k; k=i-1; i=temp;
1092         temp=l; l=j+1; j=temp;
1093         type=pair[S1[i]][S2[j]];
1094         *dG+=E_IntLoop(i-k-1, l-j-1, type2, rtype[type],SS1[k+1], SS2[l-1], SS1[i-1], SS2[j+1],P);
1095         *dGx+=di1;
1096         *dGy+=dj1;
1097         i=k;
1098         j=l;
1099         state=1;
1100         traced=1;
1101         break;
1102       }
1103     case 6:
1104       if(iny[i][j]==(iny[i][j+1]+iext_ass+dj1)) {
1105         j++;
1106         *dGplex+=iext_ass;
1107         *dGx+=dj1;
1108         state=6;
1109         traced=1;
1110         break;
1111       }
1112       type2=pair[S2[j+1]][S1[i-1]];
1113       if(type2 && iny[i][j]==(c[i-1][j+1]+P->mismatch1nI[type2][SS2[j]][SS1[i]]+iopen+iext_s+di1+dj1)){
1114         *dGplex+=P->mismatch1nI[type2][SS2[j]][SS1[i]]+iopen+iext_s;
1115         int temp; temp=k; k=i-1; i=temp;
1116         temp=l; l=j+1; j=temp;
1117         type=pair[S1[i]][S2[j]];
1118         *dG+=E_IntLoop(i-k-1, l-j-1, type2, rtype[type],SS1[k+1], SS2[l-1], SS1[i-1], SS2[j+1],P);
1119         *dGx+=di1;
1120         *dGy+=dj1;
1121         i=k;
1122         j=l;
1123         state=1;
1124         traced=1;
1125         break;
1126       }
1127     }
1128   }
1129   if (!traced) {
1130     idiff=i_pos-(n3-10-i);
1131     jdiff=j_pos+j-11;
1132     E=c[i][j];
1133     /**
1134     *** if (i>1) {E -= P->dangle5[type][SS1[i-1]]; *dG+=P->dangle5[type][SS1[i-1]];*dGplex+=P->dangle5[type][SS1[i-1]];}
1135     *** if (j<n4){E -= P->dangle3[type][SS2[j+1]]; *dG+=P->dangle3[type][SS2[j+1]];*dGplex+=P->dangle3[type][SS2[j+1]];}
1136     *** if (type>2) {E -= P->TerminalAU; *dG+=P->TerminalAU;*dGplex+=P->TerminalAU;}
1137     **/
1138     int correction;
1139     correction = E_ExtLoop(type, (i>1) ? SS1[i-1] : -1, (j<n4) ? SS2[j+1] : -1, P);
1140     *dG+=correction;
1141     *dGplex+=correction;
1142     E-=correction;
1143
1144     if (E != P->DuplexInit+access_s1[1][idiff]+access_s2[1][jdiff]) {
1145       nrerror("backtrack failed in second fold duplex");
1146     }
1147     else{
1148       *dG+=P->DuplexInit;
1149       *dGplex+=P->DuplexInit;
1150       *dGx+=access_s1[1][idiff];
1151       *dGy+=access_s2[1][jdiff];
1152       st1[i-1]='(';
1153       st2[j-1]=')';
1154     }
1155   }
1156   if (i>11)  i--;
1157   if (j<n4-10) j++;
1158   struc = (char *) space(i0-i+1+j-j0+1+2);
1159   for (k=MAX2(i,1); k<=i0; k++) if (!st1[k-1]) st1[k-1] = '.';
1160   for (k=j0; k<=j; k++) if (!st2[k-1]) st2[k-1] = '.';
1161   strcpy(struc, st1+MAX2(i-1,0));
1162   strcat(struc, "&");
1163   strcat(struc, st2+j0-1);
1164   /* printf("%s %3d,%-3d : %3d,%-3d\n", struc, i,i0,j0,j);  */
1165   free(st1); free(st2);
1166   return struc;
1167 }
1168
1169
1170 duplexT ** Lduplexfold_XS(const char *s1, const char *s2, const int **access_s1, const int **access_s2, const int threshold, const int alignment_length, const int delta, const int fast, const int il_a, const int il_b, const int b_a, const int b_b)
1171 {
1172   /**
1173   *** See variable definition in fduplexfold_XS
1174   **/
1175   int i, j;
1176   int bopen=b_b;
1177   int bext=b_a;
1178   int iopen=il_b;
1179   int iext_s=2*il_a;
1180   int iext_ass=50+il_a;
1181   int min_colonne=INF;
1182   int i_length;
1183   int max_pos;
1184   int max_pos_j;
1185   int min_j_colonne;
1186   int max=INF;
1187   int *position;
1188   int *position_j;
1189   int maxPenalty[4];
1190   int **DJ;
1191   /**
1192   *** 1D array corresponding to the standard 2d recursion matrix
1193   *** Makes the computation 20% faster
1194   **/
1195   int *SA;
1196   /**
1197   *** variable initialization
1198   **/
1199   n1 = (int) strlen(s1);
1200   n2 = (int) strlen(s2);
1201   /**
1202   *** Sequence encoding
1203   **/
1204
1205   if ((!P) || (fabs(P->temperature - temperature)>1e-6)) {
1206     update_dfold_params();  if(P) free(P); P = scale_parameters();
1207     make_pair_matrix();
1208   }
1209   encode_seqs(s1,s2);
1210   /**
1211   *** Position of the high score on the target and query sequence
1212   **/
1213   position = (int *) space ((delta+n1+3+delta) * sizeof(int));
1214   position_j= (int *) space ((delta+n1+3+delta) * sizeof(int));
1215   /**
1216   *** extension penalty, computed only once, further reduce the computation time
1217   **/
1218   maxPenalty[0]=(int) -1*P->stack[2][2]/2;
1219   maxPenalty[1]=(int) -1*P->stack[2][2];
1220   maxPenalty[2]=(int) -3*P->stack[2][2]/2;
1221   maxPenalty[3]=(int) -2*P->stack[2][2];
1222
1223   DJ=(int **) space(4*sizeof(int*));
1224   DJ[0]=(int *) space(n2*sizeof(int));
1225   DJ[1]=(int *) space(n2*sizeof(int));
1226   DJ[2]=(int *) space(n2*sizeof(int));
1227   DJ[3]=(int *) space(n2*sizeof(int));
1228   j=n2-9;
1229   while(--j>10){
1230     DJ[0][j] = access_s2[5][j+4] - access_s2[4][j+4]           ;
1231     DJ[1][j] = access_s2[5][j+5] - access_s2[4][j+5] + DJ[0][j];
1232     DJ[2][j] = access_s2[5][j+6] - access_s2[4][j+6] + DJ[1][j];
1233     DJ[3][j] = access_s2[5][j+7] - access_s2[4][j+7] + DJ[2][j];
1234     DJ[0][j] = MIN2(DJ[0][j],maxPenalty[0]);
1235     DJ[1][j] = MIN2(DJ[1][j],maxPenalty[1]);
1236     DJ[2][j] = MIN2(DJ[2][j],maxPenalty[2]);
1237     DJ[3][j] = MIN2(DJ[3][j],maxPenalty[3]);
1238   }
1239   /**
1240   *** instead of having 4 2-dim arrays we use a unique 1-dim array
1241   *** The mapping 2d -> 1D is done based ont the macro
1242   *** LCI(i,j,l)      ((i     )*l + j)
1243   *** LINI(i,j,l)     ((i +  5)*l + j)
1244   *** LBXI(i,j,l)     ((i + 10)*l + j)
1245   *** LBYI(i,j,l)     ((i + 15)*l + j)
1246   *** LINIX(i,j,l)    ((i + 20)*l + j)
1247   *** LINIY(i,j,l)    ((i + 25)*l + j)
1248   ***
1249   *** SA has a length of 5 (number of columns we look back) *
1250   ***                  * 6 (number of structures we look at) *
1251   ***                  * length of the sequence
1252   **/
1253
1254   SA=(int *) space(sizeof(int)*5*6*(n2+5));
1255   for(j=n2+4;j>=0;j--) {
1256     SA[(j*30)   ]=SA[(j*30)+1   ]=SA[(j*30)+2   ]=SA[(j*30)+3   ]=SA[(j*30)+4   ]=INF;
1257     SA[(j*30)+5 ]=SA[(j*30)+1+5 ]=SA[(j*30)+2+5 ]=SA[(j*30)+3+5 ]=SA[(j*30)+4+5 ]=INF;
1258     SA[(j*30)+10]=SA[(j*30)+1+10]=SA[(j*30)+2+10]=SA[(j*30)+3+10]=SA[(j*30)+4+10]=INF;
1259     SA[(j*30)+15]=SA[(j*30)+1+15]=SA[(j*30)+2+15]=SA[(j*30)+3+15]=SA[(j*30)+4+15]=INF;
1260     SA[(j*30)+20]=SA[(j*30)+1+20]=SA[(j*30)+2+20]=SA[(j*30)+3+20]=SA[(j*30)+4+20]=INF;
1261     SA[(j*30)+25]=SA[(j*30)+1+25]=SA[(j*30)+2+25]=SA[(j*30)+3+25]=SA[(j*30)+4+25]=INF;
1262   }
1263
1264   i=10 ;
1265   i_length= n1 - 9  ;
1266   while(i < i_length) {
1267     int di1,di2,di3,di4;
1268     int idx=i%5;
1269     int idx_1=(i-1)%5;
1270     int idx_2=(i-2)%5;
1271     int idx_3=(i-3)%5;
1272     int idx_4=(i-4)%5;
1273     di1 = access_s1[5][i]   - access_s1[4][i-1];
1274     di2 = access_s1[5][i-1] - access_s1[4][i-2] + di1;
1275     di3 = access_s1[5][i-2] - access_s1[4][i-3] + di2;
1276     di4 = access_s1[5][i-3] - access_s1[4][i-4] + di3;
1277     di1=MIN2(di1,maxPenalty[0]);
1278     di2=MIN2(di2,maxPenalty[1]);
1279     di3=MIN2(di3,maxPenalty[2]);
1280     di4=MIN2(di4,maxPenalty[3]);
1281     j=n2 - 9;
1282     while (--j > 9) {
1283       int dj1,dj2,dj3,dj4;
1284       dj1=DJ[0][j];
1285       dj2=DJ[1][j];
1286       dj3=DJ[2][j];
1287       dj4=DJ[3][j];
1288       int type2, type,temp;
1289       type  = pair[S1[i]][S2[j]];
1290       /**
1291       *** Start duplex
1292       **/
1293       SA[LCI(idx,j,n2)] = type ? P->DuplexInit + access_s1[1][i] + access_s2[1][j] : INF;
1294       /**
1295       *** update lin bx by linx liny matrix
1296       **/
1297       type2=pair[S2[j+1]][S1[i-1]];
1298       /**
1299       *** start/extend interior loop
1300       **/
1301       SA[LINI(idx,j,n2)]=MIN2(SA[LCI(idx_1,j+1,n2)]+P->mismatchI[type2][SS2[j]][SS1[i]]+di1+dj1+iopen+iext_s,
1302                               SA[LINI(idx_1,j,n2)]+iext_ass + di1);
1303
1304       /**
1305       *** start/extend nx1 target
1306       *** use same type2 as for in
1307       **/
1308       SA[LINIX(idx,j,n2)]=MIN2(SA[LCI(idx_1,j+1,n2)]+P->mismatch1nI[type2][SS2[j]][SS1[i]]+di1+dj1+iopen+iext_s,
1309                                SA[LINIX(idx_1,j,n2)]+iext_ass + di1);
1310       /**
1311       *** start/extend 1xn target
1312       *** use same type2 as for in
1313       **/
1314       SA[LINIY(idx,j,n2)]=MIN2(SA[LCI(idx_1,j+1,n2)]+P->mismatch1nI[type2][SS2[j]][SS1[i]]+di1+dj1+iopen+iext_s,
1315                                SA[LINIY(idx,j+1,n2)]+iext_ass + dj1);
1316       /**
1317       *** extend interior loop
1318       **/
1319       SA[LINI(idx,j,n2)]=MIN2(SA[LINI(idx,j,n2)],SA[LINI(idx,j+1,n2)]+iext_ass + dj1);
1320       SA[LINI(idx,j,n2)]=MIN2(SA[LINI(idx,j,n2)],SA[LINI(idx_1,j+1,n2)]+iext_s + di1 + dj1);
1321       /**
1322       *** start/extend bulge target
1323       **/
1324       type2=pair[S2[j]][S1[i-1]];
1325       SA[LBXI(idx,j,n2)]=MIN2(SA[LBXI(idx_1,j,n2)]+bext + di1, SA[LCI(idx_1,j,n2)]+bopen+bext+(type2>2?P->TerminalAU:0) + di1);
1326       /**
1327       *** start/extend bulge query
1328       **/
1329       type2=pair[S2[j+1]][S1[i]];
1330       SA[LBYI(idx,j,n2)]=MIN2(SA[LBYI(idx,j+1,n2)]+bext + dj1 , SA[LCI(idx,j+1,n2)]+bopen+bext+(type2>2?P->TerminalAU:0)+ dj1);
1331       /**
1332       ***end update recursion
1333       **/
1334       if(!type){continue;}
1335       /**
1336       *** stack extension
1337       **/
1338       SA[LCI(idx,j,n2)]+= E_ExtLoop(type, SS1[i-1] , SS2[j+1], P);
1339       /**
1340       *** stack extension
1341       **/
1342       if((type2=pair[S1[i-1]][S2[j+1]]))
1343         SA[LCI(idx,j,n2)]=MIN2(SA[LCI(idx_1,j+1,n2)]+P->stack[rtype[type]][type2]+di1+dj1, SA[LCI(idx,j,n2)]);
1344       /**
1345       *** 1x0 / 0x1 stack extension
1346       **/
1347       if((type2=pair[S1[i-1]][S2[j+2]]))
1348         SA[LCI(idx,j,n2)]=MIN2(SA[LCI(idx_1,j+2,n2)]+P->bulge[1]+P->stack[rtype[type]][type2]+di1+dj2, SA[LCI(idx,j,n2)]);
1349       if((type2=pair[S1[i-2]][S2[j+1]]))
1350         SA[LCI(idx,j,n2)]=MIN2(SA[LCI(idx_2,j+1,n2)]+P->bulge[1]+P->stack[type2][rtype[type]]+di2+dj1, SA[LCI(idx,j,n2)]);
1351       /**
1352       *** 1x1 / 2x2 stack extension
1353       **/
1354       if((type2=pair[S1[i-2]][S2[j+2]]))
1355         SA[LCI(idx,j,n2)]=MIN2(SA[LCI(idx_2,j+2,n2)]+P->int11[type2][rtype[type]][SS1[i-1]][SS2[j+1]]+di2+dj2, SA[LCI(idx,j,n2)]);
1356       if((type2 = pair[S1[i-3]][S2[j+3]]))
1357         SA[LCI(idx,j,n2)]=MIN2(SA[LCI(idx_3,j+3,n2)]+P->int22[type2][rtype[type]][SS1[i-2]][SS1[i-1]][SS2[j+1]][SS2[j+2]]+di3+dj3, SA[LCI(idx,j,n2)]);
1358       /**
1359       *** 1x2 / 2x1 stack extension
1360       *** E_IntLoop(1,2,type2, rtype[type],SS1[i-1], SS2[j+2], SS1[i-1], SS2[j+1], P) corresponds to
1361       *** P->int21[rtype[type]][type2][SS2[j+2]][SS1[i-1]][SS1[i-1]]
1362       **/
1363       if((type2 = pair[S1[i-3]][S2[j+2]]))
1364         SA[LCI(idx,j,n2)]=MIN2(SA[LCI(idx_3,j+2,n2)]+P->int21[rtype[type]][type2][SS2[j+1]][SS1[i-2]][SS1[i-1]]+di3+dj2, SA[LCI(idx,j,n2)]);
1365       if((type2 = pair[S1[i-2]][S2[j+3]]))
1366         SA[LCI(idx,j,n2)]=MIN2(SA[LCI(idx_2,j+3,n2)]+P->int21[type2][rtype[type]][SS1[i-1]][SS2[j+1]][SS2[j+2]]+di2+dj3, SA[LCI(idx,j,n2)]);
1367       /**
1368       *** 2x3 / 3x2 stack extension
1369       **/
1370       if((type2 = pair[S1[i-4]][S2[j+3]]))
1371         SA[LCI(idx,j,n2)]=MIN2(SA[LCI(idx_4,j+3,n2)]+P->internal_loop[5]+P->ninio[2]+
1372                                P->mismatch23I[type2][SS1[i-3]][SS2[j+2]]+P->mismatch23I[rtype[type]][SS2[j+1]][SS1[i-1]]+di4+dj3, SA[LCI(idx,j,n2)]);
1373       if((type2 = pair[S1[i-3]][S2[j+4]]))
1374         SA[LCI(idx,j,n2)]=MIN2(SA[LCI(idx_3,j+4,n2)]+P->internal_loop[5]+P->ninio[2]+
1375                                P->mismatch23I[type2][SS1[i-2]][SS2[j+3]]+P->mismatch23I[rtype[type]][SS2[j+1]][SS1[i-1]]+di3+dj4, SA[LCI(idx,j,n2)]);
1376       /**
1377       *** So now we have to handle 1x3, 3x1, 3x3, and mxn m,n > 3
1378       **/
1379       /**
1380       *** 3x3 or more
1381       **/
1382       SA[LCI(idx,j,n2)]=MIN2(SA[LINI(idx_3,j+3,n2)]+P->mismatchI[rtype[type]][SS1[i-1]][SS2[j+1]]+2*iext_s+di3+dj3,SA[LCI(idx,j,n2)]);
1383       /**
1384       *** 2xn or more
1385       **/
1386       SA[LCI(idx,j,n2)]=MIN2(SA[LINI(idx_4,j+2,n2)]+P->mismatchI[rtype[type]][SS1[i-1]][SS2[j+1]]+iext_s+2*iext_ass+di4+dj2, SA[LCI(idx,j,n2)]);
1387       /**
1388       *** nx2 or more
1389       **/
1390       SA[LCI(idx,j,n2)]=MIN2(SA[LINI(idx_2,j+4,n2)]+P->mismatchI[rtype[type]][SS1[i-1]][SS2[j+1]]+iext_s+2*iext_ass+di2+dj4, SA[LCI(idx,j,n2)]);
1391       /**
1392       *** nx1 n>2
1393       **/
1394       SA[LCI(idx,j,n2)]=MIN2(SA[LINIX(idx_3,j+1,n2)]+P->mismatch1nI[rtype[type]][SS1[i-1]][SS2[j+1]]+iext_ass+iext_ass+di3+dj1, SA[LCI(idx,j,n2)]);
1395       /**
1396       *** 1xn n>2
1397       **/
1398       SA[LCI(idx,j,n2)]=MIN2(SA[LINIY(idx_1,j+3,n2)]+P->mismatch1nI[rtype[type]][SS1[i-1]][SS2[j+1]]+iext_ass+iext_ass+dj3+di1, SA[LCI(idx,j,n2)]);
1399       /**
1400       *** nx0 n>1
1401       **/
1402       int bAU;
1403       bAU=(type>2?P->TerminalAU:0);
1404       SA[LCI(idx,j,n2)]=MIN2(SA[LBXI(idx_2,j+1,n2)]+di2+dj1+bext+bAU,SA[LCI(idx,j,n2)]);
1405       /**
1406       *** 0xn n>1
1407       **/
1408       SA[LCI(idx,j,n2)]=MIN2(SA[LBYI(idx_1,j+2,n2)]+di1+dj2+bext+bAU,SA[LCI(idx,j,n2)]);
1409       temp=min_colonne;
1410       /**
1411       *** (type>2?P->TerminalAU:0)+
1412       *** P->dangle3[rtype[type]][SS1[i+1]]+
1413       *** P->dangle5[rtype[type]][SS2[j-1]],
1414       **/
1415       min_colonne=MIN2(SA[LCI(idx,j,n2)]+E_ExtLoop(rtype[type], SS2[j-1] , SS1[i+1] , P), min_colonne);
1416
1417       if(temp>min_colonne){
1418         min_j_colonne=j;
1419       }
1420
1421       /* ---------------------------------------------------------------------end update */
1422     }
1423     if(max>=min_colonne){
1424       max=min_colonne;
1425       max_pos=i;
1426       max_pos_j=min_j_colonne;
1427     }
1428     position[i+delta]=min_colonne;min_colonne=INF;
1429     position_j[i+delta]=min_j_colonne;
1430     i++;
1431   }
1432   /* printf("MAX: %d",max); */
1433   free(S1); free(S2); free(SS1); free(SS2);free(SA);
1434   if(max<threshold){
1435     find_max_XS(position, position_j, delta, threshold, alignment_length, s1, s2, access_s1, access_s2, fast,il_a, il_b,b_a, b_b);
1436   }
1437   if(max<INF){
1438     plot_max_XS(max, max_pos, max_pos_j, alignment_length, s1, s2, access_s1, access_s2,fast, il_a, il_b, b_a, b_b);
1439   }
1440   for (i=0; i<=3; i++) {
1441     free(DJ[i]);
1442   }
1443   free(DJ);
1444   free(position);
1445   free(position_j);
1446   return NULL;
1447 }
1448
1449 PRIVATE void find_max_XS(const int *position, const int *position_j,const int delta, const int threshold, const int alignment_length, const char *s1, const char *s2, const int **access_s1, const int **access_s2, const int fast,const int il_a, const int il_b, const int b_a, const int b_b){
1450
1451   int n_max;n_max=16;
1452   int pos=n1-9;
1453   if(fast==1){
1454     while(10 < pos--){
1455       printf("%d %d\n",pos,position[pos+delta]);
1456       int temp_min=0;
1457       if(position[pos+delta]<(threshold)){
1458         int search_range;
1459         search_range=delta+1;
1460         while(--search_range){
1461           if(position[pos+delta-search_range]<=position[pos+delta-temp_min]){
1462             temp_min=search_range;
1463           }
1464         }
1465         pos-=temp_min;
1466         int max_pos_j;
1467         max_pos_j=position_j[pos+delta];
1468         int max;
1469         max=position[pos+delta];
1470         printf("target upper bound %d: query lower bound %d  (%5.2f) \n", pos-10, max_pos_j-10, ((double)max)/100);
1471         pos=MAX2(10,pos-delta);
1472       }
1473     }
1474   }
1475   else if(fast==2){
1476     pos=n1-9;
1477     while(10 < pos--){
1478       int temp_min=0;
1479       if(position[pos+delta]<(threshold)){
1480         int search_range;
1481         search_range=delta+1;
1482         while(--search_range){
1483           if(position[pos+delta-search_range]<=position[pos+delta-temp_min]){
1484             temp_min=search_range;
1485           }
1486         }
1487         pos-=temp_min;
1488         int max_pos_j;
1489         max_pos_j=position_j[pos+delta];
1490         /* max_pos_j und pos entsprechen die realen position
1491             in der erweiterten sequenz.
1492            pos=1 -> position 1 in the sequence (and not 0 like in C)
1493            max_pos_j -> position 1 in the sequence ( not 0 like in C)
1494         */
1495         int begin_t=MAX2(11, pos-alignment_length+1);/* 10 */
1496         int end_t  =MIN2(n1-10, pos+1);
1497         int begin_q=MAX2(11, max_pos_j-1); /* 10 */
1498         int end_q  =MIN2(n2-10, max_pos_j+alignment_length-1);
1499         char *s3 = (char*) space(sizeof(char)*(end_t - begin_t +2 + 20));
1500         char *s4 = (char*) space(sizeof(char)*(end_q - begin_q +2 + 20));
1501         strcpy(s3,"NNNNNNNNNN");strcpy(s4,"NNNNNNNNNN");
1502         strncat(s3, (s1+begin_t-1),  end_t - begin_t +1);
1503         strncat(s4, (s2+begin_q-1) , end_q - begin_q +1);
1504         strcat(s3,"NNNNNNNNNN");strcat(s4,"NNNNNNNNNN");
1505         s3[end_t -begin_t +1 +20 ]='\0';
1506         s4[end_q -begin_q +1 +20]='\0';
1507         duplexT test;
1508         test = fduplexfold_XS(s3, s4, access_s1, access_s2, end_t, begin_q,threshold, il_a, il_b, b_a, b_b);
1509         if(test.energy * 100 < threshold){
1510           int l1=strchr(test.structure, '&')-test.structure;
1511           printf(" %s %3d,%-3d : %3d,%-3d (%5.2f = %5.2f + %5.2f + %5.2f) [%5.2f = %5.2f + %5.2f + %5.2f]  %d \n", test.structure,
1512                  begin_t-10+test.i-l1-10,
1513                  begin_t-10+test.i-1-10,
1514                  begin_q-10 + test.j-1-10 ,
1515                  (begin_q -11) + test.j + strlen(test.structure)-l1-2-10,
1516                  test.ddG, test.energy, test.opening_backtrack_x, test.opening_backtrack_y,
1517                  test.energy_backtrack + test.dG1 + test.dG2, test.energy_backtrack , test.dG1 , test.dG2,
1518                  position[pos+delta]);
1519           pos=MAX2(10,pos-delta);
1520           free(s3);free(s4);
1521           free(test.structure);
1522         }
1523         else{
1524           free(s3);free(s4);
1525         }
1526       }
1527     }
1528   }
1529   else{
1530     pos=n1-9;
1531     while( pos-- > 10 ){
1532       int temp_min=0;
1533       if(position[pos+delta]<(threshold)){
1534         int search_range;
1535         search_range=delta+1;
1536         while(--search_range){
1537           if(position[pos+delta-search_range]<=position[pos+delta-temp_min]){
1538             temp_min=search_range;
1539           }
1540         }
1541         pos-=temp_min; /* position on i */
1542         int max_pos_j;
1543         max_pos_j=position_j[pos+delta]; /* position on j */
1544         int begin_t=MAX2(11,pos-alignment_length);
1545         int end_t  =MIN2(n1-10, pos+1);
1546         int begin_q=MAX2(11,max_pos_j-1);
1547         int end_q  =MIN2(n2-10,max_pos_j+alignment_length-1);
1548         int i_flag;
1549         int j_flag;
1550         i_flag = (end_t   ==  pos+1?1:0);
1551         j_flag = (begin_q == max_pos_j-1?1:0);
1552         char *s3 = (char*) space(sizeof(char)*(end_t - begin_t +2));
1553         char *s4 = (char*) space(sizeof(char)*(end_q - begin_q +2));
1554
1555         strncpy(s3, (s1+begin_t),  end_t - begin_t+1);
1556         strncpy(s4, (s2+begin_q) , end_q - begin_q+1);
1557         s3[end_t -begin_t +1 ]='\0';
1558         s4[end_q -begin_q +1 ]='\0';
1559         duplexT test;
1560         test = duplexfold_XS(s3,s4,access_s1,access_s2,pos, max_pos_j,threshold,i_flag,j_flag);
1561         if(test.energy * 100 < threshold){
1562           printf("%s %3d,%-3d : %3d,%-3d (%5.2f = %5.2f + %5.2f + %5.2f)\n", test.structure,
1563                  test.tb,test.te,test.qb,test.qe, test.ddG, test.energy, test.dG1, test.dG2);
1564           pos=MAX2(10,pos-delta);
1565         }
1566         free(s3);free(s4);
1567         free(test.structure);
1568       }
1569     }
1570   }
1571 }
1572
1573 #if 0
1574 PRIVATE int compare(const void *sub1, const void *sub2) {
1575   int d;
1576   if (((duplexT *) sub1)->ddG > ((duplexT *) sub2)->ddG)
1577     return 1;
1578   if (((duplexT *) sub1)->ddG < ((duplexT *) sub2)->ddG)
1579     return -1;
1580   d = ((duplexT *) sub1)->i - ((duplexT *) sub2)->i;
1581   if (d!=0) return d;
1582   return  ((duplexT *) sub1)->j - ((duplexT *) sub2)->j;
1583 }
1584 #endif
1585
1586 PRIVATE void plot_max_XS(const int max, const int max_pos, const int max_pos_j, const int alignment_length, const char *s1, const char *s2, const int ** access_s1, const int ** access_s2, const int fast,const int il_a, const int il_b, const int b_a, const int b_b)
1587 {
1588   if(fast==1){
1589     printf("target upper bound %d: query lower bound %d (%5.2f)\n", max_pos-3, max_pos_j, ((double)max)/100);
1590   }
1591   else if(fast==2){
1592     int begin_t=MAX2(11, max_pos-alignment_length+1);/* 10 */
1593     int end_t  =MIN2(n1-10, max_pos+1);
1594     int begin_q=MAX2(11, max_pos_j-1); /* 10 */
1595     int end_q  =MIN2(n2-10, max_pos_j+alignment_length-1);
1596     char *s3 = (char*) space(sizeof(char)*(end_t - begin_t +2 + 20));
1597     char *s4 = (char*) space(sizeof(char)*(end_q - begin_q +2 + 20));
1598     strcpy(s3,"NNNNNNNNNN");strcpy(s4,"NNNNNNNNNN");
1599     strncat(s3, (s1+begin_t-1),  end_t - begin_t +1);
1600     strncat(s4, (s2+begin_q-1) , end_q - begin_q +1);
1601     strcat(s3,"NNNNNNNNNN");strcat(s4,"NNNNNNNNNN");
1602     s3[end_t -begin_t +1 +20 ]='\0';
1603     s4[end_q -begin_q +1 +20]='\0';
1604     duplexT test;
1605     test = fduplexfold_XS(s3, s4, access_s1, access_s2, end_t, begin_q, INF, il_a, il_b, b_a, b_b);
1606     int l1=strchr(test.structure, '&')-test.structure;
1607     printf("%s %3d,%-3d : %3d,%-3d (%5.2f = %5.2f + %5.2f + %5.2f) [%5.2f = %5.2f + %5.2f + %5.2f] \n", test.structure,
1608            begin_t-10+test.i-l1-10,
1609            begin_t-10+test.i-1-10,
1610            begin_q-10 + test.j-1-10 ,
1611            (begin_q -11) + test.j + strlen(test.structure)-l1-2-10,
1612            test.ddG, test.energy, test.opening_backtrack_x, test.opening_backtrack_y,
1613            test.energy_backtrack + test.dG1 + test.dG2, test.energy_backtrack , test.dG1 , test.dG2);
1614     free(s3);free(s4);
1615     free(test.structure);
1616   }
1617   else{
1618     int begin_t=MAX2(11,max_pos-alignment_length);
1619     int end_t  =MIN2(n1-10, max_pos+1);
1620     int begin_q=MAX2(11, max_pos_j-1);
1621     int end_q  =MIN2(n2-10,max_pos_j+alignment_length-1);
1622     int i_flag;
1623     int j_flag;
1624     i_flag = (end_t   == max_pos+1?1:0);
1625     j_flag = (begin_q == max_pos_j-1?1:0);
1626     char *s3 = (char*) space(sizeof(char)*(end_t - begin_t +2)); /* +1 for \0 +1 for distance */
1627     char *s4 = (char*) space(sizeof(char)*(end_q - begin_q +2));
1628
1629     strncpy(s3, (s1+begin_t-1),  end_t - begin_t+1);/* -1 to go from  */
1630     strncpy(s4, (s2+begin_q-1) , end_q - begin_q+1);/* -1 to go from  */
1631     s3[end_t -begin_t +1 ]='\0';/*  */
1632     s4[end_q -begin_q +1 ]='\0';
1633     duplexT test;
1634     test = duplexfold_XS(s3,s4,access_s1,access_s2,max_pos, max_pos_j,INF,i_flag,j_flag);
1635     printf("%s %3d,%-3d : %3d,%-3d (%5.2f = %5.2f + %5.2f + %5.2f)\n", test.structure,
1636            test.tb,test.te,test.qb,test.qe, test.ddG, test.energy, test.dG1, test.dG2);
1637     free(s3);free(s4);free(test.structure);
1638   }
1639 }
1640
1641
1642 /*---------------------------------------------------------duplexfold----------------------------------------------------------------------------------*/
1643
1644
1645 PRIVATE duplexT duplexfold(const char *s1, const char *s2, const int extension_cost) {
1646   int i, j, l1, Emin=INF, i_min=0, j_min=0;
1647   char *struc;
1648   duplexT mfe;
1649
1650   n3 = (int) strlen(s1);
1651   n4 = (int) strlen(s2);
1652
1653   if ((!P) || (fabs(P->temperature - temperature)>1e-6)) {
1654     update_fold_params();  if(P) free(P); P = scale_parameters();
1655     make_pair_matrix();
1656   }
1657
1658   c = (int **) space(sizeof(int *) * (n3+1));
1659   for (i=0; i<=n3; i++) c[i] = (int *) space(sizeof(int) * (n4+1));
1660   encode_seqs(s1, s2);
1661   for (i=1; i<=n3; i++) {
1662     for (j=n4; j>0; j--) {
1663       int type, type2, E, k,l;
1664       type = pair[S1[i]][S2[j]];
1665       c[i][j] = type ? P->DuplexInit +2 * extension_cost: INF;
1666       if (!type) continue;
1667       /**
1668       ***       if (i>1)  c[i][j] += P->dangle5[type][SS1[i-1]]+ extension_cost;
1669       ***       if (j<n4) c[i][j] += P->dangle3[type][SS2[j+1]]+ extension_cost;
1670       ***       if (type>2) c[i][j] += P->TerminalAU;
1671       **/
1672       c[i][j] += E_ExtLoop(type, (i>1) ? SS1[i-1] : -1, (j<n4) ? SS2[j+1] : -1, P);
1673       for (k=i-1; k>0 && k>i-MAXLOOP-2; k--) {
1674         for (l=j+1; l<=n4; l++) {
1675           if (i-k+l-j-2>MAXLOOP) break;
1676           type2 = pair[S1[k]][S2[l]];
1677           if (!type2) continue;
1678           E = E_IntLoop(i-k-1, l-j-1, type2, rtype[type],
1679                         SS1[k+1], SS2[l-1], SS1[i-1], SS2[j+1],P)+(i-k+l-j)*extension_cost;
1680           c[i][j] = MIN2(c[i][j], c[k][l]+E);
1681         }
1682       }
1683       E = c[i][j];
1684       /**
1685       ***      if (i<n3) E += P->dangle3[rtype[type]][SS1[i+1]]+extension_cost;
1686       ***      if (j>1)  E += P->dangle5[rtype[type]][SS2[j-1]]+extension_cost;
1687       ***      if (type>2) E += P->TerminalAU;
1688       ***
1689       **/
1690       E += E_ExtLoop(rtype[type], (j > 1) ? SS2[j-1] : -1, (i<n3) ? SS1[i+1] : -1, P);
1691       if (E<Emin) {
1692         Emin=E; i_min=i; j_min=j;
1693       }
1694     }
1695   }
1696   struc = backtrack(i_min, j_min, extension_cost);
1697   if (i_min<n3) i_min++;
1698   if (j_min>1 ) j_min--;
1699   l1 = strchr(struc, '&')-struc;
1700   int size;
1701   size=strlen(struc)-1;
1702   Emin-= size * (extension_cost);
1703   mfe.i = i_min;
1704   mfe.j = j_min;
1705   mfe.energy = (double) Emin/100.;
1706   mfe.structure = struc;
1707   for (i=0; i<=n3; i++) free(c[i]);
1708   free(c);
1709   free(S1); free(S2); free(SS1); free(SS2);
1710   return mfe;
1711 }
1712
1713 PRIVATE char *backtrack(int i, int j, const int extension_cost) {
1714   /* backtrack structure going backwards from i, and forwards from j
1715      return structure in bracket notation with & as separator */
1716   int k, l, type, type2, E, traced, i0, j0;
1717   char *st1, *st2, *struc;
1718
1719   st1 = (char *) space(sizeof(char)*(n3+1));
1720   st2 = (char *) space(sizeof(char)*(n4+1));
1721
1722   i0=MIN2(i+1,n3); j0=MAX2(j-1,1);
1723
1724   while (i>0 && j<=n4) {
1725     E = c[i][j]; traced=0;
1726     st1[i-1] = '(';
1727     st2[j-1] = ')';
1728     type = pair[S1[i]][S2[j]];
1729     if (!type) nrerror("backtrack failed in fold duplex");
1730     for (k=i-1; k>0 && k>i-MAXLOOP-2; k--) {
1731       for (l=j+1; l<=n4; l++) {
1732         int LE;
1733         if (i-k+l-j-2>MAXLOOP) break;
1734         type2 = pair[S1[k]][S2[l]];
1735         if (!type2) continue;
1736         LE = E_IntLoop(i-k-1, l-j-1, type2, rtype[type],
1737                        SS1[k+1], SS2[l-1], SS1[i-1], SS2[j+1],P)+(i-k+l-j)*extension_cost;
1738         if (E == c[k][l]+LE) {
1739           traced=1;
1740           i=k; j=l;
1741           break;
1742         }
1743       }
1744       if (traced) break;
1745     }
1746     if (!traced) {
1747
1748       E -=  E_ExtLoop(type, (i>1) ? SS1[i-1] : -1, (j<n4) ? SS2[j+1] : -1, P);
1749       /**
1750       ***      if (i>1) E -= P->dangle5[type][SS1[i-1]]+extension_cost;
1751       ***      if (j<n4) E -= P->dangle3[type][SS2[j+1]]+extension_cost;
1752       ***      if (type>2) E -= P->TerminalAU;
1753       **/
1754       if (E != P->DuplexInit+2*extension_cost) {
1755         nrerror("backtrack failed in fold duplex");
1756       } else break;
1757     }
1758   }
1759   if (i>1)  i--;
1760   if (j<n4) j++;
1761
1762   struc = (char *) space(i0-i+1+j-j0+1+2);
1763   for (k=MAX2(i,1); k<=i0; k++) if (!st1[k-1]) st1[k-1] = '.';
1764   for (k=j0; k<=j; k++) if (!st2[k-1]) st2[k-1] = '.';
1765   strcpy(struc, st1+MAX2(i-1,0));
1766   strcat(struc, "&");
1767   strcat(struc, st2+j0-1);
1768   /* printf("%s %3d,%-3d : %3d,%-3d\n", struc, i,i0,j0,j);  */
1769   free(st1); free(st2);
1770   return struc;
1771 }
1772
1773 PRIVATE duplexT fduplexfold(const char *s1, const char *s2, const int extension_cost,const int il_a, const int il_b, const int b_a, const int b_b) {
1774   int i, j, Emin, i_min, j_min,l1;
1775   duplexT mfe;
1776   char *struc;
1777   int bopen=b_b;
1778   int bext=b_a+extension_cost;
1779   int iopen=il_b;
1780   int iext_s=2*(il_a+extension_cost);/* iext_s 2 nt nucleotide extension of interior loop, on i and j side */
1781   int iext_ass=50+il_a+extension_cost;/* iext_ass assymetric extension of interior loop, either on i or on j side. */
1782   int min_colonne=INF; /* enthaelt das maximum einer kolonne */
1783   int i_length;
1784   int max_pos;/* get position of the best hit */
1785   int max_pos_j;
1786   int temp=INF;
1787   int min_j_colonne;
1788   int max=INF;
1789   /* FOLLOWING NEXT 4 LINE DEFINES AN ARRAY CONTAINING POSITION OF THE SUBOPT IN S1 */
1790
1791   n3 = (int) strlen(s1);
1792   n4 = (int) strlen(s2);
1793   /* delta_check is the minimal distance allowed for two hits to be accepted */
1794   /* if both hits are closer, reject the smaller ( in term of position)  hits  */
1795   /* i want to implement a function that, given a position in a long sequence and a small sequence, */
1796   /* duplexfold them at this position and report the result at the command line */
1797   /* for this i first need to rewrite backtrack in order to remove the printf functio */
1798   /* END OF DEFINITION FOR NEEDED SUBOPT DATA  */
1799   if ((!P) || (fabs(P->temperature - temperature)>1e-6)) {
1800     update_fold_params();  if(P) free(P); P = scale_parameters();
1801     make_pair_matrix();
1802   }
1803   /*local c array initialization---------------------------------------------*/
1804   c  = (int**)  space(sizeof(int *) * (n3+1));
1805   in = (int**) space (sizeof(int *) * (n3+1));
1806   bx = (int**) space (sizeof(int *) * (n3+1));
1807   by = (int**) space (sizeof(int *) * (n3+1));
1808   inx= (int**) space (sizeof(int *) * (n3+1));
1809   iny= (int**) space (sizeof(int *) * (n3+1));
1810   for (i=0; i<=n3; i++){
1811     c[i]  = (int *) space(sizeof(int) * (n4+1));
1812     in[i] = (int *) space(sizeof(int) * (n4+1));
1813     bx[i] = (int *) space(sizeof(int) * (n4+1));
1814     by[i] = (int *) space(sizeof(int) * (n4+1));
1815     inx[i] = (int *) space(sizeof(int) * (n4+1));
1816     iny[i] = (int *) space(sizeof(int) * (n4+1));
1817   }
1818   /*-------------------------------------------------------------------------*/
1819   /*end of array initialisation----------------------------------*/
1820   /*maybe int *** would be better*/
1821   encode_seqs(s1,s2);
1822   /* ------------------------------------------matrix initialisierung */
1823   for(i=0; i<n3; i++){
1824     for(j=0; j<n4; j++){
1825       in[i][j]=INF;/* no in before  1 */
1826       c[i][j] =INF; /* no bulge and no in before n2 */
1827       bx[i][j]=INF;/* no bulge before 1 */
1828       by[i][j]=INF;
1829       inx[i][j]=INF;/* no bulge before 1 */
1830       iny[i][j]=INF;
1831     }
1832   }
1833
1834   /*--------------------------------------------------------local array*/
1835
1836
1837   /* -------------------------------------------------------------matrix initialisierung */
1838   i=11;
1839   i_length=n3-9;
1840   while(i < i_length) {
1841     j=n4-9;
1842     min_colonne=INF;
1843     while (10 < --j) {
1844       int type, type2;
1845       type = pair[S1[i]][S2[j]];
1846       /**
1847       *** Start duplex
1848       **/
1849       c[i][j]=type ? P->DuplexInit + 2*extension_cost : INF;
1850       /**
1851       *** update lin bx by linx liny matrix
1852       **/
1853       type2=pair[S2[j+1]][S1[i-1]];
1854       /**
1855       *** start/extend interior loop
1856       **/
1857       in[i][j]=MIN2(c[i - 1][j+1]+P->mismatchI[type2][SS2[j]][SS1[i]]+iopen+iext_s, in[i - 1][j]+iext_ass);
1858       /**
1859       *** start/extend nx1 target
1860       *** use same type2 as for in
1861       **/
1862       inx[i][j]=MIN2(c[i-1][j+1]+P->mismatch1nI[type2][SS2[j]][SS1[i]]+iopen+iext_s,
1863                      inx[i-1][j]+iext_ass);
1864       /**
1865       *** start/extend 1xn target
1866       *** use same type2 as for in
1867       **/
1868       iny[i][j]=MIN2(c[i-1][j+1]+P->mismatch1nI[type2][SS2[j]][SS1[i]]+iopen+iext_s,
1869                      iny[i][j+1]+iext_ass);
1870       /**
1871       *** extend interior loop
1872       **/
1873       in[i][j]=MIN2(in[i][j],in[i][j+1]+iext_ass);
1874       in[i][j]=MIN2(in[i][j],in[i - 1][j+1]+iext_s);
1875       /**
1876       *** start/extend bulge target
1877       **/
1878       type2=pair[S2[j]][S1[i-1]];
1879       bx[i][j]=MIN2(bx[i - 1][j]+bext, c[i - 1][j]+bopen+bext+(type2>2?P->TerminalAU:0));
1880       /**
1881       *** start/extend bulge query
1882       **/
1883       type2=pair[S2[j+1]][S1[i]];
1884       by[i][j]=MIN2(by[i][j+1]+bext, c[i][j+1]+bopen+bext+(type2>2?P->TerminalAU:0));
1885       /**
1886       ***end update recursion
1887       ***######################## Start stack extension##############################
1888       **/
1889       if(!type){continue;}
1890       c[i][j]+=E_ExtLoop(type, SS1[i-1] , SS2[j+1], P) + 2*extension_cost;
1891       /**
1892       *** stack extension
1893       **/
1894       if((type2=pair[S1[i-1]][S2[j+1]]))
1895         c[i][j]=MIN2(c[i - 1][j+1]+P->stack[rtype[type]][type2]+2*extension_cost, c[i][j]);
1896       /**
1897       *** 1x0 / 0x1 stack extension
1898       **/
1899       type2=pair[S1[i-1]][S2[j+2]];
1900       c[i][j]=MIN2(c[i - 1][j+2]+P->bulge[1]+P->stack[rtype[type]][type2]+3*extension_cost,c[i][j]);
1901       type2=pair[S1[i-2]][S2[j+1]];
1902       c[i][j]=MIN2(c[i - 2][j+1]+P->bulge[1]+P->stack[type2][rtype[type]]+3*extension_cost,c[i][j]);
1903       /**
1904       *** 1x1 / 2x2 stack extension
1905       **/
1906       type2=pair[S1[i-2]][S2[j+2]];
1907       c[i][j]=MIN2(c[i - 2][j+2]+P->int11[type2][rtype[type]][SS1[i-1]][SS2[j+1]]+4*extension_cost, c[i][j]);
1908       type2 = pair[S1[i-3]][S2[j+3]];
1909       c[i][j]=MIN2(c[i - 3][j+3]+P->int22[type2][rtype[type]][SS1[i-2]][SS1[i-1]][SS2[j+1]][SS2[j+2]]+6*extension_cost,c[i][j]);
1910       /**
1911       *** 1x2 / 2x1 stack extension
1912       *** E_IntLoop(1,2,type2, rtype[type],SS1[i-1], SS2[j+2], SS1[i-1], SS2[j+1], P) corresponds to
1913       *** P->int21[rtype[type]][type2][SS2[j+2]][SS1[i-1]][SS1[i-1]]
1914       **/
1915       type2 = pair[S1[i-3]][S2[j+2]];
1916       c[i][j]=MIN2(c[i - 3][j+2]+P->int21[rtype[type]][type2][SS2[j+1]][SS1[i-2]][SS1[i-1]]+5*extension_cost, c[i][j]);
1917       type2 = pair[S1[i-2]][S2[j+3]];
1918       c[i][j]=MIN2(c[i - 2][j+3]+P->int21[type2][rtype[type]][SS1[i-1]][SS2[j+1]][SS2[j+2]]+5*extension_cost, c[i][j]);
1919
1920       /**
1921       *** 2x3 / 3x2 stack extension
1922       **/
1923       if((type2 = pair[S1[i-4]][S2[j+3]]))
1924         c[i][j]=MIN2(c[i - 4][j+3]+P->internal_loop[5]+P->ninio[2]+
1925                      P->mismatch23I[type2][SS1[i-3]][SS2[j+2]]+P->mismatch23I[rtype[type]][SS2[j+1]][SS1[i-1]]+7*extension_cost, c[i][j]);
1926       if((type2 = pair[S1[i-3]][S2[j+4]]))
1927         c[i][j]=MIN2(c[i - 3][j+4]+P->internal_loop[5]+P->ninio[2]+
1928                      P->mismatch23I[type2][SS1[i-2]][SS2[j+3]]+P->mismatch23I[rtype[type]][SS2[j+1]][SS1[i-1]]+7*extension_cost, c[i][j]);
1929       /**
1930       *** So now we have to handle 1x3, 3x1, 3x3, and mxn m,n > 3
1931       **/
1932       /**
1933       *** 3x3 or more
1934       **/
1935       c[i][j]=MIN2(in[i - 3][j+3]+P->mismatchI[rtype[type]][SS1[i-1]][SS2[j+1]]+2*iext_s+2*extension_cost, c[i][j]);
1936       /**
1937       *** 2xn or more
1938       **/
1939       c[i][j]=MIN2(in[i - 4][j+2]+P->mismatchI[rtype[type]][SS1[i-1]][SS2[j+1]]+iext_s+2*iext_ass+2*extension_cost, c[i][j]);
1940       /**
1941       *** nx2 or more
1942       **/
1943       c[i][j]=MIN2(in[i - 2][j+4]+P->mismatchI[rtype[type]][SS1[i-1]][SS2[j+1]]+iext_s+2*iext_ass+2*extension_cost, c[i][j]);
1944       /**
1945       *** nx1 n>2
1946       **/
1947       c[i][j]=MIN2(inx[i - 3][j+1]+P->mismatch1nI[rtype[type]][SS1[i-1]][SS2[j+1]]+iext_ass+iext_ass+2*extension_cost, c[i][j]);
1948       /**
1949       *** 1xn n>2
1950       **/
1951       c[i][j]=MIN2(iny[i - 1][j+3]+P->mismatch1nI[rtype[type]][SS1[i-1]][SS2[j+1]]+iext_ass+iext_ass+2*extension_cost, c[i][j]);
1952       /**
1953       *** nx0 n>1
1954       **/
1955       int bAU;
1956       bAU=(type>2?P->TerminalAU:0);
1957       c[i][j]=MIN2(bx[i - 2][j+1]+2*extension_cost+bext+bAU, c[i][j]);
1958       /**
1959       *** 0xn n>1
1960       **/
1961       c[i][j]=MIN2(by[i - 1][j+2]+2*extension_cost+bext+bAU, c[i][j]);
1962       temp=min_colonne;
1963       min_colonne=MIN2(c[i][j]+E_ExtLoop(rtype[type], SS2[j-1] , SS1[i+1] , P) + 2*extension_cost, min_colonne);
1964       if(temp>min_colonne){
1965         min_j_colonne=j;
1966       }
1967       /* ---------------------------------------------------------------------end update */
1968     }
1969     if(max>=min_colonne){
1970       max=min_colonne;
1971       max_pos=i;
1972       max_pos_j=min_j_colonne;
1973     }
1974     i++;
1975   }
1976   Emin=max;
1977   i_min=max_pos;
1978   j_min=max_pos_j;
1979   int dGe;
1980   dGe=0;
1981   struc = fbacktrack(i_min, j_min, extension_cost, il_a, il_b, b_a, b_b,&dGe);
1982   if (i_min<n3-10) i_min++;
1983   if (j_min>11 ) j_min--;
1984   l1 = strchr(struc, '&')-struc;
1985   int size;
1986   size=strlen(struc)-1;
1987   Emin-= size * (extension_cost);
1988   mfe.i = i_min;
1989   mfe.j = j_min;
1990   mfe.energy = (double) Emin/100.;
1991   mfe.energy_backtrack = (double) dGe/100.;
1992   mfe.structure = struc;
1993   free(S1); free(S2); free(SS1); free(SS2);
1994   for (i=0; i<=n3; i++) {
1995     free(c[i]);
1996     free(in[i]);
1997     free(bx[i]);
1998     free(by[i]);
1999     free(inx[i]);
2000     free(iny[i]);
2001   }
2002   free(c);free(in);free(bx);free(by);free(inx);free(iny);
2003   return mfe;
2004 }
2005
2006
2007 PRIVATE char *fbacktrack(int i, int j, const int extension_cost,const int il_a, const int il_b, const int b_a, const int b_b, int *dG) {
2008   /* backtrack structure going backwards from i, and forwards from j
2009      return structure in bracket notation with & as separator */
2010   int k, l, type, type2, E, traced, i0, j0;
2011   char *st1, *st2, *struc;
2012   int bopen=b_b;
2013   int bext=b_a+extension_cost;
2014   int iopen=il_b;
2015   int iext_s=2*(il_a+extension_cost);/* iext_s 2 nt nucleotide extension of interior loop, on i and j side */
2016   int iext_ass=50+il_a+extension_cost;/* iext_ass assymetric extension of interior loop, either on i or on j side. */
2017   st1 = (char *) space(sizeof(char)*(n3+1));
2018   st2 = (char *) space(sizeof(char)*(n4+1));
2019   i0=MIN2(i+1,n3-10); j0=MAX2(j-1,11);
2020   int state;
2021   state=1; /* we start backtracking from a a pair , i.e. c-matrix */
2022   /* state 1 -> base pair, c
2023      state 2 -> interior loop, in
2024      state 3 -> bx loop, bx
2025      state 4 -> by loop, by
2026   */
2027   traced=1;
2028   k=i;l=j;
2029   type=pair[S1[i]][S2[j]];
2030   *dG+=E_ExtLoop(rtype[type], SS2[j-1] , SS1[i+1] , P);
2031   /*     (type>2?P->TerminalAU:0)+P->dangle3[rtype[type]][SS1[i+1]]+P->dangle5[rtype[type]][SS2[j-1]]; */
2032   while (i>10 && j<=n4-9 && traced) {
2033     traced=0;
2034     switch(state){
2035     case 1:
2036       type = pair[S1[i]][S2[j]];
2037       int bAU;
2038       bAU=(type>2?P->TerminalAU:0);
2039       if(!type) nrerror("backtrack failed in fold duplex");
2040       type2=pair[S1[i-1]][S2[j+1]];
2041       if(type2 && c[i][j]== (c[i - 1][j+1]+P->stack[rtype[type]][type2]+2*extension_cost)){
2042         k=i-1;
2043         l=j+1;
2044         (*dG)+=E_IntLoop(i-k-1, l-j-1, type2, rtype[type],SS1[k+1], SS2[l-1], SS1[i-1], SS2[j+1],P);
2045         st1[i-1] = '(';
2046         st2[j-1] = ')';
2047         i=k;
2048         j=l;
2049         state=1;
2050         traced=1;
2051         break;
2052       }
2053       type2=pair[S1[i-1]][S2[j+2]];
2054       if(type2 && c[i][j]==(c[i - 1][j+2]+P->bulge[1]+P->stack[rtype[type]][type2]+3*extension_cost)){
2055         k=i-1;
2056         l=j+2;
2057         *dG+=E_IntLoop(i-k-1, l-j-1, type2, rtype[type],SS1[k+1], SS2[l-1], SS1[i-1], SS2[j+1],P);
2058         st1[i-1] = '(';
2059         st2[j-1] = ')';
2060         i=k;
2061         j=l;
2062         state=1;
2063         traced=1;
2064         break;
2065       }
2066       type2=pair[S1[i-2]][S2[j+1]];
2067       if(type2 && c[i][j]==(c[i - 2][j+1]+P->bulge[1]+P->stack[type2][rtype[type]]+3*extension_cost)){
2068         k=i-2;
2069         l=j+1;
2070         *dG+=E_IntLoop(i-k-1, l-j-1, type2, rtype[type],SS1[k+1], SS2[l-1], SS1[i-1], SS2[j+1],P);
2071         st1[i-1] = '(';
2072         st2[j-1] = ')';
2073         i=k;
2074         j=l;
2075         state=1;
2076         traced=1;
2077         break;
2078       }
2079       type2=pair[S1[i-2]][S2[j+2]];
2080       if(type2 && c[i][j]==(c[i - 2][j+2]+P->int11[type2][rtype[type]][SS1[i-1]][SS2[j+1]]+4*extension_cost)){
2081         k=i-2;
2082         l=j+2;
2083         *dG+=E_IntLoop(i-k-1, l-j-1, type2, rtype[type],SS1[k+1], SS2[l-1], SS1[i-1], SS2[j+1],P);
2084         st1[i-1] = '(';
2085         st2[j-1] = ')';
2086         i=k;
2087         j=l;
2088         state=1;
2089         traced=1;
2090         break;
2091       }
2092       type2 = pair[S1[i-3]][S2[j+3]];
2093       if(type2 && c[i][j]==(c[i - 3][j+3]+P->int22[type2][rtype[type]][SS1[i-2]][SS1[i-1]][SS2[j+1]][SS2[j+2]]+6*extension_cost)){
2094         k=i-3;
2095         l=j+3;
2096         *dG+=E_IntLoop(i-k-1, l-j-1, type2, rtype[type],SS1[k+1], SS2[l-1], SS1[i-1], SS2[j+1],P);
2097         st1[i-1] = '(';
2098         st2[j-1] = ')';
2099         i=k;
2100         j=l;
2101         state=1;
2102         traced=1;
2103         break;
2104       }
2105       type2 = pair[S1[i-3]][S2[j+2]];
2106       if(type2 && c[i][j]==(c[i - 3][j+2]+P->int21[rtype[type]][type2][SS2[j+1]][SS1[i-2]][SS1[i-1]]+5*extension_cost)){
2107         k=i-3;
2108         l=j+2;
2109         *dG+=E_IntLoop(i-k-1, l-j-1, type2, rtype[type],SS1[k+1], SS2[l-1], SS1[i-1], SS2[j+1],P);
2110         st1[i-1] = '(';
2111         st2[j-1] = ')';
2112         i=k;
2113         j=l;
2114         state=1;
2115         traced=1;
2116         break;
2117       }
2118       type2 = pair[S1[i-2]][S2[j+3]];
2119       if(type2 && c[i][j]==(c[i - 2][j+3]+P->int21[type2][rtype[type]][SS1[i-1]][SS2[j+1]][SS2[j+2]]+5*extension_cost)){
2120         k=i-2;
2121         l=j+3;
2122         *dG+=E_IntLoop(i-k-1, l-j-1, type2, rtype[type],SS1[k+1], SS2[l-1], SS1[i-1], SS2[j+1],P);
2123         st1[i-1] = '(';
2124         st2[j-1] = ')';
2125         i=k;
2126         j=l;
2127         state=1;
2128         traced=1;
2129         break;
2130       }
2131       type2 = pair[S1[i-4]][S2[j+3]];
2132       if(type2 && c[i][j]==(c[i - 4][j+3]+P->internal_loop[5]+P->ninio[2]+
2133                             P->mismatch23I[type2][SS1[i-3]][SS2[j+2]]+P->mismatch23I[rtype[type]][SS2[j+1]][SS1[i-1]]+7*extension_cost)){
2134         k=i-4;
2135         l=j+3;
2136         *dG+=E_IntLoop(i-k-1, l-j-1, type2, rtype[type],SS1[k+1], SS2[l-1], SS1[i-1], SS2[j+1],P);
2137         st1[i-1] = '(';
2138         st2[j-1] = ')';
2139         i=k;
2140         j=l;
2141         state=1;
2142         traced=1;
2143         break;
2144       }
2145       type2 = pair[S1[i-3]][S2[j+4]];
2146       if(type2 && c[i][j]==(c[i - 3][j+4]+P->internal_loop[5]+P->ninio[2]+
2147                             P->mismatch23I[type2][SS1[i-2]][SS2[j+3]]+P->mismatch23I[rtype[type]][SS2[j+1]][SS1[i-1]]+7*extension_cost)){
2148         k=i-3;
2149         l=j+4;
2150         *dG+=E_IntLoop(i-k-1, l-j-1, type2, rtype[type],SS1[k+1], SS2[l-1], SS1[i-1], SS2[j+1],P);
2151         st1[i-1] = '(';
2152         st2[j-1] = ')';
2153         i=k;
2154         j=l;
2155         state=1;
2156         traced=1;
2157         break;
2158       }
2159       if(c[i][j]==(in[i - 3][j+3]+P->mismatchI[rtype[type]][SS1[i-1]][SS2[j+1]]+2*extension_cost+2*iext_s)){
2160         k=i;
2161         l=j;
2162         st1[i-1] = '(';
2163         st2[j-1] = ')';
2164         i=i-3;
2165         j=j+3;
2166         state=2;
2167         traced=1;
2168         break;
2169       }
2170       if(c[i][j]==(in[i - 4][j+2]+P->mismatchI[rtype[type]][SS1[i-1]][SS2[j+1]]+iext_s+2*iext_ass+2*extension_cost)){
2171         k=i;
2172         l=j;
2173         st1[i-1] = '(';
2174         st2[j-1] = ')';
2175         i=i-4;
2176         j=j+2;
2177         state=2;
2178         traced=1;
2179         break;
2180       }
2181       if(c[i][j]==(in[i - 2][j+4]+P->mismatchI[rtype[type]][SS1[i-1]][SS2[j+1]]+iext_s+2*iext_ass+2*extension_cost)){
2182         k=i;
2183         l=j;
2184         st1[i-1] = '(';
2185         st2[j-1] = ')';
2186         i=i-2;
2187         j=j+4;
2188         state=2;
2189         traced=1;
2190         break;
2191       }
2192       if(c[i][j]==(inx[i - 3][j+1]+P->mismatch1nI[rtype[type]][SS1[i-1]][SS2[j+1]]+iext_ass+iext_ass+2*extension_cost)){
2193         k=i;
2194         l=j;
2195         st1[i-1] = '(';
2196         st2[j-1] = ')';
2197         i=i-3;
2198         j=j+1;
2199         state=5;
2200         traced=1;
2201         break;
2202       }
2203       if(c[i][j]==(iny[i - 1][j+3]+P->mismatch1nI[rtype[type]][SS1[i-1]][SS2[j+1]]+iext_ass+iext_ass+2*extension_cost)){
2204         k=i;
2205         l=j;
2206         st1[i-1] = '(';
2207         st2[j-1] = ')';
2208         i=i-1;
2209         j=j+3;
2210         state=6;
2211         traced=1;
2212         break;
2213       }
2214       if(c[i][j]==(bx[i - 2][j+1]+2*extension_cost+bext+bAU)){
2215         k=i;
2216         l=j;
2217         st1[i-1] = '(';
2218         st2[j-1] = ')';
2219         i=i-2;
2220         j=j+1;
2221         state=3;
2222         traced=1;
2223         break;
2224       }
2225       if(c[i][j]==(by[i - 1][j+2]+2*extension_cost+bext+bAU)){
2226         k=i;
2227         l=j;
2228         st1[i-1] = '(';
2229         st2[j-1] = ')';
2230         i=i-1;
2231         j=j+2;
2232         state=4;
2233         traced=1;
2234         break;
2235       }
2236       break;
2237     case 2:
2238       if(in[i][j]==(in[i - 1][j+1]+iext_s)){
2239         i--;
2240         j++;
2241         state=2;
2242         traced=1;
2243         break;
2244       }
2245       if(in[i][j]==(in[i - 1][j]+iext_ass)){
2246         i=i-1;
2247         state=2;
2248         traced=1;
2249         break;
2250       }
2251       if(in[i][j]==(in[i][j+1]+iext_ass)){
2252         j++;
2253         state=2;
2254         traced=1;
2255         break;
2256       }
2257       type2=pair[S2[j+1]][S1[i-1]];
2258       if(type2 && in[i][j]==(c[i - 1][j+1]+P->mismatchI[type2][SS2[j]][SS1[i]]+iopen+iext_s)){
2259         int temp; temp=k; k=i-1; i=temp;
2260         temp=l; l=j+1; j=temp;
2261         type=pair[S1[i]][S2[j]];
2262         *dG+=E_IntLoop(i-k-1, l-j-1, type2, rtype[type],SS1[k+1], SS2[l-1], SS1[i-1], SS2[j+1],P);
2263         i=k;
2264         j=l;
2265         state=1;
2266         traced=1;
2267         break;
2268       }
2269     case 3:
2270       if(bx[i][j]==(bx[i - 1][j]+bext)){
2271         i--;
2272         state=3;
2273         traced=1;
2274         break;
2275       }
2276       type2=pair[S2[j]][S1[i-1]];
2277       if(type2 && bx[i][j]==(c[i - 1][j]+bopen+bext+(type2>2?P->TerminalAU:0))){
2278         int temp; temp=k; k=i-1; i=temp;
2279         temp=l; l=j; j=temp;
2280         type=pair[S1[i]][S2[j]];
2281         *dG+=E_IntLoop(i-k-1, l-j-1, type2, rtype[type],SS1[k+1], SS2[l-1], SS1[i-1], SS2[j+1],P);
2282         i=k;
2283         j=l;
2284         state=1;
2285         traced=1;
2286         break;
2287       }
2288     case 4:
2289       if(by[i][j]==(by[i][j+1] + bext)){
2290         j++;
2291
2292         state=4;
2293         traced=1;
2294         break;
2295       }
2296       type2=pair[S2[j+1]][S1[i]];
2297       if(type2 && by[i][j]==(c[i][j+1]+bopen+bext+(type2>2?P->TerminalAU:0))){
2298         int temp; temp=k; k=i; i=temp;
2299         temp=l; l=j+1; j=temp;
2300         type=pair[S1[i]][S2[j]];
2301         *dG+=E_IntLoop(i-k-1, l-j-1, type2, rtype[type],SS1[k+1], SS2[l-1], SS1[i-1], SS2[j+1],P);
2302         i=k;
2303         j=l;
2304         state=1;
2305         traced=1;
2306         break;
2307       }
2308     case 5:
2309       if(inx[i][j]==(inx[i-1][j]+iext_ass)) {
2310         i--;
2311         state=5;
2312         traced=1;
2313         break;
2314       }
2315       type2=pair[S2[j+1]][S1[i-1]];
2316       if(type2 && inx[i][j]==(c[i-1][j+1]+P->mismatch1nI[type2][SS2[j]][SS1[i]]+iopen+iext_s)){
2317         int temp; temp=k; k=i-1; i=temp;
2318         temp=l; l=j+1; j=temp;
2319         type=pair[S1[i]][S2[j]];
2320         *dG+=E_IntLoop(i-k-1, l-j-1, type2, rtype[type],SS1[k+1], SS2[l-1], SS1[i-1], SS2[j+1],P);
2321         i=k;
2322         j=l;
2323         state=1;
2324         traced=1;
2325         break;
2326       }
2327     case 6:
2328       if(iny[i][j]==(iny[i][j+1]+iext_ass)) {
2329         j++;
2330         state=6;
2331         traced=1;
2332         break;
2333       }
2334       type2=pair[S2[j+1]][S1[i-1]];
2335       if(type2 && iny[i][j]==(c[i-1][j+1]+P->mismatch1nI[type2][SS2[j]][SS1[i]]+iopen+iext_s)){
2336         int temp; temp=k; k=i-1; i=temp;
2337         temp=l; l=j+1; j=temp;
2338         type=pair[S1[i]][S2[j]];
2339         *dG+=E_IntLoop(i-k-1, l-j-1, type2, rtype[type],SS1[k+1], SS2[l-1], SS1[i-1], SS2[j+1],P);
2340         i=k;
2341         j=l;
2342         state=1;
2343         traced=1;
2344         break;
2345       }
2346     }
2347   }
2348   if (!traced) {
2349     E=c[i][j];
2350     /**
2351     ***    if (i>1) {E -= P->dangle5[type][SS1[i-1]]+extension_cost; *dG+=P->dangle5[type][SS1[i-1]];}
2352     ***    if (j<n4){E -= P->dangle3[type][SS2[j+1]]+extension_cost; *dG+=P->dangle3[type][SS2[j+1]];}
2353     ***    if (type>2) {E -= P->TerminalAU; *dG+=P->TerminalAU;}
2354     **/
2355     int correction;
2356     correction = E_ExtLoop(type, (i>1) ? SS1[i-1] : -1, (j<n4) ? SS2[j+1] : -1, P);
2357     *dG+=correction;
2358     E-=correction+2*extension_cost;
2359     if (E != P->DuplexInit+2*extension_cost) {
2360       nrerror("backtrack failed in second fold duplex");
2361     }
2362     else{
2363       *dG+=P->DuplexInit;
2364       st1[i-1]='(';
2365       st2[j-1]=')';
2366     }
2367   }
2368   if (i>11)  i--;
2369   if (j<n4-10) j++;
2370   struc = (char *) space(i0-i+1+j-j0+1+2);
2371   for (k=MAX2(i,1); k<=i0; k++) if (!st1[k-1]) st1[k-1] = '.';
2372   for (k=j0; k<=j; k++) if (!st2[k-1]) st2[k-1] = '.';
2373   strcpy(struc, st1+MAX2(i-1,0));
2374   strcat(struc, "&");
2375   strcat(struc, st2+j0-1);
2376   /* printf("%s %3d,%-3d : %3d,%-3d\n", struc, i,i0,j0,j);  */
2377   free(st1); free(st2);
2378   return struc;
2379 }
2380
2381
2382 duplexT ** Lduplexfold(const char *s1, const char *s2, const int threshold, const int extension_cost, const int alignment_length, const int delta, const int fast, const int il_a, const int il_b, const int b_a, const int b_b)
2383 {
2384   /**
2385   *** See variable definition in fduplexfold_XS
2386   **/
2387   int i, j;
2388   int bopen=b_b;
2389   int bext=b_a+extension_cost;
2390   int iopen=il_b;
2391   int iext_s=2*(il_a+extension_cost);/* iext_s 2 nt nucleotide extension of interior loop, on i and j side */
2392   int iext_ass=50+il_a+extension_cost;/* iext_ass assymetric extension of interior loop, either on i or on j side. */
2393   int min_colonne=INF; /* enthaelt das maximum einer kolonne */
2394   int i_length;
2395   int max_pos;/* get position of the best hit */
2396   int max_pos_j;
2397   int temp=INF;
2398   int min_j_colonne;
2399   int max=INF;
2400   int *position; /* contains the position of the hits with energy > E */
2401   int *position_j;
2402   /**
2403   *** 1D array corresponding to the standard 2d recursion matrix
2404   *** Makes the computation 20% faster
2405   **/
2406   int *SA;
2407   /**
2408   *** variable initialization
2409   **/
2410   n1 = (int) strlen(s1);
2411   n2 = (int) strlen(s2);
2412   /**
2413   *** Sequence encoding
2414   **/
2415   if ((!P) || (fabs(P->temperature - temperature)>1e-6)) {
2416     update_fold_params();  if(P) free(P); P = scale_parameters();
2417     make_pair_matrix();
2418   }
2419   encode_seqs(s1,s2);
2420   /**
2421   *** Position of the high score on the target and query sequence
2422   **/
2423   position = (int *) space ((delta+n1+3+delta) * sizeof(int));
2424   position_j= (int *) space ((delta+n1+3+delta) * sizeof(int));
2425   /**
2426   *** instead of having 4 2-dim arrays we use a unique 1-dim array
2427   *** The mapping 2d -> 1D is done based ont the macro
2428   *** LCI(i,j,l)      ((i     )*l + j)
2429   *** LINI(i,j,l)     ((i +  5)*l + j)
2430   *** LBXI(i,j,l)     ((i + 10)*l + j)
2431   *** LBYI(i,j,l)     ((i + 15)*l + j)
2432   *** LINIX(i,j,l)    ((i + 20)*l + j)
2433   *** LINIY(i,j,l)    ((i + 25)*l + j)
2434   ***
2435   *** SA has a length of 5 (number of columns we look back) *
2436   ***                  * 6 (number of structures we look at) *
2437   ***                  * length of the sequence
2438   **/
2439   SA=(int *) space(sizeof(int)*5*6*(n2+5));
2440   for(j=n2+4;j>=0;j--) {
2441     SA[(j*30)   ]=SA[(j*30)+1   ]=SA[(j*30)+2   ]=SA[(j*30)+3   ]=SA[(j*30)+4   ]=INF;
2442     SA[(j*30)+5 ]=SA[(j*30)+1+5 ]=SA[(j*30)+2+5 ]=SA[(j*30)+3+5 ]=SA[(j*30)+4+5 ]=INF;
2443     SA[(j*30)+10]=SA[(j*30)+1+10]=SA[(j*30)+2+10]=SA[(j*30)+3+10]=SA[(j*30)+4+10]=INF;
2444     SA[(j*30)+15]=SA[(j*30)+1+15]=SA[(j*30)+2+15]=SA[(j*30)+3+15]=SA[(j*30)+4+15]=INF;
2445     SA[(j*30)+20]=SA[(j*30)+1+20]=SA[(j*30)+2+20]=SA[(j*30)+3+20]=SA[(j*30)+4+20]=INF;
2446     SA[(j*30)+25]=SA[(j*30)+1+25]=SA[(j*30)+2+25]=SA[(j*30)+3+25]=SA[(j*30)+4+25]=INF;
2447   }
2448   i=10;
2449   i_length= n1 - 9 ;
2450   while(i < i_length) {
2451     int idx=i%5;
2452     int idx_1=(i-1)%5;
2453     int idx_2=(i-2)%5;
2454     int idx_3=(i-3)%5;
2455     int idx_4=(i-4)%5;
2456     j=n2-9;
2457     while (9 < --j) {
2458       int type, type2;
2459       type = pair[S1[i]][S2[j]];
2460       /**
2461       *** Start duplex
2462       **/
2463       SA[LCI(idx,j,n2)]=type ? P->DuplexInit + 2*extension_cost : INF;
2464       /**
2465       *** update lin bx by linx liny matrix
2466       **/
2467       type2=pair[S2[j+1]][S1[i-1]];
2468       /**
2469       *** start/extend interior loop
2470       **/
2471       SA[LINI(idx,j,n2)]=MIN2(SA[LCI(idx_1,j+1,n2)]+P->mismatchI[type2][SS2[j]][SS1[i]]+iopen+iext_s,
2472                               SA[LINI(idx_1,j,n2)]+iext_ass);
2473       /**
2474       *** start/extend nx1 target
2475       *** use same type2 as for in
2476       **/
2477       SA[LINIX(idx,j,n2)]=MIN2(SA[LCI(idx_1,j+1,n2)]+P->mismatch1nI[type2][SS2[j]][SS1[i]]+iopen+iext_s,
2478                                SA[LINIX(idx_1,j,n2)]+iext_ass);
2479       /**
2480       *** start/extend 1xn target
2481       *** use same type2 as for in
2482       **/
2483       SA[LINIY(idx,j,n2)]=MIN2(SA[LCI(idx_1,j+1,n2)]+P->mismatch1nI[type2][SS2[j]][SS1[i]]+iopen+iext_s,
2484                                SA[LINIY(idx,j+1,n2)]+iext_ass);
2485       /**
2486       *** extend interior loop
2487       **/
2488       SA[LINI(idx,j,n2)]=MIN2(SA[LINI(idx,j,n2)],SA[LINI(idx,j+1,n2)]+iext_ass);
2489       SA[LINI(idx,j,n2)]=MIN2(SA[LINI(idx,j,n2)],SA[LINI(idx_1,j+1,n2)]+iext_s);
2490       /**
2491       *** start/extend bulge target
2492       **/
2493       type2=pair[S2[j]][S1[i-1]];
2494       SA[LBXI(idx,j,n2)]=MIN2(SA[LBXI(idx_1,j,n2)]+bext, SA[LCI(idx_1,j,n2)]+bopen+bext+(type2>2?P->TerminalAU:0));
2495       /**
2496       *** start/extend bulge query
2497       **/
2498       type2=pair[S2[j+1]][S1[i]];
2499       SA[LBYI(idx,j,n2)]=MIN2(SA[LBYI(idx,j+1,n2)]+bext, SA[LCI(idx,j+1,n2)]+bopen+bext+(type2>2?P->TerminalAU:0));
2500       /**
2501       ***end update recursion
2502       ***##################### Start stack extension ######################
2503       **/
2504       if(!type){continue;}
2505       /**
2506       *** stack extension
2507       **/
2508       SA[LCI(idx,j,n2)]+= E_ExtLoop(type, SS1[i-1] , SS2[j+1], P) + 2*extension_cost;
2509       /**
2510       *** stack extension
2511       **/
2512       if((type2=pair[S1[i-1]][S2[j+1]]))
2513         SA[LCI(idx,j,n2)]=MIN2(SA[LCI(idx_1,j+1,n2)]+P->stack[rtype[type]][type2]+2*extension_cost, SA[LCI(idx,j,n2)]);
2514       /**
2515       *** 1x0 / 0x1 stack extension
2516       **/
2517       if((type2=pair[S1[i-1]][S2[j+2]]))
2518         SA[LCI(idx,j,n2)]=MIN2(SA[LCI(idx_1,j+2,n2)]+P->bulge[1]+P->stack[rtype[type]][type2]+3*extension_cost,SA[LCI(idx,j,n2)]);
2519       if((type2=pair[S1[i-2]][S2[j+1]]))
2520         SA[LCI(idx,j,n2)]=MIN2(SA[LCI(idx_2,j+1,n2)]+P->bulge[1]+P->stack[type2][rtype[type]]+3*extension_cost,SA[LCI(idx,j,n2)]);
2521       /**
2522       *** 1x1 / 2x2 stack extension
2523       **/
2524       if((type2=pair[S1[i-2]][S2[j+2]]))
2525         SA[LCI(idx,j,n2)]=MIN2(SA[LCI(idx_2,j+2,n2)]+P->int11[type2][rtype[type]][SS1[i-1]][SS2[j+1]]+4*extension_cost, SA[LCI(idx,j,n2)]);
2526       if((type2 = pair[S1[i-3]][S2[j+3]]))
2527         SA[LCI(idx,j,n2)]=MIN2(SA[LCI(idx_3,j+3,n2)]+P->int22[type2][rtype[type]][SS1[i-2]][SS1[i-1]][SS2[j+1]][SS2[j+2]]+6*extension_cost,SA[LCI(idx,j,n2)]);
2528       /**
2529       *** 1x2 / 2x1 stack extension
2530       *** E_IntLoop(1,2,type2, rtype[type],SS1[i-1], SS2[j+2], SS1[i-1], SS2[j+1], P) corresponds to
2531       *** P->int21[rtype[type]][type2][SS2[j+2]][SS1[i-1]][SS1[i-1]]
2532       **/
2533       if((type2 = pair[S1[i-3]][S2[j+2]]))
2534         SA[LCI(idx,j,n2)]=MIN2(SA[LCI(idx_3,j+2,n2)]+P->int21[rtype[type]][type2][SS2[j+1]][SS1[i-2]][SS1[i-1]]+5*extension_cost, SA[LCI(idx,j,n2)]);
2535       if((type2 = pair[S1[i-2]][S2[j+3]]))
2536         SA[LCI(idx,j,n2)]=MIN2(SA[LCI(idx_2,j+3,n2)]+P->int21[type2][rtype[type]][SS1[i-1]][SS2[j+1]][SS2[j+2]]+5*extension_cost, SA[LCI(idx,j,n2)]);
2537       /**
2538       *** 2x3 / 3x2 stack extension
2539       **/
2540       if((type2 = pair[S1[i-4]][S2[j+3]]))
2541         SA[LCI(idx,j,n2)]=MIN2(SA[LCI(idx_4,j+3,n2)]+P->internal_loop[5]+P->ninio[2]+
2542                                P->mismatch23I[type2][SS1[i-3]][SS2[j+2]]+P->mismatch23I[rtype[type]][SS2[j+1]][SS1[i-1]]+7*extension_cost, SA[LCI(idx,j,n2)]);
2543       if((type2 = pair[S1[i-3]][S2[j+4]]))
2544         SA[LCI(idx,j,n2)]=MIN2(SA[LCI(idx_3,j+4,n2)]+P->internal_loop[5]+P->ninio[2]+
2545                                P->mismatch23I[type2][SS1[i-2]][SS2[j+3]]+P->mismatch23I[rtype[type]][SS2[j+1]][SS1[i-1]]+7*extension_cost, SA[LCI(idx,j,n2)]);
2546       /**
2547       *** So now we have to handle 1x3, 3x1, 3x3, and mxn m,n > 3
2548       **/
2549       /**
2550       *** 3x3 or more
2551       **/
2552       SA[LCI(idx,j,n2)]=MIN2(SA[LINI(idx_3,j+3,n2)]+P->mismatchI[rtype[type]][SS1[i-1]][SS2[j+1]]+2*iext_s+2*extension_cost,SA[LCI(idx,j,n2)]);
2553       /**
2554       *** 2xn or more
2555       **/
2556       SA[LCI(idx,j,n2)]=MIN2(SA[LINI(idx_4,j+2,n2)]+P->mismatchI[rtype[type]][SS1[i-1]][SS2[j+1]]+iext_s+2*iext_ass+2*extension_cost, SA[LCI(idx,j,n2)]);
2557       /**
2558       *** nx2 or more
2559       **/
2560       SA[LCI(idx,j,n2)]=MIN2(SA[LINI(idx_2,j+4,n2)]+P->mismatchI[rtype[type]][SS1[i-1]][SS2[j+1]]+iext_s+2*iext_ass+2*extension_cost, SA[LCI(idx,j,n2)]);
2561       /**
2562       *** nx1 n>2
2563       **/
2564       SA[LCI(idx,j,n2)]=MIN2(SA[LINIX(idx_3,j+1,n2)]+P->mismatch1nI[rtype[type]][SS1[i-1]][SS2[j+1]]+iext_ass+iext_ass+2*extension_cost, SA[LCI(idx,j,n2)]);
2565       /**
2566       *** 1xn n>2
2567       **/
2568       SA[LCI(idx,j,n2)]=MIN2(SA[LINIY(idx_1,j+3,n2)]+P->mismatch1nI[rtype[type]][SS1[i-1]][SS2[j+1]]+iext_ass+iext_ass+2*extension_cost, SA[LCI(idx,j,n2)]);
2569       /**
2570       *** nx0 n>1
2571       **/
2572       int bAU;
2573       bAU=(type>2?P->TerminalAU:0);
2574       SA[LCI(idx,j,n2)]=MIN2(SA[LBXI(idx_2,j+1,n2)]+2*extension_cost+bext+bAU,SA[LCI(idx,j,n2)]);
2575       /**
2576       *** 0xn n>1
2577       **/
2578       SA[LCI(idx,j,n2)]=MIN2(SA[LBYI(idx_1,j+2,n2)]+2*extension_cost+bext+bAU,SA[LCI(idx,j,n2)]);
2579       temp=min_colonne;
2580
2581       min_colonne=MIN2(SA[LCI(idx,j,n2)]+E_ExtLoop(rtype[type], SS2[j-1] , SS1[i+1] , P) + 2*extension_cost, min_colonne);
2582       if(temp>min_colonne){
2583         min_j_colonne=j;
2584       }
2585     }
2586     if(max>=min_colonne){
2587       max=min_colonne;
2588       max_pos=i;
2589       max_pos_j=min_j_colonne;
2590     }
2591     position[i+delta]=min_colonne;min_colonne=INF;
2592     position_j[i+delta]=min_j_colonne;
2593     i++;
2594   }
2595   /* printf("MAX: %d",max); */
2596   free(S1); free(S2); free(SS1); free(SS2);
2597   if(max<threshold){
2598     find_max(position, position_j, delta, threshold, alignment_length, s1, s2, extension_cost, fast, il_a, il_b, b_a, b_b);
2599   }
2600   if(max<INF){
2601     plot_max(max, max_pos, max_pos_j,alignment_length, s1, s2, extension_cost,fast, il_a, il_b, b_a, b_b);
2602   }
2603   free(SA);
2604   free(position);
2605   free(position_j);
2606   return NULL;
2607 }
2608
2609
2610
2611
2612 PRIVATE void find_max(const int *position, const int *position_j,const int delta, const int threshold, const int alignment_length, const char *s1, const char *s2, const int extension_cost, const int fast,const int il_a, const int il_b, const int b_a, const int b_b){
2613   int pos=n1-9;
2614   if(fast==1){
2615     while(10 < pos--){
2616       int temp_min=0;
2617       if(position[pos+delta]<(threshold)){
2618         int search_range;
2619         search_range=delta+1;
2620         while(--search_range){
2621           if(position[pos+delta-search_range]<=position[pos+delta-temp_min]){
2622             temp_min=search_range;
2623           }
2624         }
2625         pos-=temp_min;
2626         int max_pos_j;
2627         max_pos_j=position_j[pos+delta];
2628         int max;
2629         max=position[pos+delta];
2630         printf("target upper bound %d: query lower bound %d  (%5.2f) \n", pos-10, max_pos_j-10, ((double)max)/100);
2631         pos=MAX2(10,pos-alignment_length);
2632       }
2633     }
2634   }
2635   else if(fast==2){
2636     pos=n1-9;
2637     while(10 < pos--){
2638       int temp_min=0;
2639       if(position[pos+delta]<(threshold)){
2640         int search_range;
2641         search_range=delta+1;
2642         while(--search_range){
2643           if(position[pos+delta-search_range]<=position[pos+delta-temp_min]){
2644             temp_min=search_range;
2645           }
2646         }
2647         pos-=temp_min;
2648         int max_pos_j;
2649         max_pos_j=position_j[pos+delta];
2650         /* max_pos_j und pos entsprechen die realen position
2651            in der erweiterten sequenz.
2652            pos=1 -> position 1 in the sequence (and not 0 like in C)
2653            max_pos_j -> position 1 in the sequence ( not 0 like in C)
2654         */
2655         int begin_t=MAX2(11, pos-alignment_length+1);/* 10 */
2656         int end_t  =MIN2(n1-10, pos+1);
2657         int begin_q=MAX2(11, max_pos_j-1); /* 10 */
2658         int end_q  =MIN2(n2-10, max_pos_j+alignment_length-1);
2659         char *s3 = (char*) space(sizeof(char)*(end_t - begin_t +2 + 20));
2660         char *s4 = (char*) space(sizeof(char)*(end_q - begin_q +2 + 20));
2661         strcpy(s3,"NNNNNNNNNN");strcpy(s4,"NNNNNNNNNN");
2662         strncat(s3, (s1+begin_t-1),  end_t - begin_t +1);
2663         strncat(s4, (s2+begin_q-1) , end_q - begin_q +1);
2664         strcat(s3,"NNNNNNNNNN");strcat(s4,"NNNNNNNNNN");
2665         s3[end_t -begin_t +1 +20 ]='\0';
2666         s4[end_q -begin_q +1 +20]='\0';
2667         duplexT test;
2668         test = fduplexfold(s3, s4, extension_cost,il_a, il_b, b_a, b_b);
2669         if(test.energy * 100 < threshold){
2670           int l1=strchr(test.structure, '&')-test.structure;
2671           printf("%s %3d,%-3d : %3d,%-3d (%5.2f) [%5.2f]\n", test.structure,
2672                  begin_t-10+test.i-l1-10,
2673                  begin_t-10+test.i-1-10,
2674                  begin_q-10 + test.j-1-10 ,
2675                  (begin_q -11) + test.j + strlen(test.structure)-l1-2-10,
2676                  test.energy,test.energy_backtrack);
2677           pos=MAX2(10,pos-delta);
2678         }
2679         free(s3);free(s4);
2680         free(test.structure);
2681       }
2682     }
2683   }
2684   else{
2685     pos=n1-9;
2686     while(10 < pos--){
2687       int temp_min=0;
2688       if(position[pos+delta]<(threshold)){
2689         int search_range;
2690         search_range=delta+1;
2691         while(--search_range){
2692           if(position[pos+delta-search_range]<=position[pos+delta-temp_min]){
2693             temp_min=search_range;
2694           }
2695         }
2696         pos-=temp_min;
2697         int max_pos_j;
2698         max_pos_j=position_j[pos+delta];
2699         /* max_pos_j und pos entsprechen die realen position
2700            in der erweiterten sequenz.
2701            pos=1 -> position 1 in the sequence (and not 0 like in C)
2702            max_pos_j -> position 1 in the sequence ( not 0 like in C)
2703         */
2704         int begin_t=MAX2(11, pos-alignment_length+1);/* 10 */
2705         int end_t  =MIN2(n1-10, pos+1);
2706         int begin_q=MAX2(11, max_pos_j-1); /* 10 */
2707         int end_q  =MIN2(n2-10, max_pos_j+alignment_length-1);
2708         char *s3 = (char*) space(sizeof(char)*(end_t - begin_t +2));
2709         char *s4 = (char*) space(sizeof(char)*(end_q - begin_q +2));
2710         strncpy(s3, (s1+begin_t-1),  end_t - begin_t +1);
2711         strncpy(s4, (s2+begin_q-1) , end_q - begin_q +1);
2712         s3[end_t -begin_t +1 ]='\0';
2713         s4[end_q -begin_q +1 ]='\0';
2714         duplexT test;
2715         test = duplexfold(s3, s4, extension_cost);
2716         if(test.energy * 100 < threshold){
2717           int l1=strchr(test.structure, '&')-test.structure;
2718           printf("%s %3d,%-3d : %3d,%-3d (%5.2f)\n", test.structure,
2719                  begin_t-10+test.i-l1,
2720                  begin_t-10+test.i-1,
2721                  begin_q-10 + test.j-1 ,
2722                  (begin_q -11) + test.j + strlen(test.structure)-l1-2,
2723                  test.energy);
2724           pos=MAX2(10,pos-delta);
2725         }
2726         free(s3);free(s4);
2727         free(test.structure);
2728       }
2729     }
2730   }
2731 }
2732 PRIVATE void plot_max(const int max, const int max_pos, const int max_pos_j, const int alignment_length, const char *s1, const char *s2, const int extension_cost, const int fast,const int il_a, const int il_b, const int b_a, const int b_b)
2733 {
2734   if(fast==1){
2735     printf("target upper bound %d: query lower bound %d (%5.2f)\n", max_pos-10, max_pos_j-10, ((double)max)/100);
2736   }
2737   else if(fast==2){
2738     int begin_t=MAX2(11, max_pos-alignment_length+1);/* 10 */
2739     int end_t  =MIN2(n1-10, max_pos+1);
2740     int begin_q=MAX2(11, max_pos_j-1); /* 10 */
2741     int end_q  =MIN2(n2-10, max_pos_j+alignment_length-1);
2742     char *s3 = (char*) space(sizeof(char)*(end_t - begin_t +2 + 20));
2743     char *s4 = (char*) space(sizeof(char)*(end_q - begin_q +2 + 20));
2744     strcpy(s3,"NNNNNNNNNN");strcpy(s4,"NNNNNNNNNN");
2745     strncat(s3, (s1+begin_t-1),  end_t - begin_t +1);
2746     strncat(s4, (s2+begin_q-1) , end_q - begin_q +1);
2747     strcat(s3,"NNNNNNNNNN");strcat(s4,"NNNNNNNNNN");
2748     s3[end_t -begin_t +1 +20 ]='\0';
2749     s4[end_q -begin_q +1 +20]='\0';
2750     duplexT test;
2751     test = fduplexfold(s3, s4, extension_cost,il_a, il_b, b_a, b_b);
2752     int l1=strchr(test.structure, '&')-test.structure;
2753     printf("%s %3d,%-3d : %3d,%-3d (%5.2f) [%5.2f]\n", test.structure,
2754            begin_t-10+test.i-l1-10,
2755            begin_t-10+test.i-1-10,
2756            begin_q-10 + test.j-1-10 ,
2757            (begin_q -11) + test.j + strlen(test.structure)-l1-2-10,
2758            test.energy, test.energy_backtrack);
2759     free(s3);free(s4);free(test.structure);
2760   }
2761   else{
2762     duplexT test;
2763     int begin_t=MAX2(11, max_pos-alignment_length+1);
2764     int end_t  =MIN2(n1-10, max_pos+1);
2765     int begin_q=MAX2(11, max_pos_j-1);
2766     int end_q  =MIN2(n2-10, max_pos_j+alignment_length-1);
2767     char *s3 = (char*) space(sizeof(char)*(end_t - begin_t +2));
2768     char *s4 = (char*) space(sizeof(char)*(end_q - begin_q +2));
2769     strncpy(s3, (s1+begin_t-1),  end_t - begin_t + 1);
2770     strncpy(s4, (s2+begin_q-1) , end_q - begin_q +1 );
2771     s3[end_t -begin_t +1 ]='\0';
2772     s4[end_q -begin_q +1 ]='\0';
2773     test = duplexfold(s3, s4, extension_cost);
2774     int l1=strchr(test.structure, '&')-test.structure;
2775     printf("%s %3d,%-3d : %3d,%-3d (%5.2f)\n", test.structure,
2776            begin_t-10+test.i-l1,
2777            begin_t-10+test.i-1,
2778            begin_q-10 +test.j-1 ,
2779            (begin_q -11) + test.j + strlen(test.structure)-l1-2,
2780            test.energy);
2781     free(s3);free(s4);free(test.structure);
2782   }
2783 }
2784
2785
2786 PRIVATE void update_dfold_params(void)
2787 {
2788   if(P) free(P);
2789   P = scale_parameters();
2790   make_pair_matrix();
2791 }
2792
2793 PRIVATE void encode_seqs(const char *s1, const char *s2) {
2794   unsigned int i,l;
2795
2796   l = strlen(s1);
2797   S1 = encode_seq(s1);
2798   SS1= (short *) space(sizeof(short)*(l+1));
2799   /* SS1 exists only for the special X K and I bases and energy_set!=0 */
2800
2801   for (i=1; i<=l; i++) { /* make numerical encoding of sequence */
2802     SS1[i] = alias[S1[i]];   /* for mismatches of nostandard bases */
2803   }
2804
2805   l = strlen(s2);
2806   S2 = encode_seq(s2);
2807   SS2= (short *) space(sizeof(short)*(l+1));
2808   /* SS2 exists only for the special X K and I bases and energy_set!=0 */
2809
2810   for (i=1; i<=l; i++) { /* make numerical encoding of sequence */
2811     SS2[i] = alias[S2[i]];   /* for mismatches of nostandard bases */
2812   }
2813 }
2814
2815 PRIVATE short * encode_seq(const char *sequence) {
2816   unsigned int i,l;
2817   short *S;
2818   l = strlen(sequence);
2819   S = (short *) space(sizeof(short)*(l+2));
2820   S[0] = (short) l;
2821
2822   /* make numerical encoding of sequence */
2823   for (i=1; i<=l; i++)
2824     S[i]= (short) encode_char(toupper(sequence[i-1]));
2825
2826   /* for circular folding add first base at position n+1 */
2827   S[l+1] = S[1];
2828
2829   return S;
2830 }
2831
2832 int arraySize(duplexT** array)
2833 {
2834   int site_count=0;
2835   while(array[site_count]!=NULL){
2836     site_count++;
2837   }
2838   return site_count;
2839 }
2840
2841 void freeDuplexT(duplexT** array)
2842 {
2843   int size=arraySize(array);
2844   while(--size){
2845     free(array[size]->structure);
2846     free(array[size]);
2847   }
2848   free(array[0]->structure);
2849   free(array);
2850 }
2851
2852
2853