#include <stdio.h>
#include <stdlib.h>
#include <math.h>
// #include <stdarg.h>
#include <string.h>
// #include <unistd.h>

#include "io_lib_header.h"
#include "util_lib_header.h"
#include "define_header.h"
#include "dp_lib_header.h"
#include "fastal_lib_header.h"
#include "fast_tree_header.h"


//TODO: seq_pair2diagonal delete num points from parameters
//TODO: reuse list


//Fastal_param *param_set;


/*! \mainpage T-Coffee Index Page
 *
 * \section intro_sec Introduction
 *
 * This is the introduction.
 *
 * \section install_sec Installation
 *
 * \subsection step1 Step 1: Opening the box
 *  
 * etc...
 * \section fastal_sec Fastal
 * 
 * This program is a very fast aligner. It is capable of aligning huge sets of sequences because it keeps as much as necessary on hard disk.
 */







/*!
 *	\file fastal.c
 *	\brief Source code for the fastal algorithm
 */


/**
 * \brief Calculates scores for diagonal segments.
 * 
 * \param seq1 Sequence 1
 * \param seq2 Sequence 2
 * \param *diagonals The diagonals. Three consecutive entries belong togehter. 1. pos in \a seq1 , 2. pos in \a seq2 and 3. length of diagonal
 * \param num_diagonals Number of diagonals
 * \param s1_length Length of \a seq1
 * \param list length of list.
 * \param list An 2-dim array to save the scores in.
 * \return new list
 */
int **
diag2pair_list(char* seq1,
			   char* seq2,
			   int *diagonals,
			   int num_diagonals,
			   int ***list_in,
			   int *current_length,
			   int *current_num_points,
			   int additional_needed,
			   Fastal_param *param_set)
{
	int **mat = param_set->M;
	int i, j, diag_length, pos1, pos2;
	int **list = list_in[0];

// 	printf("NUM: %i\n",num_diagonals);

	int l1 = strlen(seq1), l2 = strlen(seq2);
	int x = *current_num_points;


	int s1_length = strlen(seq1);
	int mini;
	for (i = 0; i < num_diagonals; ++i)
	{
		pos1 = diagonals[i*3];
		pos2 = diagonals[i*3+1];
		diag_length = diagonals[i*3+2];
		mini = MIN(pos1, pos2);
		pos1 -= mini;
		pos2 -= mini;
		while ((pos1 < l1) && (pos2 < l2))
		{
			if (x==*current_length)
			{
				*current_length+=1000;
				list=vrealloc (list,(*current_length)*sizeof(int*));
			}
			if (!list[x])
				list[x]=vcalloc (7, sizeof (int));
			list[x][0] = pos1+1;
			list[x][1] = pos2+1;
			list[x][2] = mat[toupper(seq1[pos1])-'A'][toupper(seq2[pos2])-'A'];

			++x;
			++pos1;
			++pos2;
		}
	}
	*current_num_points = x;
	list_in[0]=list;
}

void
guessalignment(Fastal_profile prf)
{
	
}

int 
fastal_compare (const void * a, const void * b)
{
	return (*(int*)a - *(int*)b);
}

int **
diagonals2int(int *diagonals,
			  int num_diagonals,
			  char *seq1,
			  char *seq2,
			  int *num_points,
			  Fastal_param *param_set)
{
	int l1 = strlen(seq1);
	int l2 = strlen(seq2);
	int gep = param_set->gep;
	
	int current_size = l2+l1;

	int **list = vcalloc(current_size, sizeof(int*));
	int *diags = vcalloc(num_diagonals, sizeof(int));
	int i;
// 	printf("SEQ: %s\nSEQ:%s\n",seq1, seq2);
// 	printf("X: %i\n",num_diagonals);
	for (i = 0; i < num_diagonals; ++i)
	{
		diags[i] = l1 - diagonals[i*3] + diagonals[i*3+1];
	}

	qsort (diags, num_diagonals, sizeof(int), fastal_compare);

	int *diagx = vcalloc(num_diagonals, sizeof(int));
	int *diagy = vcalloc(num_diagonals, sizeof(int));
	int *old_pos = vcalloc(num_diagonals, sizeof(int));

	//+1 because diagonals start here at position 1, like in "real" dynamic programming
	int a = 0, b = -1;
	for (i = 0; i < num_diagonals; ++i)
	{
		
		if (diags[i] < l1)
		{
			diagx[i] = l1 - diags[i];
			diagy[i] = 0;
			a= i;
		}
		else
			break;
	}
	++a;
	b=a-1;
	for (; i < num_diagonals; ++i)
	{
		diagx[i] = 0;
		diagy[i] = diags[i]-l1;
		b = i;
	}

	int tmpy_pos;
	int tmpy_value;
	int **M = param_set->M;
	int *last_y = vcalloc(l2+1, sizeof(int));
	int *last_x = vcalloc(l1+1, sizeof(int));
	last_y[0] = 0;
	
	last_x[0] = 0;
	list[0] = vcalloc(6, sizeof(int));

	int list_pos = 1;
	int dig_num = l1;
	int tmp_l2 = l2 + 1;
	
	//left border
	for (; list_pos < tmp_l2; ++list_pos)
	{
		list[list_pos] = vcalloc(6, sizeof(int));
		list[list_pos][0] = 0;
		list[list_pos][1] = list_pos;
		last_y[list_pos] = list_pos;
		list[list_pos][2] = list_pos*gep;
		list[list_pos][4] = list_pos-1;
	}

	int pos_x = 0;
	int diags_old = l2;
	
	int tmp = l1;
	int y;
	int tmp_l1 = l1-1;
	while (pos_x < tmp_l1)
	{
		if (list_pos + num_diagonals+2 > current_size)
		{
			current_size += num_diagonals*1000;
			list = vrealloc(list, current_size * sizeof(int*));
		}
		//upper border
		list[list_pos] = vcalloc(6, sizeof(int));
		list[list_pos][0] = ++pos_x;
		list[list_pos][1] = 0;
		list[list_pos][2] = pos_x * gep;
		list[list_pos][3] = last_y[0];
		tmpy_value = list_pos;
		tmpy_pos = 0;
		last_x[pos_x] = list_pos;
		++list_pos;

		//diagonals
		for (i = a; i <= b; ++i)
		{
			list[list_pos] = vcalloc(6, sizeof(int));
			
			list[list_pos][0] = ++diagx[i];
			
			list[list_pos][1] = ++diagy[i];
			list[list_pos][3] = last_y[diagy[i]];
			list[list_pos][4] = list_pos-1;
			list[list_pos][5] = last_y[diagy[i]-1];
			list[list_pos][2] = M[toupper(seq1[diagx[i]-1])-'A'][toupper(seq2[diagy[i]-1])-'A'];
			last_y[tmpy_pos] = tmpy_value;
			tmpy_value = list_pos;
			tmpy_pos = diagy[i];
			
			++list_pos;
		}
		last_y[tmpy_pos] = tmpy_value;

		
		//lower border
		if (list[list_pos-1][1] != l2)
		{
			list[list_pos] = vcalloc(6, sizeof(int));
			list[list_pos][0] = pos_x;
			list[list_pos][1] = l2;
			list[list_pos][3] = last_y[l2];
			
			list[list_pos][2] = -1000;
			list[list_pos][4] = list_pos-1;
			if (pos_x > l2)
				list[list_pos][5] = last_x[pos_x-l2];
			else
				list[list_pos][5] = l2-pos_x;
			last_y[l2] = list_pos;
			++list_pos;
			
		}


		if ((b >= 0) && (diagy[b] == l2))
			--b;
		
		if ((a >0) && (diagx[a-1] == pos_x))
			--a;
	}


	dig_num = -1;
	if (list_pos + l2+2 > current_size)
	{
		current_size += list_pos + l2 + 2;
		list = vrealloc(list, current_size * sizeof(int*));
	}
	
	
// 	right border	
	list[list_pos] = vcalloc(6, sizeof(int));
	list[list_pos][0] = l1;
	list[list_pos][1] = 0;
	list[list_pos][3] = last_x[l1-1];
	list[list_pos][2] = -1000;
	++list_pos;
	
	

	for (i = 1; i <= l2; ++i)
	{
		list[list_pos] = vcalloc(6, sizeof(int));
		list[list_pos][0] = l1;
		list[list_pos][1] = i;
		list[list_pos][3] = last_y[i];
		list[list_pos][4] = list_pos-1;
		y = last_y[i-1];
		if ((list[y][0] == l1-1) && (list[y][1] == i-1))
		{
			list[list_pos][5] = y;
			list[list_pos][2] = M[toupper(seq1[l1-1])-'A'][toupper(seq2[i-1])-'A'];
		}
		else
		{
			if (i <= l1)
			{
				list[list_pos][5] = last_x[l1-i];
			}
			else
			{
				list[list_pos][5] = i-l1;
			}
			list[list_pos][2] = -1000;
		}
		++list_pos;
	}
	
	list[list_pos - l2][2] = -1000;

	*num_points = list_pos;
	
	
// 	int blb;
// 	for (blb = 0; blb <list_pos; ++blb)
// 	{
//  	printf("LIST_A: %i %i %i %i %i %i %i %i\n",blb, list[blb][0],list[blb][1],list[blb][3],list[blb][2], list[blb][4], list[blb][5], list[blb][6]);
// 	}
	return list;
}
 
/**
 * \brief Makes a sorted list out of diagonals.
 * 
 * \param diagonals A list of diagonals to use during dynamic programming.
 * \param num_diagonals Number of diagonals.
 * \param seq1 Sequence 1.
 * \param seq2 Sequence 2.
 * \param gep cost for gap extension.
 * \param num_points Number of points in the list
 * \return A 2-dim array which contains all points needed for the sparse dynamic programming algorithm.
 */
int **
diagonals2int5(int *diagonals,
			  int num_diagonals,
			  char *seq1,
			  char *seq2,
			  int *num_points,
			  Fastal_profile *prf1,
			  Fastal_profile *prf2,
			  char *pos2char,
			  Fastal_param *param_set)
{

	int l1 = strlen(seq1);
	int l2 = strlen(seq2);

	int gep = param_set->gep;

	int current_size = l2+l1;
	int **list = vcalloc(current_size, sizeof(int*));
	int *diags = vcalloc(num_diagonals, sizeof(int));
	int i;
	for (i = 0; i < num_diagonals; ++i)
	{
		diags[i] = l1 - diagonals[i*3] + diagonals[i*3+1];
	
	}

 	qsort (diags, num_diagonals, sizeof(int), fastal_compare);

	int *diagx = vcalloc(num_diagonals, sizeof(int));
	int *diagy = vcalloc(num_diagonals, sizeof(int));
	int *old_pos = vcalloc(num_diagonals, sizeof(int));

	//+1 because diagonals start here at position 1, like in "real" dynamic programming
	int a = 0, b = -1;
	for (i = 0; i < num_diagonals; ++i)
	{
		
		if (diags[i] < l1)
		{
			
			diagx[i] = l1 - diags[i];
			diagy[i] = 0;

			a= i;
		}
		else
			break;
	}
	++a;
	b=a-1;
	for (; i < num_diagonals; ++i)
	{
			diagx[i] = 0;
			diagy[i] = diags[i]-l1;
			b = i;

	}

	int tmpy_pos;
	int tmpy_value;
	int **M = param_set->M;

	int *last_y = vcalloc(l2+1, sizeof(int));
	int *last_x = vcalloc(l1+1, sizeof(int));
	last_y[0] = 0;
	
	last_x[0] = 0;
	list[0] = vcalloc(6, sizeof(int));
// 	list[0][3] = l1;
	int list_pos = 1;
	int dig_num = l1;
	int tmp_l2 = l2 + 1;
	
	//left border
	for (; list_pos < tmp_l2; ++list_pos)
	{
		list[list_pos] = vcalloc(6, sizeof(int));
		list[list_pos][0] = 0;
		list[list_pos][1] = list_pos;
		last_y[list_pos] = list_pos;
		list[list_pos][2] = list_pos*gep;
		list[list_pos][3] = ++dig_num;
		list[list_pos][5] = list_pos-1;
	}

	int pos_x = 0;
	int diags_old = l2;
	
	int bla;
	int bla2, bla3, tmp_x;
	
	int tmp = l1;
	int y;
 	int tmp_l1 = l1-1;
	while (pos_x < tmp_l1)
	{
		if (list_pos + num_diagonals+2 > current_size)
		{
			current_size += num_diagonals*50;
			list = vrealloc(list, current_size * sizeof(int*));
		}
		//upper border
		list[list_pos] = vcalloc(6, sizeof(int));
		list[list_pos][0] = ++pos_x;
		list[list_pos][1] = 0;
		list[list_pos][2] = pos_x * gep;
		list[list_pos][3] = --tmp;
		list[list_pos][4] = last_y[0];
		tmpy_value = list_pos;
		tmpy_pos = 0;
		last_x[pos_x] = list_pos;
		++list_pos;

		//diagonals
		for (i = a; i <= b; ++i)
		{
			list[list_pos] = vcalloc(6, sizeof(int));
			list[list_pos][0] = ++diagx[i];
			list[list_pos][1] = ++diagy[i];
			list[list_pos][3] = diags[i];

			list[list_pos][4] = last_y[diagy[i]];
			list[list_pos][5] = list_pos-1;
			list[list_pos][6] = last_y[diagy[i]-1];
			
			list[list_pos][2] = 0;
			
			bla3 = 0;
			bla2 = 0;
			tmp_x = 0;
			for (bla = 0; bla<10; ++bla)
			{
				
				for (bla2 = 0; bla2<10; ++bla2)
				{
					bla3  += prf2->prf[bla2][diagy[i]-1] * prf1->prf[bla][diagx[i]-1];
					tmp_x += prf2->prf[bla2][diagy[i]-1] * prf1->prf[bla][diagx[i]-1] * M[pos2char[bla]-'A'][pos2char[bla2] -'A'];

				}
			}
			list[list_pos][2] = (int)tmp_x / bla3;
			
// 			for (bla = 0; bla<10; ++bla)
// 				bla2 += prf2->prf[bla][diagy[i]-1];
// 			bla2 = bla2/prf2->num_sequences;
// 			
// 			for (bla = 0; bla<10; ++bla)
// 				bla3 += prf1->prf[bla][diagy[i]-1];
//
// 			bla3 = bla3/prf1->num_sequences;
//
//
// 			if ((bla2 > 0.7) && (bla3 > 0.7))
// 				list[list_pos][2] = M[toupper(seq1[diagx[i]-1])-'A'][toupper(seq2[diagy[i]-1])-'A'];
// 			else if ((bla< 0.7) && (bla3 < 0.7))
// 				list[list_pos][2] = M[toupper(seq1[diagx[i]-1])-'A'][toupper(seq2[diagy[i]-1])-'A'] = 3;
// 			else
// 				list[list_pos][2] = M[toupper(seq1[diagx[i]-1])-'A'][toupper(seq2[diagy[i]-1])-'A'] * ((bla< 0.7) && (bla3 < 0.7));
//  			list[list_pos][2] = M[toupper(seq1[diagx[i]-1])-'A'][toupper(seq2[diagy[i]-1])-'A'];//* ((bla2+bla3)/2);
			last_y[tmpy_pos] = tmpy_value;
			tmpy_value = list_pos;
			tmpy_pos = diagy[i];

			++list_pos;
		}
		last_y[tmpy_pos] = tmpy_value;


		//lower border
		if (list[list_pos-1][1] != l2)
		{
			list[list_pos] = vcalloc(6, sizeof(int));
			list[list_pos][0] = pos_x;
			list[list_pos][1] = l2;
			list[list_pos][4] = last_y[l2];
			
			list[list_pos][2] = -1000;
			list[list_pos][3] = l1 - pos_x + l2;
			list[list_pos][5] = list_pos-1;
			if (pos_x > l2)
				list[list_pos][6] = last_x[pos_x-l2];
			else
				list[list_pos][6] = l2-pos_x;
			last_y[l2] = list_pos;
			++list_pos;
		}


		if ((b >= 0) && (diagy[b] == l2))
			--b;
		
		if ((a >0) && (diagx[a-1] == pos_x))
			--a;
	}


	dig_num = -1;
	if (list_pos + l2+2 > current_size)
	{
		current_size += list_pos + l2 + 2;
		list = vrealloc(list, current_size * sizeof(int*));
	}


// 	right border	
	list[list_pos] = vcalloc(6, sizeof(int));
	list[list_pos][0] = l1;
	list[list_pos][1] = 0;
	list[list_pos][3] = ++dig_num;
	list[list_pos][4] = last_x[l1-1];
	list[list_pos][2] = -1000;
	++list_pos;

	for (i = 1; i <= l2; ++i)
	{
		list[list_pos] = vcalloc(6, sizeof(int));
		list[list_pos][0] = l1;
		list[list_pos][1] = i;
		list[list_pos][3] = ++dig_num;
		list[list_pos][4] = last_y[i];
		list[list_pos][5] = list_pos-1;
		y = last_y[i-1];
		if ((list[y][0] == l1-1) && (list[y][1] == i-1))
		{
			list[list_pos][6] = y;
			list[list_pos][2] = M[toupper(seq1[l1-1])-'A'][toupper(seq2[i-1])-'A'];
		}
		else
		{
			if (i <= l1)
			{
				list[list_pos][6] = last_x[l1-i];
			}
			else
			{
				list[list_pos][6] = i-l1;
			}
			list[list_pos][2] = -1000;
		}
		++list_pos;
	}
	
	list[list_pos - l2][2] = -1000;

	*num_points = list_pos;

	return list;
}
 



//**************************   sparse dynamic aligning **********************************************************


void
combine_profiles2file(int **prf1,
					  int **prf2,
					  int pos1,
					  int pos2,
					  Fastal_param *param_set,
					  FILE *prof_f,
					  char state)
{
	int alphabet_size = param_set->alphabet_size;
	char *pos2aa = &(param_set->pos2char[0]);	
	int i;
	int x = 0;
	if (state == 'M')
	{
		for (i = 0; i < alphabet_size; ++i)
			if (prf1[i][pos1] + prf2[i][pos2] > 0)
			{
				if (x)
					fprintf(prof_f," %c%i", pos2aa[i],prf1[i][pos1]+prf2[i][pos2]);
				else
					fprintf(prof_f,"%c%i", pos2aa[i],prf1[i][pos1]+prf2[i][pos2]);
				x = 1;
			}
		fprintf(prof_f,"\n");
	}
	else if (state == 'D')
	{
		for (i = 0; i < alphabet_size; ++i)
			if (prf2[i][pos2] > 0)
		{
			if (x)
				fprintf(prof_f," %c%i", pos2aa[i],prf2[i][pos2]);
			else
				fprintf(prof_f,"%c%i", pos2aa[i],prf2[i][pos2]);
			x = 1;
		}
		fprintf(prof_f,"\n");
	}
	else
	{
		for (i = 0; i < alphabet_size; ++i)
			if (prf1[i][pos1] > 0)
		{
			if (x)
				fprintf(prof_f," %c%i", pos2aa[i],prf1[i][pos1]);
			else
				fprintf(prof_f,"%c%i", pos2aa[i],prf1[i][pos1]);
			x = 1;
		}
		fprintf(prof_f,"\n");
	}
}



#define LIN(a,b,c) a[b*5+c]
/**
 * Calculates a fast and sparse dynamic programming matrix
 * 
 * \param prf1 Profile of first sequence.
 * \param prf2 Profile of second sequence.
 * \param param_set The parameter for the alignment.
 * \param list The list of diagonals.
 * \param n number of dots.
 * \param edit_f File to save the edit information.
 * \param prof_f File to save the profile.
 * \param node_number Number of the new profile.
 */
int
list2linked_pair_wise_fastal(Fastal_profile *prf1,
							 Fastal_profile *prf2,
							 Fastal_param *param_set,
							 int **list,
							 int n,
							 FILE *edit_f,
							 FILE *prof_f,
							 int node_number)
{
	int a,b,c, i, j, LEN=0, start_trace;
	int pi, pj,ij, delta_i, delta_j, prev_i, prev_j;
	static int **slist;
	static long *MI, *MJ, *MM,*MT2;
	static int *sortseq;
	static int max_size;
	int gop, gep, igop, igep;
	int l1, l2, l, ls;
	char **al;
	char **aln,*char_buf;
	int ni=0, nj=0;
	long score;
	int nomatch = param_set->nomatch;
 
	l1=prf1->length;
	l2=prf2->length;

	al=declare_char (2,l1+l2+1);


	
	igop=param_set->gop;
	gep=igep=param_set->gep;
	if (n>max_size)
	{
		max_size=n;

		vfree (MI);vfree (MJ); vfree (MM);
		free_int (slist, -1);

		slist=declare_int (n,3);

		MI=vcalloc (5*n, sizeof (long));
		MJ=vcalloc (5*n, sizeof (long));
		MM=vcalloc (5*n, sizeof (long));

	}
	else
	{
		for (a=0; a<n; a++)
			for (b=0; b<5; b++)LIN(MI,a,b)=LIN(MJ,a,b)=LIN(MJ,a,b)=-1000000;
	}
  
	for (a=0; a<n; a++)
	{
		i=list[a][0];
		j=list[a][1];


		if (i==l1 || j==l2)gop=0;
		else gop=igop;

		if (i==l1 && j==l2)start_trace=a;
		else if ( i==0 || j==0)
		{
			LIN(MM,a,0)=-1000000;
			if (j==0)
			{
				LIN(MJ,a,0)=-10000000;
				LIN(MI,a,0)=gep*i;
			}
			else if (i==0)
			{
				LIN(MI,a,0)=-10000000;
				LIN(MJ,a,0)=gep*j;
			}

			LIN(MI,a,1)=LIN(MJ,a,1)=-1;
			LIN(MI,a,2)=LIN(MJ,a,2)=i;
			LIN(MI,a,3)=LIN(MJ,a,3)=j;
			continue;
		}

		pi = list[a][3];
		pj = list[a][4];
		ij = list[a][5];

		prev_i=list[pi][0];
		prev_j=list[pj][1];

		delta_i=list[a][0]-list[pi][0];
		delta_j=list[a][1]-list[pj][1];

		/*Linear Notation*/
		LIN(MI,a,0)=MAX(LIN(MI,pi,0),(LIN(MM,pi,0)+gop))+delta_i*gep;
		LIN(MI,a,1)=pi;
		LIN(MI,a,2)=delta_i;
		LIN(MI,a,3)=0;
		LIN(MI,a,4)=(LIN(MI,pi,0) =(LIN(MM,pi,0)+gop))?'i':'m';

		LIN(MJ,a,0)=MAX(LIN(MJ,pj,0),(LIN(MM,pj,0)+gop))+delta_j*gep;
		LIN(MJ,a,1)=pj;
		LIN(MJ,a,2)=0;
		LIN(MJ,a,3)=delta_j;

		LIN(MJ,a,4)=(LIN(MJ,pj,0) =LIN(MM,pj,0)+gop)?'j':'m';

		if (a>1 && (ls=list[a][0]-list[ij][0])==(list[a][1]-list[ij][1]))
		{
			LIN(MM,a,0)=MAX3(LIN(MM,ij,0),LIN(MI,ij,0),LIN(MJ,ij,0))+list[a][2]-(ls*nomatch);

			LIN(MM,a,1)=ij;
			LIN(MM,a,2)=ls;
			LIN(MM,a,3)=ls;
			if ( LIN(MM,ij,0)>=LIN(MI,ij,0) && LIN(MM,ij,0)>=LIN(MJ,ij,0))LIN(MM,a,4)='m';
			else if ( LIN(MI,ij,0) >= LIN(MJ,ij,0))LIN(MM,a,4)='i';
			else LIN(MM,a,4)='j';
	  
		}
		else
		{
			LIN(MM,a,0)=UNDEFINED;
			LIN(MM,a,1)=-1;
		}  
	}
  
	a=start_trace;
	if (LIN(MM,a,0)>=LIN(MI,a,0) && LIN(MM,a,0) >=LIN(MJ,a,0))MT2=MM;
	else if ( LIN(MI,a,0)>=LIN(MJ,a,0))MT2=MI;
	else MT2=MJ;

	score=MAX3(LIN(MM,a,0), LIN(MI,a,0), LIN(MJ,a,0));
  
	i=l1;
	j=l2;
  
  
	while (!(i==0 &&j==0))
	{
		int next_a;
		l=MAX(LIN(MT2,a,2),LIN(MT2,a,3));
      // HERE ("%c from %c %d %d SCORE=%d [%d %d] [%2d %2d]", T2[a][5],T2[a][4], T2[a][2], T2[a][3], T2[a][0], gop, gep, i, j);
		if (i==0)
		{
			while ( j>0)
			{
				al[0][LEN]=0;
				al[1][LEN]=1;
				j--; LEN++;
			}
		}
		else if (j==0)
		{
			while ( i>0)
			{
				al[0][LEN]=1;
				al[1][LEN]=0;
				i--; LEN++;
			}
		}
      
		else if (l==0) {HERE ("L=0 i=%d j=%d",l, i, j);exit (0);}
		else 
		{
			for (b=0; b<l; b++, LEN++)
			{
				if (LIN(MT2,a,2)){al[0][LEN]=1;i--;ni++;}
				else al[0][LEN]=0;
	      
				if (LIN(MT2,a,3)){al[1][LEN]=1;j--;nj++;}
				else al[1][LEN]=0;
			}
	  
			next_a=LIN(MT2,a,1);
			if (LIN(MT2,a,4)=='m')MT2=MM;
			else if (LIN(MT2,a,4)=='i')MT2=MI;
			else if (LIN(MT2,a,4)=='j')MT2=MJ;
			a=next_a;
		}
	}

	invert_list_char ( al[0], LEN);
	invert_list_char ( al[1], LEN);

	fprintf(edit_f, "%i\n%i\n%i\n%i\n",prf1->prf_number, prf2->prf_number, prf1->is_leaf, prf2->is_leaf);
	fprintf(prof_f, "%i\n0\n%i\n1\n", node_number,LEN);

	char statec[] = {'M','D','I'};
	int num = 0;
	int state = 0;
	i = 0;
	j = 0;

	for ( b=0; b< LEN; b++)
	{
		if ((al[0][b]==1) && (al[1][b]==1))
		{

			combine_profiles2file(prf1->prf, prf2->prf, i, j, param_set, prof_f, 'M');
			++i;
			++j;
			if (state != 0)
			{
				fprintf(edit_f, "%c%i\n",statec[state], num);
				num =1;
				state = 0;
			}
			else
				++num;
		}
		else if (al[0][b]==1)
		{
// 			prf1->prf[param_set->alphabet_size-1] += prf2->num_sequences;
			combine_profiles2file(prf1->prf, prf2->prf, i, j, param_set, prof_f, 'I');
			++i;
			if (state != 2)
			{
				fprintf(edit_f, "%c%i\n",statec[state], num);
				num =1;
				state = 2;
			}
			else
				++num;
		} 
		else if (al[1][b]==1)
		{
//			prf2->prf[param_set->alphabet_size-1] += prf1->num_sequences;
			combine_profiles2file(prf1->prf, prf2->prf, i, j, param_set, prof_f, 'D');
			++j;
			if (state != 1)
			{
				fprintf(edit_f, "%c%i\n",statec[state], num);
				num =1;
				state = 1;
			}
			else
				++num;
		}
	}
	

	fprintf(edit_f, "%c%i\n",statec[state], num);
	num =1;
	state = 1;
	
	
	fprintf(edit_f,"*\n");
	fprintf(prof_f,"*\n");
	free_char (al, -1);
//    exit(0);
	return LEN;
}






/**
 * \brief Tuns a profile into a consensus sequence.
 * 
 * The character with the highest number of occurences is used as consensus. Gaps are not included. For example: 10 '-' and one 'A' would give 'A' as consensus.
 * \param profile The profile.
 * \param file_name Name of the file to save the consensus sequence in.
 * \param param_set The parameter of the fastal algorithm.
 * \return the sequence
 */
char*
profile2consensus(Fastal_profile *profile, char *file_name, Fastal_param *param_set)
{
	FILE *cons_f = fopen(file_name,"w");
	fprintf(cons_f, ">%i\n", profile->prf_number);
	char* seq = vcalloc(profile->length+1, sizeof(char));
	int i, j;
	int most_pos, most;
	int alphabet_size = param_set->alphabet_size;
	int **prf = profile->prf;
	char *pos2char = param_set->pos2char;
	for (i = 0; i < profile->length; ++i)
	{
		most = -1;
		for (j = 0; j < alphabet_size; ++j)
		{
			if (prf[j][i] > most)
			{
				most = prf[j][i];
				most_pos = j;
			}
		}
		seq[i] = pos2char[most_pos];
		fprintf(cons_f, "%c",pos2char[most_pos]);
	}
	seq[i] = '\0';
	fprintf( cons_f, "\n");
	fclose(cons_f);
	return seq;
}




/**
 * \brief Calculates the diagonals between two sequences.
 * 
 * Uses bl2seq to calculate the diagonals.
 * \param seq_file1 File with sequence 1.
 * \param seq_file2 File with sequence 2.
 * \param diagonals An array where the diagonal points will be stored.
 * \param dig_length length of \a diagonals .
 * \param num_points Number of points in all diagonals.
 * \return number of diagonals;
 */
int
seq_pair2blast_diagonal(char *seq_file_name1,
						char *seq_file_name2,
						int **diagonals,
						int *dig_length,
						int l1,
						int l2,
						int is_dna)
{
	int *diag = vcalloc(l1 + l2, sizeof(int));
	char *out_file = vtmpnam(NULL);
	char blast_command[600];

	if (is_dna)
		sprintf(blast_command, "bl2seq -p blastn -i %s -j %s -D 1 -g F -o %s", seq_file_name1, seq_file_name2, out_file);
	else
 		sprintf(blast_command, "bl2seq -p blastp -i %s -j %s -D 1 -g F -o %s", seq_file_name1, seq_file_name2, out_file);
	system(blast_command);

	int *diags = diagonals[0];
	FILE *diag_f = fopen(out_file,"r");
	char line[300];
	fgets(line, 300, diag_f);
	fgets(line, 300, diag_f);
	fgets(line, 300, diag_f);
	
	
	char delims[] = "\t";
	char *result = NULL;
	int length, pos_q, pos_d, i;
	int current_pos = 0;
	while (fgets(line, 300, diag_f) != NULL)
	{
		strtok(line, delims);
		strtok(NULL, delims);
		strtok(NULL, delims);
		length =  atoi(strtok(NULL, delims));
		strtok(NULL, delims);
		strtok(NULL, delims);
		pos_q = atoi(strtok(NULL, delims))-1;
		strtok(NULL, delims);
		pos_d = atoi(strtok(NULL, delims))-1;

		if (current_pos >= *dig_length)
		{
			(*dig_length) += 90;
			diags = vrealloc(diags, sizeof(int)*(*dig_length));
		}
		if (diag[l1-pos_q+pos_d] == 0)
		{
 			diag[l1-pos_q+pos_d] =1;
			diags[current_pos++] = pos_q;
			diags[current_pos++] = pos_d;
			diags[current_pos++] = length;
		}
	}
	vfree(diag);
	fclose(diag_f);
	diagonals[0] = diags;
	return current_pos/3;
}




//******************************* OTHER STUFF ***********************

/**
 *	\brief Reads the sequence from a given position in a fasta file and turns it into a profile.
 * 
 * \param seq_file The file where the sequence is stored.
 * \param off_set The off_set from the beginning of the file to the position of the sequence name.
 * \param profile The profile where the sequence will be stored into.
 * \param prf_number The number of this profile.
 */
void
file_pos2profile(FILE *seq_file,			//File with sequences
				 long off_set,				//offset of sequence from the beginning of file point to the sequence name, not to the sequence itself
				 Fastal_profile *profile,	//profile to save into
				 int prf_number,			//number of the profile
				 Fastal_param *param_set)			
{
	int alphabet_size = param_set->alphabet_size;
	profile->is_leaf = 1;
	int *aa2pos = &(param_set->char2pos[0]);
	const int LINE_LENGTH = 500;
 	char line[LINE_LENGTH];
	profile->num_sequences = 1;
	profile->prf_number = prf_number;
	fseek (seq_file , off_set , SEEK_SET );
	
	fgets (line, LINE_LENGTH , seq_file);
	int seq_length = 0;
	int i, j;

	while(fgets(line, LINE_LENGTH, seq_file)!=NULL)
	{
		if (line[0] != '>')
		{

			line[LINE_LENGTH-1] = '\n';
			if (seq_length + LINE_LENGTH >= profile->allocated_memory)
			{
				for (i = 0; i < alphabet_size; ++i)
				{
					profile->prf[i] = vrealloc(profile->prf[i], (profile->allocated_memory+PROFILE_ENLARGEMENT)*sizeof(int));
				}
				profile->allocated_memory += PROFILE_ENLARGEMENT;
			}

			i = 0;
			while (line[i] != '\n')
			{
				for(j = 0; j<alphabet_size; ++j )
					profile->prf[j][seq_length+i] = 0;
				profile->prf[aa2pos[toupper(line[i])-'A']][seq_length+i] = 1;
				++i;
			}
			seq_length += i;

		}
		else
			break;
	}
	profile->length = seq_length;
}



/**
*	constructs index of fasta_file
*/
int
make_index_of_file(char *file_name, 		//file with sequences
				   long **file_positions)	//array to save the positions
{
	const int LINE_LENGTH = 150;
	(*file_positions) = vcalloc(ENLARGEMENT_PER_STEP,  sizeof(long));

	int current_size = ENLARGEMENT_PER_STEP;
	int current_pos = 0;
	
	FILE *file = fopen(file_name,"r");
	
	char *sequence = vcalloc(3*LINE_LENGTH,sizeof(char));
	int seq_length=0;
	int allocated_length=3*LINE_LENGTH;
	char line[LINE_LENGTH];

	int num_of_sequences = 0;
	int mem_for_pos = ENLARGEMENT_PER_STEP;
	
	if (file == NULL)
	{
		printf("FILE NOT FOUND\n");
		exit(1);
	}
	else
	{
		(*file_positions)[num_of_sequences] = ftell(file);
		while(fgets(line, LINE_LENGTH , file)!=NULL)
		{
			int length = strlen(line);
			if (line[0] == '>')
			{
				++num_of_sequences;
				
				if (num_of_sequences == mem_for_pos)
				{
					(*file_positions) = vrealloc((*file_positions),(ENLARGEMENT_PER_STEP+mem_for_pos) * sizeof(long));
					mem_for_pos += ENLARGEMENT_PER_STEP;
				}
			}
			(*file_positions)[num_of_sequences] = ftell(file);
		}
	}
	fclose(file);
	return num_of_sequences;
}


/**
*	reads a profile from a profile_file
*/
profile_file2profile(Fastal_profile *prof,	//structure to save the profile in
					 FILE *profile_f,		//file where the profile is stored
					 long position,			//position in profile_f where the profile is stored
					 Fastal_param *param_set)
{
	
	int alphabet_size = param_set->alphabet_size;
			
	int *aa2pos = &(param_set->char2pos[0]);


	fseek(profile_f,position,SEEK_SET);
	const int LINE_LENGTH = 500;
	char line[500];
	
	fgets(line, LINE_LENGTH, profile_f);

	prof->prf_number = atoi(line);
// 	fgets(line, LINE_LENGTH, profile_f);
// 	prof->num_sequences = atoi(line);
// 	fgets(line, LINE_LENGTH, profile_f); //is-dna is already known
	fgets(line, LINE_LENGTH, profile_f);
	prof->is_leaf = atoi(line);

	fgets(line, LINE_LENGTH, profile_f);
	prof->length = atoi(line);
	fgets(line, LINE_LENGTH, profile_f);
	prof->weight = atoi(line);
	int i,j;
	if (prof->length > prof->allocated_memory)
		for (i = 0;i < alphabet_size; ++i)
	{
		prof->prf[i] = vrealloc(prof->prf[i],prof->length*sizeof(int));
	}
	
	char delims[] = " ";
	char *result = NULL;
	char *result_num = NULL;
	
	int length = prof->length;

	for (i = 0; i < length; ++i)
	{
		for(j = 0; j<alphabet_size; ++j )
			prof->prf[j][i] = 0;
		fgets(line, LINE_LENGTH , profile_f);
		result = strtok( line, delims );
		
		while( result != NULL)
		{
			result_num = &result[1];
			prof->prf[aa2pos[result[0]-'A']][i] = atoi(result_num);
			result = strtok( NULL, delims );
		}
	}
}



/**
*	writes a profile into a file
*/
void 
profile2file(Fastal_profile *profile,	//the profile to save
			 FILE* file,				//file to save in
			 Fastal_param *param_set)
{
	int alphabet_size = param_set->alphabet_size;

	char *pos2aa = &(param_set->pos2char[0]);

	fseek(file,0,SEEK_SET);
	
	fprintf(file,"%i\n", profile->prf_number);
// 	fprintf(file,"%i\n", profile->num_sequences);
	
	fprintf(file,"%i\n", profile->is_leaf);
	fprintf(file,"%i\n", profile->length);
	fprintf(file,"%i\n", profile->weight);
	int i = 0, j = 0;
	int max = profile->length;
	int x= 0;
	--alphabet_size;
	while (i < max)
	{
		for (j = 0; j < alphabet_size; ++j)
			if (profile->prf[j][i] > 0)
			{
				if (x)
					fprintf(file," %c%i", pos2aa[j],profile->prf[j][i]);
				else
					fprintf(file,"%c%i", pos2aa[j],profile->prf[j][i]);
					x = 1;
			}
		if (profile->prf[j][i] > 0)
			if (x)
				fprintf(file," %c%i", pos2aa[j],profile->prf[j][i]);
			else
				fprintf(file,"%c%i", pos2aa[j],profile->prf[j][i]);
			x = 1;
		x = 0;
		fprintf(file,"\n");
		++i;
	}
	fprintf(file,"*\n");
}



/**
*	Reads the profile out of an alignment
*/
void
file2profile(FILE* profile_f,			//file to read the profile of
			 Fastal_profile *prof,		//profile saved in here
			 int prf_number,			//number of the profile
			 Fastal_param *param_set)
{
	int alphabet_size = param_set->alphabet_size;

	int *aa2pos =  &(param_set->char2pos[0]);


	fseek(profile_f,0,SEEK_SET);
	const int LINE_LENGTH = 500;
	char line[500];
	
	fgets(line, LINE_LENGTH, profile_f);
	prof->prf_number = atoi(line);
// 	fgets(line, LINE_LENGTH, profile_f); //is-dna is already known
	fgets(line, LINE_LENGTH, profile_f);
	prof->is_leaf = atoi(line);

	fgets(line, LINE_LENGTH, profile_f);
	prof->length = atoi(line);
	
	fgets(line, LINE_LENGTH, profile_f);
	prof->weight = atoi(line);
	int i,j;
	if (prof->length > prof->allocated_memory)
		for (i = 0;i < alphabet_size; ++i)
		{
			prof->prf[i] = vrealloc(prof->prf[i],prof->length*sizeof(int));
		}
	
	char delims[] = " ";
	char *result = NULL;
	char *result_num = NULL;
	
	int length = prof->length;

	for (i = 0; i < length; ++i)
	{
		for(j = 0; j<alphabet_size; ++j )
			prof->prf[j][i] = 0;
		fgets(line, LINE_LENGTH , profile_f);
		result = strtok( line, delims );
		
		while( result != NULL)
		{
			result_num = &result[1];
			prof->prf[aa2pos[result[0]-'A']][i] = atoi(result_num);
			result = strtok( NULL, delims );
		}
	}
}



/**
*	This method takes a profile and turns it into a sumed up version of same size.
*/
int**
sumup_profile(Fastal_profile *profile,	//profile to sum-up
			  int **sumup,
			  Fastal_param *param_set)	//summed_up_profile
{
	
	char *pos2aa = &(param_set->pos2char[0]);
	int alphabet_size = param_set->alphabet_size;
	int **M = param_set->M;
	int prof_length = profile->length;
	
	int i,j,k;

	for (i = 0; i < prof_length; ++i)
	{
		sumup[alphabet_size][i] = 0;
		for (k = 0; k < alphabet_size; ++k)
		{
			sumup[k][i] = 0;
			sumup[alphabet_size][i] += profile->prf[k][i];
			for (j = 0; j < alphabet_size; ++j)
			{
				sumup[k][i] += profile->weight * profile->prf[j][i] * M[pos2aa[j]-'A'][pos2aa[k]-'A'];
			}
		}
	}

	return sumup;
}



/**
*	Turns the dynamic programming matrix into a editfile and calculates the new profile
*/
int
nw_matrix2edit_file(double **prog_matrix,	//dynamic programming matrix
					Fastal_profile *prf1,	//profile of dim1
					Fastal_profile *prf2,	//profile of dim2
	 				FILE *edit_f,			//file to safe the edit in
	  				int **prf_field,		//space to safe the new profile
					int *field_length,
					Fastal_param *param_set)		//length of prf_field
{
	int **M = param_set->M;
	int alphabet_size = param_set->alphabet_size;
	double gap_cost = param_set -> gop;
	fprintf(edit_f, "%i\n%i\n%i\n%i\n",prf1->prf_number, prf2->prf_number, prf1->is_leaf, prf2->is_leaf);
	int sum[] = {0,0,0};
	char sumc[] = {'M','I','D'};
	int last = 0;
	int n = 0;
	int m = 0;
	int field_pos = 0;
	int i;
	int prf1_length = prf1->length;
	int prf2_length = prf2->length;
	while ((n < prf1_length) && (m < prf2_length))
	{
		//if necesarry allocate more memory for result
		if ((*field_length)-alphabet_size < field_pos)
		{
			(*field_length) += ENLARGEMENT_PER_STEP;
			
			for (i = 0; i <alphabet_size+1; ++i)
			{
				prf_field[i] = vrealloc(prf_field[i], (*field_length)*sizeof(int));
			}
		}
		
		if (prog_matrix[n][m] == (prog_matrix[n+1][m] +gap_cost))
		{
			for (i = 0; i<alphabet_size; ++i)
			{
				prf_field[i][field_pos] = prf1->prf[i][n];
			}
			++n;
			++ field_pos;
			
			if (last != 1)
			{
				fprintf(edit_f,"%c%i\n",sumc[last],sum[last]);
				sum[last] = 0;
			}
			last = 1;
			++sum[last];
		}
		else if (prog_matrix[n][m] == (prog_matrix[n][m+1] +gap_cost))
		{
			
			for (i = 0; i<alphabet_size; ++i)
			{
				prf_field[i][field_pos] = prf2->prf[i][m];
			}
			++m;
			++ field_pos;
			if (last != 2)
			{
				fprintf(edit_f,"%c%i\n",sumc[last],sum[last]);
				sum[last] = 0;
			}
			last = 2;
			++sum[last];
		}
		else 
		{
			for (i = 0; i<alphabet_size; ++i)
			{
				prf_field[i][field_pos] = prf1->prf[i][n] + prf2->prf[i][m];
			}
			++n;
			++m;
			++ field_pos;
			if (last != 0)
			{
				fprintf(edit_f,"%c%i\n",sumc[last],sum[last]);
				sum[last] = 0;
			}
			last = 0;
			++sum[last];
		}
	}
	fprintf(edit_f,"%c%i\n",sumc[last],sum[last]);
	
	//gaps in prf2
	last = 0;
	while (n < prf1_length)
	{
		for (i = 0; i<alphabet_size; ++i)
		{
			prf_field[i][field_pos] = prf1->prf[i][n];
		}
		++n;
		++ field_pos;
		++last;
	}
	if (last > 0)
		fprintf(edit_f,"I%i\n",last);
	
	//gaps in prf1
	last = 0;
	while (m < prf2_length)
	{
		for (i = 0; i<alphabet_size; ++i)
		{
			prf_field[i][field_pos] = prf2->prf[i][m];
		}
		++m;
		++ field_pos;
		++last;
	}
	if (last > 0)
		fprintf(edit_f,"D%i\n",last);
	fprintf(edit_f,"*\n");
	return field_pos;
}




/**
 * \brief Pairwise alignments of profile is done here.
 *
 * \param profile1 Profile of sequence 1
 * \param profile2 Profile of sequence 2
 * \param prog_matrix Matrix for dynamic programming
 * \param edit_file_name The edit_file_name
 * \param sumup_prf The sumup version of profile 1, which later contains the aligned profile.
 * \param sumup_length Contains the length of the aligned profile.
 * \return length of the aligned profile
 */
int
prf_nw(Fastal_profile *profile1,	//profile of sequence 1
	   Fastal_profile *profile2,	//profile of sequence 2
	   double **prog_matrix,		//matrix for dynamic programming (at least as long as necessary for alignment)
	   FILE *edit_file_name,		//name of edit file
	   int **sumup_prf,				//sum_up
	   int *sumup_length, 
	   Fastal_param *param_set)			//sum_up length
{
	int alphabet_size = param_set->alphabet_size;
	double gap_cost = param_set->gop;
	
	int i;
	if (*sumup_length < profile1->length)
	{
		for (i = 0; i < alphabet_size+1; ++i)
		{
			sumup_prf[i] = vrealloc(sumup_prf[i], profile1->length*sizeof(int));
		}
		*sumup_length = profile1->length;
	}
	sumup_prf = sumup_profile(profile1, sumup_prf, param_set);
	
	

	int j,k;
	int prof1_length = profile1->length;
	int prof2_length = profile2->length;

	int** M = param_set->M;
	double match_score;
	int amino_counter;
	int residue_pairs = 0;

	for (i = prof2_length; i > 0; --i)
	{
		prog_matrix[prof1_length][i] = gap_cost * (prof2_length-i);
	}

	i = prof1_length-1;
	prog_matrix[prof1_length][prof2_length] = 0.0;
	while (i >=0)
	{
		j = prof2_length-1;

		prog_matrix[i][prof2_length] = gap_cost*(prof1_length-i);
		while (j >=0)
		{
			match_score = 0.0;
			residue_pairs = 0;
			for (k = 0; k < alphabet_size; ++k)
			{
				residue_pairs += profile2->prf[k][j];
				match_score += (profile2->prf[k][j] * sumup_prf[k][i]);
			}
			match_score /= (residue_pairs * sumup_prf[alphabet_size][i]);
			prog_matrix[i][j] = MAX3(prog_matrix[i+1][j+1]+match_score, prog_matrix[i+1][j]+gap_cost, prog_matrix[i][j+1]+gap_cost);
			
			--j;
		}
		--i;
	}
	return nw_matrix2edit_file(prog_matrix, profile1, profile2, edit_file_name, sumup_prf, sumup_length, param_set);
}



/**
 * \brief Writes the sequence into the alignment_file.
 * 
 * \param aligned_sequence Pattern of aligned sequence.
 * \param sequence_file File with sequences.
 * \param sequence_position Positions of sequences in \a sequence_file. 
 * \param alignment_file The file to write the sequence into.
 * 
*/
void
edit_seq2aligned_seq(char *aligned_sequence,	//pattern for aligned sequence
					 FILE *sequence_file,		//file with all the sequences
					 long sequence_position,	//position in sequence file with the correct sequence
					 FILE *alignment_file)		//file to write the alignment into
{
	fseek(sequence_file, sequence_position, SEEK_SET);
	const int LINE_LENGTH = 300;
	char line[LINE_LENGTH];
	fgets (line, LINE_LENGTH , sequence_file);
	fprintf(alignment_file,"%s", line);	//writing of sequence name
	int pos = 0;
	int i = 0;
	while(fgets(line, LINE_LENGTH, sequence_file)!=NULL)
	{
		if (line[0] != '>')
		{
			line[LINE_LENGTH-1] = '\n';
			i = 0;
			while (line[i] != '\n')
			{
				while (aligned_sequence[pos] == '-')
				{
					fprintf(alignment_file,"-");
					++pos;
				}
				fprintf(alignment_file,"%c",line[i]);
				++i;
				++pos;
			}
		}
		else
			break;
	}
	while (aligned_sequence[pos] != '\n')
	{
		fprintf(alignment_file,"-");
		++pos;
	}
	fprintf(alignment_file,"\n");
}



/**
 * \brief Recursive function to turn the edit_file into the alignment.
 * 
 * \param sequence_file File with all sequences.
 * \param sequence_position The array of sequence positions in \a sequence_file
 * \param edit_file File to safe the edit profiles in.
 * \param edit_positions Array saving the coorespondence between edit profile and position in \a edit_file
 * \param node_number The current node.
 * \param number_of_sequences The number of sequences.
 * \param aligned_sequence The sequence that is edited.
 * \param alignment_length The length of the alignment.
 * \param edit_seq_file File that saves the edited_sequences of the internal nodes.
 * \param offset Saves the size of the edited_sequences.
 * \param alignment_file File where the alignment is saved.
 * 
 */
void
edit2alignment(FILE *sequence_file,		//sequence file
			   long *seq_positions,		//sequence positions
			   FILE *edit_file,			//file saving the edit profiles
			   long *edit_positions,	//array saving the correspondence between edit profile and position in edit_file
			   int node_number,			//the current node
			   int number_of_sequences,	//number of sequences
			   char *aligned_sequence,	//the sequence that is edited
			   int alignment_length,	//length of the alignment - and thus of aligned_sequence
			   FILE *edit_seq_file,		//file saving the edited_sequences of the internal nodes
			   int offset,				//saves the size of the edited_sequence
			   FILE* alignment_file)	//file saving the alignments
{
	fseek(edit_file, edit_positions[node_number-number_of_sequences], SEEK_SET);
	const LINE_LENGTH = 50;
	char line[LINE_LENGTH];
	fgets(line, LINE_LENGTH , edit_file);
	int child1 = atoi(line);
	fgets(line, LINE_LENGTH , edit_file);
	int child2 = atoi(line);
	fgets(line, LINE_LENGTH , edit_file);
	int is_leaf1 = atoi(line);
	fgets(line, LINE_LENGTH , edit_file);
	int is_leaf2 = atoi(line);
	
	static char seq_line[10];

	char x;
	int number;
	int pos = 0;
	
	//first child
	while(fgets(line, LINE_LENGTH , edit_file)!=NULL)
	{
		x = line[0];
		if (x == '*')
			break;
		number = atoi(&line[1]);
		if (x == 'M')
		{
			while (number > 0)
			{
				if (aligned_sequence[pos] == 'X')
					--number;
				++pos;
			}
		}
		else if (x == 'I')
		{
			while (number > 0)
			{
				if (aligned_sequence[pos] == 'X')
					--number;
				++pos;
			}
		}
		else if (x == 'D')
		{
			while (number > 0)
			{
				if (aligned_sequence[pos] == 'X')
				{
					aligned_sequence[pos] = '-';
					--number;
				}
				++pos;
			}
		}
	}

	if (is_leaf1)
	{
		edit_seq2aligned_seq(aligned_sequence, sequence_file, seq_positions[child1], alignment_file);
	}
	else
	{
		fprintf(edit_seq_file, "%s", aligned_sequence);
		edit2alignment(sequence_file, seq_positions, edit_file, edit_positions, child1, number_of_sequences, aligned_sequence, alignment_length, edit_seq_file, offset, alignment_file);
	}
	
	//second child
	fseek(edit_seq_file, offset, SEEK_CUR);
	fgets(aligned_sequence, alignment_length+3, edit_seq_file);
	fseek(edit_seq_file, offset, SEEK_CUR);
	
	pos = 0;
	fseek(edit_file, edit_positions[node_number-number_of_sequences], SEEK_SET);
	while(fgets(line, LINE_LENGTH , edit_file)!=NULL)
	{
		x = line[0];
		if (x == '*')
			break;
		number = atoi(&line[1]);
		if (x == 'M')
		{
			while (number > 0)
			{
				if (aligned_sequence[pos] == 'X')
					--number;
				++pos;
			}
		}
		else if (x == 'I')
		{
			while (number > 0)
			{
				if (aligned_sequence[pos] == 'X')
				{
					aligned_sequence[pos] = '-';
					--number;
				}
				++pos;
			}
		}
		else if (x == 'D')
		{
			while (number > 0)
			{
				if (aligned_sequence[pos] == 'X')
					--number;
				++pos;
			}
		}
	}

	if (is_leaf2)
	{
		edit_seq2aligned_seq(aligned_sequence, sequence_file, seq_positions[child2], alignment_file);
	}
	else
	{
		fprintf(edit_seq_file, "%s", aligned_sequence);
		edit2alignment(sequence_file, seq_positions, edit_file, edit_positions, child2, number_of_sequences, aligned_sequence, alignment_length, edit_seq_file, offset, alignment_file);
	}
}




//  * The file has the follwing format (# and text behind are only comments and not included into the file):
// 		 * 1		# Number of profile.
// 		 * 1		# is DNA or not.
// 		 * 5		# Number of columns in the profile.
// 		 * 4A 1C	# In this column are 4 'A' and 1 'C'
// 		 * 3G		# In this column are 3 'G'
// 		 * 5A		# In this column are 5 'A'
// 		 * 2A 3C	# In this column are 2 'A' and 3 'C'
// 		 * 5C		# In this column are 5 'C'
// 		 * *		# Marks the end of this profile



/**
 * \brief Writes a profile to a file.
 * 
 * \param sumup_prf The profile array, not a real profile.
 * \param length The length of the profile.
 * \param file The FILE object to write the the profile into.
 * \param is_dna The type of sequence.
 * \param number The number of the profile.
 */
void
write2file(int **sumup_prf,
		   int length, 
		   FILE *file,
		   int number,
		   Fastal_param *param_set)
{
	char *pos2aa = &(param_set->pos2char[0]);
	fprintf(file,"%i\n0\n%i\n1\n",number, length );
	int i, j;
	int alphabet_size = param_set->alphabet_size;
	
	i = 0;
	int x = 0;
	while (i < length)
	{
		for (j = 0; j < alphabet_size; ++j)
			if (sumup_prf[j][i] > 0)
			{
				if (x)
					fprintf(file," %c%i", pos2aa[j],sumup_prf[j][i]);
				else
					fprintf(file,"%c%i", pos2aa[j],sumup_prf[j][i]);
				x = 1;
			}
// 		x = 1;
		x = 0;
		fprintf(file,"\n");
		++i;
	}
	fprintf(file,"*\n");
}









/**
*	main of the fastal algorithm
*/
int
fastal(int argc,	//number of arguments
	   char **argv)	//arguments first = fastal, second = tree
{
	
	int test;
	for (test = 0; test < argc; ++test)
	{
		printf("%s\n",argv[test]);
	}
	
	struct fastal_arguments arguments;

	arguments.output_file = "out.aln";
	arguments.tree_file = NULL;
	arguments.gep = -1;
	arguments.gop = -10;
	arguments.method = "fast";
	
// 	argp_parse (&argp, argc, argv, 0, 0, &arguments);

	Fastal_param *param_set = vcalloc(1,sizeof(Fastal_param));
	fill_parameters(arguments.is_dna, param_set, arguments.method);
	param_set->gep = arguments.gep;
	param_set->gop = arguments.gop;


	int alphabet_size = param_set->alphabet_size;

	//sequence file management
	char **seq_name;
	long *file_positions = NULL;
	long **tmp = &file_positions;
	int number_of_sequences = make_index_of_file(arguments.sequence_file, tmp);
	FILE *seq_file = fopen(arguments.sequence_file,"r");


	//edit file management
	FILE *edit_file = fopen("edit_tmp","w+");
	long current_edit_pos;
	long *edit_positions = vcalloc(number_of_sequences,sizeof(long));


	//profile management
	Fastal_profile **profiles = vcalloc(3,sizeof(Fastal_profile*));
	initiate_profiles(profiles, param_set);
	FILE * prof_file = fopen("prf_tmp","w+");
	long* profile_positions = vcalloc(4,sizeof(long*));
	int max_prof = 4;
	int saved_prof = 0;
	
	
	//dynamic programming matrix
	double ** dyn_matrix = vcalloc(1,sizeof(double*));
	dyn_matrix[0] = vcalloc(1,sizeof(double));
	int *length1 = vcalloc(1,sizeof(int));
	int *length2 = vcalloc(1,sizeof(int));
	*length1 = 1;
	*length2 = 2;
	int i;
	int **sumup_prf = vcalloc(alphabet_size+1,sizeof(int*));
	for (i = 0; i < alphabet_size+1; ++i)
		sumup_prf[i] = vcalloc(1,sizeof(int));
	int *sumup_length = vcalloc(1,sizeof(int));
	*sumup_length = 1;



	if (arguments.tree_file == NULL)
	{
		arguments.tree_file = "HUMAN.tree";
		printf("CONSTRUCT TREE\n");
		make_partTree(arguments.sequence_file, arguments.tree_file, 4, 20);
	}


	printf("CONSTRUCT ALIGNMENT\n");
	//tree file management
	FILE *tree_file = fopen(arguments.tree_file,"r");
	const int LINE_LENGTH = 100;
	char line[LINE_LENGTH];
	char delims[] = " ";
	int node[3];
	char *result = NULL;
	int j;
	int alignment_length;


	//memory for sparse dynamic
	int *diagonals = vcalloc(3,sizeof(int));
	int *dig_length = vcalloc(1,sizeof(int));
	*dig_length = 3;
	int **list = NULL;//vcalloc(1,sizeof(int*));
// 	list[0] = vcalloc(7,sizeof(int));
	int *list_length = vcalloc(1,sizeof(int));

	*list_length = 0;
	int ***list_p = vcalloc(1,sizeof(int**));



	//bottom-up traversal
	while(fgets(line, LINE_LENGTH, tree_file)!=NULL)
	{
		//read profiles
		node[0] = atoi(strtok(line,delims));
		node[1] = atoi(strtok(NULL,delims));
		node[2] = atoi(strtok(NULL,delims));
		//getting profile of second child
		if (node[1] < number_of_sequences)
		{
			file_pos2profile(seq_file, file_positions[node[1]], profiles[1], node[1], param_set);	//profile to save into
		}
		else
		{
			profile_file2profile(profiles[1], prof_file, profile_positions[--saved_prof], param_set);
			fseek (prof_file , profile_positions[saved_prof] , SEEK_SET);
		}

		//getting profile of first child
		if (node[0] < number_of_sequences)
		{
			file_pos2profile(seq_file, file_positions[node[0]], profiles[0], node[0], param_set);	//profile to save into
		}
		else
		{
			profile_file2profile(profiles[0], prof_file, profile_positions[--saved_prof], param_set);
			fseek (prof_file , profile_positions[saved_prof] , SEEK_SET);
		}
		if (saved_prof == max_prof)
		{
			max_prof += 5;
			profile_positions = vrealloc(profile_positions, max_prof*sizeof(long));
		}
		edit_positions[node[2]-number_of_sequences] = ftell(edit_file);
		profile_positions[saved_prof] = ftell(prof_file);
		++saved_prof;
		if (!strcmp(param_set->method,"nw"))
		{
			dyn_matrix = resize_dyn_matrix(dyn_matrix, length1, length2, profiles[0]->length+1, profiles[1]->length+1);
			alignment_length = prf_nw(profiles[0], profiles[1], dyn_matrix, edit_file, sumup_prf, sumup_length, param_set);
			write2file(sumup_prf, alignment_length, prof_file, node[2], param_set);
		}
		else if (!strcmp(param_set->method, "fast"))
		{
			char *file_name1 = vtmpnam(NULL);
			char *file_name2 = vtmpnam(NULL);
			char *seq1 = profile2consensus(profiles[0], file_name1, param_set);
			char *seq2 = profile2consensus(profiles[1], file_name2, param_set);
			int **diagonals_p = &diagonals;
			int num_diagonals = seq_pair2blast_diagonal(file_name1, file_name2, diagonals_p, dig_length, strlen(seq1),strlen(seq2), arguments.is_dna);
			diagonals = diagonals_p[0];
			char *p = &param_set->pos2char[0];
			list = diagonals2int(diagonals, num_diagonals, seq1, seq2, list_length, param_set);//, profiles[0], profiles[1], p);
			alignment_length = list2linked_pair_wise_fastal(profiles[0], profiles[1], param_set, list, *list_length, edit_file, prof_file, node[2]);
			int x;

			for (x = 0; x < *list_length; ++x)
			{
				vfree(list[x]);
			}
			vfree(list);
			list = NULL;
			vfree(seq1);
			vfree(seq2);
		}
	}

	//free_memory & close files
	vfree(diagonals);
	fclose(tree_file);
	fclose(prof_file);
	free_fastal_profile(profiles[0], alphabet_size);
	free_fastal_profile(profiles[1], alphabet_size);
	vfree(profiles);
	vfree(profile_positions);
	free_dyn_matrix(*length1,dyn_matrix);
	for (i = 0; i <= alphabet_size; ++i)
	{
		vfree(sumup_prf[i]);
	}
	vfree(sumup_prf);
	vfree(param_set);

	//bottom-down traversal (edit_files --> alignment)
	char file_name[FILENAMELEN];
	sprintf(file_name,arguments.output_file);

	FILE *alignment_file = fopen(file_name, "w");
	FILE *edit_seq_file = fopen("edit_seq.tmp","w+");

	char *aligned_sequence = vcalloc(alignment_length+3, sizeof(char));

	
	long offset = ftell(edit_seq_file);
	for (i = 0; i < alignment_length; ++i)
	{
		fprintf(edit_seq_file, "X");
		aligned_sequence[i] = 'X';
	}
	aligned_sequence[i]= '\n';
	aligned_sequence[i+1]= '\0';
	fprintf(edit_seq_file, "\n");
	offset = (ftell(edit_seq_file) - offset)*-1;


	edit2alignment(seq_file, file_positions, edit_file, edit_positions, node[2], number_of_sequences, aligned_sequence, alignment_length, edit_seq_file, offset, alignment_file);


	//free_memory & close files

	vfree(edit_positions);
	fclose(edit_file);
	fclose(seq_file);

	return 0;
}




//******************   toolbox   ***************************


/**
*	enlargement of the dynamic programming matrix in case it is to small.
*/
double**
resize_dyn_matrix(double **dyn_matrix,	//the dynamic programming matrix
				  int *old_length1,		//old length of dimension 1
				  int *old_length2,		//old length of dimension 2
				  int length1,			//new minimum length of dimension 1
				  int length2)			//new maximum length of dimension 2
{
	int i;
	if (*old_length1 < length1)
	{
		dyn_matrix = vrealloc(dyn_matrix,length1*sizeof(double*));
		for (i = *old_length1; i < length1; ++i)
			dyn_matrix[i] = vcalloc(*old_length2,sizeof(double));
		*old_length1 = length1;
	}
	if (*old_length2 < length2)
	{
		for (i = 0;i<*old_length1; ++i)
			dyn_matrix[i] = vrealloc(dyn_matrix[i], length2*sizeof(double));
		*old_length2 = length2;
	}
	return dyn_matrix;
}



/**
*	frees the memory of a dynamic programming matrix
*/
void
free_dyn_matrix(int length1,			//length of first dimension
				double **dyn_matrix)	//dynamic matrix
{
	int i = 0;
	for (; i<length1; ++i)
		vfree(dyn_matrix[i]);
	vfree(dyn_matrix);
}



/**
*	initialises the profiles.
*/
void
initiate_profiles(Fastal_profile **profiles,	//profiles pointer
				  Fastal_param *param_set)
{
	int alphabet_size = param_set->alphabet_size;
	int i,j;
	for (i =0; i < 3; ++i)
	{
		profiles[i] = vcalloc(1,sizeof(Fastal_profile));
		profiles[i]->weight = 1;
		profiles[i]->is_leaf = 1;
		profiles[i]->prf = vcalloc(alphabet_size, sizeof(int*));
		for (j = 0; j < alphabet_size; ++j)
		{
			profiles[i]->prf[j] = vcalloc(PROFILE_ENLARGEMENT, sizeof(int));
		}
		profiles[i]->allocated_memory = PROFILE_ENLARGEMENT;
	}
}


/**
*	initalises the files where the profiles are temporarly stored.
*/
void
initiate_profile_files(FILE **profile_files)
{
	char names[10];
	int i = 0;
	for (;i < 4; ++i)
	{
		sprintf(names,"tmp_prf_%i",i);
		profile_files[i] = fopen(names,"w+");
	}
}



/**
 *	frees all memory occupied by the profile
 */
void
free_fastal_profile(Fastal_profile* profile, int alphabet_size)
{
	--alphabet_size;
	for (;alphabet_size >= 0; --alphabet_size)
		vfree(profile->prf[alphabet_size]);
	vfree(profile->prf);
}


/**
*	initialize the parameters
*/
void
fill_parameters(int is_dna, Fastal_param *param_set, char *method)
{
	sprintf(param_set->method,"%s",method);
	int i;
	if (is_dna)
	{
		param_set->alphabet_size = 10;
		char tmp1[] = {'A','C','G','T','N','R','Y','D','M','W'};
		int  tmp2[] = { 0, -1,  1,  7, -1, -1, 2, -1, -1, -1, -1, -1, 8, 4, -1, -1, -1, 5, -1, 3, -1, -1, 9, -1, 6, -1};
		for (i = 0; i<param_set->alphabet_size; ++i)
			param_set->pos2char[i] = tmp1[i];
		for (i = 0; i<26; ++i)
			param_set->char2pos[i] = tmp2[i];
		param_set->M = read_matrice("dna_idmat");
	}
	else
	{
		param_set->alphabet_size = 24;
		char tmp1[] = {'A','C','G','T','F','D','H','I','K','L','M','N','P','Q','R','S','E','V','W','Y','B','J','X','Z'};
		int tmp2[] = { 0, 20,  1,  5, 16,  4,  2,  6,  7, 21,  8,  9,  10, 11, -1, 12, 13, 14, 15, 3, -1, 17,  18,22, 19,23};
		for (i = 0; i<param_set->alphabet_size; ++i)
			param_set->pos2char[i] = tmp1[i];
		for (i = 0; i<26; ++i)
			param_set->char2pos[i] = tmp2[i];
		param_set->M = read_matrice("blosum62mt");
	}
}
/*********************************COPYRIGHT NOTICE**********************************/
/*© Centro de Regulacio Genomica */
/*and */
/*Cedric Notredame */
/*Tue Oct 27 10:12:26 WEST 2009. */
/*All rights reserved.*/
/*This file is part of T-COFFEE.*/
/**/
/*    T-COFFEE is free software; you can redistribute it and/or modify*/
/*    it under the terms of the GNU General Public License as published by*/
/*    the Free Software Foundation; either version 2 of the License, or*/
/*    (at your option) any later version.*/
/**/
/*    T-COFFEE is distributed in the hope that it will be useful,*/
/*    but WITHOUT ANY WARRANTY; without even the implied warranty of*/
/*    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the*/
/*    GNU General Public License for more details.*/
/**/
/*    You should have received a copy of the GNU General Public License*/
/*    along with Foobar; if not, write to the Free Software*/
/*    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA*/
/*...............................................                                                                                      |*/
/*  If you need some more information*/
/*  cedric.notredame@europe.com*/
/*...............................................                                                                                                                                     |*/
/**/
/**/
/*	*/
/*********************************COPYRIGHT NOTICE**********************************/