/*
 * Author : DaeGon Kim
 * Date   : Oct 29, 2007
 * Desc   : A sequential implementation of Needleman-Wunsch algorithm
 */

#include <stdlib.h>
#include <stdio.h>
#include <assert.h>
#include "timer.h"

#define SMatrix(i,j)    sMatrix[(i)*alphabets+(j)]

#define MAX(a,b)        ((a)>(b) ? (a) : (b))
#define MIN(a,b)        ((a)>(b) ? (b) : (a))

#define MAX3(a,b,c)     ((c)>MAX(a,b) ? (c) : MAX(a,b))

int main(int argc, char **argv) {

   int alphabets, N, M, b;
   int penalty = -3;
   
   int *sMatrix;
   int *seqA, *seqB;

   int *prev, *cur, *temp;

   int size, i, j;

   // Time
   double time;

   alphabets = 4;
   N         = 100;
   M         = 100;
   b         = 100;
   
   if ( argc > 1 ) alphabets = atoi(argv[1]);
   if ( argc > 2 ) N = atoi(argv[2]);
   if ( argc > 3 ) M = atoi(argv[3]);
   if ( argc > 4 ) b = atoi(argv[4]);
   
   printf("alphabets = %d, N = %d, M = %d, b = %d\n", alphabets, N, M, b);

   // Similarity Matrix - Identity
   size    = alphabets * alphabets * sizeof(int);
   sMatrix = (int *)malloc(size);
   assert(sMatrix != NULL);
   for ( i=0 ; i < alphabets ; i++ ) {
      for (j=0 ; j < alphabets ; j++ ) {
         SMatrix(i,j) = 0;
      }
      SMatrix(i,i) = 1;
   }
   // END of Similarity Matrix - Identity

   // Input Sequence generation
   size    = N * sizeof(int);
   seqA    = (int *)malloc(size);
   assert(seqA != NULL);

   size    = M * sizeof(int);
   seqB    = (int *)malloc(size);
   assert(seqB != NULL);

   srand(7);
   for ( i=0 ; i < N ; i++ ) {
      seqA[i] = rand() % alphabets;
   }

   srand(13);
   for ( i=0 ; i < M ; i++ ) {
      seqB[i] = rand() % alphabets;
   }
   #ifdef DEBUG
      printf("SeqA of size %d\n", N);
      for ( i=0 ; i < N ; i++ ) {
         printf("SeqA[%d] = %d\n", i, seqA[i]);
      }
      printf("\n\n\nSeqB of size %d\n", M);
      for ( i=0 ; i < M ; i++ ) {
         printf("SeqB[%d] = %d\n", i, seqB[i]);
      }
      printf("\n\n\n");
   #endif
   // END of Input Sequence generation

   size    = (N) * sizeof(int);
   prev    = (int *)malloc(size);
   cur     = (int *)malloc(size);
   assert(prev != NULL && cur  != NULL);
   
   initialize_timer ();
   start_timer();

   for ( i=1 ; i <= M ; i++ ) {
      for ( j=1 ; j <= N ; j++ ) {
         if ( i==1 && j==1 ) {
            cur[j-1] = MAX3(penalty,          penalty,           SMatrix(seqB[i],seqA[j]));
         } else if ( i==1 && j>1) {
            cur[j-1] = MAX3(cur[j-2]+penalty, penalty*j        , penalty*(j-1) + SMatrix(seqB[i-1],seqA[j-1]));
         } else if ( i>1  && j == 1 ) {
            cur[j-1] = MAX3(penalty*i,        prev[j-1]+penalty, penalty*(i-1) + SMatrix(seqB[i-1],seqA[j-1]));
         } else {
            cur[j-1] = MAX3(cur[j-2]+penalty, prev[j-1]+penalty, prev[j-2]     + SMatrix(seqB[i-1],seqA[j-1]));
         }
      }
      #ifdef DEBUG
         printf("i= %d\n", i);
         for ( j=0 ; j < N ; j++ ) {
            printf("cur[%d] = %d\n", j, cur[j]);
         }
         printf("\n");
      #endif
      temp = prev; prev = cur; cur = temp;
   }

   stop_timer();
   time = elapsed_time ();

   printf("The optimal value is %d, Time taken : %lf\n", prev[N-1], time);

   return 0;
}