From d9ea0205f466bc690b3fc6071227f1a930225f68 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Fri, 2 Jun 2023 05:40:18 -0700 Subject: [PATCH 1/2] format --- docs/src/solvers/NonlinearSystemSolvers.md | 6 +- src/NonlinearSolve.jl | 31 ++- src/jacobian.jl | 2 - src/trustRegion.jl | 284 +++++++++++---------- src/utils.jl | 2 +- test/basictests.jl | 80 +++--- test/runtests.jl | 6 +- 7 files changed, 221 insertions(+), 190 deletions(-) diff --git a/docs/src/solvers/NonlinearSystemSolvers.md b/docs/src/solvers/NonlinearSystemSolvers.md index a19a3f8a7..b1781a435 100644 --- a/docs/src/solvers/NonlinearSystemSolvers.md +++ b/docs/src/solvers/NonlinearSystemSolvers.md @@ -21,13 +21,13 @@ is small, or the eigenvalues of the Jacobian are within a few orders of magnitud then `NLSolveJL`'s `:anderson` can be a good choice. !!! note - + `TrustRegion` and `SimpleTrustRegion` are still in development. ## Full List of Methods !!! note - + For the full details on the capabilities and constructors of the different solvers, see the Detailed Solver APIs section! @@ -60,7 +60,7 @@ methods excel at small problems and problems defined with static arrays. large-scale nonlinear systems of equations. !!! note - + When used with certain types for the states `u` such as a `Number` or `StaticArray`, these solvers are very efficient and non-allocating. These implementations are thus well-suited for small systems of equations. diff --git a/src/NonlinearSolve.jl b/src/NonlinearSolve.jl index a83660705..6b4b6274e 100644 --- a/src/NonlinearSolve.jl +++ b/src/NonlinearSolve.jl @@ -41,24 +41,27 @@ include("ad.jl") import PrecompileTools -PrecompileTools.@compile_workload begin for T in (Float32, Float64) - prob = NonlinearProblem{false}((u, p) -> u .* u .- p, T(0.1), T(2)) +PrecompileTools.@compile_workload begin + for T in (Float32, Float64) + prob = NonlinearProblem{false}((u, p) -> u .* u .- p, T(0.1), T(2)) - precompile_algs = if VERSION >= v"1.7" - (NewtonRaphson(), TrustRegion(), LevenbergMarquardt()) - else - (NewtonRaphson(),) - end + precompile_algs = if VERSION >= v"1.7" + (NewtonRaphson(), TrustRegion(), LevenbergMarquardt()) + else + (NewtonRaphson(),) + end - for alg in precompile_algs - solve(prob, alg, abstol = T(1e-2)) - end + for alg in precompile_algs + solve(prob, alg, abstol = T(1e-2)) + end - prob = NonlinearProblem{true}((du, u, p) -> du[1] = u[1] * u[1] - p[1], T[0.1], T[2]) - for alg in precompile_algs - solve(prob, alg, abstol = T(1e-2)) + prob = NonlinearProblem{true}((du, u, p) -> du[1] = u[1] * u[1] - p[1], T[0.1], + T[2]) + for alg in precompile_algs + solve(prob, alg, abstol = T(1e-2)) + end end -end end +end export RadiusUpdateSchemes diff --git a/src/jacobian.jl b/src/jacobian.jl index 9d8e3bfc5..7661d8e55 100644 --- a/src/jacobian.jl +++ b/src/jacobian.jl @@ -144,5 +144,3 @@ function jacobian_autodiff(f, x::AbstractArray, nonlinfun, alg) jac_prototype = jac_prototype, chunksize = chunk_size), num_of_chunks) end - - diff --git a/src/trustRegion.jl b/src/trustRegion.jl index b3af32fbe..06ede64ba 100644 --- a/src/trustRegion.jl +++ b/src/trustRegion.jl @@ -177,15 +177,19 @@ mutable struct TrustRegionCache{iip, fType, algType, uType, resType, pType, jac_config::JC, iter::Int, force_stop::Bool, maxiters::Int, internalnorm::INType, retcode::SciMLBase.ReturnCode.T, abstol::tolType, - prob::probType, radius_update_scheme::RadiusUpdateSchemes.T, trust_r::trustType, + prob::probType, + radius_update_scheme::RadiusUpdateSchemes.T, + trust_r::trustType, max_trust_r::trustType, step_threshold::suType, shrink_threshold::trustType, expand_threshold::trustType, shrink_factor::trustType, expand_factor::trustType, loss::floatType, loss_new::floatType, H::jType, g::resType, shrink_counter::Int, step_size::su2Type, u_tmp::tmpType, fu_new::resType, make_new_J::Bool, - r::floatType, p1::floatType, p2::floatType, p3::floatType, - p4::floatType, ϵ::floatType) where {iip, fType, algType, uType, + r::floatType, p1::floatType, p2::floatType, + p3::floatType, + p4::floatType, + ϵ::floatType) where {iip, fType, algType, uType, resType, pType, INType, tolType, probType, ufType, L, jType, JC, floatType, trustType, @@ -195,7 +199,8 @@ mutable struct TrustRegionCache{iip, fType, algType, uType, resType, pType, trustType, suType, su2Type, tmpType}(f, alg, u, fu, p, uf, linsolve, J, jac_config, iter, force_stop, maxiters, internalnorm, retcode, - abstol, prob, radius_update_scheme, trust_r, max_trust_r, + abstol, prob, radius_update_scheme, + trust_r, max_trust_r, step_threshold, shrink_threshold, expand_threshold, shrink_factor, expand_factor, loss, @@ -287,47 +292,47 @@ function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg::TrustRegion, p4 = convert(eltype(u), 0.0) ϵ = convert(eltype(u), 1.0e-8) if radius_update_scheme === RadiusUpdateSchemes.Hei - step_threshold = convert(eltype(u), 0.0) - shrink_threshold = convert(eltype(u), 0.25) - expand_threshold = convert(eltype(u), 0.25) - p1 = convert(eltype(u), 5.0) # M - p2 = convert(eltype(u), 0.1) # β - p3 = convert(eltype(u), 0.15) # γ1 - p4 = convert(eltype(u), 0.15) # γ2 - initial_trust_radius = convert(eltype(u), 1.0) + step_threshold = convert(eltype(u), 0.0) + shrink_threshold = convert(eltype(u), 0.25) + expand_threshold = convert(eltype(u), 0.25) + p1 = convert(eltype(u), 5.0) # M + p2 = convert(eltype(u), 0.1) # β + p3 = convert(eltype(u), 0.15) # γ1 + p4 = convert(eltype(u), 0.15) # γ2 + initial_trust_radius = convert(eltype(u), 1.0) elseif radius_update_scheme === RadiusUpdateSchemes.Yuan - step_threshold = convert(eltype(u), 0.0001) - shrink_threshold = convert(eltype(u), 0.25) - expand_threshold = convert(eltype(u), 0.25) - p1 = convert(eltype(u), 2.0) # μ - p2 = convert(eltype(u), 1/6) # c5 - p3 = convert(eltype(u), 6.0) # c6 - p4 = convert(eltype(u), 0.0) - if iip - auto_jacvec!(g, (fu, x) -> f(fu, x, p), u, fu) - else - if isa(u, Number) - g = ForwardDiff.derivative(x -> f(x, p), u) + step_threshold = convert(eltype(u), 0.0001) + shrink_threshold = convert(eltype(u), 0.25) + expand_threshold = convert(eltype(u), 0.25) + p1 = convert(eltype(u), 2.0) # μ + p2 = convert(eltype(u), 1 / 6) # c5 + p3 = convert(eltype(u), 6.0) # c6 + p4 = convert(eltype(u), 0.0) + if iip + auto_jacvec!(g, (fu, x) -> f(fu, x, p), u, fu) else - g = auto_jacvec(x -> f(x, p), u, fu) + if isa(u, Number) + g = ForwardDiff.derivative(x -> f(x, p), u) + else + g = auto_jacvec(x -> f(x, p), u, fu) + end end - end - initial_trust_radius = convert(eltype(u), p1 * norm(g)) + initial_trust_radius = convert(eltype(u), p1 * norm(g)) elseif radius_update_scheme === RadiusUpdateSchemes.Fan - step_threshold = convert(eltype(u), 0.0001) - shrink_threshold = convert(eltype(u), 0.25) - expand_threshold = convert(eltype(u), 0.75) - p1 = convert(eltype(u), 0.1) # μ - p2 = convert(eltype(u), 1/4) # c5 - p3 = convert(eltype(u), 12) # c6 - p4 = convert(eltype(u), 1.0e18) # M - initial_trust_radius = convert(eltype(u), p1 * (norm(fu)^0.99)) + step_threshold = convert(eltype(u), 0.0001) + shrink_threshold = convert(eltype(u), 0.25) + expand_threshold = convert(eltype(u), 0.75) + p1 = convert(eltype(u), 0.1) # μ + p2 = convert(eltype(u), 1 / 4) # c5 + p3 = convert(eltype(u), 12) # c6 + p4 = convert(eltype(u), 1.0e18) # M + initial_trust_radius = convert(eltype(u), p1 * (norm(fu)^0.99)) end - return TrustRegionCache{iip}(f, alg, u, fu, p, uf, linsolve, J, jac_config, 1, false, maxiters, internalnorm, - ReturnCode.Default, abstol, prob, radius_update_scheme, initial_trust_radius, + ReturnCode.Default, abstol, prob, radius_update_scheme, + initial_trust_radius, max_trust_radius, step_threshold, shrink_threshold, expand_threshold, shrink_factor, expand_factor, loss, loss_new, H, g, shrink_counter, step_size, u_tmp, fu_new, @@ -385,100 +390,105 @@ function trust_region_step!(cache::TrustRegionCache) cache.r = -(loss - cache.loss_new) / (step_size' * g + step_size' * H * step_size / 2) @unpack r = cache - if radius_update_scheme === RadiusUpdateSchemes.Simple - # Update the trust region radius. - if r < cache.shrink_threshold - cache.trust_r *= cache.shrink_factor - cache.shrink_counter += 1 - else - cache.shrink_counter = 0 - end - if r > cache.step_threshold - take_step!(cache) - cache.loss = cache.loss_new - - # Update the trust region radius. - if r > cache.expand_threshold - cache.trust_r = min(cache.expand_factor * cache.trust_r, max_trust_r) - end - - cache.make_new_J = true - else - # No need to make a new J, no step was taken, so we try again with a smaller trust_r - cache.make_new_J = false - end - - if iszero(cache.fu) || cache.internalnorm(cache.fu) < cache.abstol - cache.force_stop = true - end - + if radius_update_scheme === RadiusUpdateSchemes.Simple + # Update the trust region radius. + if r < cache.shrink_threshold + cache.trust_r *= cache.shrink_factor + cache.shrink_counter += 1 + else + cache.shrink_counter = 0 + end + if r > cache.step_threshold + take_step!(cache) + cache.loss = cache.loss_new + + # Update the trust region radius. + if r > cache.expand_threshold + cache.trust_r = min(cache.expand_factor * cache.trust_r, max_trust_r) + end + + cache.make_new_J = true + else + # No need to make a new J, no step was taken, so we try again with a smaller trust_r + cache.make_new_J = false + end + + if iszero(cache.fu) || cache.internalnorm(cache.fu) < cache.abstol + cache.force_stop = true + end + elseif radius_update_scheme === RadiusUpdateSchemes.Hei - if r > cache.step_threshold - take_step!(cache) - cache.loss = cache.loss_new - cache.make_new_J = true - else - cache.make_new_J = false - end - # Hei's radius update scheme - @unpack shrink_threshold, p1, p2, p3, p4 = cache - if rfunc(r, shrink_threshold, p1, p3, p4, p2) * cache.internalnorm(step_size) < cache.trust_r - cache.shrink_counter += 1 - else - cache.shrink_counter = 0 - end - cache.trust_r = rfunc(r, shrink_threshold, p1, p3, p4, p2) * cache.internalnorm(step_size) - - if iszero(cache.fu) || cache.internalnorm(cache.fu) < cache.abstol || cache.internalnorm(g) < cache.ϵ - cache.force_stop = true - end + if r > cache.step_threshold + take_step!(cache) + cache.loss = cache.loss_new + cache.make_new_J = true + else + cache.make_new_J = false + end + # Hei's radius update scheme + @unpack shrink_threshold, p1, p2, p3, p4 = cache + if rfunc(r, shrink_threshold, p1, p3, p4, p2) * cache.internalnorm(step_size) < + cache.trust_r + cache.shrink_counter += 1 + else + cache.shrink_counter = 0 + end + cache.trust_r = rfunc(r, shrink_threshold, p1, p3, p4, p2) * + cache.internalnorm(step_size) + if iszero(cache.fu) || cache.internalnorm(cache.fu) < cache.abstol || + cache.internalnorm(g) < cache.ϵ + cache.force_stop = true + end elseif radius_update_scheme === RadiusUpdateSchemes.Yuan - if r < cache.shrink_threshold - cache.p1 = cache.p2 * cache.p1 - cache.shrink_counter += 1 - elseif r >= cache.expand_threshold && cache.internalnorm(step_size) > cache.trust_r / 2 - cache.p1 = cache.p3 * cache.p1 - cache.shrink_counter = 0 - end - - if r > cache.step_threshold - take_step!(cache) - cache.loss = cache.loss_new - cache.make_new_J = true - else - cache.make_new_J = false - end - - @unpack p1= cache - cache.trust_r = p1 * cache.internalnorm(jvp!(cache)) - if iszero(cache.fu) || cache.internalnorm(cache.fu) < cache.abstol || cache.internalnorm(g) < cache.ϵ - cache.force_stop = true - end - #Fan's update scheme + if r < cache.shrink_threshold + cache.p1 = cache.p2 * cache.p1 + cache.shrink_counter += 1 + elseif r >= cache.expand_threshold && + cache.internalnorm(step_size) > cache.trust_r / 2 + cache.p1 = cache.p3 * cache.p1 + cache.shrink_counter = 0 + end + + if r > cache.step_threshold + take_step!(cache) + cache.loss = cache.loss_new + cache.make_new_J = true + else + cache.make_new_J = false + end + + @unpack p1 = cache + cache.trust_r = p1 * cache.internalnorm(jvp!(cache)) + if iszero(cache.fu) || cache.internalnorm(cache.fu) < cache.abstol || + cache.internalnorm(g) < cache.ϵ + cache.force_stop = true + end + #Fan's update scheme elseif radius_update_scheme === RadiusUpdateSchemes.Fan - if r < cache.shrink_threshold - cache.p1 *= cache.p2 - cache.shrink_counter += 1 - elseif r > cache.expand_threshold - cache.p1 = min(cache.p1*cache.p3, cache.p4) - cache.shrink_counter = 0 - end - - if r > cache.step_threshold - take_step!(cache) - cache.loss = cache.loss_new - cache.make_new_J = true - else - cache.make_new_J = false - end - - @unpack p1 = cache - cache.trust_r = p1 * (cache.internalnorm(cache.fu)^0.99) - if iszero(cache.fu) || cache.internalnorm(cache.fu) < cache.abstol || cache.internalnorm(g) < cache.ϵ - cache.force_stop = true - end + if r < cache.shrink_threshold + cache.p1 *= cache.p2 + cache.shrink_counter += 1 + elseif r > cache.expand_threshold + cache.p1 = min(cache.p1 * cache.p3, cache.p4) + cache.shrink_counter = 0 + end + + if r > cache.step_threshold + take_step!(cache) + cache.loss = cache.loss_new + cache.make_new_J = true + else + cache.make_new_J = false + end + + @unpack p1 = cache + cache.trust_r = p1 * (cache.internalnorm(cache.fu)^0.99) + if iszero(cache.fu) || cache.internalnorm(cache.fu) < cache.abstol || + cache.internalnorm(g) < cache.ϵ + cache.force_stop = true + end end end @@ -520,20 +530,20 @@ function take_step!(cache::TrustRegionCache{false}) end function jvp!(cache::TrustRegionCache{false}) - @unpack f, u, fu, p = cache - if isa(u, Number) - return value_derivative(x -> f(x, p), u) - end - return auto_jacvec(x -> f(x, p), u, fu) + @unpack f, u, fu, p = cache + if isa(u, Number) + return value_derivative(x -> f(x, p), u) + end + return auto_jacvec(x -> f(x, p), u, fu) end function jvp!(cache::TrustRegionCache{true}) - @unpack g, f, u, fu, p = cache - if isa(u, Number) - return value_derivative(x -> f(x, p), u) - end - auto_jacvec!(g, (fu, x) -> f(fu, x, p), u, fu) - g + @unpack g, f, u, fu, p = cache + if isa(u, Number) + return value_derivative(x -> f(x, p), u) + end + auto_jacvec!(g, (fu, x) -> f(fu, x, p), u, fu) + g end function SciMLBase.solve!(cache::TrustRegionCache) diff --git a/src/utils.jl b/src/utils.jl index f8fd04188..3e2fb413e 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -86,7 +86,7 @@ function dolinsolve(precs::P, linsolve; A = nothing, linu = nothing, b = nothing (linsolve.Pr isa Diagonal ? linsolve.Pr.diag : linsolve.Pr.inner.diag) : weight Pl, Pr = wrapprecs(_Pl, _Pr, _weight) - linsolve.Pl = Pl + linsolve.Pl = Pl linsolve.Pr = Pr end diff --git a/test/basictests.jl b/test/basictests.jl index 4dc696162..9dce5d45e 100644 --- a/test/basictests.jl +++ b/test/basictests.jl @@ -115,7 +115,7 @@ end f = (u, p) -> u * u - p g = function (p_range) probN = NonlinearProblem{false}(f, 0.5, p_range[begin]) - cache = init(probN, NewtonRaphson(); maxiters = 100, abstol=1e-10) + cache = init(probN, NewtonRaphson(); maxiters = 100, abstol = 1e-10) sols = zeros(length(p_range)) for (i, p) in enumerate(p_range) reinit!(cache, cache.u; p = p) @@ -130,7 +130,7 @@ p = range(0.01, 2, length = 200) f = (res, u, p) -> (res[begin] = u[1] * u[1] - p) g = function (p_range) probN = NonlinearProblem{true}(f, [0.5], p_range[begin]) - cache = init(probN, NewtonRaphson(); maxiters = 100, abstol=1e-10) + cache = init(probN, NewtonRaphson(); maxiters = 100, abstol = 1e-10) sols = zeros(length(p_range)) for (i, p) in enumerate(p_range) reinit!(cache, [cache.u[1]]; p = p) @@ -165,32 +165,36 @@ end function benchmark_immutable(f, u0, radius_update_scheme) probN = NonlinearProblem{false}(f, u0) - solver = init(probN, TrustRegion(radius_update_scheme = radius_update_scheme), abstol = 1e-9) + solver = init(probN, TrustRegion(radius_update_scheme = radius_update_scheme), + abstol = 1e-9) sol = solve!(solver) end function benchmark_mutable(f, u0, radius_update_scheme) probN = NonlinearProblem{false}(f, u0) - solver = init(probN, TrustRegion(radius_update_scheme = radius_update_scheme), abstol = 1e-9) + solver = init(probN, TrustRegion(radius_update_scheme = radius_update_scheme), + abstol = 1e-9) sol = solve!(solver) end function benchmark_scalar(f, u0, radius_update_scheme) probN = NonlinearProblem{false}(f, u0) - sol = (solve(probN, TrustRegion(radius_update_scheme = radius_update_scheme), abstol = 1e-9)) + sol = (solve(probN, TrustRegion(radius_update_scheme = radius_update_scheme), + abstol = 1e-9)) end -function ff(u, p=nothing) +function ff(u, p = nothing) u .* u .- 2 end -function sf(u, p=nothing) +function sf(u, p = nothing) u * u - 2 end u0 = [1.0, 1.0] -radius_update_schemes = [RadiusUpdateSchemes.Simple, RadiusUpdateSchemes.Hei, RadiusUpdateSchemes.Yuan, - RadiusUpdateSchemes.Fan] +radius_update_schemes = [RadiusUpdateSchemes.Simple, RadiusUpdateSchemes.Hei, + RadiusUpdateSchemes.Yuan, + RadiusUpdateSchemes.Fan] for radius_update_scheme in radius_update_schemes sol = benchmark_immutable(ff, cu0, radius_update_scheme) @@ -204,14 +208,14 @@ for radius_update_scheme in radius_update_schemes @test abs(sol.u * sol.u - 2) < 1e-9 end - function benchmark_inplace(f, u0, radius_update_scheme) probN = NonlinearProblem{true}(f, u0) - solver = init(probN, TrustRegion(radius_update_scheme = radius_update_scheme), abstol = 1e-9) + solver = init(probN, TrustRegion(radius_update_scheme = radius_update_scheme), + abstol = 1e-9) sol = solve!(solver) end -function ffiip(du, u, p=nothing) +function ffiip(du, u, p = nothing) du .= u .* u .- 2 end u0 = [1.0, 1.0] @@ -224,7 +228,8 @@ end for radius_update_scheme in radius_update_schemes probN = NonlinearProblem{true}(ffiip, u0) - solver = init(probN, TrustRegion(radius_update_scheme = radius_update_scheme), abstol = 1e-9) + solver = init(probN, TrustRegion(radius_update_scheme = radius_update_scheme), + abstol = 1e-9) @test (@ballocated solve!(solver)) < 200 end @@ -247,7 +252,8 @@ end g = function (p) probN = NonlinearProblem{false}(f, csu0, p) - sol = solve(probN, TrustRegion(radius_update_scheme = RadiusUpdateSchemes.Hei), abstol = 1e-9) + sol = solve(probN, TrustRegion(radius_update_scheme = RadiusUpdateSchemes.Hei), + abstol = 1e-9) return sol.u[end] end @@ -258,7 +264,8 @@ end g = function (p) probN = NonlinearProblem{false}(f, csu0, p) - sol = solve(probN, TrustRegion(radius_update_scheme = RadiusUpdateSchemes.Yuan), abstol = 1e-9) + sol = solve(probN, TrustRegion(radius_update_scheme = RadiusUpdateSchemes.Yuan), + abstol = 1e-9) return sol.u[end] end @@ -269,7 +276,8 @@ end g = function (p) probN = NonlinearProblem{false}(f, csu0, p) - sol = solve(probN, TrustRegion(radius_update_scheme = RadiusUpdateSchemes.Fan), abstol = 1e-9) + sol = solve(probN, TrustRegion(radius_update_scheme = RadiusUpdateSchemes.Fan), + abstol = 1e-9) return sol.u[end] end @@ -296,7 +304,8 @@ end g = function (p) probN = NonlinearProblem{false}(f, oftype(p, u0), p) - sol = solve(probN, TrustRegion(radius_update_scheme = RadiusUpdateSchemes.Hei), abstol = 1e-10) + sol = solve(probN, TrustRegion(radius_update_scheme = RadiusUpdateSchemes.Hei), + abstol = 1e-10) return sol.u end @@ -309,7 +318,8 @@ end g = function (p) probN = NonlinearProblem{false}(f, oftype(p, u0), p) - sol = solve(probN, TrustRegion(radius_update_scheme = RadiusUpdateSchemes.Yuan), abstol = 1e-10) + sol = solve(probN, TrustRegion(radius_update_scheme = RadiusUpdateSchemes.Yuan), + abstol = 1e-10) return sol.u end @@ -322,7 +332,8 @@ end g = function (p) probN = NonlinearProblem{false}(f, oftype(p, u0), p) - sol = solve(probN, TrustRegion(radius_update_scheme = RadiusUpdateSchemes.Fan), abstol = 1e-10) + sol = solve(probN, TrustRegion(radius_update_scheme = RadiusUpdateSchemes.Fan), + abstol = 1e-10) return sol.u end @@ -372,7 +383,7 @@ end f = (u, p) -> u * u - p g = function (p_range) probN = NonlinearProblem{false}(f, 0.5, p_range[begin]) - cache = init(probN, TrustRegion(); maxiters = 100, abstol=1e-10) + cache = init(probN, TrustRegion(); maxiters = 100, abstol = 1e-10) sols = zeros(length(p_range)) for (i, p) in enumerate(p_range) reinit!(cache, cache.u; p = p) @@ -387,7 +398,7 @@ p = range(0.01, 2, length = 200) f = (res, u, p) -> (res[begin] = u[1] * u[1] - p) g = function (p_range) probN = NonlinearProblem{true}(f, [0.5], p_range[begin]) - cache = init(probN, TrustRegion(); maxiters = 100, abstol=1e-10) + cache = init(probN, TrustRegion(); maxiters = 100, abstol = 1e-10) sols = zeros(length(p_range)) for (i, p) in enumerate(p_range) reinit!(cache, [cache.u[1]]; p = p) @@ -406,14 +417,20 @@ probN = NonlinearProblem(f, u0) @test solve(probN, TrustRegion()).u[end] ≈ sqrt(2.0) @test solve(probN, TrustRegion(; autodiff = false)).u[end] ≈ sqrt(2.0) -@test solve(probN, TrustRegion(radius_update_scheme = RadiusUpdateSchemes.Hei)).u[end] ≈ sqrt(2.0) -@test solve(probN, TrustRegion(; radius_update_scheme = RadiusUpdateSchemes.Hei, autodiff = false)).u[end] ≈ sqrt(2.0) +@test solve(probN, TrustRegion(radius_update_scheme = RadiusUpdateSchemes.Hei)).u[end] ≈ + sqrt(2.0) +@test solve(probN, TrustRegion(; radius_update_scheme = RadiusUpdateSchemes.Hei, autodiff = false)).u[end] ≈ + sqrt(2.0) -@test solve(probN, TrustRegion(radius_update_scheme = RadiusUpdateSchemes.Yuan)).u[end] ≈ sqrt(2.0) -@test solve(probN, TrustRegion(; radius_update_scheme = RadiusUpdateSchemes.Yuan, autodiff = false)).u[end] ≈ sqrt(2.0) +@test solve(probN, TrustRegion(radius_update_scheme = RadiusUpdateSchemes.Yuan)).u[end] ≈ + sqrt(2.0) +@test solve(probN, TrustRegion(; radius_update_scheme = RadiusUpdateSchemes.Yuan, autodiff = false)).u[end] ≈ + sqrt(2.0) -@test solve(probN, TrustRegion(radius_update_scheme = RadiusUpdateSchemes.Fan)).u[end] ≈ sqrt(2.0) -@test solve(probN, TrustRegion(; radius_update_scheme = RadiusUpdateSchemes.Fan, autodiff = false)).u[end] ≈ sqrt(2.0) +@test solve(probN, TrustRegion(radius_update_scheme = RadiusUpdateSchemes.Fan)).u[end] ≈ + sqrt(2.0) +@test solve(probN, TrustRegion(; radius_update_scheme = RadiusUpdateSchemes.Fan, autodiff = false)).u[end] ≈ + sqrt(2.0) for u0 in [1.0, [1, 1.0]] local f, probN, sol @@ -449,7 +466,8 @@ f(u, p) g = function (p) probN = NonlinearProblem{false}(f, u0, p) - sol = solve(probN, TrustRegion(radius_update_scheme = RadiusUpdateSchemes.Fan), abstol = 1e-10) + sol = solve(probN, TrustRegion(radius_update_scheme = RadiusUpdateSchemes.Fan), + abstol = 1e-10) return sol.u end p = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] @@ -492,11 +510,13 @@ maxiterations = [2, 3, 4, 5] u0 = [1.0, 1.0] function iip_oop(f, fip, u0, radius_update_scheme, maxiters) prob_iip = NonlinearProblem{true}(fip, u0) - solver = init(prob_iip, TrustRegion(radius_update_scheme = radius_update_scheme), abstol = 1e-9, maxiters = maxiters) + solver = init(prob_iip, TrustRegion(radius_update_scheme = radius_update_scheme), + abstol = 1e-9, maxiters = maxiters) sol_iip = solve!(solver) prob_oop = NonlinearProblem{false}(f, u0) - solver = init(prob_oop, TrustRegion(radius_update_scheme = radius_update_scheme), abstol = 1e-9, maxiters = maxiters) + solver = init(prob_oop, TrustRegion(radius_update_scheme = radius_update_scheme), + abstol = 1e-9, maxiters = maxiters) sol_oop = solve!(solver) return sol_iip.u[end], sol_oop.u[end] diff --git a/test/runtests.jl b/test/runtests.jl index 621900996..a84fc3cb1 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -13,12 +13,12 @@ end @time begin if GROUP == "All" || GROUP == "Core" - @time @safetestset "Basic Tests + Some AD" begin include("basictests.jl") end - @time @safetestset "Sparsity Tests" begin include("sparse.jl") end + @time @safetestset "Basic Tests + Some AD" include("basictests.jl") + @time @safetestset "Sparsity Tests" include("sparse.jl") end if GROUP == "GPU" activate_downstream_env() - @time @safetestset "GPU Tests" begin include("gpu.jl") end + @time @safetestset "GPU Tests" include("gpu.jl") end end From d0dcc4c51e4bec4558898b67ae3f8ce3d150d818 Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Tue, 25 Apr 2023 09:35:07 -0400 Subject: [PATCH 2/2] use NLSolveSafeTerminationResult --- src/raphson.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/raphson.jl b/src/raphson.jl index 08a156190..96b64e602 100644 --- a/src/raphson.jl +++ b/src/raphson.jl @@ -89,6 +89,7 @@ mutable struct NewtonRaphsonCache{iip, fType, algType, uType, duType, resType, p jac_config::JC, iter::Int, force_stop::Bool, maxiters::Int, internalnorm::INType, retcode::SciMLBase.ReturnCode.T, abstol::tolType, + reltol::tolType, termination_condition::TC, prob::probType) where { iip, fType, algType, uType, duType, resType, pType, INType,