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

initialize.c
P. Clote

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

#ifndef __SPECTRUM_H__
#include "spectrum.h"
#endif

void initializeExistBP(MatrixType M, TripletMatrixType S, short *a, short ExistBP[NN][NN]){
  //RNA sequence a1,...,an with n length of sequence
  int i,j,k,k0,k1,d,max,num,r,n,x,y,index;
  //Check if there exist base pairs in region i...j

  n = a[0]; //length of RNA sequence a[1],...,a[n]

  //First, set ALL values
  for (i=0;i<=n;i++)
    for (j=0;j<=n;j++){
      ExistBP[i][j]=0;
      for (k=0;k<KK;k++){
        BP(i,j,k)  = 0;
        SS(i,j,k)[0]  = 0;
        SS(i,j,k)[1]  = 0;
        SS(i,j,k)[2]  = 0;
        }//for k
     }//for j

  for (d=THETA+1;d<=n-1;d++){
    for (i=1;i<=n-THETA-1;i++){
      j = i+d;
      if (j>n) break;
      if (basePair(a[i],a[j])||ExistBP[i][j-1]||ExistBP[i+1][j]) 
        ExistBP[i][j]=1;
    }//endfor i
  }//end{for d
}





void initialize(MatrixType M, TripletMatrixType S, short *a,short ExistBP[NN][NN]){
  //For this to work, need to previously call initializeExistBP
  int i,j,k,k0,k1,d,max,num,r,n,x,y,index;
  int K0=0,K1=0;
  //Check if there exist base pairs in region i...j

  n = a[0]; //length of RNA sequence a[1],...,a[n]
  for (d=THETA+1;d<=n-1;d++)
    for (i=1;i<=n-THETA-1;i++){
      j = i+d;
      if (j>n) break;
      for (k=1;k<=FF(i,j);k++){       
        //CASE 1
        max   = BP(i,j-1,k);
        index = 0;
        //CASE 2
        if ( basePair(a[i],a[j]) ){
          num = BP(i+1,j-1,k);
          if (k==1&&num==0 || num>0){
            num += 1; //extra 1 for basepair (i,j)
            if (num>max){
              max   = num;
              index = i;
              }
            }
          }
        //CASE 3
        for (r=i+1;r<=j-THETA-1;r++){
          if ( basePair(a[r],a[j]) ){
            //Case  3: k = k0+k1 and r,j basepair
            //Case 3a: k0=0, k1=k and no bp in region i..r-1
//            if (!ExistBP[i][r-1]){
//+++
              num = BP(r+1,j-1,k); 
              if (k==1&&num==0 || num>0){
                num += 1; //add extra bp for (r,j)
                if (max<num){
                  max   = num;
                  index = r;
                  K0    = 0;
                  K1    = k;
                  }
              } 
//            }//endif
            for (k0=1;k0<k;k0++){ 
              k1=k-k0;
              x = BP(i,r-1,k0); 
              y = BP(r+1,j-1,k1); 
              if (k1==1&&x>0&&y==0 || x>0&&y>0)
                num = x+y+1; //extra 1 for basepair (r,j)
                if (max<num){
                   max=num;
                   index = r;
                   K0    = k0;
                   K1    = k1;
                   }
              } //end for k0
            }//endif
          }//endfor r
        BP(i,j,k) = max;
        SS(i,j,k)[0] = index;
        SS(i,j,k)[1] = K0;
        SS(i,j,k)[2] = K1;
        }//endfor k
    }//endfor i
}


int backtrack(int i, int j, int k, MatrixType M, TripletMatrixType S,
   short ExistBP[NN][NN], char *rna, char paren[NN]){
   int k0,k1,r;
   
   if (k>0 && j-i>THETA){ 
     r  = S[i][j][k][0];
     k0 = S[i][j][k][1];
     k1 = S[i][j][k][2];
     if (r==0) //Case 1: j unpaired in [i,j]
       backtrack(i,j-1,k,M,S,ExistBP,rna,paren);
     else{// r>0 so (r,j) is base pair in structure with k hairpins
       if (DEBUG) printf("Base pair r=%d,j=%d %c %c\n",r,j,rna[r-1],rna[j-1]);
       //WARNING rna is 0-indexed, a[NN] and paren[NN] are 1-indexed
       paren[r] = '('; 
       paren[j] = ')';
       if (r==i) { //Case 2: r=i, no left side
         if (M[i+1][j-1][k]>0)
           backtrack(r+1,j-1,k,M,S,ExistBP,rna,paren);
         }
       else {
         backtrack(i,r-1,k0,M,S,ExistBP,rna,paren);
         backtrack(r,j,k1,M,S,ExistBP,rna,paren);
         } 
     } 
  }  
}// endBacktrack