Improve support for relational algebra

[odow]

Motivated by https://jump.dev/2023/07/20/gams-blog/

To make this syntax a bit nicer

[odow]

The alternative is to parse SparseAxisArray into a data frame: #3441 (comment)

Ideally, we should parse the condition of the JuMP.Containers.NestedIterator then it would be an expression tree of conditions with && and ||. Then we put it into conjunctive normal form c1 && c2 && ... && cm. Then for each condition, we should check: which indices does it involve ? is this an equality lhs == rhs
If it involves a subset of the indices then we should filter with this condition earlier and not in the most nested level of the for loops.
If it is an equality, we could do the indexing ourself (as in https://discourse.julialang.org/t/is-it-possible-to-incorporate-the-relational-algebra-technique-into-jump-in-the-future-to-expedite-model-generation/101343/4?u=blegat) or use DataFrames but that's kind of an heavy dependency for just this.
Now what's a bit tricky is that you could have [i = I, j = J, k = K, l = L; f(i, j, k) == g(j, k, l)]. Then, it means you need to build the set of (i, j, k) and the set of (j, k, l) and then index them with the Dict to take the intersection. One could wonder if it wouldn't be better to just to the 4-nested for loops instead of having to do twice a 3-nested for loops. I think however than unless K has only 2 elements, it is always better to do two 3-nested for-loops than a 4-nested for loops so doing this intersection with dictionaries is always better.

In some cases, it might be better to do the nested for loops in an order that is different from the order of the indices used by the user but that can be left as future work since I don't think this is easy to find that. But detecting the two things I mentioned above and improving the iteration shouldn't be too complicated and it would solve the issue in the GAMS example without the user to do anything. The SparseAxisArray would automatically build itself efficiently.

[odow]

This tutorial is relevant for any changes in the future looking at this:

https://jump.dev/JuMP.jl/stable/tutorials/linear/multi_commodity_network/

[odow]

See also #3212

[odow]

Perhaps I do like #3441

julia> using JuMP, DataFrames

julia> function JuMP.Containers.container(
           f::Function,
           indices,
           ::Type{DataFrames.DataFrame},
           names,
       )
           rows = vec(collect(indices))
           df = DataFrames.DataFrame(NamedTuple{tuple(names...)}(arg) for arg in rows)
           df.value = [f(arg...) for arg in rows]
           return df
       end

julia> model = Model();

julia> @variables(model, begin
           x[i in 1:2, j in i:2], (container = DataFrames.DataFrame)
           y[i in 1:2, 1:2], (container = DataFrames.DataFrame)
           z[i in 1:2, j in 1:2; isodd(i+j)], (container = DataFrames.DataFrame)
       end);

julia> x
3×3 DataFrame
 Row │ i      j      value     
     │ Int64  Int64  GenericV… 
─────┼─────────────────────────
   1 │     1      1  x[1,1]
   2 │     1      2  x[1,2]
   3 │     2      2  x[2,2]

julia> y
4×3 DataFrame
 Row │ i      ##288  value     
     │ Int64  Int64  GenericV… 
─────┼─────────────────────────
   1 │     1      1  y[1,1]
   2 │     2      1  y[2,1]
   3 │     1      2  y[1,2]
   4 │     2      2  y[2,2]

julia> z
2×3 DataFrame
 Row │ i      j      value     
     │ Int64  Int64  GenericV… 
─────┼─────────────────────────
   1 │     1      2  z[1,2]
   2 │     2      1  z[2,1]

[odow]

If we supported that syntax, then improvements to how we parsed and constructed the sparse index sets would just be a performance optimization.