#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <ctype.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"


// #include <omp.h>
// #define CHUNKSIZE 100
// #define N     1000


//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
 */



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


void
fill_arguments_sparse(Sparse_dynamic_param* method_arguments_p)
{
	method_arguments_p->diagonals = vcalloc(3,sizeof(Diagonal));
	method_arguments_p->dig_length = vcalloc(1,sizeof(int));
	*method_arguments_p->dig_length = 3;
	method_arguments_p->list = NULL;
	method_arguments_p->list_length = vcalloc(1,sizeof(int));
	*method_arguments_p->list_length = 0;
 	method_arguments_p->file_name1 = vtmpnam(NULL);
 	method_arguments_p->file_name2 = vtmpnam(NULL);
}

void
free_sparse(Sparse_dynamic_param* method_arguments_p)
{
	vfree(method_arguments_p->diagonals);
	vfree(method_arguments_p->dig_length);
	vfree(method_arguments_p->list_length);
}


/**
 * \brief One run of sparse dynamic programming.
 *
 * \param profiles The profiles.
 * \param param_set The fastal parameters.
 * \param method_arguments_p The method arguments.
 * \param is_dna Sequences are DNA (\a is_dna = 1) or protein.
 * \param edit_file The edit file.
 * \param prof_file the profile file.
 * \param number Number of the parent node.
 * \return The length of the alignment.
 */
int
sparse_dyn(Fastal_profile **profiles,
		   Fastal_param *param_set,
		   void *method_arguments_p,
		   int is_dna,
		   FILE *edit_file,
		   FILE *prof_file,
		   int number)
{
// 	printf("WHAT THE HELL ARE YOU DOING HERE?\n");
	Sparse_dynamic_param *arguments = (Sparse_dynamic_param*)method_arguments_p;
// 	static char *file_name1 = vtmpnam(NULL);
// 	static char *file_name2 = vtmpnam(NULL);
	char *file_name1 = arguments->file_name1;
	char *file_name2 = arguments->file_name2;
	char *seq1, *seq2;
	Fastal_profile *tmp1 = profiles[0];
	Fastal_profile *tmp2 = profiles[1];

	seq1 = profile2consensus(tmp1, param_set);
	seq2 = profile2consensus(tmp2, param_set);


	int **diagonals_p = &(arguments->diagonals);
	int num_diagonals = -1;
	if (!strcmp(param_set->diag_method, "blastz"))
	{
		FILE *cons_f = fopen(file_name1,"w");
		fprintf(cons_f, ">%i\n", tmp1->prf_number);
		fprintf(cons_f, "%s", seq1);
		fprintf( cons_f, "\n");
		fclose(cons_f);
		cons_f = fopen(file_name2,"w");
		fprintf(cons_f, ">%i\n", tmp2->prf_number);
		fprintf(cons_f, "%s", seq2);
		fprintf( cons_f, "\n");
		fclose(cons_f);
		num_diagonals = seq_pair2blastz_diagonal(file_name1, file_name2, diagonals_p, arguments->dig_length, strlen(seq1),strlen(seq2), is_dna);
	}
	else if (!strcmp(param_set->diag_method, "blast"))
	{
		FILE *cons_f = fopen(file_name1,"w");
		fprintf(cons_f, ">%i\n", tmp1->prf_number);
		fprintf(cons_f, "%s", seq1);
		fprintf( cons_f, "\n");
		fclose(cons_f);
		cons_f = fopen(file_name2,"w");
		fprintf(cons_f, ">%i\n", tmp2->prf_number);
		fprintf(cons_f, "%s", seq2);
		fprintf( cons_f, "\n");
		fclose(cons_f);
		int l1 = strlen(seq1);
		int l2 = strlen(seq2);
		num_diagonals = seq_pair2blast_diagonal(file_name1, file_name2, diagonals_p, arguments->dig_length, l1, l2, is_dna);
// 		int *num_p = &num_diagonals;
// 		Segment* seg = extend_diagonals(*diagonals_p, num_p, l1, l2);
// 		printf("A: %i\n", num_diagonals);
	}
	else if (!strcmp(param_set->diag_method, "ktup"))
	{
		num_diagonals = seq_pair2diagonal_own(seq1, seq2, diagonals_p, arguments->dig_length, strlen(seq1),strlen(seq2), is_dna, 3);
//	num_diagonals = seq_pair2diagonal_swift(seq1, seq2, diagonals_p, arguments->dig_length, strlen(seq1),strlen(seq2), is_dna, 3);
	}


// 	arguments->diagonals = diagonals_p[0];
// 	arguments->list = segments2int(seg, *num_p, seq1, seq2, profiles[0], profiles[1], arguments->list_length, param_set);

// t ** segments2int_gap(Segment *diagonals, int num_diagonals, char *seq1, char *seq2, Fastal_profile *profile1, Fastal_profile *profile2, int *num_points, Fastal_param *param_set);

//  arguments->list = diagonals2int_dot(arguments->diagonals, num_diagonals, seq1, seq2, profiles[0], profiles[1], arguments->list_length, param_set);
// 	arguments->list = diagonals2int_euclidf(arguments->diagonals, num_diagonals, seq1, seq2, profiles[0], profiles[1], arguments->list_length, param_set);
	arguments->list = diagonals2int_gap_test(arguments->diagonals, num_diagonals, seq1, seq2, profiles[0], profiles[1], arguments->list_length, param_set);
// 	arguments->list = diagonals2int(arguments->diagonals, num_diagonals, seq1, seq2, arguments->list_length, param_set);
	int alignment_length = list2linked_pair_wise_fastal(profiles[0], profiles[1], param_set, arguments->list, *arguments->list_length, edit_file, prof_file, number);
	int x;

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


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


/**
 * \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 num_points Number of points in the list
 * \param param_set Fastal parameters.
 * \return A 2-dim array which contains all points needed for the sparse dynamic programming algorithm.
 */
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 dig_length;
	if (seq1 > seq2)
		dig_length = l1;
	else
		dig_length = l2;

	int current_size = num_diagonals*dig_length + l1 +l2;

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


	//+1 because diagonals start here at position 1, like in "real" dynamic programming
	int a = -1, 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;
	}

	vfree(diags);
	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 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;
	vfree(diagx);
	vfree(diagy);


	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 num_points Number of points in the list
 * \param param_set Fastal parameters.
 * \return A 2-dim array which contains all points needed for the sparse dynamic programming algorithm.
 */
int **
diagonals2int_gap_test(int *diagonals, int num_diagonals, char *seq1, char *seq2, Fastal_profile *profile1, Fastal_profile *profile2, int *num_points, Fastal_param *param_set)
{
	int alphabet_size = param_set->alphabet_size;
	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 = -1, 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;
	}

	vfree(diags);
	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];




//SIMPLEGAP
			int num_seq = profile1->number_of_sequences + profile2->number_of_sequences;
			double gap_num = 0;
			int char_c;
			for (char_c = 0; char_c < alphabet_size; ++char_c)
			{

				gap_num += profile1->prf[char_c][diagx[i]-1] + profile2->prf[char_c][diagy[i]-1];
			}

			gap_num /= num_seq;

			list[list_pos][2] = M[toupper(seq1[diagx[i]-1])-'A'][toupper(seq2[diagy[i]-1])-'A'] * gap_num;

// CLUSTAL
// 			int num_seq = profile1->number_of_sequences + profile2->number_of_sequences;
// 			double gap_num = 0;
// 			int char_c, char_c2;
// 			for (char_c = 0; char_c < alphabet_size; ++char_c)
// 				for (char_c2 = 0; char_c2 < alphabet_size; ++char_c2)
// 			{
// 				gap_num += (profile1->prf[char_c][diagx[i]-1]/profile1->number_of_sequences) * (profile2->prf[char_c][diagy[i]-1]/profile2->number_of_sequences) * M[param_set->pos2char[char_c]-'A'][param_set->pos2char[char_c2]-'A'];
// 			}
// 			list[list_pos][2] = gap_num;


			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;
			int num_seq = profile1->number_of_sequences + profile2->number_of_sequences;
			double gap_num = 0;
			int char_c;
			for (char_c = 0; char_c < alphabet_size; ++char_c)
			{
				gap_num += profile1->prf[char_c][l1-1] + profile2->prf[char_c][i-1];
			}

			gap_num /= num_seq;

			list[list_pos][2] = M[toupper(seq1[l1-1])-'A'][toupper(seq2[i-1])-'A'] * gap_num;
		}
		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;
	vfree(diagx);
	vfree(diagy);
	vfree(old_pos);

	return list;
}


int **
diagonals2int_euclidf(int *diagonals, int num_diagonals, char *seq1, char *seq2, Fastal_profile *profile1, Fastal_profile *profile2, int *num_points, Fastal_param *param_set)
{
	int alphabet_size = param_set->alphabet_size;
	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 = -1, 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;
	}

	vfree(diags);
	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];
			int char_c;
			double tmp_score = 0;
			double freq1, freq2;
			for (char_c = 0; char_c < alphabet_size; ++char_c)
			{
				freq1 = (double)profile1->prf[char_c][diagx[i]-1] / profile1->number_of_sequences;

				freq2 = (double)profile2->prf[char_c][diagy[i]-1] / profile2->number_of_sequences;

				tmp_score += ( freq1 - freq2) * (freq1 - freq2);
			}

			list[list_pos][2] = 10 - sqrt(tmp_score);

			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;
			int char_c;
			int tmp_score = 0;
			double freq1, freq2;
			for (char_c = 0; char_c < alphabet_size; ++char_c)
			{
				freq1 = profile1->prf[char_c][l1-1] / profile1->number_of_sequences;
				freq2 = profile2->prf[char_c][i-1] / profile2->number_of_sequences;
				tmp_score += ( freq1 - freq2) * (freq2 - freq1);
			}
			list[list_pos][2] = 10 - sqrt(tmp_score);
// 			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;
	vfree(diagx);
	vfree(diagy);
	vfree(old_pos);

	return list;
}

int **
diagonals2int_dot(int *diagonals, int num_diagonals, char *seq1, char *seq2, Fastal_profile *profile1, Fastal_profile *profile2, int *num_points, Fastal_param *param_set)
{
	int alphabet_size = param_set->alphabet_size;
	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 = -1, 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;
	}

	vfree(diags);
	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];
			int char_c;
			double tmp_score = 0;
			double freq1, freq2;
			for (char_c = 0; char_c < alphabet_size; ++char_c)
			{
				freq1 = (double)profile1->prf[char_c][diagx[i]-1] / profile1->number_of_sequences;

				freq2 = (double)profile2->prf[char_c][diagy[i]-1] / profile2->number_of_sequences;

				tmp_score += freq1 * freq2;
			}

			list[list_pos][2] = tmp_score * 10;

			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;
			int char_c;
			int tmp_score = 0;
			double freq1, freq2;
			for (char_c = 0; char_c < alphabet_size; ++char_c)
			{
				freq1 = profile1->prf[char_c][l1-1] / profile1->number_of_sequences;
				freq2 = profile2->prf[char_c][i-1] / profile2->number_of_sequences;
				tmp_score += freq2 * freq1;
			}
			list[list_pos][2] = tmp_score * 10;
// 			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;
	vfree(diagx);
	vfree(diagy);
	vfree(old_pos);

	return list;
}


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, i, j, LEN=0, start_trace = -1;
	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;
	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);

		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%i\n", node_number,LEN, prf1->number_of_sequences+prf2->number_of_sequences);

	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 Turns 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, 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 = -1, 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]);
	}
	return seq;
}


int
diag_compare (const void * a, const void * b)
{
	return (((Diagonal_counter*)b)->count - ((Diagonal_counter*)a)->count);
}

/**
 * \brief Calculates the diagonals between two sequences.
 *
 * Uses  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_pair2diagonal_own(char *seq1,
					  char *seq2,
					  int **diagonals,
					  int *dig_length,
					  int l1,
					  int l2,
					  int is_dna,
					  int word_length)
{
// 	word_length = 7;
	int word_number, i;
	int ng;
	if (is_dna)
	{
		word_number = (int)pow(5, word_length);
		ng = 4;
	}
	else
	{
		word_number = (int)pow(24, word_length);
		ng = 24;
	}
	int **word_index = vcalloc(word_number, sizeof(int*));
	for (i = 0 ; i < word_number; ++i)
	{
		word_index[i] = vcalloc(20, sizeof(int));
		word_index[i][0] = 2;
		word_index[i][1] = 20;
	}


	//making of k-tup index of seq1

	int *prod=vcalloc (word_length, sizeof(int));
	for ( i=0; i<word_length; i++)
	{
		prod[word_length-i-1]=(int)pow(ng,i);
	}

	int aa[256];
	if (is_dna)
	{
		aa['A'] = 0;
		aa['C'] = 1;
		aa['G'] = 2;
		aa['T'] = 3;
		aa['U'] = 3;
	}
	else
	{
		aa['A'] = 0;
 		aa['B'] = 20;
		aa['C'] = 1;
		aa['D'] = 2;
		aa['E'] = 3;
		aa['F'] = 4;
		aa['G'] = 5;
		aa['H'] = 6;
		aa['I'] = 7;
		aa['J'] = 20;
		aa['K'] = 8;
		aa['L'] = 9;
		aa['M'] = 10;
		aa['N'] = 11;
		aa['P'] = 12;
		aa['Q'] = 13;
		aa['R'] = 14;
		aa['S'] = 15;
		aa['T'] = 16;
		aa['V'] = 17;
		aa['W'] = 18;
 		aa['X'] = 20;
		aa['Y'] = 19;
		aa['X'] = 20;
	}
	int index = 0;
	for (i = 0; i < word_length; ++i)
	{
		index += aa[(short)seq1[i]] *prod[i];
	}
	word_index[index][2] = 0;
	word_index[index][0] = 3;
	int z = -1;
	int *tmp;
	for (; i < l1; ++i)
	{
		index -= aa[(short)seq1[++z]] * prod[0];
		index *= ng;
		index += aa[(short)seq1[i]];
		tmp = word_index[index];
		if (tmp[0] == tmp[1])
		{
			tmp[1] += 25;
			tmp = vrealloc(tmp, word_index[index][1] *sizeof(int));
			word_index[index] = tmp;
		}
		tmp[tmp[0]++] = i;
	}



	//counting diagonals
	const int window_length = 14;

	Diagonal_counter *diag_index = vcalloc(l1+l2, sizeof(Diagonal_counter));
	int num = l1+l2;
	for (i = 0; i < num; ++i)
	{
		diag_index[i].diagonal = i;
		diag_index[i].count = 0;
	}
	index = 0;

	int j;
	for (i = 0; i < word_length; ++i)
	{
		index += aa[(short)seq2[i]] *prod[i];
		for (j = 2; j < word_index[index][0]; ++j)
		{
			++(diag_index[i - word_index[index][j] + l1].count);
		}
	}

	z = -1;
	int i2 = i-1;
	int second_index = index;
	for (; i < window_length; ++i)
	{
		index -= aa[(short)seq2[++z]] * prod[0];
		index *= ng;
		index += aa[(short)seq2[i]];
		tmp = word_index[index];
		for (j = 2; j < tmp[0]; ++j)
		{
			++(diag_index[i - tmp[j] + l1].count);
		}
	}
	int z2 = -1;
	for (; i < l2; ++i)
	{
		index -= aa[(short)seq2[++z]] * prod[0];
		index *= ng;
		index += aa[(short)seq2[i]];
		second_index -= aa[(short)seq2[++z2]] * prod[0];
		second_index *= ng;
		second_index += aa[(short)seq2[++i2]];

		tmp = word_index[index];
		for (j = 2; j < tmp[0]; ++j)
		{
			++(diag_index[i - tmp[j] + l1].count);
		}


		tmp = word_index[second_index];
		for (j = 2; j < tmp[0]; ++j)
		{
			if (diag_index[i2 - tmp[j] + l1].count > window_length-3)
				diag_index[i2 - tmp[j] + l1].count = window_length+100;
			else
				--diag_index[i2 - tmp[j] + l1].count;
		}
	}


	//choose diagonals
	int *diags = diagonals[0];
	int current_pos = 0;


	qsort (diag_index, num, sizeof(Diagonal_counter*), diag_compare);

	i = 0;
	int y, x;
	while (diag_index[i].count > window_length+10)
	{
		if (current_pos > (*dig_length)-3)
		{
			(*dig_length) += 30;
			diags = vrealloc(diags, sizeof(int)*(*dig_length));
		}


		y = diag_index[i].diagonal - l1;
		if (y < 0)
		{
			x = y * (-1);
			y = 0;
		}
		else
		{
			x = 0;
		}
		diags[current_pos++] = x;
		diags[current_pos++] = y;
		diags[current_pos++] = 200;
		++i;
	}

	vfree(diag_index);
	for (i = 0; i < word_number; ++i)
		vfree(word_index[i]);
	vfree(word_index);
	diagonals[0] = diags;
	return current_pos/3;
}



int
seq_pair2diagonal_swift(char *seq1,
						char *seq2,
						int **diagonals,
						int *dig_length,
						int l1,
						int l2,
						int is_dna,
						int word_length)
{
	int word_number, i;
	int ng;
	if (is_dna)
	{
		word_number = (int)pow(5, word_length);
		ng = 5;
	}
	else
	{
		word_number = (int)pow(24, word_length);
		ng = 24;
	}
	int **word_index = vcalloc(word_number, sizeof(int*));
	for (i = 0 ; i < word_number; ++i)
	{
		word_index[i] = vcalloc(20, sizeof(int));
		word_index[i][0] = 2;
		word_index[i][1] = 20;
	}


	//making of k-tup index of seq1

	int *prod=vcalloc (word_length, sizeof(int));
	for ( i=0; i<word_length; i++)
	{
		prod[word_length-i-1]=(int)pow(ng,i);
	}

	int aa[256];
	aa['A'] = 0;
	aa['C'] = 1;
	aa['G'] = 2;
	aa['T'] = 3;
	aa['U'] = 4;
	int index = 0;
	for (i = 0; i < word_length; ++i)
	{
		index += aa[(short)seq1[i]] *prod[i];
	}
	word_index[index][2] = 0;
	word_index[index][0] = 3;
	int z = -1;
	int *tmp;
	for (; i < l1; ++i)
	{
		index -= aa[(short)seq1[++z]] * prod[0];
		index *= ng;
		index += aa[(short)seq1[i]];
		tmp = word_index[index];
		if (tmp[0] == tmp[1])
		{
			tmp[1] += 25;
			tmp = vrealloc(tmp, word_index[index][1] *sizeof(int));
			word_index[index] = tmp;
		}
		tmp[tmp[0]++] = i;
	}


	//counting diagonals
	const int window_length = 14;
	const int threshold = 12;

	Swift_diagonal *diag_index = vcalloc(l1+l2, sizeof(Swift_diagonal));
	int num = l1+l2;
	for (i = 0; i < num; ++i)
	{
		diag_index[i].diagonal = i;
		diag_index[i].count = 0;
		diag_index[i].start = -99999;
		diag_index[i].end = -99999;
	}

	index = 0;

	int j;
	for (i = 0; i < word_length; ++i)
	{
		index += aa[(short)seq2[i]] *prod[i];
		for (j = 2; j < word_index[index][0]; ++j)
		{
			++(diag_index[i - word_index[index][j] + l1].count);
		}
	}

	z = -1;
	int tmp_index;
	for (; i < l2; ++i)
	{
		index -= aa[(short)seq2[++z]] * prod[0];
		index *= ng;
		index += aa[(short)seq2[i]];
		tmp = word_index[index];
		for (j = 2; j < tmp[0]; ++j)
		{
			tmp_index = i - tmp[j] + l1;
			if (i - diag_index[tmp_index].end > window_length)
			{
				if (diag_index[tmp_index].count < threshold)
				{
					diag_index[tmp_index].count = 0;
					diag_index[tmp_index].start = i;
					diag_index[tmp_index].end = i + word_length;
				}

			}
			else
			{
				++(diag_index[tmp_index].count);
			}
		}

	}



	// choose diagonals
	int *diags = diagonals[0];
	int current_pos = 0;
	int x, y;
	for (i = 0; i < num; ++i)
	{
		if (diag_index[i].count > threshold)
		{
			if (current_pos > (*dig_length)-3)
			{
				(*dig_length) += 30;
				diags = vrealloc(diags, sizeof(int)*(*dig_length));
			}
			y = diag_index[i].diagonal - l1;
			if (y < 0)
			{
				x = y * (-1);
				y = 0;
			}
			else
			{
				x = 0;
			}
			diags[current_pos++] = x;
			diags[current_pos++] = y;
			diags[current_pos++] = 200;
		}
	}

	vfree(diag_index);
	for (i = 0; i < word_number; ++i)
		vfree(word_index[i]);
	vfree(word_index);
	diagonals[0] = diags;

	return current_pos/3;
}



/**
 * \brief Calculates the diagonals between two sequences.
 *
 * Uses a k-tup index to choose 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)
{
// 	static int blast_measure[12]={0,0,0,0,0,0,0,0,0,0,0,0};
	int *diag = vcalloc(l1 + l2, sizeof(int));
	char *out_file = vtmpnam(NULL);
	char blast_command[200];
// 	char blast_command2[200];
// 	if (x)
// 	{
// 		int i;
// 		printf("BLAST-Types:\n");
// 		for (i = 0; i < 11; ++i)
// 		{
// 			printf("Type %i: %i\n", i, blast_measure[i]);
// 		}
// 		return 0;
// 	}
// 	char blast_command2[600];
// 	sprintf(blast_command2, "less %s", out_file);

	if (is_dna)
	{
		sprintf(blast_command, "bl2seq -p blastn -i %s -j %s -D 1 -g F -o %s -S 1 -F F", 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 -F F -S 1", 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";
	int length, pos_q, pos_d;
	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-20)
		{
			(*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;
		}

	}
	fclose(diag_f);
	int round = 0;
	int e_threshold = 10;
	while ((current_pos == 0) && (round < 10))
	{
		if (is_dna)
		{
			sprintf(blast_command, "bl2seq -p blastn -i %s -j %s -D 1 -g F -o %s -S 1 -F F -W 6 -e %i", seq_file_name1, seq_file_name2, out_file, e_threshold);
		}
		else
		{
		  sprintf(blast_command, "bl2seq -p blastp -i %s -j %s -D 1 -g F -o %s -F F -S 1 -e %i", seq_file_name1, seq_file_name2, out_file, e_threshold);
		}
		system(blast_command);
		e_threshold *= 10;
		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";
		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-20)
			{
				(*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;
			}
		}
		fclose(diag_f);
		++round;
		if (current_pos < 27)
			current_pos = 0;
	}
// 	++blast_measure[round];

	if (current_pos == 0)
	{
		printf("BLAST NOT SUCCESFULL\n");
		if (l1 < l2)
		{
			int i;
			int diff = l2 - l1 + 10;
			for (i = diff; i > 0; --i)
			{
				if (current_pos >= *dig_length-20)
				{
					(*dig_length) += 90;
					diags = vrealloc(diags, sizeof(int)*(*dig_length));
				}
				diags[current_pos++] = 0;
				diags[current_pos++] = i;
				diags[current_pos++] = 100;
			}
			diff = 10;
// 			printf("A: %i\n", diff);
			for (i = 0; i < diff; ++i)
			{
				if (current_pos >= *dig_length-20)
				{
					(*dig_length) += 90;
					diags = vrealloc(diags, sizeof(int)*(*dig_length));
				}
				diags[current_pos++] = i;
				diags[current_pos++] = 0;
				diags[current_pos++] = 100;
			}
		}
		else
		{
			int i;
			int diff = 10;
// 			printf("A: %i\n", diff);
			for (i = diff; i > 0; --i)
			{
				if (current_pos >= *dig_length-20)
				{
					(*dig_length) += 90;
					diags = vrealloc(diags, sizeof(int)*(*dig_length));
				}
				diags[current_pos++] = 0;
				diags[current_pos++] = i;
				diags[current_pos++] = 100;
			}
			diff = l1 - l2 + 10;
			for (i = 0; i < diff; ++i)
			{
				if (current_pos >= *dig_length-20)
				{
					(*dig_length) += 90;
					diags = vrealloc(diags, sizeof(int)*(*dig_length));
				}
				diags[current_pos++] = i;
				diags[current_pos++] = 0;
				diags[current_pos++] = 100;
			}

		}
	}


	vfree(diag);

	diagonals[0] = diags;
	return current_pos/3;
}



int
seq_pair2blat_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[200];
// 	char blast_command2[200];
// 	char blast_command2[600];
// 	sprintf(blast_command2, "less %s", out_file);

	if (is_dna)
	{
		sprintf(blast_command, "blat %s %s %s -out=blast8 -q=dna -t=dna -maxGap=0 >/dev/null 2>/dev/null", seq_file_name2, seq_file_name1, out_file);
	}
	else
	{
 		sprintf(blast_command, "blat %s %s %s -out=blast8 -prot -maxGap=0 >/dev/null 2>/dev/null", seq_file_name2, seq_file_name1, 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";
	int length, pos_q, pos_d;
	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-20)
		{
			(*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;
		}
	}
	if (current_pos == 0)
	{
		printf("BLAT NOT SUCCESFULL\n");
		if (l1 < l2)
		{
			int i;
			int diff = l2 - l1 + 10;
			for (i = diff; i > 0; --i)
			{
				if (current_pos >= *dig_length-20)
				{
					(*dig_length) += 90;
					diags = vrealloc(diags, sizeof(int)*(*dig_length));
				}
				diags[current_pos++] = 0;
				diags[current_pos++] = i;
				diags[current_pos++] = 100;
			}
			diff = 10;
// 			printf("A: %i\n", diff);
			for (i = 0; i < diff; ++i)
			{
				if (current_pos >= *dig_length-20)
				{
					(*dig_length) += 90;
					diags = vrealloc(diags, sizeof(int)*(*dig_length));
				}
				diags[current_pos++] = i;
				diags[current_pos++] = 0;
				diags[current_pos++] = 100;
			}
		}
		else
		{
			int i;
			int diff = 10;
// 			printf("A: %i\n", diff);
			for (i = diff; i > 0; --i)
			{
				if (current_pos >= *dig_length-20)
				{
					(*dig_length) += 90;
					diags = vrealloc(diags, sizeof(int)*(*dig_length));
				}
				diags[current_pos++] = 0;
				diags[current_pos++] = i;
				diags[current_pos++] = 100;
			}
			diff = l1 - l2 + 10;
			for (i = 0; i < diff; ++i)
			{
				if (current_pos >= *dig_length-20)
				{
					(*dig_length) += 90;
					diags = vrealloc(diags, sizeof(int)*(*dig_length));
				}
				diags[current_pos++] = i;
				diags[current_pos++] = 0;
				diags[current_pos++] = 100;
			}

		}
	}

// 	printf("END\n");
	vfree(diag);
	fclose(diag_f);
	diagonals[0] = diags;
	return current_pos/3;
}



/**
 * \brief Calculates the diagonals between two sequences.
 *
 * Uses blastz 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_pair2blastz_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[200];
// 	char blast_command2[200];
// 	char blast_command2[600];
// 	sprintf(blast_command2, "less %s", out_file);

	if (is_dna)
	{
		sprintf(blast_command, "~/Download/blastz-source/blastz %s %s B=0 K=10000> %s", seq_file_name1, seq_file_name2, out_file);
	}
	else
	{
		printf("SORRY - no BLASTZ with amino accid\n");
		exit(0);
	}
	system(blast_command);

	int *diags = diagonals[0];
	FILE *diag_f = fopen(out_file,"r");
	char line[300];
	char delims[] = " ";
// 	char *result = NULL;
	int length, pos_q, pos_d;
	int current_pos = 0;
	while (fgets(line, 300, diag_f) != NULL)
	{
		if (line[0] == 'a')
		{
			char *line_tmp;
			while (fgets(line, 300, diag_f) != NULL)
			{
				if (line[0] == '}')
					break;

				if (line[2] == 'l')
				{
					line_tmp = &line[4];
					if (current_pos >= *dig_length-20)
					{
						(*dig_length) += 90;
						diags = vrealloc(diags, sizeof(int)*(*dig_length));
					}
					pos_q = atoi(strtok(line_tmp, delims));
					pos_d = atoi(strtok(NULL, delims));
					length = atoi(strtok(NULL, delims) - pos_q);
					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;
					}
				}
			}
		}
	}


	if (current_pos == 0)
	{
		printf("BLASTZ NOT SUCCESFULL\n");
		if (l1 < l2)
		{
			int i;
			int diff = l2 - l1 + 10;

			for (i = diff; i > 0; --i)
			{
				if (current_pos >= *dig_length-20)
				{
					(*dig_length) += 90;
					diags = vrealloc(diags, sizeof(int)*(*dig_length));
				}
				diags[current_pos++] = 0;
				diags[current_pos++] = i;
				diags[current_pos++] = 100;
			}
			diff = 10;

			for (i = 0; i < diff; ++i)
			{
				if (current_pos >= *dig_length-20)
				{
					(*dig_length) += 90;
					diags = vrealloc(diags, sizeof(int)*(*dig_length));
				}
				diags[current_pos++] = i;
				diags[current_pos++] = 0;
				diags[current_pos++] = 100;
			}
		}
		else
		{
			int i;
			int diff = 10;

			for (i = diff; i > 0; --i)
			{
				if (current_pos >= *dig_length-20)
				{
					(*dig_length) += 90;
					diags = vrealloc(diags, sizeof(int)*(*dig_length));
				}
				diags[current_pos++] = 0;
				diags[current_pos++] = i;
				diags[current_pos++] = 100;
			}
			diff = l1 - l2 + 10;
			for (i = 0; i < diff; ++i)
			{
				if (current_pos >= *dig_length-20)
				{
					(*dig_length) += 90;
					diags = vrealloc(diags, sizeof(int)*(*dig_length));
				}
				diags[current_pos++] = i;
				diags[current_pos++] = 0;
				diags[current_pos++] = 100;
			}

		}
	}

	vfree(diag);
	fclose(diag_f);
	diagonals[0] = diags;
	return current_pos/3;
}



//**************************   needleman-wunsch aligning **********************************************************


void
fill_arguments_nw(Nw_param* method_arguments_p, int alphabet_size)
{
	method_arguments_p-> dyn_matrix = vcalloc(1,sizeof(double*));
	method_arguments_p->dyn_matrix[0] = vcalloc(1,sizeof(double));
	method_arguments_p->length1 = vcalloc(1,sizeof(int));
	method_arguments_p->length2 = vcalloc(1,sizeof(int));
	*method_arguments_p->length1 = 1;
	*method_arguments_p->length2 = 1;
	method_arguments_p->sumup_prf = vcalloc(alphabet_size+1,sizeof(int*));
	int i;
	for (i = 0; i < alphabet_size+1; ++i)
		method_arguments_p->sumup_prf[i] = vcalloc(1,sizeof(int));
	method_arguments_p->sumup_length = vcalloc(1,sizeof(int));
	*method_arguments_p->sumup_length = 1;
}


void
free_nw(Nw_param* method_arguments_p, int alphabet_size)
{
	free_dyn_matrix(*method_arguments_p->length1,method_arguments_p->dyn_matrix);
	int i;
	for (i = 0; i <= alphabet_size; ++i)
	{
		vfree(method_arguments_p->sumup_prf[i]);
	}
	vfree(method_arguments_p->sumup_prf);
	vfree(method_arguments_p->length1);
	vfree(method_arguments_p->length2);
	vfree(method_arguments_p->sumup_length);
}


/**
 * \brief One run of needleman-wunsch dynamic programming.
 *
 * \param profiles The profiles.
 * \param param_set The fastal parameters.
 * \param method_arguments_p The method arguments.
 * \param is_dna Sequences are DNA (\a is_dna = 1) or protein.
 * \param edit_file The edit file.
 * \param prof_file the profile file.
 * \param number Number of the parent node.
 * \return The length of the alignment.
 */
int
nw_dyn(Fastal_profile **profiles, Fastal_param *param_set, void *method_arguments_p, int is_dna, FILE *edit_file, FILE *prof_file, int number)
{
	Nw_param *arguments = (Nw_param*)method_arguments_p;
// 	int old_length1 = *arguments->length1;
// 	int old_length2 = *arguments->length2;
	arguments->dyn_matrix = resize_dyn_matrix(arguments->dyn_matrix, *arguments->length1, *arguments->length2, profiles[0]->length+1, profiles[1]->length+1);
	*arguments->length1 = profiles[0]->length+1;
	*arguments->length2 = profiles[1]->length+1;
	int alignment_length = prf_nw(profiles[0], profiles[1], arguments->dyn_matrix, edit_file, arguments->sumup_prf, arguments->sumup_length, param_set);
	write2file(arguments->sumup_prf, alignment_length, prof_file, number,profiles[0]->number_of_sequences + profiles[1]->number_of_sequences, param_set);
	return alignment_length;
}


/**
 * \brief This method takes a profile and turns it into a sumed up version.
 *
 * Required for NW-algorithm.
 * \param profile The profile to sum up.
 * \param sumup A field where the result will be stored.
 * \param param_set Parameters for the fastal algorithm.
 * \return The new \a sumup.
 */
int**
sumup_profile(Fastal_profile *profile,
			  int **sumup,
			  Fastal_param *param_set)
{

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


/**
 * \brief Turns the dynamic programming matrix into a editfile and calculates the new profile.
 *
 * Required for NW-algorithm.
 * \param prog_matrix The dynamic programming matrix.
 * \param prf1 The first profile.
 * \param prf2 The second profile.
 * \param edit_f A File object (already opened) to write the edit sequence into.
 * \param prf_field A 2-dim array to save the new profile into.
 * \param field_length Length of the new profile.
 * \param param_set Parameters of the Fastal-Algorithm
 */
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,
	   Fastal_profile *profile2,
	   double **prog_matrix,
	   FILE *edit_file_name,
	   int **sumup_prf,
	   int *sumup_length,
	   Fastal_param *param_set)
{
	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);
}


/************** GOTOH ***********************/


void
fill_arguments_gotoh(Gotoh_param* method_arguments_p, int alphabet_size)
{
	method_arguments_p->m_matrix = vcalloc(1,sizeof(double*));
	method_arguments_p->m_matrix[0] = vcalloc(1,sizeof(double));
	method_arguments_p->d_matrix = vcalloc(1,sizeof(double*));
	method_arguments_p->d_matrix[0] = vcalloc(1,sizeof(double));
	method_arguments_p->i_matrix = vcalloc(1,sizeof(double*));
	method_arguments_p->i_matrix[0] = vcalloc(1,sizeof(double));
	method_arguments_p->length1 = vcalloc(1,sizeof(int));
	method_arguments_p->length2 = vcalloc(1,sizeof(int));
	method_arguments_p->log_saver = vcalloc(alphabet_size+1, sizeof(double*));
	*method_arguments_p->length1 = 1;
	*method_arguments_p->length2 = 1;
	method_arguments_p->sumup_prf = vcalloc(alphabet_size+1,sizeof(int*));
	int i;
	for (i = 0; i < alphabet_size+1; ++i)
	{
		method_arguments_p->sumup_prf[i] = vcalloc(1,sizeof(int));
		method_arguments_p->log_saver[i] = vcalloc(1, sizeof(double));
	}
	method_arguments_p->sumup_length = vcalloc(1,sizeof(int));
	*method_arguments_p->sumup_length = 1;
}


void
free_gotoh(Gotoh_param* method_arguments_p, int alphabet_size)
{
	free_dyn_matrix(*method_arguments_p->length1,method_arguments_p->m_matrix);
	free_dyn_matrix(*method_arguments_p->length1,method_arguments_p->i_matrix);
	free_dyn_matrix(*method_arguments_p->length1,method_arguments_p->d_matrix);

	int i;
	for (i = 0; i <= alphabet_size; ++i)
	{
		vfree(method_arguments_p->sumup_prf[i]);
	}
	vfree(method_arguments_p->sumup_prf);
	vfree(method_arguments_p->length1);
	vfree(method_arguments_p->length2);
	vfree(method_arguments_p->sumup_length);
}


int
gotoh_dyn(Fastal_profile **profiles, Fastal_param *param_set, void *method_arguments_p, int is_dna, FILE *edit_file, FILE *prof_file, int number)
{
	Gotoh_param *arguments = (Gotoh_param*)method_arguments_p;
	arguments->m_matrix = resize_dyn_matrix(arguments->m_matrix, *arguments->length1, *arguments->length2, profiles[0]->length+1, profiles[1]->length+1);
	arguments->i_matrix = resize_dyn_matrix(arguments->i_matrix, *arguments->length1, *arguments->length2, profiles[0]->length+1, profiles[1]->length+1);
	arguments->d_matrix = resize_dyn_matrix(arguments->d_matrix, *arguments->length1, *arguments->length2, profiles[0]->length+1, profiles[1]->length+1);
	int i;
	if (profiles[1]->length > *arguments->length2-1)
	{
		for (i = 0; i < param_set->alphabet_size; ++i)
		{
			arguments->log_saver[i] = vrealloc(arguments->log_saver[i], (profiles[1]->length)*sizeof(double*));
		}
	}
	*arguments->length1 = profiles[0]->length+1;
	*arguments->length2 = profiles[1]->length+1;
	int alignment_length = prf_gotoh(profiles[0], profiles[1], edit_file, arguments, param_set);
	write2file(arguments->sumup_prf, alignment_length, prof_file, number, profiles[0]->number_of_sequences + profiles[1]->number_of_sequences, param_set);
	return alignment_length;
}


int
gotoh_matrix2edit_file(double **m_matrix,		//dynamic programming matrix
					   double **v_matrix,		//dynamic programming matrix
					   double **h_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
{
	double comp_num = log((double)prf1->number_of_sequences) + log((double)prf2->number_of_sequences);
	int** M = param_set->M;
	int alphabet_size = param_set->alphabet_size;
	double gep = param_set -> gep;
	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;
	int current_mode = 0;
	//determine start mode
	char *pos2char = param_set->pos2char;
	if (h_matrix[n][m] == m_matrix[n][m])
	{
		current_mode = 2;
	}
	else
	{

		if (v_matrix[n][m] == m_matrix[n][m])
		{
			current_mode = 1;
		}
		else
		{
			current_mode = 0;
		}
	}
// 	printf("%f %f %f - %i\n",h_matrix[n][m],v_matrix[n][m],m_matrix[n][m], current_mode);
	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 (current_mode == 2)
		{
			for (i = 0; i<alphabet_size; ++i)
			{
				prf_field[i][field_pos] = prf2->prf[i][m];
			}
			if (h_matrix[n][m] != (h_matrix[n][m+1]+gep))
			{
				current_mode = 0;
			}
			++m;
			++ field_pos;
			if (last != 2)
			{
				fprintf(edit_f,"%c%i\n",sumc[last],sum[last]);
				sum[last] = 0;
			}
			last = 2;
			++sum[last];
		}
		else
		{
			if (current_mode== 1)
			{
				for (i = 0; i<alphabet_size; ++i)
				{
					prf_field[i][field_pos] = prf1->prf[i][n];
				}
				if (v_matrix[n][m] != (v_matrix[n+1][m]+gep))
				{
					current_mode = 0;
				}
				++n;
				++ field_pos;

				if (last != 1)
				{
					fprintf(edit_f,"%c%i\n",sumc[last],sum[last]);
					sum[last] = 0;
				}
				last = 1;
				++sum[last];
			}
			else
			{
				double match_score = 0.0;
				int char_c, char_c2;
				for (char_c = 0; char_c < alphabet_size; ++char_c)
				{
					for (char_c2 = 0; char_c2 < alphabet_size; ++char_c2)
					{

						if ((log(prf1->prf[char_c][n]) != -1) && ( log(prf2->prf[char_c2][m]) != -1))
						{
							match_score += exp(log((double)prf1->prf[char_c][n]) + log((double)prf2->prf[char_c2][m])-comp_num) * M[pos2char[char_c]-'A'][pos2char[char_c2]-'A'];
						}
					}
				}
				if (m_matrix[n+1][m+1] + match_score != m_matrix[n][m])
				{
					if (m_matrix[n][m] == v_matrix[n][m])
					{
						current_mode = 1;
						continue;
					}
					if (m_matrix[n][m] == h_matrix[n][m])
					{
						current_mode = 2;
						continue;
					}
				}
				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]);

	int needed = MAX(prf1_length -n, prf2_length -m);

	if ((*field_length) - needed -10 < field_pos)
	{
		(*field_length) += needed +10;

		for (i = 0; i <alphabet_size+1; ++i)
		{
			prf_field[i] = vrealloc(prf_field[i], (*field_length)*sizeof(int));
		}
	}
	//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 The gotoh dynamic programming algorithm.
 *
 * \param profile1 The first profile.
 */
int
prf_gotoh(Fastal_profile *profile1,
		  Fastal_profile *profile2,
		  FILE *edit_file_name,
		  Gotoh_param *arguments,
		  Fastal_param *param_set)
{

// printf("I AM HERE - again\n");
	int **sumup_prf   = arguments->sumup_prf;
	int *sumup_length = arguments->sumup_length;
	int alphabet_size = param_set->alphabet_size;
	double gop = param_set->gop;
	double gep = param_set->gep;

	const int INF = -999999;
	int i;

	double **m_matrix = arguments->m_matrix;
	double **h_matrix = arguments->i_matrix;
	double **v_matrix = arguments->d_matrix;

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

	int** M = param_set->M;
	double match_score;
	for (i = prof2_length; i >= 0; --i)
	{
		m_matrix[prof1_length][i] = gop + gep * (prof2_length-i);
		v_matrix[prof1_length][i] = INF;
		h_matrix[prof1_length][i] = m_matrix[prof1_length][i];
	}

	m_matrix[prof1_length][prof2_length] = 0.0;
	h_matrix[prof1_length][prof2_length] = INF;
	v_matrix[prof1_length][prof2_length] = INF;
	int l;
	double comp_num = log((double)profile1->number_of_sequences) + log((double)profile2->number_of_sequences);
	static double *log_test = NULL;
	if (!log_test)
		log_test = vcalloc(alphabet_size, sizeof(double));
// 	int k;
	int **prf1 = profile1->prf;
	int **prf2 = profile2->prf;
	double **log_test2 = arguments->log_saver;
	for (l = 0; l < alphabet_size; ++l)
	{
		for (i = 0; i < profile2->length; ++i)
		{
			if (prf2[l][i] > 0)
			{
				log_test2[l][i] = log(prf2[l][i]);
			}
			else
				log_test2[l][i] = -1;
		}
	}

	char *pos2char = param_set->pos2char;
	i = prof1_length-1;
	while (i >=0)
	{
		j = prof2_length-1;

		for (l = 0; l < alphabet_size; ++l)
		{
			if (prf1[l][i] > 0)
				log_test[l] = log((double)prf1[l][i]);
			else
				log_test[l] = -1;
		}
		m_matrix[i][prof2_length] = gop + gep *(prof1_length-i);
		v_matrix[i][prof2_length] = m_matrix[i][prof2_length];
		h_matrix[i][prof2_length] = INF;
		while (j >=0)
		{

			match_score = 0.0;
			v_matrix[i][j] = (MAX(v_matrix[i+1][j], m_matrix[i+1][j] + gop) + gep);
			h_matrix[i][j] = (MAX(h_matrix[i][j+1], m_matrix[i][j+1] + gop) + gep);

			int char_c, char_c2;
			int num = 0;
			for (char_c = 0; char_c < alphabet_size; ++char_c)
			{
				for (char_c2 = 0; char_c2 < alphabet_size; ++char_c2)
				{

					if ((log_test[char_c] != -1) && (log_test2[char_c2][j] != -1))
					{
						match_score += exp(log_test[char_c] + log_test2[char_c2][j]-comp_num) * M[pos2char[char_c]-'A'][pos2char[char_c2]-'A'];
					}
				}
			}

			m_matrix[i][j] = m_matrix[i+1][j+1]+match_score;

			if (m_matrix[i][j] < v_matrix[i][j])
			{
				m_matrix[i][j] = v_matrix[i][j];
			}
			if (m_matrix[i][j] < h_matrix[i][j])
			{
				m_matrix[i][j] = h_matrix[i][j];
			}

			--j;
		}
		--i;
	}
	return gotoh_matrix2edit_file(m_matrix, v_matrix, h_matrix, profile1, profile2, edit_file_name, sumup_prf, sumup_length, param_set);
}


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

/**
 * \brief Writes the alignment into the profile file and the edit file.
 *
 * \param profiles The two profiles to combine.
 * \param alignment The alinment information.
 * \param alignment The length of the alignment.
 * \param edit_f The edit file.
 * \param prof_f The profile file.
 * \param node_number the new node number.
 */
void
alignment2files(Fastal_profile **profiles,
				Fastal_param *param_set,
				int **alignment,
				int alignment_length,
				FILE *edit_f,
				FILE *prof_f,
				int node_number)
{
	fprintf(edit_f, "%i\n%i\n%i\n%i\n",profiles[0]->prf_number, profiles[1]->prf_number, profiles[0]->is_leaf, profiles[1]->is_leaf);
	fprintf(prof_f, "%i\n0\n%i\n1\n", node_number, alignment_length);

	int **prf1 = profiles[0]->prf;
	int **prf2 = profiles[1]->prf;
	int i = 0;
	int pos = 0;
	int pos1, pos2;

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

	while (i < alignment_length)
	{

		pos1 = alignment[0][pos];
		pos2 = alignment[1][pos];
		// match
		if ((pos1 != -1) && (pos2 != -1))
		{

			combine_profiles2file(prf1, prf2, pos1, pos2, param_set, prof_f, 'M');
			if (state != 0)
			{
				fprintf(edit_f, "%c%i\n",statec[state], num);
				num =1;
				state = 0;
			}
			else
				++num;
			++i;
		}
		// insertion in seq 1
		else if (pos1 != -1)
		{
			combine_profiles2file(prf1, prf2, pos1, pos2, param_set, prof_f, 'I');
			if (state != 2)
			{
				fprintf(edit_f, "%c%i\n",statec[state], num);
				num =1;
				state = 2;
			}
			else
				++num;
			++i;
		}
		// deletion in seq 1
		else if (pos2 != -1)
		{
			combine_profiles2file(prf1, prf2, pos1, pos2, param_set, prof_f, 'D');
			if (state != 1)
			{
				fprintf(edit_f, "%c%i\n",statec[state], num);
				num =1;
				state = 1;
			}
			else
				++num;
			++i;
		}
		++pos;
	}
	fprintf(edit_f, "%c%i\n",statec[state], num);

	fprintf(edit_f,"*\n");
	fprintf(prof_f,"*\n");

}


//******************************* 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;
	profile->number_of_sequences = 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, x;

	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;
			x = 0;
			while ((line[i] != '\n') && (line[i] != '\0'))
			{
				if (line[i] != '-')
				{
					for(j = 0; j<alphabet_size; ++j )
						profile->prf[j][seq_length+x] = 0;
					profile->prf[aa2pos[toupper((short)line[i])]][seq_length+x] = 1;
					++x;
				}
				++i;
			}
			seq_length += x;

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

}



/**
 * \brief Constructs index of fasta_file.
 *
 * The index is of length n (n= number of sequences in the given multi fasta file.). In the order of appearance in the file the position of each sequence in the file is stored.
 * \param file_name The file with the sequences.
 * \param file_positions Array to save the positions in.
 * \return The number of sequences in \a file_name.
 */
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));

	FILE *file = fopen(file_name,"r");

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



/**
 * \brief Reads a profile from a profile file.
 *
 * \param prof A Fastal_profile object to save the profile in.
 * \param profile_f file where the profile is stored.
 * \param position Position of the profile in \a profile_f.
 * \param param_set The parameter set for Fastal
 */

void
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->is_leaf = atoi(line);

	fgets(line, LINE_LENGTH, profile_f);
	prof->length = atoi(line);
	fgets(line, LINE_LENGTH, profile_f);
	prof->weight = atoi(line);
	fgets(line, LINE_LENGTH, profile_f);
	prof->number_of_sequences = 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));
		}
	prof->allocated_memory = prof->length;
	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[(short)result[0]]][i] = atoi(result_num);
			result = strtok( NULL, delims );
		}
	}
}



/**
 * \brief Writes a profile into a file
 *
 * \param profile Pointer to the profile which has to be saved.
 * \param file A File object (already opened) to write the profile to.
 * \param param_set The parameters for the fastal algorithm.
 */
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->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 (NOT IN USE)
*/
// 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);
// 	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[(short)result[0]]][i] = atoi(result_num);
// 			result = strtok( NULL, delims );
// 		}
// 	}
// }






/**
 * \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') && (line[i] != '\0'))
			{
				while (aligned_sequence[pos] == '-')
				{

					fprintf(alignment_file,"-");
					++pos;
				}
				if (line[i] != '-')
				{

					fprintf(alignment_file,"%c",line[i]);
					++pos;
				}
				++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 int 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];
// 	printf("SO EINE VERDAMMTE SCHEISE ABER AUCH\n");
	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)
	{
// 		printf("LEAF\n");
		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 leaf.
// 	* 5		# Number of columns in the profile.
// 	* 4A 1C	# In the first column are 4 'A' and 1 'C'
// 	* 3G	# In the second column are 3 'G'
// 	* 5A	# In the third column are 5 'A'
// 	* 2A 3C	# In the fourth column are 2 'A' and 3 'C'
// 	* 5C	# In the fifth 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. The format can be seen in ./test.txt
 * \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,
		   int num_sequences,
		   Fastal_param *param_set)
{
	char *pos2aa = &(param_set->pos2char[0]);
	fprintf(file,"%i\n0\n%i\n1\n%i\n",number, length, num_sequences );
	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 = 0;
		fprintf(file,"\n");
		++i;
	}
	fprintf(file,"*\n");
}



//*************************************   main function of the fasta algorithm ***********************************************

/**
*	\brief main of the fastal algorithm
*/
int
fastal_main(int argc,		//number of arguments
			char **argv)	//arguments first = fastal, second = tree
{

	int i;
	//pointer to arguments
	void * method_arguments_p;
	int (*alignment_function)(Fastal_profile **profiles, Fastal_param *param_set, void *method_arguments_p, int is_dna, FILE *edit_file, FILE *prof_file, int number);

	struct Fastal_arguments arguments;

	arg_parse (argc, argv, &arguments);

	Fastal_param *param_set = vcalloc(1,sizeof(Fastal_param));

	fill_parameters(arguments.is_dna, param_set, arguments.method, arguments.diag_method, arguments.mat);
	param_set->gep = arguments.gep;
	param_set->gop = arguments.gop;

// 	printf("%s",arguments.mat);
	if (arguments.evaluate)
	{
		printf("Calculate Sum of pairs Score\n");
		printf("Score: %f\n", calculate_sum_of_pairs_score_affine(arguments.sequence_file, param_set->M, param_set->gop, param_set->gep));
		vfree(param_set);
		exit(0);
	}


	if (arguments.agreement_score)
	{
		complete_agreement_score(arguments.aln2test, arguments.aln_ref);
		return 0;
	}


	if (arguments.num_ref_aln)
	{
		compute_ref_alignments(arguments.sequence_file, arguments.aln_ref, arguments.num_ref_aln, arguments.num_seq_in_ref);
		return 0;
	}



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



	//edit file management

// 	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(vtmpnam(NULL),"w+");
	long* profile_positions = vcalloc(4,sizeof(long*));
	int max_prof = 4;
	int saved_prof = 0;


	printf("METHOD: %s\n",param_set->method);
	if (strcmp(param_set->method, "fast") == 0)
	{
		method_arguments_p = vcalloc(1,sizeof(Sparse_dynamic_param));
		fill_arguments_sparse((Sparse_dynamic_param*)method_arguments_p);
		alignment_function = sparse_dyn;
	}
	else if (strcmp(param_set->method, "nw") == 0)
	{
		method_arguments_p = vcalloc(1,sizeof(Nw_param));
		fill_arguments_nw((Nw_param*)method_arguments_p, alphabet_size);
		alignment_function = nw_dyn;
	}
	else if (strcmp(param_set->method, "gotoh") == 0)
	{
		method_arguments_p = vcalloc(1,sizeof(Gotoh_param));
		fill_arguments_gotoh((Gotoh_param*)method_arguments_p, alphabet_size);
		alignment_function = gotoh_dyn;
	}
	else if (strcmp(param_set->method, "udisc") == 0)
	{
		method_arguments_p = vcalloc(1,sizeof(Udisc_param));
		fill_arguments_gotoh((Gotoh_param*)method_arguments_p, alphabet_size);
		alignment_function = gotoh_dyn;

	}
	else
	{
		printf("ERROR - METHOD");
		exit(1);
	}


	if (arguments.gap_iterate)
	{
		iterate(param_set, method_arguments_p, arguments.sequence_file, arguments.output_file, arguments.gap_iterate);
		return 0;
	}

	if (arguments.tree_file == NULL)
	{
		arguments.tree_file = vtmpnam(NULL);
		printf("CONSTRUCT TREE\n");
		if (strcmp(arguments.tree_method, "parttree")==0)
		{
			make_partTree(arguments.sequence_file, arguments.tree_file, arguments.tree_param1, arguments.tree_param2, arguments.is_dna, 0);
		}
		else if (strcmp(arguments.tree_method, "oligotree") == 0)
		{
			compute_oligomer_distance_tree(arguments.sequence_file, param_set->char2pos, arguments.tree_file, arguments.tree_param2, arguments.tree_param1, param_set->alphabet_size);
		}

		if (arguments.tree_only == 1)
			return 0;
	}


	if (arguments.tree_out == 1)
	{
		char tree_out_file_name[500];
		sprintf(tree_out_file_name, "%s.tree",arguments.output_file);
		char const LINE_LENGTH = 50;
		char line[LINE_LENGTH];

		FILE* in = fopen(arguments.tree_file, "r");
		FILE* out = fopen(tree_out_file_name, "w");
		while( (fgets(line, LINE_LENGTH, in)) != NULL)
			fprintf(out, "%s", line);
		fclose(in);
		fclose(out);
	}





	FILE *seq_file = fopen(arguments.sequence_file,"r");
// 	FILE *edit_file = fopen(vtmpnam(NULL),"w+");
	FILE *edit_file = fopen("aha","w+");

	printf("CONSTRUCT ALIGNMENT\n");
	FILE *tree_file = fopen(arguments.tree_file,"r");
	const int LINE_LENGTH = 100;
	char line[LINE_LENGTH];
	char delims[] = " ";
	int node[3];
	int alignment_length = -1;
	node[2] = -1;


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

		//aligning the sequences
		alignment_length = alignment_function(profiles, param_set, method_arguments_p, arguments.is_dna, edit_file, prof_file, node[2]);
	}


//	bottom-down traversal (edit_files --> alignment)
//	tmp_out_file_name = vtmpnam(NULL);

// 	FILE *alignment_file = fopen(tmp_out_file_name, "w");
	FILE *alignment_file = fopen(arguments.output_file, "w");
	FILE *edit_seq_file = fopen(vtmpnam(NULL),"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);
	fclose(alignment_file);
	fclose(tree_file);
	fclose(edit_file);
	fclose(seq_file);
	fclose(edit_seq_file);

		//set stuff for the next cycle
// 		arguments.sequence_file	= tmp_out_file_name;


// // 		//DEBUG
// // 		char copy_command[500];
// // 		sprintf(copy_command, "cp %s %s_%i", tmp_out_file_name, arguments.output_file, cycle);
// // 		system(copy_command);
// 		++cycle;
// 	}

// 	printf("HERE_COPY\n");
// 	char copy_command[2000];
// 	sprintf(copy_command, "mv %s %s", tmp_out_file_name, arguments.output_file);
// 	printf("%s\n", copy_command);
// 	int error = system(copy_command);
// 	printf("ERROR %i\n", error);


	//free_memory & close files
	fclose(prof_file);
	free_fastal_profile(profiles[0], alphabet_size);
	free_fastal_profile(profiles[1], alphabet_size);
	vfree(profiles);
	vfree(profile_positions);



// number_of_sequences


	if (arguments.score)
	{
		printf("Calculate Score\n");
		double aln_score = calculate_sum_of_pairs_score_affine(arguments.output_file, param_set->M, param_set->gop, param_set->gep);
		printf("SCORE: %f\n", aln_score);
	}





	if (!strcmp(param_set->method, "fast"))
	{
		free_sparse((Sparse_dynamic_param*)method_arguments_p);
	}
	else if (!strcmp(param_set->method, "nw"))
	{
		free_nw((Nw_param*)method_arguments_p, alphabet_size);
	}
	else if (!strcmp(param_set->method, "gotoh"))
	{
		free_gotoh((Gotoh_param*)method_arguments_p, alphabet_size);
	}

	vfree(param_set);

	//free_memory & close files

	vfree(edit_positions);


	return 0;
}




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


/**
 * \brief Enlargement of the dynamic programming matrix in case it is to small.
 *
 * \param dyn_matrix The dynamic programming matrix.
 * \param old_length1 Current size of dimension 1.
 * \param old_length2 Current size of dimension 2.
 * \param length1 New size of dimension 1.
 * \param length2 New size of dimension 2.
 * \return Pointer to the new array.
 */
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;
}



/**
 * \brief Frees the memory of a dynamic programming matrix.
 *
 * \param length1 Size of the first dimension of the matrix.
 * \param dyn_matrix The 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);
}



/**
 * \brief Initialises the profiles with basic values.
 *
 * \param profiles Array of 3 profiles.
 * \param param_set The fastal parameters
 */
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;
	}
}




/**
 * \brief frees all memory occupied by the profile
 *
 * \param profile The profile to free.
 * \param alphabet_size The alphabet_size.
 */
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);
}


/**
 * \brief Initialize the Fastal parameter set.
 *
 * \param is_dna 1 when sequences are dna, 0 when not
 * \param param_set The fastal parameter set.
 * \param method The method to use in Fastal.
*/
void
fill_parameters(int is_dna, Fastal_param *param_set, char *method, char *diag_method, char *mat)
{
	sprintf(param_set->method,"%s",method);
	sprintf(param_set->diag_method,"%s",diag_method);
	int i;
	printf("%s",mat);
	param_set->M = read_matrice(mat);
	if (is_dna)
	{
		param_set->alphabet_size = 4;
		char tmp1[] = {'A','C','G','T'};

// 		int  tmp2[] = { 0, 1,  1,  0, -1, -1, 2, 0, -1, -1, 3, -1, 1, 0, -1, -1, -1, 0, 1, 3, 4, -1, 3, -1, 1, -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->char2pos['A'] = 0;
		param_set->char2pos['B'] = 1;
		param_set->char2pos['C'] = 1;
		param_set->char2pos['D'] = 0;
		param_set->char2pos['G'] = 2;
		param_set->char2pos['H'] = 0;
		param_set->char2pos['K'] = 3;
		param_set->char2pos['M'] = 1;
		param_set->char2pos['N'] = 0;
		param_set->char2pos['R'] = 0;
		param_set->char2pos['S'] = 1;
		param_set->char2pos['T'] = 3;
		param_set->char2pos['U'] = 3;
		param_set->char2pos['W'] = 3;
		param_set->char2pos['Y'] = 1;
// 		param_set->M[0][3] = param_set->M[3][0] = -10;
// 		param_set->M[1][2] = param_set->M[2][1] = -10;
// 		param_set->M[0][1] = param_set->M[0][2] = param_set->M[1][0] = param_set->M[2][0] = -10;
// 		param_set->M[3][1] = param_set->M[3][2] = param_set->M[1][3] = param_set->M[2][3] = -10;
	}
	else
	{
		param_set->alphabet_size = 21;
		char tmp1[] = {'A','C','D','E','F','G','H','I','K','L','M','N','P','Q','R','S','T','V','W','Y','X'};
// 		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->char2pos['A'] = 0;
 		param_set->char2pos['B'] = 20;
		param_set->char2pos['C'] = 1;
		param_set->char2pos['D'] = 2;
		param_set->char2pos['E'] = 3;
		param_set->char2pos['F'] = 4;
		param_set->char2pos['G'] = 5;
		param_set->char2pos['H'] = 6;
		param_set->char2pos['I'] = 7;
		param_set->char2pos['J'] = 20;
		param_set->char2pos['K'] = 8;
		param_set->char2pos['L'] = 9;
		param_set->char2pos['M'] = 10;
		param_set->char2pos['N'] = 11;
		param_set->char2pos['P'] = 12;
		param_set->char2pos['Q'] = 13;
		param_set->char2pos['R'] = 14;
		param_set->char2pos['S'] = 15;
		param_set->char2pos['T'] = 16;
		param_set->char2pos['V'] = 17;
		param_set->char2pos['W'] = 18;
 		param_set->char2pos['X'] = 20;
		param_set->char2pos['Y'] = 19;
		param_set->char2pos['X'] = 20;
	}
}



/******************************COPYRIGHT NOTICE*******************************/
/*© Centro de Regulacio Genomica */
/*and */
/*Cedric Notredame */
/*Fri Feb 18 08:27:45 CET 2011 - Revision 596. */
/*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*******************************/