#include <iostream>
#include <fstream>
#include <cstdlib>
#include <algorithm>
#include <cassert>
#include <ctime>
#include <cmath>
#include <vector>

using namespace std;

long int **read_matrix(ifstream &input);
double ran01(long int *idum);
vector <long int> generate_random_vector();
unsigned long int compute_evaluation_function();
void first_2_opt_symmetric ();
void swap(long int pos_1, long int pos_2, vector <long int>);
void perturbation(long int pert_size);
vector <long int> p;//the vector generated from generate_random_vector()
long int **d;//one of the matrices read from file
long int **f;//one of the matrices read from file
long int n;
long int actual_solution_value;
long int new_eval;
long int current_eval;

#define IA 16807
#define IM 2147483647
#define AM (1.0/IM)
#define IQ 127773
#define IR 2836

int main()
{
    string name;
    long int optimum;
    ifstream infile;
    long int initial_eval;
    vector <long int> ILS_best_p;
    long int new_vs_current;
    new_vs_current = INT_MAX;
    long int current_best_eval;
    current_best_eval = INT_MAX;
    //long int array[n];
    long int new_current_eval;
    long int best_eval;
    long int improvement;
    long int i;
    long int j;
    long int k;

    infile.open("nug12.dat");
    infile >> name;
    cout << "Name: " << name << "\n";

    infile >> n;
    cout << "Rows and Columns: " << n << "\n";

    infile >> optimum;
    cout << "Optimum solution: " << optimum << "\n";

    d=read_matrix(infile);
    f=read_matrix(infile);
    cout << endl;

    p=generate_random_vector();
    for (i = 0; i < n; i++)
    {
    	cout << p[i] << " ";
	}

    compute_evaluation_function();
    cout << endl;
    first_2_opt_symmetric ();
    cout << endl;
    initial_eval = compute_evaluation_function();
    cout << endl << endl;

    for (i = 0; i < 10000; i++) //ILS loop
    {
        if (new_eval < current_eval)
        {
            new_vs_current = new_eval;
            cout << "iteration 1st = " << i+1 << " " << endl;
        }
        else if (current_eval < new_eval)
        {
            new_vs_current = current_eval;
            cout << "iteration 2nd = " << i+1 << " " << endl;
        }
        cout << "iteration = " << i+1 << " " << endl;
        cout << " Vector before perturbation = ";
        for (j = 0; j < n; j++)
        {
            cout << p[j] << " ";
        }
        cout << endl << endl;
        cout << endl << "current best eval = " << current_best_eval << endl;
        perturbation(3);
        cout << endl << endl;
        cout << endl << "current best eval = " << current_best_eval << endl;

        cout << " Vector after perturbation = ";
        for (j = 0; j < n; j++)
        {
            cout << p[j] << " ";
        }
        cout << endl << endl;

        cout << endl << "current best eval = " << current_best_eval << endl;

        first_2_opt_symmetric ();

        new_current_eval = compute_evaluation_function();
        cout << endl << "current best eval = " << current_best_eval << endl;

        cout << endl;

        cout << " Vector after local search = ";
        for (j = 0; j < n; j++)
        {
            cout << p[j] << " ";
        }
        cout << endl << endl;
        cout << endl << "current best eval = " << current_best_eval << endl;

        cout << " Vector p = ";
        for (j = 0; j < n; j++)
        {
            cout << p[j] << " ";
        }
        cout << endl << endl;

        if (new_current_eval < current_best_eval)
        {
            new_eval = new_current_eval;
            cout << "iteration 3rd = " << i+1 << " " << endl;
        }
        else if (current_best_eval <= new_current_eval)
        {
            new_eval = current_best_eval;
            cout << "iteration 4th = " << i+1 << " " << endl;
        }

        if (new_eval < new_vs_current)
        {
            current_best_eval = new_eval;
            cout << "iteration 5th = " << i+1 << " " << endl;
            //for(k = 0; k < n; k++)
            //{
            	ILS_best_p = p;
			//}
            //for(k = 0; k < n; k++)
            //{
            //	cout << p[k];
			//}
            //for(k = 0; k < n; k++)
            //{
            //	cout << array[k];
			//}
		//	for (k = 0; k < n; k++)
		//	{
		//		ILS_best_p[k] = array[k];
		//	}
        }
        else if (new_vs_current < new_eval)
        {
            current_best_eval = new_vs_current;
            cout << "iteration 6th = " << i+1 << " " << endl;
            for(k = 0; k < n; k++)
            {
				ILS_best_p = p;
			}
            //for(k = 0; k < n; k++)
            //{
            //	cout << p[k];
			//}
            //for(k = 0; k < n; k++)
            //{
            //	cout << array[k];
			//}
			//for (k = 0; k < n; k++)
			//{
			//	ILS_best_p[k] = array[k];
			//}
        }

        cout << endl << "current best eval = " << current_best_eval << endl;
    }
    cout << endl << "current best eval = " << current_best_eval << endl;

    cout << "ILS best vector = ";
    for (i = 0; i < n; i++)
    {
        cout << ILS_best_p[i] << " ";
    }

    return 0;
}

long int **read_matrix(ifstream &input)
{

    /*
    FUNCTION:       read instance matrices
    INPUT:          instance file name, instance size
    OUTPUT:         d-matrix and f-matrix
    (SIDE)EFFECTS:  none
    COMMENTS:       read the distance matrix and flow matrix
    */


    long int i, j;
    long int **matrix = new long int *[n];
    for (i = 0; i < n; i++)
    {
        matrix[i] = new long int[n];
        for (j = 0; j < n; j++)
        {
            if( !(input >> matrix[i][j]) )
            {
                cerr << "Error reading at " << i << j << endl;
            	exit(1);
            }
        }
    }
    return matrix;
}

double ran01(long int *idum)
{

    /*
    FUNCTION:      returns a pseudo-random number
    INPUT:         a pointer to the seed variable
    OUTPUT:        a pseudo-random number uniformly distributed in [0,1]
    (SIDE)EFFECTS: changes the value of seed
    */


    long k;
    double ans;

    k =(*idum)/IQ;
    *idum = IA * (*idum - k * IQ) - IR * k;
    if (*idum < 0 )
    {
        *idum += IM;
    }
    ans = AM * (*idum);
    return ans;
}

vector <long int> generate_random_vector()
{

    /*
    FUNCTION:      generates a random vector
    INPUT:         vector dimension
    OUTPUT:        returns pointer to vector, free memory when vector is not needed anymore
    (SIDE)EFFECTS: none
    */


    long int  i, j, help;
    vector <long int> v;
    srand(time(0));
    long int seed=rand();

    //v = new long int[ n ];

    for ( i = 0 ; i < n; i++ )
    {
        v[i] = i;
    }

    for ( i = 0 ; i < n-1 ; i++)
    {
        j = (long int) ( ran01( &seed ) * (n - i));
        assert( i + j < n );
        help = v[i];
        v[i] = v[i+j];
        v[i+j] = help;
    }
    return v;
    //delete []v;
}

unsigned long int compute_evaluation_function()
{

    /*
    FUNCTION:      compute evaluation function
    INPUT:         pointer to solution
    OUTPUT:        evaluation function
    (SIDE)EFFECTS: none
    COMMENTS:      none
    */

    long int i, j;
    unsigned long obj_f_value;

    obj_f_value = 0;
    for(i = 0; i < n; i++)
    {
        for(j = 0; j < n; j++)
        {
            obj_f_value += d[i][j] * f[p[i]][p[j]];
        }
    }
    cout << obj_f_value;
    return obj_f_value;
}

void first_2_opt_symmetric ()
{
    /*
    FUNCTION:      first improvement 2-opt local search for symmetric instances
    INPUT:         pointer to some initial assignment
    OUTPUT:        none
    (SIDE)EFFECTS: initial assignment is locally optimized, dlb is used to increase the search speed, the dlb value is changed to false if item change it's location during perturbation
    COMMENT:       neighborhood is scanned in random order, use of register for faster computation
    */

    long int improvement = 1;
    long int improve_item = 0;
    long int u, v, i, j, k;
    long int tmp;
    vector <long int> x;

    improvement = 1;
    x = generate_random_vector();

    while (improvement)
    {
        improvement = 0;
        for ( i = 0 ; i < n ; i++ )
        {
            u = x[i];
            improve_item = 0;
            for ( j = 0 ; j < n ; j++ )
            {
                v = x[j];
                if (u == v)
                    continue;
                tmp = 0;
                for ( k = 0 ; k < n ; k++ )
                {
                    if ( (k != u) && (k != v) )
                    {
                        tmp += (d[u][k] - d[v][k]) * (f[p[v]][p[k]] - f[p[u]][p[k]]);
                    }
                }
                if (tmp < 0)
                {
                    improvement = 1;
                    improve_item = 1;
                    swap (u, v, p);
                }
            }
        }
    }
    //delete []x;
}

void swap(long int pos_1, long int pos_2, vector <long int> p)
{

    /*
    FUNCTION:      swap position of two items in a vector
    INPUT:         position 1, position 2, pointer to vector
    OUTPUT:        none
    (SIDE)EFFECTS: none
    COMMENTS:      none
    */

    long int tmp;

    tmp = p[pos_1];
    p[pos_1] = p[pos_2];
    p[pos_2] = tmp;
}

void perturbation(long int pert_size)
{
    /*
    FUNCTION:         apply random perturbation
    INPUT:            pointer to solution, perturbation scheme, perturbation strength, change in perturbation, end perturbation size
    OUTPUT:           none
    (SIDE)EFFECTS:    new solution formed from old solution
    COMMENTS:         none
    */

    vector <long int> hold_value;         //allocate memory for items to be chosen in move/*
    int chosen[pert_size];          //*allocate temporary memory to determine items to be moved */*
    long int i;

    for(i = 0; i < n; i++)
    {
        hold_value[i] = p[i];               //initialize hold_value with the same value from local search/*
    }
    int j = n - 1;
    int rand_number;
    int rand_no;

    for(i = 1; i <= pert_size; i++)
    {
        rand_number = rand();
        rand_no = rand_number % j;              //choose random number from 1 to size - 1/
        chosen[i - 1] = rand_no + i;            //copy the value to chosen[i]
        j--;
    }
    long int temp;

    for(i = 0; i < pert_size; i++)
    {
        temp = chosen[i];
        swap(i, temp, hold_value);              //swap chosen[i] and hold_value[i]; n first item in hold_value[i] is the index array that will be included in perturbation/
    }

    long int temp1;
    long int temp2;
    temp1 = p[hold_value[1]];//
    temp2 = p[hold_value[0]];

    for(i = 0; i < pert_size - 1; i++)
    {
        p[hold_value[i + 1]] = temp2;
        temp2 = temp1;
        temp1 = p[hold_value[i + 2]];
    }

    p[hold_value[0]] = temp2;

    actual_solution_value = compute_evaluation_function();
    new_eval = actual_solution_value;
    current_eval = new_eval;
}

And this is the file I read from
Nug12
12
578
0 1 2 3 1 2 3 4 2 3 4 5
1 0 1 2 2 1 2 3 3 2 3 4
2 1 0 1 3 2 1 2 4 3 2 3
3 2 1 0 4 3 2 1 5 4 3 2
1 2 3 4 0 1 2 3 1 2 3 4
2 1 2 3 1 0 1 2 2 1 2 3
3 2 1 2 2 1 0 1 3 2 1 2
4 3 2 1 3 2 1 0 4 3 2 1
2 3 4 5 1 2 3 4 0 1 2 3
3 2 3 4 2 1 2 3 1 0 1 2
4 3 2 3 3 2 1 2 2 1 0 1
5 4 3 2 4 3 2 1 3 2 1 0

0  5  2  4  1  0  0  6  2  1  1  1
5  0  3  0  2  2  2  0  4  5  0  0
2  3  0  0  0  0  0  5  5  2  2  2
4  0  0  0  5  2  2 10  0  0  5  5
1  2  0  5  0 10  0  0  0  5  1  1
0  2  0  2 10  0  5  1  1  5  4  0
0  2  0  2  0  5  0 10  5  2  3  3
6  0  5 10  0  1 10  0  0  0  5  0
2  4  5  0  0  1  5  0  0  0 10 10
1  5  2  0  5  5  2  0  0  0  5  0
1  0  2  5  1  4  3  5 10  5  0  2
1  0  2  5  1  0  3  0 10  0  2  0