GitHub

Strategies for dealing with many indices

This tutorial was generated using Literate.jl. Download the source as a .jl file.

This tutorial demonstrates two techniques—keyword indexing and typed indices —for avoiding common index-ordering bugs in JuMP models with many sets. Such bugs are easy to introduce in models with variables like x[p, f, m, t] and hard to diagnose because the wrong ordering is syntactically valid.

Learning intentions:

  • Understand how index-ordering bugs arise in models with many sets, and why they are hard to detect when all orderings are syntactically valid
  • Use keyword indexing with DenseAxisArray to make the role of each index explicit and catch wrong-order access at runtime
  • Define Julia structs as typed index sets so that passing indices in the wrong order is caught by the type system at the point of misuse

Required packages

This tutorial uses the following packages:

julia> using JuMP

The problem

Consider a transshipment model of a single product from factories to markets. One way to model the flow variable is:

julia> F, M = 3, 3(3, 3)
julia> model = Model()A JuMP Model
├ solver: none
├ objective_sense: FEASIBILITY_SENSE
├ num_variables: 0
├ num_constraints: 0
└ Names registered in the model: none
julia> @variable(model, x[1:F, 1:M] >= 0)3×3 Matrix{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]

Now the question is: what does x[1, 2] represent? Was it x[f, m] or x[m, f]? In this simple example, it's probably easy to remember the flow from a factory to a market, but it gets harder if there are more sets. Was time the first index or the last?

Mis-typing the order of two sets is a common error, and it can result in bugs that are surprisingly hard to find and debug. If you didn't already know, would you be able to spot the difference in these two constraints, especially if the constraint expression was large or you were looking through many constraints trying to find out why the solution wasn't what you expected?

julia> @constraint(model, [f in 1:F], sum(x[f, m] for m in 1:M) == 1)3-element Vector{ConstraintRef{Model, MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64}, MathOptInterface.EqualTo{Float64}}, ScalarShape}}:
 x[1,1] + x[1,2] + x[1,3] = 1
 x[2,1] + x[2,2] + x[2,3] = 1
 x[3,1] + x[3,2] + x[3,3] = 1
julia> @constraint(model, [f in 1:F], sum(x[m, f] for m in 1:M) == 1)3-element Vector{ConstraintRef{Model, MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64}, MathOptInterface.EqualTo{Float64}}, ScalarShape}}:
 x[1,1] + x[2,1] + x[3,1] = 1
 x[1,2] + x[2,2] + x[3,2] = 1
 x[1,3] + x[2,3] + x[3,3] = 1

There are two approaches you can use to avoid bugs like x[m, f]: keyword indexing and typed indices.

Keyword indexing

The first approach is to use keyword indexing. Because keyword indexing does not work for Array containers, we must explicitly choose the DenseAxisArray container.

julia> F, M = 3, 3(3, 3)
julia> model = Model()A JuMP Model
├ solver: none
├ objective_sense: FEASIBILITY_SENSE
├ num_variables: 0
├ num_constraints: 0
└ Names registered in the model: none
julia> @variable(model, x[f in 1:F, m in 1:M] >= 0, container = DenseAxisArray)2-dimensional DenseAxisArray{VariableRef,2,...} with index sets:
    Dimension 1, Base.OneTo(3)
    Dimension 2, Base.OneTo(3)
And data, a 3×3 Matrix{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]

Accessing x in the correct order works:

julia> x[f=1, m=2]x[1,2]

But the wrong order errors:

julia>     x[m=2, f=1]Invalid index m in position 1.

When using keyword indexing, the indices must match the exact name and order used when creating the container.

Check the index names and their order in the container definition.

Keyword indexing means we can write constraints like this:

julia> @constraint(model, [f in 1:F], sum(x[f=f, m=m] for m in 1:M) == 1)3-element Vector{ConstraintRef{Model, MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64}, MathOptInterface.EqualTo{Float64}}, ScalarShape}}:
 x[1,1] + x[1,2] + x[1,3] = 1
 x[2,1] + x[2,2] + x[2,3] = 1
 x[3,1] + x[3,2] + x[3,3] = 1

Typed indices

A second approach is to define new types for each index:

julia> struct Factory
           x::Int
       end
julia> struct Market
           x::Int
       end

Then, we define each axis as a vector of these types:

julia> factories = Factory.(1:F)3-element Vector{Main.Factory}:
 Main.Factory(1)
 Main.Factory(2)
 Main.Factory(3)
julia> markets = Market.(1:M)3-element Vector{Main.Market}:
 Main.Market(1)
 Main.Market(2)
 Main.Market(3)
julia> model = Model()A JuMP Model
├ solver: none
├ objective_sense: FEASIBILITY_SENSE
├ num_variables: 0
├ num_constraints: 0
└ Names registered in the model: none
julia> @variable(model, x[factories, markets] >= 0)2-dimensional DenseAxisArray{VariableRef,2,...} with index sets:
    Dimension 1, Main.Factory[Main.Factory(1), Main.Factory(2), Main.Factory(3)]
    Dimension 2, Main.Market[Main.Market(1), Main.Market(2), Main.Market(3)]
And data, a 3×3 Matrix{VariableRef}:
 x[Main.Factory(1),Main.Market(1)]  …  x[Main.Factory(1),Main.Market(3)]
 x[Main.Factory(2),Main.Market(1)]     x[Main.Factory(2),Main.Market(3)]
 x[Main.Factory(3),Main.Market(1)]     x[Main.Factory(3),Main.Market(3)]

Accessing x in the correct order works:

julia> x[Factory(1), Market(2)]x[Main.Factory(1),Main.Market(2)]

But the wrong order errors:

julia>     x[Market(2), Factory(1)]KeyError: key Main.Factory(1) not found

Typed indices means we can write constraints like this:

julia> @constraint(model, [f in 1:F], sum(x[Factory(f), Market(m)] for m in 1:M) == 1)3-element Vector{ConstraintRef{Model, MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64}, MathOptInterface.EqualTo{Float64}}, ScalarShape}}:
 x[Main.Factory(1),Main.Market(1)] + x[Main.Factory(1),Main.Market(2)] + x[Main.Factory(1),Main.Market(3)] = 1
 x[Main.Factory(2),Main.Market(1)] + x[Main.Factory(2),Main.Market(2)] + x[Main.Factory(2),Main.Market(3)] = 1
 x[Main.Factory(3),Main.Market(1)] + x[Main.Factory(3),Main.Market(2)] + x[Main.Factory(3),Main.Market(3)] = 1

or like this:

julia> @constraint(model, [f in factories], sum(x[f, m] for m in markets) == 1)1-dimensional DenseAxisArray{ConstraintRef{Model, MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64}, MathOptInterface.EqualTo{Float64}}, ScalarShape},1,...} with index sets:
    Dimension 1, Main.Factory[Main.Factory(1), Main.Factory(2), Main.Factory(3)]
And data, a 3-element Vector{ConstraintRef{Model, MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64}, MathOptInterface.EqualTo{Float64}}, ScalarShape}}:
 x[Main.Factory(1),Main.Market(1)] + x[Main.Factory(1),Main.Market(2)] + x[Main.Factory(1),Main.Market(3)] = 1
 x[Main.Factory(2),Main.Market(1)] + x[Main.Factory(2),Main.Market(2)] + x[Main.Factory(2),Main.Market(3)] = 1
 x[Main.Factory(3),Main.Market(1)] + x[Main.Factory(3),Main.Market(2)] + x[Main.Factory(3),Main.Market(3)] = 1