JuMP, moi wrapper, interface problem

zhenweilin · October 26, 2024, 6:27am

JuMP, moi wrapper, interface problem.

I am currently developing an SOCP solver designed to address the following optimization problem.

min_{x=(x_1,x_2),l \leq x_1 \leq u, x_2\in\mathcal{K}}c^\top x \text{ s.t.}Gx-h\in \mathcal{K}_G, With the capabilities of the SCS solver, we were able to adapt this formulation by making minor adjustments to its jump interface. However, we encountered the following issues during its application. When we activate VectorOfVariables with the following code,

# # customized
function MOI.supports_constraint(
    ::Optimizer,
    ::Type{MOI.VectorOfVariables},
    ::Type{
        <:Union{
            MOI.SecondOrderCone,
            MOI.RotatedSecondOrderCone,
            # may be not supported for two exponential cones
            MOI.ExponentialCone,
            MOI.DualExponentialCone,
        }
    }
)
    return true
end

Our variable may be subject to both a box-bound constraint and a cone constraint, which may not align with our intended modeling approach. Alternatively, if x is simultaneously modeled in two cones, we aim to retain only one cone constraint, incorporating the other directly into the constraint matrix G. I present the following code to demonstrate the issues encountered. My modeling code is as follows:

using Pkg
Pkg.activate("test_cpu_env")
include("../../src/rpdhg_clp_cpu/RPDHG_CLP_CPU.jl")
using .RPDHG_CLP_CPU
using Arpack, LinearAlgebra
using JuMP, MosekTools
using Profile
using Random, SparseArrays
using SCS

model = JuMP.Model(RPDHG_CLP_CPU.Optimizer);

set_optimizer_attribute(model, "use_scaling", true)
# set_optimizer_attribute(model, "rescaling_method", :ruiz)
set_optimizer_attribute(model, "use_adaptive_restart", true)
set_optimizer_attribute(model, "use_adaptive_step_size_weight", true)
set_optimizer_attribute(model, "use_resolving", true)
set_optimizer_attribute(model, "use_accelerated", false)
set_optimizer_attribute(model, "time_limit_secs", 10000.0)
set_optimizer_attribute(model, "verbose", 2)
n = 10
k_box = 2
k_soc1 = 4
k_soc2 = 4
c = rand(n)
c = abs.(c)
@variable(model, y[1:n] >= 1);
@objective(model, Min, c' * y);
x_soc1 = @constraint(model, [y[k_box + 1]; y[k_box + 2: k_box + k_soc1]] in SecondOrderCone());
x_soc2 = @constraint(model, [y[k_box + 1]; y[k_box + 2: k_box + k_soc1 + 1]]in SecondOrderCone());
x_soc3 = @constraint(model, [2 * y[k_box + 1]; y[k_box + 2: k_box + k_soc1 + 1]]in SecondOrderCone());
optimize!(model)

In this example, I model the variable yyy, with y[4] serving as a representative case. This variable is subject to a box constraint and three cone constraints: x_soc1, x_soc2, and x_soc3. Here, I intend to incorporate all three cone constraints into the constraint matrix G.

However, after enabling VectorOfVariables in the SCS jump interface, not all cone constraints are incorporated into G, leading to the following output:

This indicates that two of the cone constraints (x_soc1, x_soc2) are retained.

Below is the full implementation of my MOI_wrapper code.

import MathOptInterface as MOI
"""
    struct ScaledPSDCone <: MOI.AbstractVectorSet
        side_dimension::Int
    end

Similar to `MOI.Scaled{MOI.PositiveSemidefiniteConeTriangle}` but it the
vectorization is the lower triangular part column-wise (or the upper triangular
part row-wise).
"""
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,F} <: MOI.Bridges.Constraint.SetMapBridge{
    T,
    ScaledPSDCone,
    MOI.Scaled{MOI.PositiveSemidefiniteConeTriangle},
    F,
    F,
}
    constraint::MOI.ConstraintIndex{F,ScaledPSDCone}
end

function MOI.Bridges.Constraint.concrete_bridge_type(
    ::Type{ScaledPSDConeBridge{T}},
    ::Type{F},
    ::Type{MOI.Scaled{MOI.PositiveSemidefiniteConeTriangle}},
) where {T,F<:MOI.AbstractVectorFunction}
    return ScaledPSDConeBridge{T,F}
end

function MOI.Bridges.map_set(
    ::Type{<:ScaledPSDConeBridge},
    set::MOI.Scaled{MOI.PositiveSemidefiniteConeTriangle},
)
    return ScaledPSDCone(MOI.side_dimension(set))
end

function MOI.Bridges.inverse_map_set(
    ::Type{<:ScaledPSDConeBridge},
    set::ScaledPSDCone,
)
    return MOI.Scaled(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_to_scs::Bool)
    scalars = MOI.Utilities.eachscalar(func)
    d = length(scalars)
    upper_to_lower, lower_to_upper = _upper_to_lower_triangular_permutation(d)
    if moi_to_scs
        return scalars[lower_to_upper]
    else
        return scalars[upper_to_lower]
    end
end

# Map ConstraintFunction from MOI -> SCS
function MOI.Bridges.map_function(::Type{<:ScaledPSDConeBridge}, f)
    return _transform_function(f, true)
end

# Used to map the ConstraintPrimal from SCS -> MOI
function MOI.Bridges.inverse_map_function(::Type{<:ScaledPSDConeBridge}, f)
    return _transform_function(f, false)
end

# Used to map the ConstraintDual from SCS -> MOI
function MOI.Bridges.adjoint_map_function(::Type{<:ScaledPSDConeBridge}, f)
    return _transform_function(f, false)
end

# Used to set ConstraintDualStart
function MOI.Bridges.inverse_adjoint_map_function(
    ::Type{<:ScaledPSDConeBridge},
    f,
)
    return _transform_function(f, true)
end

MOI.Utilities.@product_of_sets(
    _Cones,
    MOI.Zeros,
    MOI.Nonnegatives,
    MOI.SecondOrderCone,
    MOI.RotatedSecondOrderCone,
    MOI.ExponentialCone,
    MOI.DualExponentialCone,
    # MOI.PowerCone{T},
    # MOI.DualPowerCone{T},
)

struct _SetConstants{T}
    b::Vector{T}
    power_coefficients::Dict{Int,T}
    _SetConstants{T}() where {T} = new{T}(T[], Dict{Int,T}())
end

function Base.empty!(x::_SetConstants)
    empty!(x.b)
    empty!(x.power_coefficients)
    return x
end

Base.resize!(x::_SetConstants, n) = resize!(x.b, n)

function MOI.Utilities.load_constants(x::_SetConstants, offset, f)
    MOI.Utilities.load_constants(x.b, offset, f)
    return
end

function MOI.Utilities.load_constants(
    x::_SetConstants{T},
    offset,
    set::Union{MOI.PowerCone{T},MOI.DualPowerCone{T}},
) where {T}
    x.power_coefficients[offset+1] = set.exponent
    return
end

function MOI.Utilities.function_constants(x::_SetConstants, rows)
    return MOI.Utilities.function_constants(x.b, rows)
end

function MOI.Utilities.set_from_constants(x::_SetConstants, S, rows)
    return MOI.Utilities.set_from_constants(x.b, S, rows)
end

function MOI.Utilities.set_from_constants(
    x::_SetConstants{T},
    ::Type{S},
    rows,
) where {T,S<:Union{MOI.PowerCone{T},MOI.DualPowerCone{T}}}
    @assert length(rows) == 3
    return S(x.power_coefficients[first(rows)])
end

const OptimizerCache{T} = MOI.Utilities.GenericModel{
    Float64,
    MOI.Utilities.ObjectiveContainer{Float64},
    MOI.Utilities.VariablesContainer{Float64},
    MOI.Utilities.MatrixOfConstraints{
        Float64,
        MOI.Utilities.MutableSparseMatrixCSC{
            Float64,
            T,
            MOI.Utilities.OneBasedIndexing,
        },
        _SetConstants{Float64},
        _Cones{Float64},
    },
}

function _to_sparse(
    ::Type{T},
    A::MOI.Utilities.MutableSparseMatrixCSC{
        Float64,
        T,
        MOI.Utilities.OneBasedIndexing,
    },
) where {T}
    return -A.nzval, A.rowval, A.colptr
end

mutable struct MOISolution
    primal::Vector{Float64}
    dual::Vector{Float64}
    slack::Vector{Float64}
    exit_code::Int
    exit_status::Symbol
    objective_value::Float64
    dual_objective_value::Float64
    solve_time_sec::Float64
    iterations::Int
    objective_constant::Float64
    function MOISolution(;primal = Float64[], dual = Float64[], slack = Float64[], exit_code = 0, exit_status = :unknown, objective_value = NaN, dual_objective_value = NaN, solve_time_sec = NaN, iterations = 0, objective_constant = 0.0)
        new(primal, dual, slack, exit_code, exit_status, objective_value, dual_objective_value, solve_time_sec, iterations, objective_constant)
    end
end

"""
    Optimizer()

Create a new PDHG-CLP optimizer.
"""
mutable struct Optimizer <: MOI.AbstractOptimizer
    cones::Union{Nothing,_Cones{Float64}}
    sol::MOISolution
    silent::Bool
    options::Dict{Symbol,Any}

    Optimizer() = new(nothing, MOISolution(), false, Dict{Symbol,Any}())
end

function MOI.get(::Optimizer, ::MOI.Bridges.ListOfNonstandardBridges)
    return [ScaledPSDConeBridge{Float64}]
end

MOI.get(::Optimizer, ::MOI.SolverName) = "PDHG-CLP"

MOI.is_empty(optimizer::Optimizer) = optimizer.cones === nothing

function MOI.empty!(optimizer::Optimizer)
    optimizer.cones = nothing
    optimizer.sol = MOISolution()
    return
end

###
### MOI.RawOptimizerAttribute
###

MOI.supports(::Optimizer, ::MOI.RawOptimizerAttribute) = true

function MOI.set(optimizer::Optimizer, param::MOI.RawOptimizerAttribute, value)
    return optimizer.options[Symbol(param.name)] = value
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)
    optimizer.silent = value
    return
end

MOI.get(optimizer::Optimizer, ::MOI.Silent) = optimizer.silent

###
### MOI.TimeLimitSec
###

MOI.supports(::Optimizer, ::MOI.TimeLimitSec) = true

function MOI.set(optimizer::Optimizer, ::MOI.TimeLimitSec, time_limit::Real)
    optimizer.options[:time_limit_secs] = convert(Float64, time_limit)
    return
end

function MOI.set(optimizer::Optimizer, ::MOI.TimeLimitSec, ::Nothing)
    delete!(optimizer.options, :time_limit_secs)
    return
end

function MOI.get(optimizer::Optimizer, ::MOI.TimeLimitSec)
    value = get(optimizer.options, :time_limit_secs, nothing)
    return value::Union{Float64,Nothing}
end

###
### MOI.AbstractModelAttribute
###

function MOI.supports(
    ::Optimizer,
    ::Union{
        MOI.ObjectiveSense,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}},
        # MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}},
    },
)
    return true
end

###
### MOI.AbstractVariableAttribute
###

function MOI.supports(
    ::Optimizer,
    ::MOI.VariablePrimalStart,
    ::Type{MOI.VariableIndex},
)
    return true
end

###
### MOI.AbstractConstraintAttribute
###

function MOI.supports(
    ::Optimizer,
    ::Union{MOI.ConstraintPrimalStart,MOI.ConstraintDualStart},
    ::Type{<:MOI.ConstraintIndex},
)
    return true
end

function MOI.supports_constraint(
    ::Optimizer,
    ::Type{MOI.VectorAffineFunction{Float64}},
    ::Type{
        <:Union{
            MOI.Zeros,
            MOI.Nonnegatives,
            MOI.SecondOrderCone,
            MOI.RotatedSecondOrderCone,
            MOI.ExponentialCone,
            MOI.DualExponentialCone,
            # MOI.PowerCone{Float64},
            # MOI.DualPowerCone{Float64},
        },
    },
)
    return true
end

# customized
function MOI.supports_constraint(
    ::Optimizer,
    ::Type{MOI.VariableIndex},
    ::Type{
        <:Union{
            MOI.EqualTo{Float64},
            MOI.LessThan{Float64},
            MOI.GreaterThan{Float64},
            MOI.Interval{Float64},
        }
    }
)
    return true
end

# customized
function MOI.supports_constraint(
    ::Optimizer,
    ::Type{MOI.VectorOfVariables},
    ::Type{
        <:Union{
            MOI.SecondOrderCone,
            MOI.RotatedSecondOrderCone,
            # may be not supported for two exponential cones
            MOI.ExponentialCone,
            MOI.DualExponentialCone,
        }
    }
)
    return true
end


function _map_sets(f, ::Type{T}, sets, ::Type{S}) where {T,S}
    F = MOI.VectorAffineFunction{Float64}
    cis = MOI.get(sets, MOI.ListOfConstraintIndices{F,S}())
    return T[f(MOI.get(sets, MOI.ConstraintSet(), ci)) for ci in cis]
end

function MOI.optimize!(
    dest::Optimizer,
    src::MOI.Utilities.UniversalFallback{OptimizerCache{T}},
) where {T}
    # The real stuff starts here.
    MOI.empty!(dest)
    index_map = MOI.Utilities.identity_index_map(src)
    Ab = src.model.constraints
    A = Ab.coefficients

    for (F, S) in keys(src.constraints) # src.constraints is very important
        # throw(MOI.UnsupportedConstraint{F,S}())
        if F == MOI.VectorOfVariables
            if S == MOI.SecondOrderCone
                println("SecondOrderCone")
            elseif S == MOI.RotatedSecondOrderCone
                println("RotatedSecondOrderCone")
            elseif S == MOI.ExponentialCone
                println("ExponentialCone")
            elseif S == MOI.DualExponentialCone
                println("DualExponentialCone")
            else
                throw(MOI.UnsupportedConstraint{F,S}())
            end
        else
            throw(MOI.UnsupportedConstraint{F,S}())
        end
    end
    # key: src.model.variables is also very important
    model_attributes = MOI.get(src, MOI.ListOfModelAttributesSet())
    max_sense = false
    obj_attr = nothing
    for attr in model_attributes
        if attr == MOI.ObjectiveSense()
            max_sense = MOI.get(src, attr) == MOI.MAX_SENSE
        elseif attr == MOI.Name()
            continue  # This can be skipped without consequence
        elseif attr isa MOI.ObjectiveFunction
            obj_attr = attr
        else
            throw(MOI.UnsupportedAttribute(attr))
        end
    end
    objective_constant, c = 0.0, zeros(A.n)
    if obj_attr == MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}()
        obj = MOI.get(src, obj_attr)
        objective_constant = MOI.constant(obj)
        for term in obj.terms
            c[term.variable.value] += (max_sense ? -1 : 1) * term.coefficient
        end
    elseif obj_attr !== nothing
        throw(MOI.UnsupportedAttribute(obj_attr))
    end

    # `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
    if A.m != length(dest.sol.dual)
        resize!(dest.sol.dual, A.m)
        fill!(dest.sol.dual, 0.0)
    end
    if A.m != length(dest.sol.slack)
        resize!(dest.sol.slack, A.m)
        fill!(dest.sol.slack, 0.0)
    end
    # Set starting values and throw error for other variable attributes
    has_warm_start = false
    vis_src = MOI.get(src, MOI.ListOfVariableIndices())
    for attr in MOI.get(src, MOI.ListOfVariableAttributesSet())
        if attr == MOI.VariableName()
            # Skip
        elseif attr == MOI.VariablePrimalStart()
            has_warm_start = true
            for (i, x) in enumerate(vis_src)
                dest.sol.primal[i] = something(MOI.get(src, attr, x), 0.0)
            end
        else
            throw(MOI.UnsupportedAttribute(attr))
        end
    end
    # 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}())
        for attr in MOI.get(src, MOI.ListOfConstraintAttributesSet{F,S}())
            if attr == MOI.ConstraintName()
                # Skip
            elseif attr == MOI.ConstraintPrimalStart()
                has_warm_start = true
                for ci in cis_src
                    start = MOI.get(src, attr, ci)
                    if start !== nothing
                        rows = MOI.Utilities.rows(Ab, ci)
                        dest.sol.slack[rows] .= start
                    end
                end
            elseif attr == MOI.ConstraintDualStart()
                has_warm_start = true
                for ci in cis_src
                    start = MOI.get(src, attr, ci)
                    if start !== nothing
                        rows = MOI.Utilities.rows(Ab, ci)
                        dest.sol.dual[rows] .= start
                    end
                end
            else
                throw(MOI.UnsupportedAttribute(attr))
            end
        end
    end
    options = copy(dest.options)
    if dest.silent
        options[:verbose] = 0
    end
    dest.cones = deepcopy(Ab.sets)
    mGzero = sum(_map_sets(MOI.dimension, T, Ab, MOI.Zeros))
    mGnonnegative = sum(_map_sets(MOI.dimension, T, Ab, MOI.Nonnegatives))
    G = SparseMatrixCSC{Float64, Int64}(A.m, A.n, Int64.(A.colptr), Int64.(A.rowval), A.nzval)
    bl = fill(-Inf, A.n)
    bl .= src.model.variables.lower
    bu = fill(Inf, A.n)
    bu .= src.model.variables.upper
    if isempty(G)
        error("The constraint matrix G is empty, PDHG solver cannot be applied.")
    end
    if !haskey(options, :verbose)
        options[:verbose] = 0
    end
    if !haskey(options, :use_scaling)
        options[:use_scaling] = true
        if options[:verbose] > 0
            println("use_scaling is not set, using default value: true")
        end
    end
    if !haskey(options, :rescaling_method)
        options[:rescaling_method] = :ruiz_pock_chambolle
        if options[:verbose] > 0
            println("rescaling_method is not set, using default value: :ruiz_pock_chambolle. Note that if :use_scaling is false, this option is ignored.")
        end
    end
    if !haskey(options, :use_adaptive_restart)
        options[:use_adaptive_restart] = true
        if options[:verbose] > 0
            println("use_adaptive_restart is not set, using default value: true")
        end
    end
    if !haskey(options, :use_adaptive_step_size_weight)
        options[:use_adaptive_step_size_weight] = true
        if options[:verbose] > 0
            println("use_adaptive_step_size_weight is not set, using default value: true")
        end
    end
    if !haskey(options, :use_resolving)
        options[:use_resolving] = false
        if options[:verbose] > 0
            println("use_resolving is not set, using default value: true")
        end
    end
    if !haskey(options, :use_accelerated)
        options[:use_accelerated] = true
        if options[:verbose] > 0
            println("use_accelerated is not set, using default value: true")
        end
    end
    if !haskey(options, :rel_tol)
        options[:rel_tol] = 1e-6
    end
    if !haskey(options, :abs_tol)
        options[:abs_tol] = 1e-6
    end
    ##### SOLVER CODE
  ########
        dest.sol = MOISolution(
            primal = sol_res.x.primal_sol.x,
            dual = sol_res.y.dual_sol.y,
            slack = sol_res.y.slack.primal_sol.x,
            exit_code = sol_res.info.exit_code,
            exit_status = sol_res.info.exit_status,
            objective_value = (max_sense ? -1 : 1) * sol_res.info.pObj,
            dual_objective_value = (max_sense ? -1 : 1) * sol_res.info.dObj,
            solve_time_sec = sol_res.info.time,
            iterations = sol_res.info.iter,
            objective_constant = objective_constant,
        )
    end
    return index_map, false
end

function MOI.optimize!(dest::Optimizer, src::MOI.ModelLike)
    cache = MOI.Utilities.UniversalFallback(OptimizerCache{Integer}())
    index_map = MOI.copy_to(cache, src)
    MOI.optimize!(dest, cache)
    return index_map, false
end

function MOI.get(optimizer::Optimizer, ::MOI.SolveTimeSec)
    return optimizer.sol.solve_time_sec
end

function MOI.get(optimizer::Optimizer, ::MOI.RawStatusString)
    return optimizer.sol.raw_status
end

"""
    PDHGIterations()

The number of PDHG iterations completed during the solve.
"""
struct PDHGIterations <: MOI.AbstractModelAttribute end

MOI.is_set_by_optimize(::PDHGIterations) = true

function MOI.get(optimizer::Optimizer, ::PDHGIterations)
    return optimizer.sol.iterations
end

"""
exit_status:
    :optimal 0
    :max_iter 1
    :primal_infeasible_low_acc 2
    :primal_infeasible_high_acc 3
    :dual_infeasible_low_acc 4
    :dual_infeasible_high_acc 5
    :time_limit 6   
    :continue 7
    :numerical_error 8
"""
function MOI.get(optimizer::Optimizer, ::MOI.TerminationStatus)
    s = optimizer.sol.exit_code
    @assert 0 <= s <= 6
    if s == 0
        return MOI.OPTIMAL
    elseif s == 1
        return MOI.ITERATION_LIMIT
    elseif s == 2
        return MOI.INFEASIBLE
    elseif s == 3
        return MOI.DUAL_INFEASIBLE
    elseif s == 4
        return MOI.INFEASIBLE
    elseif s == 5
        return MOI.DUAL_INFEASIBLE
    elseif s == 6
        return MOI.TIME_LIMIT
    elseif s == 7
        return MOI.UNFINISHED
    elseif s == 8
        return MOI.NUMERICAL_ERROR
    else
        throw(MOI.UnsupportedAttribute(MOI.TerminationStatus()))
        if occursin("reached time_limit_secs", optimizer.sol.raw_status)
            return MOI.TIME_LIMIT
        elseif occursin("reached max_iters", optimizer.sol.raw_status)
            return MOI.ITERATION_LIMIT
        else
            return MOI.ALMOST_OPTIMAL
        end
    # elseif s == -5
    #     return MOI.INTERRUPTED
    # elseif s == -4
    #     return MOI.INVALID_MODEL
    # elseif s == -3
    #     return MOI.SLOW_PROGRESS
    # elseif s == -2
    #     return MOI.INFEASIBLE
    # elseif s == -1
    #     return MOI.DUAL_INFEASIBLE
    # elseif s == 1
    #     return MOI.OPTIMAL
    # else
    #     @assert s == 0
    #     return MOI.OPTIMIZE_NOT_CALLED
    end
end

function MOI.get(optimizer::Optimizer, attr::MOI.ObjectiveValue)
    MOI.check_result_index_bounds(optimizer, attr)
    value = optimizer.sol.objective_value
    if !MOI.Utilities.is_ray(MOI.get(optimizer, MOI.PrimalStatus()))
        value += optimizer.sol.objective_constant
    end
    return value
end

function MOI.get(optimizer::Optimizer, attr::MOI.DualObjectiveValue)
    MOI.check_result_index_bounds(optimizer, attr)
    value = optimizer.sol.dual_objective_value
    if !MOI.Utilities.is_ray(MOI.get(optimizer, MOI.DualStatus()))
        value += optimizer.sol.objective_constant
    end
    return value
end

function MOI.get(optimizer::Optimizer, attr::MOI.PrimalStatus)
    if attr.result_index > MOI.get(optimizer, MOI.ResultCount())
        return MOI.NO_SOLUTION
    elseif optimizer.sol.exit_code in (0)
        return MOI.FEASIBLE_POINT
    elseif optimizer.sol.exit_code in (2, 3, 4, 5, 6)
        return MOI.INFEASIBILITY_CERTIFICATE
    elseif optimizer.sol.exit_code == (1, 7)
        return MOI.NO_SOLUTION
    end
    return MOI.INFEASIBLE_POINT
end

function MOI.get(
    optimizer::Optimizer,
    attr::MOI.VariablePrimal,
    vi::MOI.VariableIndex,
)
    MOI.check_result_index_bounds(optimizer, attr)
    return optimizer.sol.primal[vi.value]
end

function MOI.get(
    optimizer::Optimizer,
    attr::MOI.ConstraintPrimal,
    ci::MOI.ConstraintIndex{F,S},
) where {F,S}
    MOI.check_result_index_bounds(optimizer, attr)
    return optimizer.sol.slack[MOI.Utilities.rows(optimizer.cones, ci)]
end

function MOI.get(optimizer::Optimizer, attr::MOI.DualStatus)
    if attr.result_index > MOI.get(optimizer, MOI.ResultCount())
        return MOI.NO_SOLUTION
    elseif optimizer.sol.exit_code in (0)
        return MOI.FEASIBLE_POINT
    elseif optimizer.sol.exit_code in (2, 3, 4, 5, 6)
        return MOI.INFEASIBILITY_CERTIFICATE
    elseif optimizer.sol.exit_code == (1, 7)
        return MOI.NO_SOLUTION
    end
    return MOI.INFEASIBLE_POINT
end

function MOI.get(
    optimizer::Optimizer,
    attr::MOI.ConstraintDual,
    ci::MOI.ConstraintIndex{F,S},
) where {F,S}
    MOI.check_result_index_bounds(optimizer, attr)
    return optimizer.sol.dual[MOI.Utilities.rows(optimizer.cones, ci)]
end

MOI.get(optimizer::Optimizer, ::MOI.ResultCount) = 1

I understand that the original SCS code does not include VectorOfVariables; it incorporates all cone constraints for each variable directly into the constraints. However, this approach may not align with our modeling objectives. Therefore, I would like to know how to adapt this MOI_wrapper to achieve the following:

If all variables have box constraints, model all cone constraints within the constraint matrix.
If a single variable appears in multiple cone constraints, as in the previous example, retain one active cone constraint for the variable and place the remaining cone constraints within the constraint matrix.

Many thanks to the Julia community, and I would greatly appreciate insights from experienced members on this issue.

odow · October 28, 2024, 8:18am

Hi @zhenweilin, welcome to the forum

I’ve moved your question to the “Optimization (Mathematical)” section. I didn’t see it when you first posted because it was in the general category.

You need to use MOI.add_constrained_variable and MOI.supports_add_constrained_variable, but I don’t know how well these are supported with product of sets. It’s late here, so I will take a look tomorrow.

@blegat may have some ideas.

odow · October 28, 2024, 8:58pm

One option would just be:

Greedily add variables to bound constraints, and then if a VectorOfVariables appears in a cone and it already has a bound constraint, then add the cone as an affine constraint instead.

This is what we do for the CBF writer:

github.com

jump-dev/MathOptInterface.jl/blob/master/src/FileFormats/CBF/write.jl

# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

"""
    _CBFDataStructure

A helpful datastructure for converting between MOI and CBF.
"""
mutable struct _CBFDataStructure
    num_rows::Int
    power_cones::Int
    dual_power_cones::Int
    psd_side_dims::Vector{Int}
    cones::Vector{Tuple{String,Int}}
    acoord::Vector{Tuple{Int,Int,Float64}}
    bcoord::Vector{Tuple{Int,Float64}}
    hcoord::Vector{Tuple{Int,Int,Int,Int,Float64}}

This file has been truncated. show original

zhenweilin · October 29, 2024, 1:50am

Thank you for your answer! For testing the operation that write model into cbf file can help me build the right model. I try to write the model into .cbf file and read it again. I use the following code to test:

using Pkg
Pkg.activate("test_cpu_env")
using Arpack, LinearAlgebra
using JuMP, MosekTools
using Profile
using Random, SparseArrays
using SCS, COPT

module my_moi_wrapper
import MathOptInterface as MOI
using Random, SparseArrays, LinearAlgebra, KrylovKit
using Printf
using CSV
using DataFrames
using FilePathsBase
using Match
using DataStructures
using Base.Threads
using JuMP
using MosekTools
using Polynomials
using Statistics
"""
    struct ScaledPSDCone <: MOI.AbstractVectorSet
        side_dimension::Int
    end

Similar to `MOI.Scaled{MOI.PositiveSemidefiniteConeTriangle}` but it the
vectorization is the lower triangular part column-wise (or the upper triangular
part row-wise).
"""
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,F} <: MOI.Bridges.Constraint.SetMapBridge{
    T,
    ScaledPSDCone,
    MOI.Scaled{MOI.PositiveSemidefiniteConeTriangle},
    F,
    F,
}
    constraint::MOI.ConstraintIndex{F,ScaledPSDCone}
end

function MOI.Bridges.Constraint.concrete_bridge_type(
    ::Type{ScaledPSDConeBridge{T}},
    ::Type{F},
    ::Type{MOI.Scaled{MOI.PositiveSemidefiniteConeTriangle}},
) where {T,F<:MOI.AbstractVectorFunction}
    return ScaledPSDConeBridge{T,F}
end

function MOI.Bridges.map_set(
    ::Type{<:ScaledPSDConeBridge},
    set::MOI.Scaled{MOI.PositiveSemidefiniteConeTriangle},
)
    return ScaledPSDCone(MOI.side_dimension(set))
end

function MOI.Bridges.inverse_map_set(
    ::Type{<:ScaledPSDConeBridge},
    set::ScaledPSDCone,
)
    return MOI.Scaled(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_to_pdhgclp::Bool)
    scalars = MOI.Utilities.eachscalar(func)
    d = length(scalars)
    upper_to_lower, lower_to_upper = _upper_to_lower_triangular_permutation(d)
    if moi_to_pdhgclp
        return scalars[lower_to_upper]
    else
        return scalars[upper_to_lower]
    end
end

# Map ConstraintFunction from MOI -> PDHGCLP
function MOI.Bridges.map_function(::Type{<:ScaledPSDConeBridge}, f)
    return _transform_function(f, true)
end

# Used to map the ConstraintPrimal from PDHGCLP -> MOI
function MOI.Bridges.inverse_map_function(::Type{<:ScaledPSDConeBridge}, f)
    return _transform_function(f, false)
end

# Used to map the ConstraintDual from PDHGCLP -> MOI
function MOI.Bridges.adjoint_map_function(::Type{<:ScaledPSDConeBridge}, f)
    return _transform_function(f, false)
end

# Used to set ConstraintDualStart
function MOI.Bridges.inverse_adjoint_map_function(
    ::Type{<:ScaledPSDConeBridge},
    f,
)
    return _transform_function(f, true)
end

MOI.Utilities.@product_of_sets(
    _Cones,
    MOI.Zeros,
    MOI.Nonnegatives,
    MOI.SecondOrderCone,
    MOI.RotatedSecondOrderCone,
    MOI.ExponentialCone,
    MOI.DualExponentialCone,
    # MOI.PowerCone{T},
    # MOI.DualPowerCone{T},
)

struct _SetConstants{T}
    b::Vector{T}
    power_coefficients::Dict{Int,T}
    _SetConstants{T}() where {T} = new{T}(T[], Dict{Int,T}())
end

function Base.empty!(x::_SetConstants)
    empty!(x.b)
    empty!(x.power_coefficients)
    return x
end

Base.resize!(x::_SetConstants, n) = resize!(x.b, n)

function MOI.Utilities.load_constants(x::_SetConstants, offset, f)
    MOI.Utilities.load_constants(x.b, offset, f)
    return
end

function MOI.Utilities.load_constants(
    x::_SetConstants{T},
    offset,
    set::Union{MOI.PowerCone{T},MOI.DualPowerCone{T}},
) where {T}
    x.power_coefficients[offset+1] = set.exponent
    return
end

function MOI.Utilities.function_constants(x::_SetConstants, rows)
    return MOI.Utilities.function_constants(x.b, rows)
end

function MOI.Utilities.set_from_constants(x::_SetConstants, S, rows)
    return MOI.Utilities.set_from_constants(x.b, S, rows)
end

function MOI.Utilities.set_from_constants(
    x::_SetConstants{T},
    ::Type{S},
    rows,
) where {T,S<:Union{MOI.PowerCone{T},MOI.DualPowerCone{T}}}
    @assert length(rows) == 3
    return S(x.power_coefficients[first(rows)])
end

const OptimizerCache{T} = MOI.Utilities.GenericModel{
    Float64,
    MOI.Utilities.ObjectiveContainer{Float64},
    MOI.Utilities.VariablesContainer{Float64},
    MOI.Utilities.MatrixOfConstraints{
        Float64,
        MOI.Utilities.MutableSparseMatrixCSC{
            Float64,
            T,
            MOI.Utilities.OneBasedIndexing,
        },
        _SetConstants{Float64},
        _Cones{Float64},
    },
}

function _to_sparse(
    ::Type{T},
    A::MOI.Utilities.MutableSparseMatrixCSC{
        Float64,
        T,
        MOI.Utilities.OneBasedIndexing,
    },
) where {T}
    return -A.nzval, A.rowval, A.colptr
end

mutable struct MOISolution
    primal::Vector{Float64}
    dual::Vector{Float64}
    slack::Vector{Float64}
    exit_code::Int
    exit_status::Symbol
    objective_value::Float64
    dual_objective_value::Float64
    solve_time_sec::Float64
    iterations::Int
    objective_constant::Float64
    function MOISolution(;primal = Float64[], dual = Float64[], slack = Float64[], exit_code = 0, exit_status = :unknown, objective_value = NaN, dual_objective_value = NaN, solve_time_sec = NaN, iterations = 0, objective_constant = 0.0)
        new(primal, dual, slack, exit_code, exit_status, objective_value, dual_objective_value, solve_time_sec, iterations, objective_constant)
    end
end

"""
    Optimizer()

Create a new PDHG-CLP optimizer.
"""
mutable struct Optimizer <: MOI.AbstractOptimizer
    cones::Union{Nothing,_Cones{Float64}}
    sol::MOISolution
    silent::Bool
    options::Dict{Symbol,Any}

    Optimizer() = new(nothing, MOISolution(), false, Dict{Symbol,Any}())
end

function MOI.get(::Optimizer, ::MOI.Bridges.ListOfNonstandardBridges)
    return [ScaledPSDConeBridge{Float64}]
end

MOI.get(::Optimizer, ::MOI.SolverName) = "PDHG-CLP"

MOI.is_empty(optimizer::Optimizer) = optimizer.cones === nothing

function MOI.empty!(optimizer::Optimizer)
    optimizer.cones = nothing
    optimizer.sol = MOISolution()
    return
end

###
### MOI.RawOptimizerAttribute
###

MOI.supports(::Optimizer, ::MOI.RawOptimizerAttribute) = true

function MOI.set(optimizer::Optimizer, param::MOI.RawOptimizerAttribute, value)
    return optimizer.options[Symbol(param.name)] = value
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)
    optimizer.silent = value
    return
end

MOI.get(optimizer::Optimizer, ::MOI.Silent) = optimizer.silent

###
### MOI.TimeLimitSec
###

MOI.supports(::Optimizer, ::MOI.TimeLimitSec) = true

function MOI.set(optimizer::Optimizer, ::MOI.TimeLimitSec, time_limit::Real)
    optimizer.options[:time_limit_secs] = convert(Float64, time_limit)
    return
end

function MOI.set(optimizer::Optimizer, ::MOI.TimeLimitSec, ::Nothing)
    delete!(optimizer.options, :time_limit_secs)
    return
end

function MOI.get(optimizer::Optimizer, ::MOI.TimeLimitSec)
    value = get(optimizer.options, :time_limit_secs, nothing)
    return value::Union{Float64,Nothing}
end

###
### MOI.AbstractModelAttribute
###

function MOI.supports(
    ::Optimizer,
    ::Union{
        MOI.ObjectiveSense,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}},
        # MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}},
    },
)
    return true
end

###
### MOI.AbstractVariableAttribute
###

function MOI.supports(
    ::Optimizer,
    ::MOI.VariablePrimalStart,
    ::Type{MOI.VariableIndex},
)
    return true
end

###
### MOI.AbstractConstraintAttribute
###

function MOI.supports(
    ::Optimizer,
    ::Union{MOI.ConstraintPrimalStart,MOI.ConstraintDualStart},
    ::Type{<:MOI.ConstraintIndex},
)
    return true
end

function MOI.supports_constraint(
    ::Optimizer,
    ::Type{MOI.VectorAffineFunction{Float64}},
    ::Type{
        <:Union{
            MOI.Zeros,
            MOI.Nonnegatives,
            MOI.SecondOrderCone,
            MOI.RotatedSecondOrderCone,
            MOI.ExponentialCone,
            MOI.DualExponentialCone,
            # MOI.PowerCone{Float64},
            # MOI.DualPowerCone{Float64},
        },
    },
)
    return true
end

# customized
function MOI.supports_constraint(
    ::Optimizer,
    ::Type{MOI.VariableIndex},
    ::Type{
        <:Union{
            MOI.EqualTo{Float64},
            MOI.LessThan{Float64},
            MOI.GreaterThan{Float64},
            MOI.Interval{Float64},
        }
    }
)
    return true
end

# # customized
# function MOI.supports_constraint(
#     ::Optimizer,
#     ::Type{MOI.VectorOfVariables},
#     ::Type{
#         <:Union{
#             MOI.SecondOrderCone,
#             MOI.RotatedSecondOrderCone,
#             # may be not supported for two exponential cones
#             MOI.ExponentialCone,
#             MOI.DualExponentialCone,
#         }
#     }
# )
#     return true
# end


function _map_sets(f, ::Type{T}, sets, ::Type{S}) where {T,S}
    F = MOI.VectorAffineFunction{Float64}
    cis = MOI.get(sets, MOI.ListOfConstraintIndices{F,S}())
    return T[f(MOI.get(sets, MOI.ConstraintSet(), ci)) for ci in cis]
end

function MOI.optimize!(
    dest::Optimizer,
    src::MOI.Utilities.UniversalFallback{OptimizerCache{T}},
) where {T}
    # The real stuff starts here.
    MOI.empty!(dest)
    index_map = MOI.Utilities.identity_index_map(src)
    Ab = src.model.constraints
    A = Ab.coefficients

    for (F, S) in keys(src.constraints) # src.constraints is very important
        throw(MOI.UnsupportedConstraint{F,S}())
        # if F == MOI.VectorOfVariables
        #     if S == MOI.SecondOrderCone
        #         println("SecondOrderCone")
        #     elseif S == MOI.RotatedSecondOrderCone
        #         println("RotatedSecondOrderCone")
        #     elseif S == MOI.ExponentialCone
        #         println("ExponentialCone")
        #     elseif S == MOI.DualExponentialCone
        #         println("DualExponentialCone")
        #     else
        #         throw(MOI.UnsupportedConstraint{F,S}())
        #     end
        # else
        #     throw(MOI.UnsupportedConstraint{F,S}())
        # end
    end
    # key: src.model.variables is also very important
    model_attributes = MOI.get(src, MOI.ListOfModelAttributesSet())
    max_sense = false
    obj_attr = nothing
    for attr in model_attributes
        if attr == MOI.ObjectiveSense()
            max_sense = MOI.get(src, attr) == MOI.MAX_SENSE
        elseif attr == MOI.Name()
            continue  # This can be skipped without consequence
        elseif attr isa MOI.ObjectiveFunction
            obj_attr = attr
        else
            throw(MOI.UnsupportedAttribute(attr))
        end
    end
    objective_constant, c = 0.0, zeros(A.n)
    if obj_attr == MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}()
        obj = MOI.get(src, obj_attr)
        objective_constant = MOI.constant(obj)
        for term in obj.terms
            c[term.variable.value] += (max_sense ? -1 : 1) * term.coefficient
        end
    elseif obj_attr !== nothing
        throw(MOI.UnsupportedAttribute(obj_attr))
    end

    # `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
    if A.m != length(dest.sol.dual)
        resize!(dest.sol.dual, A.m)
        fill!(dest.sol.dual, 0.0)
    end
    if A.m != length(dest.sol.slack)
        resize!(dest.sol.slack, A.m)
        fill!(dest.sol.slack, 0.0)
    end
    # Set starting values and throw error for other variable attributes
    has_warm_start = false
    vis_src = MOI.get(src, MOI.ListOfVariableIndices())
    for attr in MOI.get(src, MOI.ListOfVariableAttributesSet())
        if attr == MOI.VariableName()
            # Skip
        elseif attr == MOI.VariablePrimalStart()
            has_warm_start = true
            for (i, x) in enumerate(vis_src)
                dest.sol.primal[i] = something(MOI.get(src, attr, x), 0.0)
            end
        else
            throw(MOI.UnsupportedAttribute(attr))
        end
    end
    # 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}())
        for attr in MOI.get(src, MOI.ListOfConstraintAttributesSet{F,S}())
            if attr == MOI.ConstraintName()
                # Skip
            elseif attr == MOI.ConstraintPrimalStart()
                has_warm_start = true
                for ci in cis_src
                    start = MOI.get(src, attr, ci)
                    if start !== nothing
                        rows = MOI.Utilities.rows(Ab, ci)
                        dest.sol.slack[rows] .= start
                    end
                end
            elseif attr == MOI.ConstraintDualStart()
                has_warm_start = true
                for ci in cis_src
                    start = MOI.get(src, attr, ci)
                    if start !== nothing
                        rows = MOI.Utilities.rows(Ab, ci)
                        dest.sol.dual[rows] .= start
                    end
                end
            else
                throw(MOI.UnsupportedAttribute(attr))
            end
        end
    end
    options = copy(dest.options)
    if dest.silent
        options[:verbose] = 0
    end
    dest.cones = deepcopy(Ab.sets)
    mGzero = sum(_map_sets(MOI.dimension, T, Ab, MOI.Zeros))
    mGnonnegative = sum(_map_sets(MOI.dimension, T, Ab, MOI.Nonnegatives))
    G = SparseMatrixCSC{Float64, Int64}(A.m, A.n, Int64.(A.colptr), Int64.(A.rowval), A.nzval)
    bl = fill(-Inf, A.n)
    bl .= src.model.variables.lower
    bu = fill(Inf, A.n)
    bu .= src.model.variables.upper
    if isempty(G)
        error("The constraint matrix G is empty, PDHG solver cannot be applied.")
    end
    if !haskey(options, :verbose)
        options[:verbose] = 0
    end
    if !haskey(options, :use_scaling)
        options[:use_scaling] = true
        if options[:verbose] > 0
            println("use_scaling is not set, using default value: true")
        end
    end
    if !haskey(options, :rescaling_method)
        options[:rescaling_method] = :ruiz_pock_chambolle
        if options[:verbose] > 0
            println("rescaling_method is not set, using default value: :ruiz_pock_chambolle. Note that if :use_scaling is false, this option is ignored.")
        end
    end
    if !haskey(options, :use_adaptive_restart)
        options[:use_adaptive_restart] = true
        if options[:verbose] > 0
            println("use_adaptive_restart is not set, using default value: true")
        end
    end
    if !haskey(options, :use_adaptive_step_size_weight)
        options[:use_adaptive_step_size_weight] = true
        if options[:verbose] > 0
            println("use_adaptive_step_size_weight is not set, using default value: true")
        end
    end
    if !haskey(options, :use_resolving)
        options[:use_resolving] = false
        if options[:verbose] > 0
            println("use_resolving is not set, using default value: true")
        end
    end
    if !haskey(options, :print_freq)
        options[:print_freq] = 2000
        if options[:verbose] > 0
            println("print_freq is not set, using default value: 2000")
        end
    end
    if !haskey(options, :use_accelerated)
        options[:use_accelerated] = true
        if options[:verbose] > 0
            println("use_accelerated is not set, using default value: true")
        end
    end
    if !haskey(options, :rel_tol)
        options[:rel_tol] = 1e-6
    end
    if !haskey(options, :abs_tol)
        options[:abs_tol] = 1e-6
    end
    if options[:use_scaling]
        return true
    else
        return false
    end
    return index_map, false
end

function MOI.optimize!(dest::Optimizer, src::MOI.ModelLike)
    cache = MOI.Utilities.UniversalFallback(OptimizerCache{Integer}())
    index_map = MOI.copy_to(cache, src)
    MOI.optimize!(dest, cache)
    return index_map, false
end

function MOI.get(optimizer::Optimizer, ::MOI.SolveTimeSec)
    return optimizer.sol.solve_time_sec
end

function MOI.get(optimizer::Optimizer, ::MOI.RawStatusString)
    return optimizer.sol.raw_status
end

"""
    PDHGIterations()

The number of PDHG iterations completed during the solve.
"""
struct PDHGIterations <: MOI.AbstractModelAttribute end

MOI.is_set_by_optimize(::PDHGIterations) = true

function MOI.get(optimizer::Optimizer, ::PDHGIterations)
    return optimizer.sol.iterations
end

"""
exit_status:
    :optimal 0
    :max_iter 1
    :primal_infeasible_low_acc 2
    :primal_infeasible_high_acc 3
    :dual_infeasible_low_acc 4
    :dual_infeasible_high_acc 5
    :time_limit 6   
    :continue 7
    :numerical_error 8
"""
function MOI.get(optimizer::Optimizer, ::MOI.TerminationStatus)
    s = optimizer.sol.exit_code
    @assert 0 <= s <= 8
    if s == 0
        return MOI.OPTIMAL
    elseif s == 1
        return MOI.ITERATION_LIMIT
    elseif s == 2
        return MOI.INFEASIBLE
    elseif s == 3
        return MOI.DUAL_INFEASIBLE
    elseif s == 4
        return MOI.INFEASIBLE
    elseif s == 5
        return MOI.DUAL_INFEASIBLE
    elseif s == 6
        return MOI.TIME_LIMIT
    elseif s == 7
        return MOI.UNFINISHED
    elseif s == 8
        return MOI.NUMERICAL_ERROR
    else
        throw(MOI.UnsupportedAttribute(MOI.TerminationStatus()))
        if occursin("reached time_limit_secs", optimizer.sol.raw_status)
            return MOI.TIME_LIMIT
        elseif occursin("reached max_iters", optimizer.sol.raw_status)
            return MOI.ITERATION_LIMIT
        else
            return MOI.ALMOST_OPTIMAL
        end
    # elseif s == -5
    #     return MOI.INTERRUPTED
    # elseif s == -4
    #     return MOI.INVALID_MODEL
    # elseif s == -3
    #     return MOI.SLOW_PROGRESS
    # elseif s == -2
    #     return MOI.INFEASIBLE
    # elseif s == -1
    #     return MOI.DUAL_INFEASIBLE
    # elseif s == 1
    #     return MOI.OPTIMAL
    # else
    #     @assert s == 0
    #     return MOI.OPTIMIZE_NOT_CALLED
    end
end

function MOI.get(optimizer::Optimizer, attr::MOI.ObjectiveValue)
    MOI.check_result_index_bounds(optimizer, attr)
    value = optimizer.sol.objective_value
    if !MOI.Utilities.is_ray(MOI.get(optimizer, MOI.PrimalStatus()))
        value += optimizer.sol.objective_constant
    end
    return value
end

function MOI.get(optimizer::Optimizer, attr::MOI.DualObjectiveValue)
    MOI.check_result_index_bounds(optimizer, attr)
    value = optimizer.sol.dual_objective_value
    if !MOI.Utilities.is_ray(MOI.get(optimizer, MOI.DualStatus()))
        value += optimizer.sol.objective_constant
    end
    return value
end

function MOI.get(optimizer::Optimizer, attr::MOI.PrimalStatus)
    if attr.result_index > MOI.get(optimizer, MOI.ResultCount())
        return MOI.NO_SOLUTION
    elseif optimizer.sol.exit_code in (0)
        return MOI.FEASIBLE_POINT
    elseif optimizer.sol.exit_code in (2, 3, 4, 5, 6)
        return MOI.INFEASIBILITY_CERTIFICATE
    elseif optimizer.sol.exit_code == (1, 7)
        return MOI.NO_SOLUTION
    end
    return MOI.INFEASIBLE_POINT
end

function MOI.get(
    optimizer::Optimizer,
    attr::MOI.VariablePrimal,
    vi::MOI.VariableIndex,
)
    MOI.check_result_index_bounds(optimizer, attr)
    return optimizer.sol.primal[vi.value]
end

function MOI.get(
    optimizer::Optimizer,
    attr::MOI.ConstraintPrimal,
    ci::MOI.ConstraintIndex{F,S},
) where {F,S}
    MOI.check_result_index_bounds(optimizer, attr)
    return optimizer.sol.slack[MOI.Utilities.rows(optimizer.cones, ci)]
end

function MOI.get(optimizer::Optimizer, attr::MOI.DualStatus)
    if attr.result_index > MOI.get(optimizer, MOI.ResultCount())
        return MOI.NO_SOLUTION
    elseif optimizer.sol.exit_code in (0)
        return MOI.FEASIBLE_POINT
    elseif optimizer.sol.exit_code in (2, 3, 4, 5, 6)
        return MOI.INFEASIBILITY_CERTIFICATE
    elseif optimizer.sol.exit_code == (1, 7)
        return MOI.NO_SOLUTION
    end
    return MOI.INFEASIBLE_POINT
end

function MOI.get(
    optimizer::Optimizer,
    attr::MOI.ConstraintDual,
    ci::MOI.ConstraintIndex{F,S},
) where {F,S}
    MOI.check_result_index_bounds(optimizer, attr)
    return optimizer.sol.dual[MOI.Utilities.rows(optimizer.cones, ci)]
end

MOI.get(optimizer::Optimizer, ::MOI.ResultCount) = 1
end


model = JuMP.Model(my_moi_wrapper.Optimizer);
n = 10
k_box = 2
k_soc1 = 4
k_soc2 = 4
c = rand(n)
c = abs.(c)
@variable(model, y[1:n]>=1);
@objective(model, Min, c' * y);
x_soc1 = @constraint(model, [y[k_box + 1]; y[k_box + 2: k_box + k_soc1]] in SecondOrderCone());
x_soc2 = @constraint(model, [y[k_box]; y[k_box + 2: k_box + k_soc1 + 1]]in SecondOrderCone());
x_soc3 = @constraint(model, [2 * y[k_box + 1]; y[k_box + 2: k_box + k_soc1 + 1]]in SecondOrderCone());

MOI.write_to_file(backend(model), "debug.cbf")
println("model: ", model)
optimize!(model)


moi_model = MOI.Utilities.Model{Float64}()
MOI.read_from_file(moi_model, "debug.cbf")
model = JuMP.Model(my_moi_wrapper.Optimizer);
MOI.copy_to(model, moi_model)
println("model: ", model)
optimize!(model)

And the result is

model: Min 0.10493433193272206 y[1] + 0.32840368388224184 y[2] + 0.33018114859384595 y[3] + 0.7065883184332403 y[4] + 0.7824328700290073 y[5] + 0.4232458504565929 y[6] + 0.6553639898868799 y[7] + 0.8133423400201902 y[8] + 0.33104566306610916 y[9] + 0.1490942880929491 y[10]
Subject to
 [y[3], y[4], y[5], y[6]] ∈ MathOptInterface.SecondOrderCone(4)
 [y[2], y[4], y[5], y[6], y[7]] ∈ MathOptInterface.SecondOrderCone(5)
 [2 y[3], y[4], y[5], y[6], y[7]] ∈ MathOptInterface.SecondOrderCone(5)
 y[1] ≥ 1
 y[2] ≥ 1
 y[3] ≥ 1
 y[4] ≥ 1
 y[5] ≥ 1
 y[6] ≥ 1
 y[7] ≥ 1
 y[8] ≥ 1
 y[9] ≥ 1
 y[10] ≥ 1

model: Min 0.10493433193272206 _[1] + 0.32840368388224184 _[2] + 0.33018114859384595 _[3] + 0.7065883184332403 _[4] + 0.7824328700290073 _[5] + 0.4232458504565929 _[6] + 0.6553639898868799 _[7] + 0.8133423400201902 _[8] + 0.33104566306610916 _[9] + 0.1490942880929491 _[10]
Subject to
 [_[3], _[4], _[5], _[6]] ∈ MathOptInterface.SecondOrderCone(4)
 [_[2], _[4], _[5], _[6], _[7]] ∈ MathOptInterface.SecondOrderCone(5)
 [2 _[3], _[4], _[5], _[6], _[7]] ∈ MathOptInterface.SecondOrderCone(5)

It seems that the >=1 constraint is disappear? This is a bug? Or what did I miss?

odow · October 30, 2024, 1:18am

This is a slight mis-match because you used MOI.write_to_file(backend(model), "debug.cbf") instead of JuMP.write_to_file(model, "debug.cbf").

I need to look into the details of why:

julia> using JuMP

julia> model = JuMP.Model();

julia> @variable(model, x >= 1)
x

julia> print(model)
Feasibility
Subject to
 x ≥ 1

julia> MOI.write_to_file(backend(model), "debug.cbf")

shell> cat debug.cbf
VER
3

OBJSENSE
MIN

VAR
1 1
F 1


julia> write_to_file(model, "debug.cbf")

shell> cat debug.cbf
VER
3

OBJSENSE
MIN

VAR
1 1
L+ 1

Edit: [FileFormats.CBF] bug copying unsupported variable bounds to CBF · Issue #2570 · jump-dev/MathOptInterface.jl · GitHub

zhenweilin · October 30, 2024, 2:17am

When I use code write_to_file(model, "debug.cbf"), the model read from debug.cbf file is

model: Min 0.12681356180282743 _[1] + 0.835694412153161 _[2] + 0.6947070995655321 _[3] + 0.5122754754644325 _[4] + 0.418850346186468 _[5] + 0.20975067953051874 _[6] + 0.5074994417374445 _[7] + 0.1902321155945571 _[8] + 0.07055879275182997 _[9] + 0.6747303010006362 _[10] + 2.405528625040456
Subject to
 [_[1]] ∈ MathOptInterface.Nonnegatives(1)
 [_[2]] ∈ MathOptInterface.Nonnegatives(1)
 [_[7]] ∈ MathOptInterface.Nonnegatives(1)
 [_[8]] ∈ MathOptInterface.Nonnegatives(1)
 [_[9]] ∈ MathOptInterface.Nonnegatives(1)
 [_[10]] ∈ MathOptInterface.Nonnegatives(1)
 [_[3], _[4], _[5], _[6]] ∈ MathOptInterface.SecondOrderCone(4)
 [_[3] - 1] ∈ MathOptInterface.Nonnegatives(1)
 [_[4] - 1] ∈ MathOptInterface.Nonnegatives(1)
 [_[5] - 1] ∈ MathOptInterface.Nonnegatives(1)
 [_[6] - 1] ∈ MathOptInterface.Nonnegatives(1)
 [_[2] + 1, _[4], _[5], _[6], _[7] + 1] ∈ MathOptInterface.SecondOrderCone(5)
 [2 _[3], _[4], _[5], _[6], _[7] + 1] ∈ MathOptInterface.SecondOrderCone(5)

It appears that the model currently lacks a box constraint, which does not align with the intended requirements. The write to .cbf file seems to build a model with free variable and all constraint will be put into constraint matrix.

odow · October 30, 2024, 2:24am

What’s the purpose of writing to CBF?

The formulation doesn’t line up exactly because we have bridged it. Some of the variables y >= 1 in your problem have been rewritten to z >= 0, and y = 1 + z. That’s why you see, for example, _[2] + 1.

zhenweilin · October 30, 2024, 2:41am

I realize now that I previously misinterpreted your explanation. I initially thought that the code generating the .cbf file would directly yield the model I wanted, and that I could reference the .cbf file code. I sincerely apologize for having overlooked your earlier suggestion to use add_constrained_variable and will proceed to implement it.

odow · October 30, 2024, 3:11am

No worries. I mentioned CBF only because, if you look through the source code, you may find ideas in the CBF writer.

It also has the condition that a variable can belong to at most one cone.