diff --git a/Project.toml b/Project.toml index c8defac2..f8a77881 100644 --- a/Project.toml +++ b/Project.toml @@ -8,6 +8,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" NLPModels = "a4795742-8479-5a88-8948-cc11e1c8c1a6" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" +QuadraticModels = "f468eda6-eac5-11e8-05a5-ff9e497bcd19" SolverCore = "ff4d7338-4cf1-434d-91df-b86cb86fb843" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" @@ -19,6 +20,7 @@ NLPModels = "0.21.6" NLPModelsTest = "0.10.6" Percival = "0.7.3" Printf = "1.10" +QuadraticModels = "0.9.16" SolverCore = "0.3.9" SparseArrays = "1.10" Test = "1.10" diff --git a/src/MOI_wrapper.jl b/src/MOI_wrapper.jl index 9ec6d2ce..8818e0d7 100644 --- a/src/MOI_wrapper.jl +++ b/src/MOI_wrapper.jl @@ -4,7 +4,7 @@ mutable struct Optimizer <: MOI.AbstractOptimizer options::Dict{String, Any} silent::Bool solver - nlp::Union{Nothing, MathOptNLPModel} + nlp::Union{Nothing, AbstractNLPModel{Float64, Vector{Float64}}} stats::Union{ Nothing, SolverCore.GenericExecutionStats{Float64, Vector{Float64}, Vector{Float64}, Any}, @@ -92,7 +92,11 @@ function MOI.copy_to(dest::Optimizer, src::MOI.ModelLike) "No solver specified, use for instance `using Percival; JuMP.set_attribute(model, \"solver\", PercivalSolver)`", ) end - dest.nlp, index_map = nlp_model(src) + if is_qp_model(src) + dest.nlp, index_map = qp_model(src) + else + dest.nlp, index_map = nlp_model(src) + end dest.solver = dest.options["solver"](dest.nlp) return index_map end diff --git a/src/NLPModelsJuMP.jl b/src/NLPModelsJuMP.jl index 5e726539..5955397c 100644 --- a/src/NLPModelsJuMP.jl +++ b/src/NLPModelsJuMP.jl @@ -3,6 +3,7 @@ module NLPModelsJuMP include("utils.jl") include("moi_nlp_model.jl") include("moi_nls_model.jl") +include("moi_qp_model.jl") include("MOI_wrapper.jl") end diff --git a/src/moi_qp_model.jl b/src/moi_qp_model.jl new file mode 100644 index 00000000..5299d747 --- /dev/null +++ b/src/moi_qp_model.jl @@ -0,0 +1,195 @@ +export MathOptQPModel, qp_model, is_qp_model + +import QuadraticModels + +MOI.Utilities.@product_of_sets( + _QPProductOfSets, + MOI.EqualTo{T}, + MOI.GreaterThan{T}, + MOI.LessThan{T}, + MOI.Interval{T}, +) + +const QPOptimizerCache = MOI.Utilities.GenericModel{ + Float64, + MOI.Utilities.ObjectiveContainer{Float64}, + MOI.Utilities.VariablesContainer{Float64}, + MOI.Utilities.MatrixOfConstraints{ + Float64, + MOI.Utilities.MutableSparseMatrixCSC{Float64, Int, MOI.Utilities.OneBasedIndexing}, + MOI.Utilities.Hyperrectangle{Float64}, + _QPProductOfSets{Float64}, + }, +} + +""" + is_qp_model(moimodel::MOI.ModelLike) + +Return `true` if `moimodel` can be represented as a +`QuadraticModels.QuadraticModel`: no `NLPBlock`, no user-defined nonlinear +functions, a linear or quadratic objective, and constraints restricted to +`VariableIndex` and `ScalarAffineFunction` in scalar linear sets (`EqualTo`, +`GreaterThan`, `LessThan`, `Interval`). +""" +function is_qp_model(model::MOI.ModelLike) + nlp_block = MOI.get(model, MOI.NLPBlock()) + if nlp_block !== nothing && + (length(nlp_block.constraint_bounds) > 0 || nlp_block.has_objective) + return false + end + for attr in MOI.get(model, MOI.ListOfModelAttributesSet()) + if attr isa MOI.UserDefinedFunction + return false + end + end + F_obj = MOI.get(model, MOI.ObjectiveFunctionType()) + if !(F_obj <: Union{MOI.VariableIndex, SAF, SQF}) + return false + end + for (F, S) in MOI.get(model, MOI.ListOfConstraintTypesPresent()) + if F == MOI.VariableIndex + S <: ALS || return false + elseif F == SAF + S <: ALS || return false + else + return false + end + end + return true +end + +""" + qp_model(moimodel::MOI.ModelLike; name::String = "Generic") + +Build a `QuadraticModels.QuadraticModel` from `moimodel`. The model must +contain only linear constraints and at most a quadratic objective (see +[`is_qp_model`](@ref)). The matrices are extracted by copying `moimodel` into a +[`MOI.Utilities.MatrixOfConstraints`](@ref) cache. + +Return `(qp, index_map)`. +""" +function qp_model(moimodel::MOI.ModelLike; name::String = "Generic") + cache = MOI.Utilities.UniversalFallback(QPOptimizerCache()) + index_map = MOI.copy_to(cache, moimodel) + return _qp_from_cache(cache, name), index_map +end + +function _qp_from_cache( + cache::MOI.Utilities.UniversalFallback{QPOptimizerCache}, + name::String, +) + src = cache.model + + for (F, S) in MOI.get(cache, MOI.ListOfConstraintTypesPresent()) + if !MOI.supports_constraint(src, F, S) + throw(MOI.UnsupportedConstraint{F, S}()) + end + end + + Ab = src.constraints + A_csc = convert(SparseArrays.SparseMatrixCSC{Float64, Int}, Ab.coefficients) + m, nvar = size(A_csc) + + Arows = Int[] + Acols = Int[] + Avals = Float64[] + rowvals = SparseArrays.rowvals(A_csc) + nzvals = SparseArrays.nonzeros(A_csc) + for j = 1:nvar + for k in SparseArrays.nzrange(A_csc, j) + push!(Arows, rowvals[k]) + push!(Acols, j) + push!(Avals, nzvals[k]) + end + end + + lcon = copy(Ab.constants.lower) + ucon = copy(Ab.constants.upper) + + vc = src.variables + lvar = copy(vc.lower) + uvar = copy(vc.upper) + + x0 = zeros(Float64, nvar) + if MOI.VariablePrimalStart() in MOI.get(src, MOI.ListOfVariableAttributesSet()) + for vi in MOI.get(src, MOI.ListOfVariableIndices()) + val = MOI.get(src, MOI.VariablePrimalStart(), vi) + if val !== nothing + x0[vi.value] = val + end + end + end + + c = zeros(Float64, nvar) + c0 = 0.0 + Hrows = Int[] + Hcols = Int[] + Hvals = Float64[] + sense = MOI.get(src, MOI.ObjectiveSense()) + if sense != MOI.FEASIBILITY_SENSE + F = MOI.get(src, MOI.ObjectiveFunctionType()) + obj = MOI.get(src, MOI.ObjectiveFunction{F}()) + if F == MOI.VariableIndex + c[obj.value] = 1.0 + elseif F == SAF + c0 = obj.constant + for term in obj.terms + c[term.variable.value] += term.coefficient + end + elseif F == SQF + c0 = obj.constant + for term in obj.affine_terms + c[term.variable.value] += term.coefficient + end + for term in obj.quadratic_terms + i, j = term.variable_1.value, term.variable_2.value + if i ≥ j + push!(Hrows, i) + push!(Hcols, j) + else + push!(Hrows, j) + push!(Hcols, i) + end + push!(Hvals, term.coefficient) + end + else + error("Objective function type $F is not supported by qp_model.") + end + end + + minimize = sense != MOI.MAX_SENSE + + return QuadraticModels.QuadraticModel( + c, + Hrows, + Hcols, + Hvals; + Arows = Arows, + Acols = Acols, + Avals = Avals, + lcon = lcon, + ucon = ucon, + lvar = lvar, + uvar = uvar, + c0 = c0, + x0 = x0, + minimize = minimize, + name = name, + ) +end + +""" + MathOptQPModel(jmodel::JuMP.Model; name::String = "Generic") + MathOptQPModel(moimodel::MOI.ModelLike; name::String = "Generic") + +Construct a [`QuadraticModels.QuadraticModel`](@extref) from a JuMP or MOI +model containing only linear constraints and at most a quadratic objective. +""" +function MathOptQPModel(jmodel::JuMP.Model; kws...) + _nlp_sync!(jmodel) + return MathOptQPModel(backend(jmodel); kws...) +end + +function MathOptQPModel(moimodel::MOI.ModelLike; kws...) + return qp_model(moimodel; kws...)[1] +end diff --git a/test/Project.toml b/test/Project.toml new file mode 100644 index 00000000..e7492b91 --- /dev/null +++ b/test/Project.toml @@ -0,0 +1,11 @@ +[deps] +JuMP = "4076af6c-e467-56ae-b986-b466b2749572" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +NLPModels = "a4795742-8479-5a88-8948-cc11e1c8c1a6" +NLPModelsJuMP = "792afdf1-32c1-5681-94e0-d7bf7a5df49e" +NLPModelsTest = "7998695d-6960-4d3a-85c4-e1bceb8cd856" +Percival = "01435c0c-c90d-11e9-3788-63660f8fbccc" +Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" +QuadraticModels = "f468eda6-eac5-11e8-05a5-ff9e497bcd19" +SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"