309 lines
9.2 KiB
Go
309 lines
9.2 KiB
Go
|
// Copyright (c) 2026 Petr Balvín <opensource@petrbalvin.org> (https://petrbalvin.org)
|
|||
|
|
// SPDX-License-Identifier: MIT
|
|||
|
|
|
|||
|
|
package tensor_test
|
|||
|
|
|
|||
|
|
// The godoc examples: every flagship workflow as a runnable, checked
|
|||
|
|
// snippet. pkg.go.dev renders these beside the API, and `go test`
|
|||
|
|
// executes them, so the documentation cannot rot.
|
|||
|
|
|
|||
|
|
import (
|
|||
|
|
"fmt"
|
|||
|
|
"log"
|
|||
|
|
"math"
|
|||
|
|
|
|||
|
|
tensor "sourcedock.dev/petrbalvin/tensor"
|
|||
|
|
grad "sourcedock.dev/petrbalvin/tensor/grad"
|
|||
|
|
)
|
|||
|
|
|
|||
|
|
// A rank-1 array from literals, element-wise arithmetic, a reduction.
|
|||
|
|
func ExampleAdd() {
|
|||
|
|
a, err := tensor.FromFloats([]float64{1, 2, 3, 4}, 4)
|
|||
|
|
if err != nil {
|
|||
|
|
log.Fatal(err)
|
|||
|
|
}
|
|||
|
|
b, _ := tensor.FromFloats([]float64{10, 20, 30, 40}, 4)
|
|||
|
|
sum, _ := tensor.Add(a, b)
|
|||
|
|
mean, _ := tensor.Mean(sum)
|
|||
|
|
fmt.Println(sum, mean)
|
|||
|
|
// Output: float (4) [11, 22, 33, 44] 27.5
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// The 2-D matrix product.
|
|||
|
|
func ExampleMatMul2D() {
|
|||
|
|
a, _ := tensor.FromFloats([]float64{1, 2, 3, 4}, 2, 2)
|
|||
|
|
b, _ := tensor.FromFloats([]float64{5, 6, 7, 8}, 2, 2)
|
|||
|
|
prod, err := tensor.MatMul2D(a, b)
|
|||
|
|
if err != nil {
|
|||
|
|
log.Fatal(err)
|
|||
|
|
}
|
|||
|
|
fmt.Println(prod)
|
|||
|
|
// Output: float (2, 2) [19, 22, 43, 50]
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// Einstein summation, batched over an ellipsis axis.
|
|||
|
|
func ExampleEinsum() {
|
|||
|
|
a, _ := tensor.FromFloats([]float64{1, 2, 3, 4, 5, 6, 7, 8}, 2, 2, 2)
|
|||
|
|
b, _ := tensor.FromFloats([]float64{1, 0, 0, 1, 1, 0, 0, 1}, 2, 2, 2)
|
|||
|
|
// Batched matrix product: each batch times its identity.
|
|||
|
|
got, err := tensor.Einsum("...ij,...jk->...ik", a, b)
|
|||
|
|
if err != nil {
|
|||
|
|
log.Fatal(err)
|
|||
|
|
}
|
|||
|
|
fmt.Println(got.Shape(), got.FloatAt(0), got.FloatAt(3))
|
|||
|
|
// Output: [2 2 2] 1 4
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// Reverse-mode autograd: the gradient of Σ (x·w)² arrives exact.
|
|||
|
|
func ExampleTensor_Backward() {
|
|||
|
|
x, _ := grad.FromFloat64s([]float64{1, 2, 3}, true, 3)
|
|||
|
|
w, _ := grad.FromFloat64s([]float64{0.5, -1, 2}, false, 3)
|
|||
|
|
prod, _ := x.Mul(w)
|
|||
|
|
sq, _ := prod.Pow(2)
|
|||
|
|
loss, _ := sq.Sum()
|
|||
|
|
if err := loss.Backward(); err != nil {
|
|||
|
|
log.Fatal(err)
|
|||
|
|
}
|
|||
|
|
// dL/dx = 2·x·w², non-negative because w enters squared.
|
|||
|
|
g := x.Grad()
|
|||
|
|
fmt.Printf("%.4g %.4g %.4g\n", g.FloatAt(0), g.FloatAt(1), g.FloatAt(2))
|
|||
|
|
// Output: 0.5 4 24
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// Ordinary least squares with the full classical inference.
|
|||
|
|
func ExampleLinearRegression() {
|
|||
|
|
// y = 2 + 3x on x = 0..5, the intercept column first.
|
|||
|
|
design, _ := tensor.FromFloats([]float64{
|
|||
|
|
1, 0, 1, 1, 1, 2, 1, 3, 1, 4, 1, 5,
|
|||
|
|
}, 6, 2)
|
|||
|
|
y, _ := tensor.FromFloats([]float64{2.1, 4.9, 8.2, 11.1, 13.8, 17.2}, 6)
|
|||
|
|
res, err := tensor.LinearRegression(design, y)
|
|||
|
|
if err != nil {
|
|||
|
|
log.Fatal(err)
|
|||
|
|
}
|
|||
|
|
fmt.Printf("slope %.3f ± %.3f, R2 %.4f\n",
|
|||
|
|
res.Coefficients[1], res.StandardErrors[1], res.RSquared)
|
|||
|
|
// Output: slope 3.003 ± 0.044, R2 0.9991
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// Solving an initial value problem: the exponential decay y' = −y.
|
|||
|
|
func ExampleIntegrateODE() {
|
|||
|
|
f := func(t float64, y *tensor.Array) (*tensor.Array, error) {
|
|||
|
|
return tensor.MulF(y, -1), nil
|
|||
|
|
}
|
|||
|
|
y0, _ := tensor.FromFloats([]float64{1}, 1)
|
|||
|
|
end, err := tensor.IntegrateODE(f, 0, 1, y0, tensor.ODEOptions{})
|
|||
|
|
if err != nil {
|
|||
|
|
log.Fatal(err)
|
|||
|
|
}
|
|||
|
|
fmt.Printf("y(1) = %.6f\n", end.FloatAt(0))
|
|||
|
|
// Output: y(1) = 0.367880
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// The globally adaptive cubature over a box: the 2-D Gaussian.
|
|||
|
|
func ExampleIntegrateND() {
|
|||
|
|
got, err := tensor.IntegrateND(func(x []float64) float64 {
|
|||
|
|
return math.Exp(-x[0]*x[0] - x[1]*x[1])
|
|||
|
|
}, []float64{-3, -3}, []float64{3, 3}, tensor.CubatureOptions{Tolerance: 1e-11})
|
|||
|
|
if err != nil {
|
|||
|
|
log.Fatal(err)
|
|||
|
|
}
|
|||
|
|
fmt.Printf("%.6f\n", got)
|
|||
|
|
// Output: 3.141454
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// The Crank-Nicolson heat equation holds its eigenmode shape while it
|
|||
|
|
// decays.
|
|||
|
|
func ExampleIntegrateHeat1D() {
|
|||
|
|
const (
|
|||
|
|
n = 49
|
|||
|
|
kappa = 0.1
|
|||
|
|
dx = 1.0 / 50
|
|||
|
|
)
|
|||
|
|
u0 := make([]float64, n)
|
|||
|
|
for i := range n {
|
|||
|
|
u0[i] = math.Sin(math.Pi * float64(i+1) * dx)
|
|||
|
|
}
|
|||
|
|
u0Arr, _ := tensor.FromFloats(u0, n)
|
|||
|
|
states, err := tensor.IntegrateHeat1D(u0Arr, kappa, dx, 1.0, 0.001, 2, 0, 0)
|
|||
|
|
if err != nil {
|
|||
|
|
log.Fatal(err)
|
|||
|
|
}
|
|||
|
|
last := (states.Shape()[0] - 1) * n
|
|||
|
|
// The mode's centre started at 1 and decays by exp(−κπ²t).
|
|||
|
|
decay := math.Exp(-kappa * math.Pi * math.Pi)
|
|||
|
|
fmt.Printf("decay %.4f, centre %.4f\n", decay, states.FloatAt(last+n/2))
|
|||
|
|
// Output: decay 0.3727, centre 0.3728
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// The Lanczos eigensolver on a sparse symmetric matrix.
|
|||
|
|
func ExampleSpEigen() {
|
|||
|
|
dense, _ := tensor.FromFloats([]float64{2, 1, 1, 2}, 2, 2)
|
|||
|
|
sp, _ := tensor.SparseFrom(dense)
|
|||
|
|
vals, _, err := tensor.SpEigen(sp, 2, tensor.NewGenerator(7))
|
|||
|
|
if err != nil {
|
|||
|
|
log.Fatal(err)
|
|||
|
|
}
|
|||
|
|
fmt.Printf("%.6f %.6f\n", vals.FloatAt(0), vals.FloatAt(1))
|
|||
|
|
// Output: 3.000000 1.000000
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// The complete elliptic integral of the first kind by the AGM.
|
|||
|
|
func ExampleEllipticK() {
|
|||
|
|
m, _ := tensor.FromFloats([]float64{0, 0.5}, 2)
|
|||
|
|
k, err := tensor.EllipticK(m)
|
|||
|
|
if err != nil {
|
|||
|
|
log.Fatal(err)
|
|||
|
|
}
|
|||
|
|
fmt.Printf("%.10f %.10f\n", k.FloatAt(0), k.FloatAt(1))
|
|||
|
|
// Output: 1.5707963268 1.8540746773
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// Haar wavelets turn an off-grid step into one detail coefficient,
|
|||
|
|
// and the inverse restores the signal exactly.
|
|||
|
|
func ExampleDWT() {
|
|||
|
|
vals := make([]float64, 16)
|
|||
|
|
for i := 5; i < 16; i++ {
|
|||
|
|
vals[i] = 1
|
|||
|
|
}
|
|||
|
|
x, _ := tensor.FromFloats(vals, 16)
|
|||
|
|
c, err := tensor.DWT(x, 1)
|
|||
|
|
if err != nil {
|
|||
|
|
log.Fatal(err)
|
|||
|
|
}
|
|||
|
|
restored, _ := tensor.IDWT(c, 1)
|
|||
|
|
ok := math.Abs(restored.FloatAt(7)-1) < 1e-12 && restored.FloatAt(4) == 0
|
|||
|
|
fmt.Printf("detail[2] = %.4f, restored: %v\n", c.FloatAt(10), ok)
|
|||
|
|
// Output: detail[2] = -0.7071, restored: true
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// The autocorrelation of a periodic signal is itself periodic.
|
|||
|
|
func ExampleAutocorrelate() {
|
|||
|
|
vals := make([]float64, 64)
|
|||
|
|
for i := range vals {
|
|||
|
|
vals[i] = math.Cos(2 * math.Pi * float64(i) / 16)
|
|||
|
|
}
|
|||
|
|
x, _ := tensor.FromFloats(vals, 64)
|
|||
|
|
acf, err := tensor.Autocorrelate(x, 16)
|
|||
|
|
if err != nil {
|
|||
|
|
log.Fatal(err)
|
|||
|
|
}
|
|||
|
|
// The biased normalisation tapers the tail: 48/64 of the pairs
|
|||
|
|
// overlap at lag 16.
|
|||
|
|
fmt.Printf("acf(0) = %.4f, acf(16) = %.4f\n", acf.FloatAt(0), acf.FloatAt(16))
|
|||
|
|
// Output: acf(0) = 1.0000, acf(16) = 0.7500
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// Differential evolution crosses Rastrigin's minefield of local
|
|||
|
|
// minima to the global one.
|
|||
|
|
func ExampleMinimiseDifferentialEvolution() {
|
|||
|
|
lower, _ := tensor.FromFloats([]float64{-5.12, -5.12}, 2)
|
|||
|
|
upper, _ := tensor.FromFloats([]float64{5.12, 5.12}, 2)
|
|||
|
|
_, fv, err := tensor.MinimiseDifferentialEvolution(func(a *tensor.Array) (float64, error) {
|
|||
|
|
s := 0.0
|
|||
|
|
for i := range 2 {
|
|||
|
|
z := a.FloatAt(i)
|
|||
|
|
s += z*z - 10*math.Cos(2*math.Pi*z) + 10
|
|||
|
|
}
|
|||
|
|
return s, nil
|
|||
|
|
}, lower, upper, tensor.DifferentialEvolutionOptions{Generations: 800})
|
|||
|
|
if err != nil {
|
|||
|
|
log.Fatal(err)
|
|||
|
|
}
|
|||
|
|
fmt.Printf("%.1e\n", fv)
|
|||
|
|
// Output: 0.0e+00
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// Slicing a contiguous selection returns a read-only view sharing the
|
|||
|
|
// storage; an interior range is copied.
|
|||
|
|
func ExampleSlice() {
|
|||
|
|
m, _ := tensor.FromFloats([]float64{
|
|||
|
|
1, 2, 3, 4,
|
|||
|
|
5, 6, 7, 8,
|
|||
|
|
9, 10, 11, 12,
|
|||
|
|
}, 3, 4)
|
|||
|
|
// Whole rows: a view.
|
|||
|
|
rows, err := tensor.Slice(m, 0, 1, 3)
|
|||
|
|
if err != nil {
|
|||
|
|
log.Fatal(err)
|
|||
|
|
}
|
|||
|
|
// An interior column range: a copy.
|
|||
|
|
block, _ := tensor.Slice(rows, 1, 1, 3)
|
|||
|
|
fmt.Println(rows.Shape(), rows.FloatAt(0), block.Shape(), block.FloatAt(0))
|
|||
|
|
// Output: [2 4] 5 [2 2] 6
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// A size-1 dimension replicates over the target, and anything that
|
|||
|
|
// would not is an error rather than a silent broadcast.
|
|||
|
|
func ExampleBroadcastTo() {
|
|||
|
|
col, _ := tensor.FromFloats([]float64{1, 2, 3}, 3, 1)
|
|||
|
|
wide, err := tensor.BroadcastTo(col, 3, 4)
|
|||
|
|
if err != nil {
|
|||
|
|
log.Fatal(err)
|
|||
|
|
}
|
|||
|
|
_, bad := tensor.BroadcastTo(col, 2, 4)
|
|||
|
|
fmt.Println(wide)
|
|||
|
|
fmt.Println("refused:", bad != nil)
|
|||
|
|
// Output: float (3, 4) [1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3]
|
|||
|
|
// refused: true
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// The gamma function, including the negative half line.
|
|||
|
|
func ExampleGamma() {
|
|||
|
|
x, _ := tensor.FromFloats([]float64{-0.5, 0.5, 5}, 3)
|
|||
|
|
g, err := tensor.Gamma(x)
|
|||
|
|
if err != nil {
|
|||
|
|
log.Fatal(err)
|
|||
|
|
}
|
|||
|
|
fmt.Printf("%.7f %.7f %.1f\n", g.FloatAt(0), g.FloatAt(1), g.FloatAt(2))
|
|||
|
|
// Output: -3.5449077 1.7724539 24.0
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// Linear interpolation clamps outside the sampled range; the monotone
|
|||
|
|
// cubic passes through the same knots with a bounded slope.
|
|||
|
|
func ExampleInterpolate() {
|
|||
|
|
xs, _ := tensor.FromFloats([]float64{0, 1, 3}, 3)
|
|||
|
|
ys, _ := tensor.FromFloats([]float64{0, 2, 2.5}, 3)
|
|||
|
|
query, _ := tensor.FromFloats([]float64{0.5, 2, 5}, 3)
|
|||
|
|
lin, err := tensor.Interpolate(xs, ys, query)
|
|||
|
|
if err != nil {
|
|||
|
|
log.Fatal(err)
|
|||
|
|
}
|
|||
|
|
mono, _ := tensor.InterpolateMonotone(xs, ys, query)
|
|||
|
|
// Both stay inside the bracketing samples, and the tail clamps to
|
|||
|
|
// the last knot.
|
|||
|
|
fmt.Printf("linear %.4f %.4f %.4f\n", lin.FloatAt(0), lin.FloatAt(1), lin.FloatAt(2))
|
|||
|
|
fmt.Printf("pchip %.4f %.4f %.4f\n", mono.FloatAt(0), mono.FloatAt(1), mono.FloatAt(2))
|
|||
|
|
// Output: linear 1.0000 2.2500 2.5000
|
|||
|
|
// pchip 1.2621 2.3716 2.5000
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// The Halton sequence stratifies progressively in every base; the
|
|||
|
|
// index-zero origin is dropped, as it carries no information.
|
|||
|
|
func ExampleHaltonPoints() {
|
|||
|
|
pts, err := tensor.HaltonPoints(3, 2, 0)
|
|||
|
|
if err != nil {
|
|||
|
|
log.Fatal(err)
|
|||
|
|
}
|
|||
|
|
fmt.Println(pts.Shape())
|
|||
|
|
for i := range 3 {
|
|||
|
|
fmt.Printf("(%.4f, %.4f) ", pts.FloatAt(2*i), pts.FloatAt(2*i+1))
|
|||
|
|
}
|
|||
|
|
fmt.Println()
|
|||
|
|
// Output: [3 2]
|
|||
|
|
// (0.5000, 0.3333) (0.2500, 0.6667) (0.7500, 0.1111)
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// A fixed seed replays the same stream, which is what makes a random
|
|||
|
|
// workflow testable.
|
|||
|
|
func ExampleNewGenerator() {
|
|||
|
|
u, err := tensor.Floats(tensor.NewGenerator(7), 3)
|
|||
|
|
if err != nil {
|
|||
|
|
log.Fatal(err)
|
|||
|
|
}
|
|||
|
|
again, _ := tensor.Floats(tensor.NewGenerator(7), 3)
|
|||
|
|
fmt.Printf("%.6f %.6f %.6f, equal: %v\n",
|
|||
|
|
u.FloatAt(0), u.FloatAt(1), u.FloatAt(2), tensor.Equal(u, again))
|
|||
|
|
// Output: 0.055360 0.172116 0.717576, equal: true
|
|||
|
|
}
|