Search code examples
julia-jump

Julia JuMP - how to use variable with CartesianIndex?

In order to reduce the size of my variable, I wish to define a JuMP model with a CartesianIndex.

I have tried the following code, but I keep getting problems dealing with the CartesianIndex when using the variable:

using JuMP, Cbc, IterTools

Random.seed!(3)
N = 10
cost = rand(1:4, N, N)

# reducing size using an index
idx = findall(cost .< 2)
# size reduce from 100 to 22 variables

model = Model(with_optimizer(Cbc.Optimizer));
@variable(model, x[idx], Bin);

# constraint to have only sum(x) = 1 over each row
for i in 1:N
  # creating a second CartesianIndex to select row i, then intersect with the 22 indices of x
  cart_idx = CartesianIndex.(collect(product(i,1:N)))
  cart_idx = intersect(cart_idx, idx)
  
  @constraint(model, sum(x[cart_idx]) == 1)
  # ERROR: KeyError: key CartesianIndex{2}[CartesianIndex(1, 5), CartesianIndex(1, 6)] not found
  
end

The last line of code in the loop throws the error:

ERROR: KeyError: key CartesianIndex{2}[CartesianIndex(1, 5), CartesianIndex(1, 6)] not found

Is it possible to use CartesianIndex with a JuMP variable, and how to fix the issue above?

Solution

  • The next version of JuMP will throw an error if you attempt to construct a DenseAxisArray using CartesianIndices. They are not intended for that purpose.

    You should use an alternative approach:

    Option 1

    using JuMP, IterTools
Random.seed!(3)
N = 10
cost = rand(1:4, N, N)
model = Model()
@variable(model, x[i=1:N, j=1:N; cost[i, j] < 2], Bin)
@constraint(
    model, 
    [ii=1:N], 
    sum(x[i, j] for (i, j) in eachindex(x) if ii == i) == 1,
)

    Option 2

    using JuMP, IterTools
Random.seed!(3)
N = 10
cost = rand(1:4, N, N)
S = [(i, j) for i in 1:N, j in 1:N if cost[i, j] < 2]
model = Model()
@variable(model, x[S], Bin)
@constraint(model, [ii=1:N],  sum(x[i, j] for (i, j) in S if ii == i) == 1)

    Prior post:

    I think there are a few bugs relating to cartesian indices in JuMP:

    julia> model = Model()
A JuMP Model
Feasibility problem with:
Variables: 0
Model mode: AUTOMATIC
CachingOptimizer state: NO_OPTIMIZER
Solver name: No optimizer attached.

julia> S = CartesianIndex.([2, 4])
2-element Vector{CartesianIndex{1}}:
 CartesianIndex(2,)
 CartesianIndex(4,)

julia> @variable(model, x[S])
1-dimensional DenseAxisArray{VariableRef,1,...} with index sets:
    Dimension 1, CartesianIndex{1}[CartesianIndex(2,), CartesianIndex(4,)]
And data, a 2-element Vector{VariableRef}:
 x[CartesianIndex(2,)]
 x[CartesianIndex(4,)]

julia> x[CartesianIndex(2,)]
x[CartesianIndex(4,)]

julia> x[CartesianIndex(1,)]
x[CartesianIndex(2,)]

julia> x[CartesianIndex(4,)]
ERROR: BoundsError: attempt to access 2-element Vector{VariableRef} at index [4]
Stacktrace:
 [1] getindex
   @ ./array.jl:801 [inlined]
 [2] getindex
   @ ./multidimensional.jl:637 [inlined]
 [3] getindex(A::JuMP.Containers.DenseAxisArray{VariableRef, 1, Tuple{Vector{CartesianIndex{1}}}, Tuple{JuMP.Containers._AxisLookup{Dict{CartesianIndex{1}, Int64}}}}, idx::CartesianIndex{1})
   @ JuMP.Containers ~/.julia/packages/JuMP/2IF9U/src/Containers/DenseAxisArray.jl:323
 [4] top-level scope
   @ REPL[42]:1

    It seems that if we have cartesian indices as keys, we should look them up via the key, not via the typical cartesian index.

    Edit: I opened an issue https://github.com/jump-dev/JuMP.jl/issues/2825