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

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


int Max(int x, int y)
{ if (x >= y) 
		return x;
  else 
	  return y;
}

int Min(int x, int y)
{ if (x < y) 
		return x;
  else 
	  return y;
}

int main()
{
int *MaxLoc,*maxL,*St;
int maximum;
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 difference, similarity;
int space, dif, sim; // количество '_' в выравненных строках,
								 // dif - количество несовпадающих символов 
								 // sim - количество совпадающих символов в 
								 // выравненных строках;
int penalty; // штраф за разрыв
int previous;
	int match;  // w(S1[i],S2[j])= match, если S1[i] = S2[j]
	int mismatch;  // w(S1[i],S2[j])= mismatch, если S1[i] != S2[j]
	int max,maxI,maxJ; // max - максимальное значение
						   // maxI - значение индекса i этого элемента в массиве Н
						   // maxJ - значение индекса j этого элемента в массиве Н
						   // maximum - максимальное значение матрицы Н
	int w,*wp1,*wp2;

int Sum ; // Sum - оценка выравнивания(score)
int i,j,b,k1,k2,p,d;

//int curI,curJ;  // текущие индексы i и j
int **H;  // двумерный массив (или матрица)H используется
		  // для нахождения выравнивания с наивысшей оценкой. матрица
		  // содержит столько же строк, сколько символов в последова-
		  // тельности S1, и столько же столбцов, сколько символов в 
		  // последовательности S2. 

int **BLOSUM;	// матрица сравнений нуклеотидов,
				// в которой для каждой пары нуклеотидов 
                // есть свое значение их сравнения
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;
int Size;

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

space = 0;
dif = 0;
sim = 0;
Sum = 0;
penalty = -15;
match = 2;
mismatch = -1;
// открываем файл 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");
	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;
	}

	for(i =0;i<=M;i++)
		for(j =0;j<=N;j++)
			H[i][j]=0;

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

//			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;
Size = M;
p = 0;
if (N < M){
	p = 1;
	Size = N;
}
d = M - Size;
wp1 = (int*)malloc(M*sizeof(int));
wp2 = (int*)malloc(N*sizeof(int));

for( i = 1; i <= M; i++)
	if (S1[i-1] != '_')
			wp1[i-1] = mismatch;
	else wp1[i-1] = match;

for( i = 1; i <= N; i++)
	if (S2[i-1] != '_')
			wp2[i-1] = mismatch;
	else wp2[i-1] = match;


// заносим максимумы матрицы Н в стэк MaxLoc, при этом 
// в стэке хранится не только максимальное значение, но 
// и указатель на этот элемент в другом стэке ft, в который 
// записываем для каждого элемента матрицы Н индексы i и j и
// индексы элементов из которых мы получили этот элемент
// (согдасно построению (1))

  
  if (N < M+1)
  {
      for (j = 2; j <= N; j = j + 1)
      {
	    #pragma omp parallel for
        for (i = 1; i <= j - 1; i = i + 1)
        {
              if (S1[i - 1] == S2[(j - 1) - i])
              {
                H[i][j - i] = H[i - 1][(j - 1) - i] + 2;
                Prev[i][j - i] = 2;
              }
              else
              {
                H[i][j - i] = H[i - 1][(j - 1) - i] - 1;
                Prev[i][j - i] = 2;
              }
              if (H[i][(j - 1) - i] + wp2[(j - i) - 1] > H[i][j - i])
              {
                H[i][j - i] = H[i][(j - 1) - i] + wp2[(j - i) - 1];
                Prev[i][j - i] = 1;
              }
              if (H[i][j - i] < 0)
              {
                Prev[i][j - i] = 0;
                H[i][j - i] = 0;
              }
              if (H[i - 1][j - i] + wp1[i - 1] > H[i][j - i])
              {
                H[i][j - i] = H[i - 1][j - i] + wp1[i - 1];
                Prev[i][j - i] = 3;
              }
        }
      }
  }
  else
  {
      for (j = 2; j < M+1; j = j + 1)
      {
        #pragma omp parallel for
  	    for (i = 1; i <= j - 1; i = i + 1)
        {
              if (S1[i - 1] == S2[(j - 1) - i])
              {
                H[i][j - i] = H[i - 1][(j - 1) - i] + 2;
                Prev[i][j - i] = 2;
              }
              else
              {
                H[i][j - i] = H[i - 1][(j - 1) - i] - 1;
                Prev[i][j - i] = 2;
              }
              if (H[i][(j - 1) - i] + wp2[(j - i) - 1] > H[i][j - i])
              {
                H[i][j - i] = H[i][(j - 1) - i] + wp2[(j - i) - 1];
                Prev[i][j - i] = 1;
              }
              if (H[i][j - i] < 0)
              {
                Prev[i][j - i] = 0;
                H[i][j - i] = 0;
              }
              if (H[i - 1][j - i] + wp1[i - 1] > H[i][j - i])
              {
                H[i][j - i] = H[i - 1][j - i] + wp1[i - 1];
                Prev[i][j - i] = 3;
              }
        }
      }
      for (j = M+1; j <= N; j = j + 1)
      {
	    #pragma omp parallel for
        for (i = 1; i <= M; i = i + 1)
        {
              if (S1[i - 1] == S2[(j - 1) - i])
              {
                H[i][j - i] = H[i - 1][(j - 1) - i] + 2;
                Prev[i][j - i] = 2;
              }
              else
              {
                H[i][j - i] = H[i - 1][(j - 1) - i] - 1;
                Prev[i][j - i] = 2;
              }
              if (H[i][(j - 1) - i] + wp2[(j - i) - 1] > H[i][j - i])
              {
                H[i][j - i] = H[i][(j - 1) - i] + wp2[(j - i) - 1];
                Prev[i][j - i] = 1;
              }
              if (H[i][j - i] < 0)
              {
                Prev[i][j - i] = 0;
                H[i][j - i] = 0;
              }
              if (H[i - 1][j - i] + wp1[i - 1] > H[i][j - i])
              {
                H[i][j - i] = H[i - 1][j - i] + wp1[i - 1];
                Prev[i][j - i] = 3;
              }
        }
      }
  }
  if (N+1 > M+1)
      for (j = N+1; j <= M + N; j = j + 1)
      {
	    #pragma omp parallel for
        for (i = j - N; i <= M; i = i + 1)
        {
              if (S1[i - 1] == S2[(j - 1) - i])
              {
                H[i][j - i] = H[i - 1][(j - 1) - i] + 2;
                Prev[i][j - i] = 2;
              }
              else
              {
                H[i][j - i] = H[i - 1][(j - 1) - i] - 1;
                Prev[i][j - i] = 2;
              }
              if (H[i][(j - 1) - i] + wp2[(j - i) - 1] > H[i][j - i])
              {
                H[i][j - i] = H[i][(j - 1) - i] + wp2[(j - i) - 1];
                Prev[i][j - i] = 1;
              }
              if (H[i][j - i] < 0)
              {
                Prev[i][j - i] = 0;
                H[i][j - i] = 0;
              }
              if (H[i - 1][j - i] + wp1[i - 1] > H[i][j - i])
              {
                H[i][j - i] = H[i - 1][j - i] + wp1[i - 1];
                Prev[i][j - i] = 3;
              }
        }
      }
  else
  {
      for (j = N+1; j < M + 1; j = j + 1)
      {
	    #pragma omp parallel for
        for (i = j - N; i <= j - 1; i = i + 1)
        {
              if (S1[i - 1] == S2[(j - 1) - i])
              {
                H[i][j - i] = H[i - 1][(j - 1) - i] + 2;
                Prev[i][j - i] = 2;
              }
              else
              {
                H[i][j - i] = H[i - 1][(j - 1) - i] - 1;
                Prev[i][j - i] = 2;
              }
              if (H[i][(j - 1) - i] + wp2[(j - i) - 1] > H[i][j - i])
              {
                H[i][j - i] = H[i][(j - 1) - i] + wp2[(j - i) - 1];
                Prev[i][j - i] = 1;
              }
              if (H[i][j - i] < 0)
              {
                Prev[i][j - i] = 0;
                H[i][j - i] = 0;
              }
              if (H[i - 1][j - i] + wp1[i - 1] > H[i][j - i])
              {
                H[i][j - i] = H[i - 1][j - i] + wp1[i - 1];
                Prev[i][j - i] = 3;
              }
        }
      }
      for (j = M+1; j <= M + N; j = j + 1)
      {
	    #pragma omp parallel for
        for (i = j - N; i <= M; i = i + 1)
        {
              if (S1[i - 1] == S2[(j - 1) - i])
              {
                H[i][j - i] = H[i - 1][(j - 1) - i] + 2;
                Prev[i][j - i] = 2;
              }
              else
              {
                H[i][j - i] = H[i - 1][(j - 1) - i] - 1;
                Prev[i][j - i] = 2;
              }
              if (H[i][(j - 1) - i] + wp2[(j - i) - 1] > H[i][j - i])
              {
                H[i][j - i] = H[i][(j - 1) - i] + wp2[(j - i) - 1];
                Prev[i][j - i] = 1;
              }
              if (H[i][j - i] < 0)
              {
                Prev[i][j - i] = 0;
                H[i][j - i] = 0;
              }
              if (H[i - 1][j - i] + wp1[i - 1] > H[i][j - i])
              {
                H[i][j - i] = H[i - 1][j - i] + wp1[i - 1];
                Prev[i][j - i] = 3;
              }
        }
      }
  }
      
maximum = 0;
MaxLoc = (int*)malloc(4*sizeof(int));
MaxLoc = NULL;

for (i = 1; i<= M; i++)
	for( j = 1; j <= N; j++)
		if ( H[i][j] >= maximum ){
			maximum = H[i][j];
			printf("Maximum %i ",maximum);
			maxL = (int*)malloc(4*sizeof(int));
			maxL[0] = maximum;
			maxL[1] = i;
			maxL[2] = j;
			maxL[3] = (int)MaxLoc;
			MaxLoc = maxL;
	}

// выводим на экран построеннную матрицу Н
	printf("\n matrix H \n");
	for(i=0;i<M+1;i++){
		for(j=0;j<N+1;j++)
			printf(" %i ",H[i][j]);
		printf("\n");
	}

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


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];
}

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,ft

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);

system("pause");

    return 0;
}