binaries/src/ViennaRNA/lib/snoop.c

   1
   2
   3 /* Last changed Time-stamp: <2007-08-26 11:59:45 ivo> */
   4 /*
   5            compute the duplex structure of two RNA strands,
   6                 allowing only inter-strand base pairs.
   7          see cofold() for computing hybrid structures without
   8                              restriction.
   9
  10                              Ivo Hofacker
  11                           Vienna RNA package
  12 */
  13
  14 #include <config.h>
  15 #include <stdio.h>
  16 #include <stdlib.h>
  17 #include <math.h>
  18 #include <ctype.h>
  19 #include <string.h>
  20 #include "utils.h"
  21 #include "energy_par.h"
  22 #include "fold_vars.h"
  23 #include "snofold.h"
  24 #include "pair_mat.h"
  25 #include "params.h"
  26 #include "snoop.h"
  27 #include "PS_dot.h"
  28 /* #include "fold.h" */
  29 #include "duplex.h"
  30 #include "loop_energies.h"
  31 /*@unused@*/
  32 static char rcsid[] UNUSED = "$Id: duplex.c,v 1.8 2007/08/26 10:08:44 ivo Exp $";
  33 #define UNIT 100
  34 #define MINPSCORE -2*UNIT
  35 #define PUBLIC
  36 #define PRIVATE static
  37
  38 #define STACK_BULGE1  1   /* stacking energies for bulges of size 1 */
  39 #define NEW_NINIO     1   /* new asymetry penalty */
  40
  41
  42
  43 PRIVATE void  encode_seqs(const char *s1, const char *s2);
  44 PRIVATE short *encode_seq(const char *seq);
  45
  46 PRIVATE void find_max_snoop(const char *s1, const char *s2, const int max,
  47                             const int alignment_length, const int* position,
  48                             const int delta, const int distance,  const int penalty,
  49                             const int threshloop, const int threshLE, const int threshRE,
  50                             const int threshDE, const int threshTE, const int threshSE, const int threshD,
  51                             const int half_stem, const int max_half_stem,  const int min_s2,
  52                             const int max_s2, const int min_s1, const int  max_s1, const int min_d1, const int min_d2, const char* name, const int fullStemEnergy);
  53
  54 PRIVATE void find_max_snoop_XS(const char *s1, const char *s2, const int **access_s1, const int max,
  55                                const int alignment_length, const int* position, const int *position_j,
  56                                const int delta, const int distance,  const int penalty,
  57                                const int threshloop, const int threshLE, const int threshRE,
  58                                const int threshDE, const int threshTE, const int threshSE, const int threshD,
  59                                const int half_stem, const int max_half_stem,  const int min_s2,
  60                                const int max_s2, const int min_s1, const int  max_s1, const int min_d1, const int min_d2, const char *name, const int fullStemEnergy);
  61
  62
  63
  64
  65
  66 PRIVATE char * alisnoop_backtrack(int i, int j, const char ** s2,
  67                                   int* Duplex_El, int* Duplex_Er, int* Loop_E, int *Loop_D, int *u,
  68                                   int *pscd, int *psct, int *pscg,
  69                                   const int penalty, const int threshloop,
  70                                   const int threshLE, const int threshRE, const int threshDE, const int threshD,
  71                                   const int half_stem, const int max_half_stem,
  72                                   const int min_s2, const int max_s2, const int min_s1,
  73                                   const int max_s1, const int min_d1, const int min_d2,
  74                                   const short **S1, const short **S2);
  75
  76
  77 PRIVATE char * snoop_backtrack(int i, int j, const char* s2, int* Duplex_El, int* Duplex_Er, int* Loop_E, int *Loop_D, int *u,
  78                                const int penalty, const int threshloop, const int threshLE, const int threshRE, const int threshDE,
  79                                const int threshD,
  80                                const int half_stem, const int max_half_stem,
  81                                const int min_s2, const int max_s2, const int min_s1, const int max_s1, const int min_d1, const int min_d2);
  82
  83 PRIVATE char * snoop_backtrack_XS(int i, int j, const char* s2, int* Duplex_El, int* Duplex_Er, int* Loop_E, int *Loop_D, int *u,
  84                                const int penalty, const int threshloop, const int threshLE, const int threshRE, const int threshDE,
  85                                const int threshD,
  86                                const int half_stem, const int max_half_stem,
  87                                const int min_s2, const int max_s2, const int min_s1, const int max_s1, const int min_d1, const int min_d2);
  88
  89
  90
  91
  92 PRIVATE int compare(const void *sub1, const void *sub2);
  93 PRIVATE int covscore(const int *types, int n_seq);
  94 PRIVATE short * aliencode_seq(const char *sequence);
  95
  96 PUBLIC  int snoop_subopt_sorted=0; /* from subopt.c, default 0 */
  97
  98
  99 /*@unused@*/
 100
 101
 102
 103
 104 #define MAXLOOP_L        3
 105 #define MIN2(A, B)      ((A) < (B) ? (A) : (B))
 106 #define MAX2(A, B)      ((A) > (B) ? (A) : (B))
 107 #define ASS                1
 108 PRIVATE paramT *P = NULL;
 109
 110 PRIVATE int   **c = NULL;      /* energy array, given that i-j pair */
 111 PRIVATE int   **r = NULL;
 112 PRIVATE int   **lc = NULL;      /* energy array, given that i-j pair */
 113 PRIVATE int   **lr = NULL;
 114 PRIVATE int   **c_fill = NULL;
 115 PRIVATE int   **r_fill = NULL;
 116 PRIVATE int   **lpair = NULL;
 117
 118
 119 PRIVATE short  *S1 = NULL, *SS1 = NULL, *S2 = NULL, *SS2 = NULL;
 120 PRIVATE short *S1_fill = NULL, *SS1_fill = NULL, *S2_fill = NULL, *SS2_fill = NULL;
 121 PRIVATE int   n1,n2;    /* sequence lengths */
 122
 123 extern int cut_point;
 124
 125 PRIVATE int delay_free=0;
 126 /*--------------------------------------------------------------------------*/
 127
 128 snoopT alisnoopfold(const char **s1, const char **s2,
 129                     const int penalty, const int threshloop,
 130                     const int threshLE, const int threshRE, const int threshDE, const int threshD,
 131                     const int half_stem, const int max_half_stem,
 132                     const int min_s2, const int max_s2, const int min_s1,
 133                     const int max_s1, const int min_d1, const int min_d2) {
 134
 135   int s,n_seq;
 136   int i, j, E, l1,Emin=INF, i_min=0, j_min=0;
 137   char *struc;
 138   snoopT mfe;
 139   int *indx;
 140   int *mLoop;
 141   int *cLoop;
 142   folden  **foldlist; folden **foldlist_XS;
 143   int Duplex_El, Duplex_Er,pscd,psct,pscg;
 144   int Loop_D;
 145   int u;
 146   int Loop_E;
 147   short **Sali1,**Sali2;
 148   int *type,*type2,*type3;
 149   Duplex_El=0;Duplex_Er=0;Loop_E=0; Loop_D=0;pscd=0;psct=0;pscg=0;
 150   snoexport_fold_arrays(&indx, &mLoop, &cLoop,&foldlist, &foldlist_XS);
 151   n1 = (int) strlen(s1[0]);
 152   n2 = (int) strlen(s2[0]);
 153
 154   for (s=0; s1[s]!=NULL; s++);
 155   n_seq = s;
 156   for (s=0; s2[s]!=NULL; s++);
 157   if (n_seq != s) nrerror("unequal number of sequences in aliduplexfold()\n");
 158
 159   if ((!P) || (fabs(P->temperature - temperature)>1e-6)) {
 160     snoupdate_fold_params();  if(P) free(P); P = scale_parameters();
 161     make_pair_matrix();
 162   }
 163
 164   c = (int **) space(sizeof(int *) * (n1+1));
 165   r = (int **) space(sizeof(int *) * (n1+1));
 166   for (i=0; i<=n1; i++) {
 167           c[i] = (int *) space(sizeof(int) * (n2+1));
 168         r[i] = (int *) space(sizeof(int) * (n2+1));
 169         for(j=n2; j>-1; j--){
 170                 c[i][j]=INF;
 171                 r[i][j]=INF;
 172         }
 173   }
 174   Sali1 = (short **) space((n_seq+1)*sizeof(short *));
 175   Sali2 = (short **) space((n_seq+1)*sizeof(short *));
 176   for (s=0; s<n_seq; s++) {
 177     if (strlen(s1[s]) != n1) nrerror("uneqal seqence lengths");
 178     if (strlen(s2[s]) != n2) nrerror("uneqal seqence lengths");
 179     Sali1[s] = aliencode_seq(s1[s]);
 180     Sali2[s] = aliencode_seq(s2[s]);
 181   }
 182   type = (int *) space(n_seq*sizeof(int));
 183   type2 = (int *) space(n_seq*sizeof(int));
 184   type3 = (int *) space(n_seq*sizeof(int));
 185   /*   encode_seqs(s1, s2); */
 186   for (i=6; i<=n1-5; i++) {
 187     int U; U=0;
 188     for (s=0; s<n_seq; s++) {
 189       U+=Sali1[s][i-2];
 190     }
 191     U = (U==(n_seq)*4?1:0);
 192     for (j=n2-min_d2; j>min_d1; j--) {
 193       int type4, k,l,psc,psc2,psc3;
 194       for (s=0; s<n_seq; s++) {
 195         type[s] = pair[Sali1[s][i]][Sali2[s][j]];
 196       }
 197       psc = covscore(type, n_seq);
 198       for (s=0; s<n_seq; s++) if (type[s]==0) type[s]=7;
 199       c[i][j] = (psc>=MINPSCORE) ? (n_seq*P->DuplexInit) : INF;
 200       if (psc<MINPSCORE) continue;
 201       if(/*  pair[Sali1[i+1]][Sali2[j-1]] &&  */
 202          U && j < max_s1 && j > min_s1 &&
 203          j > n2 - max_s2 - max_half_stem &&
 204          j < n2 -min_s2 -half_stem ) { /*constraint on s2 and i*/
 205         folden *temp;
 206         temp=foldlist[j+1];
 207         while(temp->next){
 208           int k = temp->k;
 209           for (s=0; s<n_seq; s++) {
 210             type2[s]= pair[Sali1[s][i-3]][Sali2[s][k+1]];
 211             type3[s]= pair[Sali1[s][i-4]][Sali2[s][k+1]];
 212           }
 213           psc2 = covscore(type2, n_seq);
 214           psc3 = covscore(type3, n_seq);
 215           if(psc2 > MINPSCORE){
 216             r[i][j]=MIN2(r[i][j],c[i-3][k+1]+temp->energy);
 217           }
 218           if(psc3 > MINPSCORE){
 219             r[i][j]=MIN2(r[i][j],c[i-4][k+1]+temp->energy);
 220           }
 221           temp=temp->next;
 222         }
 223       }
 224       /* dangle 5'SIDE relative to the mRNA  */
 225       for (s=0; s<n_seq; s++) {
 226         c[i][j] += E_ExtLoop(type[s], Sali1[s][i-1],Sali2[s][j+1],P);
 227       }
 228       for (k=i-1; k>0 && (i-k)<MAXLOOP_L; k--) {
 229         for (l=j+1; l<=n2 ; l++) {
 230           if (i-k+l-j>2*MAXLOOP_L-2) break;
 231           if (abs(i-k-l+j) >= ASS ) continue;
 232           for (E=s=0; s<n_seq; s++) {
 233             type4 = pair[Sali1[s][k]][Sali2[s][l]];
 234             if (type4==0) type4=7;
 235             E += E_IntLoop(i-k-1, l-j-1, type4, rtype[type[s]],
 236                            Sali1[s][k+1], Sali2[s][l-1], Sali1[s][i-1], Sali2[s][j+1],P);
 237           }
 238           c[i][j] = MIN2(c[i][j], c[k][l] + E);
 239           r[i][j] = MIN2(r[i][j], r[k][l] + E);
 240         }
 241       }
 242       c[i][j]-=psc;
 243       r[i][j]-=psc;
 244       E = r[i][j];
 245       for (s=0; s<n_seq; s++) {
 246         E+= E_ExtLoop(rtype[type[s]], Sali2[s][j-1], Sali1[s][i+1], P);
 247         /**
 248         *** if (i<n1) E += P->dangle3[rtype[type[s]]][Sali1[s][i+1]];
 249         *** if (j>1)  E += P->dangle5[rtype[type[s]]][Sali2[s][j-1]];
 250         *** if (type[s]>2) E += P->TerminalAU;
 251         **/
 252       }
 253       if (E<Emin) {
 254         Emin=E; i_min=i; j_min=j;
 255       }
 256     }
 257   }
 258   if(Emin > 0){
 259           printf("no target found under the constraints chosen\n");
 260         for (i=0; i<=n1; i++) {free(r[i]);free(c[i]);}
 261         free(c);
 262         free(r);
 263         for(s=0; s<n_seq;s++){
 264           free(Sali1[s]);
 265           free(Sali2[s]);
 266         }
 267         free(Sali1); free(Sali2);
 268         free(S2); free(SS1); free(SS2);free(type);free(type2);free(type3);
 269         mfe.energy=INF;
 270         mfe.structure=NULL;
 271         return mfe;
 272   }
 273   struc = alisnoop_backtrack(i_min, j_min,(const char**) s2,
 274                              &Duplex_El, &Duplex_Er, &Loop_E,
 275                              &Loop_D, &u, &pscd, &psct, &pscg,
 276                              penalty, threshloop, threshLE,
 277                              threshRE,threshDE, threshD,
 278                              half_stem, max_half_stem, min_s2,
 279                              max_s2, min_s1, max_s1, min_d1,
 280                              min_d2,(const short**) Sali1,(const short**) Sali2);
 281   /* if (i_min<n1-5) i_min++; */
 282   /* if (j_min>6 ) j_min--; */
 283   l1 = strchr(struc, '&')-struc;
 284   mfe.i = i_min-5;
 285   mfe.j = j_min-5;
 286   mfe.u = u -5;
 287   mfe.Duplex_Er = (float) Duplex_Er/100;
 288   mfe.Duplex_El = (float) Duplex_El/100;
 289   mfe.Loop_D = (float) Loop_D/100;
 290   mfe.Loop_E = (float) Loop_E/100;
 291   mfe.energy = (float) Emin/100 ;
 292   /* mfe.fullStemEnergy = (float) fullStemEnergy/100; */
 293   mfe.pscd = pscd;
 294   mfe.psct = psct;
 295   mfe.structure = struc;
 296   for(s=0; s<n_seq;s++){
 297     free(Sali1[s]);free(Sali2[s]);
 298   }
 299   free(Sali1);free(Sali2);free(type);free(type2);free(type3);
 300
 301   if (!delay_free) {
 302     for (i=0; i<=n1; i++) {free(r[i]);free(c[i]);}
 303     free(c);
 304     free(r);
 305     free(S2); free(SS1); free(SS2);
 306   }
 307   return mfe;
 308 }
 309
 310 PUBLIC snoopT *alisnoop_subopt(const char **s1, const char **s2, int delta, int w,
 311                               const int penalty, const int threshloop,
 312                                const int threshLE, const int threshRE, const int threshDE, const int threshTE, const int threshSE, const int threshD,
 313                               const int distance, const int half_stem, const int max_half_stem,
 314                                const int min_s2, const int max_s2, const int min_s1, const int max_s1, const int min_d1, const int min_d2) {
 315
 316
 317   short **Sali1, **Sali2;
 318   /* printf("%d %d\n", min_s2, max_s2); */
 319   int i,j,s,n_seq, n1, n2, E, n_subopt=0, n_max;
 320   char *struc;
 321   snoopT mfe;
 322   snoopT *subopt;
 323   int thresh;
 324   int *type;
 325   int Duplex_El, Duplex_Er, Loop_E,pscd,psct,pscg;
 326   int Loop_D;
 327   Duplex_El=0; Duplex_Er=0; Loop_E=0;Loop_D=0;pscd=0;psct=0;pscg=0;
 328   int u;
 329   u=0;
 330   n_max=16;
 331   subopt = (snoopT *) space(n_max*sizeof(snoopT));
 332   delay_free=1;
 333   mfe = alisnoopfold(s1, s2, penalty, threshloop, threshLE, threshRE, threshDE,threshD,
 334                   half_stem, max_half_stem,
 335                      min_s2, max_s2, min_s1, max_s1, min_d1, min_d2);
 336   if(mfe.energy > 0){
 337           free(subopt);
 338         delay_free=0;
 339         return NULL;
 340   }
 341   thresh = MIN2((int) ((mfe.Duplex_Er + mfe.Duplex_El + mfe.Loop_E)*100+0.1 + 410) + delta, threshTE );
 342  /* subopt[n_subopt++]=mfe; */
 343   free(mfe.structure);
 344   n1 = strlen(s1[0]);
 345   n2 = strlen(s2[0]);
 346   for (s=0; s1[s]!=NULL; s++);
 347   n_seq = s;
 348   Sali1 = (short **) space((n_seq+1)*sizeof(short *));
 349   Sali2 = (short **) space((n_seq+1)*sizeof(short *));
 350   for (s=0; s<n_seq; s++) {
 351     if (strlen(s1[s]) != n1) nrerror("uneqal seqence lengths");
 352     if (strlen(s2[s]) != n2) nrerror("uneqal seqence lengths");
 353     Sali1[s] = aliencode_seq(s1[s]);
 354     Sali2[s] = aliencode_seq(s2[s]);
 355   }
 356   Sali1[n_seq]=NULL;  Sali2[n_seq]=NULL;
 357   type = (int *) space(n_seq*sizeof(int));
 358   for (i=n1; i>1; i--){
 359     for (j=1; j<=n2; j++) {
 360       int  ii,jj, Ed,psc,skip;
 361       for (s=0; s<n_seq; s++) {
 362         type[s] = pair[Sali2[s][j]][Sali1[s][i]];
 363       }
 364       psc = covscore(type, n_seq);
 365       for (s=0; s<n_seq; s++) if (type[s]==0) type[s]=7;
 366       if (psc<MINPSCORE) continue;
 367       E = Ed = r[i][j];
 368       for  (s=0; s<n_seq; s++) {
 369         /*         if (i<n1-5) Ed += P->dangle3[type[s]][Sali1[s][i+1]]; */
 370         /*       if (j>6)  Ed += P->dangle5[type[s]][Sali2[s][j-1]]; */
 371         if (type[s]>2) Ed += P->TerminalAU;
 372       }
 373       if (Ed>thresh) continue;
 374       /* too keep output small, remove hits that are dominated by a
 375          better one close (w) by. For simplicity we do test without
 376          adding dangles, which is slightly inaccurate.
 377       */
 378       w=1;
 379       for (skip=0, ii=MAX2(i-w,1); (ii<=MIN2(i+w,n1)) && type; ii++) {
 380         for (jj=MAX2(j-w,1); jj<=MIN2(j+w,n2); jj++)
 381           if (r[ii][jj]<E) {skip=1; break;}
 382       }
 383       if (skip){continue;}
 384       psct=0;
 385       pscg=0;
 386       struc = alisnoop_backtrack(i,j,s2, &Duplex_El,
 387                                  &Duplex_Er, &Loop_E, &Loop_D, &u, &pscd, &psct,&pscg,
 388                                  penalty, threshloop,threshLE,threshRE,threshDE, threshD,
 389                                  half_stem, max_half_stem, min_s2, max_s2, min_s1, max_s1, min_d1, min_d2,(const short int**) Sali1,(const int short **) Sali2);
 390
 391       if (Duplex_Er > threshRE || Duplex_El > threshLE || Loop_D > threshD ||
 392          (Duplex_Er + Duplex_El) > threshDE ||
 393          (Duplex_Er + Duplex_El + Loop_E) > threshTE ||
 394          (Duplex_Er + Duplex_El + Loop_E + Loop_D + 410) > threshSE) {
 395                  /* printf(" Duplex_Er %d threshRE %d Duplex_El %d threshLE %d \n" */
 396                 /*        " Duplex_Er + Duplex_El %d  threshDE %d \n" */
 397                 /*        " Duplex_Er + Duplex_El + Loop_E %d  threshTE %d \n" */
 398                 /*        " Duplex_Er + Duplex_El + Loop_E + Loop_D %d  threshSE %d \n",  */
 399                 /*          Duplex_Er , threshRE , Duplex_El ,threshLE, */
 400                 /*          Duplex_Er + Duplex_El, threshDE, */
 401                 /*          Duplex_Er + Duplex_El+  Loop_E , threshTE, */
 402                 /*          Duplex_Er + Duplex_El+  Loop_E + Loop_D, threshSE);  */
 403                  Duplex_Er=0;
 404                 Duplex_El=0;
 405                 Loop_E = 0;
 406                 Loop_D = 0;
 407                 u=0,
 408                 free(struc);
 409                 continue;
 410         }
 411
 412       if (n_subopt+1>=n_max) {
 413         n_max *= 2;
 414         subopt = (snoopT *) xrealloc(subopt, n_max*sizeof(snoopT));
 415       }
 416
 417       subopt[n_subopt].i = i-5;
 418       subopt[n_subopt].j = j-5;
 419       subopt[n_subopt].u = u-5;
 420       subopt[n_subopt].Duplex_Er = Duplex_Er * 0.01;
 421       subopt[n_subopt].Duplex_El = Duplex_El * 0.01;
 422       subopt[n_subopt].Loop_E = Loop_E * 0.01;
 423       subopt[n_subopt].Loop_D = Loop_D * 0.01;
 424       subopt[n_subopt].energy = (Duplex_Er +Duplex_El + Loop_E + Loop_D + 410) * 0.01 ;
 425       subopt[n_subopt].pscd = pscd * 0.01;
 426       subopt[n_subopt].psct = -psct * 0.01;
 427       subopt[n_subopt++].structure = struc;
 428
 429       /* i=u; */
 430       Duplex_Er=0; Duplex_El=0; Loop_E=0; Loop_D=0;u=0;pscd=0;psct=0;
 431     }
 432   }
 433
 434   for (i=0; i<=n1; i++) {free(c[i]);free(r[i]);}
 435   free(c);free(r);
 436   for (s=0; s<n_seq; s++) {
 437     free(Sali1[s]); free(Sali2[s]);
 438   }
 439   free(Sali1); free(Sali2); free(type);
 440
 441   if (snoop_subopt_sorted)
 442     qsort(subopt, n_subopt, sizeof(snoopT), compare);
 443   subopt[n_subopt].i =0;
 444   subopt[n_subopt].j =0;
 445   subopt[n_subopt].structure = NULL;
 446   return subopt;
 447 }
 448
 449
 450
 451
 452
 453
 454
 455 PRIVATE char *alisnoop_backtrack(int i, int j, const char ** snoseq, int *Duplex_El,
 456                                  int *Duplex_Er, int *Loop_E, int *Loop_D, int *u,
 457                                  int *pscd, int *psct, int *pscg,
 458                                  const int penalty, const int threshloop, const int threshLE,
 459                                  const int threshRE, const int threshDE, const int threshD, const int half_stem,
 460                                  const int max_half_stem,
 461                                  const int min_s2, const int max_s2, const int min_s1,
 462                                  const int max_s1,
 463                                  const int min_d1, const int min_d2,const short **Sali1, const short **Sali2) {
 464   /* backtrack structure going backwards from i, and forwards from j
 465      return structure in bracket notation with & as separator */
 466   int k, l, *type,*type2,*type3,type4, E, traced, i0, j0,s,n_seq,psc;
 467   int traced_r=0; /* flag for following backtrack in c or r */
 468   char *st1, *st2, *struc;
 469   char *struc_loop;
 470   n1 = (int) Sali1[0][0];
 471   n2 = (int) Sali2[0][0];
 472
 473   for (s=0; Sali1[s]!=NULL; s++);
 474   n_seq = s;
 475   for (s=0; Sali2[s]!=NULL; s++);
 476   if (n_seq != s) nrerror("unequal number of sequences in alibacktrack()\n");
 477
 478   st1 = (char *) space(sizeof(char)*(n1+1));
 479   st2 = (char *) space(sizeof(char)*(n2+1));
 480   type = (int *) space(n_seq*sizeof(int));
 481   type2 = (int *) space(n_seq*sizeof(int));
 482   type3 = (int *) space(n_seq*sizeof(int));
 483   int *indx;
 484   int *mLoop;
 485   int *cLoop;
 486   folden **foldlist, **foldlist_XS;
 487   snoexport_fold_arrays(&indx, &mLoop, &cLoop,&foldlist, &foldlist_XS );
 488   i0=i; j0=j; /* MIN2(i+1,n1); j0=MAX2(j-1,1);!modified */
 489   for (s=0; s<n_seq; s++) {
 490     type[s] = pair[Sali1[s][i]][Sali2[s][j]];
 491     if(type[s]==0) type[s] = 7;
 492     *Duplex_Er += E_ExtLoop(rtype[type[s]], (j>1) ? Sali2[s][j-1] : -1, (i<n1) ? Sali1[s][i+1] : -1, P);
 493     /**
 494     *** if (i<n1)   *Duplex_Er += P->dangle3[rtype[type[s]]][Sali1[s][i+1]];
 495     *** if (j>1)    *Duplex_Er += P->dangle5[rtype[type[s]]][Sali2[s][j-1]];
 496     *** if (type[s]>2) *Duplex_Er += P->TerminalAU;
 497     **/
 498   }
 499   while (i>0 && j<=n2-min_d2 ) {
 500     if(!traced_r) {
 501       E = r[i][j]; traced=0;
 502       st1[i-1] = '<';
 503       st2[j-1] = '>';
 504       for (s=0; s<n_seq; s++) {
 505         type[s] = pair[Sali1[s][i]][Sali2[s][j]];
 506       }
 507       psc = covscore(type,n_seq);
 508       for (s=0; s<n_seq; s++) if (type[s]==0) type[s] = 7;
 509       E += psc;
 510       *pscd +=psc;
 511       for (k=i-1; k>0 && (i-k)<MAXLOOP_L; k--) {
 512         for (l=j+1; l<=n2 ; l++) {
 513           int LE;
 514           if (i-k+l-j>2*MAXLOOP_L-2) break;
 515           if (abs(i-k-l+j) >= ASS) continue;
 516           for (s=LE=0; s<n_seq; s++) {
 517             type4 = pair[Sali1[s][k]][Sali2[s][l]];
 518             if (type4==0) type4=7;
 519             LE += E_IntLoop(i-k-1, l-j-1, type4, rtype[type[s]], Sali1[s][k+1], Sali2[s][l-1], Sali1[s][i-1], Sali2[s][j+1],P);
 520           }
 521           if (E == r[k][l]+LE) {
 522             traced=1;
 523             i=k; j=l;
 524             *Duplex_Er+=LE;
 525             break;
 526           }
 527         }
 528         if (traced) break;
 529       }
 530       if(!traced){
 531         int U=0;
 532         for (s=0; s<n_seq; s++) {
 533           U+=Sali1[s][i-2];
 534         }
 535         U = (U==(n_seq)*4?1:0);
 536         if(/*  pair[Sali1[i+1]][Sali2[j-1]] && */   /* only U's are allowed */
 537             U && j < max_s1 && j > min_s1 &&
 538             j > n2 - max_s2 - max_half_stem &&
 539             j < n2 -min_s2 -half_stem ) {
 540           int min_k, max_k;
 541           max_k = MIN2(n2-min_s2,j+max_half_stem+1);
 542           min_k = MAX2(j+half_stem+1, n2-max_s2);
 543           folden * temp;
 544           temp=foldlist[j+1];
 545           while(temp->next) {
 546             int psc2, psc3;
 547             int k = temp->k;
 548             for (s=0; s<n_seq; s++) {
 549               type2[s]= pair[Sali1[s][i-3]][Sali2[s][k+1]];
 550               type3[s]= pair[Sali1[s][i-4]][Sali2[s][k+1]];
 551             }
 552             psc2 = covscore(type2, n_seq);
 553             psc3 = covscore(type3, n_seq);
 554             if(psc2>MINPSCORE /*&& pair[Sali1[i-4]][Sali2[k+2]]*/    ){  /* introduce structure from RNAfold */
 555               if(E==c[i-3][k+1]+temp->energy){
 556                 *Loop_E=temp->energy;
 557                 st1[i-3]= '|';
 558                 *u=i-2;
 559                 int a,b;
 560                 /* int fix_ij=indx[k-1+1]+j+1; */
 561                 for(a=0; a< MISMATCH ;a++){
 562                   for(b=0; b< MISMATCH ; b++){
 563                     int ij=indx[k-1-a+1]+j+1+b;
 564                     if(cLoop[ij]==temp->energy) {
 565                       /* int bla; */
 566                       struc_loop=alisnobacktrack_fold_from_pair(snoseq, j+1+b, k-a-1+1,psct);
 567                     a=INF; b=INF;
 568                     }
 569                   }
 570                 }
 571                 traced=1;
 572                 traced_r=1;
 573                 i=i-3;j=k+1;
 574                 break;
 575               }
 576             }
 577              if (psc3>MINPSCORE  /*&& pair[Sali1[i-5]][Sali2[k+2]]*/){  /* introduce structure from RNAfold */
 578               if(E==c[i-4][k+1]+temp->energy){
 579                 *Loop_E=temp->energy;
 580                 st1[i-3]= '|';
 581                 *u=i-2;
 582                 int a,b;
 583                 /* int fix_ij=indx[k-1+1]+j+1; */
 584                 for(a=0; a< MISMATCH ;a++){
 585                   for(b=0; b< MISMATCH ; b++){
 586                     int ij=indx[k-1-a+1]+j+1+b;
 587                     if(cLoop[ij]==temp->energy) {
 588                       /* int bla; */
 589                       struc_loop=alisnobacktrack_fold_from_pair(snoseq, j+1+b, k-a-1+1,psct);
 590                       a=INF; b=INF;
 591                     }
 592                   }
 593                 }
 594                 traced=1;
 595                 traced_r=1;
 596                 i=i-4;j=k+1;
 597                 break;
 598               }
 599             } /* else if */
 600             temp=temp->next;
 601           } /* while temp-> next */
 602         } /* test on j  */
 603       }/* traced? */
 604     }/* traced_r? */
 605     else{
 606       E = c[i][j]; traced=0;
 607       st1[i-1] = '<';
 608       st2[j-1] = '>';
 609       for (s=0; s<n_seq; s++) {
 610         type[s] = pair[Sali1[s][i]][Sali2[s][j]];
 611       }
 612       psc = covscore(type,n_seq);
 613       for (s=0; s<n_seq; s++) if (type[s]==0) type[s] = 7;
 614       E += psc;
 615       *pscd+=psc;
 616       if (!type) nrerror("backtrack failed in fold duplex c");
 617       for (k=i-1; (i-k)<MAXLOOP_L; k--) {
 618         for (l=j+1; l<=n2; l++) {
 619           int LE;
 620           if (i-k+l-j>2*MAXLOOP_L-2) break;
 621           if (abs(i-k-l+j) >= ASS) continue;
 622           for (s=LE=0; s<n_seq; s++) {
 623             type4 = pair[Sali1[s][k]][Sali2[s][l]];
 624             if (type4==0) type4=7;
 625             LE += E_IntLoop(i-k-1, l-j-1, type4, rtype[type[s]], Sali1[s][k+1], Sali2[s][l-1], Sali1[s][i-1], Sali2[s][j+1],P);
 626           }
 627           if (E == c[k][l]+LE) {
 628             traced=1;
 629             i=k; j=l;
 630             *Duplex_El+=LE;
 631             break;
 632           }
 633         }
 634         if (traced) break;
 635       }
 636     }
 637     if (!traced) {
 638       for (s=0; s<n_seq; s++) {
 639         int correction;
 640         correction = E_ExtLoop(type[s], (i>1) ? Sali1[s][i-1] : -1, (j<n2) ? Sali2[s][j+1] : -1, P);
 641         *Duplex_El += correction;
 642         E          -= correction;
 643         /**
 644         *** if (i>1)    {E -= P->dangle5[type[s]][Sali1[s][i-1]]; *Duplex_El +=P->dangle5[type[s]][Sali1[s][i-1]];}
 645         *** if (j<n2)   {E -= P->dangle3[type[s]][Sali2[s][j+1]]; *Duplex_El +=P->dangle3[type[s]][Sali2[s][j+1]];}
 646         *** if (type[s]>2) {E -= P->TerminalAU;                      *Duplex_El +=P->TerminalAU;}
 647         **/
 648       }
 649       if (E != n_seq * P->DuplexInit) {
 650         nrerror("backtrack failed in fold duplex end");
 651       } else break;
 652     }
 653   }
 654 /*   if (i>1)  i--;  */
 655 /*   if (j<n2) j++;  */
 656   /* struc = (char *) space(i0-i+1+j-j0+1+2); */ /* declare final duplex structure */
 657   struc = (char *) space(i0-i+1+n2-1+1+2); /* declare final duplex structure */
 658   char * struc2;
 659   struc2 = (char *) space(n2+1);
 660   /* char * struct_const; */
 661   for (k=MAX2(i,1); k<=i0; k++) if (!st1[k-1]) st1[k-1] = '.';
 662   /* for (k=j0; k<=j; k++) if (!st2[k-1]) st2[k-1] = struc_loop[k-1];*/ /* '.'; normal */
 663   /*  char * struct_const; */
 664   /*  struct_const = (char *) space(sizeof(char)*(n2+1));   */
 665   for (k=1; k<=n2; k++) {
 666     if (!st2[k-1]) st2[k-1] = struc_loop[k-1];/* '.'; */
 667     struc2[k-1] = st2[k-1];/* '.'; */
 668     /*      if (k>=j0 && k<=j){ */
 669     /*              struct_const[k-1]='x'; */
 670     /*      } */
 671     /*      else{ */
 672     /*              if(k<j0) {struct_const[k-1]='<';} */
 673     /*              if(k>j) {struct_const[k-1]='>';} */
 674     /*      } */
 675   }
 676
 677   /*   char duplexseq_1[j0+1]; */
 678   /* char duplexseq_2[n2-j+3]; */
 679   if(j<n2){
 680     char **duplexseq_1, **duplexseq_2;
 681     duplexseq_1 = (char**) space((n_seq+1) * sizeof(char*));
 682     duplexseq_2 = (char**) space((n_seq+1) * sizeof(char*));
 683     for(s=0; s<n_seq; s++){
 684       duplexseq_1[s] = (char*) space((j0)*sizeof(char)); /* modfied j0+1 */
 685       duplexseq_2[s] = (char*) space((n2-j+2)*sizeof(char)); /* modified j+3 */
 686       strncpy(duplexseq_1[s], snoseq[s], j0-1); /* modified j0 */
 687       strcpy(duplexseq_2[s], snoseq[s] + j); /* modified j-1 */
 688       duplexseq_1[s][j0-1]='\0'; /* modified j0 */
 689       duplexseq_2[s][n2-j+1]='\0';/* modified j+2 */
 690     }
 691     duplexseq_1[n_seq]=NULL;
 692     duplexseq_2[n_seq]=NULL;
 693     duplexT temp;
 694     temp=aliduplexfold((const char**)duplexseq_1, (const char**)duplexseq_2);
 695     *Loop_D =  MIN2(0,-410 + (int) 100 * temp.energy*n_seq);
 696     if(*Loop_D){
 697       int l1,ibegin, iend, jbegin, jend;
 698       l1=strchr(temp.structure, '&')-temp.structure;
 699       ibegin=temp.i-l1;
 700       iend  =temp.i-1;
 701       jbegin=temp.j;
 702       jend  =temp.j+strlen(temp.structure)-l1-2-1;
 703       for(k=ibegin+1; k<=iend+1; k++){
 704         struc2[k-1]=temp.structure[k-ibegin-1];
 705       }
 706       for(k=jbegin+j; k<=jend+j; k++){
 707         struc2[k-1]=temp.structure[l1+k-j-jbegin+1];
 708       }
 709     }
 710     for(s=0; s<n_seq; s++){
 711       free(duplexseq_1[s]);
 712       free(duplexseq_2[s]);
 713     }
 714     free(duplexseq_1);free(duplexseq_2);
 715     free(temp.structure);
 716   }
 717   strcpy(struc, st1+MAX2(i-1,0)); strcat(struc, "&");
 718   /* strcat(struc, st2); */
 719   strncat(struc, struc2+5, strlen(struc2)-10);
 720   free(struc2);
 721   free(struc_loop);
 722   free(st1); free(st2);
 723   free(type);free(type2);free(type3);
 724   /* free_arrays(); */
 725   return struc;
 726 }
 727
 728
 729
 730
 731
 732
 733
 734 void Lsnoop_subopt(const char *s1, const char *s2, int delta, int w,
 735                    const int penalty, const int threshloop,
 736                    const int threshLE, const int threshRE, const int threshDE, const int threshTE,const int threshSE,const int threshD,
 737                    const int distance,
 738                    const int half_stem, const int max_half_stem,
 739                    const int min_s2, const int max_s2, const int min_s1, const int max_s1, const int min_d1, const int min_d2, const int alignment_length, const char* name, const int fullStemEnergy)
 740 {
 741
 742  int min_colonne=INF;
 743  int max_pos;
 744  int max;max=INF;
 745  /* int temp; */
 746  /* int nsubopt=10; */
 747  n1 = (int) strlen(s1);
 748  n2 = (int) strlen(s2);
 749  int *position;
 750  position = (int*) space((n1+3)*sizeof(int));
 751
 752
 753   /* int Eminj, Emin_l; */
 754   int i, j; /* l1, Emin=INF, i_min=0, j_min=0; */
 755   /* char *struc; */
 756   /* snoopT mfe; */
 757   int *indx;
 758   int *mLoop;
 759   int *cLoop;
 760   folden **foldlist, **foldlist_XS;
 761   int Duplex_El, Duplex_Er;
 762   int Loop_D;
 763   /* int u; */
 764   int Loop_E;
 765   Duplex_El=0;Duplex_Er=0;Loop_E=0, Loop_D=0;
 766   snoexport_fold_arrays(&indx, &mLoop, &cLoop, &foldlist, &foldlist_XS);
 767   if ((!P) || (fabs(P->temperature - temperature)>1e-6)) {
 768     snoupdate_fold_params();  if(P) free(P); P = scale_parameters();
 769     make_pair_matrix();
 770   }
 771
 772   lc = (int **) space(sizeof(int *) * (5));
 773   lr = (int **) space(sizeof(int *) * (5));
 774   for (i=0; i<5; i++) {
 775           lc[i] = (int *) space(sizeof(int) * (n2+1));
 776         lr[i] = (int *) space(sizeof(int) * (n2+1));
 777         for(j=n2; j>-1; j--){
 778                 lc[i][j]=INF;
 779                 lr[i][j]=INF;
 780         }
 781   }
 782   encode_seqs(s1, s2);
 783   for (i=1; i<=n1; i++) {
 784       int idx=i%5;
 785       int idx_1=(i-1)%5;
 786       int idx_2=(i-2)%5;
 787       int idx_3=(i-3)%5;
 788       int idx_4=(i-4)%5;
 789     for (j=n2-min_d2; j>min_d1; j--) {
 790       int type, type2, k;
 791       type = pair[S1[i]][S2[j]];
 792       lc[idx][j] = (type) ? P->DuplexInit + 2*penalty : INF;
 793       lr[idx][j] = INF;
 794       if(!type) continue;
 795       if( /*pair[S1[i+1]][S2[j-1]]   &&  check that we have a solid base stack after the mLoop */
 796           j < max_s1 && j > min_s1 &&
 797           j > n2 - max_s2 - max_half_stem &&
 798           j < n2 -min_s2 -half_stem && S1[i-2]==4) { /*constraint on s2 and i*/
 799         int min_k, max_k;
 800         max_k = MIN2(n2-min_s2,j+max_half_stem+1);
 801         min_k = MAX2(j+half_stem+1, n2-max_s2);
 802         for(k=min_k; k <= max_k ; k++){
 803           if(mLoop[indx[k-1]+j+1] < 0){
 804                 }
 805           if(pair[S1[i-3]][S2[k]] /*genau zwei ungepaarte nucleotiden --NU--*/
 806              && mLoop[indx[k-1]+j+1] < threshloop){
 807             lr[idx][j]=MIN2(lr[idx][j], lc[idx_3][k]+mLoop[indx[k-1]+j+1]);
 808           }
 809           else if(pair[S1[i-4]][S2[k]] &&  mLoop[indx[k-1]+j+1] < threshloop){/*--NUN--*/
 810             lr[idx][j]=MIN2(lr[idx][j], lc[idx_4][k]+mLoop[indx[k-1]+j+1]);
 811           }
 812               }
 813       }
 814       /* dangle 5'SIDE relative to the mRNA  */
 815       lc[idx][j] += E_ExtLoop(type, (i>1) ? SS1[i-1] : -1, (j<n2) ? SS2[j+1] : -1, P);
 816       /**
 817       *** if (i>1)    lc[idx][j] += P->dangle5[type][SS1[i-1]];
 818       *** if (j<n2)   lc[idx][j] += P->dangle3[type][SS2[j+1]];
 819       *** if (type>2) lc[idx][j] += P->TerminalAU;
 820       **/
 821
 822       if(j<n2 && i>1){
 823         type2=pair[S1[i-1]][S2[j+1]];
 824               if(type2>0){
 825           lc[idx][j]=MIN2(lc[idx_1][j+1]+E_IntLoop(0,0,type2, rtype[type],SS1[i], SS2[j], SS1[i-1], SS2[j+1], P)+2*penalty, lc[idx][j]);
 826           lr[idx][j]=MIN2(lr[idx_1][j+1]+E_IntLoop(0,0,type2, rtype[type],SS1[i], SS2[j], SS1[i-1], SS2[j+1], P)+2*penalty, lr[idx][j]);
 827               }
 828       }
 829       if(j<n2-1 && i>2){
 830         type2=pair[S1[i-2]][S2[j+2]];
 831         if(type2>0 ){
 832           lc[idx][j]=MIN2(lc[idx_2][j+2]+E_IntLoop(1,1,type2, rtype[type],SS1[i-1], SS2[j+1], SS1[i-1], SS2[j+1], P)+4*penalty, lc[idx][j]);
 833           lr[idx][j]=MIN2(lr[idx_2][j+2]+E_IntLoop(1,1,type2, rtype[type],SS1[i-1], SS2[j+1], SS1[i-1], SS2[j+1], P)+4*penalty, lr[idx][j]);
 834         }
 835       }
 836       if(j<n2-2 && i>3){
 837         type2 = pair[S1[i-3]][S2[j+3]];
 838         if(type2>0){
 839           lc[idx][j]=MIN2(lc[idx_3][j+3]+E_IntLoop(2,2,type2, rtype[type],SS1[i-2], SS2[j+2], SS1[i-1], SS2[j+1], P)+6*penalty,lc[idx][j]);
 840           lr[idx][j]=MIN2(lr[idx_3][j+3]+E_IntLoop(2,2,type2, rtype[type],SS1[i-2], SS2[j+2], SS1[i-1], SS2[j+1], P)+6*penalty,lr[idx][j]);
 841               }
 842       }
 843       /**
 844       *** (type>2?P->TerminalAU:0)+(i<(n1)?P->dangle3[rtype[type]][SS1[i+1]]+penalty:0)+(j>1?P->dangle5[rtype[type]][SS2[j-1]]+penalty:0)
 845       **/
 846       min_colonne=MIN2(lr[idx][j]+E_ExtLoop(rtype[type], (j > 1) ? SS2[j-1] : -1, (i<n1) ? SS1[i+1] : -1, P), min_colonne);
 847       }
 848       position[i]=min_colonne;
 849       if(max>=min_colonne){
 850         max=min_colonne;
 851         max_pos=i;
 852       }
 853       min_colonne=INF;
 854  }
 855
 856   free(S1); free(S2); free(SS1); free(SS2);
 857   if(max<threshTE){
 858         find_max_snoop(s1, s2, max, alignment_length, position, delta,
 859                        distance, penalty, threshloop, threshLE, threshRE, threshDE,
 860                        threshTE, threshSE, threshD, half_stem, max_half_stem, min_s2, max_s2, min_s1, max_s1, min_d1, min_d2,name, fullStemEnergy);
 861    }
 862   for (i=1; i<5; i++) {free(lc[i]);free(lr[i]);}
 863   free(lc[0]);free(lr[0]);
 864   free(lc);free(lr);
 865   free(position);
 866
 867 }
 868
 869
 870
 871 void Lsnoop_subopt_list(const char *s1, const char *s2, int delta, int w,
 872                         const int penalty, const int threshloop,
 873                         const int threshLE, const int threshRE, const int threshDE, const int threshTE,const int threshSE,const int threshD,
 874                         const int distance,
 875                         const int half_stem, const int max_half_stem,
 876                         const int min_s2, const int max_s2, const int min_s1, const int max_s1, const int min_d1, const int min_d2, const int alignment_length,const char *name,const int fullStemEnergy)
 877 {
 878
 879  int min_colonne=INF;
 880  int max_pos;
 881  int max;max=INF;
 882  /* int temp; */
 883  /* int nsubopt=10; */
 884  n1 = (int) strlen(s1);
 885  n2 = (int) strlen(s2);
 886  int *position;
 887  position = (int*) space((n1+3)*sizeof(int));
 888
 889
 890  /* int Eminj, Emin_l; */
 891   int i, j;/*  l1, Emin=INF, i_min=0, j_min=0; */
 892   /* char *struc; */
 893   /* snoopT mfe; */
 894   int *indx;
 895   int *mLoop;
 896   int *cLoop;
 897   folden **foldlist, **foldlist_XS;
 898   int Duplex_El, Duplex_Er;
 899   int Loop_D;
 900   /* int u; */
 901   int Loop_E;
 902   Duplex_El=0;Duplex_Er=0;Loop_E=0, Loop_D=0;
 903   snoexport_fold_arrays(&indx, &mLoop, &cLoop, &foldlist,  &foldlist_XS);
 904   if ((!P) || (fabs(P->temperature - temperature)>1e-6)) {
 905     snoupdate_fold_params();  if(P) free(P); P = scale_parameters();
 906     make_pair_matrix();
 907   }
 908
 909   lpair = (int **) space(sizeof(int *) * (6));
 910   lc    = (int **) space(sizeof(int *) * (6));
 911   lr    = (int **) space(sizeof(int *) * (6));
 912   for (i=0; i<6; i++) {
 913           lc[i] = (int *) space(sizeof(int) * (n2+1));
 914         lr[i] = (int *) space(sizeof(int) * (n2+1));
 915         lpair[i] = (int *) space(sizeof(int) * (n2+1));
 916         for(j=n2; j>-1; j--){
 917                 lc[i][j]=INF;
 918                 lr[i][j]=INF;
 919                 lpair[i][j]=0;
 920         }
 921   }
 922   encode_seqs(s1, s2);
 923   int lim_maxj=n2-min_d2 ;
 924   int lim_minj=min_d1;
 925   int lim_maxi=n1;
 926   for (i=5; i<=lim_maxi; i++) {
 927     int idx=i%5;
 928     int idx_1=(i-1)%5;
 929     int idx_2=(i-2)%5;
 930     int idx_3=(i-3)%5;
 931     int idx_4=(i-4)%5;
 932
 933     for (j=lim_maxj; j>lim_minj; j--) {
 934       int type, type2;/*  E, k,l; */
 935       type = pair[S1[i]][S2[j]];
 936       lpair[idx][j] = type;
 937       lc[idx][j] = (type) ? P->DuplexInit + 2*penalty : INF;
 938       lr[idx][j] = INF;
 939       if(!type) continue;
 940       if( /*pair[S1[i+1]][S2[j-1]] && Be sure it binds*/
 941           j < max_s1 && j > min_s1 &&
 942           j > n2 - max_s2 - max_half_stem &&
 943           j < n2 -min_s2 -half_stem && S1[i-2]==4 ) { /*constraint on s2 and i*/
 944         int min_k, max_k;
 945         max_k = MIN2(n2-min_s2,j+max_half_stem+1);
 946         min_k = MAX2(j+half_stem+1, n2-max_s2);
 947         folden * temp;
 948         temp=foldlist[j+1];
 949         while(temp->next){
 950           int k = temp->k;
 951           /* if(k >= min_k-1 && k < max_k){ comment to recover normal behaviour */
 952           if(lpair[idx_3][k+1] /*&& lpair[idx_4][k+2]*/){
 953             lr[idx][j]=MIN2(lr[idx][j], lc[idx_3][k+1]+temp->energy);/*--NU--*/
 954           }
 955           /*else*/ if(lpair[idx_4][k+1]){/*--NUN--*/
 956             lr[idx][j]=MIN2(lr[idx][j], lc[idx_4][k+1]+temp->energy);
 957           }
 958             /*  } */
 959           temp=temp->next;
 960         }
 961       }
 962       /* dangle 5'SIDE relative to the mRNA  */
 963       lc[idx][j] += E_ExtLoop(type,  SS1[i-1] , SS2[j+1] , P);
 964       /**
 965       ***      lc[idx][j] += P->dangle5[type][SS1[i-1]];
 966       ***      lc[idx][j] += P->dangle3[type][SS2[j+1]];
 967       ***      if (type>2) lc[idx][j] += P->TerminalAU;
 968       **/
 969       /*       if(j<n2 && i>1){ */
 970       /* type2=pair[S1[i-1]][S2[j+1]]; */
 971         type2=lpair[idx_1][j+1];
 972         if(type2>0 ){
 973           lc[idx][j]=MIN2(lc[idx_1][j+1]+E_IntLoop(0,0,type2, rtype[type],SS1[i], SS2[j], SS1[i-1], SS2[j+1], P)+2*penalty, lc[idx][j]);
 974           lr[idx][j]=MIN2(lr[idx_1][j+1]+E_IntLoop(0,0,type2, rtype[type],SS1[i], SS2[j], SS1[i-1], SS2[j+1], P)+2*penalty, lr[idx][j]);
 975         }
 976         /* } */
 977         /*       if(j<n2-1 && i>2){ */
 978         /* type2=pair[S1[i-2]][S2[j+2]]; */
 979         type2=lpair[idx_2][j+2];
 980         if(type2>0 ){
 981           lc[idx][j]=MIN2(lc[idx_2][j+2]+E_IntLoop(1,1,type2, rtype[type],SS1[i-1], SS2[j+1], SS1[i-1], SS2[j+1], P), lc[idx][j]);
 982           lr[idx][j]=MIN2(lr[idx_2][j+2]+E_IntLoop(1,1,type2, rtype[type],SS1[i-1], SS2[j+1], SS1[i-1], SS2[j+1], P), lr[idx][j]);
 983           /*         } */
 984       }
 985         /*       if(j<n2-2 && i>3){ */
 986         /* type2 = pair[S1[i-3]][S2[j+3]]; */
 987         type2 =lpair[idx_3][j+3];
 988         if(type2>0 ){
 989           lc[idx][j]=MIN2(lc[idx_3][j+3]+E_IntLoop(2,2,type2, rtype[type],SS1[i-2], SS2[j+2], SS1[i-1], SS2[j+1], P)+6*penalty,lc[idx][j]);
 990           lr[idx][j]=MIN2(lr[idx_3][j+3]+E_IntLoop(2,2,type2, rtype[type],SS1[i-2], SS2[j+2], SS1[i-1], SS2[j+1], P)+6*penalty,lr[idx][j]);
 991           /*         } */
 992       }
 993       /* min_colonne=MIN2(lr[idx][j]+(type>2?P->TerminalAU:0)+P->dangle3[rtype[type]][SS1[i+1]]+P->dangle5[rtype[type]][SS2[j-1]], min_colonne); */
 994       int bla;
 995       bla=lr[idx][j]+E_ExtLoop(rtype[type], SS2[j-1] , SS1[i+1], P)+2*penalty;
 996       min_colonne=MIN2(bla, min_colonne);
 997     }
 998     position[i]=min_colonne;
 999     if(max>=min_colonne){
1000       max=min_colonne;
1001       max_pos=i;
1002       }
1003     min_colonne=INF;
1004  }
1005
1006   free(S1); free(S2); free(SS1); free(SS2);
1007   if(max<threshTE){
1008       find_max_snoop(s1, s2, max, alignment_length, position,
1009                      delta, distance, penalty, threshloop,
1010                      threshLE, threshRE, threshDE, threshTE, threshSE, threshD,
1011                      half_stem, max_half_stem, min_s2, max_s2, min_s1, max_s1, min_d1, min_d2,name, fullStemEnergy);
1012    }
1013   for (i=1; i<6; i++) {free(lc[i]);free(lr[i]);free(lpair[i]);}
1014   free(lc[0]);free(lr[0]);free(lpair[0]);
1015   free(lc);free(lr);free(lpair);
1016   free(position);
1017 }
1018
1019
1020 PRIVATE void find_max_snoop(const char *s1, const char *s2,const int max,  const int alignment_length, const int* position, const int delta,
1021                             const int distance, const int penalty, const int threshloop,  const int threshLE, const int threshRE,
1022                             const int threshDE, const int threshTE, const int threshSE, const int threshD,
1023                             const int half_stem, const int max_half_stem, const int min_s2, const int max_s2, const int min_s1, const int max_s1, const int min_d1, const int min_d2, const char* name, const int fullStemEnergy)
1024 {
1025   int count=0;
1026   int pos=n1+1;
1027   int threshold = MIN2(threshTE , max + delta );
1028   /* printf("threshTE %d max %d\n", threshTE, max); */
1029   /* #pragma omp parallel for */
1030   /*   for(pos=n1+1;pos>distance;pos--){ */
1031   while(pos-- > 5){
1032     int temp_min=0;
1033     if(position[pos]<(threshold)){
1034       int search_range;
1035       search_range=distance+1;
1036       while(--search_range){
1037         if(position[pos-search_range]<=position[pos-temp_min]){
1038           temp_min=search_range;
1039         }
1040       }
1041       pos-=temp_min;
1042       int begin=MAX2(6, pos-alignment_length+1);
1043       char *s3 = (char*) space(sizeof(char)*(pos-begin+3+12));
1044       strcpy(s3, "NNNNN");
1045       strncat(s3, (s1+begin-1), pos-begin+2);
1046       strcat(s3,"NNNNN\0");
1047       /* printf("%s s3\n", s3);  */
1048       snoopT test;
1049       test = snoopfold(s3, s2, penalty, threshloop, threshLE, threshRE, threshDE, threshD, half_stem, max_half_stem, min_s2, max_s2, min_s1,
1050                        max_s1, min_d1, min_d2,fullStemEnergy);
1051       if(test.energy==INF){
1052         free(s3);
1053         continue;
1054       }
1055       if(test.Duplex_El > threshLE * 0.01 || test.Duplex_Er > threshRE * 0.01  ||
1056          test.Loop_D > threshD * 0.01 || (test.Duplex_Er + test.Duplex_El) > threshDE * 0.01 ||
1057          (test.Duplex_Er + test.Duplex_El + test.Loop_E + test.Loop_D + 410) > threshSE*0.01) {
1058         free(test.structure);free(s3);
1059         continue;
1060       }
1061       int l1;
1062       l1 = strchr(test.structure, '&')-test.structure;
1063
1064       int shift=0;
1065       if(test.i > strlen(s3)-10){
1066         test.i--;
1067         l1--;
1068       }
1069       if(test.i-l1<0){
1070         l1--;
1071         shift++;
1072       }
1073       char *target_struct =  (char*) space(sizeof(char) * (strlen(test.structure)+1));
1074       strncpy(target_struct, test.structure+shift, l1);
1075       strncat(target_struct, test.structure + (strchr(test.structure, '&')-
1076                                                test.structure), strlen(test.structure) - (strchr(test.structure, '&')-
1077                                                                                           test.structure));
1078       strcat(target_struct,"\0");
1079       char *target;
1080       target = (char *) space(l1+1);
1081       strncpy(target, (s3+test.i+5-l1), l1);
1082       target[l1]='\0';
1083       char *s4;
1084       s4 = (char*) space(sizeof(char) *(strlen(s2)-9));
1085       strncpy(s4, s2+5, strlen(s2)-10);
1086       s4[strlen(s2)-10]='\0';
1087       printf("%s %3d,%-3d;%3d : %3d,%-3d (%5.2f = %5.2f + %5.2f + %5.2f + %5.2f + 4.1 ) (%5.2f) \n%s&%s\n",
1088              target_struct,begin + test.i-5-l1,begin + test.i -6 , begin + test.u -6,
1089              test.j+1, test.j + (strrchr(test.structure,'>') - strchr(test.structure,'>'))+1 ,
1090              test.Loop_D + test.Duplex_El + test.Duplex_Er + test.Loop_E + 4.10, test.Duplex_El,
1091              test.Duplex_Er, test.Loop_E, test.Loop_D,test.fullStemEnergy, target,s4);
1092       if(name){
1093         char *temp_seq;
1094         char *temp_struc;
1095         char psoutput[100];
1096         temp_seq = (char*) space(sizeof(char)*(l1+n2-9));
1097         temp_struc = (char*) space(sizeof(char)*(l1+n2-9));
1098         strcpy(temp_seq, target);
1099         strcat(temp_seq, s4);
1100         strncpy(temp_struc, target_struct, l1);
1101         strcat(temp_struc, target_struct+l1+1);
1102         temp_seq[n2+l1-10]='\0';
1103         temp_struc[n2+l1-10]='\0';
1104         cut_point = l1+1;
1105         char str[16];char upos[16];
1106         strcpy(psoutput,"sno_");
1107         sprintf(str,"%d",count);
1108         strcat(psoutput,str);
1109         sprintf(upos,"%d",begin + test.u - 6);
1110         strcat(psoutput,"_u_");
1111         strcat(psoutput,upos);
1112         strcat(psoutput,"_");
1113         strcat(psoutput,name);
1114         strcat(psoutput,".ps\0");
1115         PS_rna_plot_snoop_a(temp_seq, temp_struc, psoutput, NULL, NULL);
1116         cut_point = -1;
1117         free(temp_seq);
1118         free(temp_struc);
1119         count++;
1120         /* free(psoutput); */
1121       }
1122       free(s4);
1123       free(test.structure);
1124       free(target_struct);
1125       free(target);
1126       free(s3);
1127     }
1128   }
1129
1130 }
1131
1132
1133
1134
1135
1136
1137
1138
1139 snoopT snoopfold(const char *s1, const char *s2,
1140                  const int penalty, const int threshloop, const int threshLE, const int threshRE, const int threshDE,
1141                  const int threshD,
1142                  const int half_stem, const int max_half_stem,
1143                  const int min_s2, const int max_s2, const int min_s1, const int max_s1, const int min_d1, const int min_d2, const int fullStemEnergy) {
1144   /* int Eminj, Emin_l; */
1145   int i, j, l1, Emin=INF, i_min=0, j_min=0;
1146   char *struc;
1147   snoopT mfe;
1148   int *indx;
1149   int *mLoop;
1150   int *cLoop;
1151   folden** foldlist, **foldlist_XS;
1152   int Duplex_El, Duplex_Er;
1153   int Loop_D;
1154   int u;
1155   int Loop_E;
1156   Duplex_El=0;Duplex_Er=0;Loop_E=0, Loop_D=0;
1157   snoexport_fold_arrays(&indx, &mLoop, &cLoop,&foldlist, &foldlist_XS );
1158   n1 = (int) strlen(s1);
1159   n2 = (int) strlen(s2);
1160
1161   if ((!P) || (fabs(P->temperature - temperature)>1e-6)) {
1162     snoupdate_fold_params();  if(P) free(P); P = scale_parameters();
1163     make_pair_matrix();
1164   }
1165
1166   c = (int **) space(sizeof(int *) * (n1+1));
1167   r = (int **) space(sizeof(int *) * (n1+1));
1168   for (i=0; i<=n1; i++) {
1169           c[i] = (int *) space(sizeof(int) * (n2+1));
1170         r[i] = (int *) space(sizeof(int) * (n2+1));
1171         for(j=n2; j>-1; j--){
1172                 c[i][j]=INF;
1173                 r[i][j]=INF;
1174         }
1175   }
1176   encode_seqs(s1, s2);
1177   for (i=6; i<=n1-5; i++) {
1178     for (j=n2-min_d2; j>min_d1; j--) {
1179       int type, type2, E, k,l;
1180       type = pair[S1[i]][S2[j]];
1181       c[i][j] = (type ) ? P->DuplexInit : INF;
1182       if(!type) continue;
1183       if(/*  pair[S1[i+1]][S2[j-1]] &&  */
1184          j < max_s1 && j > min_s1 &&
1185          j > n2 - max_s2 - max_half_stem &&
1186          j < n2 -min_s2 -half_stem && S1[i-2]==4  ) { /*constraint on s2 and i*/
1187         int min_k, max_k;
1188         max_k = MIN2(n2-min_s2,j+max_half_stem);
1189         min_k = MAX2(j+half_stem, n2-max_s2);
1190         folden * temp;
1191         temp=foldlist[j+1];
1192         while(temp->next){
1193           int k = temp->k;
1194           /* if(k >= min_k-1 && k < max_k){ uncomment to recovernormal behaviour */
1195           if(pair[S1[i-3]][S2[k+1]] /*&& pair[S1[i-4]][S2[k+2]]*/ ){
1196             r[i][j]=MIN2(r[i][j], c[i-3][k+1]+temp->energy);
1197           }
1198           /*else*/ if(pair[S1[i-4]][S2[k+1]] /*&& pair[S1[i-5]][S2[k+2]]*/ ){
1199             r[i][j]=MIN2(r[i][j], c[i-4][k+1]+temp->energy);
1200           }
1201           /* } */
1202           temp=temp->next;
1203         }
1204       }
1205       /* dangle 5'SIDE relative to the mRNA  */
1206       /**
1207       *** c[i][j] += P->dangle5[type][SS1[i-1]];
1208       *** c[i][j] += P->dangle3[type][SS2[j+1]];
1209       *** if (type>2) c[i][j] += P->TerminalAU;
1210       **/
1211       c[i][j]+=E_ExtLoop(type, SS1[i-1] , SS2[j+1], P);
1212       for (k=i-1; k>0 && (i-k)<MAXLOOP_L; k--) {
1213         for (l=j+1; l<=n2 ; l++) {
1214           if (i-k+l-j>2*MAXLOOP_L-2) break;
1215           if (abs(i-k-l+j) >= ASS ) continue;
1216           type2 = pair[S1[k]][S2[l]];
1217           if (!type2) continue;
1218           E = E_IntLoop(i-k-1, l-j-1, type2, rtype[type],
1219                         SS1[k+1], SS2[l-1], SS1[i-1], SS2[j+1],P);
1220           c[i][j] = MIN2(c[i][j], c[k][l]+E+(i-k+l-j)*penalty);
1221           r[i][j] = MIN2(r[i][j], r[k][l]+E+(i-k+l-j)*penalty);
1222         }
1223       }
1224       E = r[i][j];
1225       /**
1226       *** if (i<n1) E += P->dangle3[rtype[type]][SS1[i+1]];
1227       *** if (j>1)  E += P->dangle5[rtype[type]][SS2[j-1]];
1228       *** f (type>2) E += P->TerminalAU;
1229       **/
1230       E+=E_ExtLoop(rtype[type], (j > 1) ? SS2[j-1] : -1, (i<n1) ? SS1[i+1] : -1, P);
1231       if (E<Emin) {
1232         Emin=E; i_min=i; j_min=j;
1233       }
1234     }
1235   }
1236   if(Emin > 0){
1237           printf("no target found under the constraints chosen\n");
1238         for (i=0; i<=n1; i++) {free(r[i]);free(c[i]);}
1239         free(c);
1240         free(r);
1241         free(S1); free(S2); free(SS1); free(SS2);
1242         mfe.energy=INF;
1243         return mfe;
1244   }
1245   struc = snoop_backtrack(i_min, j_min,s2, &Duplex_El, &Duplex_Er, &Loop_E, &Loop_D,
1246                           &u, penalty, threshloop, threshLE, threshRE,threshDE, threshD,
1247                             half_stem, max_half_stem, min_s2, max_s2, min_s1, max_s1, min_d1, min_d2);
1248 /*   if (i_min<n1-5) i_min++; */
1249 /*   if (j_min>1 ) j_min--; */
1250   l1 = strchr(struc, '&')-struc;
1251   mfe.i = i_min-5;
1252   mfe.j = j_min-5;
1253   mfe.u = u -5;
1254   mfe.Duplex_Er = (float) Duplex_Er/100;
1255   mfe.Duplex_El = (float) Duplex_El/100;
1256   mfe.Loop_D = (float) Loop_D/100;
1257   mfe.Loop_E = (float) Loop_E/100;
1258   mfe.energy = (float) Emin/100 ;
1259   mfe.fullStemEnergy = (float) fullStemEnergy/100;
1260   mfe.structure = struc;
1261   if (!delay_free) {
1262     for (i=0; i<=n1; i++) {free(r[i]);free(c[i]);}
1263     free(c);
1264     free(r);
1265     free(S1); free(S2); free(SS1); free(SS2);
1266   }
1267   return mfe;
1268 }
1269
1270 PRIVATE int snoopfold_XS_fill(const char *s1, const char *s2, const int **access_s1,
1271                  const int penalty, const int threshloop, const int threshLE, const int threshRE, const int threshDE,
1272                  const int threshD,
1273                  const int half_stem, const int max_half_stem,
1274                  const int min_s2, const int max_s2, const int min_s1, const int max_s1, const int min_d1, const int min_d2) {
1275   /* int Eminj, Emin_l; */
1276   int i, j, Emin=INF, i_min=0, j_min=0;
1277   /* char *struc; */
1278   /* snoopT mfe; */
1279   int *indx;
1280   int *mLoop;
1281   int *cLoop;
1282   folden** foldlist, **foldlist_XS;
1283   int Duplex_El, Duplex_Er;
1284   int Loop_D;
1285   /* int u; */
1286   int Loop_E;
1287   Duplex_El=0;Duplex_Er=0;Loop_E=0, Loop_D=0;
1288   snoexport_fold_arrays(&indx, &mLoop, &cLoop,&foldlist, &foldlist_XS );
1289   n1 = (int) strlen(s1);
1290   n2 = (int) strlen(s2);
1291
1292   if ((!P) || (fabs(P->temperature - temperature)>1e-6)) {
1293     snoupdate_fold_params();  if(P) free(P); P = scale_parameters();
1294     make_pair_matrix();
1295   }
1296
1297   c_fill = (int **) space(sizeof(int *) * (n1+1));
1298   r_fill = (int **) space(sizeof(int *) * (n1+1));
1299   for (i=0; i<=n1; i++) {
1300           c_fill[i] = (int *) space(sizeof(int) * (n2+1));
1301         r_fill[i] = (int *) space(sizeof(int) * (n2+1));
1302         for(j=n2; j>-1; j--){
1303                 c_fill[i][j]=INF;
1304                 r_fill[i][j]=INF;
1305         }
1306   }
1307   encode_seqs(s1, s2);
1308
1309   int di[5];
1310   di[0]=0;
1311   for (i=6; i<=n1-5; i++) {
1312     di[1]=access_s1[5][i]   - access_s1[4][i-1];
1313     di[2]=access_s1[5][i-1] - access_s1[4][i-2] + di[1];
1314     di[3]=access_s1[5][i-2] - access_s1[4][i-3] + di[2];
1315     di[4]=access_s1[5][i-3] - access_s1[4][i-4] + di[3];
1316     di[1]=MIN2(di[1],165);
1317     di[2]=MIN2(di[2],330);
1318     di[3]=MIN2(di[3],495);
1319     di[4]=MIN2(di[4],660);
1320     for (j=n2-min_d2; j>min_d1; j--) {
1321       int type, type2, E, k,l;
1322       type = pair[S1[i]][S2[j]];
1323       c_fill[i][j] = (type ) ? P->DuplexInit : INF;
1324       if(!type) continue;
1325       if(/*  pair[S1[i+1]][S2[j-1]] &&  */
1326          j < max_s1 && j > min_s1 &&
1327          j > n2 - max_s2 - max_half_stem &&
1328          j < n2 -min_s2 -half_stem && S1[i-2]==4  ) { /*constraint on s2 and i*/
1329         int min_k, max_k;
1330         max_k = MIN2(n2-min_s2,j+max_half_stem);
1331         min_k = MAX2(j+half_stem, n2-max_s2);
1332         folden * temp;
1333         temp=foldlist[j+1];
1334         while(temp->next){
1335           int k = temp->k;
1336           /* if(k >= min_k-1 && k < max_k){ uncomment to recovernormal behaviour */
1337           if(pair[S1[i-3]][S2[k+1]] /*&& pair[S1[i-4]][S2[k+2]]*/ ){
1338             r_fill[i][j]=MIN2(r_fill[i][j], c_fill[i-3][k+1]+temp->energy+ di[3]);
1339           }
1340           /*else*/ if(pair[S1[i-4]][S2[k+1]] /*&& pair[S1[i-5]][S2[k+2]]*/ ){
1341             r_fill[i][j]=MIN2(r_fill[i][j], c_fill[i-4][k+1]+temp->energy + di[4]);
1342           }
1343           /* } */
1344           temp=temp->next;
1345         }
1346       }
1347       /* dangle 5'SIDE relative to the mRNA  */
1348       /**
1349       *** c_fill[i][j] += P->dangle5[type][SS1[i-1]];
1350       *** c_fill[i][j] += P->dangle3[type][SS2[j+1]];
1351       *** if (type>2) c_fill[i][j] += P->TerminalAU;
1352       **/
1353       c_fill[i][j]+= E_ExtLoop(type, SS1[i-1], SS2[j+1],P);
1354       for (k=i-1; k>0 && (i-k)<MAXLOOP_L; k--) {
1355         for (l=j+1; l<=n2 ; l++) {
1356           if (i-k+l-j>2*MAXLOOP_L-2) break;
1357           if (abs(i-k-l+j) >= ASS ) continue;
1358           type2 = pair[S1[k]][S2[l]];
1359           if (!type2) continue;
1360           E = E_IntLoop(i-k-1, l-j-1, type2, rtype[type],
1361                         SS1[k+1], SS2[l-1], SS1[i-1], SS2[j+1],P);
1362           c_fill[i][j] = MIN2(c_fill[i][j], c_fill[k][l]+E+di[i-k]);
1363           r_fill[i][j] = MIN2(r_fill[i][j], r_fill[k][l]+E+di[i-k]);
1364         }
1365       }
1366       E = r_fill[i][j];
1367       /**
1368       ***  if (i<n1) E += P->dangle3[rtype[type]][SS1[i+1]];
1369       ***  if (j>1)  E += P->dangle5[rtype[type]][SS2[j-1]];
1370       *** if (type>2) E += P->TerminalAU;
1371       **/
1372       E+= E_ExtLoop(rtype[type], (j > 1) ? SS2[j-1] : -1, (i<n1) ? SS1[i+1] : -1, P);
1373       if (E<Emin) {
1374         Emin=E; i_min=i; j_min=j;
1375       }
1376     }
1377   }
1378   return Emin;
1379 }
1380
1381
1382
1383 PUBLIC snoopT *snoop_subopt(const char *s1, const char *s2, int delta, int w,
1384                             const int penalty, const int threshloop,
1385                             const int threshLE, const int threshRE, const int threshDE, const int threshTE, const int threshSE, const int threshD,
1386                             const int distance, const int half_stem, const int max_half_stem,
1387                             const int min_s2, const int max_s2, const int min_s1, const int max_s1, const int min_d1, const int min_d2, const int fullStemEnergy) {
1388
1389
1390
1391   /* printf("%d %d\n", min_s2, max_s2); */
1392   int i,j, n1, n2, E, n_subopt=0, n_max;
1393   char *struc;
1394   snoopT mfe;
1395   snoopT *subopt;
1396   int thresh;
1397
1398   int Duplex_El, Duplex_Er, Loop_E;
1399   int Loop_D;
1400   Duplex_El=0; Duplex_Er=0; Loop_E=0;Loop_D=0;
1401   int u;
1402   u=0;
1403   n_max=16;
1404   subopt = (snoopT *) space(n_max*sizeof(snoopT));
1405   delay_free=1;
1406   mfe = snoopfold(s1, s2, penalty, threshloop, threshLE, threshRE, threshDE,threshD,
1407                   half_stem, max_half_stem,
1408                   min_s2, max_s2, min_s1, max_s1, min_d1, min_d2, fullStemEnergy);
1409
1410
1411
1412   if(mfe.energy > 0){
1413           free(subopt);
1414         delay_free=0;
1415         return NULL;
1416   }
1417   thresh = MIN2((int) ((mfe.Duplex_Er + mfe.Duplex_El + mfe.Loop_E)*100+0.1 + 410) + delta, threshTE );
1418  /* subopt[n_subopt++]=mfe; */
1419   free(mfe.structure);
1420
1421   n1 = strlen(s1); n2=strlen(s2);
1422   for (i=n1; i>0; i--) {
1423     for (j=1; j<=n2; j++) {
1424       int type, Ed;
1425       type = pair[S2[j]][S1[i]];
1426       if (!type) continue;
1427       E = Ed = r[i][j];
1428       /**
1429       *** if (i<n1) Ed += P->dangle3[type][SS1[i+1]];
1430       *** if (j>1)  Ed += P->dangle5[type][SS2[j-1]];
1431       *** if (type>2) Ed += P->TerminalAU;
1432       **/
1433       Ed+= E_ExtLoop(type, (j > 1) ? SS2[j-1] : -1, (i<n1) ? SS1[i+1] : -1, P);
1434       if (Ed>thresh) continue;
1435       /* too keep output small, remove hits that are dominated by a
1436          better one close (w) by. For simplicity we do test without
1437          adding dangles, which is slightly inaccurate.
1438       */
1439       /* w=1; */
1440 /*       for (ii=MAX2(i-w,1); (ii<=MIN2(i+w,n1)) && type; ii++) {  */
1441 /*         for (jj=MAX2(j-w,1); jj<=MIN2(j+w,n2); jj++) */
1442 /*           if (r[ii][jj]<E) {type=0; break;} */
1443 /*       } */
1444       if (!type) continue;
1445
1446       struc = snoop_backtrack(i,j,s2, &Duplex_El, &Duplex_Er, &Loop_E, &Loop_D, &u, penalty, threshloop,threshLE,threshRE,threshDE,threshD,
1447                         half_stem, max_half_stem, min_s2, max_s2, min_s1, max_s1, min_d1, min_d2);
1448       if (Duplex_Er > threshRE || Duplex_El > threshLE || Loop_D > threshD ||
1449          (Duplex_Er + Duplex_El) > threshDE ||
1450          (Duplex_Er + Duplex_El + Loop_E) > threshTE ||
1451          (Duplex_Er + Duplex_El + Loop_E + Loop_D + 410) > threshSE) {
1452                  /* printf(" Duplex_Er %d threshRE %d Duplex_El %d threshLE %d \n" */
1453                 /*        " Duplex_Er + Duplex_El %d  threshDE %d \n" */
1454                 /*        " Duplex_Er + Duplex_El + Loop_E %d  threshTE %d \n" */
1455                 /*        " Duplex_Er + Duplex_El + Loop_E + Loop_D %d  threshSE %d \n",  */
1456                 /*          Duplex_Er , threshRE , Duplex_El ,threshLE, */
1457                 /*          Duplex_Er + Duplex_El, threshDE, */
1458                 /*          Duplex_Er + Duplex_El+  Loop_E , threshTE, */
1459                 /*          Duplex_Er + Duplex_El+  Loop_E + Loop_D, threshSE);  */
1460                  Duplex_Er=0;
1461                 Duplex_El=0;
1462                 Loop_E = 0;
1463                 Loop_D = 0;
1464                 u=0,
1465                 free(struc);
1466                 continue;
1467         }
1468
1469       if (n_subopt+1>=n_max) {
1470         n_max *= 2;
1471         subopt = (snoopT *) xrealloc(subopt, n_max*sizeof(snoopT));
1472       }
1473       subopt[n_subopt].i = i-5;
1474       subopt[n_subopt].j = j-5;
1475       subopt[n_subopt].u = u-5;
1476       subopt[n_subopt].Duplex_Er = Duplex_Er * 0.01;
1477       subopt[n_subopt].Duplex_El = Duplex_El * 0.01;
1478       subopt[n_subopt].Loop_E = Loop_E * 0.01;
1479       subopt[n_subopt].Loop_D = Loop_D * 0.01;
1480       subopt[n_subopt].energy = (Duplex_Er +Duplex_El + Loop_E + Loop_D + 410) * 0.01 ;
1481       subopt[n_subopt].fullStemEnergy = (float) fullStemEnergy * 0.01;
1482       subopt[n_subopt++].structure = struc;
1483
1484       Duplex_Er=0; Duplex_El=0; Loop_E=0; Loop_D=0;u=0;
1485     }
1486   }
1487
1488   for (i=0; i<=n1; i++) {free(c[i]);free(r[i]);}
1489   free(c);free(r);
1490   free(S1); free(S2); free(SS1); free(SS2);
1491   delay_free=0;
1492
1493   if (snoop_subopt_sorted) qsort(subopt, n_subopt, sizeof(snoopT), compare);
1494   subopt[n_subopt].i =0;
1495   subopt[n_subopt].j =0;
1496   subopt[n_subopt].structure = NULL;
1497   return subopt;
1498 }
1499
1500 PUBLIC void snoop_subopt_XS(const char *s1, const char *s2, const int **access_s1, int delta, int w,
1501                             const int penalty, const int threshloop,
1502                             const int threshLE, const int threshRE, const int threshDE, const int threshTE, const int threshSE, const int threshD,
1503                             const int distance, const int half_stem, const int max_half_stem,
1504                             const int min_s2, const int max_s2, const int min_s1, const int max_s1, const int min_d1, const int min_d2, const int alignment_length, const char *name, const int fullStemEnergy) {
1505
1506
1507
1508   /* printf("%d %d\n", min_s2, max_s2); */
1509   int i,j, E,  n_max;
1510   /* char *struc; */
1511   /* snoopT mfe; */
1512
1513   int thresh;
1514
1515   int Duplex_El, Duplex_Er, Loop_E;
1516   int Loop_D;
1517   Duplex_El=0; Duplex_Er=0; Loop_E=0;Loop_D=0;
1518   int u;
1519   u=0;
1520   n_max=16;
1521   delay_free=1;
1522   int Emin = snoopfold_XS_fill(s1, s2, access_s1,penalty, threshloop, threshLE, threshRE, threshDE,threshD,
1523                                half_stem, max_half_stem,
1524                                min_s2, max_s2, min_s1, max_s1, min_d1, min_d2);
1525   if(Emin > 0){
1526     delay_free=0;
1527   }
1528   thresh = MIN2(-100, threshTE +alignment_length*30);
1529   /*   n1=strlen(s1);  */
1530   /*   n2=strlen(s2); */
1531
1532   int n3=strlen(s1);
1533   int n4=strlen(s2);
1534   S1_fill = (short*)space(sizeof(short)*(n3+2));
1535   S2_fill = (short*)space(sizeof(short)*(n4+2));
1536   SS1_fill = (short*)space(sizeof(short)*(n3+1));
1537   SS2_fill = (short*)space(sizeof(short)*(n4+1));
1538   memcpy(S1_fill, S1, sizeof(short)*n3+2);
1539   memcpy(S2_fill, S2, sizeof(short)*n4+2);
1540   memcpy(SS1_fill, SS1, sizeof(short)*n3+1);
1541   memcpy(SS2_fill, SS2, sizeof(short)*n4+1);
1542   free(S1);free(S2);free(SS1);free(SS2);
1543   int count=0;
1544   for (i=n3-5; i>0; i--) {
1545     for (j=1; j<=n4; j++) {
1546       int type, Ed;
1547       type = pair[S2_fill[j]][S1_fill[i]];
1548       if (!type) continue;
1549       E = Ed = r_fill[i][j];
1550       /**
1551       ***if (i<n3) Ed += P->dangle3[type][SS1_fill[i+1]];
1552       ***if (j>1)  Ed += P->dangle5[type][SS2_fill[j-1]];
1553       ***if (type>2) Ed += P->TerminalAU;
1554       **/
1555       Ed+=E_ExtLoop(type, (j > 1) ? SS2[j-1] : -1, (i<n3) ? SS1[i+1] : -1, P);
1556       if (Ed>thresh) continue;
1557
1558       /* to keep output small, remove hits that are dominated by a
1559          better one close (w) by. For simplicity we do test without
1560          adding dangles, which is slightly inaccurate.
1561       */
1562 /*       w=10;  */
1563 /*       for (ii=MAX2(i-w,1); (ii<=MIN2(i+w,n3-5)) && type; ii++) {   */
1564 /*         for (jj=MAX2(j-w,1); jj<=MIN2(j+w,n4-5); jj++)  */
1565 /*           if (r_fill[ii][jj]<E) {type=0; break;}  */
1566 /*       }  */
1567 /*       i=ii;j=jj; */
1568       if (!type) continue;
1569       int begin=MAX2(5, i-alignment_length);
1570       int end  =MIN2(n3-5, i-1);
1571       char *s3 = (char*) space(sizeof(char)*(end-begin+2)+5);
1572       strncpy(s3, (s1+begin), end - begin +1);
1573       strcat(s3,"NNNNN\0");
1574       int n5 = strlen(s3);
1575       snoopT test = snoopfold_XS(s3, s2, access_s1, i, j ,penalty,
1576                           threshloop, threshLE, threshRE,
1577                           threshDE, threshD, half_stem,
1578                           max_half_stem, min_s2, max_s2, min_s1,
1579                                  max_s1, min_d1, min_d2,fullStemEnergy);
1580       if(test.energy==INF){
1581         free(s3);
1582         continue;
1583       }
1584       if( test.Duplex_El > threshLE * 0.01 ||test.Duplex_Er > threshRE * 0.01  ||
1585           test.Loop_D > threshD * 0.01 || (test.Duplex_Er + test.Duplex_El) > threshDE * 0.01 ||
1586           (test.Duplex_Er + test.Duplex_El + test.Loop_E) > threshTE*0.01 || (test.Duplex_Er + test.Duplex_El + test.Loop_E + test.Loop_D + 410) > threshSE*0.01)
1587         {
1588           free(test.structure);free(s3);
1589           continue;
1590         }
1591       char *s4;
1592       s4 = (char*) space(sizeof(char) *(n4-9));
1593       strncpy(s4, s2+5, n4-10);
1594       s4[n4-10]='\0';
1595
1596       char *s5 = space(sizeof(char) * n5-test.i+2-5);
1597       strncpy(s5,s3+test.i-1,n5-test.i+1-5);
1598       s5[n5-test.i+1-5]='\0';
1599       float dE = ((float) (access_s1[n5-test.i+1-5][i]))*0.01;
1600       printf("%s %3d,%-3d;%3d : %3d,%-3d (%5.2f = %5.2f + %5.2f + %5.2f + %5.2f + %5.2f + 4.10)  (%5.2f)\n%s&%s\n" ,
1601              test.structure, i  -  (n5 - test.i) ,i - 5, i - (n5 - test.u ),
1602              j-5, j-5 + (strrchr(test.structure,'>') - strchr(test.structure,'>')),
1603              test.Loop_D + test.Duplex_El + test.Duplex_Er + test.Loop_E + 4.10+dE, test.Duplex_El,
1604              test.Duplex_Er, test.Loop_E, test.Loop_D,dE , test.fullStemEnergy, s5,s4);
1605       if(name){
1606         int begin_t, end_t, begin_q, end_q, and, pipe,k;
1607         char psoutput[100];
1608         begin_q=0;
1609         end_q=n4-10;
1610         begin_t=0;
1611         end_t=n5-test.i+ 1-5;
1612         and=end_t+1;
1613         pipe=test.u -test.i + 1;
1614         cut_point = end_t +1 ;
1615         char *catseq, *catstruct;/*  *fname;  */
1616         catseq = (char*) space(n5 + end_q -begin_q +2);
1617         catstruct = (char*) space(n5 + end_q -begin_q +2);
1618         strcpy(catseq, s5);
1619         strncpy(catstruct, test.structure, end_t);
1620         strcat(catseq, s4);
1621         strncat(catstruct, test.structure+end_t+1, end_q-begin_q+1);
1622         catstruct[end_t - begin_t + end_q-begin_q+2]='\0';
1623         catseq[end_t - begin_t + end_q-begin_q+2]='\0';
1624         int *relative_access;
1625         relative_access = space(sizeof(int)*strlen(s5));
1626         relative_access[0] = access_s1[1][i - (n5  - test.i) + 5];
1627         for(k=1;k<strlen(s5);k++){
1628           relative_access[k] =  access_s1[k+1][i - (n5  - test.i) + k + 5] -  access_s1[k][i - (n5  - test.i) + k + 4];
1629         }
1630         char str[16];char upos[16];
1631         strcpy(psoutput,"sno_XS_");
1632         sprintf(str,"%d",count);
1633         strcat(psoutput,str);
1634         sprintf(upos,"%d",i - (n5 - test.u ));
1635         strcat(psoutput,"_u_");
1636         strcat(psoutput,upos);
1637         strcat(psoutput,"_");
1638         strcat(psoutput,name);
1639         strcat(psoutput,".ps\0");
1640         PS_rna_plot_snoop_a(catseq, catstruct, psoutput,relative_access,NULL);
1641         free(catseq);free(catstruct);free(relative_access);
1642         count++;
1643       }
1644       free(s3);free(s4);free(s5);free(test.structure);
1645     }
1646   }
1647   for (i=0; i<=n3; i++) {free(c_fill[i]);free(r_fill[i]);}
1648   free(c_fill);free(r_fill);
1649   free(S1_fill); free(S2_fill); free(SS1_fill); free(SS2_fill);
1650   delay_free=0;
1651 }
1652
1653
1654
1655
1656 PRIVATE char *snoop_backtrack(int i, int j, const char* snoseq,
1657                               int *Duplex_El, int *Duplex_Er,
1658                               int *Loop_E, int *Loop_D, int *u,
1659                               const int penalty, const int threshloop,
1660                               const int threshLE, const int threshRE, const int threshDE, const int threshD,
1661                               const int half_stem, const int max_half_stem,
1662                               const int min_s2, const int max_s2, const int min_s1,
1663                               const int max_s1, const int min_d1, const int min_d2) {
1664   /* backtrack structure going backwards from i, and forwards from j
1665      return structure in bracket notation with & as separator */
1666   int k, l, type, type2, E, traced, i0, j0;
1667   int traced_r=0; /* flag for following backtrack in c or r */
1668   char *st1, *st2, *struc;
1669   char *struc_loop;
1670
1671   st1 = (char *) space(sizeof(char)*(n1+1));
1672   st2 = (char *) space(sizeof(char)*(n2+1));
1673   int *indx;
1674   int *mLoop;
1675   int *cLoop;
1676   folden **foldlist, **foldlist_XS;
1677   type=pair[S1[i]][S2[j]];
1678   snoexport_fold_arrays(&indx, &mLoop, &cLoop,&foldlist, &foldlist_XS );
1679   i0=i; j0=j;
1680   /**
1681   *** if (i<n1)   *Duplex_Er += P->dangle3[rtype[type]][SS1[i+1]];
1682   *** if (j>1)    *Duplex_Er += P->dangle5[rtype[type]][SS2[j-1]];
1683   *** if (type>2) *Duplex_Er += P->TerminalAU;
1684   **/
1685   *Duplex_Er += E_ExtLoop(rtype[type], (j > 1) ? SS2[j-1] : -1, (i<n1) ? SS1[i+1] : -1, P);
1686   while (i>0 && j<=n2-min_d2 ) {
1687     if(!traced_r) {
1688       E = r[i][j]; traced=0;
1689       st1[i-1] = '<';
1690       st2[j-1] = '>';
1691       type = pair[S1[i]][S2[j]];
1692       if (!type) nrerror("backtrack failed in fold duplex r");
1693       for (k=i-1; k>0 && (i-k)<MAXLOOP_L; k--) {
1694         for (l=j+1; l<=n2 ; l++) {
1695           int LE;
1696           if (i-k+l-j>2*MAXLOOP_L-2) break;
1697           if (abs(i-k-l+j) >= ASS) continue;
1698
1699           type2 = pair[S1[k]][S2[l]];
1700           if (!type2) continue;
1701           LE = E_IntLoop(i-k-1, l-j-1, type2, rtype[type],
1702                          SS1[k+1], SS2[l-1], SS1[i-1], SS2[j+1],P);
1703           if (E == r[k][l]+LE+(i-k+l-j)*penalty) {
1704             traced=1;
1705             i=k; j=l;
1706             *Duplex_Er+=LE;
1707             break;
1708           }
1709         }
1710         if (traced) break;
1711       }
1712       if(!traced){
1713         if(/*  pair[S1[i+1]][S2[j-1]] && */
1714             j < max_s1 && j > min_s1 &&
1715             j > n2 - max_s2 - max_half_stem &&
1716             j < n2 -min_s2 -half_stem &&
1717             S1[i-2]==4) {
1718           int min_k, max_k;
1719           max_k = MIN2(n2-min_s2,j+max_half_stem+1);
1720           min_k = MAX2(j+half_stem+1, n2-max_s2);
1721           folden * temp;
1722           temp=foldlist[j+1];
1723           while(temp->next) {
1724             int k = temp->k;
1725             if(pair[S1[i-3]][S2[k+1]] /*&& pair[S1[i-4]][S2[k+2]]*/    ){  /* introduce structure from RNAfold */
1726               if(E==c[i-3][k+1]+temp->energy){
1727                 *Loop_E=temp->energy;
1728                 st1[i-3]= '|';
1729                 *u=i-2;
1730                 int a,b;
1731                 /* int fix_ij=indx[k-1+1]+j+1; */
1732                 for(a=0; a< MISMATCH ;a++){
1733                   for(b=0; b< MISMATCH ; b++){
1734                     int ij=indx[k-1-a+1]+j+1+b;
1735                     if(cLoop[ij]==temp->energy) {
1736                       struc_loop=snobacktrack_fold_from_pair(snoseq, j+1+b, k-a-1+1);
1737                     a=INF; b=INF;
1738                     }
1739                   }
1740                 }
1741                 traced=1;
1742                 traced_r=1;
1743                 i=i-3;j=k+1;
1744                 break;
1745               }
1746             }
1747             /*else*/ if (pair[S1[i-4]][S2[k+1]] /*&& pair[S1[i-5]][S2[k+2]]*/){  /* introduce structure from RNAfold */
1748               if(E==c[i-4][k+1]+temp->energy){
1749                 *Loop_E=temp->energy;
1750                 st1[i-3]= '|';
1751                 *u=i-2;
1752                 int a,b;
1753                 /* int fix_ij=indx[k-1+1]+j+1; */
1754                 for(a=0; a< MISMATCH ;a++){
1755                   for(b=0; b< MISMATCH ; b++){
1756                     int ij=indx[k-1-a+1]+j+1+b;
1757                     if(cLoop[ij]==temp->energy) {
1758                       struc_loop=snobacktrack_fold_from_pair(snoseq, j+1+b, k-a-1+1);
1759                       a=INF; b=INF;
1760                     }
1761                   }
1762                 }
1763                 traced=1;
1764                 traced_r=1;
1765                 i=i-4;j=k+1;
1766                 break;
1767               }
1768             } /* else if */
1769             temp=temp->next;
1770           } /* while temp-> next */
1771         } /* test on j  */
1772       }/* traced? */
1773     }/* traced_r? */
1774     else{
1775       E = c[i][j]; traced=0;
1776       st1[i-1] = '<';
1777       st2[j-1] = '>';
1778       type = pair[S1[i]][S2[j]];
1779       if (!type) nrerror("backtrack failed in fold duplex c");
1780       for (k=i-1; (i-k)<MAXLOOP_L; k--) {
1781         for (l=j+1; l<=n2; l++) {
1782           int LE;
1783           if (i-k+l-j>2*MAXLOOP_L-2) break;
1784           if (abs(i-k-l+j) >= ASS) continue;
1785           type2 = pair[S1[k]][S2[l]];
1786           if (!type2) continue;
1787           LE = E_IntLoop(i-k-1, l-j-1, type2, rtype[type],
1788                          SS1[k+1], SS2[l-1], SS1[i-1], SS2[j+1],P);
1789           if (E == c[k][l]+LE+(i-k+l-j)*penalty) {
1790             traced=1;
1791             i=k; j=l;
1792             *Duplex_El+=LE;
1793             break;
1794           }
1795         }
1796         if (traced) break;
1797       }
1798     }
1799     if (!traced) {
1800       int correction;
1801       correction = E_ExtLoop(type, (i>1) ? SS1[i-1] : -1, (j<n2) ? SS2[j+1] : -1, P);
1802       E-=correction;
1803       *Duplex_El+=correction;
1804       /**
1805       *** if (i>1)    {E -= P->dangle5[type][SS1[i-1]]; *Duplex_El +=P->dangle5[type][SS1[i-1]];}
1806       *** if (j<n2)   {E -= P->dangle3[type][SS2[j+1]]; *Duplex_El +=P->dangle3[type][SS2[j+1]];}
1807       *** if (type>2) {E -= P->TerminalAU;                    *Duplex_El +=P->TerminalAU;}
1808       **/
1809       if (E != P->DuplexInit) {
1810         nrerror("backtrack failed in fold duplex end");
1811       } else break;
1812     }
1813   }
1814 /*   if (i>1)  i--; */
1815 /*   if (j<n2) j++; */
1816   /* struc = (char *) space(i0-i+1+j-j0+1+2); */ /* declare final duplex structure */
1817   struc = (char *) space(i0-i+1+n2-1+1+2); /* declare final duplex structure */
1818   char * struc2;
1819   struc2 = (char *) space(n2+1);
1820   /* char * struct_const; */
1821   for (k=MAX2(i,1); k<=i0; k++) if (!st1[k-1]) st1[k-1] = '.';
1822   /* for (k=j0; k<=j; k++) if (!st2[k-1]) st2[k-1] = struc_loop[k-1];*/ /* '.'; normal */
1823   /*  char * struct_const; */
1824   /*  struct_const = (char *) space(sizeof(char)*(n2+1));   */
1825   for (k=1; k<=n2; k++) {
1826     if (!st2[k-1]) st2[k-1] = struc_loop[k-1];/* '.'; */
1827     struc2[k-1] = st2[k-1];/* '.'; */
1828     /*      if (k>=j0 && k<=j){ */
1829     /*              struct_const[k-1]='x'; */
1830     /*      } */
1831     /*      else{ */
1832     /*              if(k<j0) {struct_const[k-1]='<';} */
1833     /*              if(k>j) {struct_const[k-1]='>';} */
1834     /*      } */
1835   }
1836   char duplexseq_1[j0];
1837   char duplexseq_2[n2-j+2];
1838   if(j<n2){
1839     strncpy(duplexseq_1, snoseq, j0-1);
1840     strcpy(duplexseq_2, snoseq + j);
1841     duplexseq_1[j0-1]='\0';
1842     duplexseq_2[n2-j+1]='\0';
1843     duplexT temp;
1844     temp=duplexfold(duplexseq_1, duplexseq_2);
1845     *Loop_D =  MIN2(0,-410 + (int) 100 * temp.energy);
1846     if(*Loop_D){
1847       int l1,ibegin, iend, jbegin, jend;
1848       l1=strchr(temp.structure, '&')-temp.structure;
1849       ibegin=temp.i-l1;
1850       iend  =temp.i-1;
1851       jbegin=temp.j;
1852       jend  =temp.j+strlen(temp.structure)-l1-2-1;
1853       for(k=ibegin+1; k<=iend+1; k++){
1854         struc2[k-1]=temp.structure[k-ibegin-1];
1855       }
1856       for(k=jbegin+j; k<=jend+j; k++){
1857         struc2[k-1]=temp.structure[l1+k-j-jbegin+1];
1858       }
1859     }
1860     free(temp.structure);
1861   }
1862   strcpy(struc, st1+MAX2(i-1,0)); strcat(struc, "&");
1863   /* strcat(struc, st2); */
1864   strncat(struc, struc2+5, strlen(struc2)-10);
1865   free(struc2);
1866   free(struc_loop);
1867   free(st1); free(st2);
1868   /* free_arrays(); */
1869   return struc;
1870 }
1871
1872 void Lsnoop_subopt_list_XS(const char *s1, const char *s2,  const int **access_s1, int delta, int w,
1873                         const int penalty, const int threshloop,
1874                         const int threshLE, const int threshRE, const int threshDE, const int threshTE,const int threshSE,const int threshD,
1875                         const int distance,
1876                         const int half_stem, const int max_half_stem,
1877                            const int min_s2, const int max_s2, const int min_s1, const int max_s1, const int min_d1, const int min_d2, const int alignment_length, const char *name, const int fullStemEnergy)
1878 {
1879
1880  int min_colonne=INF;
1881  int max_pos;
1882  int max;max=INF;
1883  /* int temp; */
1884  /* int nsubopt=10; */
1885  n1 = (int) strlen(s1);
1886  n2 = (int) strlen(s2);
1887  int *position;
1888  int *position_j;
1889  int min_j_colonne;
1890  int max_pos_j=INF;
1891  position = (int*) space((n1+3)*sizeof(int));
1892  position_j = (int*) space((n1+3)*sizeof(int));
1893
1894  /* int Eminj, Emin_l; */
1895   int i, j;/*  l1, Emin=INF, i_min=0, j_min=0; */
1896   /* char *struc; */
1897   /* snoopT mfe; */
1898   int *indx;
1899   int *mLoop;
1900   int *cLoop;
1901   folden **foldlist, **foldlist_XS;
1902   int Duplex_El, Duplex_Er;
1903   int Loop_D;
1904   /* int u; */
1905   int Loop_E;
1906   Duplex_El=0;Duplex_Er=0;Loop_E=0, Loop_D=0;
1907   snoexport_fold_arrays(&indx, &mLoop, &cLoop, &foldlist, &foldlist_XS);
1908   if ((!P) || (fabs(P->temperature - temperature)>1e-6)) {
1909     snoupdate_fold_params();  if(P) free(P); P = scale_parameters();
1910     make_pair_matrix();
1911   }
1912
1913   lpair = (int **) space(sizeof(int *) * (6));
1914   lc    = (int **) space(sizeof(int *) * (6));
1915   lr    = (int **) space(sizeof(int *) * (6));
1916   for (i=0; i<6; i++) {
1917           lc[i] = (int *) space(sizeof(int) * (n2+1));
1918         lr[i] = (int *) space(sizeof(int) * (n2+1));
1919         lpair[i] = (int *) space(sizeof(int) * (n2+1));
1920         for(j=n2; j>-1; j--){
1921                 lc[i][j]=INF;
1922                 lr[i][j]=INF;
1923                 lpair[i][j]=0;
1924         }
1925   }
1926   encode_seqs(s1, s2);
1927   int lim_maxj=n2-min_d2 ;
1928   int lim_minj=min_d1;
1929   int lim_maxi=n1-5;
1930   for (i=5; i<=lim_maxi; i++) {
1931     int idx=i%5;
1932     int idx_1=(i-1)%5;
1933     int idx_2=(i-2)%5;
1934     int idx_3=(i-3)%5;
1935     int idx_4=(i-4)%5;
1936     int di1, di2, di3, di4;
1937     di1 = access_s1[5][i]   - access_s1[4][i-1];
1938     di2 =access_s1[5][i-1] - access_s1[4][i-2] + di1;
1939     di3 = access_s1[5][i-2] - access_s1[4][i-3] + di2;
1940     di4 = access_s1[5][i-3] - access_s1[4][i-4] + di3;
1941     di1=MIN2(di1,165);
1942     di2=MIN2(di2,330);
1943     di3=MIN2(di3,495);
1944     di4=MIN2(di4,660);
1945     for (j=lim_maxj; j>lim_minj; j--) {
1946       int type, type2;/*  E, k,l; */
1947       type = pair[S1[i]][S2[j]];
1948       lpair[idx][j] = type;
1949       lc[idx][j] = (type) ? P->DuplexInit + access_s1[1][i] : INF;
1950       lr[idx][j] = INF;
1951       if(!type) continue;
1952       if( /*pair[S1[i+1]][S2[j-1]] && Be sure it binds*/
1953           j < max_s1 && j > min_s1 &&
1954           j > n2 - max_s2 - max_half_stem &&
1955           j < n2 -min_s2 -half_stem && S1[i-2]==4 ) { /*constraint on s2 and i*/
1956         int min_k, max_k;
1957         max_k = MIN2(n2-min_s2,j+max_half_stem+1);
1958         min_k = MAX2(j+half_stem+1, n2-max_s2);
1959         folden * temp;
1960         temp=foldlist[j+1];
1961         while(temp->next){
1962           int k = temp->k;
1963           /* if(k >= min_k-1 && k < max_k){ comment to recover normal behaviour */
1964           if(lpair[idx_3][k+1] && lc[idx_3][k+1] /*+di3*/ < 411 /*&& lpair[idx_4][k+2]*/){ /*  remove second condition */
1965             lr[idx][j]=MIN2(lr[idx][j], di3 + lc[idx_3][k+1]+temp->energy);/*--NU--*/
1966           }
1967           /*else*/ if(lpair[idx_4][k+1] && /*di4 +*/ lc[idx_4][k+1] < 411  ){/*--NUN--*/ /*  remove second condition  */
1968             lr[idx][j]=MIN2(lr[idx][j], di4 + lc[idx_4][k+1]+temp->energy);
1969           }
1970             /*  } */
1971           temp=temp->next;
1972         }
1973       }
1974       /* dangle 5'SIDE relative to the mRNA  */
1975       /**
1976       *** lc[idx][j] += P->dangle5[type][SS1[i-1]];
1977       *** lc[idx][j] += P->dangle3[type][SS2[j+1]];
1978       *** if (type>2) lc[idx][j] += P->TerminalAU;
1979       **/
1980       lc[idx][j]+=E_ExtLoop(type, SS1[i-1] ,  SS2[j+1] , P);
1981       /*       if(j<n2 && i>1){ */
1982       /* type2=pair[S1[i-1]][S2[j+1]]; */
1983         type2=lpair[idx_1][j+1];
1984         if(type2>0 ){
1985           lc[idx][j]=MIN2(lc[idx_1][j+1]+E_IntLoop(0,0,type2, rtype[type],SS1[i], SS2[j], SS1[i-1], SS2[j+1], P)+di1, lc[idx][j]);
1986           lr[idx][j]=MIN2(lr[idx_1][j+1]+E_IntLoop(0,0,type2, rtype[type],SS1[i], SS2[j], SS1[i-1], SS2[j+1], P)+di1, lr[idx][j]);
1987         }
1988         type2=lpair[idx_2][j+2];
1989         if(type2>0 ){
1990           lc[idx][j]=MIN2(lc[idx_2][j+2]+E_IntLoop(1,1,type2, rtype[type],SS1[i-1], SS2[j+1], SS1[i-1], SS2[j+1], P)+di2, lc[idx][j]);
1991           lr[idx][j]=MIN2(lr[idx_2][j+2]+E_IntLoop(1,1,type2, rtype[type],SS1[i-1], SS2[j+1], SS1[i-1], SS2[j+1], P)+di2, lr[idx][j]);
1992
1993       }
1994         type2 =lpair[idx_3][j+3];
1995         if(type2>0 ){
1996           lc[idx][j]=MIN2(lc[idx_3][j+3]+E_IntLoop(2,2,type2, rtype[type],SS1[i-2], SS2[j+2], SS1[i-1], SS2[j+1], P)+di3,lc[idx][j]);
1997           lr[idx][j]=MIN2(lr[idx_3][j+3]+E_IntLoop(2,2,type2, rtype[type],SS1[i-2], SS2[j+2], SS1[i-1], SS2[j+1], P)+di3,lr[idx][j]);
1998
1999       }
2000       int bla;
2001       int temp2;
2002       temp2=min_colonne;
2003       bla=lr[idx][j] + E_ExtLoop(rtype[type], SS2[j-1], SS1[i+1] , P);
2004         /**
2005         *** (type>2?P->TerminalAU:0)+P->dangle3[rtype[type]][SS1[i+1]]+P->dangle5[rtype[type]][SS2[j-1]];
2006         **/
2007       min_colonne=MIN2(bla, min_colonne);
2008       if(temp2>min_colonne){
2009         min_j_colonne=j;
2010       }
2011     }
2012     position[i]=min_colonne;
2013     if(max>=min_colonne){
2014       max=min_colonne;
2015       max_pos=i;
2016       max_pos_j=min_j_colonne;
2017       }
2018     position_j[i]=min_j_colonne;
2019     min_colonne=INF;
2020  }
2021   free(S1); free(S2); free(SS1); free(SS2);
2022
2023   if(max<threshTE + 30*alignment_length){
2024     find_max_snoop_XS(s1, s2, access_s1,max,alignment_length, position, position_j,
2025                      delta, distance, penalty, threshloop,
2026                      threshLE, threshRE, threshDE, threshTE, threshSE, threshD,
2027                       half_stem, max_half_stem, min_s2, max_s2, min_s1, max_s1, min_d1, min_d2,name,fullStemEnergy);
2028    }
2029   for (i=1; i<6; i++) {free(lc[i]);free(lr[i]);free(lpair[i]);}
2030   free(lc[0]);free(lr[0]);free(lpair[0]);
2031   free(lc);free(lr);free(lpair);
2032   free(position);free(position_j);
2033 }
2034
2035
2036 PRIVATE void find_max_snoop_XS(const char *s1, const char *s2, const int **access_s1,
2037                                const int max,  const int alignment_length,
2038                                const int* position, const int* position_j, const int delta,
2039                                const int distance, const int penalty, const int threshloop,  const int threshLE, const int threshRE,
2040                                const int threshDE, const int threshTE, const int threshSE, const int threshD,
2041                                const int half_stem, const int max_half_stem, const int min_s2, const int max_s2, const int min_s1, const int max_s1, const int min_d1, const int min_d2, const char *name, const int fullStemEnergy){
2042   int count=0;
2043   int n3=strlen(s1);
2044   int n4=strlen(s2);
2045   int pos=n1-4;
2046   int max_pos_j;
2047   int threshold = MIN2(threshTE + alignment_length*30, -100);
2048   /* printf("threshTE %d max %d\n", threshTE, max); */
2049   /* #pragma omp parallel for */
2050   /*   for(pos=n1+1;pos>distance;pos--){ */
2051   while(pos-->5){
2052     int temp_min=0;
2053     if(position[pos]<(threshold)){
2054       int search_range;
2055       search_range=distance+1;
2056       while(--search_range){
2057         if(position[pos-search_range]<=position[pos-temp_min]){
2058           temp_min=search_range;
2059         }
2060       }
2061       pos-=temp_min;
2062       max_pos_j=position_j[pos];
2063       int begin=MAX2(5, pos-alignment_length);
2064       int end  =MIN2(n3-5, pos-1);
2065       char *s3 = (char*) space(sizeof(char)*(end-begin+2)+5);
2066       strncpy(s3, (s1+begin), end - begin +1);
2067       strcat(s3,"NNNNN\0");
2068
2069       int n5 = strlen(s3);
2070       snoopT test;
2071       test = snoopfold_XS(s3, s2, access_s1, pos, max_pos_j ,penalty,
2072                           threshloop, threshLE, threshRE,
2073                           threshDE, threshD, half_stem,
2074                           max_half_stem, min_s2, max_s2, min_s1,
2075                           max_s1, min_d1, min_d2, fullStemEnergy);
2076       if(test.energy==INF){
2077         free(s3);
2078         continue;
2079       }
2080       if( test.Duplex_El > threshLE * 0.01 ||test.Duplex_Er > threshRE * 0.01  ||
2081          test.Loop_D > threshD * 0.01 || (test.Duplex_Er + test.Duplex_El) > threshDE * 0.01 ||
2082          (test.Duplex_Er + test.Duplex_El + test.Loop_E) > threshTE*0.01 || (test.Duplex_Er + test.Duplex_El + test.Loop_E + test.Loop_D + 410) > threshSE*0.01) {
2083         free(test.structure);free(s3);
2084         continue;
2085       }
2086
2087       char *s4;
2088       s4 = (char*) space(sizeof(char) *(n4-9));
2089       strncpy(s4, s2+5, n4-10);
2090       s4[n4-10]='\0';
2091
2092       char *s5 = space(sizeof(char) * n5-test.i+2-5);
2093       strncpy(s5,s3+test.i-1,n5-test.i+1-5);
2094       s5[n5-test.i+1-5]='\0';
2095       float dE = ((float) (access_s1[n5-test.i+1-5][pos]))*0.01;
2096       printf("%s %3d,%-3d;%3d : %3d,%-3d (%5.2f = %5.2f + %5.2f + %5.2f + %5.2f + %5.2f + 4.10) (%5.2f)\n%s&%s\n" ,
2097              test.structure, pos  -  (n5 - test.i) ,pos - 5, pos - (n5 - test.u ),
2098              max_pos_j-5, max_pos_j-5 + (strrchr(test.structure,'>') - strchr(test.structure,'>')),
2099              test.Loop_D + test.Duplex_El + test.Duplex_Er + test.Loop_E + 4.10+dE, test.Duplex_El,
2100              test.Duplex_Er, test.Loop_E, test.Loop_D,dE ,test.fullStemEnergy, s5,s4);
2101       if(name){
2102         int begin_t, end_t, begin_q, end_q, and, pipe, i;
2103         char psoutput[100];
2104         begin_q=0;
2105         end_q=n4-10;
2106         begin_t=0;
2107         end_t=n5-test.i+ 1-5;
2108         and=end_t+1;
2109         pipe=test.u -test.i + 1;
2110         cut_point = end_t +1 ;
2111         char *catseq, *catstruct;/*  *fname;  */
2112         catseq = (char*) space(n5 + end_q -begin_q +2);
2113         catstruct = (char*) space(n5 + end_q -begin_q +2);
2114         strcpy(catseq, s5);
2115         strncpy(catstruct, test.structure, end_t);
2116         strcat(catseq, s4);
2117         strncat(catstruct, test.structure+end_t+1, end_q-begin_q+1);
2118         catstruct[end_t - begin_t + end_q-begin_q+2]='\0';
2119         catseq[end_t - begin_t + end_q-begin_q+2]='\0';
2120         int *relative_access;
2121         relative_access = space(sizeof(int)*strlen(s5));
2122
2123         relative_access[0] = access_s1[1][pos - (n5  - test.i) + 5];
2124         for(i=1;i<strlen(s5);i++){
2125           relative_access[i] =  access_s1[i+1][pos - (n5  - test.i) + i + 5] -  access_s1[i][pos - (n5  - test.i) + i + 4];
2126         }
2127         char str[16];
2128         char upos[16];
2129         strcpy(psoutput,"sno_XS_");
2130         sprintf(str,"%d",count);
2131         strcat(psoutput,str);
2132         sprintf(upos,"%d",pos - (n5 - test.u ));
2133         strcat(psoutput,"_u_");
2134         strcat(psoutput,upos);
2135         strcat(psoutput,"_");
2136         strcat(psoutput,name);
2137         strcat(psoutput,".ps\0");
2138         PS_rna_plot_snoop_a(catseq, catstruct, psoutput,relative_access,NULL);
2139         free(catseq);free(catstruct);free(relative_access);
2140         count++;
2141       }
2142       free(s3);free(s4);free(s5);free(test.structure);
2143     }
2144   }
2145 }
2146
2147 snoopT snoopfold_XS(const char *s1, const char *s2, const int **access_s1, const int pos_i, const int pos_j,
2148                  const int penalty, const int threshloop, const int threshLE, const int threshRE, const int threshDE,
2149                  const int threshD,
2150                  const int half_stem, const int max_half_stem,
2151                     const int min_s2, const int max_s2, const int min_s1, const int max_s1, const int min_d1, const int min_d2, const int fullStemEnergy) {
2152   /*   int Eminj, Emin_l; */
2153   int a,b,i, j, Emin=INF, a_min=0, b_min=0;
2154   char *struc;
2155   snoopT mfe;
2156   int *indx;
2157   int *mLoop;
2158   int *cLoop;
2159   folden** foldlist, **foldlist_XS;
2160   int Duplex_El, Duplex_Er;
2161   int Loop_D;
2162   int u;
2163   int Loop_E;
2164
2165   Duplex_El=0;Duplex_Er=0;Loop_E=0, Loop_D=0;
2166   snoexport_fold_arrays(&indx, &mLoop, &cLoop,&foldlist, &foldlist_XS );
2167   n1 = (int) strlen(s1);
2168   n2 = (int) strlen(s2);
2169
2170   if ((!P) || (fabs(P->temperature - temperature)>1e-6)) {
2171     snoupdate_fold_params(); if(P) free(P); P = scale_parameters();
2172     make_pair_matrix();
2173   }
2174
2175   c = (int **) space(sizeof(int *) * (n1+1));
2176   r = (int **) space(sizeof(int *) * (n1+1));
2177   for (i=0; i<=n1; i++) {
2178           c[i] = (int *) space(sizeof(int) * (n2+1));
2179         r[i] = (int *) space(sizeof(int) * (n2+1));
2180         for(j=n2; j>-1; j--){
2181                 c[i][j]=INF;
2182                 r[i][j]=INF;
2183         }
2184   }
2185   encode_seqs(s1, s2);
2186   i=n1-5;
2187   j=pos_j;
2188   /* printf("tar: %s\nsno: %s\n ", s1, s2); */
2189   /* printf("pos_i %d pos_j %d\n", pos_i, pos_j); */
2190   /* printf("type %d n1 %d n2 %d S1[n1] %d S2[n2] %d", pair[S1[i]][S2[j]], n1, n2, S1[n1], S2[n2]); */
2191   int type, type2, E, p,q;
2192   r[i][j] = P->DuplexInit;
2193   /* r[i][j] += P->dangle3[rtype[type]][SS1[i+1]] + P->dangle5[rtype[type]][SS2[j-1]];  */
2194
2195   if(pair[S1[i]][S2[j]]>2) r[i][j]+=P->TerminalAU;
2196   for(a=i-1; a>0;a--){ /* i-1 */
2197     r[a+1][0]=INF;
2198     for(b=j+1; b<=n2-min_d2;b++){ /* j+1 */
2199       r[a][b]=INF;
2200       type = pair[S1[a]][S2[b]];
2201        if(!type) continue;
2202        if(S1[a+1]==4){
2203          folden * temp;
2204          temp=foldlist_XS[b-1];
2205          while(temp->next) {
2206            int k = temp->k;
2207            if(pair[S1[a+3]][S2[k-1]] && k< max_s1 && k > min_s1 && k > n2 - max_s2 - max_half_stem &&  k < n2 -min_s2 -half_stem /*&& r[a+3][k-1] + access_s1[i-(a+3)+1][pos_i] < 411*/) { /* remove last condition last condition is to check if the interaction is stable enough */
2208              c[a][b]=MIN2(c[a][b], r[a+3][k-1]+temp->energy);
2209            }
2210            temp=temp->next;
2211          }
2212        }
2213        if(S1[a+2]==4){
2214          folden * temp;
2215          temp=foldlist_XS[b-1];
2216          while(temp->next){
2217            int k = temp->k;
2218            if(pair[S1[a+4]][S2[k-1]] &&  k< max_s1 && k > min_s1 && k > n2 - max_s2 - max_half_stem &&  k < n2 -min_s2 -half_stem /*&& r[a+4][k-1] + access_s1[i-(a+4)+1][pos_i] < 411 */ ) { /* remove last condition  */
2219              c[a][b]=MIN2(c[a][b], r[a+4][k-1]+temp->energy);
2220            }
2221            temp=temp->next;
2222          }
2223        }
2224        for(p=a+1; p<n1 && (p-a) < MAXLOOP_L; p++) { /* p < n1 */
2225          for (q=b-1; q > 1  ; q--) {  /* q > 1 */
2226            if (p-a+b-q>2*MAXLOOP_L-2) break;
2227            if (abs((p-a)-(b-q)) >= ASS ) continue;
2228            type2 = pair[S1[p]][S2[q]];
2229            if (!type2) continue;
2230            E = E_IntLoop(p-a-1, b-q-1, type2, rtype[type],SS1[a+1], SS2[b-1], SS1[p-1], SS2[q+1],P);
2231            c[a][b] = MIN2(c[a][b], c[p][q]+E);
2232            r[a][b] = MIN2(r[a][b], r[p][q]+E);
2233          }
2234        }
2235        E = c[a][b];
2236        if (type>2) E += P->TerminalAU;
2237        /*        E +=P->dangle5[rtype[type]][SS1[i+1]]; */
2238        /* E +=P->dangle5[rtype[type]][SS2[j-1]];  */
2239        E+=access_s1[i-a+1][pos_i];
2240        if (E<Emin) {
2241          Emin=E; a_min=a; b_min=b;
2242        }
2243     }
2244   }
2245   if(Emin > 0){
2246     printf("no target found under the constraints chosen\n");
2247     for (i=0; i<=n1; i++) {free(r[i]);free(c[i]);}
2248     free(c);
2249     free(r);
2250     free(S1); free(S2); free(SS1); free(SS2);
2251     mfe.energy=INF;
2252     return mfe;
2253   }
2254   type2=pair[S1[a_min]][S2[b_min]];
2255   if(type2>2) Emin +=P->TerminalAU;
2256   mfe.energy = ((float) (Emin))/100;
2257   struc = snoop_backtrack_XS(a_min, b_min,s2, &Duplex_El, &Duplex_Er, &Loop_E, &Loop_D,
2258                              &u, penalty, threshloop, threshLE, threshRE,threshDE, threshD,
2259                              half_stem, max_half_stem, min_s2, max_s2, min_s1, max_s1, min_d1, min_d2);
2260
2261   mfe.i = a_min;
2262   mfe.j = b_min;
2263   mfe.u = u;
2264   mfe.Duplex_Er = (float) Duplex_Er/100;
2265   mfe.Duplex_El = (float) Duplex_El/100;
2266   mfe.Loop_D = (float) Loop_D/100;
2267   mfe.Loop_E = (float) Loop_E/100;
2268   mfe.energy = (float) Emin/100 ;
2269   mfe.fullStemEnergy = (float) fullStemEnergy/100;
2270   mfe.structure = struc;
2271   return mfe;
2272 }
2273
2274 PRIVATE char *snoop_backtrack_XS(int i, int j, const char* snoseq,
2275                               int *Duplex_El, int *Duplex_Er,
2276                               int *Loop_E, int *Loop_D, int *u,
2277                               const int penalty, const int threshloop,
2278                               const int threshLE, const int threshRE, const int threshDE, const int threshD,
2279                               const int half_stem, const int max_half_stem,
2280                               const int min_s2, const int max_s2, const int min_s1,
2281                               const int max_s1, const int min_d1, const int min_d2) {
2282   /* backtrack structure going backwards from i, and forwards from j
2283      return structure in bracket notation with & as separator */
2284   int k, l, type, type2, E, traced, i0, j0;
2285   int traced_c=0; /* flag for following backtrack in c or r */
2286   char *st1, *st2, *struc;
2287   char *struc_loop;
2288
2289   st1 = (char *) space(sizeof(char)*(n1+1));
2290   st2 = (char *) space(sizeof(char)*(n2+1));
2291   int *indx;
2292   int *mLoop;
2293   int *cLoop;
2294   folden **foldlist, **foldlist_XS;
2295   type=pair[S1[i]][S2[j]];
2296   snoexport_fold_arrays(&indx, &mLoop, &cLoop,&foldlist, &foldlist_XS );
2297   i0=i;j0=j;
2298   /*   i0=MAX2(i,1); j0=MIN2(j+1,n2); */
2299   while (i<=n1 && j>=1 ) {
2300     if(!traced_c) {
2301       E = c[i][j]; traced=0;
2302       st1[i] = '<';
2303       st2[j-1] = '>';
2304       type = pair[S1[i]][S2[j]];
2305       if (!type) nrerror("backtrack failed in fold duplex c");
2306       for (k=i+1; k>0 && (k-i)<MAXLOOP_L; k++) {
2307         for (l=j-1; l>=1 ; l--) {
2308           int LE;
2309           if (k-i+j-l>2*MAXLOOP_L-2) break;
2310           if (abs(k-i-j+l) >= ASS) continue;
2311           type2 = pair[S1[k]][S2[l]];
2312           if (!type2) continue;
2313           LE = E_IntLoop(k-i-1, j-l-1, type2, rtype[type],
2314                          SS1[i+1], SS2[j-1], SS1[k-1], SS2[l+1],P);
2315           if (E == c[k][l]+LE) {
2316             traced=1;
2317             i=k; j=l;
2318             *Duplex_El+=LE;
2319             break;
2320           }
2321         }
2322         if (traced) break;
2323       }
2324       if(!traced){
2325         if(S1[i+1]==4) {
2326           folden * temp;
2327           temp=foldlist_XS[j-1];
2328           while(temp->next) {
2329             int k = temp->k;
2330             if(pair[S1[i+3]][S2[k-1]] && k< max_s1 && k > min_s1 && k > n2 - max_s2 - max_half_stem &&  k < n2 -min_s2 -half_stem ) {
2331               if(E==r[i+3][k-1]+temp->energy){
2332                 *Loop_E=temp->energy;
2333                 st1[i+1]= '|';
2334                 st1[i+2]='.';
2335                 *u=i+1;
2336                 int a,b;
2337                 for(a=0; a< MISMATCH ;a++){
2338                   for(b=0; b< MISMATCH ; b++){
2339                     int ij=indx[j-1-a]+k+b;
2340                     if(cLoop[ij]==temp->energy) {
2341                       struc_loop=snobacktrack_fold_from_pair(snoseq, k+b, j-1-a);
2342                       a=INF; b=INF;
2343                     }
2344                   }
2345                 }
2346                 traced=1;
2347                 traced_c=1;
2348                 i=i+3;j=k-1;
2349                 break;
2350               }
2351             }
2352             temp=temp->next;
2353           }
2354         }
2355         if (S1[i+2]==4){  /* introduce structure from RNAfold */
2356           folden * temp;
2357           temp=foldlist_XS[j-1];
2358           while(temp->next) {
2359             int k = temp->k;
2360             if(pair[S1[i+4]][S2[k-1]] && k< max_s1 && k > min_s1 && k > n2 - max_s2 - max_half_stem &&  k < n2 -min_s2 -half_stem ) {
2361               if(E==r[i+4][k-1]+temp->energy){
2362                 *Loop_E=temp->energy;
2363                 st1[i+2]= '|';
2364                 st1[i+1]=st1[i+3]='.';
2365                 *u=i+2;
2366                 int a,b;
2367                 for(a=0; a< MISMATCH ;a++){
2368                   for(b=0; b< MISMATCH ; b++){
2369                     int ij=indx[j-1-a]+k+b;
2370                     if(cLoop[ij]==temp->energy) {
2371                       struc_loop=snobacktrack_fold_from_pair(snoseq, k+b, j-a-1);
2372                       a=INF; b=INF;
2373                     }
2374                   }
2375                 }
2376                 traced=1;
2377                 traced_c=1;
2378                 i=i+4;j=k-1;
2379                 break;
2380               }
2381             }
2382             temp=temp->next;
2383           }
2384         }
2385       }/* traced? */
2386     }/* traced_r? */
2387     else{
2388       E = r[i][j]; traced=0;
2389       st1[i] = '<';
2390       st2[j-1] = '>';
2391       type = pair[S1[i]][S2[j]];
2392       if (!type) nrerror("backtrack failed in fold duplex r");
2393       for (k=i+1; k>0 && (k-i)<MAXLOOP_L; k++) {
2394         for (l=j-1; l>=1 ; l--) {
2395           int LE;
2396           if (k-i+j-l>2*MAXLOOP_L-2) break;
2397           if (abs(k-i-j+l) >= ASS) continue;
2398           type2 = pair[S1[k]][S2[l]];
2399           if (!type2) continue;
2400           LE = E_IntLoop(k-i-1, j-l-1, type2, rtype[type],
2401                          SS1[i+1], SS2[j-1], SS1[k-1], SS2[l+1],P);
2402           if (E == r[k][l]+LE) {
2403             traced=1;
2404             i=k; j=l;
2405             *Duplex_Er+=LE;
2406             break;
2407           }
2408         }
2409         if (traced) break;
2410       }
2411     }
2412     if (!traced) {
2413 /*       if (i>1)    {E -= P->dangle5[type][SS1[i-1]]; *Duplex_El +=P->dangle5[type][SS1[i-1]];} */
2414 /*       if (j<n2)   {E -= P->dangle3[type][SS2[j+1]]; *Duplex_El +=P->dangle3[type][SS2[j+1]];} */
2415       if (type>2) {E -= P->TerminalAU;        *Duplex_Er +=P->TerminalAU;}
2416       if (E != P->DuplexInit) {
2417         nrerror("backtrack failed in fold duplex end");
2418       } else break;
2419     }
2420   }
2421
2422
2423   /* struc = (char *) space(i0-i+1+j-j0+1+2); */ /* declare final duplex structure */
2424   struc = (char *) space(i-i0+1+n2); /* declare final duplex structure */
2425   char * struc2;
2426   struc2 = (char *) space(n2+1);
2427   /* char * struct_const; */
2428
2429   for (k=MIN2(i0,1); k<=i; k++) if (!st1[k-1]) st1[k-1] = '.';
2430   /* for (k=j0; k<=j; k++) if (!st2[k-1]) st2[k-1] = struc_loop[k-1];*/ /* '.'; normal */
2431   /*  char * struct_const; */
2432   /*  struct_const = (char *) space(sizeof(char)*(n2+1));   */
2433   for (k=1; k<=n2; k++) {
2434     if (!st2[k-1]) st2[k-1] = struc_loop[k-1];/* '.'; */
2435     struc2[k-1] = st2[k-1];/* '.'; */
2436     /*      if (k>=j0 && k<=j){ */
2437     /*              struct_const[k-1]='x'; */
2438     /*      } */
2439     /*      else{ */
2440     /*              if(k<j0) {struct_const[k-1]='<';} */
2441     /*              if(k>j) {struct_const[k-1]='>';} */
2442     /*      } */
2443   }
2444   char duplexseq_1[j];
2445   char duplexseq_2[n2-j0+2];
2446   if(j0<n2){
2447     strncpy(duplexseq_1, snoseq, j-1);
2448     strcpy(duplexseq_2, snoseq + j0);
2449     duplexseq_1[j-1]='\0';
2450     duplexseq_2[n2-j0+1]='\0';
2451     duplexT temp;
2452     temp=duplexfold(duplexseq_1, duplexseq_2);
2453     *Loop_D =  MIN2(0,-410 + (int) 100 * temp.energy);
2454     if(*Loop_D){
2455       int l1,ibegin, iend, jbegin, jend;
2456       l1=strchr(temp.structure, '&')-temp.structure;
2457       ibegin=temp.i-l1;
2458       iend  =temp.i-1;
2459       jbegin=temp.j;
2460       jend  =temp.j+strlen(temp.structure)-l1-2-1;
2461       for(k=ibegin+1; k<=iend+1; k++){
2462         struc2[k-1]=temp.structure[k-ibegin-1];
2463       }
2464       for(k=jbegin+j0; k<=jend+j0; k++){
2465         struc2[k-1]=temp.structure[l1+k-j0-jbegin+1];
2466       }
2467     }
2468     free(temp.structure);
2469   }
2470   strcpy(struc, st1+MAX2(i0,1)); strcat(struc, "&");
2471   /* strcat(struc, st2); */
2472   strncat(struc, struc2+5, strlen(struc2)-10);
2473   free(struc2);
2474   free(struc_loop);
2475   free(st1); free(st2);
2476
2477     for (i=0; i<=n1; i++) {free(r[i]);free(c[i]);}
2478     free(c);
2479     free(r);
2480     free(S1);free(S2);free(SS1);free(SS2);
2481     /* free_arrays(); */
2482   return struc;
2483 }
2484
2485 PRIVATE int covscore(const int *types, int n_seq) {
2486   /* calculate co-variance bonus for a pair depending on  */
2487   /* compensatory/consistent mutations and incompatible seqs */
2488   /* should be 0 for conserved pairs, >0 for good pairs      */
2489 #define NONE -10000 /* score for forbidden pairs */
2490   int k,l,s,score, pscore;
2491   int dm[7][7]={{0,0,0,0,0,0,0}, /* hamming distance between pairs */
2492                 {0,0,2,2,1,2,2} /* CG */,
2493                 {0,2,0,1,2,2,2} /* GC */,
2494                 {0,2,1,0,2,1,2} /* GU */,
2495                 {0,1,2,2,0,2,1} /* UG */,
2496                 {0,2,2,1,2,0,2} /* AU */,
2497                 {0,2,2,2,1,2,0} /* UA */};
2498
2499   int pfreq[8]={0,0,0,0,0,0,0,0};
2500   for (s=0; s<n_seq; s++)
2501     pfreq[types[s]]++;
2502
2503   if (pfreq[0]*2>n_seq) return NONE;
2504   for (k=1,score=0; k<=6; k++) /* ignore pairtype 7 (gap-gap) */
2505     for (l=k+1; l<=6; l++)
2506       /* scores for replacements between pairtypes    */
2507       /* consistent or compensatory mutations score 1 or 2  */
2508       score += pfreq[k]*pfreq[l]*dm[k][l];
2509
2510   /* counter examples score -1, gap-gap scores -0.25   */
2511   pscore = cv_fact *
2512     ((UNIT*score)/n_seq - nc_fact*UNIT*(pfreq[0] + pfreq[7]*0.25));
2513   return pscore;
2514 }
2515
2516 /*---------------------------------------------------------------------------*/
2517
2518 PRIVATE short * aliencode_seq(const char *sequence) {
2519   unsigned int i,l;
2520   short *Stemp;
2521   l = strlen(sequence);
2522   Stemp = (short *) space(sizeof(short)*(l+2));
2523   Stemp[0] = (short) l;
2524
2525   /* make numerical encoding of sequence */
2526   for (i=1; i<=l; i++)
2527     Stemp[i]= (short) encode_char(toupper(sequence[i-1]));
2528
2529   /* for circular folding add first base at position n+1 */
2530   /* Stemp[l+1] = Stemp[1]; */
2531
2532   return Stemp;
2533 }
2534
2535 PRIVATE short * encode_seq(const char *sequence) {
2536   unsigned int i,l;
2537   short *S;
2538   l = strlen(sequence);
2539 extern double nc_fact;
2540   S = (short *) space(sizeof(short)*(l+2));
2541   S[0] = (short) l;
2542
2543   /* make numerical encoding of sequence */
2544   for (i=1; i<=l; i++)
2545     S[i]= (short) encode_char(toupper(sequence[i-1]));
2546   /* for circular folding add first base at position n+1 */
2547   S[l+1] = S[1];
2548
2549   return S;
2550 }
2551
2552 PRIVATE void encode_seqs(const char *s1, const char *s2) {
2553   unsigned int i,l;
2554
2555   l = strlen(s1);
2556   S1 = encode_seq(s1);
2557   SS1= (short *) space(sizeof(short)*(l+1));
2558   /* SS1 exists only for the special X K and I bases and energy_set!=0 */
2559
2560   for (i=1; i<=l; i++) { /* make numerical encoding of sequence */
2561     SS1[i] = alias[S1[i]];   /* for mismatches of nostandard bases */
2562   }
2563
2564   l = strlen(s2);
2565   S2 = encode_seq(s2);
2566   SS2= (short *) space(sizeof(short)*(l+1));
2567   /* SS2 exists only for the special X K and I bases and energy_set!=0 */
2568
2569   for (i=1; i<=l; i++) { /* make numerical encoding of sequence */
2570     SS2[i] = alias[S2[i]];   /* for mismatches of nostandard bases */
2571   }
2572 }
2573
2574 PRIVATE int compare(const void *sub1, const void *sub2) {
2575   int d;
2576   if (((snoopT *) sub1)->energy > ((snoopT *) sub2)->energy)
2577     return 1;
2578   if (((snoopT *) sub1)->energy < ((snoopT *) sub2)->energy)
2579     return -1;
2580   d = ((snoopT *) sub1)->i - ((snoopT *) sub2)->i;
2581   if (d!=0) return d;
2582   return  ((snoopT *) sub1)->j - ((snoopT *) sub2)->j;
2583 }
2584
2585
2586