Merge pull request #82 from psrenergy/px/indicator

Add support for indicator variables

src/compiler/constraints.jl

@@ -6,7 +6,12 @@ function constraints!(model::Virtual.Model, arch::AbstractArchitecture)
     return nothing
 end
 
-function constraints!(model::Virtual.Model, ::Type{F}, ::Type{S}, arch::AbstractArchitecture) where {F,S}
+function constraints!(
+    model::Virtual.Model,
+    ::Type{F},
+    ::Type{S},
+    arch::AbstractArchitecture,
+) where {F,S}
     for ci in MOI.get(model, MOI.ListOfConstraintIndices{F,S}())
         f = MOI.get(model, MOI.ConstraintFunction(), ci)
         s = MOI.get(model, MOI.ConstraintSet(), ci)
@@ -170,6 +175,20 @@ function constraint(
     return g^2
 end
 
+@doc raw"""
+
+"""
+function constraint(
+    model::Virtual.Model{T},
+    ci::CI,
+    f::SAF{T},
+    s::MOI.Interval{T},
+    arch::AbstractArchitecture,
+) where {T}
+    return constraint(model, ci, f, LT{T}(s.upper), arch) +
+           constraint(model, ci, f, GT{T}(s.lower), arch)
+end
+
 @doc raw"""
     constraint(
         model::Virtual.Model{T}, 
@@ -468,7 +487,8 @@ function constraint(
             for (ω, _) in Virtual.expansion(vi)
                 g[ω] = one(T)
             end
-        elseif Virtual.encoding(vi) isa Encoding.OneHot || Virtual.encoding(vi) isa Encoding.DomainWall
+        elseif Virtual.encoding(vi) isa Encoding.OneHot ||
+               Virtual.encoding(vi) isa Encoding.DomainWall
             ξ = Virtual.expansion(vi)
             a = ξ[nothing]
 
@@ -509,3 +529,74 @@ function constraint(
 
     return g^2 + h
 end
+
+function constraint(
+    model::Virtual.Model{T},
+    ci::CI,
+    f::MOI.VectorAffineFunction{T},
+    s::MOI.Indicator{A,S},
+    arch::AbstractArchitecture,
+) where {T,A,S}
+    # Indicator Constraint: y = 0|1 => {g(x)}
+
+    xi = first(f.terms).scalar_term.variable # Indicator Variable
+    vi = model.source[xi]
+
+    @assert Virtual.encoding(vi) isa Encoding.Mirror
+
+    yi = only(Virtual.target(vi))
+
+    g = MOI.ScalarAffineFunction{T}(
+        SAT{T}[f.terms[i].scalar_term for i = 2:length(f.terms)],
+        sum(f.constants[i] for i = 2:length(f.constants)),
+    )
+
+    # Tell the compiler that quadratization is necessary
+    MOI.set(model, Attributes.Quadratize(), true)
+
+    if A === MOI.ACTIVATE_ON_ONE
+        return PBO.PBF{VI,T}(yi) * constraint(model, ci, g, s.set, arch)
+    elseif A === MOI.ACTIVATE_ON_ZERO
+        return (one(T) - PBO.PBF{VI,T}(yi)) * constraint(model, ci, g, s.set, arch)
+    else
+        error("Indicator constraint activation type $(A) not supported")
+    end
+
+    return nothing
+end
+
+function constraint(
+    model::Virtual.Model{T},
+    ci::CI,
+    f::MOI.VectorQuadraticFunction{T},
+    s::MOI.Indicator{A,S},
+    arch::AbstractArchitecture,
+) where {T,A,S}
+    # Indicator Constraint: y = 0|1 => {g(x)}
+
+    xi = first(f.affine_terms).scalar_term.variable # Indicator Variable
+    vi = model.source[xi]
+
+    @assert Virtual.encoding(vi) isa Encoding.Mirror
+
+    yi = only(Virtual.target(vi))
+
+    g = MOI.ScalarQuadraticFunction{T}(
+        SQT{T}[f.quadratic_terms[i].scalar_term for i = 2:length(f.quadratic_terms)],
+        SAT{T}[f.affine_terms[i].scalar_term for i = 2:length(f.affine_terms)],
+        sum(f.constants[i] for i = 2:length(f.constants)),
+    )
+
+    # Tell the compiler that quadratization is necessary
+    MOI.set(model, Attributes.Quadratize(), true)
+
+    if A === MOI.ACTIVATE_ON_ONE
+        return PBO.PBF{VI,T}(yi) * constraint(model, ci, g, s.set, arch)
+    elseif A === MOI.ACTIVATE_ON_ZERO
+        return (one(T) - PBO.PBF{VI,T}(yi)) * constraint(model, ci, g, s.set, arch)
+    else
+        error("Indicator constraint activation type $(A) not supported")
+    end
+
+    return nothing
+end

src/compiler/variables.jl

@@ -98,15 +98,19 @@ function variable_𝔹!(model::Virtual.Model{T}, i::Union{VI,CI}) where {T}
     return Encoding.encode!(model, i, Encoding.Mirror{T}())
 end
 
-function variable_ℤ!(model::Virtual.Model{T}, vi::VI, (a, b)::Tuple{T,T}) where {T}
-    if isnothing(a) || isnothing(b)
-        error("Unbounded variable $(vi) ∈ ℤ")
-    else
+function variable_ℤ!(model::Virtual.Model{T}, vi::VI, (a, b)::Tuple{A,B}) where {T,A<:Union{T,Nothing},B<:Union{T,Nothing}}
+    if !isnothing(a) && !isnothing(b)
         let e = Attributes.variable_encoding_method(model, vi)
             S = (a, b)
 
             return Encoding.encode!(model, vi, e, S)
         end
+    elseif !isnothing(b)
+        error("Unbounded variable $(vi) ∈ (-∞, $(b)] ⊂ ℤ ")
+    elseif !isnothing(a)
+        error("Unbounded variable $(vi) ∈ [$(a), +∞) ⊂ ℤ")
+    else
+        error("Unbounded variable $(vi) ∈ ℤ")
     end
 end
 
@@ -122,10 +126,8 @@ function variable_ℤ!(model::Virtual.Model{T}, ci::CI, (a, b)::Tuple{T,T}) wher
     end
 end
 
-function variable_ℝ!(model::Virtual.Model{T}, vi::VI, (a, b)::Tuple{T,T}) where {T}
-    if isnothing(a) || isnothing(b)
-        error("Unbounded variable $(vi) ∈ ℝ")
-    else
+function variable_ℝ!(model::Virtual.Model{T}, vi::VI, (a, b)::Tuple{A,B}) where {T,A<:Union{T,Nothing},B<:Union{T,Nothing}}
+    if !isnothing(a) && !isnothing(b)
         # TODO: Solve this bit-guessing magic??? (DONE)
         # IDEA: 
         #     Let x̂ ~ U[a, b], K = 2ᴺ, γ = [a, b]
@@ -151,6 +153,12 @@ function variable_ℝ!(model::Virtual.Model{T}, vi::VI, (a, b)::Tuple{T,T}) wher
                 return Encoding.encode!(model, vi, e, S; tol)
             end
         end
+    elseif !isnothing(b)
+        error("Unbounded variable $(vi) ∈ (-∞, $(b)]")
+    elseif !isnothing(a)
+        error("Unbounded variable $(vi) ∈ [$(a), +∞)")
+    else
+        error("Unbounded variable $(vi) ∈ ℝ")
     end
 end
 

src/encoding/variables/interval/arithmetic.jl

@@ -31,28 +31,25 @@ Arithmetic() = Arithmetic{Float64}()
     encode(var::Function, ::Arithmetic{T}, S::Tuple{T,T}) where {T}
 """
 function encode(var::Function, e::Arithmetic{T}, S::Tuple{T,T}; tol::Union{T,Nothing} = nothing) where {T}
-    isnothing(tol) || return encode(var, e, S, nothing; tol)
+    !isnothing(tol) && return encode(var, e, S, nothing; tol)
 
     a, b = integer_interval(S)
 
-    @assert b > a
+    if a == b
+        return (VI[], PBO.PBF{VI,T}(a), nothing)
+    end
 
     M = trunc(Int, b - a)
     N = ceil(Int, (sqrt(1 + 8M) - 1) / 2)
 
-    if N == 0
-        y = VI[]
-        ξ = PBO.PBF{VI,T}((a + b) / 2)
-    else
-        y = var(N)::Vector{VI}
-        ξ = PBO.PBF{VI,T}(
-            [
-                a
-                [y[i] => i for i = 1:N-1]
-                y[N] => M - N * (N - 1) / 2
-            ],
-        )
-    end
+    y = var(N)::Vector{VI}
+    ξ = PBO.PBF{VI,T}(
+        [
+            a
+            [y[i] => i for i = 1:N-1]
+            y[N] => M - N * (N - 1) / 2
+        ],
+    )
 
     return (y, ξ, nothing) # No penalty function
 end

src/encoding/variables/interval/binary.jl

@@ -33,7 +33,7 @@ Binary() = Binary{Float64}()
 
 # Integer
 function encode(var::Function, e::Binary{T}, S::Tuple{T,T}; tol::Union{T,Nothing} = nothing) where {T}
-    isnothing(tol) || return encode(var, e, S, nothing; tol)
+    !isnothing(tol) && return encode(var, e, S, nothing; tol)
 
     a, b = integer_interval(S)
 
@@ -47,19 +47,14 @@ function encode(var::Function, e::Binary{T}, S::Tuple{T,T}; tol::Union{T,Nothing
     M = trunc(Int, b - a)
     N = ceil(Int, log2(M + 1))
 
-    if N == 0
-        y = VI[]
-        ξ = PBO.PBF{VI,T}((a + b) / 2)
-    else
-        y = var(N)::Vector{VI}
-        ξ = PBO.PBF{VI,T}(
-            [
-                a
-                [y[i] => 2^(i - 1) for i = 1:N-1]
-                y[N] => M - 2^(N - 1) + 1
-            ],
-        )
-    end
+    y = var(N)::Vector{VI}
+    ξ = PBO.PBF{VI,T}(
+        [
+            a
+            [y[i] => 2^(i - 1) for i = 1:N-1]
+            y[N] => M - 2^(N - 1) + 1
+        ],
+    )
 
     return (y, ξ, nothing)
 end

src/encoding/variables/interval/bounded.jl

@@ -37,10 +37,14 @@ function encode(
     S::Tuple{T,T};
     tol::Union{T,Nothing} = nothing,
 ) where {E<:IntervalVariableEncodingMethod,T}
-    isnothing(tol) || return encode(var, e, S, nothing; tol)
+    !isnothing(tol) && return encode(var, e, S, nothing; tol)
 
     a, b = integer_interval(S)
 
+    if a == b
+        return (VI[], PBO.PBF{VI,T}(a), nothing)
+    end
+
     m = floor(Int, e.μ)
     n = trunc(Int, b - a)
     r = n ÷ m

src/encoding/variables/interval/unary.jl

@@ -39,20 +39,18 @@ Given ``S = [a, b] \subset \mathbb{Z}``, ``a < b``, let ``n = b - a`` and ``\mat
 ```
 """
 function encode(var::Function, e::Unary{T}, S::Tuple{T,T}; tol::Union{T,Nothing} = nothing) where {T}
-    isnothing(tol) || return encode(var, e, S, nothing; tol)
+    !isnothing(tol) && return encode(var, e, S, nothing; tol)
 
     a, b = integer_interval(S)
 
-    @assert b > a
+    if a == b
+        return (VI[], PBO.PBF{VI,T}(a), nothing)
+    end
 
     N = trunc(Int, b - a)
 
     y = var(N)::Vector{VI}
-    ξ = if N == 0
-        PBO.PBF{VI,T}((a + b) / 2)
-    else
-        PBO.PBF{VI,T}([a; [y[i] => one(T) for i = 1:N]])
-    end
+    ξ = PBO.PBF{VI,T}([a; [y[i] => one(T) for i = 1:N]])
 
     return (y, ξ, nothing) # No penalty function
 end

src/model/prequbo.jl

@@ -1,13 +1,34 @@
+const INDICATOR_EQ_ONE{T}        = MOI.Indicator{MOI.ACTIVATE_ON_ONE,EQ{T}}
+const INDICATOR_LT_ONE{T}        = MOI.Indicator{MOI.ACTIVATE_ON_ONE,LT{T}}
+const INDICATOR_GT_ONE{T}        = MOI.Indicator{MOI.ACTIVATE_ON_ONE,GT{T}}
+const INDICATOR_Interval_ONE{T}  = MOI.Indicator{MOI.ACTIVATE_ON_ONE,MOI.Interval{T}}
+const INDICATOR_EQ_ZERO{T}       = MOI.Indicator{MOI.ACTIVATE_ON_ZERO,EQ{T}}
+const INDICATOR_LT_ZERO{T}       = MOI.Indicator{MOI.ACTIVATE_ON_ZERO,LT{T}}
+const INDICATOR_GT_ZERO{T}       = MOI.Indicator{MOI.ACTIVATE_ON_ZERO,GT{T}}
+const INDICATOR_Interval_ZERO{T} = MOI.Indicator{MOI.ACTIVATE_ON_ZERO,MOI.Interval{T}}
+
 MOIU.@model(
-    PreQUBOModel,       # Name of model
-    (MOI.Integer, MOI.ZeroOne), # untyped scalar sets
-    (EQ, LT, GT),               #   typed scalar sets
-    (),                         # untyped vector sets
-    (MOI.SOS1,),                #   typed vector sets
-    (VI,),                      # untyped scalar functions
-    (SAF, SQF),                 #   typed scalar functions
-    (MOI.VectorOfVariables,),   # untyped vector functions
-    (),                         #   typed vector functions
-    false,                      # is optimizer?
+    PreQUBOModel,                # Name of model
+    (MOI.Integer, MOI.ZeroOne),  # untyped scalar sets
+    (EQ, LT, GT),                #   typed scalar sets
+    (),                          # untyped vector sets
+    (                            #   typed vector sets
+        MOI.SOS1,
+        INDICATOR_EQ_ONE,
+        INDICATOR_LT_ONE,
+        INDICATOR_GT_ONE,
+        INDICATOR_Interval_ONE,
+        INDICATOR_EQ_ZERO,
+        INDICATOR_LT_ZERO,
+        INDICATOR_GT_ZERO,
+        INDICATOR_Interval_ZERO,
+    ),
+    (VI,),                       # untyped scalar functions
+    (SAF, SQF),                  #   typed scalar functions
+    (MOI.VectorOfVariables,),    # untyped vector functions
+    (                            #   typed vector functions
+        MOI.VectorAffineFunction,
+        MOI.VectorQuadraticFunction,
+    ),
+    false,                       # is optimizer?
 )
-

src/wrapper.jl

@@ -127,38 +127,25 @@ function MOI.copy_to(model::Optimizer{T}, source::MOI.ModelLike) where {T}
 end
 
 # Objective Function Support
-MOI.supports(
-    ::Optimizer{T},
-    ::MOI.ObjectiveFunction{<:Union{VI,SAF{T},SQF{T}}},
-) where {T} = true
+function MOI.supports(model::Optimizer, f::MOI.ObjectiveFunction{F}) where {F}
+    return MOI.supports(model.source_model, f)
+end
 
 # Constraint Support
-MOI.supports_constraint(
-    ::Optimizer{T},
-    ::Type{VI},
-    ::Type{
-        <:Union{MOI.ZeroOne,MOI.Integer,MOI.Interval{T},MOI.LessThan{T},MOI.GreaterThan{T}},
-    },
-) where {T} = true
-
-MOI.supports_constraint(
-    ::Optimizer{T},
-    ::Type{<:Union{SAF{T},SQF{T}}},
-    ::Type{<:Union{MOI.EqualTo{T},MOI.LessThan{T},MOI.GreaterThan{T}}},
-) where {T} = true
-
-MOI.supports_constraint(
-    ::Optimizer{T},
-    ::Type{<:MOI.VectorOfVariables},
-    ::Type{<:MOI.SOS1},
-) where {T} = true
-
-MOI.supports_add_constrained_variable(
-    ::Optimizer{T},
-    ::Type{
-        <:Union{MOI.ZeroOne,MOI.Integer,MOI.Interval{T},MOI.LessThan{T},MOI.GreaterThan{T}},
-    },
-) where {T} = true
+function MOI.supports_constraint(
+    model::Optimizer,
+    ::Type{F},
+    ::Type{S},
+) where {F<:MOI.AbstractFunction,S<:MOI.AbstractSet}
+    return MOI.supports_constraint(model.source_model, F, S)
+end
+
+function MOI.supports_add_constrained_variable(
+    model::Optimizer,
+    ::Type{S},
+) where {S<:MOI.AbstractScalarSet}
+    return MOI.supports_add_constrained_variable(model.source_model, S)
+end
 
 function Base.show(io::IO, model::Optimizer)
     print(