Iterating SparseAxisArray does not preserve order

[odow]

x-ref https://discourse.julialang.org/t/unexpected-jump-variable-order-when-indexed-by-tuple/110262

julia> using JuMP

julia> model = Model();

julia> N = 1:2;

julia> M = [[1.0, 11.0, 100.0], [0.0, 17.0, 104.0]];

julia> @variable(model, foo[n in N, M[n]])
JuMP.Containers.SparseAxisArray{VariableRef, 2, Tuple{Int64, Float64}} with 6 entries:
  [1, 1.0  ]  =  foo[1,1.0]
  [1, 100.0]  =  foo[1,100.0]
  [1, 11.0 ]  =  foo[1,11.0]
  [2, 0.0  ]  =  foo[2,0.0]
  [2, 104.0]  =  foo[2,104.0]
  [2, 17.0 ]  =  foo[2,17.0]

julia> foo.data
Dict{Tuple{Int64, Float64}, VariableRef} with 6 entries:
  (2, 17.0)  => foo[2,17.0]
  (1, 11.0)  => foo[1,11.0]
  (1, 1.0)   => foo[1,1.0]
  (1, 100.0) => foo[1,100.0]
  (2, 0.0)   => foo[2,0.0]
  (2, 104.0) => foo[2,104.0]

julia> @constraint(model, foo[1, :] in SecondOrderCone())
[foo[1,11.0], foo[1,1.0], foo[1,100.0]] ∈ MathOptInterface.SecondOrderCone(3)

User could reasonably expect constraint function to be [foo[1,1.0], foo[1,11.0], foo[1,100.0]].

This is because we store the data in a Dict.

We could swap to an OrderedDict, or throw an error or a warning when we convert a SparseAxisArray to a vector.

[odow]

This is actually pretty nasty.

We fallback to AbstractVector methods:

julia> using JuMP

julia> A = [[1, 2, 10], [2, 3, 30]]
2-element Vector{Vector{Int64}}:
 [1, 2, 10]
 [2, 3, 30]

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

julia> @variable(model, x[i in 1:2, j in A[i]])
JuMP.Containers.SparseAxisArray{VariableRef, 2, Tuple{Int64, Int64}} with 6 entries:
  [1, 1 ]  =  x[1,1]
  [1, 10]  =  x[1,10]
  [1, 2 ]  =  x[1,2]
  [2, 2 ]  =  x[2,2]
  [2, 3 ]  =  x[2,3]
  [2, 30]  =  x[2,30]

julia> @constraint(model, x[1, :] in SecondOrderCone())
[x[1,2], x[1,10], x[1,1]] ∈ MathOptInterface.SecondOrderCone(3)

julia> x[1, :] isa AbstractVector{<:AbstractJuMPScalar}
true

entry point is:

https://github.com/jump-dev/JuMP.jl/blob/a21e616b28e101b8ba4382fb23e2a8d5b3370653/src/macros/%40constraint.jl#L842-L848

which calls:

https://github.com/jump-dev/JuMP.jl/blob/a21e616b28e101b8ba4382fb23e2a8d5b3370653/src/constraints.jl#L652-L659

And eachindex isn't ordered:

julia> y = x[1, :]
JuMP.Containers.SparseAxisArray{VariableRef, 1, Tuple{Int64}} with 3 entries:
  [1 ]  =  x[1,1]
  [10]  =  x[1,10]
  [2 ]  =  x[1,2]

julia> [y[k] for k in eachindex(y)]
3-element Vector{VariableRef}:
 x[1,2]
 x[1,10]
 x[1,1]

I don't know what a good answer is.

We can't throw an error because that will break existing code like this:

julia> @constraint(model, x[1, :] >= 0)
[x[1,2], x[1,10], x[1,1]] ∈ MathOptInterface.Nonnegatives(3)

The issue is that some vector sets are ordered, and some are not. Here's another set that is problematic:

julia> @constraint(model, x[1, :] in SecondOrderCone())
[x[1,2], x[1,10], x[1,1]] ∈ MathOptInterface.SecondOrderCone(3)

[blegat]

The user could still do dot(a, x[1, :]) where a is some vector of constant so I don't think this is about vector constraints. Maybe we could use an OrderedDict

[odow]

Hmm:

julia> using JuMP

julia> A = [[1, 2, 10], [2, 3, 30]]
2-element Vector{Vector{Int64}}:
 [1, 2, 10]
 [2, 3, 30]

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

julia> @variable(model, x[i in 1:2, j in A[i]])
JuMP.Containers.SparseAxisArray{VariableRef, 2, Tuple{Int64, Int64}} with 6 entries:
  [1, 1 ]  =  x[1,1]
  [1, 10]  =  x[1,10]
  [1, 2 ]  =  x[1,2]
  [2, 2 ]  =  x[2,2]
  [2, 3 ]  =  x[2,3]
  [2, 30]  =  x[2,30]

julia> y = x[1, :]
JuMP.Containers.SparseAxisArray{VariableRef, 1, Tuple{Int64}} with 3 entries:
  [1 ]  =  x[1,1]
  [10]  =  x[1,10]
  [2 ]  =  x[1,2]

julia> z = rand(3)
3-element Vector{Float64}:
 0.508529766285344
 0.32511007041899975
 0.6771958164053773

julia> using LinearAlgebra

julia> dot(z, y)
0.508529766285344 x[1,2] + 0.32511007041899975 x[1,10] + 0.6771958164053773 x[1,1]

This is more than just a slicing issue:

julia> model = Model();

julia> @variable(model, x[i in 1:2], container = SparseAxisArray)
JuMP.Containers.SparseAxisArray{VariableRef, 1, Tuple{Int64}} with 2 entries:
  [1]  =  x[1]
  [2]  =  x[2]

julia> dot(1:2, x)
x[2] + 2 x[1]

[odow]

So perhaps we do need an OrderedDict