257 lines
8.3 KiB
Go
257 lines
8.3 KiB
Go
|
// Copyright (c) 2026 Petr Balvín <opensource@petrbalvin.org> (https://petrbalvin.org)
|
|||
|
|
// SPDX-License-Identifier: MIT
|
|||
|
|
|
|||
|
|
package signal
|
|||
|
|
|
|||
|
|
import (
|
|||
|
|
"math"
|
|||
|
|
"testing"
|
|||
|
|
|
|||
|
|
"sourcedock.dev/petrbalvin/tensor/internal/core"
|
|||
|
|
)
|
|||
|
|
|
|||
|
|
// TestSumKahanRecoversSmallTerms checks the compensation: small terms
|
|||
|
|
// between huge ones survive, where a plain left-to-right sum loses
|
|||
|
|
// them entirely (the naive sum of this sequence is 0).
|
|||
|
|
func TestSumKahanRecoversSmallTerms(t *testing.T) {
|
|||
|
|
vals := mustFloats(t, []float64{1e16, 1, 2, -1e16})
|
|||
|
|
got, err := SumKahan(vals)
|
|||
|
|
if err != nil {
|
|||
|
|
t.Fatalf("SumKahan: %v", err)
|
|||
|
|
}
|
|||
|
|
if got != 4 {
|
|||
|
|
t.Fatalf("SumKahan = %g, want 4", got)
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// TestGradient1DDifferentiatesLinesExactly: both the central stencil
|
|||
|
|
// and the one-sided endpoints are exact on linear signals.
|
|||
|
|
func TestGradient1DDifferentiatesLinesExactly(t *testing.T) {
|
|||
|
|
const n = 16
|
|||
|
|
const dx = 0.25
|
|||
|
|
y := make([]float64, n)
|
|||
|
|
for i := range n {
|
|||
|
|
y[i] = 3*float64(i)*dx + 7
|
|||
|
|
}
|
|||
|
|
got, err := Gradient1D(mustFloats(t, y, n), dx)
|
|||
|
|
if err != nil {
|
|||
|
|
t.Fatalf("Gradient1D: %v", err)
|
|||
|
|
}
|
|||
|
|
for i := range n {
|
|||
|
|
if math.Abs(got.FloatAt(i)-3) > 1e-12 {
|
|||
|
|
t.Fatalf("gradient[%d] = %.14g, want 3", i, got.FloatAt(i))
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// TestGradient1DKnownDifferentiatesQuadratics: the central interior
|
|||
|
|
// stencil is exact on quadratics; the first-order endpoint stencils
|
|||
|
|
// carry an O(dx) offset.
|
|||
|
|
func TestGradient1DKnownDifferentiatesQuadratics(t *testing.T) {
|
|||
|
|
const n = 16
|
|||
|
|
const dx = 0.25
|
|||
|
|
y := make([]float64, n)
|
|||
|
|
for i := range n {
|
|||
|
|
x := float64(i) * dx
|
|||
|
|
y[i] = x * x
|
|||
|
|
}
|
|||
|
|
got, err := Gradient1D(mustFloats(t, y, n), dx)
|
|||
|
|
if err != nil {
|
|||
|
|
t.Fatalf("Gradient1D: %v", err)
|
|||
|
|
}
|
|||
|
|
for i := 1; i < n-1; i++ {
|
|||
|
|
want := 2 * float64(i) * dx
|
|||
|
|
if math.Abs(got.FloatAt(i)-want) > 1e-12 {
|
|||
|
|
t.Fatalf("gradient[%d] = %.14g, want %.14g", i, got.FloatAt(i), want)
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
if math.Abs(got.FloatAt(0)-dx) > 1e-12 {
|
|||
|
|
t.Fatalf("left endpoint = %.14g, want the one-sided 2x+h = %.14g", got.FloatAt(0), dx)
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// TestGradient1DErrors pins the validation contract.
|
|||
|
|
func TestGradient1DErrors(t *testing.T) {
|
|||
|
|
if _, err := Gradient1D(mustFloats(t, []float64{1}), 0.1); err == nil {
|
|||
|
|
t.Fatal("expected an error for a single-point signal")
|
|||
|
|
}
|
|||
|
|
if _, err := Gradient1D(mustFloats(t, []float64{1, 2, 3}), 0); err == nil {
|
|||
|
|
t.Fatal("expected an error for a zero spacing")
|
|||
|
|
}
|
|||
|
|
rank2, _ := core.FromFloats([]float64{1, 2, 3, 4}, 2, 2)
|
|||
|
|
if _, err := Gradient1D(rank2, 0.1); err == nil {
|
|||
|
|
t.Fatal("expected an error for a rank-2 signal")
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// TestLaplacian1DKnownDifferentiatesSines checks the 3-point stencil
|
|||
|
|
// against −sin on a full period and the boundary copy at both ends.
|
|||
|
|
func TestLaplacian1DKnownDifferentiatesSines(t *testing.T) {
|
|||
|
|
const n = 64
|
|||
|
|
const dx = 2 * math.Pi / n
|
|||
|
|
y := make([]float64, n)
|
|||
|
|
for i := range n {
|
|||
|
|
y[i] = math.Sin(float64(i) * dx)
|
|||
|
|
}
|
|||
|
|
got, err := Laplacian(mustFloats(t, y, n), dx)
|
|||
|
|
if err != nil {
|
|||
|
|
t.Fatalf("Laplacian: %v", err)
|
|||
|
|
}
|
|||
|
|
// The interior matches −sin; the boundary copies the neighbour, so
|
|||
|
|
// it is checked separately below.
|
|||
|
|
for i := 1; i < n-1; i++ {
|
|||
|
|
want := -math.Sin(float64(i) * dx)
|
|||
|
|
if math.Abs(got.FloatAt(i)-want) > 1e-3 {
|
|||
|
|
t.Fatalf("laplacian[%d] = %.6g, want %.6g", i, got.FloatAt(i), want)
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
// The boundaries copy their nearest interior value.
|
|||
|
|
if got.FloatAt(0) != got.FloatAt(1) || got.FloatAt(n-1) != got.FloatAt(n-2) {
|
|||
|
|
t.Fatal("the 1-D boundaries did not copy their interior values")
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// TestLaplacian2DKnownDifferentiatesSines checks the 5-point stencil
|
|||
|
|
// against −2·sin x·cos y and the boundary copy on every edge.
|
|||
|
|
func TestLaplacian2DKnownDifferentiatesSines(t *testing.T) {
|
|||
|
|
const n = 32
|
|||
|
|
const d = 2 * math.Pi / n
|
|||
|
|
grid := make([]float64, n*n)
|
|||
|
|
for r := range n {
|
|||
|
|
for c := range n {
|
|||
|
|
grid[r*n+c] = math.Sin(float64(c)*d) * math.Cos(float64(r)*d)
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
got, err := Laplacian(mustFloats(t, grid, n, n), d, d)
|
|||
|
|
if err != nil {
|
|||
|
|
t.Fatalf("Laplacian: %v", err)
|
|||
|
|
}
|
|||
|
|
// The interior matches −2·sin x·cos y to the stencil's own
|
|||
|
|
// truncation error 2h²/12 ≈ 0.0064 at h = 2π/32; the boundary
|
|||
|
|
// copies its nearest interior value, checked separately below.
|
|||
|
|
for r := 1; r < n-1; r++ {
|
|||
|
|
for c := 1; c < n-1; c++ {
|
|||
|
|
want := -2 * math.Sin(float64(c)*d) * math.Cos(float64(r)*d)
|
|||
|
|
if math.Abs(got.FloatAt(r*n+c)-want) > 8e-3 {
|
|||
|
|
t.Fatalf("laplacian[%d,%d] = %.6g, want %.6g", r, c, got.FloatAt(r*n+c), want)
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
// Every boundary point must equal its nearest interior neighbour.
|
|||
|
|
for i := range n {
|
|||
|
|
if got.FloatAt(i) != got.FloatAt(n+i) {
|
|||
|
|
t.Fatalf("top row %d did not copy the row below", i)
|
|||
|
|
}
|
|||
|
|
if got.FloatAt((n-1)*n+i) != got.FloatAt((n-2)*n+i) {
|
|||
|
|
t.Fatalf("bottom row %d did not copy the row above", i)
|
|||
|
|
}
|
|||
|
|
if got.FloatAt(i*n) != got.FloatAt(i*n+1) {
|
|||
|
|
t.Fatalf("left column %d did not copy its interior neighbour", i)
|
|||
|
|
}
|
|||
|
|
if got.FloatAt(i*n+n-1) != got.FloatAt(i*n+n-2) {
|
|||
|
|
t.Fatalf("right column %d did not copy its interior neighbour", i)
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// TestLaplacian3DKnownDifferentiatesSines checks the 7-point stencil
|
|||
|
|
// against −3·sin x·sin y·sin z, and pins the boundary behaviour: all
|
|||
|
|
// six faces copy their nearest interior plane (the doc's promise).
|
|||
|
|
func TestLaplacian3DKnownDifferentiatesSines(t *testing.T) {
|
|||
|
|
const n = 16
|
|||
|
|
const d = 2 * math.Pi / n
|
|||
|
|
at := func(z, y, x int) float64 {
|
|||
|
|
return math.Sin(float64(x)*d) * math.Sin(float64(y)*d) * math.Sin(float64(z)*d)
|
|||
|
|
}
|
|||
|
|
grid := make([]float64, n*n*n)
|
|||
|
|
for k := range n {
|
|||
|
|
for j := range n {
|
|||
|
|
for i := range n {
|
|||
|
|
grid[(k*n+j)*n+i] = at(k, j, i)
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
got, err := Laplacian(mustFloats(t, grid, n, n, n), d, d, d)
|
|||
|
|
if err != nil {
|
|||
|
|
t.Fatalf("Laplacian: %v", err)
|
|||
|
|
}
|
|||
|
|
// Interior accuracy first; the six boundary faces are copies, so
|
|||
|
|
// they are excluded from the truncation-error bound and checked
|
|||
|
|
// against their interior neighbours below.
|
|||
|
|
worst := 0.0
|
|||
|
|
for k := 1; k < n-1; k++ {
|
|||
|
|
for j := 1; j < n-1; j++ {
|
|||
|
|
for i := 1; i < n-1; i++ {
|
|||
|
|
want := -3 * at(k, j, i)
|
|||
|
|
d := math.Abs(got.FloatAt((k*n+j)*n+i) - want)
|
|||
|
|
if d > worst {
|
|||
|
|
worst = d
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
if worst > 5e-2 {
|
|||
|
|
// The bound is the stencil's own truncation error 3h²/12 ≈ 0.039
|
|||
|
|
// at h = 2π/16; a wrong stencil is orders of magnitude worse.
|
|||
|
|
t.Fatalf("7-point stencil error %.3g, want under 5e-2", worst)
|
|||
|
|
}
|
|||
|
|
// All six boundary faces copy their nearest interior plane.
|
|||
|
|
for j := range n {
|
|||
|
|
for i := range n {
|
|||
|
|
if got.FloatAt(j*n+i) != got.FloatAt((n+j)*n+i) {
|
|||
|
|
t.Fatalf("front plane (%d,%d) did not copy the plane behind it", j, i)
|
|||
|
|
}
|
|||
|
|
if got.FloatAt(((n-1)*n+j)*n+i) != got.FloatAt(((n-2)*n+j)*n+i) {
|
|||
|
|
t.Fatalf("back plane (%d,%d) did not copy the plane before it", j, i)
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
for k := range n {
|
|||
|
|
for i := range n {
|
|||
|
|
if got.FloatAt((k*n)*n+i) != got.FloatAt((k*n+1)*n+i) {
|
|||
|
|
t.Fatalf("top face (%d,%d) did not copy the row below", k, i)
|
|||
|
|
}
|
|||
|
|
if got.FloatAt(((k+1)*n-1)*n+i) != got.FloatAt(((k+1)*n-2)*n+i) {
|
|||
|
|
t.Fatalf("bottom face (%d,%d) did not copy the row above", k, i)
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
for j := range n {
|
|||
|
|
if got.FloatAt((k*n+j)*n) != got.FloatAt((k*n+j)*n+1) {
|
|||
|
|
t.Fatalf("left face (%d,%d) did not copy its interior neighbour", k, j)
|
|||
|
|
}
|
|||
|
|
if got.FloatAt((k*n+j)*n+n-1) != got.FloatAt((k*n+j)*n+n-2) {
|
|||
|
|
t.Fatalf("right face (%d,%d) did not copy its interior neighbour", k, j)
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// TestLaplacianDegenerateErrors pins the extent contract: every rank-1
|
|||
|
|
// signal needs at least 2 points and every axis of a rank-2 or rank-3
|
|||
|
|
// grid at least 3, so the stencil and the boundary copies have
|
|||
|
|
// somewhere to stand.
|
|||
|
|
func TestLaplacianDegenerateErrors(t *testing.T) {
|
|||
|
|
if _, err := Laplacian(mustFloats(t, []float64{1}), 0.1); err == nil {
|
|||
|
|
t.Fatal("expected an error for a single-point 1-D signal")
|
|||
|
|
}
|
|||
|
|
tiny2D, _ := core.FromFloats(make([]float64, 16), 2, 8)
|
|||
|
|
if _, err := Laplacian(tiny2D, 0.1, 0.1); err == nil {
|
|||
|
|
t.Fatal("expected an error for a 2-D grid with an extent of 2")
|
|||
|
|
}
|
|||
|
|
tiny3D, _ := core.FromFloats(make([]float64, 32), 2, 4, 4)
|
|||
|
|
if _, err := Laplacian(tiny3D, 0.1, 0.1, 0.1); err == nil {
|
|||
|
|
t.Fatal("expected an error for a 3-D grid with an extent of 2")
|
|||
|
|
}
|
|||
|
|
twoD, _ := core.FromFloats(make([]float64, 9), 3, 3)
|
|||
|
|
if _, err := Laplacian(twoD, 0.1); err == nil {
|
|||
|
|
t.Fatal("expected an error for a missing spacing")
|
|||
|
|
}
|
|||
|
|
if _, err := Laplacian(mustFloats(t, []float64{1, 2, 3}), 0); err == nil {
|
|||
|
|
t.Fatal("expected an error for a zero spacing")
|
|||
|
|
}
|
|||
|
|
rank4, _ := core.FromFloats(make([]float64, 16), 2, 2, 2, 2)
|
|||
|
|
if _, err := Laplacian(rank4, 0.1, 0.1, 0.1, 0.1); err == nil {
|
|||
|
|
t.Fatal("expected an error for a rank-4 array")
|
|||
|
|
}
|
|||
|
|
}
|