#include <iostream.h>
#include <fstream.h>
#include <math.h>
#include <time.h>
#include <conio.h>
#include <stdlib.h>
#include <iomanip.h>
//using namespace std;
#define N 51
#define M 30
double inittao=1;
double tao[N][N];
double detatao[N][N];
double distance[N][N];
double yita[N][N];
int tabu[M][N];
int route[M][N];
double solution[M];
int BestRoute[N];
double BestSolution=10000000000;
double alfa,beta,rou,Q;
int NcMax;
void initparameter();
double EvalueSolution(int* a);
void InCityXY(double x[],double y[],char *infile);

void main()
{
	int NC=0;
	initparameter();
	double x[N],y[N];
	InCityXY(x,y, "EIL51.TXT");
	for(int i = 0; i<N; i++)
		for(int j=i+1; j<N; j++)
		{
			distance[j][i] = sqrt((x[i]-x[j])*(x[i]-x[j])+(y[i]-y[j])*(y[i]-y[j]));
			distance[i][j] = distance[j][i];
		}

	//calculate the heuristic parameters
	for(i=0; i<N; i++)
		for(int j=0; j<N; j++)
		{
			tao[i][j] = inittao;
			if(j != i)
				yita[i][j] = 100/(distance[i][j] + 50 - int(distance[i][j])%50);
		}
	for(int k=0; k<M; k++)
		for(i=0; i<N; i++)
			route[k][i]=-1;
	for(k=0; k<M; k++)
	{
		route[k][0]=(k+N) % N;
		tabu[k][route[k][0]]=1;
	}

	//each ant try to find the optimal path
	srand( time(NULL));
	do
	{
		int s=1;
		double partsum;
		double pper;
		double drand;

		//ant choose one whole path
		while(s<N)
		{
			for(int k=0; k<M; k++)
			{
				int jrand = rand()%3000;
				drand = double(jrand) / 3001;
				partsum = 0;
				pper = 0;
				for(int j=0; j<N; j++)
				{
					if(tabu[k][j] == 0)
						partsum += pow(tao[route[k][s-1]][j], alfa)* pow(yita[route[k][s-1]][j], beta);
				}
				for(j=0; j<N; j++)
				{
					if(tabu[k][j] == 0)
						pper += pow(tao[route[k][s-1]][j], alfa)* pow(yita[route[k][s-1]][j], beta) / partsum;
					if(pper>drand) break;
				}

				tabu[k][j]=1;
				route[k][s]=j;
			}
			s++;
		}

		//update pheromone
		for(i=0; i<N; i++)
			for(int j=0; j<N; j++)
			{
				detatao[i][j]=0;
			}
		for(int k=0; k<M; k++)
		{
			solution[k] = EvalueSolution(route[k]);
			if(solution[k] < BestSolution)
			{
				BestSolution = solution[k];
				for(s=0; s<N; s++)
					BestRoute[s] = route[k][s];
			}
		}
		for(k=0; k<M; k++)
		{
			for(s=0; s<N-1; s++)
				detatao[route[k][s]][route[k][s+1]] += Q/solution[k];
			detatao[route[k][N-1]][route[k][0]] += Q/solution[k];
		}
		for(i=0; i<N; i++)
			for(int j=0; j<N; j++)
			{
				tao[i][j]=rou*tao[i][j] + detatao[i][j];
				if(tao[i][j] < 0.00001)
					tao[i][j] = 0.00001;
				if(tao[i][j] > 20)
					tao[i][j] =20;
			}
		for(k=0; k<M; k++)
			for(int j=1; j<N; j++)
			{
				tabu[k][route[k][j]]=0;
				route[k][j]=-1;
			}
		NC++;
	}while(NC<NcMax);

	//output the calculating results
	fstream result;
	result.open("Optimal_Result.dat",ios::app);
	if(!result)
	{
		cout<<"Can't open the 'Optimal_Result' file\n";
		abort();
	}

	result<<"*****************************************"<<endl;
	result<<"The initialized parameter of ACA are as follows:"<<endl;
	result<<"alfa="<<alfa<<", beta="<<beta<<endl;
	result<<"Q="<<Q<<", The maximum iteration number of ACA is :"<<NcMax<<endl;
	result<<"The evaporaion parameter of ACA is:"<<rou<<endl;
	result<<"*****************************************"<<endl;
	result<<"The shortest path is :"<<BestSolution<<endl;
	for(i=0; i<N; i++)
		result<<BestRoute[i]<<" ->";
	result<<endl;
	result.close();

	cout<<"The shortest path is :"<<BestSolution<<endl;
}

void initparameter()
{
	alfa=1;
	beta=5;
	rou=0.9;
	Q=100;
	NcMax=200;
}
double EvalueSolution(int* a)
{
	double dist=0;
	for(int i=0; i<N-1; i++)
		dist += distance[a[i]][a[i+1]];
	dist += distance[a[i]][a[0]];
	return dist;
}
void InCityXY(double x[],double y[], char *infile)
{
	fstream inxyfile( infile, ios::in | ios::nocreate);
	if(!inxyfile)
	{
		cout<<"Can't open the <"<<infile<<">file!\n";
		exit(0);
	}
	int i=0;
	while( !inxyfile.eof())
	{
		inxyfile>>x[i]>>y[i];
		if( ++i >= N) break;
	}
}