ceph/src/boost/boost/math/interpolators/detail/cardinal_cubic_b_spline_detail.hpp

   1 // Copyright Nick Thompson, 2017
   2 // Use, modification and distribution are subject to the
   3 // Boost Software License, Version 1.0.
   4 // (See accompanying file LICENSE_1_0.txt
   5 // or copy at http://www.boost.org/LICENSE_1_0.txt)
   6
   7 #ifndef BOOST_MATH_INTERPOLATORS_CARDINAL_CUBIC_B_SPLINE_DETAIL_HPP
   8 #define BOOST_MATH_INTERPOLATORS_CARDINAL_CUBIC_B_SPLINE_DETAIL_HPP
   9
  10 #include <limits>
  11 #include <cmath>
  12 #include <vector>
  13 #include <memory>
  14 #include <boost/math/constants/constants.hpp>
  15 #include <boost/math/special_functions/fpclassify.hpp>
  16 #include <boost/math/special_functions/trunc.hpp>
  17
  18 namespace boost{ namespace math{ namespace interpolators{ namespace detail{
  19
  20
  21 template <class Real>
  22 class cardinal_cubic_b_spline_imp
  23 {
  24 public:
  25     // If you don't know the value of the derivative at the endpoints, leave them as nans and the routine will estimate them.
  26     // f[0] = f(a), f[length -1] = b, step_size = (b - a)/(length -1).
  27     template <class BidiIterator>
  28     cardinal_cubic_b_spline_imp(BidiIterator f, BidiIterator end_p, Real left_endpoint, Real step_size,
  29                        Real left_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN(),
  30                        Real right_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN());
  31
  32     Real operator()(Real x) const;
  33
  34     Real prime(Real x) const;
  35
  36     Real double_prime(Real x) const;
  37
  38 private:
  39     std::vector<Real> m_beta;
  40     Real m_h_inv;
  41     Real m_a;
  42     Real m_avg;
  43 };
  44
  45
  46
  47 template <class Real>
  48 Real b3_spline(Real x)
  49 {
  50     using std::abs;
  51     Real absx = abs(x);
  52     if (absx < 1)
  53     {
  54         Real y = 2 - absx;
  55         Real z = 1 - absx;
  56         return boost::math::constants::sixth<Real>()*(y*y*y - 4*z*z*z);
  57     }
  58     if (absx < 2)
  59     {
  60         Real y = 2 - absx;
  61         return boost::math::constants::sixth<Real>()*y*y*y;
  62     }
  63     return static_cast<Real>(0);
  64 }
  65
  66 template<class Real>
  67 Real b3_spline_prime(Real x)
  68 {
  69     if (x < 0)
  70     {
  71         return -b3_spline_prime(-x);
  72     }
  73
  74     if (x < 1)
  75     {
  76         return x*(3*boost::math::constants::half<Real>()*x - 2);
  77     }
  78     if (x < 2)
  79     {
  80         return -boost::math::constants::half<Real>()*(2 - x)*(2 - x);
  81     }
  82     return static_cast<Real>(0);
  83 }
  84
  85 template<class Real>
  86 Real b3_spline_double_prime(Real x)
  87 {
  88     if (x < 0)
  89     {
  90         return b3_spline_double_prime(-x);
  91     }
  92
  93     if (x < 1)
  94     {
  95         return 3*x - 2;
  96     }
  97     if (x < 2)
  98     {
  99         return (2 - x);
 100     }
 101     return static_cast<Real>(0);
 102 }
 103
 104
 105 template <class Real>
 106 template <class BidiIterator>
 107 cardinal_cubic_b_spline_imp<Real>::cardinal_cubic_b_spline_imp(BidiIterator f, BidiIterator end_p, Real left_endpoint, Real step_size,
 108                                              Real left_endpoint_derivative, Real right_endpoint_derivative) : m_a(left_endpoint), m_avg(0)
 109 {
 110     using boost::math::constants::third;
 111
 112     std::size_t length = end_p - f;
 113
 114     if (length < 5)
 115     {
 116         if (boost::math::isnan(left_endpoint_derivative) || boost::math::isnan(right_endpoint_derivative))
 117         {
 118             throw std::logic_error("Interpolation using a cubic b spline with derivatives estimated at the endpoints requires at least 5 points.\n");
 119         }
 120         if (length < 3)
 121         {
 122             throw std::logic_error("Interpolation using a cubic b spline requires at least 3 points.\n");
 123         }
 124     }
 125
 126     if (boost::math::isnan(left_endpoint))
 127     {
 128         throw std::logic_error("Left endpoint is NAN; this is disallowed.\n");
 129     }
 130     if (left_endpoint + length*step_size >= (std::numeric_limits<Real>::max)())
 131     {
 132         throw std::logic_error("Right endpoint overflows the maximum representable number of the specified precision.\n");
 133     }
 134     if (step_size <= 0)
 135     {
 136         throw std::logic_error("The step size must be strictly > 0.\n");
 137     }
 138     // Storing the inverse of the stepsize does provide a measurable speedup.
 139     // It's not huge, but nonetheless worthwhile.
 140     m_h_inv = 1/step_size;
 141
 142     // Following Kress's notation, s'(a) = a1, s'(b) = b1
 143     Real a1 = left_endpoint_derivative;
 144     // See the finite-difference table on Wikipedia for reference on how
 145     // to construct high-order estimates for one-sided derivatives:
 146     // https://en.wikipedia.org/wiki/Finite_difference_coefficient#Forward_and_backward_finite_difference
 147     // Here, we estimate then to O(h^4), as that is the maximum accuracy we could obtain from this method.
 148     if (boost::math::isnan(a1))
 149     {
 150         // For simple functions (linear, quadratic, so on)
 151         // almost all the error comes from derivative estimation.
 152         // This does pairwise summation which gives us another digit of accuracy over naive summation.
 153         Real t0 = 4*(f[1] + third<Real>()*f[3]);
 154         Real t1 = -(25*third<Real>()*f[0] + f[4])/4  - 3*f[2];
 155         a1 = m_h_inv*(t0 + t1);
 156     }
 157
 158     Real b1 = right_endpoint_derivative;
 159     if (boost::math::isnan(b1))
 160     {
 161         size_t n = length - 1;
 162         Real t0 = 4*(f[n-3] + third<Real>()*f[n - 1]);
 163         Real t1 = -(25*third<Real>()*f[n - 4] + f[n])/4  - 3*f[n - 2];
 164
 165         b1 = m_h_inv*(t0 + t1);
 166     }
 167
 168     // s(x) = \sum \alpha_i B_{3}( (x- x_i - a)/h )
 169     // Of course we must reindex from Kress's notation, since he uses negative indices which make C++ unhappy.
 170     m_beta.resize(length + 2, std::numeric_limits<Real>::quiet_NaN());
 171
 172     // Since the splines have compact support, they decay to zero very fast outside the endpoints.
 173     // This is often very annoying; we'd like to evaluate the interpolant a little bit outside the
 174     // boundary [a,b] without massive error.
 175     // A simple way to deal with this is just to subtract the DC component off the signal, so we need the average.
 176     // This algorithm for computing the average is recommended in
 177     // http://www.heikohoffmann.de/htmlthesis/node134.html
 178     Real t = 1;
 179     for (size_t i = 0; i < length; ++i)
 180     {
 181         if (boost::math::isnan(f[i]))
 182         {
 183             std::string err = "This function you are trying to interpolate is a nan at index " + std::to_string(i) + "\n";
 184             throw std::logic_error(err);
 185         }
 186         m_avg += (f[i] - m_avg) / t;
 187         t += 1;
 188     }
 189
 190
 191     // Now we must solve an almost-tridiagonal system, which requires O(N) operations.
 192     // There are, in fact 5 diagonals, but they only differ from zero on the first and last row,
 193     // so we can patch up the tridiagonal row reduction algorithm to deal with two special rows.
 194     // See Kress, equations 8.41
 195     // The the "tridiagonal" matrix is:
 196     // 1  0 -1
 197     // 1  4  1
 198     //    1  4  1
 199     //       1  4  1
 200     //          ....
 201     //          1  4  1
 202     //          1  0 -1
 203     // Numerical estimate indicate that as N->Infinity, cond(A) -> 6.9, so this matrix is good.
 204     std::vector<Real> rhs(length + 2, std::numeric_limits<Real>::quiet_NaN());
 205     std::vector<Real> super_diagonal(length + 2, std::numeric_limits<Real>::quiet_NaN());
 206
 207     rhs[0] = -2*step_size*a1;
 208     rhs[rhs.size() - 1] = -2*step_size*b1;
 209
 210     super_diagonal[0] = 0;
 211
 212     for(size_t i = 1; i < rhs.size() - 1; ++i)
 213     {
 214         rhs[i] = 6*(f[i - 1] - m_avg);
 215         super_diagonal[i] = 1;
 216     }
 217
 218
 219     // One step of row reduction on the first row to patch up the 5-diagonal problem:
 220     // 1 0 -1 | r0
 221     // 1 4 1  | r1
 222     // mapsto:
 223     // 1 0 -1 | r0
 224     // 0 4 2  | r1 - r0
 225     // mapsto
 226     // 1 0 -1 | r0
 227     // 0 1 1/2| (r1 - r0)/4
 228     super_diagonal[1] = 0.5;
 229     rhs[1] = (rhs[1] - rhs[0])/4;
 230
 231     // Now do a tridiagonal row reduction the standard way, until just before the last row:
 232     for (size_t i = 2; i < rhs.size() - 1; ++i)
 233     {
 234         Real diagonal = 4 - super_diagonal[i - 1];
 235         rhs[i] = (rhs[i] - rhs[i - 1])/diagonal;
 236         super_diagonal[i] /= diagonal;
 237     }
 238
 239     // Now the last row, which is in the form
 240     // 1 sd[n-3] 0      | rhs[n-3]
 241     // 0  1     sd[n-2] | rhs[n-2]
 242     // 1  0     -1      | rhs[n-1]
 243     Real final_subdiag = -super_diagonal[rhs.size() - 3];
 244     rhs[rhs.size() - 1] = (rhs[rhs.size() - 1] - rhs[rhs.size() - 3])/final_subdiag;
 245     Real final_diag = -1/final_subdiag;
 246     // Now we're here:
 247     // 1 sd[n-3] 0         | rhs[n-3]
 248     // 0  1     sd[n-2]    | rhs[n-2]
 249     // 0  1     final_diag | (rhs[n-1] - rhs[n-3])/diag
 250
 251     final_diag = final_diag - super_diagonal[rhs.size() - 2];
 252     rhs[rhs.size() - 1] = rhs[rhs.size() - 1] - rhs[rhs.size() - 2];
 253
 254
 255     // Back substitutions:
 256     m_beta[rhs.size() - 1] = rhs[rhs.size() - 1]/final_diag;
 257     for(size_t i = rhs.size() - 2; i > 0; --i)
 258     {
 259         m_beta[i] = rhs[i] - super_diagonal[i]*m_beta[i + 1];
 260     }
 261     m_beta[0] = m_beta[2] + rhs[0];
 262 }
 263
 264 template<class Real>
 265 Real cardinal_cubic_b_spline_imp<Real>::operator()(Real x) const
 266 {
 267     // See Kress, 8.40: Since B3 has compact support, we don't have to sum over all terms,
 268     // just the (at most 5) whose support overlaps the argument.
 269     Real z = m_avg;
 270     Real t = m_h_inv*(x - m_a) + 1;
 271
 272     using std::max;
 273     using std::min;
 274     using std::ceil;
 275     using std::floor;
 276
 277     size_t k_min = static_cast<size_t>((max)(static_cast<long>(0), boost::math::ltrunc(ceil(t - 2))));
 278     size_t k_max = static_cast<size_t>((max)((min)(static_cast<long>(m_beta.size() - 1), boost::math::ltrunc(floor(t + 2))), 0l));
 279
 280     for (size_t k = k_min; k <= k_max; ++k)
 281     {
 282         z += m_beta[k]*b3_spline(t - k);
 283     }
 284
 285     return z;
 286 }
 287
 288 template<class Real>
 289 Real cardinal_cubic_b_spline_imp<Real>::prime(Real x) const
 290 {
 291     Real z = 0;
 292     Real t = m_h_inv*(x - m_a) + 1;
 293
 294     using std::max;
 295     using std::min;
 296     using std::ceil;
 297     using std::floor;
 298
 299     size_t k_min = static_cast<size_t>((max)(static_cast<long>(0), boost::math::ltrunc(ceil(t - 2))));
 300     size_t k_max = static_cast<size_t>((min)(static_cast<long>(m_beta.size() - 1), boost::math::ltrunc(floor(t + 2))));
 301
 302     for (size_t k = k_min; k <= k_max; ++k)
 303     {
 304         z += m_beta[k]*b3_spline_prime(t - k);
 305     }
 306     return z*m_h_inv;
 307 }
 308
 309 template<class Real>
 310 Real cardinal_cubic_b_spline_imp<Real>::double_prime(Real x) const
 311 {
 312     Real z = 0;
 313     Real t = m_h_inv*(x - m_a) + 1;
 314
 315     using std::max;
 316     using std::min;
 317     using std::ceil;
 318     using std::floor;
 319
 320     size_t k_min = static_cast<size_t>((max)(static_cast<long>(0), boost::math::ltrunc(ceil(t - 2))));
 321     size_t k_max = static_cast<size_t>((min)(static_cast<long>(m_beta.size() - 1), boost::math::ltrunc(floor(t + 2))));
 322
 323     for (size_t k = k_min; k <= k_max; ++k)
 324     {
 325         z += m_beta[k]*b3_spline_double_prime(t - k);
 326     }
 327     return z*m_h_inv*m_h_inv;
 328 }
 329
 330 }}}}
 331 #endif