Add solution inspection to example (#24)

* Add solution inspection to example

* Fix format

* Finish README

* Fix format

* Add tests

* Fix format

Author: blegat

README.md

@@ -43,6 +43,8 @@ adapted to exploit the low-rank constraints.
 julia> include(joinpath(dirname(dirname(pathof(LowRankOpt))), "examples", "maxcut.jl"))
 maxcut (generic function with 2 methods)
 
+julia> weights = [0 5 7 6; 5 0 0 1; 7 0 0 1; 6 1 1 0];
+
 julia> model = maxcut(weights, SDPLR.Optimizer);
 
 julia> optimize!(model)
@@ -68,4 +70,71 @@ DIMACS error measures: 5.86e-06 0.00e+00 0.00e+00 0.00e+00 2.25e-05 9.84e-06
 
 julia> objective_value(model)
 17.99998881724702
+
+julia> con_ref = VariableInSetRef(model[:dot_prod_set]);
+```
+
+Use `LRO.InnerAttribute` to request the `SDPLR.Factor` attribute.
+```julia-repl
+julia> MOI.get(model, LRO.InnerAttribute(SDPLR.Factor()), VariableInSetRef(model[:dot_prod_set]))
+4×3 Matrix{Float64}:
+  0.949505   0.31101    0.0414433
+ -0.950269  -0.308646  -0.0415714
+ -0.949503  -0.311095  -0.0406689
+ -0.948855  -0.312898  -0.0420402
+```
+We can see that SDPLR decided to search for a solution of rank at most 3.
+To check if SDPLR found an optimal solution, we need to check whether the dual solution is feasible.
+This can be achieved as follows:
+```julia-repl
+julia> dual_set = MOI.dual_set(constraint_object(con_ref).set);
+
+julia> MOI.Utilities.distance_to_set(dual(con_ref), dual_set)
+1.602881211577939e-8
+```
+
+For the MAX-CUT problem, we know there exists a rank-1 solution where
+the entries of the factor are `-1` or `1` depending on the side of the cut
+the nodes are on. Let's now be greedy and search for a solution of rank-1.
+```julia
+julia> set_attribute(model, "maxrank", (m, n) -> 1)
+
+julia> optimize!(model)
+
+            ***   SDPLR 1.03-beta   ***
+
+===================================================
+ major   minor        val        infeas      time  
+---------------------------------------------------
+    1        0   2.68869170e-01  8.7e-01       0
+    2        7   8.10081717e+01  1.0e+01       0
+    3       10  -1.81587378e+01  2.3e-01       0
+    4       11  -1.80062673e+01  1.3e-01       0
+    5       13  -1.79951904e+01  5.0e-02       0
+    6       15  -1.79981240e+01  1.4e-02       0
+    7       16  -1.79998759e+01  2.2e-03       0
+    8       17  -1.79999980e+01  1.9e-04       0
+    9       18  -1.80000000e+01  1.4e-05       0
+   10       19  -1.80000000e+01  7.1e-07       0
+===================================================
+
+DIMACS error measures: 7.12e-07 0.00e+00 0.00e+00 1.31e-05 1.59e-06 -7.81e-06
+
+
+julia> MOI.get(model, LRO.InnerAttribute(SDPLR.Factor()), VariableInSetRef(model[:dot_prod_set]))
+4×1 Matrix{Float64}:
+ -0.9999998713637199
+  0.9999995531365233
+  0.9999997036204064
+  0.9999995498147483
+```
+
+Note that even though a solution of rank-1 exists, SDPLR is more likely to
+converge to a spurious local minimum if we use a lower-rank, so we should
+be careful before claiming that we found the optimal solution.
+Luckily, the following shows that the dual is feasible which gives a certificate
+of primal optimality.
+```julia-repl
+julia> MOI.Utilities.distance_to_set(dual(con_ref), dual_set)
+0.0
 ```

examples/maxcut.jl

@@ -1,4 +1,11 @@
 using LinearAlgebra, JuMP, LowRankOpt
+import LowRankOpt as LRO
+
+function e_i(i, n)
+    ei = zeros(n)
+    ei[i] = 1
+    return ei
+end
 
 function maxcut(weights, solver)
     N = LinearAlgebra.checksquare(weights)

src/Bridges/Constraint/Constraint.jl

@@ -12,6 +12,14 @@ for filename in readdir(joinpath(@__DIR__, "bridges"); join = true)
     include(filename)
 end
 
+function MOI.get(
+    model::MOI.ModelLike,
+    attr::LRO.InnerAttribute,
+    bridge::MOI.Bridges.Constraint.SetMapBridge,
+)
+    return MOI.get(model, attr, bridge.constraint)
+end
+
 const ConversionBridge{W,T} = MOI.Bridges.Constraint.SetConversionBridge{
     T,
     LRO.LinearCombinationInSet{

src/Bridges/Variable/Variable.jl

@@ -12,6 +12,14 @@ for filename in readdir(joinpath(@__DIR__, "bridges"); join = true)
     include(filename)
 end
 
+function MOI.get(
+    model::MOI.ModelLike,
+    attr::LRO.InnerAttribute,
+    bridge::MOI.Bridges.Variable.SetMapBridge,
+)
+    return MOI.get(model, attr, bridge.constraint)
+end
+
 const ConversionBridge{W,T} = SetConversionBridge{
     T,
     LRO.SetDotProducts{

src/LowRankOpt.jl

@@ -11,6 +11,8 @@ import MathOptInterface as MOI
 import FillArrays
 
 include("sets.jl")
+include("attributes.jl")
+include("distance_to_set.jl")
 include("Test/Test.jl")
 include("Bridges/Bridges.jl")
 

src/attributes.jl

@@ -0,0 +1,35 @@
+"""
+    struct InnerAttribute{A<:MOI.AbstractConstraintAttribute} <: MOI.AbstractConstraintAttribute
+        inner::A
+    end
+
+This attributes goes through all `Bridges.Variable.SetMapBridge` and `Bridges.Constraint.SetMapBridge`, ignoring the corresponding linear transformation.
+"""
+struct InnerAttribute{A<:MOI.AbstractConstraintAttribute} <:
+       MOI.AbstractConstraintAttribute
+    inner::A
+end
+
+function MOI.is_copyable(a::InnerAttribute)
+    return MOI.is_copyable(a.inner)
+end
+
+function MOI.is_set_by_optimize(a::InnerAttribute)
+    return MOI.is_set_by_optimize(a.inner)
+end
+
+function MOI.get(
+    model::MOI.Utilities.UniversalFallback,
+    attr::InnerAttribute,
+    ci::MOI.ConstraintIndex,
+)
+    return MOI.get(model, attr.inner, ci)
+end
+
+function MOI.get_fallback(
+    model::MOI.ModelLike,
+    attr::InnerAttribute,
+    ci::MOI.ConstraintIndex,
+)
+    return MOI.get(model, attr.inner, ci)
+end

src/distance_to_set.jl

@@ -0,0 +1,28 @@
+function MOI.Utilities.distance_to_set(
+    d::MOI.Utilities.ProjectionUpperBoundDistance,
+    x::AbstractVector{T},
+    set::SetDotProducts{WITH_SET},
+) where {T,W}
+    MOI.Utilities._check_dimension(x, set)
+    n = length(set.vectors)
+    vec = x[(n+1):end]
+    init = MOI.Utilities.distance_to_set(d, vec, set.set)^2
+    return √sum(1:n; init) do i
+        return (x[i] - MOI.Utilities.set_dot(set.vectors[i], vec, set.set))^2
+    end
+end
+
+function MOI.Utilities.distance_to_set(
+    d::MOI.Utilities.ProjectionUpperBoundDistance,
+    x::AbstractVector{T},
+    set::LinearCombinationInSet{W},
+) where {T,W}
+    MOI.Utilities._check_dimension(x, set)
+    if W == WITH_SET
+        init = x[(length(set.vectors)+1):end]
+    else
+        init = zeros(T, MOI.dimension(set.set))
+    end
+    y = sum(x[i] * set.vectors[i] for i in eachindex(set.vectors); init)
+    return MOI.Utilities.distance_to_set(d, y, set.set)
+end

src/sets.jl

@@ -120,6 +120,26 @@ function MOI.dual_set_type(::Type{LinearCombinationInSet{W,S,V}}) where {W,S,V}
     return SetDotProducts{W,S,V}
 end
 
+function MOI.Utilities.set_dot(
+    x::AbstractVector,
+    y::AbstractVector,
+    set::Union{SetDotProducts{WITH_SET},LinearCombinationInSet{WITH_SET}},
+)
+    n = length(set.vectors)
+    return LinearAlgebra.dot(x[1:n], y[1:n]) +
+           MOI.Utilities.set_dot(x[(n+1):end], y[(n+1):end], set.set)
+end
+
+function MOI.Utilities.dot_coefficients(
+    x::AbstractVector,
+    set::Union{SetDotProducts{WITH_SET},LinearCombinationInSet{WITH_SET}},
+)
+    c = copy(x)
+    n = length(set.vectors)
+    c[(n+1):end] = MOI.Utilities.dot_coefficients(x[(n+1):end], set.set)
+    return c
+end
+
 abstract type AbstractFactorization{T,F} <: AbstractMatrix{T} end
 
 function Base.size(m::AbstractFactorization)

test/Bridges/Constraint/LinearCombinationBridge.jl

@@ -21,27 +21,30 @@ function runtests()
     return
 end
 
+function _model(T, model)
+    _, cx = MOI.add_constrained_variables(
+        model,
+        LRO.LinearCombinationInSet{LRO.WITH_SET}(
+            MOI.PositiveSemidefiniteConeTriangle(2),
+            LRO.TriangleVectorization.([
+                T[
+                    1 2
+                    2 3
+                ],
+                T[
+                    4 5
+                    5 6
+                ],
+            ]),
+        ),
+    )
+    return cx
+end
+
 function test_psd(T::Type)
     MOI.Bridges.runtests(
         LRO.Bridges.Constraint.LinearCombinationBridge,
-        model -> begin
-            x, _ = MOI.add_constrained_variables(
-                model,
-                LRO.LinearCombinationInSet{LRO.WITH_SET}(
-                    MOI.PositiveSemidefiniteConeTriangle(2),
-                    LRO.TriangleVectorization.([
-                        T[
-                            1 2
-                            2 3
-                        ],
-                        T[
-                            4 5
-                            5 6
-                        ],
-                    ]),
-                ),
-            )
-        end,
+        Base.Fix1(_model, T),
         model -> begin
             x = MOI.add_variables(model, 5)
             MOI.add_constraint(
@@ -61,6 +64,25 @@ function test_psd(T::Type)
     return
 end
 
+struct Custom <: MOI.AbstractConstraintAttribute end
+MOI.is_set_by_optimize(::Custom) = false
+
+function test_attribute(T::Type)
+    inner = MOI.Utilities.UniversalFallback(MOI.Utilities.Model{T}())
+    model = MOI.Bridges._bridged_model(
+        LRO.Bridges.Constraint.LinearCombinationBridge{T},
+        inner,
+    )
+    cx = _model(T, model)
+    F = MOI.VectorAffineFunction{T}
+    S = MOI.PositiveSemidefiniteConeTriangle
+    ci = only(MOI.get(inner, MOI.ListOfConstraintIndices{F,S}()))
+    MOI.set(inner, Custom(), ci, "test")
+    @test MOI.get(inner, Custom(), ci) == "test"
+    @test MOI.get(inner, LRO.InnerAttribute(Custom()), ci) == "test"
+    @test MOI.get(model, LRO.InnerAttribute(Custom()), cx) == "test"
+end
+
 end  # module
 
 TestConstraintLinearCombination.runtests()

test/Bridges/Variable/set_dot.jl

@@ -21,29 +21,32 @@ function runtests()
     return
 end
 
+function _model(T, model)
+    x, cx = MOI.add_constrained_variables(
+        model,
+        LRO.SetDotProducts{LRO.WITH_SET}(
+            MOI.PositiveSemidefiniteConeTriangle(2),
+            LRO.TriangleVectorization.([
+                T[
+                    1 2
+                    2 3
+                ],
+                T[
+                    4 5
+                    5 6
+                ],
+            ]),
+        ),
+    )
+    MOI.add_constraint(model, one(T) * x[1], MOI.EqualTo(zero(T)))
+    MOI.add_constraint(model, one(T) * x[2], MOI.LessThan(zero(T)))
+    return cx
+end
+
 function test_psd(T::Type)
     MOI.Bridges.runtests(
         LRO.Bridges.Variable.DotProductsBridge,
-        model -> begin
-            x, _ = MOI.add_constrained_variables(
-                model,
-                LRO.SetDotProducts{LRO.WITH_SET}(
-                    MOI.PositiveSemidefiniteConeTriangle(2),
-                    LRO.TriangleVectorization.([
-                        T[
-                            1 2
-                            2 3
-                        ],
-                        T[
-                            4 5
-                            5 6
-                        ],
-                    ]),
-                ),
-            )
-            MOI.add_constraint(model, one(T) * x[1], MOI.EqualTo(zero(T)))
-            MOI.add_constraint(model, one(T) * x[2], MOI.LessThan(zero(T)))
-        end,
+        Base.Fix1(_model, T),
         model -> begin
             Q, _ = MOI.add_constrained_variables(
                 model,
@@ -66,6 +69,43 @@ function test_psd(T::Type)
     return
 end
 
+struct Custom <: MOI.AbstractConstraintAttribute
+    is_copyable::Bool
+    is_set_by_optimize::Bool
+end
+MOI.is_copyable(c::Custom) = c.is_copyable
+MOI.is_set_by_optimize(c::Custom) = c.is_set_by_optimize
+
+function test_attribute(T::Type)
+    inner = MOI.Utilities.UniversalFallback(MOI.Utilities.Model{T}())
+    model = MOI.Bridges._bridged_model(
+        LRO.Bridges.Variable.DotProductsBridge{T},
+        inner,
+    )
+    cx = _model(T, model)
+    F = MOI.VectorOfVariables
+    S = MOI.PositiveSemidefiniteConeTriangle
+    ci = only(MOI.get(inner, MOI.ListOfConstraintIndices{F,S}()))
+    attr = Custom(true, false)
+    MOI.set(inner, attr, ci, "test")
+    @test MOI.get(inner, attr, ci) == "test"
+    attr = LRO.InnerAttribute(attr)
+    @test MOI.get(inner, attr, ci) == "test"
+    @test MOI.get(model, attr, cx) == "test"
+    for is_copyable in [false, true]
+        for is_set_by_optimize in [false, true]
+            attr = LRO.InnerAttribute(Custom(is_copyable, is_set_by_optimize))
+            @test MOI.is_copyable(attr) == is_copyable
+            @test MOI.is_set_by_optimize(attr) == is_set_by_optimize
+        end
+    end
+    @test_throws MOI.GetAttributeNotAllowed{Custom} MOI.get(
+        inner.model,
+        attr,
+        ci,
+    )
+end
+
 end  # module
 
 TestVariableDotProducts.runtests()