/*---------------------------
  rnaSecStr.c

  P. Clote
---------------------------*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <limits.h>
#include <string.h>
#include <time.h>

#define DEBUG 0
#define NMAX 1000
#define INF INT_MAX/2
#define GUbond -1
#define AUbond -1
#define CGbond -1
#define THRESHOLD 3
        // for i and j to basepair, |i-j| >= THRESHOLD

/*---------prototypes ----------*/
void error(char *);
int guPair(char, char);
int auPair(char, char);
int cgPair(char, char);
int a(int, int, char *);
int basePair(char, char);
int watsonCrick(char, char);
void initialEnergyMatrix(int[NMAX][NMAX], int);
void printMatrix(int[NMAX][NMAX],int);
void computeEnergyMatrix(int[NMAX][NMAX],char[NMAX],int);
int backtrack(int,int,int[NMAX][NMAX],char[NMAX],char[NMAX]);


void error(char *s){
   fprintf(stderr,"%s\n",s);
   exit(1);
}

int guPair(char x, char y){
   return ( x=='G'&&y=='U' || x=='U'&&y=='G' );
}

int auPair(char x, char y){
   return ( x=='A'&&y=='U' || x=='U'&&y=='A' );
}

int cgPair(char x, char y){
   return ( x=='C'&&y=='G' || x=='G'&&y=='C' );
}



int a(int i, int j, char *s){
   if (i<strlen(s) && j<strlen(s)) {
      if (auPair(s[i],s[j]))
         return AUbond;
      else if (cgPair(s[i],s[j]))
         return CGbond;
      else if (guPair(s[i],s[j]))
         return GUbond;
      else
         return INF;
      }
    else
         error("Indices out of bounds for RNA sequence\n");
}




int basePair(char x, char y){
   return( auPair(x,y) || cgPair(x,y) || guPair(x,y) );
}


int watsonCrick(char x, char y){
   return( auPair(x,y) || cgPair(x,y) );
}


void initialEnergyMatrix(int energy[NMAX][NMAX], int n){
   int i,j;
   for (i=0;i<n;i++) 
      for (j=i;j<n;j++){
         energy[j][i]=-1;
         energy[i][j]=0;
         }

}

void printMatrix(int a[NMAX][NMAX], int n){
   int i,j;
   for (i=0;i<n;i++) {
      for (j=0;j<n;j++) 
         printf("%d  ",a[i][j]);
      printf("\n");
   }
}




main(int argc, char *argv[]){
   char S[NMAX],paren[NMAX]; 
   int i;
       /*-------------------------------
        S is RNA sequence
        paren is parenthesis expression for 
        optimal RNA secondary structure
       -------------------------------*/

   int E[NMAX][NMAX];
   int n;
   
   /*----- error checking for input --*/ 
   if(argc < 2){ 
      printf("%s RNAnucleotideString\n",argv[0]);
      exit(1);
   }
    
   
   /*----- initialization  --*/ 
   strcpy(S,argv[1]);
   n = strlen(argv[1]);
   initialEnergyMatrix(E,n);
   for (i=0;i<n;i++) paren[i]='.';
   for (i=n;i<NMAX;i++) paren[i]=' ';
   paren[n]='\0'; // paren is string of same length as S
   
   computeEnergyMatrix(E,S,n);
   backtrack(0,n-1,E,S,paren);
   printf("%s\n",S);
   printf("%s\n",paren);
   if (DEBUG){
     printf("Number of nucleotides = %d\n",n);
     printMatrix(E,n);
     }


}// endmain



void computeEnergyMatrix(int E[NMAX][NMAX],char S[NMAX],int n){   
   int i,j,d,k,min,val,index;

   for (d = THRESHOLD+1; d < n; d++){
     for(i=0; i< n; i++){
        // printf("d=%d\n",d); 
        j=i+d;
        if (j < n){ 
          min =0;
          index=n; 
           /*-------------------------------------
           if index<n at end of for-loop, then this
           means that index and j form a base pair,
           and this is noted by E[j][i]=index.
           if index=n at end of for-loop, then this
           means that j is not base paired.
           -------------------------------------*/

          if ( E[i][j-1]<min ){
             min = E[i][j-1];
             index = n;
             // j not basepaired with some k such that i<k<j
          }
          if ( a(i,j,S) + E[i+1][j-1] < min ){
             min = a(i,j,S) + E[i+1][j-1];
             index=i;
          }
          for(k=i+1; k<j-THRESHOLD; k++){ 
             val = a(k,j,S) + E[i][k-1] + E[k+1][j-1];
             if (val < min) {
                min = val;
                index=k;
             }
          }
          E[i][j]=min;
          if (index<n) 
             E[j][i]=index;
          else
             E[j][i]=-1;
       } //endif
     } // endfor
   } // endfor

}// endComputeEnergyMatrix



int backtrack(int i, int j, int E[NMAX][NMAX], 
              char S[NMAX], char paren[NMAX]){
   int k;
   if (j-i>THRESHOLD){ 
     k = E[j][i];
     if(k != -1)  {
       if (DEBUG) printf("Base pair k=%d,j=%d   %c %c\n",k,j, S[k], S[j] );
       paren[k] = '('; 
       paren[j] = ')';
       if( THRESHOLD <= (j-1)-(k+1))
         backtrack(k+1,j-1,E,S,paren);
       if (THRESHOLD <= k-1-i  )
         backtrack(i,k-1,E,S,paren);
     } 
     else{ // k==-1
       if (DEBUG) printf("k=-1\n");
       if( THRESHOLD <= j-1-i ){
          // printf("i=%d j=%d\n",i,j-1);         
          backtrack(i,j-1,E,S,paren);
          }
       else 
          return 0;
     }
  }  
   
}// endBacktrack