Integrate with JuMP

Author: blegat

Project.toml

@@ -1,14 +1,16 @@
 name = "ComplexOptInterface"
 uuid = "25c3070e-1275-41b5-afef-2f982c87090a"
-repo = "https://github.com/jump-dev/ComplexOptInterface.jl.git"
 authors = ["Benoît Legat <[email protected]>"]
+repo = "https://github.com/jump-dev/ComplexOptInterface.jl.git"
 version = "0.0.3"
 
 [deps]
+JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
 MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
 
 [compat]
+JuMP = "0.23"
 MathOptInterface = "1"
 MutableArithmetics = "1"
 julia = "1.6"

README.md

@@ -8,6 +8,23 @@
 | [![Codecov branch][codecov-img]][codecov-url] | [<img src="https://upload.wikimedia.org/wikipedia/en/a/af/Discourse_logo.png" width="64">][discourse-url] |
 
 An extension of MathOptInterface to complex sets.
+
+## How to use
+
+Add to the JuMP model the bridges of this package with `add_all_bridges` and then you can create complex equality constraints and complex hermitian matrices:
+```julia
+using JuMP
+import ComplexOptInterface
+const COI = ComplexOptInterface
+
+model = Model()
+COI.add_all_bridges(model)
+@variable(model, x[1:2, 1:2] in COI.HermitianPSDCone())
+@constraint(model, x[1, 1] + x[2, 2] * im == 1 + 2im)
+```
+
+## Design considerations
+
 There are two types of complex sets:
 1) Some sets have complex entries, i.e. `MOI.EqualTo{Complex{Float64}}`,
 2) Some sets have real entries even if they model a complex mathematical set, i.e.

src/ComplexOptInterface.jl

@@ -9,11 +9,14 @@ const MOI = MathOptInterface
 struct HermitianPositiveSemidefiniteConeTriangle <: MOI.AbstractVectorSet
     side_dimension::Int
 end
+
 function MOI.dimension(set::HermitianPositiveSemidefiniteConeTriangle)
     return MOI.dimension(MOI.PositiveSemidefiniteConeTriangle(set.side_dimension)) +
            MOI.dimension(MOI.PositiveSemidefiniteConeTriangle(set.side_dimension - 1))
 end
 
+Base.copy(set::HermitianPositiveSemidefiniteConeTriangle) = set
+
 function MOI.Utilities.set_dot(
     x::AbstractVector,
     y::AbstractVector,
@@ -28,5 +31,6 @@ function MOI.Utilities.set_dot(
 end
 
 include("Bridges/Bridges.jl")
+include("jump.jl")
 
 end # module

src/jump.jl

@@ -0,0 +1,103 @@
+import LinearAlgebra
+import JuMP
+
+function add_all_bridges(model::JuMP.Model)
+    JuMP.add_bridge(model, Bridges.Variable.HermitianToSymmetricPSDBridge)
+    JuMP.add_bridge(model, Bridges.Constraint.SplitEqualToBridge)
+    JuMP.add_bridge(model, Bridges.Constraint.SplitZeroBridge)
+end
+
+struct HermitianPSDCone end
+
+struct HermitianMatrixShape <: JuMP.AbstractShape
+    side_dimension::Int
+end
+
+function JuMP.vectorize(matrix::Matrix, ::HermitianMatrixShape)
+    n = LinearAlgebra.checksquare(matrix)
+    return vcat(
+        JuMP.vectorize(_real.(matrix), JuMP.SymmetricMatrixShape(n)),
+        JuMP.vectorize(_imag.(matrix[1:(end-1),2:end]), JuMP.SymmetricMatrixShape(n - 1)),
+    )
+end
+
+function JuMP.reshape_vector(v::Vector{T}, shape::HermitianMatrixShape) where T
+    NewType = MA.promote_operation(MA.add_mul, T, Complex{Bool}, T)
+    n = shape.side_dimension
+    matrix = Matrix{NewType}(undef, n, n)
+    real_k = 0
+    imag_k = MOI.dimension(MOI.PositiveSemidefiniteConeTriangle(n))
+    for j in 1:n
+        for i in 1:(j-1)
+            real_k += 1
+            imag_k += 1
+            matrix[i, j] = v[real_k] + im * v[imag_k]
+            matrix[j, i] = v[real_k] - im * v[imag_k]
+        end
+        real_k += 1
+        matrix[j, j] = v[real_k]
+    end
+    return matrix
+end
+
+function _mapinfo(f::Function, v::JuMP.ScalarVariable)
+    info = v.info
+    return JuMP.ScalarVariable(JuMP.VariableInfo(
+        info.has_lb,
+        f(info.lower_bound),
+        info.has_ub,
+        f(info.upper_bound),
+        info.has_fix,
+        f(info.fixed_value),
+        info.has_start,
+        f(info.start),
+        info.binary,
+        info.integer,
+    ))
+end
+
+_real(s::String) = s
+_imag(s::String) = s
+
+_real(v::JuMP.ScalarVariable) = _mapinfo(real, v)
+_imag(v::JuMP.ScalarVariable) = _mapinfo(imag, v)
+_conj(v::JuMP.ScalarVariable) = _mapinfo(conj, v)
+function _isreal(v::JuMP.ScalarVariable)
+    info = v.info
+    return isreal(info.lower_bound) && isreal(info.upper_bound) && isreal(info.fixed_value) && isreal(info.start)
+end
+
+function _vectorize_variables(_error::Function, matrix::Matrix)
+    n = LinearAlgebra.checksquare(matrix)
+    for j in 1:n
+        if !_isreal(matrix[j, j])
+            _error(
+                "Non-real bounds or starting values for diagonal of hermitian variable.",
+            )
+        end
+        for i in 1:j
+            if matrix[i, j] != _conj(matrix[j, i])
+                _error(
+                    "Non-conjugate bounds, integrality or starting values for hermitian variable.",
+                )
+            end
+        end
+    end
+    JuMP.vectorize(matrix, HermitianMatrixShape(n))
+end
+
+function JuMP.build_variable(
+    _error::Function,
+    variables::Matrix{<:JuMP.AbstractVariable},
+    ::HermitianPSDCone,
+)
+    n = JuMP._square_side(_error, variables)
+    set = HermitianPositiveSemidefiniteConeTriangle(n)
+    shape = HermitianMatrixShape(n)
+    return JuMP.VariablesConstrainedOnCreation(
+        _vectorize_variables(_error, variables),
+        set,
+        shape,
+    )
+end
+

test/jump.jl

@@ -0,0 +1,69 @@
+module TestJuMP
+
+using JuMP
+import GLPK
+import ComplexOptInterface
+const COI = ComplexOptInterface
+using Test
+
+function test_readme_example()
+    model = Model()
+    COI.add_all_bridges(model)
+    @variable(model, x[1:2, 1:2] in COI.HermitianPSDCone())
+    @test num_variables(model) == 4
+    v = all_variables(model)
+    @test x[1, 1] == 1v[1]
+    @test x[1, 2] == v[2] + v[4] * im
+    @test x[2, 2] == 1v[3]
+    # FIXME needs https://github.com/jump-dev/JuMP.jl/pull/2899
+    #@test x[2, 1] == conj(x[1, 2])
+    @constraint(model, x[1, 1] + x[2, 2] * im == 1 + 2im)
+end
+
+function test_simple_equality()
+    model = Model(GLPK.Optimizer)
+    add_bridge(model, COI.Bridges.Constraint.SplitEqualToBridge)
+
+    @variable(model, 0 <= x[1:3] <= 1)
+
+    @objective(model, Max, x[1] + 2x[2] - x[3])
+
+    # FIXME needs https://github.com/jump-dev/JuMP.jl/pull/2895
+    #@constraint(model, (1 + im) * x[1] + (1 + im) * x[2] + (1 - 3im) * x[3] == 1.0)
+    aff = (1 + im) * x[1] + (1 + im) * x[2] + (1 - 3im) * x[3]
+    @constraint(model, aff == 1.0)
+
+    # Optimize the model:
+
+    optimize!(model)
+
+    # Check the termination and primal status to see if we have a solution:
+
+    println("Termination status : ", termination_status(model))
+    println("Primal status      : ", primal_status(model))
+
+    # Print the solution:
+
+    obj_value = objective_value(model)
+    x_values = value.(x)
+
+    println("Objective value : ", obj_value)
+    println("x values        : ", x_values)
+
+    @test obj_value ≈ 1.25
+    @test x_values ≈ [0.0, 0.75, 0.25]
+end
+
+function runtests()
+    for name in names(@__MODULE__; all = true)
+        if startswith("$name", "test_")
+            @testset "$(name)" begin
+                getfield(@__MODULE__, name)()
+            end
+        end
+    end
+end
+
+end
+
+TestJuMP.runtests()