/* Last changed Time-stamp: <2007-08-26 11:59:45 ivo> */
/*                
	   compute the duplex structure of two RNA strands,
		allowing only inter-strand base pairs.
	 see cofold() for computing hybrid structures without
			     restriction.

			     Ivo Hofacker
			  Vienna RNA package
*/

#include <config.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <ctype.h>
#include <string.h>
#include "utils.h"
#include "energy_par.h"
#include "fold_vars.h"
#include "snofold.h"
#include "pair_mat.h"
#include "params.h"
#include "snoop.h"
#include "PS_dot.h"
//#include "fold.h"
#include "duplex.h"
#include "loop_energies.h"
/*@unused@*/
static char rcsid[] UNUSED = "$Id: duplex.c,v 1.8 2007/08/26 10:08:44 ivo Exp $";
#define UNIT 100
#define MINPSCORE -2*UNIT
#define PUBLIC
#define PRIVATE static

#define STACK_BULGE1  1   /* stacking energies for bulges of size 1 */
#define NEW_NINIO     1   /* new asymetry penalty */



PRIVATE void  encode_seqs(const char *s1, const char *s2);
PRIVATE short *encode_seq(const char *seq);

PRIVATE void find_max_snoop(const char *s1, const char *s2, const int max, 
			    const int alignment_length, const int* position, 
			    const int delta, const int distance,  const int penalty, 
			    const int threshloop, const int threshLE, const int threshRE, 
			    const int threshDE, const int threshTE, const int threshSE, const int threshD,
			    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);

PRIVATE void find_max_snoop_XS(const char *s1, const char *s2, const int **access_s1, const int max, 
			       const int alignment_length, const int* position, const int *position_j,
			       const int delta, const int distance,  const int penalty, 
			       const int threshloop, const int threshLE, const int threshRE, 
			       const int threshDE, const int threshTE, const int threshSE, const int threshD,
			       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);





PRIVATE char * alisnoop_backtrack(int i, int j, const char ** s2, 
				  int* Duplex_El, int* Duplex_Er, int* Loop_E, int *Loop_D, int *u, 
				  int *pscd, int *psct, int *pscg,
				  const int penalty, const int threshloop, 
				  const int threshLE, const int threshRE, const int threshDE, const int threshD,
				  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 short **S1, const short **S2);

   
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, 
			       const int penalty, const int threshloop, const int threshLE, const int threshRE, const int threshDE, 
			       const int threshD,
			       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);

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, 
			       const int penalty, const int threshloop, const int threshLE, const int threshRE, const int threshDE, 
			       const int threshD,
			       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);




PRIVATE int compare(const void *sub1, const void *sub2);
PRIVATE int covscore(const int *types, int n_seq);
PRIVATE short * aliencode_seq(const char *sequence);

PUBLIC  int snoop_subopt_sorted=0; /* from subopt.c, default 0 */


/*@unused@*/




#define MAXLOOP_L	3
#define MIN2(A, B)      ((A) < (B) ? (A) : (B))
#define MAX2(A, B)      ((A) > (B) ? (A) : (B))
#define ASS		1
PRIVATE paramT *P = NULL;

PRIVATE int   **c = NULL;      /* energy array, given that i-j pair */
PRIVATE int   **r = NULL;
PRIVATE int   **lc = NULL;      /* energy array, given that i-j pair */
PRIVATE int   **lr = NULL;
PRIVATE int   **c_fill = NULL;
PRIVATE int   **r_fill = NULL;
PRIVATE int   **lpair = NULL;


PRIVATE short  *S1 = NULL, *SS1 = NULL, *S2 = NULL, *SS2 = NULL;
PRIVATE short *S1_fill = NULL, *SS1_fill = NULL, *S2_fill = NULL, *SS2_fill = NULL;
PRIVATE int   n1,n2;    /* sequence lengths */

extern int cut_point;

PRIVATE int delay_free=0;
/*--------------------------------------------------------------------------*/

snoopT alisnoopfold(const char **s1, const char **s2, 
		    const int penalty, const int threshloop, 
		    const int threshLE, const int threshRE, const int threshDE, const int threshD,
		    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) {
  
  int s,n_seq;
  int i, j, E, l1,Emin=INF, i_min=0, j_min=0;
  char *struc;
  snoopT mfe;
  int *indx;
  int *mLoop;
  int *cLoop;
  folden  **foldlist; folden **foldlist_XS;
  int Duplex_El, Duplex_Er,pscd,psct,pscg;
  int Loop_D;
  int u;
  int Loop_E;
  short **Sali1,**Sali2;
  int *type,*type2,*type3;
  Duplex_El=0;Duplex_Er=0;Loop_E=0; Loop_D=0;pscd=0;psct=0;pscg=0;
  snoexport_fold_arrays(&indx, &mLoop, &cLoop,&foldlist, &foldlist_XS); 
  n1 = (int) strlen(s1[0]);
  n2 = (int) strlen(s2[0]);
  
  for (s=0; s1[s]!=NULL; s++);
  n_seq = s;
  for (s=0; s2[s]!=NULL; s++);
  if (n_seq != s) nrerror("unequal number of sequences in aliduplexfold()\n");
  
  if ((!P) || (fabs(P->temperature - temperature)>1e-6)) {
    snoupdate_fold_params();  if(P) free(P); P = scale_parameters();
    make_pair_matrix();
  }
  
  c = (int **) space(sizeof(int *) * (n1+1));
  r = (int **) space(sizeof(int *) * (n1+1));
  for (i=0; i<=n1; i++) {
  	c[i] = (int *) space(sizeof(int) * (n2+1));
	r[i] = (int *) space(sizeof(int) * (n2+1));
	for(j=n2; j>-1; j--){
		c[i][j]=INF;
		r[i][j]=INF;
	}
  }
  Sali1 = (short **) space((n_seq+1)*sizeof(short *));
  Sali2 = (short **) space((n_seq+1)*sizeof(short *));
  for (s=0; s<n_seq; s++) {
    if (strlen(s1[s]) != n1) nrerror("uneqal seqence lengths");
    if (strlen(s2[s]) != n2) nrerror("uneqal seqence lengths");
    Sali1[s] = aliencode_seq(s1[s]);
    Sali2[s] = aliencode_seq(s2[s]);
  }
  type = (int *) space(n_seq*sizeof(int));
  type2 = (int *) space(n_seq*sizeof(int));
  type3 = (int *) space(n_seq*sizeof(int));
  //  encode_seqs(s1, s2);
  for (i=6; i<=n1-5; i++) {
    int U; U=0;
    for (s=0; s<n_seq; s++) {
      U+=Sali1[s][i-2];
    }
    U = (U==(n_seq)*4?1:0);
    for (j=n2-min_d2; j>min_d1; j--) {
      int type4, k,l,psc,psc2,psc3;
      for (s=0; s<n_seq; s++) {
	type[s] = pair[Sali1[s][i]][Sali2[s][j]];
      }
      psc = covscore(type, n_seq);
      for (s=0; s<n_seq; s++) if (type[s]==0) type[s]=7;
      c[i][j] = (psc>=MINPSCORE) ? (n_seq*P->DuplexInit) : INF;
      if (psc<MINPSCORE) continue;
      if(/*  pair[Sali1[i+1]][Sali2[j-1]] &&  */
	 U && j < max_s1 && j > min_s1 &&  
	 j > n2 - max_s2 - max_half_stem && 
	 j < n2 -min_s2 -half_stem ) { /*constraint on s2 and i*/
	folden *temp;
	temp=foldlist[j+1];
	while(temp->next){
	  int k = temp->k;
	  for (s=0; s<n_seq; s++) {
	    type2[s]= pair[Sali1[s][i-3]][Sali2[s][k+1]];
	    type3[s]= pair[Sali1[s][i-4]][Sali2[s][k+1]];
	  }
	  psc2 = covscore(type2, n_seq);
	  psc3 = covscore(type3, n_seq);
	  if(psc2 > MINPSCORE){
	    r[i][j]=MIN2(r[i][j],c[i-3][k+1]+temp->energy);
	  }
	  if(psc3 > MINPSCORE){
	    r[i][j]=MIN2(r[i][j],c[i-4][k+1]+temp->energy);
	  }
	  temp=temp->next;
	}
      }
      //dangle 5'SIDE relative to the mRNA 
      for (s=0; s<n_seq; s++) {
	c[i][j] += E_ExtLoop(type[s], Sali1[s][i-1],Sali2[s][j+1],P);
      }
      for (k=i-1; k>0 && (i-k)<MAXLOOP_L; k--) {
	for (l=j+1; l<=n2 ; l++) {
	  if (i-k+l-j>2*MAXLOOP_L-2) break;
	  if (abs(i-k-l+j) >= ASS ) continue;
	  for (E=s=0; s<n_seq; s++) { 
	    type4 = pair[Sali1[s][k]][Sali2[s][l]];
	    if (type4==0) type4=7;
	    E += 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);
	  }
	  c[i][j] = MIN2(c[i][j], c[k][l] + E);
	  r[i][j] = MIN2(r[i][j], r[k][l] + E);
	}
      }
      c[i][j]-=psc;
      r[i][j]-=psc;
      E = r[i][j]; 
      for (s=0; s<n_seq; s++) {
	E+= E_ExtLoop(rtype[type[s]], Sali2[s][j-1], Sali1[s][i+1], P);
	/**
	*** if (i<n1) E += P->dangle3[rtype[type[s]]][Sali1[s][i+1]];
	*** if (j>1)  E += P->dangle5[rtype[type[s]]][Sali2[s][j-1]];
	*** if (type[s]>2) E += P->TerminalAU;
	**/
      }
      if (E<Emin) {
	Emin=E; i_min=i; j_min=j;
      } 
    }
  }
  if(Emin > 0){
  	printf("no target found under the constraints chosen\n");
	for (i=0; i<=n1; i++) {free(r[i]);free(c[i]);}
	free(c);
	free(r);
	for(s=0; s<n_seq;s++){
	  free(Sali1[s]);
	  free(Sali2[s]);
	}
	free(Sali1); free(Sali2);
	free(S2); free(SS1); free(SS2);free(type);free(type2);free(type3);
	mfe.energy=INF;
	mfe.structure=NULL;
	return mfe;
  }
  struc = alisnoop_backtrack(i_min, j_min,(const char**) s2, 
			     &Duplex_El, &Duplex_Er, &Loop_E, 
			     &Loop_D, &u, &pscd, &psct, &pscg,
			     penalty, threshloop, threshLE, 
			     threshRE,threshDE, threshD,
			     half_stem, max_half_stem, min_s2, 
			     max_s2, min_s1, max_s1, min_d1, 
			     min_d2,(const short**) Sali1,(const short**) Sali2);
  //if (i_min<n1-5) i_min++;
  //if (j_min>6 ) j_min--;
  l1 = strchr(struc, '&')-struc;
  mfe.i = i_min-5;
  mfe.j = j_min-5;
  mfe.u = u -5;
  mfe.Duplex_Er = (float) Duplex_Er/100;
  mfe.Duplex_El = (float) Duplex_El/100;
  mfe.Loop_D = (float) Loop_D/100;
  mfe.Loop_E = (float) Loop_E/100;
  mfe.energy = (float) Emin/100 ;
  mfe.pscd = pscd;
  mfe.psct = psct;
  mfe.structure = struc;
  for(s=0; s<n_seq;s++){
    free(Sali1[s]);free(Sali2[s]);
  }
  free(Sali1);free(Sali2);free(type);free(type2);free(type3);

  if (!delay_free) {
    for (i=0; i<=n1; i++) {free(r[i]);free(c[i]);}
    free(c);
    free(r);
    free(S2); free(SS1); free(SS2);
  }
  return mfe;
}

PUBLIC snoopT *alisnoop_subopt(const char **s1, const char **s2, int delta, int w, 
			      const int penalty, const int threshloop, 
			       const int threshLE, const int threshRE, const int threshDE, const int threshTE, const int threshSE, const int threshD,
			      const int distance, 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) {


  short **Sali1, **Sali2;
  //printf("%d %d\n", min_s2, max_s2);
  int i,j,s,n_seq, n1, n2, E, n_subopt=0, n_max;
  char *struc;
  snoopT mfe;
  snoopT *subopt;
  int thresh;
  int *type;
  int Duplex_El, Duplex_Er, Loop_E,pscd,psct,pscg;
  int Loop_D;
  Duplex_El=0; Duplex_Er=0; Loop_E=0;Loop_D=0;pscd=0;psct=0;pscg=0;
  int u;
  u=0;
  n_max=16;
  subopt = (snoopT *) space(n_max*sizeof(snoopT));
  delay_free=1;
  mfe = alisnoopfold(s1, s2, penalty, threshloop, threshLE, threshRE, threshDE,threshD,
                  half_stem, max_half_stem,
		  min_s2, max_s2, min_s1, max_s1, min_d1, min_d2);
  if(mfe.energy > 0){
  	free(subopt);
	delay_free=0;
	return NULL;
  }
  thresh = MIN2((int) ((mfe.Duplex_Er + mfe.Duplex_El + mfe.Loop_E)*100+0.1 + 410) + delta, threshTE );
 //subopt[n_subopt++]=mfe;
  free(mfe.structure);
  n1 = strlen(s1[0]);
  n2 = strlen(s2[0]);
  for (s=0; s1[s]!=NULL; s++);
  n_seq = s;
  Sali1 = (short **) space((n_seq+1)*sizeof(short *));
  Sali2 = (short **) space((n_seq+1)*sizeof(short *));
  for (s=0; s<n_seq; s++) {
    if (strlen(s1[s]) != n1) nrerror("uneqal seqence lengths");
    if (strlen(s2[s]) != n2) nrerror("uneqal seqence lengths");
    Sali1[s] = aliencode_seq(s1[s]);
    Sali2[s] = aliencode_seq(s2[s]);
  }
  Sali1[n_seq]=NULL;  Sali2[n_seq]=NULL;
  type = (int *) space(n_seq*sizeof(int));
  for (i=n1; i>1; i--){
    for (j=1; j<=n2; j++) {
      int  ii,jj, Ed,psc,skip;
      for (s=0; s<n_seq; s++) {
        type[s] = pair[Sali2[s][j]][Sali1[s][i]];
      }
      psc = covscore(type, n_seq);
      for (s=0; s<n_seq; s++) if (type[s]==0) type[s]=7;
      if (psc<MINPSCORE) continue;
      E = Ed = r[i][j];
      for  (s=0; s<n_seq; s++) {
	//	if (i<n1-5) Ed += P->dangle3[type[s]][Sali1[s][i+1]];
	//      if (j>6)  Ed += P->dangle5[type[s]][Sali2[s][j-1]];
	if (type[s]>2) Ed += P->TerminalAU;
      }
      if (Ed>thresh) continue;
      /* too keep output small, remove hits that are dominated by a
	 better one close (w) by. For simplicity we do test without
	 adding dangles, which is slightly inaccurate. 
      */ 
      w=1;
      for (skip=0, ii=MAX2(i-w,1); (ii<=MIN2(i+w,n1)) && type; ii++) { 
        for (jj=MAX2(j-w,1); jj<=MIN2(j+w,n2); jj++)
          if (r[ii][jj]<E) {skip=1; break;}
      }
      if (skip){continue;}
      psct=0;
      pscg=0;
      struc = alisnoop_backtrack(i,j,s2, &Duplex_El, 
				 &Duplex_Er, &Loop_E, &Loop_D, &u, &pscd, &psct,&pscg, 
				 penalty, threshloop,threshLE,threshRE,threshDE, threshD,
				 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);
	      
      if (Duplex_Er > threshRE || Duplex_El > threshLE || Loop_D > threshD ||
         (Duplex_Er + Duplex_El) > threshDE || 
	 (Duplex_Er + Duplex_El + Loop_E) > threshTE ||
	 (Duplex_Er + Duplex_El + Loop_E + Loop_D + 410) > threshSE) {
	 	//printf(" Duplex_Er %d threshRE %d Duplex_El %d threshLE %d \n"
		//       " Duplex_Er + Duplex_El %d  threshDE %d \n"
		//       " Duplex_Er + Duplex_El + Loop_E %d  threshTE %d \n"
		//       " Duplex_Er + Duplex_El + Loop_E + Loop_D %d  threshSE %d \n", 
		//	 Duplex_Er , threshRE , Duplex_El ,threshLE,
	        //         Duplex_Er + Duplex_El, threshDE,
        	//         Duplex_Er + Duplex_El+  Loop_E , threshTE,
                //	 Duplex_Er + Duplex_El+  Loop_E + Loop_D, threshSE); 
	 	Duplex_Er=0; 
		Duplex_El=0;
		Loop_E = 0;
		Loop_D = 0;
		u=0,
		free(struc);
		continue;
	}

      if (n_subopt+1>=n_max) {
	n_max *= 2;
	subopt = (snoopT *) xrealloc(subopt, n_max*sizeof(snoopT));
      }
      
      subopt[n_subopt].i = i-5;
      subopt[n_subopt].j = j-5;
      subopt[n_subopt].u = u-5;
      subopt[n_subopt].Duplex_Er = Duplex_Er * 0.01;
      subopt[n_subopt].Duplex_El = Duplex_El * 0.01;
      subopt[n_subopt].Loop_E = Loop_E * 0.01;
      subopt[n_subopt].Loop_D = Loop_D * 0.01;
      subopt[n_subopt].energy = (Duplex_Er +Duplex_El + Loop_E + Loop_D + 410) * 0.01 ;
      subopt[n_subopt].pscd = pscd * 0.01;
      subopt[n_subopt].psct = -psct * 0.01;
      subopt[n_subopt++].structure = struc;
      
      //i=u;
      Duplex_Er=0; Duplex_El=0; Loop_E=0; Loop_D=0;u=0;pscd=0;psct=0;
    }
  }
  
  for (i=0; i<=n1; i++) {free(c[i]);free(r[i]);}
  free(c);free(r);
  for (s=0; s<n_seq; s++) {
    free(Sali1[s]); free(Sali2[s]);
  }
  free(Sali1); free(Sali2); free(type);
  
  if (snoop_subopt_sorted) 
    qsort(subopt, n_subopt, sizeof(snoopT), compare);
  subopt[n_subopt].i =0;
  subopt[n_subopt].j =0;
  subopt[n_subopt].structure = NULL;
  return subopt;
}







PRIVATE char *alisnoop_backtrack(int i, int j, const char ** snoseq, int *Duplex_El, 
				 int *Duplex_Er, int *Loop_E, int *Loop_D, int *u, 
				 int *pscd, int *psct, int *pscg,
				 const int penalty, const int threshloop, const int threshLE, 
				 const int threshRE, const int threshDE, const int threshD, 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 short **Sali1, const short **Sali2) {
  /* backtrack structure going backwards from i, and forwards from j 
     return structure in bracket notation with & as separator */
  int k, l, *type,*type2,*type3,type4, E, traced, i0, j0,s,n_seq,psc;
  int traced_r=0; //flag for following backtrack in c or r
  char *st1, *st2, *struc;
  char *struc_loop;
  n1 = (int) Sali1[0][0];
  n2 = (int) Sali2[0][0];
  
  for (s=0; Sali1[s]!=NULL; s++);
  n_seq = s;
  for (s=0; Sali2[s]!=NULL; s++);
  if (n_seq != s) nrerror("unequal number of sequences in alibacktrack()\n");
  
  st1 = (char *) space(sizeof(char)*(n1+1));
  st2 = (char *) space(sizeof(char)*(n2+1));
  type = (int *) space(n_seq*sizeof(int));
  type2 = (int *) space(n_seq*sizeof(int));
  type3 = (int *) space(n_seq*sizeof(int));
  int *indx;
  int *mLoop;
  int *cLoop;
  folden **foldlist, **foldlist_XS;
  snoexport_fold_arrays(&indx, &mLoop, &cLoop,&foldlist, &foldlist_XS ); 
  i0=i; j0=j; //MIN2(i+1,n1); j0=MAX2(j-1,1);!modified
  for (s=0; s<n_seq; s++) {
    type[s] = pair[Sali1[s][i]][Sali2[s][j]];
    if(type[s]==0) type[s] = 7;
    *Duplex_Er += E_ExtLoop(rtype[type[s]], (j>1) ? Sali2[s][j-1] : -1, (i<n1) ? Sali1[s][i+1] : -1, P);
    /**
    *** if (i<n1)   *Duplex_Er += P->dangle3[rtype[type[s]]][Sali1[s][i+1]];
    *** if (j>1)    *Duplex_Er += P->dangle5[rtype[type[s]]][Sali2[s][j-1]];
    *** if (type[s]>2) *Duplex_Er += P->TerminalAU;
    **/
  }
  while (i>0 && j<=n2-min_d2 ) {
    if(!traced_r) {
      E = r[i][j]; traced=0;
      st1[i-1] = '<';
      st2[j-1] = '>'; 
      for (s=0; s<n_seq; s++) {
	type[s] = pair[Sali1[s][i]][Sali2[s][j]];
      }
      psc = covscore(type,n_seq);
      for (s=0; s<n_seq; s++) if (type[s]==0) type[s] = 7;
      E += psc;
      *pscd +=psc;
      for (k=i-1; k>0 && (i-k)<MAXLOOP_L; k--) {
	for (l=j+1; l<=n2 ; l++) {
	  int LE;
	  if (i-k+l-j>2*MAXLOOP_L-2) break;
	  if (abs(i-k-l+j) >= ASS) continue;
	  for (s=LE=0; s<n_seq; s++) {
	    type4 = pair[Sali1[s][k]][Sali2[s][l]];
	    if (type4==0) type4=7;
	    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);
	  }
	  if (E == r[k][l]+LE) {
	    traced=1; 
	    i=k; j=l;
	    *Duplex_Er+=LE;
	    break;
	  }
	}
	if (traced) break;
      }
      if(!traced){
	int U=0;
	for (s=0; s<n_seq; s++) {
	  U+=Sali1[s][i-2];
	}
	U = (U==(n_seq)*4?1:0);
	if(/*  pair[Sali1[i+1]][Sali2[j-1]] && */   //only U's are allowed
	    U && j < max_s1 && j > min_s1 && 
	    j > n2 - max_s2 - max_half_stem && 
	    j < n2 -min_s2 -half_stem ) {
	  int min_k, max_k;
	  max_k = MIN2(n2-min_s2,j+max_half_stem+1);
	  min_k = MAX2(j+half_stem+1, n2-max_s2);
	  folden * temp;
	  temp=foldlist[j+1];
	  while(temp->next) {
	    int psc2, psc3;
	    int k = temp->k;
	    for (s=0; s<n_seq; s++) {
	      type2[s]= pair[Sali1[s][i-3]][Sali2[s][k+1]];
	      type3[s]= pair[Sali1[s][i-4]][Sali2[s][k+1]];
	    }
	    psc2 = covscore(type2, n_seq);
	    psc3 = covscore(type3, n_seq);
	    if(psc2>MINPSCORE /*&& pair[Sali1[i-4]][Sali2[k+2]]*/    ){  //introduce structure from RNAfold
	      if(E==c[i-3][k+1]+temp->energy){
		*Loop_E=temp->energy;
		st1[i-3]= '|';
		*u=i-2;
		int a,b;
		//int fix_ij=indx[k-1+1]+j+1;
		for(a=0; a< MISMATCH ;a++){
		  for(b=0; b< MISMATCH ; b++){
		    int ij=indx[k-1-a+1]+j+1+b;
		    if(cLoop[ij]==temp->energy) {
		      //int bla;
		      struc_loop=alisnobacktrack_fold_from_pair(snoseq, j+1+b, k-a-1+1,psct);
		    a=INF; b=INF;	
		    }
		  }
		}
		traced=1;
		traced_r=1;
		i=i-3;j=k+1;
		break;
	      }
	    }
	     if (psc3>MINPSCORE  /*&& pair[Sali1[i-5]][Sali2[k+2]]*/){  //introduce structure from RNAfold
	      if(E==c[i-4][k+1]+temp->energy){
		*Loop_E=temp->energy;
		st1[i-3]= '|';
		*u=i-2;
		int a,b;
		//int fix_ij=indx[k-1+1]+j+1;
		for(a=0; a< MISMATCH ;a++){
		  for(b=0; b< MISMATCH ; b++){
		    int ij=indx[k-1-a+1]+j+1+b;
		    if(cLoop[ij]==temp->energy) {
		      //int bla;
		      struc_loop=alisnobacktrack_fold_from_pair(snoseq, j+1+b, k-a-1+1,psct);
		      a=INF; b=INF;	
		    }
		  }
		}
		traced=1;
		traced_r=1;
		i=i-4;j=k+1;
		break;
	      }
	    } //else if
	    temp=temp->next;
	  } //while temp-> next
	} //test on j 
      }//traced?
    }//traced_r?
    else{
      E = c[i][j]; traced=0;
      st1[i-1] = '<';
      st2[j-1] = '>'; 
      for (s=0; s<n_seq; s++) {
	type[s] = pair[Sali1[s][i]][Sali2[s][j]];
      }
      psc = covscore(type,n_seq);
      for (s=0; s<n_seq; s++) if (type[s]==0) type[s] = 7;
      E += psc;
      *pscd+=psc;
      if (!type) nrerror("backtrack failed in fold duplex c");
      for (k=i-1; (i-k)<MAXLOOP_L; k--) {
	for (l=j+1; l<=n2; l++) {
	  int LE;
	  if (i-k+l-j>2*MAXLOOP_L-2) break;
	  if (abs(i-k-l+j) >= ASS) continue;
	  for (s=LE=0; s<n_seq; s++) {
	    type4 = pair[Sali1[s][k]][Sali2[s][l]];
	    if (type4==0) type4=7;
	    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);
	  }
	  if (E == c[k][l]+LE) {
	    traced=1; 
	    i=k; j=l;
	    *Duplex_El+=LE;
	    break;
	  }
	}
	if (traced) break;
      }
    }
    if (!traced) { 
      for (s=0; s<n_seq; s++) {
	int correction;
	correction = E_ExtLoop(type[s], (i>1) ? Sali1[s][i-1] : -1, (j<n2) ? Sali2[s][j+1] : -1, P);
	*Duplex_El += correction;
	E          -= correction;
	/**
	*** if (i>1)    {E -= P->dangle5[type[s]][Sali1[s][i-1]]; *Duplex_El +=P->dangle5[type[s]][Sali1[s][i-1]];}
	*** if (j<n2)   {E -= P->dangle3[type[s]][Sali2[s][j+1]]; *Duplex_El +=P->dangle3[type[s]][Sali2[s][j+1]];}
	*** if (type[s]>2) {E -= P->TerminalAU;		      *Duplex_El +=P->TerminalAU;}
	**/
      }
      if (E != n_seq * P->DuplexInit) {
	nrerror("backtrack failed in fold duplex end");
      } else break;
    }
  }
/*   if (i>1)  i--;  */
/*   if (j<n2) j++;  */
  ////struc = (char *) space(i0-i+1+j-j0+1+2); //declare final duplex structure
  struc = (char *) space(i0-i+1+n2-1+1+2); //declare final duplex structure
  char * struc2;
  struc2 = (char *) space(n2+1);
  ////char * struct_const;
  for (k=MAX2(i,1); k<=i0; k++) if (!st1[k-1]) st1[k-1] = '.';
  //for (k=j0; k<=j; k++) if (!st2[k-1]) st2[k-1] = struc_loop[k-1];//'.'; normal
  // char * struct_const;
  // struct_const = (char *) space(sizeof(char)*(n2+1));  
  for (k=1; k<=n2; k++) {
    if (!st2[k-1]) st2[k-1] = struc_loop[k-1];//'.';
    struc2[k-1] = st2[k-1];//'.';
    //     if (k>=j0 && k<=j){
    //     	struct_const[k-1]='x';
    //     }
    //     else{
    //     	if(k<j0) {struct_const[k-1]='<';}
    //     	if(k>j) {struct_const[k-1]='>';}
    //     }
  }

  //  char duplexseq_1[j0+1];
  //char duplexseq_2[n2-j+3];
  if(j<n2){
    char **duplexseq_1, **duplexseq_2;
    duplexseq_1 = (char**) space((n_seq+1) * sizeof(char*));
    duplexseq_2 = (char**) space((n_seq+1) * sizeof(char*));
    for(s=0; s<n_seq; s++){
      duplexseq_1[s] = (char*) space((j0)*sizeof(char)); //modfied j0+1
      duplexseq_2[s] = (char*) space((n2-j+2)*sizeof(char)); //modified j+3
      strncpy(duplexseq_1[s], snoseq[s], j0-1); //modified j0
      strcpy(duplexseq_2[s], snoseq[s] + j); //modified j-1
      duplexseq_1[s][j0-1]='\0'; //modified j0
      duplexseq_2[s][n2-j+1]='\0';//modified j+2
    }
    duplexseq_1[n_seq]=NULL;
    duplexseq_2[n_seq]=NULL;
    duplexT temp;
    temp=aliduplexfold((const char**)duplexseq_1, (const char**)duplexseq_2);
    *Loop_D =  MIN2(0,-410 + (int) 100 * temp.energy*n_seq);
    if(*Loop_D){
      int l1,ibegin, iend, jbegin, jend;
      l1=strchr(temp.structure, '&')-temp.structure;
      ibegin=temp.i-l1;
      iend  =temp.i-1;
      jbegin=temp.j;
      jend  =temp.j+strlen(temp.structure)-l1-2-1;
      for(k=ibegin+1; k<=iend+1; k++){
	struc2[k-1]=temp.structure[k-ibegin-1];
      }
      for(k=jbegin+j; k<=jend+j; k++){
	struc2[k-1]=temp.structure[l1+k-j-jbegin+1];
      }
    }
    for(s=0; s<n_seq; s++){
      free(duplexseq_1[s]);
      free(duplexseq_2[s]);
    }
    free(duplexseq_1);free(duplexseq_2);
    free(temp.structure);
  }
  strcpy(struc, st1+MAX2(i-1,0)); strcat(struc, "&"); 
  //strcat(struc, st2);
  strncat(struc, struc2+5, strlen(struc2)-10);
  free(struc2);
  free(struc_loop);
  free(st1); free(st2);
  free(type);free(type2);free(type3);
  //free_arrays();
  return struc;
}







void Lsnoop_subopt(const char *s1, const char *s2, int delta, int w, 
		   const int penalty, const int threshloop, 
		   const int threshLE, const int threshRE, const int threshDE, const int threshTE,const int threshSE,const int threshD,
		   const int distance,
		   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 int alignment_length, const char* name)
{

 int min_colonne=INF;
 int max_pos;
 int max;max=INF;
 //int temp;
 //int nsubopt=10;
 n1 = (int) strlen(s1);
 n2 = (int) strlen(s2);
 int *position;
 position = (int*) space((n1+3)*sizeof(int));


  //int Eminj, Emin_l;
  int i, j; //l1, Emin=INF, i_min=0, j_min=0;
  //char *struc;
  //snoopT mfe;
  int *indx;
  int *mLoop;
  int *cLoop;
  folden **foldlist, **foldlist_XS;
  int Duplex_El, Duplex_Er;
  int Loop_D;
  //int u;
  int Loop_E;
  Duplex_El=0;Duplex_Er=0;Loop_E=0, Loop_D=0;
  snoexport_fold_arrays(&indx, &mLoop, &cLoop, &foldlist, &foldlist_XS); 
  if ((!P) || (fabs(P->temperature - temperature)>1e-6)) {
    snoupdate_fold_params();  if(P) free(P); P = scale_parameters();
    make_pair_matrix();
  }
  
  lc = (int **) space(sizeof(int *) * (5));
  lr = (int **) space(sizeof(int *) * (5));
  for (i=0; i<5; i++) {
  	lc[i] = (int *) space(sizeof(int) * (n2+1));
	lr[i] = (int *) space(sizeof(int) * (n2+1));
	for(j=n2; j>-1; j--){
		lc[i][j]=INF;
		lr[i][j]=INF;
	}
  }
  encode_seqs(s1, s2);
  for (i=1; i<=n1; i++) {
      int idx=i%5;
      int idx_1=(i-1)%5;
      int idx_2=(i-2)%5;
      int idx_3=(i-3)%5;
      int idx_4=(i-4)%5;
    for (j=n2-min_d2; j>min_d1; j--) {
      int type, type2, k;
      type = pair[S1[i]][S2[j]];
      lc[idx][j] = (type) ? P->DuplexInit + 2*penalty : INF;
      lr[idx][j] = INF;
      if(!type) continue;
      if( /*pair[S1[i+1]][S2[j-1]]   &&  check that we have a solid base stack after the mLoop */
	  j < max_s1 && j > min_s1 &&  
	  j > n2 - max_s2 - max_half_stem && 
	  j < n2 -min_s2 -half_stem && S1[i-2]==4) { /*constraint on s2 and i*/
	int min_k, max_k;
	max_k = MIN2(n2-min_s2,j+max_half_stem+1);
	min_k = MAX2(j+half_stem+1, n2-max_s2);
	for(k=min_k; k <= max_k ; k++){	 
	  if(mLoop[indx[k-1]+j+1] < 0){
		}
	  if(pair[S1[i-3]][S2[k]] /*genau zwei ungepaarte nucleotiden --NU--*/
	     && mLoop[indx[k-1]+j+1] < threshloop){  
	    lr[idx][j]=MIN2(lr[idx][j], lc[idx_3][k]+mLoop[indx[k-1]+j+1]);
	  }
	  else if(pair[S1[i-4]][S2[k]] &&  mLoop[indx[k-1]+j+1] < threshloop){/*--NUN--*/
	    lr[idx][j]=MIN2(lr[idx][j], lc[idx_4][k]+mLoop[indx[k-1]+j+1]);
	  }
      	}
      }
      //dangle 5'SIDE relative to the mRNA 
      lc[idx][j] += E_ExtLoop(type, (i>1) ? SS1[i-1] : -1, (j<n2) ? SS2[j+1] : -1, P);
      /**
      *** if (i>1)    lc[idx][j] += P->dangle5[type][SS1[i-1]];
      *** if (j<n2)   lc[idx][j] += P->dangle3[type][SS2[j+1]];
      *** if (type>2) lc[idx][j] += P->TerminalAU;
      **/
      
      if(j<n2 && i>1){
	type2=pair[S1[i-1]][S2[j+1]];
      	if(type2>0){
	  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]);
	  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]);
      	}
      }
      if(j<n2-1 && i>2){
	type2=pair[S1[i-2]][S2[j+2]];
	if(type2>0 ){
	  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]);
	  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]);
	}
      }
      if(j<n2-2 && i>3){
	type2 = pair[S1[i-3]][S2[j+3]];
	if(type2>0){
	  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]);
	  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]);
      	}
      }
      /**
      *** (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)
      **/
      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);
      }
      position[i]=min_colonne;
      if(max>=min_colonne){
	max=min_colonne;
	max_pos=i;
      }
      min_colonne=INF;
 }
 
  free(S1); free(S2); free(SS1); free(SS2);
  if(max<threshTE){
	find_max_snoop(s1, s2, max, alignment_length, position, delta, 
		       distance, penalty, threshloop, threshLE, threshRE, threshDE, 
		       threshTE, threshSE, threshD, half_stem, max_half_stem, min_s2, max_s2, min_s1, max_s1, min_d1, min_d2,name);
   }
  for (i=1; i<5; i++) {free(lc[i]);free(lr[i]);}
  free(lc[0]);free(lr[0]);
  free(lc);free(lr);
  free(position);
  
}  



void Lsnoop_subopt_list(const char *s1, const char *s2, int delta, int w, 
			const int penalty, const int threshloop, 
			const int threshLE, const int threshRE, const int threshDE, const int threshTE,const int threshSE,const int threshD,
			const int distance,
			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 int alignment_length,const char *name)
{
 
 int min_colonne=INF;
 int max_pos;
 int max;max=INF;
 //int temp;
 //int nsubopt=10;
 n1 = (int) strlen(s1);
 n2 = (int) strlen(s2);
 int *position;
 position = (int*) space((n1+3)*sizeof(int));


 //int Eminj, Emin_l;
  int i, j;// l1, Emin=INF, i_min=0, j_min=0;
  //char *struc;
  //snoopT mfe;
  int *indx;
  int *mLoop;
  int *cLoop;
  folden **foldlist, **foldlist_XS;
  int Duplex_El, Duplex_Er;
  int Loop_D;
  //int u;
  int Loop_E;
  Duplex_El=0;Duplex_Er=0;Loop_E=0, Loop_D=0;
  snoexport_fold_arrays(&indx, &mLoop, &cLoop, &foldlist,  &foldlist_XS); 
  if ((!P) || (fabs(P->temperature - temperature)>1e-6)) {
    snoupdate_fold_params();  if(P) free(P); P = scale_parameters();
    make_pair_matrix();
  }
  
  lpair = (int **) space(sizeof(int *) * (6));
  lc    = (int **) space(sizeof(int *) * (6));
  lr    = (int **) space(sizeof(int *) * (6));
  for (i=0; i<6; i++) {
  	lc[i] = (int *) space(sizeof(int) * (n2+1));
	lr[i] = (int *) space(sizeof(int) * (n2+1));
	lpair[i] = (int *) space(sizeof(int) * (n2+1));
	for(j=n2; j>-1; j--){
		lc[i][j]=INF;
		lr[i][j]=INF;
		lpair[i][j]=0;
	}
  }
  encode_seqs(s1, s2);
  int lim_maxj=n2-min_d2 ;
  int lim_minj=min_d1;
  int lim_maxi=n1;
  for (i=5; i<=lim_maxi; i++) {
    int idx=i%5;
    int idx_1=(i-1)%5;
    int idx_2=(i-2)%5;
    int idx_3=(i-3)%5;
    int idx_4=(i-4)%5;

    for (j=lim_maxj; j>lim_minj; j--) {
      int type, type2;// E, k,l;
      type = pair[S1[i]][S2[j]];
      lpair[idx][j] = type;
      lc[idx][j] = (type) ? P->DuplexInit + 2*penalty : INF;
      lr[idx][j] = INF;
      if(!type) continue;
      if( /*pair[S1[i+1]][S2[j-1]] && Be sure it binds*/
	  j < max_s1 && j > min_s1 &&  
	  j > n2 - max_s2 - max_half_stem && 
	  j < n2 -min_s2 -half_stem && S1[i-2]==4 ) { /*constraint on s2 and i*/
	int min_k, max_k;
	max_k = MIN2(n2-min_s2,j+max_half_stem+1);
	min_k = MAX2(j+half_stem+1, n2-max_s2);
	folden * temp;
	temp=foldlist[j+1];
	while(temp->next){
	  int k = temp->k;
	  //if(k >= min_k-1 && k < max_k){ comment to recover normal behaviour
	  if(lpair[idx_3][k+1] /*&& lpair[idx_4][k+2]*/){
	    lr[idx][j]=MIN2(lr[idx][j], lc[idx_3][k+1]+temp->energy);/*--NU--*/
	  }
	  /*else*/ if(lpair[idx_4][k+1]){/*--NUN--*/
	    lr[idx][j]=MIN2(lr[idx][j], lc[idx_4][k+1]+temp->energy);
	  }
	    // }
	  temp=temp->next;
	}
      }
      //dangle 5'SIDE relative to the mRNA 
      lc[idx][j] += E_ExtLoop(type,  SS1[i-1] , SS2[j+1] , P);
      /**
      ***      lc[idx][j] += P->dangle5[type][SS1[i-1]];
      ***      lc[idx][j] += P->dangle3[type][SS2[j+1]];
      ***      if (type>2) lc[idx][j] += P->TerminalAU;
      **/
      //      if(j<n2 && i>1){
      //type2=pair[S1[i-1]][S2[j+1]];
	type2=lpair[idx_1][j+1];
	if(type2>0 ){
	  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]);
	  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]);
	}
	//}
	//      if(j<n2-1 && i>2){
	//type2=pair[S1[i-2]][S2[j+2]];
	type2=lpair[idx_2][j+2];
	if(type2>0 ){
	  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]);
	  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]);
	  //	}
      }
	//      if(j<n2-2 && i>3){
	//type2 = pair[S1[i-3]][S2[j+3]];
	type2 =lpair[idx_3][j+3];
	if(type2>0 ){
	  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]);
	  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]);
	  //	}
      }
      //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);
      int bla;
      bla=lr[idx][j]+E_ExtLoop(rtype[type], SS2[j-1] , SS1[i+1], P)+2*penalty;
      min_colonne=MIN2(bla, min_colonne);
    }
    position[i]=min_colonne;
    if(max>=min_colonne){
      max=min_colonne;
      max_pos=i;
      }
    min_colonne=INF;
 }
 
  free(S1); free(S2); free(SS1); free(SS2);
  if(max<threshTE){
      find_max_snoop(s1, s2, max, alignment_length, position, 
		     delta, distance, penalty, threshloop, 
		     threshLE, threshRE, threshDE, threshTE, threshSE, threshD,
		     half_stem, max_half_stem, min_s2, max_s2, min_s1, max_s1, min_d1, min_d2,name);
   }
  for (i=1; i<6; i++) {free(lc[i]);free(lr[i]);free(lpair[i]);}
  free(lc[0]);free(lr[0]);free(lpair[0]);
  free(lc);free(lr);free(lpair);
  free(position);
}  


PRIVATE void find_max_snoop(const char *s1, const char *s2,const int max,  const int alignment_length, const int* position, const int delta, 
			    const int distance, const int penalty, const int threshloop,  const int threshLE, const int threshRE, 
			    const int threshDE, const int threshTE, const int threshSE, const int threshD, 
			    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)
{
  int count=0;
  int pos=n1+1;
  int threshold = MIN2(threshTE , max + delta );
  //printf("threshTE %d max %d\n", threshTE, max);
  //#pragma omp parallel for
  //  for(pos=n1+1;pos>distance;pos--){
  while(pos-- > 5){
    int temp_min=0;
    if(position[pos]<(threshold)){
      int search_range;
      search_range=distance+1;
      while(--search_range){
	if(position[pos-search_range]<=position[pos-temp_min]){
	  temp_min=search_range;
	}
      }
      pos-=temp_min;
      int begin=MAX2(6, pos-alignment_length+1);
      char *s3 = (char*) space(sizeof(char)*(pos-begin+3+12));
      strcpy(s3, "NNNNN");
      strncat(s3, (s1+begin-1), pos-begin+2);
      strcat(s3,"NNNNN\0");
      //printf("%s s3\n", s3); 
      snoopT test;
      test = snoopfold(s3, s2, penalty, threshloop, threshLE, threshRE, threshDE, threshD, half_stem, max_half_stem, min_s2, max_s2, min_s1, 
		       max_s1, min_d1, min_d2);
      if(test.energy==INF){
	free(s3);
	continue;
      }
      if(test.Duplex_El > threshLE * 0.01 || test.Duplex_Er > threshRE * 0.01  ||
	 test.Loop_D > threshD * 0.01 || (test.Duplex_Er + test.Duplex_El) > threshDE * 0.01 ||
	 (test.Duplex_Er + test.Duplex_El + test.Loop_E + test.Loop_D + 410) > threshSE*0.01) {
	free(test.structure);free(s3);
	continue;
      }
      int l1;
      l1 = strchr(test.structure, '&')-test.structure;
      
      int shift=0;
      if(test.i > strlen(s3)-10){
	test.i--;
	l1--; 
      }
      if(test.i-l1<0){
	l1--;
	shift++;
      }
      char *target_struct =  (char*) space(sizeof(char) * (strlen(test.structure)+1));
      strncpy(target_struct, test.structure+shift, l1);
      strncat(target_struct, test.structure + (strchr(test.structure, '&')-
					       test.structure), strlen(test.structure) - (strchr(test.structure, '&')-
											  test.structure));
      strcat(target_struct,"\0");
      char *target; 
      target = (char *) space(l1+1);
      strncpy(target, (s3+test.i+5-l1), l1);
      target[l1]='\0';
      char *s4;
      s4 = (char*) space(sizeof(char) *(strlen(s2)-9));
      strncpy(s4, s2+5, strlen(s2)-10);
      s4[strlen(s2)-10]='\0';
      printf("%s %3d,%-3d;%3d : %3d,%-3d (%5.2f = %5.2f + %5.2f + %5.2f + %5.2f + 4.1 )\n%s&%s\n", 
	     target_struct,begin + test.i-5-l1,begin + test.i -6 , begin + test.u -6, 
	     test.j+1, test.j + (strrchr(test.structure,'>') - strchr(test.structure,'>'))+1 ,
	     test.Loop_D + test.Duplex_El + test.Duplex_Er + test.Loop_E + 4.10, test.Duplex_El,
	     test.Duplex_Er, test.Loop_E, test.Loop_D,target,s4);
      if(name){
	char *temp_seq;
	char *temp_struc;
	char psoutput[100];
	temp_seq = (char*) space(sizeof(char)*(l1+n2-9));
	temp_struc = (char*) space(sizeof(char)*(l1+n2-9));
	strcpy(temp_seq, target);
	strcat(temp_seq, s4);
	strncpy(temp_struc, target_struct, l1);
	strcat(temp_struc, target_struct+l1+1);
	temp_seq[n2+l1-10]='\0';
	temp_struc[n2+l1-10]='\0';
	cut_point = l1+1;
	char str[16];char upos[16];
	strcpy(psoutput,"sno_");
	sprintf(str,"%d",count);
	strcat(psoutput,str);
	sprintf(upos,"%d",begin + test.u - 6);
	strcat(psoutput,"_u_");
	strcat(psoutput,upos);
	strcat(psoutput,"_");
	strcat(psoutput,name);
	strcat(psoutput,".ps\0");
	PS_rna_plot_snoop_a(temp_seq, temp_struc, psoutput, NULL, NULL);
	cut_point = -1;
	free(temp_seq);
	free(temp_struc);
	count++;
	//free(psoutput);
      }
      free(s4);
      free(test.structure);
      free(target_struct);
      free(target);
      free(s3);
    }
  }
  
}








snoopT snoopfold(const char *s1, const char *s2, 
		 const int penalty, const int threshloop, const int threshLE, const int threshRE, const int threshDE,
		 const int threshD,
		 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) {
  //int Eminj, Emin_l;
  int i, j, l1, Emin=INF, i_min=0, j_min=0;
  char *struc;
  snoopT mfe;
  int *indx;
  int *mLoop;
  int *cLoop;
  folden** foldlist, **foldlist_XS;
  int Duplex_El, Duplex_Er;
  int Loop_D;
  int u;
  int Loop_E;
  Duplex_El=0;Duplex_Er=0;Loop_E=0, Loop_D=0;
  snoexport_fold_arrays(&indx, &mLoop, &cLoop,&foldlist, &foldlist_XS ); 
  n1 = (int) strlen(s1);
  n2 = (int) strlen(s2);
  
  if ((!P) || (fabs(P->temperature - temperature)>1e-6)) {
    snoupdate_fold_params();  if(P) free(P); P = scale_parameters();
    make_pair_matrix();
  }
  
  c = (int **) space(sizeof(int *) * (n1+1));
  r = (int **) space(sizeof(int *) * (n1+1));
  for (i=0; i<=n1; i++) {
  	c[i] = (int *) space(sizeof(int) * (n2+1));
	r[i] = (int *) space(sizeof(int) * (n2+1));
	for(j=n2; j>-1; j--){
		c[i][j]=INF;
		r[i][j]=INF;
	}
  }
  encode_seqs(s1, s2);
  for (i=6; i<=n1-5; i++) {
    for (j=n2-min_d2; j>min_d1; j--) {
      int type, type2, E, k,l;
      type = pair[S1[i]][S2[j]];
      c[i][j] = (type ) ? P->DuplexInit : INF;
      if(!type) continue;
      if(/*  pair[S1[i+1]][S2[j-1]] &&  */
	 j < max_s1 && j > min_s1 &&  
	 j > n2 - max_s2 - max_half_stem && 
	 j < n2 -min_s2 -half_stem && S1[i-2]==4  ) { /*constraint on s2 and i*/
	int min_k, max_k;
	max_k = MIN2(n2-min_s2,j+max_half_stem);
	min_k = MAX2(j+half_stem, n2-max_s2);
	folden * temp;
	temp=foldlist[j+1];
	while(temp->next){
	  int k = temp->k;
	  //if(k >= min_k-1 && k < max_k){ uncomment to recovernormal behaviour
	  if(pair[S1[i-3]][S2[k+1]] /*&& pair[S1[i-4]][S2[k+2]]*/ ){
	    r[i][j]=MIN2(r[i][j], c[i-3][k+1]+temp->energy);
	  }
	  /*else*/ if(pair[S1[i-4]][S2[k+1]] /*&& pair[S1[i-5]][S2[k+2]]*/ ){
	    r[i][j]=MIN2(r[i][j], c[i-4][k+1]+temp->energy);
	  }
	  //}
	  temp=temp->next;
	}
      }
      //dangle 5'SIDE relative to the mRNA 
      /**
      *** c[i][j] += P->dangle5[type][SS1[i-1]];
      *** c[i][j] += P->dangle3[type][SS2[j+1]];
      *** if (type>2) c[i][j] += P->TerminalAU;
      **/
      c[i][j]+=E_ExtLoop(type, SS1[i-1] , SS2[j+1], P);
      for (k=i-1; k>0 && (i-k)<MAXLOOP_L; k--) {
	for (l=j+1; l<=n2 ; l++) {
	  if (i-k+l-j>2*MAXLOOP_L-2) break;
	  if (abs(i-k-l+j) >= ASS ) continue;
	  type2 = pair[S1[k]][S2[l]];
	  if (!type2) continue;
	  E = E_IntLoop(i-k-1, l-j-1, type2, rtype[type],
			SS1[k+1], SS2[l-1], SS1[i-1], SS2[j+1],P);
	  c[i][j] = MIN2(c[i][j], c[k][l]+E+(i-k+l-j)*penalty);
	  r[i][j] = MIN2(r[i][j], r[k][l]+E+(i-k+l-j)*penalty);
	}
      }
      E = r[i][j]; 
      /**
      *** if (i<n1) E += P->dangle3[rtype[type]][SS1[i+1]];              
      *** if (j>1)  E += P->dangle5[rtype[type]][SS2[j-1]]; 
      *** f (type>2) E += P->TerminalAU;
      **/
      E+=E_ExtLoop(rtype[type], (j > 1) ? SS2[j-1] : -1, (i<n1) ? SS1[i+1] : -1, P);
      if (E<Emin) {
	Emin=E; i_min=i; j_min=j;
      } 
    }
  }
  if(Emin > 0){
  	printf("no target found under the constraints chosen\n");
	for (i=0; i<=n1; i++) {free(r[i]);free(c[i]);}
	free(c);
	free(r);
	free(S1); free(S2); free(SS1); free(SS2);
	mfe.energy=INF;
	return mfe;
  }
  struc = snoop_backtrack(i_min, j_min,s2, &Duplex_El, &Duplex_Er, &Loop_E, &Loop_D, 
			  &u, penalty, threshloop, threshLE, threshRE,threshDE, threshD,
  		          half_stem, max_half_stem, min_s2, max_s2, min_s1, max_s1, min_d1, min_d2);
/*   if (i_min<n1-5) i_min++; */
/*   if (j_min>1 ) j_min--; */
  l1 = strchr(struc, '&')-struc;
  mfe.i = i_min-5;
  mfe.j = j_min-5;
  mfe.u = u -5;
  mfe.Duplex_Er = (float) Duplex_Er/100;
  mfe.Duplex_El = (float) Duplex_El/100;
  mfe.Loop_D = (float) Loop_D/100;
  mfe.Loop_E = (float) Loop_E/100;
  mfe.energy = (float) Emin/100 ;
  mfe.structure = struc;
  if (!delay_free) {
    for (i=0; i<=n1; i++) {free(r[i]);free(c[i]);}
    free(c);
    free(r);
    free(S1); free(S2); free(SS1); free(SS2);
  }
  return mfe;
}

PRIVATE int snoopfold_XS_fill(const char *s1, const char *s2, const int **access_s1,
		 const int penalty, const int threshloop, const int threshLE, const int threshRE, const int threshDE,
		 const int threshD,
		 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) {
  //int Eminj, Emin_l;
  int i, j, Emin=INF, i_min=0, j_min=0;
  //char *struc;
  //snoopT mfe;
  int *indx;
  int *mLoop;
  int *cLoop;
  folden** foldlist, **foldlist_XS;
  int Duplex_El, Duplex_Er;
  int Loop_D;
  //int u;
  int Loop_E;
  Duplex_El=0;Duplex_Er=0;Loop_E=0, Loop_D=0;
  snoexport_fold_arrays(&indx, &mLoop, &cLoop,&foldlist, &foldlist_XS ); 
  n1 = (int) strlen(s1);
  n2 = (int) strlen(s2);
  
  if ((!P) || (fabs(P->temperature - temperature)>1e-6)) {
    snoupdate_fold_params();  if(P) free(P); P = scale_parameters();
    make_pair_matrix();
  }
  
  c_fill = (int **) space(sizeof(int *) * (n1+1));
  r_fill = (int **) space(sizeof(int *) * (n1+1));
  for (i=0; i<=n1; i++) {
  	c_fill[i] = (int *) space(sizeof(int) * (n2+1));
	r_fill[i] = (int *) space(sizeof(int) * (n2+1));
	for(j=n2; j>-1; j--){
		c_fill[i][j]=INF;
		r_fill[i][j]=INF;
	}
  }
  encode_seqs(s1, s2);

  int di[5];
  di[0]=0;  
  for (i=6; i<=n1-5; i++) {
    di[1]=access_s1[5][i]   - access_s1[4][i-1];	   
    di[2]=access_s1[5][i-1] - access_s1[4][i-2] + di[1];
    di[3]=access_s1[5][i-2] - access_s1[4][i-3] + di[2];
    di[4]=access_s1[5][i-3] - access_s1[4][i-4] + di[3];
    di[1]=MIN2(di[1],165);
    di[2]=MIN2(di[2],330);
    di[3]=MIN2(di[3],495);
    di[4]=MIN2(di[4],660);
    for (j=n2-min_d2; j>min_d1; j--) {
      int type, type2, E, k,l;
      type = pair[S1[i]][S2[j]];
      c_fill[i][j] = (type ) ? P->DuplexInit : INF;
      if(!type) continue;
      if(/*  pair[S1[i+1]][S2[j-1]] &&  */
	 j < max_s1 && j > min_s1 &&  
	 j > n2 - max_s2 - max_half_stem && 
	 j < n2 -min_s2 -half_stem && S1[i-2]==4  ) { /*constraint on s2 and i*/
	int min_k, max_k;
	max_k = MIN2(n2-min_s2,j+max_half_stem);
	min_k = MAX2(j+half_stem, n2-max_s2);
	folden * temp;
	temp=foldlist[j+1];
	while(temp->next){
	  int k = temp->k;
	  //if(k >= min_k-1 && k < max_k){ uncomment to recovernormal behaviour
	  if(pair[S1[i-3]][S2[k+1]] /*&& pair[S1[i-4]][S2[k+2]]*/ ){
	    r_fill[i][j]=MIN2(r_fill[i][j], c_fill[i-3][k+1]+temp->energy+ di[3]);
	  }
	  /*else*/ if(pair[S1[i-4]][S2[k+1]] /*&& pair[S1[i-5]][S2[k+2]]*/ ){
	    r_fill[i][j]=MIN2(r_fill[i][j], c_fill[i-4][k+1]+temp->energy + di[4]);
	  }
	  //}
	  temp=temp->next;
	}
      }
      //dangle 5'SIDE relative to the mRNA 
      /**
      *** c_fill[i][j] += P->dangle5[type][SS1[i-1]];
      *** c_fill[i][j] += P->dangle3[type][SS2[j+1]];
      *** if (type>2) c_fill[i][j] += P->TerminalAU;
      **/
      c_fill[i][j]+= E_ExtLoop(type, SS1[i-1], SS2[j+1],P);
      for (k=i-1; k>0 && (i-k)<MAXLOOP_L; k--) {
	for (l=j+1; l<=n2 ; l++) {
	  if (i-k+l-j>2*MAXLOOP_L-2) break;
	  if (abs(i-k-l+j) >= ASS ) continue;
	  type2 = pair[S1[k]][S2[l]];
	  if (!type2) continue;
	  E = E_IntLoop(i-k-1, l-j-1, type2, rtype[type],
			SS1[k+1], SS2[l-1], SS1[i-1], SS2[j+1],P);
	  c_fill[i][j] = MIN2(c_fill[i][j], c_fill[k][l]+E+di[i-k]);
	  r_fill[i][j] = MIN2(r_fill[i][j], r_fill[k][l]+E+di[i-k]);
	}
      }
      E = r_fill[i][j]; 
      /**
      ***  if (i<n1) E += P->dangle3[rtype[type]][SS1[i+1]];              
      ***  if (j>1)  E += P->dangle5[rtype[type]][SS2[j-1]]; 
      *** if (type>2) E += P->TerminalAU;
      **/
      E+= E_ExtLoop(rtype[type], (j > 1) ? SS2[j-1] : -1, (i<n1) ? SS1[i+1] : -1, P);
      if (E<Emin) {
	Emin=E; i_min=i; j_min=j;
      } 
    }
  }
  return Emin;
}



PUBLIC snoopT *snoop_subopt(const char *s1, const char *s2, int delta, int w, 
			    const int penalty, const int threshloop, 
			    const int threshLE, const int threshRE, const int threshDE, const int threshTE, const int threshSE, const int threshD,
			    const int distance, 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) {



  //printf("%d %d\n", min_s2, max_s2);
  int i,j, n1, n2, E, n_subopt=0, n_max;
  char *struc;
  snoopT mfe;
  snoopT *subopt;
  int thresh;

  int Duplex_El, Duplex_Er, Loop_E;
  int Loop_D;
  Duplex_El=0; Duplex_Er=0; Loop_E=0;Loop_D=0;
  int u;
  u=0;
  n_max=16;
  subopt = (snoopT *) space(n_max*sizeof(snoopT));
  delay_free=1;
  mfe = snoopfold(s1, s2, penalty, threshloop, threshLE, threshRE, threshDE,threshD,
                  half_stem, max_half_stem,
		  min_s2, max_s2, min_s1, max_s1, min_d1, min_d2);



  if(mfe.energy > 0){
  	free(subopt);
	delay_free=0;
	return NULL;
  }
  thresh = MIN2((int) ((mfe.Duplex_Er + mfe.Duplex_El + mfe.Loop_E)*100+0.1 + 410) + delta, threshTE );
 //subopt[n_subopt++]=mfe;
  free(mfe.structure);
  
  n1 = strlen(s1); n2=strlen(s2);
  for (i=n1; i>0; i--) {
    for (j=1; j<=n2; j++) {
      int type, Ed;
      type = pair[S2[j]][S1[i]];
      if (!type) continue;
      E = Ed = r[i][j];
      /**
      *** if (i<n1) Ed += P->dangle3[type][SS1[i+1]]; 
      *** if (j>1)  Ed += P->dangle5[type][SS2[j-1]]; 
      *** if (type>2) Ed += P->TerminalAU;
      **/
      Ed+= E_ExtLoop(type, (j > 1) ? SS2[j-1] : -1, (i<n1) ? SS1[i+1] : -1, P);
      if (Ed>thresh) continue;
      /* too keep output small, remove hits that are dominated by a
	 better one close (w) by. For simplicity we do test without
	 adding dangles, which is slightly inaccurate. 
      */ 
      /* w=1; */
/*       for (ii=MAX2(i-w,1); (ii<=MIN2(i+w,n1)) && type; ii++) {  */
/*         for (jj=MAX2(j-w,1); jj<=MIN2(j+w,n2); jj++) */
/*           if (r[ii][jj]<E) {type=0; break;} */
/*       } */
      if (!type) continue;

      struc = snoop_backtrack(i,j,s2, &Duplex_El, &Duplex_Er, &Loop_E, &Loop_D, &u, penalty, threshloop,threshLE,threshRE,threshDE,threshD, 
                        half_stem, max_half_stem, min_s2, max_s2, min_s1, max_s1, min_d1, min_d2);
      if (Duplex_Er > threshRE || Duplex_El > threshLE || Loop_D > threshD ||
         (Duplex_Er + Duplex_El) > threshDE || 
	 (Duplex_Er + Duplex_El + Loop_E) > threshTE ||
	 (Duplex_Er + Duplex_El + Loop_E + Loop_D + 410) > threshSE) {
	 	//printf(" Duplex_Er %d threshRE %d Duplex_El %d threshLE %d \n"
		//       " Duplex_Er + Duplex_El %d  threshDE %d \n"
		//       " Duplex_Er + Duplex_El + Loop_E %d  threshTE %d \n"
		//       " Duplex_Er + Duplex_El + Loop_E + Loop_D %d  threshSE %d \n", 
		//	 Duplex_Er , threshRE , Duplex_El ,threshLE,
	        //         Duplex_Er + Duplex_El, threshDE,
        	//         Duplex_Er + Duplex_El+  Loop_E , threshTE,
                //	 Duplex_Er + Duplex_El+  Loop_E + Loop_D, threshSE); 
	 	Duplex_Er=0; 
		Duplex_El=0;
		Loop_E = 0;
		Loop_D = 0;
		u=0,
		free(struc);
		continue;
	}

      if (n_subopt+1>=n_max) {
	n_max *= 2;
	subopt = (snoopT *) xrealloc(subopt, n_max*sizeof(snoopT));
      }
      subopt[n_subopt].i = i-5;
      subopt[n_subopt].j = j-5;
      subopt[n_subopt].u = u-5;
      subopt[n_subopt].Duplex_Er = Duplex_Er * 0.01;
      subopt[n_subopt].Duplex_El = Duplex_El * 0.01;
      subopt[n_subopt].Loop_E = Loop_E * 0.01;
      subopt[n_subopt].Loop_D = Loop_D * 0.01;
      subopt[n_subopt].energy = (Duplex_Er +Duplex_El + Loop_E + Loop_D + 410) * 0.01 ;
      subopt[n_subopt++].structure = struc;
      Duplex_Er=0; Duplex_El=0; Loop_E=0; Loop_D=0;u=0;
    }
  }
  
  for (i=0; i<=n1; i++) {free(c[i]);free(r[i]);}
  free(c);free(r);
  free(S1); free(S2); free(SS1); free(SS2);
  delay_free=0;

  if (snoop_subopt_sorted) qsort(subopt, n_subopt, sizeof(snoopT), compare);
  subopt[n_subopt].i =0;
  subopt[n_subopt].j =0;
  subopt[n_subopt].structure = NULL;
  return subopt;
}

PUBLIC void snoop_subopt_XS(const char *s1, const char *s2, const int **access_s1, int delta, int w, 
			    const int penalty, const int threshloop, 
			    const int threshLE, const int threshRE, const int threshDE, const int threshTE, const int threshSE, const int threshD,
			    const int distance, 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 int alignment_length, const char *name) {



  //printf("%d %d\n", min_s2, max_s2);
  int i,j, E,  n_max;
  //char *struc;
  //snoopT mfe;

  int thresh;

  int Duplex_El, Duplex_Er, Loop_E;
  int Loop_D;
  Duplex_El=0; Duplex_Er=0; Loop_E=0;Loop_D=0;
  int u;
  u=0;
  n_max=16;
  delay_free=1;
  int Emin = snoopfold_XS_fill(s1, s2, access_s1,penalty, threshloop, threshLE, threshRE, threshDE,threshD,
			       half_stem, max_half_stem,
			       min_s2, max_s2, min_s1, max_s1, min_d1, min_d2);
  if(Emin > 0){
    delay_free=0;
  }
  thresh = MIN2(-100, threshTE +alignment_length*30);  
  //  n1=strlen(s1); 
  //  n2=strlen(s2);
  
  int n3=strlen(s1);
  int n4=strlen(s2);
  S1_fill = (short*)space(sizeof(short)*(n3+2));
  S2_fill = (short*)space(sizeof(short)*(n4+2));
  SS1_fill = (short*)space(sizeof(short)*(n3+1));
  SS2_fill = (short*)space(sizeof(short)*(n4+1));
  memcpy(S1_fill, S1, sizeof(short)*n3+2);
  memcpy(S2_fill, S2, sizeof(short)*n4+2);
  memcpy(SS1_fill, SS1, sizeof(short)*n3+1);
  memcpy(SS2_fill, SS2, sizeof(short)*n4+1);
  free(S1);free(S2);free(SS1);free(SS2);
  int count=0;
  for (i=n3-5; i>0; i--) {
    for (j=1; j<=n4; j++) {
      int type, Ed;
      type = pair[S2_fill[j]][S1_fill[i]];
      if (!type) continue;
      E = Ed = r_fill[i][j];
      /**
      ***if (i<n3) Ed += P->dangle3[type][SS1_fill[i+1]]; 
      ***if (j>1)  Ed += P->dangle5[type][SS2_fill[j-1]]; 
      ***if (type>2) Ed += P->TerminalAU;
      **/
      Ed+=E_ExtLoop(type, (j > 1) ? SS2[j-1] : -1, (i<n3) ? SS1[i+1] : -1, P);
      if (Ed>thresh) continue;
      
      /* to keep output small, remove hits that are dominated by a
	 better one close (w) by. For simplicity we do test without
	 adding dangles, which is slightly inaccurate. 
      */ 
/*       w=10;  */
/*       for (ii=MAX2(i-w,1); (ii<=MIN2(i+w,n3-5)) && type; ii++) {   */
/* 	for (jj=MAX2(j-w,1); jj<=MIN2(j+w,n4-5); jj++)  */
/* 	  if (r_fill[ii][jj]<E) {type=0; break;}  */
/*       }  */
/*       i=ii;j=jj; */
      if (!type) continue;
      int begin=MAX2(5, i-alignment_length);
      int end  =MIN2(n3-5, i-1); 
      char *s3 = (char*) space(sizeof(char)*(end-begin+2)+5);
      strncpy(s3, (s1+begin), end - begin +1);
      strcat(s3,"NNNNN\0");
      int n5 = strlen(s3);
      snoopT test = snoopfold_XS(s3, s2, access_s1, i, j ,penalty, 
			  threshloop, threshLE, threshRE, 
			  threshDE, threshD, half_stem, 
			  max_half_stem, min_s2, max_s2, min_s1, 
			  max_s1, min_d1, min_d2);
      if(test.energy==INF){
	free(s3);
	continue;
      }
      if( test.Duplex_El > threshLE * 0.01 ||test.Duplex_Er > threshRE * 0.01  || 
	  test.Loop_D > threshD * 0.01 || (test.Duplex_Er + test.Duplex_El) > threshDE * 0.01 || 
	  (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) 
	{ 
	  free(test.structure);free(s3); 
	  continue; 
	}
      char *s4; 
      s4 = (char*) space(sizeof(char) *(n4-9)); 
      strncpy(s4, s2+5, n4-10); 
      s4[n4-10]='\0';
      
      char *s5 = space(sizeof(char) * n5-test.i+2-5);
      strncpy(s5,s3+test.i-1,n5-test.i+1-5);
      s5[n5-test.i+1-5]='\0';
      float dE = ((float) (access_s1[n5-test.i+1-5][i]))*0.01;
      printf("%s %3d,%-3d;%3d : %3d,%-3d (%5.2f = %5.2f + %5.2f + %5.2f + %5.2f + %5.2f + 4.10)\n%s&%s\n" ,  
	     test.structure, i  -  (n5 - test.i) ,i - 5, i - (n5 - test.u ),
             j-5, j-5 + (strrchr(test.structure,'>') - strchr(test.structure,'>')), 
	     test.Loop_D + test.Duplex_El + test.Duplex_Er + test.Loop_E + 4.10+dE, test.Duplex_El, 
	     test.Duplex_Er, test.Loop_E, test.Loop_D,dE , s5,s4);
      if(name){
	int begin_t, end_t, begin_q, end_q, and, pipe,k; 
	char psoutput[100];
	begin_q=0;
	end_q=n4-10;
	begin_t=0;
	end_t=n5-test.i+ 1-5;
	and=end_t+1;
	pipe=test.u -test.i + 1;
	cut_point = end_t +1 ;
	char *catseq, *catstruct;// *fname; 
	catseq = (char*) space(n5 + end_q -begin_q +2);
	catstruct = (char*) space(n5 + end_q -begin_q +2);
	strcpy(catseq, s5);
	strncpy(catstruct, test.structure, end_t);
	strcat(catseq, s4);
	strncat(catstruct, test.structure+end_t+1, end_q-begin_q+1);
	catstruct[end_t - begin_t + end_q-begin_q+2]='\0';
	catseq[end_t - begin_t + end_q-begin_q+2]='\0';
	int *relative_access;
	relative_access = space(sizeof(int)*strlen(s5));
	relative_access[0] = access_s1[1][i - (n5  - test.i) + 5];
	for(k=1;k<strlen(s5);k++){
	  relative_access[k] =  access_s1[k+1][i - (n5  - test.i) + k + 5] -  access_s1[k][i - (n5  - test.i) + k + 4];
	}
	char str[16];char upos[16];
	strcpy(psoutput,"sno_XS_");
	sprintf(str,"%d",count);
	strcat(psoutput,str);
	sprintf(upos,"%d",i - (n5 - test.u ));
	strcat(psoutput,"_u_");
	strcat(psoutput,upos);
	strcat(psoutput,"_");
	strcat(psoutput,name);
	strcat(psoutput,".ps\0");
	PS_rna_plot_snoop_a(catseq, catstruct, psoutput,relative_access,NULL);
	free(catseq);free(catstruct);free(relative_access);
	count++;
      }
      free(s3);free(s4);free(s5);free(test.structure);
    }
  }  
  for (i=0; i<=n3; i++) {free(c_fill[i]);free(r_fill[i]);}
  free(c_fill);free(r_fill);
  free(S1_fill); free(S2_fill); free(SS1_fill); free(SS2_fill);
  delay_free=0;
}




PRIVATE char *snoop_backtrack(int i, int j, const char* snoseq, 
			      int *Duplex_El, int *Duplex_Er, 
			      int *Loop_E, int *Loop_D, int *u, 
			      const int penalty, const int threshloop, 
			      const int threshLE, const int threshRE, const int threshDE, const int threshD,
			      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) {
  /* backtrack structure going backwards from i, and forwards from j 
     return structure in bracket notation with & as separator */
  int k, l, type, type2, E, traced, i0, j0;
  int traced_r=0; //flag for following backtrack in c or r
  char *st1, *st2, *struc;
  char *struc_loop;

  st1 = (char *) space(sizeof(char)*(n1+1));
  st2 = (char *) space(sizeof(char)*(n2+1));
  int *indx;
  int *mLoop;
  int *cLoop;
  folden **foldlist, **foldlist_XS;
  type=pair[S1[i]][S2[j]];
  snoexport_fold_arrays(&indx, &mLoop, &cLoop,&foldlist, &foldlist_XS ); 
  i0=i; j0=j;
  /**
  *** if (i<n1)   *Duplex_Er += P->dangle3[rtype[type]][SS1[i+1]];
  *** if (j>1)    *Duplex_Er += P->dangle5[rtype[type]][SS2[j-1]];
  *** if (type>2) *Duplex_Er += P->TerminalAU;
  **/
  *Duplex_Er += E_ExtLoop(rtype[type], (j > 1) ? SS2[j-1] : -1, (i<n1) ? SS1[i+1] : -1, P);
  while (i>0 && j<=n2-min_d2 ) {
    if(!traced_r) {
      E = r[i][j]; traced=0;
      st1[i-1] = '<';
      st2[j-1] = '>'; 
      type = pair[S1[i]][S2[j]];
      if (!type) nrerror("backtrack failed in fold duplex r");
      for (k=i-1; k>0 && (i-k)<MAXLOOP_L; k--) {
	for (l=j+1; l<=n2 ; l++) {
	  int LE;
	  if (i-k+l-j>2*MAXLOOP_L-2) break;
	  if (abs(i-k-l+j) >= ASS) continue;
	  
	  type2 = pair[S1[k]][S2[l]];
	  if (!type2) continue;
	  LE = E_IntLoop(i-k-1, l-j-1, type2, rtype[type],
			 SS1[k+1], SS2[l-1], SS1[i-1], SS2[j+1],P);
	  if (E == r[k][l]+LE+(i-k+l-j)*penalty) {
	    traced=1; 
	    i=k; j=l;
	    *Duplex_Er+=LE;
	    break;
	  }
	}
	if (traced) break;
      }
      if(!traced){
	if(/*  pair[S1[i+1]][S2[j-1]] && */  
	    j < max_s1 && j > min_s1 && 
	    j > n2 - max_s2 - max_half_stem && 
	    j < n2 -min_s2 -half_stem && 
	    S1[i-2]==4) {
	  int min_k, max_k;
	  max_k = MIN2(n2-min_s2,j+max_half_stem+1);
	  min_k = MAX2(j+half_stem+1, n2-max_s2);
	  folden * temp;
	  temp=foldlist[j+1];
	  while(temp->next) {
	    int k = temp->k;
	    if(pair[S1[i-3]][S2[k+1]] /*&& pair[S1[i-4]][S2[k+2]]*/    ){  //introduce structure from RNAfold
	      if(E==c[i-3][k+1]+temp->energy){
		*Loop_E=temp->energy;
		st1[i-3]= '|';
		*u=i-2;
		int a,b;
		//int fix_ij=indx[k-1+1]+j+1;
		for(a=0; a< MISMATCH ;a++){
		  for(b=0; b< MISMATCH ; b++){
		    int ij=indx[k-1-a+1]+j+1+b;
		    if(cLoop[ij]==temp->energy) {
		      struc_loop=snobacktrack_fold_from_pair(snoseq, j+1+b, k-a-1+1);
		    a=INF; b=INF;	
		    }
		  }
		}
		traced=1;
		traced_r=1;
		i=i-3;j=k+1;
		break;
	      }
	    }
	    /*else*/ if (pair[S1[i-4]][S2[k+1]] /*&& pair[S1[i-5]][S2[k+2]]*/){  //introduce structure from RNAfold
	      if(E==c[i-4][k+1]+temp->energy){
		*Loop_E=temp->energy;
		st1[i-3]= '|';
		*u=i-2;
		int a,b;
		//int fix_ij=indx[k-1+1]+j+1;
		for(a=0; a< MISMATCH ;a++){
		  for(b=0; b< MISMATCH ; b++){
		    int ij=indx[k-1-a+1]+j+1+b;
		    if(cLoop[ij]==temp->energy) {
		      struc_loop=snobacktrack_fold_from_pair(snoseq, j+1+b, k-a-1+1);
		      a=INF; b=INF;	
		    }
		  }
		}
		traced=1;
		traced_r=1;
		i=i-4;j=k+1;
		break;
	      }
	    } //else if
	    temp=temp->next;
	  } //while temp-> next
	} //test on j 
      }//traced?
    }//traced_r?
    else{
      E = c[i][j]; traced=0;
      st1[i-1] = '<';
      st2[j-1] = '>'; 
      type = pair[S1[i]][S2[j]];
      if (!type) nrerror("backtrack failed in fold duplex c");
      for (k=i-1; (i-k)<MAXLOOP_L; k--) {
	for (l=j+1; l<=n2; l++) {
	  int LE;
	  if (i-k+l-j>2*MAXLOOP_L-2) break;
	  if (abs(i-k-l+j) >= ASS) continue;
	  type2 = pair[S1[k]][S2[l]];
	  if (!type2) continue;
	  LE = E_IntLoop(i-k-1, l-j-1, type2, rtype[type],
			 SS1[k+1], SS2[l-1], SS1[i-1], SS2[j+1],P);
	  if (E == c[k][l]+LE+(i-k+l-j)*penalty) {
	    traced=1; 
	    i=k; j=l;
	    *Duplex_El+=LE;
	    break;
	  }
	}
	if (traced) break;
      }
    }
    if (!traced) { 
      int correction;
      correction = E_ExtLoop(type, (i>1) ? SS1[i-1] : -1, (j<n2) ? SS2[j+1] : -1, P);
      E-=correction;
      *Duplex_El+=correction;
      /**
      *** if (i>1)    {E -= P->dangle5[type][SS1[i-1]]; *Duplex_El +=P->dangle5[type][SS1[i-1]];}
      *** if (j<n2)   {E -= P->dangle3[type][SS2[j+1]]; *Duplex_El +=P->dangle3[type][SS2[j+1]];}
      *** if (type>2) {E -= P->TerminalAU;		    *Duplex_El +=P->TerminalAU;}
      **/
      if (E != P->DuplexInit) {
	nrerror("backtrack failed in fold duplex end");
      } else break;
    }
  }
/*   if (i>1)  i--; */
/*   if (j<n2) j++; */  
  ////struc = (char *) space(i0-i+1+j-j0+1+2); //declare final duplex structure
  struc = (char *) space(i0-i+1+n2-1+1+2); //declare final duplex structure
  char * struc2;
  struc2 = (char *) space(n2+1);
  ////char * struct_const;
  for (k=MAX2(i,1); k<=i0; k++) if (!st1[k-1]) st1[k-1] = '.';
  //for (k=j0; k<=j; k++) if (!st2[k-1]) st2[k-1] = struc_loop[k-1];//'.'; normal
  // char * struct_const;
  // struct_const = (char *) space(sizeof(char)*(n2+1));  
  for (k=1; k<=n2; k++) {
    if (!st2[k-1]) st2[k-1] = struc_loop[k-1];//'.';
    struc2[k-1] = st2[k-1];//'.';
    //     if (k>=j0 && k<=j){
    //     	struct_const[k-1]='x';
    //     }
    //     else{
    //     	if(k<j0) {struct_const[k-1]='<';}
    //     	if(k>j) {struct_const[k-1]='>';}
    //     }
  }
  char duplexseq_1[j0];
  char duplexseq_2[n2-j+2];
  if(j<n2){
    strncpy(duplexseq_1, snoseq, j0-1);
    strcpy(duplexseq_2, snoseq + j);
    duplexseq_1[j0-1]='\0';
    duplexseq_2[n2-j+1]='\0';
    duplexT temp;
    temp=duplexfold(duplexseq_1, duplexseq_2);
    *Loop_D =  MIN2(0,-410 + (int) 100 * temp.energy);
    if(*Loop_D){
      int l1,ibegin, iend, jbegin, jend;
      l1=strchr(temp.structure, '&')-temp.structure;
      ibegin=temp.i-l1;
      iend  =temp.i-1;
      jbegin=temp.j;
      jend  =temp.j+strlen(temp.structure)-l1-2-1;
      for(k=ibegin+1; k<=iend+1; k++){
	struc2[k-1]=temp.structure[k-ibegin-1];
      }
      for(k=jbegin+j; k<=jend+j; k++){
	struc2[k-1]=temp.structure[l1+k-j-jbegin+1];
      } 
    }
    free(temp.structure);
  } 
  strcpy(struc, st1+MAX2(i-1,0)); strcat(struc, "&"); 
  //strcat(struc, st2);
  strncat(struc, struc2+5, strlen(struc2)-10);
  free(struc2);
  free(struc_loop);
  free(st1); free(st2);
  //free_arrays();
  return struc;
}

void Lsnoop_subopt_list_XS(const char *s1, const char *s2,  const int **access_s1, int delta, int w, 
			const int penalty, const int threshloop, 
			const int threshLE, const int threshRE, const int threshDE, const int threshTE,const int threshSE,const int threshD,
			const int distance,
			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 int alignment_length, const char *name)
{
 
 int min_colonne=INF;
 int max_pos;
 int max;max=INF;
 //int temp;
 //int nsubopt=10;
 n1 = (int) strlen(s1);
 n2 = (int) strlen(s2);
 int *position;
 int *position_j;
 int min_j_colonne;
 int max_pos_j=INF; 
 position = (int*) space((n1+3)*sizeof(int));
 position_j = (int*) space((n1+3)*sizeof(int));

 //int Eminj, Emin_l;
  int i, j;// l1, Emin=INF, i_min=0, j_min=0;
  //char *struc;
  //snoopT mfe;
  int *indx;
  int *mLoop;
  int *cLoop;
  folden **foldlist, **foldlist_XS;
  int Duplex_El, Duplex_Er;
  int Loop_D;
  //int u;
  int Loop_E;
  Duplex_El=0;Duplex_Er=0;Loop_E=0, Loop_D=0;
  snoexport_fold_arrays(&indx, &mLoop, &cLoop, &foldlist, &foldlist_XS); 
  if ((!P) || (fabs(P->temperature - temperature)>1e-6)) {
    snoupdate_fold_params();  if(P) free(P); P = scale_parameters();
    make_pair_matrix();
  }
  
  lpair = (int **) space(sizeof(int *) * (6));
  lc    = (int **) space(sizeof(int *) * (6));
  lr    = (int **) space(sizeof(int *) * (6));
  for (i=0; i<6; i++) {
  	lc[i] = (int *) space(sizeof(int) * (n2+1));
	lr[i] = (int *) space(sizeof(int) * (n2+1));
	lpair[i] = (int *) space(sizeof(int) * (n2+1));
	for(j=n2; j>-1; j--){
		lc[i][j]=INF;
		lr[i][j]=INF;
		lpair[i][j]=0;
	}
  }
  encode_seqs(s1, s2);
  int lim_maxj=n2-min_d2 ;
  int lim_minj=min_d1;
  int lim_maxi=n1-5;
  for (i=5; i<=lim_maxi; i++) {
    int idx=i%5;
    int idx_1=(i-1)%5;
    int idx_2=(i-2)%5;
    int idx_3=(i-3)%5;
    int idx_4=(i-4)%5;
    int di1, di2, di3, di4;
    di1 = access_s1[5][i]   - access_s1[4][i-1];	   
    di2 =access_s1[5][i-1] - access_s1[4][i-2] + di1;
    di3 = access_s1[5][i-2] - access_s1[4][i-3] + di2;
    di4 = access_s1[5][i-3] - access_s1[4][i-4] + di3;
    di1=MIN2(di1,165);
    di2=MIN2(di2,330);
    di3=MIN2(di3,495);
    di4=MIN2(di4,660);
    for (j=lim_maxj; j>lim_minj; j--) {
      int type, type2;// E, k,l;
      type = pair[S1[i]][S2[j]];
      lpair[idx][j] = type;
      lc[idx][j] = (type) ? P->DuplexInit + access_s1[1][i] : INF;
      lr[idx][j] = INF;
      if(!type) continue;
      if( /*pair[S1[i+1]][S2[j-1]] && Be sure it binds*/
	  j < max_s1 && j > min_s1 &&  
	  j > n2 - max_s2 - max_half_stem && 
	  j < n2 -min_s2 -half_stem && S1[i-2]==4 ) { /*constraint on s2 and i*/
	int min_k, max_k;
	max_k = MIN2(n2-min_s2,j+max_half_stem+1);
	min_k = MAX2(j+half_stem+1, n2-max_s2);
	folden * temp;
	temp=foldlist[j+1];
	while(temp->next){
	  int k = temp->k;
	  //if(k >= min_k-1 && k < max_k){ comment to recover normal behaviour
	  if(lpair[idx_3][k+1] && lc[idx_3][k+1] /*+di3*/ < 411 /*&& lpair[idx_4][k+2]*/){ // remove second condition
	    lr[idx][j]=MIN2(lr[idx][j], di3 + lc[idx_3][k+1]+temp->energy);/*--NU--*/
	  }
	  /*else*/ if(lpair[idx_4][k+1] && /*di4 +*/ lc[idx_4][k+1] < 411  ){/*--NUN--*/ // remove second condition 
	    lr[idx][j]=MIN2(lr[idx][j], di4 + lc[idx_4][k+1]+temp->energy);
	  }
	    // }
	  temp=temp->next;
	}
      }
      //dangle 5'SIDE relative to the mRNA 
      /**
      *** lc[idx][j] += P->dangle5[type][SS1[i-1]];
      *** lc[idx][j] += P->dangle3[type][SS2[j+1]];
      *** if (type>2) lc[idx][j] += P->TerminalAU;
      **/
      lc[idx][j]+=E_ExtLoop(type, SS1[i-1] ,  SS2[j+1] , P);
      //      if(j<n2 && i>1){
      //type2=pair[S1[i-1]][S2[j+1]];
	type2=lpair[idx_1][j+1];
	if(type2>0 ){
	  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]);
	  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]);
	}
	type2=lpair[idx_2][j+2];
	if(type2>0 ){
	  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]);
	  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]);
	 
      }
	type2 =lpair[idx_3][j+3];
	if(type2>0 ){
	  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]);
	  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]);

      }
      int bla;
      int temp2;
      temp2=min_colonne;
      bla=lr[idx][j] + E_ExtLoop(rtype[type], SS2[j-1], SS1[i+1] , P);
	/**
	*** (type>2?P->TerminalAU:0)+P->dangle3[rtype[type]][SS1[i+1]]+P->dangle5[rtype[type]][SS2[j-1]];
	**/
      min_colonne=MIN2(bla, min_colonne);
      if(temp2>min_colonne){
	min_j_colonne=j;
      }
    }
    position[i]=min_colonne;
    if(max>=min_colonne){
      max=min_colonne;
      max_pos=i;
      max_pos_j=min_j_colonne;
      }
    position_j[i]=min_j_colonne;
    min_colonne=INF;
 }
  free(S1); free(S2); free(SS1); free(SS2);

  if(max<threshTE + 30*alignment_length){
    find_max_snoop_XS(s1, s2, access_s1,max,alignment_length, position, position_j,
		     delta, distance, penalty, threshloop, 
		     threshLE, threshRE, threshDE, threshTE, threshSE, threshD,
		      half_stem, max_half_stem, min_s2, max_s2, min_s1, max_s1, min_d1, min_d2,name);
   }
  for (i=1; i<6; i++) {free(lc[i]);free(lr[i]);free(lpair[i]);}
  free(lc[0]);free(lr[0]);free(lpair[0]);
  free(lc);free(lr);free(lpair);
  free(position);free(position_j);
}  


PRIVATE void find_max_snoop_XS(const char *s1, const char *s2, const int **access_s1, 
			       const int max,  const int alignment_length, 
			       const int* position, const int* position_j, const int delta, 
			       const int distance, const int penalty, const int threshloop,  const int threshLE, const int threshRE, 
			       const int threshDE, const int threshTE, const int threshSE, const int threshD, 
			       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){
  int count=0;
  int n3=strlen(s1);
  int n4=strlen(s2);
  int pos=n1-4;
  int max_pos_j;
  int threshold = MIN2(threshTE + alignment_length*30, -100);
  //printf("threshTE %d max %d\n", threshTE, max);
  //#pragma omp parallel for
  //  for(pos=n1+1;pos>distance;pos--){
  while(pos-->5){
    int temp_min=0;
    if(position[pos]<(threshold)){
      int search_range;
      search_range=distance+1;
      while(--search_range){
	if(position[pos-search_range]<=position[pos-temp_min]){
	  temp_min=search_range;
	}
      }
      pos-=temp_min;
      max_pos_j=position_j[pos];
      int begin=MAX2(5, pos-alignment_length);
      int end  =MIN2(n3-5, pos-1); 
      char *s3 = (char*) space(sizeof(char)*(end-begin+2)+5);
      strncpy(s3, (s1+begin), end - begin +1);
      strcat(s3,"NNNNN\0");

      int n5 = strlen(s3);
      snoopT test;
      test = snoopfold_XS(s3, s2, access_s1, pos, max_pos_j ,penalty, 
			  threshloop, threshLE, threshRE, 
			  threshDE, threshD, half_stem, 
			  max_half_stem, min_s2, max_s2, min_s1, 
			  max_s1, min_d1, min_d2);
      if(test.energy==INF){
	free(s3);
	continue;
      }
      if( test.Duplex_El > threshLE * 0.01 ||test.Duplex_Er > threshRE * 0.01  || 
	 test.Loop_D > threshD * 0.01 || (test.Duplex_Er + test.Duplex_El) > threshDE * 0.01 || 
	 (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) { 
        free(test.structure);free(s3); 
        continue; 
      }
      
      char *s4; 
      s4 = (char*) space(sizeof(char) *(n4-9)); 
      strncpy(s4, s2+5, n4-10); 
      s4[n4-10]='\0';

      char *s5 = space(sizeof(char) * n5-test.i+2-5);
      strncpy(s5,s3+test.i-1,n5-test.i+1-5);
      s5[n5-test.i+1-5]='\0';
      float dE = ((float) (access_s1[n5-test.i+1-5][pos]))*0.01;
      printf("%s %3d,%-3d;%3d : %3d,%-3d (%5.2f = %5.2f + %5.2f + %5.2f + %5.2f + %5.2f + 4.10)\n%s&%s\n" ,  
	     test.structure, pos  -  (n5 - test.i) ,pos - 5, pos - (n5 - test.u ),
             max_pos_j-5, max_pos_j-5 + (strrchr(test.structure,'>') - strchr(test.structure,'>')), 
	     test.Loop_D + test.Duplex_El + test.Duplex_Er + test.Loop_E + 4.10+dE, test.Duplex_El, 
	     test.Duplex_Er, test.Loop_E, test.Loop_D,dE , s5,s4);
      if(name){
	int begin_t, end_t, begin_q, end_q, and, pipe, i; 
	char psoutput[100];
	begin_q=0;
	end_q=n4-10;
	begin_t=0;
	end_t=n5-test.i+ 1-5;
	and=end_t+1;
	pipe=test.u -test.i + 1;
	cut_point = end_t +1 ;
	char *catseq, *catstruct;// *fname; 
	catseq = (char*) space(n5 + end_q -begin_q +2);
	catstruct = (char*) space(n5 + end_q -begin_q +2);
	strcpy(catseq, s5);
	strncpy(catstruct, test.structure, end_t);
	strcat(catseq, s4);
	strncat(catstruct, test.structure+end_t+1, end_q-begin_q+1);
	catstruct[end_t - begin_t + end_q-begin_q+2]='\0';
	catseq[end_t - begin_t + end_q-begin_q+2]='\0';
	int *relative_access;
	relative_access = space(sizeof(int)*strlen(s5));

	relative_access[0] = access_s1[1][pos - (n5  - test.i) + 5];
	for(i=1;i<strlen(s5);i++){
	  relative_access[i] =  access_s1[i+1][pos - (n5  - test.i) + i + 5] -  access_s1[i][pos - (n5  - test.i) + i + 4];
	}
	char str[16];
	char upos[16];
	strcpy(psoutput,"sno_XS_");
	sprintf(str,"%d",count);
	strcat(psoutput,str);
	sprintf(upos,"%d",pos - (n5 - test.u ));
	strcat(psoutput,"_u_");
	strcat(psoutput,upos);
	strcat(psoutput,"_");
	strcat(psoutput,name);
	strcat(psoutput,".ps\0");
	PS_rna_plot_snoop_a(catseq, catstruct, psoutput,relative_access,NULL);
	free(catseq);free(catstruct);free(relative_access);
	count++;
      }
      free(s3);free(s4);free(s5);free(test.structure);
    }
  }
}

snoopT snoopfold_XS(const char *s1, const char *s2, const int **access_s1, const int pos_i, const int pos_j,
		 const int penalty, const int threshloop, const int threshLE, const int threshRE, const int threshDE,
		 const int threshD,
		 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) {
  //  int Eminj, Emin_l;
  int a,b,i, j, Emin=INF, a_min=0, b_min=0;
  char *struc;
  snoopT mfe;
  int *indx;
  int *mLoop;
  int *cLoop;
  folden** foldlist, **foldlist_XS;
  int Duplex_El, Duplex_Er;
  int Loop_D;
  int u;
  int Loop_E;

  Duplex_El=0;Duplex_Er=0;Loop_E=0, Loop_D=0;
  snoexport_fold_arrays(&indx, &mLoop, &cLoop,&foldlist, &foldlist_XS ); 
  n1 = (int) strlen(s1);
  n2 = (int) strlen(s2);
  
  if ((!P) || (fabs(P->temperature - temperature)>1e-6)) {
    snoupdate_fold_params(); if(P) free(P); P = scale_parameters();
    make_pair_matrix();
  }
  
  c = (int **) space(sizeof(int *) * (n1+1));
  r = (int **) space(sizeof(int *) * (n1+1));
  for (i=0; i<=n1; i++) {
  	c[i] = (int *) space(sizeof(int) * (n2+1));
	r[i] = (int *) space(sizeof(int) * (n2+1));
	for(j=n2; j>-1; j--){
		c[i][j]=INF;
		r[i][j]=INF;
	}
  }
  encode_seqs(s1, s2);
  i=n1-5;
  j=pos_j;
  //printf("tar: %s\nsno: %s\n ", s1, s2);
  //printf("pos_i %d pos_j %d\n", pos_i, pos_j);
  //printf("type %d n1 %d n2 %d S1[n1] %d S2[n2] %d", pair[S1[i]][S2[j]], n1, n2, S1[n1], S2[n2]);
  int type, type2, E, p,q;    
  r[i][j] = P->DuplexInit; 
  //r[i][j] += P->dangle3[rtype[type]][SS1[i+1]] + P->dangle5[rtype[type]][SS2[j-1]]; 
  
  if(pair[S1[i]][S2[j]]>2) r[i][j]+=P->TerminalAU;
  for(a=i-1; a>0;a--){ //i-1
    r[a+1][0]=INF;
    for(b=j+1; b<=n2-min_d2;b++){ //j+1
      r[a][b]=INF;
      type = pair[S1[a]][S2[b]]; 
       if(!type) continue; 
       if(S1[a+1]==4){ 
	 folden * temp; 
	 temp=foldlist_XS[b-1]; 
	 while(temp->next) {     
	   int k = temp->k; 
	   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
	     c[a][b]=MIN2(c[a][b], r[a+3][k-1]+temp->energy); 
	   }
	   temp=temp->next;
	 }
       }
       if(S1[a+2]==4){
	 folden * temp; 
	 temp=foldlist_XS[b-1]; 
	 while(temp->next){ 
	   int k = temp->k; 
	   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 
	     c[a][b]=MIN2(c[a][b], r[a+4][k-1]+temp->energy); 
	   }
	   temp=temp->next;
	 }
       }
       for(p=a+1; p<n1 && (p-a) < MAXLOOP_L; p++) { //p < n1
	 for (q=b-1; q > 1  ; q--) {  //q > 1
	   if (p-a+b-q>2*MAXLOOP_L-2) break; 
	   if (abs((p-a)-(b-q)) >= ASS ) continue; 
	   type2 = pair[S1[p]][S2[q]]; 
	   if (!type2) continue; 
	   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); 
	   c[a][b] = MIN2(c[a][b], c[p][q]+E); 
	   r[a][b] = MIN2(r[a][b], r[p][q]+E); 
	 }
       }
       E = c[a][b];  
       if (type>2) E += P->TerminalAU;  
       //       E +=P->dangle5[rtype[type]][SS1[i+1]];
       //E +=P->dangle5[rtype[type]][SS2[j-1]]; 
       E+=access_s1[i-a+1][pos_i]; 
       if (E<Emin) { 
	 Emin=E; a_min=a; b_min=b; 
       }  
    }
  }
  if(Emin > 0){ 
    printf("no target found under the constraints chosen\n");
    for (i=0; i<=n1; i++) {free(r[i]);free(c[i]);}
    free(c);
    free(r); 
    free(S1); free(S2); free(SS1); free(SS2);
    mfe.energy=INF;
    return mfe;
  }  
  type2=pair[S1[a_min]][S2[b_min]];
  if(type2>2) Emin +=P->TerminalAU;
  mfe.energy = ((float) (Emin))/100;
  struc = snoop_backtrack_XS(a_min, b_min,s2, &Duplex_El, &Duplex_Er, &Loop_E, &Loop_D, 
			     &u, penalty, threshloop, threshLE, threshRE,threshDE, threshD, 
			     half_stem, max_half_stem, min_s2, max_s2, min_s1, max_s1, min_d1, min_d2); 

  mfe.i = a_min;
  mfe.j = b_min;
  mfe.u = u;
  mfe.Duplex_Er = (float) Duplex_Er/100;
  mfe.Duplex_El = (float) Duplex_El/100;
  mfe.Loop_D = (float) Loop_D/100;
  mfe.Loop_E = (float) Loop_E/100;
  mfe.energy = (float) Emin/100 ;
  mfe.structure = struc;
  return mfe;
}

PRIVATE char *snoop_backtrack_XS(int i, int j, const char* snoseq, 
			      int *Duplex_El, int *Duplex_Er, 
			      int *Loop_E, int *Loop_D, int *u, 
			      const int penalty, const int threshloop, 
			      const int threshLE, const int threshRE, const int threshDE, const int threshD,
			      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) {
  /* backtrack structure going backwards from i, and forwards from j 
     return structure in bracket notation with & as separator */
  int k, l, type, type2, E, traced, i0, j0;
  int traced_c=0; //flag for following backtrack in c or r
  char *st1, *st2, *struc;
  char *struc_loop;

  st1 = (char *) space(sizeof(char)*(n1+1));
  st2 = (char *) space(sizeof(char)*(n2+1));
  int *indx;
  int *mLoop;
  int *cLoop;
  folden **foldlist, **foldlist_XS;
  type=pair[S1[i]][S2[j]];
  snoexport_fold_arrays(&indx, &mLoop, &cLoop,&foldlist, &foldlist_XS ); 
  i0=i;j0=j;
  //  i0=MAX2(i,1); j0=MIN2(j+1,n2);
  while (i<=n1 && j>=1 ) {
    if(!traced_c) {
      E = c[i][j]; traced=0;
      st1[i] = '<';
      st2[j-1] = '>'; 
      type = pair[S1[i]][S2[j]];
      if (!type) nrerror("backtrack failed in fold duplex c");
      for (k=i+1; k>0 && (k-i)<MAXLOOP_L; k++) {
	for (l=j-1; l>=1 ; l--) {
	  int LE;
	  if (k-i+j-l>2*MAXLOOP_L-2) break;
	  if (abs(k-i-j+l) >= ASS) continue;
	  type2 = pair[S1[k]][S2[l]];
	  if (!type2) continue;
	  LE = E_IntLoop(k-i-1, j-l-1, type2, rtype[type],
			 SS1[i+1], SS2[j-1], SS1[k-1], SS2[l+1],P);
	  if (E == c[k][l]+LE) {
	    traced=1; 
	    i=k; j=l;
	    *Duplex_El+=LE;
	    break;
	  }
	}
	if (traced) break;
      }
      if(!traced){
	if(S1[i+1]==4) {
	  folden * temp;
	  temp=foldlist_XS[j-1];
	  while(temp->next) {
	    int k = temp->k;
	    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 ) {
	      if(E==r[i+3][k-1]+temp->energy){
		*Loop_E=temp->energy;
		st1[i+1]= '|';
		st1[i+2]='.';
		*u=i+1;
		int a,b;
		for(a=0; a< MISMATCH ;a++){
		  for(b=0; b< MISMATCH ; b++){
		    int ij=indx[j-1-a]+k+b;
		    if(cLoop[ij]==temp->energy) {
		      struc_loop=snobacktrack_fold_from_pair(snoseq, k+b, j-1-a); 
		      a=INF; b=INF;	
		    }
		  }
		}
		traced=1;
		traced_c=1;
		i=i+3;j=k-1;
		break;
	      }
	    }
	    temp=temp->next;
	  }
	}
	if (S1[i+2]==4){  //introduce structure from RNAfold
	  folden * temp;
	  temp=foldlist_XS[j-1];
	  while(temp->next) {
	    int k = temp->k;
	    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 ) {
	      if(E==r[i+4][k-1]+temp->energy){
		*Loop_E=temp->energy;
		st1[i+2]= '|';
		st1[i+1]=st1[i+3]='.';
		*u=i+2;
		int a,b;
		for(a=0; a< MISMATCH ;a++){
		  for(b=0; b< MISMATCH ; b++){
		    int ij=indx[j-1-a]+k+b;
		    if(cLoop[ij]==temp->energy) {
		      struc_loop=snobacktrack_fold_from_pair(snoseq, k+b, j-a-1);
		      a=INF; b=INF;	
		    }
		  }
		}
		traced=1;
		traced_c=1;
		i=i+4;j=k-1;
		break;
	      }
	    }
	    temp=temp->next;
	  }
	}
      }//traced?
    }//traced_r?
    else{
      E = r[i][j]; traced=0;
      st1[i] = '<';
      st2[j-1] = '>'; 
      type = pair[S1[i]][S2[j]];
      if (!type) nrerror("backtrack failed in fold duplex r");
      for (k=i+1; k>0 && (k-i)<MAXLOOP_L; k++) {
	for (l=j-1; l>=1 ; l--) {
	  int LE;
	  if (k-i+j-l>2*MAXLOOP_L-2) break;
	  if (abs(k-i-j+l) >= ASS) continue;
	  type2 = pair[S1[k]][S2[l]];
	  if (!type2) continue;
	  LE = E_IntLoop(k-i-1, j-l-1, type2, rtype[type],
			 SS1[i+1], SS2[j-1], SS1[k-1], SS2[l+1],P);
	  if (E == r[k][l]+LE) {
	    traced=1; 
	    i=k; j=l;
	    *Duplex_Er+=LE;
	    break;
	  }
	}
	if (traced) break;
      }
    }
    if (!traced) { 
/*       if (i>1)    {E -= P->dangle5[type][SS1[i-1]]; *Duplex_El +=P->dangle5[type][SS1[i-1]];} */
/*       if (j<n2)   {E -= P->dangle3[type][SS2[j+1]]; *Duplex_El +=P->dangle3[type][SS2[j+1]];} */
      if (type>2) {E -= P->TerminalAU;	*Duplex_Er +=P->TerminalAU;}
      if (E != P->DuplexInit) {
	nrerror("backtrack failed in fold duplex end");
      } else break;
    }
  }

  
  ////struc = (char *) space(i0-i+1+j-j0+1+2); //declare final duplex structure
  struc = (char *) space(i-i0+1+n2); //declare final duplex structure
  char * struc2;
  struc2 = (char *) space(n2+1);
  ////char * struct_const;

  for (k=MIN2(i0,1); k<=i; k++) if (!st1[k-1]) st1[k-1] = '.';
  //for (k=j0; k<=j; k++) if (!st2[k-1]) st2[k-1] = struc_loop[k-1];//'.'; normal
  // char * struct_const;
  // struct_const = (char *) space(sizeof(char)*(n2+1));  
  for (k=1; k<=n2; k++) {
    if (!st2[k-1]) st2[k-1] = struc_loop[k-1];//'.';
    struc2[k-1] = st2[k-1];//'.';
    //     if (k>=j0 && k<=j){
    //     	struct_const[k-1]='x';
    //     }
    //     else{
    //     	if(k<j0) {struct_const[k-1]='<';}
    //     	if(k>j) {struct_const[k-1]='>';}
    //     }
  }
  char duplexseq_1[j];
  char duplexseq_2[n2-j0+2];
  if(j0<n2){
    strncpy(duplexseq_1, snoseq, j-1);
    strcpy(duplexseq_2, snoseq + j0);
    duplexseq_1[j-1]='\0';
    duplexseq_2[n2-j0+1]='\0';
    duplexT temp;
    temp=duplexfold(duplexseq_1, duplexseq_2);
    *Loop_D =  MIN2(0,-410 + (int) 100 * temp.energy);
    if(*Loop_D){
      int l1,ibegin, iend, jbegin, jend;
      l1=strchr(temp.structure, '&')-temp.structure;
      ibegin=temp.i-l1;
      iend  =temp.i-1;
      jbegin=temp.j;
      jend  =temp.j+strlen(temp.structure)-l1-2-1;
      for(k=ibegin+1; k<=iend+1; k++){
	struc2[k-1]=temp.structure[k-ibegin-1];
      }
      for(k=jbegin+j0; k<=jend+j0; k++){
	struc2[k-1]=temp.structure[l1+k-j0-jbegin+1];
      } 
    }
    free(temp.structure);
  } 
  strcpy(struc, st1+MAX2(i0,1)); strcat(struc, "&"); 
  //strcat(struc, st2);
  strncat(struc, struc2+5, strlen(struc2)-10);
  free(struc2);
  free(struc_loop);
  free(st1); free(st2);
  
    for (i=0; i<=n1; i++) {free(r[i]);free(c[i]);}
    free(c);
    free(r);
    free(S1);free(S2);free(SS1);free(SS2);
    //free_arrays();
  return struc;
}

PRIVATE int covscore(const int *types, int n_seq) {
  /* calculate co-variance bonus for a pair depending on  */
  /* compensatory/consistent mutations and incompatible seqs */
  /* should be 0 for conserved pairs, >0 for good pairs      */
#define NONE -10000 /* score for forbidden pairs */
  int k,l,s,score, pscore;
  int dm[7][7]={{0,0,0,0,0,0,0}, /* hamming distance between pairs */
                {0,0,2,2,1,2,2} /* CG */,
                {0,2,0,1,2,2,2} /* GC */,
                {0,2,1,0,2,1,2} /* GU */,
                {0,1,2,2,0,2,1} /* UG */,
                {0,2,2,1,2,0,2} /* AU */,
                {0,2,2,2,1,2,0} /* UA */};
  
  int pfreq[8]={0,0,0,0,0,0,0,0};
  for (s=0; s<n_seq; s++)
    pfreq[types[s]]++;

  if (pfreq[0]*2>n_seq) return NONE;
  for (k=1,score=0; k<=6; k++) /* ignore pairtype 7 (gap-gap) */
    for (l=k+1; l<=6; l++)
      /* scores for replacements between pairtypes    */
      /* consistent or compensatory mutations score 1 or 2  */
      score += pfreq[k]*pfreq[l]*dm[k][l];
  
  /* counter examples score -1, gap-gap scores -0.25   */
  pscore = cv_fact * 
    ((UNIT*score)/n_seq - nc_fact*UNIT*(pfreq[0] + pfreq[7]*0.25));
  return pscore;
}

/*---------------------------------------------------------------------------*/

PRIVATE short * aliencode_seq(const char *sequence) {
  unsigned int i,l;
  short *Stemp;
  l = strlen(sequence);
  Stemp = (short *) space(sizeof(short)*(l+2));
  Stemp[0] = (short) l;

  /* make numerical encoding of sequence */
  for (i=1; i<=l; i++)
    Stemp[i]= (short) encode_char(toupper(sequence[i-1]));

  /* for circular folding add first base at position n+1 */
  //Stemp[l+1] = Stemp[1];

  return Stemp;
}

PRIVATE short * encode_seq(const char *sequence) {
  unsigned int i,l;
  short *S;
  l = strlen(sequence);
extern double nc_fact;
  S = (short *) space(sizeof(short)*(l+2));
  S[0] = (short) l;

  /* make numerical encoding of sequence */
  for (i=1; i<=l; i++)
    S[i]= (short) encode_char(toupper(sequence[i-1]));
  /* for circular folding add first base at position n+1 */
  S[l+1] = S[1];

  return S;
}

PRIVATE void encode_seqs(const char *s1, const char *s2) {
  unsigned int i,l;

  l = strlen(s1);
  S1 = encode_seq(s1);
  SS1= (short *) space(sizeof(short)*(l+1));
  /* SS1 exists only for the special X K and I bases and energy_set!=0 */
  
  for (i=1; i<=l; i++) { /* make numerical encoding of sequence */
    SS1[i] = alias[S1[i]];   /* for mismatches of nostandard bases */
  }

  l = strlen(s2);
  S2 = encode_seq(s2);
  SS2= (short *) space(sizeof(short)*(l+1));
  /* SS2 exists only for the special X K and I bases and energy_set!=0 */
  
  for (i=1; i<=l; i++) { /* make numerical encoding of sequence */
    SS2[i] = alias[S2[i]];   /* for mismatches of nostandard bases */
  }
}

PRIVATE int compare(const void *sub1, const void *sub2) {
  int d;
  if (((snoopT *) sub1)->energy > ((snoopT *) sub2)->energy)
    return 1;
  if (((snoopT *) sub1)->energy < ((snoopT *) sub2)->energy)
    return -1;
  d = ((snoopT *) sub1)->i - ((snoopT *) sub2)->i;
  if (d!=0) return d;
  return  ((snoopT *) sub1)->j - ((snoopT *) sub2)->j;
}