BSpline.hpp

#pragma once
#include "LinAlg/LinAlg.hpp"
#include "qip/Vector.hpp"
#include <gsl/gsl_bspline.h>
#include <gsl/gsl_version.h>
#include <iostream>
#include <numeric>
#include <vector>
 
// This should make code work with old versions of GSL
// See below for documentation of old GSL function
#ifdef GSL_MAJOR_VERSION
#if GSL_MAJOR_VERSION == 1
#define GSL_VERSION_1
#endif
// GSL 2.8 introduces substantial changes to bspline library.
#if GSL_MAJOR_VERSION <= 2 && GSL_MINOR_VERSION < 8
#define GSL_VERSION_PRIOR_2_8
#endif
#endif
 
class BSpline {
 
  // K is _order_ (degree is K-1)
  std::size_t m_K;
  // Total number of splines (i=0,1,...,N-1)
  std::size_t m_N;
  // Splines defined over range [0,x_max]
  double m_xmax;
  // Worskspace used by GSL library
  gsl_bspline_workspace *gsl_bspl_work{nullptr};
 
#ifdef GSL_VERSION_1
  // Worskspace used by old version of GSL library
  gsl_bspline_deriv_workspace *gsl_bspl_deriv_work{nullptr};
#endif
 
public:
  enum class KnotDistro { logarithmic, linear, loglinear };
  //============================================================================
  BSpline(std::size_t n, std::size_t k, double x0, double xmax,
          KnotDistro kd = KnotDistro::logarithmic)
      : m_K(k), m_N(n), m_xmax(xmax) {
    assert(m_N > m_K && "Require N>K for B-splines");
    assert(xmax > x0 && x0 > 0.0 && xmax > 0.0 && "xmax>x0 and both positive");
 
    // Allocate memory for the GSL workspace:
    const std::size_t n_break = m_N - m_K + 2;
    gsl_bspl_work = gsl_bspline_alloc(m_K, n_break);
 
#ifdef GSL_VERSION_1
    // Worskspace used by old version of GSL library
    gsl_bspl_deriv_work = gsl_bspline_deriv_alloc(m_K);
#endif
 
    set_knots(x0, xmax, n_break, kd);
 
    // Consistancy check:
#ifdef GSL_VERSION_PRIOR_2_8
    // gsl_bspline_ncoeffs will be deprecated in future GSL
    assert(gsl_bspline_ncoeffs(gsl_bspl_work) == m_N);
#else
    assert(gsl_bspline_ncontrol(gsl_bspl_work) == m_N);
#endif
    assert(gsl_bspline_order(gsl_bspl_work) == m_K);
    assert(gsl_bspline_nbreak(gsl_bspl_work) == n_break);
  }
 
  // delete copy assignment (no simple way to copy gsl_bspl_work)
  BSpline &operator=(const BSpline &) = delete;
  // delete copy constructor (no simple way to copy gsl_bspl_work)
  BSpline(const BSpline &) = delete;
 
  // Desctructor
  ~BSpline() {
    gsl_bspline_free(gsl_bspl_work);
#ifdef GSL_VERSION_1
    gsl_bspline_deriv_free(gsl_bspl_deriv_work);
#endif
  }
 
  //============================================================================
  std::size_t K() { return m_K; }
  std::size_t d() { return m_K - 1; }
  std::size_t N() { return m_N; }
 
  //============================================================================
  std::size_t find_i0(double x) {
    for (std::size_t i = 0; i < m_N; ++i) {
      if (gsl_vector_get(gsl_bspl_work->knots, i + m_K) > x) {
        return i;
      }
    }
    return m_N;
  }
 
  //============================================================================
  std::vector<double> get(double x) {
    std::vector<double> b(m_N);
    gsl_vector_view b_gsl = gsl_vector_view_array(b.data(), b.size());
    if (x >= 0.0 && x <= m_xmax) {
 
#ifdef GSL_VERSION_PRIOR_2_8
      // gsl_bspline_eval will be deprecated in future GSL
      gsl_bspline_eval(x, &b_gsl.vector, gsl_bspl_work);
#else
      gsl_bspline_eval_basis(x, &b_gsl.vector, gsl_bspl_work);
#endif
    }
    return b;
  }
 
  //============================================================================
 
  std::pair<std::size_t, LinAlg::Matrix<double>>
  get_nonzero(double x, std::size_t n_deriv) {
    std::pair<std::size_t, LinAlg::Matrix<double>> out{0, {m_K, n_deriv + 1}};
 
    if (x >= 0.0 && x <= m_xmax) {
      // outside this range, spline set to zero by default. Invalid to call
      // gsl_bspline_deriv_eval_nonzero() outside spline range
 
      auto &i0 = out.first;
      auto &bij = out.second;
      gsl_matrix_view b_gsl = bij.as_gsl_view();
 
#ifdef GSL_VERSION_1
      size_t i_end{};
      gsl_bspline_deriv_eval_nonzero(x, n_deriv, &b_gsl.matrix, &i0, &i_end,
                                     gsl_bspl_work, gsl_bspl_deriv_work);
#elif defined GSL_VERSION_PRIOR_2_8
      size_t i_end{};
      gsl_bspline_deriv_eval_nonzero(x, n_deriv, &b_gsl.matrix, &i0, &i_end,
                                     gsl_bspl_work);
#else
      gsl_bspline_basis_deriv(x, n_deriv, &b_gsl.matrix, &i0, gsl_bspl_work);
#endif
    }
 
    return out;
  }
 
  //============================================================================
  void print_knots() const {
 
    std::size_t n_cols_to_print = 4;
    std::cout << "Knots:\n ";
    std::size_t count = 0;
    while (count < m_N + m_K) {
      auto t = gsl_bspl_work->knots->data[count];
      printf("%.5e, ", t);
      ++count;
      if (count % n_cols_to_print == 0)
        std::cout << "\n ";
    }
    if (count % n_cols_to_print != 0)
      std::cout << "\n";
  }
 
private:
  //============================================================================
  void set_knots(double x0, double xmax, std::size_t n_break, KnotDistro kd) {
    // Make version that takes custom (internal) knot sequence?
    // Option for linear/log-spaced knots?
 
    // n_break break points range from [0,xmax] in n_break steps
    // We define the _internal_ (non-zero) breakpoints
    // these run from [x0, xmax] in n_break - 1 steps
    // Knots are the same as breakpoints, but the first (0) and last (xmax)
    // points are repeated K times
    auto breaks = (kd == KnotDistro::logarithmic) ?
                      qip::logarithmic_range(x0, xmax, n_break - 1) :
                  (kd == KnotDistro::linear) ?
                      qip::uniform_range(x0, xmax, n_break - 1) :
                  (kd == KnotDistro::loglinear) ?
                      qip::loglinear_range(x0, xmax, 0.5 * xmax, n_break - 1) :
                      std::vector<double>{};
 
    breaks.insert(breaks.begin(), 0.0);
 
    auto gsl_break_vec = gsl_vector_view_array(breaks.data(), breaks.size());
 
#ifdef GSL_VERSION_PRIOR_2_8
    // gsl_bspline_knots will be deprecated in future GSL
    gsl_bspline_knots(&gsl_break_vec.vector, gsl_bspl_work);
#else
    gsl_bspline_init_augment(&gsl_break_vec.vector, gsl_bspl_work);
#endif
  }
};
 
// Documentation from OLD GSL v:1.x
//  -- Function: gsl_bspline_deriv_workspace * gsl_bspline_deriv_alloc
//           (const size_t K)
//      This function allocates a workspace for computing the derivatives
//      of a B-spline basis function of order K.  The size of the workspace
//      is O(2k^2).
//
//  -- Function: int gsl_bspline_deriv_eval (const double X, const size_t
//           NDERIV, gsl_matrix * DB, gsl_bspline_workspace * W,
//           gsl_bspline_deriv_workspace * DW)
//      This function evaluates all B-spline basis function derivatives of
//      orders 0 through nderiv (inclusive) at the position X and stores
//      them in the matrix DB.  The (i,j)-th element of DB is
//      d^jB_i(x)/dx^j.  The matrix DB must be of size n = nbreak + k - 2
//      by nderiv + 1.  The value n may also be obtained by calling
//      `gsl_bspline_ncoeffs'.  Note that function evaluations are
//      included as the zeroth order derivatives in DB.  Computing all the
//      basis function derivatives at once is more efficient than
//      computing them individually, due to the nature of the defining
//      recurrence relation.
//
//  -- Function: int gsl_bspline_deriv_eval_nonzero (const double X, const
//           size_t NDERIV, gsl_matrix * DB, size_t * ISTART, size_t *
//           IEND, gsl_bspline_workspace * W, gsl_bspline_deriv_workspace
//           * DW)
//      This function evaluates all potentially nonzero B-spline basis
//      function derivatives of orders 0 through nderiv (inclusive) at the
//      position X and stores them in the matrix DB.  The (i,j)-th element
//      of DB is d^j/dx^j B_(istart+i)(x).  The last row of DB contains
//      d^j/dx^j B_(iend)(x).  The matrix DB must be of size k by at least
//      nderiv + 1.  Note that function evaluations are included as the
//      zeroth order derivatives in DB.  By returning only the nonzero
//      basis functions, this function allows quantities involving linear
//      combinations of the B_i(x) and their derivatives to be computed
//      without unnecessary terms.