/** \file
    \brief R's spline function wrapped into template class that can be used with TMB.
*/

// Atomic functions required by splines
namespace atomic {

TMB_ATOMIC_VECTOR_FUNCTION(
                           // ATOMIC_NAME
                           findInt
                           ,
                           // OUTPUT_DIM
                           tx.size() - CppAD::Integer(tx[0]) - 1
                           ,
                           // ATOMIC_DOUBLE
                           int nx = (int) tx[0];
                           int nu = tx.size() - nx - 1;
                           int n_1 = nx - 1;
                           const double* x = &tx[1];
                           const double* v = x + nx;
                           int i = 0;
                           for(int l = 0; l < nu; l++) {
                             double ul = v[l];
                             if(ul < x[i] || (i < n_1 && x[i+1] < ul)) {
                               /* reset i  such that  x[i] <= ul <= x[i+1] : */
                               i = 0;
                               int j = nx;
                               do {
                                 int k = (i+j)/2;
                                 if(ul < x[k]) j = k;
                                 else i = k;
                               }
                               while(j > i+1);
                             }
                             ty[l] = i;
                           }
                           ,
                           // ATOMIC_REVERSE
                           for (size_t i=0; i<px.size(); i++) px[i] = 0;
                           )

template<class Type>
vector<Type> findInterval(vector<Type> xt, vector<Type> x) {
  CppAD::vector<Type> arg(1 + x.size() + xt.size());
  arg[0] = x.size();
  for(int i=0; i<x.size(); i++) { arg[1+i] = x(i); }
  for(int i=0; i<xt.size(); i++){ arg[1+x.size()+i] = xt(i);}
  return findInt(arg);
}

void order_work(const CppAD::vector<double> &tx,
                CppAD::vector<double> &ty) CSKIP({
  size_t n = tx.size();
  std::vector<std::pair<double, size_t> > y(n);
  for (size_t i=0; i<n; i++) {
    y[i].first = tx[i];
    y[i].second = i;
  }
  std::sort(y.begin(), y.end());
  for (size_t i=0; i<n; i++) {
    ty[i] = y[i].second;
  }
})

TMB_ATOMIC_VECTOR_FUNCTION(
                           // ATOMIC_NAME
                           order
                           ,
                           // OUTPUT_DIM
                           tx.size()
                           ,
                           // ATOMIC_DOUBLE
                           order_work(tx, ty);
                           ,
                           // ATOMIC_REVERSE
                           for (size_t i=0; i<px.size(); i++) px[i] = 0;
                           )

template<class Type>
vector<Type> order(vector<Type> x) {
  CppAD::vector<Type> arg(x.size());
  for(int i=0; i<x.size(); i++) { arg[i] = x(i); }
  return order(arg);
}

// subset(x, i) := x[i]  (linear sparse operator y = A * x)
// Y = subset(X, i) *OR* Y += X[i]
template <bool adjoint=false>
CppAD::vector<double> subset_work(const CppAD::vector<double> &X) {
  int ni = (int) X[0];
  int nx = (int) X[1];
  int ny = (adjoint ? nx : ni);
  CppAD::vector<double> Y(ny);
  const double* v = &X[2];
  const double* x = v + ni;
  if (!adjoint) {
    for (int i=0; i<ni; i++) {
      int i_ = (int) v[i];
      bool ok = (i_ >= 0) && (i_ < nx);
      if (ok)
        Y[i] = x[i_];
      else
        Y[i] = NA_REAL;
    }
  } else {
    for (int i=0; i<nx; i++)
      Y[i] = 0;
    for (int i=0; i<ni; i++) {
      int i_ = (int) v[i];
      bool ok = (i_ >= 0) && (i_ < nx);
      if (ok) Y[i_] += x[i];
    }
  }
  return Y;
}

// subset:     c(length(i), length(x), i, x)   -> length(i)
// subset_adj: c(length(i), length(x), i, py)  -> length(x)
template<class Type>
CppAD::vector<Type> arg_adj(const CppAD::vector<Type> &X, const CppAD::vector<Type> &Y) {
  int ni = CppAD::Integer(X[0]);
  int nkeep = 2 + ni;
  int ny = Y.size();
  CppAD::vector<Type> ans(nkeep + ny);
  for (int i=0; i<nkeep; i++) ans[i] = X[i];
  for (int i=0; i<ny; i++) ans[nkeep + i] = Y[i];
  return ans;
}

template<class Type>
void tail_set(CppAD::vector<Type> &x, const CppAD::vector<Type> &y) {
  size_t offset = x.size() - y.size();
  for (size_t i=0; i<offset; i++) {
    x[i] = 0;
  }
  for (size_t i=0; i<y.size(); i++) {
    x[offset + i] = y[i];
  }
}

TMB_ATOMIC_VECTOR_FUNCTION_DECLARE(subset)     // So can be used by subset_adj
TMB_ATOMIC_VECTOR_FUNCTION_DECLARE(subset_adj) // So can be used by subset
TMB_ATOMIC_VECTOR_FUNCTION_DEFINE( subset,     CppAD::Integer(tx[0]), ty = subset_work<0>(tx), tail_set(px, subset_adj(arg_adj(tx, py))))
TMB_ATOMIC_VECTOR_FUNCTION_DEFINE( subset_adj, CppAD::Integer(tx[1]), ty = subset_work<1>(tx), tail_set(px, subset    (arg_adj(tx, py))))

template<class Type>
vector<Type> subset(vector<Type> x, vector<Type> i) {
  bool i_constant = true;
  for (int k = 0; k < i.size(); k++)
    i_constant = i_constant && !CppAD::Variable(i[k]);
  if (i_constant) {
    vector<Type> ans(i.size());
    for (int k = 0; k < i.size(); k++) {
      int i_ = (int) asDouble(i[k]);
      bool ok = (i_ >= 0) && (i_ < x.size());
      if (ok)
        ans[k] = x[i_];
      else
        ans[k] = NA_REAL;
    }
    return ans;
  }
  CppAD::vector<Type> arg(2 + x.size() + i.size());
  arg[0] = i.size();
  arg[1] = x.size();
  for(int k=0; k<i.size(); k++) { arg[2+k] = i(k); }
  for(int k=0; k<x.size(); k++) { arg[2+i.size()+k] = x(k); }
  return subset(arg);
}

template<class Type>
vector<Type> sort(vector<Type> x) {
  return subset(x, order(x));
}

TMB_ATOMIC_STATIC_FUNCTION(
                           // ATOMIC_NAME
                           fmod
                           ,
                           // INPUT_DIM
                           2
                           ,
                           // ATOMIC_DOUBLE
                           ty[0] = std::fmod(tx[0], tx[1]);
                           ,
                           // ATOMIC_REVERSE
                           px[0] = py[0];
                           px[1] = ( (ty[0] - tx[0]) / tx[1] ) * py[0];
                           )

template<class Type>
Type fmod(Type x, Type y) {
  CppAD::vector<Type> arg(2);
  arg[0] = x; arg[1] = y;
  return fmod(arg)[0];
}

}

namespace tmbutils {

// Copyright (C) 2013-2015 Kasper Kristensen
// License: GPL-2

/*
 *  R : A Computer Language for Statistical Data Analysis
 *  Copyright (C) 1995, 1996  Robert Gentleman and Ross Ihaka
 *  Copyright (C) 1998--2001  Robert Gentleman, Ross Ihaka and the
 *                            R Development Core Team
 *
 *  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, a copy is available at
 *  http://www.r-project.org/Licenses/
 */

/*	Spline Interpolation
 *	--------------------
 *	C code to perform spline fitting and interpolation.
 *	There is code here for:
 *
 *	1. Natural splines.
 *
 *	2. Periodic splines
 *
 *	3. Splines with end-conditions determined by fitting
 *	   cubics in the start and end intervals (Forsythe et al).
 *
 *
 *	Computational Techniques
 *	------------------------
 *	A special LU decomposition for symmetric tridiagonal matrices
 *	is used for all computations, except for periodic splines where
 *	Choleski is more efficient.
 */

/**  \brief Spline Interpolation

     R's spline function wrapped into template class that can be used with TMB.

     \note Spline evaluation involves parameter dependent branching on the x-values when these are not constant. The present implementation accounts for this, but for efficiency it's important to use vectorized evaluation.
*/
template <class Type>
class splinefun
{
private:
  /* Input to R's "spline_coef" */
  int method[1];
  int n[1];
  Type *x;
  Type *y;
  Type *b;
  Type *c;
  Type *d;
  Type *e;

  bool constant_knots;
  /* Not used */
  //int errno, EDOM;

  /* Memory management helpers */
  void clear() {
    if (*n != 0) {
      delete [] x;
      delete [] y;
      delete [] b;
      delete [] c;
      delete [] d;
      delete [] e;
    }
  }
  void alloc(int n_) {
    *n = n_;
    x = new Type[*n];
    y = new Type[*n];
    b = new Type[*n];
    c = new Type[*n];
    d = new Type[*n];
    e = new Type[*n];
  }
  void copy_from(const splinefun &fun) {
    *method = fun.method[0];
    *n = fun.n[0];
    alloc(*n);
    for(int i = 0; i < *n; i++) {
      x[i] = fun.x[i];
      y[i] = fun.y[i];
      b[i] = fun.b[i];
      c[i] = fun.c[i];
      d[i] = fun.d[i];
      e[i] = fun.e[i];
    }
    constant_knots = fun.constant_knots;
  }

public:
  /* Construct spline */
  splinefun() {
    x = y = b = c = d = e = NULL;
    *method = *n = 0;
  };
  splinefun(const splinefun &fun){
    copy_from(fun);
  }
  splinefun& operator=(const splinefun &x) {
    if (this != &x) {
      (*this).clear();
      (*this).copy_from(x);
    }
    return *this;
  }
  ~splinefun() {
    clear();
  };
  /** \brief Construct spline function object
      \param x_ x values values
      \param y_ y values - same length as x.
      \param method Integer determining the spline type:
      1. Natural splines.
      2. Periodic splines
      3. Splines with end-conditions determined by fitting
         cubics in the start and end intervals (Forsythe et al).
   */
  splinefun(vector<Type> x_,
            vector<Type> y_,
            int method_ = 3) {
    method[0] = method_;
    n[0] = x_.size();
    alloc( x_.size() );
    // sort x_
    vector<Type> ord = atomic::order(x_);
    x_ = atomic::subset(x_, ord);
    y_ = atomic::subset(y_, ord);
    for(int i=0; i < *n; i++) {
      x[i] = x_[i];
      y[i] = y_[i];
    }
    constant_knots = !any_variable(x, *n);
    spline_coef(method, n, x, y, b, c, d, e);
  }

  /** \brief Evaluate spline - scalar argument case
      \note If `x` (or knots) are *variables* use vectorized evaluaton for autodiff efficiency
  */
  Type operator()(const Type &x_) {
    Type u[1];
    Type v[1];
    int nu[1];
    u[0] = x_;
    nu[0] = 1;
    spline_eval(method, nu, u, v,
		n, x, y, b, c, d);
    return v[0];
  }
  /** \brief Evaluate spline - vector argument case
      \note If `x` (or knots) are *variables* use vectorized evaluaton for autodiff efficiency
  */
  vector<Type> operator() (const vector<Type> &x_) {
    vector<Type> u(x_); // input
    vector<Type> v(x_.size()); // output
    int nu[1]; *nu = x_.size();
    spline_eval(method, nu, u.data(), v.data(),
		n, x, y, b, c, d);
    return v;
  }

  /* ------------------------------------------------------------------ 
   * The following is copy-pasted from R-package "stats" where "double"
   * has been replaced by "Type".
   *-------------------------------------------------------------------
   */

  /*
   *	Natural Splines
   *	---------------
   *	Here the end-conditions are determined by setting the second
   *	derivative of the spline at the end-points to equal to zero.
   *
   *	There are n-2 unknowns (y[i]'' at x[2], ..., x[n-1]) and n-2
   *	equations to determine them.  Either Choleski or Gaussian
   *	elimination could be used.
   */
  
  void natural_spline(int n, Type *x, Type *y, Type *b, Type *c, Type *d)
  {
    int nm1, i;
    Type t;
    
    x--; y--; b--; c--; d--;
    
    if(n < 2) {
      //errno = EDOM;
      return;
    }
    
    if(n < 3) {
      t = (y[2] - y[1]);
      b[1] = t / (x[2]-x[1]);
      b[2] = b[1];
      c[1] = c[2] = d[1] = d[2] = 0.0;
      return;
    }

    nm1 = n - 1;

    /* Set up the tridiagonal system */
    /* b = diagonal, d = offdiagonal, c = right hand side */
    
    d[1] = x[2] - x[1];
    c[2] = (y[2] - y[1])/d[1];
    for( i=2 ; i<n ; i++) {
      d[i] = x[i+1] - x[i];
      b[i] = 2.0 * (d[i-1] + d[i]);
      c[i+1] = (y[i+1] - y[i])/d[i];
      c[i] = c[i+1] - c[i];
    }
    
    /* Gaussian elimination */
    
    for(i=3 ; i<n ; i++) {
      t = d[i-1]/b[i-1];
      b[i] = b[i] - t*d[i-1];
      c[i] = c[i] - t*c[i-1];
    }
    
    /* Backward substitution */
    
    c[nm1] = c[nm1]/b[nm1];
    for(i=n-2 ; i>1 ; i--)
      c[i] = (c[i]-d[i]*c[i+1])/b[i];
    
    /* End conditions */
    
    c[1] = c[n] = 0.0;
    
    /* Get cubic coefficients */
    
    b[1] = (y[2] - y[1])/d[1] - d[i] * c[2];
    c[1] = 0.0;
    d[1] = c[2]/d[1];
    b[n] = (y[n] - y[nm1])/d[nm1] + d[nm1] * c[nm1];
    for(i=2 ; i<n ; i++) {
      b[i] = (y[i+1]-y[i])/d[i] - d[i]*(c[i+1]+2.0*c[i]);
      d[i] = (c[i+1]-c[i])/d[i];
      c[i] = 3.0*c[i];
    }
    c[n] = 0.0;
    d[n] = 0.0;
    
    return;
  }
  
  /*
   *	Splines a la Forsythe Malcolm and Moler
   *	---------------------------------------
   *	In this case the end-conditions are determined by fitting
   *	cubic polynomials to the first and last 4 points and matching
   *	the third derivitives of the spline at the end-points to the
   *	third derivatives of these cubics at the end-points.
   */
  
  void fmm_spline(int n, Type *x, Type *y, Type *b, Type *c, Type *d)
  {
    int nm1, i;
    Type t;
    
    /* Adjustment for 1-based arrays */
    
    x--; y--; b--; c--; d--;
    
    if(n < 2) {
      //errno = EDOM;
      return;
    }
    
    if(n < 3) {
      t = (y[2] - y[1]);
      b[1] = t / (x[2]-x[1]);
      b[2] = b[1];
      c[1] = c[2] = d[1] = d[2] = 0.0;
      return;
    }
    
    nm1 = n - 1;
    
    /* Set up tridiagonal system */
    /* b = diagonal, d = offdiagonal, c = right hand side */
    
    d[1] = x[2] - x[1];
    c[2] = (y[2] - y[1])/d[1];/* = +/- Inf	for x[1]=x[2] -- problem? */
    for(i=2 ; i<n ; i++) {
      d[i] = x[i+1] - x[i];
      b[i] = 2.0 * (d[i-1] + d[i]);
      c[i+1] = (y[i+1] - y[i])/d[i];
      c[i] = c[i+1] - c[i];
    }
    
    /* End conditions. */
    /* Third derivatives at x[0] and x[n-1] obtained */
    /* from divided differences */
    
    b[1] = -d[1];
    b[n] = -d[nm1];
    c[1] = c[n] = 0.0;
    if(n > 3) {
      c[1] = c[3]/(x[4]-x[2]) - c[2]/(x[3]-x[1]);
      c[n] = c[nm1]/(x[n] - x[n-2]) - c[n-2]/(x[nm1]-x[n-3]);
      c[1] = c[1]*d[1]*d[1]/(x[4]-x[1]);
      c[n] = -c[n]*d[nm1]*d[nm1]/(x[n]-x[n-3]);
    }
    
    /* Gaussian elimination */
    
    for(i=2 ; i<=n ; i++) {
      t = d[i-1]/b[i-1];
      b[i] = b[i] - t*d[i-1];
      c[i] = c[i] - t*c[i-1];
    }
    
    /* Backward substitution */
    
    c[n] = c[n]/b[n];
    for(i=nm1 ; i>=1 ; i--)
      c[i] = (c[i]-d[i]*c[i+1])/b[i];
    
    /* c[i] is now the sigma[i-1] of the text */
    /* Compute polynomial coefficients */
    
    b[n] = (y[n] - y[n-1])/d[n-1] + d[n-1]*(c[n-1]+ 2.0*c[n]);
    for(i=1 ; i<=nm1 ; i++) {
      b[i] = (y[i+1]-y[i])/d[i] - d[i]*(c[i+1]+2.0*c[i]);
      d[i] = (c[i+1]-c[i])/d[i];
      c[i] = 3.0*c[i];
    }
    c[n] = 3.0*c[n];
    d[n] = d[nm1];
    return;
  }
  
  
  /*
   *	Periodic Spline
   *	---------------
   *	The end conditions here match spline (and its derivatives)
   *	at x[1] and x[n].
   *
   *	Note: There is an explicit check that the user has supplied
   *	data with y[1] equal to y[n].
   */
  
  void periodic_spline(int n, Type *x, Type *y,
		       Type *b, Type *c, Type *d, Type *e)
  {
    Type s;
    int i, nm1;
    
    /* Adjustment for 1-based arrays */
    
    x--; y--; b--; c--; d--; e--;
    
    if(n < 2 || y[1] != y[n]) {
      errno = EDOM;
      return;
    }
    
    if(n == 2) {
      b[1] = b[2] = c[1] = c[2] = d[1] = d[2] = 0.0;
      return;
    } else if (n == 3) {
      b[1] = b[2] = b[3] = -(y[1] - y[2])*(x[1] - 2*x[2] + x[3])/(x[3]-x[2])/(x[2]-x[1]);
      c[1] = -3*(y[1]-y[2])/(x[3]-x[2])/(x[2]-x[1]);
      c[2] = -c[1];
      c[3] = c[1];
      d[1] = -2*c[1]/3/(x[2]-x[1]);
      d[2] = -d[1]*(x[2]-x[1])/(x[3]-x[2]);
      d[3] = d[1];
      return;
    }
    
    /* else --------- n >= 4 --------- */
    nm1 = n-1;
    
    /* Set up the matrix system */
    /* A = diagonal	 B = off-diagonal  C = rhs */
    
#define A	b
#define B	d
#define C	c
    
    B[1]  = x[2] - x[1];
    B[nm1]= x[n] - x[nm1];
    A[1] = 2.0 * (B[1] + B[nm1]);
    C[1] = (y[2] - y[1])/B[1] - (y[n] - y[nm1])/B[nm1];
    
    for(i = 2; i < n; i++) {
      B[i] = x[i+1] - x[i];
      A[i] = 2.0 * (B[i] + B[i-1]);
      C[i] = (y[i+1] - y[i])/B[i] - (y[i] - y[i-1])/B[i-1];
    }
    
    /* Choleski decomposition */
    
#define L	b
#define M	d
#define E	e
    
    L[1] = sqrt(A[1]);
    E[1] = (x[n] - x[nm1])/L[1];
    s = 0.0;
    for(i = 1; i <= nm1 - 2; i++) {
      M[i] = B[i]/L[i];
      if(i != 1) E[i] = -E[i-1] * M[i-1] / L[i];
      L[i+1] = sqrt(A[i+1]-M[i]*M[i]);
      s = s + E[i] * E[i];
    }
    M[nm1-1] = (B[nm1-1] - E[nm1-2] * M[nm1-2])/L[nm1-1];
    L[nm1] = sqrt(A[nm1] - M[nm1-1]*M[nm1-1] - s);
    
    /* Forward Elimination */
    
#define Y	c
#define D	c
    
    Y[1] = D[1]/L[1];
    s = 0.0;
    for(i=2 ; i<=nm1-1 ; i++) {
      Y[i] = (D[i] - M[i-1]*Y[i-1])/L[i];
      s = s + E[i-1] * Y[i-1];
    }
    Y[nm1] = (D[nm1] - M[nm1-1] * Y[nm1-1] - s) / L[nm1];
    
#define X	c
    
    X[nm1] = Y[nm1]/L[nm1];
    X[nm1-1] = (Y[nm1-1] - M[nm1-1] * X[nm1])/L[nm1-1];
    for(i=nm1-2 ; i>=1 ; i--)
      X[i] = (Y[i] - M[i] * X[i+1] - E[i] * X[nm1])/L[i];
    
    /* Wrap around */
    
    X[n] = X[1];
    
    /* Compute polynomial coefficients */
    
    for(i=1 ; i<=nm1 ; i++) {
      s = x[i+1] - x[i];
      b[i] = (y[i+1]-y[i])/s - s*(c[i+1]+2.0*c[i]);
      d[i] = (c[i+1]-c[i])/s;
      c[i] = 3.0*c[i];
    }
    b[n] = b[1];
    c[n] = c[1];
    d[n] = d[1];
    return;
  }
#undef A
#undef B
#undef C
#undef L
#undef M
#undef E
#undef Y
#undef D
#undef X
  
  void spline_coef(int *method, int *n, Type *x, Type *y,
		   Type *b, Type *c, Type *d, Type *e)
  {
    switch(*method) {
    case 1:
      periodic_spline(*n, x, y, b, c, d, e);	break;
      
    case 2:
      natural_spline(*n, x, y, b, c, d);	break;
      
    case 3:
      fmm_spline(*n, x, y, b, c, d);	break;
    }
  }

  bool any_variable(Type* x, int n) {
    bool ans = false;
    for (int i=0; i<n; i++) ans = ans || CppAD::Variable(x[i]);
    return ans;
  }
  typedef Eigen::Map<Eigen::Array<Type, -1, 1> > VMap;
  vector<Type> taped_subset(Type* x, vector<Type> i) {
    vector<Type> x_(VMap(x, *n));
    return atomic::subset(x_, i);
  }
  
  void spline_eval(int *method, int *nu, Type *u, Type *v,
		   int *n, Type *x, Type *y, Type *b, Type *c, Type *d)
  {
    /* Evaluate  v[l] := spline(u[l], ...),	    l = 1,..,nu, i.e. 0:(nu-1)
     * Nodes x[i], coef (y[i]; b[i],c[i],d[i]); i = 1,..,n , i.e. 0:(*n-1)
     */
    const int n_1 = *n - 1;
    int i, j, k, l;
    Type ul, dx, tmp;
    
    if(*method == 1 && *n > 1) { /* periodic */
      dx = x[n_1] - x[0];
      for(l = 0; l < *nu; l++) {
        v[l] = atomic::fmod(u[l]-x[0], dx);
        v[l] += CppAD::CondExpLt(v[l], Type(0), dx, Type(0));
        v[l] += x[0];
      }
    }
    else {
      for(l = 0; l < *nu; l++)
	v[l] = u[l];
    }

    if (any_variable(v, *nu) || !constant_knots) {
      // Taped spline eval
      vector<Type> v_(VMap(v, *nu));
      vector<Type> x_(VMap(x, *n));
      vector<Type> i_ = atomic::findInterval(v_, x_);
      vector<Type> xi_ = taped_subset(x, i_);
      vector<Type> yi_ = taped_subset(y, i_);
      vector<Type> bi_ = taped_subset(b, i_);
      vector<Type> ci_ = taped_subset(c, i_);
      vector<Type> di_ = taped_subset(d, i_);
      vector<Type> dxi_ = v_ - xi_;
      if (*method == 2) {
        for (l = 0; l<di_.size(); l++) {
          di_[l] = CppAD::CondExpLt(v[l], x[0], Type(0), di_[l]);
        }
      }
      VMap(v, *nu) = yi_ + dxi_ * (bi_ + dxi_ * (ci_ + dxi_ * di_));
      return;
    }
    i = 0;
    for(l = 0; l < *nu; l++) {
      ul = v[l];
      if(ul < x[i] || (i < n_1 && x[i+1] < ul)) {
	/* reset i  such that  x[i] <= ul <= x[i+1] : */
	i = 0;
	j = *n;
	do {
	  k = (i+j)/2;
	  if(ul < x[k]) j = k;
	  else i = k;
	}
	while(j > i+1);
      }
      dx = ul - x[i];
      /* for natural splines extrapolate linearly left */
      tmp = (*method == 2 && ul < x[0]) ? 0.0 : d[i];
      
      v[l] = y[i] + dx*(b[i] + dx*(c[i] + dx*tmp));
    }
  }
  
};

}