interp: Add Fritsch-Butland piecewise cubic interpolation method
Also adds example of fit output for piecewise constant, piecewise linear, Akima and Fritsch-Butland methods.
MODIFY
interp/interp.go
interp/interp.go
@@ -86,15 +86,7 @@ func (pl *PiecewiseLinear) Fit(xs, ys []float64) error {
if n < 2 {
panic(tooFewPoints)
}
- m := n - 1
- pl.slopes = make([]float64, m)
- for i := 0; i < m; i++ {
- dx := xs[i+1] - xs[i]
- if dx <= 0 {
- panic(xsNotStrictlyIncreasing)
- }
- pl.slopes[i] = (ys[i+1] - ys[i]) / dx
- }
+ pl.slopes = calculateSlopes(xs, ys)
pl.xs = make([]float64, n)
pl.ys = make([]float64, n)
copy(pl.xs, xs)
@@ -365,3 +357,129 @@ func akimaWeights(slopes []float64, i int) (float64, float64) {
wRight := math.Abs(slopes[i+1] - slopes[i])
return wLeft, wRight
}
+
+// FritschButland is a piecewise cubic 1-dimensional interpolator with
+// continuous value and first derivative, which can be fitted to (X, Y)
+// value pairs without providing derivatives.
+// It is monotone, local and produces a C^1 curve. Its downside is that
+// exhibits high tension, flattening out unnaturally the interpolated
+// curve between the nodes.
+// See Fritsch, F. N. and Butland, J., "A method for constructing local
+// monotone piecewise cubic interpolants" (1984), SIAM J. Sci. Statist.
+// Comput., 5(2), pp. 300-304.
+type FritschButland struct {
+ cubic PiecewiseCubic
+}
+
+// Predict returns the interpolation value at x.
+func (fb *FritschButland) Predict(x float64) float64 {
+ return fb.cubic.Predict(x)
+}
+
+// PredictDerivative returns the predicted derivative at x.
+func (fb *FritschButland) PredictDerivative(x float64) float64 {
+ return fb.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). Always returns nil.
+func (fb *FritschButland) Fit(xs, ys []float64) error {
+ n := len(xs)
+ if n < 2 {
+ panic(tooFewPoints)
+ }
+ if len(ys) != n {
+ panic(differentLengths)
+ }
+ dydxs := make([]float64, n)
+
+ if n == 2 {
+ dx := xs[1] - xs[0]
+ slope := (ys[1] - ys[0]) / dx
+ dydxs[0] = slope
+ dydxs[1] = slope
+ fb.cubic.FitWithDerivatives(xs, ys, dydxs)
+ return nil
+ }
+ slopes := calculateSlopes(xs, ys)
+ m := len(slopes)
+ prevSlope := slopes[0]
+ for i := 1; i < m; i++ {
+ slope := slopes[i]
+ if slope*prevSlope > 0 {
+ dydxs[i] = 3 * (xs[i+1] - xs[i-1]) / ((2*xs[i+1]-xs[i-1]-xs[i])/slopes[i-1] +
+ (xs[i+1]+xs[i]-2*xs[i-1])/slopes[i])
+ } else {
+ dydxs[i] = 0
+ }
+ prevSlope = slope
+ }
+ dydxs[0] = fritschButlandEdgeDerivative(xs, ys, slopes, true)
+ dydxs[m] = fritschButlandEdgeDerivative(xs, ys, slopes, false)
+ fb.cubic.FitWithDerivatives(xs, ys, dydxs)
+ return nil
+}
+
+// fritschButlandEdgeDerivative calculates dy/dx approximation for the
+// Fritsch-Butland method for the left or right edge node.
+func fritschButlandEdgeDerivative(xs, ys, slopes []float64, leftEdge bool) float64 {
+ n := len(xs)
+ var dE, dI, h, hE, f float64
+ if leftEdge {
+ dE = slopes[0]
+ dI = slopes[1]
+ xE := xs[0]
+ xM := xs[1]
+ xI := xs[2]
+ hE = xM - xE
+ h = xI - xE
+ f = xM + xI - 2*xE
+ } else {
+ dE = slopes[n-2]
+ dI = slopes[n-3]
+ xE := xs[n-1]
+ xM := xs[n-2]
+ xI := xs[n-3]
+ hE = xE - xM
+ h = xE - xI
+ f = 2*xE - xI - xM
+ }
+ g := (f*dE - hE*dI) / h
+ if g*dE <= 0 {
+ return 0
+ }
+ if dE*dI <= 0 && math.Abs(g) > 3*math.Abs(dE) {
+ return 3 * dE
+ }
+ return g
+}
+
+// calculateSlopes calculates slopes (ys[i+1] - ys[i]) / (xs[i+1] - xs[i]).
+// It panics if len(xs) < 2, elements of xs are not strictly increasing
+// or len(xs) != len(ys).
+func calculateSlopes(xs, ys []float64) []float64 {
+ n := len(xs)
+ if n <= 2 {
+ panic(tooFewPoints)
+ }
+ if len(ys) != n {
+ panic(differentLengths)
+ }
+ m := n - 1
+ slopes := make([]float64, m)
+ prevX := xs[0]
+ prevY := ys[0]
+ for i := 0; i < m; i++ {
+ x := xs[i+1]
+ y := ys[i+1]
+ dx := x - prevX
+ if dx <= 0 {
+ panic(xsNotStrictlyIncreasing)
+ }
+ slopes[i] = (y - prevY) / dx
+ prevX = x
+ prevY = y
+ }
+ return slopes
+}

MODIFY
interp/interp_test.go
interp/interp_test.go
@@ -659,6 +659,234 @@ func TestAkimaWeights(t *testing.T) {
}
}
+func TestFritschButland(t *testing.T) {
+ t.Parallel()
+ const (
+ tol = 1e-14
+ nPts = 100
+ )
+ for k, test := range []struct {
+ xs, ys []float64
+ }{
+ {
+ xs: []float64{0, 2},
+ ys: []float64{0, 0.5},
+ },
+ {
+ xs: []float64{0, 2},
+ ys: []float64{0, -0.5},
+ },
+ {
+ xs: []float64{0, 2},
+ ys: []float64{0, 0},
+ },
+ {
+ xs: []float64{0, 2, 3, 4},
+ ys: []float64{0, 1, 2, 2.5},
+ },
+ {
+ xs: []float64{0, 2, 3, 4},
+ ys: []float64{0, 1.5, 1.5, 2.5},
+ },
+ {
+ xs: []float64{0, 2, 3, 4},
+ ys: []float64{0, 1.5, 1.5, 1},
+ },
+ {
+ xs: []float64{0, 2, 3, 4},
+ ys: []float64{0, 2.5, 1.5, 1},
+ },
+ {
+ xs: []float64{0, 2, 3, 4},
+ ys: []float64{0, 2.5, 1.5, 2},
+ },
+ {
+ xs: []float64{0, 2, 3, 4},
+ ys: []float64{4, 3, 2, 1},
+ },
+ {
+ xs: []float64{0, 2, 3, 4},
+ ys: []float64{4, 3, 2, 2},
+ },
+ {
+ xs: []float64{0, 2, 3, 4},
+ ys: []float64{4, 3, 2, 5},
+ },
+ {
+ xs: []float64{0, 2, 3, 4, 5, 6},
+ ys: []float64{0, 1, 0.5, 0.5, 1.5, 1.5},
+ },
+ {
+ xs: []float64{0, 2, 3, 4, 5, 6},
+ ys: []float64{0, 1, 1.5, 2.5, 1.5, 1},
+ },
+ {
+ xs: []float64{0, 2, 3, 4, 5, 6},
+ ys: []float64{0, -1, -1.5, -2.5, -1.5, -1},
+ },
+ {
+ xs: []float64{0, 2, 3, 4, 5, 6},
+ ys: []float64{0, 1, 0.5, 1.5, 1, 2},
+ },
+ {
+ xs: []float64{0, 2, 3, 4, 5, 6},
+ ys: []float64{0, 1, 1.5, 2.5, 3, 4},
+ },
+ {
+ xs: []float64{0, 2, 3, 4, 5, 6},
+ ys: []float64{0, 0.0001, -1.5, -2.5, -0.0001, 0},
+ },
+ } {
+ var fb FritschButland
+ err := fb.Fit(test.xs, test.ys)
+ if err != nil {
+ t.Errorf("Error when fitting FritschButland in test case %d: %v", k, err)
+ }
+ n := len(test.xs)
+ for i := 0; i < n; i++ {
+ got := fb.Predict(test.xs[i])
+ want := test.ys[i]
+ if got != want {
+ t.Errorf("Mismatch in interpolated value for node %d in test case %d: got %v, want %g", i, k, got, want)
+ }
+ }
+ if n == 2 {
+ h := test.xs[1] - test.xs[0]
+ want := (test.ys[1] - test.ys[0]) / h
+ for i := 0; i < 2; i++ {
+ got := fb.PredictDerivative(test.xs[i])
+ if !scalar.EqualWithinAbs(got, want, tol) {
+ t.Errorf("Mismatch in approximated derivative for node %d in 2-node test case %d: got %v, want %g", i, k, got, want)
+ }
+ }
+ dx := h / (nPts + 1)
+ for i := 1; i < nPts; i++ {
+ x := test.xs[0] + float64(i)*dx
+ got := fb.PredictDerivative(x)
+ if !scalar.EqualWithinAbs(got, want, tol) {
+ t.Errorf("Mismatch in interpolated derivative for x == %g in 2-node test case %d: got %v, want %g", x, k, got, want)
+ }
+ }
+ } else {
+ m := n - 1
+ for i := 1; i < m; i++ {
+ got := fb.PredictDerivative(test.xs[i])
+ slope := (test.ys[i+1] - test.ys[i]) / (test.xs[i+1] - test.xs[i])
+ prevSlope := (test.ys[i] - test.ys[i-1]) / (test.xs[i] - test.xs[i-1])
+ if slope*prevSlope > 0 {
+ if got == 0 {
+ t.Errorf("Approximated derivative is zero for node %d in test case %d: %g", i, k, got)
+ } else if math.Signbit(slope) != math.Signbit(got) {
+ t.Errorf("Approximated derivative has wrong sign for node %d in test case %d: got %g, want %g", i, k, math.Copysign(1, got), math.Copysign(1, slope))
+ }
+ } else {
+ if got != 0 {
+ t.Errorf("Approximated derivative is not zero for node %d in test case %d: %g", i, k, got)
+ }
+ }
+ }
+ for i := 0; i < m; i++ {
+ yL := test.ys[i]
+ yR := test.ys[i+1]
+ xL := test.xs[i]
+ dx := (test.xs[i+1] - xL) / (nPts + 1)
+ if yL == yR {
+ for j := 1; j < nPts; j++ {
+ x := xL + float64(j)*dx
+ got := fb.Predict(x)
+ if got != yL {
+ t.Errorf("Mismatch in interpolated value for x == %g in test case %d: got %v, want %g", x, k, got, yL)
+ }
+ got = fb.PredictDerivative(x)
+ if got != 0 {
+ t.Errorf("Interpolated derivative not zero for x == %g in test case %d: got %v", x, k, got)
+ }
+ }
+ } else {
+ minY := math.Min(yL, yR)
+ maxY := math.Max(yL, yR)
+ for j := 1; j < nPts; j++ {
+ x := xL + float64(j)*dx
+ got := fb.Predict(x)
+ if got < minY || got > maxY {
+ t.Errorf("Interpolated value out of [%g, %g] bounds for x == %g in test case %d: got %v", minY, maxY, x, k, got)
+ }
+ got = fb.PredictDerivative(x)
+ dy := yR - yL
+ if got*dy < 0 {
+ t.Errorf("Interpolated derivative has wrong sign for x == %g in test case %d: want %g, got %g", x, k, math.Copysign(1, dy), math.Copysign(1, got))
+ }
+ }
+ }
+ }
+ }
+ }
+}
+
+func TestFritschButlandErrors(t *testing.T) {
+ t.Parallel()
+ for _, test := range []struct {
+ xs, ys []float64
+ }{
+ {
+ xs: []float64{0},
+ ys: []float64{0},
+ },
+ {
+ xs: []float64{0, 1, 2},
+ ys: []float64{0, 1}},
+ {
+ xs: []float64{0, 0, 1},
+ ys: []float64{0, 0, 0},
+ },
+ {
+ xs: []float64{0, 1, 0},
+ ys: []float64{0, 0, 0},
+ },
+ } {
+ var fb FritschButland
+ if !panics(func() { _ = fb.Fit(test.xs, test.ys) }) {
+ t.Errorf("expected panic for xs: %v and ys: %v", test.xs, test.ys)
+ }
+ }
+}
+
+func TestCalculateSlopesErrors(t *testing.T) {
+ t.Parallel()
+ for _, test := range []struct {
+ xs, ys []float64
+ }{
+ {
+ xs: []float64{0},
+ ys: []float64{0},
+ },
+ {
+ xs: []float64{0, 1, 2},
+ ys: []float64{0, 1}},
+ {
+ xs: []float64{0, 0, 1},
+ ys: []float64{0, 0, 0},
+ },
+ {
+ xs: []float64{0, 1, 0},
+ ys: []float64{0, 0, 0},
+ },
+ } {
+ if !panics(func() { calculateSlopes(test.xs, test.ys) }) {
+ t.Errorf("expected panic for xs: %v and ys: %v", test.xs, test.ys)
+ }
+ }
+}
+
+func TestCalculateSlopes(t *testing.T) {
+ t.Parallel()
+ got := calculateSlopes([]float64{0, 2, 3, 5}, []float64{0, 1, 1, -1})
+ want := []float64{0.5, 0, -1}
+ if !floats.EqualApprox(got, want, 1e-14) {
+ t.Errorf("Mismatch in calculated slops: got %v, want %v", got, want)
+ }
+}
+
func applyFunc(xs []float64, f func(x float64) float64) []float64 {
ys := make([]float64, len(xs))
for i, x := range xs {

interp: Add Fritsch-Butland piecewise cubic interpolation method 

Also adds example of fit output for piecewise constant, piecewise linear, Akima and Fritsch-Butland methods.