//-----------------------------------------------------------
// File and Version Information:
// $Id$
//
// Description:
//      Implementation of class BiCubSpline
//      see BiCubSpline.h for details
//
// Environment:
//      Software developed for the PANDA Detector at FAIR.
//
// Author List:
//      Felix Boehmer        TUM            (original author)
//
//
//-----------------------------------------------------------

// Panda Headers ----------------------

// This Class' Header ------------------
#include "BiCubSpline.h"

// C/C++ Headers ----------------------

// Collaborating Class Headers --------
#include "TError.h"
#include "TF2.h"
#include "TArrayD.h"
#include "SplineTF2Interface.h"
#include <sstream>
#include <string>




// Class Member definitions -----------



BiCubSpline::BiCubSpline() 
:  
  _func(NULL),
  _ifc(NULL),
  _M(NULL),
  _N(NULL)
{
  _xoutofrange=false;        //TODO: implement some kind of "initialized" flag
  _youtofrange=false;	
  _xLength = 0;
  _yLength = 0;
    
}  

BiCubSpline::BiCubSpline(const std::vector<double>& kx,
			 const std::vector<double>& ky,
			 const std::vector<std::vector<double>*>* c)
:
  _func(NULL),
  _ifc(NULL)   
{
  _xoutofrange=false;
  _youtofrange=false;
  
  _kx = kx;
  _ky = ky;
  _xLength = _kx.size();
  _yLength = _ky.size();
  
  if(_xLength<5 || _yLength<5)
    Fatal("BiCubSpline::BiCubSpline()", "not enough knots, aborting...");

  _coeffs.ResizeTo(_yLength-4,_xLength-4);
  if(c!=NULL) {
    for(unsigned int iy=0; iy<_yLength-4; iy++)
      for(unsigned int ix=0; ix<_xLength-4; ix++)
	_coeffs[iy][ix] = (*(*c)[iy])[ix];
  }
  else {    //build standard coeffs...
    for(unsigned int iy=0; iy<_yLength-4; iy++)
      for(unsigned int ix=0; ix<_xLength-4; ix++)
	_coeffs[iy][ix] = 1.;
  }          
  //first 4 knots in x and y dir. have no attached Splines...
  if(_coeffs.GetNrows()+4==_yLength && _coeffs.GetNcols()+4==_xLength)
    ;
  else
    Fatal("BiCubSpline::BiCubSpline()", "coeffs-knots inconsistent! Aborting");
  //build B-Spline basis
  _M = new TClonesArray("BSpline");
  _N = new TClonesArray("BSpline");
  for(int i=0; i<_xLength; ++i) {
    if(i>3) 
      new ((*_M)[_M->GetEntriesFast()]) BSpline(_kx,i);
  }
  for(int j=0; j<_yLength; ++j) {
    if(j>3) 
      new ((*_N)[_N->GetEntriesFast()]) BSpline(_ky,j);
  }
}

BiCubSpline::BiCubSpline(const TMatrixT<double>* knots,
			 const TMatrixT<double>* c)
:
  _func(NULL),
  _ifc(NULL)   
{
  _xoutofrange=false;
  _youtofrange=false;
  
  std::vector<double> kx;
  for(unsigned int ik=0; ik<knots->GetNcols(); ik++)
    kx.push_back((*knots)[0][ik]);
  std::vector<double> ky;
  for(unsigned int ik=0; ik<knots->GetNrows(); ik++)
    ky.push_back((*knots)[ik][0]);
      
  _kx = kx;
  _ky = ky;
  _xLength = _kx.size();
  _yLength = _ky.size();
  
  if(_xLength<5 || _yLength<5)
    Fatal("BiCubSpline::BiCubSpline()", "not enough knots, aborting...");

  _coeffs.ResizeTo(_yLength-4,_xLength-4);
  if(c!=NULL) {
    for(unsigned int iy=0; iy<_yLength-4; iy++)
      for(unsigned int ix=0; ix<_xLength-4; ix++)
	_coeffs[iy][ix] = (*c)[iy][ix];
  }
  else {    //build standard coeffs...
    for(unsigned int iy=0; iy<_yLength-4; iy++)
      for(unsigned int ix=0; ix<_xLength-4; ix++)
	_coeffs[iy][ix] = 1.;
  }          
  //first 4 knots in x and y dir. have no attached Splines...
  if(_coeffs.GetNrows()+4==_yLength && _coeffs.GetNcols()+4==_xLength)
    ;
  else
    Fatal("BiCubSpline::BiCubSpline()", "coeffs-knots inconsistent! Aborting");
  //build B-Spline basis
  _M = new TClonesArray("BSpline");
  _N = new TClonesArray("BSpline");
  for(int i=0; i<_xLength; ++i) {
    if(i>3) 
      new ((*_M)[_M->GetEntriesFast()]) BSpline(_kx,i);
  }
  for(int j=0; j<_yLength; ++j) {
    if(j>3) 
      new ((*_N)[_N->GetEntriesFast()]) BSpline(_ky,j);
  }
}


BiCubSpline::~BiCubSpline() {
  if(_func!=NULL)
    delete _func;
  if(_ifc!=NULL)
    delete _ifc;
  if(_M!=NULL)
    delete _M;
  if(_N!=NULL)
    delete _N;
}
  
void
BiCubSpline::setCoeffsByArray(double* ca) {   //TO DO: checks. (impossible!)
  //have to follow TF2 interface restrictions and take a double* !
  //user is responsible for proper array format!
  for(int y=0; y<_yLength-4; y++) {
    for(int x=0; x<_xLength-4; x++) 
      _coeffs[y][x] = ca[y*(_xLength-4) + x];
  }
}  


double 
BiCubSpline::eval(double x, double y) {
  //search for the right BSplines:
  int i,j;
  for(i=0; i<_xLength; ++i) {
    if(_kx[i]>x)
      break;
  }
  for(j=0; j<_yLength; ++j) {
    if(_ky[j]>y)
      break;
  }
  
  //catch out of bound cases
  if(i>=_xLength-4 || i<4) {  //x was out of range
    if(_xoutofrange==false) {
      Warning("BiCubSpline::eval()", "x was out of range!");
      _xoutofrange=true;
    }
    return 0;
  }
  if(j>=_yLength-4 || j<4) {  //y was out of range
	if(_youtofrange==false) {
	  Warning("BiCubSpline::eval()", "y was out of range!");
	  _youtofrange=true;
	}
    return 0;
  }

  double result=0;

  for(int k=i; k<i+4; k++) {
    for(int l=j; l<j+4; l++) 
      result+=((BSpline*)_M->At(k-4))->eval(x) * ((BSpline*)_N->At(l-4))->eval(y) * _coeffs[l-4][k-4];
  }

  return result;
}

TF2*
BiCubSpline::getTF2(double xMin, double xMax, double yMin, double yMax) {
  
  
  if(_func != NULL) 
    return _func;

  bool fail = false;
  
  //check for consistency:
  if(xMin < _kx[0] || xMax > _kx[_kx.size()-1])
    fail = true;
  if(yMin < _ky[0] || yMax > _ky[_ky.size()-1])
    fail = true;
  if(fail) {
    std::cerr<<"\nBiCubSpline::getTF2: invalid range! returning NULL"<<std::endl;
    return NULL;
  }
  
  //need a unique name for the TF2, will use stream cast of the this pointer
  std::stringstream ss;
  ss<<this;
  std::string name = ss.str();
  
  _ifc =  new SplineTF2Interface(this);
  _ifc->setReadOnly(true);
  
  _func = new TF2(name.c_str(),_ifc,&SplineTF2Interface::eval,xMin,xMax,
		  yMin,yMax,(_xLength+4)*(_yLength+4),"Function", "eval");
  
  //dummy parameters, we are in read-only mode
  TArrayD blub = TArrayD((_xLength+4)*(_yLength+4));
  double* dummy = blub.GetArray();
  _func->SetParameters(dummy);
  _func->SetNpx(50);
  _func->SetNpy(50);

  return _func;
}
  

  
  
  


ClassImp(BiCubSpline)