/*
	Algorithm Name:	AS-QAP
	Version:		0.9 ---- need verification
	Author:			Chen Ye, Sichuan University, China
	E-mail:			arrowcy@163.com
	Reference:		Thomas Stutzle, Macro Dorigo.
					ACO Algorithms for Quadratic Assignment Problem,
					It is Chapter 3 of the book New Ideas in Optimization edited by
					D.Corne, M.Dorigo & F.Glover.
*/

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

#define ERROR			-1
#define SUCCESS			0

#define FILE_NAME		"chr12a.dat"

#define	MAX_ITER		500
#define PROBLEM_SIZE	12
#define NUM_ANTS		20
#define	INIT_PHEROMONE	0.01	///???
#define	ALPHA			1
#define	BETA			1
#define Q				100
#define	RHO				0.9

int   problem_size;

float heuristic	[PROBLEM_SIZE][PROBLEM_SIZE];
int   distance	[PROBLEM_SIZE][PROBLEM_SIZE];
int   dsum		[PROBLEM_SIZE];
int   flow		[PROBLEM_SIZE][PROBLEM_SIZE];
int   fsum		[PROBLEM_SIZE];
float pheromone	[PROBLEM_SIZE][PROBLEM_SIZE];
int   tabu		[NUM_ANTS][PROBLEM_SIZE];
float fit		[NUM_ANTS];	///////integer works properly.

int   global_best_tour[PROBLEM_SIZE];
float global_best_fit;

int   ant		[NUM_ANTS][PROBLEM_SIZE];

int   assignment_order[PROBLEM_SIZE];
int   display_order   [PROBLEM_SIZE];

int read_data()
{
	int i,j;
	FILE *fp;

	if(!(fp=fopen(FILE_NAME,"r")))
	{
		printf("Error opening problem description file.\n");
		return ERROR;
	}

	fscanf(fp,"%d",&problem_size);
	if (problem_size>PROBLEM_SIZE)
	{
		printf("PROBLEM_SIZE is too small.\n");
		return ERROR;
	}

	//no error detection will be take below
	for (i=0;i<problem_size;i++)
		for (j=0;j<problem_size;j++)
			fscanf(fp,"%d",&flow[i][j]);
	for (i=0;i<problem_size;i++)
		for (j=0;j<problem_size;j++)
			fscanf(fp,"%d",&distance[i][j]);
	fclose(fp);

	return SUCCESS;
}

void calc_heuristic()
{
	int i,j;
	for (i=0;i<problem_size;i++)
	{
		dsum[i]=0;
		fsum[i]=0;
		for (j=0;j<problem_size;j++)
		{
			if(i!=j)
			{
				dsum[i]=dsum[i]+distance[i][j];
				fsum[i]=fsum[i]+flow[i][j];
			}
		}
	}
	for (i=0;i<problem_size;i++)
		for (j=0;j<problem_size;j++)
			heuristic[i][j]=1.0/((float)fsum[i]*dsum[j]);
}

void calc_assignment_order()
{
	int i,j,bestf,bestp;
	int assigned[PROBLEM_SIZE];

	for (i=0;i<problem_size;i++)
	{
		assignment_order[i]=-1;
		assigned[i]=-1;
	}
	for (i=0;i<problem_size;i++)
	{
		bestf=0;
		bestp=-1;
		for (j=0;j<problem_size;j++)
		{
			if(bestf<fsum[j] && assigned[j]==-1)
			{
				bestf=fsum[j];
				bestp=j;
			}
		}
		if(bestp==-1)
			printf("Unknown Error!\n");

		assignment_order[i]=bestp;
		assigned[bestp]=i;
		display_order[i]=bestp;
	}
}

void init_pheromone()
{
	int i,j;
	for (i=0;i<problem_size;i++)
		for (j=0;j<problem_size;j++)
			pheromone[i][j]=INIT_PHEROMONE;
}

void solution_construction()
{
	int i,j,k;
	float sum;
	int cur_facility;
	float prob,sumprob;
	float t[PROBLEM_SIZE];

	//////// clear tabu
	for (i=0;i<NUM_ANTS;i++)
		for (j=0;j<problem_size;j++)
			tabu[i][j]=0;

	for (i=0;i<problem_size;i++)
	{
		//// each ant will take a step in this loop
		cur_facility=assignment_order[i];
		for (j=0;j<NUM_ANTS;j++)
		{
			sum=0;
			for (k=0;k<problem_size;k++)
			{
				t[k]=0;		//possibly not needed
				if(tabu[j][k]==0)
					t[k]=pow(pheromone[cur_facility][k],ALPHA)*pow(heuristic[cur_facility][k],BETA);
				sum=sum+t[k];
			}
			prob=sum*(float)rand()/RAND_MAX;
			sumprob=0;
			for (k=0;k<problem_size;k++)
			{
				if (tabu[j][k]==0)
				{
					sumprob=sumprob+t[k];
					if(sumprob>prob)
					{
						/////////// a location for the current facility has been found
						ant[j][i]=k;
						tabu[j][k]=1;
						break;
					}
				}
			}
			if(k==problem_size)
			{
				ant[j][i]=problem_size-1;
				tabu[j][problem_size-1]=1;
			}
		}
	}
}

void calc_fit()
{
	int i,j,k,f1,f2;
	float s;

	for (i=0;i<NUM_ANTS;i++)
	{
		s=0;
		for (j=0;j<problem_size;j++)
		{
			f1=assignment_order[j];
			for (k=0;k<problem_size;k++)
			{
				f2=assignment_order[k];
				s=s+flow[f1][f2]*distance[ant[i][j]][ant[i][k]];
			}
		}
		fit[i]=s;
	}
}

void pheromone_update()
{
	//float dph;
	int i,j;

	for (i=0;i<problem_size;i++)
		for(j=0;j<problem_size;j++)
			pheromone[i][j]=pheromone[i][j]*RHO;

	for (i=0;i<NUM_ANTS;i++)
		for (j=0;j<problem_size;j++)
			pheromone[assignment_order[j]][ant[i][j]] =
				pheromone[assignment_order[j]][ant[i][j]] + Q/fit[i];

}

void find_best(int first_time)
{
	int i;
	float local_best_fit;
	int   local_best_ant;

	local_best_ant=0;
	local_best_fit=fit[0];

	for (i=1;i<NUM_ANTS;i++)
	{
		if(fit[i]<local_best_fit)
		{
			local_best_ant=i;
			local_best_fit=fit[i];
		}
	}

	if(first_time || (!first_time && global_best_fit>local_best_fit))
	{
		global_best_fit=local_best_fit;
		for (i=0;i<problem_size;i++)
			global_best_tour[i]=ant[local_best_ant][i];
	}

	printf("best fit: %f",global_best_fit);
	printf(" the tour is: ");
	for (i=0;i<problem_size;i++)
		printf(" %d",global_best_tour[i]);
	printf("\n");
}

int main()
{
	int iter, i;
	unsigned int seed=time(NULL);

	srand(seed);

	if(read_data()==ERROR)
	{
		printf("Program cannot countinue.\n");
		return -1;
	}

	calc_heuristic();
	init_pheromone();
	calc_assignment_order();

	for (iter=0; iter<MAX_ITER; iter++)
	{
		printf("Cur Iter=%d\n",iter+1);
		solution_construction();
		calc_fit();
		pheromone_update();
		find_best(iter==0);
	}

	printf("Seed=%d\n",seed);
	printf("Best fit=%f\n",global_best_fit);
	printf("Best tour=");
	for (i=0;i<problem_size;i++)
		printf(" %d",global_best_tour[display_order[i]]);
	printf("\n");

	return 0;
}