Indexing sparse and dense axis arrays

Specific Domains Optimization (Mathematical)
1

I have a JuMP model with a mix of sparse and dense axis arrays. I find the behaviour of the indexes to be inconsistent in the sense that:

eachindex(X)

will give string-based indexes if X is a sparse axis array but cartesian ones for a dense array. Is there a way to convert the cartesian arrays to string versions?

As an example, I have code like this, which works if all arrays are sparse but breaks if one is dense.

for key in eachindex(model[:βm_t])
p,i = key;
fix(model[:βm_t][p,i], model[:βm][p,i,m], force=true);
end

  • PS how do I format the code?
2

You can format code with backticks

```julia
1+1
```

Can you share a reproducible example of what you are trying to do? I find that i rarely, if ever, need to use explicit indexing like this.

3

Hi @odow here is a minimal example. I have found a workaround, but not sure if this is the ideal approach.

Preferably there would be one way of performing this indexing operation, that works for both sparse and dense arrays…

## Sparse
model = Model()
abc = ["a","b","c"]
@expression(model, Xdata[i in abc, j in abc, m in 1:12; i==j], 1.0);
@variable(model, 1 <= month <= 12, (start=1));
@variable(model, X[i in abc, j in abc; i==j], start=0.0);

m = 1
for key in eachindex(X)
    fix(X[key], Xdata[key...,m])
end

## Dense
model = Model()
@expression(model, Xdata[i in abc, j in abc, m in 1:12], 1.0);
@variable(model, 1 <= month <= 12, (start=1));
@variable(model, X[i in abc, j in abc], start=0.0);

# This fails when array is dense
m = 1
for key in eachindex(X)
    fix(X[key], Xdata[key...,m])
end

# This works!!
for key in keys(X)
    i,j = key.I
    fix(X[i,j], Xdata[i,j,m])
end
4

Preferably there would be one way of performing this indexing operation, that works for both sparse and dense arrays

In this case there isn’t really, because you are doing the additional operation that , m. The keys of X are not compatible with the keys of Xdata.

Note that SparseAxisArrays are backed by a Dict, and DenseAxisArrays are backed by an Array. We’ve tried to paper over the differences, but they are fundamentally different objects.

This shows up in a bunch of places, for example:

julia> axes(dense_data)
(["a", "b", "c"], ["a", "b", "c"], Base.OneTo(12))

julia> axes(sparse_data)
ERROR: `Base.size` is not implemented for `SparseAxisArray` because although it is a subtype of `AbstractArray`, it is conceptually closer to a dictionary with `N`-dimensional keys. If you encounter this error and you didn't call `size` explicitly, it is because you called a method that is unsupported for `SparseAxisArray`s. Consult the JuMP documentation for a list of supported operations.
Stacktrace:
[...]

If the arrays have compatible dimensions, one approach is to use broadcasting and let JuMP figure out the details:

julia> using JuMP

julia> begin
           abc = ["a", "b", "c"]
           sparse_model = Model()
           @expression(sparse_model, sparse_data[i in abc, j in abc, m in 1:12; i == j], 1.0)
           @variable(sparse_model, sparse_x[i in abc, j in abc; i == j], start = 0.0)

           dense_model = Model()
           @expression(dense_model, dense_data[i in abc, j in abc, m in 1:12], 1.0)
           @variable(dense_model, dense_x[i in abc, j in abc], start = 0.0)
       end
2-dimensional DenseAxisArray{VariableRef,2,...} with index sets:
    Dimension 1, ["a", "b", "c"]
    Dimension 2, ["a", "b", "c"]
And data, a 3×3 Matrix{VariableRef}:
 dense_x[a,a]  dense_x[a,b]  dense_x[a,c]
 dense_x[b,a]  dense_x[b,b]  dense_x[b,c]
 dense_x[c,a]  dense_x[c,b]  dense_x[c,c]

julia> m = 1
1

julia> fix.(sparse_x, sparse_data[:, :, m])
JuMP.Containers.SparseAxisArray{Nothing, 2, Tuple{String, String}} with 3 entries:
  [a, a]  =  nothing
  [b, b]  =  nothing
  [c, c]  =  nothing

julia> fix.(dense_x, dense_data[:, :, m])
2-dimensional DenseAxisArray{Nothing,2,...} with index sets:
    Dimension 1, ["a", "b", "c"]
    Dimension 2, ["a", "b", "c"]
And data, a 3×3 Matrix{Nothing}:
 nothing  nothing  nothing
 nothing  nothing  nothing
 nothing  nothing  nothing

julia> print(sparse_model)
Feasibility
Subject to
 sparse_x[a,a] = 1
 sparse_x[b,b] = 1
 sparse_x[c,c] = 1

julia> print(dense_model)
Feasibility
Subject to
 dense_x[a,a] = 1
 dense_x[b,a] = 1
 dense_x[c,a] = 1
 dense_x[a,b] = 1
 dense_x[b,b] = 1
 dense_x[c,b] = 1
 dense_x[a,c] = 1
 dense_x[b,c] = 1
 dense_x[c,c] = 1

Another approach would be to use Julia’s multiple dispatch:

function my_fix(
    x::JuMP.Containers.SparseAxisArray,
    data::JuMP.Containers.SparseAxisArray,
    m,
)
    for key in eachindex(x)
        fix(x[key], data[key..., m])
    end
    return
end

function my_fix(
    x::JuMP.Containers.DenseAxisArray,
    data::JuMP.Containers.DenseAxisArray,
    m,
)
    for i in axes(x, 1), j in axes(x, 2)
        fix(x[i, j], data[i, j, m])
    end
    return
end
m = 1
my_fix(sparse_x, sparse_data, m)
my_fix(dense_x, dense_data, m)

But I would encourage you to think about different representations of the data that do not involve SparseAxisArray

julia> using JuMP

julia> begin
           abc = ["a", "b", "c"]
           pairs = [(i, j) for i in abc, j in abc if i == j]
           sparse_model = Model()
           @expression(sparse_model, sparse_data[pairs, 1:12], 1.0)
           @variable(sparse_model, sparse_x[pairs], start = 0.0)
           for ij in pairs
               fix(sparse_x[ij], sparse_data[ij, m])
           end
           print(sparse_model)
       end
Feasibility
Subject to
 sparse_x[("a", "a")] = 1
 sparse_x[("b", "b")] = 1
 sparse_x[("c", "c")] = 1

julia> begin
           pairs = [(i, j) for i in abc, j in abc]
           dense_model = Model()
           @expression(dense_model, dense_data[pairs, m in 1:12], 1.0)
           @variable(dense_model, dense_x[pairs], start = 0.0)
           for ij in pairs
               fix(dense_x[ij], dense_data[ij, m])
           end
           print(dense_model)
       end
Feasibility
Subject to
 dense_x[("a", "a")] = 1
 dense_x[("b", "a")] = 1
 dense_x[("c", "a")] = 1
 dense_x[("a", "b")] = 1
 dense_x[("b", "b")] = 1
 dense_x[("c", "b")] = 1
 dense_x[("a", "c")] = 1
 dense_x[("b", "c")] = 1
 dense_x[("c", "c")] = 1

Since you’re likely coming from GAMS, I’ll point you to:

The biggest lesson is that you do not need to use JuMP’s built-in data structures. They’re made to be convenient, but they may not be the best choice for your problem.

5

Wow thanks for that very fast and comprehensive reply.

Your suggestion of using sparse arrays with a single tuple index is very neat. I wonder if there are any performance gains from this approach vs using sparse arrays…

6

Yes, sparse arrays have known performance issues, in that, it is easy to write slow code that uses them. (See the sum-if formulation tutorial. It’s not exactly the same, but general rules apply. JuMP does not treat sparse indices in the same way that GAMS does, so writing code like you would in GAMS is not performant.)

7

Thanks, looking at the sum-if example it seems the run time issues apply at model creation time, given the way the macros are executed.

I have noticed issues in my own code very similar to this, but I guess my question was more: would sparse vs dense arrays affect the model solution time.

8

And the answer is: it depends how you write the code. Maybe no. Maybe yes. But it’s common for badly performing code to involve mishandled sparsity.