From e8b377ecc3f04aea99fba4176c74235893246d63 Mon Sep 17 00:00:00 2001 From: Avik Pal Date: Tue, 29 Oct 2024 17:34:00 -0400 Subject: [PATCH] test: move tests for QuasiNewton solvers --- .../CI_NonlinearSolveQuasiNewton.yml | 71 ++++++ common/nlls_problem_workloads.jl | 13 +- common/nonlinear_problem_workloads.jl | 8 +- .../ext/NonlinearSolveBaseLinearSolveExt.jl | 6 +- lib/NonlinearSolveBase/src/tracing.jl | 11 +- lib/NonlinearSolveBase/src/utils.jl | 18 +- lib/NonlinearSolveQuasiNewton/Project.toml | 22 +- lib/NonlinearSolveQuasiNewton/src/broyden.jl | 14 +- .../src/initialization.jl | 12 +- lib/NonlinearSolveQuasiNewton/src/klement.jl | 2 +- lib/NonlinearSolveQuasiNewton/src/solve.jl | 53 ++--- .../test/core_tests.jl | 219 ++++++++++++++++++ .../test/qa_tests.jl | 22 ++ .../test/runtests.jl | 22 ++ .../test/core_tests.jl | 10 +- test/core/rootfind_tests.jl | 213 ----------------- 16 files changed, 440 insertions(+), 276 deletions(-) create mode 100644 .github/workflows/CI_NonlinearSolveQuasiNewton.yml create mode 100644 lib/NonlinearSolveQuasiNewton/test/core_tests.jl create mode 100644 lib/NonlinearSolveQuasiNewton/test/qa_tests.jl diff --git a/.github/workflows/CI_NonlinearSolveQuasiNewton.yml b/.github/workflows/CI_NonlinearSolveQuasiNewton.yml new file mode 100644 index 000000000..5d6024b06 --- /dev/null +++ b/.github/workflows/CI_NonlinearSolveQuasiNewton.yml @@ -0,0 +1,71 @@ +name: CI (NonlinearSolveQuasiNewton) + +on: + pull_request: + branches: + - master + paths: + - "lib/NonlinearSolveQuasiNewton/**" + - ".github/workflows/CI_NonlinearSolveQuasiNewton.yml" + - "lib/NonlinearSolveBase/**" + - "lib/SciMLJacobianOperators/**" + push: + branches: + - master + +concurrency: + # Skip intermediate builds: always. + # Cancel intermediate builds: only if it is a pull request build. + group: ${{ github.workflow }}-${{ github.ref }} + cancel-in-progress: ${{ startsWith(github.ref, 'refs/pull/') }} + +jobs: + test: + runs-on: ${{ matrix.os }} + strategy: + fail-fast: false + matrix: + version: + - "lts" + - "1" + os: + - ubuntu-latest + - macos-latest + - windows-latest + steps: + - uses: actions/checkout@v4 + - uses: julia-actions/setup-julia@v2 + with: + version: ${{ matrix.version }} + - uses: actions/cache@v4 + env: + cache-name: cache-artifacts + with: + path: ~/.julia/artifacts + key: ${{ runner.os }}-test-${{ env.cache-name }}-${{ hashFiles('**/Project.toml') }} + restore-keys: | + ${{ runner.os }}-test-${{ env.cache-name }}- + ${{ runner.os }}-test- + ${{ runner.os }}- + - name: "Install Dependencies and Run Tests" + run: | + import Pkg + Pkg.Registry.update() + # Install packages present in subdirectories + dev_pks = Pkg.PackageSpec[] + for path in ("lib/SciMLJacobianOperators", "lib/NonlinearSolveBase") + push!(dev_pks, Pkg.PackageSpec(; path)) + end + Pkg.develop(dev_pks) + Pkg.instantiate() + Pkg.test(; coverage="user") + shell: julia --color=yes --code-coverage=user --depwarn=yes --project=lib/NonlinearSolveQuasiNewton {0} + - uses: julia-actions/julia-processcoverage@v1 + with: + directories: lib/NonlinearSolveQuasiNewton/src,lib/NonlinearSolveBase/src,lib/NonlinearSolveBase/ext,lib/SciMLJacobianOperators/src + - uses: codecov/codecov-action@v4 + with: + file: lcov.info + token: ${{ secrets.CODECOV_TOKEN }} + verbose: true + fail_ci_if_error: true diff --git a/common/nlls_problem_workloads.jl b/common/nlls_problem_workloads.jl index 4c49e5fa1..0d7c0f7ef 100644 --- a/common/nlls_problem_workloads.jl +++ b/common/nlls_problem_workloads.jl @@ -1,16 +1,19 @@ -using SciMLBase: NonlinearLeastSquaresProblem, NonlinearFunction +using SciMLBase: NonlinearLeastSquaresProblem, NonlinearFunction, NoSpecialize nonlinear_functions = ( - (NonlinearFunction{false}((u, p) -> (u .^ 2 .- p)[1:1]), [0.1, 0.0]), - (NonlinearFunction{false}((u, p) -> vcat(u .* u .- p, u .* u .- p)), [0.1, 0.1]), + (NonlinearFunction{false, NoSpecialize}((u, p) -> (u .^ 2 .- p)[1:1]), [0.1, 0.0]), ( - NonlinearFunction{true}( + NonlinearFunction{false, NoSpecialize}((u, p) -> vcat(u .* u .- p, u .* u .- p)), + [0.1, 0.1] + ), + ( + NonlinearFunction{true, NoSpecialize}( (du, u, p) -> du[1] = u[1] * u[1] - p, resid_prototype = zeros(1) ), [0.1, 0.0] ), ( - NonlinearFunction{true}( + NonlinearFunction{true, NoSpecialize}( (du, u, p) -> du .= vcat(u .* u .- p, u .* u .- p), resid_prototype = zeros(4) ), [0.1, 0.1] diff --git a/common/nonlinear_problem_workloads.jl b/common/nonlinear_problem_workloads.jl index c97592f90..43d0de29a 100644 --- a/common/nonlinear_problem_workloads.jl +++ b/common/nonlinear_problem_workloads.jl @@ -1,9 +1,9 @@ -using SciMLBase: NonlinearProblem, NonlinearFunction +using SciMLBase: NonlinearProblem, NonlinearFunction, NoSpecialize nonlinear_functions = ( - (NonlinearFunction{false}((u, p) -> u .* u .- p), 0.1), - (NonlinearFunction{false}((u, p) -> u .* u .- p), [0.1]), - (NonlinearFunction{true}((du, u, p) -> du .= u .* u .- p), [0.1]) + (NonlinearFunction{false, NoSpecialize}((u, p) -> u .* u .- p), 0.1), + (NonlinearFunction{false, NoSpecialize}((u, p) -> u .* u .- p), [0.1]), + (NonlinearFunction{true, NoSpecialize}((du, u, p) -> du .= u .* u .- p), [0.1]) ) nonlinear_problems = NonlinearProblem[] diff --git a/lib/NonlinearSolveBase/ext/NonlinearSolveBaseLinearSolveExt.jl b/lib/NonlinearSolveBase/ext/NonlinearSolveBaseLinearSolveExt.jl index 28b8b1937..d1829804d 100644 --- a/lib/NonlinearSolveBase/ext/NonlinearSolveBaseLinearSolveExt.jl +++ b/lib/NonlinearSolveBase/ext/NonlinearSolveBaseLinearSolveExt.jl @@ -26,8 +26,10 @@ function (cache::LinearSolveJLCache)(; if cache.precs === nothing Pl, Pr = nothing, nothing else - Pl, Pr = cache.precs(cache.lincache.A, du, linu, p, nothing, - A !== nothing, Plprev, Prprev, cachedata) + Pl, Pr = cache.precs( + cache.lincache.A, du, linu, p, nothing, + A !== nothing, Plprev, Prprev, cachedata + ) end if Pl !== nothing || Pr !== nothing diff --git a/lib/NonlinearSolveBase/src/tracing.jl b/lib/NonlinearSolveBase/src/tracing.jl index e818e05e8..d090c7c40 100644 --- a/lib/NonlinearSolveBase/src/tracing.jl +++ b/lib/NonlinearSolveBase/src/tracing.jl @@ -203,10 +203,13 @@ function update_trace!( if show_now || store_now entry = if trace.trace_level.trace_mode isa Val{:minimal} NonlinearSolveTraceEntry(trace.prob, iter, fu, δu .* α, missing, missing) - elseif trace.trace_level.trace_mode isa Val{:condition_number} - NonlinearSolveTraceEntry(trace.prob, iter, fu, δu .* α, J, missing) else - NonlinearSolveTraceEntry(trace.prob, iter, fu, δu .* α, J, u) + J = convert(AbstractArray, J) + if trace.trace_level.trace_mode isa Val{:condition_number} + NonlinearSolveTraceEntry(trace.prob, iter, fu, δu .* α, J, missing) + else + NonlinearSolveTraceEntry(trace.prob, iter, fu, δu .* α, J, u) + end end show_now && show(stdout, MIME"text/plain"(), entry) store_now && push!(trace.history, entry) @@ -224,7 +227,7 @@ function update_trace!(cache, α = true; uses_jac_inverse = Val(false)) trace, cache.nsteps + 1, get_u(cache), get_fu(cache), nothing, cache.du, α ) else - J = uses_jac_inverse isa Val{true} ? pinv(J) : J + J = uses_jac_inverse isa Val{true} ? Utils.Pinv(cache.J) : cache.J update_trace!(trace, cache.nsteps + 1, get_u(cache), get_fu(cache), J, cache.du, α) end end diff --git a/lib/NonlinearSolveBase/src/utils.jl b/lib/NonlinearSolveBase/src/utils.jl index 5e4374e00..5b0bc2484 100644 --- a/lib/NonlinearSolveBase/src/utils.jl +++ b/lib/NonlinearSolveBase/src/utils.jl @@ -1,6 +1,7 @@ module Utils using ArrayInterface: ArrayInterface +using ConcreteStructs: @concrete using FastClosures: @closure using LinearAlgebra: LinearAlgebra, Diagonal, Symmetric, norm, dot, cond, diagind, pinv using MaybeInplace: @bb @@ -15,6 +16,17 @@ is_extension_loaded(::Val) = false fast_scalar_indexing(xs...) = all(ArrayInterface.fast_scalar_indexing, xs) +@concrete struct Pinv + J +end + +function Base.convert(::Type{AbstractArray}, A::Pinv) + hasmethod(pinv, Tuple{typeof(A.J)}) && return pinv(A.J) + @warn "`pinv` not defined for $(typeof(A.J)). Jacobian will not be inverted when \ + tracing." maxlog=1 + return A.J +end + function nonallocating_isapprox(x::Number, y::Number; atol = false, rtol = atol > 0 ? false : sqrt(eps(promote_type(typeof(x), typeof(y))))) return isapprox(x, y; atol, rtol) @@ -223,7 +235,11 @@ end make_identity!!(::T, α) where {T <: Number} = T(α) function make_identity!!(A::AbstractVector{T}, α) where {T} - @bb @. A = T(α) + if ArrayInterface.can_setindex(A) + @. A = α + else + A = one.(A) .* α + end return A end function make_identity!!(::SMatrix{S1, S2, T, L}, α) where {S1, S2, T, L} diff --git a/lib/NonlinearSolveQuasiNewton/Project.toml b/lib/NonlinearSolveQuasiNewton/Project.toml index ae49ced0c..dee429de8 100644 --- a/lib/NonlinearSolveQuasiNewton/Project.toml +++ b/lib/NonlinearSolveQuasiNewton/Project.toml @@ -8,7 +8,6 @@ ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9" CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2" ConcreteStructs = "2569d6c7-a4a2-43d3-a901-331e8e4be471" DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" -LineSearch = "87fe0de2-c867-4266-b59a-2f0a94fc965b" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" MaybeInplace = "bb5d69b7-63fc-4a16-80bd-7e42200c7bdb" @@ -20,16 +19,22 @@ SciMLOperators = "c0aeaf25-5076-4817-a8d5-81caf7dfa961" StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c" [compat] +ADTypes = "1.9.0" Aqua = "0.8" ArrayInterface = "7.16.0" +BenchmarkTools = "1.5.0" CommonSolve = "0.2.4" ConcreteStructs = "0.2.3" DiffEqBase = "6.155.3" +Enzyme = "0.13.12" ExplicitImports = "1.5" +FiniteDiff = "2.26.0" +ForwardDiff = "0.10.36" Hwloc = "3" InteractiveUtils = "<0.0.1, 1" LineSearch = "0.1.4" -LinearAlgebra = "1.11.0" +LineSearches = "7.3.0" +LinearAlgebra = "1.10" LinearSolve = "2.36.1" MaybeInplace = "0.1.4" NonlinearProblemLibrary = "0.1.2" @@ -41,20 +46,31 @@ Reexport = "1" SciMLBase = "2.54" SciMLOperators = "0.3.11" StableRNGs = "1" +StaticArrays = "1.9.8" StaticArraysCore = "1.4.3" Test = "1.10" +Zygote = "0.6.72" julia = "1.10" [extras] +ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" +BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" +Enzyme = "7da242da-08ed-463a-9acd-ee780be4f1d9" ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7" +FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41" +ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" Hwloc = "0e44f5e4-bd66-52a0-8798-143a42290a1d" InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" +LineSearch = "87fe0de2-c867-4266-b59a-2f0a94fc965b" +LineSearches = "d3d80556-e9d4-5f37-9878-2ab0fcc64255" NonlinearProblemLibrary = "b7050fa9-e91f-4b37-bcee-a89a063da141" Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" ReTestItems = "817f1d60-ba6b-4fd5-9520-3cf149f6a823" StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" +StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" [targets] -test = ["Aqua", "ExplicitImports", "Hwloc", "InteractiveUtils", "NonlinearProblemLibrary", "Pkg", "ReTestItems", "StableRNGs", "Test"] +test = ["ADTypes", "Aqua", "BenchmarkTools", "Enzyme", "ExplicitImports", "FiniteDiff", "ForwardDiff", "Hwloc", "InteractiveUtils", "LineSearch", "LineSearches", "NonlinearProblemLibrary", "Pkg", "ReTestItems", "StableRNGs", "StaticArrays", "Test", "Zygote"] diff --git a/lib/NonlinearSolveQuasiNewton/src/broyden.jl b/lib/NonlinearSolveQuasiNewton/src/broyden.jl index 08570a735..ed487bf94 100644 --- a/lib/NonlinearSolveQuasiNewton/src/broyden.jl +++ b/lib/NonlinearSolveQuasiNewton/src/broyden.jl @@ -148,16 +148,14 @@ function InternalAPI.solve!( @bb @. cache.dfu = fu - cache.dfu J⁻¹_diag = Utils.restructure(cache.dfu, diag(J⁻¹)) if cache.rule isa GoodBroydenUpdateRule - @bb @. J⁻¹_diag = J⁻¹_diag * cache.dfu * du - denom = sum(J⁻¹_diag) - @bb @. J⁻¹_diag = J⁻¹_diag + - (du - J⁻¹_diag * cache.dfu) * du * J⁻¹_diag / - ifelse(iszero(denom), T(1e-5), denom) + @bb @. cache.J⁻¹dfu = J⁻¹_diag * cache.dfu * du + denom = sum(cache.J⁻¹dfu) + @bb @. J⁻¹_diag += (du - cache.J⁻¹dfu) * du * J⁻¹_diag / + ifelse(iszero(denom), T(1e-5), denom) else denom = cache.internalnorm(cache.dfu)^2 - @bb @. J⁻¹_diag = J⁻¹_diag + - (du - J⁻¹_diag * cache.dfu) * cache.dfu / - ifelse(iszero(denom), T(1e-5), denom) + @bb @. J⁻¹_diag += (du - J⁻¹_diag * cache.dfu) * cache.dfu / + ifelse(iszero(denom), T(1e-5), denom) end @bb copyto!(cache.dfu, fu) return Diagonal(vec(J⁻¹_diag)) diff --git a/lib/NonlinearSolveQuasiNewton/src/initialization.jl b/lib/NonlinearSolveQuasiNewton/src/initialization.jl index 6c22e0240..4af1fb921 100644 --- a/lib/NonlinearSolveQuasiNewton/src/initialization.jl +++ b/lib/NonlinearSolveQuasiNewton/src/initialization.jl @@ -166,6 +166,10 @@ function InternalAPI.init( J = BroydenLowRankJacobian(fu, u; alg.threshold, alpha = α) else threshold = min(Utils.unwrap_val(alg.threshold), maxiters) + if threshold > length(u) + @warn "`threshold` is larger than the size of the state, which may cause \ + numerical instability. Consider reducing `threshold`." + end J = BroydenLowRankJacobian(fu, u; threshold, alpha = α) end return InitializedApproximateJacobianCache( @@ -240,9 +244,9 @@ function LinearAlgebra.mul!(y::AbstractVector, J::BroydenLowRankJacobian, x::Abs @. y = -J.alpha * x return y end - _, U, Vᵀ = get_components(J) + cache, U, Vᵀ = get_components(J) @bb cache = Vᵀ × x - mul!(y, U, cache) + LinearAlgebra.mul!(y, U, cache) @bb @. y -= J.alpha * x return y end @@ -258,9 +262,9 @@ function LinearAlgebra.mul!(y::AbstractVector, x::AbstractVector, J::BroydenLowR @. y = -J.alpha * x return y end - _, U, Vᵀ = get_components(J) + cache, U, Vᵀ = get_components(J) @bb cache = transpose(U) × x - mul!(y, transpose(Vᵀ), cache) + LinearAlgebra.mul!(y, transpose(Vᵀ), cache) @bb @. y -= J.alpha * x return y end diff --git a/lib/NonlinearSolveQuasiNewton/src/klement.jl b/lib/NonlinearSolveQuasiNewton/src/klement.jl index 7031b7cec..4bd09f45e 100644 --- a/lib/NonlinearSolveQuasiNewton/src/klement.jl +++ b/lib/NonlinearSolveQuasiNewton/src/klement.jl @@ -112,7 +112,7 @@ function InternalAPI.solve!( T = eltype(u) J = Utils.restructure(u, diag(J)) @bb @. cache.Jdu = (J^2) * (du^2) - @bb @. J += ((fu - cache.fu_cache - cache.Jdu) / + @bb @. J += ((fu - cache.fu_cache - J * du) / ifelse(iszero(cache.Jdu), T(1e-5), cache.Jdu)) * du * (J^2) @bb copyto!(cache.fu_cache, fu) return Diagonal(vec(J)) diff --git a/lib/NonlinearSolveQuasiNewton/src/solve.jl b/lib/NonlinearSolveQuasiNewton/src/solve.jl index 5b987a10d..c52a425ae 100644 --- a/lib/NonlinearSolveQuasiNewton/src/solve.jl +++ b/lib/NonlinearSolveQuasiNewton/src/solve.jl @@ -95,34 +95,31 @@ end kwargs end -# XXX: Implement -# function __reinit_internal!(cache::QuasiNewtonCache{INV, GB, iip}, -# args...; p = cache.p, u0 = cache.u, alias_u0::Bool = false, -# maxiters = 1000, maxtime = nothing, kwargs...) where {INV, GB, iip} -# if iip -# recursivecopy!(cache.u, u0) -# cache.prob.f(cache.fu, cache.u, p) -# else -# cache.u = __maybe_unaliased(u0, alias_u0) -# set_fu!(cache, cache.prob.f(cache.u, p)) -# end -# cache.p = p - -# __reinit_internal!(cache.stats) -# cache.nsteps = 0 -# cache.nresets = 0 -# cache.steps_since_last_reset = 0 -# cache.maxiters = maxiters -# cache.maxtime = maxtime -# cache.total_time = 0.0 -# cache.force_stop = false -# cache.force_reinit = false -# cache.retcode = ReturnCode.Default - -# reset!(cache.trace) -# reinit!(cache.termination_cache, get_fu(cache), get_u(cache); kwargs...) -# reset_timer!(cache.timer) -# end +function InternalAPI.reinit_self!( + cache::QuasiNewtonCache, args...; p = cache.p, u0 = cache.u, + alias_u0::Bool = false, maxiters = 1000, maxtime = nothing, kwargs... +) + Utils.reinit_common!(cache, u0, p, alias_u0) + + InternalAPI.reinit!(cache.stats) + cache.nsteps = 0 + cache.nresets = 0 + cache.steps_since_last_reset = 0 + cache.maxiters = maxiters + cache.maxtime = maxtime + cache.total_time = 0.0 + cache.force_stop = false + cache.force_reinit = false + cache.retcode = ReturnCode.Default + + NonlinearSolveBase.reset!(cache.trace) + SciMLBase.reinit!( + cache.termination_cache, NonlinearSolveBase.get_fu(cache), + NonlinearSolveBase.get_u(cache); kwargs... + ) + NonlinearSolveBase.reset_timer!(cache.timer) + return +end NonlinearSolveBase.@internal_caches(QuasiNewtonCache, :initialization_cache, :descent_cache, :linesearch_cache, :trustregion_cache, diff --git a/lib/NonlinearSolveQuasiNewton/test/core_tests.jl b/lib/NonlinearSolveQuasiNewton/test/core_tests.jl new file mode 100644 index 000000000..a50558867 --- /dev/null +++ b/lib/NonlinearSolveQuasiNewton/test/core_tests.jl @@ -0,0 +1,219 @@ +@testsetup module CoreRootfindTesting + +include("../../../common/common_core_testing.jl") + +end + +@testitem "Broyden" setup=[CoreRootfindTesting] tags=[:core] begin + using ADTypes, LineSearch + using LineSearches: LineSearches + using BenchmarkTools: @ballocated + using StaticArrays: @SVector + using Zygote, Enzyme, ForwardDiff, FiniteDiff + + u0s = ([1.0, 1.0], @SVector[1.0, 1.0], 1.0) + + @testset for ad in (AutoForwardDiff(), AutoZygote(), AutoFiniteDiff(), AutoEnzyme()) + @testset "$(nameof(typeof(linesearch)))" for linesearch in ( + # LineSearchesJL(; method = LineSearches.Static(), autodiff = ad), + # LineSearchesJL(; method = LineSearches.BackTracking(), autodiff = ad), + # LineSearchesJL(; method = LineSearches.MoreThuente(), autodiff = ad), + # LineSearchesJL(; method = LineSearches.StrongWolfe(), autodiff = ad), + LineSearchesJL(; method = LineSearches.HagerZhang(), autodiff = ad), + BackTracking(; autodiff = ad), + LiFukushimaLineSearch() + ) + @testset for init_jacobian in (Val(:identity), Val(:true_jacobian)), + update_rule in (Val(:good_broyden), Val(:bad_broyden), Val(:diagonal)) + + @testset "[OOP] u0: $(typeof(u0))" for u0 in ( + [1.0, 1.0], @SVector[1.0, 1.0], 1.0 + ) + solver = Broyden(; linesearch, init_jacobian, update_rule) + sol = solve_oop(quadratic_f, u0; solver) + @test SciMLBase.successful_retcode(sol) + err = maximum(abs, quadratic_f(sol.u, 2.0)) + @test err < 1e-9 + + cache = init( + NonlinearProblem{false}(quadratic_f, u0, 2.0), solver, abstol = 1e-9 + ) + @test (@ballocated solve!($cache)) < 200 + end + + @testset "[IIP] u0: $(typeof(u0))" for u0 in ([1.0, 1.0],) + ad isa AutoZygote && continue + + solver = Broyden(; linesearch, init_jacobian, update_rule) + sol = solve_iip(quadratic_f!, u0; solver) + @test SciMLBase.successful_retcode(sol) + err = maximum(abs, quadratic_f(sol.u, 2.0)) + @test err < 1e-9 + + cache = init( + NonlinearProblem{true}(quadratic_f!, u0, 2.0), solver, abstol = 1e-9 + ) + @test (@ballocated solve!($cache)) ≤ 64 + end + end + end + end +end + +@testitem "Broyden: Iterator Interface" setup=[CoreRootfindTesting] tags=[:core] begin + p = range(0.01, 2, length = 200) + @test nlprob_iterator_interface(quadratic_f, p, false, Broyden()) ≈ sqrt.(p) + @test nlprob_iterator_interface(quadratic_f!, p, true, Broyden()) ≈ sqrt.(p) +end + +@testitem "Broyden Termination Conditions" setup=[CoreRootfindTesting] tags=[:core] begin + using StaticArrays: @SVector + + @testset "TC: $(nameof(typeof(termination_condition)))" for termination_condition in TERMINATION_CONDITIONS + @testset "u0: $(typeof(u0))" for u0 in ([1.0, 1.0], 1.0, @SVector([1.0, 1.0])) + probN = NonlinearProblem(quadratic_f, u0, 2.0) + sol = solve(probN, Broyden(); termination_condition) + @test all(abs.(quadratic_f(sol.u, 2.0)) .< 1e-9) + end + end +end + +@testitem "Klement" setup=[CoreRootfindTesting] tags=[:core] begin + using ADTypes, LineSearch + using LineSearches: LineSearches + using BenchmarkTools: @ballocated + using StaticArrays: @SVector + using Zygote, Enzyme, ForwardDiff, FiniteDiff + + @testset for ad in (AutoForwardDiff(), AutoZygote(), AutoFiniteDiff(), AutoEnzyme()) + @testset "$(nameof(typeof(linesearch)))" for linesearch in ( + # LineSearchesJL(; method = LineSearches.Static(), autodiff = ad), + # LineSearchesJL(; method = LineSearches.BackTracking(), autodiff = ad), + # LineSearchesJL(; method = LineSearches.MoreThuente(), autodiff = ad), + # LineSearchesJL(; method = LineSearches.StrongWolfe(), autodiff = ad), + LineSearchesJL(; method = LineSearches.HagerZhang(), autodiff = ad), + BackTracking(; autodiff = ad), + LiFukushimaLineSearch() + ) + @testset for init_jacobian in ( + Val(:identity), Val(:true_jacobian), Val(:true_jacobian_diagonal)) + @testset "[OOP] u0: $(typeof(u0))" for u0 in ( + [1.0, 1.0], @SVector[1.0, 1.0], 1.0 + ) + solver = Klement(; linesearch, init_jacobian) + sol = solve_oop(quadratic_f, u0; solver) + # XXX: some tests are failing by a margin + # @test SciMLBase.successful_retcode(sol) + err = maximum(abs, quadratic_f(sol.u, 2.0)) + @test err < 1e-9 + + cache = init( + NonlinearProblem{false}(quadratic_f, u0, 2.0), solver, abstol = 1e-9 + ) + @test (@ballocated solve!($cache)) < 200 + end + + @testset "[IIP] u0: $(typeof(u0))" for u0 in ([1.0, 1.0],) + ad isa AutoZygote && continue + + solver = Klement(; linesearch, init_jacobian) + sol = solve_iip(quadratic_f!, u0; solver) + @test SciMLBase.successful_retcode(sol) + err = maximum(abs, quadratic_f(sol.u, 2.0)) + @test err < 1e-9 + + cache = init( + NonlinearProblem{true}(quadratic_f!, u0, 2.0), solver, abstol = 1e-9 + ) + @test (@ballocated solve!($cache)) ≤ 64 + end + end + end + end +end + +@testitem "Klement: Iterator Interface" setup=[CoreRootfindTesting] tags=[:core] begin + p = range(0.01, 2, length = 200) + @test nlprob_iterator_interface(quadratic_f, p, false, Klement()) ≈ sqrt.(p) + @test nlprob_iterator_interface(quadratic_f!, p, true, Klement()) ≈ sqrt.(p) +end + +@testitem "Klement Termination Conditions" setup=[CoreRootfindTesting] tags=[:core] begin + using StaticArrays: @SVector + + @testset "TC: $(nameof(typeof(termination_condition)))" for termination_condition in TERMINATION_CONDITIONS + @testset "u0: $(typeof(u0))" for u0 in ([1.0, 1.0], 1.0, @SVector([1.0, 1.0])) + probN = NonlinearProblem(quadratic_f, u0, 2.0) + sol = solve(probN, Klement(); termination_condition) + @test all(abs.(quadratic_f(sol.u, 2.0)) .< 1e-9) + end + end +end + +@testitem "LimitedMemoryBroyden" setup=[CoreRootfindTesting] tags=[:core] begin + using ADTypes, LineSearch + using LineSearches: LineSearches + using BenchmarkTools: @ballocated + using StaticArrays: @SVector + using Zygote, Enzyme, ForwardDiff, FiniteDiff + + @testset for ad in (AutoForwardDiff(), AutoZygote(), AutoFiniteDiff(), AutoEnzyme()) + @testset "$(nameof(typeof(linesearch)))" for linesearch in ( + # LineSearchesJL(; method = LineSearches.Static(), autodiff = ad), + # LineSearchesJL(; method = LineSearches.BackTracking(), autodiff = ad), + # LineSearchesJL(; method = LineSearches.MoreThuente(), autodiff = ad), + # LineSearchesJL(; method = LineSearches.StrongWolfe(), autodiff = ad), + LineSearchesJL(; method = LineSearches.HagerZhang(), autodiff = ad), + BackTracking(; autodiff = ad), + LiFukushimaLineSearch() + ) + @testset "[OOP] u0: $(typeof(u0))" for u0 in (ones(32), @SVector(ones(2)), 1.0) + solver = LimitedMemoryBroyden(; linesearch) + sol = solve_oop(quadratic_f, u0; solver) + @test SciMLBase.successful_retcode(sol) + err = maximum(abs, quadratic_f(sol.u, 2.0)) + @test err < 1e-9 + + cache = init( + NonlinearProblem{false}(quadratic_f, u0, 2.0), solver, abstol = 1e-9 + ) + @test (@ballocated solve!($cache)) ≤ 320 + end + + @testset "[IIP] u0: $(typeof(u0))" for u0 in (ones(32),) + ad isa AutoZygote && continue + + solver = LimitedMemoryBroyden(; linesearch) + sol = solve_iip(quadratic_f!, u0; solver) + @test SciMLBase.successful_retcode(sol) + err = maximum(abs, quadratic_f(sol.u, 2.0)) + @test err < 1e-9 + + cache = init( + NonlinearProblem{true}(quadratic_f!, u0, 2.0), solver, abstol = 1e-9 + ) + @test (@ballocated solve!($cache)) ≤ 64 + end + end + end +end + +@testitem "LimitedMemoryBroyden: Iterator Interface" setup=[CoreRootfindTesting] tags=[:core] begin + p = range(0.01, 2, length = 200) + @test nlprob_iterator_interface(quadratic_f, p, false, LimitedMemoryBroyden()) ≈ + sqrt.(p) + @test nlprob_iterator_interface(quadratic_f!, p, true, LimitedMemoryBroyden()) ≈ + sqrt.(p) +end + +@testitem "LimitedMemoryBroyden Termination Conditions" setup=[CoreRootfindTesting] tags=[:core] begin + using StaticArrays: @SVector + + @testset "TC: $(nameof(typeof(termination_condition)))" for termination_condition in TERMINATION_CONDITIONS + @testset "u0: $(typeof(u0))" for u0 in ([1.0, 1.0], 1.0, @SVector([1.0, 1.0])) + probN = NonlinearProblem(quadratic_f, u0, 2.0) + sol = solve(probN, LimitedMemoryBroyden(); termination_condition) + @test all(abs.(quadratic_f(sol.u, 2.0)) .< 1e-9) + end + end +end diff --git a/lib/NonlinearSolveQuasiNewton/test/qa_tests.jl b/lib/NonlinearSolveQuasiNewton/test/qa_tests.jl new file mode 100644 index 000000000..badeab3ef --- /dev/null +++ b/lib/NonlinearSolveQuasiNewton/test/qa_tests.jl @@ -0,0 +1,22 @@ +@testitem "Aqua" tags=[:core] begin + using Aqua, NonlinearSolveQuasiNewton + + Aqua.test_all( + NonlinearSolveQuasiNewton; + piracies = false, ambiguities = false, stale_deps = false, deps_compat = false + ) + Aqua.test_stale_deps(NonlinearSolveQuasiNewton; ignore = [:SciMLJacobianOperators]) + Aqua.test_deps_compat(NonlinearSolveQuasiNewton; ignore = [:SciMLJacobianOperators]) + Aqua.test_piracies(NonlinearSolveQuasiNewton) + Aqua.test_ambiguities(NonlinearSolveQuasiNewton; recursive = false) +end + +@testitem "Explicit Imports" tags=[:core] begin + using ExplicitImports, NonlinearSolveQuasiNewton + + @test check_no_implicit_imports( + NonlinearSolveQuasiNewton; skip = (Base, Core, SciMLBase) + ) === nothing + @test check_no_stale_explicit_imports(NonlinearSolveQuasiNewton) === nothing + @test check_all_qualified_accesses_via_owners(NonlinearSolveQuasiNewton) === nothing +end diff --git a/lib/NonlinearSolveQuasiNewton/test/runtests.jl b/lib/NonlinearSolveQuasiNewton/test/runtests.jl index 8b1378917..807441bbc 100644 --- a/lib/NonlinearSolveQuasiNewton/test/runtests.jl +++ b/lib/NonlinearSolveQuasiNewton/test/runtests.jl @@ -1 +1,23 @@ +using ReTestItems, NonlinearSolveQuasiNewton, Hwloc, InteractiveUtils, Pkg +@info sprint(InteractiveUtils.versioninfo) + +const GROUP = lowercase(get(ENV, "GROUP", "All")) + +const RETESTITEMS_NWORKERS = parse( + Int, get(ENV, "RETESTITEMS_NWORKERS", string(min(Hwloc.num_physical_cores(), 4))) +) +const RETESTITEMS_NWORKER_THREADS = parse(Int, + get( + ENV, "RETESTITEMS_NWORKER_THREADS", + string(max(Hwloc.num_virtual_cores() ÷ RETESTITEMS_NWORKERS, 1)) + ) +) + +@info "Running tests for group: $(GROUP) with $(RETESTITEMS_NWORKERS) workers" + +ReTestItems.runtests( + NonlinearSolveQuasiNewton; tags = (GROUP == "all" ? nothing : [Symbol(GROUP)]), + nworkers = RETESTITEMS_NWORKERS, nworker_threads = RETESTITEMS_NWORKER_THREADS, + testitem_timeout = 3600 +) diff --git a/lib/NonlinearSolveSpectralMethods/test/core_tests.jl b/lib/NonlinearSolveSpectralMethods/test/core_tests.jl index 7f9a411e4..217086dc4 100644 --- a/lib/NonlinearSolveSpectralMethods/test/core_tests.jl +++ b/lib/NonlinearSolveSpectralMethods/test/core_tests.jl @@ -13,7 +13,8 @@ end @testset "[OOP] u0: $(typeof(u0))" for u0 in u0s sol = solve_oop(quadratic_f, u0; solver = DFSane()) @test SciMLBase.successful_retcode(sol) - @test all(abs.(sol.u .* sol.u .- 2) .< 1e-9) + err = maximum(abs, quadratic_f(sol.u, 2.0)) + @test err < 1e-9 cache = init(NonlinearProblem{false}(quadratic_f, u0, 2.0), DFSane(), abstol = 1e-9) @test (@ballocated solve!($cache)) < 200 @@ -22,7 +23,8 @@ end @testset "[IIP] u0: $(typeof(u0))" for u0 in ([1.0, 1.0],) sol = solve_iip(quadratic_f!, u0; solver = DFSane()) @test SciMLBase.successful_retcode(sol) - @test all(abs.(sol.u .* sol.u .- 2) .< 1e-9) + err = maximum(abs, quadratic_f(sol.u, 2.0)) + @test err < 1e-9 cache = init(NonlinearProblem{true}(quadratic_f!, u0, 2.0), DFSane(), abstol = 1e-9) @test (@ballocated solve!($cache)) ≤ 64 @@ -73,8 +75,10 @@ end end @testitem "DFSane Termination Conditions" setup=[CoreRootfindTesting] tags=[:core] begin + using StaticArrays: @SVector + @testset "TC: $(nameof(typeof(termination_condition)))" for termination_condition in TERMINATION_CONDITIONS - @testset "u0: $(typeof(u0))" for u0 in ([1.0, 1.0], 1.0) + @testset "u0: $(typeof(u0))" for u0 in ([1.0, 1.0], 1.0, @SVector([1.0, 1.0])) probN = NonlinearProblem(quadratic_f, u0, 2.0) sol = solve(probN, DFSane(); termination_condition) @test all(abs.(quadratic_f(sol.u, 2.0)) .< 1e-10) diff --git a/test/core/rootfind_tests.jl b/test/core/rootfind_tests.jl index 275502d77..ad09ffd2f 100644 --- a/test/core/rootfind_tests.jl +++ b/test/core/rootfind_tests.jl @@ -345,52 +345,6 @@ end end end -# --- DFSane tests --- - -@testitem "DFSane" setup=[CoreRootfindTesting] tags=[:core] begin - # Test that `DFSane` passes a test that `NewtonRaphson` fails on. - @testset "Newton Raphson Fails" begin - 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; solver = DFSane()) - @test SciMLBase.successful_retcode(sol) - @test all(abs.(newton_fails(sol.u, p)) .< 1e-9) - end - - # Test kwargs in `DFSane` - @testset "Keyword Arguments" begin - σ_min = [1e-10, 1e-5, 1e-4] - σ_max = [1e10, 1e5, 1e4] - σ_1 = [1.0, 0.5, 2.0] - M = [10, 1, 100] - γ = [1e-4, 1e-3, 1e-5] - τ_min = [0.1, 0.2, 0.3] - τ_max = [0.5, 0.8, 0.9] - nexp = [2, 1, 2] - η_strategy = [(f_1, k, x, F) -> f_1 / k^2, (f_1, k, x, F) -> f_1 / k^3, - (f_1, k, x, F) -> f_1 / k^4] - - list_of_options = zip(σ_min, σ_max, σ_1, M, γ, τ_min, τ_max, nexp, η_strategy) - for options in list_of_options - local probN, sol, alg - alg = DFSane(σ_min = options[1], σ_max = options[2], σ_1 = options[3], - M = options[4], γ = options[5], τ_min = options[6], - τ_max = options[7], n_exp = options[8], η_strategy = options[9]) - - probN = NonlinearProblem{false}(quadratic_f, [1.0, 1.0], 2.0) - sol = solve(probN, alg, abstol = 1e-11) - @test all(abs.(quadratic_f(sol.u, 2.0)) .< 1e-6) - end - end - - @testset "Termination condition: $(_nameof(termination_condition)) u0: $(_nameof(u0))" for termination_condition in TERMINATION_CONDITIONS, - u0 in (1.0, [1.0, 1.0]) - - probN = NonlinearProblem(quadratic_f, u0, 2.0) - @test all(solve(probN, DFSane(); termination_condition).u .≈ sqrt(2.0)) - end -end - # --- PseudoTransient tests --- @testitem "PseudoTransient" setup=[CoreRootfindTesting] tags=[:core] begin @@ -460,173 +414,6 @@ end end end -# --- Broyden tests --- - -@testitem "Broyden" setup=[CoreRootfindTesting] tags=[:core] begin - @testset "LineSearch: $(_nameof(linesearch)) LineSearch AD: $(_nameof(ad)) Init Jacobian: $(init_jacobian) Update Rule: $(update_rule)" for ad in ( - AutoForwardDiff(), AutoZygote(), AutoFiniteDiff() - ), - linesearch in ( - LineSearchesStatic(; autodiff = ad), LineSearchesStrongWolfe(; autodiff = ad), - LineSearchesBackTracking(; autodiff = ad), BackTracking(; autodiff = ad), - LineSearchesHagerZhang(; autodiff = ad), - LineSearchesMoreThuente(; autodiff = ad) - ), - init_jacobian in (Val(:identity), Val(:true_jacobian)), - update_rule in (Val(:good_broyden), Val(:bad_broyden), Val(:diagonal)) - - u0s = ([1.0, 1.0], @SVector[1.0, 1.0], 1.0) - - @testset "[OOP] u0: $(typeof(u0))" for u0 in u0s - solver = Broyden(; linesearch, init_jacobian, update_rule) - sol = benchmark_nlsolve_oop(quadratic_f, u0; solver) - @test SciMLBase.successful_retcode(sol) - @test all(abs.(sol.u .* sol.u .- 2) .< 1e-9) - - cache = init(NonlinearProblem{false}(quadratic_f, u0, 2.0), - Broyden(; linesearch, update_rule, init_jacobian), abstol = 1e-9) - @test (@ballocated solve!($cache)) < 200 - end - - @testset "[IIP] u0: $(typeof(u0))" for u0 in ([1.0, 1.0],) - ad isa AutoZygote && continue - solver = Broyden(; linesearch, init_jacobian, update_rule) - sol = benchmark_nlsolve_iip(quadratic_f!, u0; solver) - @test SciMLBase.successful_retcode(sol) - @test all(abs.(sol.u .* sol.u .- 2) .< 1e-9) - - cache = init(NonlinearProblem{true}(quadratic_f!, u0, 2.0), - Broyden(; linesearch, update_rule, init_jacobian), abstol = 1e-9) - @test (@ballocated solve!($cache)) ≤ 64 - end - end - - # Iterator interface - p = range(0.01, 2, length = 200) - @test nlprob_iterator_interface(quadratic_f, p, Val(false), Broyden()) ≈ sqrt.(p) - @test nlprob_iterator_interface(quadratic_f!, p, Val(true), Broyden()) ≈ sqrt.(p) - - @testset "Termination condition: $(_nameof(termination_condition)) u0: $(_nameof(u0))" for termination_condition in TERMINATION_CONDITIONS, - u0 in (1.0, [1.0, 1.0]) - - probN = NonlinearProblem(quadratic_f, u0, 2.0) - @test all(solve(probN, Broyden(); termination_condition).u .≈ sqrt(2.0)) - end -end - -# --- Klement tests --- - -@testitem "Klement" setup=[CoreRootfindTesting] tags=[:core] begin - @testset "LineSearch: $(_nameof(linesearch)) LineSearch AD: $(_nameof(ad)) Init Jacobian: $(init_jacobian)" for ad in ( - AutoForwardDiff(), AutoZygote(), AutoFiniteDiff() - ), - linesearch in ( - LineSearchesStatic(; autodiff = ad), LineSearchesStrongWolfe(; autodiff = ad), - LineSearchesBackTracking(; autodiff = ad), BackTracking(; autodiff = ad), - LineSearchesHagerZhang(; autodiff = ad), - LineSearchesMoreThuente(; autodiff = ad) - ), - init_jacobian in (Val(:identity), Val(:true_jacobian), Val(:true_jacobian_diagonal)) - - u0s = ([1.0, 1.0], @SVector[1.0, 1.0], 1.0) - - @testset "[OOP] u0: $(typeof(u0))" for u0 in u0s - solver = Klement(; linesearch, init_jacobian) - sol = benchmark_nlsolve_oop(quadratic_f, u0; solver) - # Some are failing by a margin - # @test SciMLBase.successful_retcode(sol) - @test all(abs.(sol.u .* sol.u .- 2) .< 3e-9) - - cache = init(NonlinearProblem{false}(quadratic_f, u0, 2.0), - Klement(; linesearch), abstol = 1e-9) - @test (@ballocated solve!($cache)) < 200 - end - - @testset "[IIP] u0: $(typeof(u0))" for u0 in ([1.0, 1.0],) - ad isa AutoZygote && continue - solver = Klement(; linesearch, init_jacobian) - sol = benchmark_nlsolve_iip(quadratic_f!, u0; solver) - @test SciMLBase.successful_retcode(sol) - @test all(abs.(sol.u .* sol.u .- 2) .< 1e-9) - - cache = init(NonlinearProblem{true}(quadratic_f!, u0, 2.0), - Klement(; linesearch), abstol = 1e-9) - @test (@ballocated solve!($cache)) ≤ 64 - end - end - - # Iterator interface - p = range(0.01, 2, length = 200) - @test nlprob_iterator_interface(quadratic_f, p, Val(false), Klement()) ≈ sqrt.(p) - @test nlprob_iterator_interface(quadratic_f!, p, Val(true), Klement()) ≈ sqrt.(p) - - @testset "Termination condition: $(_nameof(termination_condition)) u0: $(_nameof(u0))" for termination_condition in TERMINATION_CONDITIONS, - u0 in (1.0, [1.0, 1.0]) - - probN = NonlinearProblem(quadratic_f, u0, 2.0) - @test all(solve(probN, Klement(); termination_condition).u .≈ sqrt(2.0)) - end -end - -# --- LimitedMemoryBroyden tests --- - -@testitem "LimitedMemoryBroyden" setup=[CoreRootfindTesting] tags=[:core] begin - @testset "LineSearch: $(_nameof(linesearch)) LineSearch AD: $(_nameof(ad))" for ad in ( - AutoForwardDiff(), AutoZygote(), AutoFiniteDiff() - ), - linesearch in ( - LineSearchesStatic(; autodiff = ad), LineSearchesStrongWolfe(; autodiff = ad), - LineSearchesBackTracking(; autodiff = ad), BackTracking(; autodiff = ad), - LineSearchesHagerZhang(; autodiff = ad), - LineSearchesMoreThuente(; autodiff = ad), LiFukushimaLineSearch() - ) - - u0s = ([1.0, 1.0], @SVector[1.0, 1.0], 1.0) - - @testset "[OOP] u0: $(typeof(u0))" for u0 in u0s - broken = linesearch isa BackTracking && ad isa AutoFiniteDiff && u0 isa Vector - - solver = LimitedMemoryBroyden(; linesearch) - sol = benchmark_nlsolve_oop(quadratic_f, u0; solver) - @test SciMLBase.successful_retcode(sol) broken=broken - @test all(abs.(sol.u .* sol.u .- 2) .< 1e-9) broken=broken - - cache = init(NonlinearProblem{false}(quadratic_f, u0, 2.0), - LimitedMemoryBroyden(; linesearch), abstol = 1e-9) - @test (@ballocated solve!($cache)) < 200 - end - - @testset "[IIP] u0: $(typeof(u0))" for u0 in ([1.0, 1.0],) - broken = linesearch isa BackTracking && ad isa AutoFiniteDiff && u0 isa Vector - ad isa AutoZygote && continue - - solver = LimitedMemoryBroyden(; linesearch) - sol = benchmark_nlsolve_iip(quadratic_f!, u0; solver) - @test SciMLBase.successful_retcode(sol) broken=broken - @test all(abs.(sol.u .* sol.u .- 2) .< 1e-9) broken=broken - - cache = init(NonlinearProblem{true}(quadratic_f!, u0, 2.0), - LimitedMemoryBroyden(; linesearch), abstol = 1e-9) - @test (@ballocated solve!($cache)) ≤ 64 - end - end - - # Iterator interface - p = range(0.01, 2, length = 200) - @test nlprob_iterator_interface( - quadratic_f, p, Val(false), LimitedMemoryBroyden())≈sqrt.(p) atol=1e-2 - @test nlprob_iterator_interface( - quadratic_f!, p, Val(true), LimitedMemoryBroyden())≈sqrt.(p) atol=1e-2 - - @testset "Termination condition: $(_nameof(termination_condition)) u0: $(_nameof(u0))" for termination_condition in TERMINATION_CONDITIONS, - u0 in (1.0, [1.0, 1.0]) - - probN = NonlinearProblem(quadratic_f, u0, 2.0) - @test all(solve(probN, LimitedMemoryBroyden(); termination_condition).u .≈ - sqrt(2.0)) - end -end - # Miscellaneous Tests @testitem "Custom JVP" setup=[CoreRootfindTesting] tags=[:core] begin function F(u::Vector{Float64}, p::Vector{Float64})