bug fix

Author: Deltadahl

src/SimpleNonlinearSolve.jl

@@ -31,7 +31,7 @@ SnoopPrecompile.@precompile_all_calls begin for T in (Float32, Float64)
         solve(prob_no_brack, alg(), tol = T(1e-2))
     end
 
-    for alg in (TrustRegion(1.0),)
+    for alg in (TrustRegion(10.0),)
         solve(prob_no_brack, alg, tol = T(1e-2))
     end
 

src/trustRegion.jl

@@ -30,9 +30,9 @@ solver
 """
 struct TrustRegion{CS, AD, FDT} <: AbstractNewtonAlgorithm{CS, AD, FDT}
     max_trust_radius::Number
-    function TrustRegion(max_turst_radius::Number; chunk_size = Val{0}(),
-                                            autodiff = Val{true}(),
-                                            diff_type = Val{:forward})
+    function TrustRegion(max_trust_radius::Number; chunk_size = Val{0}(),
+                         autodiff = Val{true}(),
+                         diff_type = Val{:forward})
         new{SciMLBase._unwrap_val(chunk_size), SciMLBase._unwrap_val(autodiff),
             SciMLBase._unwrap_val(diff_type)}(max_trust_radius)
     end
@@ -46,8 +46,8 @@ function SciMLBase.solve(prob::NonlinearProblem,
     x = float(prob.u0)
     T = typeof(x)
     Δₘₐₓ = float(alg.max_trust_radius)  # The maximum trust region radius.
-    Δ = Δₘₐₓ / 5  # Initial trust region radius.
-    η₁ = 0.1   # Threshold for taking a step.
+    Δ = Δₘₐₓ / 11  # Initial trust region radius.
+    η₁ = 0.0   # Threshold for taking a step.
     η₂ = 0.25  # Threshold for shrinking the trust region.
     η₃ = 0.75  # Threshold for expanding the trust region.
     t₁ = 0.25  # Factor to shrink the trust region with.
@@ -88,38 +88,44 @@ function SciMLBase.solve(prob::NonlinearProblem,
 
         # Update the trust region radius.
         if r < η₂
-            Δ *= t₁
-        if r > η₁
-            if isapprox(x̂, x, atol = atol, rtol = rtol)
+            Δ = t₁ * Δ
+
+            if Δ < 1e-10
                 return SciMLBase.build_solution(prob, alg, x, F;
                                                 retcode = ReturnCode.Success)
             end
-
+        end
+        if r > η₁
+            if isapprox(xₖ₊₁, x, atol = atol, rtol = rtol)
+                return SciMLBase.build_solution(prob, alg, xₖ₊₁, Fₖ₊₁;
+                                                retcode = ReturnCode.Success)
+            end
+            # Take the step.
             x = xₖ₊₁
             F = Fₖ₊₁
             if alg_autodiff(alg)
                 F, ∇f = value_derivative(f, x)
             elseif x isa AbstractArray
                 ∇f = FiniteDiff.finite_difference_jacobian(f, x, diff_type(alg), eltype(x),
-                                                            F)
+                                                           F)
             else
                 ∇f = FiniteDiff.finite_difference_derivative(f, x, diff_type(alg),
-                                                              eltype(x),
-                                                              F)
+                                                             eltype(x),
+                                                             F)
             end
 
             iszero(F) &&
                 return SciMLBase.build_solution(prob, alg, x, F;
                                                 retcode = ReturnCode.Success)
+
             # Update the trust region radius.
             if r > η₃ && norm(δ) ≈ Δ
                 Δ = min(t₂ * Δ, Δₘₐₓ)
             end
-            fₖ = f̂
+            fₖ = fₖ₊₁
             H = ∇f * ∇f
             g = ∇f * F
         end
     end
-
     return SciMLBase.build_solution(prob, alg, x, F; retcode = ReturnCode.MaxIters)
 end

test/basictests.jl

@@ -46,10 +46,10 @@ sol = benchmark_scalar(sf, csu0)
 @test sol.u * sol.u - 2 < 1e-9
 @test (@ballocated benchmark_scalar(sf, csu0)) == 0
 
-# SimpleNewtonRaphsonTrustRegion
+# TrustRegion
 function benchmark_scalar(f, u0)
     probN = NonlinearProblem{false}(f, u0)
-    sol = (solve(probN, SimpleNewtonRaphsonTrustRegion(1.0)))
+    sol = (solve(probN, TrustRegion(10.0)))
 end
 
 sol = benchmark_scalar(sf, csu0)
@@ -63,7 +63,7 @@ using ForwardDiff
 f, u0 = (u, p) -> u .* u .- p, @SVector[1.0, 1.0]
 
 for alg in [SimpleNewtonRaphson(), Broyden(), Klement(),
-            SimpleNewtonRaphsonTrustRegion(10.0)]
+            TrustRegion(10.0)]
     g = function (p)
         probN = NonlinearProblem{false}(f, csu0, p)
         sol = solve(probN, alg, tol = 1e-9)
@@ -79,7 +79,7 @@ end
 # Scalar
 f, u0 = (u, p) -> u * u - p, 1.0
 for alg in [SimpleNewtonRaphson(), Broyden(), Klement(),
-            SimpleNewtonRaphsonTrustRegion(10.0)]
+            TrustRegion(10.0)]
     g = function (p)
         probN = NonlinearProblem{false}(f, oftype(p, u0), p)
         sol = solve(probN, alg)
@@ -121,7 +121,7 @@ for alg in [Bisection(), Falsi()]
 end
 
 for alg in [SimpleNewtonRaphson(), Broyden(), Klement(),
-            SimpleNewtonRaphsonTrustRegion(1.0)]
+            TrustRegion(10.0)]
     global g, p
     g = function (p)
         probN = NonlinearProblem{false}(f, 0.5, p)
@@ -141,10 +141,10 @@ probN = NonlinearProblem(f, u0)
 @test solve(probN, SimpleNewtonRaphson(; autodiff = false)).u[end] ≈ sqrt(2.0)
 @test solve(probN, SimpleNewtonRaphson(; autodiff = false)).u[end] ≈ sqrt(2.0)
 
-@test solve(probN, SimpleNewtonRaphsonTrustRegion(1.0)).u[end] ≈ sqrt(2.0)
-@test solve(probN, SimpleNewtonRaphsonTrustRegion(1.0); immutable = false).u[end] ≈ sqrt(2.0)
-@test solve(probN, SimpleNewtonRaphsonTrustRegion(1.0; autodiff = false)).u[end] ≈ sqrt(2.0)
-@test solve(probN, SimpleNewtonRaphsonTrustRegion(1.0; autodiff = false)).u[end] ≈ sqrt(2.0)
+@test solve(probN, TrustRegion(10.0)).u[end] ≈ sqrt(2.0)
+@test solve(probN, TrustRegion(10.0); immutable = false).u[end] ≈ sqrt(2.0)
+@test solve(probN, TrustRegion(10.0; autodiff = false)).u[end] ≈ sqrt(2.0)
+@test solve(probN, TrustRegion(10.0; autodiff = false)).u[end] ≈ sqrt(2.0)
 
 @test solve(probN, Broyden()).u[end] ≈ sqrt(2.0)
 @test solve(probN, Broyden(); immutable = false).u[end] ≈ sqrt(2.0)
@@ -162,9 +162,9 @@ for u0 in [1.0, [1, 1.0]]
     @test solve(probN, SimpleNewtonRaphson()).u ≈ sol
     @test solve(probN, SimpleNewtonRaphson(; autodiff = false)).u ≈ sol
 
-    @test solve(probN, SimpleNewtonRaphsonTrustRegion(1.0)).u ≈ sol
-    @test solve(probN, SimpleNewtonRaphsonTrustRegion(1.0)).u ≈ sol
-    @test solve(probN, SimpleNewtonRaphsonTrustRegion(1.0; autodiff = false)).u ≈ sol
+    @test solve(probN, TrustRegion(10.0)).u ≈ sol
+    @test solve(probN, TrustRegion(10.0)).u ≈ sol
+    @test solve(probN, TrustRegion(10.0; autodiff = false)).u ≈ sol
 
     @test solve(probN, Broyden()).u ≈ sol