Extracting JuMP results

It seems like matrix multiplication does not work in the usual way with DenseAxisArray,
If I use, ones(5,15)*x using the code in the MWE it throws an error. Interestingly,
array slicing ones(5,15)[1,:]'*x works fine, but if I loop through it I get an error.
The error message points to aff_expr.jl:319, so after a bit of battle with JuMP, I managed to fix it by defining an affine expression. However, just to multiply the matrices I am using this loop in an inefficient manner as iterations are through the rows.

One =ones(5,15)
Affine=[AffExpr(0) for i = 1:5]
    for i in eachindex(Affine)
        Affine[i] = One[i,:]'*x

Extracting the results into a normal array/vector after solving the problem is proving too difficult.
I have two questions:
1. How can I extract optimization results into a normal array/vector or a dataframe? The standard method for extracting results do not seem to work with DenseAxisArray.
2. I don’t mind using a loop to multiply matrices, but is there a more efficient way of multiplying matrices than that in the code above?

using JuMP, Clp, DataFrames
function get_data()
    Relationship = DataFrame(ID = ["A1","A2","A4"], LINK_1 = ["A5","A6","A10"],
    LINK_2 = ["A3","A7",""], LINK_3 = ["A9","A11","A12"])
    ID_t =  collect(1:15)
    ID = Vector{String}(undef,15)
    for i in 1:15

    PRICE = rand(10.0:100.0,15)
    Master = DataFrame(ID = ID,PRICE = PRICE)
    return Relationship, Master

    Relationship, Master = get_data()
    price = Master.PRICE
    model = Model(Clp.Optimizer)
    @variable(model, 0<=x[Master.ID]<=1)

    for (i, row) in enumerate(eachrow(Relationship))
    for key in values(row)
        if !isempty(key)
            @constraint(model, t[i] == x[key])
@objective(model,Min,sum(price*x[id] for (id, price) in zip(Master.ID, Master.PRICE))
@constraint(model,con1, sum(price*x[id] for (id, price) in zip(Master.ID, Master.PRICE)) >=200)
@constraint(model,con2, sum(price*x[id] for (id, price) in zip(Master.ID, Master.PRICE)) <=300)
Res = value.(x)

This part throws an error if I try to multiply x by any matrix (apart from vectors) before solving the model:

One = ones(5,15)
# This can be fixed by
#One =ones(5,15)
#Affine=[AffExpr(0) for i = 1:5]
 #   for i in eachindex(Affine)
  #     Affine[i] = One[i,:]'*x
  #  end

I have managed to find an answer myself to the first question. Second question remains unanswered.
This is how I would extract results:

Normal_results = Vector{Float64}(undef,15)
for i in eachindex(Res)
    Normal_results[i]= Res[i]

In the near-term, you can get the underlying array object with x.data.

In the longer-term, this is something we should think about, because it has come up quite a bit.

The key issue is that

model = Model()
S = ["A", "B"]
@variable(model, x[S])
A = rand(2, 2)
c = rand(2)

c' * x  # Works
A * x  # Doesn't work.

One problem is that it isn’t obvious what type A * x should return. A vector? A DenseAxisArray?

Perhaps what we need to do is define * between two DenseAxisArrays, which checks if the named dimensions are compatible, and potentially include a conversion from Array to a DenseAxisArray with OneTo dimensions. cc @blegat, because he has a better idea of the design here.


Thank you. A long term solution will take some time to implement. My suggestion will be adding a one liner to the documentation on extracting results with this container, and some guidance on when to use expressions. It could have saved me several hours .

1 Like

We’re in the process of overhauling much of our documentation, but many hands make light work.

You can edit the documentation in the browser by choosing a file here:
https://github.com/jump-dev/JuMP.jl/tree/master/docs/src (e.g., https://github.com/jump-dev/JuMP.jl/blob/master/docs/src/variables.md) and clicking the Pencil icon in the top right.

Make the changes, then scroll to the bottom and follow the instructions to open a pull request.

It’s sometimes hard for core developers to know what parts need documenting because we forget which parts are not obvious to new users.


In regards to Question 1., this is the code I use to convert optimisation results to DataFrames:

Returns a `DataFrame` with the values of the variables from the JuMP container `var`.
The column names of the `DataFrame` can be specified for the indexing columns in `dim_names`,
and the name of the data value column by a Symbol `value_col` e.g. :Value
function convert_jump_container_to_df(var::JuMP.Containers.DenseAxisArray;

    if isempty(var)
        return DataFrame()

    if length(dim_names) == 0
        dim_names = [Symbol("dim$i") for i in 1:length(var.axes)]

    if length(dim_names) != length(var.axes)
        throw(ArgumentError("Length of given name list does not fit the number of variable dimensions"))

    tup_dim = (dim_names...,)

    # With a product over all axis sets of size M, form an Mx1 Array of all indices to the JuMP container `var`
    ind = reshape([collect(k[i] for i in 1:length(dim_names)) for k in Base.Iterators.product(var.axes...)],:,1)

    var_val  = value.(var)

    df = DataFrame([merge(NamedTuple{tup_dim}(ind[i]), NamedTuple{(value_col,)}(var_val[(ind[i]...,)...])) for i in 1:length(ind)])

    return df