/* Pairwise Alignment - Linear Gap Penalty
 * 
 * Used to align a pair of aminoacid sequences given in the file 'nizovi.txt'.
 * This program implements Needleman-Wunsch algorithm and uses the BLOSUM50
 * matrix.
 *
 * Copyright (c) 2006 Filip Niksic
 * 
 * Permission is hereby granted, free of charge, to any person obtaining
 * a copy of this software and associated documentation files (the "Software"),
 * to deal in the Software without restriction, including without limitation
 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
 * and/or sell copies of the Software, and to permit persons to whom the
 * Software is furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included
 * in all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 *
 ***/

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

#define maxn 1000

int B[20][20] = {
    {5, -2, -1, -2, -1, -1, -1, 0, -2, -1, -2, -1, -1, -3, -1,  1,  0, -3, -2,  0},
    {-2, 7,-1, -2, -4,  1,  0, -3,  0, -4, -3, 3, -2, -3, -3, -1, -1, -3, -1, -3},
    {-1, -1, 7, 2, -2,  0,  0,  0,  1, -3, -4, 0, -2, -4, -2,  1,  0, -4, -2, -3},
    {-2, -2, 2, 8, -4,  0,  2, -1, -1, -4, -4, -1, -4, -5, -1,  0, -1, -5, -3, -4},
    {-1, -4, -2, -4, 13, -3, -3, -3, -3, -2, -2, -3, -2, -2, -4, -1, -1, -5, -3, -1},
    {-1, 1, 0, 0, -3, 7, 2, -2,  1, -3, -2,  2, 0, -4, -1,  0, -1, -1, -1, -3 },
    {-1, 0, 0, 2, -3, 2, 6, -3,  0, -4, -3,  1, -2, -3, -1, -1, -1, -3, -2, -3},
    {0, -3, 0, -1, -3, -2, -3, 8, -2, -4, -4, -2, -3, -4, -2,  0, -2, -3, -3, -4},
    {-2, 0, 1, -1, -3, 1, 0, -2, 10, -4, -3, 0, -1, -1, -2, -1, -2, -3,  2, -4},
    {-1, -4, -3, -4, -2, -3, -4, -4, -4, 5, 2, -3,  2, 0, -3, -3, -1, -3, -1,  4},
    {-2, -3, -4, -4, -2, -2, -3, -4, -3, 2, 5, -3,  3,  1, -4, -3, -1, -2, -1,  1},
    {-1, 3, 0, -1, -3, 2, 1, -2, 0, -3, -3, 6, -2, -4, -1,  0, -1, -3, -2, -3},
    {-1, -2, -2, -4, -2, 0, -2, -3, -1, 2, 3, -2, 7, 0, -3, -2, -1, -1,  0,  1},
    {-3, -3, -4, -5, -2, -4, -3, -4, -1, 0, 1, -4, 0, 8, -4, -3, -2,  1,  4, -1},
    {-1, -3, -2, -1, -4, -1, -1, -2, -2, -3, -4, -1, -3, -4, 10, -1, -1, -4, -3, -3},
    {1, -1, 1, 0, -1, 0, -1, 0, -1, -3, -3, 0, -2, -3, -1, 5, 2, -4, -2, -2},
    {0, -1, 0, -1, -1, -1, -1, -2, -2, -1, -1, -1, -1, -2, -1, 2, 5, -3, -2,  0},
    {-3, -3, -4, -5, -5, -1, -3, -3, -3, -3, -2, -3, -1, 1, -4, -4, -3, 15, 2, -3},
    {-2, -1, -2, -3, -3, -1, -2, -3, 2, -1, -1, -2, 0, 4, -3, -2, -2, 2, 8, -1},
    {0, -3, -3, -4, -1, -3, -3, -4, -4, 4, 1, -3, 1, -1, -3, -2, 0, -3, -1, 5}
};

// Inicijaliziramo matricu tipa mxn
int **matrica(int m,int n) {
    int **A,i;

    A=(int **)malloc(m*sizeof(int *));
    A[0]=(int *)calloc(m*n,sizeof(int));
    for (i=1;i<m;i++) A[i]=A[i-1]+n;

    return A;
}

// Oslobodimo alociranu matricu
void brisi_matricu(int **A) {
    free(A[0]);
    free(A);
}

// Pridruzivanje aminokiselina -> broj
int broj(char c){
    switch (c) {
	case 'A': return 0; break;
	case 'R': return 1; break;
	case 'N': return 2; break;
	case 'D': return 3; break;
	case 'B': return 3; break;
	case 'C': return 4; break;
	case 'Q': return 5; break;
	case 'E': return 6; break;
	case 'Z': return 6; break;
	case 'G': return 7; break;
	case 'H': return 8; break;
	case 'I': return 9; break;
	case 'L': return 10; break;
	case 'K': return 11; break;
	case 'M': return 12; break;
	case 'F': return 13; break;
	case 'P': return 14; break;
	case 'S': return 15; break;
	case 'T': return 16; break;
	case 'W': return 17; break;
	case 'Y': return 18; break;
	case 'V': return 19; break;
	default: return -1;
    }
}

int main(void) {
    FILE *fp;
    char buffer[maxn];
    char *x,*y,*xp,*yp,t;
    
    int d=8;
    int m,n,i,j,l,a,b,c;

    int **F,**P;

    /* Prva faza je ucitavanje nizova */
    
    fp=fopen("nizovi.txt","r");
    
    fgets(buffer,maxn,fp);
    // Prvi znak u nizu koji ne predstavlja aminokiselinu
    // zamijenimo s \0.
    for (i=0;broj(buffer[i])!=-1;i++);
    buffer[i]='\0';
    
    n=strlen(buffer);
    x=(char *)malloc((n+1)*sizeof(char));
    strcpy(x,buffer);

    fgets(buffer,maxn,fp);
    for (i=0;broj(buffer[i])!=-1;i++);
    buffer[i]='\0';
    
    m=strlen(buffer);
    y=(char *)malloc((m+1)*sizeof(char));
    strcpy(y,buffer);

    fclose(fp);

    /* Imamo nizove, sad inicijaliziramo matrice.
     * Pointeri u matrici P: 1 znaci <-, 2 znaci |, 4 znaci \ */
    F=matrica(m+1,n+1);
    P=matrica(m+1,n+1);

    for (i=0;i<=m;i++) {
	F[i][0]=-i*d;
	P[i][0]=2;
    }
    for (j=0;j<=n;j++) {
	F[0][j]=-j*d;
	P[0][j]=1;
    }

    /* Krecemo s popunjavanjem matrica */
    for (i=1;i<=m;i++)
	for (j=1;j<=n;j++) {
	    a=F[i-1][j-1]+B[broj(y[i-1])][broj(x[j-1])];
	    b=F[i-1][j]-d;
	    c=F[i][j-1]-d;
	    if (a>=b && a>=c) {
		P[i][j]=4;
		F[i][j]=a;
	    }
	    else if (b>=a && b>=c) {
		P[i][j]=2;
		F[i][j]=b;
	    }
	    else {
		P[i][j]=1;
		F[i][j]=c;
	    }
	}

    /* Sad treba napraviti poravnavanje */
    xp=(char *)malloc((m+n+1)*sizeof(char));
    yp=(char *)malloc((m+n+1)*sizeof(char));

    l=0;i=m;j=n;
    while (i>0 || j>0) {
	if (P[i][j]==1) {
	    xp[l]=x[j-1];
	    yp[l]='-';
	    j--;
	}
	else if (P[i][j]==2) {
	    xp[l]='-';
	    yp[l]=y[i-1];
	    i--;
	}
	else {
	    xp[l]=x[j-1];
	    yp[l]=y[i-1];
	    i--;j--;
	}
	l++;
    }
    xp[l]=yp[l]='\0';

    /* Dobili smo obrnute nizove, treba ih okrenuti */
    for (i=0;i<l/2;i++) {
	t=xp[i];xp[i]=xp[l-i-1];xp[l-i-1]=t;
	t=yp[i];yp[i]=yp[l-i-1];yp[l-i-1]=t;
    }

    /* Zapisemo rezultat, oslobodimo memoriju */
    fp=fopen("poravnavanje.txt","w");
    fprintf(fp,"%s\n%s\n",xp,yp);
    fclose(fp);

    brisi_matricu(F);
    brisi_matricu(P);
    free(x);
    free(y);
    free(xp);
    free(yp);

    return 0;
}