From 2cfa3f19c4c7e9823af4f483063e67aad07c51a5 Mon Sep 17 00:00:00 2001 From: Yash Raj Singh Date: Fri, 24 Mar 2023 04:13:12 +0100 Subject: [PATCH 01/11] jvp for yuan's method --- src/jacobian.jl | 11 +++++++++++ src/trustRegion.jl | 9 +++++---- 2 files changed, 16 insertions(+), 4 deletions(-) diff --git a/src/jacobian.jl b/src/jacobian.jl index 7661d8e55..555e10a7b 100644 --- a/src/jacobian.jl +++ b/src/jacobian.jl @@ -144,3 +144,14 @@ function jacobian_autodiff(f, x::AbstractArray, nonlinfun, alg) jac_prototype = jac_prototype, chunksize = chunk_size), num_of_chunks) end + +function jvp(cache::TrustRegionCache{false}) + @unpack f, u, fu = cache + auto_jacvec(f, u, fu) +end + +function jvp(cache::TrustRegionCache{true}) + @unpack g, f, u, fu = cache + auto_jacvec!(g, f, u, fu) + g +end diff --git a/src/trustRegion.jl b/src/trustRegion.jl index fc1662fc0..3e0df5ffa 100644 --- a/src/trustRegion.jl +++ b/src/trustRegion.jl @@ -417,9 +417,7 @@ function trust_region_step!(cache::TrustRegionCache) elseif r >= cache.expand_threshold && cache.internalnorm(step_size) > cache.trust_r / 2 cache.p1 = cache.p3 * cache.p1 end - @unpack p1, fu, f, J = cache - #cache.trust_r = p1 * cache.internalnorm(jacobian!(J, cache) * fu) # we need the gradient at the new (k+1)th point WILL THIS BECOME ALLOCATING? - + if r > cache.step_threshold take_step!(cache) cache.loss = cache.loss_new @@ -427,7 +425,10 @@ function trust_region_step!(cache::TrustRegionCache) else cache.make_new_J = false end - + + @unpack p1= cache + cache.trust_r = p1 * cache.internalnorm(jvp(cache)) # we need the gradient at the new (k+1)th point WILL THIS BECOME ALLOCATING? + if iszero(cache.fu) || cache.internalnorm(cache.fu) < cache.abstol || cache.internalnorm(g) < cache.ϵ # parameters to be defined cache.force_stop = true end From 47fef1e4d6a91067ccfc8bf599209f9886f77f5c Mon Sep 17 00:00:00 2001 From: Yash Raj Singh Date: Sat, 25 Mar 2023 03:40:16 +0100 Subject: [PATCH 02/11] jvp modifications --- src/NonlinearSolve.jl | 1 + src/jacobian.jl | 9 --------- src/trustRegion.jl | 20 +++++++++++++++----- 3 files changed, 16 insertions(+), 14 deletions(-) diff --git a/src/NonlinearSolve.jl b/src/NonlinearSolve.jl index d6b538059..0eac4b767 100644 --- a/src/NonlinearSolve.jl +++ b/src/NonlinearSolve.jl @@ -61,6 +61,7 @@ SnoopPrecompile.@precompile_all_calls begin for T in (Float32, Float64) end end export RadiusUpdateSchemes +export auto_jacvec, auto_jacvec! export NewtonRaphson, TrustRegion, LevenbergMarquardt diff --git a/src/jacobian.jl b/src/jacobian.jl index 555e10a7b..9d8e3bfc5 100644 --- a/src/jacobian.jl +++ b/src/jacobian.jl @@ -145,13 +145,4 @@ function jacobian_autodiff(f, x::AbstractArray, nonlinfun, alg) num_of_chunks) end -function jvp(cache::TrustRegionCache{false}) - @unpack f, u, fu = cache - auto_jacvec(f, u, fu) -end -function jvp(cache::TrustRegionCache{true}) - @unpack g, f, u, fu = cache - auto_jacvec!(g, f, u, fu) - g -end diff --git a/src/trustRegion.jl b/src/trustRegion.jl index 3e0df5ffa..730b9db5e 100644 --- a/src/trustRegion.jl +++ b/src/trustRegion.jl @@ -403,9 +403,9 @@ function trust_region_step!(cache::TrustRegionCache) if rfunc(r, shrink_threshold, p1, p3, p4, p2) * cache.internalnorm(step_size) < cache.trust_r cache.shrink_counter += 1 end - cache.trust_r = rfunc(r, shrink_threshold, p1, p3, p4, p2) * cache.internalnorm(step_size) # parameters to be defined + 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.ϵ # parameters to be defined + if iszero(cache.fu) || cache.internalnorm(cache.fu) < cache.abstol || cache.internalnorm(g) < cache.ϵ cache.force_stop = true end @@ -427,9 +427,8 @@ function trust_region_step!(cache::TrustRegionCache) end @unpack p1= cache - cache.trust_r = p1 * cache.internalnorm(jvp(cache)) # we need the gradient at the new (k+1)th point WILL THIS BECOME ALLOCATING? - - if iszero(cache.fu) || cache.internalnorm(cache.fu) < cache.abstol || cache.internalnorm(g) < cache.ϵ # parameters to be defined + 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 @@ -474,6 +473,17 @@ function take_step!(cache::TrustRegionCache{false}) cache.fu = cache.fu_new end +function jvp(cache::TrustRegionCache{false}) + @unpack f, u, fu = cache + auto_jacvec(f, u, fu) +end + +function jvp(cache::TrustRegionCache{true}) + @unpack g, f, u, fu = cache + auto_jacvec!(g, f, u, fu) + g +end + function SciMLBase.solve!(cache::TrustRegionCache) while !cache.force_stop && cache.iter < cache.maxiters && cache.shrink_counter < cache.alg.max_shrink_times From 12b4655a4dc72856426e053a24a9bf14f6f525f6 Mon Sep 17 00:00:00 2001 From: Yash Raj Singh Date: Sun, 26 Mar 2023 05:09:35 +0200 Subject: [PATCH 03/11] fixes --- src/NonlinearSolve.jl | 1 - src/trustRegion.jl | 2 +- 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/src/NonlinearSolve.jl b/src/NonlinearSolve.jl index 0eac4b767..d6b538059 100644 --- a/src/NonlinearSolve.jl +++ b/src/NonlinearSolve.jl @@ -61,7 +61,6 @@ SnoopPrecompile.@precompile_all_calls begin for T in (Float32, Float64) end end export RadiusUpdateSchemes -export auto_jacvec, auto_jacvec! export NewtonRaphson, TrustRegion, LevenbergMarquardt diff --git a/src/trustRegion.jl b/src/trustRegion.jl index 730b9db5e..280f54b31 100644 --- a/src/trustRegion.jl +++ b/src/trustRegion.jl @@ -301,7 +301,7 @@ function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg::TrustRegion, p2 = convert(eltype(u), 1/6) # c5 p3 = convert(eltype(u), 6.0) # c6 p4 = convert(eltype(u), 0.0) - end + end return TrustRegionCache{iip}(f, alg, u, fu, p, uf, linsolve, J, jac_config, 1, false, maxiters, internalnorm, From e251104f671fce6b787d0291d618e35cefd05446 Mon Sep 17 00:00:00 2001 From: Yash Raj Singh Date: Sun, 26 Mar 2023 05:34:37 +0200 Subject: [PATCH 04/11] added some simple tests --- test/basictests.jl | 39 ++++++++++++++++++++++----------------- 1 file changed, 22 insertions(+), 17 deletions(-) diff --git a/test/basictests.jl b/test/basictests.jl index 947975959..7c6e38920 100644 --- a/test/basictests.jl +++ b/test/basictests.jl @@ -132,41 +132,46 @@ end # --- TrustRegion tests --- -function benchmark_immutable(f, u0) +function benchmark_immutable(f, u0, radius_update_scheme) probN = NonlinearProblem{false}(f, u0) - solver = init(probN, TrustRegion(), 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) +function benchmark_mutable(f, u0, radius_update_scheme) probN = NonlinearProblem{false}(f, u0) - solver = init(probN, TrustRegion(), 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) +function benchmark_scalar(f, u0, radius_update_scheme) probN = NonlinearProblem{false}(f, u0) - sol = (solve(probN, TrustRegion(), abstol = 1e-9)) + sol = (solve(probN, TrustRegion(radius_update_scheme = radius_update_scheme), abstol = 1e-9)) end -function ff(u, p) +function ff(u, p=nothing) u .* u .- 2 end -function sf(u, p) +function sf(u, p=nothing) u * u - 2 end u0 = [1.0, 1.0] -sol = benchmark_immutable(ff, cu0) -@test sol.retcode === ReturnCode.Success -@test all(abs.(sol.u .* sol.u .- 2) .< 1e-9) -sol = benchmark_mutable(ff, u0) -@test sol.retcode === ReturnCode.Success -@test all(abs.(sol.u .* sol.u .- 2) .< 1e-9) -sol = benchmark_scalar(sf, csu0) -@test sol.retcode === ReturnCode.Success -@test abs(sol.u * sol.u - 2) < 1e-9 +radius_update_schemes = [RadiusUpdateSchemes.Simple, RadiusUpdateSchemes.Hei, RadiusUpdateSchemes.Yuan] + +for radius_update_scheme in radius_update_schemes + sol = benchmark_immutable(ff, cu0, radius_update_scheme) + @test sol.retcode === ReturnCode.Success + @test all(abs.(sol.u .* sol.u .- 2) .< 1e-9) + sol = benchmark_mutable(ff, u0, radius_update_scheme) + @test sol.retcode === ReturnCode.Success + @test all(abs.(sol.u .* sol.u .- 2) .< 1e-9) + sol = benchmark_scalar(sf, csu0, radius_update_scheme) + @test sol.retcode === ReturnCode.Success + @test abs(sol.u * sol.u - 2) < 1e-9 +end + function benchmark_inplace(f, u0) probN = NonlinearProblem{true}(f, u0) From 7d0e31d83a5ad6bbb7bed8f81ef4f04851cbcbdb Mon Sep 17 00:00:00 2001 From: Yash Raj Singh Date: Mon, 27 Mar 2023 00:52:51 +0200 Subject: [PATCH 05/11] more tests and jvp modifications --- src/trustRegion.jl | 12 +++++++++--- test/basictests.jl | 27 +++++++++++++++------------ 2 files changed, 24 insertions(+), 15 deletions(-) diff --git a/src/trustRegion.jl b/src/trustRegion.jl index 280f54b31..ca13b5525 100644 --- a/src/trustRegion.jl +++ b/src/trustRegion.jl @@ -428,7 +428,7 @@ function trust_region_step!(cache::TrustRegionCache) @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.ϵ + if iszero(cache.fu) || cache.internalnorm(cache.fu) < cache.abstol # || cache.internalnorm(g) < cache.ϵ cache.force_stop = true end @@ -475,12 +475,18 @@ end function jvp(cache::TrustRegionCache{false}) @unpack f, u, fu = cache - auto_jacvec(f, u, fu) + if isa(u, Number) + return value_derivative(f, u) + end + return auto_jacvec(f, u, fu) end function jvp(cache::TrustRegionCache{true}) @unpack g, f, u, fu = cache - auto_jacvec!(g, f, u, fu) + if isa(u, Number) + return value_derivative(f, u) + end + return auto_jacvec!(g, f, u, fu) g end diff --git a/test/basictests.jl b/test/basictests.jl index 7c6e38920..7f8fb2bb8 100644 --- a/test/basictests.jl +++ b/test/basictests.jl @@ -156,9 +156,9 @@ end function sf(u, p=nothing) u * u - 2 end -u0 = [1.0, 1.0] -radius_update_schemes = [RadiusUpdateSchemes.Simple, RadiusUpdateSchemes.Hei, RadiusUpdateSchemes.Yuan] +u0 = [1.0, 1.0] +radius_update_schemes = [RadiusUpdateSchemes.Simple, RadiusUpdateSchemes.Hei] for radius_update_scheme in radius_update_schemes sol = benchmark_immutable(ff, cu0, radius_update_scheme) @@ -173,25 +173,28 @@ for radius_update_scheme in radius_update_schemes end -function benchmark_inplace(f, u0) +function benchmark_inplace(f, u0, radius_update_scheme) probN = NonlinearProblem{true}(f, u0) - solver = init(probN, TrustRegion(), 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) +function ffiip(du, u, p=nothing) du .= u .* u .- 2 end u0 = [1.0, 1.0] -sol = benchmark_inplace(ffiip, u0) -@test sol.retcode === ReturnCode.Success -@test all(abs.(sol.u .* sol.u .- 2) .< 1e-9) +for radius_update_scheme in radius_update_schemes + sol = benchmark_inplace(ffiip, u0, radius_update_scheme) + @test sol.retcode === ReturnCode.Success + @test all(abs.(sol.u .* sol.u .- 2) .< 1e-9) +end -u0 = [1.0, 1.0] -probN = NonlinearProblem{true}(ffiip, u0) -solver = init(probN, TrustRegion(), abstol = 1e-9) -@test (@ballocated solve!(solver)) < 200 +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) + @test (@ballocated solve!(solver)) < 200 +end # AD Tests using ForwardDiff From df859fff396966c001a435e52d6450e838c297aa Mon Sep 17 00:00:00 2001 From: Yash Raj Singh Date: Wed, 29 Mar 2023 23:23:27 +0200 Subject: [PATCH 06/11] more tests --- test/basictests.jl | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/test/basictests.jl b/test/basictests.jl index 7f8fb2bb8..2f47dead4 100644 --- a/test/basictests.jl +++ b/test/basictests.jl @@ -213,6 +213,28 @@ for p in 1.1:0.1:100.0 @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 + # Scalar f, u0 = (u, p) -> u * u - p, 1.0 From 08e34dc6de5e3f68c5284ad327ddf1e99f37b645 Mon Sep 17 00:00:00 2001 From: Yash Raj Singh Date: Wed, 29 Mar 2023 23:45:57 +0200 Subject: [PATCH 07/11] fan's method --- src/trustRegion.jl | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/trustRegion.jl b/src/trustRegion.jl index ca13b5525..39f78a1b1 100644 --- a/src/trustRegion.jl +++ b/src/trustRegion.jl @@ -85,6 +85,7 @@ EnumX.@enumx RadiusUpdateSchemes begin Hei Yuan Bastin + Fan end struct TrustRegion{CS, AD, FDT, L, P, ST, CJ, MTR} <: @@ -402,6 +403,8 @@ function trust_region_step!(cache::TrustRegionCache) @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) @@ -416,6 +419,7 @@ function trust_region_step!(cache::TrustRegionCache) 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 @@ -428,7 +432,7 @@ function trust_region_step!(cache::TrustRegionCache) @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.ϵ + if iszero(cache.fu) || cache.internalnorm(cache.fu) < cache.abstol || cache.internalnorm(g) < cache.ϵ cache.force_stop = true end From 230ccb1656f40518c97c2b9cdbe9bb9fe2908a21 Mon Sep 17 00:00:00 2001 From: Yash Raj Singh Date: Fri, 31 Mar 2023 05:37:31 +0200 Subject: [PATCH 08/11] iip issue resolved --- src/trustRegion.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/trustRegion.jl b/src/trustRegion.jl index 39f78a1b1..df7daec53 100644 --- a/src/trustRegion.jl +++ b/src/trustRegion.jl @@ -445,7 +445,7 @@ function dogleg!(cache::TrustRegionCache) # Test if the full step is within the trust region. if norm(u_tmp) ≤ trust_r - cache.step_size = u_tmp + cache.step_size = deepcopy(u_tmp) return end From e76281b36af084637111d47c8b226ac9630807a7 Mon Sep 17 00:00:00 2001 From: Yash Raj Singh Date: Fri, 31 Mar 2023 05:37:51 +0200 Subject: [PATCH 09/11] more tests --- test/basictests.jl | 59 +++++++++++++++++++++++++++++++++++++++------- 1 file changed, 50 insertions(+), 9 deletions(-) diff --git a/test/basictests.jl b/test/basictests.jl index 2f47dead4..ffa42d7a5 100644 --- a/test/basictests.jl +++ b/test/basictests.jl @@ -158,7 +158,7 @@ function sf(u, p=nothing) end u0 = [1.0, 1.0] -radius_update_schemes = [RadiusUpdateSchemes.Simple, RadiusUpdateSchemes.Hei] +radius_update_schemes = [RadiusUpdateSchemes.Simple, RadiusUpdateSchemes.Hei, RadiusUpdateSchemes.Yuan] for radius_update_scheme in radius_update_schemes sol = benchmark_immutable(ff, cu0, radius_update_scheme) @@ -224,23 +224,37 @@ for p in 1.1:0.1:100.0 @test ForwardDiff.derivative(g, p) ≈ 1 / (2 * sqrt(p)) end +## FAIL BECAUSE JVP CANNOT ACCEPT PARAMETERS IN FUNCTIONS +# 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 + +# Scalar +f, u0 = (u, p) -> u * u - p, 1.0 + 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] + 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 -# 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) + sol = solve(probN, TrustRegion(radius_update_scheme = RadiusUpdateSchemes.Hei), abstol = 1e-10) return sol.u end @@ -251,6 +265,19 @@ for p in 1.1:0.1:100.0 @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 + f = (u, p) -> p[1] * u * u - p[2] t = (p) -> [sqrt(p[2] / p[1])] p = [0.9, 50.0] @@ -262,13 +289,27 @@ 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) + # Error Checks -f, u0 = (u, p) -> u .* u .- 2.0, @SVector[1.0, 1.0] +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) + +@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) + for u0 in [1.0, [1, 1.0]] local f, probN, sol f = (u, p) -> u .* u .- 2.0 From d80286dcb7960b11bec64740ac8593fde8fd1b36 Mon Sep 17 00:00:00 2001 From: Yash Raj Singh Date: Sun, 2 Apr 2023 04:20:56 +0200 Subject: [PATCH 10/11] fixed parameter issue in jvp --- src/trustRegion.jl | 18 +++++++++--------- test/basictests.jl | 38 +++++++++++++++++++------------------- 2 files changed, 28 insertions(+), 28 deletions(-) diff --git a/src/trustRegion.jl b/src/trustRegion.jl index df7daec53..1feea3f72 100644 --- a/src/trustRegion.jl +++ b/src/trustRegion.jl @@ -431,7 +431,7 @@ function trust_region_step!(cache::TrustRegionCache) end @unpack p1= cache - cache.trust_r = p1 * cache.internalnorm(jvp(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 @@ -477,20 +477,20 @@ function take_step!(cache::TrustRegionCache{false}) cache.fu = cache.fu_new end -function jvp(cache::TrustRegionCache{false}) - @unpack f, u, fu = cache +function jvp!(cache::TrustRegionCache{false}) + @unpack f, u, fu, p = cache if isa(u, Number) - return value_derivative(f, u) + return value_derivative(x -> f(x, p), u) end - return auto_jacvec(f, u, fu) + return auto_jacvec(x -> f(x, p), u, fu) end -function jvp(cache::TrustRegionCache{true}) - @unpack g, f, u, fu = cache +function jvp!(cache::TrustRegionCache{true}) + @unpack g, f, u, fu, p = cache if isa(u, Number) - return value_derivative(f, u) + return value_derivative(x -> f(x, p), u) end - return auto_jacvec!(g, f, u, fu) + return auto_jacvec!(g, (fu, x) -> f(fu, x, p), u, fu) g end diff --git a/test/basictests.jl b/test/basictests.jl index ffa42d7a5..e1e15fade 100644 --- a/test/basictests.jl +++ b/test/basictests.jl @@ -225,16 +225,16 @@ for p in 1.1:0.1:100.0 end ## FAIL BECAUSE JVP CANNOT ACCEPT PARAMETERS IN FUNCTIONS -# 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 +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 +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 @@ -265,18 +265,18 @@ for p in 1.1:0.1:100.0 @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 +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)) +@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 +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])] From 50f266d7ade9a1f79cf80be022c0c63f69e58fd8 Mon Sep 17 00:00:00 2001 From: Yash Raj Singh <34894852+yash2798@users.noreply.github.com> Date: Sun, 2 Apr 2023 19:31:25 +0200 Subject: [PATCH 11/11] tolerance decreased --- src/trustRegion.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/trustRegion.jl b/src/trustRegion.jl index 693396bf2..3e157f11a 100644 --- a/src/trustRegion.jl +++ b/src/trustRegion.jl @@ -235,7 +235,7 @@ function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg::TrustRegion, args...; alias_u0 = false, maxiters = 1000, - abstol = 1e-6, + abstol = 1e-8, internalnorm = DEFAULT_NORM, kwargs...) where {uType, iip} if alias_u0