//******************************************************
//  Copyright 1999 University of Paderborn
//
//  FILE : SA_Problem.cc
//  
//  includes         : Problem.h, parrandom.h
//                            
//
//  This is an example to show how to use the interface of 
//  the ParSA-Library. 
//  The target is to find the (absolute) minimum of an 
//  instance of the quadratic assignment problem
//  For the ParSA-Library that means:
//     class PROBLEM represents the QAP
//     class SOLUTION represents a certain permutation.
//     The COST-FUNCTION is the value of 
//       ___   ___
//       \     \ 
//       /     /    a     * b
//       ---   ---    ij      p(i)p(j)
//        i     j
//
//     The NEIGHBORHOOD-OPERATION is the random choice of
//       two indices which permutation values will be swapped
//
//  authors:
//    Karsten Klohs 
//  contact:
//    parsa@uni-paderborn.de
//  last update:
//    17.10.1999
//*********************************************************
#include <parrandom.h>     	                 // randomnumbergenerator
#include "SA_QAPProblem.h"
// -- EXAMPLE
// -- for reading the initializationfile
//#include <fstream.h>
//#include <stdlib.h>
//#include <iomanip.h>
#include <fstream>
#include <iomanip>
//#include <constants.h>

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

// ===============> CLASS SOLUTION <=================

SA_Solution::SA_Solution() {
  p = NULL;
}
// -----------------------------------------------------
// SA_Solution::Copy
//
//         Copies the given SA_Solution in current instance
//
// -----------------------------------------------------
void SA_Solution::Copy(SA_Solution & source) 
{ 
  n = source.n;
  if (p == NULL) {
    p = new int[n];
  }
  for (int k = 0; k < n; k++) {     // copy the permutation array 
    p[k] = source.p[k];
  }
  i = source.i;
  j = source.j;
  oldcost = source.oldcost;
  newcost = source.newcost;
}

std::ostream & operator <<(std::ostream &out, SA_Solution &sol) 
{
  out.width(5);
  out << setiosflags(std::ios::right) << sol.n << " " << sol.newcost << std::endl;
  for (int k = 0; k < sol.n; k++) {
    out.width(5);
    out << setiosflags(std::ios::right) << sol.p[k];
    if ((k+1)%20 == 0) {
      out << std::endl;
    }
  }
  return(out); 
}


// ===========> CLASS MOVE (optional) <==============
SA_Move::SA_Move(){}
SA_Move::~SA_Move(){}
SA_Move::SA_Move(SA_Problem &prb){}

// ================> CLASS PROBLEM <=================


//******************************************************
// class PROBLEM: 
//
//       constructor / destructor   
//
//****************************************************** 

// -----------------------------------------------------
// constructor
// -----------------------------------------------------
// -- EXAMPLE
// -- The default constructor initializes the instance 
// -- with n = 0 and empty matrix arrays
SA_Problem::SA_Problem() 
{ 
  n = 0;
//  A = NULL; B = NULL;
  globaln = 0;
  localn = 0;
  void randomize();
}

// -----------------------------------------------------
// destructor
// -----------------------------------------------------
SA_Problem::~SA_Problem()
{
}

void SA_Problem::EvaluateOptions(int *argc,char ***argv) {

}

//******************************************************
// class PROBLEM
//
//       methods for SEQUENTIAL SA
//  
//****************************************************** 



// -----------------------------------------------------
// SA_Problem::ReadProblemData
//    
//        param
//           - filename: name of the input file
//        return 
//           TRUE / FALSE  
//    
// -----------------------------------------------------
// EXAMPLE
// Reads the problem data from file "filename".
int SA_Problem::ReadProblemData(char *filename)
{   
  std::ifstream initfile(filename);

  std::cout << "Reading datafile " << filename << std::endl;
  if (initfile) {
    initfile >> n;
    std::cout << "Matrix size: " << n << std::endl;
  } else {
    std::cerr << "Error reading matrix size!" << std::endl;
    exit(1);
  }
  // Memory allocation
  // A = (int***) new MELEMTYPE[n*n];
//  A = malloc(n*n*sizeof(MELEMTYPE));
//  if (A == NULL) {
//    std::cout << "Allocated " << n*n << " * " << sizeof(MELEMTYPE) << " = " << n*n*sizeof(MELEMTYPE) << " bytes of memory for matrix A" << std::endl;
//  } else {
//    std::cerr << "Memory allocation for matrix A failed!" << std::endl;
//    exit(1);
//  }  
//  B = (int***) new MELEMTYPE[n*n];
//  if (B == NULL) {
//    std::cout << "Allocated " << n*n << " * " << sizeof(MELEMTYPE) << " = " << n*n*sizeof(MELEMTYPE) << " bytes of memory for matrix B" << std::endl;
//  } else {
//    std::cerr << "Memory allocation for matrix B failed!" << std::endl;
//    exit(1);
//  }
  
  MELEMTYPE inelem;
  for (int k = 0; k < n; k++) {
    for (int l = 0; l < n; l++) {
      if (initfile >> inelem) {
        A[k][l] = inelem;
      } else {
        std::cerr << "Error while reading element a" << k << "," << l << " !" << std::endl;
        exit(1); 
      }
    }
  } 
  for (int k = 0; k < n; k++) {
    for (int l = 0; l < n; l++) {
      if (initfile >> inelem) {
        B[k][l] = inelem;
      } else {
        std::cerr << "Error while reading element b" << k << "," << l << " !" << std::endl;
        exit(1); 
      }
    }
  }

  // computing number of feasible solutions ( n! )
  // compute them rather exactly if n < 30 otherwise
  // set global n to 1E30
  if (n <= 30) {
    globaln = 1;
    for (int k = 1; k <= n; k++) {
      globaln *= globaln;
    }
  } else {
    globaln = 1E30;
  }
  // computing number of neighbors ( n over 2 )
  localn = (n*n - n) / 2;

  return 1;
}
	
// -----------------------------------------------------
// SA_Problem::CreateSolution()
//
//        return
//           a new solution instance
//
// -----------------------------------------------------

SA_Solution * SA_Problem::CreateSolution()
{ 
   SA_Solution *newsol = new SA_Solution;
   
   newsol->n = n;
   newsol->p = new MELEMTYPE[n];
   return(newsol);
}

// -----------------------------------------------------
// SA_Problem::GetInitialSolution()
//
//        calculates an initial solution for the
//        SA-Algorithm (deterministic)
// 
// -----------------------------------------------------
// -- EXAMPLE
// -- Choose the identity function
void SA_Problem::GetInitialSolution(SA_Solution & sol)
{ 
  for (int k = 0; k < n; k++) {
    sol.p[k] = k;
  }
  sol.i = 0; sol.j = 0;
  sol.oldcost = 0;
  for (int k = 0; k < n; k++) {
    for (int l = 0; l < n; l++) {
      sol.oldcost += A[k][l] * B[sol.p[k]][sol.p[l]]; 
    }
  }
  sol.newcost = sol.oldcost;

}


// -----------------------------------------------------
// SA_Problem::GetRandomSolution()
//
//        generates a random initial solution for the
//        SA-Algorithm (non-deterministic)
// 
// -----------------------------------------------------
// -- EXAMPLE
// -- not implemented, using the identity function
void SA_Problem::GetRandomSolution(SA_Solution & sol)
{
  GetInitialSolution(sol);
}

// -----------------------------------------------------
// SA_Problem::GetCost()
//
//        param
//           solution: solution which cost has to be 
//                     calculated 
//        return
//           cost of given solution
// 
// -----------------------------------------------------
// -- EXAMPLE
// -- as the cost calculation requires O(n^2) 
// -- multiplikations it is always done once when a 
// -- neighbor is created

float SA_Problem::GetCost(SA_Solution & sol)
{
  return((float)sol.newcost);
}

// -----------------------------------------------------
// SA_Problem::GetNeighbor()
//
//        generates a reversible random change in the given 
//        solution. The range of change can be modified 
//        depending on rel_temp.
//
//        param
//           rel_temp: relative temperature - ranges from
//                     1 (=high) to 0 (=low)  
// 
// -----------------------------------------------------
// -- EXAMPLE
// -- two indices are choosen which permutation values are 
// -- changed
// -- Cost calculation is done by investigating only elements
// -- which have on or both of the new choosen swap indices.
// -- Old products are subtracted and new ones are added to 
// -- the old cost. To save some multiplications this is done
// -- a little bit more tricky than the straight forward 
// -- implementation.

void SA_Problem::GetNeighbor(SA_Solution & sol, float rel_temp)
{
  // first of all, forget the predecessor of current solution

  //  sol.oldcost = sol.newcost;
  //  sol.i = sol.j;

  long randomnumber;
  randomnumber = rand();

  int i = randomnumber % n;
  randomnumber = rand();
  int j = randomnumber %n;
  int *p = sol.p;
  sol.i = i;
  sol.j = j;

  // updating the permutation
  swap(p[i],p[j]);
  sol.newcost = 0;
  for (int k = 0; k < n; k++) {
    for (int l = 0; l < n; l++) {
      sol.newcost += A[k][l]*B[p[k]][p[l]];
    }
  }

  // determining the new cost value
  //  sol.newcost = sol.oldcost;
  //  for (int k = 0; k < n; k++) {      // for all indices
  //    // k represents the "other" index in a pair which second index is one of i or j
  //    if ((k != i) && (k != j)) {         
  //      sol.newcost -= (A[k][i] - A[k][j]) * (B[p[k]][p[j]] - B[p[k]][p[i]]); // All elements with k in first position
  //      sol.newcost -= (A[i][k] - A[j][k]) * (B[p[j]][p[k]] - B[p[i]][p[k]]); // All elements with k in second position
  //    }
  //  }
  //  sol.newcost -= (A[i,i] - A[j,j]) * (B[p[j]][p[j]] - B[p[i]][p[i]]);
  //  sol.newcost -= (A[i,j] - A[j,i]) * (B[p[j]][p[i]] - B[p[i]][p[j]]);
}

// -----------------------------------------------------
// SA_Problem::ResetSolution()
//
//        resets the change of solution made by 
//        GetNeighbor()
//
// -----------------------------------------------------
// -- EXAMPLE
// -- restore the old permutation, forget the new cost value
// -- and set i = j

void SA_Problem::ResetSolution(SA_Solution & sol)
{
  swap(sol.p[sol.i],sol.p[sol.j]);
  sol.j = sol.i;
  sol.newcost = sol.oldcost;
}

// -----------------------------------------------------
// SA_Problem::UpdateSolution()
//
//        makes the change solution made by GetNeighbor
//        permanent
//
// -----------------------------------------------------
// -- EXAMPLE
// -- makes the new permutation permanent by forgetting the
// -- choosen swapindices and the old cost value
	
void SA_Problem::UpdateSolution(SA_Solution & sol)
{
  sol.i = sol.j;
  sol.oldcost = sol.newcost;
}


// -----------------------------------------------------
// SA_Problem::GetLocalN
//    
//        return  
//          number of solution neighbours (approximately)  
//    
// -----------------------------------------------------

float SA_Problem::GetLocalN()
{ 
  return((float) localn);
}

// -----------------------------------------------------
// SA_Problem::GetGlobalN
//    
//        return  
//          number of siutable solutions  
//    
// -----------------------------------------------------

float SA_Problem::GetGlobalN()
{ 
  return((float) globaln);
}

// -----------------------------------------------------
// SA_Problem::CopySolution()
//
//        copies a solution from source to dest
// 
// -----------------------------------------------------

void SA_Problem::Copy(SA_Solution & source, SA_Solution & dest)
{
  dest.n = source.n;
  dest.i = source.i;
  dest.j = source.j;
  dest.newcost = source.newcost;
  dest.oldcost = source.oldcost;
  for (int k = 0; k < n; k++) {
    dest.p[k] = source.p[k];
  }
}	


//******************************************************
// class PROBLEM: 
//
//       methods for PARRALLEL SA   
//
//****************************************************** 

// -----------------------------------------------------
// SA_Problem::OutputSolution()      (STREAM COMM)
//
//        writes a solution in stream
// 
// -----------------------------------------------------

void SA_Problem::OutputSolution(std::ostream& out,SA_Solution & sol)
{ 
  out << sol.n << std::endl;
  for (int k = 0; k < n; k++) {
    out << sol.p[k] << std::endl;
  }
  out << sol.i << std::endl;
  out << sol.j << std::endl;
  out << sol.oldcost << std::endl;
  out << sol.newcost << std::endl;
}

// -----------------------------------------------------
// SA_Problem::OutputSolution()      (STREAM COMM)
//
//        reads a solution from stream
// 
// -----------------------------------------------------

int SA_Problem::InputSolution(std::istream& in,SA_Solution & sol)
{
  in >> sol.n;
  for (int k = 0; k < n; k++) {
    in >> sol.p[k];
  }
  in >> sol.i; 
  in >> sol.j;
  in >> sol.oldcost;
  in >> sol.newcost;
}


#ifdef MPI
// -----------------------------------------------------
// SA_Problem::BcastSolution()        (MPI COMM)
//
//        broadcasts a solution ???
// 
// -----------------------------------------------------

int  SA_Problem::BcastSolution(SA_Solution &sol, int root, MPI_Comm comm, int me_comm)
{ 

}

#endif

//******************************************************
// class PROBLEM: 
//
//       ADDITIONAL methods  
//
//****************************************************** 


// -----------------------------------------------------
// SA_Problem::CreateMove()           (MOVE_USE)         
//
//        return
//           new instance of SA_Move
// 
// -----------------------------------------------------

SA_Move * SA_Problem::CreateMove()
{

}


// -----------------------------------------------------
// SA_Problem::ExtractMove()           (MOVE_USE)         
//
//        extracts the made move from sol and stores it 
//        in move m.
// 
// -----------------------------------------------------

int	SA_Problem::ExtractMove(SA_Solution & sol, SA_Move & m)
{
}


// -----------------------------------------------------
// SA_Problem::ApplyMove()           (MOVE_USE)         
//
//      changes Solution sol according to move m
//
// -----------------------------------------------------

int SA_Problem::ApplyMove(SA_Solution & sol, SA_Move & m)	
{
}


// -----------------------------------------------------
// SA_Problem::OutputMove()           (MOVE_USE, STREAM COMM)         
//
//      writes a move in stream
//
// -----------------------------------------------------

void SA_Problem::OutputMove(std::ostream & out, SA_Move & m)
{
}

// -----------------------------------------------------
// SA_Problem::InputApplyMoves()      (MOVE_USE, STREAM COMM)         
//
//      reads a number of moves from stream and changes 
//      solution accordingly
//
//      return
//         number of steps
//
// -----------------------------------------------------

int SA_Problem::InputApplyMoves(std::istream & in, SA_Solution & sol)												
{
}

#ifdef MPI
// -----------------------------------------------------
// SA_Problem::BcastApplyMove()      (MOVE_USE, MPI COMM)         
//
//      ???
//
//      return
//         number of steps
//
// -----------------------------------------------------

// Unter der Voraussetzung, dass alle Loesungen im comm gleich bis auf die von root sind, 
// und die letzte sich durch genau einen Uebergang unterscheidet,
// bringe alle Loesungen auf den gleichen Stand;
// Nicht vergessen : UpdateSolution() auf die Loesung von root anzuwenden 

int	SA_Problem::BcastApplyMove(SA_Solution &sol, int root, MPI_Comm comm, int me_comm)
{
}

#endif

// -----------------------------------------------------
// SA_Problem::OutputSolution()      
//
//      writes a solution in a file
//
// -----------------------------------------------------


int SA_Problem::OutputSolution(char * filename,SA_Solution & sol)
{
  std::fstream outfile(filename,std::ios::out);

  if (!outfile) {
    std::cerr << "Error opening file: " << filename << std::endl;
    exit(1);
  } else {
    outfile << sol.n << "  ";
    outfile << sol.newcost;
    outfile << "\n";
    for (int k = 0; k < n; k++) {
      outfile << sol.p[k] << " ";
    }
    outfile.close();
  }
}

// -----------------------------------------------------
// SA_Problem::InputSolution()      
//
//      reads a solution from a file
//
// -----------------------------------------------------
	
int SA_Problem::InputSolution(char * filename,SA_Solution & sol)
{
}
											

// -----------------------------------------------------
// SA_Problem::Postprocessing()      
//
//      ???
//
// -----------------------------------------------------

// Methode Postprocessing
// wird nach dem Annealing mit der besten gefundenen Loesung zwecks
// Nachbearbeitung aufgerufen;

void	SA_Problem::Postprocessing(SA_Solution & sol)
{
}

// -----------------------------------------------------
// SA_Problem::Reoptimize()      
//
//      ???
//
// -----------------------------------------------------

// Methode Reoptimize
// wird nach dem Annealing mit der besten gefundenen Loesung zwecks
// Nachbearbeitung aufgerufen;

int	SA_Problem::Reoptimize(SA_Solution & sol)
{
}