From c5c6d4bffb45edf0e7b6d8f056f387d4bec689e1 Mon Sep 17 00:00:00 2001 From: ErikQQY <2283984853@qq.com> Date: Mon, 26 Jun 2023 17:18:21 +0800 Subject: [PATCH 1/2] Fix TrustRegion docstring Signed-off-by: ErikQQY <2283984853@qq.com> --- docs/make.jl | 36 ++++----- docs/pages.jl | 26 +++--- docs/src/basics/FAQ.md | 4 +- src/NonlinearSolve.jl | 36 ++++----- src/ad.jl | 22 +++--- src/jacobian.jl | 48 +++++------ src/levenberg.jl | 154 ++++++++++++++++++------------------ src/raphson.jl | 66 ++++++++-------- src/trustRegion.jl | 167 ++++++++++++++++++++------------------- src/utils.jl | 78 +++++++++--------- test/23_test_cases.jl | 56 ++++++------- test/basictests.jl | 62 +++++++-------- test/convergencetests.jl | 4 +- test/sparse.jl | 4 +- 14 files changed, 382 insertions(+), 381 deletions(-) diff --git a/docs/make.jl b/docs/make.jl index d2815c502..340624ecb 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,5 +1,5 @@ using Documenter, NonlinearSolve, SimpleNonlinearSolve, Sundials, SciMLNLSolve, - NonlinearSolveMINPACK, SteadyStateDiffEq + NonlinearSolveMINPACK, SteadyStateDiffEq cp("./docs/Manifest.toml", "./docs/src/assets/Manifest.toml", force = true) cp("./docs/Project.toml", "./docs/src/assets/Project.toml", force = true) @@ -7,22 +7,22 @@ cp("./docs/Project.toml", "./docs/src/assets/Project.toml", force = true) include("pages.jl") makedocs(sitename = "NonlinearSolve.jl", - authors = "Chris Rackauckas", - modules = [NonlinearSolve, NonlinearSolve.SciMLBase, NonlinearSolve.DiffEqBase, - SimpleNonlinearSolve, Sundials, SciMLNLSolve, NonlinearSolveMINPACK, - SteadyStateDiffEq], - clean = true, doctest = false, - strict = [ - :doctest, - :linkcheck, - :parse_error, - :example_block, - # Other available options are - # :autodocs_block, :cross_references, :docs_block, :eval_block, :example_block, :footnote, :meta_block, :missing_docs, :setup_block - ], - format = Documenter.HTML(assets = ["assets/favicon.ico"], - canonical = "https://docs.sciml.ai/NonlinearSolve/stable/"), - pages = pages) + authors = "Chris Rackauckas", + modules = [NonlinearSolve, NonlinearSolve.SciMLBase, NonlinearSolve.DiffEqBase, + SimpleNonlinearSolve, Sundials, SciMLNLSolve, NonlinearSolveMINPACK, + SteadyStateDiffEq], + clean = true, doctest = false, + strict = [ + :doctest, + :linkcheck, + :parse_error, + :example_block, + # Other available options are + # :autodocs_block, :cross_references, :docs_block, :eval_block, :example_block, :footnote, :meta_block, :missing_docs, :setup_block + ], + format = Documenter.HTML(assets = ["assets/favicon.ico"], + canonical = "https://docs.sciml.ai/NonlinearSolve/stable/"), + pages = pages) deploydocs(repo = "github.com/SciML/NonlinearSolve.jl.git"; - push_preview = true) + push_preview = true) diff --git a/docs/pages.jl b/docs/pages.jl index be8d610bc..6d9b5d0e2 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -2,20 +2,20 @@ pages = ["index.md", "Tutorials" => Any["tutorials/nonlinear.md", - "tutorials/iterator_interface.md"], + "tutorials/iterator_interface.md"], "Basics" => Any["basics/NonlinearProblem.md", - "basics/NonlinearFunctions.md", - "basics/solve.md", - "basics/NonlinearSolution.md", - "basics/TerminationCondition.md", - "basics/FAQ.md"], + "basics/NonlinearFunctions.md", + "basics/solve.md", + "basics/NonlinearSolution.md", + "basics/TerminationCondition.md", + "basics/FAQ.md"], "Solver Summaries and Recommendations" => Any["solvers/NonlinearSystemSolvers.md", - "solvers/BracketingSolvers.md", - "solvers/SteadyStateSolvers.md"], + "solvers/BracketingSolvers.md", + "solvers/SteadyStateSolvers.md"], "Detailed Solver APIs" => Any["api/nonlinearsolve.md", - "api/simplenonlinearsolve.md", - "api/minpack.md", - "api/nlsolve.md", - "api/sundials.md", - "api/steadystatediffeq.md"], + "api/simplenonlinearsolve.md", + "api/minpack.md", + "api/nlsolve.md", + "api/sundials.md", + "api/steadystatediffeq.md"], ] diff --git a/docs/src/basics/FAQ.md b/docs/src/basics/FAQ.md index 7e6fbd5b5..e9bb4dd4e 100644 --- a/docs/src/basics/FAQ.md +++ b/docs/src/basics/FAQ.md @@ -16,14 +16,14 @@ myfun(x, lv) = x * sin(x) - lv function f(out, levels, u0) for i in 1:N out[i] = solve(IntervalNonlinearProblem{false}(IntervalNonlinearFunction{false}(myfun), - u0, levels[i]), Falsi()).u + u0, levels[i]), Falsi()).u end end function f2(out, levels, u0) for i in 1:N out[i] = solve(NonlinearProblem{false}(NonlinearFunction{false}(myfun), - u0, levels[i]), SimpleNewtonRaphson()).u + u0, levels[i]), SimpleNewtonRaphson()).u end end diff --git a/src/NonlinearSolve.jl b/src/NonlinearSolve.jl index 6a93a63e4..448d2800f 100644 --- a/src/NonlinearSolve.jl +++ b/src/NonlinearSolve.jl @@ -27,8 +27,8 @@ abstract type AbstractNewtonAlgorithm{CS, AD, FDT, ST, CJ} <: AbstractNonlinearSolveAlgorithm end function SciMLBase.__solve(prob::NonlinearProblem, - alg::AbstractNonlinearSolveAlgorithm, args...; - kwargs...) + alg::AbstractNonlinearSolveAlgorithm, args...; + kwargs...) cache = init(prob, alg, args...; kwargs...) sol = solve!(cache) end @@ -42,27 +42,25 @@ include("ad.jl") import PrecompileTools -PrecompileTools.@compile_workload begin - for T in (Float32, Float64) - prob = NonlinearProblem{false}((u, p) -> u .* u .- p, T(0.1), T(2)) +PrecompileTools.@compile_workload begin for T in (Float32, Float64) + prob = NonlinearProblem{false}((u, p) -> u .* u .- p, T(0.1), T(2)) - precompile_algs = if VERSION >= v"1.7" - (NewtonRaphson(), TrustRegion(), LevenbergMarquardt()) - else - (NewtonRaphson(),) - end + precompile_algs = if VERSION >= v"1.7" + (NewtonRaphson(), TrustRegion(), LevenbergMarquardt()) + else + (NewtonRaphson(),) + end - for alg in precompile_algs - solve(prob, alg, abstol = T(1e-2)) - end + for alg in precompile_algs + solve(prob, alg, abstol = T(1e-2)) + end - prob = NonlinearProblem{true}((du, u, p) -> du[1] = u[1] * u[1] - p[1], T[0.1], - T[2]) - for alg in precompile_algs - solve(prob, alg, abstol = T(1e-2)) - end + prob = NonlinearProblem{true}((du, u, p) -> du[1] = u[1] * u[1] - p[1], T[0.1], + T[2]) + for alg in precompile_algs + solve(prob, alg, abstol = T(1e-2)) end -end +end end export RadiusUpdateSchemes diff --git a/src/ad.jl b/src/ad.jl index 0dad74c56..eb3013e47 100644 --- a/src/ad.jl +++ b/src/ad.jl @@ -24,21 +24,21 @@ function scalar_nlsolve_ad(prob, alg, args...; kwargs...) end function SciMLBase.solve(prob::NonlinearProblem{<:Union{Number, StaticArraysCore.SVector}, - iip, - <:Dual{T, V, P}}, - alg::AbstractNewtonAlgorithm, - args...; kwargs...) where {iip, T, V, P} + iip, + <:Dual{T, V, P}}, + alg::AbstractNewtonAlgorithm, + args...; kwargs...) where {iip, T, V, P} sol, partials = scalar_nlsolve_ad(prob, alg, args...; kwargs...) return SciMLBase.build_solution(prob, alg, Dual{T, V, P}(sol.u, partials), sol.resid; - retcode = sol.retcode) + retcode = sol.retcode) end function SciMLBase.solve(prob::NonlinearProblem{<:Union{Number, StaticArraysCore.SVector}, - iip, - <:AbstractArray{<:Dual{T, V, P}}}, - alg::AbstractNewtonAlgorithm, - args...; - kwargs...) where {iip, T, V, P} + iip, + <:AbstractArray{<:Dual{T, V, P}}}, + alg::AbstractNewtonAlgorithm, + args...; + kwargs...) where {iip, T, V, P} sol, partials = scalar_nlsolve_ad(prob, alg, args...; kwargs...) return SciMLBase.build_solution(prob, alg, Dual{T, V, P}(sol.u, partials), sol.resid; - retcode = sol.retcode) + retcode = sol.retcode) end diff --git a/src/jacobian.jl b/src/jacobian.jl index cb696400e..eb212dfad 100644 --- a/src/jacobian.jl +++ b/src/jacobian.jl @@ -17,11 +17,11 @@ end function jacobian_finitediff_forward!(J, f, x, jac_config, forwardcache, cache) (FiniteDiff.finite_difference_jacobian!(J, f, x, jac_config, forwardcache); - maximum(jac_config.colorvec)) + maximum(jac_config.colorvec)) end function jacobian_finitediff!(J, f, x, jac_config, cache) (FiniteDiff.finite_difference_jacobian!(J, f, x, jac_config); - 2 * maximum(jac_config.colorvec)) + 2 * maximum(jac_config.colorvec)) end function jacobian!(J::AbstractMatrix{<:Number}, cache) @@ -43,7 +43,7 @@ function jacobian!(J::AbstractMatrix{<:Number}, cache) uf(fx, x) #cache.destats.nf += 1 tmp = jacobian_finitediff_forward!(J, uf, x, jac_config, fx, - cache) + cache) else # not forward difference tmp = jacobian_finitediff!(J, uf, x, jac_config, cache) end @@ -71,18 +71,18 @@ function build_jac_config(alg, f::F1, uf::F2, du1, u, tmp, du2) where {F1, F2} typeof(ForwardDiff.Tag(uf, eltype(u))) end jac_config = ForwardColorJacCache(uf, u, _chunksize; colorvec = colorvec, - sparsity = sparsity, tag = T) + sparsity = sparsity, tag = T) else if alg_difftype(alg) !== Val{:complex} jac_config = FiniteDiff.JacobianCache(tmp, du1, du2, alg_difftype(alg), - colorvec = colorvec, - sparsity = sparsity) + colorvec = colorvec, + sparsity = sparsity) else jac_config = FiniteDiff.JacobianCache(Complex{eltype(tmp)}.(tmp), - Complex{eltype(du1)}.(du1), nothing, - alg_difftype(alg), eltype(u), - colorvec = colorvec, - sparsity = sparsity) + Complex{eltype(du1)}.(du1), nothing, + alg_difftype(alg), eltype(u), + colorvec = colorvec, + sparsity = sparsity) end end else @@ -92,27 +92,27 @@ function build_jac_config(alg, f::F1, uf::F2, du1, u, tmp, du2) where {F1, F2} end function get_chunksize(jac_config::ForwardDiff.JacobianConfig{ - T, - V, - N, - D, -}) where {T, V, N, D -} + T, + V, + N, + D + }) where {T, V, N, D + } Val(N) end # don't degrade compile time information to runtime information function jacobian_finitediff(f, x, ::Type{diff_type}, dir, colorvec, sparsity, - jac_prototype) where {diff_type} + jac_prototype) where {diff_type} (FiniteDiff.finite_difference_derivative(f, x, diff_type, eltype(x), dir = dir), 2) end function jacobian_finitediff(f, x::AbstractArray, ::Type{diff_type}, dir, colorvec, - sparsity, jac_prototype) where {diff_type} + sparsity, jac_prototype) where {diff_type} f_in = diff_type === Val{:forward} ? f(x) : similar(x) ret_eltype = eltype(f_in) J = FiniteDiff.finite_difference_jacobian(f, x, diff_type, ret_eltype, f_in, - dir = dir, colorvec = colorvec, - sparsity = sparsity, - jac_prototype = jac_prototype) + dir = dir, colorvec = colorvec, + sparsity = sparsity, + jac_prototype = jac_prototype) return J, _nfcount(maximum(colorvec), diff_type) end function jacobian(cache, f::F) where {F} @@ -130,7 +130,7 @@ function jacobian(cache, f::F) where {F} sparsity, colorvec = sparsity_colorvec(cache.f, x) dir = true J, tmp = jacobian_finitediff(uf, x, alg_difftype(alg), dir, colorvec, sparsity, - jac_prototype) + jac_prototype) end J end @@ -146,6 +146,6 @@ function jacobian_autodiff(f, x::AbstractArray, nonlinfun, alg) SparseDiffTools.getsize(ForwardDiff.pickchunksize(maxcolor)))) : Int(ceil(maxcolor / _unwrap_val(chunk_size))) (forwarddiff_color_jacobian(f, x, colorvec = colorvec, sparsity = sparsity, - jac_prototype = jac_prototype, chunksize = chunk_size), - num_of_chunks) + jac_prototype = jac_prototype, chunksize = chunk_size), + num_of_chunks) end diff --git a/src/levenberg.jl b/src/levenberg.jl index 92dbe5f29..eb907dd6b 100644 --- a/src/levenberg.jl +++ b/src/levenberg.jl @@ -118,36 +118,36 @@ struct LevenbergMarquardt{CS, AD, FDT, L, P, ST, CJ, T} <: end function LevenbergMarquardt(; chunk_size = Val{0}(), - autodiff = Val{true}(), - standardtag = Val{true}(), - concrete_jac = nothing, - diff_type = Val{:forward}, - 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) + autodiff = Val{true}(), + standardtag = Val{true}(), + concrete_jac = nothing, + diff_type = Val{:forward}, + 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) LevenbergMarquardt{_unwrap_val(chunk_size), _unwrap_val(autodiff), diff_type, - typeof(linsolve), typeof(precs), _unwrap_val(standardtag), - _unwrap_val(concrete_jac), - typeof(min_damping_D)}(linsolve, precs, - damping_initial, - damping_increase_factor, - damping_decrease_factor, - finite_diff_step_geodesic, - α_geodesic, - b_uphill, - min_damping_D) + typeof(linsolve), typeof(precs), _unwrap_val(standardtag), + _unwrap_val(concrete_jac), + typeof(min_damping_D)}(linsolve, precs, + damping_initial, + damping_increase_factor, + damping_decrease_factor, + finite_diff_step_geodesic, + α_geodesic, + b_uphill, + min_damping_D) end mutable struct LevenbergMarquardtCache{iip, fType, algType, uType, duType, resType, pType, - INType, tolType, probType, ufType, L, jType, JC, - DᵀDType, λType, lossType, -} + INType, tolType, probType, ufType, L, jType, JC, + DᵀDType, λType, lossType + } f::fType alg::algType u::uType @@ -187,42 +187,42 @@ mutable struct LevenbergMarquardtCache{iip, fType, algType, uType, duType, resTy stats::NLStats function LevenbergMarquardtCache{iip}(f::fType, alg::algType, u::uType, fu::resType, - p::pType, uf::ufType, linsolve::L, J::jType, - du_tmp::duType, jac_config::JC, - force_stop::Bool, maxiters::Int, - internalnorm::INType, - retcode::SciMLBase.ReturnCode.T, abstol::tolType, - prob::probType, DᵀD::DᵀDType, JᵀJ::jType, - λ::λType, λ_factor::λType, - damping_increase_factor::λType, - damping_decrease_factor::λType, h::λType, - α_geodesic::λType, b_uphill::λType, - min_damping_D::λType, v::uType, - a::uType, tmp_vec::uType, v_old::uType, - norm_v_old::lossType, δ::uType, - loss_old::lossType, make_new_J::Bool, - fu_tmp::resType, - mat_tmp::jType, - stats::NLStats) where { - iip, fType, algType, - uType, duType, resType, - pType, INType, tolType, - probType, ufType, L, - jType, JC, DᵀDType, - λType, lossType, - } + p::pType, uf::ufType, linsolve::L, J::jType, + du_tmp::duType, jac_config::JC, + force_stop::Bool, maxiters::Int, + internalnorm::INType, + retcode::SciMLBase.ReturnCode.T, abstol::tolType, + prob::probType, DᵀD::DᵀDType, JᵀJ::jType, + λ::λType, λ_factor::λType, + damping_increase_factor::λType, + damping_decrease_factor::λType, h::λType, + α_geodesic::λType, b_uphill::λType, + min_damping_D::λType, v::uType, + a::uType, tmp_vec::uType, v_old::uType, + norm_v_old::lossType, δ::uType, + loss_old::lossType, make_new_J::Bool, + fu_tmp::resType, + mat_tmp::jType, + stats::NLStats) where { + iip, fType, algType, + uType, duType, resType, + pType, INType, tolType, + probType, ufType, L, + jType, JC, DᵀDType, + λType, lossType + } new{iip, fType, algType, uType, duType, resType, pType, INType, tolType, probType, ufType, L, jType, JC, DᵀDType, λType, lossType}(f, alg, u, fu, p, uf, linsolve, J, du_tmp, - jac_config, force_stop, maxiters, - internalnorm, retcode, abstol, 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, - norm_v_old, δ, loss_old, make_new_J, - fu_tmp, mat_tmp, stats) + jac_config, force_stop, maxiters, + internalnorm, retcode, abstol, 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, + norm_v_old, δ, loss_old, make_new_J, + fu_tmp, mat_tmp, stats) end end @@ -235,9 +235,9 @@ function jacobian_caches(alg::LevenbergMarquardt, f, u, p, ::Val{true}) recursivefill!(weight, false) Pl, Pr = wrapprecs(alg.precs(J, nothing, u, p, nothing, nothing, nothing, nothing, - nothing)..., weight) + nothing)..., weight) linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true, - Pl = Pl, Pr = Pr) + Pl = Pl, Pr = Pr) du1 = zero(u) du2 = zero(u) @@ -252,12 +252,12 @@ function jacobian_caches(alg::LevenbergMarquardt, f, u, p, ::Val{false}) end function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg::LevenbergMarquardt, - args...; - alias_u0 = false, - maxiters = 1000, - abstol = 1e-6, - internalnorm = DEFAULT_NORM, - kwargs...) where {uType, iip} + args...; + alias_u0 = false, + maxiters = 1000, + abstol = 1e-6, + internalnorm = DEFAULT_NORM, + kwargs...) where {uType, iip} if alias_u0 u = prob.u0 else @@ -302,13 +302,13 @@ function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg::LevenbergMarq mat_tmp = zero(J) return LevenbergMarquardtCache{iip}(f, alg, u, fu, p, uf, linsolve, J, du_tmp, - jac_config, false, maxiters, internalnorm, - ReturnCode.Default, abstol, 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, mat_tmp, NLStats(1, 0, 0, 0, 0)) + jac_config, false, maxiters, internalnorm, + ReturnCode.Default, abstol, 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, mat_tmp, NLStats(1, 0, 0, 0, 0)) end function perform_step!(cache::LevenbergMarquardtCache{true}) @unpack fu, f, make_new_J = cache @@ -330,7 +330,7 @@ function perform_step!(cache::LevenbergMarquardtCache{true}) mul!(cache.fu_tmp, J', fu) @. cache.mat_tmp = JᵀJ + λ * DᵀD linres = dolinsolve(alg.precs, linsolve, A = cache.mat_tmp, b = _vec(cache.fu_tmp), - linu = _vec(cache.du_tmp), p = p, reltol = cache.abstol) + linu = _vec(cache.du_tmp), p = p, reltol = cache.abstol) cache.linsolve = linres.cache @. cache.v = -cache.du_tmp @@ -342,7 +342,7 @@ function perform_step!(cache::LevenbergMarquardtCache{true}) mul!(cache.du_tmp, J, v) @. cache.fu_tmp = (2 / h) * ((cache.fu_tmp - fu) / h - cache.du_tmp) linres = dolinsolve(alg.precs, linsolve, A = J, b = _vec(cache.fu_tmp), - linu = _vec(cache.du_tmp), p = p, reltol = cache.abstol) + linu = _vec(cache.du_tmp), p = p, reltol = cache.abstol) cache.linsolve = linres.cache @. cache.a = -cache.du_tmp cache.stats.nsolve += 2 @@ -451,5 +451,5 @@ function SciMLBase.solve!(cache::LevenbergMarquardtCache) end SciMLBase.build_solution(cache.prob, cache.alg, cache.u, cache.fu; - retcode = cache.retcode, stats = cache.stats) + retcode = cache.retcode, stats = cache.stats) end diff --git a/src/raphson.jl b/src/raphson.jl index f2b2842fc..480e759e7 100644 --- a/src/raphson.jl +++ b/src/raphson.jl @@ -55,17 +55,17 @@ struct NewtonRaphson{CS, AD, FDT, L, P, ST, CJ} <: end function NewtonRaphson(; chunk_size = Val{0}(), autodiff = Val{true}(), - standardtag = Val{true}(), concrete_jac = nothing, - diff_type = Val{:forward}, linsolve = nothing, precs = DEFAULT_PRECS) + standardtag = Val{true}(), concrete_jac = nothing, + diff_type = Val{:forward}, linsolve = nothing, precs = DEFAULT_PRECS) NewtonRaphson{_unwrap_val(chunk_size), _unwrap_val(autodiff), diff_type, - typeof(linsolve), typeof(precs), _unwrap_val(standardtag), - _unwrap_val(concrete_jac)}(linsolve, - precs) + typeof(linsolve), typeof(precs), _unwrap_val(standardtag), + _unwrap_val(concrete_jac)}(linsolve, + precs) end mutable struct NewtonRaphsonCache{iip, fType, algType, uType, duType, resType, pType, - INType, tolType, - probType, ufType, L, jType, JC} + INType, tolType, + probType, ufType, L, jType, JC} f::fType alg::algType u::uType @@ -85,20 +85,22 @@ mutable struct NewtonRaphsonCache{iip, fType, algType, uType, duType, resType, p stats::NLStats function NewtonRaphsonCache{iip}(f::fType, alg::algType, u::uType, fu::resType, - p::pType, uf::ufType, linsolve::L, J::jType, du1::duType, - jac_config::JC, force_stop::Bool, maxiters::Int, internalnorm::INType, - retcode::SciMLBase.ReturnCode.T, abstol::tolType, - prob::probType, - stats::NLStats) where { - iip, fType, algType, uType, - duType, resType, pType, INType, - tolType, - probType, ufType, L, jType, JC} + p::pType, uf::ufType, linsolve::L, J::jType, + du1::duType, + jac_config::JC, force_stop::Bool, maxiters::Int, + internalnorm::INType, + retcode::SciMLBase.ReturnCode.T, abstol::tolType, + prob::probType, + stats::NLStats) where { + iip, fType, algType, uType, + duType, resType, pType, INType, + tolType, + probType, ufType, L, jType, JC} new{iip, fType, algType, uType, duType, resType, pType, INType, tolType, probType, ufType, L, jType, JC}(f, alg, u, fu, p, - uf, linsolve, J, du1, jac_config, - force_stop, maxiters, internalnorm, - retcode, abstol, prob, stats) + uf, linsolve, J, du1, jac_config, + force_stop, maxiters, internalnorm, + retcode, abstol, prob, stats) end end @@ -111,9 +113,9 @@ function jacobian_caches(alg::NewtonRaphson, f, u, p, ::Val{true}) recursivefill!(weight, false) Pl, Pr = wrapprecs(alg.precs(J, nothing, u, p, nothing, nothing, nothing, nothing, - nothing)..., weight) + nothing)..., weight) linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true, - Pl = Pl, Pr = Pr) + Pl = Pl, Pr = Pr) du1 = zero(u) du2 = zero(u) @@ -128,12 +130,12 @@ function jacobian_caches(alg::NewtonRaphson, f, u, p, ::Val{false}) end function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg::NewtonRaphson, - args...; - alias_u0 = false, - maxiters = 1000, - abstol = 1e-6, - internalnorm = DEFAULT_NORM, - kwargs...) where {uType, iip} + args...; + alias_u0 = false, + maxiters = 1000, + abstol = 1e-6, + internalnorm = DEFAULT_NORM, + kwargs...) where {uType, iip} if alias_u0 u = prob.u0 else @@ -150,8 +152,8 @@ function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg::NewtonRaphson uf, linsolve, J, du1, jac_config = jacobian_caches(alg, f, u, p, Val(iip)) return NewtonRaphsonCache{iip}(f, alg, u, fu, p, uf, linsolve, J, du1, jac_config, - false, maxiters, internalnorm, - ReturnCode.Default, abstol, prob, NLStats(1, 0, 0, 0, 0)) + false, maxiters, internalnorm, + ReturnCode.Default, abstol, prob, NLStats(1, 0, 0, 0, 0)) end function perform_step!(cache::NewtonRaphsonCache{true}) @@ -161,7 +163,7 @@ function perform_step!(cache::NewtonRaphsonCache{true}) # u = u - J \ fu linres = dolinsolve(alg.precs, linsolve, A = J, b = _vec(fu), linu = _vec(du1), - p = p, reltol = cache.abstol) + p = p, reltol = cache.abstol) cache.linsolve = linres.cache @. u = u - du1 f(fu, u, p) @@ -204,11 +206,11 @@ function SciMLBase.solve!(cache::NewtonRaphsonCache) end SciMLBase.build_solution(cache.prob, cache.alg, cache.u, cache.fu; - retcode = cache.retcode, stats = cache.stats) + retcode = cache.retcode, stats = cache.stats) end function SciMLBase.reinit!(cache::NewtonRaphsonCache{iip}, u0 = cache.u; p = cache.p, - abstol = cache.abstol, maxiters = cache.maxiters) where {iip} + abstol = cache.abstol, maxiters = cache.maxiters) where {iip} cache.p = p if iip recursivecopy!(cache.u, u0) diff --git a/src/trustRegion.jl b/src/trustRegion.jl index bb0735560..2a332a8c0 100644 --- a/src/trustRegion.jl +++ b/src/trustRegion.jl @@ -1,3 +1,11 @@ +EnumX.@enumx RadiusUpdateSchemes begin + Simple + Hei + Yuan + Bastin + Fan +end + """ ```julia TrustRegion(; chunk_size = Val{0}(), autodiff = Val{true}(), @@ -80,14 +88,6 @@ for large-scale and numerically-difficult nonlinear systems. Currently, the linear solver and chunk size choice only applies to in-place defined `NonlinearProblem`s. That is expected to change in the future. """ -EnumX.@enumx RadiusUpdateSchemes begin - Simple - Hei - Yuan - Bastin - Fan -end - struct TrustRegion{CS, AD, FDT, L, P, ST, CJ, MTR} <: AbstractNewtonAlgorithm{CS, AD, FDT, ST, CJ} linsolve::L @@ -104,34 +104,34 @@ struct TrustRegion{CS, AD, FDT, L, P, ST, CJ, MTR} <: end function TrustRegion(; chunk_size = Val{0}(), - autodiff = Val{true}(), - standardtag = Val{true}(), concrete_jac = nothing, - diff_type = Val{:forward}, 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 // 10, - 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) + autodiff = Val{true}(), + standardtag = Val{true}(), concrete_jac = nothing, + diff_type = Val{:forward}, 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 // 10, + 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) TrustRegion{_unwrap_val(chunk_size), _unwrap_val(autodiff), diff_type, - typeof(linsolve), typeof(precs), _unwrap_val(standardtag), - _unwrap_val(concrete_jac), typeof(max_trust_radius), - }(linsolve, precs, radius_update_scheme, max_trust_radius, - initial_trust_radius, - step_threshold, - shrink_threshold, - expand_threshold, - shrink_factor, - expand_factor, - max_shrink_times) + typeof(linsolve), typeof(precs), _unwrap_val(standardtag), + _unwrap_val(concrete_jac), typeof(max_trust_radius) + }(linsolve, precs, radius_update_scheme, max_trust_radius, + initial_trust_radius, + step_threshold, + shrink_threshold, + expand_threshold, + shrink_factor, + expand_factor, + max_shrink_times) end mutable struct TrustRegionCache{iip, fType, algType, uType, resType, pType, - INType, tolType, probType, ufType, L, jType, JC, floatType, - trustType, suType, su2Type, tmpType} + INType, tolType, probType, ufType, L, jType, JC, floatType, + trustType, suType, su2Type, tmpType} f::fType alg::algType u_prev::uType @@ -175,39 +175,40 @@ mutable struct TrustRegionCache{iip, fType, algType, uType, resType, pType, stats::NLStats function TrustRegionCache{iip}(f::fType, alg::algType, u_prev::uType, u::uType, - fu_prev::resType, fu::resType, p::pType, - uf::ufType, linsolve::L, J::jType, jac_config::JC, - force_stop::Bool, maxiters::Int, internalnorm::INType, - retcode::SciMLBase.ReturnCode.T, abstol::tolType, - prob::probType, - radius_update_scheme::RadiusUpdateSchemes.T, - trust_r::trustType, - max_trust_r::trustType, step_threshold::suType, - shrink_threshold::trustType, expand_threshold::trustType, - shrink_factor::trustType, expand_factor::trustType, - loss::floatType, loss_new::floatType, H::jType, - g::resType, shrink_counter::Int, step_size::su2Type, - u_tmp::tmpType, fu_new::resType, make_new_J::Bool, - r::floatType, p1::floatType, p2::floatType, - p3::floatType, p4::floatType, ϵ::floatType, - stats::NLStats) where {iip, fType, algType, uType, - resType, pType, INType, - tolType, probType, ufType, L, - jType, JC, floatType, trustType, - suType, su2Type, tmpType} + fu_prev::resType, fu::resType, p::pType, + uf::ufType, linsolve::L, J::jType, jac_config::JC, + force_stop::Bool, maxiters::Int, internalnorm::INType, + retcode::SciMLBase.ReturnCode.T, abstol::tolType, + prob::probType, + radius_update_scheme::RadiusUpdateSchemes.T, + trust_r::trustType, + max_trust_r::trustType, step_threshold::suType, + shrink_threshold::trustType, expand_threshold::trustType, + shrink_factor::trustType, expand_factor::trustType, + loss::floatType, loss_new::floatType, H::jType, + g::resType, shrink_counter::Int, step_size::su2Type, + u_tmp::tmpType, fu_new::resType, make_new_J::Bool, + r::floatType, p1::floatType, p2::floatType, + p3::floatType, p4::floatType, ϵ::floatType, + stats::NLStats) where {iip, fType, algType, uType, + resType, pType, INType, + tolType, probType, ufType, L, + jType, JC, floatType, trustType, + suType, su2Type, tmpType} new{iip, fType, algType, uType, resType, pType, INType, tolType, probType, ufType, L, jType, JC, floatType, - trustType, suType, su2Type, tmpType}(f, alg, u_prev, u, fu_prev, fu, p, uf, linsolve, J, - jac_config, force_stop, - maxiters, internalnorm, retcode, - abstol, prob, radius_update_scheme, - trust_r, max_trust_r, - step_threshold, shrink_threshold, - expand_threshold, shrink_factor, - expand_factor, loss, - loss_new, H, g, shrink_counter, - step_size, u_tmp, fu_new, - make_new_J, r, p1, p2, p3, p4, ϵ, stats) + trustType, suType, su2Type, tmpType}(f, alg, u_prev, u, fu_prev, fu, p, uf, + linsolve, J, + jac_config, force_stop, + maxiters, internalnorm, retcode, + abstol, prob, radius_update_scheme, + trust_r, max_trust_r, + step_threshold, shrink_threshold, + expand_threshold, shrink_factor, + expand_factor, loss, + loss_new, H, g, shrink_counter, + step_size, u_tmp, fu_new, + make_new_J, r, p1, p2, p3, p4, ϵ, stats) end end @@ -220,9 +221,9 @@ function jacobian_caches(alg::TrustRegion, f, u, p, ::Val{true}) recursivefill!(weight, false) Pl, Pr = wrapprecs(alg.precs(J, nothing, u, p, nothing, nothing, nothing, nothing, - nothing)..., weight) + nothing)..., weight) linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true, - Pl = Pl, Pr = Pr) + Pl = Pl, Pr = Pr) du1 = zero(u) du2 = zero(u) @@ -238,12 +239,12 @@ function jacobian_caches(alg::TrustRegion, f, u, p, ::Val{false}) end function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg::TrustRegion, - args...; - alias_u0 = false, - maxiters = 1000, - abstol = 1e-8, - internalnorm = DEFAULT_NORM, - kwargs...) where {uType, iip} + args...; + alias_u0 = false, + maxiters = 1000, + abstol = 1e-8, + internalnorm = DEFAULT_NORM, + kwargs...) where {uType, iip} if alias_u0 u = prob.u0 else @@ -342,14 +343,14 @@ function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg::TrustRegion, end return TrustRegionCache{iip}(f, alg, u_prev, u, fu_prev, fu, p, uf, linsolve, J, - jac_config, - false, maxiters, internalnorm, - ReturnCode.Default, abstol, prob, radius_update_scheme, - initial_trust_radius, - max_trust_radius, step_threshold, shrink_threshold, - expand_threshold, shrink_factor, expand_factor, loss, - loss_new, H, g, shrink_counter, step_size, u_tmp, fu_new, - make_new_J, r, p1, p2, p3, p4, ϵ, NLStats(1, 0, 0, 0, 0)) + jac_config, + false, maxiters, internalnorm, + ReturnCode.Default, abstol, prob, radius_update_scheme, + initial_trust_radius, + max_trust_radius, step_threshold, shrink_threshold, + expand_threshold, shrink_factor, expand_factor, loss, + loss_new, H, g, shrink_counter, step_size, u_tmp, fu_new, + make_new_J, r, p1, p2, p3, p4, ϵ, NLStats(1, 0, 0, 0, 0)) end function perform_step!(cache::TrustRegionCache{true}) @@ -362,8 +363,8 @@ function perform_step!(cache::TrustRegionCache{true}) end linres = dolinsolve(alg.precs, linsolve, A = cache.H, b = _vec(cache.g), - linu = _vec(u_tmp), - p = p, reltol = cache.abstol) + linu = _vec(u_tmp), + p = p, reltol = cache.abstol) cache.linsolve = linres.cache cache.u_tmp .= -1 .* u_tmp dogleg!(cache) @@ -626,11 +627,11 @@ function SciMLBase.solve!(cache::TrustRegionCache) end SciMLBase.build_solution(cache.prob, cache.alg, cache.u, cache.fu; - retcode = cache.retcode, stats = cache.stats) + retcode = cache.retcode, stats = cache.stats) end function SciMLBase.reinit!(cache::TrustRegionCache{iip}, u0 = cache.u; p = cache.p, - abstol = cache.abstol, maxiters = cache.maxiters) where {iip} + abstol = cache.abstol, maxiters = cache.maxiters) where {iip} cache.p = p if iip recursivecopy!(cache.u, u0) diff --git a/src/utils.jl b/src/utils.jl index 9eb7cc6f1..64d51a6ea 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -5,11 +5,11 @@ sqrt(real(sum(abs2, u)) / length(u)) end @inline function DEFAULT_NORM(u::StaticArraysCore.StaticArray{ - T, -}) where { - T <: Union{ - AbstractFloat, - Complex}} + T + }) where { + T <: Union{ + AbstractFloat, + Complex}} sqrt(real(sum(abs2, u)) / length(u)) end @inline function DEFAULT_NORM(u::RecursiveArrayTools.AbstractVectorOfArray) @@ -46,56 +46,56 @@ _vec(v::Number) = v _vec(v::AbstractVector) = v function alg_difftype(alg::AbstractNewtonAlgorithm{ - CS, - AD, - FDT, - ST, - CJ, -}) where {CS, AD, FDT, - ST, CJ} + CS, + AD, + FDT, + ST, + CJ + }) where {CS, AD, FDT, + ST, CJ} FDT end function concrete_jac(alg::AbstractNewtonAlgorithm{ - CS, - AD, - FDT, - ST, - CJ, -}) where {CS, AD, FDT, - ST, CJ} + CS, + AD, + FDT, + ST, + CJ + }) where {CS, AD, FDT, + ST, CJ} CJ end function get_chunksize(alg::AbstractNewtonAlgorithm{ - CS, - AD, - FDT, - ST, - CJ, -}) where {CS, AD, - FDT, - ST, CJ} + CS, + AD, + FDT, + ST, + CJ + }) where {CS, AD, + FDT, + ST, CJ} Val(CS) end function standardtag(alg::AbstractNewtonAlgorithm{ - CS, - AD, - FDT, - ST, - CJ, -}) where {CS, AD, FDT, - ST, CJ} + CS, + AD, + FDT, + ST, + CJ + }) where {CS, AD, FDT, + ST, CJ} ST 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) @@ -106,7 +106,7 @@ function dolinsolve(precs::P, linsolve; A = nothing, linu = nothing, b = nothing linsolve.Pr _Pl, _Pr = precs(linsolve.A, du, u, p, nothing, A !== nothing, Plprev, Prprev, - cachedata) + cachedata) if (_Pl !== nothing || _Pr !== nothing) _weight = weight === nothing ? (linsolve.Pr isa Diagonal ? linsolve.Pr.diag : linsolve.Pr.inner.diag) : @@ -128,7 +128,7 @@ end function wrapprecs(_Pl, _Pr, weight) if _Pl !== nothing Pl = LinearSolve.ComposePreconditioner(LinearSolve.InvPreconditioner(Diagonal(_vec(weight))), - _Pl) + _Pl) else Pl = LinearSolve.InvPreconditioner(Diagonal(_vec(weight))) end diff --git a/test/23_test_cases.jl b/test/23_test_cases.jl index 3cb0eb310..3a325c936 100644 --- a/test/23_test_cases.jl +++ b/test/23_test_cases.jl @@ -16,7 +16,7 @@ x_sol = ones(n) x_start = ones(n) x_start[1] = -1.2 p1_dict = Dict("n" => n, "start" => x_start, "sol" => x_sol, - "title" => "Generalized Rosenbrock function") + "title" => "Generalized Rosenbrock function") # ------------------------------------- Problem 2 ------------------------------------------ function p2_f!(out, x, p = nothing) @@ -31,7 +31,7 @@ n = 4 x_sol = zeros(n) x_start = [3.0, -1.0, 0.0, 1.0] p2_dict = Dict("n" => n, "start" => x_start, "sol" => x_sol, - "title" => "Powell singular function") + "title" => "Powell singular function") # ------------------------------------- Problem 3 ------------------------------------------ function p3_f!(out, x, p = nothing) @@ -44,7 +44,7 @@ n = 2 x_sol = [1.098159e-05, 9.106146] x_start = [0.0, 1.0] p3_dict = Dict("n" => n, "start" => x_start, "sol" => x_sol, - "title" => "Powell badly scaled function") + "title" => "Powell badly scaled function") # ------------------------------------- Problem 4 ------------------------------------------ function p4_f!(out, x, p = nothing) @@ -62,7 +62,7 @@ n = 4 x_sol = ones(n) x_start = [-3.0, -1.0, -3.0, -1.0] p4_dict = Dict("n" => n, "start" => x_start, "sol" => x_sol, - "title" => "Wood function") + "title" => "Wood function") # ------------------------------------- Problem 5 ------------------------------------------ function p5_f!(out, x, p = nothing) @@ -84,7 +84,7 @@ n = 3 x_sol = [1.0, 0.0, 0.0] x_start = [-1.0, 0.0, 0.0] p5_dict = Dict("n" => n, "start" => x_start, "sol" => x_sol, - "title" => "Helical valley function") + "title" => "Helical valley function") # ------------------------------------- Problem 6 ------------------------------------------ function p6_f!(out, x, p = nothing) @@ -121,7 +121,7 @@ n = 2 x_sol = [] x_start = zeros(n) p6_dict = Dict("n" => n, "start" => x_start, "sol" => x_sol, - "title" => "Watson function") + "title" => "Watson function") # ------------------------------------- Problem 7 ------------------------------------------ function p7_f!(out, x, p = nothing) @@ -156,7 +156,7 @@ for i in 1:n x_start[i] = (2 * i - n) / (n + 1) end p7_dict = Dict("n" => n, "start" => x_start, "sol" => x_sol, - "title" => "Chebyquad function") + "title" => "Chebyquad function") # ------------------------------------- Problem 8 ------------------------------------------ function p8_f!(out, x, p = nothing) @@ -170,7 +170,7 @@ n = 10 x_sol = ones(n) x_start = ones(n) ./ 2 p8_dict = Dict("n" => n, "start" => x_start, "sol" => x_sol, - "title" => "Brown almost linear function") + "title" => "Brown almost linear function") # ------------------------------------- Problem 9 ------------------------------------------ function p9_f!(out, x, p = nothing) @@ -195,7 +195,7 @@ for i in 1:n x_start[i] = (i * (i - n - 1)) / (n + 1)^2 end p9_dict = Dict("n" => n, "start" => x_start, "sol" => x_sol, - "title" => "Discrete boundary value function") + "title" => "Discrete boundary value function") # ------------------------------------- Problem 10 ----------------------------------------- function p10_f!(out, x, p = nothing) @@ -226,7 +226,7 @@ for i in 1:n x_start[i] = (i * (i - n - 1)) / (n + 1)^2 end p10_dict = Dict("n" => n, "start" => x_start, "sol" => x_sol, - "title" => "Discrete integral equation function") + "title" => "Discrete integral equation function") # ------------------------------------- Problem 11 ----------------------------------------- function p11_f!(out, x, p = nothing) @@ -242,7 +242,7 @@ n = 10 x_sol = [] x_start = ones(n) / n p11_dict = Dict("n" => n, "start" => x_start, "sol" => x_sol, - "title" => "Trigonometric function") + "title" => "Trigonometric function") # ------------------------------------- Problem 12 ----------------------------------------- function p12_f!(out, x, p = nothing) @@ -264,7 +264,7 @@ for i in 1:n x_start[i] = 1.0 - i / n end p12_dict = Dict("n" => n, "start" => x_start, "sol" => x_sol, - "title" => "Variably dimensioned function") + "title" => "Variably dimensioned function") # ------------------------------------- Problem 13 ----------------------------------------- function p13_f!(out, x, p = nothing) @@ -285,7 +285,7 @@ n = 10 x_sol = [] x_start = ones(n) .* (-1.0) p13_dict = Dict("n" => n, "start" => x_start, "sol" => x_sol, - "title" => "Broyden tridiagonal function") + "title" => "Broyden tridiagonal function") # ------------------------------------- Problem 14 ----------------------------------------- function p14_f!(out, x, p = nothing) @@ -311,7 +311,7 @@ n = 10 x_sol = [] x_start = ones(n) .* (-1.0) p14_dict = Dict("n" => n, "start" => x_start, "sol" => x_sol, - "title" => "Broyden banded function") + "title" => "Broyden banded function") # ------------------------------------- Problem 15 ----------------------------------------- function p15_f!(out, x, p = nothing) @@ -326,7 +326,7 @@ n = 4 x_sol = [0.01, 50.0, 0.0, 0.01] x_start = [1.0, 0.0, 0.0, 1.0] p15_dict = Dict("n" => n, "start" => x_start, "sol" => x_sol, - "title" => "Hammarling 2 by 2 matrix square root problem") + "title" => "Hammarling 2 by 2 matrix square root problem") # ------------------------------------- Problem 16 ----------------------------------------- function p16_f!(out, x, p = nothing) @@ -346,7 +346,7 @@ n = 9 x_sol = [0.01, 50.0, 0.0, 0.0, 0.01, 0.0, 0.0, 0.0, 0.01] x_start = [1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0] p16_dict = Dict("n" => n, "start" => x_start, "sol" => x_sol, - "title" => "Hammarling 3 by 3 matrix square root problem") + "title" => "Hammarling 3 by 3 matrix square root problem") # ------------------------------------- Problem 17 ----------------------------------------- function p17_f!(out, x, p = nothing) @@ -359,7 +359,7 @@ n = 2 x_sol = [0.0, 3.0] x_start = [1.0, 5.0] p17_dict = Dict("n" => n, "start" => x_start, "sol" => x_sol, - "title" => "Dennis and Schnabel 2 by 2 example") + "title" => "Dennis and Schnabel 2 by 2 example") # ------------------------------------- Problem 18 ----------------------------------------- function p18_f!(out, x, p = nothing) @@ -380,7 +380,7 @@ n = 2 x_sol = zeros(n) x_start = [2.0, 2.0] p18_dict = Dict("n" => n, "start" => x_start, "sol" => x_sol, - "title" => "Sample problem 18") + "title" => "Sample problem 18") # ------------------------------------- Problem 19 ----------------------------------------- function p19_f!(out, x, p = nothing) @@ -393,7 +393,7 @@ n = 2 x_sol = zeros(n) x_start = [3.0, 3.0] p19_dict = Dict("n" => n, "start" => x_start, "sol" => x_sol, - "title" => "Sample problem 19") + "title" => "Sample problem 19") # ------------------------------------- Problem 20 ----------------------------------------- function p20_f!(out, x, p = nothing) @@ -405,7 +405,7 @@ n = 1 x_sol = [5.0] # OR [0.0]... x_start = [1.0] p20_dict = Dict("n" => n, "start" => x_start, "sol" => x_sol, - "title" => "Scalar problem f(x) = x(x - 5)^2") + "title" => "Scalar problem f(x) = x(x - 5)^2") # ------------------------------------- Problem 21 ----------------------------------------- function p21_f!(out, x, p = nothing) @@ -418,7 +418,7 @@ n = 2 x_sol = [5.0, 4.0] x_start = [0.5, -2.0] p21_dict = Dict("n" => n, "start" => x_start, "sol" => x_sol, - "title" => "Freudenstein-Roth function") + "title" => "Freudenstein-Roth function") # ------------------------------------- Problem 22 ----------------------------------------- function p22_f!(out, x, p = nothing) @@ -431,7 +431,7 @@ n = 2 x_sol = [0.0, 1.0] x_start = [1.0, 0.0] p22_dict = Dict("n" => n, "start" => x_start, "sol" => x_sol, - "title" => "Boggs function") + "title" => "Boggs function") # ------------------------------------- Problem 23 ----------------------------------------- function p23_f!(out, x, p = nothing) @@ -456,15 +456,15 @@ n = 10 x_sol = [] x_start = ones(n) p23_dict = Dict("n" => n, "start" => x_start, "sol" => x_sol, - "title" => "Chandrasekhar function") + "title" => "Chandrasekhar function") # ----------------------------------- Solve problems --------------------------------------- problems = (p1_f!, p2_f!, p3_f!, p4_f!, p5_f!, p6_f!, p7_f!, p8_f!, p9_f!, p10_f!, p11_f!, - p12_f!, p13_f!, p14_f!, p15_f!, p16_f!, p17_f!, p18_f!, p19_f!, p20_f!, p21_f!, - p22_f!, p23_f!) + p12_f!, p13_f!, p14_f!, p15_f!, p16_f!, p17_f!, p18_f!, p19_f!, p20_f!, p21_f!, + p22_f!, p23_f!) dicts = (p1_dict, p2_dict, p3_dict, p4_dict, p5_dict, p6_dict, p7_dict, p8_dict, p9_dict, - p10_dict, p11_dict, p12_dict, p13_dict, p14_dict, p15_dict, p16_dict, p17_dict, - p18_dict, p19_dict, p20_dict, p21_dict, p22_dict, p23_dict) + p10_dict, p11_dict, p12_dict, p13_dict, p14_dict, p15_dict, p16_dict, p17_dict, + p18_dict, p19_dict, p20_dict, p21_dict, p22_dict, p23_dict) algs = (NewtonRaphson(), TrustRegion(), LevenbergMarquardt()) names = ("NewtonRaphson", "TrustRegion", "LevenbergMarquardt") @@ -475,7 +475,7 @@ for (problem, dict) in zip(problems, dicts) local out = similar(x) try problem(out, - solve(nlprob, alg, abstol = 1e-15, reltol = 1e-15).u, nothing) + solve(nlprob, alg, abstol = 1e-15, reltol = 1e-15).u, nothing) dict["error_" * name] = "" catch # println("Error in $name") diff --git a/test/basictests.jl b/test/basictests.jl index 2fe3d3d64..5e6392053 100644 --- a/test/basictests.jl +++ b/test/basictests.jl @@ -166,21 +166,21 @@ end function benchmark_immutable(f, u0, radius_update_scheme) probN = NonlinearProblem{false}(f, u0) solver = init(probN, TrustRegion(radius_update_scheme = radius_update_scheme), - abstol = 1e-9) + abstol = 1e-9) sol = solve!(solver) end function benchmark_mutable(f, u0, radius_update_scheme) probN = NonlinearProblem{false}(f, u0) solver = init(probN, TrustRegion(radius_update_scheme = radius_update_scheme), - abstol = 1e-9) + abstol = 1e-9) sol = solve!(solver) end function benchmark_scalar(f, u0, radius_update_scheme) probN = NonlinearProblem{false}(f, u0) sol = (solve(probN, TrustRegion(radius_update_scheme = radius_update_scheme), - abstol = 1e-9)) + abstol = 1e-9)) end function ff(u, p = nothing) @@ -210,7 +210,7 @@ end function benchmark_inplace(f, u0, radius_update_scheme) probN = NonlinearProblem{true}(f, u0) solver = init(probN, TrustRegion(radius_update_scheme = radius_update_scheme), - abstol = 1e-9) + abstol = 1e-9) sol = solve!(solver) end @@ -228,7 +228,7 @@ end for radius_update_scheme in radius_update_schemes probN = NonlinearProblem{true}(ffiip, u0) solver = init(probN, TrustRegion(radius_update_scheme = radius_update_scheme), - abstol = 1e-9) + abstol = 1e-9) @test (@ballocated solve!(solver)) < 200 end @@ -252,7 +252,7 @@ end g = function (p) probN = NonlinearProblem{false}(f, csu0, p) sol = solve(probN, TrustRegion(radius_update_scheme = RadiusUpdateSchemes.Hei), - abstol = 1e-9) + abstol = 1e-9) return sol.u[end] end @@ -264,7 +264,7 @@ end g = function (p) probN = NonlinearProblem{false}(f, csu0, p) sol = solve(probN, TrustRegion(radius_update_scheme = RadiusUpdateSchemes.Yuan), - abstol = 1e-9) + abstol = 1e-9) return sol.u[end] end @@ -276,7 +276,7 @@ end g = function (p) probN = NonlinearProblem{false}(f, csu0, p) sol = solve(probN, TrustRegion(radius_update_scheme = RadiusUpdateSchemes.Fan), - abstol = 1e-9) + abstol = 1e-9) return sol.u[end] end @@ -316,7 +316,7 @@ end g = function (p) probN = NonlinearProblem{false}(f, oftype(p, u0), p) sol = solve(probN, TrustRegion(radius_update_scheme = RadiusUpdateSchemes.Hei), - abstol = 1e-10) + abstol = 1e-10) return sol.u end @@ -330,7 +330,7 @@ end g = function (p) probN = NonlinearProblem{false}(f, oftype(p, u0), p) sol = solve(probN, TrustRegion(radius_update_scheme = RadiusUpdateSchemes.Yuan), - abstol = 1e-10) + abstol = 1e-10) return sol.u end @@ -344,7 +344,7 @@ end g = function (p) probN = NonlinearProblem{false}(f, oftype(p, u0), p) sol = solve(probN, TrustRegion(radius_update_scheme = RadiusUpdateSchemes.Fan), - abstol = 1e-10) + abstol = 1e-10) return sol.u end @@ -505,7 +505,7 @@ f(u, p) g = function (p) probN = NonlinearProblem{false}(f, u0, p) sol = solve(probN, TrustRegion(radius_update_scheme = RadiusUpdateSchemes.Fan), - abstol = 1e-10) + abstol = 1e-10) return sol.u end p = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] @@ -535,18 +535,18 @@ expand_factor = [1.5, 2.0, 3.0] max_shrink_times = [10, 20, 30] list_of_options = zip(max_trust_radius, initial_trust_radius, step_threshold, - shrink_threshold, expand_threshold, shrink_factor, - expand_factor, max_shrink_times) + shrink_threshold, expand_threshold, shrink_factor, + expand_factor, max_shrink_times) for options in list_of_options local probN, sol, alg alg = TrustRegion(max_trust_radius = options[1], - initial_trust_radius = options[2], - step_threshold = options[3], - shrink_threshold = options[4], - expand_threshold = options[5], - shrink_factor = options[6], - expand_factor = options[7], - max_shrink_times = options[8]) + initial_trust_radius = options[2], + step_threshold = options[3], + shrink_threshold = options[4], + expand_threshold = options[5], + shrink_factor = options[6], + expand_factor = options[7], + max_shrink_times = options[8]) probN = NonlinearProblem{false}(f, u0, p) sol = solve(probN, alg, abstol = 1e-10) @@ -560,12 +560,12 @@ u0 = [1.0, 1.0] function iip_oop(f, fip, u0, radius_update_scheme, maxiters) prob_iip = NonlinearProblem{true}(fip, u0) solver = init(prob_iip, TrustRegion(radius_update_scheme = radius_update_scheme), - abstol = 1e-9, maxiters = maxiters) + abstol = 1e-9, maxiters = maxiters) sol_iip = solve!(solver) prob_oop = NonlinearProblem{false}(f, u0) solver = init(prob_oop, TrustRegion(radius_update_scheme = radius_update_scheme), - abstol = 1e-9, maxiters = maxiters) + abstol = 1e-9, maxiters = maxiters) sol_oop = solve!(solver) return sol_iip.u[end], sol_oop.u[end] @@ -747,17 +747,17 @@ b_uphill = [0, 1, 2] min_damping_D = [1e-12, 1e-9, 1e-4] list_of_options = zip(damping_initial, damping_increase_factor, damping_decrease_factor, - finite_diff_step_geodesic, α_geodesic, b_uphill, - min_damping_D) + finite_diff_step_geodesic, α_geodesic, b_uphill, + min_damping_D) for options in list_of_options local probN, sol, alg alg = LevenbergMarquardt(damping_initial = options[1], - damping_increase_factor = options[2], - damping_decrease_factor = options[3], - finite_diff_step_geodesic = options[4], - α_geodesic = options[5], - b_uphill = options[6], - min_damping_D = options[7]) + damping_increase_factor = options[2], + damping_decrease_factor = options[3], + finite_diff_step_geodesic = options[4], + α_geodesic = options[5], + b_uphill = options[6], + min_damping_D = options[7]) probN = NonlinearProblem{false}(f, u0, p) sol = solve(probN, alg, abstol = 1e-10) diff --git a/test/convergencetests.jl b/test/convergencetests.jl index 751948522..c2d4e0056 100644 --- a/test/convergencetests.jl +++ b/test/convergencetests.jl @@ -29,8 +29,8 @@ tol = 1e-9 function convergence_test_oop(f, u0, p, radius_update_scheme) prob = NonlinearProblem{false}(f, oftype(p, u0), p) cache = init(prob, - TrustRegion(radius_update_scheme = radius_update_scheme), - abstol = 1e-9) + TrustRegion(radius_update_scheme = radius_update_scheme), + abstol = 1e-9) sol = solve!(cache) return cache.internalnorm(cache.u_prev - cache.u), cache.iter, sol.retcode end diff --git a/test/sparse.jl b/test/sparse.jl index 577a1c7be..6ec40509f 100644 --- a/test/sparse.jl +++ b/test/sparse.jl @@ -11,7 +11,7 @@ function brusselator_2d_loop(du, u, p) i, j = Tuple(I) x, y = xyd_brusselator[I[1]], xyd_brusselator[I[2]] ip1, im1, jp1, jm1 = limit(i + 1, N), limit(i - 1, N), limit(j + 1, N), - limit(j - 1, N) + limit(j - 1, N) du[i, j, 1] = alpha * (u[im1, j, 1] + u[ip1, j, 1] + u[i, jp1, 1] + u[i, jm1, 1] - 4u[i, j, 1]) + B + u[i, j, 1]^2 * u[i, j, 2] - (A + 1) * u[i, j, 1] + @@ -40,7 +40,7 @@ sol = solve(prob_brusselator_2d, NewtonRaphson()) du0 = copy(u0) jac_sparsity = Symbolics.jacobian_sparsity((du, u) -> brusselator_2d_loop(du, u, p), du0, - u0) + u0) ff = NonlinearFunction(brusselator_2d_loop; jac_prototype = float.(jac_sparsity)) prob_brusselator_2d = NonlinearProblem(ff, u0, p) From 8f43450f92fdeb9030219a7e7bb59d7c788ec05d Mon Sep 17 00:00:00 2001 From: ErikQQY <2283984853@qq.com> Date: Mon, 26 Jun 2023 20:41:14 +0800 Subject: [PATCH 2/2] Use new format rule Signed-off-by: ErikQQY <2283984853@qq.com> --- docs/make.jl | 36 ++++----- docs/pages.jl | 26 +++---- docs/src/basics/FAQ.md | 4 +- src/NonlinearSolve.jl | 36 ++++----- src/ad.jl | 22 +++--- src/jacobian.jl | 48 ++++++------ src/levenberg.jl | 154 +++++++++++++++++++-------------------- src/raphson.jl | 68 ++++++++--------- src/trustRegion.jl | 150 +++++++++++++++++++------------------- src/utils.jl | 78 ++++++++++---------- test/23_test_cases.jl | 56 +++++++------- test/basictests.jl | 68 ++++++++--------- test/convergencetests.jl | 4 +- test/sparse.jl | 4 +- 14 files changed, 378 insertions(+), 376 deletions(-) diff --git a/docs/make.jl b/docs/make.jl index 340624ecb..d2815c502 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,5 +1,5 @@ using Documenter, NonlinearSolve, SimpleNonlinearSolve, Sundials, SciMLNLSolve, - NonlinearSolveMINPACK, SteadyStateDiffEq + NonlinearSolveMINPACK, SteadyStateDiffEq cp("./docs/Manifest.toml", "./docs/src/assets/Manifest.toml", force = true) cp("./docs/Project.toml", "./docs/src/assets/Project.toml", force = true) @@ -7,22 +7,22 @@ cp("./docs/Project.toml", "./docs/src/assets/Project.toml", force = true) include("pages.jl") makedocs(sitename = "NonlinearSolve.jl", - authors = "Chris Rackauckas", - modules = [NonlinearSolve, NonlinearSolve.SciMLBase, NonlinearSolve.DiffEqBase, - SimpleNonlinearSolve, Sundials, SciMLNLSolve, NonlinearSolveMINPACK, - SteadyStateDiffEq], - clean = true, doctest = false, - strict = [ - :doctest, - :linkcheck, - :parse_error, - :example_block, - # Other available options are - # :autodocs_block, :cross_references, :docs_block, :eval_block, :example_block, :footnote, :meta_block, :missing_docs, :setup_block - ], - format = Documenter.HTML(assets = ["assets/favicon.ico"], - canonical = "https://docs.sciml.ai/NonlinearSolve/stable/"), - pages = pages) + authors = "Chris Rackauckas", + modules = [NonlinearSolve, NonlinearSolve.SciMLBase, NonlinearSolve.DiffEqBase, + SimpleNonlinearSolve, Sundials, SciMLNLSolve, NonlinearSolveMINPACK, + SteadyStateDiffEq], + clean = true, doctest = false, + strict = [ + :doctest, + :linkcheck, + :parse_error, + :example_block, + # Other available options are + # :autodocs_block, :cross_references, :docs_block, :eval_block, :example_block, :footnote, :meta_block, :missing_docs, :setup_block + ], + format = Documenter.HTML(assets = ["assets/favicon.ico"], + canonical = "https://docs.sciml.ai/NonlinearSolve/stable/"), + pages = pages) deploydocs(repo = "github.com/SciML/NonlinearSolve.jl.git"; - push_preview = true) + push_preview = true) diff --git a/docs/pages.jl b/docs/pages.jl index 6d9b5d0e2..be8d610bc 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -2,20 +2,20 @@ pages = ["index.md", "Tutorials" => Any["tutorials/nonlinear.md", - "tutorials/iterator_interface.md"], + "tutorials/iterator_interface.md"], "Basics" => Any["basics/NonlinearProblem.md", - "basics/NonlinearFunctions.md", - "basics/solve.md", - "basics/NonlinearSolution.md", - "basics/TerminationCondition.md", - "basics/FAQ.md"], + "basics/NonlinearFunctions.md", + "basics/solve.md", + "basics/NonlinearSolution.md", + "basics/TerminationCondition.md", + "basics/FAQ.md"], "Solver Summaries and Recommendations" => Any["solvers/NonlinearSystemSolvers.md", - "solvers/BracketingSolvers.md", - "solvers/SteadyStateSolvers.md"], + "solvers/BracketingSolvers.md", + "solvers/SteadyStateSolvers.md"], "Detailed Solver APIs" => Any["api/nonlinearsolve.md", - "api/simplenonlinearsolve.md", - "api/minpack.md", - "api/nlsolve.md", - "api/sundials.md", - "api/steadystatediffeq.md"], + "api/simplenonlinearsolve.md", + "api/minpack.md", + "api/nlsolve.md", + "api/sundials.md", + "api/steadystatediffeq.md"], ] diff --git a/docs/src/basics/FAQ.md b/docs/src/basics/FAQ.md index e9bb4dd4e..7e6fbd5b5 100644 --- a/docs/src/basics/FAQ.md +++ b/docs/src/basics/FAQ.md @@ -16,14 +16,14 @@ myfun(x, lv) = x * sin(x) - lv function f(out, levels, u0) for i in 1:N out[i] = solve(IntervalNonlinearProblem{false}(IntervalNonlinearFunction{false}(myfun), - u0, levels[i]), Falsi()).u + u0, levels[i]), Falsi()).u end end function f2(out, levels, u0) for i in 1:N out[i] = solve(NonlinearProblem{false}(NonlinearFunction{false}(myfun), - u0, levels[i]), SimpleNewtonRaphson()).u + u0, levels[i]), SimpleNewtonRaphson()).u end end diff --git a/src/NonlinearSolve.jl b/src/NonlinearSolve.jl index 448d2800f..6a93a63e4 100644 --- a/src/NonlinearSolve.jl +++ b/src/NonlinearSolve.jl @@ -27,8 +27,8 @@ abstract type AbstractNewtonAlgorithm{CS, AD, FDT, ST, CJ} <: AbstractNonlinearSolveAlgorithm end function SciMLBase.__solve(prob::NonlinearProblem, - alg::AbstractNonlinearSolveAlgorithm, args...; - kwargs...) + alg::AbstractNonlinearSolveAlgorithm, args...; + kwargs...) cache = init(prob, alg, args...; kwargs...) sol = solve!(cache) end @@ -42,25 +42,27 @@ include("ad.jl") import PrecompileTools -PrecompileTools.@compile_workload begin for T in (Float32, Float64) - prob = NonlinearProblem{false}((u, p) -> u .* u .- p, T(0.1), T(2)) +PrecompileTools.@compile_workload begin + for T in (Float32, Float64) + prob = NonlinearProblem{false}((u, p) -> u .* u .- p, T(0.1), T(2)) - precompile_algs = if VERSION >= v"1.7" - (NewtonRaphson(), TrustRegion(), LevenbergMarquardt()) - else - (NewtonRaphson(),) - end + precompile_algs = if VERSION >= v"1.7" + (NewtonRaphson(), TrustRegion(), LevenbergMarquardt()) + else + (NewtonRaphson(),) + end - for alg in precompile_algs - solve(prob, alg, abstol = T(1e-2)) - end + for alg in precompile_algs + solve(prob, alg, abstol = T(1e-2)) + end - prob = NonlinearProblem{true}((du, u, p) -> du[1] = u[1] * u[1] - p[1], T[0.1], - T[2]) - for alg in precompile_algs - solve(prob, alg, abstol = T(1e-2)) + prob = NonlinearProblem{true}((du, u, p) -> du[1] = u[1] * u[1] - p[1], T[0.1], + T[2]) + for alg in precompile_algs + solve(prob, alg, abstol = T(1e-2)) + end end -end end +end export RadiusUpdateSchemes diff --git a/src/ad.jl b/src/ad.jl index eb3013e47..0dad74c56 100644 --- a/src/ad.jl +++ b/src/ad.jl @@ -24,21 +24,21 @@ function scalar_nlsolve_ad(prob, alg, args...; kwargs...) end function SciMLBase.solve(prob::NonlinearProblem{<:Union{Number, StaticArraysCore.SVector}, - iip, - <:Dual{T, V, P}}, - alg::AbstractNewtonAlgorithm, - args...; kwargs...) where {iip, T, V, P} + iip, + <:Dual{T, V, P}}, + alg::AbstractNewtonAlgorithm, + args...; kwargs...) where {iip, T, V, P} sol, partials = scalar_nlsolve_ad(prob, alg, args...; kwargs...) return SciMLBase.build_solution(prob, alg, Dual{T, V, P}(sol.u, partials), sol.resid; - retcode = sol.retcode) + retcode = sol.retcode) end function SciMLBase.solve(prob::NonlinearProblem{<:Union{Number, StaticArraysCore.SVector}, - iip, - <:AbstractArray{<:Dual{T, V, P}}}, - alg::AbstractNewtonAlgorithm, - args...; - kwargs...) where {iip, T, V, P} + iip, + <:AbstractArray{<:Dual{T, V, P}}}, + alg::AbstractNewtonAlgorithm, + args...; + kwargs...) where {iip, T, V, P} sol, partials = scalar_nlsolve_ad(prob, alg, args...; kwargs...) return SciMLBase.build_solution(prob, alg, Dual{T, V, P}(sol.u, partials), sol.resid; - retcode = sol.retcode) + retcode = sol.retcode) end diff --git a/src/jacobian.jl b/src/jacobian.jl index eb212dfad..cb696400e 100644 --- a/src/jacobian.jl +++ b/src/jacobian.jl @@ -17,11 +17,11 @@ end function jacobian_finitediff_forward!(J, f, x, jac_config, forwardcache, cache) (FiniteDiff.finite_difference_jacobian!(J, f, x, jac_config, forwardcache); - maximum(jac_config.colorvec)) + maximum(jac_config.colorvec)) end function jacobian_finitediff!(J, f, x, jac_config, cache) (FiniteDiff.finite_difference_jacobian!(J, f, x, jac_config); - 2 * maximum(jac_config.colorvec)) + 2 * maximum(jac_config.colorvec)) end function jacobian!(J::AbstractMatrix{<:Number}, cache) @@ -43,7 +43,7 @@ function jacobian!(J::AbstractMatrix{<:Number}, cache) uf(fx, x) #cache.destats.nf += 1 tmp = jacobian_finitediff_forward!(J, uf, x, jac_config, fx, - cache) + cache) else # not forward difference tmp = jacobian_finitediff!(J, uf, x, jac_config, cache) end @@ -71,18 +71,18 @@ function build_jac_config(alg, f::F1, uf::F2, du1, u, tmp, du2) where {F1, F2} typeof(ForwardDiff.Tag(uf, eltype(u))) end jac_config = ForwardColorJacCache(uf, u, _chunksize; colorvec = colorvec, - sparsity = sparsity, tag = T) + sparsity = sparsity, tag = T) else if alg_difftype(alg) !== Val{:complex} jac_config = FiniteDiff.JacobianCache(tmp, du1, du2, alg_difftype(alg), - colorvec = colorvec, - sparsity = sparsity) + colorvec = colorvec, + sparsity = sparsity) else jac_config = FiniteDiff.JacobianCache(Complex{eltype(tmp)}.(tmp), - Complex{eltype(du1)}.(du1), nothing, - alg_difftype(alg), eltype(u), - colorvec = colorvec, - sparsity = sparsity) + Complex{eltype(du1)}.(du1), nothing, + alg_difftype(alg), eltype(u), + colorvec = colorvec, + sparsity = sparsity) end end else @@ -92,27 +92,27 @@ function build_jac_config(alg, f::F1, uf::F2, du1, u, tmp, du2) where {F1, F2} end function get_chunksize(jac_config::ForwardDiff.JacobianConfig{ - T, - V, - N, - D - }) where {T, V, N, D - } + T, + V, + N, + D, +}) where {T, V, N, D +} Val(N) end # don't degrade compile time information to runtime information function jacobian_finitediff(f, x, ::Type{diff_type}, dir, colorvec, sparsity, - jac_prototype) where {diff_type} + jac_prototype) where {diff_type} (FiniteDiff.finite_difference_derivative(f, x, diff_type, eltype(x), dir = dir), 2) end function jacobian_finitediff(f, x::AbstractArray, ::Type{diff_type}, dir, colorvec, - sparsity, jac_prototype) where {diff_type} + sparsity, jac_prototype) where {diff_type} f_in = diff_type === Val{:forward} ? f(x) : similar(x) ret_eltype = eltype(f_in) J = FiniteDiff.finite_difference_jacobian(f, x, diff_type, ret_eltype, f_in, - dir = dir, colorvec = colorvec, - sparsity = sparsity, - jac_prototype = jac_prototype) + dir = dir, colorvec = colorvec, + sparsity = sparsity, + jac_prototype = jac_prototype) return J, _nfcount(maximum(colorvec), diff_type) end function jacobian(cache, f::F) where {F} @@ -130,7 +130,7 @@ function jacobian(cache, f::F) where {F} sparsity, colorvec = sparsity_colorvec(cache.f, x) dir = true J, tmp = jacobian_finitediff(uf, x, alg_difftype(alg), dir, colorvec, sparsity, - jac_prototype) + jac_prototype) end J end @@ -146,6 +146,6 @@ function jacobian_autodiff(f, x::AbstractArray, nonlinfun, alg) SparseDiffTools.getsize(ForwardDiff.pickchunksize(maxcolor)))) : Int(ceil(maxcolor / _unwrap_val(chunk_size))) (forwarddiff_color_jacobian(f, x, colorvec = colorvec, sparsity = sparsity, - jac_prototype = jac_prototype, chunksize = chunk_size), - num_of_chunks) + jac_prototype = jac_prototype, chunksize = chunk_size), + num_of_chunks) end diff --git a/src/levenberg.jl b/src/levenberg.jl index eb907dd6b..92dbe5f29 100644 --- a/src/levenberg.jl +++ b/src/levenberg.jl @@ -118,36 +118,36 @@ struct LevenbergMarquardt{CS, AD, FDT, L, P, ST, CJ, T} <: end function LevenbergMarquardt(; chunk_size = Val{0}(), - autodiff = Val{true}(), - standardtag = Val{true}(), - concrete_jac = nothing, - diff_type = Val{:forward}, - 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) + autodiff = Val{true}(), + standardtag = Val{true}(), + concrete_jac = nothing, + diff_type = Val{:forward}, + 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) LevenbergMarquardt{_unwrap_val(chunk_size), _unwrap_val(autodiff), diff_type, - typeof(linsolve), typeof(precs), _unwrap_val(standardtag), - _unwrap_val(concrete_jac), - typeof(min_damping_D)}(linsolve, precs, - damping_initial, - damping_increase_factor, - damping_decrease_factor, - finite_diff_step_geodesic, - α_geodesic, - b_uphill, - min_damping_D) + typeof(linsolve), typeof(precs), _unwrap_val(standardtag), + _unwrap_val(concrete_jac), + typeof(min_damping_D)}(linsolve, precs, + damping_initial, + damping_increase_factor, + damping_decrease_factor, + finite_diff_step_geodesic, + α_geodesic, + b_uphill, + min_damping_D) end mutable struct LevenbergMarquardtCache{iip, fType, algType, uType, duType, resType, pType, - INType, tolType, probType, ufType, L, jType, JC, - DᵀDType, λType, lossType - } + INType, tolType, probType, ufType, L, jType, JC, + DᵀDType, λType, lossType, +} f::fType alg::algType u::uType @@ -187,42 +187,42 @@ mutable struct LevenbergMarquardtCache{iip, fType, algType, uType, duType, resTy stats::NLStats function LevenbergMarquardtCache{iip}(f::fType, alg::algType, u::uType, fu::resType, - p::pType, uf::ufType, linsolve::L, J::jType, - du_tmp::duType, jac_config::JC, - force_stop::Bool, maxiters::Int, - internalnorm::INType, - retcode::SciMLBase.ReturnCode.T, abstol::tolType, - prob::probType, DᵀD::DᵀDType, JᵀJ::jType, - λ::λType, λ_factor::λType, - damping_increase_factor::λType, - damping_decrease_factor::λType, h::λType, - α_geodesic::λType, b_uphill::λType, - min_damping_D::λType, v::uType, - a::uType, tmp_vec::uType, v_old::uType, - norm_v_old::lossType, δ::uType, - loss_old::lossType, make_new_J::Bool, - fu_tmp::resType, - mat_tmp::jType, - stats::NLStats) where { - iip, fType, algType, - uType, duType, resType, - pType, INType, tolType, - probType, ufType, L, - jType, JC, DᵀDType, - λType, lossType - } + p::pType, uf::ufType, linsolve::L, J::jType, + du_tmp::duType, jac_config::JC, + force_stop::Bool, maxiters::Int, + internalnorm::INType, + retcode::SciMLBase.ReturnCode.T, abstol::tolType, + prob::probType, DᵀD::DᵀDType, JᵀJ::jType, + λ::λType, λ_factor::λType, + damping_increase_factor::λType, + damping_decrease_factor::λType, h::λType, + α_geodesic::λType, b_uphill::λType, + min_damping_D::λType, v::uType, + a::uType, tmp_vec::uType, v_old::uType, + norm_v_old::lossType, δ::uType, + loss_old::lossType, make_new_J::Bool, + fu_tmp::resType, + mat_tmp::jType, + stats::NLStats) where { + iip, fType, algType, + uType, duType, resType, + pType, INType, tolType, + probType, ufType, L, + jType, JC, DᵀDType, + λType, lossType, + } new{iip, fType, algType, uType, duType, resType, pType, INType, tolType, probType, ufType, L, jType, JC, DᵀDType, λType, lossType}(f, alg, u, fu, p, uf, linsolve, J, du_tmp, - jac_config, force_stop, maxiters, - internalnorm, retcode, abstol, 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, - norm_v_old, δ, loss_old, make_new_J, - fu_tmp, mat_tmp, stats) + jac_config, force_stop, maxiters, + internalnorm, retcode, abstol, 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, + norm_v_old, δ, loss_old, make_new_J, + fu_tmp, mat_tmp, stats) end end @@ -235,9 +235,9 @@ function jacobian_caches(alg::LevenbergMarquardt, f, u, p, ::Val{true}) recursivefill!(weight, false) Pl, Pr = wrapprecs(alg.precs(J, nothing, u, p, nothing, nothing, nothing, nothing, - nothing)..., weight) + nothing)..., weight) linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true, - Pl = Pl, Pr = Pr) + Pl = Pl, Pr = Pr) du1 = zero(u) du2 = zero(u) @@ -252,12 +252,12 @@ function jacobian_caches(alg::LevenbergMarquardt, f, u, p, ::Val{false}) end function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg::LevenbergMarquardt, - args...; - alias_u0 = false, - maxiters = 1000, - abstol = 1e-6, - internalnorm = DEFAULT_NORM, - kwargs...) where {uType, iip} + args...; + alias_u0 = false, + maxiters = 1000, + abstol = 1e-6, + internalnorm = DEFAULT_NORM, + kwargs...) where {uType, iip} if alias_u0 u = prob.u0 else @@ -302,13 +302,13 @@ function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg::LevenbergMarq mat_tmp = zero(J) return LevenbergMarquardtCache{iip}(f, alg, u, fu, p, uf, linsolve, J, du_tmp, - jac_config, false, maxiters, internalnorm, - ReturnCode.Default, abstol, 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, mat_tmp, NLStats(1, 0, 0, 0, 0)) + jac_config, false, maxiters, internalnorm, + ReturnCode.Default, abstol, 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, mat_tmp, NLStats(1, 0, 0, 0, 0)) end function perform_step!(cache::LevenbergMarquardtCache{true}) @unpack fu, f, make_new_J = cache @@ -330,7 +330,7 @@ function perform_step!(cache::LevenbergMarquardtCache{true}) mul!(cache.fu_tmp, J', fu) @. cache.mat_tmp = JᵀJ + λ * DᵀD linres = dolinsolve(alg.precs, linsolve, A = cache.mat_tmp, b = _vec(cache.fu_tmp), - linu = _vec(cache.du_tmp), p = p, reltol = cache.abstol) + linu = _vec(cache.du_tmp), p = p, reltol = cache.abstol) cache.linsolve = linres.cache @. cache.v = -cache.du_tmp @@ -342,7 +342,7 @@ function perform_step!(cache::LevenbergMarquardtCache{true}) mul!(cache.du_tmp, J, v) @. cache.fu_tmp = (2 / h) * ((cache.fu_tmp - fu) / h - cache.du_tmp) linres = dolinsolve(alg.precs, linsolve, A = J, b = _vec(cache.fu_tmp), - linu = _vec(cache.du_tmp), p = p, reltol = cache.abstol) + linu = _vec(cache.du_tmp), p = p, reltol = cache.abstol) cache.linsolve = linres.cache @. cache.a = -cache.du_tmp cache.stats.nsolve += 2 @@ -451,5 +451,5 @@ function SciMLBase.solve!(cache::LevenbergMarquardtCache) end SciMLBase.build_solution(cache.prob, cache.alg, cache.u, cache.fu; - retcode = cache.retcode, stats = cache.stats) + retcode = cache.retcode, stats = cache.stats) end diff --git a/src/raphson.jl b/src/raphson.jl index 480e759e7..437af65c3 100644 --- a/src/raphson.jl +++ b/src/raphson.jl @@ -55,17 +55,17 @@ struct NewtonRaphson{CS, AD, FDT, L, P, ST, CJ} <: end function NewtonRaphson(; chunk_size = Val{0}(), autodiff = Val{true}(), - standardtag = Val{true}(), concrete_jac = nothing, - diff_type = Val{:forward}, linsolve = nothing, precs = DEFAULT_PRECS) + standardtag = Val{true}(), concrete_jac = nothing, + diff_type = Val{:forward}, linsolve = nothing, precs = DEFAULT_PRECS) NewtonRaphson{_unwrap_val(chunk_size), _unwrap_val(autodiff), diff_type, - typeof(linsolve), typeof(precs), _unwrap_val(standardtag), - _unwrap_val(concrete_jac)}(linsolve, - precs) + typeof(linsolve), typeof(precs), _unwrap_val(standardtag), + _unwrap_val(concrete_jac)}(linsolve, + precs) end mutable struct NewtonRaphsonCache{iip, fType, algType, uType, duType, resType, pType, - INType, tolType, - probType, ufType, L, jType, JC} + INType, tolType, + probType, ufType, L, jType, JC} f::fType alg::algType u::uType @@ -85,22 +85,22 @@ mutable struct NewtonRaphsonCache{iip, fType, algType, uType, duType, resType, p stats::NLStats function NewtonRaphsonCache{iip}(f::fType, alg::algType, u::uType, fu::resType, - p::pType, uf::ufType, linsolve::L, J::jType, - du1::duType, - jac_config::JC, force_stop::Bool, maxiters::Int, - internalnorm::INType, - retcode::SciMLBase.ReturnCode.T, abstol::tolType, - prob::probType, - stats::NLStats) where { - iip, fType, algType, uType, - duType, resType, pType, INType, - tolType, - probType, ufType, L, jType, JC} + p::pType, uf::ufType, linsolve::L, J::jType, + du1::duType, + jac_config::JC, force_stop::Bool, maxiters::Int, + internalnorm::INType, + retcode::SciMLBase.ReturnCode.T, abstol::tolType, + prob::probType, + stats::NLStats) where { + iip, fType, algType, uType, + duType, resType, pType, INType, + tolType, + probType, ufType, L, jType, JC} new{iip, fType, algType, uType, duType, resType, pType, INType, tolType, probType, ufType, L, jType, JC}(f, alg, u, fu, p, - uf, linsolve, J, du1, jac_config, - force_stop, maxiters, internalnorm, - retcode, abstol, prob, stats) + uf, linsolve, J, du1, jac_config, + force_stop, maxiters, internalnorm, + retcode, abstol, prob, stats) end end @@ -113,9 +113,9 @@ function jacobian_caches(alg::NewtonRaphson, f, u, p, ::Val{true}) recursivefill!(weight, false) Pl, Pr = wrapprecs(alg.precs(J, nothing, u, p, nothing, nothing, nothing, nothing, - nothing)..., weight) + nothing)..., weight) linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true, - Pl = Pl, Pr = Pr) + Pl = Pl, Pr = Pr) du1 = zero(u) du2 = zero(u) @@ -130,12 +130,12 @@ function jacobian_caches(alg::NewtonRaphson, f, u, p, ::Val{false}) end function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg::NewtonRaphson, - args...; - alias_u0 = false, - maxiters = 1000, - abstol = 1e-6, - internalnorm = DEFAULT_NORM, - kwargs...) where {uType, iip} + args...; + alias_u0 = false, + maxiters = 1000, + abstol = 1e-6, + internalnorm = DEFAULT_NORM, + kwargs...) where {uType, iip} if alias_u0 u = prob.u0 else @@ -152,8 +152,8 @@ function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg::NewtonRaphson uf, linsolve, J, du1, jac_config = jacobian_caches(alg, f, u, p, Val(iip)) return NewtonRaphsonCache{iip}(f, alg, u, fu, p, uf, linsolve, J, du1, jac_config, - false, maxiters, internalnorm, - ReturnCode.Default, abstol, prob, NLStats(1, 0, 0, 0, 0)) + false, maxiters, internalnorm, + ReturnCode.Default, abstol, prob, NLStats(1, 0, 0, 0, 0)) end function perform_step!(cache::NewtonRaphsonCache{true}) @@ -163,7 +163,7 @@ function perform_step!(cache::NewtonRaphsonCache{true}) # u = u - J \ fu linres = dolinsolve(alg.precs, linsolve, A = J, b = _vec(fu), linu = _vec(du1), - p = p, reltol = cache.abstol) + p = p, reltol = cache.abstol) cache.linsolve = linres.cache @. u = u - du1 f(fu, u, p) @@ -206,11 +206,11 @@ function SciMLBase.solve!(cache::NewtonRaphsonCache) end SciMLBase.build_solution(cache.prob, cache.alg, cache.u, cache.fu; - retcode = cache.retcode, stats = cache.stats) + retcode = cache.retcode, stats = cache.stats) end function SciMLBase.reinit!(cache::NewtonRaphsonCache{iip}, u0 = cache.u; p = cache.p, - abstol = cache.abstol, maxiters = cache.maxiters) where {iip} + abstol = cache.abstol, maxiters = cache.maxiters) where {iip} cache.p = p if iip recursivecopy!(cache.u, u0) diff --git a/src/trustRegion.jl b/src/trustRegion.jl index 2a332a8c0..995463a92 100644 --- a/src/trustRegion.jl +++ b/src/trustRegion.jl @@ -104,34 +104,34 @@ struct TrustRegion{CS, AD, FDT, L, P, ST, CJ, MTR} <: end function TrustRegion(; chunk_size = Val{0}(), - autodiff = Val{true}(), - standardtag = Val{true}(), concrete_jac = nothing, - diff_type = Val{:forward}, 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 // 10, - 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) + autodiff = Val{true}(), + standardtag = Val{true}(), concrete_jac = nothing, + diff_type = Val{:forward}, 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 // 10, + 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) TrustRegion{_unwrap_val(chunk_size), _unwrap_val(autodiff), diff_type, - typeof(linsolve), typeof(precs), _unwrap_val(standardtag), - _unwrap_val(concrete_jac), typeof(max_trust_radius) - }(linsolve, precs, radius_update_scheme, max_trust_radius, - initial_trust_radius, - step_threshold, - shrink_threshold, - expand_threshold, - shrink_factor, - expand_factor, - max_shrink_times) + typeof(linsolve), typeof(precs), _unwrap_val(standardtag), + _unwrap_val(concrete_jac), typeof(max_trust_radius), + }(linsolve, precs, radius_update_scheme, max_trust_radius, + initial_trust_radius, + step_threshold, + shrink_threshold, + expand_threshold, + shrink_factor, + expand_factor, + max_shrink_times) end mutable struct TrustRegionCache{iip, fType, algType, uType, resType, pType, - INType, tolType, probType, ufType, L, jType, JC, floatType, - trustType, suType, su2Type, tmpType} + INType, tolType, probType, ufType, L, jType, JC, floatType, + trustType, suType, su2Type, tmpType} f::fType alg::algType u_prev::uType @@ -175,40 +175,40 @@ mutable struct TrustRegionCache{iip, fType, algType, uType, resType, pType, stats::NLStats function TrustRegionCache{iip}(f::fType, alg::algType, u_prev::uType, u::uType, - fu_prev::resType, fu::resType, p::pType, - uf::ufType, linsolve::L, J::jType, jac_config::JC, - force_stop::Bool, maxiters::Int, internalnorm::INType, - retcode::SciMLBase.ReturnCode.T, abstol::tolType, - prob::probType, - radius_update_scheme::RadiusUpdateSchemes.T, - trust_r::trustType, - max_trust_r::trustType, step_threshold::suType, - shrink_threshold::trustType, expand_threshold::trustType, - shrink_factor::trustType, expand_factor::trustType, - loss::floatType, loss_new::floatType, H::jType, - g::resType, shrink_counter::Int, step_size::su2Type, - u_tmp::tmpType, fu_new::resType, make_new_J::Bool, - r::floatType, p1::floatType, p2::floatType, - p3::floatType, p4::floatType, ϵ::floatType, - stats::NLStats) where {iip, fType, algType, uType, - resType, pType, INType, - tolType, probType, ufType, L, - jType, JC, floatType, trustType, - suType, su2Type, tmpType} + fu_prev::resType, fu::resType, p::pType, + uf::ufType, linsolve::L, J::jType, jac_config::JC, + force_stop::Bool, maxiters::Int, internalnorm::INType, + retcode::SciMLBase.ReturnCode.T, abstol::tolType, + prob::probType, + radius_update_scheme::RadiusUpdateSchemes.T, + trust_r::trustType, + max_trust_r::trustType, step_threshold::suType, + shrink_threshold::trustType, expand_threshold::trustType, + shrink_factor::trustType, expand_factor::trustType, + loss::floatType, loss_new::floatType, H::jType, + g::resType, shrink_counter::Int, step_size::su2Type, + u_tmp::tmpType, fu_new::resType, make_new_J::Bool, + r::floatType, p1::floatType, p2::floatType, + p3::floatType, p4::floatType, ϵ::floatType, + stats::NLStats) where {iip, fType, algType, uType, + resType, pType, INType, + tolType, probType, ufType, L, + jType, JC, floatType, trustType, + suType, su2Type, tmpType} new{iip, fType, algType, uType, resType, pType, INType, tolType, probType, ufType, L, jType, JC, floatType, trustType, suType, su2Type, tmpType}(f, alg, u_prev, u, fu_prev, fu, p, uf, - linsolve, J, - jac_config, force_stop, - maxiters, internalnorm, retcode, - abstol, prob, radius_update_scheme, - trust_r, max_trust_r, - step_threshold, shrink_threshold, - expand_threshold, shrink_factor, - expand_factor, loss, - loss_new, H, g, shrink_counter, - step_size, u_tmp, fu_new, - make_new_J, r, p1, p2, p3, p4, ϵ, stats) + linsolve, J, + jac_config, force_stop, + maxiters, internalnorm, retcode, + abstol, prob, radius_update_scheme, + trust_r, max_trust_r, + step_threshold, shrink_threshold, + expand_threshold, shrink_factor, + expand_factor, loss, + loss_new, H, g, shrink_counter, + step_size, u_tmp, fu_new, + make_new_J, r, p1, p2, p3, p4, ϵ, stats) end end @@ -221,9 +221,9 @@ function jacobian_caches(alg::TrustRegion, f, u, p, ::Val{true}) recursivefill!(weight, false) Pl, Pr = wrapprecs(alg.precs(J, nothing, u, p, nothing, nothing, nothing, nothing, - nothing)..., weight) + nothing)..., weight) linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true, - Pl = Pl, Pr = Pr) + Pl = Pl, Pr = Pr) du1 = zero(u) du2 = zero(u) @@ -239,12 +239,12 @@ function jacobian_caches(alg::TrustRegion, f, u, p, ::Val{false}) end function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg::TrustRegion, - args...; - alias_u0 = false, - maxiters = 1000, - abstol = 1e-8, - internalnorm = DEFAULT_NORM, - kwargs...) where {uType, iip} + args...; + alias_u0 = false, + maxiters = 1000, + abstol = 1e-8, + internalnorm = DEFAULT_NORM, + kwargs...) where {uType, iip} if alias_u0 u = prob.u0 else @@ -343,14 +343,14 @@ function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg::TrustRegion, end return TrustRegionCache{iip}(f, alg, u_prev, u, fu_prev, fu, p, uf, linsolve, J, - jac_config, - false, maxiters, internalnorm, - ReturnCode.Default, abstol, prob, radius_update_scheme, - initial_trust_radius, - max_trust_radius, step_threshold, shrink_threshold, - expand_threshold, shrink_factor, expand_factor, loss, - loss_new, H, g, shrink_counter, step_size, u_tmp, fu_new, - make_new_J, r, p1, p2, p3, p4, ϵ, NLStats(1, 0, 0, 0, 0)) + jac_config, + false, maxiters, internalnorm, + ReturnCode.Default, abstol, prob, radius_update_scheme, + initial_trust_radius, + max_trust_radius, step_threshold, shrink_threshold, + expand_threshold, shrink_factor, expand_factor, loss, + loss_new, H, g, shrink_counter, step_size, u_tmp, fu_new, + make_new_J, r, p1, p2, p3, p4, ϵ, NLStats(1, 0, 0, 0, 0)) end function perform_step!(cache::TrustRegionCache{true}) @@ -363,8 +363,8 @@ function perform_step!(cache::TrustRegionCache{true}) end linres = dolinsolve(alg.precs, linsolve, A = cache.H, b = _vec(cache.g), - linu = _vec(u_tmp), - p = p, reltol = cache.abstol) + linu = _vec(u_tmp), + p = p, reltol = cache.abstol) cache.linsolve = linres.cache cache.u_tmp .= -1 .* u_tmp dogleg!(cache) @@ -627,11 +627,11 @@ function SciMLBase.solve!(cache::TrustRegionCache) end SciMLBase.build_solution(cache.prob, cache.alg, cache.u, cache.fu; - retcode = cache.retcode, stats = cache.stats) + retcode = cache.retcode, stats = cache.stats) end function SciMLBase.reinit!(cache::TrustRegionCache{iip}, u0 = cache.u; p = cache.p, - abstol = cache.abstol, maxiters = cache.maxiters) where {iip} + abstol = cache.abstol, maxiters = cache.maxiters) where {iip} cache.p = p if iip recursivecopy!(cache.u, u0) diff --git a/src/utils.jl b/src/utils.jl index 64d51a6ea..9eb7cc6f1 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -5,11 +5,11 @@ sqrt(real(sum(abs2, u)) / length(u)) end @inline function DEFAULT_NORM(u::StaticArraysCore.StaticArray{ - T - }) where { - T <: Union{ - AbstractFloat, - Complex}} + T, +}) where { + T <: Union{ + AbstractFloat, + Complex}} sqrt(real(sum(abs2, u)) / length(u)) end @inline function DEFAULT_NORM(u::RecursiveArrayTools.AbstractVectorOfArray) @@ -46,56 +46,56 @@ _vec(v::Number) = v _vec(v::AbstractVector) = v function alg_difftype(alg::AbstractNewtonAlgorithm{ - CS, - AD, - FDT, - ST, - CJ - }) where {CS, AD, FDT, - ST, CJ} + CS, + AD, + FDT, + ST, + CJ, +}) where {CS, AD, FDT, + ST, CJ} FDT end function concrete_jac(alg::AbstractNewtonAlgorithm{ - CS, - AD, - FDT, - ST, - CJ - }) where {CS, AD, FDT, - ST, CJ} + CS, + AD, + FDT, + ST, + CJ, +}) where {CS, AD, FDT, + ST, CJ} CJ end function get_chunksize(alg::AbstractNewtonAlgorithm{ - CS, - AD, - FDT, - ST, - CJ - }) where {CS, AD, - FDT, - ST, CJ} + CS, + AD, + FDT, + ST, + CJ, +}) where {CS, AD, + FDT, + ST, CJ} Val(CS) end function standardtag(alg::AbstractNewtonAlgorithm{ - CS, - AD, - FDT, - ST, - CJ - }) where {CS, AD, FDT, - ST, CJ} + CS, + AD, + FDT, + ST, + CJ, +}) where {CS, AD, FDT, + ST, CJ} ST 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) @@ -106,7 +106,7 @@ function dolinsolve(precs::P, linsolve; A = nothing, linu = nothing, b = nothing linsolve.Pr _Pl, _Pr = precs(linsolve.A, du, u, p, nothing, A !== nothing, Plprev, Prprev, - cachedata) + cachedata) if (_Pl !== nothing || _Pr !== nothing) _weight = weight === nothing ? (linsolve.Pr isa Diagonal ? linsolve.Pr.diag : linsolve.Pr.inner.diag) : @@ -128,7 +128,7 @@ end function wrapprecs(_Pl, _Pr, weight) if _Pl !== nothing Pl = LinearSolve.ComposePreconditioner(LinearSolve.InvPreconditioner(Diagonal(_vec(weight))), - _Pl) + _Pl) else Pl = LinearSolve.InvPreconditioner(Diagonal(_vec(weight))) end diff --git a/test/23_test_cases.jl b/test/23_test_cases.jl index 3a325c936..3cb0eb310 100644 --- a/test/23_test_cases.jl +++ b/test/23_test_cases.jl @@ -16,7 +16,7 @@ x_sol = ones(n) x_start = ones(n) x_start[1] = -1.2 p1_dict = Dict("n" => n, "start" => x_start, "sol" => x_sol, - "title" => "Generalized Rosenbrock function") + "title" => "Generalized Rosenbrock function") # ------------------------------------- Problem 2 ------------------------------------------ function p2_f!(out, x, p = nothing) @@ -31,7 +31,7 @@ n = 4 x_sol = zeros(n) x_start = [3.0, -1.0, 0.0, 1.0] p2_dict = Dict("n" => n, "start" => x_start, "sol" => x_sol, - "title" => "Powell singular function") + "title" => "Powell singular function") # ------------------------------------- Problem 3 ------------------------------------------ function p3_f!(out, x, p = nothing) @@ -44,7 +44,7 @@ n = 2 x_sol = [1.098159e-05, 9.106146] x_start = [0.0, 1.0] p3_dict = Dict("n" => n, "start" => x_start, "sol" => x_sol, - "title" => "Powell badly scaled function") + "title" => "Powell badly scaled function") # ------------------------------------- Problem 4 ------------------------------------------ function p4_f!(out, x, p = nothing) @@ -62,7 +62,7 @@ n = 4 x_sol = ones(n) x_start = [-3.0, -1.0, -3.0, -1.0] p4_dict = Dict("n" => n, "start" => x_start, "sol" => x_sol, - "title" => "Wood function") + "title" => "Wood function") # ------------------------------------- Problem 5 ------------------------------------------ function p5_f!(out, x, p = nothing) @@ -84,7 +84,7 @@ n = 3 x_sol = [1.0, 0.0, 0.0] x_start = [-1.0, 0.0, 0.0] p5_dict = Dict("n" => n, "start" => x_start, "sol" => x_sol, - "title" => "Helical valley function") + "title" => "Helical valley function") # ------------------------------------- Problem 6 ------------------------------------------ function p6_f!(out, x, p = nothing) @@ -121,7 +121,7 @@ n = 2 x_sol = [] x_start = zeros(n) p6_dict = Dict("n" => n, "start" => x_start, "sol" => x_sol, - "title" => "Watson function") + "title" => "Watson function") # ------------------------------------- Problem 7 ------------------------------------------ function p7_f!(out, x, p = nothing) @@ -156,7 +156,7 @@ for i in 1:n x_start[i] = (2 * i - n) / (n + 1) end p7_dict = Dict("n" => n, "start" => x_start, "sol" => x_sol, - "title" => "Chebyquad function") + "title" => "Chebyquad function") # ------------------------------------- Problem 8 ------------------------------------------ function p8_f!(out, x, p = nothing) @@ -170,7 +170,7 @@ n = 10 x_sol = ones(n) x_start = ones(n) ./ 2 p8_dict = Dict("n" => n, "start" => x_start, "sol" => x_sol, - "title" => "Brown almost linear function") + "title" => "Brown almost linear function") # ------------------------------------- Problem 9 ------------------------------------------ function p9_f!(out, x, p = nothing) @@ -195,7 +195,7 @@ for i in 1:n x_start[i] = (i * (i - n - 1)) / (n + 1)^2 end p9_dict = Dict("n" => n, "start" => x_start, "sol" => x_sol, - "title" => "Discrete boundary value function") + "title" => "Discrete boundary value function") # ------------------------------------- Problem 10 ----------------------------------------- function p10_f!(out, x, p = nothing) @@ -226,7 +226,7 @@ for i in 1:n x_start[i] = (i * (i - n - 1)) / (n + 1)^2 end p10_dict = Dict("n" => n, "start" => x_start, "sol" => x_sol, - "title" => "Discrete integral equation function") + "title" => "Discrete integral equation function") # ------------------------------------- Problem 11 ----------------------------------------- function p11_f!(out, x, p = nothing) @@ -242,7 +242,7 @@ n = 10 x_sol = [] x_start = ones(n) / n p11_dict = Dict("n" => n, "start" => x_start, "sol" => x_sol, - "title" => "Trigonometric function") + "title" => "Trigonometric function") # ------------------------------------- Problem 12 ----------------------------------------- function p12_f!(out, x, p = nothing) @@ -264,7 +264,7 @@ for i in 1:n x_start[i] = 1.0 - i / n end p12_dict = Dict("n" => n, "start" => x_start, "sol" => x_sol, - "title" => "Variably dimensioned function") + "title" => "Variably dimensioned function") # ------------------------------------- Problem 13 ----------------------------------------- function p13_f!(out, x, p = nothing) @@ -285,7 +285,7 @@ n = 10 x_sol = [] x_start = ones(n) .* (-1.0) p13_dict = Dict("n" => n, "start" => x_start, "sol" => x_sol, - "title" => "Broyden tridiagonal function") + "title" => "Broyden tridiagonal function") # ------------------------------------- Problem 14 ----------------------------------------- function p14_f!(out, x, p = nothing) @@ -311,7 +311,7 @@ n = 10 x_sol = [] x_start = ones(n) .* (-1.0) p14_dict = Dict("n" => n, "start" => x_start, "sol" => x_sol, - "title" => "Broyden banded function") + "title" => "Broyden banded function") # ------------------------------------- Problem 15 ----------------------------------------- function p15_f!(out, x, p = nothing) @@ -326,7 +326,7 @@ n = 4 x_sol = [0.01, 50.0, 0.0, 0.01] x_start = [1.0, 0.0, 0.0, 1.0] p15_dict = Dict("n" => n, "start" => x_start, "sol" => x_sol, - "title" => "Hammarling 2 by 2 matrix square root problem") + "title" => "Hammarling 2 by 2 matrix square root problem") # ------------------------------------- Problem 16 ----------------------------------------- function p16_f!(out, x, p = nothing) @@ -346,7 +346,7 @@ n = 9 x_sol = [0.01, 50.0, 0.0, 0.0, 0.01, 0.0, 0.0, 0.0, 0.01] x_start = [1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0] p16_dict = Dict("n" => n, "start" => x_start, "sol" => x_sol, - "title" => "Hammarling 3 by 3 matrix square root problem") + "title" => "Hammarling 3 by 3 matrix square root problem") # ------------------------------------- Problem 17 ----------------------------------------- function p17_f!(out, x, p = nothing) @@ -359,7 +359,7 @@ n = 2 x_sol = [0.0, 3.0] x_start = [1.0, 5.0] p17_dict = Dict("n" => n, "start" => x_start, "sol" => x_sol, - "title" => "Dennis and Schnabel 2 by 2 example") + "title" => "Dennis and Schnabel 2 by 2 example") # ------------------------------------- Problem 18 ----------------------------------------- function p18_f!(out, x, p = nothing) @@ -380,7 +380,7 @@ n = 2 x_sol = zeros(n) x_start = [2.0, 2.0] p18_dict = Dict("n" => n, "start" => x_start, "sol" => x_sol, - "title" => "Sample problem 18") + "title" => "Sample problem 18") # ------------------------------------- Problem 19 ----------------------------------------- function p19_f!(out, x, p = nothing) @@ -393,7 +393,7 @@ n = 2 x_sol = zeros(n) x_start = [3.0, 3.0] p19_dict = Dict("n" => n, "start" => x_start, "sol" => x_sol, - "title" => "Sample problem 19") + "title" => "Sample problem 19") # ------------------------------------- Problem 20 ----------------------------------------- function p20_f!(out, x, p = nothing) @@ -405,7 +405,7 @@ n = 1 x_sol = [5.0] # OR [0.0]... x_start = [1.0] p20_dict = Dict("n" => n, "start" => x_start, "sol" => x_sol, - "title" => "Scalar problem f(x) = x(x - 5)^2") + "title" => "Scalar problem f(x) = x(x - 5)^2") # ------------------------------------- Problem 21 ----------------------------------------- function p21_f!(out, x, p = nothing) @@ -418,7 +418,7 @@ n = 2 x_sol = [5.0, 4.0] x_start = [0.5, -2.0] p21_dict = Dict("n" => n, "start" => x_start, "sol" => x_sol, - "title" => "Freudenstein-Roth function") + "title" => "Freudenstein-Roth function") # ------------------------------------- Problem 22 ----------------------------------------- function p22_f!(out, x, p = nothing) @@ -431,7 +431,7 @@ n = 2 x_sol = [0.0, 1.0] x_start = [1.0, 0.0] p22_dict = Dict("n" => n, "start" => x_start, "sol" => x_sol, - "title" => "Boggs function") + "title" => "Boggs function") # ------------------------------------- Problem 23 ----------------------------------------- function p23_f!(out, x, p = nothing) @@ -456,15 +456,15 @@ n = 10 x_sol = [] x_start = ones(n) p23_dict = Dict("n" => n, "start" => x_start, "sol" => x_sol, - "title" => "Chandrasekhar function") + "title" => "Chandrasekhar function") # ----------------------------------- Solve problems --------------------------------------- problems = (p1_f!, p2_f!, p3_f!, p4_f!, p5_f!, p6_f!, p7_f!, p8_f!, p9_f!, p10_f!, p11_f!, - p12_f!, p13_f!, p14_f!, p15_f!, p16_f!, p17_f!, p18_f!, p19_f!, p20_f!, p21_f!, - p22_f!, p23_f!) + p12_f!, p13_f!, p14_f!, p15_f!, p16_f!, p17_f!, p18_f!, p19_f!, p20_f!, p21_f!, + p22_f!, p23_f!) dicts = (p1_dict, p2_dict, p3_dict, p4_dict, p5_dict, p6_dict, p7_dict, p8_dict, p9_dict, - p10_dict, p11_dict, p12_dict, p13_dict, p14_dict, p15_dict, p16_dict, p17_dict, - p18_dict, p19_dict, p20_dict, p21_dict, p22_dict, p23_dict) + p10_dict, p11_dict, p12_dict, p13_dict, p14_dict, p15_dict, p16_dict, p17_dict, + p18_dict, p19_dict, p20_dict, p21_dict, p22_dict, p23_dict) algs = (NewtonRaphson(), TrustRegion(), LevenbergMarquardt()) names = ("NewtonRaphson", "TrustRegion", "LevenbergMarquardt") @@ -475,7 +475,7 @@ for (problem, dict) in zip(problems, dicts) local out = similar(x) try problem(out, - solve(nlprob, alg, abstol = 1e-15, reltol = 1e-15).u, nothing) + solve(nlprob, alg, abstol = 1e-15, reltol = 1e-15).u, nothing) dict["error_" * name] = "" catch # println("Error in $name") diff --git a/test/basictests.jl b/test/basictests.jl index 5e6392053..ea792dd4f 100644 --- a/test/basictests.jl +++ b/test/basictests.jl @@ -166,21 +166,21 @@ end function benchmark_immutable(f, u0, radius_update_scheme) probN = NonlinearProblem{false}(f, u0) solver = init(probN, TrustRegion(radius_update_scheme = radius_update_scheme), - abstol = 1e-9) + abstol = 1e-9) sol = solve!(solver) end function benchmark_mutable(f, u0, radius_update_scheme) probN = NonlinearProblem{false}(f, u0) solver = init(probN, TrustRegion(radius_update_scheme = radius_update_scheme), - abstol = 1e-9) + abstol = 1e-9) sol = solve!(solver) end function benchmark_scalar(f, u0, radius_update_scheme) probN = NonlinearProblem{false}(f, u0) sol = (solve(probN, TrustRegion(radius_update_scheme = radius_update_scheme), - abstol = 1e-9)) + abstol = 1e-9)) end function ff(u, p = nothing) @@ -210,7 +210,7 @@ end function benchmark_inplace(f, u0, radius_update_scheme) probN = NonlinearProblem{true}(f, u0) solver = init(probN, TrustRegion(radius_update_scheme = radius_update_scheme), - abstol = 1e-9) + abstol = 1e-9) sol = solve!(solver) end @@ -228,7 +228,7 @@ end for radius_update_scheme in radius_update_schemes probN = NonlinearProblem{true}(ffiip, u0) solver = init(probN, TrustRegion(radius_update_scheme = radius_update_scheme), - abstol = 1e-9) + abstol = 1e-9) @test (@ballocated solve!(solver)) < 200 end @@ -252,7 +252,7 @@ end g = function (p) probN = NonlinearProblem{false}(f, csu0, p) sol = solve(probN, TrustRegion(radius_update_scheme = RadiusUpdateSchemes.Hei), - abstol = 1e-9) + abstol = 1e-9) return sol.u[end] end @@ -264,7 +264,7 @@ end g = function (p) probN = NonlinearProblem{false}(f, csu0, p) sol = solve(probN, TrustRegion(radius_update_scheme = RadiusUpdateSchemes.Yuan), - abstol = 1e-9) + abstol = 1e-9) return sol.u[end] end @@ -276,7 +276,7 @@ end g = function (p) probN = NonlinearProblem{false}(f, csu0, p) sol = solve(probN, TrustRegion(radius_update_scheme = RadiusUpdateSchemes.Fan), - abstol = 1e-9) + abstol = 1e-9) return sol.u[end] end @@ -288,7 +288,7 @@ end g = function (p) probN = NonlinearProblem{false}(f, csu0, p) sol = solve(probN, TrustRegion(radius_update_scheme = RadiusUpdateSchemes.Bastin), - abstol = 1e-9) + abstol = 1e-9) return sol.u[end] end @@ -316,7 +316,7 @@ end g = function (p) probN = NonlinearProblem{false}(f, oftype(p, u0), p) sol = solve(probN, TrustRegion(radius_update_scheme = RadiusUpdateSchemes.Hei), - abstol = 1e-10) + abstol = 1e-10) return sol.u end @@ -330,7 +330,7 @@ end g = function (p) probN = NonlinearProblem{false}(f, oftype(p, u0), p) sol = solve(probN, TrustRegion(radius_update_scheme = RadiusUpdateSchemes.Yuan), - abstol = 1e-10) + abstol = 1e-10) return sol.u end @@ -344,7 +344,7 @@ end g = function (p) probN = NonlinearProblem{false}(f, oftype(p, u0), p) sol = solve(probN, TrustRegion(radius_update_scheme = RadiusUpdateSchemes.Fan), - abstol = 1e-10) + abstol = 1e-10) return sol.u end @@ -358,7 +358,7 @@ end g = function (p) probN = NonlinearProblem{false}(f, oftype(p, u0), p) sol = solve(probN, TrustRegion(radius_update_scheme = RadiusUpdateSchemes.Bastin), - abstol = 1e-10) + abstol = 1e-10) return sol.u end @@ -505,7 +505,7 @@ f(u, p) g = function (p) probN = NonlinearProblem{false}(f, u0, p) sol = solve(probN, TrustRegion(radius_update_scheme = RadiusUpdateSchemes.Fan), - abstol = 1e-10) + abstol = 1e-10) return sol.u end p = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] @@ -516,7 +516,7 @@ f(u, p) g = function (p) probN = NonlinearProblem{false}(f, u0, p) sol = solve(probN, TrustRegion(radius_update_scheme = RadiusUpdateSchemes.Bastin), - abstol = 1e-10) + abstol = 1e-10) return sol.u end p = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] @@ -535,18 +535,18 @@ expand_factor = [1.5, 2.0, 3.0] max_shrink_times = [10, 20, 30] list_of_options = zip(max_trust_radius, initial_trust_radius, step_threshold, - shrink_threshold, expand_threshold, shrink_factor, - expand_factor, max_shrink_times) + shrink_threshold, expand_threshold, shrink_factor, + expand_factor, max_shrink_times) for options in list_of_options local probN, sol, alg alg = TrustRegion(max_trust_radius = options[1], - initial_trust_radius = options[2], - step_threshold = options[3], - shrink_threshold = options[4], - expand_threshold = options[5], - shrink_factor = options[6], - expand_factor = options[7], - max_shrink_times = options[8]) + initial_trust_radius = options[2], + step_threshold = options[3], + shrink_threshold = options[4], + expand_threshold = options[5], + shrink_factor = options[6], + expand_factor = options[7], + max_shrink_times = options[8]) probN = NonlinearProblem{false}(f, u0, p) sol = solve(probN, alg, abstol = 1e-10) @@ -560,12 +560,12 @@ u0 = [1.0, 1.0] function iip_oop(f, fip, u0, radius_update_scheme, maxiters) prob_iip = NonlinearProblem{true}(fip, u0) solver = init(prob_iip, TrustRegion(radius_update_scheme = radius_update_scheme), - abstol = 1e-9, maxiters = maxiters) + abstol = 1e-9, maxiters = maxiters) sol_iip = solve!(solver) prob_oop = NonlinearProblem{false}(f, u0) solver = init(prob_oop, TrustRegion(radius_update_scheme = radius_update_scheme), - abstol = 1e-9, maxiters = maxiters) + abstol = 1e-9, maxiters = maxiters) sol_oop = solve!(solver) return sol_iip.u[end], sol_oop.u[end] @@ -747,17 +747,17 @@ b_uphill = [0, 1, 2] min_damping_D = [1e-12, 1e-9, 1e-4] list_of_options = zip(damping_initial, damping_increase_factor, damping_decrease_factor, - finite_diff_step_geodesic, α_geodesic, b_uphill, - min_damping_D) + finite_diff_step_geodesic, α_geodesic, b_uphill, + min_damping_D) for options in list_of_options local probN, sol, alg alg = LevenbergMarquardt(damping_initial = options[1], - damping_increase_factor = options[2], - damping_decrease_factor = options[3], - finite_diff_step_geodesic = options[4], - α_geodesic = options[5], - b_uphill = options[6], - min_damping_D = options[7]) + damping_increase_factor = options[2], + damping_decrease_factor = options[3], + finite_diff_step_geodesic = options[4], + α_geodesic = options[5], + b_uphill = options[6], + min_damping_D = options[7]) probN = NonlinearProblem{false}(f, u0, p) sol = solve(probN, alg, abstol = 1e-10) diff --git a/test/convergencetests.jl b/test/convergencetests.jl index c2d4e0056..751948522 100644 --- a/test/convergencetests.jl +++ b/test/convergencetests.jl @@ -29,8 +29,8 @@ tol = 1e-9 function convergence_test_oop(f, u0, p, radius_update_scheme) prob = NonlinearProblem{false}(f, oftype(p, u0), p) cache = init(prob, - TrustRegion(radius_update_scheme = radius_update_scheme), - abstol = 1e-9) + TrustRegion(radius_update_scheme = radius_update_scheme), + abstol = 1e-9) sol = solve!(cache) return cache.internalnorm(cache.u_prev - cache.u), cache.iter, sol.retcode end diff --git a/test/sparse.jl b/test/sparse.jl index 6ec40509f..577a1c7be 100644 --- a/test/sparse.jl +++ b/test/sparse.jl @@ -11,7 +11,7 @@ function brusselator_2d_loop(du, u, p) i, j = Tuple(I) x, y = xyd_brusselator[I[1]], xyd_brusselator[I[2]] ip1, im1, jp1, jm1 = limit(i + 1, N), limit(i - 1, N), limit(j + 1, N), - limit(j - 1, N) + limit(j - 1, N) du[i, j, 1] = alpha * (u[im1, j, 1] + u[ip1, j, 1] + u[i, jp1, 1] + u[i, jm1, 1] - 4u[i, j, 1]) + B + u[i, j, 1]^2 * u[i, j, 2] - (A + 1) * u[i, j, 1] + @@ -40,7 +40,7 @@ sol = solve(prob_brusselator_2d, NewtonRaphson()) du0 = copy(u0) jac_sparsity = Symbolics.jacobian_sparsity((du, u) -> brusselator_2d_loop(du, u, p), du0, - u0) + u0) ff = NonlinearFunction(brusselator_2d_loop; jac_prototype = float.(jac_sparsity)) prob_brusselator_2d = NonlinearProblem(ff, u0, p)