diff --git a/src/Tearing/Dispersion/GrowthRateExtraction.jl b/src/Tearing/Dispersion/GrowthRateExtraction.jl index 45f5bc74dd..d16667889c 100644 --- a/src/Tearing/Dispersion/GrowthRateExtraction.jl +++ b/src/Tearing/Dispersion/GrowthRateExtraction.jl @@ -112,7 +112,10 @@ single-surface scans; `mc.surfaces[mc.ref_idx].tauk` for coupled scans). reports the raw marching-squares estimate, which is sensitive to grid and floating-point rounding. - `polish_maxit` -- max polish iterations per root (default 20). Each costs a - handful of `f` evaluations; failures fall back to the unpolished point. + handful of `f` evaluations (two for the real 2×2 Jacobian, one trial step, + at most three backtracks); failures fall back to the unpolished point. The + polish does not assume `f` is holomorphic, so it also converges on the SLAYER + residual, which is a function of conj(Q). - `validity_rtol` -- with `residual` supplied, a polished root is kept only if `|residual| < validity_rtol · median(|Δ|)` (default 1e-3). Intersections that don't polish to a true zero (contour near-misses, or surfaces whose @@ -500,9 +503,12 @@ end # Bounded Gauss-Newton minimization of |f(Q)| within a disk of radius `R` # around `Q0`. Returns (Q_best, n_evals, improved). `improved == false` ⇒ -# Q0 is returned unchanged (hard fallback). Secant derivative + trust-region -# projection + |f|-decrease backtracking keep it strictly "improve-or-no-op": -# it cannot diverge, jump to another root, or be pulled into a pole. +# Q0 is returned unchanged (hard fallback). Trust-region projection + +# |f|-decrease backtracking keep it strictly "improve-or-no-op": it cannot +# diverge, jump to another root, or be pulled into a pole. +# +# The local model is the real 2×2 forward-difference Jacobian of (Re f, Im f), exact for holomorphic +# and anti-holomorphic f alike; SLAYER's residual is a function of conj(Q) (see Riccati.jl). function _polish_root(f, Q0::ComplexF64, R::Float64; maxit::Int=8, tol_step::Float64=1e-12, tol_f::Float64=1e-8) (R > 0) || return (Q0, 0, false) @@ -512,17 +518,21 @@ function _polish_root(f, Q0::ComplexF64, R::Float64; maxit::Int=8, a0 = abs(f0) a0 == 0 && return (Q0, nev, false) - # Second point for the secant slope: a small real offset scaled to R. + h = max(R * 1e-3, eps(Float64)) # finite-difference offset, scaled to R Qk, fk = Q0, f0 - Qp = Q0 + complex(max(R * 1e-3, eps(Float64)), 0.0) - fp = ComplexF64(f(Qp)) - nev += 1 Qbest, abest = Q0, a0 for _ in 1:maxit - df = (Qp == Qk) ? ComplexF64(0) : (fp - fk) / (Qp - Qk) - (isfinite(df) && abs(df) > 0) || break - step = -fk / df # zero of local linear model + fx = ComplexF64(f(Qk + h)) + fy = ComplexF64(f(Qk + im * h)) + nev += 2 + (isfinite(fx) && isfinite(fy)) || break + j11, j12 = real(fx - fk) / h, real(fy - fk) / h + j21, j22 = imag(fx - fk) / h, imag(fy - fk) / h + jdet = j11 * j22 - j12 * j21 + (isfinite(jdet) && jdet != 0) || break + # zero of the local linear model: J * [dRe, dIm] = -[Re f, Im f] + step = complex((-real(fk) * j22 + j12 * imag(fk)) / jdet, (-j11 * imag(fk) + j21 * real(fk)) / jdet) Qt = Qk + step if abs(Qt - Q0) > R # project into trust region Qt = Q0 + R * (Qt - Q0) / abs(Qt - Q0) @@ -538,7 +548,6 @@ function _polish_root(f, Q0::ComplexF64, R::Float64; maxit::Int=8, bt += 1 end if isfinite(ft) && abs(ft) < abest - Qp, fp = Qk, fk Qk, fk = Qt, ft Qbest, abest = Qt, abs(ft) (abs(step) < tol_step * (abs(Qk) + tol_step) || abest < tol_f * a0) && break diff --git a/test/runtests_dispersion_polish.jl b/test/runtests_dispersion_polish.jl index 8bc4c97025..43d06363b2 100644 --- a/test/runtests_dispersion_polish.jl +++ b/test/runtests_dispersion_polish.jl @@ -14,7 +14,7 @@ @test improved @test abs(Qp - r) < 1e-10 @test abs(f(Qp)) < 1e-10 - @test nev <= 3 * 8 + 2 # bounded eval budget + @test nev <= 1 + 6 # one iteration: 2 differences, the trial step, ≤ 3 backtracks end @testset "_polish_root: nonlinear residual converges" begin @@ -25,6 +25,43 @@ @test abs(Qp - r) < 1e-7 end + @testset "_polish_root: residual that is neither holomorphic nor anti-holomorphic converges" begin + r = 0.15 - 0.35im + f(Q) = (ComplexF64(Q) - r) + 0.5 * conj(ComplexF64(Q) - r) # zero only at r + Qp, _, improved = D._polish_root(f, r + 0.01 + 0.012im, 0.1) + @test improved + @test abs(Qp - r) < 1e-10 + end + + # The SLAYER residual is a function of conj(Q) (solve_inner evaluates + # conj(Δ̂(ĝ))), so the polish must not assume a holomorphic f. + @testset "_polish_root: anti-holomorphic residual converges" begin + r = 0.42 + 0.27im + f(Q) = (2 - 1im) * conj(ComplexF64(Q) - r) + Qp, _, improved = D._polish_root(f, r + (0.013 - 0.009im), 0.1) + @test improved + @test abs(Qp - r) < 1e-10 + @test abs(f(Qp)) < 1e-9 + end + + @testset "_polish_root: anti-holomorphic root offset along Im Q (mirrored complex step)" begin + # A complex-slope step on f = a·conj(Q - r) maps an offset i·y to 2i·y and + # every halving backtrack stays worse, so a complex-secant polish stalls here. + r = -0.75 - 0.009im + f(Q) = 1e3 * conj(ComplexF64(Q) - r) + Qp, _, improved = D._polish_root(f, r + 0.0017im, 0.01) + @test improved + @test abs(Qp - r) < 1e-10 + end + + @testset "_polish_root: steep nonlinear anti-holomorphic residual converges" begin + r = -0.75 - 0.009im + f(Q) = 1e3 * conj((ComplexF64(Q) - r) * (1 + 0.4 * (ComplexF64(Q) - r))) + Qp, _, improved = D._polish_root(f, r + (0.004 + 0.003im), 0.02) + @test improved + @test abs(Qp - r) < 1e-9 + end + @testset "_polish_root: no-op-or-better invariant (never worse)" begin # The polish must never return a point with larger |f| than the start. for r in (0.1 + 0.2im, -0.5 + 0.0im, 0.3 - 0.4im) @@ -181,4 +218,17 @@ @test !(:spurious in gr3.warning_flags) @test !isnan(real(gr3.Q_root)) end + + @testset "find_growth_rates: steep anti-holomorphic root passes the validity gate" begin + # A steep residual leaves the raw contour intersection far above the gate + # threshold, so the root is kept only if the polish reaches the true zero. + r = 0.42 + 0.27im + f(Q) = 1e3 * conj((ComplexF64(Q) - r) * (1 + 0.4 * (ComplexF64(Q) - r))) + scan = brute_force_scan(f, (0.0, 1.0), (0.0, 0.6); nre=9, nim=7, threaded=false) # raw |f| above the gate + # the pole classifier's threshold is absolute: scale it with f, as SLAYER's adaptive one does + gr = find_growth_rates(scan, 1.0; residual=f, validity_rtol=1e-3, pole_threshold=1e3 * 10.0) + @test !(:spurious in gr.warning_flags) + @test !isnan(real(gr.Q_root)) + @test abs(gr.Q_root - r) < 1e-8 + end end