diff --git a/base/deprecated.jl b/base/deprecated.jl index 151709861bb00..667bd10c5c4e7 100644 --- a/base/deprecated.jl +++ b/base/deprecated.jl @@ -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). +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 diff --git a/base/sparse/sparsematrix.jl b/base/sparse/sparsematrix.jl index d8d3fa47d0f42..87cb03485d304 100644 --- a/base/sparse/sparsematrix.jl +++ b/base/sparse/sparsematrix.jl @@ -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 @@ -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)) diff --git a/test/sparsedir/sparse.jl b/test/sparsedir/sparse.jl index 4cfc2b248097d..98c10af32dd02 100644 --- a/test/sparsedir/sparse.jl +++ b/test/sparsedir/sparse.jl @@ -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 @@ -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 @@ -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 @@ -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