Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

De-eval-ify vectorized unary functions over SparseMatrixCSCs, and then transition them to compact broadcast syntax #17265

Merged
merged 2 commits into from
Aug 27, 2016
Merged

De-eval-ify vectorized unary functions over SparseMatrixCSCs, and then transition them to compact broadcast syntax #17265

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
4ea088a
Rewrite vectorized unary functions over SparseMatrixCSCs, leveraging …
Sacha0 Jul 2, 2016
bd7da66
Transition vectorized unary functions over `SparseMatrixCSC`s to comp…
Sacha0 Jul 3, 2016
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
13 changes: 13 additions & 0 deletions base/deprecated.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -852,4 +852,17 @@ function convert(::Type{UpperTriangular}, A::Bidiagonal)
end
end

# Deprecate vectorized unary functions over sparse matrices in favor of compact broadcast syntax (#17265).
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

we won't be deprecating these in 0.5, move after the "to be deprecated in 0.6"

for f in (:sin, :sinh, :sind, :asin, :asinh, :asind,
:tan, :tanh, :tand, :atan, :atanh, :atand,
:sinpi, :cosc, :ceil, :floor, :trunc, :round, :real, :imag,
:log1p, :expm1, :abs, :abs2, :conj,
:log, :log2, :log10, :exp, :exp2, :exp10, :sinc, :cospi,
:cos, :cosh, :cosd, :acos, :acosd,
:cot, :coth, :cotd, :acot, :acotd,
:sec, :sech, :secd, :asech,
:csc, :csch, :cscd, :acsch)
@eval @deprecate $f(A::SparseMatrixCSC) $f.(A)
end

# End deprecations scheduled for 0.6
220 changes: 102 additions & 118 deletions base/sparse/sparsematrix.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1305,134 +1305,118 @@ end

## Unary arithmetic and boolean operators

macro _unary_op_nz2z_z2z(op,A,Tv,Ti)
esc(quote
nfilledA = nnz($A)
colptrB = Array{$Ti}($A.n+1)
rowvalB = Array{$Ti}(nfilledA)
nzvalB = Array{$Tv}(nfilledA)

nzvalA = $A.nzval
colptrA = $A.colptr
rowvalA = $A.rowval

k = 0 # number of additional zeros introduced by op(A)
@inbounds for i = 1 : $A.n
colptrB[i] = colptrA[i] - k
for j = colptrA[i] : colptrA[i+1]-1
opAj = $(op)(nzvalA[j])
if opAj == 0
k += 1
else
rowvalB[j - k] = rowvalA[j]
nzvalB[j - k] = opAj
end
end
end
colptrB[end] = $A.colptr[end] - k
deleteat!(rowvalB, colptrB[end]:nfilledA)
deleteat!(nzvalB, colptrB[end]:nfilledA)
return SparseMatrixCSC($A.m, $A.n, colptrB, rowvalB, nzvalB)
end) # quote
end

# Operations that may map nonzeros to zero, and zero to zero
# Result is sparse
for op in (:ceil, :floor, :trunc, :round,
:sin, :tan, :asin, :atan,
:sinh, :tanh, :asinh, :atanh,
:sinpi, :cosc,
:sind, :tand, :asind, :atand)
@eval begin
$(op){Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}) = @_unary_op_nz2z_z2z($op,A,Tv,Ti)
end # quote
end # macro

for op in (:real, :imag)
@eval begin
($op){Tv<:Complex,Ti}(A::SparseMatrixCSC{Tv,Ti}) = @_unary_op_nz2z_z2z($op,A,Tv.parameters[1],Ti)
end # quote
end # macro
real{Tv<:Number,Ti}(A::SparseMatrixCSC{Tv,Ti}) = copy(A)
imag{Tv<:Number,Ti}(A::SparseMatrixCSC{Tv,Ti}) = spzeros(Tv, Ti, A.m, A.n)

for op in (:ceil, :floor, :trunc, :round)
@eval begin
($op){T,Tv,Ti}(::Type{T},A::SparseMatrixCSC{Tv,Ti}) = @_unary_op_nz2z_z2z($op,A,T,Ti)
end # quote
end # macro


# Operations that map nonzeros to nonzeros, and zeros to zeros
# Result is sparse
for op in (:-, :log1p, :expm1)
@eval begin

function ($op)(A::SparseMatrixCSC)
B = similar(A)
nzvalB = B.nzval
nzvalA = A.nzval
@simd for i=1:length(nzvalB)
@inbounds nzvalB[i] = ($op)(nzvalA[i])
end
return B
end

end
end

function abs{Tv<:Complex,Ti}(A::SparseMatrixCSC{Tv,Ti})
T = Tv.parameters[1]
(T <: Integer) && (T = (T <: BigInt) ? BigFloat : Float64)
@_unary_op_nz2z_z2z(abs,A,T,Ti)
end
abs2{Tv<:Complex,Ti}(A::SparseMatrixCSC{Tv,Ti}) = @_unary_op_nz2z_z2z(abs2,A,Tv.parameters[1],Ti)
for op in (:abs, :abs2)
@eval begin
function ($op){Tv<:Number,Ti}(A::SparseMatrixCSC{Tv,Ti})
B = similar(A)
nzvalB = B.nzval
nzvalA = A.nzval
@simd for i=1:length(nzvalB)
@inbounds nzvalB[i] = ($op)(nzvalA[i])
"""
Helper macro for the unary broadcast definitions below. Takes parent method `fp` and a set
of desired child methods `fcs`, and builds an expression defining each of the child methods
such that `broadcast(::typeof(fc), A::SparseMatrixCSC) = fp(fc, A)`.
"""
macro _enumerate_childmethods(fp, fcs...)
fcexps = Expr(:block)
for fc in fcs
push!(fcexps.args, :( broadcast(::typeof($(esc(fc))), A::SparseMatrixCSC) = $(esc(fp))($(esc(fc)), A) ) )
end
return fcexps
end

# Operations that map zeros to zeros and may map nonzeros to zeros, yielding a sparse matrix
"""
Takes unary function `f` that maps zeros to zeros and may map nonzeros to zeros, and returns
a new `SparseMatrixCSC{TiA,TvB}` `B` generated by applying `f` to each nonzero entry in
`A` and retaining only the resulting nonzeros.
"""
function _broadcast_unary_nz2z_z2z_T{TvA,TiA,TvB}(f::Function, A::SparseMatrixCSC{TvA,TiA}, ::Type{TvB})
Bcolptr = Array{TiA}(A.n + 1)
Browval = Array{TiA}(nnz(A))
Bnzval = Array{TvB}(nnz(A))
Bk = 1
@inbounds for j in 1:A.n
Bcolptr[j] = Bk
for Ak in nzrange(A, j)
x = f(A.nzval[Ak])
if x != 0
Browval[Bk] = A.rowval[Ak]
Bnzval[Bk] = x
Bk += 1
end
return B
end
end
end

Bcolptr[A.n + 1] = Bk
resize!(Browval, Bk - 1)
resize!(Bnzval, Bk - 1)
return SparseMatrixCSC(A.m, A.n, Bcolptr, Browval, Bnzval)
end
function _broadcast_unary_nz2z_z2z{Tv}(f::Function, A::SparseMatrixCSC{Tv})
_broadcast_unary_nz2z_z2z_T(f, A, Tv)
end
@_enumerate_childmethods(_broadcast_unary_nz2z_z2z,
sin, sinh, sind, asin, asinh, asind,
tan, tanh, tand, atan, atanh, atand,
sinpi, cosc, ceil, floor, trunc, round)
broadcast(::typeof(real), A::SparseMatrixCSC) = copy(A)
broadcast{Tv,Ti}(::typeof(imag), A::SparseMatrixCSC{Tv,Ti}) = spzeros(Tv, Ti, A.m, A.n)
broadcast{TTv}(::typeof(real), A::SparseMatrixCSC{Complex{TTv}}) = _broadcast_unary_nz2z_z2z_T(real, A, TTv)
broadcast{TTv}(::typeof(imag), A::SparseMatrixCSC{Complex{TTv}}) = _broadcast_unary_nz2z_z2z_T(imag, A, TTv)
ceil{To}(::Type{To}, A::SparseMatrixCSC) = _broadcast_unary_nz2z_z2z_T(ceil, A, To)
floor{To}(::Type{To}, A::SparseMatrixCSC) = _broadcast_unary_nz2z_z2z_T(floor, A, To)
trunc{To}(::Type{To}, A::SparseMatrixCSC) = _broadcast_unary_nz2z_z2z_T(trunc, A, To)
round{To}(::Type{To}, A::SparseMatrixCSC) = _broadcast_unary_nz2z_z2z_T(round, A, To)

# Operations that map zeros to zeros and map nonzeros to nonzeros, yielding a sparse matrix
"""
Takes unary function `f` that maps zeros to zeros and nonzeros to nonzeros, and returns a
new `SparseMatrixCSC{TiA,TvB}` `B` generated by applying `f` to each nonzero entry in `A`.
"""
function _broadcast_unary_nz2nz_z2z_T{TvA,TiA,TvB}(f::Function, A::SparseMatrixCSC{TvA,TiA}, ::Type{TvB})
Bcolptr = Vector{TiA}(A.n + 1)
Browval = Vector{TiA}(nnz(A))
Bnzval = Vector{TvB}(nnz(A))
copy!(Bcolptr, 1, A.colptr, 1, A.n + 1)
copy!(Browval, 1, A.rowval, 1, nnz(A))
@inbounds @simd for k in 1:nnz(A)
Bnzval[k] = f(A.nzval[k])
end
return SparseMatrixCSC(A.m, A.n, Bcolptr, Browval, Bnzval)
end
function _broadcast_unary_nz2nz_z2z{Tv}(f::Function, A::SparseMatrixCSC{Tv})
_broadcast_unary_nz2nz_z2z_T(f, A, Tv)
end
@_enumerate_childmethods(_broadcast_unary_nz2nz_z2z,
log1p, expm1, abs, abs2, conj)
broadcast{TTv}(::typeof(abs2), A::SparseMatrixCSC{Complex{TTv}}) = _broadcast_unary_nz2nz_z2z_T(abs2, A, TTv)
broadcast{TTv}(::typeof(abs), A::SparseMatrixCSC{Complex{TTv}}) = _broadcast_unary_nz2nz_z2z_T(abs, A, TTv)
broadcast{TTv<:Integer}(::typeof(abs), A::SparseMatrixCSC{Complex{TTv}}) = _broadcast_unary_nz2nz_z2z_T(abs, A, Float64)
broadcast{TTv<:BigInt}(::typeof(abs), A::SparseMatrixCSC{Complex{TTv}}) = _broadcast_unary_nz2nz_z2z_T(abs, A, BigFloat)
function conj!(A::SparseMatrixCSC)
nzvalA = A.nzval
@simd for i=1:length(nzvalA)
@inbounds nzvalA[i] = conj(nzvalA[i])
@inbounds @simd for k in 1:nnz(A)
A.nzval[k] = conj(A.nzval[k])
end
return A
end

conj(A::SparseMatrixCSC) = conj!(copy(A))

# Operations that map nonzeros to nonzeros, and zeros to nonzeros
# Result is dense
for op in (:cos, :cosh, :acos, :sec, :csc, :cot, :acot, :sech,
:csch, :coth, :asech, :acsch, :cospi, :sinc, :cosd,
:cotd, :cscd, :secd, :acosd, :acotd, :log, :log2, :log10,
:exp, :exp2, :exp10)
@eval begin

function ($op){Tv}(A::SparseMatrixCSC{Tv})
B = fill($(op)(zero(Tv)), size(A))
@inbounds for col = 1 : A.n
for j = A.colptr[col] : A.colptr[col+1]-1
row = A.rowval[j]
nz = A.nzval[j]
B[row,col] = $(op)(nz)
end
end
return B
# Operations that map both zeros and nonzeros to zeros, yielding a dense matrix
"""
Takes unary function `f` that maps both zeros and nonzeros to nonzeros, constructs a new
`Matrix{TvB}` `B`, populates all entries of `B` with the result of a single `f(one(zero(Tv)))`
call, and then, for each stored entry in `A`, calls `f` on the entry's value and stores the
result in the corresponding location in `B`.
"""
function _broadcast_unary_nz2nz_z2nz{Tv}(f::Function, A::SparseMatrixCSC{Tv})
B = fill(f(zero(Tv)), size(A))
@inbounds for j in 1:A.n
for k in nzrange(A, j)
i = A.rowval[k]
x = A.nzval[k]
B[i,j] = f(x)
end

end
return B
end
@_enumerate_childmethods(_broadcast_unary_nz2nz_z2nz,
log, log2, log10, exp, exp2, exp10, sinc, cospi,
cos, cosh, cosd, acos, acosd,
cot, coth, cotd, acot, acotd,
sec, sech, secd, asech,
csc, csch, cscd, acsch)


## Broadcasting kernels specialized for returning a SparseMatrixCSC
Expand Down Expand Up @@ -1734,7 +1718,7 @@ end # macro
(.^)(A::SparseMatrixCSC, B::Number) =
B==0 ? sparse(ones(typeof(one(eltype(A)).^B), A.m, A.n)) :
SparseMatrixCSC(A.m, A.n, copy(A.colptr), copy(A.rowval), A.nzval .^ B)
(.^)(::Irrational{:e}, B::SparseMatrixCSC) = exp(B)
(.^)(::Irrational{:e}, B::SparseMatrixCSC) = exp.(B)
(.^)(A::Number, B::SparseMatrixCSC) = (.^)(A, full(B))
(.^)(A::SparseMatrixCSC, B::Array) = (.^)(full(A), B)
(.^)(A::Array, B::SparseMatrixCSC) = (.^)(A, full(B))
Expand Down
67 changes: 40 additions & 27 deletions test/sparsedir/sparse.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -286,7 +286,7 @@ end

# conj
cA = sprandn(5,5,0.2) + im*sprandn(5,5,0.2)
@test full(conj(cA)) == conj(full(cA))
@test full(conj.(cA)) == conj(full(cA))

# Test SparseMatrixCSC [c]transpose[!] and permute[!] methods
let smalldim = 5, largedim = 10, nzprob = 0.4
Expand Down Expand Up @@ -461,22 +461,47 @@ end
@test maximum(sparse(-ones(3,3))) == -1
@test minimum(sparse(ones(3,3))) == 1

# Unary functions
a = sprand(5,15, 0.5)
afull = full(a)
for op in (:sin, :cos, :tan, :ceil, :floor, :abs, :abs2)
@eval begin
@test ($op)(afull) == full($(op)(a))
end
end

for op in (:ceil, :floor)
@eval begin
@test ($op)(Int,afull) == full($(op)(Int,a))
# Test unary functions with specialized broadcast over SparseMatrixCSCs
let
A = sprand(5, 15, 0.5)
C = A + im*A
Afull = full(A)
Cfull = full(C)
# Test representatives of [unary functions that map zeros to zeros and may map nonzeros to zeros]
@test sin.(Afull) == full(sin.(A))
@test tan.(Afull) == full(tan.(A)) # should be redundant with sin test
@test ceil.(Afull) == full(ceil.(A))
@test floor.(Afull) == full(floor.(A)) # should be redundant with ceil test
@test real.(Afull) == full(real.(A))
@test imag.(Afull) == full(imag.(A))
@test real.(Cfull) == full(real.(C))
@test imag.(Cfull) == full(imag.(C))
# Test representatives of [unary functions that map zeros to zeros and nonzeros to nonzeros]
@test expm1.(Afull) == full(expm1.(A))
@test abs.(Afull) == full(abs.(A))
@test abs2.(Afull) == full(abs2.(A))
@test abs.(Cfull) == full(abs.(C))
@test abs2.(Cfull) == full(abs2.(C))
# Test representatives of [unary functions that map both zeros and nonzeros to nonzeros]
@test cos.(Afull) == full(cos.(A))
# Test representatives of remaining vectorized-nonbroadcast unary functions
@test ceil(Int, Afull) == full(ceil(Int, A))
@test floor(Int, Afull) == full(floor(Int, A))
# Tests of real, imag, abs, and abs2 for SparseMatrixCSC{Int,X}s previously elsewhere
for T in (Int, Float16, Float32, Float64, BigInt, BigFloat)
R = rand(T[1:100;], 2, 2)
I = rand(T[1:100;], 2, 2)
D = R + I*im
S = sparse(D)
@test R == real.(S)
@test I == imag.(S)
@test real.(sparse(R)) == R
@test nnz(imag.(sparse(R))) == 0
@test abs.(S) == abs(D)
@test abs2.(S) == abs2(D)
end
end


# getindex tests
ni = 23
nj = 32
Expand Down Expand Up @@ -872,7 +897,7 @@ end
@test_throws ArgumentError sparsevec(Dict(-1=>1,1=>2))

# issue #8976
@test conj(sparse([1im])) == sparse(conj([1im]))
@test conj.(sparse([1im])) == sparse(conj([1im]))
@test conj!(sparse([1im])) == sparse(conj!([1im]))

# issue #9525
Expand Down Expand Up @@ -1038,18 +1063,6 @@ end
x = speye(100)
@test_throws BoundsError x[-10:10]

for T in (Int, Float16, Float32, Float64, BigInt, BigFloat)
let R=rand(T[1:100;],2,2), I=rand(T[1:100;],2,2)
D = R + I*im
S = sparse(D)
@test R == real(S)
@test I == imag(S)
@test real(sparse(R)) == R
@test nnz(imag(sparse(R))) == 0
@test abs(S) == abs(D)
@test abs2(S) == abs2(D)
end
end

# issue #10407
@test maximum(spzeros(5, 5)) == 0.0
Expand Down