Skip to content

Commit da4cdfc

Browse files
ChrisRackauckasChrisRackauckas
committed
Add MINPACK extension for nlls
1 parent 07d4e25 commit da4cdfc

File tree

11 files changed

+225
-6
lines changed

11 files changed

+225
-6
lines changed

‎Project.toml‎

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -29,10 +29,12 @@ UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
2929
BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
3030
FastLevenbergMarquardt = "7a0df574-e128-4d35-8cbd-3d84502bf7ce"
3131
LeastSquaresOptim = "0fc2ff8b-aaa3-5acd-a817-1944a5e08891"
32+
MINPACK = "4854310b-de5a-5eb6-a2a5-c1dee2bd17f9"
3233

3334
[extensions]
3435
NonlinearSolveBandedMatricesExt = "BandedMatrices"
3536
NonlinearSolveFastLevenbergMarquardtExt = "FastLevenbergMarquardt"
37+
NonlinearSolveMINPACKExt = "MINPACK"
3638
NonlinearSolveLeastSquaresOptimExt = "LeastSquaresOptim"
3739

3840
[compat]
@@ -74,6 +76,7 @@ ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
7476
LeastSquaresOptim = "0fc2ff8b-aaa3-5acd-a817-1944a5e08891"
7577
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
7678
LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
79+
MINPACK = "4854310b-de5a-5eb6-a2a5-c1dee2bd17f9"
7780
NaNMath = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
7881
NonlinearProblemLibrary = "b7050fa9-e91f-4b37-bcee-a89a063da141"
7982
Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
@@ -86,4 +89,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
8689
Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
8790

8891
[targets]
89-
test = ["Enzyme", "BenchmarkTools", "SafeTestsets", "Pkg", "Test", "ForwardDiff", "StaticArrays", "Symbolics", "LinearSolve", "Random", "LinearAlgebra", "Zygote", "SparseDiffTools", "NonlinearProblemLibrary", "LeastSquaresOptim", "FastLevenbergMarquardt", "NaNMath", "BandedMatrices", "DiffEqBase"]
92+
test = ["Enzyme", "BenchmarkTools", "SafeTestsets", "Pkg", "Test", "ForwardDiff", "StaticArrays", "Symbolics", "LinearSolve", "Random", "LinearAlgebra", "Zygote", "SparseDiffTools", "NonlinearProblemLibrary", "LeastSquaresOptim", "FastLevenbergMarquardt", "NaNMath", "BandedMatrices", "DiffEqBase", "MINPACK"]

‎docs/Project.toml‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -4,11 +4,13 @@ ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
44
BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
55
DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
66
Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
7+
FastLevenbergMarquardt = "7a0df574-e128-4d35-8cbd-3d84502bf7ce"
78
IncompleteLU = "40713840-3770-5561-ab4c-a76e7d0d7895"
9+
LeastSquaresOptim = "0fc2ff8b-aaa3-5acd-a817-1944a5e08891"
810
LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
911
ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
1012
NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
11-
NonlinearSolveMINPACK = "c100e077-885d-495a-a2ea-599e143bf69d"
13+
MINPACK = "4854310b-de5a-5eb6-a2a5-c1dee2bd17f9"
1214
SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
1315
SciMLNLSolve = "e9a6253c-8580-4d32-9898-8661bb511710"
1416
SimpleNonlinearSolve = "727e6d20-b764-4bd8-a329-72de5adea6c7"

‎docs/pages.jl‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -20,6 +20,8 @@ pages = ["index.md",
2020
"solvers/LineSearch.md"],
2121
"Detailed Solver APIs" => Any["api/nonlinearsolve.md",
2222
"api/simplenonlinearsolve.md",
23+
"api/fastlevenbergmarquardt.md",
24+
"api/leastsquaresoptim.md",
2325
"api/minpack.md",
2426
"api/nlsolve.md",
2527
"api/sundials.md",

‎docs/src/api/fastlevenbergmarquardt.md‎

Lines changed: 19 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,19 @@
1+
# FastLevenbergMarquardt.jl
2+
3+
This is a wrapper package for importing solvers from FastLevenbergMarquardt.jl into the SciML interface.
4+
Note that these solvers do not come by default, and thus one needs to install
5+
the package before using these solvers:
6+
7+
```julia
8+
using Pkg
9+
Pkg.add("FastLevenbergMarquardt")
10+
using FastLevenbergMarquardt
11+
```
12+
13+
These methods can be used independently of the rest of NonlinearSolve.jl
14+
15+
## Solver API
16+
17+
```@docs
18+
FastLevenbergMarquardtJL
19+
```

‎docs/src/api/leastsquaresoptim.md‎

Lines changed: 19 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,19 @@
1+
# LeastSquaresOptim.jl
2+
3+
This is a wrapper package for importing solvers from LeastSquaresOptim.jl into the SciML interface.
4+
Note that these solvers do not come by default, and thus one needs to install
5+
the package before using these solvers:
6+
7+
```julia
8+
using Pkg
9+
Pkg.add("LeastSquaresOptim")
10+
using LeastSquaresOptim
11+
```
12+
13+
These methods can be used independently of the rest of NonlinearSolve.jl
14+
15+
## Solver API
16+
17+
```@docs
18+
LeastSquaresOptimJL
19+
```

‎docs/src/api/minpack.md‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1,13 +1,13 @@
11
# MINPACK.jl
22

3-
This is a wrapper package for importing solvers from Sundials into the SciML interface.
3+
This is a wrapper package for importing solvers from MINPACK into the SciML interface.
44
Note that these solvers do not come by default, and thus one needs to install
55
the package before using these solvers:
66

77
```julia
88
using Pkg
9-
Pkg.add("NonlinearSolveMINPACK")
10-
using NonlinearSolveMINPACK
9+
Pkg.add("MINPACK")
10+
using MINPACK
1111
```
1212

1313
These methods can be used independently of the rest of NonlinearSolve.jl

‎docs/src/solvers/NonlinearLeastSquaresSolvers.md‎

Lines changed: 17 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -57,7 +57,23 @@ And the choices for `linsolve` are:
5757

5858
### MINPACK.jl
5959

60+
MINPACK.jl methods are fine for medium-sized nonlinear solves. They are the FORTRAN
61+
standard methods which are used in many places, such as SciPy. However, our benchmarks
62+
demonstrate that these methods are not robust or stable. In addition, they are slower
63+
than the standard methods and do not scale due to lack of sparse Jacobian support.
64+
Thus they are only recommended for benchmarking and testing code conversions.
65+
66+
- `CMINPACK()`: A wrapper for using the classic MINPACK method through [MINPACK.jl](https://github.com/sglyon/MINPACK.jl)
67+
68+
Submethod choices for this algorithm include:
69+
70+
- `:hybr`: Modified version of Powell's algorithm.
71+
- `:lm`: Levenberg-Marquardt.
72+
- `:lmdif`: Advanced Levenberg-Marquardt
73+
- `:hybrd`: Advanced modified version of Powell's algorithm
74+
6075
### Optimization.jl
6176

62-
`NonlinearLeastSquaresProblem`s can be converted into an `OptimizationProblem` and used with any solver of
77+
`NonlinearLeastSquaresProblem`s can be converted into an `OptimizationProblem`
78+
and used with any solver of
6379
[Optimization.jl](https://github.com/SciML/Optimization.jl).

‎ext/NonlinearSolveMINPACKExt.jl‎

Lines changed: 82 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,82 @@
1+
module NonlinearSolveMINPACKExt
2+
3+
using NonlinearSolve, SciMLBase
4+
using MINPACK
5+
6+
function SciMLBase.solve(prob::Union{SciMLBase.NonlinearProblem{uType, isinplace},
7+
SciMLBase.NonlinearLeastSquaresProblem{uType, isinplace}},
8+
alg::CMINPACK,
9+
reltol = 1e-3,
10+
abstol = 1e-6,
11+
maxiters = 100000,
12+
timeseries = [],
13+
ts = [],
14+
ks = [], ;
15+
kwargs...) where {uType, isinplace}
16+
if prob.u0 isa Number
17+
u0 = [prob.u0]
18+
else
19+
u0 = deepcopy(prob.u0)
20+
end
21+
22+
sizeu = size(prob.u0)
23+
p = prob.p
24+
25+
# unwrapping alg params
26+
method = alg.method
27+
show_trace = alg.show_trace
28+
tracing = alg.tracing
29+
io = alg.io
30+
31+
if !isinplace && prob.u0 isa Number
32+
f! = (du, u) -> (du .= prob.f(first(u), p); Cint(0))
33+
elseif !isinplace && prob.u0 isa Vector{Float64}
34+
f! = (du, u) -> (du .= prob.f(u, p); Cint(0))
35+
elseif !isinplace && prob.u0 isa AbstractArray
36+
f! = (du, u) -> (du .= vec(prob.f(reshape(u, sizeu), p)); Cint(0))
37+
elseif prob.u0 isa Vector{Float64}
38+
f! = (du, u) -> prob.f(du, u, p)
39+
else # Then it's an in-place function on an abstract array
40+
f! = (du, u) -> (prob.f(reshape(du, sizeu), reshape(u, sizeu), p);
41+
du = vec(du);
42+
0)
43+
end
44+
45+
u = zero(u0)
46+
resid = similar(u0)
47+
48+
m = prob.f.resid_prototype === nothing ? length(u0) : length(prob.f.resid_prototype)
49+
50+
if SciMLBase.has_jac(prob.f)
51+
if !isinplace && prob.u0 isa Number
52+
g! = (du, u) -> (du .= prob.jac(first(u), p); Cint(0))
53+
elseif !isinplace && prob.u0 isa Vector{Float64}
54+
g! = (du, u) -> (du .= prob.jac(u, p); Cint(0))
55+
elseif !isinplace && prob.u0 isa AbstractArray
56+
g! = (du, u) -> (du .= vec(prob.jac(reshape(u, sizeu), p)); Cint(0))
57+
elseif prob.u0 isa Vector{Float64}
58+
g! = (du, u) -> prob.jac(du, u, p)
59+
else # Then it's an in-place function on an abstract array
60+
g! = (du, u) -> (prob.jac(reshape(du, sizeu), reshape(u, sizeu), p);
61+
du = vec(du);
62+
0)
63+
end
64+
original = MINPACK.fsolve(f!, g!, u0, m;
65+
tol = abstol,
66+
show_trace, tracing, method,
67+
iterations = maxiters, io, kwargs...)
68+
else
69+
original = MINPACK.fsolve(f!, u0, m;
70+
tol = abstol,
71+
show_trace, tracing, method,
72+
iterations = maxiters, io, kwargs...)
73+
end
74+
75+
u = reshape(original.x, size(u))
76+
resid = original.f
77+
retcode = original.converged ? ReturnCode.Success : ReturnCode.Failure
78+
SciMLBase.build_solution(prob, alg, u, resid; retcode = retcode,
79+
original = original)
80+
end
81+
82+
end

‎src/extension_algs.jl‎

Lines changed: 42 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -77,3 +77,45 @@ function FastLevenbergMarquardtJL(linsolve::Symbol = :cholesky; factor = 1e-6,
7777
return FastLevenbergMarquardtJL{linsolve}(factor, factoraccept, factorreject,
7878
factorupdate, minscale, maxscale, minfactor, maxfactor)
7979
end
80+
81+
abstract type MINPACKAlgorithm <: SciMLBase.AbstractNonlinearAlgorithm end
82+
83+
"""
84+
```julia
85+
CMINPACK(;show_trace::Bool=false, tracing::Bool=false, method::Symbol=:hybr,
86+
io::IO=stdout)
87+
```
88+
89+
### Keyword Arguments
90+
91+
- `show_trace`: whether to show the trace.
92+
- `tracing`: who the hell knows what this does. If you find out, please open an issue/PR.
93+
- `method`: the choice of method for the solver.
94+
- `io`: the IO to print any tracing output to.
95+
96+
### Method Choices
97+
98+
The keyword argument `method` can take on different value depending on which method of `fsolve` you are calling. The standard choices of `method` are:
99+
100+
- `:hybr`: Modified version of Powell's algorithm. Uses MINPACK routine [`hybrd1`](https://github.com/devernay/cminpack/blob/d1f5f5a273862ca1bbcf58394e4ac060d9e22c76/hybrd1.c)
101+
- `:lm`: Levenberg-Marquardt. Uses MINPACK routine [`lmdif1`](https://github.com/devernay/cminpack/blob/d1f5f5a273862ca1bbcf58394e4ac060d9e22c76/lmdif1.c)
102+
- `: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
103+
- `: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
104+
105+
If a Jacobian is supplied as part of the [`NonlinearFunction`](@ref nonlinearfunctions),
106+
then the following methods are allowed:
107+
108+
- `: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
109+
- `: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
110+
"""
111+
struct CMINPACK <: MINPACKAlgorithm
112+
show_trace::Bool
113+
tracing::Bool
114+
method::Symbol
115+
io::IO
116+
end
117+
118+
function CMINPACK(; show_trace::Bool = false, tracing::Bool = false, method::Symbol = :hybr,
119+
io::IO = stdout)
120+
CMINPACK(show_trace, tracing, method, io)
121+
end

‎test/minpack.jl‎

Lines changed: 33 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,33 @@
1+
using NonlinearSolve, MINPACK, Test
2+
3+
function f_iip(du, u, p)
4+
du[1] = 2 - 2u[1]
5+
du[2] = u[1] - 4u[2]
6+
end
7+
u0 = zeros(2)
8+
prob_iip = NonlinearProblem{true}(f_iip, u0)
9+
abstol = 1e-8
10+
for alg in [CMINPACK()]
11+
local sol
12+
sol = solve(prob_iip, alg)
13+
@test sol.retcode == ReturnCode.Success
14+
p = nothing
15+
16+
du = zeros(2)
17+
f_iip(du, sol.u, nothing)
18+
@test maximum(du) < 1e-6
19+
end
20+
21+
# OOP Tests
22+
f_oop(u, p) = [2 - 2u[1], u[1] - 4u[2]]
23+
u0 = zeros(2)
24+
prob_oop = NonlinearProblem{false}(f_oop, u0)
25+
for alg in [CMINPACK()]
26+
local sol
27+
sol = solve(prob_oop, alg)
28+
@test sol.retcode == ReturnCode.Success
29+
30+
du = zeros(2)
31+
du = f_oop(sol.u, nothing)
32+
@test maximum(du) < 1e-6
33+
end

0 commit comments

Comments
 (0)