#ifndef __FILE_LINEAR_INTERPOLATOR_H_SEEN__
#define __FILE_LINEAR_INTERPOLATOR_H_SEEN__

/* ----------------------------------------------------------------------------

Copyright (C) 2018

A. Ronald Gallant
Post Office Box 659
Chapel Hill NC 27514-0659
USA

This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.

-------------------------------------------------------------------------------

This header defines a linear interpolator that interpolates between
points stored in realmat x (abcissae) with corresponding values in
realmat y (ordinates).  When used in interpolator calls, x and y will
get sorted.  The interpolator is implemented as a vector of linear
functions. Typical usage is
  
  linear_interpolator f(x,y);
  REAL xval = 10.0;
  REAL yval = f(xval);
  REAL dydx = f.derivative(xval);

Equivalently

  linear_interpolator f;
  f.update(x,y);
  REAL xval = 10.0;
  REAL yval = f(xval);
  REAL dydx = f.derivative(xval); 

Functionality is as follows:

class linear_function 
  linear_function ()
  void initialize(REAL intercept, REAL slope, REAL origin)
  REAL operator()(REAL x) const
  REAL intercept() const
  REAL slope() const
  REAL origin() const

class linear_interpolator
  linear_interpolator()
  linear_interpolator(scl::realmat& x, scl::realmat& y)
  void update(scl::realmat& x, scl::realmat& y)
  REAL operator()(REAL x) const 
  linear_function operator[](size_t i) 
  size_t size() 

-----------------------------------------------------------------------------*/

#include "realmat.h"

namespace scl {

  class linear_function {
  private:
    REAL a;
    REAL b;
    REAL x0;
  public:
    void initialize(REAL intercept, REAL slope, REAL origin)
      { a = intercept; b = slope; x0 = origin; }
    REAL operator()(REAL x) const { return a + b*(x - x0); }
    REAL intercept() const { return a; }
    REAL slope() const { return b; }
    REAL origin() const { return x0; }
  };
  
  class linear_interpolator { // Note: input x & y are sorted by update 
  private:                    // and hence by the two argument constructor.
    std::vector<linear_function> funcs;
    typedef std::vector<linear_function>::size_type lfst;
    REAL xmin;
    REAL xmax;
    lfst N;
    lfst hash(REAL x) const { return lfst( REAL(N-2)*(x-xmin)/(xmax-xmin) ); }
  public:
    linear_interpolator() 
    {
      scl::realmat grid(2,1) ; grid[1] = 0.0; grid[2] = 1.0;
      scl::realmat vals(2,1) ; vals[1] = 0.0; vals[2] = 1.0;
      update(grid,vals);
    }
    linear_interpolator(scl::realmat& x, scl::realmat& y) { update(x,y); };
    void update(scl::realmat& x, scl::realmat& y) 
    { 
      INTEGER n = x.size(); N = lfst(n);  
      funcs.clear(); funcs.reserve(N);
      if (n<2) 
        scl::error("Error, linear_interpolator, x.size() < 2");
      if (x.ncol() != 1 || y.ncol() != 1) 
        scl::error("Error, linear_interpolator, x or y not a vector");
      if (n != y.size()) 
        scl::error("Error, linear_interpolator, x and y sizes differ");
      scl::intvec rank = x.sort(); 
      y = y(rank,"");
      xmin = x[1]; xmax = x[n];
      linear_function f;
      for (INTEGER i=1; i<n; ++i) {
        f.initialize(y[i], (y[i+1]-y[i])/(x[i+1]-x[i]), x[i]);
        funcs.push_back(f);
      }
      funcs.push_back(f);
    }
    REAL operator()(REAL x) const 
    {
       if (x <= funcs[0].origin()) return funcs[0](x);
       if (x >= funcs[N-1].origin()) return funcs[N-1](x);
       lfst i = hash(x);
       if (x < funcs[i].origin()) while(x < funcs[--i].origin());
       else if (x >= funcs[i+1].origin()) while(x >= funcs[++i+1].origin());
       return funcs[i](x);
    }
    linear_function operator[](lfst i) {return funcs[i];}
    lfst size() { return funcs.size(); }
  
  };

}

#endif
