add FunctionUniform

sgreben · Dec 9, 2023 · ee7fd38 · ee7fd38
1 parent 6bce2e9
commit ee7fd38
Show file tree

Hide file tree

Showing 4 changed files with 203 additions and 19 deletions.
diff --git 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,18 +27,48 @@ 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),
 		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
+	// 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
@@ -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
@@ -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
@@ -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
 }