diff --git a/Project.toml b/Project.toml index df4cc71b..a035454e 100644 --- a/Project.toml +++ b/Project.toml @@ -15,7 +15,7 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" [compat] BinaryProvider = "0.3, 0.4, 0.5" -MathOptInterface = "0.9.7" +MathOptInterface = "0.10.2" Requires = "1" SCS_GPU_jll = "2.1" SCS_jll = "2.1" diff --git a/README.md b/README.md index 4babe099..b9ea75c9 100644 --- a/README.md +++ b/README.md @@ -185,7 +185,7 @@ mutable struct Solution{T} x::Vector{Float64} y::Vector{Float64} s::Vector{Float64} - info::SCSInfo{T} + info::ScsInfo{T} ret_val::T end ``` @@ -193,7 +193,7 @@ where `x` stores the optimal value of the primal variable, `y` stores the optimal value of the dual variable, `s` is the slack variable, and `info` contains various information about the solve step. -`SCS.raw_status(::SCSInfo)::String` describes the status, e.g. 'Solved', +`SCS.raw_status(::ScsInfo)::String` describes the status, e.g. 'Solved', 'Indeterminate', 'Infeasible/Inaccurate', etc. ## Custom Installation diff --git a/src/MOI_wrapper/MOI_wrapper.jl b/src/MOI_wrapper/MOI_wrapper.jl index 92c5b4ca..65aeec50 100644 --- a/src/MOI_wrapper/MOI_wrapper.jl +++ b/src/MOI_wrapper/MOI_wrapper.jl @@ -1,19 +1,51 @@ -import MathOptInterface +using MathOptInterface const MOI = MathOptInterface -include("matrix.jl") -include("geometric.jl") +include("scaled_psd_cone_bridge.jl") -""" - ADMMIterations() - -The number of ADMM iterations completed during the solve. -""" -struct ADMMIterations <: MOI.AbstractModelAttribute end +MOI.Utilities.@product_of_sets( + Cones, + MOI.Zeros, + MOI.Nonnegatives, + MOI.SecondOrderCone, + ScaledPSDCone, + MOI.ExponentialCone, + MOI.DualExponentialCone, + MOI.PowerCone{T}, + MOI.DualPowerCone{T}, +) -MOI.is_set_by_optimize(::ADMMIterations) = true +const OptimizerCache{T} = MOI.Utilities.GenericModel{ + Cdouble, + MOI.Utilities.ObjectiveContainer{Cdouble}, + MOI.Utilities.VariablesContainer{Cdouble}, + MOI.Utilities.MatrixOfConstraints{ + Cdouble, + MOI.Utilities.MutableSparseMatrixCSC{ + Cdouble, + T, + MOI.Utilities.ZeroBasedIndexing, + }, + Vector{Cdouble}, + Cones{Cdouble}, + }, +} + +function _to_sparse( + ::Type{T}, + A::MOI.Utilities.MutableSparseMatrixCSC{ + Cdouble, + T, + MOI.Utilities.ZeroBasedIndexing, + }, +) where {T} + return -A.nzval, A.rowval, A.colptr +end mutable struct MOISolution + primal::Vector{Float64} + dual::Vector{Float64} + slack::Vector{Float64} ret_val::Int raw_status::String objective_value::Float64 @@ -25,6 +57,9 @@ end function MOISolution() return MOISolution( + Float64[], + Float64[], + Float64[], 0, # SCS_UNFINISHED "", NaN, @@ -35,92 +70,62 @@ function MOISolution() ) end -# This is tied to SCS's internal representation -struct ConeData - qa::Vector{Int} # array of second-order cone constraints - sa::Vector{Int} # array of semi-definite constraints - p::Vector{Float64} # array of power cone params - ConeData() = new(Int[], Int[], Float64[]) -end - -const CONE_TYPES = ( - MOI.Zeros, - MOI.Nonnegatives, - MOI.SecondOrderCone, - MOI.PositiveSemidefiniteConeTriangle, - MOI.ExponentialCone, - MOI.DualExponentialCone, - MOI.PowerCone{Float64}, - MOI.DualPowerCone{Float64}, -) - -const Form{T} = - GeometricConicForm{Float64,SparseMatrixCSRtoCSC{T},typeof(CONE_TYPES)} - mutable struct Optimizer <: MOI.AbstractOptimizer - cone::ConeData - data::Union{Form{Int},Form{Int32}} + cones::Union{Nothing,Cones{Cdouble}} sol::MOISolution silent::Bool options::Dict{Symbol,Any} function Optimizer(; kwargs...) - optimizer = new( - ConeData(), - GeometricConicForm{Float64,SparseMatrixCSRtoCSC{Int}}(CONE_TYPES), - MOISolution(), - false, - Dict{Symbol,Any}(), - ) + optimizer = new(nothing, MOISolution(), false, Dict{Symbol,Any}()) for (key, value) in kwargs - MOI.set(optimizer, MOI.RawParameter(String(key)), value) + MOI.set(optimizer, MOI.RawOptimizerAttribute(String(key)), value) end return optimizer end end +function MOI.get(::Optimizer, ::MOI.Bridges.ListOfNonstandardBridges) + return [ScaledPSDConeBridge{Cdouble}] +end + MOI.get(::Optimizer, ::MOI.SolverName) = "SCS" -function MOI.set(optimizer::Optimizer, param::MOI.RawParameter, value) - # TODO(odow): remove warning in future version. - if !(param.name isa String) - Base.depwarn( - "passing `$(param.name)` to `MOI.RawParameter` as type " * - "`$(typeof(param.name))` is deprecated. Use a string instead.", - Symbol("MOI.set"), - ) - end +MOI.is_empty(optimizer::Optimizer) = optimizer.cones === nothing + +function MOI.empty!(optimizer::Optimizer) + optimizer.cones = nothing + optimizer.sol = MOISolution() + return +end + +### +### MOI.RawOptimizerAttribute +### + +function MOI.set(optimizer::Optimizer, param::MOI.RawOptimizerAttribute, value) return optimizer.options[Symbol(param.name)] = value end -function MOI.get(optimizer::Optimizer, param::MOI.RawParameter) - # TODO(odow): remove warning in future version. - if !(param.name isa String) - Base.depwarn( - "passing $(param.name) to `MOI.RawParameter` as type " * - "$(typeof(param.name)) is deprecated. Use a string instead.", - Symbol("MOI.get"), - ) - end +function MOI.get(optimizer::Optimizer, param::MOI.RawOptimizerAttribute) return optimizer.options[Symbol(param.name)] end +### +### MOI.Silent +### + MOI.supports(::Optimizer, ::MOI.Silent) = true function MOI.set(optimizer::Optimizer, ::MOI.Silent, value::Bool) - return optimizer.silent = value + optimizer.silent = value + return end MOI.get(optimizer::Optimizer, ::MOI.Silent) = optimizer.silent -MOI.is_empty(optimizer::Optimizer) = MOI.is_empty(optimizer.data) - -function MOI.empty!(optimizer::Optimizer) - empty!(optimizer.cone.qa) - empty!(optimizer.cone.sa) - empty!(optimizer.cone.p) - MOI.empty!(optimizer.data) - return optimizer.sol.ret_val = 0 -end +### +### MOI.AbstractModelAttribute +### function MOI.supports( ::Optimizer, @@ -132,231 +137,176 @@ function MOI.supports( return true end -function MOI.supports( - ::Optimizer, - ::MOI.VariablePrimalStart, - ::Type{MOI.VariableIndex}, -) - return true -end +### +### MOI.AbstractVariableAttribute +### + +# TODO(odow): support variable primal starts +# function MOI.supports( +# ::Optimizer, +# ::MOI.VariablePrimalStart, +# ::Type{MOI.VariableIndex}, +# ) +# return true +# end + +### +### MOI.AbstractConstraintAttribute +### + +# TODO(odow): support constraint primal and dual starts +# function MOI.supports( +# ::Optimizer, +# ::Union{MOI.ConstraintPrimalStart,MOI.ConstraintDualStart}, +# ::Type{<:MOI.ConstraintIndex}, +# ) +# return true +# end -function MOI.supports( +function MOI.supports_constraint( ::Optimizer, - ::Union{MOI.ConstraintPrimalStart,MOI.ConstraintDualStart}, - ::Type{<:MOI.ConstraintIndex}, + ::Type{MOI.VectorAffineFunction{Cdouble}}, + ::Type{ + <:Union{ + MOI.Zeros, + MOI.Nonnegatives, + MOI.SecondOrderCone, + ScaledPSDCone, + MOI.ExponentialCone, + MOI.DualExponentialCone, + MOI.PowerCone{Cdouble}, + MOI.DualPowerCone{Cdouble}, + }, + }, ) return true end -function MOI.supports_constraint( - optimizer::Optimizer, - F::Type{MOI.VectorAffineFunction{Float64}}, - S::Type{<:MOI.AbstractVectorSet}, -) - return MOI.supports_constraint(optimizer.data, F, S) +function _map_sets(f, ::Type{T}, sets, ::Type{S}) where {T,S} + F = MOI.VectorAffineFunction{Cdouble} + cis = MOI.get(sets, MOI.ListOfConstraintIndices{F,S}()) + return T[f(MOI.get(sets, MOI.ConstraintSet(), ci)) for ci in cis] end -function MOI.copy_to(dest::Optimizer, src::MOI.ModelLike; kws...) +function MOI.optimize!( + dest::Optimizer, + src::MOI.Utilities.UniversalFallback{OptimizerCache{T}}, +) where {T} linear_solver = get(dest.options, :linear_solver, DirectSolver) - T = scsint_t(linear_solver) - if !(dest.data isa Form{T}) - dest.data = - GeometricConicForm{Float64,SparseMatrixCSRtoCSC{T}}(CONE_TYPES) + if T != scsint_t(linear_solver) + cache = MOI.Utilities.UniversalFallback( + OptimizerCache{scsint_t(linear_solver)}(), + ) + MOI.copy_to(cache, src) + return MOI.optimize!(dest, cache) end - function preprocess_constraint(func, set) - _store_cone_data(dest.cone, set) - return _preprocess_function(func, set) + # The real stuff starts here. + MOI.empty!(dest) + index_map = MOI.Utilities.identity_index_map(src) + Ab = src.model.constraints + A = Ab.coefficients + + model_attributes = MOI.get(src, MOI.ListOfModelAttributesSet()) + objective_function_attr = + MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Cdouble}}() + for attr in model_attributes + if attr != MOI.ObjectiveSense() && attr != objective_function_attr + throw(MOI.UnsupportedAttribute(attr)) + end end - return MOI.copy_to( - dest.data, - src; - preprocess = preprocess_constraint, - kws..., - ) -end - -_store_cone_data(::ConeData, ::MOI.Zeros) = nothing - -_store_cone_data(::ConeData, ::MOI.Nonnegatives) = nothing - -function _store_cone_data(cone::ConeData, set::MOI.SecondOrderCone) - return push!(cone.qa, set.dimension) -end - -function _store_cone_data( - cone::ConeData, - set::MOI.PositiveSemidefiniteConeTriangle, -) - return push!(cone.sa, set.side_dimension) -end - -_store_cone_data(::ConeData, ::MOI.ExponentialCone) = nothing - -_store_cone_data(::ConeData, ::MOI.DualExponentialCone) = nothing - -function _store_cone_data(cone::ConeData, set::MOI.PowerCone{Float64}) - return push!(cone.p, set.exponent) -end - -function _store_cone_data(cone::ConeData, set::MOI.DualPowerCone{Float64}) - # SCS' convention: dual cones have a negative exponent. - return push!(cone.p, -set.exponent) -end - -# Vectorized length for matrix dimension n -sympackedlen(n) = div(n * (n + 1), 2) - -# Matrix dimension for vectorized length n -sympackeddim(n) = div(isqrt(1 + 8n) - 1, 2) - -trimap(i::Integer, j::Integer) = i < j ? trimap(j, i) : div((i - 1) * i, 2) + j - -function trimapL(i::Integer, j::Integer, n::Integer) - return i < j ? trimapL(j, i, n) : i + div((2n - j) * (j - 1), 2) -end - -function _sympackedto(x, n, mapfrom, mapto) - @assert length(x) == sympackedlen(n) - y = similar(x) - for i in 1:n, j in 1:i - y[mapto(i, j)] = x[mapfrom(i, j)] + max_sense = false + if MOI.ObjectiveSense() in model_attributes + max_sense = MOI.get(src, MOI.ObjectiveSense()) == MOI.MAX_SENSE end - return y -end - -function sympackedLtoU(x, n = sympackeddim(length(x))) - return _sympackedto(x, n, (i, j) -> trimapL(i, j, n), trimap) -end - -sympackedUtoL(x, n) = _sympackedto(x, n, trimap, (i, j) -> trimapL(i, j, n)) - -# Scale coefficients depending on rows index -# rows: List of row indices -# coef: List of corresponding coefficients -# d: dimension of set -# rev: if true, we unscale instead (e.g. divide by √2 instead of multiply for -# PSD cone) -function _scalecoef( - rows::AbstractVector{<:Integer}, - coef::Vector{Float64}, - d::Integer, - rev::Bool, -) - scaling = rev ? 1 / √2 : 1 * √2 - output = copy(coef) - for i in 1:length(output) - if !MOI.Utilities.is_diagonal_vectorized_index(rows[i]) - output[i] *= scaling + objective_constant = 0.0 + c = zeros(A.n) + if objective_function_attr in model_attributes + obj = MOI.get(src, objective_function_attr) + objective_constant = MOI.constant(obj) + for term in obj.terms + c[term.variable.value] += (max_sense ? -1 : 1) * term.coefficient end end - return output -end - -# Unscale the coefficients in `coef` with respective rows in `rows` for a set `s` -scalecoef(rows, coef, s) = _scalecoef(rows, coef, MOI.dimension(s), false) - -# Unscale the coefficients in `coef` with respective rows in `rows` for a set of -# type `S` with dimension `d` -function unscalecoef(coef) - return _scalecoef(eachindex(coef), coef, sympackeddim(length(coef)), true) -end - -function _scale(i, coef) - if MOI.Utilities.is_diagonal_vectorized_index(i) - return coef - else - return coef * √2 + # `model.primal` contains the result of the previous optimization. + # It is used as a warm start if its length is the same, e.g. + # probably because no variable and/or constraint has been added. + if A.n != length(dest.sol.primal) + resize!(dest.sol.primal, A.n) + fill!(dest.sol.primal, 0.0) end -end - -function _preprocess_function(func, set::MOI.PositiveSemidefiniteConeTriangle) - n = set.side_dimension - LtoU_map = sympackedLtoU(1:sympackedlen(n), n) - function map_term(t::MOI.VectorAffineTerm) - return MOI.VectorAffineTerm( - LtoU_map[t.output_index], - MOI.ScalarAffineTerm( - _scale(t.output_index, t.scalar_term.coefficient), - t.scalar_term.variable_index, - ), - ) + if A.m != length(dest.sol.dual) + resize!(dest.sol.dual, A.m) + fill!(dest.sol.dual, 0.0) end - UtoL_map = sympackedUtoL(1:sympackedlen(n), n) - function constant(row) - i = UtoL_map[row] - return _scale(i, func.constants[i]) + if A.m != length(dest.sol.slack) + resize!(dest.sol.slack, A.m) + fill!(dest.sol.slack, 0.0) end - new_func = MOI.VectorAffineFunction{Float64}( - MOI.VectorAffineTerm{Float64}[map_term(t) for t in func.terms], - constant.(eachindex(func.constants)), - ) - # The rows have been reordered in `map_term` so we need to re-canonicalize - # to reorder the rows. - MOI.Utilities.canonicalize!(new_func) - return new_func -end - -_preprocess_function(func, set) = func - -function MOI.optimize!(optimizer::Optimizer) - data = optimizer.data - m = data.A.m - n = data.A.n - b = data.b - objective_constant = data.objective_constant - c = data.c - if data.sense == MOI.MAX_SENSE - c = -c + # Set starting values and throw error for other variable attributes + vis_src = MOI.get(src, MOI.ListOfVariableIndices()) + MOI.Utilities.pass_attributes(dest, src, index_map, vis_src) + # Set starting values and throw error for other constraint attributes + for (F, S) in MOI.get(src, MOI.ListOfConstraintTypesPresent()) + cis_src = MOI.get(src, MOI.ListOfConstraintIndices{F,S}()) + MOI.Utilities.pass_attributes(dest, src, index_map, cis_src) end - - options = optimizer.options - if optimizer.silent - options = copy(options) + options = copy(dest.options) + if dest.silent options[:verbose] = 0 end - linear_solver = get(optimizer.options, :linear_solver, DirectSolver) - cone = optimizer.cone - sol = SCS_solve( + dest.cones = deepcopy(Ab.sets) + sol = scs_solve( linear_solver, - m, - n, - data.A, - b, + A.m, + A.n, + A, + Ab.constants, c, - data.num_rows[1], - data.num_rows[2], - cone.qa, - cone.sa, - div(data.num_rows[5], 3), - div(data.num_rows[6], 3), - cone.p, - data.primal, - data.dual, - data.slack; + Ab.sets.num_rows[1], + Ab.sets.num_rows[2] - Ab.sets.num_rows[1], + _map_sets(MOI.dimension, T, Ab, MOI.SecondOrderCone), + _map_sets(MOI.side_dimension, T, Ab, ScaledPSDCone), + div(Ab.sets.num_rows[5] - Ab.sets.num_rows[4], 3), + div(Ab.sets.num_rows[6] - Ab.sets.num_rows[5], 3), + vcat( + _map_sets(set -> set.exponent, Float64, Ab, MOI.PowerCone{Cdouble}), + _map_sets( + set -> -set.exponent, + Float64, + Ab, + MOI.DualPowerCone{Cdouble}, + ), + ), + dest.sol.primal, + dest.sol.dual, + dest.sol.slack; options..., ) - - data.primal = sol.x - data.dual = sol.y - data.slack = sol.s - ret_val = sol.ret_val - objective_value = (data.sense == MOI.MAX_SENSE ? -1 : 1) * sol.info.pobj - dual_objective_value = - (data.sense == MOI.MAX_SENSE ? -1 : 1) * sol.info.dobj - solve_time = (sol.info.setupTime + sol.info.solveTime) / 1000 - optimizer.sol = MOISolution( - ret_val, + dest.sol = MOISolution( + sol.x, + sol.y, + sol.s, + sol.ret_val, raw_status(sol.info), - objective_value, - dual_objective_value, + (max_sense ? -1 : 1) * sol.info.pobj, + (max_sense ? -1 : 1) * sol.info.dobj, objective_constant, - solve_time, + (sol.info.setupTime + sol.info.solveTime) / 1000, sol.info.iter, ) - return + return index_map, false +end + +function MOI.optimize!(dest::Optimizer, src::MOI.ModelLike) + linear_solver = get(dest.options, :linear_solver, DirectSolver) + T = scsint_t(linear_solver) + cache = MOI.Utilities.UniversalFallback(OptimizerCache{T}()) + index_map = MOI.copy_to(cache, src) + MOI.optimize!(dest, cache) + return index_map, false end -function MOI.get(optimizer::Optimizer, ::MOI.SolveTime) +function MOI.get(optimizer::Optimizer, ::MOI.SolveTimeSec) return optimizer.sol.solve_time_sec end @@ -364,6 +314,15 @@ function MOI.get(optimizer::Optimizer, ::MOI.RawStatusString) return optimizer.sol.raw_status end +""" + ADMMIterations() + +The number of ADMM iterations completed during the solve. +""" +struct ADMMIterations <: MOI.AbstractModelAttribute end + +MOI.is_set_by_optimize(::ADMMIterations) = true + function MOI.get(optimizer::Optimizer, ::ADMMIterations) return optimizer.sol.iterations end @@ -426,17 +385,14 @@ function MOI.get(optimizer::Optimizer, attr::MOI.DualObjectiveValue) end function MOI.get(optimizer::Optimizer, attr::MOI.PrimalStatus) - if attr.N > MOI.get(optimizer, MOI.ResultCount()) + if attr.result_index > MOI.get(optimizer, MOI.ResultCount()) return MOI.NO_SOLUTION + elseif optimizer.sol.ret_val in (-3, 1, 2) + return MOI.FEASIBLE_POINT + elseif optimizer.sol.ret_val in (-6, -1) + return MOI.INFEASIBILITY_CERTIFICATE end - s = optimizer.sol.ret_val - if s in (-3, 1, 2) - MOI.FEASIBLE_POINT - elseif s in (-6, -1) - MOI.INFEASIBILITY_CERTIFICATE - else - MOI.INFEASIBLE_POINT - end + return MOI.INFEASIBLE_POINT end function MOI.get( @@ -445,47 +401,27 @@ function MOI.get( vi::MOI.VariableIndex, ) MOI.check_result_index_bounds(optimizer, attr) - return optimizer.data.primal[vi.value] -end - -function MOI.get( - optimizer::Optimizer, - attr::MOI.VariablePrimal, - vi::Vector{MOI.VariableIndex}, -) - return MOI.get.(optimizer, attr, vi) -end - -function post_process_result(x, ::Type{MOI.PositiveSemidefiniteConeTriangle}) - return unscalecoef(sympackedLtoU(x)) + return optimizer.sol.primal[vi.value] end -post_process_result(x, ::Type) = x - function MOI.get( optimizer::Optimizer, attr::MOI.ConstraintPrimal, ci::MOI.ConstraintIndex{F,S}, ) where {F,S} MOI.check_result_index_bounds(optimizer, attr) - return post_process_result( - optimizer.data.slack[rows(optimizer.data, ci)], - S, - ) + return optimizer.sol.slack[MOI.Utilities.rows(optimizer.cones, ci)] end function MOI.get(optimizer::Optimizer, attr::MOI.DualStatus) - if attr.N > MOI.get(optimizer, MOI.ResultCount()) + if attr.result_index > MOI.get(optimizer, MOI.ResultCount()) return MOI.NO_SOLUTION + elseif optimizer.sol.ret_val in (-3, 1, 2) + return MOI.FEASIBLE_POINT + elseif optimizer.sol.ret_val in (-7, -2) + return MOI.INFEASIBILITY_CERTIFICATE end - s = optimizer.sol.ret_val - if s in (-3, 1, 2) - MOI.FEASIBLE_POINT - elseif s in (-7, -2) - MOI.INFEASIBILITY_CERTIFICATE - else - MOI.INFEASIBLE_POINT - end + return MOI.INFEASIBLE_POINT end function MOI.get( @@ -494,7 +430,7 @@ function MOI.get( ci::MOI.ConstraintIndex{F,S}, ) where {F,S} MOI.check_result_index_bounds(optimizer, attr) - return post_process_result(optimizer.data.dual[rows(optimizer.data, ci)], S) + return optimizer.sol.dual[MOI.Utilities.rows(optimizer.cones, ci)] end MOI.get(optimizer::Optimizer, ::MOI.ResultCount) = 1 diff --git a/src/MOI_wrapper/geometric.jl b/src/MOI_wrapper/geometric.jl deleted file mode 100644 index cb057d09..00000000 --- a/src/MOI_wrapper/geometric.jl +++ /dev/null @@ -1,305 +0,0 @@ -mutable struct GeometricConicForm{T,AT,C} <: MOI.ModelLike - cone_types::C - cone_types_dict::Dict{DataType,Int} - num_rows::Vector{Int} - dimension::Dict{Int,Int} - A::Union{Nothing,AT} # The constraints are - b::Vector{T} # `b - Ax in cones` - sense::MOI.OptimizationSense - objective_constant::T # The objective is - c::Vector{T} # `sense c'x + objective_constant` - primal::Vector{T} - slack::Vector{T} - dual::Vector{T} - function GeometricConicForm{T,AT}(cone_types) where {T,AT} - model = new{T,AT,typeof(cone_types)}() - model.cone_types = cone_types - model.cone_types_dict = - Dict{DataType,Int}(s => i for (i, s) in enumerate(cone_types)) - model.num_rows = zeros(Int, length(cone_types)) - model.dimension = Dict{Int,Int}() - model.A = nothing - model.primal = T[] - model.slack = T[] - model.dual = T[] - return model - end -end - -MOI.is_empty(model::GeometricConicForm) = model.A === nothing - -function MOI.empty!(model::GeometricConicForm{T}) where {T} - empty!(model.dimension) - fill!(model.num_rows, 0) - model.A = nothing - model.sense = MOI.FEASIBILITY_SENSE - model.objective_constant = zero(T) - return -end - -function MOI.supports_constraint( - model::GeometricConicForm, - ::Type{MOI.VectorAffineFunction{Float64}}, - ::Type{S}, -) where {S<:MOI.AbstractVectorSet} - return haskey(model.cone_types_dict, S) -end - -function _allocate_variables( - model::GeometricConicForm{T,AT}, - vis_src, - idxmap, -) where {T,AT} - model.A = AT(length(vis_src)) - for (i, vi) in enumerate(vis_src) - idxmap[vi] = MOI.VariableIndex(i) - end - return -end - -function rows( - model::GeometricConicForm, - ci::MOI.ConstraintIndex{MOI.VectorAffineFunction{Float64}}, -) - return ci.value .+ (1:model.dimension[ci.value]) -end - -function MOI.supports( - ::GeometricConicForm, - ::MOI.VariablePrimalStart, - ::Type{MOI.VariableIndex}, -) - return true -end - -function MOI.set( - ::GeometricConicForm, - ::MOI.VariablePrimalStart, - ::MOI.VariableIndex, - ::Nothing, -) - return -end - -function MOI.set( - model::GeometricConicForm, - ::MOI.VariablePrimalStart, - vi::MOI.VariableIndex, - value::Float64, -) - model.primal[vi.value] = value - return -end - -function MOI.supports( - ::GeometricConicForm, - ::MOI.ConstraintPrimalStart, - ::Type{<:MOI.ConstraintIndex}, -) - return true -end - -function MOI.set( - ::GeometricConicForm, - ::MOI.ConstraintPrimalStart, - ::MOI.ConstraintIndex, - ::Nothing, -) - return -end - -function MOI.set( - model::GeometricConicForm, - ::MOI.ConstraintPrimalStart, - ci::MOI.ConstraintIndex, - value, -) - # TODO(odow): why is offset not used? - offset = constroffset(model, ci) - model.slack[rows(model, ci)] .= value - return -end - -function MOI.supports( - ::GeometricConicForm, - ::MOI.ConstraintDualStart, - ::Type{<:MOI.ConstraintIndex}, -) - return true -end - -function MOI.set( - ::GeometricConicForm, - ::MOI.ConstraintDualStart, - ::MOI.ConstraintIndex, - ::Nothing, -) - return -end - -function MOI.set( - model::GeometricConicForm, - ::MOI.ConstraintDualStart, - ci::MOI.ConstraintIndex, - value, -) - # TODO(odow): why is offset not used? - offset = constroffset(model, ci) - model.dual[rows(model, ci)] .= value - return -end - -function MOI.set( - model::GeometricConicForm, - ::MOI.ObjectiveSense, - sense::MOI.OptimizationSense, -) - model.sense = sense - return -end - -variable_index_value(t::MOI.ScalarAffineTerm) = t.variable_index.value - -function variable_index_value(t::MOI.VectorAffineTerm) - return variable_index_value(t.scalar_term) -end - -function MOI.set( - model::GeometricConicForm, - ::MOI.ObjectiveFunction, - f::MOI.ScalarAffineFunction{Float64}, -) - c = Vector( - SparseArrays.sparsevec( - variable_index_value.(f.terms), - MOI.coefficient.(f.terms), - model.A.n, - ), - ) - model.objective_constant = f.constant - model.c = c - return -end - -function _allocate_constraint( - model::GeometricConicForm, - src, - indexmap, - cone_id, - ci, -) - # TODO use `CanonicalConstraintFunction` - func = MOI.get(src, MOI.ConstraintFunction(), ci) - func = - MOI.Utilities.is_canonical(func) ? func : MOI.Utilities.canonical(func) - allocate_terms(model.A, indexmap, func) - offset = model.num_rows[cone_id] - model.num_rows[cone_id] = offset + MOI.output_dimension(func) - return ci, offset, func -end - -function _allocate_constraints( - model::GeometricConicForm, - src, - indexmap, - cone_id, - ::Type{S}, -) where {S} - cis = MOI.get( - src, - MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64},S}(), - ) - return map(cis) do ci - return _allocate_constraint(model, src, indexmap, cone_id, ci) - end -end - -function _load_variables(model::GeometricConicForm, nvars::Integer) - m = sum(model.num_rows) - model.A.m = m - model.b = zeros(m) - model.c = zeros(model.A.n) - allocate_nonzeros(model.A) - # `model.primal` contains the result of the previous optimization. - # It is used as a warm start if its length is the same, e.g. - # probably because no variable and/or constraint has been added. - if length(model.primal) != nvars - model.primal = zeros(nvars) - end - @assert length(model.dual) == length(model.slack) - if length(model.dual) != m - model.dual = zeros(m) - model.slack = zeros(m) - end - return -end - -function _load_constraints( - model::GeometricConicForm, - src, - indexmap, - cone_offset, - i, - cache, - preprocess, -) - for (ci_src, offset_in_cone, func) in cache - offset = cone_offset + offset_in_cone - set = MOI.get(src, MOI.ConstraintSet(), ci_src) - new_func = preprocess(func, set) - load_terms(model.A, indexmap, new_func, offset) - copyto!(model.b, offset + 1, new_func.constants) - model.dimension[offset] = MOI.output_dimension(func) - indexmap[ci_src] = typeof(ci_src)(offset) - end -end - -preprocess(func, set) = func - -function MOI.copy_to( - dest::GeometricConicForm, - src::MOI.ModelLike; - preprocess = preprocess, - copy_names::Bool = true, -) - MOI.empty!(dest) - - vis_src = MOI.get(src, MOI.ListOfVariableIndices()) - idxmap = MOI.Utilities.IndexMap() - - has_constraints = BitSet() - for (F, S) in MOI.get(src, MOI.ListOfConstraints()) - i = get(dest.cone_types_dict, S, nothing) - if i === nothing || F != MOI.VectorAffineFunction{Float64} - throw(MOI.UnsupportedConstraint{F,S}()) - end - push!(has_constraints, i) - end - - _allocate_variables(dest, vis_src, idxmap) - - # Allocate constraints - caches = map(collect(has_constraints)) do i - return _allocate_constraints(dest, src, idxmap, i, dest.cone_types[i]) - end - - # Load variables - _load_variables(dest, length(vis_src)) - - # Set variable attributes - MOI.Utilities.pass_attributes(dest, src, copy_names, idxmap, vis_src) - - # Set model attributes - MOI.Utilities.pass_attributes(dest, src, copy_names, idxmap) - - # Load constraints - offset = 0 - for (i, cache) in zip(has_constraints, caches) - _load_constraints(dest, src, idxmap, offset, i, cache, preprocess) - offset += dest.num_rows[i] - end - - final_touch(dest.A) - - return idxmap -end diff --git a/src/MOI_wrapper/matrix.jl b/src/MOI_wrapper/matrix.jl deleted file mode 100644 index d392d255..00000000 --- a/src/MOI_wrapper/matrix.jl +++ /dev/null @@ -1,59 +0,0 @@ -mutable struct SparseMatrixCSRtoCSC{T} <: AbstractMatrix{Cdouble} - m::Int # Number of rows - n::Int # Number of columns - colptr::Vector{T} - rowval::Vector{T} - nzval::Vector{Cdouble} - function SparseMatrixCSRtoCSC{T}(n) where {T} - A = new() - A.n = n - A.colptr = zeros(T, n + 1) - return A - end -end - -function _to_sparse(::Type{T}, A::SparseMatrixCSRtoCSC{T}) where {T} - return A.nzval, A.rowval, A.colptr -end - -function allocate_nonzeros(A::SparseMatrixCSRtoCSC{T}) where {T} - for i in 3:length(A.colptr) - A.colptr[i] += A.colptr[i-1] - end - A.rowval = Vector{T}(undef, A.colptr[end]) - A.nzval = Vector{Float64}(undef, A.colptr[end]) - return -end - -function final_touch(A::SparseMatrixCSRtoCSC) - for i in length(A.colptr):-1:2 - A.colptr[i] = A.colptr[i-1] - end - A.colptr[1] = 0 - return -end - -function _allocate_terms(colptr, indexmap, terms) - for term in terms - colptr[indexmap[term.scalar_term.variable_index].value+1] += 1 - end - return -end - -function allocate_terms(A::SparseMatrixCSRtoCSC, indexmap, func) - return _allocate_terms(A.colptr, indexmap, func.terms) -end - -function _load_terms(colptr, rowval, nzval, indexmap, terms, offset) - for term in terms - ptr = colptr[indexmap[term.scalar_term.variable_index].value] += 1 - rowval[ptr] = offset + term.output_index - 1 - nzval[ptr] = -term.scalar_term.coefficient - end - return -end - -function load_terms(A::SparseMatrixCSRtoCSC, indexmap, func, offset) - _load_terms(A.colptr, A.rowval, A.nzval, indexmap, func.terms, offset) - return -end diff --git a/src/MOI_wrapper/scaled_psd_cone_bridge.jl b/src/MOI_wrapper/scaled_psd_cone_bridge.jl new file mode 100644 index 00000000..695e143d --- /dev/null +++ b/src/MOI_wrapper/scaled_psd_cone_bridge.jl @@ -0,0 +1,144 @@ +struct ScaledPSDCone <: MOI.AbstractVectorSet + side_dimension::Int +end + +function MOI.Utilities.set_with_dimension(::Type{ScaledPSDCone}, dim) + return ScaledPSDCone(div(-1 + isqrt(1 + 8 * dim), 2)) +end + +Base.copy(x::ScaledPSDCone) = ScaledPSDCone(x.side_dimension) + +MOI.side_dimension(x::ScaledPSDCone) = x.side_dimension + +function MOI.dimension(x::ScaledPSDCone) + return div(x.side_dimension * (x.side_dimension + 1), 2) +end + +struct ScaledPSDConeBridge{T} <: MOI.Bridges.Constraint.SetMapBridge{ + T, + ScaledPSDCone, + MOI.PositiveSemidefiniteConeTriangle, + MOI.VectorAffineFunction{T}, + MOI.VectorAffineFunction{T}, +} + constraint::MOI.ConstraintIndex{MOI.VectorAffineFunction{T},ScaledPSDCone} +end + +function MOI.Bridges.Constraint.concrete_bridge_type( + ::Type{ScaledPSDConeBridge{T}}, + F::Type{<:MOI.AbstractVectorFunction}, + ::Type{MOI.PositiveSemidefiniteConeTriangle}, +) where {T} + return ScaledPSDConeBridge{T} +end + +function MOI.Bridges.map_set( + ::Type{<:ScaledPSDConeBridge}, + set::MOI.PositiveSemidefiniteConeTriangle, +) + return ScaledPSDCone(set.side_dimension) +end + +function MOI.Bridges.inverse_map_set( + ::Type{<:ScaledPSDConeBridge}, + set::ScaledPSDCone, +) + return MOI.PositiveSemidefiniteConeTriangle(set.side_dimension) +end + +function _upper_to_lower_triangular_permutation(dim::Int) + side_dimension = MOI.Utilities.side_dimension_for_vectorized_dimension(dim) + permutation = zeros(Int, dim) + i = 0 + for row in 1:side_dimension + start = div(row * (row + 1), 2) + for col in row:side_dimension + i += 1 + permutation[i] = start + start += col + end + end + return sortperm(permutation), permutation +end + +function _transform_function( + func::MOI.VectorAffineFunction{T}, + scale, + moi_to_scs::Bool, +) where {T} + d = MOI.output_dimension(func) + # upper_to_lower[i] maps the i'th element of the upper matrix to the linear + # index of the lower + # lower_to_upper[i] maps the i'th element of the lower matrix to the linear + # index of the upper + upper_to_lower, lower_to_upper = _upper_to_lower_triangular_permutation(d) + # scale_factor[i] is the amount to scale the i'th linear index in the MOI + # function by. + scale_factor = fill(scale, d) + for i in 1:d + if MOI.Utilities.is_diagonal_vectorized_index(i) + scale_factor[i] = 1.0 + end + end + scaled_constants = if moi_to_scs + (func.constants.*scale_factor)[lower_to_upper] + else + func.constants[upper_to_lower] .* scale_factor + end + scaled_terms = MOI.VectorAffineTerm{T}[] + for term in func.terms + row = term.output_index + i = moi_to_scs ? upper_to_lower[row] : lower_to_upper[row] + scale_i = moi_to_scs ? scale_factor[row] : scale_factor[i] + push!( + scaled_terms, + MOI.VectorAffineTerm( + i, + MOI.ScalarAffineTerm( + term.scalar_term.coefficient * scale_i, + term.scalar_term.variable, + ), + ), + ) + end + return MOI.VectorAffineFunction(scaled_terms, scaled_constants) +end + +function _transform_function(func::Vector{T}, scale, moi_to_scs::Bool) where {T} + d = length(func) + upper_to_lower, lower_to_upper = _upper_to_lower_triangular_permutation(d) + scale_factor = fill(scale, d) + for i in 1:d + if MOI.Utilities.is_diagonal_vectorized_index(i) + scale_factor[i] = 1.0 + end + end + if moi_to_scs + return (func.*scale_factor)[lower_to_upper] + else + return func[upper_to_lower] .* scale_factor + end +end + +# Map ConstraintFunction from MOI -> SCS +function MOI.Bridges.map_function(::Type{<:ScaledPSDConeBridge}, f) + return _transform_function(f, √2, true) +end + +# Used to map the ConstraintPrimal from SCS -> MOI +function MOI.Bridges.inverse_map_function(::Type{<:ScaledPSDConeBridge}, f) + return _transform_function(f, 1 / √2, false) +end + +# Used to map the ConstraintDual from SCS -> MOI +function MOI.Bridges.adjoint_map_function(::Type{<:ScaledPSDConeBridge}, f) + return _transform_function(f, 1 / √2, false) +end + +# Used to set ConstraintDualStart +function MOI.Bridges.inverse_adjoint_map_function( + ::Type{<:ScaledPSDConeBridge}, + f, +) + return _transform_function(f, √2, true) +end diff --git a/src/SCS.jl b/src/SCS.jl index 2704d7ba..e632c0cc 100644 --- a/src/SCS.jl +++ b/src/SCS.jl @@ -42,6 +42,6 @@ include("MOI_wrapper/MOI_wrapper.jl") const available_solvers = [DirectSolver, IndirectSolver] -export SCS_solve +export scs_solve end diff --git a/src/c_wrapper.jl b/src/c_wrapper.jl index 290b0aee..4d9f8150 100644 --- a/src/c_wrapper.jl +++ b/src/c_wrapper.jl @@ -8,7 +8,7 @@ function Base.unsafe_convert(::Type{Ptr{Cvoid}}, x::AbstractSCSType) return pointer_from_objref(x) end -mutable struct SCSMatrix{T} <: AbstractSCSType +mutable struct ScsMatrix{T} <: AbstractSCSType values::Ptr{Cdouble} rowval::Ptr{T} colptr::Ptr{T} @@ -16,7 +16,7 @@ mutable struct SCSMatrix{T} <: AbstractSCSType n::T end -mutable struct SCSSettings{T} <: AbstractSCSType +mutable struct ScsSettings{T} <: AbstractSCSType normalize::T # boolean, heuristic data rescaling scale::Cdouble # if normalized, rescales by this factor rho_x::Cdouble # x equality constraint scaling @@ -28,10 +28,10 @@ mutable struct SCSSettings{T} <: AbstractSCSType warm_start::T # boolean, warm start (put initial guess in Sol struct) acceleration_lookback::T # acceleration memory parameter write_data_filename::Cstring - SCSSettings{T}() where {T} = new{T}() + ScsSettings{T}() where {T} = new{T}() end -mutable struct SCSData{T} <: AbstractSCSType +mutable struct ScsData{T} <: AbstractSCSType m::T n::T A::Ptr{Cvoid} @@ -40,13 +40,13 @@ mutable struct SCSData{T} <: AbstractSCSType stgs::Ptr{Cvoid} end -mutable struct SCSSolution <: AbstractSCSType +mutable struct ScsSolution <: AbstractSCSType x::Ptr{Cdouble} y::Ptr{Cdouble} s::Ptr{Cdouble} end -mutable struct SCSInfo{T} <: AbstractSCSType +mutable struct ScsInfo{T} <: AbstractSCSType iter::T status::NTuple{32,Cchar} statusVal::T @@ -59,10 +59,10 @@ mutable struct SCSInfo{T} <: AbstractSCSType relGap::Cdouble setupTime::Cdouble solveTime::Cdouble - SCSInfo{T}() where {T} = new{T}() + ScsInfo{T}() where {T} = new{T}() end -mutable struct SCSCone{T} <: AbstractSCSType +mutable struct ScsCone{T} <: AbstractSCSType f::T l::T q::Ptr{T} @@ -79,23 +79,23 @@ mutable struct Solution{T} x::Array{Float64,1} y::Array{Float64,1} s::Array{Float64,1} - info::SCSInfo{T} + info::ScsInfo{T} ret_val::T end """ - _SCSDataWrapper + _ScsDataWrapper A type for wrapping all data inputs to SCS to prevent them from being freed by the garbage collector during a solve. -You should not construct this manually. Call `SCS_solve` instead. +You should not construct this manually. Call `scs_solve` instead. """ -struct _SCSDataWrapper{S,T} +struct _ScsDataWrapper{S,T} linear_solver::S m::T n::T - A::SCSMatrix{T} + A::ScsMatrix{T} values::Vector{Cdouble} rowval::Vector{T} colptr::Vector{T} @@ -111,13 +111,13 @@ struct _SCSDataWrapper{S,T} primal::Vector{Cdouble} dual::Vector{Cdouble} slack::Vector{Cdouble} - settings::SCSSettings + settings::ScsSettings options::Any end function _sanitize_options(options) option_dict = Dict{Symbol,Any}() - fields = fieldnames(SCSSettings) + fields = fieldnames(ScsSettings) for (key, value) in options if key == :linear_solver continue @@ -129,7 +129,7 @@ function _sanitize_options(options) return option_dict end -function raw_status(info::SCSInfo) +function raw_status(info::ScsInfo) data = UInt8[info.status[i] for i in 1:findfirst(iszero, info.status)-1] return String(data) end @@ -145,7 +145,7 @@ function _to_sparse(::Type{T}, A::SparseArrays.SparseMatrixCSC) where {T} end """ - SCS_solve(args...) + scs_solve(args...) SCS solves a problem of the form ``` @@ -180,12 +180,18 @@ above. the dual cone) Returns a Solution object. + +!!! warning + SCS expects the semi-definite cones to be scaled by a factor of √2. That is, + the off-diagonal elements in the A matrix and b vector should be multiplied + by √2, and then the corresponding rows in the dual and slack solution + vectors should be multiplied by 1/√2. """ -function SCS_solve( +function scs_solve( linear_solver::Type{<:LinearSolver}, m::Integer, n::Integer, - A::AbstractMatrix, + A, b::Vector{Float64}, c::Vector{Float64}, f::Integer, @@ -225,11 +231,11 @@ function SCS_solve( if warm_start option_dict[:warm_start] = 1 end - model = _SCSDataWrapper( + model = _ScsDataWrapper( linear_solver, m, n, - SCSMatrix(pointer(values), pointer(rowval), pointer(colptr), m, n), + ScsMatrix(pointer(values), pointer(rowval), pointer(colptr), m, n), values, rowval, colptr, @@ -245,7 +251,7 @@ function SCS_solve( primal_sol, dual_sol, slack, - SCSSettings{T}(), + ScsSettings{T}(), option_dict, ) Base.GC.@preserve model begin @@ -253,13 +259,13 @@ function SCS_solve( end end -function _unsafe_scs_solve(model::_SCSDataWrapper{S,T}) where {S,T} - scs_solution = SCSSolution( +function _unsafe_scs_solve(model::_ScsDataWrapper{S,T}) where {S,T} + scs_solution = ScsSolution( pointer(model.primal), pointer(model.dual), pointer(model.slack), ) - scs_cone = SCSCone{T}( + scs_cone = ScsCone{T}( model.f, model.l, pointer(model.q), @@ -271,8 +277,8 @@ function _unsafe_scs_solve(model::_SCSDataWrapper{S,T}) where {S,T} pointer(model.p), length(model.p), ) - scs_info = SCSInfo{T}() - scs_data = SCSData{T}( + scs_info = ScsInfo{T}() + scs_data = ScsData{T}( model.m, model.n, pointer_from_objref(model.A), @@ -285,7 +291,7 @@ function _unsafe_scs_solve(model::_SCSDataWrapper{S,T}) where {S,T} setfield!( model.settings, key, - convert(fieldtype(SCSSettings{T}, key), value), + convert(fieldtype(ScsSettings{T}, key), value), ) end p = scs_init(model.linear_solver, scs_data, scs_cone, scs_info) diff --git a/src/linear_solvers/direct.jl b/src/linear_solvers/direct.jl index 3d53c905..226ca80e 100644 --- a/src/linear_solvers/direct.jl +++ b/src/linear_solvers/direct.jl @@ -2,7 +2,7 @@ struct DirectSolver <: LinearSolver end scsint_t(::Type{DirectSolver}) = Int -function scs_set_default_settings(::Type{DirectSolver}, data::SCSData{Int}) +function scs_set_default_settings(::Type{DirectSolver}, data::ScsData{Int}) return ccall( (:scs_set_default_settings, direct), Cvoid, @@ -13,9 +13,9 @@ end function scs_init( ::Type{DirectSolver}, - data::SCSData{Int}, - cone::SCSCone{Int}, - info::SCSInfo{Int}, + data::ScsData{Int}, + cone::ScsCone{Int}, + info::ScsInfo{Int}, ) return ccall( (:scs_init, direct), @@ -30,10 +30,10 @@ end function scs_solve( ::Type{DirectSolver}, p_work::Ptr{Cvoid}, - data::SCSData{Int}, - cone::SCSCone{Int}, - solution::SCSSolution, - info::SCSInfo{Int}, + data::ScsData{Int}, + cone::ScsCone{Int}, + solution::ScsSolution, + info::ScsInfo{Int}, ) return ccall( (:scs_solve, direct), diff --git a/src/linear_solvers/gpu_indirect.jl b/src/linear_solvers/gpu_indirect.jl index a034aa3b..60103166 100644 --- a/src/linear_solvers/gpu_indirect.jl +++ b/src/linear_solvers/gpu_indirect.jl @@ -12,7 +12,7 @@ scsint_t(::Type{GpuIndirectSolver}) = Int32 function scs_set_default_settings( ::Type{GpuIndirectSolver}, - data::SCSData{Int32}, + data::ScsData{Int32}, ) return ccall( (:scs_set_default_settings, gpuindirect), @@ -24,9 +24,9 @@ end function scs_init( ::Type{GpuIndirectSolver}, - data::SCSData{Int32}, - cone::SCSCone{Int32}, - info::SCSInfo{Int32}, + data::ScsData{Int32}, + cone::ScsCone{Int32}, + info::ScsInfo{Int32}, ) return ccall( (:scs_init, gpuindirect), @@ -41,10 +41,10 @@ end function scs_solve( ::Type{GpuIndirectSolver}, p_work::Ptr{Cvoid}, - data::SCSData{Int32}, - cone::SCSCone{Int32}, - solution::SCSSolution, - info::SCSInfo{Int32}, + data::ScsData{Int32}, + cone::ScsCone{Int32}, + solution::ScsSolution, + info::ScsInfo{Int32}, ) return ccall( (:scs_solve, gpuindirect), diff --git a/src/linear_solvers/indirect.jl b/src/linear_solvers/indirect.jl index f0ada274..d7e06622 100644 --- a/src/linear_solvers/indirect.jl +++ b/src/linear_solvers/indirect.jl @@ -2,7 +2,7 @@ struct IndirectSolver <: LinearSolver end scsint_t(::Type{IndirectSolver}) = Int -function scs_set_default_settings(::Type{IndirectSolver}, data::SCSData{Int}) +function scs_set_default_settings(::Type{IndirectSolver}, data::ScsData{Int}) return ccall( (:scs_set_default_settings, indirect), Cvoid, @@ -13,9 +13,9 @@ end function scs_init( ::Type{IndirectSolver}, - data::SCSData{Int}, - cone::SCSCone{Int}, - info::SCSInfo{Int}, + data::ScsData{Int}, + cone::ScsCone{Int}, + info::ScsInfo{Int}, ) return ccall( (:scs_init, indirect), @@ -30,10 +30,10 @@ end function scs_solve( ::Type{IndirectSolver}, p_work::Ptr{Cvoid}, - data::SCSData{Int}, - cone::SCSCone{Int}, - solution::SCSSolution, - info::SCSInfo{Int}, + data::ScsData{Int}, + cone::ScsCone{Int}, + solution::ScsSolution, + info::ScsInfo{Int}, ) return ccall( (:scs_solve, indirect), diff --git a/test/MOI_wrapper.jl b/test/MOI_wrapper.jl index fe7b33c7..de4818ae 100644 --- a/test/MOI_wrapper.jl +++ b/test/MOI_wrapper.jl @@ -1,73 +1,91 @@ -using Test +module TestSCS +using Test using MathOptInterface -const MOI = MathOptInterface -const MOIT = MOI.Test - -# UniversalFallback is needed for starting values -const CACHE = MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()) - import SCS -function moi_tests(T) - optimizer = SCS.Optimizer(linear_solver = T, eps = 1e-6) - MOI.set(optimizer, MOI.Silent(), true) - - @testset "SolverName" begin - @test MOI.get(optimizer, MOI.SolverName()) == "SCS" - end - - MOI.empty!(CACHE) - cached = MOI.Utilities.CachingOptimizer(CACHE, optimizer) - bridged = MOI.Bridges.full_bridge_optimizer(cached, Float64) - config = MOIT.TestConfig(atol = 1e-5) - - @testset "Unit" begin - MOIT.unittest( - bridged, - config, - [ - # FIXME `NumberOfThreads` not supported. - "number_threads", - # `TimeLimitSec` not supported. - "time_limit_sec", - # ArgumentError: The number of constraints in SCSModel must be greater than 0 - "solve_unbounded_model", - # Integer and ZeroOne sets are not supported - "solve_integer_edge_cases", - "solve_objbound_edge_cases", - "solve_zero_one_with_bounds_1", - "solve_zero_one_with_bounds_2", - "solve_zero_one_with_bounds_3", - ], - ) - end +const MOI = MathOptInterface - @testset "Continuous linear problems with $T" begin - MOIT.contlineartest(bridged, config) +function runtests() + for name in names(@__MODULE__; all = true) + if startswith("$(name)", "test_") + @testset "$(name)" begin + getfield(@__MODULE__, name)() + end + end end + return +end - @testset "ADMMIterations attribute with $T" begin - MOIT.linear1test(bridged, config) - @test MOI.get(bridged, SCS.ADMMIterations()) > 0 - end +test_DirectSolver() = _test_runtests(SCS.DirectSolver) - @testset "Continuous quadratic problems with $T" begin - MOIT.qcptest(bridged, config) - end +test_IndirectSolver() = _test_runtests(SCS.IndirectSolver) - @testset "Continuous conic problems with $T" begin - MOIT.contconictest(bridged, config, ["rootdets", "logdets"]) - end +function _test_runtests(linear_solver) + optimizer = SCS.Optimizer() + MOI.set( + optimizer, + MOI.RawOptimizerAttribute("linear_solver"), + linear_solver, + ) + MOI.set(optimizer, MOI.RawOptimizerAttribute("eps"), 1e-6) + MOI.set(optimizer, MOI.Silent(), true) + model = MOI.Bridges.full_bridge_optimizer( + MOI.Utilities.CachingOptimizer( + MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()), + optimizer, + ), + Float64, + ) + MOI.Test.runtests( + model, + MOI.Test.Config( + atol = 1e-4, + exclude = Any[ + MOI.ConstraintBasisStatus, + MOI.VariableBasisStatus, + MOI.ConstraintName, + MOI.VariableName, + MOI.ObjectiveBound, + ], + ), + exclude = String[ + # Expected test failures: + # ArgumentError: The number of constraints must be greater than 0 + "test_attribute_RawStatusString", + "test_attribute_SolveTimeSec", + "test_objective_ObjectiveFunction_blank", + "test_solve_TerminationStatus_DUAL_INFEASIBLE", + # Problem is a nonconvex QP + "test_basic_ScalarQuadraticFunction_EqualTo", + "test_basic_ScalarQuadraticFunction_GreaterThan", + "test_basic_ScalarQuadraticFunction_Interval", + "test_basic_VectorQuadraticFunction_", + "test_quadratic_SecondOrderCone_basic", + "test_quadratic_nonconvex_", + # power cone error, values must be in [-1,1] + "test_conic_DualPowerCone_VectorOfVariables", + "test_conic_DualPowerCone_VectorAffineFunction", + "test_conic_PowerCone_VectorAffineFunction", + "test_conic_PowerCone_VectorOfVariables", + # MathOptInterface.jl issue #1431 + "test_model_LowerBoundAlreadySet", + "test_model_UpperBoundAlreadySet", + ], + ) + return end -@testset "MOI.RawParameter" begin +function test_RawOptimizerAttribute() model = SCS.Optimizer() - # TODO(odow): remove symbol cases when deprecation is removed. - MOI.set(model, MOI.RawParameter(:eps), 1.0) - @test MOI.get(model, MOI.RawParameter(:eps)) == 1.0 - @test MOI.get(model, MOI.RawParameter("eps")) == 1.0 - MOI.set(model, MOI.RawParameter("eps"), 2.0) - @test MOI.get(model, MOI.RawParameter(:eps)) == 2.0 - @test MOI.get(model, MOI.RawParameter("eps")) == 2.0 + MOI.set(model, MOI.RawOptimizerAttribute("eps"), 1.0) + @test MOI.get(model, MOI.RawOptimizerAttribute("eps")) == 1.0 + @test MOI.get(model, MOI.RawOptimizerAttribute("eps")) == 1.0 + MOI.set(model, MOI.RawOptimizerAttribute("eps"), 2.0) + @test MOI.get(model, MOI.RawOptimizerAttribute("eps")) == 2.0 + @test MOI.get(model, MOI.RawOptimizerAttribute("eps")) == 2.0 end + +end # module + +TestSCS.runtests() diff --git a/test/runtests.jl b/test/runtests.jl index 2a695109..90304e02 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -9,12 +9,8 @@ using Test using SCS include("test_problems.jl") -include("MOI_wrapper.jl") - for s in SCS.available_solvers feasible_basic_problems(s) end -for s in SCS.available_solvers - moi_tests(s) -end +include("MOI_wrapper.jl") diff --git a/test/test_gpu.jl b/test/test_gpu.jl index f4ae8662..949e5b0e 100644 --- a/test/test_gpu.jl +++ b/test/test_gpu.jl @@ -15,5 +15,6 @@ using Test include("test_problems.jl") feasible_basic_problems(SCS.GpuIndirectSolver) -include("MOI_wrapper.jl") -moi_tests(SCS.GpuIndirectSolver) +# TODO(odow): consider re-enabling these +# include("MOI_wrapper.jl") +# moi_tests(SCS.GpuIndirectSolver) diff --git a/test/test_problems.jl b/test/test_problems.jl index 982d2a6b..64df2c5f 100644 --- a/test/test_problems.jl +++ b/test/test_problems.jl @@ -432,7 +432,7 @@ function feasible_basic_conic(T) A = SparseMatrixCSC(m, n, colptr .+ 1, rowval .+ 1, vec(values)) - sol = SCS_solve(T, m, n, A, b, c, f, l, q, s, ep, ed, Float64[]) + sol = scs_solve(T, m, n, A, b, c, f, l, q, s, ep, ed, Float64[]) @test sol.ret_val == 1 @test SCS.raw_status(sol.info) == "Solved" end @@ -556,7 +556,7 @@ function feasible_exponential_conic(T) A = SparseMatrixCSC(m, n, colptr .+ 1, rowval .+ 1, vec(values)) - sol = SCS_solve(T, m, n, A, b, c, f, l, q, s, ep, ed, Float64[]) + sol = scs_solve(T, m, n, A, b, c, f, l, q, s, ep, ed, Float64[]) @assert sol.ret_val == 1 @test SCS.raw_status(sol.info) == "Solved" end @@ -633,7 +633,7 @@ function feasible_sdp_conic(T) A = SparseMatrixCSC(m, n, colptr .+ 1, rowval .+ 1, vec(values)) - sol = SCS_solve(T, m, n, A, b, c, f, l, q, s, ep, ed, Float64[]) + sol = scs_solve(T, m, n, A, b, c, f, l, q, s, ep, ed, Float64[]) @test sol.ret_val == 1 @test SCS.raw_status(sol.info) == "Solved" end @@ -684,14 +684,14 @@ function feasible_pow_conic(T) A = SparseMatrixCSC(m, n, colptr .+ 1, rowval .+ 1, vec(values)) - sol = SCS_solve(T, m, n, A, b, c, f, l, q, s, ep, ed, p) + sol = scs_solve(T, m, n, A, b, c, f, l, q, s, ep, ed, p) @test sol.ret_val == 1 @test SCS.raw_status(sol.info) == "Solved" end function feasible_basic_problems(solver) A = reshape([1.0],(1,1)) - solution = SCS_solve(solver, 1, 1, A, [1.0], [1.0], 1, 0, Int[], Int[], 0, 0, Float64[]); + solution = scs_solve(solver, 1, 1, A, [1.0], [1.0], 1, 0, Int[], Int[], 0, 0, Float64[]); @test solution.ret_val == 1 feasible_basic_conic(solver)