/*! \file
  Tabulated Function Based on the interpolation algorithm of DIVDIF (CERNLIB)
  \author Luis Prado Jr.
  \version $Id$
  \date 26 Jan 2003
*/

#ifndef _utl_TabulatedFunction_h_
#define _utl_TabulatedFunction_h_

#include <utl/TabulatedFunctionIterators.h>
#include <vector>
#include <algorithm>
#include <numeric>
#include <string>
#include <functional>


namespace utl {

  /*!
    \class TabulatedFunction TabulatedFunction.h "utl/TabulatedFunction.h"
    \author Luis Prado Jr.
    \brief Class to hold collection (x,y) points and provide interpolation between them

    The interpolation methods in this class are based on the interpolation
    algorithm of DIVDIF (CERNLIB)

    \version $Id$
    \date 26 Jan 2003
    \ingroup math
  */

  class TabulatedFunction {

  public:
    typedef std::vector<double> Array;
    typedef Array::iterator ArrayIterator;
    typedef Array::const_iterator ArrayConstIterator;
    typedef Array::reverse_iterator ArrayReverseIterator;
    typedef Array::const_reverse_iterator ArrayConstReverseIterator;
    typedef Array::const_reference ArrayConstReference;
    typedef TabulatedFunctionIterator Iterator;
    typedef ConstTabulatedFunctionIterator ConstIterator;

    TabulatedFunction(const unsigned int interpolationOrder = 1,
                      const double boundaryTolerance = 1e-3) :
      fInterpolationOrder(interpolationOrder),
      fBoundaryTolerance(boundaryTolerance)
    { }

    TabulatedFunction(const std::vector<double>& xValues,
                      const std::vector<double>& yValues,
                      const unsigned int interpolationOrder = 1,
                      const double boundaryTolerance = 1e-3) :
      fInterpolationOrder(interpolationOrder),
      fBoundaryTolerance(boundaryTolerance)
    {
      FillTable(xValues, yValues);
    }

    //virtual ~TabulatedFunction() { }

    unsigned int GetNPoints() const { return fX.size(); }

    /// Note: cannot be used in assignment
    Pair operator[](const int idx) const { return Pair(fX[idx], fY[idx]); }

    void Clear() { fX.clear(); fY.clear(); }

    void PushBack(const double x, const double y);

    Iterator Begin()
    { return Iterator(fX.begin(), fY.begin()); }

    Iterator End()
    { return Iterator(fX.end(), fY.end()); }

    ConstIterator Begin() const
    { return ConstIterator(fX.begin(), fY.begin()); }

    ConstIterator End() const
    { return ConstIterator(fX.end(), fY.end()); }

    /// begin of array of X
    ArrayIterator XBegin() { return fX.begin(); }

    /// begin of array of X
    ArrayConstIterator XBegin() const { return fX.begin(); }

    /// begin reverse iterator for X
    ArrayReverseIterator XRBegin() { return fX.rbegin(); }

    /// begin reverse iterator for X
    ArrayConstReverseIterator XRBegin() const { return fX.rbegin(); }

    /// begin of array of Y
    ArrayIterator YBegin() { return fY.begin(); }

    /// begin of array of Y
    ArrayConstIterator YBegin() const { return fY.begin(); }

    /// begin reverse iterator for Y
    ArrayReverseIterator YRBegin() { return fY.rbegin(); }

    /// begin reverse iterator for Y
    ArrayConstReverseIterator YRBegin() const { return fY.rbegin(); }

    /// end of array of X
    ArrayIterator XEnd() { return fX.end(); }

    /// end of array of X
    ArrayConstIterator XEnd() const { return fX.end(); }

    /// end reverse iterator for X
    ArrayReverseIterator XREnd() { return fX.rend(); }

    // end reverse iterator for X
    ArrayConstReverseIterator XREnd() const { return fX.rend(); }

    /// end of array of Y
    ArrayIterator YEnd() { return fY.end(); }

    /// end of array of Y
    ArrayConstIterator YEnd() const { return fY.end(); }

    /// end reverse iterator for Y
    ArrayReverseIterator YREnd() { return fY.rend(); }

    /// end reverse iterator for Y
    ArrayConstReverseIterator YREnd() const { return fY.rend(); }

    /// read-only reference to front of array of X
    ArrayConstReference XFront() const { return fX.front(); }

    /// read-only reference to front of array of Y
    ArrayConstReference YFront() const { return fY.front(); }

    /// read-only reference to back of array of X
    ArrayConstReference XBack() const { return fX.back(); }

    /// read-only reference to back of array of Y
    ArrayConstReference YBack() const { return fY.back(); }

    Iterator Insert(const double x, const double y);

    //ConstIterator Find(const double x, const double y) const;

    ConstIterator FindX(const double x) const;

    ConstIterator FindY(const double y) const;

    /// Get or interpolate the Y value that corresponds to parameter \p x
    double Y(const double x) const
    { return InterpolateY(x, fInterpolationOrder); }

    /// Interpolate the Y value with a \p polyDegree polynomial
    double InterpolateY(const double x, const unsigned int polyDegree) const;

    // Interpolate the Y value with the polynomial degree given in the constructor
    double operator()(const double x) const { return Y(x); }

    // get/set the \p idx-th value of the Ys
    const double& GetY(const unsigned int idx) const { return fY[idx]; }

    // get/set the \p idx-th value of the Xs
    const double& GetX(const unsigned int idx) const { return fX[idx]; }

    double& GetX(const unsigned int idx) { return fX[idx]; }
    double& GetY(const unsigned int idx) { return fY[idx]; }

    /// return the sum of X values
    double SumX() const { return std::accumulate(fX.begin(), fX.end(), 0.); }

    /// return the sum of Y values
    double SumY() const { return std::accumulate(fY.begin(), fY.end(), 0.); }

    bool IsInValidRange(const double x) const
    { return fX.front() <= x && x <= fX.back(); }

    bool operator==(const TabulatedFunction& t) const
    { return fX == t.fX && fY == t.fY; }

    bool operator!=(const TabulatedFunction& t) const
    { return !operator==(t); }

    unsigned int GetInterpolationOrder() const { return fInterpolationOrder; }

    void SetInterpolationOrder(const unsigned int interpolationOrder)
    { fInterpolationOrder = interpolationOrder; }

    double GetBoundaryTolerance() const { return fBoundaryTolerance; }

    void SetBoundaryTolerance(const double tolerance) { fBoundaryTolerance = tolerance; }

    void
    ScaleX(const double s)
    {
      std::transform(
        XBegin(), XEnd(), XBegin(),
        std::bind2nd(std::multiplies<double>(), s)
      );
    }

    void
    ScaleY(const double s)
    {
      std::transform(
        YBegin(), YEnd(), YBegin(),
        std::bind2nd(std::multiplies<double>(), s)
      );
    }

    void
    Swap(TabulatedFunction& tf)
    {
      std::swap(fX, tf.fX);
      std::swap(fY, tf.fY);
      std::swap(fInterpolationOrder, tf.fInterpolationOrder);
      std::swap(fBoundaryTolerance, tf.fBoundaryTolerance);
    }

  protected:
    void FillTable(const std::vector<double>& xValues,
                   const std::vector<double>& yValues);

    Array fX;
    Array fY;

  private:
    double LinearInterpolateY(const double x) const;

    unsigned int fInterpolationOrder;
    double fBoundaryTolerance;
  };

}

namespace std {

  template<>
  inline void swap(utl::TabulatedFunction& t1, utl::TabulatedFunction& t2)
  { t1.Swap(t2); }

}


#endif

// Configure (x)emacs for this file ...
// Local Variables:
// mode: c++
// compile-command: "make -C .. -k"
// End: