From ee7fd38308ec2698daa822ed6bcdcd9d44319688 Mon Sep 17 00:00:00 2001 From: Sergey Grebenshchikov Date: Sun, 10 Dec 2023 00:51:52 +0100 Subject: [PATCH] add FunctionUniform --- README.md | 58 ++++++++++++++++++++++------ piecewiselinear.go | 77 ++++++++++++++++++++++++++++++++++++-- piecewiselinear_test.go | 83 ++++++++++++++++++++++++++++++++++++++++- span.go | 4 +- 4 files changed, 203 insertions(+), 19 deletions(-) diff --git a/README.md b/README.md index 0b9ad89..7ff5c0d 100644 --- a/README.md +++ b/README.md @@ -10,6 +10,7 @@ import "github.com/sgreben/piecewiselinear" - [Get it](#get-it) - [Use it](#use-it) + - [Fast special case](#fast-special-case) - [Benchmarks](#benchmarks) @@ -26,8 +27,35 @@ import "github.com/sgreben/piecewiselinear" func main() { f := piecewiselinear.Function{Y:[]float64{0,1,0}} // range: "hat" function - f.X = piecewiselinear.Span(0, 1, len(f.Y)) // domain: equidistant points along X axis + f.X = []float64{0, 0.5, 1} // domain: equidistant points along X axis fmt.Println( + f.At(0), // f.At(x) evaluates f at x + f.At(0.25), + f.At(0.5), + f.At(0.75), + f.At(1.0), + f.At(123.0), // outside its domain X the function is constant 0 + f.At(-123.0), // + + f.Area(), + f.AreaUpTo(0.5), + ) + // Output: + // 0 0.5 1 0.5 0 0 0 0.5 0.25 +} +``` + +### Fast special case + +If the control points are uniformly spaced, `piecewiselinear.FunctionUniform` is much faster (no search required): + +```go +import "github.com/sgreben/piecewiselinear" + +func main() { + f := piecewiselinear.FunctionUniform{Y: []float64{0, 1, 0}} // range: "hat" function + f.Xmin, f.Xmax = 0, 1 // domain: equidistant points along X axis + fmt.Println( f.At(0), // f.At(x) evaluates f at x f.At(0.25), f.At(0.5), @@ -35,9 +63,12 @@ func main() { f.At(1.0), f.At(123.0), // outside its domain X the function is constant 0 f.At(-123.0), // + + f.Area(), + f.AreaUpTo(0.5), ) - // Output: - // 0 0.5 1 0.5 0 0 0 + // Output: + // 0 0.5 1 0.5 0 0 0 0.5 0.25 } ``` @@ -48,17 +79,20 @@ On an Apple M1 Pro: - **6ns** per evaluation (`.At(x)`) for 10 control points - **320ns** per evaluation for 10 million control points. +and, for `FunctionUniform`, **2ns** per evaluation regardless of the number of control points. + ``` goos: darwin goarch: arm64 pkg: github.com/sgreben/piecewiselinear -BenchmarkAt4-10 217302646 5.461 ns/op 0 B/op 0 allocs/op -BenchmarkAt8-10 197175420 6.048 ns/op 0 B/op 0 allocs/op -BenchmarkAt10-10 188384818 6.283 ns/op 0 B/op 0 allocs/op -BenchmarkAt100-10 138276301 9.086 ns/op 0 B/op 0 allocs/op -BenchmarkAt1k-10 41258203 25.18 ns/op 0 B/op 0 allocs/op -BenchmarkAt10k-10 16852758 69.99 ns/op 0 B/op 0 allocs/op -BenchmarkAt100k-10 11745384 100.5 ns/op 0 B/op 0 allocs/op -BenchmarkAt1M-10 8501438 143.0 ns/op 0 B/op 0 allocs/op -BenchmarkAt10M-10 3659188 319.5 ns/op 0 B/op 0 allocs/op +BenchmarkAt4-10 230890022 5.499 ns/op 0 B/op 0 allocs/op +BenchmarkAt8-10 199668106 6.084 ns/op 0 B/op 0 allocs/op +BenchmarkAt10-10 192352903 6.206 ns/op 0 B/op 0 allocs/op +BenchmarkAt100-10 138742411 8.613 ns/op 0 B/op 0 allocs/op +BenchmarkAt1k-10 46360660 25.50 ns/op 0 B/op 0 allocs/op +BenchmarkAt10k-10 16649996 70.02 ns/op 0 B/op 0 allocs/op +BenchmarkAt100k-10 11696936 100.4 ns/op 0 B/op 0 allocs/op +BenchmarkAt1M-10 8512652 140.6 ns/op 0 B/op 0 allocs/op +BenchmarkAt10M-10 3769648 320.4 ns/op 0 B/op 0 allocs/op +BenchmarkUniformAt10M-10 571224222 2.185 ns/op 0 B/op 0 allocs/op ``` diff --git a/piecewiselinear.go b/piecewiselinear.go index 5082abc..1f8e804 100644 --- a/piecewiselinear.go +++ b/piecewiselinear.go @@ -1,6 +1,10 @@ // Package piecewiselinear is a tiny library for linear interpolation. package piecewiselinear +import ( + "math" +) + // Function is a piecewise-linear 1-dimensional function type Function struct { X []float64 @@ -71,9 +75,6 @@ func (f Function) At(x float64) float64 { } } if i == 0 { - if len(X) > 0 && x < X[0] { - return 0 - } if len(X) > 0 && x == X[0] { return Y[0] } @@ -85,3 +86,73 @@ func (f Function) At(x float64) float64 { w := (x - X[i-1]) / (X[i] - X[i-1]) return (1-w)*Y[i-1] + w*Y[i] } + +// Function is a piecewise-linear 1-dimensional function with uniformly spaced control points. Faster to interpolate than the equivalent `Function`. +type FunctionUniform struct { + Xmin, Xmax float64 + Y []float64 +} + +// At returns the function's value at the given point. +// Outside its domain X, the function is constant at 0. +// +// Time complexity: O(1) +// Space complexity: O(1) +func (f FunctionUniform) At(x float64) float64 { + if x < f.Xmin { + return 0 + } + if x > f.Xmax { + return 0 + } + step := (f.Xmax - f.Xmin) / float64(len(f.Y)-1) + i := math.Floor(x / step) + j := math.Ceil(x / step) + if i == j { + return f.Y[int(i)] + } + Xi := f.Xmin + (i * step) + Xj := f.Xmin + (j * step) + w := (x - Xi) / (Xj - Xi) + return (1-w)*f.Y[int(i)] + w*f.Y[int(j)] +} + +// Area returns the definite integral of the function on its domain X. +// +// Time complexity: O(N), where N is the number of points. +// Space complexity: O(1) +func (f FunctionUniform) Area() (area float64) { + Y := f.Y + step := (f.Xmax - f.Xmin) / float64(len(f.Y)-1) + for i := 1; i < len(Y); i++ { + Xi := f.Xmin + (float64(i) * step) + Xi1 := f.Xmin + (float64(i-1) * step) + area += (Xi - Xi1) * (Y[i] + Y[i-1]) / 2 + } + return area +} + +// AreaUpTo returns the definite integral of the function on its domain X intersected with [-Inf, x]. +// +// Time complexity: O(N), where N is the number of points. +// Space complexity: O(1) +func (f FunctionUniform) AreaUpTo(x float64) (area float64) { + Y := f.Y + step := (f.Xmax - f.Xmin) / float64(len(f.Y)-1) + for i := 1; i < len(Y); i++ { + Xi := f.Xmin + (float64(i) * step) + Xi1 := f.Xmin + (float64(i-1) * step) + dX := Xi - Xi1 + if x < Xi { + if x >= Xi1 { + dxX := x - Xi1 + w := dxX / dX + y := (1-w)*Y[i-1] + w*Y[i] + area += dxX * (y + Y[i-1]) / 2 + } + return area + } + area += dX * (Y[i] + Y[i-1]) / 2 + } + return area +} diff --git a/piecewiselinear_test.go b/piecewiselinear_test.go index 9498aa6..69b40fd 100644 --- a/piecewiselinear_test.go +++ b/piecewiselinear_test.go @@ -8,7 +8,7 @@ import ( func Example() { f := Function{Y: []float64{0, 1, 0}} // range: "hat" function - f.X = Span(0, 1, len(f.Y)) // domain: equidistant points along X axis + f.X = []float64{0, 0.5, 1} // domain: equidistant points along X axis fmt.Println( f.At(0), // f.At(x) evaluates f at x f.At(0.25), @@ -17,9 +17,31 @@ func Example() { f.At(1.0), f.At(123.0), // outside its domain X the function is constant 0 f.At(-123.0), // + + f.Area(), + f.AreaUpTo(0.5), ) // Output: - // 0 0.5 1 0.5 0 0 0 + // 0 0.5 1 0.5 0 0 0 0.5 0.25 +} + +func ExampleUniform() { + f := FunctionUniform{Y: []float64{0, 1, 0}} // range: "hat" function + f.Xmin, f.Xmax = 0, 1 // domain: equidistant points along X axis + fmt.Println( + f.At(0), // f.At(x) evaluates f at x + f.At(0.25), + f.At(0.5), + f.At(0.75), + f.At(1.0), + f.At(123.0), // outside its domain X the function is constant 0 + f.At(-123.0), // + + f.Area(), + f.AreaUpTo(0.5), + ) + // Output: + // 0 0.5 1 0.5 0 0 0 0.5 0.25 } func TestFunction_AreaUpTo(t *testing.T) { @@ -326,6 +348,43 @@ func TestFunction_IsInterpolatedAt(t *testing.T) { } } +func TestUniform(t *testing.T) { + var f Function + var fU FunctionUniform + ks := []int{2, 3, 4, 5, 8, 9, 16, 17, 32, 33, 64, 65, 128, 129, 256, 257} + for _, k := range ks { + f.Y = make([]float64, k) + fU.Y = f.Y + fU.Xmax = 1 + fU.Xmin = 0 + for i := range f.Y { + f.Y[i] = rand.Float64() + } + f.X = Span(0, 1, len(f.Y)) + xs := make([]float64, len(f.X)) + for i := range f.X { + xs[i] = rand.Float64() + } + for _, x := range xs { + fx := f.At(x) + fUx := fU.At(x) + if fx != fUx { + t.Errorf(".At(): %v != %v (expected)", fUx, fx) + } + Fx := f.AreaUpTo(x) + FUx := fU.AreaUpTo(x) + if Fx != FUx { + t.Errorf(".AreaUpTo(): %v != %v (expected)", FUx, Fx) + } + } + F := f.Area() + FU := fU.Area() + if F != FU { + t.Errorf(".Area(): %v != %v (expected)", F, FU) + } + } +} + func benchmarkWithKPoints(b *testing.B, k int) { var f Function f.Y = make([]float64, k) @@ -343,6 +402,22 @@ func benchmarkWithKPoints(b *testing.B, k int) { } } +func benchmarkUniformWithKPoints(b *testing.B, k int) { + var f FunctionUniform + f.Y = make([]float64, k) + for i := range f.Y { + f.Y[i] = rand.Float64() + } + xs := make([]float64, len(f.Y)) + for i := range f.Y { + xs[i] = rand.Float64() + } + b.ResetTimer() + for n := 0; n < b.N; n++ { + _ = f.At(xs[n%len(f.Y)]) + } +} + func BenchmarkAt4(b *testing.B) { benchmarkWithKPoints(b, 4) } @@ -378,3 +453,7 @@ func BenchmarkAt1M(b *testing.B) { func BenchmarkAt10M(b *testing.B) { benchmarkWithKPoints(b, 10_000_000) } + +func BenchmarkUniformAt10M(b *testing.B) { + benchmarkUniformWithKPoints(b, 10_000_000) +} diff --git a/span.go b/span.go index 8ea0f12..265f783 100644 --- a/span.go +++ b/span.go @@ -6,9 +6,9 @@ import "math" func Span(min, max float64, nPoints int) []float64 { X := make([]float64, nPoints) min, max = math.Min(max, min), math.Max(max, min) - d := max - min + step := (max - min) / float64(nPoints-1) for i := range X { - X[i] = min + d*(float64(i)/float64(nPoints-1)) + X[i] = min + float64(i)*step } return X }