Question/request about indicator and reification

[hakank]

I'm trying to implement the all_different_except_0 constraint using indicators, but am a little confused how indicators (and reifications) work. (See last for the main test function I used.)

The general idea of the (decomoposition of the) constraint is:

• for each pair x[i] and x[j (where i > j): if both x[i] != 0 and x[j] != 0 then x[i] != x[j]

The ideal way (for me) is the following, and it would be really great if that worked:

function all_different_except_0(m, a)
    len = length(a)
    for i in 2:len, j in 1:i-1
        @constraint(m, (a[i] != 0 && a[j] != 0) => {a[i] != a[j]} ) 
    end
end

However, it seems that the LHS of the indicator/reification must be a binary variable (defined with @variable).

Here's a second attempt where two binary variables are defined for the two parts of the @constraint:

function all_different_except_0(m, a)
    len = length(a)
    for i in 2:len, j in 1:i-1
            @variable(m, b1, Bin)
            @constraint(m, b1 := {(a[i] != 0 && a[j] != 0)} ) # Line 64
            @variable(m, b2, Bin)
            @constraint(m, b2 := {(a[i] != a[j])} )
            @constraint(m, b1 => {b2})
    end
end

But this don't work either inequality (!=) is not supported:

At /home/hakank/julia/constraints/all_different_except_0.jl:64: `@constraint(m, $(Expr(:(:=), :b1, :({a[i] != 0 && a[j] != 0}))))`: Constraints must be in one of the following forms:
       expr1 <= expr2
       expr1 >= expr2
       expr1 == expr2
       lb <= expr <= ub

Is it correct that '!=is not supported for indicator/reification constraints? The documentation https://wikunia.github.io/ConstraintSolver.jl/stable/supported/#Supported-constraints-1 seems to indicate that!=` is supported...

I then tried the following - using Boolean algebra - which seems to be somewhat better but - still - get some error.

function all_different_except_0(m, a)
    len = length(a)
    for i in 2:length(a), j in 1:i-1
            @constraint(m, (a[i] != c) * (a[j] != c) <= (a[i] != a[j]) )
  end
end

But this throws the following error in the optimize! function (so the parsing seems to be correct):

ERROR: LoadError: BoundsError: attempt to access 0-element Array{ConstraintSolver.Variable,1} at index [1]
Stacktrace:
 [1] getindex at ./array.jl:809 [inlined]
 [2] set_pvals!(::ConstraintSolver.ConstraintSolverModel{Float64}, ::ConstraintSolver.LinearConstraint{Float64}) at /home/hakank/.julia/packages/ConstraintSolver/yg7vZ/src/ConstraintSolver.jl:95
 [3] set_pvals!(::ConstraintSolver.Optimizer) at /home/hakank/.julia/packages/ConstraintSolver/yg7vZ/src/MOI_wrapper/constraints.jl:517
 [4] optimize!(::ConstraintSolver.Optimizer) at /home/hakank/.julia/packages/ConstraintSolver/yg7vZ/src/MOI_wrapper/MOI_wrapper.jl:121
 [5] optimize!(::MathOptInterface.Bridges.LazyBridgeOptimizer{ConstraintSolver.Optimizer}) at /home/hakank/.julia/packages/MathOptInterface/ZJFKw/src/Bridges/bridge_optimizer.jl:264
 [6] optimize!(::MathOptInterface.Utilities.CachingOptimizer{MathOptInterface.AbstractOptimizer,MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}}}) at /home/hakank/.julia/packages/MathOptInterface/ZJFKw/src/Utilities/cachingoptimizer.jl:215
 [7] optimize!(::Model, ::Nothing; bridge_constraints::Bool, ignore_optimize_hook::Bool, kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at /home/hakank/.julia/packages/JuMP/e0Uc2/src/optimizer_interface.jl:130
 [8] optimize! at /home/hakank/.julia/packages/JuMP/e0Uc2/src/optimizer_interface.jl:106 [inlined] (repeats 2 times)
 [9] all_different_except_0_test(::Int64) at /home/hakank/julia/constraints/all_different_except_0.jl:86
 [10] all_different_except_0_test() at /home/hakank/julia/constraints/all_different_except_0.jl:72
 [11] top-level scope at /home/hakank/julia/constraints/all_different_except_0.jl:109
 [12] include(::String) at ./client.jl:457
 [13] top-level scope at ./timing.jl:174 [inlined]
 [14] top-level scope at ./none:0
in expression starting at /home/hakank/julia/constraints/all_different_except_0.jl:109

I then tested the following just to try the general idea of using binary variables in a loop:

function all_different_except_0(m, a)
    len = length(a)
    for i in 2:length(a), j in 1:i-1
        @variable(m, b1, Bin)
        @variable(m, b2, Bin)
        @constraint(m, b1+b2 == 1)
  end
end

Which throws:

LoadError: An object of name b1 is already attached to this model. If this is intended, consider using the anonymous construction syntax, e.g., x = @variable(model, [1:N], ...) where the name of the object does not appear inside the macro.

So it seems that one have to prepare an array of binary variables before the loop and then use it in the loop. Here's a version that does that but it also throws an error.

function all_different_except_0(m, a)
    len = length(a)
    b_len = length([1 for i in 2:len for j in 1:i-1 for k in 1:2]) # calcuiate the number of binary variables
    @variable(m, bs[1:b_len], Bin)
    c = 1 # counter
    for i in 2:length(a), j in 1:i-1
        @constraint(m, (a[i] != 0 && a[j] != 0) => {bs[c] == 1})
        @constraint(m, bs[c] == 1 => {a[i] != a[j]})
        c += 2 # increment the counter
  end
end

This throws this rather uninformative error:

ERROR: LoadError: MethodError: no method matching _build_reified_constraint(::JuMP.var"#_error#73"{Symbol,LineNumberNode}, ::Bool, ::ScalarConstraint{GenericAffExpr{Float64,VariableRef},MathOptInterface.EqualTo{Float64}}, ::Type{ConstraintSolver.ReifiedSet{MathOptInterface.ACTIVATE_ON_ONE}})
Closest candidates are:
  _build_reified_constraint(::Function, ::AbstractVariableRef, ::ScalarConstraint, ::Type{ConstraintSolver.ReifiedSet{A}}) where A at /home/hakank/.julia/packages/ConstraintSolver/yg7vZ/src/MOI_wrapper/reified.jl:1
  _build_reified_constraint(::Function, ::AbstractVariableRef, ::VectorConstraint, ::Type{ConstraintSolver.ReifiedSet{A}}) where A at /home/hakank/.julia/packages/ConstraintSolver/yg7vZ/src/MOI_wrapper/reified.jl:9
Stacktrace:
 [1] macro expansion at /home/hakank/.julia/packages/JuMP/e0Uc2/src/macros.jl:441 [inlined]
 [2] all_different_except_c(::Model, ::Array{VariableRef,1}, ::Int64) at /home/hakank/julia/constraints/all_different_except_0.jl:61
 [3] all_different_except_0_test(::Int64) at /home/hakank/julia/constraints/all_different_except_0.jl:81
 [4] all_different_except_0_test() at /home/hakank/julia/constraints/all_different_except_0.jl:72
 [5] top-level scope at /home/hakank/julia/constraints/all_different_except_0.jl:109
 [6] include(::String) at ./client.jl:457
 [7] top-level scope at ./timing.jl:174 [inlined]
 [8] top-level scope at ./none:0
in expression starting at /home/hakank/julia/constraints/all_different_except_0.jl:109

I'm not sure if anyone of these versions should work, but - as mentioned earlier - I wonder how one work with indicator/reification in a function such as all_different_except_0.

Here's the test function I use, where the different all_different_except_0 functions can be inserted.

using ConstraintSolver, JuMP
using Cbc
const CS = ConstraintSolver

function all_different_except_0(m, a)
   # some of the functions above
end 

function all_different_except_0_test(n=4)
    println("n:$n")
    model = Model(optimizer_with_attributes(CS.Optimizer,
                                            "all_solutions"=>true,
                                            "logging"=>[],
                                            "time_limit"=>10,
                                            )
                                            )
    @variable(model, 0 <=  x[1:n] <= n, Int)

    all_different_except_0(model,x)

    println("solve")
    optimize!(model)

    status = JuMP.termination_status(model)
    println("status:$status")
    if status == MOI.OPTIMAL
        # println("x:$x")
        # x = convert.(Integer,JuMP.value.(x))
        # println(x)
        num_sols = MOI.get(model, MOI.ResultCount())
        println("\nnum_sols:$num_sols\n")

        for sol in 1:num_sols
            println("solution #$sol")
            xx = convert.(Integer,JuMP.value.(x,result=sol))
            println("$xx\n")

        end

    end

end

all_different_except_0_test()

``'

[Wikunia]

Thanks for the issue.

I'm not aware that things like (a[i] != 0 && a[j] != 0) work in JuMP in general. I definitely don't have && implemented.

For the != inside an indicator or reified constraint: This works in general so:

function all_different_except_0(m, a)
    len = length(a)
    for i in 2:len, j in 1:i-1
            @variable(m, b1, Bin)
            @constraint(m, b1 := {(a[i] != 0 && a[j] != 0)} ) 
            @variable(m, b2, Bin)
            @constraint(m, b2 := {(a[i] != a[j])} ) # Line 66
            @constraint(m, b1 => {b2})
    end
end

line 66 should work but the error is in line 64 where && fails.

Have you tried something like the following?

b1 := {a[i] != 0}
b2 := {a[j] != 0}
b3 := {b1 + b2 != 2}
b3 => {a[i] != a[j]}

And yes b{1,2,3} would need to be arrays of binary variables to distinguish their names.

[hakank]

Thanks for your comments, Ole.

Regarding your comment about &&:

@constraint(m, b1 := {(a[i] != 0 && a[j] != 0)} ) 

I tried this instead:

@constraint(m, bs[c] := {(a[i] >= 1) + (a[j] >= 1) == 2} )

But got the error:

LoadError: Unexpected comparison in expression `a[i] >= 1`.

So it seems that one cannot use these kind of conditions in the RHS.

Well, I'll tested further with your suggestion

b1 := {a[i] != 0}
b2 := {a[j] != 0}
b3 := {b1 + b2 != 2}  # <-- See comment 
b3 => {a[i] != a[j]}

However, I changed the line to b3 := {b1 + b2 != 2} to ... == 2 since it's only when both a[i] != 0 and a[j] != 0 that we should require that a[i] != a[j]. If one of them is 0 we'll ignore this constraint.

This is implemented as following. I had to move the all different except 0 constraint in the main function to be able to print the bs array, otherwise JuMP didn't recognize it.

function all_different_except_0_test(n=2)
    println("n:$n")
    model = Model(optimizer_with_attributes(CS.Optimizer,
                                            "all_solutions"=>true,
                                            "logging"=>[],
                                            "time_limit"=>10,
                                            )
                                            )
    @variable(model, 0 <=  x[1:n] <= n, Int)
    b_len = length([1 for i in 2:n for j in 1:i-1 for k in 1:3])
    println("b_len:$b_len")
    @variable(model, bs[1:b_len], Bin)
    c = 1
    for i in 2:n, j in 1:i-1
        # b1: bs[c]
        # b2: bs[c+1]
        # b3: bs[c+2]
        @constraint(model, bs[c]   := {x[i] != 0})
        @constraint(model, bs[c+1] := {x[j] != 0})
        @constraint(model, bs[c+2] := {bs[c] + bs[c+1] == 2}) # Note: I changed this to == 2
        @constraint(model, bs[c+2] => {x[i] != x[j]})
        c += 3
    end

    optimize!(model)

    status = JuMP.termination_status(model)
    println("status:$status")
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        println("\nnum_sols:$num_sols\n")

        for sol in 1:num_sols
            println("solution #$sol")
            xx = convert.(Integer,JuMP.value.(x,result=sol))
            bss = convert.(Integer,JuMP.value.(bs,result=sol))
            println("x:$xx")
            println("bs:$bss\n")

        end
    end
end

all_different_except_0_test()

It don't yield any error (which is nice!), but for n=2 the model show 27 solutions which is not correct (it should be these 7 solutions: [0,0],[0,1],[0,2],[1,0],[1,2],[2,0],[2,1]). Especially, the model yield incorrectly x = [1,1] and [2,2] multiple times. Can you spot some error in the model?

solution #1
x:[0, 0]
bs:[0, 0, 0]

solution #2
x:[0, 1]
bs:[0, 0, 0]

solution #3
x:[0, 2]
bs:[0, 0, 0]

solution #4
x:[1, 0]
bs:[0, 0, 0]

solution #5
x:[1, 1]
bs:[0, 0, 0]

solution #6
x:[1, 2]
bs:[0, 0, 0]

solution #7
x:[2, 0]
bs:[0, 0, 0]

solution #8
x:[2, 1]
bs:[0, 0, 0]

solution #9
x:[2, 2]
bs:[0, 0, 0]

solution #10
x:[1, 0]
bs:[0, 1, 0]

solution #11
x:[1, 1]
bs:[0, 1, 0]

solution #12
x:[1, 2]
bs:[0, 1, 0]

solution #13
x:[2, 0]
bs:[0, 1, 0]

solution #14
x:[2, 1]
bs:[0, 1, 0]

solution #15
x:[2, 2]
bs:[0, 1, 0]

solution #16
x:[0, 1]
bs:[1, 0, 0]

solution #17
x:[1, 1]
bs:[1, 0, 0]

solution #18
x:[2, 1]
bs:[1, 0, 0]

solution #19
x:[1, 1]
bs:[1, 1, 0]

solution #20
x:[2, 1]
bs:[1, 1, 0]

solution #21
x:[0, 2]
bs:[1, 0, 0]

solution #22
x:[1, 2]
bs:[1, 0, 0]

solution #23
x:[2, 2]
bs:[1, 0, 0]

solution #24
x:[1, 2]
bs:[1, 1, 0]

solution #25
x:[2, 2]
bs:[1, 1, 0]

solution #26
x:[1, 2]
bs:[1, 1, 1]

solution #27
x:[2, 1]
bs:[1, 1, 1]

[Wikunia]

Thanks yeah it should be == 2. Seems like there is an issue with the reified constraint for !=. Will have a look.

Here:

solution #2
x:[0, 1]
bs:[0, 0, 0]

one bs should be 1, right? As x[2] != 0.

[hakank]

Yes, bs should be [0,1,0}.

[Wikunia]

Should be solved in #204

If interested you should be able to test this by using: ] add ConstraintSolver#bugfix-202-reified

[hakank]

Wow, that was quick!

I tested the fix and it work as expected. Thanks! (And I'm impressed by the packages feature of Julia. It's really grown on me. :-))

Should we keep this issue still open so you remember the other comments/wishes?

An aside: Since I now understand indicators/reification better, I could implemented the Who killed Agatha problem without much problem: http://hakank.org/julia/constraints/who_killed_agatha.jl (It shows 8 occurrences of agatha since the problem is underconstrained). I've started to publish my ConstraintSolver.jl models at http://hakank.org/julia/constraints/

[Wikunia]

Great I'll create a new version soon then. Currently checking out #200 . I think it might make sense to create a new issue for

function all_different_except_0(m, a)
    len = length(a)
    for i in 2:len, j in 1:i-1
        @constraint(m, (a[i] != 0 && a[j] != 0) => {a[i] != a[j]} ) 
    end
end

and maybe some other parts to tackle them individually.

[Wikunia]

BTW just a small comment on: http://hakank.org/julia/constraints/to_num.jl

 # Why can't I use this first and later print all solutions?
        # This give MethodError: no method matching value(::Int64; result=1)
        #
        # x = convert.(Integer,JuMP.value.(x))
        # v = convert.(Integer,JuMP.value.(v))

You try to convert a vector of variables x into a vector of integers x by using the same name x you loose the ability to call JuMP.value(x; result=sol) afterwards.

[hakank]

BTW just a small comment on: http://hakank.org/julia/constraints/to_num.jl

 # Why can't I use this first and later print all solutions?
        # This give MethodError: no method matching value(::Int64; result=1)
        #
        # x = convert.(Integer,JuMP.value.(x))
        # v = convert.(Integer,JuMP.value.(v))

You try to convert a vector of variables x into a vector of integers x by using the same name x you loose the ability to call JuMP.value(x; result=sol) afterwards.

Thanks for this. Yes, it works if I rename the variable (as done in the loop some lines below):

        xx = convert.(Integer,JuMP.value.(x))
        vv = convert.(Integer,JuMP.value.(v))
``

I'll open another issue about the indicator/reification.