Strange error (decomposition of matrix element)

[hakank]

The following model throws a strange/uninformative error when activating one of the two commented constraints in the matrix_element function. The error is "no method matching _build_indicator_constraint". See below for the full stack trace.

The principle of the constraint is to fix the row (via the bi binary array) and the column (via bj array) and selecting the row/column in the matrix x.

What I can see, this decomposition should work:

• I tested the same approach in Picat and it works
• The values of bi and bj reification arrays are correct
  • When ii = 2 then bi is [0,1,0
  • When jj = 3 then bj is [0,0,1]

However, when uncommenting one of the two "connection" constraints in matrix_element the model don't work, i.e.

• @constraint(model, val == x[k,l] => {bi[k] + bj[l] == 2} )
• @constraint(model, (bi[k] + bj[l] == 2) => { val == x[k,l] } )

(I'm testing both to debug and they throw the same error. Also changing to reification instead give the same error.)

Is there an error in this model or is it some kind of bug in handling indicators?

using ConstraintSolver, JuMP
using GLPK
const CS = ConstraintSolver

#
#  matrix_element(model,x,i,j,val)
# 
# Ensure that val = x{i,j]
# 
function matrix_element(model, x,i,j,val)
    n,_ = size(x)
    bi = @variable(model,[1:n], Bin) # row
    bj = @variable(model,[1:n], Bin) # column
    @constraint(model, sum(bi) == 1)
    @constraint(model, sum(bj) == 1)
    for k in 1:n, l in 1:n
        @constraint(model,bi[k] := { k == i }) # fix the row
        @constraint(model,bj[l] := { l == j })  # fix the column

        # Connection constraint
        # Both these constraints yield the error "no method matching _build_indicator_constraint"
        # @constraint(model, val == x[k,l]  => {bi[k] + bj[l] == 2} )
        # @constraint(model, (bi[k] + bj[l] == 2) => { val == x[k,l] } )
    end
    return bi,bj # return for debugging
end

function matrix_element_test(n=4,print_solutions=true,all_solutions=true)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,   "all_solutions"=> all_solutions,
                                                            "logging"=>[],

                                                            "traverse_strategy"=>:BFS,
                                                            "branch_split"=>:InHalf,
                                                            "time_limit"=>3,

                                                            "lp_optimizer" => glpk_optimizer,
                                        ))

    @variable(model, 1 <= x[1:n,1:n] <= n*n, Int)

    @constraint(model, x[:] in CS.AllDifferentSet())

    # Indices and value x[ii,jj] = val
    @variable(model, 1 <= ii <= n, Int)
    @variable(model, 1 <= jj <= n, Int)
    @variable(model, 1 <= val <= n*n, Int)

    # Fix some values:
    @constraint(model, ii == 2)
    @constraint(model, jj == 3)
    # @constraint(model, val == 6)

    bi, bj = matrix_element(model,x,ii,jj,val)

    optimize!(model)

    status = JuMP.termination_status(model)
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        println("num_sols:$num_sols\n")
        if print_solutions
            for sol in 1:num_sols
                println("solution #$sol")
                xx = convert.(Integer,JuMP.value.(x; result=sol))
                iix = convert.(Integer,JuMP.value.(ii; result=sol))
                jjx = convert.(Integer,JuMP.value.(jj; result=sol))
                valx = convert.(Integer,JuMP.value.(val; result=sol))

                bix = convert.(Integer,JuMP.value.(bi; result=sol))
                bjx = convert.(Integer,JuMP.value.(bj; result=sol))
                display(xx)
                println("ii:$iix jj:$jjx val:$valx")
                println("bi:$bix")
                println("bj:$bjx")
                println()
            end
        end
    else
        println("status:$status")
    end

    return status
end

print_solutions=true
all_solutions=false
@time matrix_element_test(3,print_solutions,all_solutions)

As it stand, the program prints:

n:3
num_sols:1

solution #1
3?3 Array{Int64,2}:
 9  6  3
 8  5  2
 7  4  1
ii:2 jj:3 val:9
bi:[0, 1, 0]
bj:[0, 0, 1]

Here's the full stacktrace when the first of the two connection constraints is uncommented

n:3
x:VariableRef[x[1,1] x[1,2] x[1,3]; x[2,1] x[2,2] x[2,3]; x[3,1] x[3,2] x[3,3]] i:ii j:jj val:val
k:1 l:1 x[k,l]: x[1,1]
ERROR: LoadError: MethodError: no method matching _build_indicator_constraint(::JuMP.var"#_error#73"{Symbol,LineNumberNode}, ::Bool, ::ScalarConstraint{GenericAffExpr{Float64,VariableRef},MathOptInterface.EqualTo{Float64}}, ::Type{MathOptInterface.IndicatorSet{MathOptInterface.ACTIVATE_ON_ONE,S} where S<:MathOptInterface.AbstractScalarSet})
Closest candidates are:
  _build_indicator_constraint(::Function, ::AbstractVariableRef, ::ScalarConstraint, ::Type{MathOptInterface.IndicatorSet{A,S} where S<:MathOptInterface.AbstractScalarSet}) where A at /home/hakank/.julia/packages/JuMP/e0Uc2/src/indicator.jl:9
  _build_indicator_constraint(::Function, ::AbstractVariableRef, ::VectorConstraint, ::Type{MathOptInterface.IndicatorSet{A,S} where S<:MathOptInterface.AbstractScalarSet}) where A at /home/hakank/.julia/packages/ConstraintSolver/AIwl4/src/MOI_wrapper/constraints.jl:9
Stacktrace:
 [1] macro expansion at /home/hakank/.julia/packages/JuMP/e0Uc2/src/macros.jl:441 [inlined]
 [2] matrix_element(::Model, ::Array{VariableRef,2}, ::VariableRef, ::VariableRef, ::VariableRef) at /home/hakank/julia/constraints/matrix_element.jl:31
 [3] matrix_element_test(::Int64, ::Bool, ::Bool) at /home/hakank/julia/constraints/matrix_element.jl:76
 [4] top-level scope at ./timing.jl:174
 [5] include(::String) at ./client.jl:457
 [6] top-level scope at ./timing.jl:174 [inlined]
 [7] top-level scope at ./none:0
in expression starting at /home/hakank/julia/constraints/matrix_element.jl:116

[Wikunia]

This is a bit similar to #205 . Currently the indicator or reified variable has to be a binary variable.

In your example:

@constraint(model, val == x[k,l]  => {bi[k] + bj[l] == 2} )

Needs to be written as:

@constraint(model, binary_variable := {val == x[k,l]})
@constraint(model, binary_variable  => {bi[k] + bj[l] == 2} )

I'll have a look at supporting this.

[hakank]

Ah, of course! I should have realized that.

Adding a new variable bij[1:n,1:n] seems to work. I have to test it further but it seems promising...

function matrix_element(model, x,i,j,val)
    n,_ = size(x)
    bi = @variable(model,[1:n], Bin)
    bj = @variable(model,[1:n], Bin)
    bij = @variable(model, [1:n,1:n], Bin) # Added
    @constraint(model, sum(bi) == 1)
    @constraint(model, sum(bj) == 1)
    @constraint(model, sum(bij[:]) == 1) # Added
    for k in 1:n, l in 1:n
        @constraint(model,bi[k] := { k == i })
        @constraint(model,bj[l] := { l == j })
        @constraint(model,bij[k,l] := { bi[k] + bj[l] == 2 })  # Added
        @constraint(model, bij[k,l] => { val == x[k,l] } )
    end
    return bi,bj,bij
end