source: Sophya/trunk/SophyaLib/NTools/slininterp.cc@ 4051

/*    
   --- SOPHYA software - NTools module ---
   --- Class SLinInterp1D : Simple linear 1D interpolation class
   R. Ansari , C. Magneville   2000-2010 
   (C) UPS+LAL IN2P3/CNRS     (C) IRFU-SPP/CEA 
*/

#include "slininterp.h"
#include <fstream>

namespace SOPHYA {

/*! 
  \class SLinInterp1D
  \ingroup NTools
  \brief Simple linear interpolation class 

  \code
  #include "slininterp.h"
  #include "srandgen.h"

  ...
  vector<double> ys;
  double xmin = 0.5;
  double xmax = 0.;
  for(int i=0; i<=12; i++) {
    xmax = xmin+i*0.1;
    yreg.push_back(sin(xmax)*cos(2.2*xmax));
  }
  SLinInterp1D interpYR(xmin, xmax, yreg);
  cout << interpYR
  for(int i=0; i<=12; i++) {
    double x = drand01()*2.;
    cout << " Interpol result for X=" << x << " -> " << interpYR(x) << " ?= " << sin(x)*cos(2.2*x) << endl;
  }
  \endcode 
*/

/* --Methode-- */
SLinInterp1D::SLinInterp1D()
  : xmin_(0.), xmax_(1.), dx_(1.), ksmx_(1), npoints_(0)
{
  xs_.push_back(0.);
  xs_.push_back(1.); 
  ys_.push_back(0.);
  ys_.push_back(1.); 
}

/*!
  \brief Constructor from regularly spaced points in X with Y values defined by yreg 

  \b xmin and \b xmax are the two extreme points in X corresponding to yreg.
  Example: xmin=1, xmax=10, with yreg.size()=10 , 
  yreg[0]=Y(1) , yreg[1]=Y(2)  ...  yreg[9]=Y(10) 
*/
SLinInterp1D::SLinInterp1D(double xmin, double xmax, vector<double>& yreg)
  : xmin_(0.), xmax_(1.), dx_(1.), ksmx_(1), npoints_(0)
{
  DefinePoints(xmin, xmax, yreg);
}


/*!
  \brief Constructor from a set of \b (xs[i],ys[i]) pairs.

  if (npt > 0), interpolates to a finer regularly spaced grid, from \b xmin to \b xmax 
  with npt points.  use \b (xs[0],xs[xs.size()-1]) as limits if \b xmax<xmin 
*/
SLinInterp1D::SLinInterp1D(vector<double>& xs, vector<double>& ys, double xmin, double xmax, size_t npt)
  : xmin_(0.), xmax_(1.), dx_(1.), ksmx_(1), npoints_(0)
{
  DefinePoints(xs, ys, xmin, xmax, npt);
}

/* --Methode-- */
double SLinInterp1D::YInterp(double x)  const
{
  if (npoints_>0) {  // on utilise les points regulierement espace 
    long i = (long)((x-xmin_)/dx_);  
    if (i<0) return ( yreg_[0]+(x-xmin_)*(yreg_[1]-yreg_[0])/dx_ );
    if (i>=npoints_) return ( yreg_[npoints_]+(x-xmax_)*(yreg_[npoints_]-yreg_[npoints_-1])/dx_ );
    return (yreg_[i]+(x-X(i))/dx_*(yreg_[i+1]-yreg_[i]));
  }
  else {  // On utilise les points xs_,ys_ directement 
    if (x<=xs_[0]) return ( ys_[0]+(x-xs_[0])*(ys_[1]-ys_[0])/(xs_[1]-xs_[0]) );
    if (x>=xs_[ksmx_]) return ( ys_[ksmx_]+(x-xs_[ksmx_])*(ys_[ksmx_]-ys_[ksmx_-1])/(xs_[ksmx_]-xs_[ksmx_-1]) );

    size_t k=1;
    while(x>xs_[k]) k++; 
    if (k>=xs_.size()) {  // ne devrait pas arriver ...
      string emsg = " SLinInterp1D::YInterp()  out of range k -> BUG in code ";
      throw out_of_range(emsg);
    }

    double rv=ys_[k-1]+(x-xs_[k-1])*(ys_[k]-ys_[k-1])/(xs_[k]-xs_[k-1]);
    //    cout << " DBG- x=" << x << " k=" << k << " xs[k]=" << xs_[k] << " ys[k]" << ys_[k] 
    //   << " rv=" << rv << endl;
    return rv;
  }
}

/*!
  \brief Defines the interpolation points from regularly spaced points in X with Y values defined by yreg 

  \b xmin and \b xmax are the two extreme points in X corresponding to yreg.
  Example: xmin=1, xmax=10, with yreg.size()=10 , 
  yreg[0]=Y(1) , yreg[1]=Y(2)  ...  yreg[9]=Y(10) 
*/
void SLinInterp1D::DefinePoints(double xmin, double xmax, vector<double>& yreg)
{
  if (yreg.size()<2)  {
    string emsg = "SLinInterp1D::DefinePoints(xmin,xmax,yreg) Bad parameters yreg.size()<2 ";
    throw range_error(emsg);
  }
  xmin_ = xmin; 
  xmax_ = xmax;
  npoints_ = yreg.size()-1;
  dx_ = (xmax_-xmin_)/(double)npoints_;
  yreg_ = yreg;
}


/*!
  \brief Define the interpolation points from a set of \b (xs[i],ys[i]) pairs.

  if (npt > 0), interpolates to a finer regularly spaced grid, from \b xmin to \b xmax 
  with npt points.  use \b (xs[0],xs[xs.size()-1]) as limits if \b xmax<xmin 
*/
void SLinInterp1D::DefinePoints(vector<double>& xs, vector<double>& ys, double xmin, double xmax, size_t npt)
{
  if ((xs.size() != ys.size())||(xs.size()<2))  {
    string emsg = "SLinInterp1D::DefinePoints() Bad parameters (xs.size() != ys.size())||(xs.size()<2) ";
    throw range_error(emsg);
  }
  for(size_t k=1; k<xs.size(); k++) {
    if (xs[k-1]>=xs[k])  { 
       string emsg =  "SLinInterp1D::DefinePoints()  unsorted xs ";
      throw range_error(emsg);
    }
  }
  xs_=xs;
  ys_=ys;
  ksmx_=xs_.size()-1;
  npoints_ = npt;
  if (xmin>=xmax)  {
    xmin_ = xs_[0]; 
    xmax_ = xs_[ksmx_];
  }
  else {
    xmin_ = xmin; 
    xmax_ = xmax;
  }
  if (npoints_<1) {
    dx_=(xmax_-xmin_)/(double)(xs_.size()-1);
    return;
  }
  dx_ = (xmax_-xmin_)/(double)npoints_;
  yreg_.resize(npoints_+1);
  
  // Compute the the y values for regularly spaced x xmin <= x <= xmax 
  // and keep values in the yreg vector
  yreg_[0] = ys_[0];
  yreg_[npoints_] = ys_[ksmx_];
  size_t k=1;
  for(size_t i=0; i<npoints_; i++)  {
    double x = X(i);
    while(x>xs_[k]) k++;
    if (k>=xs_.size())  k=xs_.size()-1;   
    yreg_[i] = ys_[k-1]+(x-xs_[k-1])*(ys_[k]-ys_[k-1])/(xs_[k]-xs_[k-1]);
    //DBG cout << " DBG* i=" << i << " X(i)=" << X(i) << " yreg_[i]= " << yreg_[i] << " X^2= " << X(i)*X(i) 
    //DBG << " k=" << k << " xs[k]=" << xs_[k] << endl;
  }
  return;
}


/*!
  \brief Read  Y values from the file \b filename

  Read Y values ( one/line) for regularly spaced X's from file \b filename and call 
  DefinePoints(xmin, xmax, yreg). Return the number of Y values read.
  \param filename : input file name 
  \param xmin,xamx : X range limits 
  \param clm : comment character, lines starting with \b clm are ignored 
*/
size_t SLinInterp1D::ReadYFromFile(string const& filename, double xmin, double xmax, char clm)
{
  ifstream inputFile;
  inputFile.open(filename.c_str(), ifstream::in);  
#ifndef __DECCXX
// ifstream.is_open() ne passe pas avec OSF-cxx
  if(! inputFile.is_open()) {
    string emsg = "  SLinInterp1D::ReadYFromFile() problem opening file ";
    emsg += filename;
    throw runtime_error(emsg);
  }
#endif

  vector<double> xsv, ysv;
  size_t cnt=0;
  double cola;
  string eline;
  while(!inputFile.eof())  { 
    inputFile.clear(); 
    if (inputFile.peek() == (int)clm)  {  // skip comment lines 
      getline(inputFile, eline);
      continue;
    }
    inputFile >> cola; 
    if ( (!inputFile.good()) || inputFile.eof())   break;
    inputFile.ignore(1024,'\n');  // make sure we go to the next line
    //cout << cola<< "    "<<colb<<endl;
    ysv.push_back(cola);
    cnt++;
  }
  inputFile.close();
  cout << " SLinInterp1D::ReadYFromFile()/Info: " << cnt << " Y-values read from file " << filename << endl;  
  DefinePoints(xmin, xmax, ysv);
  return cnt;
}

/*!
  \brief Read pairs of  ( X Y ) values from file \b filename 

  Read pairs of (X Y) values, one pair / line from the specified file and call DefinePoints(xs, ys ...).
  One pair of space or tab separated numbers on each line. Return the number of Y values read.
  \param filename : input file name 
  \param xmin,xamx : X range limits. use the X limits from the file if xmax<xmin
  \param npt : number of points for regularly spaced interpolation points
  \param clm : comment character, lines starting with \b clm are ignored 
*/
size_t SLinInterp1D::ReadXYFromFile(string const& filename, double xmin, double xmax, size_t npt, char clm)
{
  ifstream inputFile;
  inputFile.open(filename.c_str(), ifstream::in); 
#ifndef __DECCXX 
// ifstream.is_open() ne passe pas avec OSF-cxx
  if(! inputFile.is_open()) {
    string emsg = "  SLinInterp1D::ReadXYFromFile() problem opening file ";
    emsg += filename;
    throw runtime_error(emsg);
  }
#endif

  vector<double> xsv, ysv;
  size_t cnt=0;
  double cola, colb;
  string eline;
  while(!inputFile.eof())  {  
    inputFile.clear();
    if (inputFile.peek() == (int)clm)  {  // skip comment lines 
      getline(inputFile, eline);
      continue;
    }
    inputFile >> cola >> colb; 
    if ( (!inputFile.good()) || inputFile.eof())   break;
    inputFile.ignore(1024,'\n');  // make sure we go to the next line
    //    cout << " DEBUG-GCount " << inputFile.gcount() << endl;
    //    cout << " DEBUG - cnt=" << cnt << " x=" << cola << "  y= "<<colb<<endl;
    xsv.push_back(cola);
    ysv.push_back(colb);
    cnt++;
  }
  inputFile.close();
  cout << " SLinInterp1D::ReadXYFromFile()/Info: " << cnt << " (x,y) pairs read from file " << filename << endl;  
  DefinePoints(xsv, ysv, xmin, xmax, npt);
  return cnt;
}


/* --Methode-- */
ostream& SLinInterp1D::Print(ostream& os, int lev)  const
{
  os << " ---- SLinInterp1D::Print() XMin=" << XMin() << " XMax=" << XMax() << " NPoints=" << npoints_ << endl;
  os << "  xs_.size()= " << xs_.size() << "  ys_.size()= " << ys_.size() << "  yreg_.size()= " << yreg_.size() << endl;
  if ((lev>0)&&(xs_.size()>0)) {
    for(size_t i=0; i<xs_.size(); i++) 
      os << " xs[" << i << " ]=" << xs_[i] << " -> ys[" << i << "]=" << ys_[i] << endl;
  } 
  if ((lev>0)&&(yreg_.size()>0)) {
    for(size_t i=0; i<yreg_.size(); i++) 
      os << " Regularly Spaced X(" << i << " )=" << X(i) << " -> yreg_[" << i << "]=" << yreg_[i] << endl;
  } 
  os << " -----------------------------------------------------------------------" << endl;
}

}  // End namespace SOPHYA