//******************************************************
//  Copyright 2007 University of Connecticut
//
//  FILE : SA_MIProblem.h
//  
//  includes         : mpi.h, iostream, TModel.h, TBitmap.h
//  is included by   : SA_main.cc
//
//  description      : 
//    This SA_Problem represents the two-wave and three-wave inverse
//    Michelson interferometry problem.  The problem begins with a 2d
//    intensity image seen in the sensor plane of the interferometer.
//    The goal is to find a model for the surface profiles of the front
//    and back surfaces of the film, as encoded in the observed interference
//    fringe pattern in the intensity image.  A by-product of the solution
//    is a model for the amplitude profiles in the signal and reference arms
//    of the interferometer.  The goal of the optimization is to minimize
//    the chi_square between the intensity image computed on the basis of
//    a parameterized model and the observed image.  That is,
//                    ___   ___
//                    \     \                       2
//    chi_square =    /     /    ( model  -  data  )
//                    ---   ---         ij       ij
//                     i     j
//    where i and j are the row and column indices of the interference
//    image, ranging from i=0 to WIDTH-1 and j=0 to HEIGHT-1 where WIDTH
//    and HEIGHT are the dimensions of the rectangular intensity image in
//    pixels. The model for a Michelson interferometer is given by
//
//                  2              2              2
//    model   = refa (x ,y ) + sig1 (x ,y ) + sig2 (x ,y ) +
//         ij          i  j           i  j           i  j
//
//            + 2 refa(x ,y ) sig1(x ,y ) cos( sig1ph(x ,y ) ) +
//                      i  j        i  j               i  j 
//
//            + 2 refa(x ,y ) sig2(x ,y ) cos( sig2ph(x ,y ) ) +
//                      i  j        i  j               i  j
//
//            + 2 sig1a(x ,y ) sig2(x ,y ) cos( sig1ph(x ,y ) - sig2ph(x ,y ) )
//                       i  j        i  j               i  j            i  j
//
//    where refa(x,y) represents the amplitude profile in the reference beam,
//    sig1(x,y) represents that of the signal beam, and sig2(x,y) is a copy
//    of sig1(x,y) multiplied by an attenuation factor to account for loss
//    beam intensity in reflections from the back surface of the film relative
//    to the front surface arising from attenuation and transmission factors.
//    Functions sig1ph(x,y) and sig2ph(x,y) are the surface profiles of the
//    front and back surfaces of the film, respectively, in units of 1/k
//    where k is the wavenumber of the light beams in the interferometer.
//    The model assumes the reference mirror to be optically flat, so that
//    it may serve to define the phase reference in the problem.  Coordinates
//       x(i) = i*(2.0/(WIDTH-1))-1
//       y(j) = j*(2.0/(HEIGHT-1))-1
//    serve to map the image pixel coordinates onto the real number range
//    [-1,1] in x and y, so that the model functions are defined on the
//    square area with corners at (-1,-1) and (1,1).
//
//    The parameterization of the model functions is encapsulated within
//    the TModel class, which exposes a two-dimensional subscript operator
//    model[m][n] referencing the parameters.  The indices m and n both
//    range from 0 to ndim-1, where ndim is specified as the dimension of
//    the model space passed to the TModel constructor, so that the total
//    number of parameters per model is ndim*ndim.  A full 3-wave
//    interference problem involves the four independent models refa, sig1a,
//    sig1ph, sig2ph so the dimension of the solution space is 4*ndim*ndim.
//
//  authors:
//    Richard Jones
//  contact:
//    jonesrt@phys.uconn.edu
//  last update:
//    July 26, 2007
//
//******************************************************


#ifndef SA_Problem_H
#define SA_Problem_H

#include <iostream>

#ifdef MPI
  #include <mpi.h>
#endif

#include <TModel.h>
#include <TBitmap.h>

class SA_Problem;
class SA_Move;


//******************************************************
// class SA_Solution :
//	- contains a solution 
//	- can copy a solution into itself
//
// class MUST be implemented by the user 
//****************************************************** 
class SA_Solution
{
  friend class SA_Problem;

  private:
   int fNdim;
   double fCost, fNewcost;
   TModel* fRefa;
   TModel* fSig1a;
   TModel* fSig2a;
   TModel* fSig1ph;
   TModel* fSig2ph;
   TBitmap* fRefaMap;
   TBitmap* fSig1aMap;
   TBitmap* fSig2aMap;
   TBitmap* fSig1phMap;
   TBitmap* fSig2phMap;
   SA_Move* fLastmove;

  private:  
   SA_Solution(); 
   SA_Solution(const unsigned int ndim); 

  public:
   ~SA_Solution(); 
   void Copy(SA_Solution& source);
   int GetNdim();

  friend std::ostream& operator<<(std::ostream& out, SA_Solution& sol);
  friend std::istream& operator>>(std::istream& in, SA_Solution*& sol);
};

//******************************************************
// class SA_Move :
//      - contains a move 
//
// class CAN be implemented by the user
//****************************************************** 
class SA_Move 
{ 
  friend class SA_Problem;

  private:
   int fModel_set;
   TModel* fModel;
   TBitmap* fBitmap;

  public: 
   SA_Move();
   ~SA_Move();
   SA_Move(SA_Problem& prob);
   SA_Move(const SA_Move& move);

  friend std::ostream& operator<<(std::ostream& out, SA_Move& move);
  friend std::istream& operator>>(std::istream& in, SA_Move*& move);
};

//******************************************************
// class SA_Problem :
//      - representation of the problem
//      - generates initial solution
//      - neighborhood operation 
//      - cost function
//
//  class MUST be implemented by user
//****************************************************** 
class SA_Problem
{
  public: 
   static const int kRefa=0;
   static const int kSig1a=1;
   static const int kSig1ph=2;
   static const int kSig2ph=3;

  private:
   TBitmap* fDataImage;    // data interference bitmap being fitted
   TBitmap* fModelImage;   // model interference bitmap of latest solution
   int fNdim;              // declared order of basis used by models
   int fNvar[4];           // variational order of basis used by models
   double fWeight[4];      // relative weights for variation of each model
   double fWeightSum;      // sum of elements of fWeight, for normalization
   std::string fInitSol;   // name of file containing initial solution

// The following factor fixes the amplitude of the reflection from the
// back surface of the sample film relative to that from the front surface.
// It is computed based on value of 2.417 for the diamond refractive index.
   static const double fSig2a_ratio = -0.83;

// The following factors set the scale of the solution neighborhood by
// which solutions differ from their neighbors in each model parameter
// set.  All parameters in a set vary on the same scale in a single step.
   static const double fRefa_step = 0.2;
   static const double fSig1a_step = 0.1;
   static const double fSig1ph_step = 0.3;
   static const double fSig2ph_step = 0.3;

// The following factors set the scale of the entire solution space which
// should be searched for the optimum.  There is no rigorous mathematical
// bound on the coefficients in the model space, but there is a practical
// limit on the number of distinguishable solutions that can be represented
// by a finite resolution bitmap interference image.
   static const double fRefa_scale = 10;
   static const double fSig1a_scale = 4;
   static const double fSig1ph_scale = 100;
   static const double fSig2ph_scale = 100;

  public:
   SA_Problem(int ndim=5);
   ~SA_Problem();
   void EvaluateOptions(int* argc, char*** argv);
   int ReadProblemData(char* file);
   SA_Solution* CreateSolution();
   void GetInitialSolution(SA_Solution& sol);	
   void	GetRandomSolution(SA_Solution& sol);
   float GetCost(SA_Solution& sol); 
   void GetNeighbor(SA_Solution& sol, float rel_temp);
   void ResetSolution(SA_Solution& sol);
   void UpdateSolution(SA_Solution& sol);
   int GetNdim();
   int GetWidth();
   int GetHeight();
   float GetLocalN();
   float GetGlobalN();	
   int GetNvar(int model);
   void SetNvar(int model, int order);
   double GetWeight(int model);
   void SetWeight(int model, double weight);
   void Copy(SA_Solution& source, SA_Solution& dest);		
   void OutputSolution(std::ostream& out, SA_Solution& sol);
   int InputSolution(std::istream& in, SA_Solution& sol);
			
#ifdef MPI								
   int BcastSolution(SA_Solution& sol, int root, MPI_Comm comm, int me_comm);
#endif

   SA_Move* CreateMove();
   int ExtractMove(SA_Solution& sol, SA_Move& mov);   
   int ApplyMove(SA_Solution& sol, SA_Move& mov);
   void OutputMove(std::ostream& out, SA_Move& mov);
   int InputApplyMoves(std::istream& in, SA_Solution& sol);	

#ifdef MPI
   int BcastApplyMove(SA_Solution& sol, int root, MPI_Comm comm, int me_comm);
#endif

   int OutputSolution(char* filename, SA_Solution& sol);
   int InputSolution(char* filename, SA_Solution& sol);
   void Postprocessing(SA_Solution& sol);
   int Reoptimize(SA_Solution& sol);
   void ControlOutput(std::ostream& out = std::cout);

  private:
   double ComputeCost(SA_Solution& sol); 
};

inline int SA_Solution::GetNdim()
{
   return fNdim;
}

inline int SA_Problem::GetNdim()
{
   return fNdim;
}

inline int SA_Problem::GetWidth()
{
   if (fDataImage)
      return fDataImage->GetWidth();
   else
      return 0;
}

inline int SA_Problem::GetHeight()
{
   if (fDataImage)
      return fDataImage->GetHeight();
   else
      return 0;
}

inline int SA_Problem::GetNvar(int model)
{
   if (model >= 0 && model < 5)
      return fNvar[model];
   else
      return -1;
}

inline void SA_Problem::SetNvar(int model, int order)
{
   if (model >= 0 && model < 5)
      fNvar[model] = order;
}

inline double SA_Problem::GetWeight(int model)
{
   if (model >= 0 && model < 5)
      return fWeight[model];
   else
      return 0;
}

inline void SA_Problem::SetWeight(int model, double weight)
{
   if (model >= 0 && model < 5)
      fWeight[model] = weight;
   fWeightSum = 0;
   for (int i=0; i<5; i++) {
      fWeightSum += fWeight[i] = (fNvar[i] == 0)? 0 : fWeight[i];
   }
}

#endif