/// Evaluates the FOM over the parameter space
/// (slices of FOM surface can be plotted by FitPlot tools)
///
/// J Amey <jjtamey@googlemail.com>
///
/// 08/12 - New File
////////////////////////////////////////////////////////

#include <RAT/OptimisedComponent.hh>
#include <RAT/Surface.hh>
#include <RAT/DS/FitResult.hh>
#include <RAT/Log.hh>
#include <TFile.h>
#include <TTree.h>
#include <iostream>
#include <fstream>
#include <sstream>
#include <math.h>
using namespace RAT::Optimisers;

#include <CLHEP/Units/SystemOfUnits.h>
#include <iostream>
using namespace CLHEP;
using namespace std;

void Surface::Initialise( const std::string& param )
{
  // Convert input string to fDiv, else make fDiv take default val
  if ( param.empty() ){
    fDiv=40;
  } else {
    stringstream strStream(param);
    strStream >> fDiv;
    if ( fDiv % 2 ) warn << "Surface::Initialise: No. of divisions must be even"
                         << newline;
  }

}

double
Surface::Minimise()
{
  Evaluate();
  return 1.0;
}

double
Surface::Maximise()
{
  Evaluate();
  return 1.0;
}

void Surface::Evaluate()
{
  // Get initial seed parameters
  vector<double> params = fComponent->GetParams();
  // Use the positive position errors since only 1 set of them is used below
  vector<double> errors = fComponent->GetPositiveErrors();
  std::string methodName = "TEMP, fix in #436"; //fComponent->GetName();
  fTotalDim=params.size();

  // Initialise starting point of mixed space grid and increments
  fStartp = params;
  fIncs = params;

  // Range to be evaluated is from (seed-error) to (seed+error)
  for (int ii=0; ii<fTotalDim; ii++) {
    fStartp[ii] -=errors[ii];
    fIncs[ii] =  ( 2.0 * errors[ii] ) / (double) fDiv  ;
  }

  stringstream fileTitle;
  fileTitle << methodName << "fomsurface.root";
  TFile* f = TFile::Open(fileTitle.str().c_str(),"RECREATE");

  ft = new TTree("Tsurf","Tree with surface");
  ft->Branch("data", &fData);

  TTree* i = new TTree("Tincs","Increments");
  i->Branch("incs", &fIncs);

  TTree* s = new TTree("Tseeds","Seeds");
  s->Branch("seeds", &params);

  i->Fill();
  s->Fill();

  // Call recursive Scan member function
  Scan( params, 0 );

  f->Write();
  f->Close();
}

void Surface::Scan( std::vector<double> coords, int dimension )
{
  for( int spacePosition=0 ; spacePosition<fDiv; spacePosition++  )
  {
    coords[dimension] = fStartp[dimension] + (spacePosition * fIncs[dimension]);
    if( dimension != fTotalDim - 1 ) {
      Scan( coords, dimension + 1 );
    } else {
          fData.clear();
          for (int q=0; q<fTotalDim; q++) fData.push_back(coords[q]);
          double fom = (*fComponent)( coords );
          fData.push_back(fom);
          //cout << "DataRow: " << fData[0] << " " << fData[1] << " " << fData[2] << " " << fData[3] << " " << fData[4] << endl;
          ft->Fill();
      }
  }

}