#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#define MAXLEN 100
#define UNPAIRED (-1)
#define UNDEF (-2)
#define NEG (-10*MAXLEN)

#define SEQtMit "TAGGATGGGGTGTGATAGGTGGCACGGAGAATTTTGGATTCTCAGGGATGGGTTCGATTCTCATAGTCCTAG"
#define SEQtC17 "GGTTCCATGGTGTAATGGTTAGCACTCTGGACTCTGAATCCAGCGATCCGAGTTCAAATCTCGGTGGAACCT"
#define SEQ SEQtC17


// Global: the matrices W and K
int W[MAXLEN][MAXLEN];
int K[MAXLEN][MAXLEN];


// Score function
inline int nussinov_score(char a, char b)
{ // A-T/U=2, C-G=3, G-T/U=1.
  a=toupper(a);
  b=toupper(b);
  switch(a)
    {
    case 'A': 
      switch(b)
        {
        case 'T': case 'U': return 2; break;
        }
      break;
    case 'C': 
      switch(b)
        {
        case 'G': return 3; break;
        }
      break;
    case 'G': 
      switch(b)
        {
        case 'C': return 3; break;
        case 'T': case 'U': return 1; break;
        }
      break;
    case 'T': case 'U':
      switch(b)
        {
        case 'A': return 2; break;
        case 'G': return 1; break;
        }
      break;
    }
  return(NEG);
}


// Nussinov - Algorithm
int nussinov(const char* seq, int d)
{
  // seq: Sequence (Treat T like U)
  // d  : minimal Hairpin loop length
  //
  // Computes global matrices
  // W  : W[i][j] = best score on seq[i..j]
  // K  : K[i][j] = optimal pair of j on seq[i..j]
  int l,L;  // l iterates over interval length of [i..j]
  int i,j;  // j:=i+l-1
  int w,bw;
  int k,bk;

  L=strlen(seq);
  if(L<d+2) d=L-2;

  // Initialize
  // Set all to undef
  for(i=0; i<L; i++) for(j=0; j<L; j++) 
    { W[i][j]=UNDEF; K[i][j]=UNDEF; }
  // Iniitiialize diagonal belonging to length 0
  for(i=1; i<L; i++) 
    { W[i][i-1]=0; K[i][i-1]=UNPAIRED; }
  // Initialize all diagonals belonging to length 1..d-1
  for(l=1; l<d+2; l++) 
    for(i=0,j=i+l-1;  j<L;  i++,j++)
      { W[i][j]=0; K[i][j]=UNPAIRED; }

  // Fill matrices diagonally, starting with interval length l=d+2
  for(l=d+2; l<=L; l++)
    for(i=0,j=i+l-1;  j<L;  i++,j++)
      {
        bw=W[i][j-1]; bk=UNPAIRED; // first attempt case I
        for(k=i; k<j-d; k++) // attempt to pair k-j for all k=i..j-d-1
          {
            w=W[i][k-1]+W[k+1][j-1]+nussinov_score(seq[k],seq[j]);
            if(w>bw) {bw=w; bk=k;}
          }
        W[i][j]=bw; K[i][j]=bk; // store best w and best k
      }

  // return score for whole sequence
  return(W[0][L-1]);
}


int print_pairs(const char* seq, int i, int j, char* brackets)
     // recursively print all basepairs in interval [i,j]
     // uses matrix K
     // returns total number of pairs in [i,j]
     // also sets '(',')' in brackets for each pair
{
  int k;
  char sk,sj;
  int p;

  if(W[i][j]<=0 || j<=i) return 0;
  switch((k=K[i][j]))
    {
    case UNPAIRED: 
      return(print_pairs(seq,i,j-1,brackets));
      break;
    default: // k in [i,j-1]
      p=1;
      p+=print_pairs(seq,i,k-1,brackets);
      sk=seq[k]; sj=seq[j];
      printf("%2d %c-%c %2d   (%2d)\n",k,sk,sj,j,nussinov_score(sk,sj));
      brackets[k]='('; brackets[j]=')';
      p+=print_pairs(seq,k+1,j-1,brackets);
    }
  return p;
}



int main(void)
{
  char seq[] = SEQ;
  char brackets[] = SEQ;
  int  L = strlen(seq);
  int score;
  int pairs;
  int i;

  // Compute Score
  score=nussinov(seq,3);
  printf("Sequence (length=%d):\n%s\n\n",strlen(seq),seq);
  printf("Score=%d\nPairs:\n",score);

  // Traceback
  for(i=0;i<L;i++) brackets[i]='.';
  pairs=print_pairs(seq,0,strlen(seq)-1,brackets);
  printf("\nTotal: %d pairs\n",pairs);
  printf("%s\n%s\n\n",seq,brackets);
  return 0;
}


//Total: 27 pairs for Mitochondiral Glutamine tRNA
//TAGGATGGGGTGTGATAGGTGGCACGGAGAATTTTGGATTCTCAGGGATGGGTTCGATTCTCATAGTCCTAG
//(((((((...(((((....)..)))(.((((((..((((((((...)).)))))))))))))))).))))).