ConstraintSolver.jl
ConstraintSolver.jl copied to clipboard
Is there an element constraint?
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()
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.
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 :smile:
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.)
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
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 :slightly_smiling_face:
Great progress, Ole!