/* Last changed Time-stamp: <2007-12-05 13:55:42 ronny> */
/*
                  Ineractive Access to folding Routines

                  c Ivo L Hofacker
                  Vienna RNA package
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <ctype.h>
#include <unistd.h>
#include <string.h>
#include <time.h>
#include "snofold.h"
#include "fold.h"
#include "energy_par.h"
#include "snoop.h"
#include "part_func.h"
#include "fold_vars.h"
#include "PS_dot.h"
#include "utils.h"
#include "aln_util.h"
#include "RNAsnoop_cmdl.h"
static void  aliprint_struc(snoopT *dup, const char **s1, const char **s2, char**,char**,int count, int);
static void  print_struc(snoopT *dup, const char *s1, const char *s2, char*, char*, int,int);
//static void  print_struc_L(snoopT const *dup, const char *s1, const char *s2, char*, char*);


static void redraw_output(char *fname, char *output, int plfold_up_flag, char *suffix, char *access, char *sname, char *tname);
static int ** read_plfold_i(char *fname, const int beg, const int end); //read plfold profiles
static int ** read_rnaup(char *fname, const int beg, const int end);
static int get_max_u(const char *s, char delim);
extern int cut_point;
/*@unused@*/
static char UNUSED rcsid[] = "$Id: RNAfold.c,v 1.23 2007/12/05 13:04:08 ivo Exp $";

#define PRIVATE static
#define MAX_NUM_NAMES    500

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

int main(int argc, char *argv[])
{
  struct        RNAsnoop_args_info args_info;  
  char *string_s, *line_s, *name_s, *temp_s; /*string for snoRNA*/
  char *string_t, *line_t, *temp_t; /*string for target RNA*/
  char  *sname=NULL,*tname=NULL/*name of the sno RNa file, mRNA file respectively*/;
  char *access;char *result_file;
  char *output_directory;
  output_directory=NULL;
  result_file=NULL;
  access=NULL;
  char *suffix;
  suffix=NULL;
  FILE *sno, *mrna;
  int fast;fast=0;
  //long int elapTicks;
  //clock_t Begin, End;
  int plfold_up_flag=0;
  int   nice,i,r,l,length_s,//, length_t, 
   penalty,         /*extension penalty*/
   threshloop,         /*energy threshold on loop*/
   threshLE,         /*energy threshold on the S2 part*/
   threshRE,        /*energy threshold on the S1 part*/
   threshDE,         /*Total duplex energy threshold*/
   threshTE,        /*Total duplex + Loop energy*/
   threshSE,        /*Sum of all energies*/
   threshD,            /*lower stem energy*/
   half_stem,         /*minimal length of stem*/
   max_half_stem,        /*maximal length of stem*/
   max_s2                /*maximale position wo stem anfangen darf in 3->5*/,
   min_s2                 /*minimale position wo stem anfangen darf in 3->5*/,
   max_s1                /*minimal position wo stem enden darf in 3->5 */,
   min_s1                /*maximale position wo stem enden darf in 3->5*/,
   distance        /*distance between two subopts*/,
   min_d1                /*minimal distance between 5' sno and first duplex*/,
   min_d2                /*minimal distance between 3' sno and second duplex*/,
   max_asymm        /*maximal asymmetry in the stem interior loop*/,
   alignment_length/*maximal target RNA alignment length*/,
   alignment       /*flag to use ali version*/,
    delta,
   redraw          /*if used (option I) allow to redraw command line output into ps files */;
  
  int  noconv=0;
  string_s=NULL;
  string_t=NULL;
  plfold_up_flag=0;
  alignment=0;
  redraw=0;
  nice=0;
  fast=0;
  snoop_subopt_sorted=0;

  if(RNAsnoop_cmdline_parser(argc,argv, &args_info)!=0) exit(1);
  
  if(args_info.produce_ps_given                      ) redraw=1;
  output_directory=strdup(args_info.output_directory_arg);
  if(args_info.direct_redraw_given                      ) nice=1;
  if(args_info.from_RNAup_given                      ) {plfold_up_flag=2;access=strdup(args_info.from_RNAup_arg);}
  suffix=strdup(args_info.suffix_arg);
  if(args_info.from_RNAplfold_given                      ) {plfold_up_flag=1;access=strdup(args_info.from_RNAplfold_arg);}
  if(args_info.alignment_mode_given                      ) alignment=1;
  if(args_info.query_given                              ) sname=strdup(args_info.query_arg);
  if(args_info.target_given                              ) tname=strdup(args_info.target_arg);
  delta=(int) (100*(args_info.energy_threshold_arg+0.1));
  fast=args_info.fast_folding_arg;
  penalty=args_info.extension_cost_arg;
  threshloop=args_info.minimal_loop_energy_arg;
  threshRE=args_info.minimal_right_duplex_arg;
  threshLE=args_info.minimal_left_duplex_arg;
  threshDE=args_info.minimal_duplex_arg;
  distance=args_info.duplex_distance_arg;
  half_stem=args_info.minimal_stem_length_arg;
  max_half_stem=args_info.maximal_stem_length_arg;
  min_s2=args_info.maximal_duplex_box_length_arg;
  max_s2=args_info.maximal_duplex_box_length_arg;
  min_s1=args_info.minimal_snoRNA_stem_loop_length_arg;
  max_s1=args_info.maximal_snoRNA_stem_loop_length_arg;
  min_d1=args_info.minimal_snoRNA_duplex_length_arg;
  min_d2=args_info.minimal_snoRNA_duplex_length_arg;
  threshTE=args_info.minimal_duplex_stem_energy_arg;
  threshSE=args_info.minimal_total_energy_arg;
  max_asymm=args_info.maximal_stem_asymmetry_arg;
  threshD=args_info.minimal_lower_stem_energy_arg;
  
  
  threshloop=MIN2(threshloop,0);

/*   if(plfold_up_flag && !fast){ */
/*     printf("Sorry, no accessibility implementation with the non fast implementation\n"); */
/*     printf("If you want to include accessibility information please run RNAsnoop with the -f option\n"); */
/*     printf("If you want to run RNAsnoop with the slow algorithm, please remove run RNAsnoop without -P\n"); */
/*     exit(0); */
/*   } */
  
  
  if(plfold_up_flag==2 && suffix==NULL){
    printf("RNAsnoop needs a suffix (-S u1-to-30) for the RNAup accessibility file in order to localize them\n");
    exit(0);
  }
  if(redraw){
    if(!alignment){
      redraw_output(result_file, output_directory, plfold_up_flag, suffix, access,NULL,NULL);
    }
    else{
      redraw_output(result_file, output_directory, plfold_up_flag, suffix, access,sname,tname);
    }
  
    //readfile
    //parse lines
    //use current ploting function to output the results
    exit(0);
  }
  //istty = isatty(fileno(stdout))&&isatty(fileno(stdin));
   
  min_s2+=5;
  max_s2+=5;
  min_s1+=5;
  max_s1+=5;
  min_d1+=5;
  min_d2+=5;
  

  if(!alignment){
    if(tname==NULL || sname==NULL) RNAsnoop_cmdline_parser_print_help();
    sno=fopen(sname, "r");
    if(sno==NULL){printf("%s: Wrong snoRNA file name\n", sname);return 0;}
    mrna=fopen(tname, "r");
    if(mrna==NULL){printf("%s: Wrong target file name\n", tname);return 0;}
    do {                                /* main loop: continue until end of file */
      if ((line_s = get_line(sno))==NULL) {
        free(line_s);
        break;
      }
      
      //  /* skip comment lines and get filenames */
      while ((*line_s=='*')||(*line_s=='\0')||(*line_s=='>')) {
        if (*line_s=='>')
          name_s = (char *) space(strlen(line_s)+1);
        (void) sscanf(line_s,"%s",name_s);
        free(line_s);
        if ((line_s = get_line(sno))==NULL) {
          free(line_s);
          break;
        }
      }
      if(name_s == NULL) {printf("Your snoRNA sequence: \n%s\nhas no header. Please update your fasta file\n", line_s); exit(0);}
      //  if ((line ==NULL) || (strcmp(line, "@") == 0)) break;
      temp_s = (char *) space(strlen(line_s)+1);
      (void) sscanf(line_s,"%s",temp_s);
      free(line_s);
      length_s = (int) strlen(temp_s);
      for (l = 0; l < length_s; l++) {
        temp_s[l] = toupper(temp_s[l]);
        if (!noconv && temp_s[l] == 'T') temp_s[l] = 'U';
      }
      string_s = (char*) space(length_s + 11);
      strcpy(string_s, "NNNNN");
      strcat(string_s+5, temp_s);
      strcat(string_s+5+length_s, "NNNNN\0");
      free(temp_s); 
      snofold(string_s, max_asymm, threshloop, min_s2, max_s2, half_stem, max_half_stem);
      do{                                /* main loop for target continue until end of file */
        snoopT mfe;
        if ((line_t = get_line(mrna))==NULL) {
          //free(line_t);
          break;
        }
        char *name_t=NULL;
        
        // //  /* skip comment lines and get filenames */
        while ((*line_t=='*')||(*line_t=='\0')||(*line_t=='>')) {
          if (*line_t=='>'){
            printf("%s\n", name_s);
            name_t = (char*) space(strlen(line_t)+1);
            (void) sscanf(line_t,"%s",name_t);
            
            printf("%s\n", name_t);
            free(line_t);
            //            free(string_t);
              } 
          if ((line_t = get_line(mrna))==NULL) {
            free(line_t);
            break;
          }
        }
        if(name_t == NULL) {printf("Your target sequence: \n%s\nhas no header. Please update your fasta file\n", line_t); exit(0);}
        //  if ((line ==NULL) || (strcmp(line, "@") == 0)) break;
        temp_t = (char *) space(strlen(line_t)+1);
        (void) sscanf(line_t,"%s",temp_t);
        int length_t;
        length_t = (int) strlen(temp_t);
        for (l = 0; l < length_t; l++) {
          temp_t[l] = toupper(temp_t[l]);
          if (!noconv && temp_t[l] == 'T') temp_t[l] = 'U';
        }
        string_t =(char *) space(length_t +11);
        strcpy(string_t, "NNNNN");
        strcat(string_t+5, temp_t);
        strcat(string_t+5+length_t, "NNNNN");
        free(temp_t);
        char *name_output;
        name_output=NULL;
        if(nice){
          name_output = (char*) space(sizeof(char)*(length_t+length_s+2));
          strcpy(name_output, name_t+1);
          strcat(name_output, "_");
          strcat(name_output, name_s+1);
          name_output[length_t + length_s +1]='\0';
        }
        if(delta>=0){
          snoopT* subopt;
          snoopT* sub;
          if(!fast && !plfold_up_flag){
            subopt = snoop_subopt(string_t, string_s, delta, 5, penalty, threshloop, 
                                  threshLE,threshRE, threshDE,threshTE,threshSE,threshD,distance,
                                  half_stem, max_half_stem, min_s2, max_s2, min_s1, max_s1, min_d1, min_d2);
            if(subopt==NULL){printf("no target found under the given constraints\n");free(subopt); free(string_t);free(line_t);continue;}
            int count=0;
            for (sub=subopt; sub->structure !=NULL; sub++) {
              print_struc(sub, string_t, string_s, name_s, name_t, count++,nice);
              free(sub->structure);
            }
            free(subopt);
            free(string_t);
            free(line_t);
          }
          else if(!plfold_up_flag){
            Lsnoop_subopt_list(string_t, string_s, delta, 5, penalty, threshloop, 
                               threshLE,threshRE,threshDE,threshTE,threshSE,threshD,distance,
                               half_stem, max_half_stem, min_s2, max_s2, min_s1, max_s1, min_d1, min_d2, alignment_length,name_output);
            free(string_t);
            free(line_t);
          }
          else {
            if(plfold_up_flag==1){
              int **access_s1; 
              char *file_s1;
              int s1_len;//k;//,j;
              s1_len=strlen(string_t); 
              file_s1 = (char *) space(sizeof(char) * (strlen(name_t+1)+strlen(access)+9));
              strcpy(file_s1, access);
              strcat(file_s1, "/");
              strcat(file_s1, name_t+1);
              strcat(file_s1, "_openen");
              access_s1 = read_plfold_i(file_s1, 1, s1_len);
              if(fast) { 
                Lsnoop_subopt_list_XS(string_t, string_s, (const int**) access_s1, delta, 5, penalty, threshloop,  
                                      threshLE,threshRE,threshDE,threshTE,threshSE,threshD,distance, 
                                      half_stem, max_half_stem, min_s2, max_s2, min_s1, max_s1, min_d1, min_d2, alignment_length,name_output); 
              }
              else{
                snoop_subopt_XS(string_t, string_s, (const int **) access_s1, delta, 5, penalty, threshloop,  
                                threshLE,threshRE,threshDE,threshTE,threshSE,threshD,distance, 
                                half_stem, max_half_stem, min_s2, max_s2, min_s1, max_s1, min_d1, min_d2, alignment_length,name_output); 
              }
/*                for(j=0; j<s1_len; j++){  */
/*                  for(k=0; k<access_s1[0][0];k++){  */
/*                    printf("%d \t",access_s1[k][j]);  */
/*                  }  */
/*                  printf("\n");  */
/*                }  */
              free(file_s1);free(string_t);free(line_t);
              int i = access_s1[0][0];
              while(--i>-1){
                free(access_s1[i]);
              }
              free(access_s1);
            }
            else if(plfold_up_flag==2){
              int **access_s1; 
              char *file_s1;
              int s1_len,k;//,j;
              s1_len=strlen(string_t); 
              file_s1 = (char *) space(sizeof(char) * (strlen(name_t+1)+strlen(suffix) + strlen(access)+3));
              strcpy(file_s1, access);
              strcat(file_s1, "/");
              strcat(file_s1, name_t+1);
              strcat(file_s1, "_");
              strcat(file_s1, suffix);
              access_s1 = read_rnaup(file_s1, 1, s1_len);
              if(fast){
                Lsnoop_subopt_list_XS(string_t, string_s, (const int**) access_s1, delta, 5, penalty, threshloop,  
                                      threshLE,threshRE,threshDE,threshTE,threshSE,threshD,distance, 
                                      half_stem, max_half_stem, min_s2, max_s2, min_s1, max_s1, min_d1, min_d2, alignment_length, name_output); 
              }
              else{
                snoop_subopt_XS(string_t, string_s, (const int**) access_s1, delta, 5, penalty, threshloop,
                                threshLE,threshRE,threshDE,threshTE,threshSE,threshD,distance,
                                half_stem, max_half_stem, min_s2, max_s2, min_s1, max_s1, min_d1, min_d2, alignment_length,name_output);
              }
                

/*                for(j=0; j<s1_len; j++){  */
/*                  for(k=0; k<access_s1[0][0];k++){  */
/*                    printf("%d \t",access_s1[k][j]);  */
/*                  }  */
/*                  printf("\n");  */
/*                }  */
              free(file_s1);free(string_t);free(line_t);
              k = access_s1[0][0];
              while(--k>-1){
                free(access_s1[k]);
              }
              free(access_s1);
            }
          }
        }
        else{
          mfe=snoopfold(string_t, string_s, 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 < INF){
            print_struc(&mfe, string_t, string_s, name_s, name_t,0,1);
            free(mfe.structure);
          }
          free(line_t);
          length_t= (int) strlen(string_t);
          free(string_t);
        }
        free(name_t);
        if(nice){
          free(name_output);
        }
      } while(1);
      rewind(mrna);
      snofree_arrays(strlen(string_s));  /* free's base_pair */
      free(string_s);
      free(name_s);
    } while (1);
  }
  else{
    if(tname==NULL || sname==NULL) RNAsnoop_cmdline_parser_print_help();
    sno=fopen(sname, "r");
    if(sno==NULL){printf("%s: Wrong snoRNA file name\n", sname);return 0;}
    mrna=fopen(tname, "r");
    if(mrna==NULL){printf("%s: Wrong target file name\n", tname);return 0;}
    char *temp1[MAX_NUM_NAMES], *temp2[MAX_NUM_NAMES],*AS1[MAX_NUM_NAMES], *AS2[MAX_NUM_NAMES], *names1[MAX_NUM_NAMES], *names2[MAX_NUM_NAMES];
    int n_seq, n_seq2;
    n_seq =  read_clustal(mrna, temp1, names1);
    n_seq2 = read_clustal(sno, temp2, names2);
    if (n_seq != n_seq2){
      for (i=0; temp1[i]; i++) {
        free(temp1[i]); 
        free(temp2[i]); 
      }  
      nrerror("unequal number of seqs in alignments");
    }
    for(i=0;temp1[i];i++){
      AS1[i] = (char*) space((strlen(temp1[i])+11)*sizeof(char));
      AS2[i] = (char*) space((strlen(temp2[i])+11)*sizeof(char));
      strcpy(AS1[i],"NNNNN");
      strcpy(AS2[i],"NNNNN");
      strcat(AS1[i],temp1[i]);
      strcat(AS2[i],temp2[i]);
      strcat(AS1[i],"NNNNN\0");
      strcat(AS2[i],"NNNNN\0");
    }
    for (i=0; temp1[i]; i++) {
      free(temp1[i]); 
      free(temp2[i]); 
    }  
    AS1[n_seq]=NULL;
    AS2[n_seq]=NULL;
    update_fold_params();
    alisnofold((const char **)AS2, max_asymm, threshloop, min_s2, max_s2, half_stem, max_half_stem);
     snoopT struc; 
     struc = alisnoopfold((const char**) AS1, (const char**) AS2,  
                  penalty, threshloop,   
                          threshLE,threshRE,threshDE,threshD, 
                  half_stem, max_half_stem,  
                  min_s2, max_s2, min_s1, 
                  max_s1, min_d1, min_d2); 
     if(!(struc.structure==NULL)){ 
       aliprint_struc(&struc, (const char**) AS1,(const char**) AS2, names1, names2,0,nice); 
       free(struc.structure); 
     }
    


    snoopT* subopt = alisnoop_subopt((const char**) AS1, (const char**) AS2, delta, 5, penalty, threshloop,                              
                                      threshLE,threshRE, threshDE, threshD,
                                      threshTE,threshSE,distance,                     
                                      half_stem, max_half_stem, min_s2, max_s2,  
                                      min_s1, max_s1, min_d1, min_d2); 
    

    snoopT *sub; 
    if(subopt==NULL){
      printf("no target found under the given constraints\n");
    }else{
      int count=1;
      for (sub=subopt; !(sub->structure ==NULL); sub++) { 
        aliprint_struc(sub, (const char**) AS1,(const char**) AS2,names1,names2,count,nice); 
        free(sub->structure); 
        count++;
      } 
    }

     alisnofree_arrays(strlen(AS2[0])); 
     free(subopt); 
    for (i=0; temp1[i]; i++) {
      free(AS1[i]); 
      free(AS2[i]); 
      free(names1[i]);
      free(names2[i]);
    }  
    
  } 
  fclose(sno);
  fclose(mrna);
  return 0;
}

static void print_struc(snoopT  *dup, const char *s1, const char *s2, char *name_s, char *name_t, int count, int nice) {
  int l1;
  l1 = strchr(dup->structure, '&')-dup->structure;
  char *target_struct;
  int shift=0, n2;
  char* psoutput;
  psoutput = (char*) space(100*sizeof(char));
/*   if(dup->i > strlen(s1)-10){ */
/*         dup->i--; */
/*         l1--; */
/*   } */
/*   if(dup->i-l1< 0){ */
/*         l1--; */
/*         shift++; */
/*   } */
  target_struct = (char*) space(sizeof(char) * (strlen(dup->structure)+1));
  strncpy(target_struct, dup->structure+shift, l1);
  strncat(target_struct, dup->structure + (strchr(dup->structure, '&')-dup->structure), strlen(dup->structure) - (strchr(dup->structure, '&')-dup->structure));
  strcat(target_struct,"\0");
  char *target;
  target = (char *) space(l1+1);
  strncpy(target, (s1+dup->i-l1+5), l1);
  target[l1]='\0';
  char *s4;
  n2 = strlen(s2);
  s4 = (char*) space(sizeof(char) *(n2-9));
  strncpy(s4, s2+5, n2-10);
  s4[n2-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, dup->i+1-l1,
         dup->i, dup->u, dup->j+1, dup->j + (strrchr(dup->structure,'>') - strchr(dup->structure, '>'))+1, 
         (dup->Loop_D + dup->Duplex_El + dup->Duplex_Er + dup->Loop_E) + 4.10, 
         dup->Duplex_El, dup->Duplex_Er, dup->Loop_E, dup->Loop_D,target,s4);
  if(nice){
    char *temp_seq;
    char *temp_struc;
    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];
    char *temp; 
    int length_name = strlen(name_t)+strlen(name_s);
    temp = (char *) space(sizeof(char)*(length_name+1));
    strcpy(temp,name_t+1);
    strcat(temp,"_");
    strcat(temp,name_s + 1);
    temp[length_name] = '\0';
    strcpy(psoutput,"sno_");
    sprintf(str,"%d",count);
    strcat(psoutput,str);
    sprintf(upos,"%d",dup->u);
    strcat(psoutput,"_u_");
    strcat(psoutput,upos);
    strcat(psoutput,"_");
    strcat(psoutput,temp);
    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);free(psoutput);free(temp);
  }
free(s4);
free(target_struct);
free(target);
}


static void aliprint_struc(snoopT  *dup, const char **s1, const char **s2, char **name_t, char **name_s,int count,int nice) {
  int l1;
  l1 = strchr(dup->structure, '&')-dup->structure;
  char *target_struct;
  int shift=0;
/*   if(dup->i > strlen(s1[0])-10){ */
/*     dup->i--; */
/*     l1--; */
/*   } */
/*   if(dup->i-l1< 0){ */
/*     l1--; */
/*     shift++; */
/*   } */
  int length_struct = strlen(dup->structure);
  target_struct = (char*) space(sizeof(char) * (length_struct+1));
  strncpy(target_struct, dup->structure+shift, l1);
  strncat(target_struct, dup->structure + (strchr(dup->structure, '&')-dup->structure), length_struct - (strchr(dup->structure, '&')-dup->structure));
  strcat(target_struct,"\0");
  //get the corresponding alignment slice
  int n_seq,s;
  for(s=0; s1[s]!=NULL; s++);
  n_seq=s;
  int n1,n2;
  n1=strlen(s1[0]);n2=strlen(s2[0]);  
  char **target;
  target = (char**) space((n_seq+1)*sizeof(char*));
  for(s=0;s<n_seq;s++){
    target[s] = (char *) space((l1 + n2-8)*sizeof(char));
    strncpy(target[s], (s1[s]+dup->i-l1+5), l1);
    strcat(target[s],"&");
    strncat(target[s], (s2[s]+5), n2-10);
    target[s][l1+n2-9]='\0';
  }
  char * consens;
  consens =  consens_mis((const char**)target);  
  consens[l1]='&';

  printf("%s %3d,%-3d;%3d : %3d,%-3d (%5.2f = %5.2f + %5.2f + %5.2f + %5.2f + 4.1; duplex cov = %5.2f; stem cov = %5.2f )\n%s\n",  
          dup->structure, dup->i+1-l1, 
          dup->i, dup->u, dup->j+1, dup->j + (strrchr(dup->structure,'>') - strchr(dup->structure, '>'))+1,
          (dup->Loop_D + dup->Duplex_El + dup->Duplex_Er + dup->Loop_E)/n_seq + 4.10,  
          dup->Duplex_El/n_seq, dup->Duplex_Er/n_seq, dup->Loop_E/n_seq, dup->Loop_D/n_seq,dup->pscd/n_seq, dup->psct/n_seq,consens ); 
  if(nice){
    char* psoutput;
    psoutput = (char*) space(100*sizeof(char));
    
    char *temp_seq, *temp_struct, **temp_target;  
    temp_seq    = (char*) space((l1 + n2 -9)*sizeof(char));
    temp_struct = (char*) space((l1 + n2 -9)*sizeof(char));
    temp_target = (char**) space((n_seq+1)*sizeof(char*));
    for(s=0;s<n_seq;s++){
      temp_target[s] = (char *) space((l1 + n2-9)*sizeof(char));
      strncpy(temp_target[s], (s1[s]+dup->i-l1+5), l1);
      strncat(temp_target[s], (s2[s]+5), n2-10);
      temp_target[s][l1+n2-10]='\0';
    }
    strncpy(temp_seq, consens, l1);
    strncpy(temp_struct, target_struct, l1);
    strcat(temp_seq, consens+l1+1);
    strcat(temp_struct, target_struct+l1+1);
    temp_seq[n2-10+l1]='\0';
    temp_struct[n2-10+l1]='\0';
    char str[16];
    char upos[16];
    char *temp; 
    int length_name = strlen(name_t[0]) + strlen(name_s[0])+1;
    temp = (char *) space(sizeof(char)*(length_name+1));
    strcpy(temp,name_t[0]);
    strcat(temp,"_");
    strcat(temp, name_s[0]);
    temp[length_name] = '\0';
    strcpy(psoutput,"snoaln_");
    sprintf(str,"%d",count);
    strcat(psoutput,str);
    sprintf(upos,"%d",dup->u);
    strcat(psoutput,"_u_");
    strcat(psoutput,upos);
    strcat(psoutput,"_");
    for(s=0; s<length_name; s++){
      if(temp[s]=='/') temp[s]='-';
    }
    strcat(psoutput,temp);
    strcat(psoutput,".ps\0");
    //  psoutput[strlen(temp)+4+strlen(str)+39]='\0';
    cut_point = l1+1;
    aliPS_color_aln(dup->structure, psoutput, (const char**) target, (const char**) name_t);  
    psoutput[1]='t';
    psoutput[2]='r';
    PS_rna_plot_snoop_a(temp_seq, temp_struct, psoutput, NULL, (const char**) temp_target);
    cut_point = -1;
    free(psoutput);
    for(s=0; s<n_seq; s++){
      free(temp_target[s]);
    } 
    free(temp);
    free(temp_seq);
    free(temp_struct);
    free(temp_target);
  }
  for(s=0; s<n_seq; s++){
    free(target[s]);
  }
  free(consens);
  free(target_struct);
  free(target);
}

static int get_max_u(const char *s, char delim){
  char * pch;
  int total;
  total=0;
  pch = strchr(s, delim);
  total++;
  while(pch!=NULL){
    pch=strchr(pch+1, delim);
    total++;
  }
  return total-2;
}



static int ** read_rnaup(char *fname, const int beg, const int end)
{
  FILE *in = fopen(fname,"r"); // open RNA_up file
        
  int i,j;
  int **access;
  int offset, temp;
  temp=0; offset=0;
  int seq_pos;
  int beg_r, end_r;
  beg_r=beg;
  end_r=end;

  if(in==NULL){
    printf("%s", fname);
    perror("RNAup File open error here, Computing next target");
 
    exit(EXIT_FAILURE);
  }
  char tmp[2048]={0x0};
  //char *ptr;


  if(fgets(tmp,sizeof(tmp),in)==0){
    perror("Empty File");
  }
  if(strchr(tmp,'>')){
    fprintf(stderr,"file %s is not in RNAup format\n",fname);
    exit(EXIT_FAILURE);
  }
  
  while(!strstr(fgets(tmp,sizeof(tmp),in), "pos")){};
 /*  (void) fgets(tmp,sizeof(tmp),in); //Datum  */
/*   (void) fgets(tmp,sizeof(tmp),in); //white line */
/*   (void) fgets(tmp,sizeof(tmp),in); //RNAup  */
/*   (void) fgets(tmp,sizeof(tmp),in); //sequence length */
/*   (void) fgets(tmp,sizeof(tmp),in); //sequence */

  int dim_x;
  dim_x = get_max_u(tmp,'S'); // get unpaired regions by conting tabs in second line
  access = (int**) space(sizeof(int *) * (dim_x+2));
  for(i=0; i< dim_x+2; i++){
    access[i] =(int *) space(sizeof(int) * (end_r-beg_r+7));
  }
  for(i=0;i<end_r -  beg_r +6;i++){
    for(j=0;j<dim_x+2;j++){
      access[j][i]=INF;
    }
  }
  access[0][0]=dim_x+2;
  while(fgets(tmp,sizeof(tmp),in)!=0 && --end_r > -1) /* read a record */
  {
    float n;
    //int i;
    //int u;
    beg_r--;
    if(beg_r<1){ 
      if(sscanf(tmp,"%d%n",&seq_pos,&temp)==1){ // read sequenz pos = 1. spalte
        offset+=temp;
        int count;
        count=1;
        while(sscanf(tmp+offset,"%f%n",&n,&temp)==1){ //read columns
          offset+=temp;
          access[count][-beg_r+5+1]=(int) rint(100*n); //seq_pos+5
          //          printf("%d %d %f\n", count, -beg_r, access[count][-beg_r]);
          count++;
        }
        offset=0;
      }
    }
  }
  fclose(in);
  return access;
}







static int ** read_plfold_i(char *fname, const int beg, const int end)
{
  FILE *in=fopen(fname,"r");
  int i,j;
  int **access;
  int offset,temp;
  temp=0;offset=0;
  int seq_pos;
  int beg_r, end_r;
  beg_r=beg;
  end_r=end;
  
  if(in==NULL){
    perror(" open error");
    exit(EXIT_FAILURE);
  }
  
  char tmp[2048]={0x0};
  //char *ptr;
  if(fgets(tmp,sizeof(tmp),in)==0){
    perror("Empty File");
  }
  if(strchr(tmp,'>')){
    fprintf(stderr,"file %s is not in RNAplfold format",fname);
    exit(EXIT_FAILURE);
  }
  if(fgets(tmp,sizeof(tmp),in)==0){
    perror("No accessibility data");
  }
  int dim_x;
  dim_x=get_max_u(tmp,'\t');
  access = (int**) space(sizeof(int *) * (dim_x+2));
  for(i=0; i< dim_x+2; i++){
    access[i] =(int *) space(sizeof(int) * (end_r-beg_r+1));
  }
  
  for(i=0;i<end_r -  beg_r +1;i++){
    for(j=0;j<dim_x+2;j++){
      access[j][i]=INF;
    }
  }
  access[0][0]=dim_x+2;
  while(fgets(tmp,sizeof(tmp),in)!=0 && --end_r > -1) /* read a record */
    {
      float n;
      //int i;
      //int u;
      beg_r--;
      if(beg_r<1){
        if(sscanf(tmp,"%d\t%n",&seq_pos,&temp)==1){
          offset+=temp;
          int count;
          count=1;
          while(sscanf(tmp+offset,"%f%n",&n,&temp)==1 ){
            offset+=temp;
            access[count][seq_pos+5]=(int)  rint( 100 * n); //round the number
            count++;
          }
          offset=0;
        }
      }
      
    }
  fclose(in);
  return access;
}



static void redraw_output(char *fname, char *output, int plfold_up_flag, char *suffix, char *access, char *sname, char *tname)
{
  char *line;
  line=NULL;
  int count;
  short two_seq = 0;
  char *results;
  char *sequence;
  char *query;
  char *target;
  char *structure;
  char *pos;
  int posi;
  int begin=0, end=0, u=0;
  if(output==NULL){
    output = (char*) space(sizeof(char)*2);
    strcpy(output,".\0");
  }
  count=0;
  if(sname==NULL && tname==NULL){
    while((line=get_line(stdin))!=NULL) {
      count++;
      if(two_seq==0 && *line =='>'){
        query=(char*) space(strlen(line)+1);
        (void) sscanf(line,"%s",query);
        free(line);
        line=NULL;
        memmove(query, query+1, strlen(query));
        two_seq++;
      }
      else if(two_seq==1 && *line =='>'){
        target=(char*) space(strlen(line)+1);
        (void) sscanf(line,"%s",target);
        free(line);
        line=NULL;
        memmove(target, target+1, strlen(target));
        two_seq++;
      }
      else if(two_seq == 2) {
        int *relative_access;
        relative_access=NULL;
        printf("%s %s\n", target,query);
        if(strchr(line, '(') && strchr(line,'&') && strchr(line, '(')&& strchr(line, ',') && strchr(line,':') && strchr(line, '-')){
          results=line;
          int length;
          pos = strchr(results, ' ');
          posi = (int)  (pos - results);
          length = posi;
          structure = (char *) space((length+1) * sizeof(char));
          sscanf(results,"%s",structure); //parse structure
          char *line2;
          if((line2=get_line(stdin))!=NULL){
            sequence=(char *) space((length+1)* sizeof(char));
            sscanf(line2,"%s",sequence);
            if(line2!=NULL){
              free(line2);
            }
          }
          else{
            printf("your file is fucked up");exit(0);
          } //parse sequence
          
          sscanf(pos, "%10d,%10d;%10d", &begin,&end,&u); //parse coordinates
          if(plfold_up_flag){
            int **access_s1;
            char *file_s1;
            //read_rnaup_file
            file_s1 = (char*) space(sizeof(char) * (strlen(target)+strlen(suffix)+strlen(access)+3));
            strcpy(file_s1, access);
            strcat(file_s1, "/");
            strcat(file_s1, target);
            strcat(file_s1, "_");
            strcat(file_s1, suffix);
            access_s1 = read_rnaup(file_s1, begin, end);
            free(file_s1);
            relative_access = (int*) space(sizeof(int)*(end-begin+2));
            relative_access[0]=access_s1[1][1+5];
            int i;
            for(i=2;i<(end-begin+2);i++){
              relative_access[i-1]=access_s1[i+1][i+5] - access_s1[i][i+4];
            }
            int l = access_s1[0][0];
            while(--l>-1){
              free(access_s1[l]);
            }
            free(access_s1);
          }
          char *catseq, *catstruct,*output_file;
          catseq = (char*) space(strlen(sequence)  *sizeof(char));
          catstruct = (char*) space(strlen(structure)  *sizeof(char));
          int l1 = strchr(structure, '&')-structure;
          strncpy(catseq, sequence, l1);
          strcat(catseq, sequence+l1+1);
          strncpy(catstruct, structure, l1);
          strcat(catstruct, structure+l1+1);
          strcat(catseq,"\0");
          strcat(catstruct,"\0");
          
          //printf("%s\n%s\n%s\n%s", catseq,sequence,catstruct,structure);
          cut_point=l1+1;
          output_file = (char*) space((strlen(output) + strlen(query)+strlen(target)+50)*sizeof(char));
          strcpy(output_file,output);
          strcat(output_file, "/");
          strcat(output_file,"sno_");
          strcat(output_file,query);
          strcat(output_file,"_");
          strcat(output_file, target);
          strcat(output_file,"_u_");
          char str[9];
          sprintf(str,"%d",u);
          strcat(output_file,str);
          strcat(output_file,"_");
          sprintf(str,"%d",count);
          strcat(output_file,str);
          strcat(output_file,".ps\0");
          PS_rna_plot_snoop_a(catseq, catstruct, output_file,relative_access,NULL);
          if(relative_access){
            free(relative_access);
          }
          free(output_file);output_file=NULL;
          free(catseq);catseq=NULL;
          free(catstruct);catstruct=NULL;
          free(structure);
          free(sequence);     
          free(line);
          line=NULL;
        }
        else if(*line=='>'){
          free(query);free(target);
          two_seq=1;
          query = (char*) space(sizeof(char)*(strlen(line)+1));
          (void) sscanf(line,"%s",query); 
          free(line);
          memmove(query, query+1, strlen(query));
        }
      }
    }
    free(target);
    free(query);
  }
  else if(!(tname==NULL && sname==NULL)) {
    FILE* sno, *mrna;
    int i;
    sno=fopen(sname, "r");
    if(sno==NULL){printf("%s: Wrong snoRNA file name\n", sname);}
    mrna=fopen(tname, "r");
    if(mrna==NULL){printf("%s: Wrong target file name\n", tname);}
    char *temp1[MAX_NUM_NAMES], *temp2[MAX_NUM_NAMES], *AS1[MAX_NUM_NAMES], *AS2[MAX_NUM_NAMES], *names1[MAX_NUM_NAMES], *names2[MAX_NUM_NAMES];
    int n_seq, n_seq2;
    n_seq =  read_clustal(mrna, temp1, names1);
    n_seq2 = read_clustal(sno, temp2, names2);
    if (n_seq != n_seq2){
      for (i=0; temp1[i]; i++) {
        free(temp1[i]); 
        free(temp2[i]); 
      }  
      nrerror("unequal number of seqs in alignments");
    }
    for(i=0;temp1[i];i++){
      AS1[i] = (char*) space((strlen(temp1[i])+11)*sizeof(char));
      AS2[i] = (char*) space((strlen(temp2[i])+11)*sizeof(char));
      strcpy(AS1[i],"NNNNN");
      strcpy(AS2[i],"NNNNN");
      strcat(AS1[i],temp1[i]);
      strcat(AS2[i],temp2[i]);
      strcat(AS1[i],"NNNNN\0");
      strcat(AS2[i],"NNNNN\0");
    }
    for (i=0; temp1[i]; i++) {
      free(temp1[i]); 
      free(temp2[i]); 
    }  
    AS1[n_seq]=NULL;
    AS2[n_seq]=NULL;
    int count=0;
    while((line=get_line(stdin))!=NULL) {
      results=line;
      if(strchr(line, '(') && strchr(line,'&') && strchr(line, '(')&& strchr(line, ',') && strchr(line,':') && strchr(line, '-')) {
        count++;
        int length;
        //int sbegin, send;
        //int energy;
        snoopT dup;
        pos = strchr(results, ' ');
        posi = (int)  (pos - results);
        length = posi;
        dup.structure = (char *) space((length+1) * sizeof(char));
        sscanf(results,"%s %10d,%10d;%10d : %10d,%10d (%10f = %10f + %10f + %10f + %10f + 4.1; duplex cov = %10f; stem cov = %10f",
               dup.structure,
               &begin,
               &(dup.i),
               &(dup.u),
               &(dup.j),
               &end,
               &(dup.energy),
               &(dup.Duplex_El),
               &(dup.Duplex_Er),
               &(dup.Loop_E),
               &(dup.Loop_D),
               &(dup.pscd),
               &(dup.psct)); //parse duplex stuff;
        dup.energy*=n_seq;
        dup.Duplex_El*=n_seq;
        dup.Duplex_Er*=n_seq;
        dup.Loop_E*=n_seq;
        dup.Loop_D*=n_seq;
        aliprint_struc(&dup, (const char**) AS1, (const char**)AS2, names1, names2,count,1);
        free(dup.structure);
        free(line);
      }
      else{
        free(line);
      }
    }
    
    for (i=0; AS1[i]; i++) {
      free(AS1[i]); 
      free(AS2[i]); 
        free(names1[i]); 
      free(names2[i]);
    }  
    fclose(mrna);
    fclose(sno);
  }
  //  free(output);
}             
