#include <iostream>
#include <fstream>
#include <cstdlib>
#include <cassert>
#include <ctime>
#include <cmath>
#include <climits>

using namespace std;

long int **read_matrix(ifstream &input, long int n);
double ran01(long int *idum);
long int *generate_random_vector(long int n);
unsigned long int compute_evaluation_function(long int n, long int **d, long int **f, long int *p);
void first_2_opt_symmetric (long int n, long int **d, long int **f, long int *p);
void swap(long int pos_1, long int pos_2, long int *);
void perturbation(long int pert_size, long int n, long int **d, long int **f, long int *p);
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 n;
    long int **d = NULL;
    long int **f = NULL;
    long int *p = NULL;
    long int initial_eval;
    long int *ILS_best_p = NULL;
    long int new_vs_current;
    new_vs_current = INT_MAX;
    long int current_best_eval;
    current_best_eval = INT_MAX;
    long int new_current_eval;
    long int count;
    count = 0;
    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";

    ILS_best_p = new long int[n];

	infile >> optimum;
    cout << "Optimum solution: " << optimum << "\n";
    d=read_matrix(infile, n);
    f=read_matrix(infile, n);
    cout << endl;

    p = generate_random_vector(n);
    compute_evaluation_function(n, d, f, p);
    cout << endl;
    first_2_opt_symmetric (n, d, f, p);
    cout << endl;
    initial_eval = compute_evaluation_function(n, d, f, p);
    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;
        perturbation(3, n, d, f, p);
        cout << endl << endl;

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

        first_2_opt_symmetric (n, d, f, p);

        new_current_eval = compute_evaluation_function(n, d, f, p);

        cout << endl;

        cout << " Vector after local search = ";
        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[k] = p[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[k] = p[k];
            }
        }

        if (current_best_eval == new_current_eval && new_current_eval == optimum)
        {
           	count += 1;
		}
    }
    cout << endl << "current best eval = " << current_best_eval << endl;

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

    cout << endl;
    cout << "Number of times achieving optimality = " << count << endl;

    for ( i = 0 ; i < n ; i++ )
	{
      delete [] d[i];
      delete [] f[i];
    }
    delete [] d;
    delete [] f;
    delete [] p;
    delete [] ILS_best_p;

    return 0;
}

long int **read_matrix(ifstream &input, long int n)
{

    /*
    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;
}

long int *generate_random_vector(long int n)
{

    /*
    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;
    long int *v = NULL;
    srand(static_cast<size_t>(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;
}

unsigned long int compute_evaluation_function(long int n, long int **d, long int **f, long int *p)
{

    /*
    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 (long int n, long int **d, long int **f, long int *p)
{
    /*
    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;
    long int *x;
    improvement = 1;
    x = generate_random_vector(n);

    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, 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, long int n, long int **d, long int **f, long int *p)
{
    /*
    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
    */

    long int *hold_value = new long int[n];         //allocate memory for items to be chosen in move/*
    int *chosen = new int[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(n, d, f, p);
    new_eval = actual_solution_value;
    current_eval = new_eval;

    delete [] hold_value;
    delete [] chosen;
}