Is there an element constraint?

[hakank]

I tried the following model which tests the x[ix] construct where x is an a array of decisions variables and ix is a decision variable, The constraint @constraint(m, x[ix] == ix) yield the error ArgumentError: invalid index: ix of type VariableRef.

Is there some way of doing an element constraint, or a plan to implement this? (In traditional constraint programming the element constraint is a quite important constraint.)

Here's a simple model testing this:

using ConstraintSolver, JuMP
const CS = ConstraintSolver

function element_test()
    n = 4
    m = Model(optimizer_with_attributes(CS.Optimizer, "all_solutions"=>true,"logging"=>[]))
    @variable(m, 1 <=  x[1:n] <= n,Int)
    @variable(m, 1 <=  ix <= n,Int)
    @constraint(m, x[ix] == ix)   # Element constraint

    optimize!(m)
    status = JuMP.termination_status(m)
    println("status:$status")
    if status == MOI.OPTIMAL
        x = convert.(Integer,JuMP.value.(x))
        println(x)
    end

end

element_test()

[Wikunia]

Thanks for opening the issue. Unfortunately there isn't yet. For your example I could create own the next week which would look like this:

@constraint(m, x in ElementSet(ix, ix)

or something like that. Writing it like you is currently out of reach as one would need to rewrite JuMP by quite a bit I think.

My way is very ugly which is somewhat the reason why I haven't implemented it yet. Additionally it would limit it quite a bit such that you would be able to write x[ix] + y[iy] >= 9 or something.

If you have any idea of a better syntax please let me know.

Something like x in Set is very easy to add. There is also a way to do something like x := constraint or x => constraint which I use for the reified and indicator constraint. Maybe that can be used for somehow defining it more user friendly.

In the long run I'll try to add:

 @constraint(m, x[ix] == ix)

but that's currently out of reach.

[Wikunia]

Let me see. Maybe I can even support the nice syntax with some of this: https://github.com/dourouc05/JuCP.jl/blob/master/src/sets.jl

Will get back to you in a week 😄

[hakank]

It would be really great if you support x[ix], but in fact this syntax is rarely suported in CP system (MiniZinc supports it for example). I especially like if one could use x[ix] + y[iy] == z[iz].

The other CP systems has variants are something like element(ix,x,val) (for x{ix] = val), so your suggestion of @constraint(m, x in ElementSet(ix, val) would work well.

(Regarding reification/indicator I will get back with some questions/comments later.)

[hakank]

FYI: Here's an decomposition of element (as my_element) using indicator/reification:

#
# my_element(model, ix, a, val)
#
# Ensure that a[ix] = val
#
function my_element(model, ix, a, val)
    n = length(a)
    b = @variable(model, [1:n], Bin)
    for i in 1:n
        @constraint(model,b[i] => {a[i] == val})
        @constraint(model,b[i] := {ix == i})
    end
end

It works - for example with this Langford model. As other decompositions, it might not scale well...

using ConstraintSolver, JuMP
using Cbc, GLPK
const CS = ConstraintSolver

function langford(k=5,symmetry_breaking=true,all_solutions=true,print_solutions=true)

    if !((k % 4 == 0) || (k % 4 == 3))
        println("There is no solutions for k ($k) unless K mod 4 == 0 or K mod 4 == 3")
        return
    end

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)

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

                                                            "traverse_strategy"=>:BFS,
                                                            # "traverse_strategy"=>:DFS, # <-
                                                            # "traverse_strategy"=>:DBFS,

                                                            "branch_split"=>:Smallest,
                                                            # "branch_split"=>:Biggest,
                                                            # "branch_split"=>:InHalf, # <-

                                                            # "simplify"=>false,
                                                            # "simplify"=>true, # default

                                                            "time_limit"=>6,

                                                            # "lp_optimizer" => cbc_optimizer,
                                                            # "lp_optimizer" => glpk_optimizer,
                                        ))
    k2 = 2*k

    @variable(model, 1 <= position[1:k2] <= k2, Int)
    @variable(model, 1 <= solution[1:k2] <= k, Int)

    @constraint(model, position in CS.AllDifferentSet())

    # symmetry breaking
    if symmetry_breaking
        @constraint(model, solution[1] <= solution[k2])
    end

    for i in 1:k
        @constraint(model, position[k+i] == position[i] + i+1)
        my_element(model, position[i], solution, i)
        my_element(model, position[k+i], solution, i)
    end

    # Solve the problem
    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")
                positionx = convert.(Integer,JuMP.value.(position; result=sol))
                solutionx = convert.(Integer,JuMP.value.(solution; result=sol))
                println("solution:$solutionx position:$positionx\n")
            end
        end
    else
        println("status:$status")
    end

    return status
end

all_solutions = true
print_solutions = true
symmetry_breaking = false
@time langford(7,symmetry_breaking,all_solutions,print_solutions)

# Show the number of solutions for n in 1:12
symmetry_breaking = true
all_solutions = true
print_solutions = false
for n in 1:12
    if n % 4 == 3  || n % 4 == 0
        println("n:$n")
        @time status = langford(n,symmetry_breaking,all_solutions,print_solutions)
        if status != MOI.OPTIMAL
            break
        end
    end
end

[Wikunia]

The first step is done with #213 but it currently only works for constant arrays so your example isn't working yet.
At least the nice syntax of z == T[y] is working 🙂

[hakank]

Great progress, Ole!