From c7e10db5a2ed73791f7bec697c0d584732a0eeac Mon Sep 17 00:00:00 2001 From: Avik Pal Date: Sat, 17 Feb 2024 22:06:30 -0500 Subject: [PATCH] Enable formatting with join_lines_based_on_source=false --- .JuliaFormatter.toml | 1 + docs/make.jl | 24 ++- docs/pages.jl | 63 ++---- docs/src/basics/diagnostics_api.md | 4 +- docs/src/basics/faq.md | 11 +- docs/src/basics/sparsity_detection.md | 13 +- docs/src/index.md | 10 +- docs/src/tutorials/getting_started.md | 4 +- docs/src/tutorials/large_systems.md | 41 ++-- docs/src/tutorials/modelingtoolkit.md | 17 +- ...NonlinearSolveFastLevenbergMarquardtExt.jl | 24 +-- ...NonlinearSolveFixedPointAccelerationExt.jl | 12 +- ext/NonlinearSolveLeastSquaresOptimExt.jl | 17 +- ext/NonlinearSolveMINPACKExt.jl | 19 +- ext/NonlinearSolveNLsolveExt.jl | 18 +- ext/NonlinearSolveSIAMFANLEquationsExt.jl | 50 ++--- ext/NonlinearSolveSpeedMappingExt.jl | 10 +- src/NonlinearSolve.jl | 19 +- src/abstract_types.jl | 70 ++++--- src/algorithms/broyden.jl | 25 ++- src/algorithms/dfsane.jl | 4 +- src/algorithms/extension_algs.jl | 48 +++-- src/algorithms/gauss_newton.jl | 4 +- src/algorithms/klement.jl | 13 +- src/algorithms/lbroyden.jl | 39 ++-- src/algorithms/levenberg_marquardt.jl | 51 ++--- src/algorithms/pseudo_transient.jl | 26 +-- src/algorithms/raphson.jl | 4 +- src/algorithms/trust_region.jl | 12 +- src/core/approximate_jacobian.jl | 95 +++++---- src/core/generalized_first_order.jl | 73 ++++--- src/core/generic.jl | 19 +- src/core/spectral_methods.jl | 29 ++- src/default.jl | 80 ++++---- src/descent/damped_newton.jl | 37 ++-- src/descent/dogleg.jl | 17 +- src/descent/geodesic_acceleration.jl | 27 ++- src/descent/newton.jl | 40 ++-- src/descent/steepest.jl | 14 +- src/globalization/line_search.jl | 46 +++-- src/globalization/trust_region.jl | 86 ++++----- src/internal/approximate_initialization.jl | 50 ++--- src/internal/forward_diff.jl | 36 ++-- src/internal/helpers.jl | 51 +++-- src/internal/jacobian.jl | 32 ++-- src/internal/linear_solve.jl | 30 +-- src/internal/operators.jl | 43 ++--- src/internal/termination.jl | 20 +- src/internal/tracing.jl | 51 ++--- test/core/23_test_problems_tests.jl | 10 +- test/core/forward_ad_tests.jl | 10 +- test/core/nlls_tests.jl | 30 ++- test/core/rootfind_tests.jl | 181 +++++++++--------- test/gpu/core_tests.jl | 7 +- test/misc/bruss_tests.jl | 22 +-- test/misc/matrix_resizing_tests.jl | 12 +- test/misc/polyalg_tests.jl | 4 +- test/misc/qa_tests.jl | 4 +- test/runtests.jl | 3 +- test/wrappers/nlls_tests.jl | 31 ++- 60 files changed, 897 insertions(+), 946 deletions(-) diff --git a/.JuliaFormatter.toml b/.JuliaFormatter.toml index 1768a1a7f..66c13bae3 100644 --- a/.JuliaFormatter.toml +++ b/.JuliaFormatter.toml @@ -2,3 +2,4 @@ style = "sciml" format_markdown = true annotate_untyped_fields_with_any = false format_docstrings = true +join_lines_based_on_source = false diff --git a/docs/make.jl b/docs/make.jl index c42e05004..01b53ce8d 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,11 +1,11 @@ using Documenter, DocumenterCitations -using NonlinearSolve, - SimpleNonlinearSolve, Sundials, SteadyStateDiffEq, SciMLBase, DiffEqBase +using NonlinearSolve, SimpleNonlinearSolve, Sundials, SteadyStateDiffEq, SciMLBase, + DiffEqBase -cp(joinpath(@__DIR__, "Manifest.toml"), joinpath(@__DIR__, "src/assets/Manifest.toml"), - force = true) -cp(joinpath(@__DIR__, "Project.toml"), joinpath(@__DIR__, "src/assets/Project.toml"), - force = true) +cp(joinpath(@__DIR__, "Manifest.toml"), + joinpath(@__DIR__, "src/assets/Manifest.toml"), force = true) +cp(joinpath(@__DIR__, "Project.toml"), + joinpath(@__DIR__, "src/assets/Project.toml"), force = true) include("pages.jl") @@ -13,11 +13,15 @@ bib = CitationBibliography(joinpath(@__DIR__, "src", "refs.bib")) makedocs(; sitename = "NonlinearSolve.jl", authors = "Chris Rackauckas", - modules = [NonlinearSolve, SimpleNonlinearSolve, SteadyStateDiffEq, Sundials, - DiffEqBase, SciMLBase], - clean = true, doctest = false, linkcheck = true, + modules = [NonlinearSolve, SimpleNonlinearSolve, + SteadyStateDiffEq, Sundials, DiffEqBase, SciMLBase], + clean = true, + doctest = false, + linkcheck = true, linkcheck_ignore = ["https://twitter.com/ChrisRackauckas/status/1544743542094020615"], - checkdocs = :exports, warnonly = [:missing_docs], plugins = [bib], + checkdocs = :exports, + warnonly = [:missing_docs], + plugins = [bib], format = Documenter.HTML(assets = ["assets/favicon.ico", "assets/citations.css"], canonical = "https://docs.sciml.ai/NonlinearSolve/stable/"), pages) diff --git a/docs/pages.jl b/docs/pages.jl index 52791bc45..0706115f7 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -1,47 +1,26 @@ # Put in a separate page so it can be used by SciMLDocs.jl -pages = [ - "index.md", +pages = ["index.md", "Getting Started with Nonlinear Rootfinding in Julia" => "tutorials/getting_started.md", - "Tutorials" => Any["tutorials/code_optimization.md", - "tutorials/large_systems.md", - "tutorials/modelingtoolkit.md", - "tutorials/small_compile.md", - "tutorials/iterator_interface.md", - "tutorials/optimizing_parameterized_ode.md"], - "Basics" => Any["basics/nonlinear_problem.md", - "basics/nonlinear_functions.md", - "basics/solve.md", - "basics/nonlinear_solution.md", - "basics/autodiff.md", - "basics/termination_condition.md", - "basics/diagnostics_api.md", - "basics/sparsity_detection.md", - "basics/faq.md"], + "Tutorials" => Any["tutorials/code_optimization.md", "tutorials/large_systems.md", + "tutorials/modelingtoolkit.md", "tutorials/small_compile.md", + "tutorials/iterator_interface.md", "tutorials/optimizing_parameterized_ode.md"], + "Basics" => Any["basics/nonlinear_problem.md", "basics/nonlinear_functions.md", + "basics/solve.md", "basics/nonlinear_solution.md", "basics/autodiff.md", + "basics/termination_condition.md", "basics/diagnostics_api.md", + "basics/sparsity_detection.md", "basics/faq.md"], "Solver Summaries and Recommendations" => Any["solvers/nonlinear_system_solvers.md", - "solvers/bracketing_solvers.md", - "solvers/steady_state_solvers.md", - "solvers/nonlinear_least_squares_solvers.md", - "solvers/fixed_point_solvers.md"], - "Native Functionalities" => Any["native/solvers.md", - "native/simplenonlinearsolve.md", - "native/steadystatediffeq.md", - "native/descent.md", - "native/globalization.md", - "native/diagnostics.md"], - "Wrapped Solver APIs" => Any["api/fastlevenbergmarquardt.md", - "api/fixedpointacceleration.md", - "api/leastsquaresoptim.md", - "api/minpack.md", - "api/nlsolve.md", - "api/siamfanlequations.md", - "api/speedmapping.md", - "api/sundials.md"], - "Development Documentation" => ["devdocs/internal_interfaces.md", - "devdocs/linear_solve.md", - "devdocs/jacobian.md", - "devdocs/operators.md", - "devdocs/algorithm_helpers.md"], + "solvers/bracketing_solvers.md", "solvers/steady_state_solvers.md", + "solvers/nonlinear_least_squares_solvers.md", "solvers/fixed_point_solvers.md"], + "Native Functionalities" => Any["native/solvers.md", "native/simplenonlinearsolve.md", + "native/steadystatediffeq.md", "native/descent.md", + "native/globalization.md", "native/diagnostics.md"], + "Wrapped Solver APIs" => Any[ + "api/fastlevenbergmarquardt.md", "api/fixedpointacceleration.md", + "api/leastsquaresoptim.md", "api/minpack.md", "api/nlsolve.md", + "api/siamfanlequations.md", "api/speedmapping.md", "api/sundials.md"], + "Development Documentation" => [ + "devdocs/internal_interfaces.md", "devdocs/linear_solve.md", + "devdocs/jacobian.md", "devdocs/operators.md", "devdocs/algorithm_helpers.md"], "Release Notes" => "release_notes.md", - "References" => "references.md" -] + "References" => "references.md"] diff --git a/docs/src/basics/diagnostics_api.md b/docs/src/basics/diagnostics_api.md index 993432a00..7da29c1f7 100644 --- a/docs/src/basics/diagnostics_api.md +++ b/docs/src/basics/diagnostics_api.md @@ -33,9 +33,7 @@ using ModelingToolkit, NonlinearSolve @parameters σ ρ β # Define a nonlinear system -eqs = [0 ~ σ * (y - x), - 0 ~ x * (ρ - z) - y, - 0 ~ x * y - β * z] +eqs = [0 ~ σ * (y - x), 0 ~ x * (ρ - z) - y, 0 ~ x * y - β * z] @named ns = NonlinearSystem(eqs, [x, y, z], [σ, ρ, β]) u0 = [x => 1.0, y => 0.0, z => 0.0] diff --git a/docs/src/basics/faq.md b/docs/src/basics/faq.md index e40b57b33..265cee264 100644 --- a/docs/src/basics/faq.md +++ b/docs/src/basics/faq.md @@ -16,8 +16,8 @@ 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]), + IntervalNonlinearProblem{false}( + IntervalNonlinearFunction{false}(myfun), u0, levels[i]), Falsi()).u end end @@ -25,8 +25,7 @@ end function f2(out, levels, u0) for i in 1:N out[i] = solve( - NonlinearProblem{false}(NonlinearFunction{false}(myfun), - u0, levels[i]), + NonlinearProblem{false}(NonlinearFunction{false}(myfun), u0, levels[i]), SimpleNewtonRaphson()).u end end @@ -75,8 +74,8 @@ be a Dual number. This causes the error. To fix it: 1. Specify the `autodiff` to be `AutoFiniteDiff` ```@example dual_error_faq -sol = solve(prob_oop, LevenbergMarquardt(; autodiff = AutoFiniteDiff()); maxiters = 10000, - abstol = 1e-8) +sol = solve(prob_oop, LevenbergMarquardt(; autodiff = AutoFiniteDiff()); + maxiters = 10000, abstol = 1e-8) ``` This worked but, Finite Differencing is not the recommended approach in any scenario. diff --git a/docs/src/basics/sparsity_detection.md b/docs/src/basics/sparsity_detection.md index 222aebe19..208fc306e 100644 --- a/docs/src/basics/sparsity_detection.md +++ b/docs/src/basics/sparsity_detection.md @@ -23,12 +23,10 @@ Now you can help the solver further by providing the color vector. This can be d ```julia prob = NonlinearProblem( - NonlinearFunction(nlfunc; sparsity = jac_prototype, - colorvec = colorvec), x0) + NonlinearFunction(nlfunc; sparsity = jac_prototype, colorvec = colorvec), x0) # OR prob = NonlinearProblem( - NonlinearFunction(nlfunc; jac_prototype = jac_prototype, - colorvec = colorvec), x0) + NonlinearFunction(nlfunc; jac_prototype = jac_prototype, colorvec = colorvec), x0) ``` If the `colorvec` is not provided, then it is computed on demand. @@ -46,12 +44,11 @@ If you don't have a Sparse Jacobian Prototype, but you know the which sparsity d algorithm you want to use, then you can create your problem as follows: ```julia -prob = NonlinearProblem(NonlinearFunction(nlfunc; sparsity = SymbolicsSparsityDetection()), - x0) # Remember to have Symbolics.jl loaded +prob = NonlinearProblem( + NonlinearFunction(nlfunc; sparsity = SymbolicsSparsityDetection()), x0) # Remember to have Symbolics.jl loaded # OR prob = NonlinearProblem( - NonlinearFunction(nlfunc; sparsity = ApproximateJacobianSparsity()), - x0) + NonlinearFunction(nlfunc; sparsity = ApproximateJacobianSparsity()), x0) ``` These Detection Algorithms are from [SparseDiffTools.jl](https://github.com/JuliaDiff/SparseDiffTools.jl), diff --git a/docs/src/index.md b/docs/src/index.md index b35f09f3a..bf1f52a43 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -84,9 +84,15 @@ using TOML using Markdown version = TOML.parse(read("../../Project.toml", String))["version"] name = TOML.parse(read("../../Project.toml", String))["name"] -link_manifest = "https://github.com/SciML/" * name * ".jl/tree/gh-pages/v" * version * +link_manifest = "https://github.com/SciML/" * + name * + ".jl/tree/gh-pages/v" * + version * "/assets/Manifest.toml" -link_project = "https://github.com/SciML/" * name * ".jl/tree/gh-pages/v" * version * +link_project = "https://github.com/SciML/" * + name * + ".jl/tree/gh-pages/v" * + version * "/assets/Project.toml" Markdown.parse("""You can also download the [manifest]($link_manifest) diff --git a/docs/src/tutorials/getting_started.md b/docs/src/tutorials/getting_started.md index 26bf9faa9..d0729eeaf 100644 --- a/docs/src/tutorials/getting_started.md +++ b/docs/src/tutorials/getting_started.md @@ -180,8 +180,8 @@ be skipped for out of place problems): ```@example 1 u0 = [0.0, 0.0] -prob = NonlinearLeastSquaresProblem(NonlinearFunction(nlls!, resid_prototype = zeros(3)), - u0) +prob = NonlinearLeastSquaresProblem( + NonlinearFunction(nlls!, resid_prototype = zeros(3)), u0) solve(prob) ``` diff --git a/docs/src/tutorials/large_systems.md b/docs/src/tutorials/large_systems.md index aedd58445..19e7493d5 100644 --- a/docs/src/tutorials/large_systems.md +++ b/docs/src/tutorials/large_systems.md @@ -78,11 +78,10 @@ function brusselator_2d_loop(du, u, p) 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] + - brusselator_f(x, y) + B + + u[i, j, 1]^2 * u[i, j, 2] - (A + 1) * u[i, j, 1] + brusselator_f(x, y) du[i, j, 2] = alpha * (u[im1, j, 2] + u[ip1, j, 2] + u[i, jp1, 2] + u[i, jm1, 2] - - 4u[i, j, 2]) + - A * u[i, j, 1] - u[i, j, 1]^2 * u[i, j, 2] + 4u[i, j, 2]) + A * u[i, j, 1] - u[i, j, 1]^2 * u[i, j, 2] end end p = (3.4, 1.0, 10.0, step(xyd_brusselator)) @@ -100,8 +99,8 @@ function init_brusselator_2d(xyd) end u0 = init_brusselator_2d(xyd_brusselator) -prob_brusselator_2d = NonlinearProblem(brusselator_2d_loop, u0, p; abstol = 1e-10, - reltol = 1e-10) +prob_brusselator_2d = NonlinearProblem( + brusselator_2d_loop, u0, p; abstol = 1e-10, reltol = 1e-10) ``` ## Choosing Jacobian Types @@ -140,11 +139,11 @@ using BenchmarkTools # for @btime @btime solve(prob_brusselator_2d, NewtonRaphson(; autodiff = AutoSparseForwardDiff(; chunksize = 32))); @btime solve(prob_brusselator_2d, - NewtonRaphson(; autodiff = AutoSparseForwardDiff(; chunksize = 32), - linsolve = KLUFactorization())); + NewtonRaphson(; + autodiff = AutoSparseForwardDiff(; chunksize = 32), linsolve = KLUFactorization())); @btime solve(prob_brusselator_2d, - NewtonRaphson(; autodiff = AutoSparseForwardDiff(; chunksize = 32), - linsolve = KrylovJL_GMRES())); + NewtonRaphson(; + autodiff = AutoSparseForwardDiff(; chunksize = 32), linsolve = KrylovJL_GMRES())); nothing # hide ``` @@ -164,8 +163,8 @@ kick out a sparse matrix with our pattern, that we can turn into our `jac_protot ```@example ill_conditioned_nlprob using Symbolics du0 = copy(u0) -jac_sparsity = Symbolics.jacobian_sparsity((du, u) -> brusselator_2d_loop(du, u, p), - du0, u0) +jac_sparsity = Symbolics.jacobian_sparsity( + (du, u) -> brusselator_2d_loop(du, u, p), du0, u0) ``` Notice that Julia gives a nice print out of the sparsity pattern. That's neat, and would be @@ -275,8 +274,8 @@ function algebraicmultigrid(W, du, u, p, t, newW, Plprev, Prprev, solverdata) end @btime solve(prob_brusselator_2d_sparse, - NewtonRaphson(linsolve = KrylovJL_GMRES(), precs = algebraicmultigrid, - concrete_jac = true)); + NewtonRaphson( + linsolve = KrylovJL_GMRES(), precs = algebraicmultigrid, concrete_jac = true)); nothing # hide ``` @@ -286,8 +285,8 @@ or with a Jacobi smoother: function algebraicmultigrid2(W, du, u, p, t, newW, Plprev, Prprev, solverdata) if newW === nothing || newW A = convert(AbstractMatrix, W) - Pl = AlgebraicMultigrid.aspreconditioner(AlgebraicMultigrid.ruge_stuben(A, - presmoother = AlgebraicMultigrid.Jacobi(rand(size(A, 1))), + Pl = AlgebraicMultigrid.aspreconditioner(AlgebraicMultigrid.ruge_stuben( + A, presmoother = AlgebraicMultigrid.Jacobi(rand(size(A, 1))), postsmoother = AlgebraicMultigrid.Jacobi(rand(size(A, 1))))) else Pl = Plprev @@ -296,8 +295,8 @@ function algebraicmultigrid2(W, du, u, p, t, newW, Plprev, Prprev, solverdata) end @btime solve(prob_brusselator_2d_sparse, - NewtonRaphson(linsolve = KrylovJL_GMRES(), precs = algebraicmultigrid2, - concrete_jac = true)); + NewtonRaphson( + linsolve = KrylovJL_GMRES(), precs = algebraicmultigrid2, concrete_jac = true)); nothing # hide ``` @@ -309,12 +308,10 @@ sparsity detection. Let's compare the two by setting the sparsity detection algo ```@example ill_conditioned_nlprob prob_brusselator_2d_exact = NonlinearProblem( - NonlinearFunction(brusselator_2d_loop; - sparsity = SymbolicsSparsityDetection()), + NonlinearFunction(brusselator_2d_loop; sparsity = SymbolicsSparsityDetection()), u0, p; abstol = 1e-10, reltol = 1e-10) prob_brusselator_2d_approx = NonlinearProblem( - NonlinearFunction(brusselator_2d_loop; - sparsity = ApproximateJacobianSparsity()), + NonlinearFunction(brusselator_2d_loop; sparsity = ApproximateJacobianSparsity()), u0, p; abstol = 1e-10, reltol = 1e-10) @btime solve(prob_brusselator_2d_exact, NewtonRaphson()); diff --git a/docs/src/tutorials/modelingtoolkit.md b/docs/src/tutorials/modelingtoolkit.md index 17dbc985e..845b727b7 100644 --- a/docs/src/tutorials/modelingtoolkit.md +++ b/docs/src/tutorials/modelingtoolkit.md @@ -12,14 +12,10 @@ using ModelingToolkit, NonlinearSolve @parameters σ ρ β # Define a nonlinear system -eqs = [0 ~ σ * (y - x), - 0 ~ x * (ρ - z) - y, - 0 ~ x * y - β * z] +eqs = [0 ~ σ * (y - x), 0 ~ x * (ρ - z) - y, 0 ~ x * y - β * z] @named ns = NonlinearSystem(eqs, [x, y, z], [σ, ρ, β]) -u0 = [x => 1.0, - y => 0.0, - z => 0.0] +u0 = [x => 1.0, y => 0.0, z => 0.0] ps = [σ => 10.0 ρ => 26.0 @@ -60,13 +56,8 @@ solved more simply. Let's take a look at a quick system: ```@example mtk @variables u1 u2 u3 u4 u5 -eqs = [ - 0 ~ u1 - sin(u5), - 0 ~ u2 - cos(u1), - 0 ~ u3 - hypot(u1, u2), - 0 ~ u4 - hypot(u2, u3), - 0 ~ u5 - hypot(u4, u1) -] +eqs = [0 ~ u1 - sin(u5), 0 ~ u2 - cos(u1), 0 ~ u3 - hypot(u1, u2), + 0 ~ u4 - hypot(u2, u3), 0 ~ u5 - hypot(u4, u1)] @named sys = NonlinearSystem(eqs, [u1, u2, u3, u4, u5], []) ``` diff --git a/ext/NonlinearSolveFastLevenbergMarquardtExt.jl b/ext/NonlinearSolveFastLevenbergMarquardtExt.jl index 2cfb98020..6221883a7 100644 --- a/ext/NonlinearSolveFastLevenbergMarquardtExt.jl +++ b/ext/NonlinearSolveFastLevenbergMarquardtExt.jl @@ -20,11 +20,11 @@ end function SciMLBase.__solve(prob::Union{NonlinearLeastSquaresProblem, NonlinearProblem}, alg::FastLevenbergMarquardtJL, args...; alias_u0 = false, abstol = nothing, reltol = nothing, maxiters = 1000, termination_condition = nothing, kwargs...) - NonlinearSolve.__test_termination_condition(termination_condition, - :FastLevenbergMarquardt) + NonlinearSolve.__test_termination_condition( + termination_condition, :FastLevenbergMarquardt) - fn, u, resid = NonlinearSolve.__construct_extension_f(prob; alias_u0, - can_handle_oop = Val(prob.u0 isa SArray)) + fn, u, resid = NonlinearSolve.__construct_extension_f( + prob; alias_u0, can_handle_oop = Val(prob.u0 isa SArray)) f = if prob.u0 isa SArray @closure (u, p) -> fn(u) else @@ -33,8 +33,8 @@ function SciMLBase.__solve(prob::Union{NonlinearLeastSquaresProblem, NonlinearPr abstol = NonlinearSolve.DEFAULT_TOLERANCE(abstol, eltype(u)) reltol = NonlinearSolve.DEFAULT_TOLERANCE(reltol, eltype(u)) - _jac_fn = NonlinearSolve.__construct_extension_jac(prob, alg, u, resid; alg.autodiff, - can_handle_oop = Val(prob.u0 isa SArray)) + _jac_fn = NonlinearSolve.__construct_extension_jac( + prob, alg, u, resid; alg.autodiff, can_handle_oop = Val(prob.u0 isa SArray)) jac_fn = if prob.u0 isa SArray @closure (u, p) -> _jac_fn(u) else @@ -42,12 +42,12 @@ function SciMLBase.__solve(prob::Union{NonlinearLeastSquaresProblem, NonlinearPr end solver_kwargs = (; xtol = reltol, ftol = reltol, gtol = abstol, maxit = maxiters, - alg.factor, alg.factoraccept, alg.factorreject, alg.minscale, alg.maxscale, - alg.factorupdate, alg.minfactor, alg.maxfactor) + alg.factor, alg.factoraccept, alg.factorreject, alg.minscale, + alg.maxscale, alg.factorupdate, alg.minfactor, alg.maxfactor) if prob.u0 isa SArray - res, fx, info, iter, nfev, njev = FastLM.lmsolve(f, jac_fn, prob.u0; - solver_kwargs...) + res, fx, info, iter, nfev, njev = FastLM.lmsolve( + f, jac_fn, prob.u0; solver_kwargs...) LM, solver = nothing, nothing else J = prob.f.jac_prototype === nothing ? similar(u, length(resid), length(u)) : @@ -55,8 +55,8 @@ function SciMLBase.__solve(prob::Union{NonlinearLeastSquaresProblem, NonlinearPr solver = _fast_lm_solver(alg, u) LM = FastLM.LMWorkspace(u, resid, J) - res, fx, info, iter, nfev, njev, LM, solver = FastLM.lmsolve!(f, jac_fn, LM; - solver, solver_kwargs...) + res, fx, info, iter, nfev, njev, LM, solver = FastLM.lmsolve!( + f, jac_fn, LM; solver, solver_kwargs...) end stats = SciMLBase.NLStats(nfev, njev, -1, -1, iter) diff --git a/ext/NonlinearSolveFixedPointAccelerationExt.jl b/ext/NonlinearSolveFixedPointAccelerationExt.jl index 0c8ff8371..9e7254886 100644 --- a/ext/NonlinearSolveFixedPointAccelerationExt.jl +++ b/ext/NonlinearSolveFixedPointAccelerationExt.jl @@ -4,13 +4,13 @@ using NonlinearSolve, FixedPointAcceleration, SciMLBase function SciMLBase.__solve(prob::NonlinearProblem, alg::FixedPointAccelerationJL, args...; abstol = nothing, maxiters = 1000, alias_u0::Bool = false, - show_trace::Val{PrintReports} = Val(false), termination_condition = nothing, - kwargs...) where {PrintReports} - NonlinearSolve.__test_termination_condition(termination_condition, - :FixedPointAccelerationJL) + show_trace::Val{PrintReports} = Val(false), + termination_condition = nothing, kwargs...) where {PrintReports} + NonlinearSolve.__test_termination_condition( + termination_condition, :FixedPointAccelerationJL) - f, u0, resid = NonlinearSolve.__construct_extension_f(prob; alias_u0, - make_fixed_point = Val(true), force_oop = Val(true)) + f, u0, resid = NonlinearSolve.__construct_extension_f( + prob; alias_u0, make_fixed_point = Val(true), force_oop = Val(true)) tol = NonlinearSolve.DEFAULT_TOLERANCE(abstol, eltype(u0)) sol = fixed_point(f, u0; Algorithm = alg.algorithm, MaxIter = maxiters, MaxM = alg.m, diff --git a/ext/NonlinearSolveLeastSquaresOptimExt.jl b/ext/NonlinearSolveLeastSquaresOptimExt.jl index 6ce6eabd4..42ff7c9d3 100644 --- a/ext/NonlinearSolveLeastSquaresOptimExt.jl +++ b/ext/NonlinearSolveLeastSquaresOptimExt.jl @@ -29,8 +29,8 @@ end function SciMLBase.__init(prob::Union{NonlinearLeastSquaresProblem, NonlinearProblem}, alg::LeastSquaresOptimJL, args...; alias_u0 = false, abstol = nothing, - show_trace::Val{ShT} = Val(false), trace_level = TraceMinimal(), reltol = nothing, - store_trace::Val{StT} = Val(false), maxiters = 1000, + show_trace::Val{ShT} = Val(false), trace_level = TraceMinimal(), + reltol = nothing, store_trace::Val{StT} = Val(false), maxiters = 1000, termination_condition = nothing, kwargs...) where {ShT, StT} NonlinearSolve.__test_termination_condition(termination_condition, :LeastSquaresOptim) @@ -43,13 +43,16 @@ function SciMLBase.__init(prob::Union{NonlinearLeastSquaresProblem, NonlinearPro J = prob.f.jac_prototype, output_length = length(resid)) else g! = NonlinearSolve.__construct_extension_jac(prob, alg, u, resid; alg.autodiff) - lsoprob = LSO.LeastSquaresProblem(; x = u, f!, y = resid, g!, - J = prob.f.jac_prototype, output_length = length(resid)) + lsoprob = LSO.LeastSquaresProblem(; + x = u, f!, y = resid, g!, J = prob.f.jac_prototype, + output_length = length(resid)) end allocated_prob = LSO.LeastSquaresProblemAllocated(lsoprob, _lso_solver(alg)) - return LeastSquaresOptimJLCache(prob, alg, allocated_prob, + return LeastSquaresOptimJLCache(prob, + alg, + allocated_prob, (; x_tol = reltol, f_tol = abstol, g_tol = abstol, iterations = maxiters, show_trace = ShT, store_trace = StT, show_every = trace_level.print_frequency)) end @@ -61,8 +64,8 @@ function SciMLBase.solve!(cache::LeastSquaresOptimJLCache) (res.iterations ≥ maxiters ? ReturnCode.MaxIters : ReturnCode.ConvergenceFailure) stats = SciMLBase.NLStats(res.f_calls, res.g_calls, -1, -1, res.iterations) - return SciMLBase.build_solution(cache.prob, cache.alg, res.minimizer, res.ssr / 2; - retcode, original = res, stats) + return SciMLBase.build_solution( + cache.prob, cache.alg, res.minimizer, res.ssr / 2; retcode, original = res, stats) end end diff --git a/ext/NonlinearSolveMINPACKExt.jl b/ext/NonlinearSolveMINPACKExt.jl index 4465051ea..8be121c13 100644 --- a/ext/NonlinearSolveMINPACKExt.jl +++ b/ext/NonlinearSolveMINPACKExt.jl @@ -3,12 +3,11 @@ module NonlinearSolveMINPACKExt using MINPACK, NonlinearSolve, SciMLBase import FastClosures: @closure -function SciMLBase.__solve(prob::Union{NonlinearLeastSquaresProblem, - NonlinearProblem}, - alg::CMINPACK, args...; abstol = nothing, maxiters = 1000, - alias_u0::Bool = false, show_trace::Val{ShT} = Val(false), - store_trace::Val{StT} = Val(false), termination_condition = nothing, - kwargs...) where {ShT, StT} +function SciMLBase.__solve( + prob::Union{NonlinearLeastSquaresProblem, NonlinearProblem}, alg::CMINPACK, + args...; abstol = nothing, maxiters = 1000, alias_u0::Bool = false, + show_trace::Val{ShT} = Val(false), store_trace::Val{StT} = Val(false), + termination_condition = nothing, kwargs...) where {ShT, StT} NonlinearSolve.__test_termination_condition(termination_condition, :CMINPACK) _f!, u0, resid = NonlinearSolve.__construct_extension_f(prob; alias_u0) @@ -23,13 +22,13 @@ function SciMLBase.__solve(prob::Union{NonlinearLeastSquaresProblem, tol = NonlinearSolve.DEFAULT_TOLERANCE(abstol, eltype(u0)) if alg.autodiff === missing && prob.f.jac === nothing - original = MINPACK.fsolve(f!, u0, m; tol, show_trace, tracing, method, - iterations = maxiters) + original = MINPACK.fsolve( + f!, u0, m; tol, show_trace, tracing, method, iterations = maxiters) else _jac! = NonlinearSolve.__construct_extension_jac(prob, alg, u0, resid; alg.autodiff) jac! = @closure (J, u) -> (_jac!(J, u); Cint(0)) - original = MINPACK.fsolve(f!, jac!, u0, m; tol, show_trace, tracing, method, - iterations = maxiters) + original = MINPACK.fsolve( + f!, jac!, u0, m; tol, show_trace, tracing, method, iterations = maxiters) end u = original.x diff --git a/ext/NonlinearSolveNLsolveExt.jl b/ext/NonlinearSolveNLsolveExt.jl index 64886c021..2af1ab5a3 100644 --- a/ext/NonlinearSolveNLsolveExt.jl +++ b/ext/NonlinearSolveNLsolveExt.jl @@ -2,11 +2,11 @@ module NonlinearSolveNLsolveExt using NonlinearSolve, NLsolve, SciMLBase -function SciMLBase.__solve(prob::NonlinearProblem, alg::NLsolveJL, args...; - abstol = nothing, maxiters = 1000, alias_u0::Bool = false, - termination_condition = nothing, store_trace::Val{StT} = Val(false), - show_trace::Val{ShT} = Val(false), trace_level = TraceMinimal(), - kwargs...) where {StT, ShT} +function SciMLBase.__solve( + prob::NonlinearProblem, alg::NLsolveJL, args...; abstol = nothing, + maxiters = 1000, alias_u0::Bool = false, termination_condition = nothing, + store_trace::Val{StT} = Val(false), show_trace::Val{ShT} = Val(false), + trace_level = TraceMinimal(), kwargs...) where {StT, ShT} NonlinearSolve.__test_termination_condition(termination_condition, :NLsolveJL) f!, u0, resid = NonlinearSolve.__construct_extension_f(prob; alias_u0) @@ -16,8 +16,8 @@ function SciMLBase.__solve(prob::NonlinearProblem, alg::NLsolveJL, args...; else jac! = NonlinearSolve.__construct_extension_jac(prob, alg, u0, resid; alg.autodiff) if prob.f.jac_prototype === nothing - J = similar(u0, promote_type(eltype(u0), eltype(resid)), length(u0), - length(resid)) + J = similar( + u0, promote_type(eltype(u0), eltype(resid)), length(u0), length(resid)) else J = zero(prob.f.jac_prototype) end @@ -30,8 +30,8 @@ function SciMLBase.__solve(prob::NonlinearProblem, alg::NLsolveJL, args...; extended_trace = !(trace_level isa TraceMinimal) || alg.extended_trace original = nlsolve(df, vec(u0); ftol = abstol, iterations = maxiters, alg.method, - store_trace, extended_trace, alg.linesearch, alg.linsolve, alg.factor, - alg.autoscale, alg.m, alg.beta, show_trace) + store_trace, extended_trace, alg.linesearch, alg.linsolve, + alg.factor, alg.autoscale, alg.m, alg.beta, show_trace) f!(vec(resid), original.zero) u = prob.u0 isa Number ? original.zero[1] : reshape(original.zero, size(prob.u0)) diff --git a/ext/NonlinearSolveSIAMFANLEquationsExt.jl b/ext/NonlinearSolveSIAMFANLEquationsExt.jl index c313477df..17fef4922 100644 --- a/ext/NonlinearSolveSIAMFANLEquationsExt.jl +++ b/ext/NonlinearSolveSIAMFANLEquationsExt.jl @@ -25,14 +25,14 @@ end # not interesting here. @inline function __siam_fanl_equations_stats_mapping(method, sol) ((method === :pseudotransient) || (method === :anderson)) && return nothing - return SciMLBase.NLStats(sum(sol.stats.ifun), sum(sol.stats.ijac), 0, 0, - sum(sol.stats.iarm)) + return SciMLBase.NLStats( + sum(sol.stats.ifun), sum(sol.stats.ijac), 0, 0, sum(sol.stats.iarm)) end function SciMLBase.__solve(prob::NonlinearProblem, alg::SIAMFANLEquationsJL, args...; - abstol = nothing, reltol = nothing, alias_u0::Bool = false, maxiters = 1000, - termination_condition = nothing, show_trace::Val{ShT} = Val(false), - kwargs...) where {ShT} + abstol = nothing, reltol = nothing, alias_u0::Bool = false, + maxiters = 1000, termination_condition = nothing, + show_trace::Val{ShT} = Val(false), kwargs...) where {ShT} NonlinearSolve.__test_termination_condition(termination_condition, :SIAMFANLEquationsJL) (; method, delta, linsolve, m, beta) = alg @@ -45,19 +45,19 @@ function SciMLBase.__solve(prob::NonlinearProblem, alg::SIAMFANLEquationsJL, arg if method == :newton sol = nsolsc(f, prob.u0; maxit = maxiters, atol, rtol, printerr = ShT) elseif method == :pseudotransient - sol = ptcsolsc(f, prob.u0; delta0 = delta, maxit = maxiters, atol, rtol, - printerr = ShT) + sol = ptcsolsc( + f, prob.u0; delta0 = delta, maxit = maxiters, atol, rtol, printerr = ShT) elseif method == :secant sol = secant(f, prob.u0; maxit = maxiters, atol, rtol, printerr = ShT) elseif method == :anderson - f_aa, u, _ = NonlinearSolve.__construct_extension_f(prob; alias_u0, - make_fixed_point = Val(true)) - sol = aasol(f_aa, u, m, __zeros_like(u, 1, 2 * m + 4); maxit = maxiters, - atol, rtol, beta) + f_aa, u, _ = NonlinearSolve.__construct_extension_f( + prob; alias_u0, make_fixed_point = Val(true)) + sol = aasol(f_aa, u, m, __zeros_like(u, 1, 2 * m + 4); + maxit = maxiters, atol, rtol, beta) end else - f, u, resid = NonlinearSolve.__construct_extension_f(prob; alias_u0, - make_fixed_point = Val(method == :anderson)) + f, u, resid = NonlinearSolve.__construct_extension_f( + prob; alias_u0, make_fixed_point = Val(method == :anderson)) N = length(u) FS = __zeros_like(u, N) @@ -67,34 +67,34 @@ function SciMLBase.__solve(prob::NonlinearProblem, alg::SIAMFANLEquationsJL, arg JVS = linsolve == :gmres ? __zeros_like(u, N, 3) : __zeros_like(u, N) linsolve_alg = String(linsolve) if method == :newton - sol = nsoli(f, u, FS, JVS; lsolver = linsolve_alg, maxit = maxiters, atol, - rtol, printerr = ShT) + sol = nsoli(f, u, FS, JVS; lsolver = linsolve_alg, + maxit = maxiters, atol, rtol, printerr = ShT) elseif method == :pseudotransient - sol = ptcsoli(f, u, FS, JVS; lsolver = linsolve_alg, maxit = maxiters, - atol, rtol, printerr = ShT) + sol = ptcsoli(f, u, FS, JVS; lsolver = linsolve_alg, + maxit = maxiters, atol, rtol, printerr = ShT) end else if prob.f.jac === nothing && alg.autodiff === missing FPS = __zeros_like(u, N, N) if method == :newton - sol = nsol(f, u, FS, FPS; sham = 1, atol, rtol, maxit = maxiters, - printerr = ShT) + sol = nsol(f, u, FS, FPS; sham = 1, atol, rtol, + maxit = maxiters, printerr = ShT) elseif method == :pseudotransient sol = ptcsol(f, u, FS, FPS; atol, rtol, maxit = maxiters, delta0 = delta, printerr = ShT) elseif method == :anderson - sol = aasol(f, u, m, zeros(T, N, 2 * m + 4); atol, rtol, - maxit = maxiters, beta) + sol = aasol( + f, u, m, zeros(T, N, 2 * m + 4); atol, rtol, maxit = maxiters, beta) end else FPS = prob.f.jac_prototype !== nothing ? zero(prob.f.jac_prototype) : __zeros_like(u, N, N) - jac = NonlinearSolve.__construct_extension_jac(prob, alg, u, resid; - alg.autodiff) + jac = NonlinearSolve.__construct_extension_jac( + prob, alg, u, resid; alg.autodiff) AJ! = @closure (J, u, x) -> jac(J, x) if method == :newton - sol = nsol(f, u, FS, FPS, AJ!; sham = 1, atol, rtol, maxit = maxiters, - printerr = ShT) + sol = nsol(f, u, FS, FPS, AJ!; sham = 1, atol, + rtol, maxit = maxiters, printerr = ShT) elseif method == :pseudotransient sol = ptcsol(f, u, FS, FPS, AJ!; atol, rtol, maxit = maxiters, delta0 = delta, printerr = ShT) diff --git a/ext/NonlinearSolveSpeedMappingExt.jl b/ext/NonlinearSolveSpeedMappingExt.jl index 23f1cba98..286b2577c 100644 --- a/ext/NonlinearSolveSpeedMappingExt.jl +++ b/ext/NonlinearSolveSpeedMappingExt.jl @@ -3,13 +3,13 @@ module NonlinearSolveSpeedMappingExt using NonlinearSolve, SciMLBase, SpeedMapping function SciMLBase.__solve(prob::NonlinearProblem, alg::SpeedMappingJL, args...; - abstol = nothing, maxiters = 1000, alias_u0::Bool = false, maxtime = nothing, - store_trace::Val{store_info} = Val(false), termination_condition = nothing, - kwargs...) where {store_info} + abstol = nothing, maxiters = 1000, alias_u0::Bool = false, + maxtime = nothing, store_trace::Val{store_info} = Val(false), + termination_condition = nothing, kwargs...) where {store_info} NonlinearSolve.__test_termination_condition(termination_condition, :SpeedMappingJL) - m!, u, resid = NonlinearSolve.__construct_extension_f(prob; alias_u0, - make_fixed_point = Val(true)) + m!, u, resid = NonlinearSolve.__construct_extension_f( + prob; alias_u0, make_fixed_point = Val(true)) tol = NonlinearSolve.DEFAULT_TOLERANCE(abstol, eltype(u)) time_limit = ifelse(maxtime === nothing, alg.time_limit, maxtime) diff --git a/src/NonlinearSolve.jl b/src/NonlinearSolve.jl index de46fe153..e30ce530c 100644 --- a/src/NonlinearSolve.jl +++ b/src/NonlinearSolve.jl @@ -90,15 +90,15 @@ include("default.jl") push!(probs_nls, NonlinearProblem(fn, T.(u0), T(2))) end - nls_algs = (NewtonRaphson(), TrustRegion(), LevenbergMarquardt(), PseudoTransient(), - Broyden(), Klement(), DFSane(), nothing) + nls_algs = (NewtonRaphson(), TrustRegion(), LevenbergMarquardt(), + PseudoTransient(), Broyden(), Klement(), DFSane(), nothing) probs_nlls = NonlinearLeastSquaresProblem[] nlfuncs = ((NonlinearFunction{false}((u, p) -> (u .^ 2 .- p)[1:1]), [0.1, 0.0]), (NonlinearFunction{false}((u, p) -> vcat(u .* u .- p, u .* u .- p)), [0.1, 0.1]), ( - NonlinearFunction{true}((du, u, p) -> du[1] = u[1] * u[1] - p, - resid_prototype = zeros(1)), + NonlinearFunction{true}( + (du, u, p) -> du[1] = u[1] * u[1] - p, resid_prototype = zeros(1)), [0.1, 0.0]), ( NonlinearFunction{true}((du, u, p) -> du .= vcat(u .* u .- p, u .* u .- p), @@ -111,8 +111,8 @@ include("default.jl") (NonlinearFunction{false}((u, p) -> vcat(u .* u .- p, u .* u .- p)), Float32[0.1, 0.1]), ( - NonlinearFunction{true}((du, u, p) -> du[1] = u[1] * u[1] - p, - resid_prototype = zeros(Float32, 1)), + NonlinearFunction{true}( + (du, u, p) -> du[1] = u[1] * u[1] - p, resid_prototype = zeros(Float32, 1)), Float32[0.1, 0.0]), ( NonlinearFunction{true}((du, u, p) -> du .= vcat(u .* u .- p, u .* u .- p), @@ -140,8 +140,8 @@ end # Core Algorithms export NewtonRaphson, PseudoTransient, Klement, Broyden, LimitedMemoryBroyden, DFSane export GaussNewton, LevenbergMarquardt, TrustRegion -export NonlinearSolvePolyAlgorithm, - RobustMultiNewton, FastShortcutNonlinearPolyalg, FastShortcutNLLSPolyalg +export NonlinearSolvePolyAlgorithm, RobustMultiNewton, FastShortcutNonlinearPolyalg, + FastShortcutNLLSPolyalg # Extension Algorithms export LeastSquaresOptimJL, FastLevenbergMarquardtJL, CMINPACK, NLsolveJL, @@ -151,8 +151,7 @@ export LeastSquaresOptimJL, FastLevenbergMarquardtJL, CMINPACK, NLsolveJL, export GeneralizedFirstOrderAlgorithm, ApproximateJacobianSolveAlgorithm, GeneralizedDFSane # Descent Algorithms -export NewtonDescent, SteepestDescent, Dogleg, DampedNewtonDescent, - GeodesicAcceleration +export NewtonDescent, SteepestDescent, Dogleg, DampedNewtonDescent, GeodesicAcceleration # Globalization ## Line Search Algorithms diff --git a/src/abstract_types.jl b/src/abstract_types.jl index 1b30f2b9f..b7127cb2c 100644 --- a/src/abstract_types.jl +++ b/src/abstract_types.jl @@ -16,15 +16,15 @@ squares solver. ### `__internal_init` specification ```julia +__internal_init(prob::NonlinearProblem{uType, iip}, alg::AbstractDescentAlgorithm, J, + fu, u; pre_inverted::Val{INV} = Val(false), linsolve_kwargs = (;), + abstol = nothing, reltol = nothing, alias_J::Bool = true, + shared::Val{N} = Val(1), kwargs...) where {INV, N, uType, iip} --> AbstractDescentCache + __internal_init( - prob::NonlinearProblem{uType, iip}, alg::AbstractDescentAlgorithm, J, fu, u; - pre_inverted::Val{INV} = Val(false), linsolve_kwargs = (;), abstol = nothing, - reltol = nothing, alias_J::Bool = true, shared::Val{N} = Val(1), - kwargs...) where {INV, N, uType, iip} --> AbstractDescentCache - -__internal_init(prob::NonlinearLeastSquaresProblem{uType, iip}, - alg::AbstractDescentAlgorithm, J, fu, u; pre_inverted::Val{INV} = Val(false), - linsolve_kwargs = (;), abstol = nothing, reltol = nothing, alias_J::Bool = true, + prob::NonlinearLeastSquaresProblem{uType, iip}, alg::AbstractDescentAlgorithm, + J, fu, u; pre_inverted::Val{INV} = Val(false), linsolve_kwargs = (;), + abstol = nothing, reltol = nothing, alias_J::Bool = true, shared::Val{N} = Val(1), kwargs...) where {INV, N, uType, iip} --> AbstractDescentCache ``` @@ -66,8 +66,8 @@ Abstract Type for all Descent Caches. ### `__internal_solve!` specification ```julia -δu, success, intermediates = __internal_solve!(cache::AbstractDescentCache, J, fu, u, - idx::Val; skip_solve::Bool = false, kwargs...) +δu, success, intermediates = __internal_solve!( + cache::AbstractDescentCache, J, fu, u, idx::Val; skip_solve::Bool = false, kwargs...) ``` - `J`: Jacobian or Inverse Jacobian (if `pre_inverted = Val(true)`). @@ -119,10 +119,10 @@ Abstract Type for all Line Search Algorithms used in NonlinearSolve.jl. ### `__internal_init` specification ```julia -__internal_init(prob::AbstractNonlinearProblem, - alg::AbstractNonlinearSolveLineSearchAlgorithm, f::F, fu, u, p, args...; - internalnorm::IN = DEFAULT_NORM, - kwargs...) where {F, IN} --> AbstractNonlinearSolveLineSearchCache +__internal_init( + prob::AbstractNonlinearProblem, alg::AbstractNonlinearSolveLineSearchAlgorithm, f::F, + fu, u, p, args...; internalnorm::IN = DEFAULT_NORM, kwargs...) where {F, IN} --> +AbstractNonlinearSolveLineSearchCache ``` """ abstract type AbstractNonlinearSolveLineSearchAlgorithm end @@ -145,8 +145,8 @@ Returns 2 values: """ abstract type AbstractNonlinearSolveLineSearchCache end -function reinit_cache!(cache::AbstractNonlinearSolveLineSearchCache, args...; p = cache.p, - kwargs...) +function reinit_cache!( + cache::AbstractNonlinearSolveLineSearchCache, args...; p = cache.p, kwargs...) cache.nf[] = 0 cache.p = p end @@ -233,10 +233,9 @@ Abstract Type for Damping Functions in DampedNewton. ### `__internal_init` specification ```julia -__internal_init( - prob::AbstractNonlinearProblem, f::AbstractDampingFunction, initial_damping, - J, fu, u, args...; internal_norm = DEFAULT_NORM, - kwargs...) --> AbstractDampingFunctionCache +__internal_init(prob::AbstractNonlinearProblem, f::AbstractDampingFunction, initial_damping, + J, fu, u, args...; internal_norm = DEFAULT_NORM, kwargs...) --> +AbstractDampingFunctionCache ``` Returns a [`AbstractDampingFunctionCache`](@ref). @@ -318,9 +317,8 @@ Abstract Type for all Jacobian Initialization Algorithms used in NonlinearSolve. ### `__internal_init` specification ```julia -__internal_init(prob::AbstractNonlinearProblem, alg::AbstractJacobianInitialization, - solver, f::F, fu, u, p; linsolve = missing, internalnorm::IN = DEFAULT_NORM, - kwargs...) +__internal_init(prob::AbstractNonlinearProblem, alg::AbstractJacobianInitialization, solver, + f::F, fu, u, p; linsolve = missing, internalnorm::IN = DEFAULT_NORM, kwargs...) ``` Returns a [`NonlinearSolve.InitializedApproximateJacobianCache`](@ref). @@ -353,10 +351,10 @@ Abstract Type for all Approximate Jacobian Update Rules used in NonlinearSolve.j ### `__internal_init` specification ```julia -__internal_init(prob::AbstractNonlinearProblem, - alg::AbstractApproximateJacobianUpdateRule, J, fu, u, du, args...; - internalnorm::F = DEFAULT_NORM, - kwargs...) where {F} --> AbstractApproximateJacobianUpdateRuleCache{INV} +__internal_init( + prob::AbstractNonlinearProblem, alg::AbstractApproximateJacobianUpdateRule, J, + fu, u, du, args...; internalnorm::F = DEFAULT_NORM, kwargs...) where {F} --> +AbstractApproximateJacobianUpdateRuleCache{INV} ``` """ abstract type AbstractApproximateJacobianUpdateRule{INV} end @@ -375,8 +373,8 @@ Abstract Type for all Approximate Jacobian Update Rule Caches used in NonlinearS ### `__internal_solve!` specification ```julia -__internal_solve!(cache::AbstractApproximateJacobianUpdateRuleCache, J, fu, u, du; - kwargs...) --> J / J⁻¹ +__internal_solve!( + cache::AbstractApproximateJacobianUpdateRuleCache, J, fu, u, du; kwargs...) --> J / J⁻¹ ``` """ abstract type AbstractApproximateJacobianUpdateRuleCache{INV} end @@ -391,8 +389,8 @@ Condition for resetting the Jacobian in Quasi-Newton's methods. ### `__internal_init` specification ```julia -__internal_init(alg::AbstractResetCondition, J, fu, u, du, args...; - kwargs...) --> ResetCache +__internal_init(alg::AbstractResetCondition, J, fu, u, du, args...; kwargs...) --> +ResetCache ``` ### `__internal_solve!` specification @@ -411,9 +409,9 @@ Abstract Type for all Trust Region Methods used in NonlinearSolve.jl. ### `__internal_init` specification ```julia -__internal_init(prob::AbstractNonlinearProblem, alg::AbstractTrustRegionMethod, - f::F, fu, u, p, args...; internalnorm::IF = DEFAULT_NORM, - kwargs...) where {F, IF} --> AbstractTrustRegionMethodCache +__internal_init(prob::AbstractNonlinearProblem, alg::AbstractTrustRegionMethod, f::F, fu, u, + p, args...; internalnorm::IF = DEFAULT_NORM, kwargs...) where {F, IF} --> +AbstractTrustRegionMethodCache ``` """ abstract type AbstractTrustRegionMethod end @@ -468,8 +466,8 @@ abstract type AbstractNonlinearSolveTraceLevel end # Default Printing for aType in (AbstractTrustRegionMethod, AbstractNonlinearSolveLineSearchAlgorithm, - AbstractResetCondition, AbstractApproximateJacobianUpdateRule, AbstractDampingFunction, - AbstractNonlinearSolveExtensionAlgorithm) + AbstractResetCondition, AbstractApproximateJacobianUpdateRule, + AbstractDampingFunction, AbstractNonlinearSolveExtensionAlgorithm) @eval function Base.show(io::IO, alg::$(aType)) print(io, "$(nameof(typeof(alg)))()") end diff --git a/src/algorithms/broyden.jl b/src/algorithms/broyden.jl index bd3b490e4..679897efa 100644 --- a/src/algorithms/broyden.jl +++ b/src/algorithms/broyden.jl @@ -28,10 +28,9 @@ search. useful for specific problems, but whether it will work may depend strongly on the problem """ -function Broyden(; - max_resets = 100, linesearch = NoLineSearch(), reset_tolerance = nothing, - init_jacobian::Val{IJ} = Val(:identity), autodiff = nothing, alpha = nothing, - update_rule::Val{UR} = Val(:good_broyden)) where {IJ, UR} +function Broyden(; max_resets = 100, linesearch = NoLineSearch(), reset_tolerance = nothing, + init_jacobian::Val{IJ} = Val(:identity), autodiff = nothing, + alpha = nothing, update_rule::Val{UR} = Val(:good_broyden)) where {IJ, UR} if IJ === :identity if UR === :diagonal initialization = IdentityInitialization(alpha, DiagonalStructure()) @@ -56,9 +55,9 @@ function Broyden(; or `:diagonal`")) end - return ApproximateJacobianSolveAlgorithm{IJ === :true_jacobian, :Broyden}(; linesearch, - descent = NewtonDescent(), update_rule, max_resets, initialization, - reinit_rule = NoChangeInStateReset(; reset_tolerance)) + return ApproximateJacobianSolveAlgorithm{IJ === :true_jacobian, :Broyden}(; + linesearch, descent = NewtonDescent(), update_rule, max_resets, + initialization, reinit_rule = NoChangeInStateReset(; reset_tolerance)) end # Checks for no significant change for `nsteps` @@ -107,8 +106,8 @@ function __internal_init(alg::NoChangeInStateReset, J, fu, u, du, args...; kwarg end T = real(eltype(u)) tol = alg.reset_tolerance === nothing ? eps(T)^(3 // 4) : T(alg.reset_tolerance) - return NoChangeInStateResetCache(dfu, tol, alg.check_du, alg.check_dfu, alg.nsteps, 0, - 0) + return NoChangeInStateResetCache( + dfu, tol, alg.check_du, alg.check_dfu, alg.nsteps, 0, 0) end function __internal_solve!(cache::NoChangeInStateResetCache, J, fu, u, du) @@ -170,8 +169,8 @@ Broyden Update Rule corresponding to "good broyden's method" [broyden1965class]( end function __internal_init(prob::AbstractNonlinearProblem, - alg::Union{GoodBroydenUpdateRule, BadBroydenUpdateRule}, J, fu, u, du, args...; - internalnorm::F = DEFAULT_NORM, kwargs...) where {F} + alg::Union{GoodBroydenUpdateRule, BadBroydenUpdateRule}, J, fu, u, + du, args...; internalnorm::F = DEFAULT_NORM, kwargs...) where {F} @bb J⁻¹dfu = similar(u) @bb dfu = copy(fu) if alg isa GoodBroydenUpdateRule || J isa Diagonal @@ -206,8 +205,8 @@ function __internal_solve!(cache::BroydenUpdateRuleCache{mode}, J⁻¹, fu, u, d return J⁻¹ end -function __internal_solve!(cache::BroydenUpdateRuleCache{mode}, J⁻¹::Diagonal, fu, u, - du) where {mode} +function __internal_solve!( + cache::BroydenUpdateRuleCache{mode}, J⁻¹::Diagonal, fu, u, du) where {mode} T = eltype(u) @bb @. cache.dfu = fu - cache.dfu J⁻¹_diag = _restructure(cache.dfu, diag(J⁻¹)) diff --git a/src/algorithms/dfsane.jl b/src/algorithms/dfsane.jl index 9122fb6a1..b42544055 100644 --- a/src/algorithms/dfsane.jl +++ b/src/algorithms/dfsane.jl @@ -20,7 +20,7 @@ function DFSane(; σ_min = 1 // 10^10, σ_max = 1e10, σ_1 = 1, M::Int = 10, γ τ_min = 1 // 10, τ_max = 1 // 2, n_exp::Int = 2, max_inner_iterations::Int = 100, η_strategy::ETA = (fn_1, n, x_n, f_n) -> fn_1 / n^2) where {ETA} linesearch = RobustNonMonotoneLineSearch(; - gamma = γ, sigma_1 = σ_1, M, tau_min = τ_min, - tau_max = τ_max, n_exp, η_strategy, maxiters = max_inner_iterations) + gamma = γ, sigma_1 = σ_1, M, tau_min = τ_min, tau_max = τ_max, + n_exp, η_strategy, maxiters = max_inner_iterations) return GeneralizedDFSane{:DFSane}(linesearch, σ_min, σ_max, nothing) end diff --git a/src/algorithms/extension_algs.jl b/src/algorithms/extension_algs.jl index eeed0ea5b..673af5314 100644 --- a/src/algorithms/extension_algs.jl +++ b/src/algorithms/extension_algs.jl @@ -80,10 +80,10 @@ see the documentation for `FastLevenbergMarquardt.jl`. maxfactor end -function FastLevenbergMarquardtJL(linsolve::Symbol = :cholesky; factor = 1e-6, - factoraccept = 13.0, factorreject = 3.0, factorupdate = :marquardt, - minscale = 1e-12, maxscale = 1e16, minfactor = 1e-28, maxfactor = 1e32, - autodiff = nothing) +function FastLevenbergMarquardtJL( + linsolve::Symbol = :cholesky; factor = 1e-6, factoraccept = 13.0, + factorreject = 3.0, factorupdate = :marquardt, minscale = 1e-12, + maxscale = 1e16, minfactor = 1e-28, maxfactor = 1e32, autodiff = nothing) @assert linsolve in (:qr, :cholesky) @assert factorupdate in (:marquardt, :nielson) @@ -148,8 +148,8 @@ NonlinearLeastSquaresProblem. autodiff end -function CMINPACK(; show_trace = missing, tracing = missing, method::Symbol = :auto, - autodiff = missing) +function CMINPACK(; show_trace = missing, tracing = missing, + method::Symbol = :auto, autodiff = missing) if Base.get_extension(@__MODULE__, :NonlinearSolveMINPACKExt) === nothing error("CMINPACK requires MINPACK.jl to be loaded") end @@ -235,8 +235,8 @@ end function NLsolveJL(; method = :trust_region, autodiff = :central, store_trace = missing, extended_trace = missing, linesearch = LineSearches.Static(), - linsolve = (x, A, b) -> copyto!(x, A \ b), factor = 1.0, autoscale = true, m = 10, - beta = one(Float64), show_trace = missing) + linsolve = (x, A, b) -> copyto!(x, A \ b), factor = 1.0, + autoscale = true, m = 10, beta = one(Float64), show_trace = missing) if Base.get_extension(@__MODULE__, :NonlinearSolveNLsolveExt) === nothing error("NLsolveJL requires NLsolve.jl to be loaded") end @@ -252,22 +252,20 @@ function NLsolveJL(; method = :trust_region, autodiff = :central, store_trace = end if store_trace !== missing - Base.depwarn( - "`store_trace` for NLsolveJL has been deprecated and will be removed \ - in v4. Use the `store_trace` keyword argument via the logging API \ - https://docs.sciml.ai/NonlinearSolve/stable/basics/Logging/ \ - instead.", + Base.depwarn("`store_trace` for NLsolveJL has been deprecated and will be removed \ + in v4. Use the `store_trace` keyword argument via the logging API \ + https://docs.sciml.ai/NonlinearSolve/stable/basics/Logging/ \ + instead.", :NLsolveJL) else store_trace = false end if extended_trace !== missing - Base.depwarn( - "`extended_trace` for NLsolveJL has been deprecated and will be \ - removed in v4. Use the `trace_level = TraceAll()` keyword argument \ - via the logging API \ - https://docs.sciml.ai/NonlinearSolve/stable/basics/Logging/ instead.", + Base.depwarn("`extended_trace` for NLsolveJL has been deprecated and will be \ + removed in v4. Use the `trace_level = TraceAll()` keyword argument \ + via the logging API \ + https://docs.sciml.ai/NonlinearSolve/stable/basics/Logging/ instead.", :NLsolveJL) else extended_trace = false @@ -277,8 +275,8 @@ function NLsolveJL(; method = :trust_region, autodiff = :central, store_trace = error("`autodiff` must be `:central` or `:forward`.") end - return NLsolveJL(method, autodiff, store_trace, extended_trace, linesearch, linsolve, - factor, autoscale, m, beta, show_trace) + return NLsolveJL(method, autodiff, store_trace, extended_trace, linesearch, + linsolve, factor, autoscale, m, beta, show_trace) end """ @@ -368,9 +366,9 @@ problems as well. condition_number_threshold end -function FixedPointAccelerationJL(; algorithm = :Anderson, m = missing, - condition_number_threshold = missing, extrapolation_period = missing, - replace_invalids = :NoAction, dampening = 1.0) +function FixedPointAccelerationJL(; + algorithm = :Anderson, m = missing, condition_number_threshold = missing, + extrapolation_period = missing, replace_invalids = :NoAction, dampening = 1.0) if Base.get_extension(@__MODULE__, :NonlinearSolveFixedPointAccelerationExt) === nothing error("FixedPointAccelerationJL requires FixedPointAcceleration.jl to be loaded") end @@ -447,8 +445,8 @@ end autodiff end -function SIAMFANLEquationsJL(; method = :newton, delta = 1e-3, linsolve = nothing, m = 0, - beta = 1.0, autodiff = missing) +function SIAMFANLEquationsJL(; method = :newton, delta = 1e-3, linsolve = nothing, + m = 0, beta = 1.0, autodiff = missing) if Base.get_extension(@__MODULE__, :NonlinearSolveSIAMFANLEquationsExt) === nothing error("SIAMFANLEquationsJL requires SIAMFANLEquations.jl to be loaded") end diff --git a/src/algorithms/gauss_newton.jl b/src/algorithms/gauss_newton.jl index 1e6384788..4e826c130 100644 --- a/src/algorithms/gauss_newton.jl +++ b/src/algorithms/gauss_newton.jl @@ -9,6 +9,6 @@ for large-scale and numerically-difficult nonlinear least squares problems. function GaussNewton(; concrete_jac = nothing, linsolve = nothing, precs = DEFAULT_PRECS, linesearch = NoLineSearch(), vjp_autodiff = nothing, autodiff = nothing) descent = NewtonDescent(; linsolve, precs) - return GeneralizedFirstOrderAlgorithm(; concrete_jac, name = :GaussNewton, - descent, jacobian_ad = autodiff, reverse_ad = vjp_autodiff) + return GeneralizedFirstOrderAlgorithm(; concrete_jac, name = :GaussNewton, descent, + jacobian_ad = autodiff, reverse_ad = vjp_autodiff) end diff --git a/src/algorithms/klement.jl b/src/algorithms/klement.jl index 8aa94fbd1..428dc382c 100644 --- a/src/algorithms/klement.jl +++ b/src/algorithms/klement.jl @@ -47,8 +47,9 @@ function Klement(; max_resets::Int = 100, linsolve = nothing, alpha = nothing, CJ = IJ === :true_jacobian || IJ === :true_jacobian_diagonal - return ApproximateJacobianSolveAlgorithm{CJ, :Klement}(; linesearch, - descent = NewtonDescent(; linsolve, precs), update_rule = KlementUpdateRule(), + return ApproximateJacobianSolveAlgorithm{CJ, :Klement}(; + linesearch, descent = NewtonDescent(; linsolve, precs), + update_rule = KlementUpdateRule(), reinit_rule = IllConditionedJacobianReset(), max_resets, initialization) end @@ -99,8 +100,8 @@ Update rule for [`Klement`](@ref). fu_cache end -function __internal_init(prob::AbstractNonlinearProblem, alg::KlementUpdateRule, J, fu, u, - du, args...; kwargs...) +function __internal_init(prob::AbstractNonlinearProblem, alg::KlementUpdateRule, + J, fu, u, du, args...; kwargs...) @bb Jdu = similar(fu) if J isa Diagonal || J isa Number J_cache, J_cache_2, Jdu_cache = nothing, nothing, nothing @@ -125,7 +126,9 @@ function __internal_solve!(cache::KlementUpdateRuleCache, J_::Diagonal, fu, u, d J = _restructure(u, diag(J_)) @bb @. cache.Jdu = (J^2) * (du^2) @bb @. J += ((fu - cache.fu_cache - J * du) / - ifelse(iszero(cache.Jdu), T(1e-5), cache.Jdu)) * du * (J^2) + ifelse(iszero(cache.Jdu), T(1e-5), cache.Jdu)) * + du * + (J^2) @bb copyto!(cache.fu_cache, fu) return Diagonal(vec(J)) end diff --git a/src/algorithms/lbroyden.jl b/src/algorithms/lbroyden.jl index 43b362819..499cc7c34 100644 --- a/src/algorithms/lbroyden.jl +++ b/src/algorithms/lbroyden.jl @@ -19,9 +19,12 @@ function LimitedMemoryBroyden(; max_resets::Int = 3, linesearch = NoLineSearch() threshold::Union{Val, Int} = Val(10), reset_tolerance = nothing, alpha = nothing) threshold isa Int && (threshold = Val(threshold)) return ApproximateJacobianSolveAlgorithm{false, :LimitedMemoryBroyden}(; linesearch, - descent = NewtonDescent(), update_rule = GoodBroydenUpdateRule(), max_resets, - initialization = BroydenLowRankInitialization{_unwrap_val(threshold)}(alpha, - threshold), reinit_rule = NoChangeInStateReset(; reset_tolerance)) + descent = NewtonDescent(), + update_rule = GoodBroydenUpdateRule(), + max_resets, + initialization = BroydenLowRankInitialization{_unwrap_val(threshold)}( + alpha, threshold), + reinit_rule = NoChangeInStateReset(; reset_tolerance)) end """ @@ -38,13 +41,13 @@ end jacobian_initialized_preinverted(::BroydenLowRankInitialization) = true -function __internal_init(prob::AbstractNonlinearProblem, - alg::BroydenLowRankInitialization{T}, solver, f::F, fu, u, p; maxiters = 1000, +function __internal_init( + prob::AbstractNonlinearProblem, alg::BroydenLowRankInitialization{T}, + solver, f::F, fu, u, p; maxiters = 1000, internalnorm::IN = DEFAULT_NORM, kwargs...) where {T, F, IN} if u isa Number # Use the standard broyden - return __internal_init(prob, IdentityInitialization(true, FullStructure()), solver, - f, fu, u, - p; maxiters, kwargs...) + return __internal_init(prob, IdentityInitialization(true, FullStructure()), + solver, f, fu, u, p; maxiters, kwargs...) end # Pay to cost of slightly more allocations to prevent type-instability for StaticArrays α = inv(__initial_alpha(alg.alpha, u, fu, internalnorm)) @@ -54,13 +57,12 @@ function __internal_init(prob::AbstractNonlinearProblem, threshold = min(_unwrap_val(alg.threshold), maxiters) J = BroydenLowRankJacobian(fu, u; threshold, alpha = α) end - return InitializedApproximateJacobianCache(J, FullStructure(), alg, nothing, true, - internalnorm) + return InitializedApproximateJacobianCache( + J, FullStructure(), alg, nothing, true, internalnorm) end function (cache::InitializedApproximateJacobianCache)( - alg::BroydenLowRankInitialization, fu, - u) + alg::BroydenLowRankInitialization, fu, u) α = __initial_alpha(alg.alpha, u, fu, cache.internalnorm) cache.J.idx = 0 cache.J.alpha = inv(α) @@ -97,14 +99,15 @@ end for op in (:adjoint, :transpose) # FIXME: adjoint might be a problem here. Fix if a complex number issue shows up @eval function Base.$(op)(operator::BroydenLowRankJacobian{T}) where {T} - return BroydenLowRankJacobian{T}(operator.Vᵀ, operator.U, - operator.idx, operator.cache, operator.alpha) + return BroydenLowRankJacobian{T}( + operator.Vᵀ, operator.U, operator.idx, operator.cache, operator.alpha) end end # Storing the transpose to ensure contiguous memory on splicing -function BroydenLowRankJacobian(fu::StaticArray{S2, T2}, u::StaticArray{S1, T1}; - alpha = true, threshold::Val{Th} = Val(10)) where {S1, S2, T1, T2, Th} +function BroydenLowRankJacobian( + fu::StaticArray{S2, T2}, u::StaticArray{S1, T1}; alpha = true, + threshold::Val{Th} = Val(10)) where {S1, S2, T1, T2, Th} T = promote_type(T1, T2) fuSize, uSize = Size(fu), Size(u) U = MArray{Tuple{prod(fuSize), Th}, T}(undef) @@ -156,8 +159,8 @@ function LinearAlgebra.mul!(y::AbstractVector, x::AbstractVector, J::BroydenLowR return y end -function LinearAlgebra.mul!(J::BroydenLowRankJacobian, u, - vᵀ::LinearAlgebra.AdjOrTransAbsVec, α::Bool, β::Bool) +function LinearAlgebra.mul!( + J::BroydenLowRankJacobian, u, vᵀ::LinearAlgebra.AdjOrTransAbsVec, α::Bool, β::Bool) @assert α & β idx_update = mod1(J.idx + 1, size(J.U, 2)) copyto!(@view(J.U[:, idx_update]), _vec(u)) diff --git a/src/algorithms/levenberg_marquardt.jl b/src/algorithms/levenberg_marquardt.jl index adb7ee9a4..510226d7d 100644 --- a/src/algorithms/levenberg_marquardt.jl +++ b/src/algorithms/levenberg_marquardt.jl @@ -30,11 +30,12 @@ nonlinear systems. For the remaining arguments, see [`GeodesicAcceleration`](@ref) and [`NonlinearSolve.LevenbergMarquardtTrustRegion`](@ref) documentations. """ -function LevenbergMarquardt(; concrete_jac = missing, linsolve = nothing, - precs = DEFAULT_PRECS, damping_initial::Real = 1.0, α_geodesic::Real = 0.75, +function LevenbergMarquardt(; + concrete_jac = missing, linsolve = nothing, precs = DEFAULT_PRECS, + damping_initial::Real = 1.0, α_geodesic::Real = 0.75, damping_increase_factor::Real = 2.0, damping_decrease_factor::Real = 3.0, - finite_diff_step_geodesic = 0.1, b_uphill::Real = 1.0, autodiff = nothing, - min_damping_D::Real = 1e-8, disable_geodesic = False) + finite_diff_step_geodesic = 0.1, b_uphill::Real = 1.0, + autodiff = nothing, min_damping_D::Real = 1e-8, disable_geodesic = False) if concrete_jac !== missing Base.depwarn("The `concrete_jac` keyword argument is deprecated and will be \ removed in v0.4. This kwarg doesn't make sense (and is currently \ @@ -43,15 +44,16 @@ function LevenbergMarquardt(; concrete_jac = missing, linsolve = nothing, :LevenbergMarquardt) end - descent = DampedNewtonDescent(; linsolve, precs, initial_damping = damping_initial, - damping_fn = LevenbergMarquardtDampingFunction(damping_increase_factor, - damping_decrease_factor, min_damping_D)) + descent = DampedNewtonDescent(; linsolve, + precs, + initial_damping = damping_initial, + damping_fn = LevenbergMarquardtDampingFunction( + damping_increase_factor, damping_decrease_factor, min_damping_D)) if disable_geodesic === False descent = GeodesicAcceleration(descent, finite_diff_step_geodesic, α_geodesic) end trustregion = LevenbergMarquardtTrustRegion(b_uphill) - return GeneralizedFirstOrderAlgorithm(; - concrete_jac = true, name = :LevenbergMarquardt, + return GeneralizedFirstOrderAlgorithm(; concrete_jac = true, name = :LevenbergMarquardt, trustregion, descent, jacobian_ad = autodiff) end @@ -87,21 +89,22 @@ function reinit_cache!(cache::LevenbergMarquardtDampingCache, args...; kwargs... cache.J_damped = cache.λ .* cache.DᵀD end -function requires_normal_form_jacobian(::Union{LevenbergMarquardtDampingFunction, - LevenbergMarquardtDampingCache}) +function requires_normal_form_jacobian(::Union{ + LevenbergMarquardtDampingFunction, LevenbergMarquardtDampingCache}) return false end -function requires_normal_form_rhs(::Union{LevenbergMarquardtDampingFunction, - LevenbergMarquardtDampingCache}) +function requires_normal_form_rhs(::Union{ + LevenbergMarquardtDampingFunction, LevenbergMarquardtDampingCache}) return false end -function returns_norm_form_damping(::Union{LevenbergMarquardtDampingFunction, - LevenbergMarquardtDampingCache}) +function returns_norm_form_damping(::Union{ + LevenbergMarquardtDampingFunction, LevenbergMarquardtDampingCache}) return true end -function __internal_init(prob::AbstractNonlinearProblem, - f::LevenbergMarquardtDampingFunction, initial_damping, J, fu, u, ::Val{NF}; +function __internal_init( + prob::AbstractNonlinearProblem, f::LevenbergMarquardtDampingFunction, + initial_damping, J, fu, u, ::Val{NF}; internalnorm::F = DEFAULT_NORM, kwargs...) where {F, NF} T = promote_type(eltype(u), eltype(fu)) DᵀD = __init_diagonal(u, T(f.min_damping)) @@ -111,15 +114,15 @@ function __internal_init(prob::AbstractNonlinearProblem, @bb J_diag_cache = similar(u) end J_damped = T(initial_damping) .* DᵀD - return LevenbergMarquardtDampingCache(T(f.increase_factor), T(f.decrease_factor), - T(f.min_damping), T(f.increase_factor), T(initial_damping), DᵀD, J_diag_cache, - J_damped, f, T(initial_damping)) + return LevenbergMarquardtDampingCache( + T(f.increase_factor), T(f.decrease_factor), T(f.min_damping), T(f.increase_factor), + T(initial_damping), DᵀD, J_diag_cache, J_damped, f, T(initial_damping)) end (damping::LevenbergMarquardtDampingCache)(::Nothing) = damping.J_damped -function __internal_solve!(damping::LevenbergMarquardtDampingCache, J, fu, ::Val{false}; - kwargs...) +function __internal_solve!( + damping::LevenbergMarquardtDampingCache, J, fu, ::Val{false}; kwargs...) if __can_setindex(damping.J_diag_cache) sum!(abs2, _vec(damping.J_diag_cache), J') elseif damping.J_diag_cache isa Number @@ -132,8 +135,8 @@ function __internal_solve!(damping::LevenbergMarquardtDampingCache, J, fu, ::Val return damping.J_damped end -function __internal_solve!(damping::LevenbergMarquardtDampingCache, JᵀJ, fu, ::Val{true}; - kwargs...) +function __internal_solve!( + damping::LevenbergMarquardtDampingCache, JᵀJ, fu, ::Val{true}; kwargs...) damping.DᵀD = __update_LM_diagonal!!(damping.DᵀD, JᵀJ) @bb @. damping.J_damped = damping.λ * damping.DᵀD return damping.J_damped diff --git a/src/algorithms/pseudo_transient.jl b/src/algorithms/pseudo_transient.jl index 957cfc904..fda39a3b9 100644 --- a/src/algorithms/pseudo_transient.jl +++ b/src/algorithms/pseudo_transient.jl @@ -20,8 +20,8 @@ function PseudoTransient(; concrete_jac = nothing, linsolve = nothing, precs = DEFAULT_PRECS, autodiff = nothing, alpha_initial = 1e-3) descent = DampedNewtonDescent(; linsolve, precs, initial_damping = alpha_initial, damping_fn = SwitchedEvolutionRelaxation()) - return GeneralizedFirstOrderAlgorithm(; concrete_jac, - name = :PseudoTransient, linesearch, descent, jacobian_ad = autodiff) + return GeneralizedFirstOrderAlgorithm(; + concrete_jac, name = :PseudoTransient, linesearch, descent, jacobian_ad = autodiff) end """ @@ -43,27 +43,27 @@ Cache for the [`SwitchedEvolutionRelaxation`](@ref) method. internalnorm end -function requires_normal_form_jacobian(cache::Union{SwitchedEvolutionRelaxation, - SwitchedEvolutionRelaxationCache}) +function requires_normal_form_jacobian(cache::Union{ + SwitchedEvolutionRelaxation, SwitchedEvolutionRelaxationCache}) return false end -function requires_normal_form_rhs(cache::Union{SwitchedEvolutionRelaxation, - SwitchedEvolutionRelaxationCache}) +function requires_normal_form_rhs(cache::Union{ + SwitchedEvolutionRelaxation, SwitchedEvolutionRelaxationCache}) return false end -function __internal_init(prob::AbstractNonlinearProblem, f::SwitchedEvolutionRelaxation, - initial_damping, J, fu, u, args...; internalnorm::F = DEFAULT_NORM, - kwargs...) where {F} +function __internal_init( + prob::AbstractNonlinearProblem, f::SwitchedEvolutionRelaxation, initial_damping, + J, fu, u, args...; internalnorm::F = DEFAULT_NORM, kwargs...) where {F} T = promote_type(eltype(u), eltype(fu)) - return SwitchedEvolutionRelaxationCache(internalnorm(fu), T(1 / initial_damping), - internalnorm) + return SwitchedEvolutionRelaxationCache( + internalnorm(fu), T(1 / initial_damping), internalnorm) end (damping::SwitchedEvolutionRelaxationCache)(::Nothing) = damping.α⁻¹ -function __internal_solve!(damping::SwitchedEvolutionRelaxationCache, J, fu, args...; - kwargs...) +function __internal_solve!( + damping::SwitchedEvolutionRelaxationCache, J, fu, args...; kwargs...) res_norm = damping.internalnorm(fu) damping.α⁻¹ *= res_norm / damping.res_norm damping.res_norm = res_norm diff --git a/src/algorithms/raphson.jl b/src/algorithms/raphson.jl index bc005f2b6..c6847f54c 100644 --- a/src/algorithms/raphson.jl +++ b/src/algorithms/raphson.jl @@ -9,6 +9,6 @@ for large-scale and numerically-difficult nonlinear systems. function NewtonRaphson(; concrete_jac = nothing, linsolve = nothing, linesearch = NoLineSearch(), precs = DEFAULT_PRECS, autodiff = nothing) descent = NewtonDescent(; linsolve, precs) - return GeneralizedFirstOrderAlgorithm(; concrete_jac, name = :NewtonRaphson, - linesearch, descent, jacobian_ad = autodiff) + return GeneralizedFirstOrderAlgorithm(; + concrete_jac, name = :NewtonRaphson, linesearch, descent, jacobian_ad = autodiff) end diff --git a/src/algorithms/trust_region.jl b/src/algorithms/trust_region.jl index 7f550a65c..ac8f42d35 100644 --- a/src/algorithms/trust_region.jl +++ b/src/algorithms/trust_region.jl @@ -35,10 +35,10 @@ function TrustRegion(; concrete_jac = nothing, linsolve = nothing, precs = DEFAU if isnothing(vjp_autodiff) && autodiff isa ADTypes.AbstractFiniteDifferencesMode vjp_autodiff = autodiff end - trustregion = GenericTrustRegionScheme(; method = radius_update_scheme, step_threshold, - shrink_threshold, expand_threshold, shrink_factor, expand_factor, - reverse_ad = vjp_autodiff, forward_ad) - return GeneralizedFirstOrderAlgorithm(; concrete_jac, name = :TrustRegion, - trustregion, descent, jacobian_ad = autodiff, reverse_ad = vjp_autodiff, - max_shrink_times) + trustregion = GenericTrustRegionScheme(; + method = radius_update_scheme, step_threshold, shrink_threshold, expand_threshold, + shrink_factor, expand_factor, reverse_ad = vjp_autodiff, forward_ad) + return GeneralizedFirstOrderAlgorithm(; + concrete_jac, name = :TrustRegion, trustregion, descent, + jacobian_ad = autodiff, reverse_ad = vjp_autodiff, max_shrink_times) end diff --git a/src/core/approximate_jacobian.jl b/src/core/approximate_jacobian.jl index ffc77cea8..d0471e99b 100644 --- a/src/core/approximate_jacobian.jl +++ b/src/core/approximate_jacobian.jl @@ -50,14 +50,14 @@ function __show_algorithm(io::IO, alg::ApproximateJacobianSolveAlgorithm, name, print(io, "$(name)(\n$(spacing)$(join(modifiers, ",\n$(spacing)"))\n$(spacing_last))") end -function ApproximateJacobianSolveAlgorithm(; concrete_jac = nothing, - name::Symbol = :unknown, kwargs...) +function ApproximateJacobianSolveAlgorithm(; + concrete_jac = nothing, name::Symbol = :unknown, kwargs...) return ApproximateJacobianSolveAlgorithm{concrete_jac, name}(; kwargs...) end -function ApproximateJacobianSolveAlgorithm{concrete_jac, name}(; linesearch = missing, - trustregion = missing, descent, update_rule, reinit_rule, initialization, - max_resets::Int = typemax(Int), +function ApproximateJacobianSolveAlgorithm{concrete_jac, name}(; + linesearch = missing, trustregion = missing, descent, update_rule, + reinit_rule, initialization, max_resets::Int = typemax(Int), max_shrink_times::Int = typemax(Int)) where {concrete_jac, name} if linesearch !== missing && !(linesearch isa AbstractNonlinearSolveLineSearchAlgorithm) Base.depwarn("Passing in a `LineSearches.jl` algorithm directly is deprecated. \ @@ -65,8 +65,9 @@ function ApproximateJacobianSolveAlgorithm{concrete_jac, name}(; linesearch = mi :GeneralizedFirstOrderAlgorithm) linesearch = LineSearchesJL(; method = linesearch) end - return ApproximateJacobianSolveAlgorithm{concrete_jac, name}(linesearch, trustregion, - descent, update_rule, reinit_rule, max_resets, max_shrink_times, initialization) + return ApproximateJacobianSolveAlgorithm{concrete_jac, name}( + linesearch, trustregion, descent, update_rule, reinit_rule, + max_resets, max_shrink_times, initialization) end @inline concrete_jac(::ApproximateJacobianSolveAlgorithm{CJ}) where {CJ} = CJ @@ -117,9 +118,9 @@ end store_inverse_jacobian(::ApproximateJacobianSolveCache{INV}) where {INV} = INV -function __reinit_internal!(cache::ApproximateJacobianSolveCache{INV, GB, iip}, args...; - p = cache.p, u0 = cache.u, alias_u0::Bool = false, maxiters = 1000, - maxtime = nothing, kwargs...) where {INV, GB, iip} +function __reinit_internal!(cache::ApproximateJacobianSolveCache{INV, GB, iip}, + args...; p = cache.p, u0 = cache.u, alias_u0::Bool = false, + maxiters = 1000, maxtime = nothing, kwargs...) where {INV, GB, iip} if iip recursivecopy!(cache.u, u0) cache.prob.f(cache.fu, cache.u, p) @@ -147,10 +148,10 @@ end @internal_caches ApproximateJacobianSolveCache :initialization_cache :descent_cache :linesearch_cache :trustregion_cache :update_rule_cache :reinit_rule_cache -function SciMLBase.__init(prob::AbstractNonlinearProblem{uType, iip}, - alg::ApproximateJacobianSolveAlgorithm, args...; alias_u0 = false, - maxtime = nothing, maxiters = 1000, abstol = nothing, reltol = nothing, - linsolve_kwargs = (;), termination_condition = nothing, +function SciMLBase.__init( + prob::AbstractNonlinearProblem{uType, iip}, alg::ApproximateJacobianSolveAlgorithm, + args...; alias_u0 = false, maxtime = nothing, maxiters = 1000, abstol = nothing, + reltol = nothing, linsolve_kwargs = (;), termination_condition = nothing, internalnorm::F = DEFAULT_NORM, kwargs...) where {uType, iip, F} timer = get_timer_output() @static_timeit timer "cache construction" begin @@ -162,18 +163,18 @@ function SciMLBase.__init(prob::AbstractNonlinearProblem{uType, iip}, INV = store_inverse_jacobian(alg.update_rule) linsolve = get_linear_solver(alg.descent) - initialization_cache = __internal_init(prob, alg.initialization, alg, f, fu, u, p; - linsolve, - maxiters, internalnorm) + initialization_cache = __internal_init( + prob, alg.initialization, alg, f, fu, u, p; linsolve, maxiters, internalnorm) - abstol, reltol, termination_cache = init_termination_cache(abstol, reltol, fu, u, - termination_condition) + abstol, reltol, termination_cache = init_termination_cache( + abstol, reltol, fu, u, termination_condition) linsolve_kwargs = merge((; abstol, reltol), linsolve_kwargs) J = initialization_cache(nothing) inv_workspace, J = INV ? __safe_inv_workspace(J) : (nothing, J) - descent_cache = __internal_init(prob, alg.descent, J, fu, u; abstol, reltol, - internalnorm, linsolve_kwargs, pre_inverted = Val(INV), timer) + descent_cache = __internal_init( + prob, alg.descent, J, fu, u; abstol, reltol, internalnorm, + linsolve_kwargs, pre_inverted = Val(INV), timer) du = get_du(descent_cache) reinit_rule_cache = __internal_init(alg.reinit_rule, J, fu, u, du) @@ -189,31 +190,31 @@ function SciMLBase.__init(prob::AbstractNonlinearProblem{uType, iip}, if alg.trustregion !== missing supports_trust_region(alg.descent) || error("Trust Region not supported by \ $(alg.descent).") - trustregion_cache = __internal_init(prob, alg.trustregion, f, fu, u, p; - internalnorm, kwargs...) + trustregion_cache = __internal_init( + prob, alg.trustregion, f, fu, u, p; internalnorm, kwargs...) GB = :TrustRegion end if alg.linesearch !== missing supports_line_search(alg.descent) || error("Line Search not supported by \ $(alg.descent).") - linesearch_cache = __internal_init(prob, alg.linesearch, f, fu, u, p; - internalnorm, - kwargs...) + linesearch_cache = __internal_init( + prob, alg.linesearch, f, fu, u, p; internalnorm, kwargs...) GB = :LineSearch end - update_rule_cache = __internal_init(prob, alg.update_rule, J, fu, u, du; - internalnorm) + update_rule_cache = __internal_init( + prob, alg.update_rule, J, fu, u, du; internalnorm) trace = init_nonlinearsolve_trace(alg, u, fu, ApplyArray(__zero, J), du; uses_jacobian_inverse = Val(INV), kwargs...) - return ApproximateJacobianSolveCache{INV, GB, iip, maxtime !== nothing}(fu, u, - u_cache, p, du, J, alg, prob, initialization_cache, descent_cache, - linesearch_cache, trustregion_cache, update_rule_cache, reinit_rule_cache, - inv_workspace, 0, 0, 0, alg.max_resets, maxiters, maxtime, alg.max_shrink_times, - 0, timer, 0.0, termination_cache, trace, ReturnCode.Default, false, false) + return ApproximateJacobianSolveCache{INV, GB, iip, maxtime !== nothing}( + fu, u, u_cache, p, du, J, alg, prob, initialization_cache, + descent_cache, linesearch_cache, trustregion_cache, + update_rule_cache, reinit_rule_cache, inv_workspace, 0, 0, 0, + alg.max_resets, maxiters, maxtime, alg.max_shrink_times, 0, timer, + 0.0, termination_cache, trace, ReturnCode.Default, false, false) end end @@ -222,10 +223,8 @@ function __step!(cache::ApproximateJacobianSolveCache{INV, GB, iip}; new_jacobian = true @static_timeit cache.timer "jacobian init/reinit" begin if get_nsteps(cache) == 0 # First Step is special ignore kwargs - J_init = __internal_solve!(cache.initialization_cache, - cache.fu, - cache.u, - Val(false)) + J_init = __internal_solve!( + cache.initialization_cache, cache.fu, cache.u, Val(false)) if INV if jacobian_initialized_preinverted(cache.initialization_cache.alg) cache.J = J_init @@ -248,8 +247,8 @@ function __step!(cache::ApproximateJacobianSolveCache{INV, GB, iip}; cache.force_reinit = false elseif recompute_jacobian === nothing # Standard Step - reinit = __internal_solve!(cache.reinit_rule_cache, cache.J, cache.fu, - cache.u, cache.du) + reinit = __internal_solve!( + cache.reinit_rule_cache, cache.J, cache.fu, cache.u, cache.du) reinit && (countable_reinit = true) elseif recompute_jacobian reinit = true # Force ReInitialization: Don't count towards resetting @@ -268,8 +267,8 @@ function __step!(cache::ApproximateJacobianSolveCache{INV, GB, iip}; end if reinit - J_init = __internal_solve!(cache.initialization_cache, cache.fu, cache.u, - Val(true)) + J_init = __internal_solve!( + cache.initialization_cache, cache.fu, cache.u, Val(true)) cache.J = INV ? __safe_inv!!(cache.inv_workspace, J_init) : J_init J = cache.J cache.steps_since_last_reset = 0 @@ -284,13 +283,11 @@ function __step!(cache::ApproximateJacobianSolveCache{INV, GB, iip}; if cache.trustregion_cache !== nothing && hasfield(typeof(cache.trustregion_cache), :trust_region) δu, descent_success, descent_intermediates = __internal_solve!( - cache.descent_cache, - J, cache.fu, cache.u; new_jacobian, + cache.descent_cache, J, cache.fu, cache.u; new_jacobian, trust_region = cache.trustregion_cache.trust_region) else δu, descent_success, descent_intermediates = __internal_solve!( - cache.descent_cache, - J, cache.fu, cache.u; new_jacobian) + cache.descent_cache, J, cache.fu, cache.u; new_jacobian) end end @@ -309,8 +306,9 @@ function __step!(cache::ApproximateJacobianSolveCache{INV, GB, iip}; end elseif GB === :TrustRegion @static_timeit cache.timer "trustregion" begin - tr_accepted, u_new, fu_new = __internal_solve!(cache.trustregion_cache, J, - cache.fu, cache.u, δu, descent_intermediates) + tr_accepted, u_new, fu_new = __internal_solve!( + cache.trustregion_cache, J, cache.fu, + cache.u, δu, descent_intermediates) if tr_accepted @bb copyto!(cache.u, u_new) @bb copyto!(cache.fu, fu_new) @@ -341,7 +339,8 @@ function __step!(cache::ApproximateJacobianSolveCache{INV, GB, iip}; update_trace!(cache, α) @bb copyto!(cache.u_cache, cache.u) - if (cache.force_stop || cache.force_reinit || + if (cache.force_stop || + cache.force_reinit || (recompute_jacobian !== nothing && !recompute_jacobian)) callback_into_cache!(cache) return nothing diff --git a/src/core/generalized_first_order.jl b/src/core/generalized_first_order.jl index 0812e7f05..753b7682f 100644 --- a/src/core/generalized_first_order.jl +++ b/src/core/generalized_first_order.jl @@ -46,14 +46,14 @@ function __show_algorithm(io::IO, alg::GeneralizedFirstOrderAlgorithm, name, ind print(io, "$(name)(\n$(spacing)$(join(modifiers, ",\n$(spacing)"))\n$(spacing_last))") end -function GeneralizedFirstOrderAlgorithm(; concrete_jac = nothing, - name::Symbol = :unknown, kwargs...) +function GeneralizedFirstOrderAlgorithm(; + concrete_jac = nothing, name::Symbol = :unknown, kwargs...) return GeneralizedFirstOrderAlgorithm{concrete_jac, name}(; kwargs...) end -function GeneralizedFirstOrderAlgorithm{concrete_jac, name}(; descent, - linesearch = missing, trustregion = missing, jacobian_ad = nothing, - forward_ad = nothing, reverse_ad = nothing, +function GeneralizedFirstOrderAlgorithm{concrete_jac, name}(; + descent, linesearch = missing, trustregion = missing, + jacobian_ad = nothing, forward_ad = nothing, reverse_ad = nothing, max_shrink_times::Int = typemax(Int)) where {concrete_jac, name} forward_ad = ifelse(forward_ad !== nothing, forward_ad, ifelse(jacobian_ad isa ADTypes.AbstractForwardMode, jacobian_ad, nothing)) @@ -67,8 +67,9 @@ function GeneralizedFirstOrderAlgorithm{concrete_jac, name}(; descent, linesearch = LineSearchesJL(; method = linesearch) end - return GeneralizedFirstOrderAlgorithm{concrete_jac, name}(linesearch, - trustregion, descent, max_shrink_times, jacobian_ad, forward_ad, reverse_ad) + return GeneralizedFirstOrderAlgorithm{concrete_jac, name}( + linesearch, trustregion, descent, max_shrink_times, + jacobian_ad, forward_ad, reverse_ad) end concrete_jac(::GeneralizedFirstOrderAlgorithm{CJ}) where {CJ} = CJ @@ -112,9 +113,9 @@ concrete_jac(::GeneralizedFirstOrderAlgorithm{CJ}) where {CJ} = CJ force_stop::Bool end -function __reinit_internal!(cache::GeneralizedFirstOrderAlgorithmCache{iip}, args...; - p = cache.p, u0 = cache.u, alias_u0::Bool = false, maxiters = 1000, - maxtime = nothing, kwargs...) where {iip} +function __reinit_internal!( + cache::GeneralizedFirstOrderAlgorithmCache{iip}, args...; p = cache.p, u0 = cache.u, + alias_u0::Bool = false, maxiters = 1000, maxtime = nothing, kwargs...) where {iip} if iip recursivecopy!(cache.u, u0) cache.prob.f(cache.fu, cache.u, p) @@ -140,11 +141,11 @@ end @internal_caches GeneralizedFirstOrderAlgorithmCache :jac_cache :descent_cache :linesearch_cache :trustregion_cache -function SciMLBase.__init(prob::AbstractNonlinearProblem{uType, iip}, - alg::GeneralizedFirstOrderAlgorithm, args...; alias_u0 = false, maxiters = 1000, - abstol = nothing, reltol = nothing, maxtime = nothing, - termination_condition = nothing, internalnorm = DEFAULT_NORM, linsolve_kwargs = (;), - kwargs...) where {uType, iip} +function SciMLBase.__init( + prob::AbstractNonlinearProblem{uType, iip}, alg::GeneralizedFirstOrderAlgorithm, + args...; alias_u0 = false, maxiters = 1000, abstol = nothing, + reltol = nothing, maxtime = nothing, termination_condition = nothing, + internalnorm = DEFAULT_NORM, linsolve_kwargs = (;), kwargs...) where {uType, iip} timer = get_timer_output() @static_timeit timer "cache construction" begin (; f, u0, p) = prob @@ -154,12 +155,13 @@ function SciMLBase.__init(prob::AbstractNonlinearProblem{uType, iip}, linsolve = get_linear_solver(alg.descent) - abstol, reltol, termination_cache = init_termination_cache(abstol, reltol, fu, u, - termination_condition) + abstol, reltol, termination_cache = init_termination_cache( + abstol, reltol, fu, u, termination_condition) linsolve_kwargs = merge((; abstol, reltol), linsolve_kwargs) - jac_cache = JacobianCache(prob, alg, f, fu, u, p; autodiff = alg.jacobian_ad, - linsolve, jvp_autodiff = alg.forward_ad, vjp_autodiff = alg.reverse_ad) + jac_cache = JacobianCache( + prob, alg, f, fu, u, p; autodiff = alg.jacobian_ad, linsolve, + jvp_autodiff = alg.forward_ad, vjp_autodiff = alg.reverse_ad) J = jac_cache(nothing) descent_cache = __internal_init(prob, alg.descent, J, fu, u; abstol, reltol, internalnorm, linsolve_kwargs, timer) @@ -176,27 +178,25 @@ function SciMLBase.__init(prob::AbstractNonlinearProblem{uType, iip}, if alg.trustregion !== missing supports_trust_region(alg.descent) || error("Trust Region not supported by \ $(alg.descent).") - trustregion_cache = __internal_init(prob, alg.trustregion, f, fu, u, p; - internalnorm, - kwargs...) + trustregion_cache = __internal_init( + prob, alg.trustregion, f, fu, u, p; internalnorm, kwargs...) GB = :TrustRegion end if alg.linesearch !== missing supports_line_search(alg.descent) || error("Line Search not supported by \ $(alg.descent).") - linesearch_cache = __internal_init(prob, alg.linesearch, f, fu, u, p; - internalnorm, - kwargs...) + linesearch_cache = __internal_init( + prob, alg.linesearch, f, fu, u, p; internalnorm, kwargs...) GB = :LineSearch end trace = init_nonlinearsolve_trace(alg, u, fu, ApplyArray(__zero, J), du; kwargs...) - return GeneralizedFirstOrderAlgorithmCache{iip, GB, maxtime !== nothing}(fu, u, - u_cache, p, du, J, alg, prob, jac_cache, descent_cache, linesearch_cache, - trustregion_cache, 0, 0, maxiters, maxtime, alg.max_shrink_times, timer, 0.0, - true, termination_cache, trace, ReturnCode.Default, false) + return GeneralizedFirstOrderAlgorithmCache{iip, GB, maxtime !== nothing}( + fu, u, u_cache, p, du, J, alg, prob, jac_cache, descent_cache, linesearch_cache, + trustregion_cache, 0, 0, maxiters, maxtime, alg.max_shrink_times, + timer, 0.0, true, termination_cache, trace, ReturnCode.Default, false) end end @@ -216,13 +216,11 @@ function __step!(cache::GeneralizedFirstOrderAlgorithmCache{iip, GB}; if cache.trustregion_cache !== nothing && hasfield(typeof(cache.trustregion_cache), :trust_region) δu, descent_success, descent_intermediates = __internal_solve!( - cache.descent_cache, - J, cache.fu, cache.u; new_jacobian, + cache.descent_cache, J, cache.fu, cache.u; new_jacobian, trust_region = cache.trustregion_cache.trust_region) else δu, descent_success, descent_intermediates = __internal_solve!( - cache.descent_cache, - J, cache.fu, cache.u; new_jacobian) + cache.descent_cache, J, cache.fu, cache.u; new_jacobian) end end @@ -230,8 +228,8 @@ function __step!(cache::GeneralizedFirstOrderAlgorithmCache{iip, GB}; cache.make_new_jacobian = true if GB === :LineSearch @static_timeit cache.timer "linesearch" begin - linesearch_failed, α = __internal_solve!(cache.linesearch_cache, - cache.u, δu) + linesearch_failed, α = __internal_solve!( + cache.linesearch_cache, cache.u, δu) end if linesearch_failed cache.retcode = ReturnCode.InternalLineSearchFailed @@ -243,8 +241,9 @@ function __step!(cache::GeneralizedFirstOrderAlgorithmCache{iip, GB}; end elseif GB === :TrustRegion @static_timeit cache.timer "trustregion" begin - tr_accepted, u_new, fu_new = __internal_solve!(cache.trustregion_cache, J, - cache.fu, cache.u, δu, descent_intermediates) + tr_accepted, u_new, fu_new = __internal_solve!( + cache.trustregion_cache, J, cache.fu, + cache.u, δu, descent_intermediates) if tr_accepted @bb copyto!(cache.u, u_new) @bb copyto!(cache.fu, fu_new) diff --git a/src/core/generic.jl b/src/core/generic.jl index 849a259f1..22fc3e9d0 100644 --- a/src/core/generic.jl +++ b/src/core/generic.jl @@ -15,20 +15,20 @@ function SciMLBase.solve!(cache::AbstractNonlinearSolveCache) # The solver might have set a different `retcode` if cache.retcode == ReturnCode.Default - cache.retcode = ifelse(get_nsteps(cache) ≥ cache.maxiters, ReturnCode.MaxIters, - ReturnCode.Success) + cache.retcode = ifelse( + get_nsteps(cache) ≥ cache.maxiters, ReturnCode.MaxIters, ReturnCode.Success) end update_from_termination_cache!(cache.termination_cache, cache) - update_trace!(cache.trace, get_nsteps(cache), get_u(cache), get_fu(cache), nothing, - nothing, nothing; last = True) + update_trace!(cache.trace, get_nsteps(cache), get_u(cache), + get_fu(cache), nothing, nothing, nothing; last = True) stats = ImmutableNLStats(get_nf(cache), get_njacs(cache), get_nfactors(cache), get_nsolve(cache), get_nsteps(cache)) - return SciMLBase.build_solution(cache.prob, cache.alg, get_u(cache), get_fu(cache); - cache.retcode, stats, cache.trace) + return SciMLBase.build_solution(cache.prob, cache.alg, get_u(cache), + get_fu(cache); cache.retcode, stats, cache.trace) end """ @@ -45,8 +45,8 @@ Performs one step of the nonlinear solver. respectively. For algorithms that don't use jacobian information, this keyword is ignored with a one-time warning. """ -function SciMLBase.step!(cache::AbstractNonlinearSolveCache{iip, timeit}, args...; - kwargs...) where {iip, timeit} +function SciMLBase.step!(cache::AbstractNonlinearSolveCache{iip, timeit}, + args...; kwargs...) where {iip, timeit} timeit && (time_start = time()) res = @static_timeit cache.timer "solve" begin __step!(cache, args...; kwargs...) @@ -55,7 +55,8 @@ function SciMLBase.step!(cache::AbstractNonlinearSolveCache{iip, timeit}, args.. if timeit cache.total_time += time() - time_start - if !cache.force_stop && cache.retcode == ReturnCode.Default && + if !cache.force_stop && + cache.retcode == ReturnCode.Default && cache.total_time ≥ cache.maxtime cache.retcode = ReturnCode.MaxTime cache.force_stop = true diff --git a/src/core/spectral_methods.jl b/src/core/spectral_methods.jl index deab95b3f..a4e824744 100644 --- a/src/core/spectral_methods.jl +++ b/src/core/spectral_methods.jl @@ -76,9 +76,9 @@ concrete_jac(::GeneralizedDFSane) = nothing force_stop::Bool end -function __reinit_internal!(cache::GeneralizedDFSaneCache{iip}, args...; p = cache.p, - u0 = cache.u, alias_u0::Bool = false, maxiters = 1000, maxtime = nothing, - kwargs...) where {iip} +function __reinit_internal!( + cache::GeneralizedDFSaneCache{iip}, args...; p = cache.p, u0 = cache.u, + alias_u0::Bool = false, maxiters = 1000, maxtime = nothing, kwargs...) where {iip} if iip recursivecopy!(cache.u, u0) cache.prob.f(cache.fu, cache.u, p) @@ -115,10 +115,10 @@ end @internal_caches GeneralizedDFSaneCache :linesearch_cache -function SciMLBase.__init(prob::AbstractNonlinearProblem, alg::GeneralizedDFSane, args...; - alias_u0 = false, maxiters = 1000, abstol = nothing, reltol = nothing, - termination_condition = nothing, internalnorm::F = DEFAULT_NORM, maxtime = nothing, - kwargs...) where {F} +function SciMLBase.__init(prob::AbstractNonlinearProblem, alg::GeneralizedDFSane, + args...; alias_u0 = false, maxiters = 1000, abstol = nothing, + reltol = nothing, termination_condition = nothing, + internalnorm::F = DEFAULT_NORM, maxtime = nothing, kwargs...) where {F} timer = get_timer_output() @static_timeit timer "cache construction" begin u = __maybe_unaliased(prob.u0, alias_u0) @@ -129,11 +129,11 @@ function SciMLBase.__init(prob::AbstractNonlinearProblem, alg::GeneralizedDFSane fu = evaluate_f(prob, u) @bb fu_cache = copy(fu) - linesearch_cache = __internal_init(prob, alg.linesearch, prob.f, fu, u, prob.p; - maxiters, internalnorm, kwargs...) + linesearch_cache = __internal_init( + prob, alg.linesearch, prob.f, fu, u, prob.p; maxiters, internalnorm, kwargs...) - abstol, reltol, tc_cache = init_termination_cache(abstol, reltol, fu, u_cache, - termination_condition) + abstol, reltol, tc_cache = init_termination_cache( + abstol, reltol, fu, u_cache, termination_condition) trace = init_nonlinearsolve_trace(alg, u, fu, nothing, du; kwargs...) if alg.σ_1 === nothing @@ -148,10 +148,9 @@ function SciMLBase.__init(prob::AbstractNonlinearProblem, alg::GeneralizedDFSane end return GeneralizedDFSaneCache{isinplace(prob), maxtime !== nothing}( - fu, fu_cache, u, - u_cache, prob.p, du, alg, prob, σ_n, T(alg.σ_min), T(alg.σ_max), - linesearch_cache, 0, 0, maxiters, maxtime, timer, 0.0, tc_cache, trace, - ReturnCode.Default, false) + fu, fu_cache, u, u_cache, prob.p, du, alg, prob, σ_n, T(alg.σ_min), + T(alg.σ_max), linesearch_cache, 0, 0, maxiters, maxtime, + timer, 0.0, tc_cache, trace, ReturnCode.Default, false) end end diff --git a/src/default.jl b/src/default.jl index dafdf7902..8b17a9393 100644 --- a/src/default.jl +++ b/src/default.jl @@ -59,24 +59,21 @@ end for (probType, pType) in ((:NonlinearProblem, :NLS), (:NonlinearLeastSquaresProblem, :NLLS)) algType = NonlinearSolvePolyAlgorithm{pType} @eval begin - function SciMLBase.__init(prob::$probType, alg::$algType{N}, args...; - kwargs...) where {N} + function SciMLBase.__init( + prob::$probType, alg::$algType{N}, args...; kwargs...) where {N} return NonlinearSolvePolyAlgorithmCache{isinplace(prob), N}( - map(solver -> SciMLBase.__init(prob, - solver, args...; kwargs...), alg.algs), + map(solver -> SciMLBase.__init(prob, solver, args...; kwargs...), alg.algs), alg, 1) end end end -@generated function SciMLBase.solve!(cache::NonlinearSolvePolyAlgorithmCache{iip, - N}) where {iip, N} - calls = [ - quote +@generated function SciMLBase.solve!(cache::NonlinearSolvePolyAlgorithmCache{ + iip, N}) where {iip, N} + calls = [quote 1 ≤ cache.current ≤ length(cache.caches) || error("Current choices shouldn't get here!") - end - ] + end] cache_syms = [gensym("cache") for i in 1:N] sol_syms = [gensym("sol") for i in 1:N] @@ -90,8 +87,9 @@ end stats = $(sol_syms[i]).stats u = $(sol_syms[i]).u fu = get_fu($(cache_syms[i])) - return SciMLBase.build_solution($(sol_syms[i]).prob, cache.alg, u, - fu; retcode = ReturnCode.Success, stats, + return SciMLBase.build_solution( + $(sol_syms[i]).prob, cache.alg, u, fu; + retcode = ReturnCode.Success, stats, original = $(sol_syms[i]), trace = $(sol_syms[i]).trace) end cache.current = $(i + 1) @@ -112,8 +110,8 @@ end stats = cache.caches[idx].stats u = cache.caches[idx].u - return SciMLBase.build_solution(cache.caches[idx].prob, cache.alg, u, - fus[idx]; retcode, stats, cache.caches[idx].trace) + return SciMLBase.build_solution(cache.caches[idx].prob, cache.alg, u, fus[idx]; + retcode, stats, cache.caches[idx].trace) end) return Expr(:block, calls...) @@ -122,19 +120,19 @@ end for (probType, pType) in ((:NonlinearProblem, :NLS), (:NonlinearLeastSquaresProblem, :NLLS)) algType = NonlinearSolvePolyAlgorithm{pType} @eval begin - @generated function SciMLBase.__solve(prob::$probType, alg::$algType{N}, args...; - kwargs...) where {N} + @generated function SciMLBase.__solve( + prob::$probType, alg::$algType{N}, args...; kwargs...) where {N} calls = [] sol_syms = [gensym("sol") for _ in 1:N] for i in 1:N cur_sol = sol_syms[i] push!(calls, quote - $(cur_sol) = SciMLBase.__solve(prob, alg.algs[$(i)], args...; - kwargs...) + $(cur_sol) = SciMLBase.__solve(prob, alg.algs[$(i)], args...; kwargs...) if SciMLBase.successful_retcode($(cur_sol)) - return SciMLBase.build_solution(prob, alg, $(cur_sol).u, - $(cur_sol).resid; $(cur_sol).retcode, $(cur_sol).stats, + return SciMLBase.build_solution( + prob, alg, $(cur_sol).u, $(cur_sol).resid; + $(cur_sol).retcode, $(cur_sol).stats, original = $(cur_sol), trace = $(cur_sol).trace) end end) @@ -145,11 +143,10 @@ for (probType, pType) in ((:NonlinearProblem, :NLS), (:NonlinearLeastSquaresProb push!(calls, :($(resid) = $(sym).resid)) end - push!(calls, - quote - resids = tuple($(Tuple(resids)...)) - minfu, idx = __findmin(DEFAULT_NORM, resids) - end) + push!(calls, quote + resids = tuple($(Tuple(resids)...)) + minfu, idx = __findmin(DEFAULT_NORM, resids) + end) for i in 1:N push!(calls, @@ -217,8 +214,9 @@ for more performance and then tries more robust techniques if the faster ones fa - `T`: The eltype of the initial guess. It is only used to check if some of the algorithms are compatible with the problem type. Defaults to `Float64`. """ -function FastShortcutNonlinearPolyalg(::Type{T} = Float64; concrete_jac = nothing, - linsolve = nothing, precs = DEFAULT_PRECS, must_use_jacobian::Val{JAC} = Val(false), +function FastShortcutNonlinearPolyalg( + ::Type{T} = Float64; concrete_jac = nothing, linsolve = nothing, + precs = DEFAULT_PRECS, must_use_jacobian::Val{JAC} = Val(false), prefer_simplenonlinearsolve::Val{SA} = Val(false), autodiff = nothing) where {T, JAC, SA} if JAC @@ -237,8 +235,7 @@ function FastShortcutNonlinearPolyalg(::Type{T} = Float64; concrete_jac = nothin # and thus are not included in the polyalgorithm if SA if __is_complex(T) - algs = (SimpleBroyden(), - Broyden(; init_jacobian = Val(:true_jacobian)), + algs = (SimpleBroyden(), Broyden(; init_jacobian = Val(:true_jacobian)), SimpleKlement(), NewtonRaphson(; concrete_jac, linsolve, precs, autodiff)) else @@ -253,8 +250,7 @@ function FastShortcutNonlinearPolyalg(::Type{T} = Float64; concrete_jac = nothin end else if __is_complex(T) - algs = (Broyden(), - Broyden(; init_jacobian = Val(:true_jacobian)), + algs = (Broyden(), Broyden(; init_jacobian = Val(:true_jacobian)), Klement(; linsolve, precs), NewtonRaphson(; concrete_jac, linsolve, precs, autodiff)) else @@ -317,8 +313,8 @@ end ## can use that! function SciMLBase.__init(prob::NonlinearProblem, ::Nothing, args...; kwargs...) must_use_jacobian = Val(prob.f.jac !== nothing) - return SciMLBase.__init(prob, - FastShortcutNonlinearPolyalg(eltype(prob.u0); must_use_jacobian), + return SciMLBase.__init( + prob, FastShortcutNonlinearPolyalg(eltype(prob.u0); must_use_jacobian), args...; kwargs...) end @@ -326,17 +322,19 @@ function SciMLBase.__solve(prob::NonlinearProblem, ::Nothing, args...; kwargs... must_use_jacobian = Val(prob.f.jac !== nothing) prefer_simplenonlinearsolve = Val(prob.u0 isa SArray) return SciMLBase.__solve(prob, - FastShortcutNonlinearPolyalg(eltype(prob.u0); must_use_jacobian, - prefer_simplenonlinearsolve), args...; kwargs...) + FastShortcutNonlinearPolyalg( + eltype(prob.u0); must_use_jacobian, prefer_simplenonlinearsolve), + args...; + kwargs...) end function SciMLBase.__init(prob::NonlinearLeastSquaresProblem, ::Nothing, args...; kwargs...) - return SciMLBase.__init(prob, FastShortcutNLLSPolyalg(eltype(prob.u0)), args...; - kwargs...) + return SciMLBase.__init( + prob, FastShortcutNLLSPolyalg(eltype(prob.u0)), args...; kwargs...) end -function SciMLBase.__solve(prob::NonlinearLeastSquaresProblem, ::Nothing, args...; - kwargs...) - return SciMLBase.__solve(prob, FastShortcutNLLSPolyalg(eltype(prob.u0)), args...; - kwargs...) +function SciMLBase.__solve( + prob::NonlinearLeastSquaresProblem, ::Nothing, args...; kwargs...) + return SciMLBase.__solve( + prob, FastShortcutNLLSPolyalg(eltype(prob.u0)), args...; kwargs...) end diff --git a/src/descent/damped_newton.jl b/src/descent/damped_newton.jl index 77ad95b54..3e21cc2d1 100644 --- a/src/descent/damped_newton.jl +++ b/src/descent/damped_newton.jl @@ -51,11 +51,10 @@ end @internal_caches DampedNewtonDescentCache :lincache :damping_fn_cache -function __internal_init( - prob::AbstractNonlinearProblem, alg::DampedNewtonDescent, J, fu, u; - pre_inverted::Val{INV} = False, linsolve_kwargs = (;), abstol = nothing, - timer = get_timer_output(), reltol = nothing, alias_J = true, - shared::Val{N} = Val(1), kwargs...) where {INV, N} +function __internal_init(prob::AbstractNonlinearProblem, alg::DampedNewtonDescent, J, + fu, u; pre_inverted::Val{INV} = False, linsolve_kwargs = (;), + abstol = nothing, timer = get_timer_output(), reltol = nothing, + alias_J = true, shared::Val{N} = Val(1), kwargs...) where {INV, N} length(fu) != length(u) && @assert !INV "Precomputed Inverse for Non-Square Jacobian doesn't make sense." @bb δu = similar(u) @@ -103,8 +102,7 @@ function __internal_init( A, b = J_cache, rhs_cache elseif mode === :simple damping_fn_cache = __internal_init( - prob, alg.damping_fn, alg.initial_damping, J, fu, - u, False; kwargs...) + prob, alg.damping_fn, alg.initial_damping, J, fu, u, False; kwargs...) J_cache = __maybe_unaliased(J, alias_J) D = damping_fn_cache(nothing) J_damped = __dampen_jacobian!!(J_cache, J, D) @@ -117,8 +115,7 @@ function __internal_init( jac_damp = requires_normal_form_jacobian(alg.damping_fn) ? JᵀJ : J rhs_damp = requires_normal_form_rhs(alg.damping_fn) ? Jᵀfu : fu damping_fn_cache = __internal_init(prob, alg.damping_fn, alg.initial_damping, - jac_damp, - rhs_damp, u, True; kwargs...) + jac_damp, rhs_damp, u, True; kwargs...) D = damping_fn_cache(nothing) @bb J_cache = similar(JᵀJ) @bb @. J_cache = 0 @@ -127,16 +124,16 @@ function __internal_init( rhs_cache = nothing end - lincache = LinearSolverCache(alg, alg.linsolve, A, b, _vec(u); abstol, reltol, - linsolve_kwargs...) + lincache = LinearSolverCache( + alg, alg.linsolve, A, b, _vec(u); abstol, reltol, linsolve_kwargs...) - return DampedNewtonDescentCache{INV, mode}(J_cache, δu, δus, lincache, JᵀJ, Jᵀfu, - rhs_cache, damping_fn_cache, timer) + return DampedNewtonDescentCache{INV, mode}( + J_cache, δu, δus, lincache, JᵀJ, Jᵀfu, rhs_cache, damping_fn_cache, timer) end -function __internal_solve!(cache::DampedNewtonDescentCache{INV, mode}, J, fu, u, - idx::Val{N} = Val(1); skip_solve::Bool = false, new_jacobian::Bool = true, - kwargs...) where {INV, N, mode} +function __internal_solve!(cache::DampedNewtonDescentCache{INV, mode}, J, fu, + u, idx::Val{N} = Val(1); skip_solve::Bool = false, + new_jacobian::Bool = true, kwargs...) where {INV, N, mode} δu = get_du(cache, idx) skip_solve && return δu, true, (;) @@ -186,8 +183,8 @@ function __internal_solve!(cache::DampedNewtonDescentCache{INV, mode}, J, fu, u, INV && (J = inv(J)) @bb cache.JᵀJ_cache = transpose(J) × J @bb cache.Jᵀfu_cache = transpose(J) × vec(fu) - D = __internal_solve!(cache.damping_fn_cache, cache.JᵀJ_cache, - cache.Jᵀfu_cache, True) + D = __internal_solve!( + cache.damping_fn_cache, cache.JᵀJ_cache, cache.Jᵀfu_cache, True) cache.J = __dampen_jacobian!!(cache.J, cache.JᵀJ_cache, D) A = __maybe_symmetric(cache.J) elseif !recompute_A @@ -203,8 +200,8 @@ function __internal_solve!(cache::DampedNewtonDescentCache{INV, mode}, J, fu, u, end @static_timeit cache.timer "linear solve" begin - δu = cache.lincache(; A, b, - reuse_A_if_factorization = !new_jacobian && !recompute_A, + δu = cache.lincache(; + A, b, reuse_A_if_factorization = !new_jacobian && !recompute_A, kwargs..., linu = _vec(δu)) δu = _restructure(get_du(cache, idx), δu) end diff --git a/src/descent/dogleg.jl b/src/descent/dogleg.jl index e1a50832f..ee725988b 100644 --- a/src/descent/dogleg.jl +++ b/src/descent/dogleg.jl @@ -33,8 +33,7 @@ function Dogleg(; linsolve = nothing, precs = DEFAULT_PRECS, damping = False, SteepestDescent(; linsolve, precs)) end -@concrete mutable struct DoglegCache{pre_inverted, normalform} <: - AbstractDescentCache +@concrete mutable struct DoglegCache{pre_inverted, normalform} <: AbstractDescentCache δu δus newton_cache @@ -49,9 +48,9 @@ end @internal_caches DoglegCache :newton_cache :cauchy_cache function __internal_init(prob::AbstractNonlinearProblem, alg::Dogleg, J, fu, u; - pre_inverted::Val{INV} = False, linsolve_kwargs = (;), abstol = nothing, - reltol = nothing, internalnorm::F = DEFAULT_NORM, shared::Val{N} = Val(1), - kwargs...) where {F, INV, N} + pre_inverted::Val{INV} = False, linsolve_kwargs = (;), + abstol = nothing, reltol = nothing, internalnorm::F = DEFAULT_NORM, + shared::Val{N} = Val(1), kwargs...) where {F, INV, N} newton_cache = __internal_init(prob, alg.newton_descent, J, fu, u; pre_inverted, linsolve_kwargs, abstol, reltol, shared, kwargs...) cauchy_cache = __internal_init(prob, alg.steepest_descent, J, fu, u; pre_inverted, @@ -82,8 +81,8 @@ function __internal_solve!(cache::DoglegCache{INV, NF}, J, fu, u, idx::Val{N} = want to use a Trust Region." δu = get_du(cache, idx) T = promote_type(eltype(u), eltype(fu)) - δu_newton, _, _ = __internal_solve!(cache.newton_cache, J, fu, u, idx; skip_solve, - kwargs...) + δu_newton, _, _ = __internal_solve!( + cache.newton_cache, J, fu, u, idx; skip_solve, kwargs...) # Newton's Step within the trust region if cache.internalnorm(δu_newton) ≤ trust_region @@ -103,8 +102,8 @@ function __internal_solve!(cache::DoglegCache{INV, NF}, J, fu, u, idx::Val{N} = @bb cache.δu_cache_mul = JᵀJ × vec(δu_cauchy) δuJᵀJδu = __dot(δu_cauchy, cache.δu_cache_mul) else - δu_cauchy, _, _ = __internal_solve!(cache.cauchy_cache, J, fu, u, idx; skip_solve, - kwargs...) + δu_cauchy, _, _ = __internal_solve!( + cache.cauchy_cache, J, fu, u, idx; skip_solve, kwargs...) J_ = INV ? inv(J) : J l_grad = cache.internalnorm(δu_cauchy) @bb cache.JᵀJ_cache = J × vec(δu_cauchy) # TODO: Rename diff --git a/src/descent/geodesic_acceleration.jl b/src/descent/geodesic_acceleration.jl index 35764783c..6cba1202e 100644 --- a/src/descent/geodesic_acceleration.jl +++ b/src/descent/geodesic_acceleration.jl @@ -31,8 +31,7 @@ for other methods are not theorectically or experimentally verified. end function Base.show(io::IO, alg::GeodesicAcceleration) - print( - io, "GeodesicAcceleration(descent = $(alg.descent), finite_diff_step_geodesic = ", + print(io, "GeodesicAcceleration(descent = $(alg.descent), finite_diff_step_geodesic = ", "$(alg.finite_diff_step_geodesic), α = $(alg.α))") end @@ -55,8 +54,8 @@ get_linear_solver(alg::GeodesicAcceleration) = get_linear_solver(alg.descent) last_step_accepted::Bool end -function __reinit_internal!(cache::GeodesicAccelerationCache, args...; p = cache.p, - kwargs...) +function __reinit_internal!( + cache::GeodesicAccelerationCache, args...; p = cache.p, kwargs...) cache.p = p cache.last_step_accepted = false end @@ -84,10 +83,10 @@ function set_acceleration!(cache::GeodesicAccelerationCache, δa, ::Val{N}) wher set_du!(cache.descent_cache, δa, Val(2N)) end -function __internal_init(prob::AbstractNonlinearProblem, alg::GeodesicAcceleration, J, fu, - u; shared::Val{N} = Val(1), pre_inverted::Val{INV} = False, linsolve_kwargs = (;), - abstol = nothing, reltol = nothing, internalnorm::F = DEFAULT_NORM, - kwargs...) where {INV, N, F} +function __internal_init(prob::AbstractNonlinearProblem, alg::GeodesicAcceleration, J, + fu, u; shared::Val{N} = Val(1), pre_inverted::Val{INV} = False, + linsolve_kwargs = (;), abstol = nothing, reltol = nothing, + internalnorm::F = DEFAULT_NORM, kwargs...) where {INV, N, F} T = promote_type(eltype(u), eltype(fu)) @bb δu = similar(u) δus = N ≤ 1 ? nothing : map(2:N) do i @@ -98,17 +97,17 @@ function __internal_init(prob::AbstractNonlinearProblem, alg::GeodesicAccelerati @bb Jv = similar(fu) @bb fu_cache = copy(fu) @bb u_cache = similar(u) - return GeodesicAccelerationCache(δu, δus, descent_cache, prob.f, prob.p, T(alg.α), - internalnorm, T(alg.finite_diff_step_geodesic), Jv, fu_cache, u_cache, false) + return GeodesicAccelerationCache( + δu, δus, descent_cache, prob.f, prob.p, T(alg.α), internalnorm, + T(alg.finite_diff_step_geodesic), Jv, fu_cache, u_cache, false) end -function __internal_solve!( - cache::GeodesicAccelerationCache, J, fu, u, idx::Val{N} = Val(1); +function __internal_solve!(cache::GeodesicAccelerationCache, J, fu, u, idx::Val{N} = Val(1); skip_solve::Bool = false, kwargs...) where {N} a, v, δu = get_acceleration(cache, idx), get_velocity(cache, idx), get_du(cache, idx) skip_solve && return δu, true, (; a, v) - v, _, _ = __internal_solve!(cache.descent_cache, J, fu, u, Val(2N - 1); skip_solve, - kwargs...) + v, _, _ = __internal_solve!( + cache.descent_cache, J, fu, u, Val(2N - 1); skip_solve, kwargs...) @bb @. cache.u_cache = u + cache.h * v cache.fu_cache = evaluate_f!!(cache.f, cache.fu_cache, cache.u_cache, cache.p) diff --git a/src/descent/newton.jl b/src/descent/newton.jl index c8ba35ed9..cfbc1c239 100644 --- a/src/descent/newton.jl +++ b/src/descent/newton.jl @@ -32,24 +32,24 @@ end @internal_caches NewtonDescentCache :lincache -function __internal_init(prob::NonlinearProblem, alg::NewtonDescent, J, fu, u; - shared::Val{N} = Val(1), pre_inverted::Val{INV} = False, linsolve_kwargs = (;), - abstol = nothing, reltol = nothing, timer = get_timer_output(), - kwargs...) where {INV, N} +function __internal_init( + prob::NonlinearProblem, alg::NewtonDescent, J, fu, u; shared::Val{N} = Val(1), + pre_inverted::Val{INV} = False, linsolve_kwargs = (;), abstol = nothing, + reltol = nothing, timer = get_timer_output(), kwargs...) where {INV, N} @bb δu = similar(u) δus = N ≤ 1 ? nothing : map(2:N) do i @bb δu_ = similar(u) end INV && return NewtonDescentCache{true, false}(δu, δus, nothing, nothing, nothing, timer) - lincache = LinearSolverCache(alg, alg.linsolve, J, _vec(fu), _vec(u); abstol, reltol, - linsolve_kwargs...) + lincache = LinearSolverCache( + alg, alg.linsolve, J, _vec(fu), _vec(u); abstol, reltol, linsolve_kwargs...) return NewtonDescentCache{false, false}(δu, δus, lincache, nothing, nothing, timer) end -function __internal_init(prob::NonlinearLeastSquaresProblem, alg::NewtonDescent, J, fu, u; - pre_inverted::Val{INV} = False, linsolve_kwargs = (;), shared::Val{N} = Val(1), - abstol = nothing, reltol = nothing, timer = get_timer_output(), - kwargs...) where {INV, N} +function __internal_init(prob::NonlinearLeastSquaresProblem, alg::NewtonDescent, J, + fu, u; pre_inverted::Val{INV} = False, linsolve_kwargs = (;), + shared::Val{N} = Val(1), abstol = nothing, reltol = nothing, + timer = get_timer_output(), kwargs...) where {INV, N} length(fu) != length(u) && @assert !INV "Precomputed Inverse for Non-Square Jacobian doesn't make sense." @@ -62,8 +62,8 @@ function __internal_init(prob::NonlinearLeastSquaresProblem, alg::NewtonDescent, JᵀJ, Jᵀfu = nothing, nothing A, b = J, _vec(fu) end - lincache = LinearSolverCache(alg, alg.linsolve, A, b, _vec(u); abstol, reltol, - linsolve_kwargs...) + lincache = LinearSolverCache( + alg, alg.linsolve, A, b, _vec(u); abstol, reltol, linsolve_kwargs...) @bb δu = similar(u) δus = N ≤ 1 ? nothing : map(2:N) do i @bb δu_ = similar(u) @@ -71,9 +71,9 @@ function __internal_init(prob::NonlinearLeastSquaresProblem, alg::NewtonDescent, return NewtonDescentCache{false, normal_form}(δu, δus, lincache, JᵀJ, Jᵀfu, timer) end -function __internal_solve!(cache::NewtonDescentCache{INV, false}, J, fu, u, - idx::Val = Val(1); skip_solve::Bool = false, new_jacobian::Bool = true, - kwargs...) where {INV} +function __internal_solve!( + cache::NewtonDescentCache{INV, false}, J, fu, u, idx::Val = Val(1); + skip_solve::Bool = false, new_jacobian::Bool = true, kwargs...) where {INV} δu = get_du(cache, idx) skip_solve && return δu, true, (;) if INV @@ -81,8 +81,9 @@ function __internal_solve!(cache::NewtonDescentCache{INV, false}, J, fu, u, @bb δu = J × vec(fu) else @static_timeit cache.timer "linear solve" begin - δu = cache.lincache(; A = J, b = _vec(fu), kwargs..., linu = _vec(δu), - du = _vec(δu), reuse_A_if_factorization = !new_jacobian || (idx !== Val(1))) + δu = cache.lincache(; + A = J, b = _vec(fu), kwargs..., linu = _vec(δu), du = _vec(δu), + reuse_A_if_factorization = !new_jacobian || (idx !== Val(1))) δu = _restructure(get_du(cache, idx), δu) end end @@ -91,8 +92,9 @@ function __internal_solve!(cache::NewtonDescentCache{INV, false}, J, fu, u, return δu, true, (;) end -function __internal_solve!(cache::NewtonDescentCache{false, true}, J, fu, u, - idx::Val = Val(1); skip_solve::Bool = false, new_jacobian::Bool = true, kwargs...) +function __internal_solve!( + cache::NewtonDescentCache{false, true}, J, fu, u, idx::Val = Val(1); + skip_solve::Bool = false, new_jacobian::Bool = true, kwargs...) δu = get_du(cache, idx) skip_solve && return δu, true, (;) if idx === Val(1) diff --git a/src/descent/steepest.jl b/src/descent/steepest.jl index d19505a86..40919cc6e 100644 --- a/src/descent/steepest.jl +++ b/src/descent/steepest.jl @@ -29,8 +29,9 @@ end @internal_caches SteepestDescentCache :lincache -@inline function __internal_init(prob::AbstractNonlinearProblem, alg::SteepestDescent, J, - fu, u; shared::Val{N} = Val(1), pre_inverted::Val{INV} = False, +@inline function __internal_init( + prob::AbstractNonlinearProblem, alg::SteepestDescent, J, fu, + u; shared::Val{N} = Val(1), pre_inverted::Val{INV} = False, linsolve_kwargs = (;), abstol = nothing, reltol = nothing, timer = get_timer_output(), kwargs...) where {INV, N} INV && @assert length(fu)==length(u) "Non-Square Jacobian Inverse doesn't make sense." @@ -39,8 +40,8 @@ end @bb δu_ = similar(u) end if INV - lincache = LinearSolverCache(alg, alg.linsolve, transpose(J), _vec(fu), _vec(u); - abstol, reltol, linsolve_kwargs...) + lincache = LinearSolverCache(alg, alg.linsolve, transpose(J), _vec(fu), + _vec(u); abstol, reltol, linsolve_kwargs...) else lincache = nothing end @@ -53,8 +54,9 @@ function __internal_solve!(cache::SteepestDescentCache{INV}, J, fu, u, idx::Val if INV A = J === nothing ? nothing : transpose(J) @static_timeit cache.timer "linear solve" begin - δu = cache.lincache(; A, b = _vec(fu), kwargs..., linu = _vec(δu), - du = _vec(δu), reuse_A_if_factorization = !new_jacobian || idx !== Val(1)) + δu = cache.lincache(; + A, b = _vec(fu), kwargs..., linu = _vec(δu), du = _vec(δu), + reuse_A_if_factorization = !new_jacobian || idx !== Val(1)) δu = _restructure(get_du(cache, idx), δu) end else diff --git a/src/globalization/line_search.jl b/src/globalization/line_search.jl index 6e13043de..9d4f52e71 100644 --- a/src/globalization/line_search.jl +++ b/src/globalization/line_search.jl @@ -9,8 +9,8 @@ struct NoLineSearch <: AbstractNonlinearSolveLineSearchAlgorithm end α end -function __internal_init(prob::AbstractNonlinearProblem, alg::NoLineSearch, f::F, fu, u, - p, args...; kwargs...) where {F} +function __internal_init(prob::AbstractNonlinearProblem, alg::NoLineSearch, + f::F, fu, u, p, args...; kwargs...) where {F} return NoLineSearchCache(promote_type(eltype(fu), eltype(u))(true)) end @@ -80,8 +80,9 @@ Base.@deprecate_binding LineSearch LineSearchesJL true nf::Base.RefValue{Int} end -function __internal_init(prob::AbstractNonlinearProblem, alg::LineSearchesJL, f::F, fu, u, - p, args...; internalnorm::IN = DEFAULT_NORM, kwargs...) where {F, IN} +function __internal_init( + prob::AbstractNonlinearProblem, alg::LineSearchesJL, f::F, fu, u, p, + args...; internalnorm::IN = DEFAULT_NORM, kwargs...) where {F, IN} T = promote_type(eltype(fu), eltype(u)) if u isa Number grad_op = @closure (u, fu, p) -> last(__value_derivative(Base.Fix2(f, p), u)) * fu @@ -94,8 +95,8 @@ function __internal_init(prob::AbstractNonlinearProblem, alg::LineSearchesJL, f: grad_op = @closure (u, fu, p) -> f.vjp(fu, u, p) end else - autodiff = get_concrete_reverse_ad(alg.autodiff, prob; - check_forward_mode = true) + autodiff = get_concrete_reverse_ad( + alg.autodiff, prob; check_forward_mode = true) vjp_op = VecJacOperator(prob, fu, u; autodiff) if isinplace(prob) g_cache = similar(u) @@ -134,17 +135,16 @@ function __internal_init(prob::AbstractNonlinearProblem, alg::LineSearchesJL, f: return obj, dot(g₀, du) end - return LineSearchesJLCache(f, p, ϕ, dϕ, ϕdϕ, alg.method, T(alg.initial_alpha), grad_op, - u_cache, fu_cache, nf) + return LineSearchesJLCache( + f, p, ϕ, dϕ, ϕdϕ, alg.method, T(alg.initial_alpha), grad_op, u_cache, fu_cache, nf) end function __internal_solve!(cache::LineSearchesJLCache, u, du; kwargs...) ϕ = @closure α -> cache.ϕ(cache.f, cache.p, u, du, α, cache.u_cache, cache.fu_cache) - dϕ = @closure α -> cache.dϕ(cache.f, cache.p, u, du, α, cache.u_cache, cache.fu_cache, - cache.grad_op) + dϕ = @closure α -> cache.dϕ( + cache.f, cache.p, u, du, α, cache.u_cache, cache.fu_cache, cache.grad_op) ϕdϕ = @closure α -> cache.ϕdϕ( - cache.f, cache.p, u, du, α, cache.u_cache, cache.fu_cache, - cache.grad_op) + cache.f, cache.p, u, du, α, cache.u_cache, cache.fu_cache, cache.grad_op) ϕ₀, dϕ₀ = ϕdϕ(zero(eltype(u))) @@ -225,8 +225,9 @@ end nf::Base.RefValue{Int} end -function __internal_init(prob::AbstractNonlinearProblem, alg::RobustNonMonotoneLineSearch, - f::F, fu, u, p, args...; internalnorm::IN = DEFAULT_NORM, kwargs...) where {F, IN} +function __internal_init( + prob::AbstractNonlinearProblem, alg::RobustNonMonotoneLineSearch, f::F, fu, + u, p, args...; internalnorm::IN = DEFAULT_NORM, kwargs...) where {F, IN} @bb u_cache = similar(u) @bb fu_cache = similar(fu) T = promote_type(eltype(fu), eltype(u)) @@ -242,9 +243,10 @@ function __internal_init(prob::AbstractNonlinearProblem, alg::RobustNonMonotoneL fn₁ = internalnorm(fu)^alg.n_exp η_strategy = @closure (n, xₙ, fₙ) -> alg.η_strategy(fn₁, n, xₙ, fₙ) - return RobustNonMonotoneLineSearchCache(f, p, ϕ, u_cache, fu_cache, internalnorm, - alg.maxiters, fill(fn₁, alg.M), T(alg.gamma), T(alg.sigma_1), alg.M, T(alg.tau_min), - T(alg.tau_max), 0, η_strategy, alg.n_exp, nf) + return RobustNonMonotoneLineSearchCache( + f, p, ϕ, u_cache, fu_cache, internalnorm, alg.maxiters, + fill(fn₁, alg.M), T(alg.gamma), T(alg.sigma_1), alg.M, + T(alg.tau_min), T(alg.tau_max), 0, η_strategy, alg.n_exp, nf) end function __internal_solve!(cache::RobustNonMonotoneLineSearchCache, u, du; kwargs...) @@ -316,8 +318,9 @@ end nf::Base.RefValue{Int} end -function __internal_init(prob::AbstractNonlinearProblem, alg::LiFukushimaLineSearch, - f::F, fu, u, p, args...; internalnorm::IN = DEFAULT_NORM, kwargs...) where {F, IN} +function __internal_init( + prob::AbstractNonlinearProblem, alg::LiFukushimaLineSearch, f::F, fu, u, + p, args...; internalnorm::IN = DEFAULT_NORM, kwargs...) where {F, IN} @bb u_cache = similar(u) @bb fu_cache = similar(fu) T = promote_type(eltype(fu), eltype(u)) @@ -330,8 +333,9 @@ function __internal_init(prob::AbstractNonlinearProblem, alg::LiFukushimaLineSea return internalnorm(fu_cache) end - return LiFukushimaLineSearchCache(ϕ, f, p, internalnorm, u_cache, fu_cache, - T(alg.lambda_0), T(alg.beta), T(alg.sigma_1), T(alg.sigma_2), T(alg.eta), + return LiFukushimaLineSearchCache( + ϕ, f, p, internalnorm, u_cache, fu_cache, T(alg.lambda_0), + T(alg.beta), T(alg.sigma_1), T(alg.sigma_2), T(alg.eta), T(alg.rho), T(true), alg.nan_max_iter, alg.maxiters, nf) end diff --git a/src/globalization/trust_region.jl b/src/globalization/trust_region.jl index 891767a4a..6004695ec 100644 --- a/src/globalization/trust_region.jl +++ b/src/globalization/trust_region.jl @@ -46,8 +46,8 @@ end nf::Int end -function reinit_cache!(cache::LevenbergMarquardtTrustRegionCache, args...; p = cache.p, - u0 = cache.v_cache, kwargs...) +function reinit_cache!(cache::LevenbergMarquardtTrustRegionCache, args...; + p = cache.p, u0 = cache.v_cache, kwargs...) cache.p = p @bb copyto!(cache.v_cache, u0) cache.loss_old = oftype(cache.loss_old, Inf) @@ -57,18 +57,18 @@ function reinit_cache!(cache::LevenbergMarquardtTrustRegionCache, args...; p = c end function __internal_init( - prob::AbstractNonlinearProblem, alg::LevenbergMarquardtTrustRegion, - f::F, fu, u, p, args...; internalnorm::IF = DEFAULT_NORM, kwargs...) where {F, IF} + prob::AbstractNonlinearProblem, alg::LevenbergMarquardtTrustRegion, f::F, + fu, u, p, args...; internalnorm::IF = DEFAULT_NORM, kwargs...) where {F, IF} T = promote_type(eltype(u), eltype(fu)) @bb v = copy(u) @bb u_cache = similar(u) @bb fu_cache = similar(fu) - return LevenbergMarquardtTrustRegionCache(f, p, T(Inf), v, T(Inf), internalnorm, - alg.β_uphill, false, u_cache, fu_cache, 0) + return LevenbergMarquardtTrustRegionCache( + f, p, T(Inf), v, T(Inf), internalnorm, alg.β_uphill, false, u_cache, fu_cache, 0) end -function __internal_solve!(cache::LevenbergMarquardtTrustRegionCache, J, fu, u, δu, - descent_stats) +function __internal_solve!( + cache::LevenbergMarquardtTrustRegionCache, J, fu, u, δu, descent_stats) # This should be true if Geodesic Acceleration is being used v = hasfield(typeof(descent_stats), :v) ? descent_stats.v : δu norm_v = cache.internalnorm(v) @@ -236,8 +236,8 @@ the value used in the respective paper. - `expand_factor`: the factor to expand the trust region radius with if `expand_threshold < r` (with `r` defined in `shrink_threshold`). Defaults to `2.0`. """ -@kwdef @concrete struct GenericTrustRegionScheme{ - M <: RadiusUpdateSchemes.AbstractRadiusUpdateScheme} +@kwdef @concrete struct GenericTrustRegionScheme{M <: + RadiusUpdateSchemes.AbstractRadiusUpdateScheme} method::M = RadiusUpdateSchemes.Simple step_threshold = nothing shrink_threshold = nothing @@ -286,14 +286,15 @@ end alg end -function reinit_cache!(cache::GenericTrustRegionSchemeCache, args...; u0 = nothing, - p = cache.p, kwargs...) +function reinit_cache!( + cache::GenericTrustRegionSchemeCache, args...; u0 = nothing, p = cache.p, kwargs...) T = eltype(cache.u_cache) cache.p = p if u0 !== nothing u0_norm = cache.internalnorm(u0) - cache.trust_region = __initial_trust_radius(cache.alg.initial_trust_radius, T, - cache.alg.method, cache.max_trust_radius, u0_norm, u0_norm) # FIXME: scheme specific + cache.trust_region = __initial_trust_radius( + cache.alg.initial_trust_radius, T, cache.alg.method, + cache.max_trust_radius, u0_norm, u0_norm) # FIXME: scheme specific end cache.last_step_accepted = false cache.shrink_counter = 0 @@ -313,14 +314,15 @@ for func in (:__max_trust_radius, :__initial_trust_radius, :__step_threshold, end @inline __max_trust_radius(::Nothing, ::Type{T}, method, u, fu_norm) where {T} = T(Inf) -@inline function __max_trust_radius(::Nothing, ::Type{T}, - ::Union{RUS.__Simple, RUS.__NocedalWright}, u, fu_norm) where {T} +@inline function __max_trust_radius( + ::Nothing, ::Type{T}, ::Union{RUS.__Simple, RUS.__NocedalWright}, + u, fu_norm) where {T} u_min, u_max = extrema(u) return max(T(fu_norm), u_max - u_min) end -@inline function __initial_trust_radius(::Nothing, ::Type{T}, method, max_tr, - u0_norm, fu_norm) where {T} +@inline function __initial_trust_radius( + ::Nothing, ::Type{T}, method, max_tr, u0_norm, fu_norm) where {T} method isa RUS.__NLsolve && return T(ifelse(u0_norm > 0, u0_norm, 1)) (method isa RUS.__Hei || method isa RUS.__Bastin) && return T(1) method isa RUS.__Fan && return T((fu_norm^0.99) / 10) @@ -365,16 +367,17 @@ end @inline __expand_factor(::Nothing, ::Type{T}, method) where {T} = T(2) -function __internal_init(prob::AbstractNonlinearProblem, alg::GenericTrustRegionScheme, - f::F, fu, u, p, args...; internalnorm::IF = DEFAULT_NORM, kwargs...) where {F, IF} +function __internal_init( + prob::AbstractNonlinearProblem, alg::GenericTrustRegionScheme, f::F, fu, + u, p, args...; internalnorm::IF = DEFAULT_NORM, kwargs...) where {F, IF} T = promote_type(eltype(u), eltype(fu)) u0_norm = internalnorm(u) fu_norm = internalnorm(fu) # Common Setup max_trust_radius = __max_trust_radius(alg.max_trust_radius, T, alg.method, u, fu_norm) - initial_trust_radius = __initial_trust_radius(alg.initial_trust_radius, T, alg.method, - max_trust_radius, u0_norm, fu_norm) + initial_trust_radius = __initial_trust_radius( + alg.initial_trust_radius, T, alg.method, max_trust_radius, u0_norm, fu_norm) step_threshold = __step_threshold(alg.step_threshold, T, alg.method) shrink_threshold = __shrink_threshold(alg.shrink_threshold, T, alg.method) expand_threshold = __expand_threshold(alg.expand_threshold, T, alg.method) @@ -412,15 +415,15 @@ function __internal_init(prob::AbstractNonlinearProblem, alg::GenericTrustRegion @bb fu_cache = similar(fu) @bb Jδu_cache = similar(fu) - return GenericTrustRegionSchemeCache(alg.method, f, p, max_trust_radius, - initial_trust_radius, initial_trust_radius, step_threshold, shrink_threshold, - expand_threshold, shrink_factor, expand_factor, p1, p2, p3, p4, ϵ, T(0), - vjp_operator, jvp_operator, Jᵀfu_cache, Jδu_cache, δu_cache, internalnorm, - u_cache, fu_cache, false, 0, 0, alg) + return GenericTrustRegionSchemeCache( + alg.method, f, p, max_trust_radius, initial_trust_radius, initial_trust_radius, + step_threshold, shrink_threshold, expand_threshold, shrink_factor, + expand_factor, p1, p2, p3, p4, ϵ, T(0), vjp_operator, jvp_operator, Jᵀfu_cache, + Jδu_cache, δu_cache, internalnorm, u_cache, fu_cache, false, 0, 0, alg) end -function __internal_solve!(cache::GenericTrustRegionSchemeCache, J, fu, u, δu, - descent_stats) +function __internal_solve!( + cache::GenericTrustRegionSchemeCache, J, fu, u, δu, descent_stats) T = promote_type(eltype(u), eltype(fu)) @bb @. cache.u_cache = u + δu cache.fu_cache = evaluate_f!!(cache.f, cache.fu_cache, cache.u_cache, cache.p) @@ -462,8 +465,8 @@ function __internal_solve!(cache::GenericTrustRegionSchemeCache, J, fu, u, δu, if cache.ρ ≥ cache.expand_threshold cache.trust_region = cache.expand_factor * cache.internalnorm(δu) elseif cache.ρ ≥ cache.p1 - cache.trust_region = max(cache.trust_region, - cache.expand_factor * cache.internalnorm(δu)) + cache.trust_region = max( + cache.trust_region, cache.expand_factor * cache.internalnorm(δu)) end end elseif cache.method isa RUS.__NocedalWright @@ -473,14 +476,14 @@ function __internal_solve!(cache::GenericTrustRegionSchemeCache, J, fu, u, δu, else cache.shrink_counter = 0 if cache.ρ > cache.expand_threshold && - abs(cache.internalnorm(δu) - cache.trust_region) < - 1e-6 * cache.trust_region + abs(cache.internalnorm(δu) - cache.trust_region) < 1e-6 * cache.trust_region cache.trust_region = cache.expand_factor * cache.trust_region end end elseif cache.method isa RUS.__Hei - tr_new = __rfunc(cache.ρ, cache.shrink_threshold, cache.p1, cache.p3, cache.p4, - cache.p2) * cache.internalnorm(δu) + tr_new = __rfunc( + cache.ρ, cache.shrink_threshold, cache.p1, cache.p3, cache.p4, cache.p2) * + cache.internalnorm(δu) if tr_new < cache.trust_region cache.shrink_counter += 1 else @@ -507,16 +510,14 @@ function __internal_solve!(cache::GenericTrustRegionSchemeCache, J, fu, u, δu, cache.shrink_counter += 1 else cache.shrink_counter = 0 - cache.ρ > cache.expand_threshold && (cache.p1 = min(cache.p1 * cache.p3, - cache.p4)) + cache.ρ > cache.expand_threshold && + (cache.p1 = min(cache.p1 * cache.p3, cache.p4)) end cache.trust_region = cache.p1 * (cache.internalnorm(cache.fu_cache)^T(0.99)) elseif cache.method isa RUS.__Bastin if cache.ρ > cache.step_threshold - jvp_op = StatefulJacobianOperator(cache.jvp_operator, cache.u_cache, - cache.p) - vjp_op = StatefulJacobianOperator(cache.vjp_operator, cache.u_cache, - cache.p) + jvp_op = StatefulJacobianOperator(cache.jvp_operator, cache.u_cache, cache.p) + vjp_op = StatefulJacobianOperator(cache.vjp_operator, cache.u_cache, cache.p) @bb cache.Jδu_cache = jvp_op × vec(δu) @bb cache.Jᵀfu_cache = vjp_op × vec(cache.fu_cache) denom_1 = dot(_vec(δu), cache.Jᵀfu_cache) @@ -541,7 +542,6 @@ end # R-function for adaptive trust region method function __rfunc(r::R, c2::R, M::R, γ1::R, γ2::R, β::R) where {R <: Real} - return ifelse(r ≥ c2, - (2 * (M - 1 - γ2) * atan(r - c2) + (1 + γ2)) / R(π), + return ifelse(r ≥ c2, (2 * (M - 1 - γ2) * atan(r - c2) + (1 + γ2)) / R(π), (1 - γ1 - β) * (exp(r - c2) + β / (1 - γ1 - β))) end diff --git a/src/internal/approximate_initialization.jl b/src/internal/approximate_initialization.jl index bb9898009..60f4a7584 100644 --- a/src/internal/approximate_initialization.jl +++ b/src/internal/approximate_initialization.jl @@ -63,17 +63,16 @@ structure as specified by `structure`. structure end -function __internal_init(prob::AbstractNonlinearProblem, alg::IdentityInitialization, - solver, - f::F, fu, u::Number, p; internalnorm::IN = DEFAULT_NORM, kwargs...) where {F, IN} +function __internal_init( + prob::AbstractNonlinearProblem, alg::IdentityInitialization, solver, f::F, + fu, u::Number, p; internalnorm::IN = DEFAULT_NORM, kwargs...) where {F, IN} α = __initial_alpha(alg.alpha, u, fu, internalnorm) - return InitializedApproximateJacobianCache(α, alg.structure, alg, nothing, true, - internalnorm) + return InitializedApproximateJacobianCache( + α, alg.structure, alg, nothing, true, internalnorm) end function __internal_init(prob::AbstractNonlinearProblem, alg::IdentityInitialization, - solver, - f::F, fu::StaticArray, u::StaticArray, p; internalnorm::IN = DEFAULT_NORM, - kwargs...) where {IN, F} + solver, f::F, fu::StaticArray, u::StaticArray, p; + internalnorm::IN = DEFAULT_NORM, kwargs...) where {IN, F} α = __initial_alpha(alg.alpha, u, fu, internalnorm) if alg.structure isa DiagonalStructure @assert length(u)==length(fu) "Diagonal Jacobian Structure must be square!" @@ -87,11 +86,12 @@ function __internal_init(prob::AbstractNonlinearProblem, alg::IdentityInitializa end J = alg.structure(J_; alias = true) end - return InitializedApproximateJacobianCache(J, alg.structure, alg, nothing, true, - internalnorm) + return InitializedApproximateJacobianCache( + J, alg.structure, alg, nothing, true, internalnorm) end -function __internal_init(prob::AbstractNonlinearProblem, alg::IdentityInitialization, - solver, f::F, fu, u, p; internalnorm::IN = DEFAULT_NORM, kwargs...) where {F, IN} +function __internal_init( + prob::AbstractNonlinearProblem, alg::IdentityInitialization, solver, + f::F, fu, u, p; internalnorm::IN = DEFAULT_NORM, kwargs...) where {F, IN} α = __initial_alpha(alg.alpha, u, fu, internalnorm) if alg.structure isa DiagonalStructure @assert length(u)==length(fu) "Diagonal Jacobian Structure must be square!" @@ -100,8 +100,8 @@ function __internal_init(prob::AbstractNonlinearProblem, alg::IdentityInitializa J_ = similar(fu, promote_type(eltype(fu), eltype(u)), length(fu), length(u)) J = alg.structure(__make_identity!!(J_, α); alias = true) end - return InitializedApproximateJacobianCache(J, alg.structure, alg, nothing, true, - internalnorm) + return InitializedApproximateJacobianCache( + J, alg.structure, alg, nothing, true, internalnorm) end @inline function __initial_alpha(α, u, fu, internalnorm::F) where {F} @@ -145,15 +145,15 @@ make a selection automatically. autodiff end -function __internal_init(prob::AbstractNonlinearProblem, alg::TrueJacobianInitialization, - solver, f::F, fu, u, p; linsolve = missing, internalnorm::IN = DEFAULT_NORM, - kwargs...) where {F, IN} - autodiff = get_concrete_forward_ad(alg.autodiff, prob; check_reverse_mode = false, - kwargs...) +function __internal_init( + prob::AbstractNonlinearProblem, alg::TrueJacobianInitialization, solver, f::F, fu, + u, p; linsolve = missing, internalnorm::IN = DEFAULT_NORM, kwargs...) where {F, IN} + autodiff = get_concrete_forward_ad( + alg.autodiff, prob; check_reverse_mode = false, kwargs...) jac_cache = JacobianCache(prob, solver, prob.f, fu, u, p; autodiff, linsolve) J = alg.structure(jac_cache(nothing)) - return InitializedApproximateJacobianCache(J, alg.structure, alg, jac_cache, false, - internalnorm) + return InitializedApproximateJacobianCache( + J, alg.structure, alg, jac_cache, false, internalnorm) end """ @@ -205,8 +205,8 @@ function (cache::InitializedApproximateJacobianCache)(::Nothing) return get_full_jacobian(cache, cache.structure, cache.J) end -function __internal_solve!(cache::InitializedApproximateJacobianCache, fu, u, - ::Val{reinit}) where {reinit} +function __internal_solve!( + cache::InitializedApproximateJacobianCache, fu, u, ::Val{reinit}) where {reinit} if reinit || !cache.initialized cache(cache.alg, fu, u) cache.initialized = true @@ -225,8 +225,8 @@ function (cache::InitializedApproximateJacobianCache)(alg::IdentityInitializatio return end -function (cache::InitializedApproximateJacobianCache)(alg::TrueJacobianInitialization, fu, - u) +function (cache::InitializedApproximateJacobianCache)( + alg::TrueJacobianInitialization, fu, u) J_new = cache.cache(u) cache.J = cache.structure(cache.J, J_new) return diff --git a/src/internal/forward_diff.jl b/src/internal/forward_diff.jl index 4555c99a6..190c80645 100644 --- a/src/internal/forward_diff.jl +++ b/src/internal/forward_diff.jl @@ -1,17 +1,17 @@ # Not part of public API but helps reduce code duplication -import SimpleNonlinearSolve: __nlsolve_ad, - __nlsolve_dual_soln, __nlsolve_∂f_∂p, __nlsolve_∂f_∂u +import SimpleNonlinearSolve: __nlsolve_ad, __nlsolve_dual_soln, __nlsolve_∂f_∂p, + __nlsolve_∂f_∂u function SciMLBase.solve( - prob::NonlinearProblem{<:Union{Number, <:AbstractArray}, - iip, <:Union{<:Dual{T, V, P}, <:AbstractArray{<:Dual{T, V, P}}}}, - alg::Union{Nothing, AbstractNonlinearAlgorithm}, args...; + prob::NonlinearProblem{<:Union{Number, <:AbstractArray}, iip, + <:Union{<:Dual{T, V, P}, <:AbstractArray{<:Dual{T, V, P}}}}, + alg::Union{Nothing, AbstractNonlinearAlgorithm}, + args...; kwargs...) where {T, V, P, iip} sol, partials = __nlsolve_ad(prob, alg, args...; kwargs...) dual_soln = __nlsolve_dual_soln(sol.u, partials, prob.p) return SciMLBase.build_solution( - prob, alg, dual_soln, sol.resid; sol.retcode, sol.stats, - sol.original) + prob, alg, dual_soln, sol.resid; sol.retcode, sol.stats, sol.original) end @concrete mutable struct NonlinearSolveForwardDiffCache @@ -25,10 +25,9 @@ end @internal_caches NonlinearSolveForwardDiffCache :cache -function reinit_cache!(cache::NonlinearSolveForwardDiffCache; p = cache.p, - u0 = get_u(cache.cache), kwargs...) - inner_cache = reinit_cache!(cache.cache; p = __value(p), u0 = __value(u0), - kwargs...) +function reinit_cache!(cache::NonlinearSolveForwardDiffCache; + p = cache.p, u0 = get_u(cache.cache), kwargs...) + inner_cache = reinit_cache!(cache.cache; p = __value(p), u0 = __value(u0), kwargs...) cache.cache = inner_cache cache.p = p cache.values_p = __value(p) @@ -37,15 +36,16 @@ function reinit_cache!(cache::NonlinearSolveForwardDiffCache; p = cache.p, end function SciMLBase.init( - prob::NonlinearProblem{<:Union{Number, <:AbstractArray}, - iip, <:Union{<:Dual{T, V, P}, <:AbstractArray{<:Dual{T, V, P}}}}, - alg::Union{Nothing, AbstractNonlinearAlgorithm}, args...; + prob::NonlinearProblem{<:Union{Number, <:AbstractArray}, iip, + <:Union{<:Dual{T, V, P}, <:AbstractArray{<:Dual{T, V, P}}}}, + alg::Union{Nothing, AbstractNonlinearAlgorithm}, + args...; kwargs...) where {T, V, P, iip} p = __value(prob.p) newprob = NonlinearProblem(prob.f, __value(prob.u0), p; prob.kwargs...) cache = init(newprob, alg, args...; kwargs...) - return NonlinearSolveForwardDiffCache(cache, newprob, alg, prob.p, p, - ForwardDiff.partials(prob.p)) + return NonlinearSolveForwardDiffCache( + cache, newprob, alg, prob.p, p, ForwardDiff.partials(prob.p)) end function SciMLBase.solve!(cache::NonlinearSolveForwardDiffCache) @@ -66,8 +66,8 @@ function SciMLBase.solve!(cache::NonlinearSolveForwardDiffCache) end dual_soln = __nlsolve_dual_soln(sol.u, partials, cache.p) - return SciMLBase.build_solution(prob, cache.alg, dual_soln, sol.resid; sol.retcode, - sol.stats, sol.original) + return SciMLBase.build_solution( + prob, cache.alg, dual_soln, sol.resid; sol.retcode, sol.stats, sol.original) end @inline __value(x) = x diff --git a/src/internal/helpers.jl b/src/internal/helpers.jl index 4f475214b..fc43543ec 100644 --- a/src/internal/helpers.jl +++ b/src/internal/helpers.jl @@ -32,29 +32,27 @@ end # AutoDiff Selection Functions struct NonlinearSolveTag end -function ForwardDiff.checktag(::Type{<:ForwardDiff.Tag{<:NonlinearSolveTag, <:T}}, f::F, - x::AbstractArray{T}) where {T, F} +function ForwardDiff.checktag(::Type{<:ForwardDiff.Tag{<:NonlinearSolveTag, <:T}}, + f::F, x::AbstractArray{T}) where {T, F} return true end function get_concrete_forward_ad( - autodiff::Union{ADTypes.AbstractForwardMode, - ADTypes.AbstractFiniteDifferencesMode}, - prob, sp::Val{test_sparse} = True, - args...; kwargs...) where {test_sparse} + autodiff::Union{ADTypes.AbstractForwardMode, ADTypes.AbstractFiniteDifferencesMode}, + prob, sp::Val{test_sparse} = True, args...; kwargs...) where {test_sparse} return autodiff end -function get_concrete_forward_ad(autodiff::ADTypes.AbstractADType, prob, - sp::Val{test_sparse} = True, args...; - check_reverse_mode = true, kwargs...) where {test_sparse} +function get_concrete_forward_ad( + autodiff::ADTypes.AbstractADType, prob, sp::Val{test_sparse} = True, + args...; check_reverse_mode = true, kwargs...) where {test_sparse} if check_reverse_mode @warn "$(autodiff)::$(typeof(autodiff)) is not a \ `Abstract(Forward/FiniteDifferences)Mode`. Use with caution." maxlog=1 end return autodiff end -function get_concrete_forward_ad(autodiff, prob, sp::Val{test_sparse} = True, args...; - kwargs...) where {test_sparse} +function get_concrete_forward_ad( + autodiff, prob, sp::Val{test_sparse} = True, args...; kwargs...) where {test_sparse} if test_sparse (; sparsity, jac_prototype) = prob.f use_sparse_ad = sparsity !== nothing || jac_prototype !== nothing @@ -71,10 +69,8 @@ function get_concrete_forward_ad(autodiff, prob, sp::Val{test_sparse} = True, ar end function get_concrete_reverse_ad( - autodiff::Union{ADTypes.AbstractReverseMode, - ADTypes.AbstractFiniteDifferencesMode}, - prob, sp::Val{test_sparse} = True, - args...; kwargs...) where {test_sparse} + autodiff::Union{ADTypes.AbstractReverseMode, ADTypes.AbstractFiniteDifferencesMode}, + prob, sp::Val{test_sparse} = True, args...; kwargs...) where {test_sparse} return autodiff end function get_concrete_reverse_ad(autodiff::Union{AutoZygote, AutoSparseZygote}, prob, @@ -88,17 +84,17 @@ function get_concrete_reverse_ad(autodiff::Union{AutoZygote, AutoSparseZygote}, end return autodiff end -function get_concrete_reverse_ad(autodiff::ADTypes.AbstractADType, prob, - sp::Val{test_sparse} = True, args...; check_reverse_mode = true, - kwargs...) where {test_sparse} +function get_concrete_reverse_ad( + autodiff::ADTypes.AbstractADType, prob, sp::Val{test_sparse} = True, + args...; check_reverse_mode = true, kwargs...) where {test_sparse} if check_reverse_mode @warn "$(autodiff)::$(typeof(autodiff)) is not a \ `Abstract(Forward/FiniteDifferences)Mode`. Use with caution." maxlog=1 end return autodiff end -function get_concrete_reverse_ad(autodiff, prob, sp::Val{test_sparse} = True, args...; - kwargs...) where {test_sparse} +function get_concrete_reverse_ad( + autodiff, prob, sp::Val{test_sparse} = True, args...; kwargs...) where {test_sparse} if test_sparse (; sparsity, jac_prototype) = prob.f use_sparse_ad = sparsity !== nothing || jac_prototype !== nothing @@ -234,19 +230,18 @@ macro internal_caches(cType, internal_cache_names...) end function __internal_caches(__source__, __module__, cType, internal_cache_names::Tuple) - fields = map(name -> :($(__query_stat)(getproperty(cache, $(name)), ST)), - internal_cache_names) + fields = map( + name -> :($(__query_stat)(getproperty(cache, $(name)), ST)), internal_cache_names) callback_caches = map( - name -> :($(callback_into_cache!)(cache, - getproperty(internalcache, $(name)), internalcache, args...)), + name -> :($(callback_into_cache!)( + cache, getproperty(internalcache, $(name)), internalcache, args...)), internal_cache_names) callbacks_self = map( - name -> :($(callback_into_cache!)(internalcache, - getproperty(internalcache, $(name)))), + name -> :($(callback_into_cache!)( + internalcache, getproperty(internalcache, $(name)))), internal_cache_names) reinit_caches = map( - name -> :($(reinit_cache!)(getproperty(cache, $(name)), - args...; kwargs...)), + name -> :($(reinit_cache!)(getproperty(cache, $(name)), args...; kwargs...)), internal_cache_names) return esc(quote function __query_stat(cache::$(cType), ST::Val{stat}) where {stat} diff --git a/src/internal/jacobian.jl b/src/internal/jacobian.jl index 6f28d5999..fc4f12f33 100644 --- a/src/internal/jacobian.jl +++ b/src/internal/jacobian.jl @@ -42,24 +42,25 @@ Construct a cache for the Jacobian of `f` w.r.t. `u`. jvp_autodiff end -function reinit_cache!(cache::JacobianCache{iip}, args...; p = cache.p, u0 = cache.u, - kwargs...) where {iip} +function reinit_cache!(cache::JacobianCache{iip}, args...; p = cache.p, + u0 = cache.u, kwargs...) where {iip} cache.njacs = 0 cache.u = u0 cache.p = p cache.uf = JacobianWrapper{iip}(cache.f, p) end -function JacobianCache(prob, alg, f::F, fu_, u, p; autodiff = nothing, - vjp_autodiff = nothing, jvp_autodiff = nothing, linsolve = missing) where {F} +function JacobianCache( + prob, alg, f::F, fu_, u, p; autodiff = nothing, vjp_autodiff = nothing, + jvp_autodiff = nothing, linsolve = missing) where {F} iip = isinplace(prob) uf = JacobianWrapper{iip}(f, p) autodiff = get_concrete_forward_ad(autodiff, prob; check_reverse_mode = false) - jvp_autodiff = get_concrete_forward_ad(jvp_autodiff, prob, Val(false); - check_reverse_mode = true) - vjp_autodiff = get_concrete_reverse_ad(vjp_autodiff, prob, Val(false); - check_forward_mode = false) + jvp_autodiff = get_concrete_forward_ad( + jvp_autodiff, prob, Val(false); check_reverse_mode = true) + vjp_autodiff = get_concrete_reverse_ad( + vjp_autodiff, prob, Val(false); check_forward_mode = false) has_analytic_jac = SciMLBase.has_jac(f) linsolve_needs_jac = concrete_jac(alg) === nothing && (linsolve === missing || @@ -72,8 +73,8 @@ function JacobianCache(prob, alg, f::F, fu_, u, p; autodiff = nothing, if !has_analytic_jac && needs_jac sd = __sparsity_detection_alg(f, autodiff) jac_cache = iip ? sparse_jacobian_cache(autodiff, sd, uf, fu, u) : - sparse_jacobian_cache(autodiff, sd, uf, __maybe_mutable(u, autodiff); - fx = fu) + sparse_jacobian_cache( + autodiff, sd, uf, __maybe_mutable(u, autodiff); fx = fu) else jac_cache = nothing end @@ -91,14 +92,13 @@ function JacobianCache(prob, alg, f::F, fu_, u, p; autodiff = nothing, end return JacobianCache{iip}( - J, f, uf, fu, u, p, jac_cache, alg, 0, autodiff, vjp_autodiff, - jvp_autodiff) + J, f, uf, fu, u, p, jac_cache, alg, 0, autodiff, vjp_autodiff, jvp_autodiff) end function JacobianCache(prob, alg, f::F, ::Number, u::Number, p; kwargs...) where {F} uf = JacobianWrapper{false}(f, p) - return JacobianCache{false}(u, f, uf, u, u, p, nothing, alg, 0, nothing, nothing, - nothing) + return JacobianCache{false}( + u, f, uf, u, u, p, nothing, alg, 0, nothing, nothing, nothing) end @inline (cache::JacobianCache)(u = cache.u) = cache(cache.J, u, cache.p) @@ -117,8 +117,8 @@ function (cache::JacobianCache)(::Number, u, p = cache.p) # Scalar return J end # Compute the Jacobian -function (cache::JacobianCache{iip})(J::Union{AbstractMatrix, Nothing}, u, - p = cache.p) where {iip} +function (cache::JacobianCache{iip})( + J::Union{AbstractMatrix, Nothing}, u, p = cache.p) where {iip} cache.njacs += 1 if iip if has_jac(cache.f) diff --git a/src/internal/linear_solve.jl b/src/internal/linear_solve.jl index 319b80005..261d9e8dc 100644 --- a/src/internal/linear_solve.jl +++ b/src/internal/linear_solve.jl @@ -16,9 +16,9 @@ handled: ### Solving the System ```julia -(cache::LinearSolverCache)(; A = nothing, b = nothing, linu = nothing, - du = nothing, p = nothing, weight = nothing, cachedata = nothing, - reuse_A_if_factorization = false, kwargs...) +(cache::LinearSolverCache)(; + A = nothing, b = nothing, linu = nothing, du = nothing, p = nothing, + weight = nothing, cachedata = nothing, reuse_A_if_factorization = false, kwargs...) ``` Returns the solution of the system `u` and stores the updated cache in `cache.lincache`. @@ -67,8 +67,10 @@ end function LinearSolverCache(alg, linsolve, A, b, u; kwargs...) u_fixed = __fix_strange_type_combination(A, b, u) - if (A isa Number && b isa Number) || (linsolve === nothing && A isa SMatrix) || - (A isa Diagonal) || (linsolve isa typeof(\)) + if (A isa Number && b isa Number) || + (linsolve === nothing && A isa SMatrix) || + (A isa Diagonal) || + (linsolve isa typeof(\)) return LinearSolverCache(nothing, nothing, A, b, nothing, 0, 0) end @bb u_ = copy(u_fixed) @@ -77,8 +79,8 @@ function LinearSolverCache(alg, linsolve, A, b, u; kwargs...) weight = __init_ones(u_fixed) if __hasfield(alg, Val(:precs)) precs = alg.precs - Pl_, Pr_ = precs(A, nothing, u, nothing, nothing, nothing, nothing, nothing, - nothing) + Pl_, Pr_ = precs( + A, nothing, u, nothing, nothing, nothing, nothing, nothing, nothing) else precs, Pl_, Pr_ = nothing, nothing, nothing end @@ -91,8 +93,8 @@ function LinearSolverCache(alg, linsolve, A, b, u; kwargs...) end # Direct Linear Solve Case without Caching -function (cache::LinearSolverCache{Nothing})(; A = nothing, b = nothing, linu = nothing, - kwargs...) +function (cache::LinearSolverCache{Nothing})(; + A = nothing, b = nothing, linu = nothing, kwargs...) cache.nsolve += 1 cache.nfactors += 1 A === nothing || (cache.A = A) @@ -107,9 +109,9 @@ function (cache::LinearSolverCache{Nothing})(; A = nothing, b = nothing, linu = return res end # Use LinearSolve.jl -function (cache::LinearSolverCache)(; A = nothing, b = nothing, linu = nothing, - du = nothing, p = nothing, weight = nothing, cachedata = nothing, - reuse_A_if_factorization = false, kwargs...) +function (cache::LinearSolverCache)(; + A = nothing, b = nothing, linu = nothing, du = nothing, p = nothing, + weight = nothing, cachedata = nothing, reuse_A_if_factorization = false, kwargs...) cache.nsolve += 1 __update_A!(cache, A, reuse_A_if_factorization) @@ -124,8 +126,8 @@ function (cache::LinearSolverCache)(; A = nothing, b = nothing, linu = nothing, if cache.precs === nothing _Pl, _Pr = nothing, nothing else - _Pl, _Pr = cache.precs(cache.lincache.A, du, linu, p, nothing, A !== nothing, - Plprev, Prprev, cachedata) + _Pl, _Pr = cache.precs(cache.lincache.A, du, linu, p, nothing, + A !== nothing, Plprev, Prprev, cachedata) end if (_Pl !== nothing || _Pr !== nothing) diff --git a/src/internal/operators.jl b/src/internal/operators.jl index d5e6f4145..a34565a9d 100644 --- a/src/internal/operators.jl +++ b/src/internal/operators.jl @@ -16,9 +16,9 @@ transposed, and `T = false` means that the Jacobian is not transposed. ### Constructor ```julia -JacobianOperator(prob::AbstractNonlinearProblem, fu, u; jvp_autodiff = nothing, - vjp_autodiff = nothing, skip_vjp::Val{NoVJP} = False, - skip_jvp::Val{NoJVP} = False) where {NoVJP, NoJVP} +JacobianOperator( + prob::AbstractNonlinearProblem, fu, u; jvp_autodiff = nothing, vjp_autodiff = nothing, + skip_vjp::Val{NoVJP} = False, skip_jvp::Val{NoJVP} = False) where {NoVJP, NoJVP} ``` See also [`NonlinearSolve.VecJacOperator`](@ref) and @@ -45,8 +45,8 @@ end for op in (:adjoint, :transpose) @eval function Base.$(op)(operator::JacobianOperator{vjp, iip, T}) where {vjp, iip, T} - return JacobianOperator{!vjp, iip, T}(operator.jvp_op, operator.vjp_op, - operator.output_cache, operator.input_cache) + return JacobianOperator{!vjp, iip, T}( + operator.jvp_op, operator.vjp_op, operator.output_cache, operator.input_cache) end end @@ -68,8 +68,8 @@ function JacobianOperator(prob::AbstractNonlinearProblem, fu, u; jvp_autodiff = @closure (v, u, p) -> FiniteDiff.finite_difference_derivative(uf, u) * v end else - vjp_autodiff = __get_nonsparse_ad(get_concrete_reverse_ad(vjp_autodiff, - prob, False)) + vjp_autodiff = __get_nonsparse_ad(get_concrete_reverse_ad( + vjp_autodiff, prob, False)) if vjp_autodiff isa AutoZygote iip && error("`AutoZygote` cannot handle inplace problems.") @closure (v, u, p) -> auto_vecjac(uf, u, v) @@ -98,17 +98,15 @@ function JacobianOperator(prob::AbstractNonlinearProblem, fu, u; jvp_autodiff = @closure (v, u, p) -> FiniteDiff.finite_difference_derivative(uf, u) * v end else - jvp_autodiff = __get_nonsparse_ad(get_concrete_forward_ad(jvp_autodiff, - prob, False)) + jvp_autodiff = __get_nonsparse_ad(get_concrete_forward_ad( + jvp_autodiff, prob, False)) if jvp_autodiff isa AutoForwardDiff || jvp_autodiff isa AutoPolyesterForwardDiff if iip # FIXME: Technically we should propagate the tag but ignoring that for now - cache1 = Dual{ - typeof(ForwardDiff.Tag(NonlinearSolveTag(), eltype(u))), eltype(u), 1 - }.(similar(u), ForwardDiff.Partials.(tuple.(u))) - cache2 = Dual{ - typeof(ForwardDiff.Tag(NonlinearSolveTag(), eltype(fu))), eltype(fu), 1 - }.(similar(fu), ForwardDiff.Partials.(tuple.(fu))) + cache1 = Dual{typeof(ForwardDiff.Tag(NonlinearSolveTag(), eltype(u))), + eltype(u), 1}.(similar(u), ForwardDiff.Partials.(tuple.(u))) + cache2 = Dual{typeof(ForwardDiff.Tag(NonlinearSolveTag(), eltype(fu))), + eltype(fu), 1}.(similar(fu), ForwardDiff.Partials.(tuple.(fu))) @closure (Jv, v, u, p) -> auto_jacvec!(Jv, uf, u, v, cache1, cache2) else @closure (v, u, p) -> auto_jacvec(uf, u, v) @@ -128,8 +126,7 @@ function JacobianOperator(prob::AbstractNonlinearProblem, fu, u; jvp_autodiff = end return JacobianOperator{false, iip, promote_type(eltype(fu), eltype(u))}( - jvp_op, vjp_op, - u, fu) + jvp_op, vjp_op, u, fu) end """ @@ -208,8 +205,8 @@ end Wrapper over a [`JacobianOperator`](@ref) which stores the input `u` and `p` and defines `mul!` and `*` for computing VJPs and JVPs. """ -@concrete struct StatefulJacobianOperator{vjp, iip, T, - J <: JacobianOperator{vjp, iip, T}} <: AbstractNonlinearSolveOperator{T} +@concrete struct StatefulJacobianOperator{ + vjp, iip, T, J <: JacobianOperator{vjp, iip, T}} <: AbstractNonlinearSolveOperator{T} jac_op::J u p @@ -230,8 +227,8 @@ function Base.:*(J_op::StatefulJacobianOperator{vjp, iip, T, J, <:Number}, return J_op.jac_op(v, J_op.u, J_op.p) end -function LinearAlgebra.mul!(Jv::AbstractArray, J::StatefulJacobianOperator, - v::AbstractArray) +function LinearAlgebra.mul!( + Jv::AbstractArray, J::StatefulJacobianOperator, v::AbstractArray) J.jac_op(Jv, v, J.u, J.p) return Jv end @@ -271,8 +268,8 @@ function Base.:*(JᵀJ::StatefulJacobianNormalFormOperator, x::AbstractArray) return JᵀJ.vjp_operator * (JᵀJ.jvp_operator * x) end -function LinearAlgebra.mul!(JᵀJx::AbstractArray, JᵀJ::StatefulJacobianNormalFormOperator, - x::AbstractArray) +function LinearAlgebra.mul!( + JᵀJx::AbstractArray, JᵀJ::StatefulJacobianNormalFormOperator, x::AbstractArray) mul!(JᵀJ.cache, JᵀJ.jvp_operator, x) mul!(JᵀJx, JᵀJ.vjp_operator, JᵀJ.cache) return JᵀJx diff --git a/src/internal/termination.jl b/src/internal/termination.jl index 59d8905f5..e55e344f7 100644 --- a/src/internal/termination.jl +++ b/src/internal/termination.jl @@ -1,6 +1,6 @@ function init_termination_cache(abstol, reltol, du, u, ::Nothing) - return init_termination_cache(abstol, reltol, du, u, - AbsSafeBestTerminationMode(; max_stalled_steps = 32)) + return init_termination_cache( + abstol, reltol, du, u, AbsSafeBestTerminationMode(; max_stalled_steps = 32)) end function init_termination_cache(abstol, reltol, du, u, tc::AbstractNonlinearTerminationMode) tc_cache = init(du, u, tc; abstol, reltol, use_deprecated_retcodes = Val(false)) @@ -12,8 +12,8 @@ function check_and_update!(cache, fu, u, uprev) end function check_and_update!(tc_cache, cache, fu, u, uprev) - return check_and_update!(tc_cache, cache, fu, u, uprev, - DiffEqBase.get_termination_mode(tc_cache)) + return check_and_update!( + tc_cache, cache, fu, u, uprev, DiffEqBase.get_termination_mode(tc_cache)) end function check_and_update!(tc_cache, cache, fu, u, uprev, mode) @@ -25,17 +25,17 @@ function check_and_update!(tc_cache, cache, fu, u, uprev, mode) end function update_from_termination_cache!(tc_cache, cache, u = get_u(cache)) - return update_from_termination_cache!(tc_cache, cache, - DiffEqBase.get_termination_mode(tc_cache), u) + return update_from_termination_cache!( + tc_cache, cache, DiffEqBase.get_termination_mode(tc_cache), u) end -function update_from_termination_cache!(tc_cache, cache, - mode::AbstractNonlinearTerminationMode, u = get_u(cache)) +function update_from_termination_cache!( + tc_cache, cache, mode::AbstractNonlinearTerminationMode, u = get_u(cache)) evaluate_f!(cache, u, cache.p) end -function update_from_termination_cache!(tc_cache, cache, - mode::AbstractSafeBestNonlinearTerminationMode, u = get_u(cache)) +function update_from_termination_cache!( + tc_cache, cache, mode::AbstractSafeBestNonlinearTerminationMode, u = get_u(cache)) if isinplace(cache) copyto!(get_u(cache), tc_cache.u) else diff --git a/src/internal/tracing.jl b/src/internal/tracing.jl index 667c6ce07..8f0f10fb1 100644 --- a/src/internal/tracing.jl +++ b/src/internal/tracing.jl @@ -90,13 +90,13 @@ function Base.show(io::IO, entry::NonlinearSolveTraceEntry) end function NonlinearSolveTraceEntry(iteration, fu, δu) - return NonlinearSolveTraceEntry(iteration, norm(fu, Inf), norm(δu, 2), nothing, - nothing, nothing, nothing, nothing) + return NonlinearSolveTraceEntry( + iteration, norm(fu, Inf), norm(δu, 2), nothing, nothing, nothing, nothing, nothing) end function NonlinearSolveTraceEntry(iteration, fu, δu, J) - return NonlinearSolveTraceEntry(iteration, norm(fu, Inf), norm(δu, 2), __cond(J), - nothing, nothing, nothing, nothing) + return NonlinearSolveTraceEntry(iteration, norm(fu, Inf), norm(δu, 2), + __cond(J), nothing, nothing, nothing, nothing) end function NonlinearSolveTraceEntry(iteration, fu, δu, J, u) @@ -104,8 +104,8 @@ function NonlinearSolveTraceEntry(iteration, fu, δu, J, u) __copy(J), __copy(u), __copy(fu), __copy(δu)) end -@concrete struct NonlinearSolveTrace{show_trace, store_trace, - Tr <: AbstractNonlinearSolveTraceLevel} +@concrete struct NonlinearSolveTrace{ + show_trace, store_trace, Tr <: AbstractNonlinearSolveTraceLevel} history trace_level::Tr end @@ -126,25 +126,26 @@ end function init_nonlinearsolve_trace(alg, u, fu, J, δu; show_trace::Val = Val(false), trace_level::AbstractNonlinearSolveTraceLevel = TraceMinimal(), store_trace::Val = Val(false), uses_jac_inverse = Val(false), kwargs...) - return init_nonlinearsolve_trace(alg, show_trace, trace_level, store_trace, u, fu, J, - δu, uses_jac_inverse) + return init_nonlinearsolve_trace( + alg, show_trace, trace_level, store_trace, u, fu, J, δu, uses_jac_inverse) end -function init_nonlinearsolve_trace(alg, ::Val{show_trace}, - trace_level::AbstractNonlinearSolveTraceLevel, ::Val{store_trace}, u, fu, J, - δu, ::Val{uses_jac_inverse}) where {show_trace, store_trace, uses_jac_inverse} +function init_nonlinearsolve_trace( + alg, ::Val{show_trace}, trace_level::AbstractNonlinearSolveTraceLevel, + ::Val{store_trace}, u, fu, J, δu, + ::Val{uses_jac_inverse}) where {show_trace, store_trace, uses_jac_inverse} if show_trace print("\nAlgorithm: ") Base.printstyled(alg, "\n\n"; color = :green, bold = true) end J_ = uses_jac_inverse ? (trace_level isa TraceMinimal ? J : __safe_inv(J)) : J - history = __init_trace_history(Val{show_trace}(), trace_level, Val{store_trace}(), u, - fu, J_, δu) + history = __init_trace_history( + Val{show_trace}(), trace_level, Val{store_trace}(), u, fu, J_, δu) return NonlinearSolveTrace{show_trace, store_trace}(history, trace_level) end -function __init_trace_history(::Val{show_trace}, trace_level, ::Val{store_trace}, u, fu, J, - δu) where {show_trace, store_trace} +function __init_trace_history(::Val{show_trace}, trace_level, ::Val{store_trace}, + u, fu, J, δu) where {show_trace, store_trace} !store_trace && !show_trace && return nothing entry = __trace_entry(trace_level, 0, u, fu, J, δu) show_trace && show(entry) @@ -167,16 +168,16 @@ function update_trace!(trace::NonlinearSolveTrace{ShT, StT}, iter, u, fu, J, δu !StT && !ShT && return nothing if L - entry = NonlinearSolveTraceEntry(-1, norm(fu, Inf), NaN32, nothing, nothing, - nothing, nothing, nothing) + entry = NonlinearSolveTraceEntry( + -1, norm(fu, Inf), NaN32, nothing, nothing, nothing, nothing, nothing) ShT && show(entry) return trace end show_now = ShT && (mod1(iter, trace.trace_level.print_frequency) == 1) store_now = StT && (mod1(iter, trace.trace_level.store_frequency) == 1) - (show_now || store_now) && (entry = __trace_entry(trace.trace_level, iter, u, fu, J, - δu, α)) + (show_now || store_now) && + (entry = __trace_entry(trace.trace_level, iter, u, fu, J, δu, α)) store_now && push!(trace.history, entry) show_now && show(entry) return trace @@ -188,13 +189,13 @@ function update_trace!(cache::AbstractNonlinearSolveCache, α = true) J = __getproperty(cache, Val(:J)) if J === nothing - update_trace!(trace, get_nsteps(cache) + 1, get_u(cache), get_fu(cache), - nothing, cache.du, α) + update_trace!( + trace, get_nsteps(cache) + 1, get_u(cache), get_fu(cache), nothing, cache.du, α) elseif cache isa ApproximateJacobianSolveCache && store_inverse_jacobian(cache) - update_trace!(trace, get_nsteps(cache) + 1, get_u(cache), get_fu(cache), - ApplyArray(__safe_inv, J), cache.du, α) + update_trace!(trace, get_nsteps(cache) + 1, get_u(cache), + get_fu(cache), ApplyArray(__safe_inv, J), cache.du, α) else - update_trace!(trace, get_nsteps(cache) + 1, get_u(cache), get_fu(cache), J, - cache.du, α) + update_trace!( + trace, get_nsteps(cache) + 1, get_u(cache), get_fu(cache), J, cache.du, α) end end diff --git a/test/core/23_test_problems_tests.jl b/test/core/23_test_problems_tests.jl index 43a361958..7924e372d 100644 --- a/test/core/23_test_problems_tests.jl +++ b/test/core/23_test_problems_tests.jl @@ -4,8 +4,8 @@ using NonlinearSolve, LinearAlgebra, LinearSolve, NonlinearProblemLibrary, Test problems = NonlinearProblemLibrary.problems dicts = NonlinearProblemLibrary.dicts -function test_on_library(problems, dicts, alg_ops, broken_tests, ϵ = 1e-4; - skip_tests = nothing) +function test_on_library( + problems, dicts, alg_ops, broken_tests, ϵ = 1e-4; skip_tests = nothing) for (idx, (problem, dict)) in enumerate(zip(problems, dicts)) x = dict["start"] res = similar(x) @@ -71,8 +71,7 @@ end @testitem "LevenbergMarquardt" setup=[RobustnessTesting] begin using LinearSolve - alg_ops = (LevenbergMarquardt(), - LevenbergMarquardt(; α_geodesic = 0.1), + alg_ops = (LevenbergMarquardt(), LevenbergMarquardt(; α_geodesic = 0.1), LevenbergMarquardt(; linsolve = CholeskyFactorization())) broken_tests = Dict(alg => Int[] for alg in alg_ops) @@ -93,8 +92,7 @@ end end @testitem "Broyden" retries=5 setup=[RobustnessTesting] begin - alg_ops = (Broyden(), - Broyden(; init_jacobian = Val(:true_jacobian)), + alg_ops = (Broyden(), Broyden(; init_jacobian = Val(:true_jacobian)), Broyden(; update_rule = Val(:bad_broyden)), Broyden(; init_jacobian = Val(:true_jacobian), update_rule = Val(:bad_broyden))) diff --git a/test/core/forward_ad_tests.jl b/test/core/forward_ad_tests.jl index 22882cbec..a13259939 100644 --- a/test/core/forward_ad_tests.jl +++ b/test/core/forward_ad_tests.jl @@ -63,10 +63,9 @@ export test_f!, test_f, jacobian_f, solve_with, __compatible end @testitem "ForwardDiff.jl Integration" setup=[ForwardADTesting] begin - for alg in (NewtonRaphson(), TrustRegion(), - LevenbergMarquardt(), PseudoTransient(; alpha_initial = 10.0), Broyden(), Klement(), - DFSane(), nothing, NLsolveJL(), CMINPACK(), - KINSOL(; globalization_strategy = :LineSearch)) + for alg in (NewtonRaphson(), TrustRegion(), LevenbergMarquardt(), + PseudoTransient(; alpha_initial = 10.0), Broyden(), Klement(), DFSane(), nothing, + NLsolveJL(), CMINPACK(), KINSOL(; globalization_strategy = :LineSearch)) us = (2.0, @SVector[1.0, 1.0], [1.0, 1.0], ones(2, 2), @SArray ones(2, 2)) @testset "Scalar AD" begin @@ -90,7 +89,8 @@ end end @testset "Jacobian" begin - for u0 in us, p in ([2.0, 1.0], [2.0 1.0; 3.0 4.0]), + for u0 in us, + p in ([2.0, 1.0], [2.0 1.0; 3.0 4.0]), mode in (:iip, :oop, :iip_cache, :oop_cache) __compatible(u0, p) || continue diff --git a/test/core/nlls_tests.jl b/test/core/nlls_tests.jl index c95fd0de0..3722ca671 100644 --- a/test/core/nlls_tests.jl +++ b/test/core/nlls_tests.jl @@ -1,7 +1,7 @@ @testsetup module CoreNLLSTesting using Reexport -@reexport using NonlinearSolve, - LinearSolve, LinearAlgebra, StableRNGs, Random, ForwardDiff, Zygote +@reexport using NonlinearSolve, LinearSolve, LinearAlgebra, StableRNGs, Random, ForwardDiff, + Zygote true_function(x, θ) = @. θ[1] * exp(θ[2] * x) * cos(θ[3] * x + θ[4]) true_function(y, x, θ) = (@. y = θ[1] * exp(θ[2] * x) * cos(θ[3] * x + θ[4])) @@ -36,16 +36,12 @@ for linsolve in [nothing, LUFactorization(), KrylovJL_GMRES(), KrylovJL_LSMR()] end end append!(solvers, - [ - LevenbergMarquardt(), - LevenbergMarquardt(; linsolve = LUFactorization()), + [LevenbergMarquardt(), LevenbergMarquardt(; linsolve = LUFactorization()), LevenbergMarquardt(; linsolve = KrylovJL_GMRES()), - LevenbergMarquardt(; linsolve = KrylovJL_LSMR()), - nothing - ]) + LevenbergMarquardt(; linsolve = KrylovJL_LSMR()), nothing]) for radius_update_scheme in [RadiusUpdateSchemes.Simple, RadiusUpdateSchemes.NocedalWright, - RadiusUpdateSchemes.NLsolve, RadiusUpdateSchemes.Hei, RadiusUpdateSchemes.Yuan, - RadiusUpdateSchemes.Fan, RadiusUpdateSchemes.Bastin] + RadiusUpdateSchemes.NLsolve, RadiusUpdateSchemes.Hei, + RadiusUpdateSchemes.Yuan, RadiusUpdateSchemes.Fan, RadiusUpdateSchemes.Bastin] push!(solvers, TrustRegion(; radius_update_scheme)) end @@ -55,8 +51,7 @@ end @testitem "General NLLS Solvers" setup=[CoreNLLSTesting] begin prob_oop = NonlinearLeastSquaresProblem{false}(loss_function, θ_init, x) prob_iip = NonlinearLeastSquaresProblem( - NonlinearFunction(loss_function; - resid_prototype = zero(y_target)), θ_init, x) + NonlinearFunction(loss_function; resid_prototype = zero(y_target)), θ_init, x) nlls_problems = [prob_oop, prob_iip] @@ -84,16 +79,15 @@ end probs = [ NonlinearLeastSquaresProblem( - NonlinearFunction{true}(loss_function; - resid_prototype = zero(y_target), vjp = vjp!), + NonlinearFunction{true}( + loss_function; resid_prototype = zero(y_target), vjp = vjp!), θ_init, x), NonlinearLeastSquaresProblem( - NonlinearFunction{false}(loss_function; - resid_prototype = zero(y_target), vjp = vjp), + NonlinearFunction{false}( + loss_function; resid_prototype = zero(y_target), vjp = vjp), θ_init, - x) - ] + x)] for prob in probs, solver in solvers sol = solve(prob, solver; maxiters = 10000, abstol = 1e-8) diff --git a/test/core/rootfind_tests.jl b/test/core/rootfind_tests.jl index 6a9402651..86dc5ff2c 100644 --- a/test/core/rootfind_tests.jl +++ b/test/core/rootfind_tests.jl @@ -25,17 +25,15 @@ function newton_fails(u, p) (0.21640425613334457 .+ 216.40425613334457 ./ (1 .+ (0.21640425613334457 .+ - 216.40425613334457 ./ - (1 .+ 0.0006250000000000001(u .^ 2.0))) .^ 2.0)) .^ 2.0) .- - 0.0011552453009332421u .- p + 216.40425613334457 ./ (1 .+ 0.0006250000000000001(u .^ 2.0))) .^ 2.0)) .^ + 2.0) .- 0.0011552453009332421u .- p end const TERMINATION_CONDITIONS = [ SteadyStateDiffEqTerminationMode(), SimpleNonlinearSolveTerminationMode(), NormTerminationMode(), RelTerminationMode(), RelNormTerminationMode(), AbsTerminationMode(), AbsNormTerminationMode(), RelSafeTerminationMode(), - AbsSafeTerminationMode(), RelSafeBestTerminationMode(), AbsSafeBestTerminationMode() -] + AbsSafeTerminationMode(), RelSafeBestTerminationMode(), AbsSafeBestTerminationMode()] function benchmark_nlsolve_oop(f, u0, p = 2.0; solver, kwargs...) prob = NonlinearProblem{false}(f, u0, p) @@ -68,8 +66,7 @@ end @testitem "NewtonRaphson" setup=[CoreRootfindTesting] begin @testset "LineSearch: $(_nameof(lsmethod)) LineSearch AD: $(_nameof(ad))" for lsmethod in ( - Static(), - StrongWolfe(), BackTracking(), HagerZhang(), MoreThuente()), + Static(), StrongWolfe(), BackTracking(), HagerZhang(), MoreThuente()), ad in (AutoFiniteDiff(), AutoZygote()) linesearch = LineSearchesJL(; method = lsmethod, autodiff = ad) @@ -81,18 +78,19 @@ end @test SciMLBase.successful_retcode(sol) @test all(abs.(sol.u .* sol.u .- 2) .< 1e-9) - cache = init(NonlinearProblem{false}(quadratic_f, u0, 2.0), NewtonRaphson(), - abstol = 1e-9) + cache = init(NonlinearProblem{false}(quadratic_f, u0, 2.0), + NewtonRaphson(), abstol = 1e-9) @test (@ballocated solve!($cache)) < 200 end - precs = [ - (u0) -> NonlinearSolve.DEFAULT_PRECS, - u0 -> ((args...) -> (Diagonal(rand!(similar(u0))), nothing)) - ] + precs = [(u0) -> NonlinearSolve.DEFAULT_PRECS, + u0 -> ((args...) -> (Diagonal(rand!(similar(u0))), nothing))] @testset "[IIP] u0: $(typeof(u0)) precs: $(_nameof(prec)) linsolve: $(_nameof(linsolve))" for u0 in ([ - 1.0, 1.0],), prec in precs, linsolve in (nothing, KrylovJL_GMRES(), \) + 1.0, 1.0],), + prec in precs, + linsolve in (nothing, KrylovJL_GMRES(), \) + ad isa AutoZygote && continue if prec === :Random prec = (args...) -> (Diagonal(randn!(similar(u0))), nothing) @@ -114,8 +112,8 @@ end @test nlprob_iterator_interface(quadratic_f!, p, Val(true), NewtonRaphson()) ≈ sqrt.(p) @testset "ADType: $(autodiff) u0: $(_nameof(u0))" for autodiff in ( - AutoSparseForwardDiff(), - AutoSparseFiniteDiff(), AutoZygote(), AutoSparseZygote(), __autosparseenzyme()), + AutoSparseForwardDiff(), AutoSparseFiniteDiff(), + AutoZygote(), AutoSparseZygote(), __autosparseenzyme()), u0 in (1.0, [1.0, 1.0]) probN = NonlinearProblem(quadratic_f, u0, 2.0) @@ -134,13 +132,14 @@ end @testitem "TrustRegion" setup=[CoreRootfindTesting] begin radius_update_schemes = [RadiusUpdateSchemes.Simple, RadiusUpdateSchemes.NocedalWright, - RadiusUpdateSchemes.NLsolve, RadiusUpdateSchemes.Hei, RadiusUpdateSchemes.Yuan, - RadiusUpdateSchemes.Fan, RadiusUpdateSchemes.Bastin] + RadiusUpdateSchemes.NLsolve, RadiusUpdateSchemes.Hei, + RadiusUpdateSchemes.Yuan, RadiusUpdateSchemes.Fan, RadiusUpdateSchemes.Bastin] u0s = ([1.0, 1.0], @SVector[1.0, 1.0], 1.0) linear_solvers = [nothing, LUFactorization(), KrylovJL_GMRES(), \] @testset "[OOP] u0: $(typeof(u0)) radius_update_scheme: $(radius_update_scheme) linear_solver: $(linsolve)" for u0 in u0s, - radius_update_scheme in radius_update_schemes, linsolve in linear_solvers + radius_update_scheme in radius_update_schemes, + linsolve in linear_solvers abstol = ifelse(linsolve isa KrylovJL, 1e-6, 1e-9) @@ -155,7 +154,8 @@ end end @testset "[IIP] u0: $(typeof(u0)) radius_update_scheme: $(radius_update_scheme) linear_solver: $(linsolve)" for u0 in ([ - 1.0, 1.0],), radius_update_scheme in radius_update_schemes, + 1.0, 1.0],), + radius_update_scheme in radius_update_schemes, linsolve in linear_solvers abstol = ifelse(linsolve isa KrylovJL, 1e-6, 1e-9) @@ -175,8 +175,8 @@ end @test nlprob_iterator_interface(quadratic_f!, p, Val(true), TrustRegion()) ≈ sqrt.(p) @testset "ADType: $(autodiff) u0: $(_nameof(u0)) radius_update_scheme: $(radius_update_scheme)" for autodiff in ( - AutoSparseForwardDiff(), - AutoSparseFiniteDiff(), AutoZygote(), AutoSparseZygote(), __autosparseenzyme()), + AutoSparseForwardDiff(), AutoSparseFiniteDiff(), + AutoZygote(), AutoSparseZygote(), __autosparseenzyme()), u0 in (1.0, [1.0, 1.0]), radius_update_scheme in radius_update_schemes @@ -207,16 +207,16 @@ end 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) + list_of_options = zip( + max_trust_radius, initial_trust_radius, step_threshold, shrink_threshold, + expand_threshold, shrink_factor, expand_factor, max_shrink_times) for options in list_of_options local probN, sol, alg - alg = TrustRegion(max_trust_radius = options[1], - initial_trust_radius = options[2], step_threshold = options[3], - shrink_threshold = options[4], expand_threshold = options[5], - shrink_factor = options[6], expand_factor = options[7], - max_shrink_times = options[8]) + alg = TrustRegion( + max_trust_radius = options[1], initial_trust_radius = options[2], + step_threshold = options[3], shrink_threshold = options[4], + expand_threshold = options[5], shrink_factor = options[6], + expand_factor = options[7], max_shrink_times = options[8]) probN = NonlinearProblem{false}(quadratic_f, [1.0, 1.0], 2.0) sol = solve(probN, alg, abstol = 1e-10) @@ -255,8 +255,8 @@ end @test SciMLBase.successful_retcode(sol) @test all(abs.(sol.u .* sol.u .- 2) .< 1e-9) - cache = init(NonlinearProblem{false}(quadratic_f, u0, 2.0), LevenbergMarquardt(), - abstol = 1e-9) + cache = init(NonlinearProblem{false}(quadratic_f, u0, 2.0), + LevenbergMarquardt(), abstol = 1e-9) @test (@ballocated solve!($cache)) < 200 end @@ -265,19 +265,20 @@ end @test SciMLBase.successful_retcode(sol) @test all(abs.(sol.u .* sol.u .- 2) .< 1e-9) - cache = init(NonlinearProblem{true}(quadratic_f!, u0, 2.0), LevenbergMarquardt(), - abstol = 1e-9) + cache = init(NonlinearProblem{true}(quadratic_f!, u0, 2.0), + LevenbergMarquardt(), abstol = 1e-9) @test (@ballocated solve!($cache)) ≤ 64 end @testset "ADType: $(autodiff) u0: $(_nameof(u0))" for autodiff in ( - AutoSparseForwardDiff(), - AutoSparseFiniteDiff(), AutoZygote(), AutoSparseZygote(), __autosparseenzyme()), + AutoSparseForwardDiff(), AutoSparseFiniteDiff(), + AutoZygote(), AutoSparseZygote(), __autosparseenzyme()), u0 in (1.0, [1.0, 1.0]) probN = NonlinearProblem(quadratic_f, u0, 2.0) - @test all(solve(probN, LevenbergMarquardt(; autodiff); abstol = 1e-9, - reltol = 1e-9).u .≈ sqrt(2.0)) + @test all(solve( + probN, LevenbergMarquardt(; autodiff); abstol = 1e-9, reltol = 1e-9).u .≈ + sqrt(2.0)) end # Test that `LevenbergMarquardt` passes a test that `NewtonRaphson` fails on. @@ -291,10 +292,10 @@ end # Iterator interface p = range(0.01, 2, length = 200) - @test abs.(nlprob_iterator_interface(quadratic_f, p, Val(false), - LevenbergMarquardt())) ≈ sqrt.(p) - @test abs.(nlprob_iterator_interface(quadratic_f!, p, Val(true), - LevenbergMarquardt())) ≈ sqrt.(p) + @test abs.(nlprob_iterator_interface( + quadratic_f, p, Val(false), LevenbergMarquardt())) ≈ sqrt.(p) + @test abs.(nlprob_iterator_interface( + quadratic_f!, p, Val(true), LevenbergMarquardt())) ≈ sqrt.(p) # Test kwargs in `LevenbergMarquardt` @testset "Keyword Arguments" begin @@ -306,13 +307,13 @@ end b_uphill = Float64[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) + list_of_options = zip( + damping_initial, damping_increase_factor, damping_decrease_factor, + 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], + 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]) @@ -341,8 +342,7 @@ end @test SciMLBase.successful_retcode(sol) @test all(abs.(sol.u .* sol.u .- 2) .< 1e-9) - cache = init(NonlinearProblem{false}(quadratic_f, u0, 2.0), DFSane(), - abstol = 1e-9) + cache = init(NonlinearProblem{false}(quadratic_f, u0, 2.0), DFSane(), abstol = 1e-9) @test (@ballocated solve!($cache)) < 200 end @@ -379,19 +379,15 @@ end τ_min = [0.1, 0.2, 0.3] τ_max = [0.5, 0.8, 0.9] nexp = [2, 1, 2] - η_strategy = [ - (f_1, k, x, F) -> f_1 / k^2, - (f_1, k, x, F) -> f_1 / k^3, - (f_1, k, x, F) -> f_1 / k^4 - ] - - list_of_options = zip(σ_min, σ_max, σ_1, M, γ, τ_min, τ_max, nexp, - η_strategy) + η_strategy = [(f_1, k, x, F) -> f_1 / k^2, (f_1, k, x, F) -> f_1 / k^3, + (f_1, k, x, F) -> f_1 / k^4] + + list_of_options = zip(σ_min, σ_max, σ_1, M, γ, τ_min, τ_max, nexp, η_strategy) for options in list_of_options local probN, sol, alg alg = DFSane(σ_min = options[1], σ_max = options[2], σ_1 = options[3], - M = options[4], γ = options[5], τ_min = options[6], τ_max = options[7], - n_exp = options[8], η_strategy = options[9]) + M = options[4], γ = options[5], τ_min = options[6], + τ_max = options[7], n_exp = options[8], η_strategy = options[9]) probN = NonlinearProblem{false}(quadratic_f, [1.0, 1.0], 2.0) sol = solve(probN, alg, abstol = 1e-11) @@ -412,8 +408,8 @@ end @testitem "PseudoTransient" setup=[CoreRootfindTesting] begin # These are tests for NewtonRaphson so we should set alpha_initial to be high so that we # converge quickly - @testset "PT: alpha_initial = 10.0 PT AD: $(ad)" for ad in (AutoFiniteDiff(), - AutoZygote()) + @testset "PT: alpha_initial = 10.0 PT AD: $(ad)" for ad in ( + AutoFiniteDiff(), AutoZygote()) u0s = ([1.0, 1.0], @SVector[1.0, 1.0], 1.0) @testset "[OOP] u0: $(typeof(u0))" for u0 in u0s @@ -431,7 +427,10 @@ end precs = [NonlinearSolve.DEFAULT_PRECS, :Random] @testset "[IIP] u0: $(typeof(u0)) precs: $(_nameof(prec)) linsolve: $(_nameof(linsolve))" for u0 in ([ - 1.0, 1.0],), prec in precs, linsolve in (nothing, KrylovJL_GMRES(), \) + 1.0, 1.0],), + prec in precs, + linsolve in (nothing, KrylovJL_GMRES(), \) + ad isa AutoZygote && continue if prec === :Random prec = (args...) -> (Diagonal(randn!(similar(u0))), nothing) @@ -441,21 +440,21 @@ end @test SciMLBase.successful_retcode(sol) @test all(abs.(sol.u .* sol.u .- 2) .< 1e-9) - cache = init(NonlinearProblem{true}(quadratic_f!, u0, 2.0), solver; - abstol = 1e-9) + cache = init( + NonlinearProblem{true}(quadratic_f!, u0, 2.0), solver; abstol = 1e-9) @test (@ballocated solve!($cache)) ≤ 64 end end p = range(0.01, 2, length = 200) - @test nlprob_iterator_interface(quadratic_f, p, Val(false), - PseudoTransient(; alpha_initial = 10.0)) ≈ sqrt.(p) - @test nlprob_iterator_interface(quadratic_f!, p, Val(true), - PseudoTransient(; alpha_initial = 10.0)) ≈ sqrt.(p) + @test nlprob_iterator_interface( + quadratic_f, p, Val(false), PseudoTransient(; alpha_initial = 10.0)) ≈ sqrt.(p) + @test nlprob_iterator_interface( + quadratic_f!, p, Val(true), PseudoTransient(; alpha_initial = 10.0)) ≈ sqrt.(p) @testset "ADType: $(autodiff) u0: $(_nameof(u0))" for autodiff in ( - AutoSparseForwardDiff(), - AutoSparseFiniteDiff(), AutoZygote(), AutoSparseZygote(), __autosparseenzyme()), + AutoSparseForwardDiff(), AutoSparseFiniteDiff(), + AutoZygote(), AutoSparseZygote(), __autosparseenzyme()), u0 in (1.0, [1.0, 1.0]) probN = NonlinearProblem(quadratic_f, u0, 2.0) @@ -467,8 +466,9 @@ end u0 in (1.0, [1.0, 1.0]) probN = NonlinearProblem(quadratic_f, u0, 2.0) - @test all(solve(probN, PseudoTransient(; alpha_initial = 10.0); - termination_condition).u .≈ sqrt(2.0)) + @test all(solve( + probN, PseudoTransient(; alpha_initial = 10.0); termination_condition).u .≈ + sqrt(2.0)) end end @@ -476,9 +476,8 @@ end @testitem "Broyden" setup=[CoreRootfindTesting] begin @testset "LineSearch: $(_nameof(lsmethod)) LineSearch AD: $(_nameof(ad)) Init Jacobian: $(init_jacobian) Update Rule: $(update_rule)" for lsmethod in ( - Static(), - StrongWolfe(), BackTracking(), HagerZhang(), MoreThuente(), - LiFukushimaLineSearch()), + Static(), StrongWolfe(), BackTracking(), HagerZhang(), + MoreThuente(), LiFukushimaLineSearch()), ad in (AutoFiniteDiff(), AutoZygote()), init_jacobian in (Val(:identity), Val(:true_jacobian)), update_rule in (Val(:good_broyden), Val(:bad_broyden), Val(:diagonal)) @@ -527,8 +526,7 @@ end @testitem "Klement" setup=[CoreRootfindTesting] begin @testset "LineSearch: $(_nameof(lsmethod)) LineSearch AD: $(_nameof(ad)) Init Jacobian: $(init_jacobian)" for lsmethod in ( - Static(), - StrongWolfe(), BackTracking(), HagerZhang(), MoreThuente()), + Static(), StrongWolfe(), BackTracking(), HagerZhang(), MoreThuente()), ad in (AutoFiniteDiff(), AutoZygote()), init_jacobian in (Val(:identity), Val(:true_jacobian), Val(:true_jacobian_diagonal)) @@ -577,9 +575,8 @@ end @testitem "LimitedMemoryBroyden" setup=[CoreRootfindTesting] begin @testset "LineSearch: $(_nameof(lsmethod)) LineSearch AD: $(_nameof(ad))" for lsmethod in ( - Static(), - StrongWolfe(), BackTracking(), HagerZhang(), MoreThuente(), - LiFukushimaLineSearch()), + Static(), StrongWolfe(), BackTracking(), HagerZhang(), + MoreThuente(), LiFukushimaLineSearch()), ad in (AutoFiniteDiff(), AutoZygote()) linesearch = LineSearchesJL(; method = lsmethod, autodiff = ad) @@ -611,17 +608,17 @@ end # Iterator interface p = range(0.01, 2, length = 200) - @test nlprob_iterator_interface(quadratic_f, p, Val(false), - LimitedMemoryBroyden())≈sqrt.(p) atol=1e-2 - @test nlprob_iterator_interface(quadratic_f!, p, Val(true), - LimitedMemoryBroyden())≈sqrt.(p) atol=1e-2 + @test nlprob_iterator_interface( + quadratic_f, p, Val(false), LimitedMemoryBroyden())≈sqrt.(p) atol=1e-2 + @test nlprob_iterator_interface( + quadratic_f!, p, Val(true), LimitedMemoryBroyden())≈sqrt.(p) atol=1e-2 @testset "Termination condition: $(termination_condition) u0: $(_nameof(u0))" for termination_condition in TERMINATION_CONDITIONS, u0 in (1.0, [1.0, 1.0]) probN = NonlinearProblem(quadratic_f, u0, 2.0) - @test all(solve(probN, LimitedMemoryBroyden(); - termination_condition).u .≈ sqrt(2.0)) + @test all(solve(probN, LimitedMemoryBroyden(); termination_condition).u .≈ + sqrt(2.0)) end end @@ -643,8 +640,8 @@ end return v + 0.1 * (u .* Δ * v + v .* Δ * u) end - function JVP!(du::Vector{Float64}, v::Vector{Float64}, u::Vector{Float64}, - p::Vector{Float64}) + function JVP!( + du::Vector{Float64}, v::Vector{Float64}, u::Vector{Float64}, p::Vector{Float64}) Δ = Tridiagonal(-ones(99), 2 * ones(100), -ones(99)) du .= v + 0.1 * (u .* Δ * v + v .* Δ * u) return nothing @@ -655,16 +652,16 @@ end prob = NonlinearProblem(NonlinearFunction{false}(F; jvp = JVP), u0, u0) sol = solve(prob, NewtonRaphson(; linsolve = KrylovJL_GMRES()); abstol = 1e-13) @test norm(sol.resid, Inf) ≤ 1e-6 - sol = solve(prob, - TrustRegion(; linsolve = KrylovJL_GMRES(), vjp_autodiff = AutoFiniteDiff()); + sol = solve( + prob, TrustRegion(; linsolve = KrylovJL_GMRES(), vjp_autodiff = AutoFiniteDiff()); abstol = 1e-13) @test norm(sol.resid, Inf) ≤ 1e-6 prob = NonlinearProblem(NonlinearFunction{true}(F!; jvp = JVP!), u0, u0) sol = solve(prob, NewtonRaphson(; linsolve = KrylovJL_GMRES()); abstol = 1e-13) @test norm(sol.resid, Inf) ≤ 1e-6 - sol = solve(prob, - TrustRegion(; linsolve = KrylovJL_GMRES(), vjp_autodiff = AutoFiniteDiff()); + sol = solve( + prob, TrustRegion(; linsolve = KrylovJL_GMRES(), vjp_autodiff = AutoFiniteDiff()); abstol = 1e-13) @test norm(sol.resid, Inf) ≤ 1e-6 end diff --git a/test/gpu/core_tests.jl b/test/gpu/core_tests.jl index 074098072..6e9806498 100644 --- a/test/gpu/core_tests.jl +++ b/test/gpu/core_tests.jl @@ -12,12 +12,11 @@ prob = NonlinearProblem(linear_f, u0) SOLVERS = (NewtonRaphson(), LevenbergMarquardt(; linsolve = QRFactorization()), - LevenbergMarquardt(; linsolve = KrylovJL_GMRES()), PseudoTransient(), Klement(), - Broyden(; linesearch = LiFukushimaLineSearch()), + LevenbergMarquardt(; linsolve = KrylovJL_GMRES()), PseudoTransient(), + Klement(), Broyden(; linesearch = LiFukushimaLineSearch()), LimitedMemoryBroyden(; threshold = 2, linesearch = LiFukushimaLineSearch()), DFSane(), TrustRegion(; linsolve = QRFactorization()), - TrustRegion(; linsolve = KrylovJL_GMRES(), - concrete_jac = true), # Needed if Zygote not loaded + TrustRegion(; linsolve = KrylovJL_GMRES(), concrete_jac = true), # Needed if Zygote not loaded nothing) @testset "[IIP] GPU Solvers" begin diff --git a/test/misc/bruss_tests.jl b/test/misc/bruss_tests.jl index a1fcd6761..08eae5436 100644 --- a/test/misc/bruss_tests.jl +++ b/test/misc/bruss_tests.jl @@ -18,12 +18,12 @@ 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] + + B + + u[i, j, 1]^2 * u[i, j, 2] - (A + 1) * u[i, j, 1] + brusselator_f(x, y) du[i, j, 2] = alpha * (u[im1, j, 2] + u[ip1, j, 2] + u[i, jp1, 2] + u[i, jm1, 2] - - 4u[i, j, 2]) + - A * u[i, j, 1] - u[i, j, 1]^2 * u[i, j, 2] + 4u[i, j, 2]) + A * u[i, j, 1] - u[i, j, 1]^2 * u[i, j, 2] end end @@ -46,17 +46,17 @@ sol = solve(prob_brusselator_2d, NewtonRaphson(); abstol = 1e-8) @test norm(sol.resid, Inf) < 1e-8 - sol = solve(prob_brusselator_2d, NewtonRaphson(autodiff = AutoSparseForwardDiff()); - abstol = 1e-8) + sol = solve(prob_brusselator_2d, + NewtonRaphson(autodiff = AutoSparseForwardDiff()); abstol = 1e-8) @test norm(sol.resid, Inf) < 1e-8 - sol = solve(prob_brusselator_2d, NewtonRaphson(autodiff = AutoSparseFiniteDiff()); - abstol = 1e-8) + sol = solve(prob_brusselator_2d, + NewtonRaphson(autodiff = AutoSparseFiniteDiff()); abstol = 1e-8) @test norm(sol.resid, Inf) < 1e-8 du0 = copy(u0) - jac_sparsity = Symbolics.jacobian_sparsity((du, u) -> brusselator_2d_loop(du, u, p), - du0, u0) + jac_sparsity = Symbolics.jacobian_sparsity( + (du, u) -> brusselator_2d_loop(du, u, p), du0, u0) jac_prototype = float.(jac_sparsity) fill!(jac_prototype, 0) @test all(iszero, jac_prototype) @@ -68,8 +68,8 @@ @test norm(sol.resid, Inf) < 1e-8 @test !all(iszero, jac_prototype) - sol = solve(prob_brusselator_2d, NewtonRaphson(autodiff = AutoSparseFiniteDiff()); - abstol = 1e-8) + sol = solve(prob_brusselator_2d, + NewtonRaphson(autodiff = AutoSparseFiniteDiff()); abstol = 1e-8) @test norm(sol.resid, Inf) < 1e-8 cache = init(prob_brusselator_2d, NewtonRaphson(; autodiff = AutoSparseForwardDiff())) diff --git a/test/misc/matrix_resizing_tests.jl b/test/misc/matrix_resizing_tests.jl index e13704c96..1d17a696a 100644 --- a/test/misc/matrix_resizing_tests.jl +++ b/test/misc/matrix_resizing_tests.jl @@ -7,9 +7,9 @@ vecprob = NonlinearProblem(ff, vec(u0), p) prob = NonlinearProblem(ff, u0, p) - for alg in (NewtonRaphson(), TrustRegion(), LevenbergMarquardt(), - PseudoTransient(), RobustMultiNewton(), FastShortcutNonlinearPolyalg(), Broyden(), - Klement(), LimitedMemoryBroyden(; threshold = 2)) + for alg in (NewtonRaphson(), TrustRegion(), LevenbergMarquardt(), PseudoTransient(), + RobustMultiNewton(), FastShortcutNonlinearPolyalg(), + Broyden(), Klement(), LimitedMemoryBroyden(; threshold = 2)) @test vec(solve(prob, alg).u) == solve(vecprob, alg).u end end @@ -23,9 +23,9 @@ end vecprob = NonlinearProblem(fiip, vec(u0), p) prob = NonlinearProblem(fiip, u0, p) - for alg in (NewtonRaphson(), TrustRegion(), LevenbergMarquardt(), - PseudoTransient(), RobustMultiNewton(), FastShortcutNonlinearPolyalg(), Broyden(), - Klement(), LimitedMemoryBroyden(; threshold = 2)) + for alg in (NewtonRaphson(), TrustRegion(), LevenbergMarquardt(), PseudoTransient(), + RobustMultiNewton(), FastShortcutNonlinearPolyalg(), + Broyden(), Klement(), LimitedMemoryBroyden(; threshold = 2)) @test vec(solve(prob, alg).u) == solve(vecprob, alg).u end end diff --git a/test/misc/polyalg_tests.jl b/test/misc/polyalg_tests.jl index 10df60429..d9433d494 100644 --- a/test/misc/polyalg_tests.jl +++ b/test/misc/polyalg_tests.jl @@ -91,8 +91,8 @@ end # This test is not meant to return success but test that all the default solvers can handle # complex valued problems @test_nowarn solve(prob; abstol = 1e-19, maxiters = 10) - @test_nowarn solve(prob, RobustMultiNewton(eltype(prob.u0)); abstol = 1e-19, - maxiters = 10) + @test_nowarn solve( + prob, RobustMultiNewton(eltype(prob.u0)); abstol = 1e-19, maxiters = 10) end @testitem "No AD" begin diff --git a/test/misc/qa_tests.jl b/test/misc/qa_tests.jl index e182ac22f..c43b9f027 100644 --- a/test/misc/qa_tests.jl +++ b/test/misc/qa_tests.jl @@ -4,8 +4,8 @@ Aqua.find_persistent_tasks_deps(NonlinearSolve) Aqua.test_ambiguities(NonlinearSolve; recursive = false) Aqua.test_deps_compat(NonlinearSolve) - Aqua.test_piracies(NonlinearSolve, - treat_as_own = [NonlinearProblem, NonlinearLeastSquaresProblem]) + Aqua.test_piracies( + NonlinearSolve, treat_as_own = [NonlinearProblem, NonlinearLeastSquaresProblem]) Aqua.test_project_extras(NonlinearSolve) # Timer Outputs needs to be enabled via Preferences Aqua.test_stale_deps(NonlinearSolve; ignore = [:TimerOutputs]) diff --git a/test/runtests.jl b/test/runtests.jl index 0c3f7d855..014e4cd0e 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,8 +3,7 @@ using ReTestItems const GROUP = get(ENV, "GROUP", "All") if GROUP == "All" || GROUP == "Core" - ReTestItems.runtests(joinpath(@__DIR__, "core/"), - joinpath(@__DIR__, "misc/"), + ReTestItems.runtests(joinpath(@__DIR__, "core/"), joinpath(@__DIR__, "misc/"), joinpath(@__DIR__, "wrappers/")) end diff --git a/test/wrappers/nlls_tests.jl b/test/wrappers/nlls_tests.jl index f578b192c..a65ab5ec9 100644 --- a/test/wrappers/nlls_tests.jl +++ b/test/wrappers/nlls_tests.jl @@ -31,8 +31,7 @@ end @testitem "LeastSquaresOptim.jl" setup=[WrapperNLLSSetup] begin prob_oop = NonlinearLeastSquaresProblem{false}(loss_function, θ_init, x) prob_iip = NonlinearLeastSquaresProblem( - NonlinearFunction(loss_function; - resid_prototype = zero(y_target)), θ_init, x) + NonlinearFunction(loss_function; resid_prototype = zero(y_target)), θ_init, x) nlls_problems = [prob_oop, prob_iip] @@ -47,8 +46,7 @@ end end end -@testitem "FastLevenbergMarquardt.jl + CMINPACK: Jacobian Provided" setup=[ - WrapperNLLSSetup] begin +@testitem "FastLevenbergMarquardt.jl + CMINPACK: Jacobian Provided" setup=[WrapperNLLSSetup] begin function jac!(J, θ, p) resid = zeros(length(p)) ForwardDiff.jacobian!(J, (resid, θ) -> loss_function(resid, θ, p), resid, θ) @@ -59,18 +57,17 @@ end probs = [ NonlinearLeastSquaresProblem( - NonlinearFunction{true}(loss_function; - resid_prototype = zero(y_target), jac = jac!), + NonlinearFunction{true}( + loss_function; resid_prototype = zero(y_target), jac = jac!), θ_init, x), NonlinearLeastSquaresProblem( - NonlinearFunction{false}(loss_function; - resid_prototype = zero(y_target), jac = jac), + NonlinearFunction{false}( + loss_function; resid_prototype = zero(y_target), jac = jac), θ_init, x), - NonlinearLeastSquaresProblem(NonlinearFunction{false}(loss_function; jac), - θ_init, x) - ] + NonlinearLeastSquaresProblem( + NonlinearFunction{false}(loss_function; jac), θ_init, x)] solvers = Any[FastLevenbergMarquardtJL(linsolve) for linsolve in (:cholesky, :qr)] push!(solvers, CMINPACK()) @@ -80,19 +77,15 @@ end end end -@testitem "FastLevenbergMarquardt.jl + CMINPACK: Jacobian Not Provided" setup=[ - WrapperNLLSSetup] begin +@testitem "FastLevenbergMarquardt.jl + CMINPACK: Jacobian Not Provided" setup=[WrapperNLLSSetup] begin probs = [ NonlinearLeastSquaresProblem( - NonlinearFunction{true}(loss_function; - resid_prototype = zero(y_target)), + NonlinearFunction{true}(loss_function; resid_prototype = zero(y_target)), θ_init, x), NonlinearLeastSquaresProblem( - NonlinearFunction{false}(loss_function; - resid_prototype = zero(y_target)), + NonlinearFunction{false}(loss_function; resid_prototype = zero(y_target)), θ_init, x), - NonlinearLeastSquaresProblem(NonlinearFunction{false}(loss_function), θ_init, x) - ] + NonlinearLeastSquaresProblem(NonlinearFunction{false}(loss_function), θ_init, x)] solvers = vec(Any[FastLevenbergMarquardtJL(linsolve; autodiff) for linsolve in (:cholesky, :qr),