/Habitats.hpp

#ifndef HABITATS_H
#define HABITATS_H

#include <math.h>
#include <fstream>

#include "Configuration.hpp"
#include "Aleatory.hpp"
#include "Specie.hpp"
#include "./src-clstr/cluster.h" 


#define MINVAL 0.01 //1% of chance to NOT link two populations. Gives 1% of chance to the closest populations not belong to the same habitat. 
#define MAXVAL 0.99 //1% of chance to link. This gives 1% of chance to the farest populations belong to the same habitat.

struct habitat {
	int h_sp_count;		//Contains how many species are in each habitat.
	int *h_sp;              //This vector contains which species are in each habitat. h_sp[SP_NUMBER]
} *H; 		                //struct that contains the habitats H[SP_NUMBER]	

unsigned short int g_habitatsSize = 0;
double ** g_distances;
Node * g_tree;
int * sp_int;
int ** sp_adj;

//=========================================================================================
// Reserva memória para as variaveis
//=========================================================================================
void initHabitats(int size)
{
	unsigned short int a = 0;
	g_distances = (double**) malloc(size*sizeof(double*));
	for(a = 0; a < size; a++)
	{
		g_distances[a] = (double*) malloc(size*sizeof(double));
	}

	H = (struct habitat*) malloc(size * sizeof(struct habitat));	
	for (a = 0; a < size; a++)
	{
		H[a].h_sp = (int*)malloc (sizeof(int)*size );
		H[a].h_sp_count = 0; 
	}

	sp_int = (int*) malloc (size * sizeof(int));
	sp_adj = (int**)malloc (size * sizeof(int*));
	for (a = 0; a < size; a++)
	{
		sp_adj[a] = (int*)malloc (size * sizeof(int));
	}
}

//=========================================================================================
// Libera a memória
//=========================================================================================
void destroyHabitats(int size)
{
	unsigned short int a = 0;
	for(a = 0; a < size; a++)
	{
		free(g_distances[a]);
	}
	free(g_distances);
	
	for (a = 0; a < size; a++)
	{
		free(H[a].h_sp);
	}
	free(H);

	free(sp_int);
	
	for (a = 0; a < size; a++)
	{
		free(sp_adj[a]);
	}
	free(sp_adj);
}

//=========================================================================================
// Calcula a matriz de distância e normaliza os valores
//=========================================================================================
void buildDistanceMatrix(int size, int dimensions, Specie species[])
{
	double euclid_max = 0;
	double euclid_sum = 0;
	double temp = 0;
	unsigned short int a = 0;
	unsigned short int b = 0;
	unsigned short int d = 0;
	Organism centroid[size];
	for(a = 0; a < size; a++)
	{
		centroid[a] = species[a].getCentroid();
	}
	for(a = 0; a < size; a++)
	{
		g_distances[a][a] = 0;
		for(b = a + 1 ; b < size; b++)
		{
			for(d = 0; d < dimensions; d++)
			{
				euclid_sum += pow(fabs(centroid[a][d] - centroid[b][d]),2.0);
			}
			g_distances[a][b] = sqrt(euclid_sum);
			g_distances[b][a] = g_distances[a][b];
			if (euclid_max < g_distances[a][b])
			{
				euclid_max = g_distances[a][b];
			}
			euclid_sum = 0.0;
		}
	}
	for(a = 0; a < size; a++)
	{
		for(int b = a +1; b < size; b++)
		{
			g_distances[a][b] = g_distances[a][b] / euclid_max;
			g_distances[b][a] = g_distances[a][b];
		}
	}
}
//=========================================================================================
// Gerar o dendograma com o algoritmo singlelinkclustering
//=========================================================================================
void buildDendogram(int size, int dimensions)
{
	double maxrange;
	double minrange;
	unsigned short int i = 0;

	g_tree = treecluster(size, size, 0, 0, 0, 0, 'e', 's', g_distances);
	
	if(!g_tree)
	{
			std::cout << " Erro ao gerar dendograma" << std::endl;
			return;
	}

	maxrange = -1;
	for (i = 0; i < size-1;i++)
	{
		if (maxrange < g_tree[i].distance)
		{
			maxrange = g_tree[i].distance;
		}
	}
	minrange = maxrange;
	for (i = 0; i < size-1; i++)
	{
		if (minrange > g_tree[i].distance)
		{
			minrange = g_tree[i].distance;
		}
	}

	for (i = 0 ;i < size-1;i++)
	{
			g_tree[i].distance=(MAXVAL-MINVAL)/(double)(maxrange-minrange)*g_tree[i].distance-(MAXVAL-MINVAL)/(double)(maxrange-minrange)*minrange+MINVAL;
	}
	
}

void printDendogram(std::string fileName, int size)
{
	std::ofstream output(fileName.c_str());

	output << "MATLAB: ==== Pairwise single linkage clustering using NORMALIZED DIST MATRIX." << std::endl;

	for(unsigned short int i=0; i< (size - 1); i++)
	{
		if (g_tree[i].left >= 0 && g_tree[i].right >= 0)
		    output << g_tree[i].left+1 << g_tree[i].right+1 << "      " << g_tree[i].distance << std::endl;
		if (g_tree[i].left >= 0 && g_tree[i].right < 0)
		    output << g_tree[i].left+1 << (-1*g_tree[i].right)+size << "      " << g_tree[i].distance << std::endl;
		if (g_tree[i].left < 0 && g_tree[i].right >= 0)
		    output << (-1*g_tree[i].left)+size << g_tree[i].right+1 << "      " << g_tree[i].distance << std::endl;
		if (g_tree[i].left < 0 && g_tree[i].right < 0)
		    output << (-1*g_tree[i].left)+size << (-1*g_tree[i].right)+size << "      "  << g_tree[i].distance << std::endl;
	}
}

void buildHabitats(int size)
{
	unsigned short int i = 0;
	unsigned short int j = 0;
	unsigned short int k = 0;
	unsigned short int sum_tabu = 0;
	int cur_node = size - 1;
	int cur_habitat = 0;
	unsigned short int cur_node_index = 0;
	double sum_left = 0.0;
	double sum_right = 0.0;
	double dist_aux = 0.0;

	Aleatory aleatory;

	g_habitatsSize = 0;

	for (i = 0; i < size; i++)
	{
		H[i].h_sp_count = 0; 
	}

	

	while (sum_tabu < size-1)
	{
		if ( aleatory.nextDouble() >= g_tree[cur_node-1].distance) //link those two itens
		{
			if (cur_node == size-1) //first iteration.
			{
				H[cur_habitat].h_sp[ H[cur_habitat].h_sp_count ] = g_tree[cur_node-1].left;
				H[cur_habitat].h_sp[ H[cur_habitat].h_sp_count+1 ] = g_tree[cur_node-1].right;
				H[cur_habitat].h_sp_count = 2;
				sum_tabu++;
			}
			else
			{
				if (cur_node == 1)
				{
					H[cur_habitat].h_sp[ cur_node_index ] = g_tree[0].left;
					H[cur_habitat].h_sp[ H[cur_habitat].h_sp_count ] = g_tree[0].right;
					H[cur_habitat].h_sp_count += 1;
					sum_tabu++;
				}
				else
				{
					H[cur_habitat].h_sp[ cur_node_index ] = g_tree[cur_node-1].left;
					H[cur_habitat].h_sp[ H[cur_habitat].h_sp_count ] = g_tree[cur_node-1].right;
					H[cur_habitat].h_sp_count += 1;
					sum_tabu++;
				}
			}
		}
		else
		{
			if (cur_node == size-1)
			{
				H[cur_habitat].h_sp[ H[cur_habitat].h_sp_count ] = g_tree[cur_node-1].left;
				g_habitatsSize++;
				H[g_habitatsSize].h_sp[ H[g_habitatsSize].h_sp_count ] = g_tree[cur_node-1].right;
				H[cur_habitat].h_sp_count += 1;
				H[g_habitatsSize].h_sp_count += 1;
				sum_tabu++;
			}
			else
			{
				if (H[cur_habitat].h_sp_count == 1 && H[cur_habitat].h_sp[0] < 0)
				{
					H[cur_habitat].h_sp[ cur_node_index ] = g_tree[cur_node-1].left;
					g_habitatsSize++;
					H[g_habitatsSize].h_sp[ H[g_habitatsSize].h_sp_count ] = g_tree[cur_node-1].right;
					H[g_habitatsSize].h_sp_count += 1;
					sum_tabu++;
				}
				if (g_tree[cur_node-1].left >= 0 && g_tree[cur_node-1].right >= 0 && H[cur_habitat].h_sp_count > 1) 
				{
					sum_left = sum_right = 0.0;
					for (i = 0; i < H[cur_habitat].h_sp_count; i++)
					{
					   if (H[cur_habitat].h_sp[ i ] >= 0)
					   {
						sum_left  += g_distances[ H[cur_habitat].h_sp[ i ] ][ g_tree[cur_node-1].left ];
						sum_right += g_distances[ H[cur_habitat].h_sp[ i ] ][ g_tree[cur_node-1].right ];
					   }
					}
					if (sum_left <= sum_right) 
					{
						H[cur_habitat].h_sp[ cur_node_index ] = g_tree[cur_node-1].left;
						g_habitatsSize++;
						H[g_habitatsSize].h_sp[ H[g_habitatsSize].h_sp_count ] = g_tree[cur_node-1].right;
					}
					else
					{
						H[cur_habitat].h_sp[ cur_node_index ] = g_tree[cur_node-1].right;
						g_habitatsSize++;
						H[g_habitatsSize].h_sp[ H[g_habitatsSize].h_sp_count ] = g_tree[cur_node-1].left;
					}
					H[g_habitatsSize].h_sp_count += 1;
					sum_tabu++;
				}
				else if (g_tree[cur_node-1].left < 0 && g_tree[cur_node-1].right < 0 && H[cur_habitat].h_sp_count > 1)
				{
					if ( g_tree[cur_node-1].left < g_tree[cur_node-1].right )
					{
						H[cur_habitat].h_sp[ cur_node_index ] = g_tree[cur_node-1].left;
						g_habitatsSize++;
						H[g_habitatsSize].h_sp[ H[g_habitatsSize].h_sp_count ] = g_tree[cur_node-1].right;
					}
					else
					{
						H[cur_habitat].h_sp[ cur_node_index ] = g_tree[cur_node-1].right;
						g_habitatsSize++;
						H[g_habitatsSize].h_sp[ H[g_habitatsSize].h_sp_count ] = g_tree[cur_node-1].left;
					}					
					H[g_habitatsSize].h_sp_count += 1;
					sum_tabu++;
				}
				else if (H[cur_habitat].h_sp_count > 1) 
				{
					if (g_tree[cur_node-1].left >= 0) 
					{
						H[cur_habitat].h_sp[ cur_node_index ] = g_tree[cur_node-1].left;
						g_habitatsSize++;
						H[g_habitatsSize].h_sp[ H[g_habitatsSize].h_sp_count ] = g_tree[cur_node-1].right;
					}else
					{
						H[cur_habitat].h_sp[ cur_node_index ] = g_tree[cur_node-1].right;
						g_habitatsSize++;
						H[g_habitatsSize].h_sp[ H[g_habitatsSize].h_sp_count ] = g_tree[cur_node-1].left;
					}
					H[g_habitatsSize].h_sp_count += 1;
					sum_tabu++;
				}
			}
		}
		i = 0;
		cur_habitat = -1;
		j = 0;
		while (i <= g_habitatsSize && cur_habitat < 0)
		{
			while (j < H[i].h_sp_count && cur_habitat < 0)
			{
				if (H[i].h_sp[j] < 0)
				{
					cur_habitat = i;
				}
				else j++;				
			}
			i++;
			j = 0;
		}
		if (cur_habitat >= 0)
		{
			i = 0;
			while (H[cur_habitat].h_sp[i] >= 0)
				i++;
			cur_node = -1*H[cur_habitat].h_sp[i];
			cur_node_index = i;
		}
	}
	g_habitatsSize++;


	for (i = 0; i < size; i++)
	{
		sp_int[i] = 0;
		for (j = 0; j < size; j++)
		{
			sp_adj[i][j] = -1;
		}
	}
	
	for ( i = 0; i < g_habitatsSize; i++)
	{
		if (H[i].h_sp_count == 1)
		{
			//update nothing
		}else if ((H[i].h_sp_count == 2))
		{
				sp_adj[ H[i].h_sp[0] ][0] = H[i].h_sp[1];
				sp_int[H[i].h_sp[0]]++;

				sp_adj[ H[i].h_sp[1] ][0] = H[i].h_sp[0];
				sp_int[H[i].h_sp[1]]++;
		}
		else
		{
			for (j=0;j<H[i].h_sp_count;j++)
			{
				k = aleatory.nextInt(H[i].h_sp_count);
				while (sp_int[H[i].h_sp[j]] == 0)
				{
					if (g_distances[ H[i].h_sp[j] ][ H[i].h_sp[k] ] == 1.0)//give a small chance to choose
					{
						dist_aux = 0.99;
					}
					else
					{
						dist_aux = g_distances[ H[i].h_sp[j] ][ H[i].h_sp[k] ];
					}

					if ( (j != k) && aleatory.nextDouble() >= dist_aux) 
					{
						sp_adj[ H[i].h_sp[j] ][sp_int[H[i].h_sp[j]]] = H[i].h_sp[k];
						sp_int[H[i].h_sp[j]]++;
						k = H[i].h_sp_count;
					}
					else 
					{
						k = aleatory.nextInt(H[i].h_sp_count);
					}
				}
			}
		}
	}

}

void printPairwiseInteractions(std::string fileName, int size)
{
	unsigned short int i = 0;
	unsigned short int j = 0;
	std::ofstream output;	
  	output.open(fileName.c_str(), std::ios_base::app);
	for(i = 0; i < size; i++)
	{
		output << "Specie :: " << i << " => ";
		for(j = 0; j < sp_int[i]; j++)
		{
			output << sp_adj[i][j] << "; ";
		}
		output << "" << std::endl;
	}
}


#endif