Sequence-matching-2

#include<bits/stdc++.h>

using namespace std;

void Alignment (char *, char *);
int Alignment_Q1 (char *A, char *B);
void Backward_Alignment (char *, char *);
void Forward_Alignment (char *, char *);
void Print(vector<char>);
int min(int ,int, int);

// Global Variables are CAPITALISED by convention
const int SIZE = 1000;
int F[SIZE], G[SIZE];
int GAP_PENALTY=2, MISMATCH_PENALTY=3;
int MATCH_PENALTY = 0;
int SCORE = 0;
vector < char > A_FINAL, B_FINAL;

int main(int argc, char **argv) {
   char a[SIZE], b[SIZE];

   cout << "s1: ";
   cin >> a;
   cout << "s2: ";
   cin >> b;


   Alignment(a,b);
   cout << "Gap Score: " << SCORE << "\n";
   Print(A_FINAL);
   Print(B_FINAL);

   return 0;
}

int Alignment_Q1 (char *A, char *B) {
   // getting length of the strings
   int m;
   for (m = 0; A[m] != '\0'; m++);

   int n;
   for (n = 0; B[n] != '\0'; n++);


   int arr[m+1][n+1];

   // Initializing Alignment Matrix
   for (int i = 0; i <= m; i++)
       arr[i][0] = i * GAP_PENALTY;
   for (int i = 0; i <= n; i++)
       arr[0][i] = i * GAP_PENALTY;

   // Store 0, mismatch_penalty depending on A[i] == B[j]
   int temp;


   for (int i = 1; i <= m; i++) {
       for (int j = 1; j <= n; j++) {
           // temp = 0 if Match-Found
           if (A[i-1] == B[j-1])
               temp = MATCH_PENALTY;

           // temp = PENALTY if Match Not Found
           else
               temp = MISMATCH_PENALTY;

           // Recurrence Relation
           arr[i][j] = min(
               temp + arr[i-1][j-1],
               GAP_PENALTY + arr[i][j-1],
               GAP_PENALTY + arr[i-1][j] );
       }
   }

   // Obtaining the Matched Sequences With Gaps
   char A_new[2*m + 1], B_new[2*n + 1];
   int i = m, j = n, a = 0, b = 0;
   while (i > 0 && j > 0)
   {
       // temp = 0 if Match-Found
       if (A[i-1] == B[j-1])
           temp = MATCH_PENALTY;

       // temp = PENALTY if Match Not Found
       else
           temp = MISMATCH_PENALTY;

       // Go diagonal => Match
       if (arr[i-1][j-1] + temp <= arr[i-1][j] + GAP_PENALTY && arr[i-1][j-1] + temp <= arr[i][j-1] + GAP_PENALTY)
       {
           A_new[a++] = A[--i];
           B_new[b++] = B[--j];
       }

       // Go Right => Anew.append(-), Bnew.append(B[j])
       else if (arr[i][j-1] + GAP_PENALTY <= arr[i-1][j-1] + temp && arr[i][j-1] +GAP_PENALTY <= arr[i-1][j] + GAP_PENALTY)
       {
           A_new[a++] = '-';
           B_new[b++] = B[--j];
       }

       // Go Down => Anew.append(A[i]), Bnew.append(-)
       else if (arr[i-1][j] + GAP_PENALTY <= arr[i-1][j-1] + temp && arr[i-1][j] + GAP_PENALTY <= arr[i][j-1] + GAP_PENALTY)
       {
           A_new[a++] = A[--i];
           B_new[b++] = '-';
       }
   }
   // Finishing Sequences: Edge Case
   while (i > 0) {
       A_new[a++] = A[--i];
       B_new[b++] = '-';
   }
   // Finishing Sequences: Edge Case
   while (j > 0) {
       A_new[a++] = '-';
       B_new[b++] = B[--j];
   }

   A_new[a++] = '\0';
   B_new[b++] = '\0';

   // Appending A i reverse order
   for (i = a - 2; i >= 0; i--)
       A_FINAL.push_back(A_new[i]);

   // Appending B in reverse order
   for (i = b - 2; i >= 0; i--)
       B_FINAL.push_back(B_new[i]);

   return arr[m][n];
}

void Forward_Alignment (char *A, char *B) {
   // Length of Sequence A
   int m;
   for (m = 0; A[m] != '\0'; m++);

   // Length of Sequnce B
   int n;
   for (n = 0; B[n] != '\0'; n++);

   // Alignment Matrix
   int arr[m+1][2];

   // Initializing Alignment Matrix
   for (int i = 0; i <= m; i++)
       arr[i][0] = i * GAP_PENALTY;

   // Store 0, mismatch_penalty depending on A[i] == B[j]
   int temp;

   // Evaluating Alignment Matrix Using Recurrence Relation
   // As per Keinberg and Tardos relation 6.16
   for (int j = 1; j <= n; j++) {
       arr[0][1] = j * GAP_PENALTY;
       for (int i = 1; i <= m; i++) {
           // temp = 0 if Match-Found
           if (A[i-1] == B[j-1])
               temp = 0;

           // temp = PENALTY if Match Not Found
           else
               temp = MISMATCH_PENALTY;

           // Recurrence Relation
           arr[i][1] = min(
               temp + arr[i-1][0],
               GAP_PENALTY + arr[i-1][1],
               GAP_PENALTY + arr[i][0] );
       }
       for (int i = 0; i <= m; i++)
           arr[i][0] = arr[i][1];
   }

   // returning f(:len_B)
   for ( int i = 0; i <= m; i++)
       F[i] = arr[i][1];
}

void Backward_Alignment (char *A, char *B) {
   // Length of Sequence A
   int m;
   for (m = 0; A[m] != '\0'; m++);

   // Length of Sequnce B
   int n;
   for (n = 0; B[n] != '\0'; n++);

   // Alignment Matrix
   int arr[m+1][2];

   // Initializing Alignment Matrix
   for (int i = 0; i <= m; i++)
       arr[i][1] = (m-i) * GAP_PENALTY;

   // Store 0, mismatch_penalty depending on A[i] == B[j]
   int temp;

   // Evaluating Alignment Matrix Using Recurrence Relation
   // As per Keinberg and Tardos relation 6.16
   for (int j = n-1; j >= 0; j--) {
       arr[m][0] = (n-j) * GAP_PENALTY;
       for (int i = m-1; i >= 0; i--) {
           // temp = 0 if Match-Found
           if (A[i] == B[j])
               temp = 0;

           // temp = PENALTY if Match Not Found
           else
               temp = MISMATCH_PENALTY;

           // Recurrence Relation
           arr[i][0] = min(
               temp + arr[i+1][1],
               GAP_PENALTY + arr[i+1][0],
               GAP_PENALTY + arr[i][1] );
       }
       for (int i = 0; i <= m; i++)
           arr[i][1] = arr[i][0];
   }

   // returning g(:len_B)
   for ( int i = 0; i <= m; i++)
       G[i] = arr[i][0];
}

// Sequencing in Linear Space-time
void Alignment(char *A, char *B) {
   // Length of Sequence A
   int m;
   for (m = 0; A[m] != '\0'; m++);

   // Length of Sequnce B
   int n;
   for (n = 0; B[n] != '\0'; n++);

   // Base Case
   if (m <= 2 || n <= 2) {
       SCORE  += Alignment_Q1(A, B);
       return ;
   }

   // Splitting B
   int i = 0;
   char B_first[n/2 + 1], B_last[n - n/2];
   for (i = 0; i <= n/2; i++)
       B_first[i] = B[i];
   B_first[i] = '\0';
   for (i = n/2 ; i < n; i++)
       B_last[i - n/2] = B[i];
   B_last[i - n/2] = '\0';

   // Dividing the problem
   Forward_Alignment(A, B_first);
   Backward_Alignment(A, B_last);

   // Conquering the problem
   // Minimising f(q, n/2) + g(q, n/2)
   int temp, index = 0, Min = F[0] + G[0];
   for (i = 1; i <= m; i++) {
       temp = F[i] + G[i];
       if (temp < Min) {
           Min = temp;
           index = i;
       }
   }

   // Splitting A
   char A_first[index + 1], A_last[m - index + 1];
   for (i = 0; i < index; i++)
       A_first[i] = A[i];
   A_first[i] = '\0';

   for (i = index; i < m; i++)
       A_last[i - index] = A[i];
   A_last[i - index] = '\0';

   // Splitting B
   for (i = 0; i < n/2; i++)
       B_first[i] = B[i];
   B_first[i] = '\0';

   for (i = n/2 ; i < n; i++)
       B_last[i - n/2] = B[i];
   B_last[i - n/2] = '\0';

   // Dividing the Problem Recursively
   Alignment(A_first, B_first);
   Alignment(A_last, B_last);
}


int min(int a, int b, int c) {
   if (a <= b && a <= c)
       return a;
   else if (b <= a && b <= c)
       return b;
   else
       return c;
}

void Print(vector<char> v) {
   for (auto it = v.begin(); it < v.end(); it++)
       cout << (*it) << "";
   cout << "\n";
}