Feature VectorAffineFuncton in AllDifferent/Table (#268)

* allow VAF in AllDifferenSet/TableSet

Author: Wikunia

CHANGELOG.md

@@ -4,7 +4,11 @@
 - Allow variables as constraint like `a || !b` instead of `a == 1 || b == 0`. [PR #267](https://github.com/Wikunia/ConstraintSolver.jl/pull/267)
   - **Attention** Does not check if variable is a binary variable
 - Support for indicator/reified in indicator/reified (without bridges) [PR #251](https://github.com/Wikunia/ConstraintSolver.jl/pull/251)
-  
+- Support for VectorAffineFunction in TableSet/AllDifferentSet 
+  - i.e `[x[i]+i for i in 1:n] in CS.AllDifferentSet()` 
+  - `[x,y,10] in CS.TableSet(...)`
+  - see [issue #235](https://github.com/Wikunia/ConstraintSolver.jl/issues/235) for in-depth examples  
+
 ## v0.6.9 (17th of July 2021)
 - set activator to false when inner violated [PR #266](https://github.com/Wikunia/ConstraintSolver.jl/pull/266)
 

docs/src/supported.md

@@ -82,7 +82,9 @@ The following list shows constraints that are implemented and those which are pl
   - [X] `>`
 - [X] All different
   - `@constraint(m, [x,y,z] in CS.AllDifferentSet())`
+  - `@constraint(m, [x,y+2,x+y] in CS.AllDifferentSet())`
 - [X] `TableSet` constraint [#130](https://github.com/Wikunia/ConstraintSolver.jl/pull/130)
+  - also with `VectorAffineFunction` i.e `[x,y+2,x+y] in CS.TableSet(...)`
 - Indicator constraints [#167](https://github.com/Wikunia/ConstraintSolver.jl/pull/167)
   - i.e `@constraint(m, b => {x + y >= 12})`
   - [X] for affine inner constraints
@@ -91,6 +93,8 @@ The following list shows constraints that are implemented and those which are pl
   - [X] Allow `||` inside the inner constraint i.e `@constraint(m, b => {x + y >= 12 || 2x + y <= 7})`
   - [-] Have `Indicator` and `Reified` inside one another
     - this is only supported for simple cases i.e bridge support is missing
+  - [ ] `VectorAffineFunction` in `AllDifferentSet` or `TableSet` 
+    - Please open an issue if needed for a problem. Probably not too hard to add
 - Reified constraints [#171](https://github.com/Wikunia/ConstraintSolver.jl/pull/171)
   - i.e `@constraint(m, b := {x + y >= 12})`
   - [X] for everything that is supported by indicator constraints

src/MOI_wrapper/constraints.jl

@@ -42,12 +42,24 @@ MOI.supports_constraint(
     ::Type{AllDifferentSetInternal},
 ) = true
 
+MOI.supports_constraint(
+    ::Optimizer,
+    ::Type{VAF{T}},
+    ::Type{AllDifferentSetInternal},
+) where {T <: Real} = true
+
 MOI.supports_constraint(
     ::Optimizer,
     ::Type{MOI.VectorOfVariables},
     ::Type{TableSetInternal},
 ) = true
 
+MOI.supports_constraint(
+    ::Optimizer,
+    ::Type{VAF{T}},
+    ::Type{TableSetInternal},
+) where {T <: Real} = true
+
 MOI.supports_constraint(
     ::Optimizer,
     ::Type{SAF{T}},
@@ -80,13 +92,17 @@ function MOI.supports_constraint(
     return A == MOI.ACTIVATE_ON_ONE || A == MOI.ACTIVATE_ON_ZERO
 end
 
+supports_inner_constraint(optimizer::Optimizer, func, set) = MOI.supports_constraint(optimizer, func, set)
+supports_inner_constraint(optimizer::Optimizer, func::Type{VAF{T}}, ::Type{AllDifferentSetInternal}) where T = false
+supports_inner_constraint(optimizer::Optimizer, func::Type{VAF{T}}, ::Type{TableSetInternal}) where T = false
+
 function MOI.supports_constraint(
     optimizer::Optimizer,
     func::Type{<:MOI.AbstractFunction},
     set::Type{OS},
 ) where {A,F,IS,OS<:CS.IndicatorSet{A,F,IS}}
     !(A == MOI.ACTIVATE_ON_ONE || A == MOI.ACTIVATE_ON_ZERO) && return false
-    return MOI.supports_constraint(optimizer, F, IS)
+    return supports_inner_constraint(optimizer, F, IS)
 end
 
 function MOI.supports_constraint(
@@ -95,7 +111,7 @@ function MOI.supports_constraint(
     set::Type{OS}
 ) where {A,F,IS,OS<:CS.ReifiedSet{A,F,IS}}
     !(A == MOI.ACTIVATE_ON_ONE || A == MOI.ACTIVATE_ON_ZERO) && return false
-    return MOI.supports_constraint(optimizer, F, IS)
+    return supports_inner_constraint(optimizer, F, IS)
 end
 
 function MOI.supports_constraint(
@@ -290,6 +306,47 @@ function MOI.add_constraint(
     return MOI.ConstraintIndex{MOI.VectorOfVariables,typeof(set)}(length(com.constraints))
 end
 
+"""
+    function MOI.add_constraint(
+        model::Optimizer,
+        vaf::VAF{T},
+        set::Union{AllDifferentSetInternal, TableSetInternal},
+    ) where T
+
+`VectorAffineFunction` constraint for `AllDifferentSet` or `TableSet` to support for things like 
+`[a+1, b+2, c+3] in CS.AllDifferentSet()` or
+`[a, b, 4] in TableSet(...)`
+"""
+function MOI.add_constraint(
+    model::Optimizer,
+    vaf::VAF{T},
+    set::Union{AllDifferentSetInternal, TableSetInternal},
+) where T
+    fs = MOIU.eachscalar(vaf)
+    variables = Vector{MOI.VariableIndex}(undef, length(fs))
+    for (i,f) in enumerate(fs)
+        # we need to create a new variable and SAF constraint when it's not a SVF
+        if !is_svf(f)
+            discrete, non_continuous_value = is_discrete_saf(f)
+            if !discrete
+                throw(DomainError(non_continuous_value, "The constant and all coefficients need to be discrete"))
+            end
+            vidx = MOI.add_variable(model)
+            variables[i] = vidx
+            min_val, max_val = get_extrema(model, f)
+            svf = MOI.SingleVariable(vidx)
+            MOI.add_constraint(model, svf, MOI.Integer())
+            MOI.add_constraint(model, svf, MOI.Interval(min_val, max_val))
+            new_constraint_fct = MOIU.operate(-, T, f, svf)
+            MOI.add_constraint(model, new_constraint_fct, MOI.EqualTo(0.0))
+        else
+            variables[i] = f.terms[1].variable_index
+        end
+    end
+    ci = MOI.add_constraint(model, MOI.VectorOfVariables(variables), set)
+    return MOI.ConstraintIndex{VAF{T},typeof(set)}(ci.value) 
+end
+
 """
     MOI.add_constraint(
         model::Optimizer,
@@ -326,10 +383,10 @@ function MOI.add_constraint(
 
     if length(func.terms) == 1
         vidx = func.terms[1].variable_index.value
-        val = convert(Int, set.value / func.terms[1].coefficient)
+        val = convert(Int, (set.value - func.constant) / func.terms[1].coefficient)
         push!(model.inner.init_fixes, (vidx, val))
         return MOI.ConstraintIndex{SAF{T},MOI.EqualTo{T}}(0)
-    elseif length(func.terms) == 2 && set.value == zero(T)
+    elseif length(func.terms) == 2 && set.value == zero(T) && func.constant == zero(T)
         if func.terms[1].coefficient == -func.terms[2].coefficient
             # we have the form a == b
             vecOfvar = MOI.VectorOfVariables([

src/MOI_wrapper/util.jl

@@ -98,4 +98,41 @@ end
 
 function get_activator_internals(A, indices)
     ActivatorConstraintInternals(A, indices[1] in indices[2:end], false, 0)
+end
+
+"""
+    saf_is_svf(saf::MOI.ScalarAffineFunction)
+
+Checks if a `ScalarAffineFunction` can be represented as a `SingleVariable`.
+This can be used for example when having a `VectorAffineFunction` in `AllDifferentSet` constraint when the decision 
+has to be made whether a new constraint + variable has to be created.
+"""
+function is_svf(saf::MOI.ScalarAffineFunction)
+    !iszero(saf.constant) && return false
+    length(saf.terms) != 1 && return false
+    !isone(saf.terms[1].coefficient) && return false
+    return true
+end
+
+function get_extrema(model::Optimizer, saf::MOI.ScalarAffineFunction{T}) where T
+    min_val = saf.constant
+    max_val = saf.constant
+    for term in saf.terms 
+        if term.coefficient < 0
+            min_val += term.coefficient*model.variable_info[term.variable_index.value].upper_bound
+            max_val += term.coefficient*model.variable_info[term.variable_index.value].lower_bound
+        else
+            min_val += term.coefficient*model.variable_info[term.variable_index.value].lower_bound
+            max_val += term.coefficient*model.variable_info[term.variable_index.value].upper_bound
+        end
+    end
+    return min_val, max_val
+end
+
+function is_discrete_saf(saf::MOI.ScalarAffineFunction{T}) where T 
+    !isapprox(saf.constant, round(saf.constant)) && return false, saf.constant
+    for term in saf.terms 
+        !isapprox(term.coefficient, round(term.coefficient)) && return false, term.coefficient
+    end
+    return true, 0
 end

test/moi.jl

@@ -52,6 +52,11 @@
         indicator_set = MOI.IndicatorSet{MOI.ACTIVATE_ON_ONE}(MOI.EqualTo(9.0))
         @test MOI.supports_constraint(optimizer, typeof(f), typeof(indicator_set))
 
+        indicator_set = CS.IndicatorSet{MOI.ACTIVATE_ON_ONE, typeof(f)}(CS.AllDifferentSetInternal(2))
+        @test !MOI.supports_constraint(optimizer, typeof(f), typeof(indicator_set))
+        indicator_set = CS.IndicatorSet{MOI.ACTIVATE_ON_ZERO, typeof(f)}(CS.TableSetInternal(2, [1 2; ]))
+        @test !MOI.supports_constraint(optimizer, typeof(f), typeof(indicator_set))
+
         @test MOI.supports_constraint(
             optimizer,
             typeof(f),

test/unit/constraints/alldifferent.jl

@@ -199,3 +199,40 @@ end
     @test CS.fix!(com, variables[2], 2; check_feasibility = false)
     @test !CS.is_constraint_violated(com, constraint, constraint.fct, constraint.set)
 end
+
+@testset "all 8queens solutions" begin 
+    n = 8
+    model = Model(optimizer_with_attributes(CS.Optimizer, "all_optimal_solutions"=>true, "logging"=>[]))
+
+    @variable(model, 1 <= x[1:n] <= n, Int)
+    @constraint(model, x in CS.AllDifferentSet())
+    @constraint(model, [x[i] + i for i in 1:n] in CS.AllDifferentSet())
+    @constraint(model, [x[i] - i for i in 1:n] in CS.AllDifferentSet())
+
+    optimize!(model)
+
+    status = JuMP.termination_status(model)
+
+    @test status == MOI.OPTIMAL
+    num_sols = MOI.get(model, MOI.ResultCount())
+    @test num_sols == 92
+    for sol in 1:num_sols
+        x_val = convert.(Integer,JuMP.value.(x; result=sol))
+        @test allunique(x_val)
+        @test allunique([x_val[i] + i for i in 1:n])
+        @test allunique([x_val[i] - i for i in 1:n])
+    end
+end
+
+@testset "domainerror due to term not discrete " begin 
+    n = 8
+    model = Model(optimizer_with_attributes(CS.Optimizer, "all_optimal_solutions"=>true, "logging"=>[]))
+
+    @variable(model, 1 <= x[1:n] <= n, Int)
+    @constraint(model, x in CS.AllDifferentSet())
+    # 1.5 not allowed
+    @constraint(model, [1.5*x[i] + i for i in 1:n] in CS.AllDifferentSet())
+    @constraint(model, [x[i] - i for i in 1:n] in CS.AllDifferentSet())
+
+    @test_throws DomainError optimize!(model)
+end

test/unit/constraints/table.jl

@@ -162,3 +162,74 @@ end
     @test CS.fix!(com, variables[constraint.indices[2]], 2; check_feasibility = false)
     @test !CS.is_constraint_violated(com, constraint, constraint.fct, constraint.set)
 end
+
+
+@testset "all modulo" begin 
+    function modulo(m, x, y, z)
+        lbx = !(x isa Integer)  ? round(Int, JuMP.lower_bound(x)) : x
+        ubx = !(x isa Integer) ? round(Int, JuMP.upper_bound(x)) : x
+        lby = !(y isa Integer) ? round(Int, JuMP.lower_bound(y)) : y
+        uby = !(y isa Integer) ? round(Int, JuMP.upper_bound(y)) : y
+    
+        table =  transpose(reduce(hcat,[ [i,j,i % j] for i in lbx:ubx, j in lby:uby if j != 0]))
+        @constraint(m, [x, y, z] in CS.TableSet(table))
+    end   
+
+    m = Model(optimizer_with_attributes(CS.Optimizer, "logging" => [], "all_solutions" => true))
+    @variable(m, 1 <= x[1:2] <= 10, Int)
+    modulo(m,x[1],2,x[2])
+    optimize!(m)
+
+    @test JuMP.termination_status(m) == MOI.OPTIMAL
+    nresults = JuMP.result_count(m)
+    results = Set{Tuple{Int,Int}}()
+    for i in 1:nresults
+        xv = convert.(Int, round.(JuMP.value.(x; result=i)))
+        @test xv[1] % 2 == xv[2] 
+        push!(results, (xv[1], xv[2]))
+    end
+
+    found_nr = 0
+    for i in 1:10, j in 1:10
+        if i % 2 == j 
+            @test (i,j) in results
+            found_nr += 1
+        end
+    end
+    @test found_nr == nresults
+end
+
+@testset "vector affine" begin 
+    function modulo_affine(m, x, y, z)
+        lbx = !(x isa Integer)  ? round(Int, JuMP.lower_bound(x)) : x
+        ubx = !(x isa Integer) ? round(Int, JuMP.upper_bound(x)) : x
+        lby = !(y isa Integer) ? round(Int, JuMP.lower_bound(y)) : y
+        uby = !(y isa Integer) ? round(Int, JuMP.upper_bound(y)) : y
+    
+        table =  transpose(reduce(hcat,[ [i,j,i % j] for i in lbx:2*ubx, j in lby:uby if j != 0]))
+        @constraint(m, [x+z, y, 2x-z] in CS.TableSet(table))
+    end   
+
+    m = Model(optimizer_with_attributes(CS.Optimizer, "logging" => [], "all_solutions" => true))
+    @variable(m, 1 <= x[1:2] <= 10, Int)
+    modulo_affine(m,x[1],2,x[2])
+    optimize!(m)
+
+    @test JuMP.termination_status(m) == MOI.OPTIMAL
+    nresults = JuMP.result_count(m)
+    results = Set{Tuple{Int,Int}}()
+    for i in 1:nresults
+        xv = convert.(Int, round.(JuMP.value.(x; result=i)))
+        @test (xv[1]+xv[2]) % 2 == 2*xv[1]-xv[2] 
+        push!(results, (xv[1], xv[2]))
+    end
+
+    found_nr = 0
+    for i in 1:10, j in 1:10
+        if (i+j) % 2 == 2*i-j 
+            @test (i,j) in results
+            found_nr += 1
+        end
+    end
+    @test found_nr == nresults
+end