/*
                               aln_util.c
               Helper functions frelated to alignments
*/

#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

#include <stdio.h>
#include <stdlib.h>
#include <errno.h>
#include <time.h>
#include <string.h>
#include <ctype.h>
#include <math.h>

#include "ViennaRNA/utils.h"
#include "ViennaRNA/fold_vars.h"
#include "ViennaRNA/pair_mat.h"
#include "ViennaRNA/model.h"
#include "ViennaRNA/ribo.h"
#include "ViennaRNA/aln_util.h"

#define MAX_NUM_NAMES    500
int read_clustal(FILE *clust, char *AlignedSeqs[], char *names[]) {
   char *line, name[100]="", *seq;
   int  n, nn=0, num_seq = 0, i;

   if ((line=vrna_read_line(clust)) == NULL) {
     vrna_message_warning("Empty CLUSTAL file"); return 0;
   }

   if ((strncmp(line,"CLUSTAL", 7) !=0) && (!strstr(line,"STOCKHOLM"))) {
     vrna_message_warning("This doesn't look like a CLUSTAL/STOCKHOLM file, sorry");
     free(line); return 0;
   }
   free(line);
   line = vrna_read_line(clust);

   while (line!=NULL) {
    if(strncmp(line, "//", 2) == 0){
      free(line);
      break;
    }

    if (((n=strlen(line))<4) || isspace((int)line[0])) {
      /* skip non-sequence line */
      free(line); line = vrna_read_line(clust);
      nn=0; /* reset seqence number */
      continue;
    }
    /* skip comments */
    if(line[0] == '#'){
      free(line);
      line = vrna_read_line(clust);
      continue;
    }

     seq = (char *) vrna_alloc( (n+1)*sizeof(char) );
     sscanf(line,"%99s %s", name, seq);

    for(i=0;i<strlen(seq);i++){
      if(seq[i] == '.') seq[i] = '-'; /* replace '.' gaps by '-' */
      /* comment the next line and think about something more difficult to deal with
         lowercase sequence letters if you really want to */
      seq[i] = toupper(seq[i]);
    }

     if (nn == num_seq) { /* first time */
       names[nn] = strdup(name);
       AlignedSeqs[nn] = strdup(seq);
     }
     else {
       if (strcmp(name, names[nn])!=0) {
         /* name doesn't match */
         vrna_message_warning("Sorry, your file is messed up (inconsitent seq-names)");
         free(line); free(seq);
         return 0;
       }
       AlignedSeqs[nn] = (char *)
         vrna_realloc(AlignedSeqs[nn], strlen(seq)+strlen(AlignedSeqs[nn])+1);
       strcat(AlignedSeqs[nn], seq);
     }
     nn++;
     if (nn>num_seq) num_seq = nn;
     free(seq);
     free(line);
     if (num_seq>=MAX_NUM_NAMES) {
       vrna_message_warning("Too many sequences in CLUSTAL/STOCKHOLM file");
       return 0;
     }

     line = vrna_read_line(clust);
   }

   AlignedSeqs[num_seq] = NULL;
   names[num_seq] = NULL;
   if (num_seq == 0) {
     vrna_message_warning("No sequences found in CLUSTAL/STOCKHOLM file");
     return 0;
   }
   n = strlen(AlignedSeqs[0]);
   for (nn=1; nn<num_seq; nn++) {
     if (strlen(AlignedSeqs[nn])!=n) {
       vrna_message_warning("Sorry, your file is messed up.\n"
                            "Unequal lengths!");
       return 0;
     }
   }

   vrna_message_info(stderr, "%d sequences; length of alignment %d.", nn, n);
   return num_seq;
}

char *consensus(const char *AS[]) {
  /* simple consensus sequence (most frequent character) */
  char *string;
  int i,n;

  string = NULL;

  if(AS){
    n = strlen(AS[0]);
    string = (char *) vrna_alloc((n+1)*sizeof(char));
    for (i=0; i<n; i++) {
      int s,c,fm, freq[8] = {0,0,0,0,0,0,0,0};
      for (s=0; AS[s]!=NULL; s++)
        freq[encode_char(AS[s][i])]++;
      for (s=c=fm=0; s<8; s++) /* find the most frequent char */
        if (freq[s]>fm) {c=s, fm=freq[c];}
      if (s>4) s++; /* skip T */
      string[i]=Law_and_Order[c];
    }
  }
  return string;
}

/* IUP nucleotide classes indexed by a bit string of the present bases */
/* A C AC G AG CG ACG U AU CU ACU GU AGU CGU ACGU */
static char IUP[17] = "-ACMGRSVUWYHKDBN";
char *consens_mis(const char*AS[]) {
  /* MIS displays the 'most informative sequence' (Freyhult et al 2004),
     elements in columns with frequency greater than the background
     frequency are projected into iupac notation. Columns where gaps are
     over-represented are in lower case. */

  char *cons;
  int i, s, n, N, c;
  int bgfreq[8] = {0,0,0,0,0,0,0,0};

  cons = NULL;

  if(AS){
    n = strlen(AS[0]);
    for (N=0; AS[N]!=NULL; N++);
    cons = (char *) vrna_alloc((n+1)*sizeof(char));

    for (i=0; i<n; i++)
      for (s=0; s<N; s++) {
        c = encode_char(AS[s][i]);
        if (c>4) c=5;
        bgfreq[c]++;
      }

    for (i=0; i<n; i++) {
      int freq[8] = {0,0,0,0,0,0,0,0};
      int code = 0;
      for (s=0; s<N; s++) {
        c = encode_char(AS[s][i]);
        if (c>4) c=5;
        freq[c]++;
      }
      for (c=4; c>0; c--) {
        code <<=1;
        if (freq[c]*n>=bgfreq[c]) code++;
      }
      cons[i] = IUP[code];
      if (freq[0]*n>bgfreq[0])
        cons[i] = tolower(IUP[code]);
    }
  }
  return cons;
}

PUBLIC char *
get_ungapped_sequence(const char *seq){

  char  *tmp_sequence, *b;
  int   i;

  tmp_sequence = strdup(seq);

  b = tmp_sequence;
  i = 0;
  do{
    if((*b=='-')||(*b=='_')||(*b=='~')||(*b=='.')) continue;
    tmp_sequence[i] = *b;
    i++;
  }while(*(++b));

  tmp_sequence = (char *)vrna_realloc(tmp_sequence, (i+1)*sizeof(char));
  tmp_sequence[i] = '\0';

  return tmp_sequence;
}

PUBLIC int
vrna_aln_mpi(const char **alignment){

  int   i, j, k, s, n_seq, n, pairnum = 0, sumident = 0;
  float ident = 0;

  if(alignment){
    n = (int)strlen(alignment[0]);
    for(s = 0; alignment[s]; s++);
    n_seq = s;

    for(j = 0; j < n_seq - 1; j++)
      for(k = j + 1; k < n_seq; k++) {
        ident = 0;
        for (i = 1; i <= n; i++){
          if(alignment[k][i] == alignment[j][i])
            ident++;
          pairnum++;
        }
        sumident+=ident;
      }

    if(pairnum > 0)
      return (int) (sumident*100/pairnum);
  }
  return 0;
}

/*---------------------------------------------------------------------------*/
PRIVATE int
compare_pinfo(const void *pi1,
              const void *pi2){

  vrna_pinfo_t *p1, *p2;
  int  i, nc1, nc2;
  p1 = (vrna_pinfo_t *)pi1;  p2 = (vrna_pinfo_t *)pi2;
  for (nc1=nc2=0, i=1; i<=6; i++) {
    if (p1->bp[i]>0) nc1++;
    if (p2->bp[i]>0) nc2++;
  }
  /* sort mostly by probability, add
     epsilon * comp_mutations/(non-compatible+1) to break ties */
  return (p1->p + 0.01*nc1/(p1->bp[0]+1.)) <
         (p2->p + 0.01*nc2/(p2->bp[0]+1.)) ? 1 : -1;
}

PUBLIC vrna_pinfo_t *
vrna_aln_pinfo( vrna_fold_compound_t *vc,
                const char *structure,
                double threshold){

  int i,j, num_p=0, max_p = 64;
  vrna_pinfo_t *pi;
  double *duck, p;
  short *ptable = NULL;

  short **S = vc->S;
  char **AS = vc->sequences;
  int n_seq = vc->n_seq;
  int n     = vc->length;
  int         *my_iindx = vc->iindx;
  FLT_OR_DBL  *probs    = vc->exp_matrices->probs;
  vrna_md_t   *md = &(vc->exp_params->model_details);

  max_p = 64; pi = vrna_alloc(max_p*sizeof(vrna_pinfo_t));
  duck =  (double *) vrna_alloc((n+1)*sizeof(double));
  if(structure)
    ptable = vrna_ptable(structure);

  for (i=1; i<n; i++)
    for (j=i+TURN+1; j<=n; j++) {
      if ((p=probs[my_iindx[i]-j])>=threshold) {
        duck[i] -=  p * log(p);
        duck[j] -=  p * log(p);

        int type, s;
        pi[num_p].i   = i;
        pi[num_p].j   = j;
        pi[num_p].p   = p;
        pi[num_p].ent = duck[i]+duck[j]-p*log(p);

        for (type=0; type<8; type++) pi[num_p].bp[type]=0;
        for (s=0; s<n_seq; s++) {
          type = md->pair[S[s][i]][S[s][j]];
          if(S[s][i]==0 && S[s][j]==0) type = 7; /* gap-gap  */
          if ((AS[s][i-1] == '-')||(AS[s][j-1] == '-')) type = 7;
          if ((AS[s][i-1] == '~')||(AS[s][j-1] == '~')) type = 7;
          pi[num_p].bp[type]++;
        }
        if(ptable)
          pi[num_p].comp = (ptable[i] == j) ? 1:0;

        num_p++;
        if (num_p>=max_p) {
          max_p *= 2;
          pi = vrna_realloc(pi, max_p * sizeof(vrna_pinfo_t));
        }
      }
    }
  free(duck);
  pi = vrna_realloc(pi, (num_p+1)*sizeof(vrna_pinfo_t));
  pi[num_p].i=0;
  qsort(pi, num_p, sizeof(vrna_pinfo_t), compare_pinfo );

  free(ptable);
  return pi;
}


PUBLIC int *
vrna_aln_pscore(const char  **alignment,
                vrna_md_t   *md){

  /* calculate co-variance bonus for each 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         i, j, k, l, s, n, n_seq, *indx, turn, max_span;
  float       **dm;
  vrna_md_t   md_default;
  int         *pscore;
  short       **S;

  int olddm[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 */};

  pscore = NULL;

  if(!md){
    vrna_md_set_default(&md_default);
    md = &md_default;
  }

  if(alignment){
    /* length of alignment */
    n = (int)strlen(alignment[0]);

    /* count number of sequences */
    for(s = 0; alignment[s]; s++);
    n_seq = s;

    /* make numeric encoding of sequences */
    S = (short **)vrna_alloc(sizeof(short *) * (n_seq + 1));
    for(s = 0; s < n_seq; s++){
      S[s] = vrna_seq_encode_simple(alignment[s], md);
    }

    indx  = vrna_idx_col_wise(n);

    turn    = md->min_loop_size;

    pscore = (int *)vrna_alloc(sizeof(int) * ((n+1)*(n+2)/2 + 2));

    if (md->ribo) {
      if (RibosumFile !=NULL) dm=readribosum(RibosumFile);
      else dm=get_ribosum(alignment, n_seq, n);
    }
    else { /*use usual matrix*/
      dm = vrna_alloc(7*sizeof(float*));
      for (i=0; i<7;i++) {
        dm[i] = vrna_alloc(7*sizeof(float));
        for (j=0; j<7; j++)
          dm[i][j] = (float) olddm[i][j];
      }
    }

    max_span = md->max_bp_span;
    if((max_span < turn+2) || (max_span > n))
      max_span = n;
    for (i=1; i<n; i++) {
      for (j=i+1; (j<i+turn+1) && (j<=n); j++)
        pscore[indx[j]+i] = NONE;
      for (j=i+turn+1; j<=n; j++) {
        int pfreq[8]={0,0,0,0,0,0,0,0};
        double score;
        for (s=0; s<n_seq; s++) {
          int type;
          if (S[s][i]==0 && S[s][j]==0) type = 7; /* gap-gap  */
          else {
            if ((alignment[s][i] == '~')||(alignment[s][j] == '~')) type = 7;
            else type = md->pair[S[s][i]][S[s][j]];
          }
          pfreq[type]++;
        }
        if (pfreq[0]*2+pfreq[7]>n_seq) { pscore[indx[j]+i] = NONE; continue;}
        for (k=1,score=0; k<=6; k++) /* ignore pairtype 7 (gap-gap) */
          for (l=k; l<=6; l++)
            score += pfreq[k]*pfreq[l]*dm[k][l];
        /* counter examples score -1, gap-gap scores -0.25   */
        pscore[indx[j]+i] = md->cv_fact *
          ((UNIT*score)/n_seq - md->nc_fact*UNIT*(pfreq[0] + pfreq[7]*0.25));

        if((j - i + 1) > max_span){
          pscore[indx[j]+i] = NONE;
        }
      }
    }

    if (md->noLP) /* remove unwanted pairs */
      for (k=1; k<n-turn-1; k++)
        for (l=1; l<=2; l++) {
          int type,ntype=0,otype=0;
          i=k; j = i+turn+l;
          type = pscore[indx[j]+i];
          while ((i>=1)&&(j<=n)) {
            if ((i>1)&&(j<n)) ntype = pscore[indx[j+1]+i-1];
            if ((otype<md->cv_fact*MINPSCORE)&&(ntype<md->cv_fact*MINPSCORE))  /* too many counterexamples */
              pscore[indx[j]+i] = NONE; /* i.j can only form isolated pairs */
            otype =  type;
            type  = ntype;
            i--; j++;
          }
        }

    /*free dm */
    for (i=0; i<7;i++) {
      free(dm[i]);
    }
    free(dm);

    for(s = 0; s < n_seq; s++){
      free(S[s]);
    }
    free(S);

    free(indx);
  }

  return pscore;
}

/*###########################################*/
/*# deprecated functions below              #*/
/*###########################################*/

PUBLIC int
get_mpi(char *Alseq[],
        int n_seq,
        int length,
        int *mini){

  int   i, j, k, pairnum = 0, sumident = 0;
  float ident = 0, minimum = 1.;

  for(j=0; j<n_seq-1; j++)
    for(k=j+1; k<n_seq; k++) {
      ident=0;
      for (i=1; i<=length; i++){
        if (Alseq[k][i]==Alseq[j][i]) ident++;
        pairnum++;
      }
      if ((ident/length)<minimum) minimum=ident/(float)length;
      sumident+=ident;
    }
  mini[0]=(int)(minimum*100.);
  if (pairnum>0)   return (int) (sumident*100/pairnum);
  else return 0;
}

PUBLIC void
alloc_sequence_arrays(const char **sequences,
                      short ***S,
                      short ***S5,
                      short ***S3,
                      unsigned short ***a2s,
                      char ***Ss,
                      int circ){

  unsigned int s, n_seq, length;
  if(sequences[0] != NULL){
    length = strlen(sequences[0]);
    for (s=0; sequences[s] != NULL; s++);
    n_seq = s;
    *S    = (short **)          vrna_alloc((n_seq+1) * sizeof(short *));
    *S5   = (short **)          vrna_alloc((n_seq+1) * sizeof(short *));
    *S3   = (short **)          vrna_alloc((n_seq+1) * sizeof(short *));
    *a2s  = (unsigned short **) vrna_alloc((n_seq+1) * sizeof(unsigned short *));
    *Ss   = (char **)           vrna_alloc((n_seq+1) * sizeof(char *));
    for (s=0; s<n_seq; s++) {
      if(strlen(sequences[s]) != length) vrna_message_error("uneqal seqence lengths");
      (*S5)[s]  = (short *)         vrna_alloc((length + 2) * sizeof(short));
      (*S3)[s]  = (short *)         vrna_alloc((length + 2) * sizeof(short));
      (*a2s)[s] = (unsigned short *)vrna_alloc((length + 2) * sizeof(unsigned short));
      (*Ss)[s]  = (char *)          vrna_alloc((length + 2) * sizeof(char));
      (*S)[s]   = (short *)         vrna_alloc((length + 2) * sizeof(short));
      encode_ali_sequence(sequences[s], (*S)[s], (*S5)[s], (*S3)[s], (*Ss)[s], (*a2s)[s], circ);
    }
    (*S5)[n_seq]  = NULL;
    (*S3)[n_seq]  = NULL;
    (*a2s)[n_seq] = NULL;
    (*Ss)[n_seq]  = NULL;
    (*S)[n_seq]   = NULL;
  }
  else vrna_message_error("alloc_sequence_arrays: no sequences in the alignment!");
}

PUBLIC void
free_sequence_arrays( unsigned int n_seq,
                      short ***S,
                      short ***S5,
                      short ***S3,
                      unsigned short ***a2s,
                      char ***Ss){

  unsigned int s;
  for (s=0; s<n_seq; s++) {
    free((*S)[s]);
    free((*S5)[s]);
    free((*S3)[s]);
    free((*a2s)[s]);
    free((*Ss)[s]);
  }
  free(*S);   *S    = NULL;
  free(*S5);  *S5   = NULL;
  free(*S3);  *S3   = NULL;
  free(*a2s); *a2s  = NULL;
  free(*Ss);  *Ss   = NULL;
}

PUBLIC void
encode_ali_sequence(const char *sequence,
                    short *S,
                    short *s5,
                    short *s3,
                    char *ss,
                    unsigned short *as,
                    int circular){

  unsigned int i,l;
  unsigned short p;
  l     = strlen(sequence);
  S[0]  = (short) l;
  s5[0] = s5[1] = 0;

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

  if (oldAliEn){
    /* use alignment sequences in all energy evaluations */
    ss[0]=sequence[0];
    for(i=1; i<l; i++){
      s5[i] = S[i-1];
      s3[i] = S[i+1];
      ss[i] = sequence[i];
      as[i] = i;
    }
    ss[l]   = sequence[l];
    as[l]   = l;
    s5[l]   = S[l-1];
    s3[l]   = 0;
    S[l+1]  = S[1];
    s5[1]   = 0;
    if (circular) {
      s5[1]   = S[l];
      s3[l]   = S[1];
      ss[l+1] = S[1];
    }
  }
  else{
    if(circular){
      for(i=l; i>0; i--){
        char c5;
        c5 = sequence[i-1];
        if ((c5=='-')||(c5=='_')||(c5=='~')||(c5=='.')) continue;
        s5[1] = S[i];
        break;
      }
      for (i=1; i<=l; i++) {
        char c3;
        c3 = sequence[i-1];
        if ((c3=='-')||(c3=='_')||(c3=='~')||(c3=='.')) continue;
        s3[l] = S[i];
        break;
      }
    }
    else  s5[1]=s3[l]=0;

    for(i=1,p=0; i<=l; i++){
      char c5;
      c5 = sequence[i-1];
      if ((c5=='-')||(c5=='_')||(c5=='~')||(c5=='.'))
        s5[i+1]=s5[i];
      else { /* no gap */
        ss[p++]=sequence[i-1]; /*start at 0!!*/
        s5[i+1]=S[i];
      }
      as[i]=p;
    }
    for (i=l; i>=1; i--) {
      char c3;
      c3 = sequence[i-1];
      if ((c3=='-')||(c3=='_')||(c3=='~')||(c3=='.'))
        s3[i-1]=s3[i];
      else
        s3[i-1]=S[i];
    }
  }
}