#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <ctype.h>
#include <string.h>
#include <sys/stat.h>
#include <dirent.h>


#include "io_lib_header.h"
#include "util_lib_header.h"
#include "fastal_lib_header.h"




/**
 * \brief Function to calculate the sum of pairs score with quasi affine gap costs.
 *
 * \param alignment_file File with the aligned sequences in fasta_format.
 * \param file_index Positions of the sequences in \a alignment_file.
 * \param number_of_sequences The number of sequences in \a alignment_file.
 * \param score_matrix Matrix containing the scores.
 * \param gop The gap opening costs.
 * \param gep The gep extension costs.
 * \return The score of the alignment.
 */
double
calculate_sum_of_pairs_score_affine(char *alignment_file_name,
									int **score_matrix,
									double gop,
									double gep)
{

	const int LINE_LENGTH = 1000;
	int number_of_sequences = 0;
	FILE *alignment_file = fopen(alignment_file_name,"r");
	if (alignment_file == NULL)
	{
		printf("FILE COULD NOT BE OPENED!\n");
		exit(1);
	}
	int alignment_length = 0;
	char line[LINE_LENGTH+1];
	line[0] = '0';

	while (line[0] != '>')
	{
		fgets(line, LINE_LENGTH, alignment_file);
	}
	fgets(line, LINE_LENGTH, alignment_file);
	while (line[0] != '>')
	{
		alignment_length += strlen(line);
		fgets(line, LINE_LENGTH, alignment_file);
	}

	int alphabet_size = 28;
	int **counting = vcalloc(alphabet_size, sizeof(int*));

	int i;
	for (i = 0; i < alphabet_size; ++i)
	{
		counting[i] = vcalloc(alignment_length, sizeof(int));
	}

	alphabet_size -= 2;
	char c;
	int was_gap = 0;
	int pos_line, pos_alignment = -1;
	fseek (alignment_file, 0, SEEK_SET);
	pos_line = 0;

	int gap_line = 0;
	int current_size = 100;
	int gap_size = 0;
	int gap_pos = 0;

	int **gaps = vcalloc(100, sizeof(int*));

	while (fgets(line, LINE_LENGTH, alignment_file) != NULL)
	{
		if (line[0] != '>')
		{
			pos_line = -1;
			while (((c = line[++pos_line]) != '\n') && (c != '\0'))
			{
				if (isalpha(c))
				{
					++counting[toupper(c)-'A'][++pos_alignment];
					if (was_gap)
					{
						gaps[gap_line][gap_pos++] = pos_alignment-1;
					}
					was_gap = 0;
				}
				else
				{
					++counting[alphabet_size][++pos_alignment];
					if (!was_gap)
					{
						++counting[alphabet_size+1][pos_alignment];
						was_gap = 1;
						if (gap_pos >= gap_size-4)
						{
							gap_size += 50;
							gaps[gap_line] = vrealloc(gaps[gap_line], gap_size * sizeof(int));
							gaps[gap_line][gap_pos++] = pos_alignment;
						}
					}
				}
			}
		}
		else
		{
			if (number_of_sequences != 0)
			{
				if ((gap_pos % 2) != 0)
				{
					gaps[gap_line][gap_pos++] = pos_alignment-1;
				}
				gaps[gap_line][gap_pos++] = alignment_length+2;
				gaps[gap_line][gap_pos] = alignment_length+2;
				++gap_line;
			}
			if (current_size == gap_line)
			{
				current_size += 50;
				gaps = vrealloc(gaps, current_size * sizeof(int*));
			}
			gap_size = 2;
			gaps[gap_line] = vcalloc(gap_size, sizeof(int));
			++number_of_sequences;
			pos_alignment = -1;
			was_gap = 0;

			gap_pos = 0;
		}
	}

	gaps[gap_line][gap_pos++] = alignment_length+2;
	gaps[gap_line][gap_pos] = alignment_length+2;


	long double score = 0;
	int j, k, tmp2;
	int non_gap;
	for (i = 0; i < alignment_length; ++i)
	{
		for (j = 0; j < alphabet_size; ++j)
		{
			if ((tmp2 = counting[j][i]) >1)
			{
				while (tmp2 > 1)
				{
					score += score_matrix[j][j] * (--tmp2);
				}
			}
			for (k = j+1; k < alphabet_size; ++k)
			{
				score += score_matrix[j][k] * counting[j][i] * counting[k][i];
			}

		}


		non_gap = number_of_sequences - counting[alphabet_size][i];
		score += counting[alphabet_size][i] * non_gap  * gep;

	}

	++gap_line;
	alignment_length += 2;
	unsigned long gap_open = 0;

	int chunk = 500;

// 	#pragma omp parallel shared(gaps, chunk, alignment_length, gap_open) private(i, j)
	{
		int alignment_length2 = alignment_length;
// 		#pragma omp for schedule(dynamic,chunk) reduction(+:gap_open) nowait
		for (i = 0; i < gap_line; ++i)
		{
			int *gaps1 = gaps[i];
			for (j = i+1; j < gap_line; ++j)
			{
				int *gaps2 = gaps[j];
				int k = 0;
				int l = 0;
				while ((gaps1[k] != alignment_length2) && (gaps2[l] != alignment_length2))
				{
					if (gaps1[k+1] < gaps2[l])
					{
						++gap_open;
						k+=2;
						continue;
					}
					if (gaps2[l+1] < gaps1[k])
					{
						++gap_open;
						l+=2;
						continue;
					}
					if (gaps1[k] < gaps2[l])
					{
						if (gaps1[k+1] < gaps2[l+1])
						{
							++gap_open;
							k+=2;
							continue;
						}
						else
						{
							l+=2;
							continue;
						}
					}
					if (gaps1[k] > gaps2[l])
					{
						if (gaps1[k+1] > gaps2[l+1])
						{
							++gap_open;
							l+=2;
							continue;
						}
						else
						{
							k+=2;
							continue;
						}
					}
					if (gaps1[k] == gaps2[l])
					{
						if (gaps1[k+1] == gaps2[l+1])
						{
							k+=2;
							l+=2;
							continue;
						}
						else
						{
							if (gaps1[k+1] <gaps2[l+1])
							{
								k+=2;
							}
							else
							{
								l+=2;
							}
						}
					}

				}
				while (gaps1[k] != alignment_length2)
				{
					++gap_open;
					k+=2;
				}
				while (gaps2[l] != alignment_length2)
				{
					++gap_open;
					l+=2;
				}
			}
		}
	}

	score += gop * gap_open;


	return score / number_of_sequences;
}


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



int
get_random_value(int a, int b)
{
	return rand() % (b - a + 1) + a;
// 	j = 1 + (int)( 10.0 * rand() / ( RAND_MAX + 1.0 ) );
}





void
compute_ref_alignments(char *seq_file_name, char* ref_directory, int num_alignments, int num_seq_in_aln)
{

	printf("%s %s %i %i\n", seq_file_name, ref_directory, num_alignments, num_seq_in_aln);
	//check if target directory exists, if not make it
	char directory[500];
	if (ref_directory[0]== '~')
	{
		sprintf(directory,"%s%s",getenv("HOME"),ref_directory+1);
	}
	else
	{
		sprintf(directory,"%s",ref_directory);
	}
	struct stat st;
	if(stat(directory, &st) != 0)
	{
		mkdir(directory,0700);
	}

	printf("made\n");

	// get sequences & produce alignments
	//make index
	long *file_positions = NULL;
	long **tmp_pos = &file_positions;
	int number_of_sequences = make_index_of_file(seq_file_name, tmp_pos);
	char *tmp_file_name = vtmpnam(NULL);
	FILE *seq_file = fopen(seq_file_name, "r");
	int *already_taken = vcalloc(num_seq_in_aln, sizeof(int));
	int i, j, k, tmp, already, num_already;
	for (i = 0; i < num_alignments; ++i)
	{
		//choose randomly the sequences

		num_already = 0;
		while (num_already < num_seq_in_aln)
		{
			already = 0;
			tmp = get_random_value(0, number_of_sequences-1);
			for (k = 0; k < num_already; ++k)
			{
				if (already_taken[k] == tmp)
				{
					already = 1;
					break;
				}
			}
			if (!already)
			{
				already_taken[num_already] = tmp;
				++num_already;
			}
		}
		//write sequences into temporary file
		const int LINE_LENGTH = 200;
		char line[LINE_LENGTH];


		FILE* tmp_file = fopen(tmp_file_name, "w");
		for (k = 0; k < num_already; ++k)
		{
			fseek(seq_file, file_positions[already_taken[k]], SEEK_SET);
			fgets(line, LINE_LENGTH, seq_file);
			fprintf(tmp_file, "%s", line);
			while( (fgets(line, LINE_LENGTH, seq_file)) != NULL)
			{
				if (line[0] == '>')
					break;
				fprintf(tmp_file, "%s", line);
			}
		}
		fclose(tmp_file);

		//calculate alignment
		char aln_command[1000];
		sprintf(aln_command,"t_coffee -in %s -outfile %s/%i.fa -output fasta_aln 2>/dev/null >/dev/null", tmp_file_name, directory, i);
		system(aln_command);
	}

	vfree(already_taken);
}



// void
// compute agreement_score(char *alignment_file_name, char *alignment_directory)
// {
//
// }




/**
 * This algorithm choses according to a tree a number of sequenes of which it computes a reference algnment.
 * \param seq_file_name The sequence file.
 * \param tree_file_name The given tree.
 * \param ref_aln_name The file where the reference alignment should be written to.
 * \param num_seq_in_ref The number of sequences to choose.
**/
void
make_ref_alignment(char *seq_file_name, char *tree_file_name, char *ref_aln_name, int num_seq_in_ref)
{
	const int LINE_LENGTH = 200;
	char line[LINE_LENGTH];

	FILE *tree_f = fopen(tree_file_name,"r");


	long *file_positions = NULL;
	long **tmp = &file_positions;
	int num_sequences = make_index_of_file(seq_file_name, tmp);


	int every_x = num_sequences/(num_seq_in_ref);

	int *seq_ids = vcalloc(num_seq_in_ref,sizeof(int));
	char delims[] = " ";

	int pos = -1;
	fseek (tree_f, 0, SEEK_SET);
	int node;
	int dist = every_x;
	while(fgets(line, LINE_LENGTH, tree_f)!=NULL)
	{
		node = atoi(strtok(line,delims));
		if (node < num_sequences)
		{
			if (dist == every_x)
			{

				seq_ids[++pos] = node;
				dist = 0;
			}
			else
			{
				++dist;
			}
		}

		node = atoi(strtok(NULL,delims));
		if (node < num_sequences)
		{
			if (dist == every_x)
			{
				seq_ids[++pos] = node;
				dist = 0;
			}
			else
			{
				++dist;
			}
		}
	}
	fclose(tree_f);


	int i = 0;
	qsort(seq_ids, num_seq_in_ref, sizeof(int), compare);

	char *seq_ref_name = vtmpnam(NULL);
	FILE *seq_ref_f = fopen(seq_ref_name, "w");;
	FILE *seq_f = fopen(seq_file_name, "r");

	for (i = 0; i < num_seq_in_ref; ++i)
	{
		fseek (seq_f, file_positions[seq_ids[i]], SEEK_SET);
		fgets(line, LINE_LENGTH, seq_f);
		fprintf(seq_ref_f, "%s", line);
		while(fgets(line, LINE_LENGTH, seq_f)!=NULL)
		{
			if (line[0] == '>')
				break;
			fprintf(seq_ref_f, "%s", line);
		}
	}
	fclose(seq_ref_f);
	fclose(seq_f);

	char command[500];
	sprintf(command, "mafft --quiet %s > %s", seq_ref_name, ref_aln_name);
	system(command);
}


/**
 * This function calculates the agreement between two alignments in a specified number of sequences.
 * \param ref_file_name The reference alignment.
 * \param aln_file_name The test alignment.
 **/
double
agreement_score(char *ref_file_name, char *aln_file_name)
{
	const int LINE_LENGTH = 200;
	char line[LINE_LENGTH];

	int number_of_sequences = 0;
	FILE *ref_f = fopen(ref_file_name,"r");
	while(fgets(line, LINE_LENGTH, ref_f)!=NULL)
	{
		if (line[0] == '>')
		{
			++number_of_sequences;
		}
	}
	fseek(ref_f, 0, SEEK_SET);
	char **seq_names = vcalloc(number_of_sequences, sizeof(char*));
	int i = 0;
	while(fgets(line, LINE_LENGTH, ref_f)!=NULL)
	{
		if (line[0] == '>')
		{
			seq_names[i] = vcalloc(LINE_LENGTH, sizeof(char));
			sprintf(seq_names[i], "%s",line);
			++i;
		}
	}


	fclose(ref_f);
	FILE *aln_f = fopen(aln_file_name,"r");
	char *tmp_name = vtmpnam(NULL);
	FILE *tmp_f = fopen(tmp_name,"w");
	int x = 0;
	long last_pos = -1;
	while(fgets(line, LINE_LENGTH, aln_f)!=NULL)
	{
		if (line[0] == '>')
		{
			for (i = 0; i < number_of_sequences; ++i)
			{
				if (!strcmp(line, seq_names[i]))
				{
					fprintf(tmp_f, "%s", line);
					last_pos = ftell(aln_f);
					while(fgets(line, LINE_LENGTH, aln_f)!=NULL)
					{
						if (line[0] == '>')
							break;
						fprintf(tmp_f, "%s", line);

						last_pos = ftell(aln_f);
					}
					++x;
					fseek(aln_f, last_pos, SEEK_SET);
					break;
				}
			}
		}
		if (x == number_of_sequences)
			break;
	}
	fclose(aln_f);
	fclose(tmp_f);
	char command[500];
// 	sprintf(command, "cp %s /users/cn/ckemena/Desktop/1.fa", tmp_name);
// 	system(command);
// 	fp = popen("ls -l", "r");
	sprintf(command, "t_coffee -other_pg aln_compare -al1 %s -al2 %s ", ref_file_name, tmp_name);


	FILE *result;
	result = popen(command, "r");
	fgets(line, LINE_LENGTH, result);
	fgets(line, LINE_LENGTH, result);
	fgets(line, LINE_LENGTH, result);

// 	sprintf(command, "cp %s ~/Destkop/1.fa", tmp_name);
// 	system(command);

	char delims[] = " ";
	char *tmp_str;
	strtok(line, delims);
	for (i = 0; i < 3; ++i)
	{
		while((tmp_str = strtok(NULL, delims))==NULL);
	}
	return(atof(tmp_str));
}


complete_agreement_score(char *aln_file_name, const char *ref_directory)
{
	struct dirent *dp;
	char *name = strrchr(aln_file_name,'/')+1;
// 	printf("%s ",name);
     // enter existing path to directory below
	char directory[500];
	if (ref_directory[0]== '~')
	{
		sprintf(directory,"%s%s",getenv("HOME"),ref_directory+1);
	}
	else
	{
		sprintf(directory,"%s",ref_directory);
	}
	DIR *dir = opendir(directory);

	char ref_file_name[200];
	while ((dp=readdir(dir)) != NULL)
	{
		if ((strcmp(dp->d_name,".")) && (strcmp(dp->d_name,"..")))
		{

			sprintf(ref_file_name, "%s/%s",directory, dp->d_name);
			printf("%s %f\n",name, agreement_score(ref_file_name, aln_file_name));
// 			printf("%f ", agreement_score(ref_file_name, aln_file_name));

		}
	}
	printf("\n");
	closedir(dir);
	return 0;

}


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