Refactor some constraint into GenericConstraint (#590)

* Refactor Constraint

* Refactor into GenericConstraint

* up

* Fix tests

* Fixes

* Fix char

* Fixes

* Fix format

* Apply suggestions from code review

* Update

* Relocate

* Update

* Add support for Zeros

* Add SecondOrderCone support

* Update

* Update changelog.md

---------

Co-authored-by: Oscar Dowson <[email protected]>
Co-authored-by: odow <[email protected]>

Author: blegat

docs/src/changelog.md

@@ -27,6 +27,9 @@ changes.
       dual variable associated with such constraints is now reversed in sign**.
       (Following the convention in MathOptInterface, the dual of `a <= b` is
       always negative, regardless of optimization sense.) (#593)
+    * The structs `LessThanConstraint`, `GreaterThanConstraint`, `EqualToConstraint`
+      `SecondOrderConeConstraint` and `PostiveSemidefinteConeConstraint` have
+      been replaced by `GenericConstraint{S}` where `S<:MOI.AbstractSet` (#590)
  * The syntaxes `dot(*)`, `dot(/)` and `dot(^)` have been removed in favor of
    explicit broadcasting (`x .* y`, `x ./ y`, and `x .^ y`). These were (mild)
    type piracy. In addition, `vecdot(x,y)` has been removed. Call

src/atoms/second_order_cone/QolElemAtom.jl

@@ -30,10 +30,8 @@ function new_conic_form!(context::Context{T}, q::QolElemAtom) where {T}
     x, y = q.children
     t = Variable(x.size)
     for i in 1:length(x)
-        add_constraint!(
-            context,
-            SecondOrderConeConstraint(y[i] + t[i], y[i] - t[i], 2 * x[i]),
-        )
+        f = vcat(y[i] + t[i], y[i] - t[i], 2 * x[i])
+        add_constraint!(context, GenericConstraint{MOI.SecondOrderCone}(f))
     end
     add_constraint!(context, y >= 0)
     return conic_form!(context, t)

src/atoms/second_order_cone/QuadOverLinAtom.jl

@@ -30,7 +30,8 @@ end
 function new_conic_form!(context::Context, q::QuadOverLinAtom)
     t = Variable()
     x, y = q.children
-    add_constraint!(context, SecondOrderConeConstraint(y + t, y - t, 2 * x))
+    f = vcat(y + t, y - t, 2 * x)
+    add_constraint!(context, GenericConstraint{MOI.SecondOrderCone}(f))
     add_constraint!(context, y >= 0)
     return conic_form!(context, t)
 end

src/constraints/EqualToConstraint.jl

@@ -0,0 +1,245 @@
+mutable struct GenericConstraint{S<:MOI.AbstractSet} <: Constraint
+    child::AbstractExpr
+    set::S
+    dual::Union{Value,Nothing}
+
+    function GenericConstraint(child::AbstractExpr, set::MOI.AbstractSet)
+        return new{typeof(set)}(child, set, nothing)
+    end
+end
+
+function GenericConstraint{S}(child::AbstractExpr) where {S<:MOI.AbstractSet}
+    return GenericConstraint(child, set_with_size(S, size(child)))
+end
+
+head(io::IO, c::GenericConstraint) = head(io, c.set)
+
+# A default fallback that skips the feasibiltiy check.
+is_feasible(f, ::MOI.AbstractSet, tol) = true
+
+AbstractTrees.children(c::GenericConstraint) = (c.child,)
+
+vexity(c::GenericConstraint) = vexity(vexity(c.child), c.set)
+
+function _add_constraint!(context::Context, c::GenericConstraint)
+    if vexity(c.child) == ConstVexity()
+        x = evaluate(c.child)
+        if !is_feasible(x, c.set, CONSTANT_CONSTRAINT_TOL[])
+            context.detected_infeasible_during_formulation[] = true
+        end
+        return
+    end
+    f = conic_form!(context, c.child)
+    context.constr_to_moi_inds[c] = MOI_add_constraint(context.model, f, c.set)
+    return
+end
+
+function populate_dual!(model::MOI.ModelLike, c::GenericConstraint, indices)
+    ret = MOI.get(model, MOI.ConstraintDual(), indices)
+    c.dual = output(reshape(ret, c.child.size))
+    return
+end
+
+function populate_dual!(
+    model::MOI.ModelLike,
+    c::GenericConstraint,
+    indices::NTuple{2},
+)
+    re = MOI.get(model, MOI.ConstraintDual(), indices[1])
+    imag = MOI.get(model, MOI.ConstraintDual(), indices[2])
+    c.dual = output(reshape(re + im * imag, c.child.size))
+    return
+end
+
+function _promote_size(lhs::AbstractExpr, rhs::AbstractExpr)
+    if lhs.size == rhs.size || lhs.size == (1, 1)
+        sz = rhs.size
+        if lhs.size == (1, 1) && rhs.size != (1, 1)
+            lhs = lhs * ones(rhs.size)
+        end
+    elseif rhs.size == (1, 1)
+        sz = lhs.size
+        if rhs.size == (1, 1) && lhs.size != (1, 1)
+            rhs = rhs * ones(lhs.size)
+        end
+    else
+        error(
+            "Cannot create constraint between expressions of size " *
+            "$(lhs.size) and $(rhs.size)",
+        )
+    end
+    return lhs, rhs
+end
+
+# ==============================================================================
+#     Nonnegatives
+# ==============================================================================
+
+function set_with_size(::Type{MOI.Nonnegatives}, sz::Tuple{Int,Int})
+    return MOI.Nonnegatives(prod(sz))
+end
+
+head(io::IO, ::MOI.Nonnegatives) = print(io, "≥")
+
+is_feasible(f, ::MOI.Nonnegatives, tol) = all(f .>= -tol)
+
+function vexity(vex, ::MOI.Nonnegatives)
+    if vex == ConvexVexity()
+        return NotDcp()
+    elseif vex == ConcaveVexity()
+        return ConvexVexity()
+    end
+    return vex
+end
+
+function Base.:>=(lhs::AbstractExpr, rhs::AbstractExpr)
+    if sign(lhs) == ComplexSign() || sign(rhs) == ComplexSign()
+        error(
+            "Cannot create constraint between expressions of sign " *
+            "$(sign(lhs)) and $(sign(rhs))",
+        )
+    end
+    lhs, rhs = _promote_size(lhs, rhs)
+    return GenericConstraint{MOI.Nonnegatives}(lhs - rhs)
+end
+
+Base.:>=(lhs::AbstractExpr, rhs::Value) = >=(lhs, constant(rhs))
+
+Base.:>=(lhs::Value, rhs::AbstractExpr) = >=(constant(lhs), rhs)
+
+# ==============================================================================
+#     Nonnpositives
+# ==============================================================================
+
+function set_with_size(::Type{MOI.Nonpositives}, sz::Tuple{Int,Int})
+    return MOI.Nonpositives(prod(sz))
+end
+
+head(io::IO, ::MOI.Nonpositives) = print(io, "≤")
+
+is_feasible(f, ::MOI.Nonpositives, tol) = all(f .<= tol)
+
+function vexity(vex, ::MOI.Nonpositives)
+    if vex == ConcaveVexity()
+        return NotDcp()
+    end
+    return vex
+end
+
+function Base.:<=(lhs::AbstractExpr, rhs::AbstractExpr)
+    if sign(lhs) == ComplexSign() || sign(rhs) == ComplexSign()
+        error(
+            "Cannot create constraint between expressions of sign " *
+            "$(sign(lhs)) and $(sign(rhs))",
+        )
+    end
+    lhs, rhs = _promote_size(lhs, rhs)
+    return GenericConstraint{MOI.Nonpositives}(lhs - rhs)
+end
+
+Base.:<=(lhs::AbstractExpr, rhs::Value) = <=(lhs, constant(rhs))
+
+Base.:<=(lhs::Value, rhs::AbstractExpr) = <=(constant(lhs), rhs)
+
+# ==============================================================================
+#     Zeros
+# ==============================================================================
+
+function set_with_size(::Type{MOI.Zeros}, sz::Tuple{Int,Int})
+    return MOI.Zeros(prod(sz))
+end
+
+head(io::IO, ::MOI.Zeros) = print(io, "==")
+
+is_feasible(f, ::MOI.Zeros, tol) = all(abs.(f) .<= tol)
+
+function vexity(vex, ::MOI.Zeros)
+    if vex == ConvexVexity() || vex == ConcaveVexity()
+        return NotDcp()
+    end
+    return vex
+end
+
+function Base.:(==)(lhs::AbstractExpr, rhs::AbstractExpr)
+    lhs, rhs = _promote_size(lhs, rhs)
+    return GenericConstraint{MOI.Zeros}(lhs - rhs)
+end
+
+Base.:(==)(lhs::AbstractExpr, rhs::Value) = ==(lhs, constant(rhs))
+
+Base.:(==)(lhs::Value, rhs::AbstractExpr) = ==(constant(lhs), rhs)
+
+# ==============================================================================
+#     PositiveSemidefiniteConeSquare
+# ==============================================================================
+
+function set_with_size(
+    ::Type{MOI.PositiveSemidefiniteConeSquare},
+    sz::Tuple{Int,Int},
+)
+    if sz[1] != sz[2]
+        error("Positive semidefinite expressions must be square")
+    end
+    return MOI.PositiveSemidefiniteConeSquare(sz[1])
+end
+
+head(io::IO, ::MOI.PositiveSemidefiniteConeSquare) = print(io, "sdp")
+
+function vexity(vex, ::MOI.PositiveSemidefiniteConeSquare)
+    if !(vex in (AffineVexity(), ConstVexity()))
+        return NotDcp()
+    end
+    return AffineVexity()
+end
+
+function is_feasible(x, ::MOI.PositiveSemidefiniteConeSquare, tol)
+    return x ≈ transpose(x) && LinearAlgebra.eigmin(x) >= -tol
+end
+
+function LinearAlgebra.isposdef(x::AbstractExpr)
+    if iscomplex(x)
+        return GenericConstraint{MOI.PositiveSemidefiniteConeSquare}(
+            [real(x) -imag(x); imag(x) real(x)],
+        )
+    end
+    return GenericConstraint{MOI.PositiveSemidefiniteConeSquare}(x)
+end
+
+⪰(x::AbstractExpr, y::AbstractExpr) = isposdef(x - y)
+
+function ⪰(x::AbstractExpr, y::Value)
+    if all(y .== 0)
+        return isposdef(x)
+    end
+    return isposdef(x - constant(y))
+end
+
+function ⪰(x::Value, y::AbstractExpr)
+    if all(x .== 0)
+        return isposdef(-y)
+    end
+    return isposdef(constant(x) - y)
+end
+
+⪯(x::AbstractExpr, y::AbstractExpr) = ⪰(y, x)
+
+⪯(x::Value, y::AbstractExpr) = ⪰(y, x)
+
+⪯(x::AbstractExpr, y::Value) = ⪰(y, x)
+
+# ==============================================================================
+#     SecondOrderCone
+# ==============================================================================
+
+function set_with_size(::Type{MOI.SecondOrderCone}, sz::Tuple{Int,Int})
+    return MOI.SecondOrderCone(prod(sz))
+end
+
+head(io::IO, ::MOI.SecondOrderCone) = print(io, "soc")
+
+function vexity(vex, ::MOI.SecondOrderCone)
+    if !(vex == ConstVexity() || vex == AffineVexity())
+        return NotDcp()
+    end
+    return ConvexVexity()
+end