610 lines
20 KiB
Go
610 lines
20 KiB
Go
|
// Copyright (c) 2026 Petr Balvín <opensource@petrbalvin.org> (https://petrbalvin.org)
|
|||
|
|
// SPDX-License-Identifier: MIT
|
|||
|
|
|
|||
|
|
package optim
|
|||
|
|
|
|||
|
|
import (
|
|||
|
|
"math"
|
|||
|
|
"math/big"
|
|||
|
|
"strings"
|
|||
|
|
"testing"
|
|||
|
|
|
|||
|
|
"sourcedock.dev/petrbalvin/tensor/internal/base"
|
|||
|
|
"sourcedock.dev/petrbalvin/tensor/internal/core"
|
|||
|
|
)
|
|||
|
|
|
|||
|
|
// linearPair returns the determined linear pair r(p) = (p₀ + 2p₁ − 3,
|
|||
|
|
// 2p₀ + p₁ − 4) with its constant Jacobian, the model whose
|
|||
|
|
// Gauss-Newton step is exact.
|
|||
|
|
func linearPair() (func(*core.Array) (*core.Array, error), func(*core.Array) (*core.Array, error)) {
|
|||
|
|
residual := func(p *core.Array) (*core.Array, error) {
|
|||
|
|
return core.FromFloats([]float64{
|
|||
|
|
p.FloatAt(0) + 2*p.FloatAt(1) - 3,
|
|||
|
|
2*p.FloatAt(0) + p.FloatAt(1) - 4,
|
|||
|
|
}, 2)
|
|||
|
|
}
|
|||
|
|
jacobian := func(*core.Array) (*core.Array, error) {
|
|||
|
|
return core.FromFloats([]float64{1, 2, 2, 1}, 2, 2)
|
|||
|
|
}
|
|||
|
|
return residual, jacobian
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// TestLevenbergMarquardtFitConverged pins the Fit surface on the model
|
|||
|
|
// the plain entry point already covers: same point, same chi2, and a
|
|||
|
|
// status that says converged.
|
|||
|
|
func TestLevenbergMarquardtFitConverged(t *testing.T) {
|
|||
|
|
residual, jacobian := linearPair()
|
|||
|
|
p0 := mustFloats(t, []float64{0, 0}, 2)
|
|||
|
|
res, err := LevenbergMarquardtFit(residual, p0, LMOptions{Jacobian: jacobian, Tolerance: 1e-14})
|
|||
|
|
if err != nil {
|
|||
|
|
t.Fatalf("LevenbergMarquardtFit: %v", err)
|
|||
|
|
}
|
|||
|
|
if res.Status != FitConverged {
|
|||
|
|
t.Fatalf("status = %d, want FitConverged", res.Status)
|
|||
|
|
}
|
|||
|
|
if math.Abs(res.Parameters.FloatAt(0)-5.0/3) > 1e-8 || math.Abs(res.Parameters.FloatAt(1)-2.0/3) > 1e-8 {
|
|||
|
|
t.Fatalf("p = (%.12g, %.12g), want (5/3, 2/3)", res.Parameters.FloatAt(0), res.Parameters.FloatAt(1))
|
|||
|
|
}
|
|||
|
|
if res.Chi2 > 1e-16 {
|
|||
|
|
t.Fatalf("chi2 = %v, want 0", res.Chi2)
|
|||
|
|
}
|
|||
|
|
if res.Covariance != nil {
|
|||
|
|
t.Fatal("a fit that did not ask for a covariance returned one")
|
|||
|
|
}
|
|||
|
|
legacyP, legacyChi2, err := LevenbergMarquardt(residual, p0, LMOptions{Jacobian: jacobian, Tolerance: 1e-14})
|
|||
|
|
if err != nil {
|
|||
|
|
t.Fatalf("LevenbergMarquardt: %v", err)
|
|||
|
|
}
|
|||
|
|
if legacyP.FloatAt(0) != res.Parameters.FloatAt(0) || legacyP.FloatAt(1) != res.Parameters.FloatAt(1) {
|
|||
|
|
t.Fatal("the two entry points returned different parameters")
|
|||
|
|
}
|
|||
|
|
if legacyChi2 != res.Chi2 {
|
|||
|
|
t.Fatalf("the two entry points reported chi2 %v and %v", legacyChi2, res.Chi2)
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// TestLevenbergMarquardtGradientExit pins the gradient test: with the
|
|||
|
|
// chi2 tolerance set far below anything the fit can reach, the only
|
|||
|
|
// exit that can fire after the first accepted step is ‖Jᵀr‖∞. The
|
|||
|
|
// residual call count proves it: the start and the one trial make
|
|||
|
|
// two, and a third call would mean the run walked into a second
|
|||
|
|
// solve, which the gradient test forbids.
|
|||
|
|
func TestLevenbergMarquardtGradientExit(t *testing.T) {
|
|||
|
|
residual, jacobian := linearPair()
|
|||
|
|
calls := 0
|
|||
|
|
counted := func(p *core.Array) (*core.Array, error) {
|
|||
|
|
calls++
|
|||
|
|
return residual(p)
|
|||
|
|
}
|
|||
|
|
p0 := mustFloats(t, []float64{0, 0}, 2)
|
|||
|
|
// Lambda 1e-300 leaves the damping factor 1+λ bit-identically 1,
|
|||
|
|
// so the first step is the exact Gauss-Newton one and the run
|
|||
|
|
// reaches the solution in one acceptance.
|
|||
|
|
res, err := LevenbergMarquardtFit(counted, p0, LMOptions{
|
|||
|
|
Jacobian: jacobian,
|
|||
|
|
Tolerance: 1e-30,
|
|||
|
|
GradTol: 1e-7,
|
|||
|
|
Lambda: 1e-300,
|
|||
|
|
})
|
|||
|
|
if err != nil {
|
|||
|
|
t.Fatalf("LevenbergMarquardtFit: %v", err)
|
|||
|
|
}
|
|||
|
|
if res.Status != FitConverged {
|
|||
|
|
t.Fatalf("status = %d, want FitConverged", res.Status)
|
|||
|
|
}
|
|||
|
|
if calls != 2 {
|
|||
|
|
t.Fatalf("the residual was evaluated %d times, want 2 (the gradient test must stop the run before its second solve)", calls)
|
|||
|
|
}
|
|||
|
|
if math.Abs(res.Parameters.FloatAt(0)-5.0/3) > 1e-8 || math.Abs(res.Parameters.FloatAt(1)-2.0/3) > 1e-8 {
|
|||
|
|
t.Fatalf("p = (%.12g, %.12g), want (5/3, 2/3)", res.Parameters.FloatAt(0), res.Parameters.FloatAt(1))
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// TestLevenbergMarquardtStepExit pins the step test: a linear
|
|||
|
|
// residual whose first accepted step is negligible against the
|
|||
|
|
// point's own scale converges on that step, with the chi2 tolerance
|
|||
|
|
// far too tight to claim the exit.
|
|||
|
|
func TestLevenbergMarquardtStepExit(t *testing.T) {
|
|||
|
|
residual := func(p *core.Array) (*core.Array, error) {
|
|||
|
|
return core.FromFloats([]float64{10 * (p.FloatAt(0) - 10)}, 1)
|
|||
|
|
}
|
|||
|
|
jacobian := func(*core.Array) (*core.Array, error) {
|
|||
|
|
return core.FromFloats([]float64{10}, 1, 1)
|
|||
|
|
}
|
|||
|
|
counted := func(count *int) func(*core.Array) (*core.Array, error) {
|
|||
|
|
return func(p *core.Array) (*core.Array, error) {
|
|||
|
|
*count++
|
|||
|
|
return residual(p)
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
stepCalls := 0
|
|||
|
|
res, err := LevenbergMarquardtFit(counted(&stepCalls), mustFloats(t, []float64{10.001}, 1), LMOptions{
|
|||
|
|
Jacobian: jacobian,
|
|||
|
|
Tolerance: 1e-15,
|
|||
|
|
StepTol: 0.5,
|
|||
|
|
})
|
|||
|
|
if err != nil {
|
|||
|
|
t.Fatalf("LevenbergMarquardtFit: %v", err)
|
|||
|
|
}
|
|||
|
|
if res.Status != FitConverged {
|
|||
|
|
t.Fatalf("status = %d, want FitConverged", res.Status)
|
|||
|
|
}
|
|||
|
|
if stepCalls != 2 {
|
|||
|
|
t.Fatalf("the step-exit run evaluated the residual %d times, want 2", stepCalls)
|
|||
|
|
}
|
|||
|
|
if math.Abs(res.Parameters.FloatAt(0)-10) > 1e-4 {
|
|||
|
|
t.Fatalf("p = %.12g, want 10", res.Parameters.FloatAt(0))
|
|||
|
|
}
|
|||
|
|
// Without the step test the same run needs the chi2 tolerance, and
|
|||
|
|
// at 1e-6 that fires one acceptance later: three evaluations.
|
|||
|
|
wideCalls := 0
|
|||
|
|
wide, err := LevenbergMarquardtFit(counted(&wideCalls), mustFloats(t, []float64{10.001}, 1), LMOptions{
|
|||
|
|
Jacobian: jacobian,
|
|||
|
|
Tolerance: 1e-6,
|
|||
|
|
})
|
|||
|
|
if err != nil {
|
|||
|
|
t.Fatalf("LevenbergMarquardtFit without StepTol: %v", err)
|
|||
|
|
}
|
|||
|
|
if wide.Status != FitConverged {
|
|||
|
|
t.Fatalf("status = %d, want FitConverged", wide.Status)
|
|||
|
|
}
|
|||
|
|
if wideCalls != 3 {
|
|||
|
|
t.Fatalf("the control run evaluated the residual %d times, want 3", wideCalls)
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// TestLevenbergMarquardtFitStalled pins the stalled contract: a
|
|||
|
|
// parameter the residual never sees makes the normal equations
|
|||
|
|
// singular on the first iteration, and the fit reports the start
|
|||
|
|
// point back with FitStalled instead of an error. The legacy entry
|
|||
|
|
// point keeps its historical error for the same run.
|
|||
|
|
func TestLevenbergMarquardtFitStalled(t *testing.T) {
|
|||
|
|
residual := func(p *core.Array) (*core.Array, error) {
|
|||
|
|
return core.FromFloats([]float64{p.FloatAt(1) - 1, p.FloatAt(1) - 1}, 2)
|
|||
|
|
}
|
|||
|
|
jacobian := func(*core.Array) (*core.Array, error) {
|
|||
|
|
return core.FromFloats([]float64{0, 1, 0, 1}, 2, 2)
|
|||
|
|
}
|
|||
|
|
p0 := mustFloats(t, []float64{5, 5}, 2)
|
|||
|
|
res, err := LevenbergMarquardtFit(residual, p0, LMOptions{Jacobian: jacobian})
|
|||
|
|
if err != nil {
|
|||
|
|
t.Fatalf("LevenbergMarquardtFit: %v", err)
|
|||
|
|
}
|
|||
|
|
if res.Status != FitStalled {
|
|||
|
|
t.Fatalf("status = %d, want FitStalled", res.Status)
|
|||
|
|
}
|
|||
|
|
if res.Parameters.FloatAt(0) != 5 || res.Parameters.FloatAt(1) != 5 {
|
|||
|
|
t.Fatalf("p = (%.12g, %.12g), want the untouched start", res.Parameters.FloatAt(0), res.Parameters.FloatAt(1))
|
|||
|
|
}
|
|||
|
|
if res.Chi2 != 32 {
|
|||
|
|
t.Fatalf("chi2 = %v, want 32", res.Chi2)
|
|||
|
|
}
|
|||
|
|
if _, _, err := LevenbergMarquardt(residual, p0, LMOptions{Jacobian: jacobian}); err == nil {
|
|||
|
|
t.Fatal("the legacy entry point must keep reporting the singular solve as an error")
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// TestLevenbergMarquardtFitBudget pins the budget report: a fit given
|
|||
|
|
// two iterations of a problem that needs more returns its last point
|
|||
|
|
// with FitBudget and no error, and the legacy entry point refuses the
|
|||
|
|
// same run unless AllowBudgetExit is set.
|
|||
|
|
func TestLevenbergMarquardtFitBudget(t *testing.T) {
|
|||
|
|
xData := []float64{0, 1, 2, 3, 4, 5, 6, 7, 8, 9}
|
|||
|
|
residual := func(p *core.Array) (*core.Array, error) {
|
|||
|
|
a, b, c := p.FloatAt(0), p.FloatAt(1), p.FloatAt(2)
|
|||
|
|
out := core.New(core.Float, len(xData))
|
|||
|
|
for i := range len(xData) {
|
|||
|
|
out.RawFloats()[i] = 3*math.Exp(-0.5*float64(xData[i])) + 0.5 -
|
|||
|
|
(a*math.Exp(-b*float64(xData[i])) + c)
|
|||
|
|
}
|
|||
|
|
return out, nil
|
|||
|
|
}
|
|||
|
|
p0 := mustFloats(t, []float64{20, 5, 20}, 3)
|
|||
|
|
res, err := LevenbergMarquardtFit(residual, p0, LMOptions{MaxIterations: 2})
|
|||
|
|
if err != nil {
|
|||
|
|
t.Fatalf("LevenbergMarquardtFit: %v", err)
|
|||
|
|
}
|
|||
|
|
if res.Status != FitBudget {
|
|||
|
|
t.Fatalf("status = %d, want FitBudget", res.Status)
|
|||
|
|
}
|
|||
|
|
if res.Parameters == nil {
|
|||
|
|
t.Fatal("a budget stop must still report its last point")
|
|||
|
|
}
|
|||
|
|
if _, _, err := LevenbergMarquardt(residual, p0, LMOptions{MaxIterations: 2}); err == nil {
|
|||
|
|
t.Fatal("the legacy entry point must refuse a budget stop without AllowBudgetExit")
|
|||
|
|
}
|
|||
|
|
if _, _, err := LevenbergMarquardt(residual, p0, LMOptions{MaxIterations: 2, AllowBudgetExit: true}); err != nil {
|
|||
|
|
t.Fatalf("LevenbergMarquardt with AllowBudgetExit: %v", err)
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// weightedModel builds the three-observation linear model r(p) = y −
|
|||
|
|
// Xp with its Jacobian, small enough for the exact rational referent.
|
|||
|
|
func weightedModel() (func(*core.Array) (*core.Array, error), func(*core.Array) (*core.Array, error)) {
|
|||
|
|
x := [][]float64{{1, 0}, {0, 1}, {1, 1}}
|
|||
|
|
y := []float64{1.2, 0.7, 2.4}
|
|||
|
|
residual := func(p *core.Array) (*core.Array, error) {
|
|||
|
|
out := core.New(core.Float, 3)
|
|||
|
|
for i := range 3 {
|
|||
|
|
out.RawFloats()[i] = y[i] - (x[i][0]*p.FloatAt(0) + x[i][1]*p.FloatAt(1))
|
|||
|
|
}
|
|||
|
|
return out, nil
|
|||
|
|
}
|
|||
|
|
jacobian := func(*core.Array) (*core.Array, error) {
|
|||
|
|
return core.FromFloats([]float64{-1, 0, 0, -1, -1, -1}, 3, 2)
|
|||
|
|
}
|
|||
|
|
return residual, jacobian
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// ratSolve solves a x = b exactly over the rationals by elimination
|
|||
|
|
// with nonzero pivots; a singular system returns nil. The inputs are
|
|||
|
|
// copied: big.Rat values are pointers, and the elimination mutates
|
|||
|
|
// every entry it touches.
|
|||
|
|
func ratSolve(a, b [][]*big.Rat) [][]*big.Rat {
|
|||
|
|
n := len(a)
|
|||
|
|
m := make([][]*big.Rat, n)
|
|||
|
|
for i := range n {
|
|||
|
|
m[i] = make([]*big.Rat, 0, len(a[i])+len(b[i]))
|
|||
|
|
for _, v := range a[i] {
|
|||
|
|
m[i] = append(m[i], new(big.Rat).Set(v))
|
|||
|
|
}
|
|||
|
|
for _, v := range b[i] {
|
|||
|
|
m[i] = append(m[i], new(big.Rat).Set(v))
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
zero := new(big.Rat)
|
|||
|
|
for col := range n {
|
|||
|
|
piv := -1
|
|||
|
|
for row := col; row < n; row++ {
|
|||
|
|
if m[row][col].Cmp(zero) != 0 {
|
|||
|
|
piv = row
|
|||
|
|
break
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
if piv < 0 {
|
|||
|
|
return nil
|
|||
|
|
}
|
|||
|
|
m[col], m[piv] = m[piv], m[col]
|
|||
|
|
inv := new(big.Rat).Inv(m[col][col])
|
|||
|
|
for j := col; j < len(m[col]); j++ {
|
|||
|
|
m[col][j].Mul(m[col][j], inv)
|
|||
|
|
}
|
|||
|
|
for row := range n {
|
|||
|
|
if row == col || m[row][col].Cmp(zero) == 0 {
|
|||
|
|
continue
|
|||
|
|
}
|
|||
|
|
f := new(big.Rat).Set(m[row][col])
|
|||
|
|
for j := col; j < len(m[row]); j++ {
|
|||
|
|
m[row][j].Sub(m[row][j], new(big.Rat).Mul(f, m[col][j]))
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
out := make([][]*big.Rat, n)
|
|||
|
|
for i := range n {
|
|||
|
|
out[i] = m[i][n:]
|
|||
|
|
}
|
|||
|
|
return out
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// ratFroms converts float columns into exact rationals.
|
|||
|
|
func ratFroms(vs []float64) []*big.Rat {
|
|||
|
|
out := make([]*big.Rat, len(vs))
|
|||
|
|
for i, v := range vs {
|
|||
|
|
out[i] = new(big.Rat).SetFloat64(v)
|
|||
|
|
}
|
|||
|
|
return out
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// ratColumn turns a vector into the one-column right-hand side
|
|||
|
|
// ratSolve takes.
|
|||
|
|
func ratColumn(vs []*big.Rat) [][]*big.Rat {
|
|||
|
|
out := make([][]*big.Rat, len(vs))
|
|||
|
|
for i, v := range vs {
|
|||
|
|
out[i] = []*big.Rat{v}
|
|||
|
|
}
|
|||
|
|
return out
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// ratMatrixInvert inverts a small square rational matrix.
|
|||
|
|
func ratMatrixInvert(a [][]*big.Rat) [][]*big.Rat {
|
|||
|
|
n := len(a)
|
|||
|
|
id := make([][]*big.Rat, n)
|
|||
|
|
for i := range n {
|
|||
|
|
row := make([]*big.Rat, n)
|
|||
|
|
for j := range n {
|
|||
|
|
row[j] = new(big.Rat)
|
|||
|
|
if i == j {
|
|||
|
|
row[j].SetInt64(1)
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
id[i] = row
|
|||
|
|
}
|
|||
|
|
return ratSolve(a, id)
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// TestLevenbergMarquardtWeights pins the weighted fit against an
|
|||
|
|
// exact rational referent: chi2 becomes rᵀC⁻¹r, the parameters
|
|||
|
|
// minimise it, and the diagonal and matrix forms of Sigma agree on
|
|||
|
|
// the same weighting.
|
|||
|
|
func TestLevenbergMarquardtWeights(t *testing.T) {
|
|||
|
|
residual, jacobian := weightedModel()
|
|||
|
|
cov := []float64{
|
|||
|
|
2, 1, 0,
|
|||
|
|
1, 3, 1,
|
|||
|
|
0, 1, 2,
|
|||
|
|
}
|
|||
|
|
sigma, err := core.FromFloats(cov, 3, 3)
|
|||
|
|
if err != nil {
|
|||
|
|
t.Fatal(err)
|
|||
|
|
}
|
|||
|
|
p0 := mustFloats(t, []float64{0, 0}, 2)
|
|||
|
|
res, err := LevenbergMarquardtFit(residual, p0, LMOptions{Jacobian: jacobian, Sigma: sigma})
|
|||
|
|
if err != nil {
|
|||
|
|
t.Fatalf("LevenbergMarquardtFit: %v", err)
|
|||
|
|
}
|
|||
|
|
if res.Status != FitConverged {
|
|||
|
|
t.Fatalf("status = %d, want FitConverged", res.Status)
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// The exact referent: r = y − Xp, so the weighted optimum solves
|
|||
|
|
// (XᵀC⁻¹X) p = XᵀC⁻¹y.
|
|||
|
|
xf := [][]float64{{1, 0}, {0, 1}, {1, 1}}
|
|||
|
|
yf := []float64{1.2, 0.7, 2.4}
|
|||
|
|
cr := make([][]*big.Rat, 3)
|
|||
|
|
for i := range 3 {
|
|||
|
|
cr[i] = ratFroms(cov[i*3 : i*3+3])
|
|||
|
|
}
|
|||
|
|
cinv := ratMatrixInvert(cr)
|
|||
|
|
if cinv == nil {
|
|||
|
|
t.Fatal("the referent covariance is singular")
|
|||
|
|
}
|
|||
|
|
xr := make([][]*big.Rat, 3)
|
|||
|
|
for i := range 3 {
|
|||
|
|
xr[i] = ratFroms(xf[i])
|
|||
|
|
}
|
|||
|
|
yr := ratFroms(yf)
|
|||
|
|
cinvX := ratSolve(cr, xr)
|
|||
|
|
normal := make([][]*big.Rat, 2)
|
|||
|
|
for i := range 2 {
|
|||
|
|
normal[i] = make([]*big.Rat, 2)
|
|||
|
|
for j := range 2 {
|
|||
|
|
s := new(big.Rat)
|
|||
|
|
for k := range 3 {
|
|||
|
|
s.Add(s, new(big.Rat).Mul(xr[k][i], cinvX[k][j]))
|
|||
|
|
}
|
|||
|
|
normal[i][j] = s
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
rhs := make([]*big.Rat, 2)
|
|||
|
|
cinvY := ratSolve(cr, ratColumn(yr))
|
|||
|
|
if cinvY == nil {
|
|||
|
|
t.Fatal("the referent solve failed")
|
|||
|
|
}
|
|||
|
|
for i := range 2 {
|
|||
|
|
s := new(big.Rat)
|
|||
|
|
for k := range 3 {
|
|||
|
|
s.Add(s, new(big.Rat).Mul(xr[k][i], cinvY[k][0]))
|
|||
|
|
}
|
|||
|
|
rhs[i] = s
|
|||
|
|
}
|
|||
|
|
wantP := ratSolve(normal, ratColumn(rhs))
|
|||
|
|
if wantP == nil {
|
|||
|
|
t.Fatal("the referent normal equations are singular")
|
|||
|
|
}
|
|||
|
|
for j := range 2 {
|
|||
|
|
got := new(big.Rat).SetFloat64(res.Parameters.FloatAt(j))
|
|||
|
|
d := new(big.Rat).Sub(got, wantP[j][0])
|
|||
|
|
d.Abs(d)
|
|||
|
|
if d.Cmp(new(big.Rat).SetFloat64(1e-10)) > 0 {
|
|||
|
|
t.Fatalf("parameter %d = %.12g, want %s", j, res.Parameters.FloatAt(j), wantP[j][0].FloatString(12))
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
// The weighted chi2: rᵀC⁻¹r at the returned point.
|
|||
|
|
rAt := make([]*big.Rat, 3)
|
|||
|
|
for k := range 3 {
|
|||
|
|
rAt[k] = new(big.Rat).SetFloat64(yf[k] - (xf[k][0]*res.Parameters.FloatAt(0) + xf[k][1]*res.Parameters.FloatAt(1)))
|
|||
|
|
}
|
|||
|
|
cinvR := make([]*big.Rat, 3)
|
|||
|
|
for k := range 3 {
|
|||
|
|
s := new(big.Rat)
|
|||
|
|
for j := range 3 {
|
|||
|
|
s.Add(s, new(big.Rat).Mul(cinv[k][j], rAt[j]))
|
|||
|
|
}
|
|||
|
|
cinvR[k] = s
|
|||
|
|
}
|
|||
|
|
chiRef := new(big.Rat)
|
|||
|
|
for k := range 3 {
|
|||
|
|
chiRef.Add(chiRef, new(big.Rat).Mul(rAt[k], cinvR[k]))
|
|||
|
|
}
|
|||
|
|
chiFloat, _ := chiRef.Float64()
|
|||
|
|
if math.Abs(res.Chi2-chiFloat) > 1e-12*math.Max(1, math.Abs(chiFloat)) {
|
|||
|
|
t.Fatalf("chi2 = %.17g, want %.17g", res.Chi2, chiFloat)
|
|||
|
|
}
|
|||
|
|
// The weighted answer must differ from the unweighted one, or the
|
|||
|
|
// test proves nothing about the weighting.
|
|||
|
|
plain, _, err := LevenbergMarquardt(residual, p0, LMOptions{Jacobian: jacobian})
|
|||
|
|
if err != nil {
|
|||
|
|
t.Fatalf("the unweighted fit: %v", err)
|
|||
|
|
}
|
|||
|
|
if plain.FloatAt(0) == res.Parameters.FloatAt(0) && plain.FloatAt(1) == res.Parameters.FloatAt(1) {
|
|||
|
|
t.Fatal("the weighted and unweighted fits returned identical parameters")
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// The diagonal form weights by the same matrix's diagonal, and the
|
|||
|
|
// legacy entry point carries Sigma too.
|
|||
|
|
variance, err := core.FromFloats([]float64{2, 3, 2}, 3)
|
|||
|
|
if err != nil {
|
|||
|
|
t.Fatal(err)
|
|||
|
|
}
|
|||
|
|
diagRes, err := LevenbergMarquardtFit(residual, p0, LMOptions{Jacobian: jacobian, Sigma: variance})
|
|||
|
|
if err != nil {
|
|||
|
|
t.Fatalf("the diagonal Sigma fit: %v", err)
|
|||
|
|
}
|
|||
|
|
if diagRes.Status != FitConverged {
|
|||
|
|
t.Fatalf("the diagonal Sigma status = %d, want FitConverged", diagRes.Status)
|
|||
|
|
}
|
|||
|
|
wp, _, err := LevenbergMarquardt(residual, p0, LMOptions{Jacobian: jacobian, Sigma: sigma})
|
|||
|
|
if err != nil {
|
|||
|
|
t.Fatalf("the legacy weighted fit: %v", err)
|
|||
|
|
}
|
|||
|
|
if wp.FloatAt(0) != res.Parameters.FloatAt(0) || wp.FloatAt(1) != res.Parameters.FloatAt(1) {
|
|||
|
|
t.Fatal("the legacy weighted fit returned different parameters")
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// TestLevenbergMarquardtWeightedCovariance pins the covariance under
|
|||
|
|
// Sigma: (JᵀC⁻¹J)⁻¹ against the exact rational inverse, entry by
|
|||
|
|
// entry.
|
|||
|
|
func TestLevenbergMarquardtWeightedCovariance(t *testing.T) {
|
|||
|
|
residual, jacobian := weightedModel()
|
|||
|
|
cov := []float64{
|
|||
|
|
2, 1, 0,
|
|||
|
|
1, 3, 1,
|
|||
|
|
0, 1, 2,
|
|||
|
|
}
|
|||
|
|
sigma, err := core.FromFloats(cov, 3, 3)
|
|||
|
|
if err != nil {
|
|||
|
|
t.Fatal(err)
|
|||
|
|
}
|
|||
|
|
xf := [][]float64{{1, 0}, {0, 1}, {1, 1}}
|
|||
|
|
cr := make([][]*big.Rat, 3)
|
|||
|
|
for i := range 3 {
|
|||
|
|
cr[i] = ratFroms(cov[i*3 : i*3+3])
|
|||
|
|
}
|
|||
|
|
xr := make([][]*big.Rat, 3)
|
|||
|
|
for i := range 3 {
|
|||
|
|
xr[i] = ratFroms(xf[i])
|
|||
|
|
}
|
|||
|
|
cinvX := ratSolve(cr, xr)
|
|||
|
|
normal := make([][]*big.Rat, 2)
|
|||
|
|
for i := range 2 {
|
|||
|
|
normal[i] = make([]*big.Rat, 2)
|
|||
|
|
for j := range 2 {
|
|||
|
|
s := new(big.Rat)
|
|||
|
|
for k := range 3 {
|
|||
|
|
s.Add(s, new(big.Rat).Mul(xr[k][i], cinvX[k][j]))
|
|||
|
|
}
|
|||
|
|
normal[i][j] = s
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
want := ratMatrixInvert(normal)
|
|||
|
|
if want == nil {
|
|||
|
|
t.Fatal("the referent normal matrix is singular")
|
|||
|
|
}
|
|||
|
|
res, err := LevenbergMarquardtFit(residual, mustFloats(t, []float64{0, 0}, 2), LMOptions{
|
|||
|
|
Jacobian: jacobian,
|
|||
|
|
Sigma: sigma,
|
|||
|
|
RequestCovariance: true,
|
|||
|
|
})
|
|||
|
|
if err != nil {
|
|||
|
|
t.Fatalf("LevenbergMarquardtFit: %v", err)
|
|||
|
|
}
|
|||
|
|
if res.Covariance == nil {
|
|||
|
|
t.Fatal("the requested covariance is missing")
|
|||
|
|
}
|
|||
|
|
if res.Covariance.NDim() != 2 || res.Covariance.Shape()[0] != 2 || res.Covariance.Shape()[1] != 2 {
|
|||
|
|
t.Fatalf("covariance shape %s, want 2×2", base.ShapeText(res.Covariance.Shape()))
|
|||
|
|
}
|
|||
|
|
for i := range 2 {
|
|||
|
|
for j := range 2 {
|
|||
|
|
got := res.Covariance.FloatAt(i*2 + j)
|
|||
|
|
wantFloat, _ := want[i][j].Float64()
|
|||
|
|
if math.Abs(got-wantFloat) > 1e-9*math.Max(1, math.Abs(wantFloat)) {
|
|||
|
|
t.Fatalf("covariance (%d, %d) = %.17g, want %.17g", i, j, got, wantFloat)
|
|||
|
|
}
|
|||
|
|
mirror := res.Covariance.FloatAt(j*2 + i)
|
|||
|
|
if got != mirror {
|
|||
|
|
t.Fatalf("covariance (%d, %d) = %.17g against (%d, %d) = %.17g", i, j, got, j, i, mirror)
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// TestLevenbergMarquardtUnweightedCovariance pins the unweighted
|
|||
|
|
// covariance (JᵀJ)⁻¹ against the exact rational inverse, including on
|
|||
|
|
// the perfect start whose fit never enters the iteration loop.
|
|||
|
|
func TestLevenbergMarquardtUnweightedCovariance(t *testing.T) {
|
|||
|
|
residual, jacobian := weightedModel()
|
|||
|
|
xf := [][]float64{{1, 0}, {0, 1}, {1, 1}}
|
|||
|
|
xr := make([][]*big.Rat, 3)
|
|||
|
|
for i := range 3 {
|
|||
|
|
xr[i] = ratFroms(xf[i])
|
|||
|
|
}
|
|||
|
|
xtX := make([][]*big.Rat, 2)
|
|||
|
|
for i := range 2 {
|
|||
|
|
xtX[i] = make([]*big.Rat, 2)
|
|||
|
|
for j := range 2 {
|
|||
|
|
s := new(big.Rat)
|
|||
|
|
for k := range 3 {
|
|||
|
|
s.Add(s, new(big.Rat).Mul(xr[k][i], xr[k][j]))
|
|||
|
|
}
|
|||
|
|
xtX[i][j] = s
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
want := ratMatrixInvert(xtX)
|
|||
|
|
if want == nil {
|
|||
|
|
t.Fatal("the referent normal matrix is singular")
|
|||
|
|
}
|
|||
|
|
starts := [][]float64{{0, 0}, {1.2, 0.7}}
|
|||
|
|
for _, s := range starts {
|
|||
|
|
res, err := LevenbergMarquardtFit(residual, mustFloats(t, s, 2), LMOptions{
|
|||
|
|
Jacobian: jacobian,
|
|||
|
|
RequestCovariance: true,
|
|||
|
|
})
|
|||
|
|
if err != nil {
|
|||
|
|
t.Fatalf("LevenbergMarquardtFit(%v): %v", s, err)
|
|||
|
|
}
|
|||
|
|
if res.Covariance == nil {
|
|||
|
|
t.Fatalf("LevenbergMarquardtFit(%v): the requested covariance is missing", s)
|
|||
|
|
}
|
|||
|
|
for i := range 2 {
|
|||
|
|
for j := range 2 {
|
|||
|
|
got := res.Covariance.FloatAt(i*2 + j)
|
|||
|
|
wantFloat, _ := want[i][j].Float64()
|
|||
|
|
if math.Abs(got-wantFloat) > 1e-9*math.Max(1, math.Abs(wantFloat)) {
|
|||
|
|
t.Fatalf("start %v: covariance (%d, %d) = %.17g, want %.17g", s, i, j, got, wantFloat)
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// TestLevenbergMarquardtSigmaErrors pins the Sigma contract: shapes,
|
|||
|
|
// positive variances, exact symmetry and positive definiteness are
|
|||
|
|
// all checked before the fit moves.
|
|||
|
|
func TestLevenbergMarquardtSigmaErrors(t *testing.T) {
|
|||
|
|
residual, jacobian := weightedModel()
|
|||
|
|
p0 := mustFloats(t, []float64{0, 0}, 2)
|
|||
|
|
cases := []struct {
|
|||
|
|
name string
|
|||
|
|
sigma *core.Array
|
|||
|
|
want string
|
|||
|
|
}{
|
|||
|
|
{"wrong length", mustFloats(t, []float64{1, 2}, 2), "one variance per residual"},
|
|||
|
|
{"wrong shape", mustFloats(t, []float64{1, 0, 0, 1, 0, 0}, 2, 3), "shape (2, 3)"},
|
|||
|
|
{"zero variance", mustFloats(t, []float64{1, 0, 1}, 3), "positive variances"},
|
|||
|
|
{"negative variance", mustFloats(t, []float64{1, -2, 1}, 3), "positive variances"},
|
|||
|
|
}
|
|||
|
|
for _, tc := range cases {
|
|||
|
|
if _, err := LevenbergMarquardtFit(residual, p0, LMOptions{Jacobian: jacobian, Sigma: tc.sigma}); err == nil {
|
|||
|
|
t.Fatalf("%s: want an error", tc.name)
|
|||
|
|
} else if !strings.Contains(err.Error(), tc.want) {
|
|||
|
|
t.Fatalf("%s: error %q, want it to mention %q", tc.name, err, tc.want)
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
asymmetric := mustFloats(t, []float64{2, 1, 0, 1, 3, 1, 2, 1, 2}, 3, 3)
|
|||
|
|
if _, err := LevenbergMarquardtFit(residual, p0, LMOptions{Jacobian: jacobian, Sigma: asymmetric}); err == nil {
|
|||
|
|
t.Fatal("an asymmetric Sigma: want an error")
|
|||
|
|
} else if !strings.Contains(err.Error(), "symmetric") {
|
|||
|
|
t.Fatalf("asymmetric Sigma: error %q, want it to mention symmetry", err)
|
|||
|
|
}
|
|||
|
|
indefinite := mustFloats(t, []float64{1, 2, 2, 1}, 2, 2)
|
|||
|
|
wrongRows := mustFloats(t, []float64{1, 0, 0, 1}, 2, 2)
|
|||
|
|
if _, err := LevenbergMarquardtFit(residual, p0, LMOptions{Jacobian: jacobian, Sigma: indefinite}); err == nil {
|
|||
|
|
t.Fatal("an indefinite Sigma: want an error")
|
|||
|
|
}
|
|||
|
|
if _, err := LevenbergMarquardtFit(residual, p0, LMOptions{Jacobian: jacobian, Sigma: wrongRows}); err == nil {
|
|||
|
|
t.Fatal("a Sigma of the wrong row count: want an error")
|
|||
|
|
}
|
|||
|
|
}
|