//段海濱教授主編的《蟻群算法原理及其應(yīng)用》附錄里的C程序代碼.  
//Basic Ant Colony Algorithm for TSP  
#include <iostream>  
#include <fstream>  
#include <math.h>  
#include <time.h>  
#include <conio.h>  
#include <stdlib.h>  
#include <iomanip>  
using namespace std;
#define N 51//city size  
#define M 31 //ant number  

double inittao=1; //初始信息量的多少.   
double tao[N][N]; //相應(yīng)路徑上的信息素 . 
double detatao[N][N]; //相應(yīng)路徑上的信息素增量. 可能有的資料里寫成:Δτ   
double Distance[N][N]; //距離矩陣  
double yita[N][N];  
int tabu[M][N];  
int route[M][N]; //route[i][j]：保存螞蟻i的經(jīng)過的路徑的數(shù)組   
double solution[M]; 
//int score[NcMax][M];//保存某次循環(huán)中路徑最優(yōu)的螞蟻號
int BestRoute[N];  //保存最佳路徑
double BestSolution=10000000000;  
double alfa,beta,rou,Q;  
 int NcMax;  
void initparameter(void); // initialize the parameters of basic ACA  
double EvalueSolution(int *a); // evaluate the solution of TSP, and calculate the length of path  
void InCityXY( double x[], double y[], char *infile ); // input the nodes' coordinates of TSP  

void initparameter(void) //初始化參數(shù)  
{  
alfa=1.5; beta=4.0; rou=0.6; Q=50;  
NcMax=200; //最大迭代次數(shù)  
}  

void main(void)  
{  
int NC=0;  
initparameter();  
double x[N];  
double y[N];  
InCityXY( x, y, "eil51.tsp" );  

//初始化距離矩陣  
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  //計(jì)算啟發(fā)函數(shù)
for( int i=0;i<N;i++)  
 for(int j=0;j<N;j++)  
 {  
 tao[i][j]=inittao;  
 if(j!=i)  
 yita[i][j]=1/Distance[i][j];  
 }  
for(int k=0;k<M;k++)  
 for(int i=0;i<N;i++)  
 route[k][i]=-1;  
 srand((unsigned)time(NULL));  
for(int k=0;k<M;k++)  
{  
 route[k][0]=k%N;  //隨機(jī)分配螞蟻的初始位置
 tabu[k][route[k][0]]=1;  
}  

//each ant try to find the optiamal path  
do  
{  
 int s=1;  
 double partsum;  
 double pper;  
 double drand;  
//static int tag=0;
// int first;
 //ant choose one whole path  
 while(s<N)  
 {  
 //開始計(jì)算第k只螞蟻的路徑  
 for(int k=0;k<M;k++)  
 {  
 int jrand=rand()%3000; //為了生成一個(gè)0到3000之間的隨機(jī)數(shù).   
  
 drand=jrand/3001.;  
 partsum=0;  
 pper=0; //（轉(zhuǎn)移概率） 
  int j;
 //根據(jù)概率函數(shù)計(jì)算螞蟻的轉(zhuǎn)移概率  
 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++;  
 }  
 //在N次循環(huán)后，所有螞蟻的禁忌表都已填滿  
  
 //計(jì)算每個(gè)螞蟻?zhàn)哌^的路徑的長度，并找到最短路徑保存，記錄此路徑并更改信息素。  
 // the pheromone is updated 更新信息素  
  
 for( int 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];
// score[NC][k]=1;
  
 }  
 }  

 for(int 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(int 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(int k=0;k<M;k++)  
 for(int j=1;j<N;j++)  
 {  
 tabu[k][route[k][j]]=0;  
 route[k][j]=-1; //將螞蟻的路徑再重新置空,為下一次循環(huán)做準(zhǔn)備,要不然每個(gè)螞蟻的路徑都已經(jīng)滿了,則沒有辦法進(jìn)行下一次迭代了.  
 }  
 NC++;  
} while(NC<NcMax);  


//output the calculating results  
fstream result;  
result.open("optimal_results.log", ios::app);  
if(!result)  
{  
 cout<<"can't open the <optimal_results.log> file!\n";  
 exit(0);  
}  
result<<"*-------------------------------------------------------------------------*"<<endl;  
//result<<"the initialized parameters of ACA are as follows:"<<endl;  
result<<"ACA的參數(shù)初始化如下:"<<endl;  

result<<"alfa="<<alfa<<", beta="<<beta<<", rou="<<rou<<", Q="<<Q<<endl;  
result<<"the maximum iteration number of ACA is:"<<NcMax<<endl;  
result<<"the shortest length of the path is:"<<BestSolution<<endl;  
result<<"the best route is:"<<endl;  
for(int i=0;i<N;i++)  
 result<<BestRoute[i]<<" ";  
result<<endl;  
result<<"*-------------------------------------------------------------------------*"<<endl<<endl;  
result.close(); 
 
}  


double EvalueSolution(int *a)  
{  
double dist=0;  int i;
for(  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 );  
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;  
}  
}