//SIE 540 Project: Metaheuristic Optimization
//Solving a symmetric Quadratic Assignment Problem
//author: Benjamin Armbruster et. al.
//date: April 11, 2003

//todo: variable problem_size, get stats (depth of valleys),
//runs with varying parameters: seem to have little effect
//10 restarts usually (~90% of time) suffices to get optimal
//maybe get stats with exhaustive search
//---------------------------------------------------------------------------
#ifdef __BORLANDC__
	#define _NO_VCL
	#include <sysdefs.h>              //this is needed to understand USEUNIT
	//USEUNIT("otherUnit1.cpp"); 		 //do USEUNIT for all other source files in
	//USEUNIT("otherUnit1.cpp"); 		//the main source file
	//USEUNIT("otherUnit1.cpp");		   //this is needed to do multi file projects
	#include <stdlib.h>
	#include <iostream.h>
	void __exit_func(void);
   void __exit_func(void)    //so it doesn't exit immediatly
      {
      cout<<endl<<"***Enter a character to exit***"<<endl;
      char dummy;
      cin>>dummy;
      }
#endif  //Borland C
//---------------------------------------------------------------------------

//System Headers
#include <iostream.h> //for input output
#include <fstream.h>  //for reading from file
#include <vector.h> //for permutations in exhaustive search
#include <deque.h>  //for tabu list
#include <math.h>	//for exp() in simulated annealing
#include <stdlib.h> //for srand() and rand()
#include <time.h> //for time()
#include <sys/timeb.h> //for ftime() getting time to 10^-3 sec

//Local Headers or Function/Object Prototypes
int main(void);
int cost(int assgnmt[]);
int delta_cost(int old_assgnmt[],int area1,int area2);
inline void swap(int& a,int& b);
bool in_tabulist(int assgnmt[]);
void init_data(void);
void local_search(void);
void new_best(void);	//checks if it is a new best solution
int best_neighbor(int& area1,int& area2); //outputs best swap in area1,area2 and returns deltacost
void create_stats(void);
void interpret_perm(int assgnmt[]);	//prints to cout an object-room matching
//void dealloc_data(void);
void rand_assgnmt(int assgnmt[]);
void run_metaheuristic(void (*heuristic)(void));
//meta heuristics
void random_search(void);
void exhaustive_search(void);
void simulated_annealing(void);
void tabu_search(void);
//tabu list helpers
void tabu_list(int assgnmt[]);
//GA helpers
void genetic_algorithm(void);
void cull_pop(void);
struct pop_member; //forward declaration
struct pop_member {
	vector<int> permutation;
        int p_cost;
        int created; //iteration created
        pop_member(int iteration);
        pop_member(pop_member &m,pop_member &f,int iteration); //crossover constructor
};
//comparison operators <
inline bool stupider(pop_member &left,pop_member &right);
inline bool older(pop_member &left,pop_member &right);

#define problem_size 10

int assgn_penalties[problem_size][problem_size]; //[item][area]
int relation[problem_size][problem_size]; //relation&dist are triangular so may be merged
int dist[problem_size][problem_size];

int assgn_penalties_weight=1;
int relation_weight=1;

int best_cost;
int best[problem_size];	//assgnmt[i] is item in area i
int curr[problem_size];

//tabu search parameters
unsigned int tabu_iterations=1000;
unsigned int tabu_size_list=7;
deque<vector<int> > tabulist;

//simulated annealing parameters
float T=10;
float r=0.95;
unsigned int num_same_temp=10;
unsigned int temp_steps=100;

//GA parameters
int pop_size=10; //number of population's solutions
float mutation_prob=0.05;
unsigned int GA_iterations=100;

//Function Definitions
int main(void) {
   #ifdef __BORLANDC__
      atexit(__exit_func);
   #endif   //Borland

   srand(time(NULL));
   init_data();

   run_metaheuristic(random_search);
   run_metaheuristic(simulated_annealing);
   run_metaheuristic(tabu_search);
   run_metaheuristic(genetic_algorithm);

   run_metaheuristic(exhaustive_search);
   create_stats(); //must be right after exhaustive search as create_stats assumes best=optimal

   //dealloc_data;   //maybe needed when go to variable problem_sizes
	return 0;
}

int cost(int assgnmt[]) {	//assgnmt[i] is item in area i
	int cost=0;
	for(int i=0;i<problem_size;i++) cost+=assgn_penalties[assgnmt[i]][i];
	for(int i=0;i<problem_size;i++)
		for(int j=0;j<i;j++) cost+=dist[i][j]*relation[assgnmt[i]][assgnmt[j]];
	return cost;
}

int delta_cost(int old_assgnmt[],int area1,int area2) {
	//optimized version O(problem_size)
	int delta=0;
	delta-=assgn_penalties[old_assgnmt[area1]][area1];  //assignment costs
	delta-=assgn_penalties[old_assgnmt[area2]][area2];
	delta+=assgn_penalties[old_assgnmt[area1]][area2];
	delta+=assgn_penalties[old_assgnmt[area2]][area1];
	//relationship costs
	for(int j=0;j<problem_size;j++){
		delta+=dist[area1][j]*(relation[old_assgnmt[area2]][old_assgnmt[j]]
			-relation[old_assgnmt[area1]][old_assgnmt[j]]);
		delta+=dist[area2][j]*(relation[old_assgnmt[area1]][old_assgnmt[j]]
			-relation[old_assgnmt[area2]][old_assgnmt[j]]);
	}
	delta+=2*dist[area1][area2]*relation[old_assgnmt[area1]][old_assgnmt[area2]];

	//unoptimized version O(problem_size^2)
	/*
	int new_assgnmt[10];
	for(int i=0;i<problem_size;i++) new_assgnmt[i]=old_assgnmt[i];
	swap(new_assgnmt[area1],new_assgnmt[area2]);
	delta=cost(new_assgnmt)-cost(old_assgnmt); */

	return delta;
}

inline
void swap(int& a,int& b) {
	int c=b;
	b=a;
	a=c;
}

void init_data(void) {
	std::string fileName;
	cout<<"In which file is the problem stored? ";
	cin>>fileName;

	cout<<"What weight should be given to assignment costs? [int] ";
	cin>>assgn_penalties_weight;
	cout<<"What weight should be given to relationships? [int] ";
	cin>>relation_weight;

	cout<<endl;
	ifstream file(fileName.c_str());
	if(file.fail()) {
		cout<<"error opening file";
		exit(1);
	}

	//cin>>problem_size;
	//assgn_penalties=new int[problem_size][problem_size];
	//relation=new int[problem_size][problem_size];
	//dist=new int[problem_size][problem_size];
	for(int i=0;i<problem_size;i++) {
		for(int j=0;j<problem_size;j++) {
			file>>assgn_penalties[i][j];
			assgn_penalties[i][j]*=assgn_penalties_weight;
		}
	}
	for(int i=0;i<problem_size;i++) {
		for(int j=0;j<i;j++) {
			file>>relation[i][j];
			relation[i][j]*=relation_weight;
			relation[j][i]=relation[i][j];
		}
		relation[i][i]=0;
	}
	for(int i=0;i<problem_size;i++){
		for(int j=0;j<i;j++) {
			file>>dist[i][j];
			dist[j][i]=dist[i][j];
		}
		dist[i][i]=0;
	}
	file.close();

        cout<<"Favorite tabu list size? [int] ";
        cin>>tabu_size_list;
        cout<<endl;

        cout<<"Starting temperature? [float] ";
        cin>>T;
        cout<<"Number of different temperature levels? [int] ";
        cin>>temp_steps;
        cout<<"Ratio between temperature levels? [float] ";
        cin>>r;
        cout<<"Number of iterations per temperature level? [int] ";
        cin>>num_same_temp;
        cout<<endl;

        cout<<"Number of iterations for GA? [int] ";
        cin>>GA_iterations;
        cout<<"Mutation probability? [float] ";
        cin>>mutation_prob;
        cout<<"Population size for GA? [int] ";
        cin>>pop_size;
        cout<<endl;
}

/*
void dealloc_data(void) {
	for(int i=0;i<problem_size;i++) {
		delete[] assgn_penalties[i];
		delete[] relation[i];
		delete[] dist[i];
	}
	delete[] assgn_penalties;
	delete[] relation;
	delete[] dist;
} */

void rand_assgnmt(int assgnmt[]) {
	//an initial permutation
	for(int i=0;i<problem_size;i++) assgnmt[i]=i;

	//this then recursively creates bigger and bigger random permutations
	//at the beginning of step i:
		//0..i-1 form a random permutation in the spots 0..i-1 (the induction hypothesis)
	//at the end of step i:
		//0..i form a random permutation in the spots 0..i
	//because step i, swaps spot i with a random spot from 0..i
	for(int i=0;i<problem_size;i++) swap(assgnmt[i],assgnmt[rand()%(i+1)]);
}

void run_metaheuristic(void (*heuristic)(void)) {
	//if one is concerned about speed one can move init stuff and polishing to heuristic

	//init stuff: random starting point, clear best incumbent
	//note: uncessary for exhaustive search
	rand_assgnmt(curr);
	for(int i=0;i<problem_size;i++) best[i]=curr[i];
	best_cost=cost(curr);

	//time_t start=time(NULL);
	timeb start,end;	//this not ANSI may need to use above but has better resolution
	::ftime(&start);	//"::" needed to resolve ambiguity

	heuristic();

	//polishing the best solution
	//note: unecessary for random_search
	//note: uncessary for exhaustive search
	for(int i=0;i<problem_size;i++) curr[i]=best[i];
	local_search();
	new_best();

	//float timediff=time(NULL)-start;	//ANSI version
	::ftime(&end);	//"::" needed to resolve ambiguity
	float timediff=end.time-start.time+0.001*(end.millitm-start.millitm);

	//presenting solution
      	cout<<"Time: "<<timediff<<" seconds"<<endl;
	cout<<"Best permutation found: [ ";
	for(int i=0;i<problem_size;i++) cout<<best[i]<<" ";
	cout<<"], with cost: "<<best_cost<<endl;
	interpret_perm(best);
	cout<<endl;
}

//assumes problem_size=10
void interpret_perm(int assgnmt[]) {	//prints to cout an object-room matching
	if(problem_size!=10) return; //this interpretation only valid for the room problem
	cout<<"This permutation matches ";
	char items[10]={'C','D','E','K','L','M','P','S','T','W'};
	char areas[10][3]={"BE","BW","DE","DN","DS","KA","UN","US","WA","WC"};
	for(int i=0;i<problem_size;i++) {
		cout<<items[assgnmt[i]]<<"-"<<areas[i]<<" ";
	}
	cout<<endl;
}

void random_search(void) {
	cout<<"Running random-restart local search..."<<endl;
	for(int i=0;i<10;i++) {

		//these statements form the core of the metaheuristic
		rand_assgnmt(curr);
		local_search();

		new_best();
	}
}

void exhaustive_search(void) {
	cout<<"Running an exhaustive search..."<<endl;
	for(int i=0;i<problem_size;i++) curr[i]=i;
	std::vector<int> vcurr(curr,curr+problem_size);
	do{
		for(int i=0;i<problem_size;i++)curr[i]=vcurr[i];
		new_best();
	}while(std::next_permutation(vcurr.begin(),vcurr.end()));
}

void create_stats(void) {
	cout<<"Creating statistics..."<<endl;
	ofstream file("stats.txt");
	file<<"Format:"<<endl;
	file<<"Cost best_neighbor_delta steps_to_local_optimum ";
	file<<"local_optimum_cost local_optimum"<<endl;

	int num_optimal=0;
	const int iterations=10000;
	for(int i=0;i<iterations;i++) {
		rand_assgnmt(curr);

		int cost_top=cost(curr);
		int area1,area2;
		file<<cost_top<<" "<<best_neighbor(area1,area2)<<" ";

		int num_steps=0;	//a version of local search that tracks steps
		while(best_neighbor(area1,area2)<0) {
			num_steps++;
			swap(curr[area1],curr[area2]);
		}
		file<<num_steps<<" ";

		int cost_bottom=cost(curr);
		file<<cost_bottom<<" ";
		for(int i=0;i<problem_size;i++) file<<curr[i]<<",";
		file<<endl;

		if(best_cost==cost_bottom) num_optimal++;
	}

	if(num_optimal>0.1*iterations) cout<<"This problem is easy: ";
	cout<<"With probability "<<((float)num_optimal/(float)iterations);
	cout<<" ("<<num_optimal<<"/"<<iterations<<") ";
	cout<<"a local search will end up at the optimum."<<endl;

	cout<<"Statistical raw data can be found in stats.txt"<<endl;

	file.close();
}

void new_best(void) {      //checking if new best
       int curr_cost=cost(curr);
       if(curr_cost<best_cost) {
		best_cost=curr_cost;
		for(int i=0;i<problem_size;i++) best[i]=curr[i];
       }
}

void local_search(void) {
	int area1,area2;
	while(best_neighbor(area1,area2)<0) swap(curr[area1],curr[area2]);
}

int best_neighbor(int& area1,int& area2) {
	int delta=0;
	for(int i=0;i<problem_size;i++) {
		for(int j=0;j<i;j++) {
		    int test=delta_cost(curr,i,j);
		    if(test<delta) {
			delta=test;
			area1=i;
			area2=j;
		    }
		}
	}
	return delta;
}

void tabu_search(void){
	cout<<"Running tabu search..."<<endl;
	int area1;
        int area2;
	for(unsigned int i=0;i<tabu_iterations;i++) {
	  int i=rand()%problem_size;
	  int j=rand()%(problem_size-i);
	  int test=delta_cost(curr,i,j);
	  if(test<0){
	     area1=i;
	     area2=j;
	     swap(curr[area1],curr[area2]);
             int curr_cost=cost(curr);
	     new_best();   //checks if a new best incumbent solution found
	     if (in_tabulist(curr) && curr_cost > best_cost){    //as this a min problem
                 swap(curr[area2],curr[area1]);
                 continue;
                 }
             else tabu_list(curr);
          }
	}
}

void tabu_list(int assignment[]){
        vector<int> temp(assignment,assignment+10);
        tabulist.push_back(temp);
        if(tabulist.size() >= tabu_size_list)  tabulist.pop_front();
        //if tabulist.size()<size_list then we let the list grow
}

bool in_tabulist(int assignment[]){
        vector<int> temp(assignment,assignment+10);
        for(unsigned int x=0;x<tabulist.size();x++) {
                if(temp==tabulist[x]) return true;
        }
        return false;
}

void simulated_annealing(void){
	cout<<"Running simulated annealing..."<<endl;
        float random_num;
	float probability;
	int area1;
	int area2;
    	for(unsigned int k=0;k<temp_steps;k++){
        for(unsigned int s=0;s<num_same_temp;s++){
	  	int i=rand()%problem_size;
	  	int j=rand()%(problem_size-i);
	  	int test=delta_cost(curr,i,j);
	  	if(test<0){
	     		area1=i;
	     		area2=j;
	     		swap(curr[area1],curr[area2]);
	     		new_best();   //checks if a new best incumbent solution found
	     	}
	  	else {
	    		random_num=(float)rand()/(float)RAND_MAX;
        		probability=exp(-test/T);
	    		if(random_num<probability)
	       			{
		 		area1=i;
		 		area2=j;
		 		swap(curr[area1],curr[area2]);
		 		new_best();   //checks if a new best incumbent solution found
                        }
	       }
	}
	T*=r;
	}
}

vector<pop_member > population(pop_size,pop_member(0));  //initializes population

pop_member::pop_member(int iteration)
	:created(iteration) {
	rand_assgnmt(curr);
        permutation=vector<int>(curr,curr+problem_size);
        p_cost=cost(curr);
        new_best(); //checks if new best
}

pop_member::pop_member(pop_member &m,pop_member &f,int iteration) //crossover constructor
	:created(iteration) {
        vector<bool> used(10,false); //keeping track of repeats
        for(int i=0;i<problem_size;i++) {
		if(rand()%2==0) curr[i]=m.permutation[i];
                else curr[i]=f.permutation[i];
                if((float)rand()/(float)RAND_MAX<mutation_prob) curr[i]=rand()%problem_size; //mutation with 5% prob.
                while(used[curr[i]]) curr[i]=rand()%problem_size;
                used[curr[i]]=true;
        }
        permutation=vector<int>(curr,curr+problem_size);
        p_cost=cost(curr);
        new_best();    //checks if new best (means we don't need to do this check in the GA loop
}

void genetic_algorithm(void){
        cout<<"Running genetic algorithm..."<<endl;
        for(unsigned int f=0;f<GA_iterations; f++){
                cull_pop();//function that assigns survival probabilities & deletes those that do not make it
                sort(population.begin(),population.end(),stupider);
                for(int i=0;i<pop_size/10;i++) {
                	population.push_back(pop_member(population[0],population[1],f));
                	population.push_back(pop_member(population[2],population[3],f));
                }
       }
}

inline bool stupider(pop_member &left,pop_member &right) {return left.p_cost>right.p_cost;}
inline bool older(pop_member &left,pop_member &right) {return left.created<right.created;}

void cull_pop(void){
            sort(population.begin(),population.end(),stupider);
            for(int i=0;i<pop_size/10;i++) population.pop_back();   //culls 10%
            sort(population.begin(),population.end(),older);
            for(int i=0;i<pop_size/10;i++) population.pop_back();  //culls 10%
           //get rid of stupidest and oldest
}