#include <stdlib.h>
#include <stdio.h>
#include <iostream>
#include <time.h>
#include <algorithm>		/*	find	*/
#include <vector>			/*	vector	*/
#include <fstream>			/*	to read and write file	*/
#include <string>       	/*	std::string, std::to_string	*/
#include <iterator>  		/*	std::begin, std::end	*/
#include <sstream>
#include <math.h>			/*	sqrt pow	*/
#include <cfloat>			/*	Limits of floating point types	*/

using namespace std;

typedef vector<int> i1d;
typedef vector<i1d> i2d;
typedef vector<double> d1d;
typedef vector<d1d> d2d;
typedef vector<bool> b1d;
typedef vector<b1d> b2d;

i1d shortest_tour;
i2d ant_position;
b2d tabu;
d2d trail;
d2d delta_trail;
d2d distanceXY;
d2d locationXY;
d2d transportation_probability;
i1d unused_position;

int edge_weight_type;
int dimension;
int sum_ant;
int Nc = 0;		//NC is the cycles counter
int NcMax;
int S;

double alpha;
double beta;
double rho;
double tau;
double Q;
double shortest_distance = DBL_MAX;

string dataset_name;

void 	init();
int 	random(int low, int up=0);
double 	random(double low, double up=0.0);
void	read_file(string filename);
void	write_file(string ,string);	//[filename, context]
void	file_context_analyse(string, int);
double	Euclidean_distance(d1d, d1d);
void	caculate_distance();
void	run();
double	get_distance(int, int);
double	get_trail(int, int);
void 	ant_system();
void 	first_trip();
void	delta_trail_update(int, int, double);
void	trail_update_and_reset_delta_trail();
void	unused_position_reset();
void	unused_position_used(int);
double	caculate_transition_probability(int, int, int);

int main(int argc, char* argv[])
{
	srand( (unsigned)time(NULL) );
	
	NcMax = 100;
	tau=0.0;
	Q=100.0;					// not sure what data to set
	alpha = 0.1;
	beta = 2;
	rho = 0.9;
	
	edge_weight_type = 2;
	dimension = 51;	dataset_name = "tsp51.txt";
	sum_ant = dimension;
	
	run();

	cout << "End of the Ant optimization." << endl;
	system("pause");
	return 0;
}

void run(){
	init();
	read_file(dataset_name);
	caculate_distance();
	ant_system();
}

void init(){
	ant_position.assign(sum_ant, i1d(dimension+1, 0));			/*need go to start at end of the trip*/
	shortest_tour.assign(dimension+1, 0);
	tabu.assign(sum_ant, b1d(dimension, true));
	trail.assign(dimension, d1d(dimension, tau));				/*	trail[i][j]	j>=i	*/
	delta_trail.assign(dimension, d1d(dimension, 0.0));			/*	delta_trail[i][j]	j>=i	*/
	locationXY.assign(dimension, d1d(edge_weight_type, 0.0));
	distanceXY.assign(dimension, d1d(dimension, 0.0));			/*	distanceXY[i][j]	j>=i*/
	transportation_probability.assign(sum_ant, d1d(dimension, 0.0));
	unused_position.assign(dimension, 0);
	for(int i=0; i<dimension; i++)	unused_position[i] = i;
}

void ant_system(){
	first_trip();
	
	do{
		trail_update_and_reset_delta_trail();
		Nc++;
		
		unused_position_reset();
		
		int this_position , next_position;
		
		for(int k=0; k<sum_ant; k++){		//k is index of ant
			unused_position_reset();
			this_position = random(dimension-1);
			
			ant_position[k][0] = this_position;
			unused_position_used(this_position);
			
			for(int s=1; s<dimension; s++){			//s is index of tabu, is number of cities
				double q;
				do{
					//if q=1.0 , r will be something wrong
					q = random(1.0);
				}while(q==1.0);
				int i = 0;
				double r = caculate_transition_probability(this_position, unused_position[i], k);;
				
				while( r<q && i<unused_position.size() ){
					i++;
					double p = caculate_transition_probability(this_position, unused_position[i], k);
					r += p;
				}
				next_position = unused_position[i];
				
				unused_position_used(next_position);
				ant_position[k][s] = next_position;
				this_position = next_position;
			}
		}
	}while(Nc<NcMax);
	
}

void first_trip(){
	int next_position;
	/*Place the k ants on the s nodes, first time all ants walk randomly without repeat city*/
	for(int k=0; k<sum_ant; k++){		/*k is index of ant*/
		unused_position_reset();
		for(int s=0; s<dimension; s++){			/*s is index of tabu, is number of cities*/
			next_position = unused_position[ random( (int)unused_position.size() - 1 ) ];
			unused_position_used(next_position);
			
			ant_position[k][s] = next_position;
		}
	}
	Nc++;
}

void trail_update_and_reset_delta_trail(){
	for(int k=0; k<sum_ant; k++){								/*k is the index of ant*/
		ant_position[k][dimension] = ant_position[k][0];		/*set last position equal to first position*/
		
		double distance=0.0;
		
		for(int s=0; s<dimension; s++)
			distance += get_distance(ant_position[k][s], ant_position[k][s+1]);
		for(int s=0; s<dimension; s++){
			delta_trail_update(ant_position[k][s], ant_position[k][s+1], distance);
		}
		
 		if(shortest_distance>distance){
			shortest_distance = distance;
			shortest_tour.assign(ant_position[k].begin(), ant_position[k].end());
			cout << "NO." << Nc << " ant= "<< k << "\tshortest_distance = " << shortest_distance << endl;
			
			/* This shown the optimal position  
			for(int i=0; i<dimension; i++){
				cout << ant_position[k][i] << " ";
			}
			cout << endl;
			/* End of shown the optimal position  */
		}
	}
	
	//
	for(int i=0; i<distanceXY.size(); i++){
		for(int j=i; j<distanceXY.size(); j++){
			trail[i][j]  = trail[i][j]*rho + delta_trail[i][j];
			delta_trail[i][j] = 0.0;				
		}
	}
}

double caculate_transition_probability(int this_position, int next_position, int k){
	double p, q=0.0;
	
	p  = pow( get_trail(this_position, next_position), alpha );
	p *= pow( 1/get_distance(this_position, next_position), beta );
	for(int j=0; j<(int)unused_position.size(); j++){
		q += ( pow( get_trail(this_position, unused_position[j]), alpha ) * pow( 1/get_distance(this_position, unused_position[j]), beta ) );
	}
	return p/q;
}

void delta_trail_update(int i, int j, double d){
	if(i>j)	delta_trail[j][i] += Q/d;
	else	delta_trail[i][j] += Q/d;
}

void unused_position_reset(){
	unused_position.assign(dimension, 0);
	for(int i=0; i<dimension; i++)	unused_position[i] = i;
}

void unused_position_used(int position){
	vector<int>::iterator pos;
	pos = find(unused_position.begin(), unused_position.end(), position);
	unused_position.erase(pos);
}

int random(int low, int up){						/*	low<=random(low,up)<=up	*/
	if(low>up)	swap(low,up);						/*		0<=random(a)<=a		*/
	return rand() % (up - low + 1) + low ;
}

double random(double low, double up){ 				/* low<=random(low,up)<=up	*/
	if(up<low)	swap(low,up);						/*		0<=random(a)<=a		*/
	return (up - low) * rand() / (RAND_MAX ) + low;
}

void read_file(string filename){
//	http://w...content-available-to-author-only...s.com/doc/tutorial/files/
	int count = 0;
	string line;
	ifstream myfile (filename.c_str());
	if (myfile.is_open()){
		while ( getline (myfile,line) ){
			file_context_analyse(line, count);
			count++;
		}
		myfile.close();
	}else{
		cout << "Read file error. Please check again." << endl;
		exit(1);
	}
}

void write_file(string filename, string context){
	ofstream myfile (filename.c_str());
	if (myfile.is_open()){
		myfile << context;
		myfile.close();
	}else {
		cout << "Unable to open file" << endl;
		exit(1);
	}
}

void file_context_analyse(string str, int count){
	istringstream ss(str);
    istream_iterator<double> begin(ss), end;
	
    //putting all the tokens in the vector
    d1d arrayTokens(begin, end);
	
	/*check the dimension is as same as file*/
    if(edge_weight_type!=arrayTokens.size()){
    	cout << "Dimension setting error, Please check again." << endl;
    	exit(1);
	}
	
    /*Save context to locationXY*/
    for(int i=0; i<arrayTokens.size(); i++){
    	locationXY[count][i] = arrayTokens[i];
	}
}

double Euclidean_distance(d1d location1, d1d location2){
	//https://e...content-available-to-author-only...a.org/wiki/Euclidean_distance
	double d=0.0;
	for(int i=0; i<location1.size(); i++){
		d += (pow(location1[i]-location2[i],2));
	}
	return sqrt(d);
}

void caculate_distance(){
	for(int i=0; i<locationXY.size(); i++){
		for(int j=i; j<locationXY.size(); j++){
			if(i==j){
				distanceXY[i][j] = 0;
			}else{
				distanceXY[i][j] = Euclidean_distance(locationXY[i], locationXY[j]);
			}
		}
	}
}

double get_distance(int i, int j){
	if(i>j)	/*swap(i, j);*/ return( distanceXY[j][i] );
	else return( distanceXY[i][j] );
}

double get_trail(int i, int j){
	if(i>j)	/*swap(i, j);*/ return( trail[j][i] );
	else return( trail[i][j] );
}