Splines

/*--------------------------------------------------------------------------*\
 |                                                                          |
 |  Copyright (C) 2016                                                      |
 |                                                                          |
 |         , __                 , __                                        |
 |        /|/  \               /|/  \                                       |
 |         | __/ _   ,_         | __/ _   ,_                                |
 |         |   \|/  /  |  |   | |   \|/  /  |  |   |                        |
 |         |(__/|__/   |_/ \_/|/|(__/|__/   |_/ \_/|/                       |
 |                           /|                   /|                        |
 |                           \|                   \|                        |
 |                                                                          |
 |      Enrico Bertolazzi                                                   |
 |      Dipartimento di Ingegneria Industriale                              |
 |      Universita` degli Studi di Trento                                   |
 |      email: enrico.bertolazzi@unitn.it                                   |
 |                                                                          |
\*--------------------------------------------------------------------------*/
 
#pragma once
 
#ifndef SPLINES_HH
#define SPLINES_HH
 
#ifdef __GNUC__
#pragma GCC diagnostic push
#endif
 
#ifdef __clang__
#pragma clang diagnostic push
#pragma clang diagnostic ignored "-Wc++98-compat"
#pragma clang diagnostic ignored "-Wglobal-constructors"
#pragma clang diagnostic ignored "-Wc++98-compat-pedantic"
#pragma clang diagnostic ignored "-Wpoison-system-directories"
#endif
 
#include "SplinesConfig.hh"
#include <fstream>
 
namespace Splines {
 
  using std::vector;
  using std::string;
  using std::exception;
  using std::runtime_error;
  using std::basic_ostream;
  using std::basic_istream;
  using std::lower_bound;
  using std::pair;
  using std::cout;
  using std::cin;
  using std::cerr;
 
  typedef double real_type; 
  typedef int    integer;   
 
  #ifndef DOXYGEN_SHOULD_SKIP_THIS
  using Malloc_real  = Utils::Malloc<real_type>;
  using ostream_type = basic_ostream<char>;
  using istream_type = basic_istream<char>;
 
  void backtrace( ostream_type & );
  #endif
 
  // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 
  using SplineType1D = enum class SplineType1D : integer {
    CONSTANT   = 0,
    LINEAR     = 1,
    CUBIC      = 2,
    AKIMA      = 3,
    BESSEL     = 4,
    PCHIP      = 5,
    QUINTIC    = 6,
    HERMITE    = 7,
    SPLINE_SET = 8,
    SPLINE_VEC = 9
  };
 
  // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 
  using SplineType2D = enum class SplineType2D : integer {
    BILINEAR  = 0,
    BICUBIC   = 1,
    BIQUINTIC = 2,
    AKIMA2D   = 3
  };
 
  #ifndef DOXYGEN_SHOULD_SKIP_THIS
  extern SplineType1D string_to_splineType1D( string const & n );
  extern SplineType2D string_to_splineType2D( string const & n );
  extern char const * to_string( SplineType2D t );
  extern char const * to_string( SplineType1D t );
  #endif
 
  using GC_namespace::GenericContainer;
  using GC_namespace::vec_real_type;
  using GC_namespace::vec_string_type;
  using GC_namespace::vector_type;
  using GC_namespace::map_type;
 
  /*
  //   _   _                     _ _
  //  | | | | ___ _ __ _ __ ___ (_) |_ ___
  //  | |_| |/ _ \ '__| '_ ` _ \| | __/ _ \
  //  |  _  |  __/ |  | | | | | | | ||  __/
  //  |_| |_|\___|_|  |_| |_| |_|_|\__\___|
  */
 
  #ifndef DOXYGEN_SHOULD_SKIP_THIS
 
  void Hermite3( real_type x, real_type H, real_type base[4] );
  void Hermite3_D( real_type x, real_type H, real_type base_D[4] );
  void Hermite3_DD( real_type x, real_type H, real_type base_DD[4] );
  void Hermite3_DDD( real_type x, real_type H, real_type base_DDD[4] );
 
  void Hermite5( real_type x, real_type H, real_type base[6] );
  void Hermite5_D( real_type x, real_type H, real_type base_D[6] );
  void Hermite5_DD( real_type x, real_type H, real_type base_DD[6] );
  void Hermite5_DDD( real_type x, real_type H, real_type base_DDD[6] );
  void Hermite5_DDDD( real_type x, real_type H, real_type base_DDDD[6] );
  void Hermite5_DDDDD( real_type x, real_type H, real_type base_DDDDD[6] );
 
  #endif
 
  static
  inline
  void
  Hermite3_to_poly(
    real_type   H,
    real_type   P0,
    real_type   P1,
    real_type   DP0,
    real_type   DP1,
    real_type & A,
    real_type & B,
    real_type & C,
    real_type & D
  ) {
    real_type H2{ H*H };
    real_type P10{ P1-P0 };
    A = (DP0+DP1-2*P10/H)/H2;
    B = (3*P10/H-(2*DP0+DP1))/H;
    C = DP0;
    D = P0;
  }
 
  static
  inline
  void
  Hermite5_to_poly(
    real_type   h,
    real_type   P0,
    real_type   P1,
    real_type   DP0,
    real_type   DP1,
    real_type   DDP0,
    real_type   DDP1,
    real_type & A,
    real_type & B,
    real_type & C,
    real_type & D,
    real_type & E,
    real_type & F
  ) {
    real_type h2{ h*h };
    real_type h3{ h*h2 };
    real_type P10{ P1-P0 };
    A = ( (DDP1-DDP0)/2+(6*P10/h-3*(DP0+DP1))/h)/h3;
    B = ( (1.5*DDP0-DDP1)+ ((8*DP0+7*DP1)-15*P10/h)/h )/h2;
    C = ( 0.5*DDP1-1.5*DDP0 + (10*P10/h -(6*DP0+4*DP1))/h )/h;
    D = DDP0/2;
    E = DP0;
    F = P0;
  }
 
  // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 
  /*
  //   ____  _ _ _
  //  | __ )(_) (_)_ __   ___  __ _ _ __
  //  |  _ \| | | | '_ \ / _ \/ _` | '__|
  //  | |_) | | | | | | |  __/ (_| | |
  //  |____/|_|_|_|_| |_|\___|\__,_|_|
  */
 
  #ifndef DOXYGEN_SHOULD_SKIP_THIS
 
  real_type
  bilinear3(
    real_type const p[4],
    real_type const M[4][4],
    real_type const q[4]
  );
 
  real_type
  bilinear5(
    real_type const p[6],
    real_type const M[6][6],
    real_type const q[6]
  );
 
  integer
  check_cubic_spline_monotonicity(
    real_type const X[],
    real_type const Y[],
    real_type const Yp[],
    integer         npts
  );
 
  #endif
 
  /*\
   |  ____                                _        _          _   _
   | |  _ \ __ _ _ __ __ _ _ __ ___   ___| |_ _ __(_)______ _| |_(_) ___  _ __
   | | |_) / _` | '__/ _` | '_ ` _ \ / _ \ __| '__| |_  / _` | __| |/ _ \| '_ \
   | |  __/ (_| | | | (_| | | | | | |  __/ |_| |  | |/ / (_| | |_| | (_) | | | |
   | |_|   \__,_|_|  \__,_|_| |_| |_|\___|\__|_|  |_/___\__,_|\__|_|\___/|_| |_|
   |
  \*/
 
  void
  uniform(
    integer         dim,
    integer         npts,
    real_type const pnts[],
    integer         ld_pnts,
    real_type       t[]
  );
 
  void
  chordal(
    integer         dim,
    integer         npts,
    real_type const pnts[],
    integer         ld_pnts,
    real_type       t[]
  );
 
  void
  centripetal(
    integer         dim,
    integer         npts,
    real_type const pnts[],
    integer         ld_pnts,
    real_type       alpha,
    real_type       t[]
  );
 
  void
  universal(
    integer         dim,
    integer         npts,
    real_type const pnts[],
    integer         ld_pnts,
    real_type       t[]
  );
 
  void
  FoleyNielsen(
    integer         dim,
    integer         npts,
    real_type const pnts[],
    integer         ld_pnts,
    real_type       t[]
  );
 
  void
  FangHung(
    integer         dim,
    integer         npts,
    real_type const pnts[],
    integer         ld_pnts,
    real_type       t[]
  );
 
  /*\
   |   ____        _ _
   |  / ___| _ __ | (_)_ __   ___
   |  \___ \| '_ \| | | '_ \ / _ \
   |   ___) | |_) | | | | | |  __/
   |  |____/| .__/|_|_|_| |_|\___|
   |        |_|
  \*/
  class Spline {
    friend class SplineSet;
  protected:
 
    #ifndef DOXYGEN_SHOULD_SKIP_THIS
 
    string m_name;
    bool   m_curve_is_closed{false};
    bool   m_curve_can_extend{true};
    bool   m_curve_extended_constant{false};
 
    integer     m_npts{0};
    integer     m_npts_reserved{0};
    real_type * m_X{nullptr}; // allocated in the derived class!
    real_type * m_Y{nullptr}; // allocated in the derived class!
 
    #ifdef SPLINES_USE_THREADS
    mutable std::mutex                                         m_last_interval_mutex;
    mutable std::map<std::thread::id,std::shared_ptr<integer>> m_last_interval;
    #else
    mutable integer m_last_interval;
    #endif
 
    void init_last_interval();
 
    Spline( Spline const & ) = delete;
    Spline const & operator = ( Spline const & ) = delete;
 
    #endif
 
  public:
 
 
    Spline( string const & name = "Spline" )
    : m_name(name)
    {
      this->init_last_interval();
    }
 
    virtual
    ~Spline()
    {}
 
 
    void search( std::pair<integer,real_type> & res ) const;
 
 
    string const & name() const { return m_name; }
 
    bool is_closed() const { return m_curve_is_closed;  }
    void make_closed() { m_curve_is_closed = true;  }
    void make_opened() { m_curve_is_closed = false; }
 
    bool is_bounded() const { return !m_curve_can_extend; }
    void make_unbounded() { m_curve_can_extend = true;  }
    void make_bounded() { m_curve_can_extend = false; }
 
    bool is_extended_constant() const { return m_curve_extended_constant;  }
    void make_extended_constant()     { m_curve_extended_constant = true;  }
    void make_extended_not_constant() { m_curve_extended_constant = false; }
 
 
 
    integer num_points() const { return m_npts; }
 
    real_type x_node( integer i ) const { return m_X[size_t(i)]; }
 
    real_type y_node( integer i ) const { return m_Y[size_t(i)]; }
 
    real_type x_begin() const { return m_X[0]; }
 
    real_type y_begin() const { return m_Y[0]; }
 
    real_type x_end() const { return m_X[size_t(m_npts-1)]; }
 
    real_type y_end() const { return m_Y[size_t(m_npts-1)]; }
 
    real_type x_min() const { return m_X[0]; }
 
    real_type x_max() const { return m_X[m_npts-1]; }
 
    real_type
    y_min() const {
      integer N{m_npts};
      if ( type() == SplineType1D::CONSTANT ) --N;
      return *std::min_element(m_Y,m_Y+N);
    }
 
    real_type
    y_max() const {
      integer N{m_npts};
      if ( type() == SplineType1D::CONSTANT ) --N;
      return *std::max_element(m_Y,m_Y+N);
    }
 
    virtual
    void
    y_min_max(
      integer   & i_min_pos,
      real_type & x_min_pos,
      real_type & y_min,
      integer   & i_max_pos,
      real_type & x_max_pos,
      real_type & y_max
    ) const;
 
    virtual
    void
    y_min_max(
      vector<integer>   & i_min_pos,
      vector<real_type> & x_min_pos,
      vector<real_type> & y_min,
      vector<integer>   & i_max_pos,
      vector<real_type> & x_max_pos,
      vector<real_type> & y_max
    ) const;
 
 
 
    void
    build( GenericContainer const & gc )
    { setup(gc); }
 
    virtual
    void
    build(
      real_type const x[], integer incx,
      real_type const y[], integer incy,
      integer n
    );
 
    inline
    void
    build(
      real_type const x[],
      real_type const y[],
      integer         n
    )
    { this->build( x, 1, y, 1, n ); }
 
    inline
    void
    build( vector<real_type> const & x, vector<real_type> const & y ) {
      integer N = integer(x.size());
      if ( N > integer(y.size()) ) N = integer(y.size());
      this->build( &x.front(), 1, &y.front(), 1, N );
    }
 
    virtual
    void build() = 0;
 
    virtual
    void setup( GenericContainer const & gc );
 
 
 
    virtual void reserve( integer npts ) = 0;
 
    void push_back( real_type x, real_type y );
 
    void drop_back() { if ( m_npts > 0 ) --m_npts; }
 
    virtual void clear() = 0;
 
 
 
    void set_origin( real_type x0 );
 
    void set_range( real_type xmin, real_type xmax );
 
 
 
    void
    dump(
      ostream_type & s,
      integer        nintervals,
      char const     header[] = "x\ty"
    ) const;
 
    void
    dump(
      char const fname[],
      integer    nintervals,
      char const header[] = "x\ty"
    ) const {
      std::ofstream file(fname);
      this->dump( file, nintervals, header );
      file.close();
    }
 
    virtual void write_to_stream( ostream_type & s ) const = 0;
 
 
 
    virtual real_type eval( real_type x ) const = 0;
 
    virtual real_type D( real_type x ) const = 0;
 
    virtual real_type DD( real_type x ) const = 0;
 
    virtual real_type DDD( real_type x ) const = 0;
 
    virtual real_type DDDD( real_type ) const { return real_type(0); }
 
    virtual real_type DDDDD( real_type ) const { return real_type(0); }
 
 
    real_type operator () ( real_type x ) const { return this->eval(x); }
    real_type eval_D( real_type x ) const { return this->D(x); }
    real_type eval_DD( real_type x ) const { return this->DD(x); }
    real_type eval_DDD( real_type x ) const { return this->DDD(x); }
    real_type eval_DDDD( real_type x ) const { return this->DDDD(x); }
    real_type eval_DDDDD( real_type x ) const { return this->DDDDD(x); }
 
 
    virtual real_type id_eval( integer ni, real_type x ) const = 0;
 
    virtual real_type id_D( integer ni, real_type x ) const = 0;
 
    virtual real_type id_DD( integer ni, real_type x ) const = 0;
 
    virtual real_type id_DDD( integer ni, real_type x ) const = 0;
 
    virtual real_type id_DDDD( integer, real_type ) const { return real_type(0); }
 
    virtual real_type id_DDDDD( integer, real_type ) const { return real_type(0); }
 
 
 
    virtual
    integer // order
    coeffs(
      real_type cfs[],
      real_type nodes[],
      bool      transpose = false
    ) const = 0;
 
    virtual integer order() const = 0;
 
    char const *
    type_name() const
    { return to_string(type()); }
 
    virtual SplineType1D type() const = 0;
 
    string info() const;
 
    void
    info( ostream_type & stream ) const
    { stream << this->info() << '\n'; }
 
 
    #ifdef SPLINES_BACK_COMPATIBILITY
    void pushBack( real_type x, real_type y ) { push_back(x,y); }
    void dropBack() { drop_back(); }
    void setOrigin( real_type x0 ) { set_origin(x0); }
    void setRange( real_type xmin, real_type xmax ) { set_range( xmin, xmax ); }
    void writeToStream( ostream_type & s ) const { write_to_stream(s); }
    #endif
 
  };
 
  real_type curvature( real_type s, Spline const & X, Spline const & Y );
 
  real_type curvature_D( real_type s, Spline const & X, Spline const & Y );
 
  real_type curvature_DD( real_type s, Spline const & X, Spline const & Y );
 
  /*\
   |    ____      _     _        ____        _ _              ____
   |   / ___|   _| |__ (_) ___  / ___| _ __ | (_)_ __   ___  | __ )  __ _ ___  ___
   |  | |  | | | | '_ \| |/ __| \___ \| '_ \| | | '_ \ / _ \ |  _ \ / _` / __|/ _ \
   |  | |__| |_| | |_) | | (__   ___) | |_) | | | | | |  __/ | |_) | (_| \__ \  __/
   |   \____\__,_|_.__/|_|\___| |____/| .__/|_|_|_| |_|\___| |____/ \__,_|___/\___|
   |                                  |_|
  \*/
  class CubicSplineBase : public Spline {
  protected:
 
    #ifndef DOXYGEN_SHOULD_SKIP_THIS
    Malloc_real m_mem_cubic;
    real_type * m_Yp{nullptr};
    bool        m_external_alloc{false};
    #endif
 
  public:
 
    #ifndef DOXYGEN_SHOULD_SKIP_THIS
    using Spline::build;
    #endif
 
    CubicSplineBase( string const & name = "CubicSplineBase" );
 
    ~CubicSplineBase() override {}
 
    void copy_spline( CubicSplineBase const & S );
 
    real_type yp_node( integer i ) const { return m_Yp[size_t(i)]; }
 
    void set_range( real_type xmin, real_type xmax );
 
    void
    reserve_external(
      integer       n,
      real_type * & p_x,
      real_type * & p_y,
      real_type * & p_dy
    );
 
    void
    y_min_max(
      integer   & i_min_pos,
      real_type & x_min_pos,
      real_type & y_min,
      integer   & i_max_pos,
      real_type & x_max_pos,
      real_type & y_max
    ) const override;
 
    void
    y_min_max(
      vector<integer>   & i_min_pos,
      vector<real_type> & x_min_pos,
      vector<real_type> & y_min,
      vector<integer>   & i_max_pos,
      vector<real_type> & x_max_pos,
      vector<real_type> & y_max
    ) const override;
 
    // --------------------------- VIRTUALS -----------------------------------
 
    real_type eval( real_type x ) const override;
    real_type D( real_type x ) const override;
    real_type DD( real_type x ) const override;
    real_type DDD( real_type x ) const override;
    real_type DDDD( real_type ) const override { return 0; }
    real_type DDDDD( real_type ) const override { return 0; }
 
    real_type id_eval( integer ni, real_type x ) const override;
    real_type id_D( integer ni, real_type x ) const override;
    real_type id_DD( integer ni, real_type x ) const override;
    real_type id_DDD( integer ni, real_type x ) const override;
    real_type id_DDDD( integer, real_type ) const override { return 0; }
    real_type id_DDDDD( integer, real_type ) const override { return 0; }
 
    void write_to_stream( ostream_type & s ) const override;
 
    // --------------------------- VIRTUALS -----------------------------------
 
 
    void
    build(
      real_type const x[],  integer incx,
      real_type const y[],  integer incy,
      real_type const yp[], integer incyp,
      integer n
    );
 
    inline
    void
    build(
      real_type const x[],
      real_type const y[],
      real_type const yp[],
      integer n
    ) {
      this->build( x, 1, y, 1, yp, 1, n );
    }
 
    void
    build(
      vector<real_type> const & x,
      vector<real_type> const & y,
      vector<real_type> const & yp
    );
 
    void reserve( integer npts ) override;
 
 
    void clear() override;
 
    integer // order
    coeffs(
      real_type cfs[],
      real_type nodes[],
      bool      transpose = false
    ) const override;
 
    integer order() const override;
 
    #ifdef SPLINES_BACK_COMPATIBILITY
    void copySpline( CubicSplineBase const & S ) { this->copy_spline(S); }
    integer numPoints() const { return m_npts; }
    real_type xNode( integer i ) const { return m_X[size_t(i)]; }
    real_type yNode( integer i ) const { return m_Y[size_t(i)]; }
    real_type ypNode( integer i ) const { return this->yp_node(i); }
    real_type xBegin() const { return m_X[0]; }
    real_type yBegin() const { return m_Y[0]; }
    real_type xEnd() const { return m_X[size_t(m_npts-1)]; }
    real_type yEnd() const { return m_Y[size_t(m_npts-1)]; }
    real_type xMin() const { return m_X[0]; }
    real_type xMax() const { return m_X[m_npts-1]; }
    real_type yMin() const { return y_min(); }
    real_type yMax() const { return y_max(); }
    #endif
 
  };
 
  /*\
   |   ____        _ _            ____              __
   |  / ___| _ __ | (_)_ __   ___/ ___| _   _ _ __ / _|
   |  \___ \| '_ \| | | '_ \ / _ \___ \| | | | '__| |_
   |   ___) | |_) | | | | | |  __/___) | |_| | |  |  _|
   |  |____/| .__/|_|_|_| |_|\___|____/ \__,_|_|  |_|
   |        |_|
  \*/
 
  class SplineSurf {
 
    SplineSurf( SplineSurf const & ) = delete; // block copy constructor
    SplineSurf const & operator = ( SplineSurf const & ) = delete; // block copy method
 
    Malloc_real m_mem;
 
  protected:
 
    #ifndef DOXYGEN_SHOULD_SKIP_THIS
 
    string const m_name;
    bool         m_x_closed{false};
    bool         m_y_closed{false};
    bool         m_x_can_extend{true};
    bool         m_y_can_extend{true};
 
    integer      m_nx{0};
    integer      m_ny{0};
 
    real_type *  m_X{nullptr};
    real_type *  m_Y{nullptr};
    real_type *  m_Z{nullptr};
 
    real_type    m_Z_min{0};
    real_type    m_Z_max{0};
 
    #ifdef SPLINES_USE_THREADS
    mutable std::mutex                                         m_last_interval_x_mutex;
    mutable std::mutex                                         m_last_interval_y_mutex;
    mutable std::map<std::thread::id,std::shared_ptr<integer>> m_last_interval_x;
    mutable std::map<std::thread::id,std::shared_ptr<integer>> m_last_interval_y;
    #else
    mutable integer m_last_interval_x;
    mutable integer m_last_interval_y;
    #endif
 
    integer search_x( real_type & x ) const;
    integer search_y( real_type & y ) const;
 
    void init_last_interval_x();
    void init_last_interval_y();
 
    integer
    ipos_C( integer i, integer j, integer ldZ ) const
    { return i*ldZ + j; }
 
    integer
    ipos_F( integer i, integer j, integer ldZ ) const
    { return i + ldZ*j; }
 
    integer
    ipos_C( integer i, integer j ) const
    { return this->ipos_C(i,j,m_ny); }
 
    integer
    ipos_F( integer i, integer j ) const
    { return this->ipos_F(i,j,m_nx); }
 
    virtual void make_spline() = 0;
 
    void
    load_Z(
      real_type const z[],
      integer         ldZ,
      bool            fortran_storage,
      bool            transposed
    );
 
    #endif
 
  public:
 
    SplineSurf( string const & name = "Spline" )
    : m_mem("SplineSurf")
    , m_name(name)
    {
      this->init_last_interval_x();
      this->init_last_interval_y();
    }
 
    virtual
    ~SplineSurf();
 
 
    bool is_x_closed() const { return m_x_closed; }
 
    void make_x_closed() { m_x_closed = true; }
 
    void make_x_opened() { m_x_closed = false; }
 
    bool is_y_closed() const { return m_y_closed; }
 
    void make_y_closed() { m_y_closed = true; }
 
    void make_y_opened() { m_y_closed = false; }
 
    bool is_x_bounded() const { return m_x_can_extend; }
 
    void make_x_unbounded() { m_x_can_extend = true; }
 
    void make_x_bounded() { m_x_can_extend = false; }
 
    bool is_y_bounded() const { return m_y_can_extend; }
 
    void make_y_unbounded() { m_y_can_extend = true; }
 
    void make_y_bounded() { m_y_can_extend = false; }
 
 
    void clear();
 
 
    string const & name() const { return m_name; }
 
    integer num_point_x() const { return m_nx; }
 
    integer num_point_y() const { return m_ny; }
 
    real_type x_node( integer i ) const { return m_X[size_t(i)]; }
 
    real_type y_node( integer i ) const { return m_Y[size_t(i)]; }
 
    real_type
    z_node( integer i, integer j ) const
    { return m_Z[size_t(this->ipos_C(i,j))]; }
 
    real_type x_min() const { return m_X[0]; }
 
    real_type x_max() const { return m_X[m_nx-1]; }
 
    real_type y_min() const { return m_Y[0]; }
 
    real_type y_max() const { return m_Y[m_ny-1]; }
 
    real_type z_min() const { return m_Z_min; }
 
    real_type z_max() const { return m_Z_max; }
 
 
 
    void
    build(
      real_type const x[], integer incx,
      real_type const y[], integer incy,
      real_type const z[], integer ldZ,
      integer nx, integer ny,
      bool fortran_storage = false,
      bool transposed      = false
    );
 
    void
    build(
      vector<real_type> const & x,
      vector<real_type> const & y,
      vector<real_type> const & z,
      bool fortran_storage = false,
      bool transposed      = false
    ) {
      bool xyswp = (fortran_storage && transposed) ||
                   (!fortran_storage && !transposed);
      this->build(
        &x.front(), 1,
        &y.front(), 1,
        &z.front(), integer(xyswp ? y.size() : x.size()),
        integer(x.size()), integer(y.size()),
        fortran_storage, transposed
      );
    }
 
    void
    build(
      real_type const z[],
      integer         ldZ,
      integer         nx,
      integer         ny,
      bool fortran_storage = false,
      bool transposed      = false
    );
 
    void
    build(
      vector<real_type> const & z,
      integer                   nx,
      integer                   ny,
      bool fortran_storage = false,
      bool transposed      = false
    ) {
      this->build( &z.front(), nx, ny, fortran_storage ? nx : ny, fortran_storage, transposed );
    }
 
    void
    setup( GenericContainer const & gc );
 
    void
    build( GenericContainer const & gc )
    { setup(gc); }
 
 
 
    virtual
    real_type
    eval( real_type x, real_type y ) const = 0;
 
    virtual
    void
    D( real_type x, real_type y, real_type d[3] ) const = 0;
 
    virtual
    real_type
    Dx( real_type x, real_type y ) const = 0;
 
    virtual
    real_type
    Dy( real_type x, real_type y ) const = 0;
 
    virtual
    void
    DD( real_type x, real_type y, real_type dd[6] ) const = 0;
 
    virtual
    real_type
    Dxx( real_type x, real_type y ) const = 0;
 
    virtual
    real_type
    Dxy( real_type x, real_type y ) const = 0;
 
    virtual
    real_type
    Dyy( real_type x, real_type y ) const = 0;
 
    real_type
    operator () ( real_type x, real_type y ) const
    { return this->eval(x,y); }
 
    real_type
    eval_D_1( real_type x, real_type y ) const
    { return this->Dx(x,y); }
 
    real_type
    eval_D_2( real_type x, real_type y ) const
    { return this->Dy(x,y); }
 
    real_type
    eval_D_1_1( real_type x, real_type y ) const
    { return this->Dxx(x,y); }
 
    real_type
    eval_D_1_2( real_type x, real_type y ) const
    { return this->Dxy(x,y); }
 
    real_type
    eval_D_2_2( real_type x, real_type y ) const
    { return this->Dyy(x,y); }
 
 
    virtual void write_to_stream( ostream_type & s ) const = 0;
 
    virtual char const * type_name() const = 0;
 
    virtual string info() const;
 
    void
    info( ostream_type & stream ) const
    { stream << this->info() << '\n'; }
 
    void dump_data( ostream_type & s ) const;
 
    #ifdef SPLINES_BACK_COMPATIBILITY
    integer numPointX() const { return m_nx; }
    integer numPointY() const { return m_ny; }
    real_type xNode( integer i ) const { return m_X[size_t(i)]; }
    real_type yNode( integer i ) const { return m_Y[size_t(i)]; }
    real_type zNode( integer i, integer j ) const { return z_node(i,j); }
    real_type xMin() const { return this->x_min(); }
    real_type xMax() const { return this->x_max(); }
    real_type yMin() const { return this->y_min(); }
    real_type yMax() const { return this->y_max(); }
    real_type zMin() const { return m_Z_min; }
    real_type zMax() const { return m_Z_max; }
    void writeToStream( ostream_type & s ) const { write_to_stream(s); }
    #endif
 
  };
 
}
 
#include "Splines/SplineAkima.hxx"
#include "Splines/SplineBessel.hxx"
#include "Splines/SplineConstant.hxx"
#include "Splines/SplineLinear.hxx"
#include "Splines/SplineCubic.hxx"
#include "Splines/SplineHermite.hxx"
#include "Splines/SplinePchip.hxx"
#include "Splines/SplineQuinticBase.hxx"
#include "Splines/SplineQuintic.hxx"
#include "Splines/SplineBilinear.hxx"
#include "Splines/SplineBiCubic.hxx"
#include "Splines/SplineAkima2D.hxx"
#include "Splines/SplineBiQuintic.hxx"
 
#include "Splines/SplineVec.hxx"
#include "Splines/SplineSet.hxx"
#include "Splines/Splines1D.hxx"
#include "Splines/Splines2D.hxx"
 
#ifndef DOXYGEN_SHOULD_SKIP_THIS
 
namespace SplinesLoad {
 
  using Splines::Spline;
  using Splines::CubicSplineBase;
  using Splines::CubicSpline;
  using Splines::AkimaSpline;
  using Splines::BesselSpline;
  using Splines::PchipSpline;
  using Splines::LinearSpline;
  using Splines::ConstantSpline;
  using Splines::QuinticSpline;
  using Splines::Spline1D;
 
  using Splines::BilinearSpline;
  using Splines::BiCubicSpline;
  using Splines::BiQuinticSpline;
  using Splines::Akima2Dspline;
  using Splines::Spline2D;
 
  using Splines::SplineVec;
  using Splines::SplineSet;
 
  using Splines::SplineType1D;
  using Splines::SplineType2D;
}
 
#endif
 
#ifdef __GNUC__
#pragma GCC diagnostic pop
#endif
#ifdef __clang__
#pragma clang diagnostic pop
#endif
 
#endif