Extracting JuMP results

Specific Domains Optimization (Mathematical)
, , ,
1

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"],
    LINK_2 = ["A3","A7",""], LINK_3 = ["A9","A11","A12"])
    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
2

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
3

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.

2 Likes
4

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
5

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
6

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
7

Thanks