getvalue method for sparse matrices of variables

[ExpandingMan]

There currently isn't one. Here's one I wrote, I'm not sure if it's sufficiently general to get added to JuMP

function getvalue{Ti <: Integer}(x::SparseMatrixCSC{JuMP.Variable, Ti})
    x_val = spzeros(size(x)...)
    rows = rowvals(x)
    vars = nonzeros(x)
    m, n = size(x)
    for j ∈ 1:n
        for idx ∈ nzrange(x, j)
            var = vars[idx]
            # this is the row number
            i = rows[idx]
            x_val[i, j] = getvalue(var)
        end
    end
    x_val
end

[mlubin]

You should use the new broadcasting syntax or map for this.

On Nov 9, 2016 12:22 PM, "ExpandingMan" notifications@github.com wrote:

There currently isn't one. Here's one I wrote, I'm not sure if it's
sufficiently general to get added to JuMP

function getvalue{Ti <: Integer}(x::SparseMatrixCSC{JuMP.Variable, Ti})
x_val = spzeros(size(x)...)
rows = rowvals(x)
vars = nonzeros(x)
m, n = size(x)
for j ∈ 1:n
for idx ∈ nzrange(x, j)
var = vars[idx]
# this is the row number
i = rows[idx]
x_val[i, j] = getvalue(var)
end
end
x_valend

—
You are receiving this because you are subscribed to this thread.
Reply to this email directly, view it on GitHub
#889, or mute the thread
https://github.com/notifications/unsubscribe-auth/ABp0M9NGxHxjWdPSPAi5H55kGC_8IgDVks5q8gFAgaJpZM4Ktxtj
.

[ExpandingMan]

The new broadcasting syntax returns a dense matrix (which certainly isn't the ideal behavior), but map (which I hadn't thought of) seems to work just fine. Thanks.

It probably wouldn't be a bad idea to still put

getvalue{Ti<:Integer}(x:SparseMatrixCSC{JuMP.Variable, Ti}) = map(getvalue, x)

into JuMP because of the bad broadcast behavior.

[ExpandingMan]

Thanks all.

[ExpandingMan]

Hello all, I feel like I'm going crazy, because I don't see how this can be, but map seems to be unbelievably slow for this. Give me a little bit to create a benchmark program...

[ExpandingMan]

Ok, see below

using JuMP

const N = 10^6
const M = 3
const Nvar = 10^4

m = Model()

x = spzeros(JuMP.Variable, N, M)

x[1:Nvar, 1] = @variable(m, [1:Nvar], lowerbound=0, basename="x")

@objective(m, Min, -sum(x))
@constraint(m, x*ones(Int64, M) .≤ 1.0)

status = solve(m)

function getvalue_map{Ti<:Integer}(x::SparseMatrixCSC{JuMP.Variable, Ti})
    map(getvalue, x)
end

function getvalue_dumbloop{Ti<:Integer}(x::SparseMatrixCSC{JuMP.Variable, Ti})
    x_val = spzeros(size(x)...)
    for j ∈ size(x, 2), i ∈ size(x, 1)
        x_val[i, j] = getvalue(x[i, j]) 
    end
    x_val
end

function getvalue_loop{Ti<:Integer}(x::SparseMatrixCSC{JuMP.Variable, Ti})
    x_val = spzeros(size(x)...)
    rows = rowvals(x)
    vars = nonzeros(x)
    m, n = size(x)
    for j ∈ 1:n
        for idx ∈ nzrange(x, j)
            var = vars[idx]
            # this is the row number
            i = rows[idx]
            x_val[i, j] = getvalue(var)
        end
    end
    x_val
end

@time x_val = getvalue_dumbloop(x)
>>>   0.073583 seconds (26.33 k allocations: 1.107 MB)
@time x_val = getvalue_map(x)
>>>   0.753596 seconds (12.04 M allocations: 596.768 MB, 37.09% gc time)
@time x_val = getvalue_loop(x)
>>>  0.013843 seconds (5.87 k allocations: 771.323 KB)

I feel like this is an issue with map, but I'm not sure what exactly I should say when I raise the issue. I'm actually not even sure whether this is an intended use of map (i.e. broadcasting returns a dense matrix, shouldn't map have the same behavior?). Regardless, getvalue should probably not be defined using map in JuMP.

[blegat]

Hey, thanks for pointing this, it is indeed surprising, we should definitely not merge #891.
It seems to be an issue in Julia Base, since I get the same performance issues with this code:

const N = 10^6
const M = 3
const Nvar = 10^4

typealias T Int

x = spzeros(T, N, M)

x[1:Nvar, 1] = 3

f = x->(x == 3 ? 1 : 0)

function getvalue_map{Ti<:Integer}(x::SparseMatrixCSC{T, Ti})
    map(f, x)
end

function getvalue_dumbloop{Ti<:Integer}(x::SparseMatrixCSC{T, Ti})
    x_val = spzeros(size(x)...)
    for j ∈ size(x, 2), i ∈ size(x, 1)
        x_val[i, j] = f(x[i, j])
    end
    x_val
end

function getvalue_loop{Ti<:Integer}(x::SparseMatrixCSC{T, Ti})
    x_val = spzeros(size(x)...)
    rows = rowvals(x)
    vars = nonzeros(x)
    m, n = size(x)
    for j ∈ 1:n
        for idx ∈ nzrange(x, j)
            var = vars[idx]
            # this is the row number
            i = rows[idx]
            x_val[i, j] = f(var)
        end
    end
    x_val
end

The fact that broadcasting returns a dense matrix isn't the expected behavior as this behavior will change soon once JuliaLang/julia#19239 is merged.

[ExpandingMan]

Ok, I have found issue JuliaLang/julia#7010 which goes into some detail on this. It is not a bug: the function argument to map is being explicitly evaluated on every element. The fact that it does this then returns a sparse (rather than dense) matrix is incredibly confusing to me, but when I think carefully about it I guess it is the intended behavior.

One partial alternative is to do map(f, nonzeros(x)), but this returns a flattened tensor, so that certainly doesn't seem appropriate in most cases. Therefore, as far as I know, the only solution at the moment is the horrible loop that I wrote above. If anybody knows of a simpler function for achieving this, please let me know.

I really don't understand why there isn't something like this in Base

for (i, j) ∈ spiterator(x)
    x[i, j] = f(x[i, j])
end

which would run over columns first, then rows.

[Sacha0]

Re. map and map! lacking specializations for sparse matrices and vectors, filed an issue: JuliaLang/julia#19363.

In the mean time,

function getvalue_map(A::SparseMatrixCSC)
    C = copy(A)
    map!(f, C.nzval) # or map!(f, nonzeros(C)) if you prefer
    return C
end

or more efficiently

getvalue_map(A::SparseMatrixCSC) =
    SparseMatrixCSC(A.m, A.n, copy(A.colptr), copy(A.rowval), map(f, A.nzval))

?

Best!

[ExpandingMan]

Thanks @Sacha0! Please note that the first solution doesn't work in this case, because the type gets converted by getvalue, however the second solution works and is quite efficient.