test_csdp.jl

#=
(soi) pkg> st
    Status `~/soi/Project.toml`
  [0a46da34] CSDP v0.5.0 #bl/moiv0.9 (https://github.com/JuliaOpt/CSDP.jl.git)
  [b8f27783] MathOptInterface v0.9.0 #master (https://github.com/JuliaOpt/MathOptInterface.jl.git)
  [f0680fed] SemidefiniteOptInterface v0.6.0 #bl/moiv0.9 (https://github.com/JuliaOpt/SemidefiniteOptInterface.jl.git)
=#

using LinearAlgebra
using MathOptInterface
const MOI = MathOptInterface
const MOIT = MOI.Test
const MOIB = MOI.Bridges
const MOIU = MOI.Utilities

#=
    Helper functions
=#
LinearAlgebra.symmetric_type(::Type{MathOptInterface.VariableIndex}) = MathOptInterface.VariableIndex
LinearAlgebra.symmetric(v::MathOptInterface.VariableIndex, ::Symbol) = v
LinearAlgebra.transpose(v::MathOptInterface.VariableIndex) = v

sympackedlen(n) = div(n*(n+1), 2)
sympackeddim(n) = div(isqrt(1+8n) - 1, 2)
function ivech!(out::AbstractMatrix{T}, v::AbstractVector{T}) where T
    n = sympackeddim(length(v))
    n1, n2 = size(out)
    @assert n == n1 == n2
    c = 0
    for j in 1:n, i in 1:j
        c += 1
        out[i,j] = v[c]
    end
    return out
end
function ivech(v::AbstractVector{T}) where T
    n = sympackeddim(length(v))
    out = zeros(n, n)
    ivech!(out, v)
    
    return out
end

#=
    CSDP
=#

# ]activate soi
using CSDP

MOIU.@model(CSDP_ModelData,
            (),
            (MOI.EqualTo, MOI.GreaterThan, MOI.LessThan),
            (MOI.Zeros, MOI.Nonnegatives, MOI.Nonpositives,
             MOI.PositiveSemidefiniteConeTriangle),
            (),
            (MOI.SingleVariable,),
            (MOI.ScalarAffineFunction,),
            (MOI.VectorOfVariables,),
            (MOI.VectorAffineFunction,))

const optimizer = 
MOIB.full_bridge_optimizer(
    MOIU.CachingOptimizer(
        CSDP_ModelData{Float64}(), CSDP.Optimizer()),#printlevel=0))
    Float64
    )

MOI.empty!(optimizer)

#=
    Problem Data
=#

n = 25
L = ones(n+1, n+1)

#=
    Optimization Problem
=#

nvars = sympackedlen(n + 1)

X = MOI.add_variables(optimizer, nvars)

#
# X is PSD
#
vov = MOI.VectorOfVariables(X)
cX = MOI.add_constraint(optimizer, vov, MOI.PositiveSemidefiniteConeTriangle(n+1))

#
# -1 <= X[i,j] <= 1
#
MOI.add_constraint(optimizer,
                    MOI.VectorAffineFunction(
                        MOI.VectorAffineTerm.(
                            collect(1:nvars), MOI.ScalarAffineTerm.(1.0, X)),
                        -ones(nvars)),
                    MOI.Nonpositives(nvars))
MOI.add_constraint(optimizer,
                    MOI.VectorAffineFunction(
                        MOI.VectorAffineTerm.(
                            collect(1:nvars), MOI.ScalarAffineTerm.(1.0, X)),
                        ones(nvars)),
                    MOI.Nonnegatives(nvars))

# define a square
Xsq = Matrix{MOI.VariableIndex}(undef, n+1,n+1)
ivech!(Xsq, X)
Xsq = Matrix(Symmetric(Xsq,:U))

#
# diag(X) .== 1
#
MOI.add_constraint(optimizer,
                    MOI.VectorAffineFunction(
                        MOI.VectorAffineTerm.(
                            collect(1:n+1), MOI.ScalarAffineTerm.(1.0, [Xsq[i,i] for i in 1:n+1])),
                        -ones(n+1)),
                    MOI.Zeros(n+1))

#
# Min sum(L[i,j], X[i,j])
#
objf_t = vec([MOI.ScalarAffineTerm(L[i,j], Xsq[i,j]) for i in 1:n+1, j in 1:n+1])
MOI.set(optimizer, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), MOI.ScalarAffineFunction(objf_t, 0.0))

MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MIN_SENSE)

MOI.optimize!(optimizer)