//******************************************************
//  Copyright 1999 University of Paderborn
//
//  FILE : SA_MIProblem.cc
//  
//  See the file SA_MIProblem.h for documentation on this class.
//
//  author:
//    Richard Jones
//  contact:
//    jonesrt@phys.uconn.edu
//  last update:
//    July 25, 2007
//*********************************************************

#include "SA_MIProblem.h"

#include <fstream>
#include <iomanip>
#include <string>
#include <sstream>
#include <values.h>

SA_Solution::SA_Solution()
 : fNdim(0), fCost(-1), fNewcost(-1), fLastmove(0),
   fRefa(0), fSig1a(0), fSig2a(0), fSig1ph(0), fSig2ph(0),
   fRefaMap(0), fSig1aMap(0), fSig2aMap(0), fSig1phMap(0), fSig2phMap(0)
{
}

SA_Solution::SA_Solution(const unsigned int ndim)
 : fNdim(ndim), fCost(-1), fNewcost(-1), fLastmove(0),
   fRefaMap(0), fSig1aMap(0), fSig2aMap(0), fSig1phMap(0), fSig2phMap(0)
{
   fRefa = new TModel(ndim);
   fSig1a = new TModel(ndim);
   fSig2a = new TModel(ndim);
   fSig1ph = new TModel(ndim);
   fSig2ph = new TModel(ndim);
}

SA_Solution::~SA_Solution()
{
   if (fRefa)
     delete fRefa;
   if (fSig1a)
      delete fSig1a;
   if (fSig2a)
      delete fSig2a;
   if (fSig1ph)
      delete fSig1ph;
   if (fSig2ph)
      delete fSig2ph;
   if (fRefaMap)
      delete fRefaMap;
   if (fSig1aMap)
      delete fSig1aMap;
   if (fSig2aMap)
      delete fSig2aMap;
   if (fSig1phMap)
      delete fSig1phMap;
   if (fSig2phMap)
      delete fSig2phMap;
   if (fLastmove)
      delete fLastmove;
} 

void SA_Solution::Copy(SA_Solution& source) 
{ 
   // Implements deep-copy sematics for SA_Solution

   if (fRefa) {
      delete fRefa;
      fRefa = 0;
   }
   if (fSig1a) {
      delete fSig1a;
      fSig1a = 0;
   }
   if (fSig2a) {
      delete fSig2a;
      fSig2a = 0;
   }
   if (fSig1ph) {
      delete fSig1ph;
      fSig1ph = 0;
   }
   if (fSig2ph) {
      delete fSig2ph;
      fSig2ph = 0;
   }
   if (fRefaMap) {
      delete fRefaMap;
      fRefaMap = 0;
   }
   if (fSig1aMap) {
      delete fSig1aMap;
      fSig1aMap = 0;
   }
   if (fSig2aMap) {
      delete fSig2aMap;
      fSig2aMap = 0;
   }
   if (fSig1phMap) {
      delete fSig1phMap;
      fSig1phMap = 0;
   }
   if (fSig2phMap) {
      delete fSig2phMap;
      fSig2phMap = 0;
   }
   if (fLastmove) {
      delete fLastmove;
      fLastmove = 0;
   }

   fNdim = source.fNdim;
   if (source.fRefa)
      fRefa = new TModel(*(source.fRefa));
   if (source.fSig1a)
      fSig1a = new TModel(*(source.fSig1a));
   if (source.fSig2a)
      fSig2a = new TModel(*(source.fSig2a));
   if (source.fSig1ph)
      fSig1ph = new TModel(*(source.fSig1ph));
   if (source.fSig2ph)
      fSig2ph = new TModel(*(source.fSig2ph));
   if (source.fRefaMap)
      fRefaMap = new TBitmap(*(source.fRefaMap));
   if (source.fSig1aMap)
      fSig1aMap = new TBitmap(*(source.fSig1aMap));
   if (source.fSig2aMap)
      fSig2aMap = new TBitmap(*(source.fSig2aMap));
   if (source.fSig1phMap)
      fSig1phMap = new TBitmap(*(source.fSig1phMap));
   if (source.fSig2phMap)
      fSig2phMap = new TBitmap(*(source.fSig2phMap));
 
   fCost = source.fCost;
   fNewcost = source.fNewcost;
   if (source.fLastmove) {
      fLastmove = new SA_Move(*(source.fLastmove));
   }
}

std::ostream& operator<<(std::ostream& out, SA_Solution& sol) 
{
   out << " " << sol.fNdim;
   out << " " << sol.fCost;
   out << " " << sol.fNewcost;
   if (sol.fNdim) {
      out << sol.fRefa;
      out << sol.fSig1a;
      out << sol.fSig1ph;
      out << sol.fSig2ph;
   }
   if (sol.fLastmove) {
      out << " 1" << *(sol.fLastmove);
   }
   else {
      out << " 0";
   }
   return out; 
}

std::istream& operator>>(std::istream& in, SA_Solution*& sol) 
{
   sol = new SA_Solution();
   in >> sol->fNdim >> sol->fCost >> sol->fNewcost;
   if (sol->fNdim) {
      in >> sol->fRefa;
      in >> sol->fSig1a;
      in >> sol->fSig1ph;
      in >> sol->fSig2ph;
   }
   int has_lastmove;
   in >> has_lastmove;
   if (has_lastmove) {
      in >> sol->fLastmove;
   }
   else {
      sol->fLastmove = 0;
   }
   return(in); 
}

SA_Move::SA_Move()
 : fModel_set(0), fModel(0), fBitmap(0)
{
}

SA_Move::~SA_Move()
{
   if (fModel)
      delete fModel;
   if (fBitmap)
      delete fBitmap;
}

SA_Move::SA_Move(SA_Problem& prob)
 : fModel_set(0)
{
   int ndim = prob.GetNdim();
   fModel = new TModel(ndim);
   fBitmap = 0;
}

SA_Move::SA_Move(const SA_Move& move)
{
   fModel_set = move.fModel_set;
   if (move.fModel) {
      fModel = new TModel(*(move.fModel));
   }
   else {
      fModel = 0;
   }
   if (move.fBitmap) {
      fBitmap = new TBitmap(*(move.fBitmap));
   }
   else {
      fBitmap = 0;
   }
}

std::ostream& operator<<(std::ostream& out, SA_Move& move) 
{
   out << " " << move.fModel_set;
   if (move.fModel) {
      out << " 1" << move.fModel;
   }
   else {
      out << " 0";
   }
   return out; 
}

std::istream& operator>>(std::istream& in, SA_Move*& move) 
{
   move = new SA_Move();
   in >> move->fModel_set;
   int has_fModel;
   in >> has_fModel;
   if (has_fModel) {
      in >> move->fModel;
   }
   else {
      move->fModel = 0;
   }
   return(in); 
}

SA_Problem::SA_Problem(int ndim) 
 : fDataImage(0), fModelImage(0), fNdim(ndim)
{ 
   fWeightSum = 0;
   for (int i=0; i<4; i++) {
      fNvar[i] = ndim;
      fWeightSum += fWeight[i] = ndim*ndim;
   }
}

SA_Problem::~SA_Problem()
{
   if (fDataImage)
      delete fDataImage;
   if (fModelImage)
      delete fModelImage;
}

void SA_Problem::EvaluateOptions(int *argc,char ***argv)
{
   // Picks through the command line for options relevant to SA_Problem

   for (int j=0; j<*argc; j++) {
      if (strstr((*argv)[j],"--ndim") == (*argv)[j]) {
         fNdim = atoi((*argv)[++j]);
      }
      else if (strstr((*argv)[j],"--nvar") == (*argv)[j]) {
         fWeightSum = 0;
         std::istringstream str((*argv)[++j]);
         for (int i=0; i<4; i++) {
            char s[100];
            str.getline(s,99,',');
            fNvar[i] = atoi(s);
            fWeightSum += fWeight[i] = fNvar[i]*fNvar[i];
         }
      }
      else if (strstr((*argv)[j],"--init") == (*argv)[j]) {
         std::istringstream str((*argv)[++j]);
         str >> fInitSol;
      }
   }
}

int SA_Problem::ReadProblemData(char *filename)
{   
   // Reads the problem data from file "filename".

   if (fDataImage) {
      delete fDataImage;
      fDataImage = 0;
   }
   if (fModelImage) {
      delete fModelImage;
      fModelImage = 0;
   }
   std::ifstream initfile(filename);
   if (initfile.good()) {
      initfile >> fDataImage;
      fModelImage = new TBitmap(fDataImage->GetWidth(),fDataImage->GetHeight());
      fModelImage->Precompute(fNdim);
      return 1;
   }
   return 0;
}
	
SA_Solution * SA_Problem::CreateSolution()
{ 
   // Create a new SA_Solution instance for this problem;
   // the user is responsible for deleting it when it is no longer needed.

   SA_Solution* newsol = new SA_Solution(fNdim);
   return newsol;
}

void SA_Problem::GetInitialSolution(SA_Solution& sol)
{ 
  // Set the SA_Solution sol to a predefined initialization value

   std::ifstream in(fInitSol.c_str());
   if (in.good()) {
      SA_Solution* s;
      in >> s;
      Copy(*s,sol);
      delete s;
      return;
   }

   if (sol.fNdim > 0) {
      UInt_t drefa = sol.fRefa->GetDegree();
      UInt_t dsig1a = sol.fSig1a->GetDegree();
      UInt_t dsig1ph = sol.fSig1ph->GetDegree();
      UInt_t dsig2ph = sol.fSig2ph->GetDegree();
      for (int i=0; i<drefa; i++)
         for (int j=0; j<drefa; j++)
            (*sol.fRefa)[i][j] = 0;
      for (int i=0; i<dsig1a; i++)
         for (int j=0; j<dsig1a; j++)
            (*sol.fSig1a)[i][j] = 0;
      for (int i=0; i<dsig1ph; i++)
         for (int j=0; j<dsig1ph; j++)
            (*sol.fSig1ph)[i][j] = 0;
      for (int i=0; i<dsig2ph; i++)
         for (int j=0; j<dsig2ph; j++)
            (*sol.fSig2ph)[i][j] = 0;
      (*sol.fRefa)[0][0] = 1;
      (*sol.fSig1a)[0][0] = 1;
      if (sol.fSig2a) {
         delete sol.fSig2a;
         sol.fSig2a = 0;
      }
      sol.fCost = ComputeCost(sol);
      sol.fNewcost = -1;
      if (sol.fLastmove) {
         delete sol.fLastmove;
         sol.fLastmove = 0;
      }
   }
}

void SA_Problem::GetRandomSolution(SA_Solution& sol)
{
  // Set the SA_Solution sol to a random value

   if (sol.fNdim > 0) {
      sol.fRefa->Randomize(fRefa_scale);
      sol.fSig1a->Randomize(fSig1a_scale);
      sol.fSig1ph->Randomize(fSig1ph_scale);
      sol.fSig2ph->Randomize(fSig2ph_scale);
      sol.fCost = ComputeCost(sol);
      sol.fNewcost = -1;
      if (sol.fSig2a) {
         delete sol.fSig2a;
         sol.fSig2a = 0;
      }
      if (sol.fLastmove) {
         delete sol.fLastmove;
         sol.fLastmove = 0;
      }
   }
}

double SA_Problem::ComputeCost(SA_Solution& sol)
{
   // Compute the cost function for this solution

   if (fDataImage == 0 || sol.fNdim != fNdim) {
      return -1;
   }
   if (sol.fSig2a == 0) {
      sol.fSig2a = new TModel(*(sol.fSig1a));
      *(sol.fSig2a) *= fSig2a_ratio;
   }
   unsigned int width = GetWidth();
   unsigned int height = GetHeight();
   if (sol.fRefaMap == 0) {
      sol.fRefaMap = new TBitmap(width,height);
      sol.fRefaMap->Precompute(sol.fNdim);
      sol.fRefaMap->Compute(*(sol.fRefa));
   }
   if (sol.fSig1aMap == 0) {
      sol.fSig1aMap = new TBitmap(width,height);
      sol.fSig1aMap->Precompute(sol.fNdim);
      sol.fSig1aMap->Compute(*(sol.fSig1a));
   }
   if (sol.fSig2aMap == 0) {
      sol.fSig2aMap = new TBitmap(width,height);
      sol.fSig2aMap->Precompute(sol.fNdim);
      sol.fSig2aMap->Compute(*(sol.fSig2a));
   }
   if (sol.fSig1phMap == 0) {
      sol.fSig1phMap = new TBitmap(width,height);
      sol.fSig1phMap->Precompute(sol.fNdim);
      sol.fSig1phMap->Compute(*(sol.fSig1ph));
   }
   if (sol.fSig2phMap == 0) {
      sol.fSig2phMap = new TBitmap(width,height);
      sol.fSig2phMap->Precompute(fNdim);
      sol.fSig2phMap->Compute(*(sol.fSig2ph));
   }
   fModelImage->Interfere(*(sol.fRefaMap),
                          *(sol.fSig1aMap),
                          *(sol.fSig2aMap),
                          *(sol.fSig1phMap),
                          *(sol.fSig2phMap));
   return fModelImage->ChiSquared(*fDataImage);
}

float SA_Problem::GetCost(SA_Solution& sol)
{
   // Return the pre-computed cost for this solution

   if (sol.fLastmove)
      return sol.fNewcost;
   else
      return sol.fCost;
}

void SA_Problem::GetNeighbor(SA_Solution& sol, float rel_temp)
{
   // Generate a reversible random change in the given 
   // solution. The range of change can be modified 
   // depending on relative temperature rel_temp
   // which ranges from 0 (=low) to 1 (=high).

   double wselect = (fWeightSum*rand()/RAND_MAX);
   double weight=0;
   int is;
   for (is = 0; is < 3; is++) {
      if (wselect < (weight += fWeight[is])) 
         break;
   }
   SA_Move* move = new SA_Move();
   move->fModel_set = is;

   TModel* orig;
   TModel* orig2;
   TBitmap* origmap;
   TBitmap* origmap2;
   double rescale = pow(rel_temp,0.08);
   switch (move->fModel_set) {
      case kRefa:
         orig = new TModel(*(sol.fRefa));
         sol.fRefa->Mutate(rescale*fRefa_step,fNvar[kRefa]);
         origmap = new TBitmap(*(sol.fRefaMap));
         sol.fRefaMap->Compute(*(sol.fRefa));
         sol.fNewcost = ComputeCost(sol);
         move->fBitmap = sol.fRefaMap;
         move->fModel = sol.fRefa;
         sol.fRefaMap = origmap;
         sol.fRefa = orig;
         break;
      case kSig1a:
         orig = new TModel(*(sol.fSig1a));
         sol.fSig1a->Mutate(rescale*fSig1a_step,fNvar[kSig1a]);
         origmap = new TBitmap(*(sol.fSig1aMap));
         sol.fSig1aMap->Compute(*(sol.fSig1a));
         orig2 = new TModel(*(sol.fSig2a));
         *sol.fSig2a = *sol.fSig1a;
         *sol.fSig2a *= fSig2a_ratio;
         origmap2 = new TBitmap(*(sol.fSig2aMap));
         *sol.fSig2aMap = *sol.fSig1aMap;
         *sol.fSig2aMap *= fSig2a_ratio;
         sol.fNewcost = ComputeCost(sol);
         move->fBitmap = sol.fSig1aMap;
         move->fModel = sol.fSig1a;
         sol.fSig1aMap = origmap;
         sol.fSig1a = orig;
         delete sol.fSig2aMap;
         sol.fSig2aMap = origmap2;
         delete sol.fSig2a;
         sol.fSig2a = orig2;
         break;
      case kSig1ph:
         orig = new TModel(*(sol.fSig1ph));
         sol.fSig1ph->Mutate(rescale*fSig1ph_step,fNvar[kSig1ph]);
         origmap = new TBitmap(*(sol.fSig1phMap));
         sol.fSig1phMap->Compute(*(sol.fSig1ph));
         sol.fNewcost = ComputeCost(sol);
         move->fBitmap = sol.fSig1phMap;
         move->fModel = sol.fSig1ph;
         sol.fSig1phMap = origmap;
         sol.fSig1ph = orig;
         break;
      case kSig2ph:
         orig = new TModel(*(sol.fSig2ph));
         sol.fSig2ph->Mutate(rescale*fSig2ph_step,fNvar[kSig2ph]);
         origmap = new TBitmap(*(sol.fSig2phMap));
         sol.fSig2phMap->Compute(*(sol.fSig2ph));
         sol.fNewcost = ComputeCost(sol);
         move->fBitmap = sol.fSig2phMap;
         move->fModel = sol.fSig2ph;
         sol.fSig2phMap = origmap;
         sol.fSig2ph = orig;
         break;
      default:
         std::cerr << "Fatal error in SA_Problem::GetNeighbor, "
                   << "this should never happen!" << std::endl;
         exit(1);
   }
   if (sol.fLastmove)
      delete sol.fLastmove;
   sol.fLastmove = move;
}

void SA_Problem::ResetSolution(SA_Solution& sol)
{
   // Discards the last change of solution made by GetNeighbor()

   if (sol.fLastmove) {
      delete sol.fLastmove;
      sol.fLastmove = 0;
   }
   sol.fNewcost = -1;
}

void SA_Problem::UpdateSolution(SA_Solution& sol)
{
   // Make the last move made by GetNeighbor() permanent

   if (sol.fLastmove == 0)
      return;

   switch (sol.fLastmove->fModel_set) {
      case kRefa:
         delete sol.fRefa;
         delete sol.fRefaMap;
         sol.fRefa = sol.fLastmove->fModel;
         sol.fRefaMap = sol.fLastmove->fBitmap;
         sol.fLastmove->fModel = 0;
         sol.fLastmove->fBitmap = 0;
         break;
      case kSig1a:
         delete sol.fSig1a;
         delete sol.fSig1aMap;
         sol.fSig1a = sol.fLastmove->fModel;
         sol.fSig1aMap = sol.fLastmove->fBitmap;
         sol.fLastmove->fModel = 0;
         sol.fLastmove->fBitmap = 0;
         delete sol.fSig2a;
         delete sol.fSig2aMap;
         sol.fSig2a = new TModel(*(sol.fSig1a));
         sol.fSig2aMap = new TBitmap(*(sol.fSig1aMap));
         *sol.fSig2a *= fSig2a_ratio;
         *sol.fSig2aMap *= fSig2a_ratio;
         break;
      case kSig1ph:
         delete sol.fSig1ph;
         delete sol.fSig1phMap;
         sol.fSig1ph = sol.fLastmove->fModel;
         sol.fSig1phMap = sol.fLastmove->fBitmap;
         sol.fLastmove->fModel = 0;
         sol.fLastmove->fBitmap = 0;
         break;
      case kSig2ph:
         delete sol.fSig2ph;
         delete sol.fSig2phMap;
         sol.fSig2ph = sol.fLastmove->fModel;
         sol.fSig2phMap = sol.fLastmove->fBitmap;
         sol.fLastmove->fModel = 0;
         sol.fLastmove->fBitmap = 0;
         break;
      default:
         std::cerr << "Fatal error in SA_Problem::GetNeighbor, "
                   << "this should never happen!" << std::endl;
   }
   sol.fCost = sol.fNewcost;
   sol.fNewcost = -1;
   delete sol.fLastmove;
   sol.fLastmove = 0;
}

float SA_Problem::GetLocalN()
{ 
   // Return  the number of solution neighbours (approximately).
   // The scaling with problem size is correct, but the magnitude
   // has been arbitrarily fixed to fit within the float range
   // for most problems.

   double logN = log(4)*(fNvar[0]*fNvar[0]
                       + fNvar[1]*fNvar[1]
                       + fNvar[2]*fNvar[2]
                       + fNvar[3]*fNvar[3]);
   return (logN < log(1e5))? exp(logN) : 1e5;
}

float SA_Problem::GetGlobalN()
{ 
   // Return the total number of distinguishable solutions (approximately).
   // The scaling with problem size is correct, but the magnitude
   // has been arbitrarily fixed to fit within the float range
   // for most problems.

   double logN = log(fRefa_scale/fRefa_step)*fNvar[0]*fNvar[0]
               + log(fSig1a_scale/fSig1a_step)*fNvar[1]*fNvar[1]
               + log(fSig1ph_scale/fSig1ph_step)*fNvar[2]*fNvar[2]
               + log(fSig2ph_scale/fSig2ph_step)*fNvar[3]*fNvar[3];
   return (logN < log(FLT_MAX))? exp(logN) : FLT_MAX;
}

void SA_Problem::Copy(SA_Solution& source, SA_Solution& dest)
{
   // Implements deep-copy semantics for SA_Solution

   dest.Copy(source);
}	

void SA_Problem::OutputSolution(std::ostream& out, SA_Solution& sol)
{ 
   // Serialize a solution into an output stream

  out << sol;
}

int SA_Problem::InputSolution(std::istream& in, SA_Solution& sol)
{
   // Reads a solution from an input stream

   SA_Solution* s;
   in >> s;
   Copy(*s,sol);
   delete s;
   return (in.gcount() > 0);
}

#ifdef MPI

int  SA_Problem::BcastSolution(SA_Solution& sol, int root,
                               MPI_Comm comm, int me_comm)
{ 
   // Broadcast a solution to all MPI slave processors

   std::cerr << "SA_Problem::BcastSolution not yet implemented!" << std::endl;
   return 0;
}

#endif

SA_Move* SA_Problem::CreateMove()
{
   // Create a new instance of SA_Move for this problem;
   // it is up to the user to delete it when it is no longer needed.

   return new SA_Move();
}

int SA_Problem::ExtractMove(SA_Solution& sol, SA_Move& mov)
{
   // Extract the current move from sol and stores it in mov

   if (sol.fLastmove) {
      mov.fModel_set = sol.fLastmove->fModel_set;
      if (mov.fModel)
         delete mov.fModel;
      if (mov.fBitmap)
         delete mov.fBitmap;
      mov.fModel = new TModel(*(sol.fLastmove->fModel));
      mov.fBitmap = new TBitmap(*(sol.fLastmove->fBitmap));
      return 1;
   }
   return 0;
}

int SA_Problem::ApplyMove(SA_Solution& sol, SA_Move& mov)	
{
   // applies the move mov to solution sol

   if (sol.fLastmove)
      delete sol.fLastmove;
   sol.fLastmove = new SA_Move(mov);
   UpdateSolution(sol);
   sol.fCost = ComputeCost(sol);
   return 1;
}

void SA_Problem::OutputMove(std::ostream& out, SA_Move& mov)
{
   // Write a move to an output stream

   out << mov;
}

int SA_Problem::InputApplyMoves(std::istream& in, SA_Solution& sol)
{
   // Read a number of moves from stream and change solution accordingly,
   // returning the number of steps successfully completed

   if (sol.fLastmove) {
      delete sol.fLastmove;
      sol.fLastmove = 0;
   }
   int nseen=0;
   while (!in.eof()) {
      in >> sol.fLastmove;
      UpdateSolution(sol);
      nseen++;
   }
   sol.fCost = ComputeCost(sol);
   return nseen;
}

#ifdef MPI

int SA_Problem::BcastApplyMove(SA_Solution& sol, int root,
                               MPI_Comm comm, int me_comm)
{
   // Broadcast a move and apply it to all slave processors

   std::cerr << "SA_Problem::BcastApplyMove not yet implemented!" << std::endl;
   return 0;
}

#endif

int SA_Problem::OutputSolution(char* filename, SA_Solution& sol)
{
   // Writes a solution to a file

   std::ofstream outfile(filename);
   if (!outfile) {
      std::cerr << "Error opening file: " << filename << std::endl;
      return 0;
   }
   else {
      outfile << sol;
      outfile.close();
   }
   return 1;
}

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

   std::ifstream infile(filename);
   if (!infile) {
      std::cerr << "Error opening file: " << filename << std::endl;
      return 0;
   }
   else {
      SA_Solution* s;
      infile >> s;
      Copy(*s,sol);
      delete s;
      infile.close();
   }
   return 1;
}

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

int SA_Problem::Reoptimize(SA_Solution& sol)
{
   return 1;
}

void SA_Problem::ControlOutput(std::ostream& out)
{
}