diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index b8e0c0c0b..3bcd2581a 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -10,6 +10,11 @@ on: - master paths-ignore: - 'docs/**' +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: ubuntu-latest diff --git a/.github/workflows/Downstream.yml b/.github/workflows/Downstream.yml index f8110ab41..b977c8409 100644 --- a/.github/workflows/Downstream.yml +++ b/.github/workflows/Downstream.yml @@ -4,7 +4,11 @@ on: branches: [master] tags: [v*] pull_request: - +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: name: ${{ matrix.package.repo }}/${{ matrix.package.group }}/${{ matrix.julia-version }} diff --git a/Project.toml b/Project.toml index 078c37f3b..af3db4c85 100644 --- a/Project.toml +++ b/Project.toml @@ -40,7 +40,7 @@ ADTypes = "0.2" ArrayInterface = "6.0.24, 7" BandedMatrices = "1" ConcreteStructs = "0.2" -DiffEqBase = "6.130" +DiffEqBase = "6.136" EnumX = "1" Enzyme = "0.11" FastBroadcast = "0.1.9, 0.2" @@ -56,7 +56,7 @@ RecursiveArrayTools = "2" Reexport = "0.2, 1" SciMLBase = "2.4" SimpleNonlinearSolve = "0.1.23" -SparseDiffTools = "2.6" +SparseDiffTools = "2.9" StaticArraysCore = "1.4" UnPack = "1.0" Zygote = "0.6" diff --git a/docs/Project.toml b/docs/Project.toml index d240bb245..9507a2966 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -2,12 +2,14 @@ AlgebraicMultigrid = "2169fc97-5a83-5252-b627-83903c6c433c" ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9" BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" +DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" IncompleteLU = "40713840-3770-5561-ab4c-a76e7d0d7895" LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec" NonlinearSolveMINPACK = "c100e077-885d-495a-a2ea-599e143bf69d" +SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" SciMLNLSolve = "e9a6253c-8580-4d32-9898-8661bb511710" SimpleNonlinearSolve = "727e6d20-b764-4bd8-a329-72de5adea6c7" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" @@ -19,12 +21,14 @@ Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" AlgebraicMultigrid = "0.5, 0.6" ArrayInterface = "6, 7" BenchmarkTools = "1" +DiffEqBase = "6.136" Documenter = "1" IncompleteLU = "0.2" LinearSolve = "2" ModelingToolkit = "8" NonlinearSolve = "1, 2" NonlinearSolveMINPACK = "0.1" +SciMLBase = "2.4" SciMLNLSolve = "0.1" SimpleNonlinearSolve = "0.1.5" StaticArrays = "1" diff --git a/docs/make.jl b/docs/make.jl index 8f895b0b4..9a063bf47 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,5 +1,5 @@ using Documenter, NonlinearSolve, SimpleNonlinearSolve, Sundials, SciMLNLSolve, - NonlinearSolveMINPACK, SteadyStateDiffEq + NonlinearSolveMINPACK, SteadyStateDiffEq, SciMLBase, DiffEqBase cp("./docs/Manifest.toml", "./docs/src/assets/Manifest.toml", force = true) cp("./docs/Project.toml", "./docs/src/assets/Project.toml", force = true) @@ -8,9 +8,8 @@ include("pages.jl") makedocs(sitename = "NonlinearSolve.jl", authors = "Chris Rackauckas", - modules = [NonlinearSolve, NonlinearSolve.SciMLBase, NonlinearSolve.DiffEqBase, - SimpleNonlinearSolve, Sundials, SciMLNLSolve, NonlinearSolveMINPACK, - SteadyStateDiffEq], + modules = [NonlinearSolve, SciMLBase, DiffEqBase, SimpleNonlinearSolve, Sundials, + SciMLNLSolve, NonlinearSolveMINPACK, SteadyStateDiffEq], clean = true, doctest = false, linkcheck = true, linkcheck_ignore = ["https://twitter.com/ChrisRackauckas/status/1544743542094020615"], warnonly = [:missing_docs, :cross_references], diff --git a/docs/pages.jl b/docs/pages.jl index f1b239bc9..a68b0f8e2 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -6,7 +6,6 @@ pages = ["index.md", "Handling Large Ill-Conditioned and Sparse Systems" => "tutorials/large_systems.md", "Symbolic System Definition and Acceleration via ModelingToolkit" => "tutorials/modelingtoolkit.md", "tutorials/small_compile.md", - "tutorials/termination_conditions.md", "tutorials/iterator_interface.md"], "Basics" => Any["basics/NonlinearProblem.md", "basics/NonlinearFunctions.md", diff --git a/docs/src/basics/TerminationCondition.md b/docs/src/basics/TerminationCondition.md index 97bec6d32..a4ff5272b 100644 --- a/docs/src/basics/TerminationCondition.md +++ b/docs/src/basics/TerminationCondition.md @@ -1,10 +1,75 @@ # [Termination Conditions](@id termination_condition) Provides a API to specify termination conditions for [`NonlinearProblem`](@ref) and -[`SteadyStateProblem`](@ref). For details on the various termination modes, i.e., -NLSolveTerminationMode, see the documentation for [`NLSolveTerminationCondition`](@ref). +[`SteadyStateProblem`](@ref). For details on the various termination modes: -## Termination Condition API +## Termination Conditions + +The termination condition is constructed as: + +```julia +cache = init(du, u, AbsNormTerminationMode(); abstol = 1e-9, reltol = 1e-9) +``` + +If `abstol` and `reltol` are not supplied, then we choose a default based on the element +types of `du` and `u`. + +We can query the `cache` using `DiffEqBase.get_termination_mode`, `DiffEqBase.get_abstol` +and `DiffEqBase.get_reltol`. + +To test for termination simply call the `cache`: + +```julia +terminated = cache(du, u, uprev) +``` + +!!! note + + The default for NonlinearSolve.jl is `AbsSafeBestTerminationMode`! + +### Absolute Tolerance + +```@docs +AbsTerminationMode +AbsNormTerminationMode +AbsSafeTerminationMode +AbsSafeBestTerminationMode +``` + +### Relative Tolerance + +```@docs +RelTerminationMode +RelNormTerminationMode +RelSafeTerminationMode +RelSafeBestTerminationMode +``` + +### Both Absolute and Relative Tolerance + +```@docs +NormTerminationMode +SteadyStateDiffEqTerminationMode +SimpleNonlinearSolveTerminationMode +``` + +### Return Codes + +```@docs +DiffEqBase.NonlinearSafeTerminationReturnCode +DiffEqBase.NonlinearSafeTerminationReturnCode.Success +DiffEqBase.NonlinearSafeTerminationReturnCode.Default +DiffEqBase.NonlinearSafeTerminationReturnCode.Failure +DiffEqBase.NonlinearSafeTerminationReturnCode.PatienceTermination +DiffEqBase.NonlinearSafeTerminationReturnCode.ProtectiveTermination +``` + +## [Deprecated] Termination Condition API + +!!! warning + + This is deprecated. Currently only parts of `SimpleNonlinearSolve` uses this API. That + will also be phased out soon! ```@docs NLSolveTerminationCondition diff --git a/docs/src/tutorials/termination_conditions.md b/docs/src/tutorials/termination_conditions.md deleted file mode 100644 index 6a7e8a3cd..000000000 --- a/docs/src/tutorials/termination_conditions.md +++ /dev/null @@ -1,3 +0,0 @@ -# [More Detailed Termination Conditions](@id termination_conditions_tutorial) - -This is a stub article to be completed soon. diff --git a/ext/NonlinearSolveFastLevenbergMarquardtExt.jl b/ext/NonlinearSolveFastLevenbergMarquardtExt.jl index 4a096747d..08d719f96 100644 --- a/ext/NonlinearSolveFastLevenbergMarquardtExt.jl +++ b/ext/NonlinearSolveFastLevenbergMarquardtExt.jl @@ -32,8 +32,8 @@ end (f::InplaceFunction{false})(fx, x, p) = (fx .= f.f(x, p)) function SciMLBase.__init(prob::NonlinearLeastSquaresProblem, - alg::FastLevenbergMarquardtJL, args...; abstol = 1e-8, reltol = 1e-8, - verbose = false, maxiters = 1000, kwargs...) + alg::FastLevenbergMarquardtJL, args...; abstol = 1e-8, reltol = 1e-8, + verbose = false, maxiters = 1000, kwargs...) iip = SciMLBase.isinplace(prob) @assert prob.f.jac!==nothing "FastLevenbergMarquardt requires a Jacobian!" diff --git a/ext/NonlinearSolveLeastSquaresOptimExt.jl b/ext/NonlinearSolveLeastSquaresOptimExt.jl index e1f3cdfb7..5cf6ab67d 100644 --- a/ext/NonlinearSolveLeastSquaresOptimExt.jl +++ b/ext/NonlinearSolveLeastSquaresOptimExt.jl @@ -33,7 +33,7 @@ end (f::FunctionWrapper{false})(du, u) = (du .= f.f(u, f.p)) function SciMLBase.__init(prob::NonlinearLeastSquaresProblem, alg::LeastSquaresOptimJL, - args...; abstol = 1e-8, reltol = 1e-8, verbose = false, maxiters = 1000, kwargs...) + args...; abstol = 1e-8, reltol = 1e-8, verbose = false, maxiters = 1000, kwargs...) iip = SciMLBase.isinplace(prob) f! = FunctionWrapper{iip}(prob.f, prob.p) diff --git a/src/NonlinearSolve.jl b/src/NonlinearSolve.jl index 0f8626314..2b6ab25e4 100644 --- a/src/NonlinearSolve.jl +++ b/src/NonlinearSolve.jl @@ -14,7 +14,7 @@ PrecompileTools.@recompile_invalidations begin import ADTypes: AbstractFiniteDifferencesMode import ArrayInterface: undefmatrix, - matrix_colors, parameterless_type, ismutable, issingular,fast_scalar_indexing + matrix_colors, parameterless_type, ismutable, issingular, fast_scalar_indexing import ConcreteStructs: @concrete import EnumX: @enumx import ForwardDiff @@ -30,6 +30,9 @@ PrecompileTools.@recompile_invalidations begin end @reexport using ADTypes, LineSearches, SciMLBase, SimpleNonlinearSolve +import DiffEqBase: AbstractNonlinearTerminationMode, + AbstractSafeNonlinearTerminationMode, AbstractSafeBestNonlinearTerminationMode, + NonlinearSafeTerminationReturnCode, get_termination_mode const AbstractSparseADType = Union{ADTypes.AbstractSparseFiniteDifferences, ADTypes.AbstractSparseForwardMode, ADTypes.AbstractSparseReverseMode} @@ -44,7 +47,7 @@ abstract type AbstractNonlinearSolveCache{iip} end isinplace(::AbstractNonlinearSolveCache{iip}) where {iip} = iip function SciMLBase.__solve(prob::Union{NonlinearProblem, NonlinearLeastSquaresProblem}, - alg::AbstractNonlinearSolveAlgorithm, args...; kwargs...) + alg::AbstractNonlinearSolveAlgorithm, args...; kwargs...) cache = init(prob, alg, args...; kwargs...) return solve!(cache) end @@ -53,6 +56,9 @@ function not_terminated(cache::AbstractNonlinearSolveCache) return !cache.force_stop && cache.stats.nsteps < cache.maxiters end get_fu(cache::AbstractNonlinearSolveCache) = cache.fu1 +set_fu!(cache::AbstractNonlinearSolveCache, fu) = (cache.fu1 = fu) +get_u(cache::AbstractNonlinearSolveCache) = cache.u +SciMLBase.set_u!(cache::AbstractNonlinearSolveCache, u) = (cache.u = u) function SciMLBase.solve!(cache::AbstractNonlinearSolveCache) while not_terminated(cache) @@ -69,7 +75,7 @@ function SciMLBase.solve!(cache::AbstractNonlinearSolveCache) end end - return SciMLBase.build_solution(cache.prob, cache.alg, cache.u, get_fu(cache); + return SciMLBase.build_solution(cache.prob, cache.alg, get_u(cache), get_fu(cache); cache.retcode, cache.stats) end @@ -96,7 +102,7 @@ PrecompileTools.@compile_workload begin NonlinearProblem{true}((du, u, p) -> du .= u .* u .- p, T[0.1], T[2])) precompile_algs = (NewtonRaphson(), TrustRegion(), LevenbergMarquardt(), - PseudoTransient(), GeneralBroyden(), GeneralKlement(), nothing) + PseudoTransient(), GeneralBroyden(), GeneralKlement(), DFSane(), nothing) for prob in probs, alg in precompile_algs solve(prob, alg, abstol = T(1e-2)) @@ -113,4 +119,10 @@ export RobustMultiNewton, FastShortcutNonlinearPolyalg export LineSearch, LiFukushimaLineSearch +# Export the termination conditions from DiffEqBase +export SteadyStateDiffEqTerminationMode, SimpleNonlinearSolveTerminationMode, + NormTerminationMode, RelTerminationMode, RelNormTerminationMode, AbsTerminationMode, + AbsNormTerminationMode, RelSafeTerminationMode, AbsSafeTerminationMode, + RelSafeBestTerminationMode, AbsSafeBestTerminationMode + end # module diff --git a/src/ad.jl b/src/ad.jl index 1af781e21..73ea45fab 100644 --- a/src/ad.jl +++ b/src/ad.jl @@ -26,16 +26,16 @@ function scalar_nlsolve_ad(prob, alg, args...; kwargs...) end function SciMLBase.solve(prob::NonlinearProblem{<:Union{Number, SVector, <:AbstractArray}, - false, <:Dual{T, V, P}}, alg::AbstractNewtonAlgorithm, args...; - kwargs...) where {T, V, P} + false, <:Dual{T, V, P}}, alg::AbstractNonlinearSolveAlgorithm, args...; + kwargs...) where {T, V, P} sol, partials = scalar_nlsolve_ad(prob, alg, args...; kwargs...) dual_soln = scalar_nlsolve_dual_soln(sol.u, partials, prob.p) return SciMLBase.build_solution(prob, alg, dual_soln, sol.resid; sol.retcode) end function SciMLBase.solve(prob::NonlinearProblem{<:Union{Number, SVector, <:AbstractArray}, - false, <:AbstractArray{<:Dual{T, V, P}}}, alg::AbstractNewtonAlgorithm, args...; - kwargs...) where {T, V, P} + false, <:AbstractArray{<:Dual{T, V, P}}}, alg::AbstractNonlinearSolveAlgorithm, + args...; kwargs...) where {T, V, P} sol, partials = scalar_nlsolve_ad(prob, alg, args...; kwargs...) dual_soln = scalar_nlsolve_dual_soln(sol.u, partials, prob.p) return SciMLBase.build_solution(prob, alg, dual_soln, sol.resid; sol.retcode) @@ -53,11 +53,11 @@ function scalar_nlsolve_∂f_∂u(f, u, p) end function scalar_nlsolve_dual_soln(u::Number, partials, - ::Union{<:AbstractArray{<:Dual{T, V, P}}, Dual{T, V, P}}) where {T, V, P} + ::Union{<:AbstractArray{<:Dual{T, V, P}}, Dual{T, V, P}}) where {T, V, P} return Dual{T, V, P}(u, partials) end function scalar_nlsolve_dual_soln(u::AbstractArray, partials, - ::Union{<:AbstractArray{<:Dual{T, V, P}}, Dual{T, V, P}}) where {T, V, P} + ::Union{<:AbstractArray{<:Dual{T, V, P}}, Dual{T, V, P}}) where {T, V, P} return map(((uᵢ, pᵢ),) -> Dual{T, V, P}(uᵢ, pᵢ), zip(u, partials)) end diff --git a/src/broyden.jl b/src/broyden.jl index 5fcbd3d51..4d3be3224 100644 --- a/src/broyden.jl +++ b/src/broyden.jl @@ -22,7 +22,7 @@ An implementation of `Broyden` with reseting and line search. end function GeneralBroyden(; max_resets = 3, linesearch = LineSearch(), - reset_tolerance = nothing) + reset_tolerance = nothing) linesearch = linesearch isa LineSearch ? linesearch : LineSearch(; method = linesearch) return GeneralBroyden(max_resets, reset_tolerance, linesearch) end @@ -52,17 +52,17 @@ end reset_check prob stats::NLStats - lscache - termination_condition - tc_storage + ls_cache + tc_cache end get_fu(cache::GeneralBroydenCache) = cache.fu +set_fu!(cache::GeneralBroydenCache, fu) = (cache.fu = fu) function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg::GeneralBroyden, args...; - alias_u0 = false, maxiters = 1000, abstol = nothing, reltol = nothing, - termination_condition = nothing, internalnorm = DEFAULT_NORM, - kwargs...) where {uType, iip} + alias_u0 = false, maxiters = 1000, abstol = nothing, reltol = nothing, + termination_condition = nothing, internalnorm::F = DEFAULT_NORM, + kwargs...) where {uType, iip, F} @unpack f, u0, p = prob u = alias_u0 ? u0 : deepcopy(u0) fu = evaluate_f(prob, u) @@ -71,37 +71,26 @@ function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg::GeneralBroyde alg.reset_tolerance reset_check = x -> abs(x) ≤ reset_tolerance - abstol, reltol, termination_condition = _init_termination_elements(abstol, - reltol, - termination_condition, - eltype(u)) - - mode = DiffEqBase.get_termination_mode(termination_condition) + abstol, reltol, tc_cache = init_termination_cache(abstol, reltol, fu, u, + termination_condition) - storage = mode ∈ DiffEqBase.SAFE_TERMINATION_MODES ? NLSolveSafeTerminationResult() : - nothing return GeneralBroydenCache{iip}(f, alg, u, zero(u), _mutable_zero(u), fu, zero(fu), zero(fu), p, J⁻¹, zero(_reshape(fu, 1, :)), _mutable_zero(u), false, 0, alg.max_resets, maxiters, internalnorm, ReturnCode.Default, abstol, reltol, - reset_tolerance, - reset_check, prob, NLStats(1, 0, 0, 0, 0), - init_linesearch_cache(alg.linesearch, f, u, p, fu, Val(iip)), termination_condition, - storage) + reset_tolerance, reset_check, prob, NLStats(1, 0, 0, 0, 0), + init_linesearch_cache(alg.linesearch, f, u, p, fu, Val(iip)), tc_cache) end function perform_step!(cache::GeneralBroydenCache{true}) - @unpack f, p, du, fu, fu2, dfu, u, u_prev, J⁻¹, J⁻¹df, J⁻¹₂, tc_storage = cache - - termination_condition = cache.termination_condition(tc_storage) + @unpack f, p, du, fu, fu2, dfu, u, u_prev, J⁻¹, J⁻¹df, J⁻¹₂ = cache T = eltype(u) - mul!(_vec(du), J⁻¹, -_vec(fu)) - α = perform_linesearch!(cache.lscache, u, du) - _axpy!(α, du, u) + mul!(_vec(du), J⁻¹, _vec(fu)) + α = perform_linesearch!(cache.ls_cache, u, du) + _axpy!(-α, du, u) f(fu2, u, p) - termination_condition(fu2, u, u_prev, cache.abstol, cache.reltol) && - (cache.force_stop = true) + check_and_update!(cache, fu2, u, u_prev) cache.stats.nf += 1 cache.force_stop && return nothing @@ -111,7 +100,7 @@ function perform_step!(cache::GeneralBroydenCache{true}) if all(cache.reset_check, du) || all(cache.reset_check, dfu) if cache.resets ≥ cache.max_resets - cache.retcode = ReturnCode.Unstable + cache.retcode = ReturnCode.ConvergenceFailure cache.force_stop = true return nothing end @@ -119,6 +108,7 @@ function perform_step!(cache::GeneralBroydenCache{true}) J⁻¹[diagind(J⁻¹)] .= T(1) cache.resets += 1 else + du .*= -1 mul!(_vec(J⁻¹df), J⁻¹, _vec(dfu)) mul!(J⁻¹₂, _vec(du)', J⁻¹) denom = dot(du, J⁻¹df) @@ -132,19 +122,16 @@ function perform_step!(cache::GeneralBroydenCache{true}) end function perform_step!(cache::GeneralBroydenCache{false}) - @unpack f, p, tc_storage = cache - - termination_condition = cache.termination_condition(tc_storage) + @unpack f, p = cache T = eltype(cache.u) - cache.du = _restructure(cache.du, cache.J⁻¹ * -_vec(cache.fu)) - α = perform_linesearch!(cache.lscache, cache.u, cache.du) - cache.u = cache.u .+ α * cache.du + cache.du = _restructure(cache.du, cache.J⁻¹ * _vec(cache.fu)) + α = perform_linesearch!(cache.ls_cache, cache.u, cache.du) + cache.u = cache.u .- α * cache.du cache.fu2 = f(cache.u, p) - termination_condition(cache.fu2, cache.u, cache.u_prev, cache.abstol, cache.reltol) && - (cache.force_stop = true) + check_and_update!(cache, cache.fu2, cache.u, cache.u_prev) cache.stats.nf += 1 cache.force_stop && return nothing @@ -153,13 +140,14 @@ function perform_step!(cache::GeneralBroydenCache{false}) cache.dfu = cache.fu2 .- cache.fu if all(cache.reset_check, cache.du) || all(cache.reset_check, cache.dfu) if cache.resets ≥ cache.max_resets - cache.retcode = ReturnCode.Unstable + cache.retcode = ReturnCode.ConvergenceFailure cache.force_stop = true return nothing end cache.J⁻¹ = __init_identity_jacobian(cache.u, cache.fu) cache.resets += 1 else + cache.du = -cache.du cache.J⁻¹df = _restructure(cache.J⁻¹df, cache.J⁻¹ * _vec(cache.dfu)) cache.J⁻¹₂ = _vec(cache.du)' * cache.J⁻¹ denom = dot(cache.du, cache.J⁻¹df) @@ -173,9 +161,8 @@ function perform_step!(cache::GeneralBroydenCache{false}) end function SciMLBase.reinit!(cache::GeneralBroydenCache{iip}, u0 = cache.u; p = cache.p, - abstol = cache.abstol, reltol = cache.reltol, - termination_condition = cache.termination_condition, - maxiters = cache.maxiters) where {iip} + abstol = cache.abstol, reltol = cache.reltol, maxiters = cache.maxiters, + termination_condition = get_termination_mode(cache.tc_cache)) where {iip} cache.p = p if iip recursivecopy!(cache.u, u0) @@ -186,12 +173,12 @@ function SciMLBase.reinit!(cache::GeneralBroydenCache{iip}, u0 = cache.u; p = ca cache.fu = cache.f(cache.u, p) end - termination_condition = _get_reinit_termination_condition(cache, abstol, reltol, + abstol, reltol, tc_cache = init_termination_cache(abstol, reltol, cache.fu, cache.u, termination_condition) cache.abstol = abstol cache.reltol = reltol - cache.termination_condition = termination_condition + cache.tc_cache = tc_cache cache.maxiters = maxiters cache.stats.nf = 1 cache.stats.nsteps = 1 diff --git a/src/default.jl b/src/default.jl index 3b4f8ef23..178d72f1c 100644 --- a/src/default.jl +++ b/src/default.jl @@ -42,7 +42,7 @@ end const Robusters = RobustMultiNewton function RobustMultiNewton(; concrete_jac = nothing, linsolve = nothing, - precs = DEFAULT_PRECS, adkwargs...) + precs = DEFAULT_PRECS, adkwargs...) return RobustMultiNewton{_unwrap_val(concrete_jac)}(adkwargs, linsolve, precs) end @@ -53,7 +53,7 @@ end end function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg::RobustMultiNewton, - args...; kwargs...) where {uType, iip} + args...; kwargs...) where {uType, iip} @unpack adkwargs, linsolve, precs = alg algs = (TrustRegion(; linsolve, precs, adkwargs...), @@ -105,7 +105,7 @@ for more performance and then tries more robust techniques if the faster ones fa end function FastShortcutNonlinearPolyalg(; concrete_jac = nothing, linsolve = nothing, - precs = DEFAULT_PRECS, adkwargs...) + precs = DEFAULT_PRECS, adkwargs...) return FastShortcutNonlinearPolyalg{_unwrap_val(concrete_jac)}(adkwargs, linsolve, precs) end @@ -118,18 +118,17 @@ end end function FastShortcutNonlinearPolyalgCache(; concrete_jac = nothing, linsolve = nothing, - precs = DEFAULT_PRECS, adkwargs...) + precs = DEFAULT_PRECS, adkwargs...) return FastShortcutNonlinearPolyalgCache{_unwrap_val(concrete_jac)}(adkwargs, linsolve, precs) end function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, - alg::FastShortcutNonlinearPolyalg, args...; kwargs...) where {uType, iip} + alg::FastShortcutNonlinearPolyalg, args...; kwargs...) where {uType, iip} @unpack adkwargs, linsolve, precs = alg - algs = ( - # Klement(), - # Broyden(), + algs = (GeneralKlement(; linsolve, precs), + GeneralBroyden(), NewtonRaphson(; linsolve, precs, adkwargs...), NewtonRaphson(; linsolve, precs, linesearch = BackTracking(), adkwargs...), TrustRegion(; linsolve, precs, adkwargs...), @@ -142,8 +141,8 @@ end # This version doesn't allocate all the caches! @generated function SciMLBase.__solve(prob::NonlinearProblem{uType, iip}, - alg::Union{FastShortcutNonlinearPolyalg, RobustMultiNewton}, args...; - kwargs...) where {uType, iip} + alg::Union{FastShortcutNonlinearPolyalg, RobustMultiNewton}, args...; + kwargs...) where {uType, iip} calls = [:(@unpack adkwargs, linsolve, precs = alg)] algs = if parameterless_type(alg) == RobustMultiNewton @@ -159,7 +158,7 @@ end ] else [ - :(GeneralKlement()), + :(GeneralKlement(; linsolve, precs)), :(GeneralBroyden()), :(NewtonRaphson(; linsolve, precs, adkwargs...)), :(NewtonRaphson(; linsolve, precs, linesearch = BackTracking(), adkwargs...)), @@ -191,7 +190,7 @@ end push!(calls, quote resids = tuple($(Tuple(resids)...)) - minfu, idx = findmin(DEFAULT_NORM, resids) + minfu, idx = __findmin(DEFAULT_NORM, resids) end) for i in 1:length(algs) @@ -211,7 +210,7 @@ end ## General shared polyalg functions @generated function SciMLBase.solve!(cache::Union{RobustMultiNewtonCache{iip, N}, - FastShortcutNonlinearPolyalgCache{iip, N}}) where {iip, N} + FastShortcutNonlinearPolyalgCache{iip, N}}) where {iip, N} calls = [ quote 1 ≤ cache.current ≤ length(cache.caches) || @@ -249,7 +248,7 @@ end retcode = ReturnCode.MaxIters fus = tuple($(Tuple(resids)...)) - minfu, idx = findmin(cache.caches[1].internalnorm, fus) + minfu, idx = __findmin(cache.caches[1].internalnorm, fus) stats = cache.caches[idx].stats u = cache.caches[idx].u @@ -261,7 +260,7 @@ end end function SciMLBase.reinit!(cache::Union{RobustMultiNewtonCache, - FastShortcutNonlinearPolyalgCache}, args...; kwargs...) + FastShortcutNonlinearPolyalgCache}, args...; kwargs...) for c in cache.caches SciMLBase.reinit!(c, args...; kwargs...) end @@ -270,11 +269,11 @@ end ## Defaults function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg::Nothing, args...; - kwargs...) where {uType, iip} + kwargs...) where {uType, iip} SciMLBase.__init(prob, FastShortcutNonlinearPolyalg(), args...; kwargs...) end function SciMLBase.__solve(prob::NonlinearProblem{uType, iip}, alg::Nothing, args...; - kwargs...) where {uType, iip} + kwargs...) where {uType, iip} SciMLBase.__solve(prob, FastShortcutNonlinearPolyalg(), args...; kwargs...) end diff --git a/src/dfsane.jl b/src/dfsane.jl index aca13c344..e158b61c6 100644 --- a/src/dfsane.jl +++ b/src/dfsane.jl @@ -10,78 +10,76 @@ see the paper: [W LaCruz, JM Martinez, and M Raydan (2006), Spectral residual ma gradient information for solving large-scale nonlinear systems of equations, Mathematics of Computation, 75, 1429-1448.](https://www.researchgate.net/publication/220576479_Spectral_Residual_Method_without_Gradient_Information_for_Solving_Large-Scale_Nonlinear_Systems_of_Equations) -See also the implementation in [SimpleNonlinearSolve.jl](https://github.com/SciML/SimpleNonlinearSolve.jl/blob/main/src/dfsane.jl) - ### Keyword Arguments -- `σ_min`: the minimum value of the spectral coefficient `σₙ` which is related to the step - size in the algorithm. Defaults to `1e-10`. -- `σ_max`: the maximum value of the spectral coefficient `σₙ` which is related to the step - size in the algorithm. Defaults to `1e10`. -- `σ_1`: the initial value of the spectral coefficient `σₙ` which is related to the step - size in the algorithm.. Defaults to `1.0`. -- `M`: The monotonicity of the algorithm is determined by a this positive integer. - A value of 1 for `M` would result in strict monotonicity in the decrease of the L2-norm - of the function `f`. However, higher values allow for more flexibility in this reduction. - Despite this, the algorithm still ensures global convergence through the use of a - non-monotone line-search algorithm that adheres to the Grippo-Lampariello-Lucidi - condition. Values in the range of 5 to 20 are usually sufficient, but some cases may call - for a higher value of `M`. The default setting is 10. -- `γ`: a parameter that influences if a proposed step will be accepted. Higher value of `γ` - will make the algorithm more restrictive in accepting steps. Defaults to `1e-4`. -- `τ_min`: if a step is rejected the new step size will get multiplied by factor, and this - parameter is the minimum value of that factor. Defaults to `0.1`. -- `τ_max`: if a step is rejected the new step size will get multiplied by factor, and this - parameter is the maximum value of that factor. Defaults to `0.5`. -- `n_exp`: the exponent of the loss, i.e. ``f_n=||F(x_n)||^{n_exp}``. The paper uses - `n_exp ∈ {1,2}`. Defaults to `2`. -- `η_strategy`: function to determine the parameter `η`, which enables growth - of ``||f_n||^2``. Called as ``η = η_strategy(fn_1, n, x_n, f_n)`` with `fn_1` initialized as - ``fn_1=||f(x_1)||^{n_exp}``, `n` is the iteration number, `x_n` is the current `x`-value and - `f_n` the current residual. Should satisfy ``η > 0`` and ``∑ₖ ηₖ < ∞``. Defaults to - ``fn_1 / n^2``. -- `max_inner_iterations`: the maximum number of iterations allowed for the inner loop of the - algorithm. Defaults to `1000`. + - `σ_min`: the minimum value of the spectral coefficient `σₙ` which is related to the step + size in the algorithm. Defaults to `1e-10`. + - `σ_max`: the maximum value of the spectral coefficient `σₙ` which is related to the step + size in the algorithm. Defaults to `1e10`. + - `σ_1`: the initial value of the spectral coefficient `σₙ` which is related to the step + size in the algorithm.. Defaults to `1.0`. + - `M`: The monotonicity of the algorithm is determined by a this positive integer. + A value of 1 for `M` would result in strict monotonicity in the decrease of the L2-norm + of the function `f`. However, higher values allow for more flexibility in this reduction. + Despite this, the algorithm still ensures global convergence through the use of a + non-monotone line-search algorithm that adheres to the Grippo-Lampariello-Lucidi + condition. Values in the range of 5 to 20 are usually sufficient, but some cases may call + for a higher value of `M`. The default setting is 10. + - `γ`: a parameter that influences if a proposed step will be accepted. Higher value of `γ` + will make the algorithm more restrictive in accepting steps. Defaults to `1e-4`. + - `τ_min`: if a step is rejected the new step size will get multiplied by factor, and this + parameter is the minimum value of that factor. Defaults to `0.1`. + - `τ_max`: if a step is rejected the new step size will get multiplied by factor, and this + parameter is the maximum value of that factor. Defaults to `0.5`. + - `n_exp`: the exponent of the loss, i.e. ``f_n=||F(x_n)||^{n_exp}``. The paper uses + `n_exp ∈ {1,2}`. Defaults to `2`. + - `η_strategy`: function to determine the parameter `η`, which enables growth + of ``||f_n||^2``. Called as ``η = η_strategy(fn_1, n, x_n, f_n)`` with `fn_1` initialized as + ``fn_1=||f(x_1)||^{n_exp}``, `n` is the iteration number, `x_n` is the current `x`-value and + `f_n` the current residual. Should satisfy ``η > 0`` and ``∑ₖ ηₖ < ∞``. Defaults to + ``fn_1 / n^2``. + - `max_inner_iterations`: the maximum number of iterations allowed for the inner loop of the + algorithm. Defaults to `1000`. """ -struct DFSane{T, F} <: AbstractNonlinearSolveAlgorithm - σ_min::T - σ_max::T - σ_1::T +@concrete struct DFSane <: AbstractNonlinearSolveAlgorithm + σ_min + σ_max + σ_1 M::Int - γ::T - τ_min::T - τ_max::T + γ + τ_min + τ_max n_exp::Int - η_strategy::F + η_strategy max_inner_iterations::Int end function DFSane(; σ_min = 1e-10, σ_max = 1e+10, σ_1 = 1.0, M = 10, γ = 1e-4, τ_min = 0.1, - τ_max = 0.5, n_exp = 2, η_strategy = (fn_1, n, x_n, f_n) -> fn_1 / n^2, - max_inner_iterations = 1000) - return DFSane{typeof(σ_min), typeof(η_strategy)}(σ_min, σ_max, σ_1, M, γ, τ_min, τ_max, - n_exp, η_strategy, max_inner_iterations) + τ_max = 0.5, n_exp = 2, η_strategy::F = (fn_1, n, x_n, f_n) -> fn_1 / n^2, + max_inner_iterations = 1000) where {F} + return DFSane(σ_min, σ_max, σ_1, M, γ, τ_min, τ_max, n_exp, η_strategy, + max_inner_iterations) end -@concrete mutable struct DFSaneCache{iip} +@concrete mutable struct DFSaneCache{iip} <: AbstractNonlinearSolveCache{iip} alg - uₙ - uₙ₋₁ - fuₙ - fuₙ₋₁ - 𝒹 - ℋ - f₍ₙₒᵣₘ₎ₙ₋₁ - f₍ₙₒᵣₘ₎₀ + u + uprev + fu + fuprev + du + history + f_norm + f_norm_0 M - σₙ - σₘᵢₙ - σₘₐₓ - α₁ + σ_n + σ_min + σ_max + α_1 γ - τₘᵢₙ - τₘₐₓ - nₑₓₚ::Int + τ_min + τ_max + n_exp::Int p force_stop::Bool maxiters::Int @@ -91,251 +89,211 @@ end reltol prob stats::NLStats - termination_condition - tc_storage + tc_cache end -function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg::DFSane, args...; - alias_u0 = false, maxiters = 1000, abstol = nothing, reltol = nothing, - termination_condition = nothing, internalnorm = DEFAULT_NORM, - kwargs...) where {uType, iip} - uₙ = alias_u0 ? prob.u0 : deepcopy(prob.u0) - - p = prob.p - T = eltype(uₙ) - σₘᵢₙ, σₘₐₓ, γ, τₘᵢₙ, τₘₐₓ = T(alg.σ_min), T(alg.σ_max), T(alg.γ), T(alg.τ_min), - T(alg.τ_max) - α₁ = one(T) - γ = T(alg.γ) - f₍ₙₒᵣₘ₎ₙ₋₁ = α₁ - σₙ = T(alg.σ_1) - M = alg.M - nₑₓₚ = alg.n_exp - 𝒹, uₙ₋₁, fuₙ, fuₙ₋₁ = copy(uₙ), copy(uₙ), copy(uₙ), copy(uₙ) - - if iip - prob.f(fuₙ₋₁, uₙ₋₁, p) - else - fuₙ₋₁ = prob.f(uₙ₋₁, p) - end +get_fu(cache::DFSaneCache) = cache.fu +set_fu!(cache::DFSaneCache, fu) = (cache.fu = fu) - f₍ₙₒᵣₘ₎ₙ₋₁ = norm(fuₙ₋₁)^nₑₓₚ - f₍ₙₒᵣₘ₎₀ = f₍ₙₒᵣₘ₎ₙ₋₁ +function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg::DFSane, args...; + alias_u0 = false, maxiters = 1000, abstol = nothing, reltol = nothing, + termination_condition = nothing, internalnorm::F = DEFAULT_NORM, + kwargs...) where {uType, iip, F} + u = alias_u0 ? prob.u0 : deepcopy(prob.u0) + T = eltype(u) - ℋ = fill(f₍ₙₒᵣₘ₎ₙ₋₁, M) + du, uprev = copy(u), copy(u) + fu = evaluate_f(prob, u) + fuprev = copy(fu) - abstol, reltol, termination_condition = _init_termination_elements(abstol, reltol, - termination_condition, T) + f_norm = internalnorm(fu)^alg.n_exp + f_norm_0 = f_norm - mode = DiffEqBase.get_termination_mode(termination_condition) + history = fill(f_norm, alg.M) - storage = mode ∈ DiffEqBase.SAFE_TERMINATION_MODES ? NLSolveSafeTerminationResult() : - nothing + abstol, reltol, tc_cache = init_termination_cache(abstol, reltol, fu, uprev, + termination_condition) - return DFSaneCache{iip}(alg, uₙ, uₙ₋₁, fuₙ, fuₙ₋₁, 𝒹, ℋ, f₍ₙₒᵣₘ₎ₙ₋₁, f₍ₙₒᵣₘ₎₀, - M, σₙ, σₘᵢₙ, σₘₐₓ, α₁, γ, τₘᵢₙ, τₘₐₓ, nₑₓₚ, p, false, maxiters, - internalnorm, ReturnCode.Default, abstol, reltol, prob, NLStats(1, 0, 0, 0, 0), - termination_condition, storage) + return DFSaneCache{iip}(alg, u, uprev, fu, fuprev, du, history, f_norm, f_norm_0, alg.M, + T(alg.σ_1), T(alg.σ_min), T(alg.σ_max), one(T), T(alg.γ), T(alg.τ_min), + T(alg.τ_max), alg.n_exp, prob.p, false, maxiters, internalnorm, ReturnCode.Default, + abstol, reltol, prob, NLStats(1, 0, 0, 0, 0), tc_cache) end function perform_step!(cache::DFSaneCache{true}) - @unpack alg, f₍ₙₒᵣₘ₎ₙ₋₁, f₍ₙₒᵣₘ₎₀, σₙ, σₘᵢₙ, σₘₐₓ, α₁, γ, τₘᵢₙ, τₘₐₓ, nₑₓₚ, M, tc_storage = cache - - termination_condition = cache.termination_condition(tc_storage) - f = (dx, x) -> cache.prob.f(dx, x, cache.p) - - T = eltype(cache.uₙ) - n = cache.stats.nsteps + @unpack alg, f_norm, σ_n, σ_min, σ_max, α_1, γ, τ_min, τ_max, n_exp, M, prob = cache + T = eltype(cache.u) + f_norm_old = f_norm # Spectral parameter range check - σₙ = sign(σₙ) * clamp(abs(σₙ), σₘᵢₙ, σₘₐₓ) + σ_n = sign(σ_n) * clamp(abs(σ_n), σ_min, σ_max) # Line search direction - @. cache.𝒹 = -σₙ * cache.fuₙ₋₁ + @. cache.du = -σ_n * cache.fuprev - η = alg.η_strategy(f₍ₙₒᵣₘ₎₀, n, cache.uₙ₋₁, cache.fuₙ₋₁) + η = alg.η_strategy(cache.f_norm_0, cache.stats.nsteps, cache.u, cache.fu) - f̄ = maximum(cache.ℋ) - α₊ = α₁ - α₋ = α₁ - @. cache.uₙ = cache.uₙ₋₁ + α₊ * cache.𝒹 + f_bar = maximum(cache.history) + α₊ = α_1 + α₋ = α_1 + _axpy!(α₊, cache.du, cache.u) - f(cache.fuₙ, cache.uₙ) - f₍ₙₒᵣₘ₎ₙ = norm(cache.fuₙ)^nₑₓₚ - for jjj in 1:(cache.alg.max_inner_iterations) - 𝒸 = f̄ + η - γ * α₊^2 * f₍ₙₒᵣₘ₎ₙ₋₁ + prob.f(cache.fu, cache.u, cache.p) + f_norm = cache.internalnorm(cache.fu)^n_exp - f₍ₙₒᵣₘ₎ₙ ≤ 𝒸 && break + # TODO: Failure mode with inner line search failed? + for _ in 1:(cache.alg.max_inner_iterations) + c = f_bar + η - γ * α₊^2 * f_norm_old - α₊ = clamp(α₊^2 * f₍ₙₒᵣₘ₎ₙ₋₁ / (f₍ₙₒᵣₘ₎ₙ + (T(2) * α₊ - T(1)) * f₍ₙₒᵣₘ₎ₙ₋₁), - τₘᵢₙ * α₊, τₘₐₓ * α₊) - @. cache.uₙ = cache.uₙ₋₁ - α₋ * cache.𝒹 + f_norm ≤ c && break - f(cache.fuₙ, cache.uₙ) - f₍ₙₒᵣₘ₎ₙ = norm(cache.fuₙ)^nₑₓₚ + α₊ = α₊ * clamp(α₊ * f_norm_old / (f_norm + (T(2) * α₊ - T(1)) * f_norm_old), + τ_min, τ_max) + @. cache.u = cache.uprev - α₋ * cache.du - f₍ₙₒᵣₘ₎ₙ .≤ 𝒸 && break + prob.f(cache.fu, cache.u, cache.p) + f_norm = cache.internalnorm(cache.fu)^n_exp - α₋ = clamp(α₋^2 * f₍ₙₒᵣₘ₎ₙ₋₁ / (f₍ₙₒᵣₘ₎ₙ + (T(2) * α₋ - T(1)) * f₍ₙₒᵣₘ₎ₙ₋₁), - τₘᵢₙ * α₋, τₘₐₓ * α₋) + f_norm ≤ c && break - @. cache.uₙ = cache.uₙ₋₁ + α₊ * cache.𝒹 - f(cache.fuₙ, cache.uₙ) - f₍ₙₒᵣₘ₎ₙ = norm(cache.fuₙ)^nₑₓₚ - end + α₋ = α₋ * clamp(α₋ * f_norm_old / (f_norm + (T(2) * α₋ - T(1)) * f_norm_old), + τ_min, τ_max) + @. cache.u = cache.uprev + α₊ * cache.du - if termination_condition(cache.fuₙ, cache.uₙ, cache.uₙ₋₁, cache.abstol, cache.reltol) - cache.force_stop = true + prob.f(cache.fu, cache.u, cache.p) + f_norm = cache.internalnorm(cache.fu)^n_exp end + check_and_update!(cache, cache.fu, cache.u, cache.uprev) + # Update spectral parameter - @. cache.uₙ₋₁ = cache.uₙ - cache.uₙ₋₁ - @. cache.fuₙ₋₁ = cache.fuₙ - cache.fuₙ₋₁ + @. cache.uprev = cache.u - cache.uprev + @. cache.fuprev = cache.fu - cache.fuprev - α₊ = sum(abs2, cache.uₙ₋₁) - @. cache.uₙ₋₁ = cache.uₙ₋₁ * cache.fuₙ₋₁ - α₋ = sum(cache.uₙ₋₁) - cache.σₙ = α₊ / α₋ + α₊ = sum(abs2, cache.uprev) + @. cache.uprev *= cache.fuprev + α₋ = sum(cache.uprev) + cache.σ_n = α₊ / α₋ # Spectral parameter bounds check - if abs(cache.σₙ) > σₘₐₓ || abs(cache.σₙ) < σₘᵢₙ - test_norm = sqrt(sum(abs2, cache.fuₙ₋₁)) - cache.σₙ = clamp(T(1) / test_norm, T(1), T(1e5)) + if !(σ_min ≤ abs(cache.σ_n) ≤ σ_max) + test_norm = sqrt(sum(abs2, cache.fuprev)) + cache.σ_n = clamp(inv(test_norm), T(1), T(1e5)) end # Take step - @. cache.uₙ₋₁ = cache.uₙ - @. cache.fuₙ₋₁ = cache.fuₙ - cache.f₍ₙₒᵣₘ₎ₙ₋₁ = f₍ₙₒᵣₘ₎ₙ + @. cache.uprev = cache.u + @. cache.fuprev = cache.fu + cache.f_norm = f_norm # Update history - cache.ℋ[n % M + 1] = f₍ₙₒᵣₘ₎ₙ + cache.history[cache.stats.nsteps % M + 1] = f_norm cache.stats.nf += 1 return nothing end function perform_step!(cache::DFSaneCache{false}) - @unpack alg, f₍ₙₒᵣₘ₎ₙ₋₁, f₍ₙₒᵣₘ₎₀, σₙ, σₘᵢₙ, σₘₐₓ, α₁, γ, τₘᵢₙ, τₘₐₓ, nₑₓₚ, M, tc_storage = cache - - termination_condition = cache.termination_condition(tc_storage) - f = x -> cache.prob.f(x, cache.p) - - T = eltype(cache.uₙ) - n = cache.stats.nsteps + @unpack alg, f_norm, σ_n, σ_min, σ_max, α_1, γ, τ_min, τ_max, n_exp, M, prob = cache + T = eltype(cache.u) + f_norm_old = f_norm # Spectral parameter range check - σₙ = sign(σₙ) * clamp(abs(σₙ), σₘᵢₙ, σₘₐₓ) + σ_n = sign(σ_n) * clamp(abs(σ_n), σ_min, σ_max) # Line search direction - cache.𝒹 = -σₙ * cache.fuₙ₋₁ + cache.du = @. -σ_n * cache.fuprev + + η = alg.η_strategy(cache.f_norm_0, cache.stats.nsteps, cache.u, cache.fu) - η = alg.η_strategy(f₍ₙₒᵣₘ₎₀, n, cache.uₙ₋₁, cache.fuₙ₋₁) + f_bar = maximum(cache.history) + α₊ = α_1 + α₋ = α_1 + cache.u = @. cache.uprev + α₊ * cache.du - f̄ = maximum(cache.ℋ) - α₊ = α₁ - α₋ = α₁ - cache.uₙ = cache.uₙ₋₁ + α₊ * cache.𝒹 + cache.fu = prob.f(cache.u, cache.p) + f_norm = cache.internalnorm(cache.fu)^n_exp - cache.fuₙ = f(cache.uₙ) - f₍ₙₒᵣₘ₎ₙ = norm(cache.fuₙ)^nₑₓₚ + # TODO: Failure mode with inner line search failed? for _ in 1:(cache.alg.max_inner_iterations) - 𝒸 = f̄ + η - γ * α₊^2 * f₍ₙₒᵣₘ₎ₙ₋₁ + c = f_bar + η - γ * α₊^2 * f_norm_old - f₍ₙₒᵣₘ₎ₙ ≤ 𝒸 && break + f_norm ≤ c && break - α₊ = clamp(α₊^2 * f₍ₙₒᵣₘ₎ₙ₋₁ / (f₍ₙₒᵣₘ₎ₙ + (T(2) * α₊ - T(1)) * f₍ₙₒᵣₘ₎ₙ₋₁), - τₘᵢₙ * α₊, τₘₐₓ * α₊) - cache.uₙ = @. cache.uₙ₋₁ - α₋ * cache.𝒹 + α₊ = α₊ * clamp(α₊ * f_norm_old / (f_norm + (T(2) * α₊ - T(1)) * f_norm_old), + τ_min, τ_max) + cache.u = @. cache.uprev - α₋ * cache.du - cache.fuₙ = f(cache.uₙ) - f₍ₙₒᵣₘ₎ₙ = norm(cache.fuₙ)^nₑₓₚ + cache.fu = prob.f(cache.u, cache.p) + f_norm = cache.internalnorm(cache.fu)^n_exp - f₍ₙₒᵣₘ₎ₙ .≤ 𝒸 && break + f_norm ≤ c && break - α₋ = clamp(α₋^2 * f₍ₙₒᵣₘ₎ₙ₋₁ / (f₍ₙₒᵣₘ₎ₙ + (T(2) * α₋ - T(1)) * f₍ₙₒᵣₘ₎ₙ₋₁), - τₘᵢₙ * α₋, τₘₐₓ * α₋) + α₋ = α₋ * clamp(α₋ * f_norm_old / (f_norm + (T(2) * α₋ - T(1)) * f_norm_old), + τ_min, τ_max) + cache.u = @. cache.uprev + α₊ * cache.du - cache.uₙ = @. cache.uₙ₋₁ + α₊ * cache.𝒹 - cache.fuₙ = f(cache.uₙ) - f₍ₙₒᵣₘ₎ₙ = norm(cache.fuₙ)^nₑₓₚ + cache.fu = prob.f(cache.u, cache.p) + f_norm = cache.internalnorm(cache.fu)^n_exp end - if termination_condition(cache.fuₙ, cache.uₙ, cache.uₙ₋₁, cache.abstol, cache.reltol) - cache.force_stop = true - end + check_and_update!(cache, cache.fu, cache.u, cache.uprev) # Update spectral parameter - cache.uₙ₋₁ = @. cache.uₙ - cache.uₙ₋₁ - cache.fuₙ₋₁ = @. cache.fuₙ - cache.fuₙ₋₁ + cache.uprev = @. cache.u - cache.uprev + cache.fuprev = @. cache.fu - cache.fuprev - α₊ = sum(abs2, cache.uₙ₋₁) - cache.uₙ₋₁ = @. cache.uₙ₋₁ * cache.fuₙ₋₁ - α₋ = sum(cache.uₙ₋₁) - cache.σₙ = α₊ / α₋ + α₊ = sum(abs2, cache.uprev) + cache.uprev = @. cache.uprev * cache.fuprev + α₋ = sum(cache.uprev) + cache.σ_n = α₊ / α₋ # Spectral parameter bounds check - if abs(cache.σₙ) > σₘₐₓ || abs(cache.σₙ) < σₘᵢₙ - test_norm = sqrt(sum(abs2, cache.fuₙ₋₁)) - cache.σₙ = clamp(T(1) / test_norm, T(1), T(1e5)) + if !(σ_min ≤ abs(cache.σ_n) ≤ σ_max) + test_norm = sqrt(sum(abs2, cache.fuprev)) + cache.σ_n = clamp(inv(test_norm), T(1), T(1e5)) end # Take step - cache.uₙ₋₁ = cache.uₙ - cache.fuₙ₋₁ = cache.fuₙ - cache.f₍ₙₒᵣₘ₎ₙ₋₁ = f₍ₙₒᵣₘ₎ₙ + cache.uprev = cache.u + cache.fuprev = cache.fu + cache.f_norm = f_norm # Update history - cache.ℋ[n % M + 1] = f₍ₙₒᵣₘ₎ₙ + cache.history[cache.stats.nsteps % M + 1] = f_norm cache.stats.nf += 1 return nothing end -function SciMLBase.solve!(cache::DFSaneCache) - while !cache.force_stop && cache.stats.nsteps < cache.maxiters - cache.stats.nsteps += 1 - perform_step!(cache) - end - - if cache.stats.nsteps == cache.maxiters - cache.retcode = ReturnCode.MaxIters - else - cache.retcode = ReturnCode.Success - end - - return SciMLBase.build_solution(cache.prob, cache.alg, cache.uₙ, cache.fuₙ; - retcode = cache.retcode, stats = cache.stats) -end - -function SciMLBase.reinit!(cache::DFSaneCache{iip}, u0 = cache.uₙ; p = cache.p, - abstol = cache.abstol, reltol = cache.reltol, - termination_condition = cache.termination_condition, - maxiters = cache.maxiters) where {iip} +function SciMLBase.reinit!(cache::DFSaneCache{iip}, u0 = cache.u; p = cache.p, + abstol = cache.abstol, reltol = cache.reltol, maxiters = cache.maxiters, + termination_condition = get_termination_mode(cache.tc_cache)) where {iip} cache.p = p if iip - recursivecopy!(cache.uₙ, u0) - recursivecopy!(cache.uₙ₋₁, u0) - cache.prob.f(cache.fuₙ, cache.uₙ, p) - cache.prob.f(cache.fuₙ₋₁, cache.uₙ, p) + recursivecopy!(cache.u, u0) + recursivecopy!(cache.uprev, u0) + cache.prob.f(cache.fu, cache.u, p) + cache.prob.f(cache.fuprev, cache.uprev, p) else - cache.uₙ = u0 - cache.uₙ₋₁ = u0 - cache.fuₙ = cache.prob.f(cache.uₙ, p) - cache.fuₙ₋₁ = cache.prob.f(cache.uₙ, p) + cache.u = u0 + cache.uprev = u0 + cache.fu = cache.prob.f(cache.u, p) + cache.fuprev = cache.prob.f(cache.uprev, p) end - cache.f₍ₙₒᵣₘ₎ₙ₋₁ = norm(cache.fuₙ₋₁)^cache.nₑₓₚ - cache.f₍ₙₒᵣₘ₎₀ = cache.f₍ₙₒᵣₘ₎ₙ₋₁ - fill!(cache.ℋ, cache.f₍ₙₒᵣₘ₎ₙ₋₁) + cache.f_norm = cache.internalnorm(cache.fu)^cache.n_exp + cache.f_norm_0 = cache.f_norm + + fill!(cache.history, cache.f_norm) - T = eltype(cache.uₙ) - cache.σₙ = T(cache.alg.σ_1) + T = eltype(cache.u) + cache.σ_n = T(cache.alg.σ_1) - termination_condition = _get_reinit_termination_condition(cache, abstol, reltol, + abstol, reltol, tc_cache = init_termination_cache(abstol, reltol, cache.fu, cache.u, termination_condition) cache.abstol = abstol cache.reltol = reltol - cache.termination_condition = termination_condition + cache.tc_cache = tc_cache cache.maxiters = maxiters cache.stats.nf = 1 cache.stats.nsteps = 1 diff --git a/src/extension_algs.jl b/src/extension_algs.jl index 6038a7623..c4e56b6a5 100644 --- a/src/extension_algs.jl +++ b/src/extension_algs.jl @@ -65,8 +65,8 @@ Wrapper over [FastLevenbergMarquardt.jl](https://github.com/kamesy/FastLevenberg end function FastLevenbergMarquardtJL(linsolve::Symbol = :cholesky; factor = 1e-6, - factoraccept = 13.0, factorreject = 3.0, factorupdate = :marquardt, - minscale = 1e-12, maxscale = 1e16, minfactor = 1e-28, maxfactor = 1e32) + factoraccept = 13.0, factorreject = 3.0, factorupdate = :marquardt, + minscale = 1e-12, maxscale = 1e16, minfactor = 1e-28, maxfactor = 1e32) @assert linsolve in (:qr, :cholesky) @assert factorupdate in (:marquardt, :nielson) diff --git a/src/gaussnewton.jl b/src/gaussnewton.jl index b066f6169..c857f2d23 100644 --- a/src/gaussnewton.jl +++ b/src/gaussnewton.jl @@ -47,11 +47,9 @@ function set_ad(alg::GaussNewton{CJ}, ad) where {CJ} end function GaussNewton(; concrete_jac = nothing, linsolve = nothing, - precs = DEFAULT_PRECS, adkwargs...) + precs = DEFAULT_PRECS, adkwargs...) ad = default_adargs_to_adtype(; adkwargs...) - return GaussNewton{_unwrap_val(concrete_jac)}(ad, - linsolve, - precs) + return GaussNewton{_unwrap_val(concrete_jac)}(ad, linsolve, precs) end @concrete mutable struct GaussNewtonCache{iip} <: AbstractNonlinearSolveCache{iip} @@ -78,27 +76,21 @@ end reltol prob stats::NLStats - tc_storage - termination_condition + tc_cache_1 + tc_cache_2 end function SciMLBase.__init(prob::NonlinearLeastSquaresProblem{uType, iip}, alg_::GaussNewton, - args...; alias_u0 = false, maxiters = 1000, abstol = nothing, reltol = nothing, - termination_condition = nothing, - internalnorm = DEFAULT_NORM, - kwargs...) where {uType, iip} + args...; alias_u0 = false, maxiters = 1000, abstol = nothing, reltol = nothing, + termination_condition = nothing, internalnorm::F = DEFAULT_NORM, + kwargs...) where {uType, iip, F} alg = get_concrete_algorithm(alg_, prob) @unpack f, u0, p = prob linsolve_with_JᵀJ = Val(_needs_square_A(alg, u0)) u = alias_u0 ? u0 : deepcopy(u0) - if iip - fu1 = f.resid_prototype === nothing ? zero(u) : f.resid_prototype - f(fu1, u, p) - else - fu1 = f(u, p) - end + fu1 = evaluate_f(prob, u) if SciMLBase._unwrap_val(linsolve_with_JᵀJ) uf, linsolve, J, fu2, jac_cache, du, JᵀJ, Jᵀf = jacobian_caches(alg, f, u, p, @@ -109,25 +101,19 @@ function SciMLBase.__init(prob::NonlinearLeastSquaresProblem{uType, iip}, alg_:: JᵀJ, Jᵀf = nothing, nothing end - abstol, reltol, termination_condition = _init_termination_elements(abstol, reltol, - termination_condition, eltype(u); mode = NLSolveTerminationMode.AbsNorm) - - mode = DiffEqBase.get_termination_mode(termination_condition) - - storage = mode ∈ DiffEqBase.SAFE_TERMINATION_MODES ? NLSolveSafeTerminationResult() : - nothing + abstol, reltol, tc_cache_1 = init_termination_cache(abstol, reltol, fu1, u, + termination_condition) + _, _, tc_cache_2 = init_termination_cache(abstol, reltol, fu1, u, termination_condition) return GaussNewtonCache{iip}(f, alg, u, copy(u), fu1, fu2, zero(fu1), du, p, uf, linsolve, J, JᵀJ, Jᵀf, jac_cache, false, maxiters, internalnorm, ReturnCode.Default, - abstol, reltol, prob, NLStats(1, 0, 0, 0, 0), storage, termination_condition) + abstol, reltol, prob, NLStats(1, 0, 0, 0, 0), tc_cache_1, tc_cache_2) end function perform_step!(cache::GaussNewtonCache{true}) - @unpack u, u_prev, fu1, f, p, alg, J, JᵀJ, Jᵀf, linsolve, du, tc_storage = cache + @unpack u, u_prev, fu1, f, p, alg, J, JᵀJ, Jᵀf, linsolve, du = cache jacobian!!(J, cache) - termination_condition = cache.termination_condition(tc_storage) - if JᵀJ !== nothing __matmul!(JᵀJ, J', J) __matmul!(Jᵀf, J', fu1) @@ -145,10 +131,11 @@ function perform_step!(cache::GaussNewtonCache{true}) @. u = u - du f(cache.fu_new, u, p) - (termination_condition(cache.fu_new .- cache.fu1, cache.u, u_prev, cache.abstol, - cache.reltol) || - termination_condition(cache.fu_new, cache.u, u_prev, cache.abstol, cache.reltol)) && - (cache.force_stop = true) + check_and_update!(cache.tc_cache_1, cache, cache.fu_new, cache.u, cache.u_prev) + if !cache.force_stop + cache.fu1 .= cache.fu_new .- cache.fu1 + check_and_update!(cache.tc_cache_2, cache, cache.fu1, cache.u, cache.u_prev) + end @. u_prev = u cache.fu1 .= cache.fu_new @@ -160,9 +147,7 @@ function perform_step!(cache::GaussNewtonCache{true}) end function perform_step!(cache::GaussNewtonCache{false}) - @unpack u, u_prev, fu1, f, p, alg, linsolve, tc_storage = cache - - termination_condition = cache.termination_condition(tc_storage) + @unpack u, u_prev, fu1, f, p, alg, linsolve = cache cache.J = jacobian!!(cache.J, cache) @@ -187,10 +172,13 @@ function perform_step!(cache::GaussNewtonCache{false}) cache.u = @. u - cache.du # `u` might not support mutation cache.fu_new = f(cache.u, p) - termination_condition(cache.fu_new, cache.u, u_prev, cache.abstol, cache.reltol) && - (cache.force_stop = true) + check_and_update!(cache.tc_cache_1, cache, cache.fu_new, cache.u, cache.u_prev) + if !cache.force_stop + cache.fu1 = cache.fu_new .- cache.fu1 + check_and_update!(cache.tc_cache_2, cache, cache.fu1, cache.u, cache.u_prev) + end - cache.u_prev = @. cache.u + cache.u_prev = cache.u cache.fu1 = cache.fu_new cache.stats.nf += 1 cache.stats.njacs += 1 @@ -200,9 +188,8 @@ function perform_step!(cache::GaussNewtonCache{false}) end function SciMLBase.reinit!(cache::GaussNewtonCache{iip}, u0 = cache.u; p = cache.p, - abstol = cache.abstol, reltol = cache.reltol, - termination_condition = cache.termination_condition, - maxiters = cache.maxiters) where {iip} + abstol = cache.abstol, reltol = cache.reltol, maxiters = cache.maxiters, + termination_condition = get_termination_mode(cache.tc_cache)) where {iip} cache.p = p if iip recursivecopy!(cache.u, u0) @@ -213,12 +200,15 @@ function SciMLBase.reinit!(cache::GaussNewtonCache{iip}, u0 = cache.u; p = cache cache.fu1 = cache.f(cache.u, p) end - termination_condition = _get_reinit_termination_condition(cache, abstol, reltol, + abstol, reltol, tc_cache_1 = init_termination_cache(abstol, reltol, cache.fu1, cache.u, + termination_condition) + _, _, tc_cache_2 = init_termination_cache(abstol, reltol, cache.fu1, cache.u, termination_condition) cache.abstol = abstol cache.reltol = reltol - cache.termination_condition = termination_condition + cache.tc_cache_1 = tc_cache_1 + cache.tc_cache_2 = tc_cache_2 cache.maxiters = maxiters cache.stats.nf = 1 cache.stats.nsteps = 1 diff --git a/src/jacobian.jl b/src/jacobian.jl index 676678bb2..41df3c092 100644 --- a/src/jacobian.jl +++ b/src/jacobian.jl @@ -49,9 +49,9 @@ end jacobian!!(::Number, cache) = last(value_derivative(cache.uf, cache.u)) # Build Jacobian Caches -function jacobian_caches(alg::AbstractNonlinearSolveAlgorithm, f, u, p, ::Val{iip}; - linsolve_kwargs = (;), lininit::Val{linsolve_init} = Val(true), - linsolve_with_JᵀJ::Val{needsJᵀJ} = Val(false)) where {iip, needsJᵀJ, linsolve_init} +function jacobian_caches(alg::AbstractNonlinearSolveAlgorithm, f::F, u, p, ::Val{iip}; + linsolve_kwargs = (;), lininit::Val{linsolve_init} = Val(true), + linsolve_with_JᵀJ::Val{needsJᵀJ} = Val(false)) where {iip, needsJᵀJ, linsolve_init, F} uf = JacobianWrapper{iip}(f, p) haslinsolve = hasfield(typeof(alg), :linsolve) @@ -135,9 +135,9 @@ __maybe_symmetric(x::StaticArray) = x __maybe_symmetric(x::SparseArrays.AbstractSparseMatrix) = x ## Special Handling for Scalars -function jacobian_caches(alg::AbstractNonlinearSolveAlgorithm, f, u::Number, p, - ::Val{false}; linsolve_with_JᵀJ::Val{needsJᵀJ} = Val(false), - kwargs...) where {needsJᵀJ} +function jacobian_caches(alg::AbstractNonlinearSolveAlgorithm, f::F, u::Number, p, + ::Val{false}; linsolve_with_JᵀJ::Val{needsJᵀJ} = Val(false), + kwargs...) where {needsJᵀJ, F} # NOTE: Scalar `u` assumes scalar output from `f` uf = JacobianWrapper{false}(f, p) needsJᵀJ && return uf, nothing, u, nothing, nothing, u, u, u diff --git a/src/klement.jl b/src/klement.jl index e60aeee9b..cdb146892 100644 --- a/src/klement.jl +++ b/src/klement.jl @@ -32,7 +32,7 @@ function set_linsolve(alg::GeneralKlement, linsolve) end function GeneralKlement(; max_resets::Int = 5, linsolve = nothing, - linesearch = LineSearch(), precs = DEFAULT_PRECS) + linesearch = LineSearch(), precs = DEFAULT_PRECS) linesearch = linesearch isa LineSearch ? linesearch : LineSearch(; method = linesearch) return GeneralKlement(max_resets, linsolve, precs, linesearch) end @@ -61,17 +61,17 @@ end reltol prob stats::NLStats - lscache - termination_condition - tc_storage + ls_cache + tc_cache end get_fu(cache::GeneralKlementCache) = cache.fu +set_fu!(cache::GeneralKlementCache, fu) = (cache.fu = fu) function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg_::GeneralKlement, args...; - alias_u0 = false, maxiters = 1000, abstol = nothing, reltol = nothing, - termination_condition = nothing, internalnorm = DEFAULT_NORM, - linsolve_kwargs = (;), kwargs...) where {uType, iip} + alias_u0 = false, maxiters = 1000, abstol = nothing, reltol = nothing, + termination_condition = nothing, internalnorm::F = DEFAULT_NORM, + linsolve_kwargs = (;), kwargs...) where {uType, iip, F} @unpack f, u0, p = prob u = alias_u0 ? u0 : deepcopy(u0) fu = evaluate_f(prob, u) @@ -89,36 +89,26 @@ function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg_::GeneralKleme linsolve = __setup_linsolve(J, _vec(fu), _vec(du), p, alg) end - abstol, reltol, termination_condition = _init_termination_elements(abstol, - reltol, - termination_condition, - eltype(u)) - - mode = DiffEqBase.get_termination_mode(termination_condition) - - storage = mode ∈ DiffEqBase.SAFE_TERMINATION_MODES ? NLSolveSafeTerminationResult() : - nothing + abstol, reltol, tc_cache = init_termination_cache(abstol, reltol, fu, u, + termination_condition) return GeneralKlementCache{iip}(f, alg, u, zero(u), fu, zero(fu), du, p, linsolve, J, zero(J), zero(J), _vec(zero(fu)), _vec(zero(fu)), 0, false, maxiters, internalnorm, ReturnCode.Default, abstol, reltol, prob, NLStats(1, 0, 0, 0, 0), - init_linesearch_cache(alg.linesearch, f, u, p, fu, Val(iip)), termination_condition, - storage) + init_linesearch_cache(alg.linesearch, f, u, p, fu, Val(iip)), tc_cache) end function perform_step!(cache::GeneralKlementCache{true}) - @unpack u, u_prev, fu, f, p, alg, J, linsolve, du, tc_storage = cache + @unpack u, u_prev, fu, f, p, alg, J, linsolve, du = cache T = eltype(J) - termination_condition = cache.termination_condition(tc_storage) - singular, fact_done = _try_factorize_and_check_singular!(linsolve, J) if singular if cache.resets == alg.max_resets cache.force_stop = true - cache.retcode = ReturnCode.Unstable + cache.retcode = ReturnCode.ConvergenceFailure return nothing end fact_done = false @@ -129,16 +119,15 @@ function perform_step!(cache::GeneralKlementCache{true}) # u = u - J \ fu linres = dolinsolve(alg.precs, linsolve; A = ifelse(fact_done, nothing, J), - b = -_vec(fu), linu = _vec(du), p, reltol = cache.abstol) + b = _vec(fu), linu = _vec(du), p, reltol = cache.abstol) cache.linsolve = linres.cache # Line Search - α = perform_linesearch!(cache.lscache, u, du) - _axpy!(α, du, u) + α = perform_linesearch!(cache.ls_cache, u, du) + _axpy!(-α, du, u) f(cache.fu2, u, p) - termination_condition(cache.fu2, u, u_prev, cache.abstol, cache.reltol) && - (cache.force_stop = true) + check_and_update!(cache, cache.fu2, cache.u, cache.u_prev) cache.stats.nf += 1 cache.stats.nsolve += 1 cache.stats.nfactors += 1 @@ -146,6 +135,7 @@ function perform_step!(cache::GeneralKlementCache{true}) cache.force_stop && return nothing # Update the Jacobian + cache.du .*= -1 cache.J_cache .= cache.J' .^ 2 cache.Jdu .= _vec(du) .^ 2 mul!(cache.Jᵀ²du, cache.J_cache, cache.Jdu) @@ -165,9 +155,7 @@ function perform_step!(cache::GeneralKlementCache{true}) end function perform_step!(cache::GeneralKlementCache{false}) - @unpack fu, f, p, alg, J, linsolve, tc_storage = cache - - termination_condition = cache.termination_condition(tc_storage) + @unpack fu, f, p, alg, J, linsolve = cache T = eltype(J) @@ -176,7 +164,7 @@ function perform_step!(cache::GeneralKlementCache{false}) if singular if cache.resets == alg.max_resets cache.force_stop = true - cache.retcode = ReturnCode.Unstable + cache.retcode = ReturnCode.ConvergenceFailure return nothing end fact_done = false @@ -186,22 +174,20 @@ function perform_step!(cache::GeneralKlementCache{false}) # u = u - J \ fu if linsolve === nothing - cache.du = -fu / cache.J + cache.du = fu / cache.J else linres = dolinsolve(alg.precs, linsolve; A = ifelse(fact_done, nothing, J), - b = -_vec(fu), linu = _vec(cache.du), p, reltol = cache.abstol) + b = _vec(fu), linu = _vec(cache.du), p, reltol = cache.abstol) cache.linsolve = linres.cache end # Line Search - α = perform_linesearch!(cache.lscache, cache.u, cache.du) - cache.u = @. cache.u + α * cache.du # `u` might not support mutation + α = perform_linesearch!(cache.ls_cache, cache.u, cache.du) + cache.u = @. cache.u - α * cache.du # `u` might not support mutation cache.fu2 = f(cache.u, p) - termination_condition(cache.fu2, cache.u, cache.u_prev, cache.abstol, cache.reltol) && - (cache.force_stop = true) - - cache.u_prev = @. cache.u + check_and_update!(cache, cache.fu2, cache.u, cache.u_prev) + cache.u_prev = cache.u cache.stats.nf += 1 cache.stats.nsolve += 1 cache.stats.nfactors += 1 @@ -209,6 +195,7 @@ function perform_step!(cache::GeneralKlementCache{false}) cache.force_stop && return nothing # Update the Jacobian + cache.du = -cache.du cache.J_cache = cache.J' .^ 2 cache.Jdu = _vec(cache.du) .^ 2 cache.Jᵀ²du = cache.J_cache * cache.Jdu @@ -225,9 +212,8 @@ function perform_step!(cache::GeneralKlementCache{false}) end function SciMLBase.reinit!(cache::GeneralKlementCache{iip}, u0 = cache.u; p = cache.p, - abstol = cache.abstol, reltol = cache.reltol, - termination_condition = cache.termination_condition, - maxiters = cache.maxiters) where {iip} + abstol = cache.abstol, reltol = cache.reltol, maxiters = cache.maxiters, + termination_condition = get_termination_mode(cache.tc_cache)) where {iip} cache.p = p if iip recursivecopy!(cache.u, u0) @@ -238,12 +224,12 @@ function SciMLBase.reinit!(cache::GeneralKlementCache{iip}, u0 = cache.u; p = ca cache.fu = cache.f(cache.u, p) end - termination_condition = _get_reinit_termination_condition(cache, abstol, reltol, + abstol, reltol, tc_cache = init_termination_cache(abstol, reltol, cache.fu, cache.u, termination_condition) cache.abstol = abstol cache.reltol = reltol - cache.termination_condition = termination_condition + cache.tc_cache = tc_cache cache.maxiters = maxiters cache.stats.nf = 1 cache.stats.nsteps = 1 diff --git a/src/lbroyden.jl b/src/lbroyden.jl index db4353b41..62fad090d 100644 --- a/src/lbroyden.jl +++ b/src/lbroyden.jl @@ -25,7 +25,7 @@ An implementation of `LimitedMemoryBroyden` with reseting and line search. end function LimitedMemoryBroyden(; max_resets::Int = 3, linesearch = LineSearch(), - threshold::Int = 10, reset_tolerance = nothing) + threshold::Int = 27, reset_tolerance = nothing) linesearch = linesearch isa LineSearch ? linesearch : LineSearch(; method = linesearch) return LimitedMemoryBroyden(max_resets, threshold, linesearch, reset_tolerance) end @@ -59,17 +59,17 @@ end reset_check prob stats::NLStats - lscache - termination_condition - tc_storage + ls_cache + tc_cache end get_fu(cache::LimitedMemoryBroydenCache) = cache.fu +set_fu!(cache::LimitedMemoryBroydenCache, fu) = (cache.fu = fu) function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg::LimitedMemoryBroyden, - args...; alias_u0 = false, maxiters = 1000, abstol = nothing, reltol = nothing, - termination_condition = nothing, internalnorm = DEFAULT_NORM, - kwargs...) where {uType, iip} + args...; alias_u0 = false, maxiters = 1000, abstol = nothing, reltol = nothing, + termination_condition = nothing, internalnorm::F = DEFAULT_NORM, + kwargs...) where {uType, iip, F} @unpack f, u0, p = prob u = alias_u0 ? u0 : deepcopy(u0) if u isa Number @@ -81,42 +81,31 @@ function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg::LimitedMemory fu = evaluate_f(prob, u) threshold = min(alg.threshold, maxiters) U, Vᵀ = __init_low_rank_jacobian(u, fu, threshold) - du = -fu + du = copy(fu) reset_tolerance = alg.reset_tolerance === nothing ? sqrt(eps(eltype(u))) : alg.reset_tolerance reset_check = x -> abs(x) ≤ reset_tolerance - abstol, reltol, termination_condition = _init_termination_elements(abstol, - reltol, - termination_condition, - eltype(u)) - - mode = DiffEqBase.get_termination_mode(termination_condition) - - storage = mode ∈ DiffEqBase.SAFE_TERMINATION_MODES ? NLSolveSafeTerminationResult() : - nothing + abstol, reltol, tc_cache = init_termination_cache(abstol, reltol, fu, u, + termination_condition) return LimitedMemoryBroydenCache{iip}(f, alg, u, zero(u), du, fu, zero(fu), zero(fu), p, U, Vᵀ, similar(u, threshold), similar(u, 1, threshold), zero(u), zero(u), false, 0, 0, alg.max_resets, maxiters, internalnorm, ReturnCode.Default, abstol, reltol, reset_tolerance, reset_check, prob, NLStats(1, 0, 0, 0, 0), - init_linesearch_cache(alg.linesearch, f, u, p, fu, Val(iip)), termination_condition, - storage) + init_linesearch_cache(alg.linesearch, f, u, p, fu, Val(iip)), tc_cache) end function perform_step!(cache::LimitedMemoryBroydenCache{true}) - @unpack f, p, du, u, tc_storage = cache + @unpack f, p, du, u = cache T = eltype(u) - termination_condition = cache.termination_condition(tc_storage) - - α = perform_linesearch!(cache.lscache, u, du) - _axpy!(α, du, u) + α = perform_linesearch!(cache.ls_cache, u, du) + _axpy!(-α, du, u) f(cache.fu2, u, p) - termination_condition(cache.fu2, cache.u, cache.u_prev, cache.abstol, cache.reltol) && - (cache.force_stop = true) + check_and_update!(cache, cache.fu2, cache.u, cache.u_prev) cache.stats.nf += 1 cache.force_stop && return nothing @@ -128,14 +117,15 @@ function perform_step!(cache::LimitedMemoryBroydenCache{true}) if cache.iterations_since_reset > size(cache.U, 1) && (all(cache.reset_check, du) || all(cache.reset_check, cache.dfu)) if cache.resets ≥ cache.max_resets - cache.retcode = ReturnCode.Unstable + cache.retcode = ReturnCode.ConvergenceFailure cache.force_stop = true return nothing end cache.iterations_since_reset = 0 cache.resets += 1 - cache.du .= -cache.fu + cache.du .= cache.fu else + cache.du .*= -1 idx = min(cache.iterations_since_reset, size(cache.U, 1)) U_part = selectdim(cache.U, 1, 1:idx) Vᵀ_part = selectdim(cache.Vᵀ, 2, 1:idx) @@ -154,7 +144,6 @@ function perform_step!(cache::LimitedMemoryBroydenCache{true}) U_part = selectdim(cache.U, 1, 1:idx) Vᵀ_part = selectdim(cache.Vᵀ, 2, 1:idx) __lbroyden_matvec!(_vec(cache.du), cache.Ux, U_part, Vᵀ_part, _vec(cache.fu2)) - cache.du .*= -1 cache.iterations_since_reset += 1 end @@ -165,18 +154,15 @@ function perform_step!(cache::LimitedMemoryBroydenCache{true}) end function perform_step!(cache::LimitedMemoryBroydenCache{false}) - @unpack f, p, tc_storage = cache - - termination_condition = cache.termination_condition(tc_storage) + @unpack f, p = cache T = eltype(cache.u) - α = perform_linesearch!(cache.lscache, cache.u, cache.du) - cache.u = cache.u .+ α * cache.du + α = perform_linesearch!(cache.ls_cache, cache.u, cache.du) + cache.u = cache.u .- α * cache.du cache.fu2 = f(cache.u, p) - termination_condition(cache.fu2, cache.u, cache.u_prev, cache.abstol, cache.reltol) && - (cache.force_stop = true) + check_and_update!(cache, cache.fu2, cache.u, cache.u_prev) cache.stats.nf += 1 cache.force_stop && return nothing @@ -188,14 +174,15 @@ function perform_step!(cache::LimitedMemoryBroydenCache{false}) if cache.iterations_since_reset > size(cache.U, 1) && (all(cache.reset_check, cache.du) || all(cache.reset_check, cache.dfu)) if cache.resets ≥ cache.max_resets - cache.retcode = ReturnCode.Unstable + cache.retcode = ReturnCode.ConvergenceFailure cache.force_stop = true return nothing end cache.iterations_since_reset = 0 cache.resets += 1 - cache.du = -cache.fu + cache.du = cache.fu else + cache.du = -cache.du idx = min(cache.iterations_since_reset, size(cache.U, 1)) U_part = selectdim(cache.U, 1, 1:idx) Vᵀ_part = selectdim(cache.Vᵀ, 2, 1:idx) @@ -215,7 +202,7 @@ function perform_step!(cache::LimitedMemoryBroydenCache{false}) U_part = selectdim(cache.U, 1, 1:idx) Vᵀ_part = selectdim(cache.Vᵀ, 2, 1:idx) cache.du = _restructure(cache.du, - -__lbroyden_matvec(U_part, Vᵀ_part, _vec(cache.fu2))) + __lbroyden_matvec(U_part, Vᵀ_part, _vec(cache.fu2))) cache.iterations_since_reset += 1 end @@ -226,7 +213,8 @@ function perform_step!(cache::LimitedMemoryBroydenCache{false}) end function SciMLBase.reinit!(cache::LimitedMemoryBroydenCache{iip}, u0 = cache.u; p = cache.p, - abstol = cache.abstol, maxiters = cache.maxiters) where {iip} + termination_condition = get_termination_mode(cache.tc_cache), + abstol = cache.abstol, reltol = cache.reltol, maxiters = cache.maxiters) where {iip} cache.p = p if iip recursivecopy!(cache.u, u0) @@ -236,7 +224,13 @@ function SciMLBase.reinit!(cache::LimitedMemoryBroydenCache{iip}, u0 = cache.u; cache.u = u0 cache.fu = cache.f(cache.u, p) end + + abstol, reltol, tc_cache = init_termination_cache(abstol, reltol, cache.fu, cache.u, + termination_condition) + cache.abstol = abstol + cache.reltol = reltol + cache.tc_cache = tc_cache cache.maxiters = maxiters cache.stats.nf = 1 cache.stats.nsteps = 1 @@ -248,7 +242,7 @@ function SciMLBase.reinit!(cache::LimitedMemoryBroydenCache{iip}, u0 = cache.u; end @views function __lbroyden_matvec!(y::AbstractVector, Ux::AbstractVector, - U::AbstractMatrix, Vᵀ::AbstractMatrix, x::AbstractVector) + U::AbstractMatrix, Vᵀ::AbstractMatrix, x::AbstractVector) # Computes Vᵀ × U × x η = size(U, 1) if η == 0 @@ -267,7 +261,7 @@ end end @views function __lbroyden_rmatvec!(y::AbstractVector, xᵀVᵀ::AbstractMatrix, - U::AbstractMatrix, Vᵀ::AbstractMatrix, x::AbstractVector) + U::AbstractMatrix, Vᵀ::AbstractMatrix, x::AbstractVector) # Computes xᵀ × Vᵀ × U η = size(U, 1) if η == 0 diff --git a/src/levenberg.jl b/src/levenberg.jl index bb9f88bd9..f181fcfcc 100644 --- a/src/levenberg.jl +++ b/src/levenberg.jl @@ -100,10 +100,11 @@ function set_ad(alg::LevenbergMarquardt{CJ}, ad) where {CJ} end function LevenbergMarquardt(; concrete_jac = nothing, linsolve = nothing, - precs = DEFAULT_PRECS, damping_initial::Real = 1.0, damping_increase_factor::Real = 2.0, - damping_decrease_factor::Real = 3.0, finite_diff_step_geodesic::Real = 0.1, - α_geodesic::Real = 0.75, b_uphill::Real = 1.0, min_damping_D::AbstractFloat = 1e-8, - adkwargs...) + precs = DEFAULT_PRECS, damping_initial::Real = 1.0, + damping_increase_factor::Real = 2.0, + damping_decrease_factor::Real = 3.0, finite_diff_step_geodesic::Real = 0.1, + α_geodesic::Real = 0.75, b_uphill::Real = 1.0, min_damping_D::AbstractFloat = 1e-8, + adkwargs...) ad = default_adargs_to_adtype(; adkwargs...) return LevenbergMarquardt{_unwrap_val(concrete_jac)}(ad, linsolve, precs, damping_initial, damping_increase_factor, damping_decrease_factor, @@ -156,16 +157,15 @@ end rhs_tmp J² stats::NLStats - termination_condition - tc_storage + tc_cache_1 + tc_cache_2 end function SciMLBase.__init(prob::Union{NonlinearProblem{uType, iip}, - NonlinearLeastSquaresProblem{uType, iip}}, alg_::LevenbergMarquardt, - args...; alias_u0 = false, maxiters = 1000, abstol = nothing, reltol = nothing, - termination_condition = nothing, - internalnorm = DEFAULT_NORM, - linsolve_kwargs = (;), kwargs...) where {uType, iip} + NonlinearLeastSquaresProblem{uType, iip}}, alg_::LevenbergMarquardt, + args...; alias_u0 = false, maxiters = 1000, abstol = nothing, reltol = nothing, + termination_condition = nothing, internalnorm::F = DEFAULT_NORM, + linsolve_kwargs = (;), kwargs...) where {uType, iip, F} alg = get_concrete_algorithm(alg_, prob) @unpack f, u0, p = prob u = alias_u0 ? u0 : deepcopy(u0) @@ -185,9 +185,6 @@ function SciMLBase.__init(prob::Union{NonlinearProblem{uType, iip}, v = similar(du) end - abstol, reltol, termination_condition = _init_termination_elements(abstol, reltol, - termination_condition, eltype(u); mode = NLSolveTerminationMode.AbsNorm) - λ = convert(eltype(u), alg.damping_initial) λ_factor = convert(eltype(u), alg.damping_increase_factor) damping_increase_factor = convert(eltype(u), alg.damping_increase_factor) @@ -213,10 +210,14 @@ function SciMLBase.__init(prob::Union{NonlinearProblem{uType, iip}, make_new_J = true fu_tmp = zero(fu1) - mode = DiffEqBase.get_termination_mode(termination_condition) - - storage = mode ∈ DiffEqBase.SAFE_TERMINATION_MODES ? NLSolveSafeTerminationResult() : - nothing + abstol, reltol, tc_cache_1 = init_termination_cache(abstol, reltol, fu1, u, + termination_condition) + if prob isa NonlinearLeastSquaresProblem + _, _, tc_cache_2 = init_termination_cache(abstol, reltol, fu1, u, + termination_condition) + else + tc_cache_2 = nothing + end if _unwrap_val(linsolve_with_JᵀJ) mat_tmp = zero(JᵀJ) @@ -231,25 +232,15 @@ function SciMLBase.__init(prob::Union{NonlinearProblem{uType, iip}, end return LevenbergMarquardtCache{iip, !_unwrap_val(linsolve_with_JᵀJ)}(f, alg, u, copy(u), - fu1, - fu2, du, p, uf, linsolve, J, - jac_cache, false, maxiters, internalnorm, ReturnCode.Default, abstol, reltol, prob, - DᵀD, - JᵀJ, λ, λ_factor, damping_increase_factor, damping_decrease_factor, h, α_geodesic, - b_uphill, min_damping_D, v, a, tmp_vec, v_old, loss, δ, loss, make_new_J, fu_tmp, - zero(u), zero(fu1), mat_tmp, rhs_tmp, J², NLStats(1, 0, 0, 0, 0), - termination_condition, storage) + fu1, fu2, du, p, uf, linsolve, J, jac_cache, false, maxiters, internalnorm, + ReturnCode.Default, abstol, reltol, prob, DᵀD, JᵀJ, λ, λ_factor, + damping_increase_factor, damping_decrease_factor, h, α_geodesic, b_uphill, + min_damping_D, v, a, tmp_vec, v_old, loss, δ, loss, make_new_J, fu_tmp, zero(u), + zero(fu1), mat_tmp, rhs_tmp, J², NLStats(1, 0, 0, 0, 0), tc_cache_1, tc_cache_2) end function perform_step!(cache::LevenbergMarquardtCache{true, fastls}) where {fastls} - @unpack fu1, f, make_new_J, tc_storage = cache - - termination_condition = cache.termination_condition(tc_storage) - - if iszero(fu1) - cache.force_stop = true - return nothing - end + @unpack fu1, f, make_new_J = cache if make_new_J jacobian!!(cache.J, cache) @@ -321,13 +312,11 @@ function perform_step!(cache::LevenbergMarquardtCache{true, fastls}) where {fast if (1 - β)^b_uphill * loss ≤ loss_old # Accept step. cache.u .+= δ - if termination_condition(cache.fu_tmp, - cache.u, - u_prev, - cache.abstol, - cache.reltol) - cache.force_stop = true - return nothing + check_and_update!(cache.tc_cache_1, cache, cache.fu_tmp, cache.u, cache.u_prev) + if !cache.force_stop && cache.tc_cache_2 !== nothing + # For NLLS Problems + cache.fu1 .= cache.fu_tmp .- cache.fu1 + check_and_update!(cache.tc_cache_2, cache, cache.fu1, cache.u, cache.u_prev) end cache.fu1 .= cache.fu_tmp _vec(cache.v_old) .= _vec(v) @@ -344,14 +333,7 @@ function perform_step!(cache::LevenbergMarquardtCache{true, fastls}) where {fast end function perform_step!(cache::LevenbergMarquardtCache{false, fastls}) where {fastls} - @unpack fu1, f, make_new_J, tc_storage = cache - - termination_condition = cache.termination_condition(tc_storage) - - if iszero(fu1) - cache.force_stop = true - return nothing - end + @unpack fu1, f, make_new_J = cache if make_new_J cache.J = jacobian!!(cache.J, cache) @@ -423,9 +405,11 @@ function perform_step!(cache::LevenbergMarquardtCache{false, fastls}) where {fas if (1 - β)^b_uphill * loss ≤ loss_old # Accept step. cache.u += δ - if termination_condition(fu_new, cache.u, u_prev, cache.abstol, cache.reltol) - cache.force_stop = true - return nothing + check_and_update!(cache.tc_cache_1, cache, fu_new, cache.u, cache.u_prev) + if !cache.force_stop && cache.tc_cache_2 !== nothing + # For NLLS Problems + cache.fu1 = fu_new .- cache.fu1 + check_and_update!(cache.tc_cache_2, cache, cache.fu1, cache.u, cache.u_prev) end cache.fu1 = fu_new cache.v_old = _restructure(cache.v_old, v) diff --git a/src/linesearch.jl b/src/linesearch.jl index 1f63d695c..760f67769 100644 --- a/src/linesearch.jl +++ b/src/linesearch.jl @@ -8,7 +8,7 @@ differentiation for fast Vector Jacobian Products. ### Arguments - - `method`: the line search algorithm to use. Defaults to `Static()`, which means that the + - `method`: the line search algorithm to use. Defaults to `nothing`, which means that the step size is fixed to the value of `alpha`. - `autodiff`: the automatic differentiation backend to use for the line search. Defaults to `AutoFiniteDiff()`, which means that finite differencing is used to compute the VJP. @@ -22,17 +22,31 @@ differentiation for fast Vector Jacobian Products. α end -function LineSearch(; method = Static(), autodiff = AutoFiniteDiff(), alpha = true) +function LineSearch(; method = nothing, autodiff = AutoFiniteDiff(), alpha = true) return LineSearch(method, autodiff, alpha) end -@inline function init_linesearch_cache(ls::LineSearch, args...) - return init_linesearch_cache(ls.method, ls, args...) +@inline function init_linesearch_cache(ls::LineSearch, f::F, u, p, fu, iip) where {F} + return init_linesearch_cache(ls.method, ls, f, u, p, fu, iip) end +@concrete struct NoLineSearchCache + α +end + +function init_linesearch_cache(::Nothing, ls::LineSearch, f::F, u, p, fu, iip) where {F} + return NoLineSearchCache(convert(eltype(u), ls.α)) +end + +perform_linesearch!(cache::NoLineSearchCache, u, du) = cache.α + # LineSearches.jl doesn't have a supertype so default to that -init_linesearch_cache(_, ls, f, u, p, fu, iip) = LineSearchesJLCache(ls, f, u, p, fu, iip) +function init_linesearch_cache(_, ls::LineSearch, f::F, u, p, fu, iip) where {F} + return LineSearchesJLCache(ls, f, u, p, fu, iip) +end +# FIXME: The closures lead to too many unnecessary runtime dispatches which leads to the +# massive increase in precompilation times. # Wrapper over LineSearches.jl algorithms @concrete mutable struct LineSearchesJLCache f @@ -43,7 +57,7 @@ init_linesearch_cache(_, ls, f, u, p, fu, iip) = LineSearchesJLCache(ls, f, u, p ls end -function LineSearchesJLCache(ls::LineSearch, f, u::Number, p, _, ::Val{false}) +function LineSearchesJLCache(ls::LineSearch, f::F, u::Number, p, _, ::Val{false}) where {F} eval_f(u, du, α) = eval_f(u - α * du) eval_f(u) = f(u, p) @@ -84,7 +98,7 @@ function LineSearchesJLCache(ls::LineSearch, f, u::Number, p, _, ::Val{false}) return LineSearchesJLCache(eval_f, ϕ, dϕ, ϕdϕ, convert(eltype(u), ls.α), ls) end -function LineSearchesJLCache(ls::LineSearch, f, u, p, fu1, IIP::Val{iip}) where {iip} +function LineSearchesJLCache(ls::LineSearch, f::F, u, p, fu1, IIP::Val{iip}) where {iip, F} fu = iip ? deepcopy(fu1) : nothing u_ = _mutable_zero(u) @@ -158,10 +172,7 @@ function perform_linesearch!(cache::LineSearchesJLCache, u, du) ϕ₀, dϕ₀ = ϕdϕ(zero(eltype(u))) - # This case is sometimes possible for large optimization problems - dϕ₀ ≥ 0 && return cache.α - - return first(cache.ls.method(ϕ, cache.dϕ(u, du), cache.ϕdϕ(u, du), cache.α, ϕ₀, dϕ₀)) + return first(cache.ls.method(ϕ, dϕ, ϕdϕ, cache.α, ϕ₀, dϕ₀)) end """ @@ -184,7 +195,7 @@ struct LiFukushimaLineSearch{T} <: AbstractNonlinearSolveLineSearchAlgorithm end function LiFukushimaLineSearch(; lambda_0 = 1.0, beta = 0.1, sigma_1 = 0.001, - sigma_2 = 0.001, eta = 0.1, rho = 0.9, nan_max_iter = 5, maxiters = 50) + sigma_2 = 0.001, eta = 0.1, rho = 0.9, nan_max_iter = 5, maxiters = 50) T = promote_type(typeof(lambda_0), typeof(beta), typeof(sigma_1), typeof(eta), typeof(rho), typeof(sigma_2)) return LiFukushimaLineSearch{T}(lambda_0, beta, sigma_1, sigma_2, eta, rho, @@ -200,8 +211,8 @@ end α end -function init_linesearch_cache(alg::LiFukushimaLineSearch, ls::LineSearch, f, _u, p, _fu, - ::Val{iip}) where {iip} +function init_linesearch_cache(alg::LiFukushimaLineSearch, ls::LineSearch, f::F, _u, p, _fu, + ::Val{iip}) where {iip, F} fu = iip ? deepcopy(_fu) : nothing u = iip ? deepcopy(_u) : nothing return LiFukushimaLineSearchCache{iip}(f, p, u, fu, alg, ls.α) @@ -224,21 +235,21 @@ function perform_linesearch!(cache::LiFukushimaLineSearchCache{iip}, u, du) wher # Early Terminate based on Eq. 2.7 if iip - cache.u_cache .= u .+ du + cache.u_cache .= u .- du cache.f(cache.fu_cache, cache.u_cache, cache.p) fxλ_norm = norm(cache.fu_cache, 2) else - fxλ_norm = norm(cache.f(u .+ du, cache.p), 2) + fxλ_norm = norm(cache.f(u .- du, cache.p), 2) end fxλ_norm ≤ ρ * fx_norm - σ₂ * norm(du, 2)^2 && return cache.α if iip - cache.u_cache .= u .+ λ₂ .* du + cache.u_cache .= u .- λ₂ .* du cache.f(cache.fu_cache, cache.u_cache, cache.p) fxλp_norm = norm(cache.fu_cache, 2) else - fxλp_norm = norm(cache.f(u .+ λ₂ .* du, cache.p), 2) + fxλp_norm = norm(cache.f(u .- λ₂ .* du, cache.p), 2) end if !isfinite(fxλp_norm) @@ -265,11 +276,11 @@ function perform_linesearch!(cache::LiFukushimaLineSearchCache{iip}, u, du) wher for _ in 1:maxiters if iip - cache.u_cache .= u .+ λ₂ .* du + cache.u_cache .= u .- λ₂ .* du cache.f(cache.fu_cache, cache.u_cache, cache.p) fxλp_norm = norm(cache.fu_cache, 2) else - fxλp_norm = norm(cache.f(u .+ λ₂ .* du, cache.p), 2) + fxλp_norm = norm(cache.f(u .- λ₂ .* du, cache.p), 2) end converged = fxλp_norm ≤ (1 + η) * fx_norm - σ₁ * λ₂^2 * norm(du, 2)^2 diff --git a/src/pseudotransient.jl b/src/pseudotransient.jl index ca7cbcdbb..5da1375d6 100644 --- a/src/pseudotransient.jl +++ b/src/pseudotransient.jl @@ -46,7 +46,7 @@ function set_ad(alg::PseudoTransient{CJ}, ad) where {CJ} end function PseudoTransient(; concrete_jac = nothing, linsolve = nothing, - precs = DEFAULT_PRECS, alpha_initial = 1e-3, adkwargs...) + precs = DEFAULT_PRECS, alpha_initial = 1e-3, adkwargs...) ad = default_adargs_to_adtype(; adkwargs...) return PseudoTransient{_unwrap_val(concrete_jac)}(ad, linsolve, precs, alpha_initial) end @@ -74,48 +74,36 @@ end reltol prob stats::NLStats - termination_condition - tc_storage + tc_cache end function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg_::PseudoTransient, - args...; alias_u0 = false, maxiters = 1000, abstol = nothing, reltol = nothing, - termination_condition = nothing, internalnorm = DEFAULT_NORM, - linsolve_kwargs = (;), kwargs...) where {uType, iip} + args...; alias_u0 = false, maxiters = 1000, abstol = nothing, reltol = nothing, + termination_condition = nothing, internalnorm = DEFAULT_NORM, + linsolve_kwargs = (;), kwargs...) where {uType, iip} alg = get_concrete_algorithm(alg_, prob) @unpack f, u0, p = prob u = alias_u0 ? u0 : deepcopy(u0) - if iip - fu1 = f.resid_prototype === nothing ? zero(u) : f.resid_prototype - f(fu1, u, p) - else - fu1 = _mutable(f(u, p)) - end + fu1 = evaluate_f(prob, u) uf, linsolve, J, fu2, jac_cache, du = jacobian_caches(alg, f, u, p, Val(iip); linsolve_kwargs) alpha = convert(eltype(u), alg.alpha_initial) res_norm = internalnorm(fu1) - abstol, reltol, termination_condition = _init_termination_elements(abstol, - reltol, termination_condition, eltype(u)) - - mode = DiffEqBase.get_termination_mode(termination_condition) - - storage = mode ∈ DiffEqBase.SAFE_TERMINATION_MODES ? NLSolveSafeTerminationResult() : - nothing + abstol, reltol, tc_cache = init_termination_cache(abstol, reltol, fu1, u, + termination_condition) return PseudoTransientCache{iip}(f, alg, u, copy(u), fu1, fu2, du, p, alpha, res_norm, - uf, - linsolve, J, jac_cache, false, maxiters, internalnorm, ReturnCode.Default, abstol, - reltol, prob, NLStats(1, 0, 0, 0, 0), termination_condition, storage) + uf, linsolve, J, jac_cache, false, maxiters, internalnorm, ReturnCode.Default, + abstol, reltol, prob, NLStats(1, 0, 0, 0, 0), tc_cache) end function perform_step!(cache::PseudoTransientCache{true}) - @unpack u, u_prev, fu1, f, p, alg, J, linsolve, du, alpha, tc_storage = cache + @unpack u, u_prev, fu1, f, p, alg, J, linsolve, du, alpha = cache jacobian!!(J, cache) - inv_alpha = inv(alpha) + inv_alpha = inv(alpha) if J isa SciMLBase.AbstractSciMLOperator J = J - inv_alpha * I else @@ -129,8 +117,6 @@ function perform_step!(cache::PseudoTransientCache{true}) end end - termination_condition = cache.termination_condition(tc_storage) - # u = u - J \ fu linres = dolinsolve(alg.precs, linsolve; A = J, b = _vec(fu1), linu = _vec(du), p, reltol = cache.abstol) @@ -142,8 +128,7 @@ function perform_step!(cache::PseudoTransientCache{true}) cache.alpha *= cache.res_norm / new_norm cache.res_norm = new_norm - termination_condition(fu1, u, u_prev, cache.abstol, cache.reltol) && - (cache.force_stop = true) + check_and_update!(cache, cache.fu1, cache.u, cache.u_prev) @. u_prev = u cache.stats.nf += 1 @@ -154,20 +139,18 @@ function perform_step!(cache::PseudoTransientCache{true}) end function perform_step!(cache::PseudoTransientCache{false}) - @unpack u, u_prev, fu1, f, p, alg, linsolve, alpha, tc_storage = cache - - tc_storage = cache.tc_storage - termination_condition = cache.termination_condition(tc_storage) + @unpack u, u_prev, fu1, f, p, alg, linsolve, alpha = cache cache.J = jacobian!!(cache.J, cache) - inv_alpha = inv(alpha) + inv_alpha = inv(alpha) cache.J = cache.J - inv_alpha * I # u = u - J \ fu if linsolve === nothing cache.du = fu1 / cache.J else - linres = dolinsolve(alg.precs, linsolve; A = cache.J,b = _vec(fu1),linu = _vec(cache.du), p, reltol = cache.abstol) + linres = dolinsolve(alg.precs, linsolve; A = cache.J, b = _vec(fu1), + linu = _vec(cache.du), p, reltol = cache.abstol) cache.linsolve = linres.cache end cache.u = @. u - cache.du # `u` might not support mutation @@ -176,9 +159,10 @@ function perform_step!(cache::PseudoTransientCache{false}) new_norm = cache.internalnorm(fu1) cache.alpha *= cache.res_norm / new_norm cache.res_norm = new_norm - termination_condition(fu1, cache.u, u_prev, cache.abstol, cache.reltol) && - (cache.force_stop = true) - cache.u_prev = @. cache.u + + check_and_update!(cache, cache.fu1, cache.u, cache.u_prev) + + cache.u_prev = cache.u cache.stats.nf += 1 cache.stats.njacs += 1 cache.stats.nsolve += 1 @@ -187,10 +171,9 @@ function perform_step!(cache::PseudoTransientCache{false}) end function SciMLBase.reinit!(cache::PseudoTransientCache{iip}, u0 = cache.u; p = cache.p, - alpha_new, - abstol = cache.abstol, reltol = cache.reltol, - termination_condition = cache.termination_condition, - maxiters = cache.maxiters) where {iip} + alpha = cache.alpha, abstol = cache.abstol, reltol = cache.reltol, + termination_condition = get_termination_mode(cache.tc_cache), + maxiters = cache.maxiters) where {iip} cache.p = p if iip recursivecopy!(cache.u, u0) @@ -201,14 +184,14 @@ function SciMLBase.reinit!(cache::PseudoTransientCache{iip}, u0 = cache.u; p = c cache.fu1 = cache.f(cache.u, p) end - termination_condition = _get_reinit_termination_condition(cache, abstol, reltol, + abstol, reltol, tc_cache = init_termination_cache(abstol, reltol, cache.fu1, cache.u, termination_condition) - cache.alpha = convert(eltype(cache.u), alpha_new) + cache.alpha = convert(eltype(cache.u), alpha) cache.res_norm = cache.internalnorm(cache.fu1) cache.abstol = abstol cache.reltol = reltol - cache.termination_condition = termination_condition + cache.tc_cache = tc_cache cache.maxiters = maxiters cache.stats.nf = 1 cache.stats.nsteps = 1 diff --git a/src/raphson.jl b/src/raphson.jl index 6e2a502bb..1b75d231e 100644 --- a/src/raphson.jl +++ b/src/raphson.jl @@ -43,13 +43,10 @@ function set_ad(alg::NewtonRaphson{CJ}, ad) where {CJ} end function NewtonRaphson(; concrete_jac = nothing, linsolve = nothing, - linesearch = LineSearch(), precs = DEFAULT_PRECS, adkwargs...) + linesearch = LineSearch(), precs = DEFAULT_PRECS, adkwargs...) ad = default_adargs_to_adtype(; adkwargs...) linesearch = linesearch isa LineSearch ? linesearch : LineSearch(; method = linesearch) - return NewtonRaphson{_unwrap_val(concrete_jac)}(ad, - linsolve, - precs, - linesearch) + return NewtonRaphson{_unwrap_val(concrete_jac)}(ad, linsolve, precs, linesearch) end @concrete mutable struct NewtonRaphsonCache{iip} <: AbstractNonlinearSolveCache{iip} @@ -73,16 +70,14 @@ end reltol prob stats::NLStats - lscache - termination_condition - tc_storage + ls_cache + tc_cache end function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg_::NewtonRaphson, args...; - alias_u0 = false, maxiters = 1000, abstol = nothing, reltol = nothing, - termination_condition = nothing, - internalnorm = DEFAULT_NORM, - linsolve_kwargs = (;), kwargs...) where {uType, iip} + alias_u0 = false, maxiters = 1000, abstol = nothing, reltol = nothing, + termination_condition = nothing, internalnorm = DEFAULT_NORM, + linsolve_kwargs = (;), kwargs...) where {uType, iip} alg = get_concrete_algorithm(alg_, prob) @unpack f, u0, p = prob u = alias_u0 ? u0 : deepcopy(u0) @@ -90,41 +85,30 @@ function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg_::NewtonRaphso uf, linsolve, J, fu2, jac_cache, du = jacobian_caches(alg, f, u, p, Val(iip); linsolve_kwargs) - abstol, reltol, termination_condition = _init_termination_elements(abstol, - reltol, - termination_condition, - eltype(u)) - - mode = DiffEqBase.get_termination_mode(termination_condition) - - storage = mode ∈ DiffEqBase.SAFE_TERMINATION_MODES ? NLSolveSafeTerminationResult() : - nothing + abstol, reltol, tc_cache = init_termination_cache(abstol, reltol, fu1, u, + termination_condition) return NewtonRaphsonCache{iip}(f, alg, u, copy(u), fu1, fu2, du, p, uf, linsolve, J, jac_cache, false, maxiters, internalnorm, ReturnCode.Default, abstol, reltol, prob, NLStats(1, 0, 0, 0, 0), - init_linesearch_cache(alg.linesearch, f, u, p, fu1, Val(iip)), - termination_condition, storage) + init_linesearch_cache(alg.linesearch, f, u, p, fu1, Val(iip)), tc_cache) end function perform_step!(cache::NewtonRaphsonCache{true}) - @unpack u, u_prev, fu1, f, p, alg, J, linsolve, du, tc_storage = cache + @unpack u, u_prev, fu1, f, p, alg, J, linsolve, du = cache jacobian!!(J, cache) - termination_condition = cache.termination_condition(tc_storage) - # u = u - J \ fu linres = dolinsolve(alg.precs, linsolve; A = J, b = _vec(fu1), linu = _vec(du), p, reltol = cache.abstol) cache.linsolve = linres.cache # Line Search - α = perform_linesearch!(cache.lscache, u, du) + α = perform_linesearch!(cache.ls_cache, u, du) _axpy!(-α, du, u) f(cache.fu1, u, p) - termination_condition(cache.fu1, u, u_prev, cache.abstol, cache.reltol) && - (cache.force_stop = true) + check_and_update!(cache, cache.fu1, cache.u, cache.u_prev) @. u_prev = u cache.stats.nf += 1 @@ -135,9 +119,7 @@ function perform_step!(cache::NewtonRaphsonCache{true}) end function perform_step!(cache::NewtonRaphsonCache{false}) - @unpack u, u_prev, fu1, f, p, alg, linsolve, tc_storage = cache - - termination_condition = cache.termination_condition(tc_storage) + @unpack u, u_prev, fu1, f, p, alg, linsolve = cache cache.J = jacobian!!(cache.J, cache) # u = u - J \ fu @@ -150,14 +132,13 @@ function perform_step!(cache::NewtonRaphsonCache{false}) end # Line Search - α = perform_linesearch!(cache.lscache, u, cache.du) + α = perform_linesearch!(cache.ls_cache, u, cache.du) cache.u = @. u - α * cache.du # `u` might not support mutation cache.fu1 = f(cache.u, p) - termination_condition(cache.fu1, cache.u, u_prev, cache.abstol, cache.reltol) && - (cache.force_stop = true) + check_and_update!(cache, cache.fu1, cache.u, cache.u_prev) - cache.u_prev = @. cache.u + cache.u_prev = cache.u cache.stats.nf += 1 cache.stats.njacs += 1 cache.stats.nsolve += 1 @@ -166,9 +147,8 @@ function perform_step!(cache::NewtonRaphsonCache{false}) end function SciMLBase.reinit!(cache::NewtonRaphsonCache{iip}, u0 = cache.u; p = cache.p, - abstol = cache.abstol, reltol = cache.reltol, - termination_condition = cache.termination_condition, - maxiters = cache.maxiters) where {iip} + abstol = cache.abstol, reltol = cache.reltol, maxiters = cache.maxiters, + termination_condition = get_termination_mode(cache.tc_cache)) where {iip} cache.p = p if iip recursivecopy!(cache.u, u0) @@ -179,12 +159,12 @@ function SciMLBase.reinit!(cache::NewtonRaphsonCache{iip}, u0 = cache.u; p = cac cache.fu1 = cache.f(cache.u, p) end - termination_condition = _get_reinit_termination_condition(cache, abstol, reltol, + abstol, reltol, tc_cache = init_termination_cache(abstol, reltol, cache.fu1, cache.u, termination_condition) cache.abstol = abstol cache.reltol = reltol - cache.termination_condition = termination_condition + cache.tc_cache = tc_cache cache.maxiters = maxiters cache.stats.nf = 1 cache.stats.nsteps = 1 diff --git a/src/trustRegion.jl b/src/trustRegion.jl index cf9f41af0..0070694ff 100644 --- a/src/trustRegion.jl +++ b/src/trustRegion.jl @@ -141,11 +141,6 @@ for large-scale and numerically-difficult nonlinear systems. `expand_threshold < r` (with `r` defined in `shrink_threshold`). Defaults to `2.0`. - `max_shrink_times`: the maximum number of times to shrink the trust region radius in a row, `max_shrink_times` is exceeded, the algorithm returns. Defaults to `32`. - -!!! warning - - `linsolve` and `precs` are used exclusively for the inplace version of the algorithm. - Support for the OOP version is planned! """ @concrete struct TrustRegion{CJ, AD, MTR} <: AbstractNewtonAlgorithm{CJ, AD} @@ -171,11 +166,11 @@ function set_ad(alg::TrustRegion{CJ}, ad) where {CJ} end function TrustRegion(; concrete_jac = nothing, linsolve = nothing, precs = DEFAULT_PRECS, - radius_update_scheme::RadiusUpdateSchemes.T = RadiusUpdateSchemes.Simple, #defaults to conventional radius update - max_trust_radius::Real = 0 // 1, initial_trust_radius::Real = 0 // 1, - step_threshold::Real = 1 // 10000, shrink_threshold::Real = 1 // 4, - expand_threshold::Real = 3 // 4, shrink_factor::Real = 1 // 4, - expand_factor::Real = 2 // 1, max_shrink_times::Int = 32, adkwargs...) + radius_update_scheme::RadiusUpdateSchemes.T = RadiusUpdateSchemes.Simple, #defaults to conventional radius update + max_trust_radius::Real = 0 // 1, initial_trust_radius::Real = 0 // 1, + step_threshold::Real = 1 // 10000, shrink_threshold::Real = 1 // 4, + expand_threshold::Real = 3 // 4, shrink_factor::Real = 1 // 4, + expand_factor::Real = 2 // 1, max_shrink_times::Int = 32, adkwargs...) ad = default_adargs_to_adtype(; adkwargs...) return TrustRegion{_unwrap_val(concrete_jac)}(ad, linsolve, precs, radius_update_scheme, max_trust_radius, initial_trust_radius, step_threshold, shrink_threshold, @@ -229,15 +224,13 @@ end p4::floatType ϵ::floatType stats::NLStats - tc_storage - termination_condition + tc_cache end function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg_::TrustRegion, args...; - alias_u0 = false, maxiters = 1000, abstol = nothing, reltol = nothing, - termination_condition = nothing, - internalnorm = DEFAULT_NORM, - linsolve_kwargs = (;), kwargs...) where {uType, iip} + alias_u0 = false, maxiters = 1000, abstol = nothing, reltol = nothing, + termination_condition = nothing, internalnorm = DEFAULT_NORM, linsolve_kwargs = (;), + kwargs...) where {uType, iip} alg = get_concrete_algorithm(alg_, prob) @unpack f, u0, p = prob u = alias_u0 ? u0 : deepcopy(u0) @@ -250,7 +243,7 @@ function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg_::TrustRegion, linsolve_kwargs) u_tmp = zero(u) u_cauchy = zero(u) - u_gauss_newton = zero(u) + u_gauss_newton = _mutable_zero(u) loss_new = loss H = zero(J' * J) @@ -338,23 +331,15 @@ function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg_::TrustRegion, initial_trust_radius = convert(trustType, 1.0) end - abstol, reltol, termination_condition = _init_termination_elements(abstol, - reltol, - termination_condition, - eltype(u)) - - mode = DiffEqBase.get_termination_mode(termination_condition) - - storage = mode ∈ DiffEqBase.SAFE_TERMINATION_MODES ? NLSolveSafeTerminationResult() : - nothing + abstol, reltol, tc_cache = init_termination_cache(abstol, reltol, fu1, u, + termination_condition) return TrustRegionCache{iip}(f, alg, u_prev, u, fu_prev, fu1, fu2, p, uf, linsolve, J, jac_cache, false, maxiters, internalnorm, ReturnCode.Default, abstol, reltol, 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, du, u_tmp, u_gauss_newton, u_cauchy, fu_new, make_new_J, r, - p1, p2, p3, p4, ϵ, - NLStats(1, 0, 0, 0, 0), storage, termination_condition) + p1, p2, p3, p4, ϵ, NLStats(1, 0, 0, 0, 0), tc_cache) end function perform_step!(cache::TrustRegionCache{true}) @@ -368,8 +353,7 @@ function perform_step!(cache::TrustRegionCache{true}) # do not use A = cache.H, b = _vec(cache.g) since it is equivalent # to A = cache.J, b = _vec(fu) as long as the Jacobian is non-singular linres = dolinsolve(alg.precs, linsolve, A = J, b = _vec(fu), - linu = _vec(u_gauss_newton), - p = p, reltol = cache.abstol) + linu = _vec(u_gauss_newton), p = p, reltol = cache.abstol) cache.linsolve = linres.cache @. cache.u_gauss_newton = -1 * u_gauss_newton end @@ -395,7 +379,18 @@ function perform_step!(cache::TrustRegionCache{false}) cache.H = J' * J cache.g = _restructure(fu, J' * _vec(fu)) cache.stats.njacs += 1 - cache.u_gauss_newton = -1 .* _restructure(cache.g, cache.H \ _vec(cache.g)) + + if cache.linsolve === nothing + # Scalar + cache.u_gauss_newton = -cache.H \ cache.g + else + # do not use A = cache.H, b = _vec(cache.g) since it is equivalent + # to A = cache.J, b = _vec(fu) as long as the Jacobian is non-singular + linres = dolinsolve(cache.alg.precs, cache.linsolve, A = cache.J, b = _vec(fu), + linu = _vec(cache.u_gauss_newton), p = p, reltol = cache.abstol) + cache.linsolve = linres.cache + @. cache.u_gauss_newton *= -1 + end end # Compute the Newton step. @@ -430,9 +425,7 @@ function retrospective_step!(cache::TrustRegionCache) end function trust_region_step!(cache::TrustRegionCache) - @unpack fu_new, du, g, H, loss, max_trust_r, radius_update_scheme, tc_storage = cache - - termination_condition = cache.termination_condition(tc_storage) + @unpack fu_new, du, g, H, loss, max_trust_r, radius_update_scheme = cache cache.loss_new = get_loss(fu_new) @@ -463,13 +456,7 @@ function trust_region_step!(cache::TrustRegionCache) # 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) || termination_condition(cache.fu, - cache.u, - cache.u_prev, - cache.abstol, - cache.reltol) - cache.force_stop = true - end + check_and_update!(cache, cache.fu, cache.u, cache.u_prev) elseif radius_update_scheme === RadiusUpdateSchemes.NLsolve # accept/reject decision @@ -491,9 +478,7 @@ function trust_region_step!(cache::TrustRegionCache) end # convergence test - if iszero(cache.fu) || cache.internalnorm(cache.fu) < cache.abstol - cache.force_stop = true - end + check_and_update!(cache, cache.fu, cache.u, cache.u_prev) elseif radius_update_scheme === RadiusUpdateSchemes.NocedalWright # accept/reject decision @@ -512,9 +497,7 @@ function trust_region_step!(cache::TrustRegionCache) end # convergence test - if iszero(cache.fu) || cache.internalnorm(cache.fu) < cache.abstol - cache.force_stop = true - end + check_and_update!(cache, cache.fu, cache.u, cache.u_prev) elseif radius_update_scheme === RadiusUpdateSchemes.Hei if r > cache.step_threshold @@ -535,15 +518,8 @@ function trust_region_step!(cache::TrustRegionCache) cache.trust_r = rfunc(r, shrink_threshold, p1, p3, p4, p2) * cache.internalnorm(du) - if iszero(cache.fu) || - termination_condition(cache.fu, - cache.u, - cache.u_prev, - cache.abstol, - cache.reltol) || - cache.internalnorm(g) < cache.ϵ - cache.force_stop = true - end + check_and_update!(cache, cache.fu, cache.u, cache.u_prev) + cache.internalnorm(g) < cache.ϵ && (cache.force_stop = true) elseif radius_update_scheme === RadiusUpdateSchemes.Yuan if r < cache.shrink_threshold @@ -565,15 +541,8 @@ function trust_region_step!(cache::TrustRegionCache) @unpack p1 = cache cache.trust_r = p1 * cache.internalnorm(jvp!(cache)) - if iszero(cache.fu) || - termination_condition(cache.fu, - cache.u, - cache.u_prev, - cache.abstol, - cache.reltol) || - cache.internalnorm(g) < cache.ϵ - cache.force_stop = true - end + check_and_update!(cache, cache.fu, cache.u, cache.u_prev) + cache.internalnorm(g) < cache.ϵ && (cache.force_stop = true) #Fan's update scheme elseif radius_update_scheme === RadiusUpdateSchemes.Fan if r < cache.shrink_threshold @@ -594,15 +563,8 @@ function trust_region_step!(cache::TrustRegionCache) @unpack p1 = cache cache.trust_r = p1 * (cache.internalnorm(cache.fu)^0.99) - if iszero(cache.fu) || - termination_condition(cache.fu, - cache.u, - cache.u_prev, - cache.abstol, - cache.reltol) || - cache.internalnorm(g) < cache.ϵ - cache.force_stop = true - end + check_and_update!(cache, cache.fu, cache.u, cache.u_prev) + cache.internalnorm(g) < cache.ϵ && (cache.force_stop = true) elseif radius_update_scheme === RadiusUpdateSchemes.Bastin if r > cache.step_threshold take_step!(cache) @@ -617,13 +579,7 @@ function trust_region_step!(cache::TrustRegionCache) cache.trust_r *= cache.p2 cache.shrink_counter += 1 end - if iszero(cache.fu) || termination_condition(cache.fu, - cache.u, - cache.u_prev, - cache.abstol, - cache.reltol) - cache.force_stop = true - end + check_and_update!(cache, cache.fu, cache.u, cache.u_prev) end end @@ -718,15 +674,23 @@ function jvp!(cache::TrustRegionCache{true}) end function not_terminated(cache::TrustRegionCache) - return !cache.force_stop && cache.stats.nsteps < cache.maxiters && - cache.shrink_counter < cache.alg.max_shrink_times + non_shrink_terminated = cache.force_stop || cache.stats.nsteps ≥ cache.maxiters + # Terminated due to convergence or maxiters + non_shrink_terminated && return false + # Terminated due to too many shrink + shrink_terminated = cache.shrink_counter ≥ cache.alg.max_shrink_times + if shrink_terminated + cache.retcode = ReturnCode.ConvergenceFailure + return false + end + return true end get_fu(cache::TrustRegionCache) = cache.fu +set_fu!(cache::TrustRegionCache, fu) = (cache.fu = fu) function SciMLBase.reinit!(cache::TrustRegionCache{iip}, u0 = cache.u; p = cache.p, - abstol = cache.abstol, reltol = cache.reltol, - termination_condition = cache.termination_condition, - maxiters = cache.maxiters) where {iip} + abstol = cache.abstol, reltol = cache.reltol, maxiters = cache.maxiters, + termination_condition = get_termination_mode(cache.tc_cache)) where {iip} cache.p = p if iip recursivecopy!(cache.u, u0) @@ -737,12 +701,12 @@ function SciMLBase.reinit!(cache::TrustRegionCache{iip}, u0 = cache.u; p = cache cache.fu = cache.f(cache.u, p) end - termination_condition = _get_reinit_termination_condition(cache, abstol, reltol, + abstol, reltol, tc_cache = init_termination_cache(abstol, reltol, cache.fu, cache.u, termination_condition) cache.abstol = abstol cache.reltol = reltol - cache.termination_condition = termination_condition + cache.tc_cache = tc_cache cache.maxiters = maxiters cache.stats.nf = 1 cache.stats.nsteps = 1 diff --git a/src/utils.jl b/src/utils.jl index 87a80e4ed..0b0447772 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -1,17 +1,19 @@ +const DEFAULT_NORM = DiffEqBase.NONLINEARSOLVE_DEFAULT_NORM -@inline UNITLESS_ABS2(x) = real(abs2(x)) -@inline DEFAULT_NORM(u::Union{AbstractFloat, Complex}) = @fastmath abs(u) -@inline function DEFAULT_NORM(u::Array{T}) where {T <: Union{AbstractFloat, Complex}} - return sqrt(real(sum(abs2, u)) / length(u)) -end -@inline function DEFAULT_NORM(u::StaticArray{<:Union{AbstractFloat, Complex}}) - return sqrt(real(sum(abs2, u)) / length(u)) +# Ignores NaN +function __findmin(f, x) + return findmin(x) do xᵢ + fx = f(xᵢ) + return isnan(fx) ? Inf : fx + end end -@inline function DEFAULT_NORM(u::AbstractVectorOfArray) - return sum(sqrt(real(sum(UNITLESS_ABS2, _u)) / length(_u)) for _u in u.u) + +struct NonlinearSolveTag end + +function ForwardDiff.checktag(::Type{<:ForwardDiff.Tag{<:NonlinearSolveTag, <:T}}, f::F, + x::AbstractArray{T}) where {T, F} + return true end -@inline DEFAULT_NORM(u::AbstractArray) = sqrt(real(sum(UNITLESS_ABS2, u)) / length(u)) -@inline DEFAULT_NORM(u) = norm(u) """ default_adargs_to_adtype(; chunk_size = Val{0}(), autodiff = Val{true}(), @@ -26,7 +28,7 @@ code. `autodiff=`. """ function default_adargs_to_adtype(; chunk_size = missing, autodiff = nothing, - standardtag = missing, diff_type = missing) + standardtag = missing, diff_type = missing) # If using the new API short circuit autodiff === nothing && return nothing autodiff isa ADTypes.AbstractADType && return autodiff @@ -45,8 +47,9 @@ function default_adargs_to_adtype(; chunk_size = missing, autodiff = nothing, autodiff === missing && (autodiff = Val{true}()) ad = _unwrap_val(autodiff) - tag = _unwrap_val(standardtag) - ad && return AutoForwardDiff{_unwrap_val(chunk_size), typeof(tag)}(tag) + # We don't really know the typeof the input yet, so we can't use the correct tag! + ad && + return AutoForwardDiff{_unwrap_val(chunk_size), NonlinearSolveTag}(NonlinearSolveTag()) return AutoFiniteDiff(; fdtype = diff_type) end @@ -79,8 +82,8 @@ end DEFAULT_PRECS(W, du, u, p, t, newW, Plprev, Prprev, cachedata) = nothing, nothing function dolinsolve(precs::P, linsolve; A = nothing, linu = nothing, b = nothing, - du = nothing, u = nothing, p = nothing, t = nothing, weight = nothing, - cachedata = nothing, reltol = nothing) where {P} + du = nothing, u = nothing, p = nothing, t = nothing, weight = nothing, + cachedata = nothing, reltol = nothing) where {P} A !== nothing && (linsolve.A = A) b !== nothing && (linsolve.b = b) linu !== nothing && (linsolve.u = linu) @@ -120,17 +123,6 @@ function wrapprecs(_Pl, _Pr, weight) return Pl, Pr end -function _nfcount(N, ::Type{diff_type}) where {diff_type} - if diff_type === Val{:complex} - tmp = N - elseif diff_type === Val{:forward} - tmp = N + 1 - else - tmp = 2N - end - return tmp -end - get_loss(fu) = norm(fu)^2 / 2 function rfunc(r::R, c2::R, M::R, γ1::R, γ2::R, β::R) where {R <: Real} # R-function for adaptive trust region method @@ -157,7 +149,7 @@ _maybe_mutable(x, _) = x # Helper function to get value of `f(u, p)` function evaluate_f(prob::Union{NonlinearProblem{uType, iip}, - NonlinearLeastSquaresProblem{uType, iip}}, u) where {uType, iip} + NonlinearLeastSquaresProblem{uType, iip}}, u) where {uType, iip} @unpack f, u0, p = prob if iip fu = f.resid_prototype === nothing ? zero(u) : f.resid_prototype @@ -205,63 +197,80 @@ function __get_concrete_algorithm(alg, prob) # Use Finite Differencing use_sparse_ad ? AutoSparseFiniteDiff() : AutoFiniteDiff() else - use_sparse_ad ? AutoSparseForwardDiff() : AutoForwardDiff{nothing, Nothing}(nothing) + (use_sparse_ad ? AutoSparseForwardDiff : AutoForwardDiff)(; + tag = NonlinearSolveTag()) end return set_ad(alg, ad) end -function _get_tolerance(η, tc_η, ::Type{T}) where {T} - fallback_η = real(oneunit(T)) * (eps(real(one(T))))^(4 // 5) - return T(ifelse(η !== nothing, η, ifelse(tc_η !== nothing, tc_η, fallback_η))) +function init_termination_cache(abstol, reltol, du, u, ::Nothing) + return init_termination_cache(abstol, reltol, du, u, AbsSafeBestTerminationMode()) +end +function init_termination_cache(abstol, reltol, du, u, tc::AbstractNonlinearTerminationMode) + tc_cache = init(du, u, tc; abstol, reltol) + return DiffEqBase.get_abstol(tc_cache), DiffEqBase.get_reltol(tc_cache), tc_cache end -function _init_termination_elements(abstol, reltol, termination_condition, - ::Type{T}; mode = NLSolveTerminationMode.AbsNorm) where {T} - if termination_condition !== nothing - if abstol !== nothing && abstol != termination_condition.abstol - error("Incompatible absolute tolerances found. The tolerances supplied as the \ - keyword argument and the one supplied in the termination condition should \ - be same.") +function check_and_update!(cache, fu, u, uprev) + return check_and_update!(cache.tc_cache, cache, fu, u, uprev) +end +function check_and_update!(tc_cache, cache, fu, u, uprev) + return check_and_update!(tc_cache, cache, fu, u, uprev, + DiffEqBase.get_termination_mode(tc_cache)) +end +function check_and_update!(tc_cache, cache, fu, u, uprev, + mode::AbstractNonlinearTerminationMode) + if tc_cache(fu, u, uprev) + # Just a sanity measure! + if isinplace(cache) + cache.prob.f(get_fu(cache), u, cache.prob.p) + else + set_fu!(cache, cache.prob.f(u, cache.prob.p)) end - if reltol !== nothing && reltol != termination_condition.reltol - error("Incompatible relative tolerances found. The tolerances supplied as the \ - keyword argument and the one supplied in the termination condition should \ - be same.") + cache.force_stop = true + end +end +function check_and_update!(tc_cache, cache, fu, u, uprev, + mode::AbstractSafeNonlinearTerminationMode) + if tc_cache(fu, u, uprev) + if tc_cache.retcode == NonlinearSafeTerminationReturnCode.Success + cache.retcode = ReturnCode.Success end - abstol = _get_tolerance(abstol, termination_condition.abstol, T) - reltol = _get_tolerance(reltol, termination_condition.reltol, T) - return abstol, reltol, termination_condition - else - abstol = _get_tolerance(abstol, nothing, T) - reltol = _get_tolerance(reltol, nothing, T) - termination_condition = NLSolveTerminationCondition(mode; abstol, reltol) - return abstol, reltol, termination_condition + if tc_cache.retcode == NonlinearSafeTerminationReturnCode.PatienceTermination + cache.retcode = ReturnCode.ConvergenceFailure + end + if tc_cache.retcode == NonlinearSafeTerminationReturnCode.ProtectiveTermination + cache.retcode = ReturnCode.Unstable + end + # Just a sanity measure! + if isinplace(cache) + cache.prob.f(get_fu(cache), u, cache.prob.p) + else + set_fu!(cache, cache.prob.f(u, cache.prob.p)) + end + cache.force_stop = true end end - -function _get_reinit_termination_condition(cache, abstol, reltol, termination_condition) - if termination_condition != cache.termination_condition - if abstol != cache.abstol && abstol != termination_condition.abstol - error("Incompatible absolute tolerances found. The tolerances supplied as the \ - keyword argument and the one supplied in the termination condition \ - should be same.") +function check_and_update!(tc_cache, cache, fu, u, uprev, + mode::AbstractSafeBestNonlinearTerminationMode) + if tc_cache(fu, u, uprev) + if tc_cache.retcode == NonlinearSafeTerminationReturnCode.Success + cache.retcode = ReturnCode.Success end - - if reltol != cache.reltol && reltol != termination_condition.reltol - error("Incompatible absolute tolerances found. The tolerances supplied as the \ - keyword argument and the one supplied in the termination condition \ - should be same.") + if tc_cache.retcode == NonlinearSafeTerminationReturnCode.PatienceTermination + cache.retcode = ReturnCode.ConvergenceFailure end - return termination_condition - else - # Build the termination_condition with new abstol and reltol - return NLSolveTerminationCondition{ - DiffEqBase.get_termination_mode(termination_condition), - eltype(abstol), - typeof(termination_condition.safe_termination_options), - }(abstol, - reltol, - termination_condition.safe_termination_options) + if tc_cache.retcode == NonlinearSafeTerminationReturnCode.ProtectiveTermination + cache.retcode = ReturnCode.Unstable + end + if isinplace(cache) + copyto!(get_u(cache), tc_cache.u) + cache.prob.f(get_fu(cache), get_u(cache), cache.prob.p) + else + set_u!(cache, tc_cache.u) + set_fu!(cache, cache.prob.f(get_u(cache), cache.prob.p)) + end + cache.force_stop = true end end diff --git a/test/23_test_problems.jl b/test/23_test_problems.jl index 77274d34a..5ec5e611f 100644 --- a/test/23_test_problems.jl +++ b/test/23_test_problems.jl @@ -11,7 +11,8 @@ function test_on_library(problems, dicts, alg_ops, broken_tests, ϵ = 1e-4) @testset "$idx: $(dict["title"])" begin for alg in alg_ops try - sol = solve(nlprob, alg) + sol = solve(nlprob, alg; + termination_condition = AbsNormTerminationMode()) problem(res, sol.u, nothing) broken = idx in broken_tests[alg] ? true : false @@ -52,7 +53,7 @@ end broken_tests[alg_ops[1]] = [6, 11, 21] broken_tests[alg_ops[2]] = [6, 11, 21] broken_tests[alg_ops[3]] = [1, 6, 11, 12, 15, 16, 21] - broken_tests[alg_ops[4]] = [1, 6, 8, 11, 16, 21, 22] + broken_tests[alg_ops[4]] = [1, 6, 8, 11, 15, 16, 21, 22] broken_tests[alg_ops[5]] = [6, 21] broken_tests[alg_ops[6]] = [6, 21] @@ -76,39 +77,27 @@ end alg_ops = (DFSane(),) broken_tests = Dict(alg => Int[] for alg in alg_ops) - broken_tests[alg_ops[1]] = [1, 2, 3, 5, 6, 21] + broken_tests[alg_ops[1]] = [1, 2, 3, 4, 5, 6, 11, 22] test_on_library(problems, dicts, alg_ops, broken_tests) end # Broyden and Klement Tests are quite flaky and failure seems to be platform dependent # needs additional investigation before we can enable them -#= @testset "GeneralBroyden 23 Test Problems" begin - alg_ops = (GeneralBroyden(), - GeneralBroyden(; linesearch = LiFukushimaLineSearch(; beta = 0.1)), - GeneralBroyden(; linesearch = BackTracking())) + alg_ops = (GeneralBroyden(),) broken_tests = Dict(alg => Int[] for alg in alg_ops) - broken_tests[alg_ops[1]] = [1, 3, 4, 5, 6, 11, 12, 13, 14, 21] - broken_tests[alg_ops[2]] = [1, 2, 3, 4, 5, 6, 9, 11, 13, 22] - broken_tests[alg_ops[3]] = [1, 2, 4, 5, 6, 11, 12, 13, 14, 21] + broken_tests[alg_ops[1]] = [1, 2, 4, 5, 6, 11, 12, 13, 14] test_on_library(problems, dicts, alg_ops, broken_tests) end @testset "GeneralKlement 23 Test Problems" begin - alg_ops = (GeneralKlement(), - GeneralKlement(; linesearch = BackTracking()), - GeneralKlement(; linesearch = HagerZhang())) + alg_ops = (GeneralKlement(),) broken_tests = Dict(alg => Int[] for alg in alg_ops) broken_tests[alg_ops[1]] = [1, 2, 4, 5, 6, 7, 11, 13, 22] - broken_tests[alg_ops[2]] = [1, 2, 4, 5, 6, 7, 11, 12, 13, 22] - broken_tests[alg_ops[3]] = [1, 2, 4, 5, 6, 8, 11, 12, 13, 22] test_on_library(problems, dicts, alg_ops, broken_tests) end -=# - -# NOTE: Not adding LimitedMemoryBroyden here since it fails on most of the preoblems diff --git a/test/basictests.jl b/test/basictests.jl index eaa48d05e..bff4bcbae 100644 --- a/test/basictests.jl +++ b/test/basictests.jl @@ -18,6 +18,13 @@ function newton_fails(u, p) 0.0011552453009332421u .- p end +const TERMINATION_CONDITIONS = [ + SteadyStateDiffEqTerminationMode(), SimpleNonlinearSolveTerminationMode(), + NormTerminationMode(), RelTerminationMode(), RelNormTerminationMode(), + AbsTerminationMode(), AbsNormTerminationMode(), RelSafeTerminationMode(), + AbsSafeTerminationMode(), RelSafeBestTerminationMode(), AbsSafeBestTerminationMode(), +] + # --- NewtonRaphson tests --- @testset "NewtonRaphson" begin @@ -27,7 +34,7 @@ end end function benchmark_nlsolve_iip(f, u0, p = 2.0; linsolve, precs, - linesearch = LineSearch()) + linesearch = LineSearch()) prob = NonlinearProblem{true}(f, u0, p) return solve(prob, NewtonRaphson(; linsolve, precs, linesearch), abstol = 1e-9) end @@ -124,16 +131,9 @@ end @test all(solve(probN, NewtonRaphson(; autodiff)).u .≈ sqrt(2.0)) end - @testset "Termination condition: $(mode) u0: $(_nameof(u0))" for mode in instances(NLSolveTerminationMode.T), + @testset "Termination condition: $(termination_condition) u0: $(_nameof(u0))" for termination_condition in TERMINATION_CONDITIONS, u0 in (1.0, [1.0, 1.0]) - if mode ∈ - (NLSolveTerminationMode.SteadyStateDefault, NLSolveTerminationMode.RelSafeBest, - NLSolveTerminationMode.AbsSafeBest) - continue - end - termination_condition = NLSolveTerminationCondition(mode; abstol = nothing, - reltol = nothing) probN = NonlinearProblem(quadratic_f, u0, 2.0) @test all(solve(probN, NewtonRaphson(); termination_condition).u .≈ sqrt(2.0)) end @@ -296,16 +296,9 @@ end end end - @testset "Termination condition: $(mode) u0: $(_nameof(u0))" for mode in instances(NLSolveTerminationMode.T), + @testset "Termination condition: $(termination_condition) u0: $(_nameof(u0))" for termination_condition in TERMINATION_CONDITIONS, u0 in (1.0, [1.0, 1.0]) - if mode ∈ - (NLSolveTerminationMode.SteadyStateDefault, NLSolveTerminationMode.RelSafeBest, - NLSolveTerminationMode.AbsSafeBest) - continue - end - termination_condition = NLSolveTerminationCondition(mode; abstol = nothing, - reltol = nothing) probN = NonlinearProblem(quadratic_f, u0, 2.0) @test all(solve(probN, TrustRegion(); termination_condition).u .≈ sqrt(2.0)) end @@ -420,16 +413,9 @@ end end end - @testset "Termination condition: $(mode) u0: $(_nameof(u0))" for mode in instances(NLSolveTerminationMode.T), + @testset "Termination condition: $(termination_condition) u0: $(_nameof(u0))" for termination_condition in TERMINATION_CONDITIONS, u0 in (1.0, [1.0, 1.0]) - if mode ∈ - (NLSolveTerminationMode.SteadyStateDefault, NLSolveTerminationMode.RelSafeBest, - NLSolveTerminationMode.AbsSafeBest) - continue - end - termination_condition = NLSolveTerminationCondition(mode; abstol = nothing, - reltol = nothing) probN = NonlinearProblem(quadratic_f, u0, 2.0) @test all(solve(probN, LevenbergMarquardt(); termination_condition).u .≈ sqrt(2.0)) end @@ -470,7 +456,6 @@ end end @testset "[OOP] [Immutable AD]" begin - broken_forwarddiff = [2.9, 3.0, 4.0, 81.0] for p in 1.1:0.1:100.0 res = abs.(benchmark_nlsolve_oop(quadratic_f, @SVector[1.0, 1.0], p).u) @@ -478,9 +463,6 @@ end @test_broken all(res .≈ sqrt(p)) @test_broken abs.(ForwardDiff.derivative(p -> benchmark_nlsolve_oop(quadratic_f, @SVector[1.0, 1.0], p).u[end], p)) ≈ 1 / (2 * sqrt(p)) - elseif p in broken_forwarddiff - @test_broken abs.(ForwardDiff.derivative(p -> benchmark_nlsolve_oop(quadratic_f, - @SVector[1.0, 1.0], p).u[end], p)) ≈ 1 / (2 * sqrt(p)) else @test all(res .≈ sqrt(p)) @test isapprox(abs.(ForwardDiff.derivative(p -> benchmark_nlsolve_oop(quadratic_f, @@ -490,7 +472,6 @@ end end @testset "[OOP] [Scalar AD]" begin - broken_forwarddiff = [1.6, 2.9, 3.0, 3.5, 4.0, 81.0] for p in 1.1:0.1:100.0 res = abs(benchmark_nlsolve_oop(quadratic_f, 1.0, p).u) @@ -498,9 +479,6 @@ end @test_broken res ≈ sqrt(p) @test_broken abs.(ForwardDiff.derivative(p -> benchmark_nlsolve_oop(quadratic_f, 1.0, p).u, p)) ≈ 1 / (2 * sqrt(p)) - elseif p in broken_forwarddiff - @test_broken abs.(ForwardDiff.derivative(p -> benchmark_nlsolve_oop(quadratic_f, - 1.0, p).u, p)) ≈ 1 / (2 * sqrt(p)) else @test res ≈ sqrt(p) @test isapprox(abs.(ForwardDiff.derivative(p -> benchmark_nlsolve_oop(quadratic_f, @@ -521,7 +499,7 @@ end cache = init(probN, DFSane(); 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) + reinit!(cache, iip ? [cache.u[1]] : cache.u; p = p) sol = solve!(cache) sols[i] = iip ? sol.u[1] : sol.u end @@ -566,21 +544,13 @@ end probN = NonlinearProblem{false}(quadratic_f, [1.0, 1.0], 2.0) sol = solve(probN, alg, abstol = 1e-11) - println(abs.(quadratic_f(sol.u, 2.0))) @test all(abs.(quadratic_f(sol.u, 2.0)) .< 1e-10) end end - @testset "Termination condition: $(mode) u0: $(_nameof(u0))" for mode in instances(NLSolveTerminationMode.T), + @testset "Termination condition: $(termination_condition) u0: $(_nameof(u0))" for termination_condition in TERMINATION_CONDITIONS, u0 in (1.0, [1.0, 1.0]) - if mode ∈ - (NLSolveTerminationMode.SteadyStateDefault, NLSolveTerminationMode.RelSafeBest, - NLSolveTerminationMode.AbsSafeBest) - continue - end - termination_condition = NLSolveTerminationCondition(mode; abstol = nothing, - reltol = nothing) probN = NonlinearProblem(quadratic_f, u0, 2.0) @test all(solve(probN, DFSane(); termination_condition).u .≈ sqrt(2.0)) end @@ -598,7 +568,7 @@ end end function benchmark_nlsolve_iip(f, u0, p = 2.0; linsolve, precs, - alpha_initial = 10.0) + alpha_initial = 10.0) prob = NonlinearProblem{true}(f, u0, p) return solve(prob, PseudoTransient(; linsolve, precs, alpha_initial), abstol = 1e-9) end @@ -669,13 +639,11 @@ end 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, - PseudoTransient(alpha_initial = 10.0); - maxiters = 100, + cache = init(probN, PseudoTransient(alpha_initial = 10.0); 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, alpha_new = 10.0) + reinit!(cache, iip ? [cache.u[1]] : cache.u; p = p, alpha = 10.0) sol = solve!(cache) sols[i] = iip ? sol.u[1] : sol.u end @@ -701,16 +669,9 @@ end @test all(abs.(newton_fails(sol.u, p)) .< 1e-10) end - @testset "Termination condition: $(mode) u0: $(_nameof(u0))" for mode in instances(NLSolveTerminationMode.T), + @testset "Termination condition: $(termination_condition) u0: $(_nameof(u0))" for termination_condition in TERMINATION_CONDITIONS, u0 in (1.0, [1.0, 1.0]) - if mode ∈ - (NLSolveTerminationMode.SteadyStateDefault, NLSolveTerminationMode.RelSafeBest, - NLSolveTerminationMode.AbsSafeBest) - continue - end - termination_condition = NLSolveTerminationCondition(mode; abstol = nothing, - reltol = nothing) probN = NonlinearProblem(quadratic_f, u0, 2.0) @test all(solve(probN, PseudoTransient(; alpha_initial = 10.0); termination_condition).u .≈ sqrt(2.0)) @@ -805,6 +766,13 @@ 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) + + @testset "Termination condition: $(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, GeneralBroyden(); termination_condition).u .≈ sqrt(2.0)) + end end # --- GeneralKlement tests --- @@ -831,7 +799,7 @@ end sol = benchmark_nlsolve_oop(quadratic_f, u0; linesearch) # Some are failing by a margin # @test SciMLBase.successful_retcode(sol) - @test all(abs.(sol.u .* sol.u .- 2) .< 1e-9) + @test all(abs.(sol.u .* sol.u .- 2) .< 3e-9) cache = init(NonlinearProblem{false}(quadratic_f, u0, 2.0), GeneralKlement(; linesearch), abstol = 1e-9) @@ -895,19 +863,30 @@ 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) + + @testset "Termination condition: $(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, GeneralKlement(); termination_condition).u .≈ sqrt(2.0)) + end end # --- LimitedMemoryBroyden tests --- @testset "LimitedMemoryBroyden" begin - function benchmark_nlsolve_oop(f, u0, p = 2.0; linesearch = LineSearch()) + function benchmark_nlsolve_oop(f, u0, p = 2.0; linesearch = LineSearch(), + termination_condition = AbsNormTerminationMode()) prob = NonlinearProblem{false}(f, u0, p) - return solve(prob, LimitedMemoryBroyden(; linesearch), abstol = 1e-9) + return solve(prob, LimitedMemoryBroyden(; linesearch); abstol = 1e-9, + termination_condition) end - function benchmark_nlsolve_iip(f, u0, p = 2.0; linesearch = LineSearch()) + function benchmark_nlsolve_iip(f, u0, p = 2.0; linesearch = LineSearch(), + termination_condition = AbsNormTerminationMode()) prob = NonlinearProblem{true}(f, u0, p) - return solve(prob, LimitedMemoryBroyden(; linesearch), abstol = 1e-9) + return solve(prob, LimitedMemoryBroyden(; linesearch); abstol = 1e-9, + termination_condition) end @testset "LineSearch: $(_nameof(lsmethod)) LineSearch AD: $(_nameof(ad))" for lsmethod in (Static(), @@ -985,4 +964,12 @@ end p = range(0.01, 2, length = 200) @test nlprob_iterator_interface(quadratic_f, p, Val(false))≈sqrt.(p) atol=1e-2 @test nlprob_iterator_interface(quadratic_f!, p, Val(true))≈sqrt.(p) atol=1e-2 + + @testset "Termination condition: $(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 diff --git a/test/gpu.jl b/test/gpu.jl index 43f12672a..daeee0c58 100644 --- a/test/gpu.jl +++ b/test/gpu.jl @@ -15,7 +15,7 @@ prob = NonlinearProblem(f, u0) # TrustRegion is broken # LimitedMemoryBroyden will diverge! for alg in (NewtonRaphson(), LevenbergMarquardt(; linsolve = QRFactorization()), - PseudoTransient(; alpha_initial = 10.0f0), GeneralKlement(), GeneralBroyden(), + PseudoTransient(; alpha_initial = 1.0f0), GeneralKlement(), GeneralBroyden(), LimitedMemoryBroyden()) @test_nowarn sol = solve(prob, alg; abstol = 1.0f-8, reltol = 1.0f-8) end @@ -27,7 +27,7 @@ prob = NonlinearProblem{false}(f, u0) # TrustRegion is broken # LimitedMemoryBroyden will diverge! for alg in (NewtonRaphson(), LevenbergMarquardt(; linsolve = QRFactorization()), - PseudoTransient(; alpha_initial = 10.0f0), GeneralKlement(), GeneralBroyden(), + PseudoTransient(; alpha_initial = 1.0f0), GeneralKlement(), GeneralBroyden(), LimitedMemoryBroyden()) @test_nowarn sol = solve(prob, alg; abstol = 1.0f-8, reltol = 1.0f-8) end diff --git a/test/infeasible.jl b/test/infeasible.jl new file mode 100644 index 000000000..001ce1f6e --- /dev/null +++ b/test/infeasible.jl @@ -0,0 +1,70 @@ +using LinearAlgebra, NonlinearSolve, StaticArrays, Test + +# this is infeasible +function f1!(out, u, p) + μ = 3.986004415e14 + x = 7000.0e3 + y = -6.970561549987071e-9 + z = -3.784706123246018e-9 + v_x = 8.550491684548064e-12 + u[1] + v_y = 6631.60076191005 + u[2] + v_z = 3600.665431405663 + u[3] + r = @SVector [x, y, z] + v = @SVector [v_x, v_y, v_z] + h = cross(r, v) + ev = cross(v, h) / μ - r / norm(r) + i = acos(h[3] / norm(h)) + e = norm(ev) + a = 1 / (2 / norm(r) - (norm(v)^2 / μ)) + out .= [a - 42.0e6, e - 1e-5, i - 1e-5] + return nothing +end + +# this is unfeasible +function f1(u, p) + μ = 3.986004415e14 + x = 7000.0e3 + y = -6.970561549987071e-9 + z = -3.784706123246018e-9 + v_x = 8.550491684548064e-12 + u[1] + v_y = 6631.60076191005 + u[2] + v_z = 3600.665431405663 + u[3] + r = @SVector [x, y, z] + v = @SVector [v_x, v_y, v_z] + h = cross(r, v) + ev = cross(v, h) / μ - r / norm(r) + i = acos(h[3] / norm(h)) + e = norm(ev) + a = 1 / (2 / norm(r) - (norm(v)^2 / μ)) + return [a - 42.0e6, e - 1e-5, i - 1e-5] +end + +@testset "[IIP] Infeasible" begin + u0 = [0.0, 0.0, 0.0] + prob = NonlinearProblem(f1!, u0) + sol = solve(prob) + + @test all(!isnan, sol.u) + @test !SciMLBase.successful_retcode(sol.retcode) +end + +@testset "[OOP] Infeasible" begin + u0 = [0.0, 0.0, 0.0] + prob = NonlinearProblem(f1, u0) + sol = solve(prob) + + @test all(!isnan, sol.u) + @test !SciMLBase.successful_retcode(sol.retcode) + + try + u0 = @SVector [0.0, 0.0, 0.0] + prob = NonlinearProblem(f1, u0) + sol = solve(prob) + + @test all(!isnan, sol.u) + @test !SciMLBase.successful_retcode(sol.retcode) + catch err + # Static Arrays has different default linearsolve which throws an error + @test err isa SingularException + end +end diff --git a/test/polyalgs.jl b/test/polyalgs.jl index 4f861c20b..4b5377159 100644 --- a/test/polyalgs.jl +++ b/test/polyalgs.jl @@ -34,7 +34,7 @@ function f(du, u, p) end prob = NonlinearProblem(f, [2.0, 2.0, 2.0], [1.0, 2.0, 2.5]) -sol = solve(prob) +sol = solve(prob; abstol = 1e-9) @test SciMLBase.successful_retcode(sol) # https://github.com/SciML/NonlinearSolve.jl/issues/187 @@ -43,11 +43,11 @@ ff(u, p) = 0.5 / 1.5 * NaNMath.log.(u ./ (1.0 .- u)) .- 2.0 * u .+ 1.0 uspan = (0.02, 0.1) prob = IntervalNonlinearProblem(ff, uspan) -sol = solve(prob) +sol = solve(prob; abstol = 1e-9) @test SciMLBase.successful_retcode(sol) u0 = 0.06 p = 2.0 prob = NonlinearProblem(ff, u0, p) -sol = solve(prob) +sol = solve(prob; abstol = 1e-9) @test SciMLBase.successful_retcode(sol) diff --git a/test/runtests.jl b/test/runtests.jl index d4f817d0a..a5365cf38 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -18,6 +18,7 @@ end @time @safetestset "Polyalgs" include("polyalgs.jl") @time @safetestset "Matrix Resizing" include("matrix_resizing.jl") @time @safetestset "Nonlinear Least Squares" include("nonlinear_least_squares.jl") + @time @safetestset "Infeasible Problems" include("infeasible.jl") end if GROUP == "All" || GROUP == "23TestProblems" diff --git a/test/sparse.jl b/test/sparse.jl index 256ca7530..729629c38 100644 --- a/test/sparse.jl +++ b/test/sparse.jl @@ -41,6 +41,13 @@ end u0 = init_brusselator_2d(xyd_brusselator) prob_brusselator_2d = NonlinearProblem(brusselator_2d_loop, u0, p) sol = solve(prob_brusselator_2d, NewtonRaphson()) +@test norm(sol.resid) < 1e-8 + +sol = solve(prob_brusselator_2d, NewtonRaphson(autodiff = AutoSparseForwardDiff())) +@test norm(sol.resid) < 1e-8 + +sol = solve(prob_brusselator_2d, NewtonRaphson(autodiff = AutoSparseFiniteDiff())) +@test norm(sol.resid) < 1e-8 du0 = copy(u0) jac_sparsity = Symbolics.jacobian_sparsity((du, u) -> brusselator_2d_loop(du, u, p), du0, @@ -52,13 +59,13 @@ fill!(jac_prototype, 0) ff = NonlinearFunction(brusselator_2d_loop; jac_prototype) prob_brusselator_2d = NonlinearProblem(ff, u0, p) -# for autodiff in [false, ] sol = solve(prob_brusselator_2d, NewtonRaphson()) @test norm(sol.resid) < 1e-8 @test !all(iszero, jac_prototype) sol = solve(prob_brusselator_2d, NewtonRaphson(autodiff = AutoSparseFiniteDiff())) -@test norm(sol.resid) < 1e-6 +@test norm(sol.resid) < 1e-8 cache = init(prob_brusselator_2d, NewtonRaphson(; autodiff = AutoSparseForwardDiff())); @test maximum(cache.jac_cache.coloring.colorvec) == 12 +@test cache.alg.ad isa AutoSparseForwardDiff