#ifndef __FILE_SPLINE_INTERPOLATOR_H_SEEN__
#define __FILE_SPLINE_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 spline 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 spline 
functions.  Typical usage is
  
  spline_interpolator f(x,y);
  REAL xval = 10.0;
  REAL yval = f(xval);
  REAL dydx = f.derivative(xval);

Equivalently

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

Functionality is as follows:

class cubic_function
  spline_function()
  void initialize
    (REAL intercept, REAL linear, REAL quadratic, REAL cubic, REAL origin)
  REAL operator()(REAL x) const
  REAL derivative(REAL x) const
  REAL intercept() const
  REAL linear() const
  REAL quadratic() const
  REAL cubic()
  REAL origin()

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

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

#include "realmat.h"

namespace scl {

  class cubic_function {
  private:
    REAL a;
    REAL b;
    REAL c;
    REAL d;
    REAL x0;
  public:
    void initialize
      (REAL intercept, REAL linear, REAL quadratic, REAL cubic, REAL origin)
      { a = intercept; b = linear; c = quadratic; d = cubic; x0 = origin; }
    REAL operator()(REAL x) const 
      { REAL z=x-x0; REAL y=d; y=y*z+c; y=y*z+b; y=y*z+a; return y; }
    REAL derivative(REAL x) const 
      { REAL z=x-x0; REAL dy=3.0*d; dy=dy*z+2.0*c; dy=dy*z+b; return dy; }
    REAL intercept() const { return a; }
    REAL linear() const { return b; }
    REAL quadratic() const { return c; }
    REAL cubic() const { return d; }
    REAL origin() const { return x0; }
  };
  
  class spline_interpolator { // Note: input x & y are sorted by update 
  private:                    // and hence by the two argument constructor.
    std::vector<cubic_function> funcs;
    typedef std::vector<cubic_function>::size_type cfst;
    REAL xmin;
    REAL xmax;
    cfst N;
    cfst hash(REAL x) const { return cfst( REAL(N-2)*(x-xmin)/(xmax-xmin) ); }
  public:
    spline_interpolator() 
    {
      scl::realmat grid(3,1) ; grid[1] = 0.0; grid[2] = 0.5; grid[3] = 1.0;
      scl::realmat vals(3,1) ; vals[1] = 0.0; vals[2] = 0.5; vals[3] = 1.0;
      update(grid,vals);
    }
    spline_interpolator(scl::realmat& x, scl::realmat& y) { update(x,y); };
    void update(scl::realmat& x, scl::realmat& y) 
    { 
      INTEGER m = x.size(); N = cfst(m-1);  
      funcs.clear(); funcs.reserve(N);
      if (m<3) 
        scl::error("Error, spline_interpolator, x.size() < 3");
      if (x.ncol() != 1 || y.ncol() != 1) 
        scl::error("Error, spline_interpolator, x or y not a vector");
      if (m != y.size()) 
        scl::error("Error, spline_interpolator, x and y sizes differ");
      scl::intvec rank = x.sort(); 
      y = y(rank,"");
      xmin = x[1]; xmax = x[m];
      // Begin from https://en.wikipedia.org/wiki/Spline_(mathematics)
      int n = m - 1;
      int i, j;
      double x0[n+1],a[n+1];
      for (i = 0; i < n+1; ++i) x0[i] = x[i+1];
      for (i = 0; i < n+1; ++i) a[i] = y[i+1];
      double b[n],d[n];
      double h[n];
      for (i = 0; i <= n-1; ++i) h[i] = x0[i+1] - x0[i];
      double A[n];
      for (i = 1; i <= n-1; ++i) {
        A[i] = 3.0*(a[i+1] - a[i])/h[i] - 3.0*(a[i] - a[i-1])/h[i-1];
      }
      double c[n+1],l[n+1],u[n+1],z[n+1];
      l[0] = 1.0; u[0] = z[0] = 0.0;
      for (i = 1; i <= n-1; ++i) {
        l[i] = 2.0*(x0[i+1] - x0[i-1]) - h[i-1]*u[i-1];
        u[i] = h[i]/l[i];
        z[i] = (A[i] - h[i-1]*z[i-1])/l[i];
      }
      l[n] = 1; z[n] = c[n] = 0;
      for (j = n-1; j >= 0; --j) {
        c[j] = z[j] - u[j]*c[j+1];
        b[j] = (a[j+1] - a[j])/h[j] - h[j]*(c[j+1] + 2.0*c[j])/3.0;
        d[j] = (c[j+1] - c[j])/(3.0*h[j]);
      }
      // End from https://en.wikipedia.org/wiki/Spline_(mathematics)
      cubic_function f;
      for (i = 0; i < n; ++i) {
        f.initialize(a[i], b[i], c[i], d[i], x0[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);
       cfst 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);
    }
    REAL derivative(REAL x) const 
    {
       if (x <= funcs[0].origin()) return funcs[0].derivative(x);
       if (x >= funcs[N-1].origin()) return funcs[N-1].derivative(x);
       cfst 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].derivative(x);
    }
    cubic_function operator[](cfst i) { return funcs[i]; }
    cfst size() { return funcs.size(); }
   };

}

#endif
