mirror of
https://github.com/gonum/gonum.git
synced 2025-11-02 23:14:01 +08:00
interp: Move out cubic spline code to a separate file
Co-authored-by: Roman Werpachowski <romanw@google.com>
This commit is contained in:
Roman Werpachowski
committed by
GitHub
GitHub
parent
973d5f0ae3
commit
caa2b8c196
299
interp/cubic.go
Normal file
299
interp/cubic.go
Normal file
@@ -0,0 +1,299 @@
|
||||
// Copyright ©2020 The Gonum Authors. All rights reserved.
|
||||
// Use of this source code is governed by a BSD-style
|
||||
// license that can be found in the LICENSE file.
|
||||
|
||||
package interp
|
||||
|
||||
import (
|
||||
"math"
|
||||
|
||||
"gonum.org/v1/gonum/mat"
|
||||
)
|
||||
|
||||
// PiecewiseCubic is a piecewise cubic 1-dimensional interpolator with
|
||||
// continuous value and first derivative.
|
||||
type PiecewiseCubic struct {
|
||||
// Interpolated X values.
|
||||
xs []float64
|
||||
|
||||
// Coefficients of interpolating cubic polynomials, with
|
||||
// len(xs) - 1 rows and 4 columns. The interpolated value
|
||||
// for xs[i] <= x < xs[i + 1] is defined as
|
||||
// sum_{k = 0}^3 coeffs.At(i, k) * (x - xs[i])^k
|
||||
// To guarantee left-continuity, coeffs.At(i, 0) == ys[i].
|
||||
coeffs mat.Dense
|
||||
|
||||
// Last interpolated Y value, corresponding to xs[len(xs) - 1].
|
||||
lastY float64
|
||||
|
||||
// Last interpolated dY/dX value, corresponding to xs[len(xs) - 1].
|
||||
lastDyDx float64
|
||||
}
|
||||
|
||||
// Predict returns the interpolation value at x.
|
||||
func (pc *PiecewiseCubic) Predict(x float64) float64 {
|
||||
i := findSegment(pc.xs, x)
|
||||
if i < 0 {
|
||||
return pc.coeffs.At(0, 0)
|
||||
}
|
||||
m := len(pc.xs) - 1
|
||||
if x == pc.xs[i] {
|
||||
if i < m {
|
||||
return pc.coeffs.At(i, 0)
|
||||
}
|
||||
return pc.lastY
|
||||
}
|
||||
if i == m {
|
||||
return pc.lastY
|
||||
}
|
||||
dx := x - pc.xs[i]
|
||||
a := pc.coeffs.RawRowView(i)
|
||||
return ((a[3]*dx+a[2])*dx+a[1])*dx + a[0]
|
||||
}
|
||||
|
||||
// PredictDerivative returns the predicted derivative at x.
|
||||
func (pc *PiecewiseCubic) PredictDerivative(x float64) float64 {
|
||||
i := findSegment(pc.xs, x)
|
||||
if i < 0 {
|
||||
return pc.coeffs.At(0, 1)
|
||||
}
|
||||
m := len(pc.xs) - 1
|
||||
if x == pc.xs[i] {
|
||||
if i < m {
|
||||
return pc.coeffs.At(i, 1)
|
||||
}
|
||||
return pc.lastDyDx
|
||||
}
|
||||
if i == m {
|
||||
return pc.lastDyDx
|
||||
}
|
||||
dx := x - pc.xs[i]
|
||||
a := pc.coeffs.RawRowView(i)
|
||||
return (3*a[3]*dx+2*a[2])*dx + a[1]
|
||||
}
|
||||
|
||||
// FitWithDerivatives fits a piecewise cubic predictor to (X, Y, dY/dX) value
|
||||
// triples provided as three slices.
|
||||
// It panics if len(xs) < 2, elements of xs are not strictly increasing,
|
||||
// len(xs) != len(ys) or len(xs) != len(dydxs).
|
||||
func (pc *PiecewiseCubic) FitWithDerivatives(xs, ys, dydxs []float64) {
|
||||
n := len(xs)
|
||||
if len(ys) != n {
|
||||
panic(differentLengths)
|
||||
}
|
||||
if len(dydxs) != n {
|
||||
panic(differentLengths)
|
||||
}
|
||||
if 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]
|
||||
// a_0
|
||||
pc.coeffs.Set(i, 0, ys[i])
|
||||
// a_1
|
||||
pc.coeffs.Set(i, 1, dydxs[i])
|
||||
// Solve a linear equation system for a_2 and a_3.
|
||||
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.lastY = ys[m]
|
||||
pc.lastDyDx = dydxs[m]
|
||||
}
|
||||
|
||||
// AkimaSpline 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.
|
||||
// See https://www.iue.tuwien.ac.at/phd/rottinger/node60.html for more details.
|
||||
type AkimaSpline struct {
|
||||
cubic PiecewiseCubic
|
||||
}
|
||||
|
||||
// Predict returns the interpolation value at x.
|
||||
func (as *AkimaSpline) Predict(x float64) float64 {
|
||||
return as.cubic.Predict(x)
|
||||
}
|
||||
|
||||
// PredictDerivative returns the predicted derivative at x.
|
||||
func (as *AkimaSpline) PredictDerivative(x float64) float64 {
|
||||
return as.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 (as *AkimaSpline) Fit(xs, ys []float64) error {
|
||||
n := len(xs)
|
||||
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
|
||||
as.cubic.FitWithDerivatives(xs, ys, dydxs)
|
||||
return nil
|
||||
}
|
||||
slopes := akimaSlopes(xs, ys)
|
||||
for i := 0; i < n; i++ {
|
||||
wLeft, wRight := akimaWeights(slopes, i)
|
||||
dydxs[i] = akimaWeightedAverage(slopes[i+1], slopes[i+2], wLeft, wRight)
|
||||
}
|
||||
as.cubic.FitWithDerivatives(xs, ys, dydxs)
|
||||
return nil
|
||||
}
|
||||
|
||||
// akimaSlopes returns slopes for Akima spline method, including the approximations
|
||||
// of slopes outside the data range (two on each side).
|
||||
// It panics if len(xs) <= 2, elements of xs are not strictly increasing
|
||||
// or len(xs) != len(ys).
|
||||
func akimaSlopes(xs, ys []float64) []float64 {
|
||||
n := len(xs)
|
||||
if n <= 2 {
|
||||
panic(tooFewPoints)
|
||||
}
|
||||
if len(ys) != n {
|
||||
panic(differentLengths)
|
||||
}
|
||||
m := n + 3
|
||||
slopes := make([]float64, m)
|
||||
for i := 2; i < m-2; i++ {
|
||||
dx := xs[i-1] - xs[i-2]
|
||||
if dx <= 0 {
|
||||
panic(xsNotStrictlyIncreasing)
|
||||
}
|
||||
slopes[i] = (ys[i-1] - ys[i-2]) / dx
|
||||
}
|
||||
slopes[0] = 3*slopes[2] - 2*slopes[3]
|
||||
slopes[1] = 2*slopes[2] - slopes[3]
|
||||
slopes[m-2] = 2*slopes[m-3] - slopes[m-4]
|
||||
slopes[m-1] = 3*slopes[m-3] - 2*slopes[m-4]
|
||||
return slopes
|
||||
}
|
||||
|
||||
// akimaWeightedAverage returns (v1 * w1 + v2 * w2) / (w1 + w2) for w1, w2 >= 0 (not checked).
|
||||
// If w1 == w2 == 0, it returns a simple average of v1 and v2.
|
||||
func akimaWeightedAverage(v1, v2, w1, w2 float64) float64 {
|
||||
w := w1 + w2
|
||||
if w > 0 {
|
||||
return (v1*w1 + v2*w2) / w
|
||||
}
|
||||
return 0.5*v1 + 0.5*v2
|
||||
}
|
||||
|
||||
// akimaWeights returns the left and right weight for approximating
|
||||
// the i-th derivative with neighbouring slopes.
|
||||
func akimaWeights(slopes []float64, i int) (float64, float64) {
|
||||
wLeft := math.Abs(slopes[i+2] - slopes[i+3])
|
||||
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
|
||||
}
|
||||
673
interp/cubic_test.go
Normal file
673
interp/cubic_test.go
Normal file
@@ -0,0 +1,673 @@
|
||||
// Copyright ©2020 The Gonum Authors. All rights reserved.
|
||||
// Use of this source code is governed by a BSD-style
|
||||
// license that can be found in the LICENSE file.
|
||||
|
||||
package interp
|
||||
|
||||
import (
|
||||
"math"
|
||||
"testing"
|
||||
|
||||
"gonum.org/v1/gonum/floats"
|
||||
"gonum.org/v1/gonum/floats/scalar"
|
||||
"gonum.org/v1/gonum/mat"
|
||||
)
|
||||
|
||||
func TestPiecewiseCubic(t *testing.T) {
|
||||
t.Parallel()
|
||||
const (
|
||||
h = 1e-8
|
||||
valueTol = 1e-13
|
||||
derivTol = 1e-6
|
||||
nPts = 100
|
||||
)
|
||||
for i, test := range []struct {
|
||||
xs []float64
|
||||
f func(float64) float64
|
||||
df func(float64) float64
|
||||
}{
|
||||
{
|
||||
xs: []float64{-1.001, 0.2, 2},
|
||||
f: func(x float64) float64 { return x * x },
|
||||
df: func(x float64) float64 { return 2 * x },
|
||||
},
|
||||
{
|
||||
xs: []float64{-1.2, -1.001, 0, 0.2, 2.01, 2.1},
|
||||
f: func(x float64) float64 { return 4*math.Pow(x, 3) - 2*x*x + 10*x - 7 },
|
||||
df: func(x float64) float64 { return 12*x*x - 4*x + 10 },
|
||||
},
|
||||
{
|
||||
xs: []float64{-1.001, 0.2, 10},
|
||||
f: func(x float64) float64 { return 1.5*x - 1 },
|
||||
df: func(x float64) float64 { return 1.5 },
|
||||
},
|
||||
{
|
||||
xs: []float64{-1.001, 0.2, 10},
|
||||
f: func(x float64) float64 { return -1 },
|
||||
df: func(x float64) float64 { return 0 },
|
||||
},
|
||||
} {
|
||||
ys := applyFunc(test.xs, test.f)
|
||||
dydxs := applyFunc(test.xs, test.df)
|
||||
var pc PiecewiseCubic
|
||||
pc.FitWithDerivatives(test.xs, ys, dydxs)
|
||||
n := len(test.xs)
|
||||
m := n - 1
|
||||
x0 := test.xs[0]
|
||||
x1 := test.xs[m]
|
||||
x := x0 - 0.1
|
||||
got := pc.Predict(x)
|
||||
want := ys[0]
|
||||
if got != want {
|
||||
t.Errorf("Mismatch in value extrapolated to the left for test case %d: got %v, want %g", i, got, want)
|
||||
}
|
||||
got = pc.PredictDerivative(x)
|
||||
want = dydxs[0]
|
||||
if got != want {
|
||||
t.Errorf("Mismatch in derivative extrapolated to the left for test case %d: got %v, want %g", i, got, want)
|
||||
}
|
||||
x = x1 + 0.1
|
||||
got = pc.Predict(x)
|
||||
want = ys[m]
|
||||
if got != want {
|
||||
t.Errorf("Mismatch in value extrapolated to the right for test case %d: got %v, want %g", i, got, want)
|
||||
}
|
||||
got = pc.PredictDerivative(x)
|
||||
want = dydxs[m]
|
||||
if got != want {
|
||||
t.Errorf("Mismatch in derivative extrapolated to the right for test case %d: got %v, want %g", i, got, want)
|
||||
}
|
||||
for j := 0; j < n; j++ {
|
||||
x := test.xs[j]
|
||||
got := pc.Predict(x)
|
||||
want := test.f(x)
|
||||
if math.Abs(got-want) > valueTol {
|
||||
t.Errorf("Mismatch in interpolated value at x == %g for test case %d: got %v, want %g", x, i, got, want)
|
||||
}
|
||||
if j < m {
|
||||
got = pc.coeffs.At(j, 0)
|
||||
if math.Abs(got-want) > valueTol {
|
||||
t.Errorf("Mismatch in 0-th order interpolation coefficient in %d-th node for test case %d: got %v, want %g", j, i, got, want)
|
||||
}
|
||||
dx := (test.xs[j+1] - x) / nPts
|
||||
for k := 1; k < nPts; k++ {
|
||||
xk := x + float64(k)*dx
|
||||
got := pc.Predict(xk)
|
||||
want := test.f(xk)
|
||||
if math.Abs(got-want) > valueTol {
|
||||
t.Errorf("Mismatch in interpolated value at x == %g for test case %d: got %v, want %g", x, i, got, want)
|
||||
}
|
||||
got = pc.PredictDerivative(xk)
|
||||
want = discrDerivPredict(&pc, x0, x1, xk, h)
|
||||
if math.Abs(got-want) > derivTol {
|
||||
t.Errorf("Mismatch in interpolated derivative at x == %g for test case %d: got %v, want %g", x, i, got, want)
|
||||
}
|
||||
}
|
||||
} else {
|
||||
got = pc.lastY
|
||||
if math.Abs(got-want) > valueTol {
|
||||
t.Errorf("Mismatch in lastY for test case %d: got %v, want %g", i, got, want)
|
||||
}
|
||||
}
|
||||
|
||||
if j > 0 {
|
||||
dx := test.xs[j] - test.xs[j-1]
|
||||
got = ((pc.coeffs.At(j-1, 3)*dx+pc.coeffs.At(j-1, 2))*dx+pc.coeffs.At(j-1, 1))*dx + pc.coeffs.At(j-1, 0)
|
||||
if math.Abs(got-want) > valueTol {
|
||||
t.Errorf("Interpolation coefficients in %d-th node produce mismatch in interpolated value at %g for test case %d: got %v, want %g", j-1, x, i, got, want)
|
||||
}
|
||||
}
|
||||
got = discrDerivPredict(&pc, x0, x1, x, h)
|
||||
want = test.df(x)
|
||||
if math.Abs(got-want) > derivTol {
|
||||
t.Errorf("Mismatch in numerical derivative of interpolated function at x == %g for test case %d: got %v, want %g", x, i, got, want)
|
||||
}
|
||||
got = pc.PredictDerivative(x)
|
||||
if math.Abs(got-want) > valueTol {
|
||||
t.Errorf("Mismatch in interpolated derivative value at x == %g for test case %d: got %v, want %g", x, i, got, want)
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
func TestPiecewiseCubicFitWithDerivatives(t *testing.T) {
|
||||
t.Parallel()
|
||||
xs := []float64{-1, 0, 1}
|
||||
ys := make([]float64, 3)
|
||||
dydxs := make([]float64, 3)
|
||||
leftPoly := func(x float64) float64 {
|
||||
return x*x - x + 1
|
||||
}
|
||||
leftPolyDerivative := func(x float64) float64 {
|
||||
return 2*x - 1
|
||||
}
|
||||
rightPoly := func(x float64) float64 {
|
||||
return x*x*x - x + 1
|
||||
}
|
||||
rightPolyDerivative := func(x float64) float64 {
|
||||
return 3*x*x - 1
|
||||
}
|
||||
ys[0] = leftPoly(xs[0])
|
||||
ys[1] = leftPoly(xs[1])
|
||||
ys[2] = rightPoly(xs[2])
|
||||
dydxs[0] = leftPolyDerivative(xs[0])
|
||||
dydxs[1] = leftPolyDerivative(xs[1])
|
||||
dydxs[2] = rightPolyDerivative(xs[2])
|
||||
var pc PiecewiseCubic
|
||||
pc.FitWithDerivatives(xs, ys, dydxs)
|
||||
lastY := rightPoly(xs[2])
|
||||
if pc.lastY != lastY {
|
||||
t.Errorf("Mismatch in lastY: got %v, want %g", pc.lastY, lastY)
|
||||
}
|
||||
lastDyDx := rightPolyDerivative(xs[2])
|
||||
if pc.lastDyDx != lastDyDx {
|
||||
t.Errorf("Mismatch in lastDxDy: got %v, want %g", pc.lastDyDx, lastDyDx)
|
||||
}
|
||||
if !floats.Equal(pc.xs, xs) {
|
||||
t.Errorf("Mismatch in xs: got %v, want %v", pc.xs, xs)
|
||||
}
|
||||
coeffs := mat.NewDense(2, 4, []float64{3, -3, 1, 0, 1, -1, 0, 1})
|
||||
if !mat.EqualApprox(&pc.coeffs, coeffs, 1e-14) {
|
||||
t.Errorf("Mismatch in coeffs: got %v, want %v", pc.coeffs, coeffs)
|
||||
}
|
||||
}
|
||||
|
||||
func TestPiecewiseCubicFitWithDerivativesErrors(t *testing.T) {
|
||||
t.Parallel()
|
||||
for _, test := range []struct {
|
||||
xs, ys, dydxs []float64
|
||||
}{
|
||||
{
|
||||
xs: []float64{0, 1, 2},
|
||||
ys: []float64{10, 20},
|
||||
dydxs: []float64{0, 0, 0},
|
||||
},
|
||||
{
|
||||
xs: []float64{0, 1, 1},
|
||||
ys: []float64{10, 20, 30},
|
||||
dydxs: []float64{0, 0, 0, 0},
|
||||
},
|
||||
{
|
||||
xs: []float64{0},
|
||||
ys: []float64{0},
|
||||
dydxs: []float64{0},
|
||||
},
|
||||
{
|
||||
xs: []float64{0, 1, 1},
|
||||
ys: []float64{10, 20, 10},
|
||||
dydxs: []float64{0, 0, 0},
|
||||
},
|
||||
} {
|
||||
var pc PiecewiseCubic
|
||||
if !panics(func() { pc.FitWithDerivatives(test.xs, test.ys, test.dydxs) }) {
|
||||
t.Errorf("expected panic for xs: %v, ys: %v and dydxs: %v", test.xs, test.ys, test.dydxs)
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
func TestAkimaSpline(t *testing.T) {
|
||||
t.Parallel()
|
||||
const (
|
||||
derivAbsTol = 1e-8
|
||||
derivRelTol = 1e-7
|
||||
h = 1e-8
|
||||
nPts = 100
|
||||
tol = 1e-14
|
||||
)
|
||||
for i, test := range []struct {
|
||||
xs []float64
|
||||
f func(float64) float64
|
||||
}{
|
||||
{
|
||||
xs: []float64{-5, -3, -2, -1.5, -1, 0.5, 1.5, 2.5, 3},
|
||||
f: func(x float64) float64 { return x * x },
|
||||
},
|
||||
{
|
||||
xs: []float64{-5, -3, -2, -1.5, -1, 0.5, 1.5, 2.5, 3},
|
||||
f: func(x float64) float64 { return math.Pow(x, 3.) - x*x + 2 },
|
||||
},
|
||||
{
|
||||
xs: []float64{-5, -3, -2, -1.5, -1, 0.5, 1.5, 2.5, 3},
|
||||
f: func(x float64) float64 { return -10 * x },
|
||||
},
|
||||
{
|
||||
xs: []float64{-5, -3, -2, -1.5, -1, 0.5, 1.5, 2.5, 3},
|
||||
f: math.Sin,
|
||||
},
|
||||
{
|
||||
xs: []float64{0, 1},
|
||||
f: math.Exp,
|
||||
},
|
||||
{
|
||||
xs: []float64{-1, 0.5},
|
||||
f: math.Cos,
|
||||
},
|
||||
} {
|
||||
var as AkimaSpline
|
||||
n := len(test.xs)
|
||||
m := n - 1
|
||||
x0 := test.xs[0]
|
||||
x1 := test.xs[m]
|
||||
ys := applyFunc(test.xs, test.f)
|
||||
err := as.Fit(test.xs, ys)
|
||||
if err != nil {
|
||||
t.Errorf("Error when fitting AkimaSpline in test case %d: %v", i, err)
|
||||
}
|
||||
for j := 0; j < n; j++ {
|
||||
x := test.xs[j]
|
||||
got := as.Predict(x)
|
||||
want := test.f(x)
|
||||
if math.Abs(got-want) > tol {
|
||||
t.Errorf("Mismatch in interpolated value at x == %g for test case %d: got %v, want %g", x, i, got, want)
|
||||
}
|
||||
if j < m {
|
||||
dx := (test.xs[j+1] - x) / nPts
|
||||
for k := 1; k < nPts; k++ {
|
||||
xk := x + float64(k)*dx
|
||||
got = as.PredictDerivative(xk)
|
||||
want = discrDerivPredict(&as, x0, x1, xk, h)
|
||||
if math.Abs(got-want) > derivRelTol*math.Abs(want)+derivAbsTol {
|
||||
t.Errorf("Mismatch in interpolated derivative at x == %g for test case %d: got %v, want %g", x, i, got, want)
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
if n == 2 {
|
||||
got := as.cubic.coeffs.At(0, 1)
|
||||
want := (ys[1] - ys[0]) / (test.xs[1] - test.xs[0])
|
||||
if math.Abs(got-want) > tol {
|
||||
t.Errorf("Mismatch in approximated slope for length-2 test case %d: got %v, want %g", i, got, want)
|
||||
}
|
||||
for j := 2; i < 4; j++ {
|
||||
got := as.cubic.coeffs.At(0, j)
|
||||
if got != 0 {
|
||||
t.Errorf("Non-zero order-%d coefficient for length-2 test case %d: got %v", j, i, got)
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
func TestAkimaSplineFitErrors(t *testing.T) {
|
||||
t.Parallel()
|
||||
for _, test := range []struct {
|
||||
xs, ys []float64
|
||||
}{
|
||||
{
|
||||
xs: []float64{0, 1, 2},
|
||||
ys: []float64{10, 20},
|
||||
},
|
||||
{
|
||||
xs: []float64{0, 1},
|
||||
ys: []float64{10, 20, 30},
|
||||
},
|
||||
{
|
||||
xs: []float64{0},
|
||||
ys: []float64{0},
|
||||
},
|
||||
{
|
||||
xs: []float64{0, 1, 1},
|
||||
ys: []float64{10, 20, 10},
|
||||
},
|
||||
{
|
||||
xs: []float64{0, 2, 1},
|
||||
ys: []float64{10, 20, 10},
|
||||
},
|
||||
{
|
||||
xs: []float64{0, 0},
|
||||
ys: []float64{-1, 2},
|
||||
},
|
||||
{
|
||||
xs: []float64{0, -1},
|
||||
ys: []float64{-1, 2},
|
||||
},
|
||||
} {
|
||||
var as AkimaSpline
|
||||
if !panics(func() { _ = as.Fit(test.xs, test.ys) }) {
|
||||
t.Errorf("expected panic for xs: %v and ys: %v", test.xs, test.ys)
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
func TestAkimaWeightedAverage(t *testing.T) {
|
||||
t.Parallel()
|
||||
for i, test := range []struct {
|
||||
v1, v2, w1, w2, want float64
|
||||
// "want" values calculated by hand.
|
||||
}{
|
||||
{
|
||||
v1: -1,
|
||||
v2: 1,
|
||||
w1: 0,
|
||||
w2: 0,
|
||||
want: 0,
|
||||
},
|
||||
{
|
||||
v1: -1,
|
||||
v2: 1,
|
||||
w1: 1e6,
|
||||
w2: 1e6,
|
||||
want: 0,
|
||||
},
|
||||
{
|
||||
v1: -1,
|
||||
v2: 1,
|
||||
w1: 1e-10,
|
||||
w2: 0,
|
||||
want: -1,
|
||||
},
|
||||
{
|
||||
v1: -1,
|
||||
v2: 1,
|
||||
w1: 0,
|
||||
w2: 1e-10,
|
||||
want: 1,
|
||||
},
|
||||
{
|
||||
v1: 0,
|
||||
v2: 1000,
|
||||
w1: 1e-13,
|
||||
w2: 3e-13,
|
||||
want: 750,
|
||||
},
|
||||
{
|
||||
v1: 0,
|
||||
v2: 1000,
|
||||
w1: 3e-13,
|
||||
w2: 1e-13,
|
||||
want: 250,
|
||||
},
|
||||
} {
|
||||
got := akimaWeightedAverage(test.v1, test.v2, test.w1, test.w2)
|
||||
if !scalar.EqualWithinAbsOrRel(got, test.want, 1e-14, 1e-14) {
|
||||
t.Errorf("Mismatch in test case %d: got %v, want %g", i, got, test.want)
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
func TestAkimaSlopes(t *testing.T) {
|
||||
t.Parallel()
|
||||
for i, test := range []struct {
|
||||
xs, ys, want []float64
|
||||
// "want" values calculated by hand.
|
||||
}{
|
||||
{
|
||||
xs: []float64{-2, 0, 1},
|
||||
ys: []float64{2, 0, 1.5},
|
||||
want: []float64{-6, -3.5, -1, 1.5, 4, 6.5},
|
||||
},
|
||||
{
|
||||
xs: []float64{-2, -0.5, 1},
|
||||
ys: []float64{-2, -0.5, 1},
|
||||
want: []float64{1, 1, 1, 1, 1, 1},
|
||||
},
|
||||
{
|
||||
xs: []float64{-2, -0.5, 1},
|
||||
ys: []float64{1, 1, 1},
|
||||
want: []float64{0, 0, 0, 0, 0, 0},
|
||||
},
|
||||
{
|
||||
xs: []float64{0, 1.5, 2, 4, 4.5, 5, 6, 7.5, 8},
|
||||
ys: []float64{-5, -4, -3.5, -3.25, -3.25, -2.5, -1.5, -1, 2},
|
||||
want: []float64{0, 1. / 3, 2. / 3, 1, 0.125, 0, 1.5, 1, 1. / 3, 6, 12 - 1./3, 18 - 2./3},
|
||||
},
|
||||
} {
|
||||
got := akimaSlopes(test.xs, test.ys)
|
||||
if !floats.EqualApprox(got, test.want, 1e-14) {
|
||||
t.Errorf("Mismatch in test case %d: got %v, want %v", i, got, test.want)
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
func TestAkimaSlopesErrors(t *testing.T) {
|
||||
t.Parallel()
|
||||
for _, test := range []struct {
|
||||
xs, ys []float64
|
||||
}{
|
||||
{
|
||||
xs: []float64{0, 1, 2},
|
||||
ys: []float64{10, 20},
|
||||
},
|
||||
{
|
||||
xs: []float64{0, 1},
|
||||
ys: []float64{10, 20, 30},
|
||||
},
|
||||
{
|
||||
xs: []float64{0, 2},
|
||||
ys: []float64{0, 1},
|
||||
},
|
||||
{
|
||||
xs: []float64{0, 1, 1},
|
||||
ys: []float64{10, 20, 10},
|
||||
},
|
||||
{
|
||||
xs: []float64{0, 2, 1},
|
||||
ys: []float64{10, 20, 10},
|
||||
},
|
||||
{
|
||||
xs: []float64{0, 0},
|
||||
ys: []float64{-1, 2},
|
||||
},
|
||||
{
|
||||
xs: []float64{0, -1},
|
||||
ys: []float64{-1, 2},
|
||||
},
|
||||
} {
|
||||
if !panics(func() { akimaSlopes(test.xs, test.ys) }) {
|
||||
t.Errorf("expected panic for xs: %v and ys: %v", test.xs, test.ys)
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
func TestAkimaWeights(t *testing.T) {
|
||||
t.Parallel()
|
||||
const tol = 1e-14
|
||||
slopes := []float64{-2, -1, -0.1, 0.2, 1.2, 2.5}
|
||||
// "want" values calculated by hand.
|
||||
want := [][]float64{
|
||||
{0.3, 1},
|
||||
{1, 0.9},
|
||||
{1.3, 0.3},
|
||||
}
|
||||
for i := 0; i < len(want); i++ {
|
||||
gotLeft, gotRight := akimaWeights(slopes, i)
|
||||
if math.Abs(gotLeft-want[i][0]) > tol {
|
||||
t.Errorf("Mismatch in left weight for node %d: got %v, want %g", i, gotLeft, want[i][0])
|
||||
}
|
||||
if math.Abs(gotRight-want[i][1]) > tol {
|
||||
t.Errorf("Mismatch in left weight for node %d: got %v, want %g", i, gotRight, want[i][1])
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
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)
|
||||
}
|
||||
}
|
||||
}
|
||||
295
interp/interp.go
295
interp/interp.go
@@ -4,12 +4,7 @@
|
||||
|
||||
package interp
|
||||
|
||||
import (
|
||||
"math"
|
||||
"sort"
|
||||
|
||||
"gonum.org/v1/gonum/mat"
|
||||
)
|
||||
import "sort"
|
||||
|
||||
const (
|
||||
differentLengths = "interp: input slices have different lengths"
|
||||
@@ -160,179 +155,6 @@ func (pc PiecewiseConstant) Predict(x float64) float64 {
|
||||
return pc.ys[i+1]
|
||||
}
|
||||
|
||||
// PiecewiseCubic is a piecewise cubic 1-dimensional interpolator with
|
||||
// continuous value and first derivative.
|
||||
type PiecewiseCubic struct {
|
||||
// Interpolated X values.
|
||||
xs []float64
|
||||
|
||||
// Coefficients of interpolating cubic polynomials, with
|
||||
// len(xs) - 1 rows and 4 columns. The interpolated value
|
||||
// for xs[i] <= x < xs[i + 1] is defined as
|
||||
// sum_{k = 0}^3 coeffs.At(i, k) * (x - xs[i])^k
|
||||
// To guarantee left-continuity, coeffs.At(i, 0) == ys[i].
|
||||
coeffs mat.Dense
|
||||
|
||||
// Last interpolated Y value, corresponding to xs[len(xs) - 1].
|
||||
lastY float64
|
||||
|
||||
// Last interpolated dY/dX value, corresponding to xs[len(xs) - 1].
|
||||
lastDyDx float64
|
||||
}
|
||||
|
||||
// Predict returns the interpolation value at x.
|
||||
func (pc *PiecewiseCubic) Predict(x float64) float64 {
|
||||
i := findSegment(pc.xs, x)
|
||||
if i < 0 {
|
||||
return pc.coeffs.At(0, 0)
|
||||
}
|
||||
m := len(pc.xs) - 1
|
||||
if x == pc.xs[i] {
|
||||
if i < m {
|
||||
return pc.coeffs.At(i, 0)
|
||||
}
|
||||
return pc.lastY
|
||||
}
|
||||
if i == m {
|
||||
return pc.lastY
|
||||
}
|
||||
dx := x - pc.xs[i]
|
||||
a := pc.coeffs.RawRowView(i)
|
||||
return ((a[3]*dx+a[2])*dx+a[1])*dx + a[0]
|
||||
}
|
||||
|
||||
// PredictDerivative returns the predicted derivative at x.
|
||||
func (pc *PiecewiseCubic) PredictDerivative(x float64) float64 {
|
||||
i := findSegment(pc.xs, x)
|
||||
if i < 0 {
|
||||
return pc.coeffs.At(0, 1)
|
||||
}
|
||||
m := len(pc.xs) - 1
|
||||
if x == pc.xs[i] {
|
||||
if i < m {
|
||||
return pc.coeffs.At(i, 1)
|
||||
}
|
||||
return pc.lastDyDx
|
||||
}
|
||||
if i == m {
|
||||
return pc.lastDyDx
|
||||
}
|
||||
dx := x - pc.xs[i]
|
||||
a := pc.coeffs.RawRowView(i)
|
||||
return (3*a[3]*dx+2*a[2])*dx + a[1]
|
||||
}
|
||||
|
||||
// FitWithDerivatives fits a piecewise cubic predictor to (X, Y, dY/dX) value
|
||||
// triples provided as three slices.
|
||||
// It panics if len(xs) < 2, elements of xs are not strictly increasing,
|
||||
// len(xs) != len(ys) or len(xs) != len(dydxs).
|
||||
func (pc *PiecewiseCubic) FitWithDerivatives(xs, ys, dydxs []float64) {
|
||||
n := len(xs)
|
||||
if len(ys) != n {
|
||||
panic(differentLengths)
|
||||
}
|
||||
if len(dydxs) != n {
|
||||
panic(differentLengths)
|
||||
}
|
||||
if 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]
|
||||
// a_0
|
||||
pc.coeffs.Set(i, 0, ys[i])
|
||||
// a_1
|
||||
pc.coeffs.Set(i, 1, dydxs[i])
|
||||
// Solve a linear equation system for a_2 and a_3.
|
||||
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.lastY = ys[m]
|
||||
pc.lastDyDx = dydxs[m]
|
||||
}
|
||||
|
||||
// AkimaSpline 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.
|
||||
// See https://www.iue.tuwien.ac.at/phd/rottinger/node60.html for more details.
|
||||
type AkimaSpline struct {
|
||||
cubic PiecewiseCubic
|
||||
}
|
||||
|
||||
// Predict returns the interpolation value at x.
|
||||
func (as *AkimaSpline) Predict(x float64) float64 {
|
||||
return as.cubic.Predict(x)
|
||||
}
|
||||
|
||||
// PredictDerivative returns the predicted derivative at x.
|
||||
func (as *AkimaSpline) PredictDerivative(x float64) float64 {
|
||||
return as.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 (as *AkimaSpline) Fit(xs, ys []float64) error {
|
||||
n := len(xs)
|
||||
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
|
||||
as.cubic.FitWithDerivatives(xs, ys, dydxs)
|
||||
return nil
|
||||
}
|
||||
slopes := akimaSlopes(xs, ys)
|
||||
for i := 0; i < n; i++ {
|
||||
wLeft, wRight := akimaWeights(slopes, i)
|
||||
dydxs[i] = akimaWeightedAverage(slopes[i+1], slopes[i+2], wLeft, wRight)
|
||||
}
|
||||
as.cubic.FitWithDerivatives(xs, ys, dydxs)
|
||||
return nil
|
||||
}
|
||||
|
||||
// akimaSlopes returns slopes for Akima spline method, including the approximations
|
||||
// of slopes outside the data range (two on each side).
|
||||
// It panics if len(xs) <= 2, elements of xs are not strictly increasing
|
||||
// or len(xs) != len(ys).
|
||||
func akimaSlopes(xs, ys []float64) []float64 {
|
||||
n := len(xs)
|
||||
if n <= 2 {
|
||||
panic(tooFewPoints)
|
||||
}
|
||||
if len(ys) != n {
|
||||
panic(differentLengths)
|
||||
}
|
||||
m := n + 3
|
||||
slopes := make([]float64, m)
|
||||
for i := 2; i < m-2; i++ {
|
||||
dx := xs[i-1] - xs[i-2]
|
||||
if dx <= 0 {
|
||||
panic(xsNotStrictlyIncreasing)
|
||||
}
|
||||
slopes[i] = (ys[i-1] - ys[i-2]) / dx
|
||||
}
|
||||
slopes[0] = 3*slopes[2] - 2*slopes[3]
|
||||
slopes[1] = 2*slopes[2] - slopes[3]
|
||||
slopes[m-2] = 2*slopes[m-3] - slopes[m-4]
|
||||
slopes[m-1] = 3*slopes[m-3] - 2*slopes[m-4]
|
||||
return slopes
|
||||
}
|
||||
|
||||
// findSegment returns 0 <= i < len(xs) such that xs[i] <= x < xs[i + 1], where xs[len(xs)]
|
||||
// is assumed to be +Inf. If no such i is found, it returns -1. It assumes that len(xs) >= 2
|
||||
// without checking.
|
||||
@@ -340,121 +162,6 @@ func findSegment(xs []float64, x float64) int {
|
||||
return sort.Search(len(xs), func(i int) bool { return xs[i] > x }) - 1
|
||||
}
|
||||
|
||||
// akimaWeightedAverage returns (v1 * w1 + v2 * w2) / (w1 + w2) for w1, w2 >= 0 (not checked).
|
||||
// If w1 == w2 == 0, it returns a simple average of v1 and v2.
|
||||
func akimaWeightedAverage(v1, v2, w1, w2 float64) float64 {
|
||||
w := w1 + w2
|
||||
if w > 0 {
|
||||
return (v1*w1 + v2*w2) / w
|
||||
}
|
||||
return 0.5*v1 + 0.5*v2
|
||||
}
|
||||
|
||||
// akimaWeights returns the left and right weight for approximating
|
||||
// the i-th derivative with neighbouring slopes.
|
||||
func akimaWeights(slopes []float64, i int) (float64, float64) {
|
||||
wLeft := math.Abs(slopes[i+2] - slopes[i+3])
|
||||
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).
|
||||
|
||||
@@ -9,8 +9,6 @@ import (
|
||||
"testing"
|
||||
|
||||
"gonum.org/v1/gonum/floats"
|
||||
"gonum.org/v1/gonum/floats/scalar"
|
||||
"gonum.org/v1/gonum/mat"
|
||||
)
|
||||
|
||||
func TestConstant(t *testing.T) {
|
||||
@@ -192,665 +190,6 @@ func TestPiecewiseConstantPredict(t *testing.T) {
|
||||
testInterpolatorPredict(t, pc, testXs, leftYs, 0)
|
||||
}
|
||||
|
||||
func TestPiecewiseCubic(t *testing.T) {
|
||||
t.Parallel()
|
||||
const (
|
||||
h = 1e-8
|
||||
valueTol = 1e-13
|
||||
derivTol = 1e-6
|
||||
nPts = 100
|
||||
)
|
||||
for i, test := range []struct {
|
||||
xs []float64
|
||||
f func(float64) float64
|
||||
df func(float64) float64
|
||||
}{
|
||||
{
|
||||
xs: []float64{-1.001, 0.2, 2},
|
||||
f: func(x float64) float64 { return x * x },
|
||||
df: func(x float64) float64 { return 2 * x },
|
||||
},
|
||||
{
|
||||
xs: []float64{-1.2, -1.001, 0, 0.2, 2.01, 2.1},
|
||||
f: func(x float64) float64 { return 4*math.Pow(x, 3) - 2*x*x + 10*x - 7 },
|
||||
df: func(x float64) float64 { return 12*x*x - 4*x + 10 },
|
||||
},
|
||||
{
|
||||
xs: []float64{-1.001, 0.2, 10},
|
||||
f: func(x float64) float64 { return 1.5*x - 1 },
|
||||
df: func(x float64) float64 { return 1.5 },
|
||||
},
|
||||
{
|
||||
xs: []float64{-1.001, 0.2, 10},
|
||||
f: func(x float64) float64 { return -1 },
|
||||
df: func(x float64) float64 { return 0 },
|
||||
},
|
||||
} {
|
||||
ys := applyFunc(test.xs, test.f)
|
||||
dydxs := applyFunc(test.xs, test.df)
|
||||
var pc PiecewiseCubic
|
||||
pc.FitWithDerivatives(test.xs, ys, dydxs)
|
||||
n := len(test.xs)
|
||||
m := n - 1
|
||||
x0 := test.xs[0]
|
||||
x1 := test.xs[m]
|
||||
x := x0 - 0.1
|
||||
got := pc.Predict(x)
|
||||
want := ys[0]
|
||||
if got != want {
|
||||
t.Errorf("Mismatch in value extrapolated to the left for test case %d: got %v, want %g", i, got, want)
|
||||
}
|
||||
got = pc.PredictDerivative(x)
|
||||
want = dydxs[0]
|
||||
if got != want {
|
||||
t.Errorf("Mismatch in derivative extrapolated to the left for test case %d: got %v, want %g", i, got, want)
|
||||
}
|
||||
x = x1 + 0.1
|
||||
got = pc.Predict(x)
|
||||
want = ys[m]
|
||||
if got != want {
|
||||
t.Errorf("Mismatch in value extrapolated to the right for test case %d: got %v, want %g", i, got, want)
|
||||
}
|
||||
got = pc.PredictDerivative(x)
|
||||
want = dydxs[m]
|
||||
if got != want {
|
||||
t.Errorf("Mismatch in derivative extrapolated to the right for test case %d: got %v, want %g", i, got, want)
|
||||
}
|
||||
for j := 0; j < n; j++ {
|
||||
x := test.xs[j]
|
||||
got := pc.Predict(x)
|
||||
want := test.f(x)
|
||||
if math.Abs(got-want) > valueTol {
|
||||
t.Errorf("Mismatch in interpolated value at x == %g for test case %d: got %v, want %g", x, i, got, want)
|
||||
}
|
||||
if j < m {
|
||||
got = pc.coeffs.At(j, 0)
|
||||
if math.Abs(got-want) > valueTol {
|
||||
t.Errorf("Mismatch in 0-th order interpolation coefficient in %d-th node for test case %d: got %v, want %g", j, i, got, want)
|
||||
}
|
||||
dx := (test.xs[j+1] - x) / nPts
|
||||
for k := 1; k < nPts; k++ {
|
||||
xk := x + float64(k)*dx
|
||||
got := pc.Predict(xk)
|
||||
want := test.f(xk)
|
||||
if math.Abs(got-want) > valueTol {
|
||||
t.Errorf("Mismatch in interpolated value at x == %g for test case %d: got %v, want %g", x, i, got, want)
|
||||
}
|
||||
got = pc.PredictDerivative(xk)
|
||||
want = discrDerivPredict(&pc, x0, x1, xk, h)
|
||||
if math.Abs(got-want) > derivTol {
|
||||
t.Errorf("Mismatch in interpolated derivative at x == %g for test case %d: got %v, want %g", x, i, got, want)
|
||||
}
|
||||
}
|
||||
} else {
|
||||
got = pc.lastY
|
||||
if math.Abs(got-want) > valueTol {
|
||||
t.Errorf("Mismatch in lastY for test case %d: got %v, want %g", i, got, want)
|
||||
}
|
||||
}
|
||||
|
||||
if j > 0 {
|
||||
dx := test.xs[j] - test.xs[j-1]
|
||||
got = ((pc.coeffs.At(j-1, 3)*dx+pc.coeffs.At(j-1, 2))*dx+pc.coeffs.At(j-1, 1))*dx + pc.coeffs.At(j-1, 0)
|
||||
if math.Abs(got-want) > valueTol {
|
||||
t.Errorf("Interpolation coefficients in %d-th node produce mismatch in interpolated value at %g for test case %d: got %v, want %g", j-1, x, i, got, want)
|
||||
}
|
||||
}
|
||||
got = discrDerivPredict(&pc, x0, x1, x, h)
|
||||
want = test.df(x)
|
||||
if math.Abs(got-want) > derivTol {
|
||||
t.Errorf("Mismatch in numerical derivative of interpolated function at x == %g for test case %d: got %v, want %g", x, i, got, want)
|
||||
}
|
||||
got = pc.PredictDerivative(x)
|
||||
if math.Abs(got-want) > valueTol {
|
||||
t.Errorf("Mismatch in interpolated derivative value at x == %g for test case %d: got %v, want %g", x, i, got, want)
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
func TestPiecewiseCubicFitWithDerivatives(t *testing.T) {
|
||||
t.Parallel()
|
||||
xs := []float64{-1, 0, 1}
|
||||
ys := make([]float64, 3)
|
||||
dydxs := make([]float64, 3)
|
||||
leftPoly := func(x float64) float64 {
|
||||
return x*x - x + 1
|
||||
}
|
||||
leftPolyDerivative := func(x float64) float64 {
|
||||
return 2*x - 1
|
||||
}
|
||||
rightPoly := func(x float64) float64 {
|
||||
return x*x*x - x + 1
|
||||
}
|
||||
rightPolyDerivative := func(x float64) float64 {
|
||||
return 3*x*x - 1
|
||||
}
|
||||
ys[0] = leftPoly(xs[0])
|
||||
ys[1] = leftPoly(xs[1])
|
||||
ys[2] = rightPoly(xs[2])
|
||||
dydxs[0] = leftPolyDerivative(xs[0])
|
||||
dydxs[1] = leftPolyDerivative(xs[1])
|
||||
dydxs[2] = rightPolyDerivative(xs[2])
|
||||
var pc PiecewiseCubic
|
||||
pc.FitWithDerivatives(xs, ys, dydxs)
|
||||
lastY := rightPoly(xs[2])
|
||||
if pc.lastY != lastY {
|
||||
t.Errorf("Mismatch in lastY: got %v, want %g", pc.lastY, lastY)
|
||||
}
|
||||
lastDyDx := rightPolyDerivative(xs[2])
|
||||
if pc.lastDyDx != lastDyDx {
|
||||
t.Errorf("Mismatch in lastDxDy: got %v, want %g", pc.lastDyDx, lastDyDx)
|
||||
}
|
||||
if !floats.Equal(pc.xs, xs) {
|
||||
t.Errorf("Mismatch in xs: got %v, want %v", pc.xs, xs)
|
||||
}
|
||||
coeffs := mat.NewDense(2, 4, []float64{3, -3, 1, 0, 1, -1, 0, 1})
|
||||
if !mat.EqualApprox(&pc.coeffs, coeffs, 1e-14) {
|
||||
t.Errorf("Mismatch in coeffs: got %v, want %v", pc.coeffs, coeffs)
|
||||
}
|
||||
}
|
||||
|
||||
func TestPiecewiseCubicFitWithDerivativesErrors(t *testing.T) {
|
||||
t.Parallel()
|
||||
for _, test := range []struct {
|
||||
xs, ys, dydxs []float64
|
||||
}{
|
||||
{
|
||||
xs: []float64{0, 1, 2},
|
||||
ys: []float64{10, 20},
|
||||
dydxs: []float64{0, 0, 0},
|
||||
},
|
||||
{
|
||||
xs: []float64{0, 1, 1},
|
||||
ys: []float64{10, 20, 30},
|
||||
dydxs: []float64{0, 0, 0, 0},
|
||||
},
|
||||
{
|
||||
xs: []float64{0},
|
||||
ys: []float64{0},
|
||||
dydxs: []float64{0},
|
||||
},
|
||||
{
|
||||
xs: []float64{0, 1, 1},
|
||||
ys: []float64{10, 20, 10},
|
||||
dydxs: []float64{0, 0, 0},
|
||||
},
|
||||
} {
|
||||
var pc PiecewiseCubic
|
||||
if !panics(func() { pc.FitWithDerivatives(test.xs, test.ys, test.dydxs) }) {
|
||||
t.Errorf("expected panic for xs: %v, ys: %v and dydxs: %v", test.xs, test.ys, test.dydxs)
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
func TestAkimaSpline(t *testing.T) {
|
||||
t.Parallel()
|
||||
const (
|
||||
derivAbsTol = 1e-8
|
||||
derivRelTol = 1e-7
|
||||
h = 1e-8
|
||||
nPts = 100
|
||||
tol = 1e-14
|
||||
)
|
||||
for i, test := range []struct {
|
||||
xs []float64
|
||||
f func(float64) float64
|
||||
}{
|
||||
{
|
||||
xs: []float64{-5, -3, -2, -1.5, -1, 0.5, 1.5, 2.5, 3},
|
||||
f: func(x float64) float64 { return x * x },
|
||||
},
|
||||
{
|
||||
xs: []float64{-5, -3, -2, -1.5, -1, 0.5, 1.5, 2.5, 3},
|
||||
f: func(x float64) float64 { return math.Pow(x, 3.) - x*x + 2 },
|
||||
},
|
||||
{
|
||||
xs: []float64{-5, -3, -2, -1.5, -1, 0.5, 1.5, 2.5, 3},
|
||||
f: func(x float64) float64 { return -10 * x },
|
||||
},
|
||||
{
|
||||
xs: []float64{-5, -3, -2, -1.5, -1, 0.5, 1.5, 2.5, 3},
|
||||
f: math.Sin,
|
||||
},
|
||||
{
|
||||
xs: []float64{0, 1},
|
||||
f: math.Exp,
|
||||
},
|
||||
{
|
||||
xs: []float64{-1, 0.5},
|
||||
f: math.Cos,
|
||||
},
|
||||
} {
|
||||
var as AkimaSpline
|
||||
n := len(test.xs)
|
||||
m := n - 1
|
||||
x0 := test.xs[0]
|
||||
x1 := test.xs[m]
|
||||
ys := applyFunc(test.xs, test.f)
|
||||
err := as.Fit(test.xs, ys)
|
||||
if err != nil {
|
||||
t.Errorf("Error when fitting AkimaSpline in test case %d: %v", i, err)
|
||||
}
|
||||
for j := 0; j < n; j++ {
|
||||
x := test.xs[j]
|
||||
got := as.Predict(x)
|
||||
want := test.f(x)
|
||||
if math.Abs(got-want) > tol {
|
||||
t.Errorf("Mismatch in interpolated value at x == %g for test case %d: got %v, want %g", x, i, got, want)
|
||||
}
|
||||
if j < m {
|
||||
dx := (test.xs[j+1] - x) / nPts
|
||||
for k := 1; k < nPts; k++ {
|
||||
xk := x + float64(k)*dx
|
||||
got = as.PredictDerivative(xk)
|
||||
want = discrDerivPredict(&as, x0, x1, xk, h)
|
||||
if math.Abs(got-want) > derivRelTol*math.Abs(want)+derivAbsTol {
|
||||
t.Errorf("Mismatch in interpolated derivative at x == %g for test case %d: got %v, want %g", x, i, got, want)
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
if n == 2 {
|
||||
got := as.cubic.coeffs.At(0, 1)
|
||||
want := (ys[1] - ys[0]) / (test.xs[1] - test.xs[0])
|
||||
if math.Abs(got-want) > tol {
|
||||
t.Errorf("Mismatch in approximated slope for length-2 test case %d: got %v, want %g", i, got, want)
|
||||
}
|
||||
for j := 2; i < 4; j++ {
|
||||
got := as.cubic.coeffs.At(0, j)
|
||||
if got != 0 {
|
||||
t.Errorf("Non-zero order-%d coefficient for length-2 test case %d: got %v", j, i, got)
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
func TestAkimaSplineFitErrors(t *testing.T) {
|
||||
t.Parallel()
|
||||
for _, test := range []struct {
|
||||
xs, ys []float64
|
||||
}{
|
||||
{
|
||||
xs: []float64{0, 1, 2},
|
||||
ys: []float64{10, 20},
|
||||
},
|
||||
{
|
||||
xs: []float64{0, 1},
|
||||
ys: []float64{10, 20, 30},
|
||||
},
|
||||
{
|
||||
xs: []float64{0},
|
||||
ys: []float64{0},
|
||||
},
|
||||
{
|
||||
xs: []float64{0, 1, 1},
|
||||
ys: []float64{10, 20, 10},
|
||||
},
|
||||
{
|
||||
xs: []float64{0, 2, 1},
|
||||
ys: []float64{10, 20, 10},
|
||||
},
|
||||
{
|
||||
xs: []float64{0, 0},
|
||||
ys: []float64{-1, 2},
|
||||
},
|
||||
{
|
||||
xs: []float64{0, -1},
|
||||
ys: []float64{-1, 2},
|
||||
},
|
||||
} {
|
||||
var as AkimaSpline
|
||||
if !panics(func() { _ = as.Fit(test.xs, test.ys) }) {
|
||||
t.Errorf("expected panic for xs: %v and ys: %v", test.xs, test.ys)
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
func TestAkimaWeightedAverage(t *testing.T) {
|
||||
t.Parallel()
|
||||
for i, test := range []struct {
|
||||
v1, v2, w1, w2, want float64
|
||||
// "want" values calculated by hand.
|
||||
}{
|
||||
{
|
||||
v1: -1,
|
||||
v2: 1,
|
||||
w1: 0,
|
||||
w2: 0,
|
||||
want: 0,
|
||||
},
|
||||
{
|
||||
v1: -1,
|
||||
v2: 1,
|
||||
w1: 1e6,
|
||||
w2: 1e6,
|
||||
want: 0,
|
||||
},
|
||||
{
|
||||
v1: -1,
|
||||
v2: 1,
|
||||
w1: 1e-10,
|
||||
w2: 0,
|
||||
want: -1,
|
||||
},
|
||||
{
|
||||
v1: -1,
|
||||
v2: 1,
|
||||
w1: 0,
|
||||
w2: 1e-10,
|
||||
want: 1,
|
||||
},
|
||||
{
|
||||
v1: 0,
|
||||
v2: 1000,
|
||||
w1: 1e-13,
|
||||
w2: 3e-13,
|
||||
want: 750,
|
||||
},
|
||||
{
|
||||
v1: 0,
|
||||
v2: 1000,
|
||||
w1: 3e-13,
|
||||
w2: 1e-13,
|
||||
want: 250,
|
||||
},
|
||||
} {
|
||||
got := akimaWeightedAverage(test.v1, test.v2, test.w1, test.w2)
|
||||
if !scalar.EqualWithinAbsOrRel(got, test.want, 1e-14, 1e-14) {
|
||||
t.Errorf("Mismatch in test case %d: got %v, want %g", i, got, test.want)
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
func TestAkimaSlopes(t *testing.T) {
|
||||
t.Parallel()
|
||||
for i, test := range []struct {
|
||||
xs, ys, want []float64
|
||||
// "want" values calculated by hand.
|
||||
}{
|
||||
{
|
||||
xs: []float64{-2, 0, 1},
|
||||
ys: []float64{2, 0, 1.5},
|
||||
want: []float64{-6, -3.5, -1, 1.5, 4, 6.5},
|
||||
},
|
||||
{
|
||||
xs: []float64{-2, -0.5, 1},
|
||||
ys: []float64{-2, -0.5, 1},
|
||||
want: []float64{1, 1, 1, 1, 1, 1},
|
||||
},
|
||||
{
|
||||
xs: []float64{-2, -0.5, 1},
|
||||
ys: []float64{1, 1, 1},
|
||||
want: []float64{0, 0, 0, 0, 0, 0},
|
||||
},
|
||||
{
|
||||
xs: []float64{0, 1.5, 2, 4, 4.5, 5, 6, 7.5, 8},
|
||||
ys: []float64{-5, -4, -3.5, -3.25, -3.25, -2.5, -1.5, -1, 2},
|
||||
want: []float64{0, 1. / 3, 2. / 3, 1, 0.125, 0, 1.5, 1, 1. / 3, 6, 12 - 1./3, 18 - 2./3},
|
||||
},
|
||||
} {
|
||||
got := akimaSlopes(test.xs, test.ys)
|
||||
if !floats.EqualApprox(got, test.want, 1e-14) {
|
||||
t.Errorf("Mismatch in test case %d: got %v, want %v", i, got, test.want)
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
func TestAkimaSlopesErrors(t *testing.T) {
|
||||
t.Parallel()
|
||||
for _, test := range []struct {
|
||||
xs, ys []float64
|
||||
}{
|
||||
{
|
||||
xs: []float64{0, 1, 2},
|
||||
ys: []float64{10, 20},
|
||||
},
|
||||
{
|
||||
xs: []float64{0, 1},
|
||||
ys: []float64{10, 20, 30},
|
||||
},
|
||||
{
|
||||
xs: []float64{0, 2},
|
||||
ys: []float64{0, 1},
|
||||
},
|
||||
{
|
||||
xs: []float64{0, 1, 1},
|
||||
ys: []float64{10, 20, 10},
|
||||
},
|
||||
{
|
||||
xs: []float64{0, 2, 1},
|
||||
ys: []float64{10, 20, 10},
|
||||
},
|
||||
{
|
||||
xs: []float64{0, 0},
|
||||
ys: []float64{-1, 2},
|
||||
},
|
||||
{
|
||||
xs: []float64{0, -1},
|
||||
ys: []float64{-1, 2},
|
||||
},
|
||||
} {
|
||||
if !panics(func() { akimaSlopes(test.xs, test.ys) }) {
|
||||
t.Errorf("expected panic for xs: %v and ys: %v", test.xs, test.ys)
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
func TestAkimaWeights(t *testing.T) {
|
||||
t.Parallel()
|
||||
const tol = 1e-14
|
||||
slopes := []float64{-2, -1, -0.1, 0.2, 1.2, 2.5}
|
||||
// "want" values calculated by hand.
|
||||
want := [][]float64{
|
||||
{0.3, 1},
|
||||
{1, 0.9},
|
||||
{1.3, 0.3},
|
||||
}
|
||||
for i := 0; i < len(want); i++ {
|
||||
gotLeft, gotRight := akimaWeights(slopes, i)
|
||||
if math.Abs(gotLeft-want[i][0]) > tol {
|
||||
t.Errorf("Mismatch in left weight for node %d: got %v, want %g", i, gotLeft, want[i][0])
|
||||
}
|
||||
if math.Abs(gotRight-want[i][1]) > tol {
|
||||
t.Errorf("Mismatch in left weight for node %d: got %v, want %g", i, gotRight, want[i][1])
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
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 {
|
||||
|
||||
Reference in New Issue
Block a user