RE: A recursive approach to construction of Bi, Tri and NcubicSpline classes.

[RGProlog]

Hi Quantlib.

Here is a proposal for stepping up from class NcubicSpline class to Npus1cubicSpline class for N >= 2 in a pseudo recursive pattern as well as
CubicSpline, BicubicSpline and TricubicSpline classes in flesh and blood.

One of the cornerstones to the proposed design is the CubicSpline class which derives from 2 base classes, Qspline, and QSplin implemented in the file QSpline.hpp below.

///////////////////////////////////////////////////////////////////////////////////////////////////////

// QSpline.hpp

// Copyright (C) 2002, 2003 Roman Gitlin

#ifndef q_spline_h
#define q_spline_h

#define NOMINMAX

#include <cmath>
#include <vector>
#include <algorithm>
#include <functional>
#include <stdexcept>
#include <iostream>

#include <config.msvc.hpp>

#define LN 0.99e30

namespace QuantLib {

    namespace Math {

template <class X, class Y> class CubicSpline;

template<class X, class Y>
class QBase :
public std::unary_function <QL_ITERATOR_TRAITS<X>::value_type,
                    QL_ITERATOR_TRAITS<Y>::value_type> {

typedef std::vector<result_type>::iterator iterator;
};

template<class X, class Y>
class QSplint : public QBase<X, Y> {
friend  class CubicSpline<X, Y>;
QSplint(const X &x,
const Y &y,
const iterator &y2,
int n)
: x_(x), y_(y), y2_(y2), i_(0), n_(n) {}
result_type result(argument_type) const;
          virtual ~QSplint() {}
// data
const X &x_;
const Y &y_;
const iterator &y2_;
const int n_;
mutable int i_;
};

template<class X, class Y>
QSplint<X, Y>::result_type
QSplint<X, Y>::result(argument_type x) const {
if(!(x >= x_[i_] && x < x_[i_ + 1]))
i_ = x >= x_[i_ + 1] ?
    std::upper_bound(x_ + i_ + 1, x_ + n_, x) - x_ - 1 :
        std::upper_bound(x_, x_ + i_ + 1, x) - x_ - 1;
argument_type h = x_[i_ + 1] - x_[i_];
result_type   a = (x_[i_ + 1] - x) / h,
          b = (x - x_[i_]) / h;
return a * y_[i_] + b * y_[i_ + 1] +
   ((a * a * a - a) * y2_[i_] +
       (b * b * b - b) * y2_[i_ + 1]) * (h * h) / 6.0;
}

template<class X, class Y>
class QSpline  :  public QBase<X, Y> {
friend  class CubicSpline<X, Y>;
QSpline(){}
QSpline(const X &,
const Y &,
const iterator &y2,
int,
bool,
result_type,
result_type);
        virtual ~QSpline() {}
};

template<class X, class Y>
QSpline<X, Y>::QSpline(const X &x,
   const Y &y,
   const std::vector<result_type>::iterator &y2,
   int n,
   bool valid_argument,
   result_type y1,
   result_type yn) {
try {
if(valid_argument) {
std::vector<result_type> bsub(n - 1);
argument_type b = 0.0, s = x[1] - x[0], t = x[2] - x[1];
result_type   w = 0.0, u = (y[1] - y[0]) / s,
  v = (y[2] - y[1]) / t;
                        y2[0] = y1 >= LN ? 0.0 : -0.5;
                        bsub[0] = y1 >= LN ? 0.0 : 3.0 * (v - y1) / t;
for(int k = 1; k < n - 1; ++k) {
v = (y[k + 1] - y[k]) / (t = x[k + 1] - x[k]);
w = s / (s + t);
y2[k] = (w - 1.0) / (b = w * y2[k - 1] + 2.0);
bsub[k] = (6.0 * (v - u) / (s + t) - w * bsub[k - 1]) / b;
s = t;
u = v;
}
                    t = yn >= LN ? 0.0 : 3.0 * (yn - v) / s;
                    s = yn >= LN ? 0.0 : 0.5;
y2[n - 1] = (t - s * bsub[n - 2]) /
(s * y2[n - 2] + 1.0);
for(int q = n - 1; q > 0; --q)
y2[q - 1] = y2[q - 1] * y2[q] + bsub[q - 1];
}
} catch(std::exception &e) {
std::cerr << "QSpline:: " << e.what() << '\n';
throw;
}
}

}

}

/////////////////////////////////////////////////////////////////////////////////////////////////////

The QSpline and QSplint classes are implemented along the lines of qspline
and qsplint procedured found in Numerical Recepies in C (pp. 115, 116)
the only difference being that the QSpline clas eliminates some of the calculational redunduncies of the qspline procedure and the QSplint clas does optimize some of the binary searches out, The difference in how the quspline procedureis implemented in the existing and the proposed design aside, the key difference between the proposed  and the existing design is that the ubicSpline class as implemented below does not hold any data types by value except for a
couple of int and a couple of bool types. The rest of its data members are
costant references. The obvious downside of this approach is that a container 
holding the solution for a tridiagonal matrix system of equations performed in CubicSpline ctor must be set up outside of the CubicSpline class's scope. This approach's upside, though, is that while to construct the current CubicSpline class object of size n 3 vectors of size n, and 6 objects of the class Array of the same size have to be created, processed, and some of them passed or returned by value during the construction, the proposed CubicSpline class's ctor creates a single local vector object of size n (see QSpline ctor above; it does not hold any data on its own).


/////////////////////////////////////////////////////////////////////////////////////////////////////

#endif
// CubicSpline.hpp

// Copyright (C) 2002, 2003 Roman Gitlin

#ifndef cubic_spline_h
#define cubic_spline_h

#include <QSpline\QSpline.hpp>

using std::vector;

namespace QuantLib {

    namespace Math {

        template <class I1, class I2>
        class CubicSpline  : public QSpline<I1, I2>, QSplint<I1, I2> {
public:
            CubicSpline(const I1 &,
                        const I2 &,
const iterator &y2_,
            int,
            bool,
            result_type = LN,
            result_type = LN);
            CubicSpline(int,
        const I1 &,
                           const I2 &,
   const iterator &y2_,
             int,
             bool);
             result_type operator()(argument_type x) const;
         virtual ~CubicSpline() {}
        private:
bool is_valid_arg(const I1 &, int) const;
// data
mutable bool is_valid_data_;
const bool allow_extrapolation_;
        };

        template <class I1, class I2>
        CubicSpline<I1,I2>::CubicSpline(const I1& x,
                                      const I2& y,
const iterator &y2,
int n,
bool allow_extrapolation,
result_type y1,
result_type yn)
        :    QSpline<I1, I2>(x, y, y2, n, is_valid_arg(x, n), y1, yn),
    QSplint<I1, I2>(x, y, y2, n),
allow_extrapolation_(allow_extrapolation) {
try {
if(n < 4)
throw std::invalid_argument("Not enough points for cubic "
"spline interpolation.");
    if(!is_valid_data_)
throw std::invalid_argument("Invalid data.");
} catch(std::exception &e) {
std::cerr << "CubicSpline<I1, I2>(I1, I2, "
         "iterator, bool, int):: "
      << e.what() << '\n';
throw;
}
        }

        template <class I1, class I2>
        CubicSpline<I1,I2>::CubicSpline(int,
                                const I1& x,
                                const I2& y,
                const iterator &y2,
int n,
bool allow_extrapolation_)
        :    QSplint<I1,I2>(x, y, y2, n),
    allow_extrapolation_(allow_extrapolation_) {
try {
if(n < 4)
throw std::invalid_argument("Not enough points for cubic "
"spline interpolation.");
    if(!is_valid_arg(x, n))
throw std::invalid_argument("Invalid data.");
} catch(std::exception &e) {
std::cerr << "CubicSpline<I1, I2>(I1, I2, "
         "iterator, bool, int):: "
      << e.what() << '\n';
throw;
}
}

        template <class I1, class I2>
        CubicSpline<I1,I2>::result_type
        CubicSpline<I1,I2>::operator()(argument_type x) const {
try{
if (x < x_[0]) {
if(!allow_extrapolation_)
throw std::range_error("CubicSpline::operator() : "
   "extrapolation not allowed "
   "[x<xMin]");
return result(x_[0]);
} else if (x > x_[n_ - 1]) {
if(!allow_extrapolation_)
throw std::range_error("CubicSpline::operator() : "
   "extrapolation not allowed "
   "[x>xMax]");
return result(x_[n_ - 1]);
} else return result(x);
} catch(std::exception &e) {
std::cerr << "CubicSpline::operator() : " << e.what() << '\n';
throw;
}
        }

        template <class I1, class I2>
bool CubicSpline<I1,I2>::is_valid_arg(const I1 &x, int n) const {
bool b = true;
for(int i = 0; i < n - 1; ++i)
if(!(b = x[i] < x[i + 1])) break;
return is_valid_data_ = b && n > 3;
}

    }

}

#endif

/////////////////////////////////////////////////////////////////////////////////////////////////////

While the proposed CubicSpline class can perfectly function on its own with the
first of its 2 ctors, it is its second ctor that  holds the key to the class being a bulding block for the BicubicSpline class, which in turn, agin by virtue of having a second  specifically modified constructor, becomes a building bolock for the TricubicSpline class, and so on, and so on. The process is implemented twice for the BicubicSpline and TricubicSpline classes in the files below and is fairly straightforward.

The thing to notice here is that the CubicSpline inherent efficiency of design,
as well as binary search optimization,  its data checking and  its extrapolation handling capacity, is passed on to and is fully adequate for  the BicubicSpline, TricubicSpline, and NcubicSpline classes.


// BicubicSpline.hpp

// Copyright (C) 2002, 2003 Roman Gitlin

#ifndef bicubicspline_hpp
#define bicubicspline_hpp

#include <CubicSpline\CubicSpline.hpp>

namespace QuantLib {

    namespace Math {

template<class X1, class X2, typename X3>
class BicubicSpline :
public std::binary_function                                      <QL_ITERATOR_TRAITS<X1>::value_type,
      QL_ITERATOR_TRAITS<X2>::value_type,
      X3>
{
public:
typedef typename
std::vector<result_type>::iterator iterator;
typedef typename
std::vector<result_type>::const_iterator const_iterator;
typedef typename
std::vector<std::vector<result_type> > data_table;
BicubicSpline(int,
  const X1 &,
                          const X2 &,
  const data_table &,
  const data_table &,
  int,
  int,
  bool);
BicubicSpline(const X1 &,
  const X2 &,
  const data_table &,
  const data_table &,
  int,
  int,
  bool);
virtual ~BicubicSpline(){}
result_type operator()(first_argument_type,
   second_argument_type) const;
private:
// data
const int m_;
const int n_;
const X1 &x1_;                                            
const X2 &x2_;                                            
const data_table &y_;                                     
const data_table &y2_;                                   
bool allow_extrapolation_;
};

template<class X1, class X2, typename X3>
BicubicSpline<X1, X2, X3>::BicubicSpline(const X1 &x1,
                    const X2 &x2,
                    const data_table &y,
                    const data_table &y2,
               int m,
               int n,
               bool extrapolation)
:   m_(m), n_(n), x1_(x1), x2_(x2), y_(y), y2_(y2),
    allow_extrapolation_(extrapolation) {
try {
for(int i = 0; i < m_; ++i)
CubicSpline<X2, const_iterator>(x2_,
                         y_[i].begin(),
          (iterator)y2_[i].begin(),
                  n_,
                          allow_extrapolation_);
}
catch(std::exception &e) {
std::cerr << "BicubicSpline<X1, X2, X3>::BicubicSpline(const "
    "X1 &, const X2 &, const std::vector<std::vector"
    << "<X3> > &, int, int, bool) : "
    <<  e.what() << std::endl;
throw;
}
}

template<class X1, class X2, typename X3>
BicubicSpline<X1, X2, X3>::BicubicSpline(int,
                                                               const X1 &x1,
                    const X2 &x2,
                    const data_table &y,
                    const data_table &y2,
                int m,
                int n,
                 bool extrapolation)
:   m_(m), n_(n), x1_(x1), x2_(x2), y_(y), y2_(y2),
    allow_extrapolation_(extrapolation) {}

template<class X1, class X2, typename X3>
BicubicSpline<X1, X2, X3>::result_type
BicubicSpline<X1, X2, X3>::operator()(first_argument_type x,
               second_argument_type y) const {
try {
std::vector<result_type> temp1(m_),
                               temp2(m_);
iterator  ya = temp1.begin(),
        y2 = temp2.begin();
for(int i = 0; i < m_; ++i)
ya[i] = CubicSpline<X2, const_iterator>(0,
                                                          x2_,
          y_[i].begin(),
         (iterator)y2_[i].begin(),
          n_,
         allow_extrapolation_)(y);
return CubicSpline<X1, const_iterator>(x1_,
                                                     ya,
               y2,
                m_,
                allow_extrapolation_)(x);
}
catch(std::exception &e) {
std::cerr << "BiCubicInterpolation<X1, X2, X3>::operator() : "
      << e.what() << std::endl;
throw;
}
}

}

}

#endif


// TricubicSpline.hpp

// Copyright (C) 2002, 2003 Roman Gitlin

#ifndef tricubicspline_hpp
#define tricubicspline_hpp

#include <BicubicSpline\BicubicSpline.hpp>

namespace QuantLib {

    namespace Math {

template<class X1, class X2, typename X3, class X4>
class TricubicSpline
{
public:
typedef typename
QL_ITERATOR_TRAITS<X1>::value_type first_argument_type;
typedef typename
QL_ITERATOR_TRAITS<X2>::value_type second_argument_type;
typedef typename
QL_ITERATOR_TRAITS<X3>::value_type third_argument_type;
typedef typename
X4 result_type;
typedef typename
std::vector<result_type>::iterator iterator;
typedef typename
std::vector<result_type>::const_iterator const_iterator;
typedef typename
std::vector<std::vector<result_type> > D2data_table;
typedef typename
std::vector<D2data_table> D3data_table;
TricubicSpline(int,
   const X1 &,
   const X2 &,
   const X3 &,
   const D3data_table &,
   const D3data_table &,
   int,
   int,
   int,
   bool);
TricubicSpline(const X1 &,
   const X2 &,
   const X3 &,
   const D3data_table &,
   const D3data_table &,
   int,
   int,
   int,
   bool);
virtual ~TricubicSpline(){}
result_type operator()(first_argument_type,
   second_argument_type,
   third_argument_type) const;
private:
// data
const int m_;
const int n_;
const int o_;
const X1 &x1_;                                            
const X2 &x2_;                                            
const X3 &x3_;                                            
const D3data_table &y_;                                     
const D3data_table &y2_;                                   
bool allow_extrapolation_;
};

template<class X1, class X2, class X3, typename X4>
TricubicSpline<X1, X2, X3, X4>::TricubicSpline(const X1 &x1,
           const X2 &x2,
           const X3 &x3,
           const D3data_table &y,
           const D3data_table &y2,
       int m,
       int n,
       int o,
       bool extrapolation)
:   m_(m), n_(n), o_(o), x1_(x1), x2_(x2), x3_(x3), y_(y), y2_(y2),
    allow_extrapolation_(extrapolation) {
try {
for(int i = 0; i < m_; ++i)
BicubicSpline<X2, X3, X4>(x2_,
                                        x3_,
                y_[i],
           y2_[i],
            n_,
            o_,
            allow_extrapolation_);
}
catch(std::exception &e) {
std::cerr << "TricubicSpline<X1, X2, X3, X4>::TricubicSpline"
"(const X1 &, const X2 &, const X3 &, const std::vector<"
     "std::vector<std::vector<X3> > >&, const std::vector<"
"std::vector<std::vector<X3> > >&,  int, int, int, bool) : "
  <<  e.what() << std::endl;
throw;
}
}

template<class X1, class X2, class X3, typename X4>
TricubicSpline<X1, X2, X3, X4>::TricubicSpline(int,
                                         const X1 &x1,
     const X2 &x2,
     const X3 &x3,
     const D3data_table &y,
     const D3data_table &y2,
int m,
int n,
int o,
bool extrapolation)
:   m_(m), n_(n), o_(o), x1_(x1), x2_(x2), x3_(x3), y_(y), y2_(y2),
    allow_extrapolation_(extrapolation) {}

template<class X1, class X2, class X3, typename X4>
TricubicSpline<X1, X2, X3, X4>::result_type
TricubicSpline<X1, X2, X3, X4>::operator()(first_argument_type x,
                  second_argument_type y,
                  third_argument_type z) const {
try {
std::vector<result_type> temp1(m_),
                     temp2(m_);
iterator ya = temp1.begin(),
     y2 = temp2.begin();
for(int i = 0; i < m_; ++i)
ya[i] = BicubicSpline<X2, X3, X4>(0,
                                                    x2_,
                                                    x3_,
                            y_[i],
                       y2_[i],
                        n_,
                        o_,
              allow_extrapolation_)(y, z);

return CubicSpline<X1, const_iterator>(x1_,
                                               ya,
         y2,
         m_,
         allow_extrapolation_)(x);
}
catch(std::exception &e) {
std::cerr << "BiCubicInterpolation<X1, X2, X3>::operator() : "
      << e.what() << std::endl;
throw;
}
}

}

}

#endif