interp: add ClampedCubic interpolator
MODIFY
interp/cubic.go
interp/cubic.go
@@ -377,8 +377,9 @@ func makeCubicSplineSecondDerivativeEquations(xs, ys []float64) (*mat.Tridiag, m
// NaturalCubic is a piecewise cubic 1-dimensional interpolator with
// continuous value, first and second derivatives, which can be fitted to (X, Y)
-// value pairs without providing derivatives. See e.g. https://www.math.drexel.edu/~tolya/cubicspline.pdf
-// for details.
+// value pairs without providing derivatives. It uses the boundary conditions
+// Y''(left end ) = Y''(right end) = 0.
+// See e.g. https://www.math.drexel.edu/~tolya/cubicspline.pdf for details.
type NaturalCubic struct {
cubic PiecewiseCubic
}
@@ -412,3 +413,48 @@ func (nc *NaturalCubic) Fit(xs, ys []float64) error {
}
return err
}
+
+// ClampedCubic is a piecewise cubic 1-dimensional interpolator with
+// continuous value, first and second derivatives, which can be fitted to (X, Y)
+// value pairs without providing derivatives. It uses the boundary conditions
+// Y'(left end ) = Y'(right end) = 0.
+type ClampedCubic struct {
+ cubic PiecewiseCubic
+}
+
+// Predict returns the interpolation value at x.
+func (cc *ClampedCubic) Predict(x float64) float64 {
+ return cc.cubic.Predict(x)
+}
+
+// PredictDerivative returns the predicted derivative at x.
+func (cc *ClampedCubic) PredictDerivative(x float64) float64 {
+ return cc.cubic.PredictDerivative(x)
+}
+
+// Fit fits a predictor to (X, Y) value pairs provided as two slices.
+// It panics if len(xs) < 2, elements of xs are not strictly increasing
+// or len(xs) != len(ys). It returns an error if solving the required system
+// of linear equations fails.
+func (cc *ClampedCubic) Fit(xs, ys []float64) error {
+ a, b := makeCubicSplineSecondDerivativeEquations(xs, ys)
+ // Add boundary conditions y''(left) = y''(right) = 0:
+ n := len(xs)
+ // Condition Y'(left end) = 0:
+ dxL := xs[1] - xs[0]
+ b.SetVec(0, (ys[1]-ys[0])/dxL)
+ a.SetBand(0, 0, dxL/3)
+ a.SetBand(0, 1, dxL/6)
+ // Condition Y'(right end) = 0:
+ m := n - 1
+ dxR := xs[m] - xs[m-1]
+ b.SetVec(m, (ys[m]-ys[m-1])/dxR)
+ a.SetBand(m, m, -dxR/3)
+ a.SetBand(m, m-1, -dxR/6)
+ x := mat.NewVecDense(n, nil)
+ err := a.SolveVecTo(x, false, b)
+ if err == nil {
+ cc.cubic.fitWithSecondDerivatives(xs, ys, x.RawVector().Data)
+ }
+ return err
+}

MODIFY
interp/cubic_test.go
interp/cubic_test.go
@@ -799,6 +799,48 @@ func TestNaturalCubicFit(t *testing.T) {
if math.Abs(got-want[i]) > 1e-14 {
t.Errorf("Mismatch in predicted Y value for x = %g: got %v, want %g", testXs[i], got, want[i])
}
+ got = nc.PredictDerivative(testXs[i])
+ want := discrDerivPredict(&nc, xs[0], xs[len(xs)-1], testXs[i], 1e-9)
+ if math.Abs(got-want) > 1e-6 {
+ t.Errorf("Mismatch in predicted dY/dX value for x = %g: got %v, want %g", testXs[i], got, want)
+ }
+ }
}
+func TestClampedCubicFit(t *testing.T) {
+ t.Parallel()
+ xs := []float64{-1, 0, 2, 3.5}
+ ys := []float64{2, 0, 2, 1.5}
+ var cc ClampedCubic
+ err := cc.Fit(xs, ys)
+ if err != nil {
+ t.Errorf("Error when fitting ClampedCubic: %v", err)
+ }
+ testXs := []float64{-1, -0.99, -0.5, 0, 0.5, 1, 1.5, 2, 2.5, 3, 3.49, 3.5}
+ // From scipy.interpolate.CubicSpline:
+ want := []float64{
+ 2.0,
+ 1.999564111111111,
+ 1.2021604938271606,
+ 0.0,
+ -0.20833333333333343,
+ 0.41975308641975284,
+ 1.337962962962962,
+ 2.0,
+ 2.026748971193416,
+ 1.7078189300411522,
+ 1.5001129711934151,
+ 1.4999999999999996,
+ }
+ for i := 0; i < len(testXs); i++ {
+ got := cc.Predict(testXs[i])
+ if math.Abs(got-want[i]) > 1e-14 {
+ t.Errorf("Mismatch in predicted Y value for x = %g: got %v, want %g", testXs[i], got, want[i])
+ }
+ got = cc.PredictDerivative(testXs[i])
+ want := discrDerivPredict(&cc, xs[0], xs[len(xs)-1], testXs[i], 1e-9)
+ if math.Abs(got-want) > 1e-6 {
+ t.Errorf("Mismatch in predicted dY/dX value for x = %g: got %v, want %g", testXs[i], got, want)
+ }
+ }
}

interp: add ClampedCubic interpolator (#1662)