diff --git a/Project.toml b/Project.toml index 46f2191da..9891f3cfb 100644 --- a/Project.toml +++ b/Project.toml @@ -33,12 +33,16 @@ UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed" BandedMatrices = "aae01518-5342-5314-be14-df237901396f" FastLevenbergMarquardt = "7a0df574-e128-4d35-8cbd-3d84502bf7ce" LeastSquaresOptim = "0fc2ff8b-aaa3-5acd-a817-1944a5e08891" +MINPACK = "4854310b-de5a-5eb6-a2a5-c1dee2bd17f9" +NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" [extensions] NonlinearSolveBandedMatricesExt = "BandedMatrices" NonlinearSolveFastLevenbergMarquardtExt = "FastLevenbergMarquardt" NonlinearSolveLeastSquaresOptimExt = "LeastSquaresOptim" +NonlinearSolveMINPACKExt = "MINPACK" +NonlinearSolveNLsolveExt = "NLsolve" NonlinearSolveZygoteExt = "Zygote" [compat] @@ -60,7 +64,9 @@ LeastSquaresOptim = "0.8" LineSearches = "7" LinearAlgebra = "<0.0.1, 1" LinearSolve = "2.12" +MINPACK = "1.2" MaybeInplace = "0.1" +NLsolve = "4.5" NaNMath = "1" NonlinearProblemLibrary = "0.1" Pkg = "1" @@ -94,17 +100,19 @@ ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" LeastSquaresOptim = "0fc2ff8b-aaa3-5acd-a817-1944a5e08891" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" +MINPACK = "4854310b-de5a-5eb6-a2a5-c1dee2bd17f9" +NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" NaNMath = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3" NonlinearProblemLibrary = "b7050fa9-e91f-4b37-bcee-a89a063da141" Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" -StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" SparseDiffTools = "47a9eef4-7e08-11e9-0b38-333d64bd3804" +StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" [targets] -test = ["Aqua", "Enzyme", "BenchmarkTools", "SafeTestsets", "Pkg", "Test", "ForwardDiff", "StaticArrays", "Symbolics", "LinearSolve", "Random", "LinearAlgebra", "Zygote", "SparseDiffTools", "NonlinearProblemLibrary", "LeastSquaresOptim", "FastLevenbergMarquardt", "NaNMath", "BandedMatrices", "DiffEqBase", "StableRNGs"] +test = ["Aqua", "Enzyme", "BenchmarkTools", "SafeTestsets", "Pkg", "Test", "ForwardDiff", "StaticArrays", "Symbolics", "LinearSolve", "Random", "LinearAlgebra", "Zygote", "SparseDiffTools", "NonlinearProblemLibrary", "LeastSquaresOptim", "FastLevenbergMarquardt", "NaNMath", "BandedMatrices", "DiffEqBase", "StableRNGs", "MINPACK", "NLsolve"] diff --git a/docs/Project.toml b/docs/Project.toml index f92656632..30583f843 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -8,10 +8,8 @@ IncompleteLU = "40713840-3770-5561-ab4c-a76e7d0d7895" LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec" -NonlinearSolveMINPACK = "c100e077-885d-495a-a2ea-599e143bf69d" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" -SciMLNLSolve = "e9a6253c-8580-4d32-9898-8661bb511710" SimpleNonlinearSolve = "727e6d20-b764-4bd8-a329-72de5adea6c7" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" SteadyStateDiffEq = "9672c7b4-1e72-59bd-8a11-6ac3964bc41f" @@ -31,7 +29,6 @@ NonlinearSolve = "3" NonlinearSolveMINPACK = "0.1" Random = "<0.0.1, 1" SciMLBase = "2.4" -SciMLNLSolve = "0.1" SimpleNonlinearSolve = "1" StaticArrays = "1" SteadyStateDiffEq = "2" diff --git a/docs/make.jl b/docs/make.jl index 9a063bf47..ed350c4bc 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,5 +1,6 @@ -using Documenter, NonlinearSolve, SimpleNonlinearSolve, Sundials, SciMLNLSolve, - NonlinearSolveMINPACK, SteadyStateDiffEq, SciMLBase, DiffEqBase +using Documenter, + NonlinearSolve, SimpleNonlinearSolve, Sundials, + SteadyStateDiffEq, SciMLBase, DiffEqBase cp("./docs/Manifest.toml", "./docs/src/assets/Manifest.toml", force = true) cp("./docs/Project.toml", "./docs/src/assets/Project.toml", force = true) @@ -9,7 +10,7 @@ include("pages.jl") makedocs(sitename = "NonlinearSolve.jl", authors = "Chris Rackauckas", modules = [NonlinearSolve, SciMLBase, DiffEqBase, SimpleNonlinearSolve, Sundials, - SciMLNLSolve, NonlinearSolveMINPACK, SteadyStateDiffEq], + SteadyStateDiffEq], clean = true, doctest = false, linkcheck = true, linkcheck_ignore = ["https://twitter.com/ChrisRackauckas/status/1544743542094020615"], warnonly = [:missing_docs, :cross_references], diff --git a/docs/src/api/fastlevenbergmarquardt.md b/docs/src/api/fastlevenbergmarquardt.md index b970cf819..8709dc303 100644 --- a/docs/src/api/fastlevenbergmarquardt.md +++ b/docs/src/api/fastlevenbergmarquardt.md @@ -1,17 +1,15 @@ # FastLevenbergMarquardt.jl -This is a wrapper package for importing solvers from FastLevenbergMarquardt.jl into the -SciML interface. Note that these solvers do not come by default, and thus one needs to -install the package before using these solvers: +This is a extension for importing solvers from FastLevenbergMarquardt.jl into the SciML +interface. Note that these solvers do not come by default, and thus one needs to install +the package before using these solvers: ```julia using Pkg Pkg.add("FastLevenbergMarquardt") -using FastLevenbergMarquardt +using FastLevenbergMarquardt, NonlinearSolve ``` -These methods can be used independently of the rest of NonlinearSolve.jl - ## Solver API ```@docs diff --git a/docs/src/api/leastsquaresoptim.md b/docs/src/api/leastsquaresoptim.md index 89c07d8b3..76850555b 100644 --- a/docs/src/api/leastsquaresoptim.md +++ b/docs/src/api/leastsquaresoptim.md @@ -1,17 +1,15 @@ # LeastSquaresOptim.jl -This is a wrapper package for importing solvers from LeastSquaresOptim.jl into the SciML +This is a extension for importing solvers from LeastSquaresOptim.jl into the SciML interface. Note that these solvers do not come by default, and thus one needs to install the package before using these solvers: ```julia using Pkg Pkg.add("LeastSquaresOptim") -using LeastSquaresOptim +using LeastSquaresOptim, NonlinearSolve ``` -These methods can be used independently of the rest of NonlinearSolve.jl - ## Solver API ```@docs diff --git a/docs/src/api/minpack.md b/docs/src/api/minpack.md index 36bec764f..d55491c0e 100644 --- a/docs/src/api/minpack.md +++ b/docs/src/api/minpack.md @@ -1,17 +1,15 @@ # MINPACK.jl -This is a wrapper package for importing solvers from Sundials into the SciML interface. -Note that these solvers do not come by default, and thus one needs to install the package -before using these solvers: +This is a extension for importing solvers from MINPACK into the SciML interface. Note that +these solvers do not come by default, and thus one needs to install the package before using +these solvers: ```julia using Pkg -Pkg.add("NonlinearSolveMINPACK") -using NonlinearSolveMINPACK +Pkg.add("MINPACK") +using MINPACK, NonlinearSolve ``` -These methods can be used independently of the rest of NonlinearSolve.jl - ## Solver API ```@docs diff --git a/docs/src/api/nlsolve.md b/docs/src/api/nlsolve.md index f0b611cc9..6de117872 100644 --- a/docs/src/api/nlsolve.md +++ b/docs/src/api/nlsolve.md @@ -1,13 +1,13 @@ # NLsolve.jl -This is a wrapper package for importing solvers from NLsolve.jl into the SciML interface. -Note that these solvers do not come by default, and thus one needs to install the package -before using these solvers: +This is a extension for importing solvers from NLsolve.jl into the SciML interface. Note +that these solvers do not come by default, and thus one needs to install the package before +using these solvers: ```julia using Pkg -Pkg.add("SciMLNLSolve") -using SciMLNLSolve +Pkg.add("NLsolve") +using NLSolve, NonlinearSolve ``` ## Solver API diff --git a/docs/src/solvers/NonlinearLeastSquaresSolvers.md b/docs/src/solvers/NonlinearLeastSquaresSolvers.md index fbd68c847..7adfd9508 100644 --- a/docs/src/solvers/NonlinearLeastSquaresSolvers.md +++ b/docs/src/solvers/NonlinearLeastSquaresSolvers.md @@ -58,6 +58,23 @@ And the choices for `linsolve` are: - `:cholesky` - `:lsmr` a conjugate gradient method (LSMR with diagonal preconditioner). +### MINPACK.jl + +MINPACK.jl methods are fine for medium-sized nonlinear solves. They are the FORTRAN +standard methods which are used in many places, such as SciPy. However, our benchmarks +demonstrate that these methods are not robust or stable. In addition, they are slower +than the standard methods and do not scale due to lack of sparse Jacobian support. +Thus they are only recommended for benchmarking and testing code conversions. + + - `CMINPACK()`: A wrapper for using the classic MINPACK method through [MINPACK.jl](https://github.com/sglyon/MINPACK.jl) + +Submethod choices for this algorithm include: + + - `:hybr`: Modified version of Powell's algorithm. + - `:lm`: Levenberg-Marquardt. + - `:lmdif`: Advanced Levenberg-Marquardt + - `:hybrd`: Advanced modified version of Powell's algorithm + ### Optimization.jl `NonlinearLeastSquaresProblem`s can be converted into an `OptimizationProblem` and used diff --git a/docs/src/solvers/NonlinearSystemSolvers.md b/docs/src/solvers/NonlinearSystemSolvers.md index 078131ff8..55e799efb 100644 --- a/docs/src/solvers/NonlinearSystemSolvers.md +++ b/docs/src/solvers/NonlinearSystemSolvers.md @@ -110,7 +110,7 @@ computationally expensive than direct methods. close to the steady state. - `SSRootfind()`: Uses a NonlinearSolve compatible solver to find the steady state. -### SciMLNLSolve.jl +### NLsolve.jl This is a wrapper package for importing solvers from NLsolve.jl into the SciML interface. diff --git a/ext/NonlinearSolveMINPACKExt.jl b/ext/NonlinearSolveMINPACKExt.jl new file mode 100644 index 000000000..588e93fa4 --- /dev/null +++ b/ext/NonlinearSolveMINPACKExt.jl @@ -0,0 +1,73 @@ +module NonlinearSolveMINPACKExt + +using NonlinearSolve, SciMLBase +using MINPACK + +function SciMLBase.__solve(prob::Union{NonlinearProblem{uType, iip}, + NonlinearLeastSquaresProblem{uType, iip}}, alg::CMINPACK, args...; + abstol = 1e-6, maxiters = 100000, alias_u0::Bool = false, + kwargs...) where {uType, iip} + if prob.u0 isa Number + u0 = [prob.u0] + else + u0 = NonlinearSolve.__maybe_unaliased(prob.u0, alias_u0) + end + + sizeu = size(prob.u0) + p = prob.p + + # unwrapping alg params + show_trace = alg.show_trace + tracing = alg.tracing + io = alg.io + + if !iip && prob.u0 isa Number + f! = (du, u) -> (du .= prob.f(first(u), p); Cint(0)) + elseif !iip && prob.u0 isa Vector{Float64} + f! = (du, u) -> (du .= prob.f(u, p); Cint(0)) + elseif !iip && prob.u0 isa AbstractArray + f! = (du, u) -> (du .= vec(prob.f(reshape(u, sizeu), p)); Cint(0)) + elseif prob.u0 isa Vector{Float64} + f! = (du, u) -> prob.f(du, u, p) + else # Then it's an in-place function on an abstract array + f! = (du, u) -> (prob.f(reshape(du, sizeu), reshape(u, sizeu), p); du = vec(du); 0) + end + + u = zero(u0) + resid = NonlinearSolve.evaluate_f(prob, u) + m = length(resid) + + method = ifelse(alg.method === :auto, + ifelse(prob isa NonlinearLeastSquaresProblem, :lm, :hydr), alg.method) + + if SciMLBase.has_jac(prob.f) + if !iip && prob.u0 isa Number + g! = (du, u) -> (du .= prob.f.jac(first(u), p); Cint(0)) + elseif !iip && prob.u0 isa Vector{Float64} + g! = (du, u) -> (du .= prob.f.jac(u, p); Cint(0)) + elseif !iip && prob.u0 isa AbstractArray + g! = (du, u) -> (du .= vec(prob.f.jac(reshape(u, sizeu), p)); Cint(0)) + elseif prob.u0 isa Vector{Float64} + g! = (du, u) -> prob.f.jac(du, u, p) + else # Then it's an in-place function on an abstract array + g! = function (du, u) + prob.f.jac(reshape(du, sizeu), reshape(u, sizeu), p) + du = vec(du) + return CInt(0) + end + end + original = MINPACK.fsolve(f!, g!, u0, m; tol = abstol, show_trace, tracing, method, + iterations = maxiters, kwargs...) + else + original = MINPACK.fsolve(f!, u0, m; tol = abstol, show_trace, tracing, method, + iterations = maxiters, kwargs...) + end + + u = reshape(original.x, size(u)) + resid = original.f + retcode = original.converged ? ReturnCode.Success : ReturnCode.Failure + + return SciMLBase.build_solution(prob, alg, u, resid; retcode, original) +end + +end diff --git a/ext/NonlinearSolveNLsolveExt.jl b/ext/NonlinearSolveNLsolveExt.jl new file mode 100644 index 000000000..a9424fdf7 --- /dev/null +++ b/ext/NonlinearSolveNLsolveExt.jl @@ -0,0 +1,3 @@ +module NonlinearSolveNLsolveExt + +end diff --git a/src/NonlinearSolve.jl b/src/NonlinearSolve.jl index 8bdfcedf5..dd0a6cc33 100644 --- a/src/NonlinearSolve.jl +++ b/src/NonlinearSolve.jl @@ -236,7 +236,7 @@ export RadiusUpdateSchemes export NewtonRaphson, TrustRegion, LevenbergMarquardt, DFSane, GaussNewton, PseudoTransient, Broyden, Klement, LimitedMemoryBroyden -export LeastSquaresOptimJL, FastLevenbergMarquardtJL +export LeastSquaresOptimJL, FastLevenbergMarquardtJL, CMINPACK, NLsolveJL export NonlinearSolvePolyAlgorithm, RobustMultiNewton, FastShortcutNonlinearPolyalg, FastShortcutNLLSPolyalg diff --git a/src/default.jl b/src/default.jl index 8b69b5a06..a4b6325ca 100644 --- a/src/default.jl +++ b/src/default.jl @@ -243,32 +243,22 @@ function FastShortcutNonlinearPolyalg(; concrete_jac = nothing, linsolve = nothi prefer_simplenonlinearsolve::Val{SA} = Val(false), autodiff = nothing) where {JAC, SA} if JAC - if SA - algs = (SimpleNewtonRaphson(; - autodiff = ifelse(autodiff === nothing, AutoForwardDiff(), autodiff)), - SimpleTrustRegion(; - autodiff = ifelse(autodiff === nothing, AutoForwardDiff(), autodiff)), - NewtonRaphson(; concrete_jac, linsolve, precs, linesearch = BackTracking(), - autodiff), - TrustRegion(; concrete_jac, linsolve, precs, - radius_update_scheme = RadiusUpdateSchemes.Bastin, autodiff)) - else - algs = (NewtonRaphson(; concrete_jac, linsolve, precs, autodiff), - NewtonRaphson(; concrete_jac, linsolve, precs, linesearch = BackTracking(), - autodiff), - TrustRegion(; concrete_jac, linsolve, precs, autodiff), - TrustRegion(; concrete_jac, linsolve, precs, - radius_update_scheme = RadiusUpdateSchemes.Bastin, autodiff)) - end + algs = (NewtonRaphson(; concrete_jac, linsolve, precs, autodiff), + NewtonRaphson(; concrete_jac, linsolve, precs, linesearch = BackTracking(), + autodiff), + TrustRegion(; concrete_jac, linsolve, precs, autodiff), + TrustRegion(; concrete_jac, linsolve, precs, + radius_update_scheme = RadiusUpdateSchemes.Bastin, autodiff)) else + # SimpleNewtonRaphson and SimpleTrustRegion are not robust to singular Jacobians + # and thus are not included in the polyalgorithm if SA algs = (SimpleBroyden(), Broyden(; init_jacobian = Val(:true_jacobian)), SimpleKlement(), - SimpleNewtonRaphson(; - autodiff = ifelse(autodiff === nothing, AutoForwardDiff(), autodiff)), - SimpleTrustRegion(; - autodiff = ifelse(autodiff === nothing, AutoForwardDiff(), autodiff)), + NewtonRaphson(; concrete_jac, linsolve, precs, autodiff), + NewtonRaphson(; concrete_jac, linsolve, precs, linesearch = BackTracking(), + autodiff), NewtonRaphson(; concrete_jac, linsolve, precs, linesearch = BackTracking(), autodiff), TrustRegion(; concrete_jac, linsolve, precs, diff --git a/src/extension_algs.jl b/src/extension_algs.jl index 3fe4b84fd..d59627a67 100644 --- a/src/extension_algs.jl +++ b/src/extension_algs.jl @@ -80,9 +80,128 @@ function FastLevenbergMarquardtJL(linsolve::Symbol = :cholesky; factor = 1e-6, autodiff isa AutoForwardDiff if Base.get_extension(@__MODULE__, :NonlinearSolveFastLevenbergMarquardtExt) === nothing - error("LeastSquaresOptimJL requires FastLevenbergMarquardt.jl to be loaded") + error("FastLevenbergMarquardtJL requires FastLevenbergMarquardt.jl to be loaded") end return FastLevenbergMarquardtJL{linsolve}(autodiff, factor, factoraccept, factorreject, factorupdate, minscale, maxscale, minfactor, maxfactor) end + +""" + CMINPACK(; show_trace::Bool=false, tracing::Bool=false, method::Symbol=:auto) + +### Keyword Arguments + + - `show_trace`: whether to show the trace. + - `tracing`: who the hell knows what this does. If you find out, please open an issue/PR. + - `method`: the choice of method for the solver. + +### Method Choices + +The keyword argument `method` can take on different value depending on which method of +`fsolve` you are calling. The standard choices of `method` are: + + - `:hybr`: Modified version of Powell's algorithm. Uses MINPACK routine + [`hybrd1`](https://github.com/devernay/cminpack/blob/d1f5f5a273862ca1bbcf58394e4ac060d9e22c76/hybrd1.c) + - `:lm`: Levenberg-Marquardt. Uses MINPACK routine + [`lmdif1`](https://github.com/devernay/cminpack/blob/d1f5f5a273862ca1bbcf58394e4ac060d9e22c76/lmdif1.c) + - `:lmdif`: Advanced Levenberg-Marquardt (more options available with `; kwargs...`). See + MINPACK routine [`lmdif`](https://github.com/devernay/cminpack/blob/d1f5f5a273862ca1bbcf58394e4ac060d9e22c76/lmdif.c) + for more information + - `:hybrd`: Advanced modified version of Powell's algorithm (more options available with + `; kwargs...`). See MINPACK routine + [`hybrd`](https://github.com/devernay/cminpack/blob/d1f5f5a273862ca1bbcf58394e4ac060d9e22c76/hybrd.c) + for more information + +If a Jacobian is supplied as part of the [`NonlinearFunction`](@ref nonlinearfunctions), +then the following methods are allowed: + + - `:hybr`: Advanced modified version of Powell's algorithm with user supplied Jacobian. + Additional arguments are available via `; kwargs...`. See MINPACK routine + [`hybrj`](https://github.com/devernay/cminpack/blob/d1f5f5a273862ca1bbcf58394e4ac060d9e22c76/hybrj.c) + for more information + - `:lm`: Advanced Levenberg-Marquardt with user supplied Jacobian. Additional arguments + are available via `;kwargs...`. See MINPACK routine + [`lmder`](https://github.com/devernay/cminpack/blob/d1f5f5a273862ca1bbcf58394e4ac060d9e22c76/lmder.c) + for more information + +The default choice of `:auto` selects `:hybr` for NonlinearProblem and `:lm` for +NonlinearLeastSquaresProblem. +""" +struct CMINPACK <: AbstractNonlinearAlgorithm + show_trace::Bool + tracing::Bool + method::Symbol +end + +function CMINPACK(; show_trace::Bool = false, tracing::Bool = false, method::Symbol = :auto) + if Base.get_extension(@__MODULE__, :NonlinearSolveMINPACKExt) === nothing + error("CMINPACK requires MINPACK.jl to be loaded") + end + + return CMINPACK(show_trace, tracing, method) +end + +""" + NLSolveJL(; method=:trust_region, autodiff=:central, store_trace=false, + extended_trace=false, linesearch=LineSearches.Static(), + linsolve=(x, A, b) -> copyto!(x, A\\b), factor = one(Float64), autoscale=true, + m=10, beta=one(Float64), show_trace=false) + +### Keyword Arguments + + - `method`: the choice of method for solving the nonlinear system. + - `autodiff`: the choice of method for generating the Jacobian. Defaults to `:central` or + central differencing via FiniteDiff.jl. The other choices are `:forward` + - `show_trace`: should a trace of the optimization algorithm's state be shown on `STDOUT`? + - `extended_trace`: should additional algorithm internals be added to the state trace? + - `linesearch`: the line search method to be used within the solver method. The choices + are line search types from + [LineSearches.jl](https://github.com/JuliaNLSolvers/LineSearches.jl). + - `linsolve`: a function `linsolve(x, A, b)` that solves `Ax = b`. + - `factor``: determines the size of the initial trust region. This size is set to the product of factor and the euclidean norm of `u0` if nonzero, or else to factor itself. + - `autoscale`: if true, then the variables will be automatically rescaled. The scaling + factors are the norms of the Jacobian columns. + - `m`: the amount of history in the Anderson method. Naive "Picard"-style iteration can be + achieved by setting m=0, but that isn't advisable for contractions whose Lipschitz + constants are close to 1. If convergence fails, though, you may consider lowering it. + - `beta`: It is also known as DIIS or Pulay mixing, this method is based on the acceleration + of the fixed-point iteration xₙ₊₁ = xₙ + beta*f(xₙ), where by default beta = 1. + - `store_trace``: should a trace of the optimization algorithm's state be stored? + +### Submethod Choice + +Choices for methods in `NLSolveJL`: + + - `:anderson`: Anderson-accelerated fixed-point iteration + - `:broyden`: Broyden's quasi-Newton method + - `:newton`: Classical Newton method with an optional line search + - `:trust_region`: Trust region Newton method (the default choice) For more information on + these arguments, consult the + [NLsolve.jl documentation](https://github.com/JuliaNLSolvers/NLsolve.jl). +""" +@concrete struct NLSolveJL <: AbstractNonlinearAlgorithm + method::Symbol + autodiff::Symbol + store_trace::Bool + extended_trace::Bool + linesearch + linsolve + factor + autoscale::Bool + m::Int + beta + show_trace::Bool +end + +function NLSolveJL(; method = :trust_region, autodiff = :central, store_trace = false, + extended_trace = false, linesearch = LineSearches.Static(), + linsolve = (x, A, b) -> copyto!(x, A \ b), factor = 1.0, autoscale = true, m = 10, + beta = one(Float64), show_trace = false) + if Base.get_extension(@__MODULE__, :NonlinearSolveNLsolveExt) === nothing + error("NLSolveJL requires NLsolve.jl to be loaded") + end + + return NLSolveJL(method, autodiff, store_trace, extended_trace, linesearch, linsolve, + factor, autoscale, m, beta, show_trace) +end diff --git a/test/minpack.jl b/test/minpack.jl new file mode 100644 index 000000000..46cc6be6e --- /dev/null +++ b/test/minpack.jl @@ -0,0 +1,35 @@ +using NonlinearSolve, MINPACK, Test + +function f_iip(du, u, p) + du[1] = 2 - 2u[1] + du[2] = u[1] - 4u[2] +end +u0 = zeros(2) +prob_iip = NonlinearProblem{true}(f_iip, u0) +abstol = 1e-8 + +for alg in [CMINPACK()] + local sol + sol = solve(prob_iip, alg) + @test sol.retcode == ReturnCode.Success + p = nothing + + du = zeros(2) + f_iip(du, sol.u, nothing) + @test maximum(du) < 1e-6 +end + +# OOP Tests +f_oop(u, p) = [2 - 2u[1], u[1] - 4u[2]] +u0 = zeros(2) +prob_oop = NonlinearProblem{false}(f_oop, u0) + +for alg in [CMINPACK()] + local sol + sol = solve(prob_oop, alg) + @test sol.retcode == ReturnCode.Success + + du = zeros(2) + du = f_oop(sol.u, nothing) + @test maximum(du) < 1e-6 +end diff --git a/test/nlsolve.jl b/test/nlsolve.jl new file mode 100644 index 000000000..e1d7714d0 --- /dev/null +++ b/test/nlsolve.jl @@ -0,0 +1,138 @@ +using NonlinearSolve, NLsolve, LinearAlgebra, Test + +# IIP Tests +function f_iip(du, u, p, t) + du[1] = 2 - 2u[1] + du[2] = u[1] - 4u[2] +end +u0 = zeros(2) +prob_iip = SteadyStateProblem(f_iip, u0) +abstol = 1e-8 + +for alg in [NLSolveJL()] + sol = solve(prob_iip, alg) + @test sol.retcode == ReturnCode.Success + p = nothing + + du = zeros(2) + f_iip(du, sol.u, nothing, 0) + @test maximum(du) < 1e-6 +end + +# OOP Tests +f_oop(u, p, t) = [2 - 2u[1], u[1] - 4u[2]] +u0 = zeros(2) +prob_oop = SteadyStateProblem(f_oop, u0) + +for alg in [NLSolveJL()] + sol = solve(prob_oop, alg) + @test sol.retcode == ReturnCode.Success + # test the solver is doing reasonable things for linear solve + # and that the stats are working properly + @test 1 <= sol.stats.nf < 10 + + du = zeros(2) + du = f_oop(sol.u, nothing, 0) + @test maximum(du) < 1e-6 +end + +# NonlinearProblem Tests + +function f_iip(du, u, p) + du[1] = 2 - 2u[1] + du[2] = u[1] - 4u[2] +end +u0 = zeros(2) +prob_iip = NonlinearProblem{true}(f_iip, u0) +abstol = 1e-8 +for alg in [NLSolveJL()] + local sol + sol = solve(prob_iip, alg) + @test sol.retcode == ReturnCode.Success + p = nothing + + du = zeros(2) + f_iip(du, sol.u, nothing) + @test maximum(du) < 1e-6 +end + +# OOP Tests +f_oop(u, p) = [2 - 2u[1], u[1] - 4u[2]] +u0 = zeros(2) +prob_oop = NonlinearProblem{false}(f_oop, u0) +for alg in [NLSolveJL()] + local sol + sol = solve(prob_oop, alg) + @test sol.retcode == ReturnCode.Success + + du = zeros(2) + du = f_oop(sol.u, nothing) + @test maximum(du) < 1e-6 +end + +# tolerance tests +f_tol(u, p) = u^2 - 2 +prob_tol = NonlinearProblem(f_tol, 1.0) +for tol in [1e-1, 1e-3, 1e-6, 1e-10, 1e-15] + sol = solve(prob_tol, NLSolveJL(), abstol = tol) + @test abs(sol.u[1] - sqrt(2)) < tol +end + +# Test the finite differencing technique +function f!(fvec, x, p) + fvec[1] = (x[1] + 3) * (x[2]^3 - 7) + 18 + fvec[2] = sin(x[2] * exp(x[1]) - 1) +end + +prob = NonlinearProblem{true}(f!, [0.1; 1.2]) +sol = solve(prob, NLSolveJL(autodiff = :central)) + +du = zeros(2) +f!(du, sol.u, nothing) +@test maximum(du) < 1e-6 + +# Test the autodiff technique +function f!(fvec, x, p) + fvec[1] = (x[1] + 3) * (x[2]^3 - 7) + 18 + fvec[2] = sin(x[2] * exp(x[1]) - 1) +end + +prob = NonlinearProblem{true}(f!, [0.1; 1.2]) +sol = solve(prob, NLSolveJL(autodiff = :forward)) + +du = zeros(2) +f!(du, sol.u, nothing) +@test maximum(du) < 1e-6 + +function problem(x, A) + return x .^ 2 - A +end + +function problemJacobian(x, A) + return diagm(2 .* x) +end + +function f!(F, u, p) + F[1:152] = problem(u, p) +end + +function j!(J, u, p) + J[1:152, 1:152] = problemJacobian(u, p) +end + +f = NonlinearFunction(f!) + +init = ones(152); +A = ones(152); +A[6] = 0.8 + +f = NonlinearFunction(f!, jac = j!) + +p = A + +ProbN = NonlinearProblem(f, init, p) +sol = solve(ProbN, NLSolveJL(), reltol = 1e-8, abstol = 1e-8) + +init = ones(Complex{Float64}, 152); +ProbN = NonlinearProblem(f, init, p) +sol = solve(ProbN, NLSolveJL(), reltol = 1e-8, abstol = 1e-8) diff --git a/test/nonlinear_least_squares.jl b/test/nonlinear_least_squares.jl index 330c2f5da..f4b8233e0 100644 --- a/test/nonlinear_least_squares.jl +++ b/test/nonlinear_least_squares.jl @@ -1,6 +1,6 @@ using NonlinearSolve, LinearSolve, LinearAlgebra, Test, StableRNGs, Random, ForwardDiff, Zygote -import FastLevenbergMarquardt, LeastSquaresOptim +import FastLevenbergMarquardt, LeastSquaresOptim, MINPACK true_function(x, θ) = @. θ[1] * exp(θ[2] * x) * cos(θ[3] * x + θ[4]) true_function(y, x, θ) = (@. y = θ[1] * exp(θ[2] * x) * cos(θ[3] * x + θ[4])) @@ -99,7 +99,8 @@ probs = [ NonlinearLeastSquaresProblem(NonlinearFunction{false}(loss_function; jac), θ_init, x), ] -solvers = [FastLevenbergMarquardtJL(linsolve) for linsolve in (:cholesky, :qr)] +solvers = Any[FastLevenbergMarquardtJL(linsolve) for linsolve in (:cholesky, :qr)] +push!(solvers, CMINPACK()) for solver in solvers, prob in probs @time sol = solve(prob, solver; maxiters = 10000, abstol = 1e-8) @@ -114,8 +115,10 @@ probs = [ NonlinearLeastSquaresProblem(NonlinearFunction{false}(loss_function), θ_init, x), ] -solvers = [FastLevenbergMarquardtJL(linsolve; autodiff) for linsolve in (:cholesky, :qr), -autodiff in (nothing, AutoForwardDiff(), AutoFiniteDiff())] +solvers = vec(Any[FastLevenbergMarquardtJL(linsolve; autodiff) + for linsolve in (:cholesky, :qr), +autodiff in (nothing, AutoForwardDiff(), AutoFiniteDiff())]) +append!(solvers, [CMINPACK(; method) for method in (:auto, :lm, :lmdif)]) for solver in solvers, prob in probs @time sol = solve(prob, solver; maxiters = 10000, abstol = 1e-8) diff --git a/test/runtests.jl b/test/runtests.jl index 465f9d7b1..855a793fc 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -14,12 +14,14 @@ end @time begin if GROUP == "All" || GROUP == "Core" @time @safetestset "Quality Assurance" include("qa.jl") - @time @safetestset "Nonlinear Least Squares" include("nonlinear_least_squares.jl") @time @safetestset "Basic Tests + Some AD" include("basictests.jl") @time @safetestset "Sparsity Tests" include("sparse.jl") @time @safetestset "Polyalgs" include("polyalgs.jl") @time @safetestset "Matrix Resizing" include("matrix_resizing.jl") @time @safetestset "Infeasible Problems" include("infeasible.jl") + @time @safetestset "Nonlinear Least Squares" include("nonlinear_least_squares.jl") + @time @safetestset "MINPACK" include("minpack.jl") + @time @safetestset "NLsolve" include("nlsolve.jl") end if GROUP == "All" || GROUP == "23TestProblems"