interp: Add NaturalCubic spline interpolator
MODIFY
interp/cubic.go
interp/cubic.go
@@ -104,8 +104,7 @@ func (pc *PiecewiseCubic) FitWithDerivatives(xs, ys, dydxs []float64) {
pc.coeffs.Set(i, 2, (3*dy-(2*dydxs[i]+dydxs[i+1])*dx)/dx/dx)
pc.coeffs.Set(i, 3, (-2*dy+(dydxs[i]+dydxs[i+1])*dx)/dx/dx/dx)
}
- pc.xs = make([]float64, n)
- copy(pc.xs, xs)
+ pc.xs = append(pc.xs[:0], xs...)
pc.lastY = ys[m]
pc.lastDyDx = dydxs[m]
}
@@ -297,3 +296,119 @@ func fritschButlandEdgeDerivative(xs, ys, slopes []float64, leftEdge bool) float
}
return g
}
+
+// fitWithSecondDerivatives fits a piecewise cubic predictor to (X, Y, d^2Y/dX^2) value
+// triples provided as three slices.
+// It panics if any of these is true:
+// - len(xs) < 2,
+// - elements of xs are not strictly increasing,
+// - len(xs) != len(ys),
+// - len(xs) != len(d2ydx2s).
+// Note that this method does not guarantee on its own the continuity of first derivatives.
+func (pc *PiecewiseCubic) fitWithSecondDerivatives(xs, ys, d2ydx2s []float64) {
+ n := len(xs)
+ switch {
+ case len(ys) != n, len(d2ydx2s) != n:
+ panic(differentLengths)
+ case n < 2:
+ panic(tooFewPoints)
+ }
+ m := n - 1
+ pc.coeffs.Reset()
+ pc.coeffs.ReuseAs(m, 4)
+ for i := 0; i < m; i++ {
+ dx := xs[i+1] - xs[i]
+ if dx <= 0 {
+ panic(xsNotStrictlyIncreasing)
+ }
+ dy := ys[i+1] - ys[i]
+ dm := d2ydx2s[i+1] - d2ydx2s[i]
+ pc.coeffs.Set(i, 0, ys[i]) // a_0
+ pc.coeffs.Set(i, 1, (dy-(d2ydx2s[i]+dm/3)*dx*dx/2)/dx) // a_1
+ pc.coeffs.Set(i, 2, d2ydx2s[i]/2) // a_2
+ pc.coeffs.Set(i, 3, dm/6/dx) // a_3
+ }
+ pc.xs = append(pc.xs[:0], xs...)
+ pc.lastY = ys[m]
+ lastDx := xs[m] - xs[m-1]
+ pc.lastDyDx = pc.coeffs.At(m-1, 1) + 2*pc.coeffs.At(m-1, 2)*lastDx + 3*pc.coeffs.At(m-1, 3)*lastDx*lastDx
+}
+
+// makeCubicSplineSecondDerivativeEquations generates the basic system of linear equations
+// which have to be satisfied by the second derivatives to make the first derivatives of a
+// cubic spline continuous. It panics if elements of xs are not strictly increasing, or
+// len(xs) != len(ys).
+// makeCubicSplineSecondDerivativeEquations returns a tri-diagonal matrix A and a vector b
+// defining a system of linear equations A*m = b for second derivatives vector m.
+func makeCubicSplineSecondDerivativeEquations(xs, ys []float64) (*mat.Tridiag, mat.MutableVector) {
+ n := len(xs)
+ if len(ys) != n {
+ panic(differentLengths)
+ }
+ b := make([]float64, n)
+ m := n - 1
+ // Diagonal of A:
+ d := make([]float64, n)
+ // Sub-diagonal of A:
+ dl := make([]float64, m)
+ // Super-diagonal of A:
+ du := make([]float64, m)
+ if n > 2 {
+ for i := 0; i < m; i++ {
+ dx := xs[i+1] - xs[i]
+ if dx <= 0 {
+ panic(xsNotStrictlyIncreasing)
+ }
+ slope := (ys[i+1] - ys[i]) / dx
+ if i > 0 {
+ b[i] += slope
+ d[i] += dx / 3
+ du[i] = dx / 6
+ }
+ if i < m-1 {
+ b[i+1] -= slope
+ d[i+1] += dx / 3
+ dl[i] = dx / 6
+ }
+ }
+ }
+ return mat.NewTridiag(n, dl, d, du), mat.NewVecDense(n, b)
+}
+
+// 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.
+type NaturalCubic struct {
+ cubic PiecewiseCubic
+}
+
+// Predict returns the interpolation value at x.
+func (nc *NaturalCubic) Predict(x float64) float64 {
+ return nc.cubic.Predict(x)
+}
+
+// PredictDerivative returns the predicted derivative at x.
+func (nc *NaturalCubic) PredictDerivative(x float64) float64 {
+ return nc.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 (nc *NaturalCubic) Fit(xs, ys []float64) error {
+ a, b := makeCubicSplineSecondDerivativeEquations(xs, ys)
+ // Add boundary conditions y''(left) = y''(right) = 0:
+ n := len(xs)
+ b.SetVec(0, 0)
+ b.SetVec(n-1, 0)
+ a.SetBand(0, 0, 1)
+ a.SetBand(n-1, n-1, 1)
+ x := mat.NewVecDense(n, nil)
+ err := a.SolveVecTo(x, false, b)
+ if err == nil {
+ nc.cubic.fitWithSecondDerivatives(xs, ys, x.RawVector().Data)
+ }
+ return err
+}

MODIFY
interp/cubic_test.go
interp/cubic_test.go
@@ -671,3 +671,134 @@ func TestFritschButlandErrors(t *testing.T) {
}
}
}
+
+func TestPiecewiseCubicFitWithSecondDerivatives(t *testing.T) {
+ t.Parallel()
+ const tol = 1e-14
+ xs := []float64{-2, 0, 3}
+ ys := []float64{2.5, 1, 2.5}
+ d2ydx2s := []float64{1, 2, 3}
+ var pc PiecewiseCubic
+ pc.fitWithSecondDerivatives(xs, ys, d2ydx2s)
+ m := len(xs) - 1
+ if pc.lastY != ys[m] {
+ t.Errorf("Mismatch in lastY: got %v, want %g", pc.lastY, ys[m])
+ }
+ if !floats.Equal(pc.xs, xs) {
+ t.Errorf("Mismatch in xs: got %v, want %v", pc.xs, xs)
+ }
+ for i := 0; i < len(xs); i++ {
+ yHat := pc.Predict(xs[i])
+ if math.Abs(yHat-ys[i]) > tol {
+ t.Errorf("Mismatch in predicted Y[%d]: got %v, want %g", i, yHat, ys[i])
+ }
+ var d2ydx2Hat float64
+ if i < m {
+ d2ydx2Hat = 2 * pc.coeffs.At(i, 2)
+ } else {
+ d2ydx2Hat = 2*pc.coeffs.At(m-1, 2) + 6*pc.coeffs.At(m-1, 3)*(xs[m]-xs[m-1])
+ }
+ if math.Abs(d2ydx2Hat-d2ydx2s[i]) > tol {
+ t.Errorf("Mismatch in predicted d2Y/dX2[%d]: got %v, want %g", i, d2ydx2Hat, d2ydx2s[i])
+ }
+ }
+ // Test pc.lastDyDx without copying verbatim the calculation from the tested method:
+ lastDyDx := pc.PredictDerivative(xs[m] - tol/1000)
+ if math.Abs(lastDyDx-pc.lastDyDx) > tol {
+ t.Errorf("Mismatch in lastDxDy: got %v, want %g", pc.lastDyDx, lastDyDx)
+ }
+}
+
+func TestPiecewiseCubicFitWithSecondDerivativesErrors(t *testing.T) {
+ t.Parallel()
+ for _, test := range []struct {
+ xs, ys, d2ydx2s []float64
+ }{
+ {
+ xs: []float64{0, 1, 2},
+ ys: []float64{10, 20},
+ d2ydx2s: []float64{0, 0, 0},
+ },
+ {
+ xs: []float64{0, 1, 1},
+ ys: []float64{10, 20, 30},
+ d2ydx2s: []float64{0, 0, 0, 0},
+ },
+ {
+ xs: []float64{0},
+ ys: []float64{0},
+ d2ydx2s: []float64{0},
+ },
+ {
+ xs: []float64{0, 1, 1},
+ ys: []float64{10, 20, 10},
+ d2ydx2s: []float64{0, 0, 0},
+ },
+ } {
+ var pc PiecewiseCubic
+ if !panics(func() { pc.fitWithSecondDerivatives(test.xs, test.ys, test.d2ydx2s) }) {
+ t.Errorf("expected panic for xs: %v, ys: %v and d2ydx2s: %v", test.xs, test.ys, test.d2ydx2s)
+ }
+ }
+}
+
+func TestMakeCubicSplineSecondDerivativeEquations(t *testing.T) {
+ t.Parallel()
+ const tol = 1e-15
+ xs := []float64{-1, 0, 2}
+ ys := []float64{2, 0, 2}
+ n := len(xs)
+ A, b := makeCubicSplineSecondDerivativeEquations(xs, ys)
+ if b.Len() != n {
+ t.Errorf("Mismatch in b size: got %v, want %d", b.Len(), n)
+ }
+ r, c := A.Dims()
+ if r != n || c != n {
+ t.Errorf("Mismatch in A size: got %d x %d, want %d x %d", r, c, n, n)
+ }
+ expectedB := mat.NewVecDense(3, []float64{0, 3, 0})
+ var diffB mat.VecDense
+ diffB.SubVec(b, expectedB)
+ if diffB.Norm(math.Inf(1)) > tol {
+ t.Errorf("Mismatch in b values: got %v, want %v", b, expectedB)
+ }
+ expectedA := mat.NewDense(3, 3, []float64{0, 0, 0, 1 / 6., 1, 2 / 6., 0, 0, 0})
+ var diffA mat.Dense
+ diffA.Sub(A, expectedA)
+ if diffA.Norm(math.Inf(1)) > tol {
+ t.Errorf("Mismatch in A values: got %v, want %v", A, expectedA)
+ }
+}
+
+func TestNaturalCubicFit(t *testing.T) {
+ t.Parallel()
+ xs := []float64{-1, 0, 2, 3.5}
+ ys := []float64{2, 0, 2, 1.5}
+ var nc NaturalCubic
+ err := nc.Fit(xs, ys)
+ if err != nil {
+ t.Errorf("Error when fitting NaturalCubic: %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.9737725526315788,
+ 0.7664473684210527,
+ 0.0,
+ 0.027960526315789477,
+ 0.6184210526315789,
+ 1.3996710526315788,
+ 2.0,
+ 2.1403508771929824,
+ 1.9122807017543857,
+ 1.508859403508772,
+ 1.5}
+ for i := 0; i < len(testXs); i++ {
+ got := nc.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])
+ }
+ }
+
+}

MODIFY
mat/band.go
mat/band.go
@@ -51,6 +51,9 @@ type RawBander interface {
// A MutableBanded can set elements of a band matrix.
type MutableBanded interface {
Banded
+
+ // SetBand sets the element at row i, column j to the value v.
+ // It panics if the location is outside the appropriate region of the matrix.
SetBand(i, j int, v float64)
}

interp: Add NaturalCubic spline interpolator (#1657)