tensor/integrate/event_boundary_pins_test.go

// Copyright (c) 2026 Petr Balvín <opensource@petrbalvin.org> (https://petrbalvin.org)
// SPDX-License-Identifier: MIT

package integrate

import (
	"math"
	"strings"
	"testing"

	"sourcedock.dev/petrbalvin/tensor/internal/core"
)

// Regression pins for the event and boundary contracts: a watch landing
// exactly on the final boundary, backward event search, the RK4 and
// Verlet refusals of non-finite states, the grid mesh origin screen and
// the cubature budget's true cost.

// TestEventExactlyOnFinalBoundary pins the hit a watch landing
// exactly on zero at the final accepted boundary produces, which the
// sign walk used to swallow.
func TestEventExactlyOnFinalBoundary(t *testing.T) {
	f := func(_ float64, y *core.Array) (*core.Array, error) {
		return mustFloats(t, []float64{1}, 1), nil
	}
	y0 := mustFloats(t, []float64{0}, 1)
	watch := ODEWatch{
		Function:  func(tt float64, _ *core.Array) (float64, error) { return tt - 1, nil },
		Direction: 1,
	}
	hits, _, err := IntegrateODEEvents(f, 0, 1, y0, []ODEWatch{watch}, ODEOptions{})
	if err != nil {
		t.Fatalf("IntegrateODEEvents: %v", err)
	}
	if len(hits) != 1 || !hits[0].Rising || math.Abs(hits[0].Time-1) > 1e-12 {
		t.Fatalf("hits = %+v, want one rising hit at t = 1", hits)
	}
}

// TestEventsBackward pins the event machinery in the backward
// direction: the watch g = t − 0.5 falls through zero at t = 0.5, and
// the backward run must report that hit with the time refined to the
// integrator's accuracy.
func TestEventsBackward(t *testing.T) {
	f := func(_ float64, y *core.Array) (*core.Array, error) {
		return mustFloats(t, []float64{1}, 1), nil
	}
	y0 := mustFloats(t, []float64{0}, 1)
	watch := ODEWatch{
		Function:  func(tt float64, _ *core.Array) (float64, error) { return tt - 0.5, nil },
		Direction: -1,
	}
	hits, _, err := IntegrateODEEvents(f, 1, 0, y0, []ODEWatch{watch}, ODEOptions{})
	if err != nil {
		t.Fatalf("IntegrateODEEvents backward: %v", err)
	}
	if len(hits) != 1 || hits[0].Rising {
		t.Fatalf("hits = %+v, want one falling hit", hits)
	}
	if math.Abs(hits[0].Time-0.5) > 1e-9 {
		t.Fatalf("hit time = %g, want 0.5", hits[0].Time)
	}
}

// TestRK4AndVerletRefuseNonFinite pins the loud refusals on
// the fixed-step integrators, which published NaN states with nil
// errors before.
func TestRK4AndVerletRefuseNonFinite(t *testing.T) {
	bad := func(_ float64, _ *core.Array) (*core.Array, error) {
		return mustFloats(t, []float64{math.NaN()}, 1), nil
	}
	y0 := mustFloats(t, []float64{0}, 1)
	if _, err := IntegrateRK4(bad, 0, 1, y0, 4); err == nil {
		t.Fatal("IntegrateRK4: expected an error for a NaN derivative")
	}
	accel := func(_ *core.Array) (*core.Array, error) {
		return mustFloats(t, []float64{math.Inf(1)}, 1), nil
	}
	q0 := mustFloats(t, []float64{0}, 1)
	p0 := mustFloats(t, []float64{1}, 1)
	if _, _, err := IntegrateVerlet(accel, 0, 1, q0, p0, 4); err == nil {
		t.Fatal("IntegrateVerlet: expected an error for an Inf acceleration")
	}
}

// TestGridMeshRejectsNonFiniteOrigin pins the origin guard.
func TestGridMeshRejectsNonFiniteOrigin(t *testing.T) {
	if _, err := GridTriangleMesh2D(math.NaN(), 0, 1, 1, 2, 2); err == nil {
		t.Fatal("GridTriangleMesh2D: expected an error for a NaN origin")
	}
	if _, err := GridTriangleMesh2D(0, math.Inf(1), 1, 1, 2, 2); err == nil {
		t.Fatal("GridTriangleMesh2D: expected an error for an Inf origin")
	}
	// An empty triangle table is refused at construction.
	v, _ := core.FromFloats([]float64{0, 0, 1, 0, 0, 1}, 3, 2)
	tri, _ := core.FromInts([]int64{}, 0, 3)
	if _, err := NewTriangleMesh2D(v, tri); err == nil {
		t.Fatal("NewTriangleMesh2D: expected an error for an empty triangle table")
	}
}

// TestCubatureBudgetAccountsTrueCost pins the true bisection
// cost 2·(5^d + 3^d): a budget that admits the root box and exactly
// one bisection must complete, and the dimension guard still refuses
// the twenties under any budget.
func TestCubatureBudgetAccountsTrueCost(t *testing.T) {
	f := func(x []float64) float64 { return x[0] * x[0] }
	lower := []float64{0}
	upper := []float64{2}
	// Root box: 5 + 3 = 8; one bisection: 2·8 = 16. A budget of 24
	// admits the box and one bisection; the old 8^1 + 6^1 = 14
	// accounting let the loop overshoot it by two evaluations.
	if _, err := IntegrateND(f, lower, upper, CubatureOptions{MaxEvals: 24}); err != nil && !strings.Contains(err.Error(), "converge") {
		t.Fatalf("budget 24: err = %v", err)
	}
	lower25 := make([]float64, 25)
	upper25 := make([]float64, 25)
	for i := range upper25 {
		upper25[i] = 1
	}
	one := func([]float64) float64 { return 1 }
	if _, err := IntegrateND(one, lower25, upper25, CubatureOptions{MaxEvals: math.MaxInt}); err == nil || !strings.Contains(err.Error(), "budget") {
		t.Fatalf("d = 25 under a MaxInt budget: err = %v", err)
	}
}