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)
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
}
|