# 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
end
``````

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"],
ID_t =  collect(1:15)
ID = Vector{String}(undef,15)
for i in 1:15
ID[i]="A".*string(ID_t[i])
end

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

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

for (i, row) in enumerate(eachrow(Relationship))
for key in values(row)
if !isempty(key)
@constraint(model, t[i] == x[key])
end
end
end
@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)
optimize!(model)
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)
One*x
# 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]
end
``````

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 `DenseAxisArray`s, 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.

2 Likes

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.

5 Likes

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;
dim_names::Vector{Symbol}=Vector{Symbol}(),
value_col::Symbol=:Value)

if isempty(var)
return DataFrame()
end

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

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

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
end
``````
2 Likes

Thanks