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

char *S1,*S2;	// исходные строки алфавита;
int *S3,*S4;


int main()
{
int match = 2;  // w(S1[i],S2[j])= match, если S1[i] = S2[j]
int mismatch = -1;  // w(S1[i],S2[j])= mismatch, если S1[i] != S2[j]
int difference, similarity;
int space = 0, dif = 0, sim = 0; // количество '_' в выравненных строках,
								 // dif - количество несовпадающих символов 
								 // sim - количество совпадающих символов в 
								 // выравненных строках;
int penalty ; // штраф за разрыв
int w,wp1,wp2,Sum = 0; // Sum - оценка выравнивания(score)
int i,j,b,c;
int max,maxI,maxJ,maximum; // max - максимальное значение
						   // maxI - значение индекса i этого элемента в массиве Н
						   // maxJ - значение индекса j этого элемента в массиве Н
						   // maximum - максимальное значение матрицы Н

int curI,curJ;  // текущие индексы i и j

int **H;  // двумерный массив (или матрица)H используется
		  // для нахождения выравнивания с наивысшей оценкой. матрица
		  // содержит столько же строк, сколько символов в последова-
		  // тельности S1, и столько же столбцов, сколько символов в 
		  // последовательности S2. 
int **Prev; // если Prev[i][j]=0, то  Н[i][j] = 0;
			// если Prev[i][j]=1, то значение Н[i][j] было получено из Н[i][j-1] 
			// согласно правилам (1), описанным ниже;
			// если Prev[i][j]=2, то значение Н[i][j] было получено из Н[i-1][j-1];
			// если Prev[i][j]=3, то значение Н[i][j] было получено из Н[i-1][j];

int previous;
int **BLOSUM;	// матрица сравнений нуклеотидов,
				// в которой для каждой пары нуклеотидов 
                // есть свое значение их сравнения
int *MaxLoc,*maxL,*St;
int N,M;	// длины строк S1 и S2;
int endI,endJ; //  указывают на позиции в исходных строках, на которых 
			   //  было закончено выравнивание;
int beginI,beginJ;// указывают позиции в строках, начиная с которых 
				  // найдено оптимальное выравнивание;
FILE *fileS1,*fileS2, *file2;// файл fileS1 содержит строку S1 длины M,
							 // файл fileS2 содержит строку S2 длины N, 
							 // файл file2 хранит матрицу сравнений 
							 // нуклеотидов BLOSUM;
FILE *file3, *file4;
int *String1;   // стэк для считывания строки S1 из файла;
int *String2;	// стэк для считывания строки S2 из файла;

int *Str,*Str1;	// элементы стэков String1 и String2;
float time;

time = clock();

S3 = (int*)malloc(2*sizeof(int)); // строка S1 после выравнивания;
S4 = (int*)malloc(2*sizeof(int)); // строка S2 после выравнивания;

// открываем файл string1.txt для считывания строки S1,
// при этом игнорируем пробелы в этой строке и разрыв 
// строки S1 на две строки файла,определяем M-длину строки S1;
fileS1=fopen("string1.txt","r");
String1 = NULL;
M = 0;
while(!feof(fileS1)){
    b = (char)fgetc(fileS1);
    if((b != ' ')&&( b != '\n')){
    Str = (int*)malloc(2*sizeof(int));
    Str[0] = b;
    Str[1] = (int) String1;
    String1 = Str;
    M++;
    }
}

M--;
fclose(fileS1);

Str = String1;
String1 = (int*)String1[1];
free(Str);
// выделяем память для строки S1 длины M;
S1 = (char*)malloc(M*sizeof(char));

// заполняем строку символами из стэка String1;
for( i = M - 1; i>=0; i--){
    S1[i]=(char)String1[0];
	Str = String1;
    String1 = (int*)String1[1];
	free(Str);
}

// открываем файл string2.txt для считывания строки S2,
// при этом игнорируем пробелы в этой строке и разрыв 
// строки S2 на две строки файла,определяем N-длину строки S2;
fileS2=fopen("string2.txt","r");
N = 0;
String2 = NULL;
while(!(feof(fileS2))){
    b = (char)fgetc(fileS2);
    if((b != ' ')&&( b != '\n')){
    Str = (int*)malloc(2*sizeof(int));
    Str[0] = b;
    Str[1] = (int) String2;
    String2= Str;
    N++;
    }
}
fclose(fileS2);
N--;
S2= (char*)malloc(N*sizeof(char));

Str1 = String2;
String2 = (int*)String2[1];
free(Str1);
// заполняем строку S2;
for( i = N - 1; i>=0; i--){
    S2[i]=(char)String2[0];
	Str1 = String2;
    String2 = (int*)String2[1];
	free(Str1);
}

printf("\n Length of First string is  %i\n",M);
printf("\n Length of Second string is  %i\n",N);
// выделяем память для двумерной матрицы BLOSUM;
BLOSUM = (int**)malloc( 5 * sizeof(int*));
	for ( i = 0; i < 5;i ++)
		BLOSUM[i] = (int*)malloc( 5 * sizeof(int));

// считываем матрицу BLOSUM из файла BLOSUM.txt; поставь пробел в матрице после всех чисел
file2 = fopen("BLOSUM.txt","r");
fscanf(file2," %i",&penalty);

	for(i = 0; i < 5; i++)
		for(j = 0; j < 5; j++){
			fscanf(file2," %i",&b);
			while((b!=' ')&&(!feof(file2))){
				BLOSUM[i][j] = b;
				b=fgetc(file2);
			}

		}

fclose(file2);
/*
// выводим матрицу сравнений BLOSUM на экран;
printf("BLOSUM\n");
for(i=0;i<5;i++){
for(j=0;j<5;j++)
printf("  %i  ",BLOSUM[i][j]);
printf("\n");
}*/
// выделяем память для матрицы Н;
    H = (int**)malloc((M+1)*sizeof(int*));
	for ( i = 0; i < M+1;i ++)
		H[i] = (int*)malloc((N+1)*sizeof(int));

	// выделяем память для матрицы Н;
    Prev = (int**)malloc((M+1)*sizeof(int*));
	for ( i = 0; i < M+1;i ++)
		Prev[i] = (int*)malloc((N+1)*sizeof(int));


// заполняем первую строку и первый столбец матрицы Н нулями;
	for(i=0;i<=M;i++){
	H[i][0]=0;
	Prev[i][0] = 0;
	}

	for(j=1;j<=N;j++){
	H[0][j]=0;
	Prev[0][j] = 0;
	}

maximum = 0;
MaxLoc = (int)malloc(4*sizeof(int));
MaxLoc = NULL;

// строим матрицу Н по следующему принципу:
// находим максимум среди элементов 

//			0,
//	Н[i][j-1] + w(S2[j],_ ), 
//	H[i - 1][j - 1] + w( S1[i], S2[j] ),   (1)
//  H[i-1][j] + w( S1[i],_ ), 

// где w(S1[i],S2[j]) = match, если S1[i] = S2[j], 
// и w(S1[i],S2[j]) = mismatch в противном случае,
// w(S1[i],_ ) = w(S2[j],_ )= - 1;

// заносим максимумы матрицы Н в стэк MaxLoc, при этом 
// в стэке хранится не только максимальное значение, но 
// и указатель на этот элемент в другом стэке ft, в который 
// записываем для каждого элемента матрицы Н индексы i и j и
// индексы элементов из которых мы получили этот элемент
// (согдасно построению (1))
for ( i = 1; i <= M;i++ )
	for ( j = 1; j <= N; j++ ){
		max = 0;
		maxI = i - 1;
		maxJ = j - 1;
		previous = 0;
		if( S1[i-1]==S2[j-1])
			w = match;
		else
			w = mismatch;

		if (S1[i-1] != '_')
			wp1 = mismatch;
		else wp1 = match;
		
		if (S2[j-1] != '_')
			wp2 = mismatch;
		else wp2 = match;

        if ( ( H[i][ j-1] + wp2 )> max ){
			max = H[i][ j-1] + wp2 ;
			maxI = i;
			maxJ = j-1;
			previous = 1;
		}
		if (( H[i-1][ j-1] + w )> max ) {
			max = H[i-1][ j-1] + w ;
			maxI = i-1;
			maxJ = j-1;
			previous = 2;
		}
		if ( ( H[i-1][ j] + wp1 )> max ) {
			max = H[i-1][ j] + wp1 ;
			maxI = i-1;
			maxJ = j;
			previous = 3;
		}

		H[i][j] = max;
		Prev[i][j] = previous;
		
		if ( H[i][j] >= maximum ){
			maximum = H[i][j];
			maxL = (int*)malloc(4*sizeof(int));
			maxL[0] = maximum;
			maxL[1] = i;
			maxL[2] = j;
			maxL[3] = (int)MaxLoc;
			MaxLoc = maxL;
		}
}

maxL = MaxLoc;
St = (int*)maxL[3];

while ( St ){
	if(St[0] != maximum ){
		maxL[3] = St[3];
		free(St);
		St = (int*)maxL[3];
	}
	else 
		maxL =  (int*)maxL[3];	
		St = (int*)maxL[3];
}

for ( j=0; j<M+1;j++)
		free (H[j]);
free (H);

file3 = fopen("newS1.txt","w");
file4 = fopen("newS2.txt","w");


// находим все возможные оптимальные выравнивания строк S1 и S2.
// выравнивание выполняется следующим образом:
//	-	после того как нашли максимальное значение в матрице Н
//      находим обратный путь от этого максимума в позицию (0,0)
//		или до первого 0 в матрице Н;
//  -   анализируем построенный путь:
//			*   если значение Н[i][j] получено из H[i-1][j-1], то 
//          записываем символы S1[i] и S2[j] в новые строки, при
//			этом символы могут либо совпадать, либо отличаться,
//			что влияет на оценку выравнивания (оценка меняется 
//			за счет значений из матрицы BLOSUM )
//			*	если значение Н[i][j] получено из H[i][j-1], то 
//			ставим в первой строке '_',а S2[j]записываем в новую 
//			строку S4
//			*	если значение Н[i][j] получено из H[i-1][j], то 
//      	ставим во второй строке '_',а S1[i] записываем в новую 
//			строку S3

while(MaxLoc){ 
Sum = 0;
space = 0;
dif = 0; 
sim = 0;

if(MaxLoc){ 
	S3 = NULL;
	S4 = NULL;
	i = MaxLoc[1];
	j = MaxLoc[2];
	endI = i;
	endJ = j;
	
	while( Prev[i][j]!=0 ){	
			beginI = i;
			beginJ = j;
			if( Prev[i][j] == 3 ){ 
                Str = (int*)malloc(2*sizeof(int));
				Str[0] = (char)'_';
				Str[1] = (int) S4;
				S4=Str;
				space++;
				Sum = Sum + penalty;
                Str1 = (int*)malloc(2*sizeof(int));
				Str1[0] = (char)S1[i-1];
				Str1[1] = (int) S3;
				S3=Str1;
				i = i - 1;
			}
			 if( Prev[i][j] == 2 ){
				Str= (int*)malloc(2*sizeof(int));
				Str[0] = (char)S2[j-1];
				Str[1] = (int)S4;
				S4=Str;
				Str1 = (int*)malloc(2*sizeof(int));

				Str1[0] = (char)S1[i-1];
				Str1[1] = (int) S3;
				S3=Str1;
				i = i - 1;
				j = j - 1;
				if (S3[0]== S4[0]){
					if(( S3[0] == 'A')||( S3[0] == 'a')) similarity = BLOSUM [0][0];
					if(( S3[0] == 'G')||( S3[0] == 'g')) similarity = BLOSUM [1][1];
					if(( S3[0] == 'C') ||( S3[0] == 'c'))similarity = BLOSUM [2][2];
					if(( S3[0] == 'T')||( S3[0] == 't')) similarity = BLOSUM [3][3];
					if( S3[0] == '_' ) similarity = BLOSUM [4][4];
					Sum = Sum + similarity;
					sim++;
				}
				else {
                    if((S3[0] == 'N')||(S4[0] == 'N')||(S3[0] == 'n')||(S4[0] == 'n'))
                        difference = penalty/2;
					if(( S3[0] == 'A')||( S3[0] == 'a')) {
						if ((S4[0] == 'G')||( S4[0] == 'g'))
							difference = BLOSUM [1][0];
						if ((S4[0] == 'C')||( S4[0] == 'c'))
							difference = BLOSUM [2][0];
						if (( S4[0] == 'T')||( S4[0] == 't'))
							difference = BLOSUM [3][0];
						if ( S4[0] == '_')
							difference = BLOSUM [4][0];
					}
					if(( S3[0] == 'G')||( S3[0] == 'g')){
						if ((S4[0] == 'A')||( S4[0] == 'a'))
							difference = BLOSUM [1][0];
						if ((S4[0] == 'C')||( S4[0] == 'c'))
							difference = BLOSUM [2][1];
						if ((S4[0] == 'T')||( S4[0] == 't'))
							difference = BLOSUM [3][1];
						if ( S4[0] == '_')
							difference = BLOSUM [4][1];
					}
					if(( S3[0] == 'C') ||( S3[0] == 'c')){
						if (( S4[0] == 'A')||( S4[0] == 'a'))
							difference = BLOSUM [2][0];
						if ((S4[0] == 'G')||( S4[0] == 'g'))
							difference = BLOSUM [2][1];
						if (( S4[0] == 'T')||( S4[0] == 't'))
							difference = BLOSUM [3][2];
						if ( S4[0] == '_')
							difference = BLOSUM [4][2];
					}
					if((S3[0] == 'T')||( S3[0] == 't')){
						if ((S4[0] == 'A')||( S4[0] == 'a'))
							difference = BLOSUM [3][0];
						if ((S4[0] == 'G')||( S4[0] == 'g'))
							difference = BLOSUM [3][1];
						if ((S4[0] == 'C')||( S4[0] == 'c'))
							difference = BLOSUM [3][2];
						if ( S4[0] == '_')
							difference = BLOSUM [4][3];
					}
					if( S3[0] == '_') {
						if ((S4[0] == 'G')||( S4[0] == 'g'))
							difference = BLOSUM [4][1];
						if ((S4[0] == 'C')||( S4[0] == 'c'))
							difference = BLOSUM [4][2];
						if (( S4[0] == 'T')||( S4[0] == 't'))
							difference = BLOSUM [4][3];
						if (( S4[0] == 'A')||( S4[0] == 'a'))
							difference = BLOSUM [4][0];
					}
					Sum = Sum + difference;
					dif++;
				}
			}

		else { 
			Sum = Sum + penalty;
            Str1 = (int*)malloc(2*sizeof(int));
            Str1[0] = (char)'_';
            Str1[1] = (int) S3;
            S3=Str1;

			Str = (int*)malloc(2*sizeof(int));
			Str[0] = (char)S2[j-1];
			Str[1] = (int)S4;
			S4=Str;
			space++;
			j = j - 1;
		}
}

Str = (int*)malloc(2*sizeof(int));
Str1 = (int*)malloc(2*sizeof(int));

// записываем полученную выравненную строку S1 в файл
// newS1.txt, каждую новую выравненную строку записываем 
// с новой строки в файле и ставим перед ней '*'
printf("\n New S1 :\n");
fprintf(file3," %c",'*');
while(S3){
	printf(" %c",S3[0]);
	fprintf(file3,"%c",S3[0]);
    Str=S3;
    S3=(int*)Str[1];
    free(Str);
}
fprintf(file3,"	  { %i  %i }",beginI,endI);
fprintf(file3,"%c",'\n');

// записываем полученную выравненную строку S2 в файл
// newS2.txt, каждую новую выравненную строку записываем 
// с новой строки в файле и ставим перед ней '*'
fprintf(file4," %c",'*'); 
printf("\n New S2 :\n");
	while(S4){
        printf(" %c",S4[0]);
		fprintf(file4,"%c",S4[0]);
        Str1=S4;
        S4=(int*)(Str1[1]);
        free(Str1);
    }
fprintf(file4,"	  { %i  %i }",beginJ,endJ);
fprintf(file4,"%c",'\n');


printf("\n Score = %i ",Sum);
printf(" \n Simularities = %i", sim);
printf("\n Differences = %i ",dif);
printf(" \n Spaces = %i \n", space);
printf(" Position of alignment at First string %i %i",beginI,endI);
printf(" \n Position of alignment at Second string %i %i \n",beginJ,endJ);
MaxLoc = (int*)MaxLoc[3];
	}
}

fclose(file3);
fclose(file4);

// очищаем память использованную для массивов BLOSUM,S1,S2
// и стэков S3,S4,MaxLoc,maxL

for ( j=0; j<M+1;j++)
		free (Prev[j]);
free (Prev);

for ( i=0; i< 4;i++)
		free (BLOSUM[i]);
free (BLOSUM);

free(S4);
free(S3);
free(S1);
free(S2);
free(MaxLoc);
free(maxL);


system("pause");

    return 0;
}