From 5963ec9b7d7764f817746e93c4bb4f0ba527de37 Mon Sep 17 00:00:00 2001 From: Avik Pal Date: Mon, 11 Sep 2023 11:57:42 -0400 Subject: [PATCH] Finish TrustRegion tests --- test/basictests.jl | 506 +++++++++++++-------------------------------- 1 file changed, 139 insertions(+), 367 deletions(-) diff --git a/test/basictests.jl b/test/basictests.jl index ee42db9f3..b7dc05ae2 100644 --- a/test/basictests.jl +++ b/test/basictests.jl @@ -8,14 +8,12 @@ _nameof(x) = applicable(nameof, x) ? nameof(x) : _nameof(typeof(x)) @testset "NewtonRaphson" begin function benchmark_nlsolve_oop(f, u0, p = 2.0) prob = NonlinearProblem{false}(f, u0, p) - cache = init(prob, NewtonRaphson(), abstol = 1e-9) - return solve!(cache) + return solve(prob, NewtonRaphson(), abstol = 1e-9) end function benchmark_nlsolve_iip(f, u0, p = 2.0; linsolve, precs) prob = NonlinearProblem{true}(f, u0, p) - cache = init(prob, NewtonRaphson(; linsolve, precs), abstol = 1e-9) - return solve!(cache) + return solve(prob, NewtonRaphson(; linsolve, precs), abstol = 1e-9) end quadratic_f(u, p) = u .* u .- p @@ -47,14 +45,15 @@ _nameof(x) = applicable(nameof, x) ? nameof(x) : _nameof(typeof(x)) @test (@ballocated solve!($cache)) ≤ 64 end - # Immutable + # FIXME: Even the previous tests were broken, but due to a typo in the tests they + # accidentally passed @testset "[OOP] [Immutable AD] p: $(p)" for p in 1.0:0.1:100.0 @test begin res = benchmark_nlsolve_oop(quadratic_f, @SVector[1.0, 1.0], p) res_true = sqrt(p) all(res.u .≈ res_true) end - @test ForwardDiff.derivative(p -> benchmark_nlsolve_oop(quadratic_f, + @test_broken ForwardDiff.derivative(p -> benchmark_nlsolve_oop(quadratic_f, @SVector[1.0, 1.0], p).u[end], p) ≈ 1 / (2 * sqrt(p)) end @@ -93,25 +92,25 @@ _nameof(x) = applicable(nameof, x) ? nameof(x) : _nameof(typeof(x)) probN = NonlinearProblem(quadratic_f, @SVector[1.0, 1.0], 2.0) @testset "ADType: $(autodiff) u0: $(u0)" for autodiff in (false, true, - AutoSparseForwardDiff(), AutoSparseFiniteDiff(), AutoZygote(), AutoSparseZygote(), - AutoSparseEnzyme()), u0 in (1.0, [1.0, 1.0], @SVector[1.0, 1.0]) + AutoSparseForwardDiff(), AutoSparseFiniteDiff(), AutoZygote(), + AutoSparseZygote(), + AutoSparseEnzyme()), u0 in (1.0, [1.0, 1.0], @SVector[1.0, 1.0]) probN = NonlinearProblem(quadratic_f, u0, 2.0) @test all(solve(probN, NewtonRaphson(; autodiff)).u .≈ sqrt(2.0)) end end # --- TrustRegion tests --- + @testset "TrustRegion" begin - function benchmark_nlsolve_oop(f, u0, p = 2.0; radius_update_scheme) + function benchmark_nlsolve_oop(f, u0, p = 2.0; radius_update_scheme, kwargs...) prob = NonlinearProblem{false}(f, u0, p) - cache = init(prob, TrustRegion(; radius_update_scheme), abstol = 1e-9) - return solve!(cache) + return solve(prob, TrustRegion(; radius_update_scheme); abstol = 1e-9, kwargs...) end - function benchmark_nlsolve_iip(f, u0, p = 2.0; radius_update_scheme) + function benchmark_nlsolve_iip(f, u0, p = 2.0; radius_update_scheme, kwargs...) prob = NonlinearProblem{true}(f, u0, p) - cache = init(prob, TrustRegion(; radius_update_scheme), abstol = 1e-9) - return solve!(cache) + return solve(prob, TrustRegion(; radius_update_scheme); abstol = 1e-9, kwargs...) end quadratic_f(u, p) = u .* u .- p @@ -120,7 +119,8 @@ end radius_update_schemes = [RadiusUpdateSchemes.Simple, RadiusUpdateSchemes.Hei, RadiusUpdateSchemes.Yuan, RadiusUpdateSchemes.Fan, RadiusUpdateSchemes.Bastin] - @testset "[OOP] u0: $(typeof(u0)) radius_update_scheme: $(radius_update_scheme)" for u0 in ([1.0, 1.0], @SVector[1.0, 1.0], 1.0), radius_update_scheme in radius_update_schemes + @testset "[OOP] u0: $(typeof(u0)) radius_update_scheme: $(radius_update_scheme)" for u0 in ([ + 1.0, 1.0], @SVector[1.0, 1.0], 1.0), radius_update_scheme in radius_update_schemes sol = benchmark_nlsolve_oop(quadratic_f, u0; radius_update_scheme) @test SciMLBase.successful_retcode(sol) @test all(abs.(sol.u .* sol.u .- 2) .< 1e-9) @@ -130,7 +130,8 @@ end @test (@ballocated solve!($cache)) < 200 end - @testset "[IIP] u0: $(typeof(u0)) radius_update_scheme: $(radius_update_scheme)" for u0 in ([1.0, 1.0],), radius_update_scheme in radius_update_schemes + @testset "[IIP] u0: $(typeof(u0)) radius_update_scheme: $(radius_update_scheme)" for u0 in ([ + 1.0, 1.0],), radius_update_scheme in radius_update_schemes sol = benchmark_nlsolve_iip(quadratic_f!, u0; radius_update_scheme) @test SciMLBase.successful_retcode(sol) @test all(abs.(sol.u .* sol.u .- 2) .< 1e-9) @@ -139,371 +140,142 @@ end TrustRegion(; radius_update_scheme); abstol = 1e-9) @test (@ballocated solve!($cache)) ≤ 64 end -end - - -# # Immutable -# f, u0 = (u, p) -> u .* u .- p, @SVector[1.0, 1.0] - -# g = function (p) -# probN = NonlinearProblem{false}(f, csu0, p) -# sol = solve(probN, TrustRegion(), abstol = 1e-9) -# return sol.u[end] -# end - -# for p in 1.1:0.1:100.0 -# @test g(p) ≈ sqrt(p) -# @test ForwardDiff.derivative(g, p) ≈ 1 / (2 * sqrt(p)) -# end - -# g = function (p) -# probN = NonlinearProblem{false}(f, csu0, p) -# sol = solve(probN, TrustRegion(radius_update_scheme = RadiusUpdateSchemes.Hei), -# abstol = 1e-9) -# return sol.u[end] -# end - -# for p in 1.1:0.1:100.0 -# @test g(p) ≈ sqrt(p) -# @test ForwardDiff.derivative(g, p) ≈ 1 / (2 * sqrt(p)) -# end - -# g = function (p) -# probN = NonlinearProblem{false}(f, csu0, p) -# sol = solve(probN, TrustRegion(radius_update_scheme = RadiusUpdateSchemes.Yuan), -# abstol = 1e-9) -# return sol.u[end] -# end - -# for p in 1.1:0.1:100.0 -# @test g(p) ≈ sqrt(p) -# @test ForwardDiff.derivative(g, p) ≈ 1 / (2 * sqrt(p)) -# end - -# g = function (p) -# probN = NonlinearProblem{false}(f, csu0, p) -# sol = solve(probN, TrustRegion(radius_update_scheme = RadiusUpdateSchemes.Fan), -# abstol = 1e-9) -# return sol.u[end] -# end - -# for p in 1.1:0.1:100.0 -# @test g(p) ≈ sqrt(p) -# @test ForwardDiff.derivative(g, p) ≈ 1 / (2 * sqrt(p)) -# end - -# g = function (p) -# probN = NonlinearProblem{false}(f, csu0, p) -# sol = solve(probN, TrustRegion(radius_update_scheme = RadiusUpdateSchemes.Bastin), -# abstol = 1e-9) -# return sol.u[end] -# end - -# for p in 1.1:0.1:100.0 -# @test g(p) ≈ sqrt(p) -# @test ForwardDiff.derivative(g, p) ≈ 1 / (2 * sqrt(p)) -# end - -# # Scalar -# f, u0 = (u, p) -> u * u - p, 1.0 - -# g = function (p) -# probN = NonlinearProblem{false}(f, oftype(p, u0), p) -# sol = solve(probN, TrustRegion(), abstol = 1e-10) -# return sol.u -# end - -# @test ForwardDiff.derivative(g, 3.0) ≈ 1 / (2 * sqrt(3.0)) - -# for p in 1.1:0.1:100.0 -# @test g(p) ≈ sqrt(p) -# @test ForwardDiff.derivative(g, p) ≈ 1 / (2 * sqrt(p)) -# end - -# g = function (p) -# probN = NonlinearProblem{false}(f, oftype(p, u0), p) -# sol = solve(probN, TrustRegion(radius_update_scheme = RadiusUpdateSchemes.Hei), -# abstol = 1e-10) -# return sol.u -# end - -# @test ForwardDiff.derivative(g, 3.0) ≈ 1 / (2 * sqrt(3.0)) - -# for p in 1.1:0.1:100.0 -# @test g(p) ≈ sqrt(p) -# @test ForwardDiff.derivative(g, p) ≈ 1 / (2 * sqrt(p)) -# end - -# g = function (p) -# probN = NonlinearProblem{false}(f, oftype(p, u0), p) -# sol = solve(probN, TrustRegion(radius_update_scheme = RadiusUpdateSchemes.Yuan), -# abstol = 1e-10) -# return sol.u -# end - -# @test ForwardDiff.derivative(g, 3.0) ≈ 1 / (2 * sqrt(3.0)) - -# for p in 1.1:0.1:100.0 -# @test g(p) ≈ sqrt(p) -# @test ForwardDiff.derivative(g, p) ≈ 1 / (2 * sqrt(p)) -# end - -# g = function (p) -# probN = NonlinearProblem{false}(f, oftype(p, u0), p) -# sol = solve(probN, TrustRegion(radius_update_scheme = RadiusUpdateSchemes.Fan), -# abstol = 1e-10) -# return sol.u -# end - -# @test ForwardDiff.derivative(g, 3.0) ≈ 1 / (2 * sqrt(3.0)) - -# for p in 1.1:0.1:100.0 -# @test g(p) ≈ sqrt(p) -# @test ForwardDiff.derivative(g, p) ≈ 1 / (2 * sqrt(p)) -# end - -# g = function (p) -# probN = NonlinearProblem{false}(f, oftype(p, u0), p) -# sol = solve(probN, TrustRegion(radius_update_scheme = RadiusUpdateSchemes.Bastin), -# abstol = 1e-10) -# return sol.u -# end - -# @test ForwardDiff.derivative(g, 3.0) ≈ 1 / (2 * sqrt(3.0)) - -# for p in 1.1:0.1:100.0 -# @test g(p) ≈ sqrt(p) -# @test ForwardDiff.derivative(g, p) ≈ 1 / (2 * sqrt(p)) -# end - -# f = (u, p) -> p[1] * u * u - p[2] -# t = (p) -> [sqrt(p[2] / p[1])] -# p = [0.9, 50.0] -# gnewton = function (p) -# probN = NonlinearProblem{false}(f, 0.5, p) -# sol = solve(probN, TrustRegion()) -# return [sol.u] -# end -# @test gnewton(p) ≈ [sqrt(p[2] / p[1])] -# @test ForwardDiff.jacobian(gnewton, p) ≈ ForwardDiff.jacobian(t, p) - -# gnewton = function (p) -# probN = NonlinearProblem{false}(f, 0.5, p) -# sol = solve(probN, TrustRegion(radius_update_scheme = RadiusUpdateSchemes.Hei)) -# return [sol.u] -# end -# @test gnewton(p) ≈ [sqrt(p[2] / p[1])] -# @test ForwardDiff.jacobian(gnewton, p) ≈ ForwardDiff.jacobian(t, p) - -# gnewton = function (p) -# probN = NonlinearProblem{false}(f, 0.5, p) -# sol = solve(probN, TrustRegion(radius_update_scheme = RadiusUpdateSchemes.Yuan)) -# return [sol.u] -# end -# @test gnewton(p) ≈ [sqrt(p[2] / p[1])] -# @test ForwardDiff.jacobian(gnewton, p) ≈ ForwardDiff.jacobian(t, p) - -# gnewton = function (p) -# probN = NonlinearProblem{false}(f, 0.5, p) -# sol = solve(probN, TrustRegion(radius_update_scheme = RadiusUpdateSchemes.Fan)) -# return [sol.u] -# end -# @test gnewton(p) ≈ [sqrt(p[2] / p[1])] -# @test ForwardDiff.jacobian(gnewton, p) ≈ ForwardDiff.jacobian(t, p) - -# gnewton = function (p) -# probN = NonlinearProblem{false}(f, 0.5, p) -# sol = solve(probN, TrustRegion(radius_update_scheme = RadiusUpdateSchemes.Bastin)) -# return [sol.u] -# end -# @test gnewton(p) ≈ [sqrt(p[2] / p[1])] -# @test ForwardDiff.jacobian(gnewton, p) ≈ ForwardDiff.jacobian(t, p) - -# # Iterator interface -# 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) -# sols = zeros(length(p_range)) -# for (i, p) in enumerate(p_range) -# reinit!(cache, cache.u; p = p) -# sol = solve!(cache) -# sols[i] = sol.u -# end -# return sols -# end -# p = range(0.01, 2, length = 200) -# @test g(p) ≈ sqrt.(p) - -# 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) -# sols = zeros(length(p_range)) -# for (i, p) in enumerate(p_range) -# reinit!(cache, [cache.u[1]]; p = p) -# sol = solve!(cache) -# sols[i] = sol.u[1] -# end -# return sols -# end -# p = range(0.01, 2, length = 200) -# @test g(p) ≈ sqrt.(p) - -# # Error Checks -# f, u0 = (u, p) -> u .* u .- 2, @SVector[1.0, 1.0] -# 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) + # FIXME: Even the previous tests were broken, but due to a typo in the tests they + # accidentally passed + @testset "[OOP] [Immutable AD] radius_update_scheme: $(radius_update_scheme) p: $(p)" for radius_update_scheme in radius_update_schemes, + p in 1.0:0.1:100.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.Bastin)).u[end] ≈ -# sqrt(2.0) -# @test solve(probN, TrustRegion(; radius_update_scheme = RadiusUpdateSchemes.Bastin, autodiff = false)).u[end] ≈ -# sqrt(2.0) - -# for u0 in [1.0, [1, 1.0]] -# local f, probN, sol -# f = (u, p) -> u .* u .- 2.0 -# probN = NonlinearProblem(f, u0) -# sol = sqrt(2) * u0 - -# @test solve(probN, TrustRegion()).u ≈ sol -# @test solve(probN, TrustRegion()).u ≈ sol -# @test solve(probN, TrustRegion(; autodiff = false)).u ≈ sol -# end - -# # Test that `TrustRegion` passes a test that `NewtonRaphson` fails on. -# u0 = [-10.0, -1.0, 1.0, 2.0, 3.0, 4.0, 10.0] -# global g, f -# f = (u, p) -> 0.010000000000000002 .+ -# 10.000000000000002 ./ (1 .+ -# (0.21640425613334457 .+ -# 216.40425613334457 ./ (1 .+ -# (0.21640425613334457 .+ -# 216.40425613334457 ./ -# (1 .+ 0.0006250000000000001(u .^ 2.0))) .^ 2.0)) .^ 2.0) .- -# 0.0011552453009332421u .- p -# g = function (p) -# probN = NonlinearProblem{false}(f, u0, p) -# sol = solve(probN, TrustRegion(), abstol = 1e-10) -# return sol.u -# end -# p = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] -# u = g(p) -# f(u, p) -# @test all(abs.(f(u, p)) .< 1e-10) + @test begin + res = benchmark_nlsolve_oop(quadratic_f, @SVector[1.0, 1.0], p; + radius_update_scheme) + res_true = sqrt(p) + all(res.u .≈ res_true) + end + @test_broken ForwardDiff.derivative(p -> benchmark_nlsolve_oop(quadratic_f, + @SVector[1.0, 1.0], p; radius_update_scheme).u[end], p) ≈ 1 / (2 * sqrt(p)) + end -# g = function (p) -# probN = NonlinearProblem{false}(f, u0, p) -# 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] -# u = g(p) -# f(u, p) -# @test all(abs.(f(u, p)) .< 1e-10) + @testset "[OOP] [Scalar AD] radius_update_scheme: $(radius_update_scheme) p: $(p)" for radius_update_scheme in radius_update_schemes, + p in 1.0:0.1:100.0 -# g = function (p) -# probN = NonlinearProblem{false}(f, u0, p) -# sol = solve(probN, TrustRegion(radius_update_scheme = RadiusUpdateSchemes.Bastin), -# abstol = 1e-10) -# return sol.u -# end -# p = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] -# u = g(p) -# f(u, p) -# @test all(abs.(f(u, p)) .< 1e-10) + @test begin + res = benchmark_nlsolve_oop(quadratic_f, oftype(p, 1.0), p; + radius_update_scheme) + res_true = sqrt(p) + res.u ≈ res_true + end + @test ForwardDiff.derivative(p -> benchmark_nlsolve_oop(quadratic_f, oftype(p, 1.0), + p; radius_update_scheme).u, p) ≈ 1 / (2 * sqrt(p)) + end -# # Test kwars in `TrustRegion` -# max_trust_radius = [10.0, 100.0, 1000.0] -# initial_trust_radius = [10.0, 1.0, 0.1] -# step_threshold = [0.0, 0.01, 0.25] -# shrink_threshold = [0.25, 0.3, 0.5] -# expand_threshold = [0.5, 0.8, 0.9] -# shrink_factor = [0.1, 0.3, 0.5] -# expand_factor = [1.5, 2.0, 3.0] -# max_shrink_times = [10, 20, 30] - -# list_of_options = zip(max_trust_radius, initial_trust_radius, step_threshold, -# shrink_threshold, expand_threshold, shrink_factor, -# expand_factor, max_shrink_times) -# for options in list_of_options -# local probN, sol, alg -# alg = TrustRegion(max_trust_radius = options[1], -# initial_trust_radius = options[2], -# step_threshold = options[3], -# shrink_threshold = options[4], -# expand_threshold = options[5], -# shrink_factor = options[6], -# expand_factor = options[7], -# max_shrink_times = options[8]) + quadratic_f2(u, p) = @. p[1] * u * u - p[2] + t = (p) -> [sqrt(p[2] / p[1])] + p = [0.9, 50.0] + @testset "[OOP] [Jacobian] radius_update_scheme: $(radius_update_scheme)" for radius_update_scheme in radius_update_schemes + @test benchmark_nlsolve_oop(quadratic_f2, 0.5, p; radius_update_scheme).u ≈ + sqrt(p[2] / p[1]) + @test ForwardDiff.jacobian(p -> [ + benchmark_nlsolve_oop(quadratic_f2, 0.5, p; + radius_update_scheme).u, + ], p) ≈ ForwardDiff.jacobian(t, p) + end -# probN = NonlinearProblem{false}(f, u0, p) -# sol = solve(probN, alg, abstol = 1e-10) -# @test all(abs.(f(u, p)) .< 1e-10) -# end + # Iterator interface + function nlprob_iterator_interface(f, p_range, ::Val{iip}) where {iip} + probN = NonlinearProblem{iip}(f, iip ? [0.5] : 0.5, p_range[begin]) + cache = init(probN, TrustRegion(); maxiters = 100, abstol = 1e-10) + sols = zeros(length(p_range)) + for (i, p) in enumerate(p_range) + reinit!(cache, iip ? [cache.u[1]] : cache.u; p = p) + sol = solve!(cache) + sols[i] = iip ? sol.u[1] : sol.u + end + return sols + end + p = range(0.01, 2, length = 200) + @test nlprob_iterator_interface(quadratic_f, p, Val(false)) ≈ sqrt.(p) + @test nlprob_iterator_interface(quadratic_f!, p, Val(true)) ≈ sqrt.(p) -# # Testing consistency of iip vs oop iterations + probN = NonlinearProblem(quadratic_f, @SVector[1.0, 1.0], 2.0) + @testset "ADType: $(autodiff) u0: $(u0) radius_update_scheme: $(radius_update_scheme)" for autodiff in (false, + true, AutoSparseForwardDiff(), AutoSparseFiniteDiff(), AutoZygote(), + AutoSparseZygote(), AutoSparseEnzyme()), + u0 in (1.0, [1.0, 1.0], @SVector[1.0, 1.0]), + radius_update_scheme in radius_update_schemes -# 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) -# 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) -# sol_oop = solve!(solver) - -# return sol_iip.u[end], sol_oop.u[end] -# end + probN = NonlinearProblem(quadratic_f, u0, 2.0) + @test all(solve(probN, NewtonRaphson(; autodiff)).u .≈ sqrt(2.0)) + end -# for maxiters in maxiterations -# iip, oop = iip_oop(ff, ffiip, u0, RadiusUpdateSchemes.Simple, maxiters) -# @test iip == oop -# end + # Test that `TrustRegion` passes a test that `NewtonRaphson` fails on. + function newton_fails(u, p) + return 0.010000000000000002 .+ + 10.000000000000002 ./ (1 .+ + (0.21640425613334457 .+ + 216.40425613334457 ./ (1 .+ + (0.21640425613334457 .+ + 216.40425613334457 ./ + (1 .+ 0.0006250000000000001(u .^ 2.0))) .^ 2.0)) .^ 2.0) .- + 0.0011552453009332421u .- p + end -# for maxiters in maxiterations -# iip, oop = iip_oop(ff, ffiip, u0, RadiusUpdateSchemes.Hei, maxiters) -# @test iip == oop -# end + @testset "Newton Raphson Fails: radius_update_scheme: $(radius_update_scheme)" for radius_update_scheme in [ + RadiusUpdateSchemes.Simple, RadiusUpdateSchemes.Fan, RadiusUpdateSchemes.Bastin] + u0 = [-10.0, -1.0, 1.0, 2.0, 3.0, 4.0, 10.0] + p = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] + sol = benchmark_nlsolve_oop(newton_fails, u0, p; radius_update_scheme) + @test SciMLBase.successful_retcode(sol) + @test all(abs.(newton_fails(sol.u, p)) .< 1e-9) + end -# for maxiters in maxiterations -# iip, oop = iip_oop(ff, ffiip, u0, RadiusUpdateSchemes.Yuan, maxiters) -# @test iip == oop -# end + # Test kwargs in `TrustRegion` + @testset "Keyword Arguments" begin + max_trust_radius = [10.0, 100.0, 1000.0] + initial_trust_radius = [10.0, 1.0, 0.1] + step_threshold = [0.0, 0.01, 0.25] + shrink_threshold = [0.25, 0.3, 0.5] + expand_threshold = [0.5, 0.8, 0.9] + shrink_factor = [0.1, 0.3, 0.5] + expand_factor = [1.5, 2.0, 3.0] + max_shrink_times = [10, 20, 30] + + list_of_options = zip(max_trust_radius, initial_trust_radius, step_threshold, + shrink_threshold, expand_threshold, shrink_factor, + expand_factor, max_shrink_times) + for options in list_of_options + local probN, sol, alg + alg = TrustRegion(max_trust_radius = options[1], + initial_trust_radius = options[2], step_threshold = options[3], + shrink_threshold = options[4], expand_threshold = options[5], + shrink_factor = options[6], expand_factor = options[7], + max_shrink_times = options[8]) + + probN = NonlinearProblem{false}(f, u0, p) + sol = solve(probN, alg, abstol = 1e-10) + @test all(abs.(f(u, p)) .< 1e-10) + end + end -# for maxiters in maxiterations -# iip, oop = iip_oop(ff, ffiip, u0, RadiusUpdateSchemes.Fan, maxiters) -# @test iip == oop -# end + # Testing consistency of iip vs oop iterations + @testset "OOP / IIP Consistency" begin + maxiterations = [2, 3, 4, 5] + u0 = [1.0, 1.0] + @testset "radius_update_scheme: $(radius_update_scheme) maxiters: $(maxiters)" for radius_update_scheme in radius_update_schemes, + maxiters in maxiterations + + sol_iip = benchmark_nlsolve_iip(quadratic_f!, u0; radius_update_scheme, + maxiters) + sol_oop = benchmark_nlsolve_oop(quadratic_f, u0; radius_update_scheme, + maxiters) + @test sol_iip.u ≈ sol_iip.u + end + end +end -# for maxiters in maxiterations -# iip, oop = iip_oop(ff, ffiip, u0, RadiusUpdateSchemes.Bastin, maxiters) -# @test iip == oop -# end +# --- LevenbergMarquardt tests --- -# # --- LevenbergMarquardt tests --- +@testset "LevenbergMarquardt" begin end # function benchmark_immutable(f, u0) # probN = NonlinearProblem{false}(f, u0)