239 lines
7.0 KiB
Go
239 lines
7.0 KiB
Go
// Copyright (c) 2026 Petr Balvín <opensource@petrbalvin.org> (https://petrbalvin.org)
|
||
// SPDX-License-Identifier: MIT
|
||
|
||
package core
|
||
|
||
import (
|
||
"math"
|
||
"strconv"
|
||
"strings"
|
||
"testing"
|
||
)
|
||
|
||
// TestHaltonRadicalInverse pins the van der Corput construction against
|
||
// the digit-reflection definition.
|
||
func TestHaltonRadicalInverse(t *testing.T) {
|
||
cases := []struct {
|
||
i, b int
|
||
want float64
|
||
}{
|
||
{1, 2, 0.5},
|
||
{2, 2, 0.25},
|
||
{3, 2, 0.75},
|
||
{4, 2, 0.125},
|
||
{5, 2, 0.625},
|
||
{1, 3, 1.0 / 3},
|
||
{2, 3, 2.0 / 3},
|
||
{3, 3, 1.0 / 9},
|
||
{7, 3, 5.0 / 9}, // 21₃ reflected is 12₃ = 1/3 + 2/9
|
||
{5, 3, 7.0 / 9}, // 12₃ reflected is 21₃ = 2/3 + 1/9
|
||
{1, 5, 0.2},
|
||
{4, 5, 0.8},
|
||
}
|
||
for _, c := range cases {
|
||
if got := radicalInverse(c.i, c.b); math.Abs(got-c.want) > 1e-15 {
|
||
t.Fatalf("radicalInverse(%d, %d) = %g, want %g", c.i, c.b, got, c.want)
|
||
}
|
||
}
|
||
}
|
||
|
||
// TestHaltonStratification pins the defining property: the first b^k
|
||
// points integrate a smooth function far better than random sampling
|
||
// would, and ∫x dx through the mean sits at 1/2 to high accuracy.
|
||
func TestHaltonStratification(t *testing.T) {
|
||
const n = 4096
|
||
pts, err := HaltonPoints(n, 2, 0)
|
||
if err != nil {
|
||
t.Fatalf("HaltonPoints: %v", err)
|
||
}
|
||
mean := 0.0
|
||
for i := range n {
|
||
mean += pts.FloatAt(i * 2)
|
||
}
|
||
mean /= float64(n)
|
||
// The base-2 marginal is a permutation of {k/8192}: its mean sits
|
||
// at 1/2 to within one point's contribution.
|
||
if math.Abs(mean-0.5) > 1.5/float64(n) {
|
||
t.Fatalf("mean of base-2 marginal = %.8f, want 1/2 within one point", mean)
|
||
}
|
||
}
|
||
|
||
// TestHaltonQuadratureAccuracy pins the Monte Carlo payoff: a Halton
|
||
// estimate of a smooth 2-D integral beats the random-sampling error
|
||
// scale by orders of magnitude.
|
||
func TestHaltonQuadratureAccuracy(t *testing.T) {
|
||
const n = 1024
|
||
pts, err := HaltonPoints(n, 2, 100)
|
||
if err != nil {
|
||
t.Fatalf("HaltonPoints: %v", err)
|
||
}
|
||
est := 0.0
|
||
for i := range n {
|
||
x := pts.FloatAt(i * 2)
|
||
y := pts.FloatAt(i*2 + 1)
|
||
est += math.Exp(-x-y*y) / float64(n)
|
||
}
|
||
// Reference: ∫₀¹∫₀¹ e^{-x-y²} = (1−e^{-1})·sqrt(π)/2·erf(1). The
|
||
// Halton error decays like O(log n / n); a thousand points sit
|
||
// around 1e-4, where random sampling would scatter at 1/sqrt(n)
|
||
// ~ 3e-2.
|
||
want := (1 - math.Exp(-1)) * math.Sqrt(math.Pi) / 2 * math.Erf(1)
|
||
if math.Abs(est-want) > 5e-4 {
|
||
t.Fatalf("Halton integral = %.8f, exact %.8f", est, want)
|
||
}
|
||
}
|
||
|
||
// TestHaltonErrors pins the input gates.
|
||
func TestHaltonErrors(t *testing.T) {
|
||
if _, err := HaltonPoints(10, 0, 0); err == nil {
|
||
t.Error("dim 0 accepted")
|
||
}
|
||
if _, err := HaltonPoints(10, 40, 0); err == nil {
|
||
t.Error("dim 40 accepted")
|
||
}
|
||
if _, err := HaltonPoints(10, 2, -1); err == nil {
|
||
t.Error("negative skip accepted")
|
||
}
|
||
// Zero points is a valid empty request.
|
||
empty, err := HaltonPoints(0, 2, 0)
|
||
if err != nil || empty.Len() != 0 {
|
||
t.Fatalf("n = 0 must give an empty array, got err %v", err)
|
||
}
|
||
}
|
||
|
||
// TestSobolCanonicalPoints pins the two dimensions every Sobol table
|
||
// agrees on: dimension 1 is Gray-coded van der Corput and the first
|
||
// four 2-D points are the canonical square-covering quartet.
|
||
func TestSobolCanonicalPoints(t *testing.T) {
|
||
want1 := []float64{0.5, 0.75, 0.25, 0.375, 0.875, 0.625, 0.125, 0.1875}
|
||
p1, err := SobolPoints(len(want1), 1, 0)
|
||
if err != nil {
|
||
t.Fatalf("SobolPoints: %v", err)
|
||
}
|
||
for i, want := range want1 {
|
||
if got := p1.FloatAt(i); got != want {
|
||
t.Fatalf("dim 1 point %d = %g, want %g", i+1, got, want)
|
||
}
|
||
}
|
||
|
||
pts, err := SobolPoints(4, 2, 0)
|
||
if err != nil {
|
||
t.Fatalf("SobolPoints: %v", err)
|
||
}
|
||
want2 := [4][2]float64{{0.5, 0.5}, {0.75, 0.25}, {0.25, 0.75}, {0.375, 0.375}}
|
||
for p, want := range want2 {
|
||
for d := range 2 {
|
||
if got := pts.FloatAt(p*2 + d); got != want[d] {
|
||
t.Fatalf("point %d dim %d = %g, want %g", p+1, d+1, got, want[d])
|
||
}
|
||
}
|
||
}
|
||
}
|
||
|
||
// TestSobolStratification pins the defining digital property across
|
||
// the whole table: in every dimension the 2^m points X_0..X_{2^m-1}
|
||
// (the origin included, since it carries the zero cell) hit each of
|
||
// the 2^m equal subintervals exactly once, which holds only if the
|
||
// direction numbers are linearly independent over GF(2).
|
||
func TestSobolStratification(t *testing.T) {
|
||
for _, b := range []int{4, 8} {
|
||
count := 1 << b
|
||
for dim := 1; dim <= len(sobolTable); dim++ {
|
||
// X_1..X_{count-1}; X_0 is the origin and fills cell 0.
|
||
pts, err := SobolPoints(count-1, dim, 0)
|
||
if err != nil {
|
||
t.Fatalf("SobolPoints dim %d: %v", dim, err)
|
||
}
|
||
for d := range dim {
|
||
seen := make([]bool, count)
|
||
seen[0] = true
|
||
for p := range count - 1 {
|
||
cell := int(pts.FloatAt(p*dim+d) * float64(count))
|
||
if cell == count {
|
||
cell = count - 1 // a value of 1.0 would be out of range
|
||
}
|
||
if seen[cell] {
|
||
t.Fatalf("dim %d, coordinate %d: subinterval %d hit twice in the first %d points", dim, d, cell, count)
|
||
}
|
||
seen[cell] = true
|
||
}
|
||
}
|
||
}
|
||
}
|
||
}
|
||
|
||
// TestSobolSkip: skipping must be equivalent to generating more points
|
||
// and discarding the leading ones, for every coordinate.
|
||
func TestSobolSkip(t *testing.T) {
|
||
full, err := SobolPoints(12, 3, 0)
|
||
if err != nil {
|
||
t.Fatalf("SobolPoints: %v", err)
|
||
}
|
||
tail, err := SobolPoints(4, 3, 8)
|
||
if err != nil {
|
||
t.Fatalf("SobolPoints: %v", err)
|
||
}
|
||
for p := range 4 {
|
||
for d := range 3 {
|
||
got := tail.FloatAt(p*3 + d)
|
||
want := full.FloatAt((p+8)*3 + d)
|
||
if got != want {
|
||
t.Fatalf("skipped point %d dim %d = %g, want %g", p, d+1, got, want)
|
||
}
|
||
}
|
||
}
|
||
}
|
||
|
||
// TestSobolQuadratureAccuracy: the same integral TestHaltonQuadrature
|
||
// pins, at the same sample count and bound. Sobol carries no
|
||
// early-point correlation to skip away, so the sample runs from the
|
||
// start, where the digital structure is densest.
|
||
func TestSobolQuadratureAccuracy(t *testing.T) {
|
||
const n = 1024
|
||
pts, err := SobolPoints(n, 2, 0)
|
||
if err != nil {
|
||
t.Fatalf("SobolPoints: %v", err)
|
||
}
|
||
est := 0.0
|
||
for i := range n {
|
||
x := pts.FloatAt(i * 2)
|
||
y := pts.FloatAt(i*2 + 1)
|
||
est += math.Exp(-x-y*y) / float64(n)
|
||
}
|
||
want := (1 - math.Exp(-1)) * math.Sqrt(math.Pi) / 2 * math.Erf(1)
|
||
if math.Abs(est-want) > 5e-4 {
|
||
t.Fatalf("Sobol integral = %.8f, exact %.8f", est, want)
|
||
}
|
||
}
|
||
|
||
// TestSobolErrors pins the input gates.
|
||
func TestSobolErrors(t *testing.T) {
|
||
if _, err := SobolPoints(10, 0, 0); err == nil {
|
||
t.Error("dim 0 accepted")
|
||
}
|
||
if _, err := SobolPoints(10, 41, 0); err == nil {
|
||
t.Error("dim 41 accepted")
|
||
}
|
||
if _, err := SobolPoints(-1, 2, 0); err == nil {
|
||
t.Error("negative n accepted")
|
||
}
|
||
if _, err := SobolPoints(10, 2, -1); err == nil {
|
||
t.Error("negative skip accepted")
|
||
}
|
||
empty, err := SobolPoints(0, 2, 0)
|
||
if err != nil || empty.Len() != 0 {
|
||
t.Fatalf("n = 0 must give an empty array, got err %v", err)
|
||
}
|
||
// The period guard needs indices a 32-bit int cannot express.
|
||
if strconv.IntSize >= 64 {
|
||
if _, err := SobolPoints(2, 1, (1<<32)-1); err == nil ||
|
||
!strings.Contains(err.Error(), "sequence period") {
|
||
t.Errorf("period wrap accepted: %v", err)
|
||
}
|
||
// The last legal index is still fine.
|
||
if _, err := SobolPoints(1, 1, (1<<32)-2); err != nil {
|
||
t.Errorf("last legal index refused: %v", err)
|
||
}
|
||
}
|
||
}
|