SciML · ChrisRackauckas · Nov 4, 2023 · Nov 4, 2023 · Nov 4, 2023
diff --git a/Project.toml b/Project.toml
@@ -29,11 +29,15 @@ 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"
 
 [extensions]
 NonlinearSolveBandedMatricesExt = "BandedMatrices"
 NonlinearSolveFastLevenbergMarquardtExt = "FastLevenbergMarquardt"
 NonlinearSolveLeastSquaresOptimExt = "LeastSquaresOptim"
+NonlinearSolveMINPACKExt = "MINPACK"
+NonlinearSolveNLsolveExt = "NLsolve"
 
 [compat]
 ADTypes = "0.2"
@@ -74,7 +78,9 @@ 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"
 NaNMath = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
+NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
 NonlinearProblemLibrary = "b7050fa9-e91f-4b37-bcee-a89a063da141"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
@@ -86,4 +92,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [targets]
-test = ["Enzyme", "BenchmarkTools", "SafeTestsets", "Pkg", "Test", "ForwardDiff", "StaticArrays", "Symbolics", "LinearSolve", "Random", "LinearAlgebra", "Zygote", "SparseDiffTools", "NonlinearProblemLibrary", "LeastSquaresOptim", "FastLevenbergMarquardt", "NaNMath", "BandedMatrices", "DiffEqBase"]
+test = ["NLsolve", "Enzyme", "BenchmarkTools", "SafeTestsets", "Pkg", "Test", "ForwardDiff", "StaticArrays", "Symbolics", "LinearSolve", "Random", "LinearAlgebra", "Zygote", "SparseDiffTools", "NonlinearProblemLibrary", "LeastSquaresOptim", "FastLevenbergMarquardt", "NaNMath", "BandedMatrices", "DiffEqBase", "MINPACK"]
diff --git a/docs/Project.toml b/docs/Project.toml
@@ -4,11 +4,13 @@ ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+FastLevenbergMarquardt = "7a0df574-e128-4d35-8cbd-3d84502bf7ce"
 IncompleteLU = "40713840-3770-5561-ab4c-a76e7d0d7895"
+LeastSquaresOptim = "0fc2ff8b-aaa3-5acd-a817-1944a5e08891"
 LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
-NonlinearSolveMINPACK = "c100e077-885d-495a-a2ea-599e143bf69d"
+MINPACK = "4854310b-de5a-5eb6-a2a5-c1dee2bd17f9"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 SciMLNLSolve = "e9a6253c-8580-4d32-9898-8661bb511710"
 SimpleNonlinearSolve = "727e6d20-b764-4bd8-a329-72de5adea6c7"
@@ -27,7 +29,6 @@ IncompleteLU = "0.2"
 LinearSolve = "2"
 ModelingToolkit = "8"
 NonlinearSolve = "1, 2"
-NonlinearSolveMINPACK = "0.1"
 SciMLBase = "2.4"
 SciMLNLSolve = "0.1"
 SimpleNonlinearSolve = "0.1.5"

diff --git a/docs/pages.jl b/docs/pages.jl
@@ -20,6 +20,8 @@ pages = ["index.md",
         "solvers/LineSearch.md"],
     "Detailed Solver APIs" => Any["api/nonlinearsolve.md",
         "api/simplenonlinearsolve.md",
+        "api/fastlevenbergmarquardt.md",
+        "api/leastsquaresoptim.md",
         "api/minpack.md",
         "api/nlsolve.md",
         "api/sundials.md",

diff --git a/docs/src/api/fastlevenbergmarquardt.md b/docs/src/api/fastlevenbergmarquardt.md
@@ -0,0 +1,19 @@
+# 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:
+
+```julia
+using Pkg
+Pkg.add("FastLevenbergMarquardt")
+using FastLevenbergMarquardt
+```
+
+These methods can be used independently of the rest of NonlinearSolve.jl
+
+## Solver API
+
+```@docs
+FastLevenbergMarquardtJL
+```
diff --git a/docs/src/api/leastsquaresoptim.md b/docs/src/api/leastsquaresoptim.md
@@ -0,0 +1,19 @@
+# LeastSquaresOptim.jl
+
+This is a wrapper package 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
+```
+
+These methods can be used independently of the rest of NonlinearSolve.jl
+
+## Solver API
+
+```@docs
+LeastSquaresOptimJL
+```
diff --git a/docs/src/api/minpack.md b/docs/src/api/minpack.md
@@ -1,13 +1,13 @@
 # MINPACK.jl
 
-This is a wrapper package for importing solvers from Sundials into the SciML interface.
+This is a wrapper package 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
 ```
 
 These methods can be used independently of the rest of NonlinearSolve.jl

diff --git a/docs/src/api/nlsolve.md b/docs/src/api/nlsolve.md
@@ -6,8 +6,8 @@ the package before using these solvers:
 
 ```julia
 using Pkg
-Pkg.add("SciMLNLSolve")
-using SciMLNLSolve
+Pkg.add("NLsolve")
+using NLsolve
 ```
 
 ## Solver API

diff --git a/docs/src/solvers/NonlinearLeastSquaresSolvers.md b/docs/src/solvers/NonlinearLeastSquaresSolvers.md
@@ -14,6 +14,8 @@ algorithm (`LevenbergMarquardt`).
 
 ## Full List of Methods
 
+### NonlinearSolve.jl
+
   - `LevenbergMarquardt()`: An advanced Levenberg-Marquardt implementation with the
     improvements suggested in the [paper](https://arxiv.org/abs/1201.5885) "Improvements to
     the Levenberg-Marquardt algorithm for nonlinear least-squares minimization". Designed for
@@ -22,5 +24,56 @@ algorithm (`LevenbergMarquardt`).
     handling of sparse matrices via colored automatic differentiation and preconditioned
     linear solvers. Designed for large-scale and numerically-difficult nonlinear least squares
     problems.
-  - `SimpleNewtonRaphson()`: Simple Gauss Newton Implementation with `QRFactorization` to
-    solve a linear least squares problem at each step!
+
+### SimpleNonlinearSolve.jl
+
+These methods are included with NonlinearSolve.jl by default, though SimpleNonlinearSolve.jl can be used 
+directly to reduce dependencies and improve load times. SimpleNonlinearSolve.jl's methods excel at small 
+problems and problems defined with static arrays.
+
+  - `SimpleGaussNewton()`: Simple Gauss Newton implementation using QR factorizations for numerical stability.
+
+### FastLevenbergMarquardt.jl
+
+A wrapper over [FastLevenbergMarquardt.jl](https://github.com/kamesy/FastLevenbergMarquardt.jl). Note that
+it is called `FastLevenbergMarquardt` since the original package is called "Fast", though benchmarks
+demonstrate `LevenbergMarquardt()` usually outperforms.
+
+  - `FastLevenbergMarquardtJL(linsolve = :cholesky)`, can also choose `linsolve = :qr`.
+
+### LeastSquaresOptim.jl
+
+A wrapper over [LeastSquaresOptim.jl](https://github.com/matthieugomez/LeastSquaresOptim.jl).
+Has a core algorithm `LeastSquaresOptimJL(alg; linsolve)` where the choices for `alg` are:
+
+- `:lm` a Levenberg-Marquardt implementation
+- `:dogleg` a trust-region dogleg Gauss-Newton
+
+And the choices for `linsolve` are:
+
+- `:qr`
+- `: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 with any solver of
+[Optimization.jl](https://github.com/SciML/Optimization.jl).
diff --git a/docs/src/solvers/NonlinearSystemSolvers.md b/docs/src/solvers/NonlinearSystemSolvers.md
@@ -113,7 +113,7 @@ computationally expensive than direct methods.
   - `DynamicSS()` : Uses an ODE solver to find the steady state. Automatically
     terminates when close to the steady state.
 
-### SciMLNLSolve.jl
+### NLsolve.jl
 
 This is a wrapper package for importing solvers from NLsolve.jl into the SciML interface.
 
@@ -127,8 +127,11 @@ Submethod choices for this algorithm include:
 
 ### MINPACK.jl
 
-MINPACK.jl methods are good for medium-sized nonlinear solves. It does not scale due to
-the lack of sparse Jacobian support, though the methods are very robust and stable.
+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)
 

diff --git a/ext/NonlinearSolveMINPACKExt.jl b/ext/NonlinearSolveMINPACKExt.jl
@@ -0,0 +1,82 @@
+module NonlinearSolveMINPACKExt
+
+using NonlinearSolve, SciMLBase
+using MINPACK
+
+function SciMLBase.solve(prob::Union{SciMLBase.NonlinearProblem{uType, isinplace},
+                                     SciMLBase.NonlinearLeastSquaresProblem{uType, isinplace}},
+                         alg::CMINPACK,
+                         reltol = 1e-3,
+                         abstol = 1e-6,
+                         maxiters = 100000,
+                         timeseries = [],
+                         ts = [],
+                         ks = [], ;
+                         kwargs...) where {uType, isinplace}
+    if prob.u0 isa Number
+        u0 = [prob.u0]
+    else
+        u0 = deepcopy(prob.u0)
+    end
+
+    sizeu = size(prob.u0)
+    p = prob.p
+
+    # unwrapping alg params
+    method = alg.method
+    show_trace = alg.show_trace
+    tracing = alg.tracing
+    io = alg.io
+
+    if !isinplace && prob.u0 isa Number
+        f! = (du, u) -> (du .= prob.f(first(u), p); Cint(0))
+    elseif !isinplace && prob.u0 isa Vector{Float64}
+        f! = (du, u) -> (du .= prob.f(u, p); Cint(0))
+    elseif !isinplace && 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 = similar(u0)
+
+    m = prob.f.resid_prototype === nothing ? length(u0) : length(prob.f.resid_prototype)
+
+    if SciMLBase.has_jac(prob.f)
+        if !isinplace && prob.u0 isa Number
+            g! = (du, u) -> (du .= prob.jac(first(u), p); Cint(0))
+        elseif !isinplace && prob.u0 isa Vector{Float64}
+            g! = (du, u) -> (du .= prob.jac(u, p); Cint(0))
+        elseif !isinplace && prob.u0 isa AbstractArray
+            g! = (du, u) -> (du .= vec(prob.jac(reshape(u, sizeu), p)); Cint(0))
+        elseif prob.u0 isa Vector{Float64}
+            g! = (du, u) -> prob.jac(du, u, p)
+        else # Then it's an in-place function on an abstract array
+            g! = (du, u) -> (prob.jac(reshape(du, sizeu), reshape(u, sizeu), p);
+                             du = vec(du);
+                             0)
+        end
+        original = MINPACK.fsolve(f!, g!, u0, m;
+                                  tol = abstol,
+                                  show_trace, tracing, method,
+                                  iterations = maxiters, io, kwargs...)
+    else
+        original = MINPACK.fsolve(f!, u0, m;
+                                  tol = abstol,
+                                  show_trace, tracing, method,
+                                  iterations = maxiters, io, kwargs...)
+    end
+
+    u = reshape(original.x, size(u))
+    resid = original.f
+    retcode = original.converged ? ReturnCode.Success : ReturnCode.Failure
+    SciMLBase.build_solution(prob, alg, u, resid; retcode = retcode,
+                             original = original)
+end
+
+end