Bug in CubicSplineInterpolator integration function

[jamiescottie1]

The way the integral in the CubicSplineInterpolator is implemented is slightly off.
We need to evaluate the integral
$$\int_0^t (1-t')y_l + t'y_r + t'(1-t')((1-t')a + t'b)\mathrm{d}t'$$
which results in (see WolframAlpha where $l\to y_l$, $r\to y_r$, $k\to t$, $t\to t'$)
$$\int_0^t\dots\mathrm{d}t' = \frac{t}{12}\left(at\left(3t^2-8t+6\right) + bt^2\left(4-3t\right) - 6ty_l + 6ty_r + 12y_l\right)$$

If we rewrite this as implemented in the current code, we get
$$\int_0^t\dots\mathrm{d}t' = a\left(\frac{1}{4}t^4 - \frac{2}{3}t^3 + \frac{1}{2}t^2\right) + b\left(\frac{1}{3}t^3 - \frac{1}{4}t^4\right) + y_l\left(t - \frac{1}{2}t^2\right) + \frac{1}{2}y_rt$$

This is almost the same as what is implemented, with the only difference being the factor $\frac{1}{2}$ in front of the quadratic term $t^2$ that gets multiplied by $a$, which is missing in the current implementation.

---

Side note: For performance reasons the various integer powers should be replaced by explicit products, which are much more efficient than the general purpose pow.

[CD3]

Thanks for bug report. I'm sure your right, but there must have been an example that caused you to notice. Can you give a similar example that we can use as a unit test?

[jamiescottie1]

Hi!

Sure, I noticed with a unit-test using an extension I implemented for my team. I extended the CubicSplineInterpolator so that the user can specify the boundary conditions, which allowed me to compare with the analytic exact integral of a third-order polynomial ($y(x) = x^3 - 2x^2 - 8x + 4$).

Since the original code is restricted to natural boundary conditions ($\mathrm{d}^2y/dx^2=0$) the cubic polynomial is not interpolated exactly, so I fell back to using scipy results as reference (using their CubicSpline implementation with natural BCs),
This is the python file used for generating the reference results:

import numpy as np
import scipy as sp

def fCub(x):
    return x**3 - 2*x**2 - 8*x + 4
    

xSetup = np.array([-2.        , -1.91687501, -1.79924508, -1.63602439, -1.57188369,
                   -1.45164833, -0.76586246, -0.68397459,  0.0585834 ,  0.19441298,
                    0.75255509,  1.18348665,  1.71673075,  1.7273539 ,  2.95432601,
                    3.87390714,  3.91713465,  4.00759876,  4.34247681,  4.51385375,
                    4.73297478,  5.])

ySetup = fCub(xSetup)
interp = sp.interpolate.CubicSpline(xSetup, ySetup, bc_type='natural')
xTest = [(0.0, 1.2), (2.1, 2.4), (1.3, 3.0)]

antiderivs = []

for x in xTest:
    antiderivs.append((interp.antiderivative()(x[0]), interp.antiderivative()(x[1])))

for d in antiderivs:
    print(f'{d[0]:.12f}, {d[1]:.12f}')

and the C++ code (an excerpt from my unit tests):

#include <boost/math/special_functions.hpp>
#include "doctest.h"

#include "libInterpolate/Interpolators/_1D/CubicSplineInterpolator.hpp"

#define boostCub(x) boost::math::pow<3>(x)
#define boostSq(x) boost::math::pow<2>(x)
#define relD(x, y) boost::math::relative_difference(x, y)

 // Common coordinates for setup
static std::array<double, 22> xSetup{ -2.        , -1.91687501, -1.79924508, -1.63602439, -1.57188369,
                                       -1.45164833, -0.76586246, -0.68397459,  0.0585834 ,  0.19441298,
                                        0.75255509,  1.18348665,  1.71673075,  1.7273539 ,  2.95432601,
                                        3.87390714,  3.91713465,  4.00759876,  4.34247681,  4.51385375,
                                        4.73297478,  5. };
static std::array<std::pair<double, double>, 3> xIntPairsTest{ std::make_pair(0.0, 1.2), std::make_pair(2.1, 2.4), std::make_pair(1.3, 3.0) };

// Third order polynomial test function
static double fCub(const double x)
{
   return boostCub(x) - 2. * boostSq(x) - 8. * x + 4.;
}

// Helper function
template<class Interpolator>
static void checkIntegrals(const Interpolator& interp, const std::array<std::pair<double, double>, 3>& arr, const double relDiff = 1E-12)
{
   const auto& p = xIntPairsTest;
   const std::array<double, 3> expected{
       (arr[0].second - arr[0].first),
       (arr[1].second - arr[1].first),
       (arr[2].second - arr[2].first)
   };

   for (int iPair = 0; iPair < expected.size(); ++iPair)
      REQUIRE(relD(interp.integral(p[iPair].first, p[iPair].second), expected[iPair]) < relDiff);
}

// The actual test
TEST_CASE("InterpolationTests_1D")
{
   constexpr double nan = std::numeric_limits<double>::quiet_NaN();

   // Compute y and dydx values for interpolator setup
   std::array<double, 22> ySetupCub, dydxSetupCub;
   std::transform(xSetup.cbegin(), xSetup.cend(), ySetupCub.begin(), fCub);

   SUBCASE("Cubic Spline (Natural) interpolator")
   {
      SUBCASE("Cubic interpolant")
      {
         _1D::CubicSplineInterpolator<double> interp(xSetup, ySetupCub);

         std::array<std::pair<double, double>, 3> aDsScipy{ std::make_pair(14.666496160949, 13.072896191366), std::make_pair(4.114519860572, 0.304890322266), std::make_pair(12.355854524369, -7.083518709314) };
         checkIntegrals(interp, aDsScipy);
      }
   }
}

This should fail with the original code and complete successfully with the fix (tested on my machine compiled with MSVC v143).