Add constraints support for NLopt

lib/OptimizationNLopt/Project.toml

@@ -6,10 +6,11 @@ version = "0.2.2"
 [deps]
 NLopt = "76087f3c-5699-56af-9a33-bf431cd00edd"
 Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 
 [compat]
-NLopt = "0.6, 1"
+NLopt = "1.0.3"
 Optimization = "3.21"
 Reexport = "1.2"
 julia = "1"

lib/OptimizationNLopt/src/OptimizationNLopt.jl

@@ -3,6 +3,7 @@ module OptimizationNLopt
 using Reexport
 @reexport using NLopt, Optimization
 using Optimization.SciMLBase
+using Optimization: deduce_retcode
 
 (f::NLopt.Algorithm)() = f
 
@@ -63,6 +64,24 @@ function SciMLBase.requiresconsjac(opt::Union{NLopt.Algorithm, NLopt.Opt}) #http
     end
 end
 
+function SciMLBase.allowsconstraints(opt::NLopt.Algorithm)
+    str_opt = string(opt)
+    if occursin("AUGLAG", str_opt) || occursin("CCSA", str_opt) || occursin("MMA", str_opt) || occursin("COBYLA", str_opt) || occursin("ISRES", str_opt) || occursin("AGS", str_opt) || occursin("ORIG_DIRECT", str_opt) || occursin("SLSQP", str_opt) 
+        return true
+    else
+        return false
+    end
+end
+
+function SciMLBase.__init(prob::SciMLBase.OptimizationProblem, opt::NLopt.Algorithm,
+; cons_tol =  1e-6,
+        callback = (args...) -> (false),
+        progress = false, kwargs...)
+    return OptimizationCache(prob, opt; cons_tol, callback, progress,
+        kwargs...)
+end
+
+
 function __map_optimizer_args!(cache::OptimizationCache, opt::NLopt.Opt;
         callback = nothing,
         maxiters::Union{Number, Nothing} = nothing,
@@ -103,7 +122,9 @@ function __map_optimizer_args!(cache::OptimizationCache, opt::NLopt.Opt;
 
     # add optimiser options from kwargs
     for j in kwargs
-        eval(Meta.parse("NLopt." * string(j.first) * "!"))(opt, j.second)
+        if j.first != :cons_tol
+            eval(Meta.parse("NLopt." * string(j.first) * "!"))(opt, j.second)
+        end
     end
 
     if cache.ub !== nothing
@@ -132,31 +153,6 @@ function __map_optimizer_args!(cache::OptimizationCache, opt::NLopt.Opt;
     return nothing
 end
 
-function __nlopt_status_to_ReturnCode(status::Symbol)
-    if status in Symbol.([
-        NLopt.SUCCESS,
-        NLopt.STOPVAL_REACHED,
-        NLopt.FTOL_REACHED,
-        NLopt.XTOL_REACHED,
-        NLopt.ROUNDOFF_LIMITED
-    ])
-        return ReturnCode.Success
-    elseif status == Symbol(NLopt.MAXEVAL_REACHED)
-        return ReturnCode.MaxIters
-    elseif status == Symbol(NLopt.MAXTIME_REACHED)
-        return ReturnCode.MaxTime
-    elseif status in Symbol.([
-        NLopt.OUT_OF_MEMORY,
-        NLopt.INVALID_ARGS,
-        NLopt.FAILURE,
-        NLopt.FORCED_STOP
-    ])
-        return ReturnCode.Failure
-    else
-        return ReturnCode.Default
-    end
-end
-
 function SciMLBase.__solve(cache::OptimizationCache{
         F,
         RC,
@@ -204,6 +200,8 @@ function SciMLBase.__solve(cache::OptimizationCache{
         return _loss(θ)
     end
 
+
+
     opt_setup = if isa(cache.opt, NLopt.Opt)
         if ndims(cache.opt) != length(cache.u0)
             error("Passed NLopt.Opt optimization dimension does not match OptimizationProblem dimension.")
@@ -219,6 +217,37 @@ function SciMLBase.__solve(cache::OptimizationCache{
         NLopt.min_objective!(opt_setup, fg!)
     end
 
+    if cache.f.cons !== nothing
+        eqinds = map((y) -> y[1]==y[2], zip(cache.lcons, cache.ucons))
+        ineqinds = map((y) -> y[1]!=y[2], zip(cache.lcons, cache.ucons))
+        if sum(ineqinds) > 0
+            ineqcons = function (res, θ, J)
+                cons_cache = zeros(eltype(res), sum(eqinds)+sum(ineqinds))
+                cache.f.cons(cons_cache, θ)
+                res .= @view(cons_cache[ineqinds])
+                if length(J) > 0
+                    Jcache = zeros(eltype(J), sum(ineqinds)+sum(eqinds), length(θ))
+                    cache.f.cons_j(Jcache, θ)
+                    J .= @view(Jcache[ineqinds, :])'
+                end
+            end
+            NLopt.inequality_constraint!(opt_setup, ineqcons, [cache.solver_args.cons_tol for i in 1:sum(ineqinds)])
+        end
+        if sum(eqinds) > 0
+            eqcons = function (res, θ, J)
+                cons_cache = zeros(eltype(res), sum(eqinds)+sum(ineqinds))
+                cache.f.cons(cons_cache, θ)
+                res .= @view(cons_cache[eqinds])
+                if length(J) > 0
+                    Jcache = zeros(eltype(res), sum(eqinds)+sum(ineqinds), length(θ))
+                    cache.f.cons_j(Jcache, θ)
+                    J .= @view(Jcache[eqinds, :])'
+                end
+            end
+            NLopt.equality_constraint!(opt_setup, eqcons, [cache.solver_args.cons_tol for i in 1:sum(eqinds)])
+        end
+    end
+
     maxiters = Optimization._check_and_convert_maxiters(cache.solver_args.maxiters)
     maxtime = Optimization._check_and_convert_maxtime(cache.solver_args.maxtime)
 
@@ -229,7 +258,7 @@ function SciMLBase.__solve(cache::OptimizationCache{
     t0 = time()
     (minf, minx, ret) = NLopt.optimize(opt_setup, cache.u0)
     t1 = time()
-    retcode = __nlopt_status_to_ReturnCode(ret)
+    retcode = deduce_retcode(ret)
 
     if retcode == ReturnCode.Failure
         @warn "NLopt failed to converge: $(ret)"

lib/OptimizationNLopt/test/runtests.jl

@@ -1,5 +1,5 @@
 using OptimizationNLopt, Optimization, Zygote
-using Test
+using Test, Random
 
 @testset "OptimizationNLopt.jl" begin
     rosenbrock(x, p) = (p[1] - x[1])^2 + p[2] * (x[2] - x[1]^2)^2
@@ -16,7 +16,7 @@ using Test
     optprob = OptimizationFunction(rosenbrock, Optimization.AutoZygote())
     prob = OptimizationProblem(optprob, x0, _p)
 
-    sol = solve(prob, NLopt.Opt(:LN_BOBYQA, 2))
+    sol = solve(prob, NLopt.Opt(:LD_LBFGS, 2))
     @test sol.retcode == ReturnCode.Success
     @test 10 * sol.objective < l1
 
@@ -26,10 +26,6 @@ using Test
     @test sol.retcode == ReturnCode.Success
     @test 10 * sol.objective < l1
 
-    sol = solve(prob, NLopt.Opt(:LD_LBFGS, 2))
-    @test sol.retcode == ReturnCode.Success
-    @test 10 * sol.objective < l1
-
     sol = solve(prob, NLopt.Opt(:G_MLSL_LDS, 2), local_method = NLopt.Opt(:LD_LBFGS, 2),
         maxiters = 10000)
     @test sol.retcode == ReturnCode.MaxIters
@@ -82,4 +78,46 @@ using Test
         #nlopt gives the last best not the one where callback stops
         @test sol.objective < 0.8
     end
+
+    @testset "constrained" begin
+        cons = (res, x, p) -> res .= [x[1]^2 + x[2]^2 - 1.0]
+        x0 = zeros(2)
+        optprob = OptimizationFunction(rosenbrock, Optimization.AutoZygote();
+            cons = cons)
+        prob = OptimizationProblem(optprob, x0, _p, lcons = [0.0], ucons = [0.0])
+        sol = solve(prob, NLopt.LN_COBYLA())
+        @test sol.retcode == ReturnCode.Success
+        @test 10 * sol.objective < l1
+
+        Random.seed!(1)
+        prob = OptimizationProblem(optprob, rand(2), _p,
+            lcons = [0.0], ucons = [0.0])
+
+        sol = solve(prob, NLopt.LD_SLSQP())
+        @test sol.retcode == ReturnCode.Success
+        @test 10 * sol.objective < l1
+
+        Random.seed!(1)
+        prob = OptimizationProblem(optprob, rand(2), _p,
+            lcons = [0.0], ucons = [0.0])
+        sol = solve(prob, NLopt.AUGLAG(), local_method = NLopt.LD_LBFGS())
+        @test sol.retcode == ReturnCode.Success
+        @test 10 * sol.objective < l1
+
+        function con2_c(res, x, p)
+            res .= [x[1]^2 + x[2]^2 - 1.0, x[2] * sin(x[1]) - x[1] - 2.0]
+        end
+
+        optprob = OptimizationFunction(rosenbrock, Optimization.AutoForwardDiff();cons = con2_c)
+        Random.seed!(1)
+        prob = OptimizationProblem(optprob, rand(2), _p, lcons = [0.0, -Inf], ucons = [0.0, 0.0])
+        sol = solve(prob, NLopt.LD_AUGLAG(), local_method = NLopt.LD_LBFGS())
+        @test sol.retcode == ReturnCode.Success
+        @test 10 * sol.objective < l1
+
+        prob = OptimizationProblem(optprob, rand(2), _p, lcons = [-Inf, -Inf], ucons = [0.0, 0.0], lb = [-1.0, -1.0], ub = [1.0, 1.0])
+        sol = solve(prob, NLopt.GN_ISRES(), maxiters = 1000)
+        @test sol.retcode == ReturnCode.MaxIters
+        @test 10 * sol.objective < l1
+    end
 end

src/utils.jl

@@ -79,7 +79,8 @@ const STOP_REASON_MAP = Dict(
     r"STOP: XTOL.TOO.SMALL" => ReturnCode.ConvergenceFailure,
     r"STOP: TERMINATION" => ReturnCode.Terminated,
     r"Optimization completed" => ReturnCode.Success,
-    r"Convergence achieved" => ReturnCode.Success
+    r"Convergence achieved" => ReturnCode.Success,
+    r"ROUNDOFF_LIMITED" => ReturnCode.Success
 )
 
 # Function to deduce ReturnCode from a stop_reason string using the dictionary