From ec25e8032020e51feff9080fc2798d7b4fb52a56 Mon Sep 17 00:00:00 2001 From: odow Date: Mon, 21 Mar 2022 20:17:33 +1300 Subject: [PATCH 1/4] Update to MathOptInterface v1.0 --- src/solve/moi.jl | 36 ++++++++++++++---------------------- test/Project.toml | 2 +- 2 files changed, 15 insertions(+), 23 deletions(-) diff --git a/src/solve/moi.jl b/src/solve/moi.jl index a32dd9ef4..2ef245094 100644 --- a/src/solve/moi.jl +++ b/src/solve/moi.jl @@ -89,23 +89,13 @@ end function __map_optimizer_args(prob::OptimizationProblem, opt::Union{MOI.AbstractOptimizer, MOI.OptimizerWithAttributes}; maxiters::Union{Number, Nothing}=nothing, maxtime::Union{Number, Nothing}=nothing, - abstol::Union{Number, Nothing}=nothing, + abstol::Union{Number, Nothing}=nothing, reltol::Union{Number, Nothing}=nothing, kwargs...) - - mapped_args = Vector{Pair{String, Any}}[] - mapped_args = [mapped_args..., [Pair(string(j.first),j.second) for j = kwargs]...] - - if isa(opt, MOI.AbstractOptimizer) - if length(mapped_args) > 0 - opt = MOI.OptimizerWithAttributes(typeof(opt), mapped_args...) - else - opt = typeof(opt) - end - end - optimizer = MOI.instantiate(opt) - + for (key, value) in kwargs + MOI.set(optimizer, MOI.RawOptimizerattribute("$(key)"), value) + end if !isnothing(maxtime) MOI.set(optimizer, MOI.TimeLimitSec(), maxtime) end @@ -121,14 +111,14 @@ function __map_optimizer_args(prob::OptimizationProblem, opt::Union{MOI.Abstract if !isnothing(maxiters) @warn "common maxiters argument is currently not used by $(optimizer). Set number of interations via optimizer specific keyword aguments." end - + return optimizer end function __solve(prob::OptimizationProblem, opt::Union{MOI.AbstractOptimizer, MOI.OptimizerWithAttributes}; maxiters::Union{Number, Nothing}=nothing, maxtime::Union{Number, Nothing}=nothing, - abstol::Union{Number, Nothing}=nothing, + abstol::Union{Number, Nothing}=nothing, reltol::Union{Number, Nothing}=nothing, kwargs...) @@ -136,25 +126,27 @@ function __solve(prob::OptimizationProblem, opt::Union{MOI.AbstractOptimizer, MO maxtime = _check_and_convert_maxtime(maxtime) opt_setup = __map_optimizer_args(prob, opt; abstol=abstol, reltol=reltol, maxiters=maxiters, maxtime=maxtime, kwargs...) - + num_variables = length(prob.u0) θ = MOI.add_variables(opt_setup, num_variables) if prob.lb !== nothing @assert eachindex(prob.lb) == Base.OneTo(num_variables) for i in 1:num_variables - MOI.add_constraint(opt_setup, MOI.SingleVariable(θ[i]), MOI.GreaterThan(prob.lb[i])) + MOI.add_constraint(opt_setup, θ[i], MOI.GreaterThan(prob.lb[i])) end end if prob.ub !== nothing @assert eachindex(prob.ub) == Base.OneTo(num_variables) for i in 1:num_variables - MOI.add_constraint(opt_setup, MOI.SingleVariable(θ[i]), MOI.LessThan(prob.ub[i])) + MOI.add_constraint(opt_setup, θ[i], MOI.LessThan(prob.ub[i])) end end @assert eachindex(prob.u0) == Base.OneTo(num_variables) - for i in 1:num_variables - MOI.set(opt_setup, MOI.VariablePrimalStart(), θ[i], prob.u0[i]) - end + if MOI.supports(opt_setup, MOI.VariablePrimalStart(), MOI.VariableIndex) + for i in 1:num_variables + MOI.set(opt_setup, MOI.VariablePrimalStart(), θ[i], prob.u0[i]) + end + end MOI.set(opt_setup, MOI.ObjectiveSense(), prob.sense === MaxSense ? MOI.MAX_SENSE : MOI.MIN_SENSE) if prob.lcons === nothing @assert prob.ucons === nothing diff --git a/test/Project.toml b/test/Project.toml index 4d3b2029e..85e381409 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -48,7 +48,7 @@ ForwardDiff = ">= 0.10.19" GCMAES = ">= 0.1.25" Ipopt = ">= 0.7.0" IterTools = ">= 1.3.0" -MathOptInterface = ">= 0.9.22" +MathOptInterface = ">= 1" Metaheuristics = ">=3.0.2" ModelingToolkit = ">= 6.4.7" MultistartOptimization = ">= 0.1.2" From 7fcb5f497d306c2e8fe8c4f12648f946bce792bd Mon Sep 17 00:00:00 2001 From: Oscar Dowson Date: Mon, 21 Mar 2022 20:26:46 +1300 Subject: [PATCH 2/4] Update src/solve/moi.jl --- src/solve/moi.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/solve/moi.jl b/src/solve/moi.jl index 2ef245094..5726a8a7b 100644 --- a/src/solve/moi.jl +++ b/src/solve/moi.jl @@ -94,7 +94,7 @@ function __map_optimizer_args(prob::OptimizationProblem, opt::Union{MOI.Abstract kwargs...) optimizer = MOI.instantiate(opt) for (key, value) in kwargs - MOI.set(optimizer, MOI.RawOptimizerattribute("$(key)"), value) + MOI.set(optimizer, MOI.RawOptimizerAttribute("$(key)"), value) end if !isnothing(maxtime) MOI.set(optimizer, MOI.TimeLimitSec(), maxtime) From c7d5d98f1387e5c0d33025270168a19552553227 Mon Sep 17 00:00:00 2001 From: odow Date: Tue, 22 Mar 2022 16:34:37 +1300 Subject: [PATCH 3/4] Variety of misc fixes --- src/solve/moi.jl | 204 +++++++++++++++++++++++++++++++---------------- 1 file changed, 134 insertions(+), 70 deletions(-) diff --git a/src/solve/moi.jl b/src/solve/moi.jl index 2ef245094..4f150f931 100644 --- a/src/solve/moi.jl +++ b/src/solve/moi.jl @@ -4,20 +4,96 @@ const MOI = MathOptInterface struct MOIOptimizationProblem{T,F<:OptimizationFunction,uType,P} <: MOI.AbstractNLPEvaluator f::F u0::uType - p::P + p::P J::Matrix{T} H::Matrix{T} cons_H::Vector{Matrix{T}} end -MOI.eval_objective(moiproblem::MOIOptimizationProblem, x) = moiproblem.f(x, moiproblem.p) +function MOIOptimizationProblem(prob::OptimizationProblem) + num_cons = prob.ucons === nothing ? 0 : length(prob.ucons) + f = instantiate_function(prob.f, prob.u0, prob.f.adtype, prob.p, num_cons) + T = eltype(prob.u0) + n = length(prob.u0) + return MOIOptimizationProblem( + f, + prob.u0, + prob.p, + zeros(T, num_cons, n), + zeros(T, n, n), + Matrix{T}[zeros(T, n, n) for i in 1:num_cons], + ) +end + +MOI.features_available(::MOIOptimizationProblem) = [:Grad, :Hess, :Jac] + +function MOI.initialize( + moiproblem::MOIOptimizationProblem, + requested_features::Vector{Symbol}, +) + available_features = MOI.features_available(moiproblem) + for feat in requested_features + if !(feat in available_features) + error("Unsupported feature $feat") + # TODO: implement Jac-vec and Hess-vec products + # for solvers that need them + end + end + return +end + +function MOI.eval_objective(moiproblem::MOIOptimizationProblem, x) + return moiproblem.f(x, moiproblem.p) +end + +function MOI.eval_constraint(moiproblem::MOIOptimizationProblem, g, x) + g .= moiproblem.f.cons(x) + return +end -MOI.eval_constraint(moiproblem::MOIOptimizationProblem, g, x) = g .= moiproblem.f.cons(x) +function MOI.eval_objective_gradient(moiproblem::MOIOptimizationProblem, G, x) + moiproblem.f.grad(G, x) + return +end -MOI.eval_objective_gradient(moiproblem::MOIOptimizationProblem, G, x) = moiproblem.f.grad(G, x) +# This structure assumes the calculation of moiproblem.J is dense. +function MOI.jacobian_structure(moiproblem::MOIOptimizationProblem) + rows, cols = size(moiproblem.J) + return Tuple{Int,Int}[(i, j) for j in 1:cols for i in 1:rows] +end -function MOI.eval_hessian_lagrangian(moiproblem::MOIOptimizationProblem{T}, h, x, σ, μ) where {T} +function MOI.eval_constraint_jacobian(moiproblem::MOIOptimizationProblem, j, x) + if isempty(j) + return + elseif moiproblem.f.cons_j === nothing + error( + "Use OptimizationFunction to pass the derivatives or " * + "automatically generate them with one of the autodiff backends", + ) + end + moiproblem.f.cons_j(moiproblem.J, x) + for i in eachindex(j) + j[i] = moiproblem.J[i] + end + return +end + +# Because the Hessian is symmetrical, we choose to store the upper-triangular +# component. We also assume that it is dense. +function MOI.hessian_lagrangian_structure(moiproblem::MOIOptimizationProblem) + num_vars = length(moiproblem.u0) + return Tuple{Int,Int}[(row, col) for col in 1:num_vars for row in 1:col] +end + +function MOI.eval_hessian_lagrangian( + moiproblem::MOIOptimizationProblem{T}, + h, + x, + σ, + μ, +) where {T} n = length(moiproblem.u0) + # TODO(odow): what is this? Why is it called? a = zeros(n, n) moiproblem.f.hess(a, x) if iszero(σ) @@ -47,107 +123,86 @@ function MOI.eval_hessian_lagrangian(moiproblem::MOIOptimizationProblem{T}, h, x return end -function MOI.eval_constraint_jacobian(moiproblem::MOIOptimizationProblem, j, x) - isempty(j) && return - moiproblem.f.cons_j === nothing && error("Use OptimizationFunction to pass the derivatives or automatically generate them with one of the autodiff backends") - n = length(moiproblem.u0) - moiproblem.f.cons_j(moiproblem.J, x) - for i in eachindex(j) - j[i] = moiproblem.J[i] - end -end - -function MOI.jacobian_structure(moiproblem::MOIOptimizationProblem) - return Tuple{Int,Int}[(con, var) for con in 1:size(moiproblem.J,1) for var in 1:size(moiproblem.J,2)] -end - -function MOI.hessian_lagrangian_structure(moiproblem::MOIOptimizationProblem) - return Tuple{Int,Int}[(row, col) for col in 1:length(moiproblem.u0) for row in 1:col] -end - -function MOI.initialize(moiproblem::MOIOptimizationProblem, requested_features::Vector{Symbol}) - for feat in requested_features - if !(feat in MOI.features_available(moiproblem)) - error("Unsupported feature $feat") - # TODO: implement Jac-vec and Hess-vec products - # for solvers that need them - end - end -end - -MOI.features_available(moiproblem::MOIOptimizationProblem) = [:Grad, :Hess, :Jac] +_create_new_optimizer(opt::MOI.AbstractOptimizer) = opt +_create_new_optimizer(opt::MOI.OptimizerWithAttributes) = MOI.instantiate(opt) -function make_moi_problem(prob::OptimizationProblem) - num_cons = prob.ucons === nothing ? 0 : length(prob.ucons) - f = instantiate_function(prob.f,prob.u0,prob.f.adtype,prob.p,num_cons) - T = eltype(prob.u0) - n = length(prob.u0) - moiproblem = MOIOptimizationProblem(f,prob.u0,prob.p,zeros(T,num_cons,n),zeros(T,n,n),Matrix{T}[zeros(T,n,n) for i in 1:num_cons]) - return moiproblem -end - -function __map_optimizer_args(prob::OptimizationProblem, opt::Union{MOI.AbstractOptimizer, MOI.OptimizerWithAttributes}; +function __map_optimizer_args( + prob::OptimizationProblem, + opt::Union{MOI.AbstractOptimizer, MOI.OptimizerWithAttributes}; maxiters::Union{Number, Nothing}=nothing, maxtime::Union{Number, Nothing}=nothing, abstol::Union{Number, Nothing}=nothing, reltol::Union{Number, Nothing}=nothing, - kwargs...) - optimizer = MOI.instantiate(opt) + kwargs..., +) + optimizer = _create_new_optimizer(opt) for (key, value) in kwargs MOI.set(optimizer, MOI.RawOptimizerattribute("$(key)"), value) end if !isnothing(maxtime) MOI.set(optimizer, MOI.TimeLimitSec(), maxtime) end - if !isnothing(reltol) @warn "common reltol argument is currently not used by $(optimizer). Set tolerances via optimizer specific keyword aguments." end - if !isnothing(abstol) @warn "common abstol argument is currently not used by $(optimizer). Set tolerances via optimizer specific keyword aguments." end - if !isnothing(maxiters) @warn "common maxiters argument is currently not used by $(optimizer). Set number of interations via optimizer specific keyword aguments." end - return optimizer end -function __solve(prob::OptimizationProblem, opt::Union{MOI.AbstractOptimizer, MOI.OptimizerWithAttributes}; +function __solve( + prob::OptimizationProblem, + opt::Union{MOI.AbstractOptimizer, MOI.OptimizerWithAttributes}; maxiters::Union{Number, Nothing}=nothing, maxtime::Union{Number, Nothing}=nothing, abstol::Union{Number, Nothing}=nothing, reltol::Union{Number, Nothing}=nothing, - kwargs...) - + kwargs..., +) maxiters = _check_and_convert_maxiters(maxiters) maxtime = _check_and_convert_maxtime(maxtime) - - opt_setup = __map_optimizer_args(prob, opt; abstol=abstol, reltol=reltol, maxiters=maxiters, maxtime=maxtime, kwargs...) - + opt_setup = __map_optimizer_args( + prob, + opt; + abstol=abstol, + reltol=reltol, + maxiters=maxiters, + maxtime=maxtime, + kwargs..., + ) num_variables = length(prob.u0) - θ = MOI.add_variables(opt_setup, num_variables) - if prob.lb !== nothing + θ = MOI.add_variables(opt_setup, num_variables) + if prob.lb !== nothing @assert eachindex(prob.lb) == Base.OneTo(num_variables) - for i in 1:num_variables - MOI.add_constraint(opt_setup, θ[i], MOI.GreaterThan(prob.lb[i])) + for i in 1:num_variables + if prob.lb[i] > -Inf + MOI.add_constraint(opt_setup, θ[i], MOI.GreaterThan(prob.lb[i])) + end end end - if prob.ub !== nothing + if prob.ub !== nothing @assert eachindex(prob.ub) == Base.OneTo(num_variables) - for i in 1:num_variables - MOI.add_constraint(opt_setup, θ[i], MOI.LessThan(prob.ub[i])) + for i in 1:num_variables + if prob.ub[i] < Inf + MOI.add_constraint(opt_setup, θ[i], MOI.LessThan(prob.ub[i])) + end end end - @assert eachindex(prob.u0) == Base.OneTo(num_variables) if MOI.supports(opt_setup, MOI.VariablePrimalStart(), MOI.VariableIndex) + @assert eachindex(prob.u0) == Base.OneTo(num_variables) for i in 1:num_variables MOI.set(opt_setup, MOI.VariablePrimalStart(), θ[i], prob.u0[i]) end end - MOI.set(opt_setup, MOI.ObjectiveSense(), prob.sense === MaxSense ? MOI.MAX_SENSE : MOI.MIN_SENSE) + MOI.set( + opt_setup, + MOI.ObjectiveSense(), + prob.sense === MaxSense ? MOI.MAX_SENSE : MOI.MIN_SENSE, + ) if prob.lcons === nothing @assert prob.ucons === nothing con_bounds = MOI.NLPBoundsPair[] @@ -155,10 +210,12 @@ function __solve(prob::OptimizationProblem, opt::Union{MOI.AbstractOptimizer, MO @assert prob.ucons !== nothing con_bounds = MOI.NLPBoundsPair.(prob.lcons, prob.ucons) end - MOI.set(opt_setup, MOI.NLPBlock(), MOI.NLPBlockData(con_bounds, make_moi_problem(prob), true)) - - MOI.optimize!(opt_setup) - + MOI.set( + opt_setup, + MOI.NLPBlock(), + MOI.NLPBlockData(con_bounds, MOIOptimizationProblem(prob), true), + ) + MOI.optimize!(opt_setup) if MOI.get(opt_setup, MOI.ResultCount()) >= 1 minimizer = MOI.get(opt_setup, MOI.VariablePrimal(), θ) minimum = MOI.get(opt_setup, MOI.ObjectiveValue()) @@ -168,5 +225,12 @@ function __solve(prob::OptimizationProblem, opt::Union{MOI.AbstractOptimizer, MO minimum = NaN opt_ret= :Default end - SciMLBase.build_solution(prob, opt, minimizer, minimum; original=opt_setup, retcode=opt_ret) + return SciMLBase.build_solution( + prob, + opt, + minimizer, + minimum; + original=opt_setup, + retcode=opt_ret, + ) end From 72155b62ce4120058c76cc6a86c21c347faaf8f5 Mon Sep 17 00:00:00 2001 From: odow Date: Tue, 22 Mar 2022 18:45:21 +1300 Subject: [PATCH 4/4] Remove unneeded hessian call --- src/solve/moi.jl | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/solve/moi.jl b/src/solve/moi.jl index 79017e3fc..3badbc4a9 100644 --- a/src/solve/moi.jl +++ b/src/solve/moi.jl @@ -93,9 +93,6 @@ function MOI.eval_hessian_lagrangian( μ, ) where {T} n = length(moiproblem.u0) - # TODO(odow): what is this? Why is it called? - a = zeros(n, n) - moiproblem.f.hess(a, x) if iszero(σ) fill!(h, zero(T)) else