Finish wrapping NLSolvers.jl

Author: avik-pal

docs/src/api/nlsolve.md

@@ -7,7 +7,7 @@ using these solvers:
 ```julia
 using Pkg
 Pkg.add("NLsolve")
-using NLSolve, NonlinearSolve
+using NLsolve, NonlinearSolve
 ```
 
 ## Solver API

docs/src/api/nlsolvers.md

@@ -0,0 +1,17 @@
+# NLSolvers.jl
+
+This is a extension for importing solvers from NLSolvers.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("NLSolvers")
+using NLSolvers, NonlinearSolve
+```
+
+## Solver API
+
+```@docs
+NLSolversJL
+```

docs/src/solvers/nonlinear_system_solvers.md

@@ -168,3 +168,12 @@ SIAMFANLEquations.jl is a wrapper for the methods in the SIAMFANLEquations.jl li
 Other solvers listed in [Fixed Point Solvers](@ref fixed_point_methods_full_list),
 [FastLevenbergMarquardt.jl](@ref fastlm_wrapper_summary) and
 [LeastSquaresOptim.jl](@ref lso_wrapper_summary) can also solve nonlinear systems.
+
+### NLSolvers.jl
+
+This is a wrapper package for importing solvers from NLSolvers.jl into the SciML interface.
+
+  - [`NLSolversJL()`](@ref): A wrapper for
+    [NLSolvers.jl](https://github.com/JuliaNLSolvers/NLSolvers.jl)
+
+For a list of possible solvers see the [NLSolvers.jl documentation](https://julianlsolvers.github.io/NLSolvers.jl/)

ext/NonlinearSolveNLSolversExt.jl

@@ -0,0 +1,84 @@
+module NonlinearSolveNLSolversExt
+
+using ADTypes, FastClosures, NonlinearSolve, NLSolvers, SciMLBase, LinearAlgebra
+using FiniteDiff, ForwardDiff
+
+function SciMLBase.__solve(prob::NonlinearProblem, alg::NLSolversJL, args...;
+        abstol = nothing, reltol = nothing, maxiters = 1000, alias_u0::Bool = false,
+        termination_condition = nothing, kwargs...)
+    NonlinearSolve.__test_termination_condition(termination_condition, :NLSolversJL)
+
+    abstol = NonlinearSolve.DEFAULT_TOLERANCE(abstol, eltype(prob.u0))
+    reltol = NonlinearSolve.DEFAULT_TOLERANCE(reltol, eltype(prob.u0))
+
+    options = NEqOptions(; maxiter = maxiters, f_abstol = abstol, f_reltol = reltol)
+
+    if prob.u0 isa Number
+        f_scalar = @closure x -> prob.f(x, prob.p)
+
+        if alg.autodiff === nothing
+            if ForwardDiff.can_dual(typeof(prob.u0))
+                autodiff_concrete = :forwarddiff
+            else
+                autodiff_concrete = :finitediff
+            end
+        else
+            if alg.autodiff isa AutoForwardDiff || alg.autodiff isa AutoPolyesterForwardDiff
+                autodiff_concrete = :forwarddiff
+            elseif alg.autodiff isa AutoFiniteDiff
+                autodiff_concrete = :finitediff
+            else
+                error("Only ForwardDiff or FiniteDiff autodiff is supported.")
+            end
+        end
+
+        if autodiff_concrete === :forwarddiff
+            fj_scalar = @closure (Jx, x) -> begin
+                T = typeof(NonlinearSolve.NonlinearSolveTag())
+                x_dual = ForwardDiff.Dual{T}(x, one(x))
+                y = prob.f(x_dual, prob.p)
+                return ForwardDiff.value(y), ForwardDiff.extract_derivative(T, y)
+            end
+        else
+            fj_scalar = @closure (Jx, x) -> begin
+                _f = Base.Fix2(prob.f, prob.p)
+                return _f(x), FiniteDiff.finite_difference_derivative(_f, x)
+            end
+        end
+
+        prob_obj = NLSolvers.ScalarObjective(; f = f_scalar, fg = fj_scalar)
+        prob_nlsolver = NEqProblem(prob_obj; inplace = false)
+        res = NLSolvers.solve(prob_nlsolver, prob.u0, alg.method, options)
+
+        retcode = ifelse(norm(res.info.best_residual, Inf) ≤ abstol, ReturnCode.Success,
+            ReturnCode.MaxIters)
+        stats = SciMLBase.NLStats(-1, -1, -1, -1, res.info.iter)
+
+        return SciMLBase.build_solution(prob, alg, res.info.solution,
+            res.info.best_residual; retcode, original = res, stats)
+    end
+
+    f!, u0, resid = NonlinearSolve.__construct_extension_f(prob; alias_u0)
+
+    jac! = NonlinearSolve.__construct_extension_jac(prob, alg, u0, resid; alg.autodiff)
+
+    FJ_vector! = @closure (Fx, Jx, x) -> begin
+        f!(Fx, x)
+        jac!(Jx, x)
+        return Fx, Jx
+    end
+
+    prob_obj = NLSolvers.VectorObjective(; F = f!, FJ = FJ_vector!)
+    prob_nlsolver = NEqProblem(prob_obj)
+
+    res = NLSolvers.solve(prob_nlsolver, u0, alg.method, options)
+
+    retcode = ifelse(norm(res.info.best_residual, Inf) ≤ abstol, ReturnCode.Success,
+        ReturnCode.MaxIters)
+    stats = SciMLBase.NLStats(-1, -1, -1, -1, res.info.iter)
+
+    return SciMLBase.build_solution(prob, alg, res.info.solution,
+        res.info.best_residual; retcode, original = res, stats)
+end
+
+end

ext/NonlinearSolveNLsolversExt.jl

@@ -281,12 +281,28 @@ function NLsolveJL(; method = :trust_region, autodiff = :central, store_trace =
         factor, autoscale, m, beta, show_trace)
 end
 
-@concrete struct NLSolversJL <: AbstractNonlinearSolveExtensionAlgorithm
-    method
-    autodiff
+"""
+    NLSolversJL(method; autodiff = nothing)
+    NLSolversJL(; method, autodiff = nothing)
+
+Wrapper over NLSolvers.jl Nonlinear Equation Solvers. We automatically construct the
+jacobian function and supply it to the solver.
+
+### Arguments
+
+  - `method`: the choice of method for solving the nonlinear system. See the documentation
+    for NLSolvers.jl for more information.
+  - `autodiff`: the choice of method for generating the Jacobian. Defaults to `nothing`
+    which means that a default is selected according to the problem specification. Can be
+    any valid ADTypes.jl autodiff type (conditional on that backend being supported in
+    NonlinearSolve.jl).
+"""
+struct NLSolversJL{M, AD} <: AbstractNonlinearSolveExtensionAlgorithm
+    method::M
+    autodiff::AD
 
     function NLSolversJL(method, autodiff)
-        if Base.get_extension(@__MODULE__, :NonlinearSolveNLsolversExt) === nothing
+        if Base.get_extension(@__MODULE__, :NonlinearSolveNLSolversExt) === nothing
             error("NLSolversJL requires NLSolvers.jl to be loaded")
         end
 

src/algorithms/extension_algs.jl

@@ -1,6 +1,9 @@
 @testsetup module WrapperRootfindImports
 using Reexport
-@reexport using LinearAlgebra, NLsolve, SIAMFANLEquations, MINPACK
+@reexport using LinearAlgebra
+import NLSolvers, NLsolve, SIAMFANLEquations, MINPACK
+
+export NLSolvers
 end
 
 @testitem "Steady State Problems" setup=[WrapperRootfindImports] begin
@@ -12,7 +15,8 @@ end
     u0 = zeros(2)
     prob_iip = SteadyStateProblem(f_iip, u0)
 
-    for alg in [NLsolveJL(), CMINPACK(), SIAMFANLEquationsJL()]
+    for alg in [
+        NLSolversJL(NLSolvers.LineSearch(NLSolvers.Newton(), NLSolvers.Backtracking())), NLsolveJL(), CMINPACK(), SIAMFANLEquationsJL()]
         sol = solve(prob_iip, alg)
         @test SciMLBase.successful_retcode(sol.retcode)
         @test maximum(abs, sol.resid) < 1e-6
@@ -23,7 +27,8 @@ end
     u0 = zeros(2)
     prob_oop = SteadyStateProblem(f_oop, u0)
 
-    for alg in [NLsolveJL(), CMINPACK(), SIAMFANLEquationsJL()]
+    for alg in [
+        NLSolversJL(NLSolvers.LineSearch(NLSolvers.Newton(), NLSolvers.Backtracking())), NLsolveJL(), CMINPACK(), SIAMFANLEquationsJL()]
         sol = solve(prob_oop, alg)
         @test SciMLBase.successful_retcode(sol.retcode)
         @test maximum(abs, sol.resid) < 1e-6
@@ -39,7 +44,8 @@ end
     u0 = zeros(2)
     prob_iip = NonlinearProblem{true}(f_iip, u0)
 
-    for alg in [NLsolveJL(), CMINPACK(), SIAMFANLEquationsJL()]
+    for alg in [
+        NLSolversJL(NLSolvers.LineSearch(NLSolvers.Newton(), NLSolvers.Backtracking())), NLsolveJL(), CMINPACK(), SIAMFANLEquationsJL()]
         local sol
         sol = solve(prob_iip, alg)
         @test SciMLBase.successful_retcode(sol.retcode)
@@ -50,7 +56,8 @@ end
     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(), CMINPACK(), SIAMFANLEquationsJL()]
+    for alg in [
+        NLSolversJL(NLSolvers.LineSearch(NLSolvers.Newton(), NLSolvers.Backtracking())), NLsolveJL(), CMINPACK(), SIAMFANLEquationsJL()]
         local sol
         sol = solve(prob_oop, alg)
         @test SciMLBase.successful_retcode(sol.retcode)
@@ -61,7 +68,10 @@ end
     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],
-        alg in [NLsolveJL(), CMINPACK(), SIAMFANLEquationsJL(; method = :newton),
+        alg in [
+            NLSolversJL(NLSolvers.LineSearch(NLSolvers.Newton(), NLSolvers.Backtracking())),
+            NLsolveJL(),
+            CMINPACK(), SIAMFANLEquationsJL(; method = :newton),
             SIAMFANLEquationsJL(; method = :pseudotransient),
             SIAMFANLEquationsJL(; method = :secant)]
 
@@ -107,6 +117,10 @@ end
 
     sol = solve(ProbN, NLsolveJL(); abstol = 1e-8)
     @test maximum(abs, sol.resid) < 1e-6
+    sol = solve(ProbN,
+        NLSolversJL(NLSolvers.LineSearch(NLSolvers.Newton(), NLSolvers.Backtracking()));
+        abstol = 1e-8)
+    @test maximum(abs, sol.resid) < 1e-6
     sol = solve(ProbN, SIAMFANLEquationsJL(; method = :newton); abstol = 1e-8)
     @test maximum(abs, sol.resid) < 1e-6
     sol = solve(ProbN, SIAMFANLEquationsJL(; method = :pseudotransient); abstol = 1e-8)