Merge pull request #13440 from JuliaLang/mb/sparsevec

RFC: SparseVectors, Take 2

Author: mbauman

NEWS.md

@@ -37,6 +37,11 @@ Library improvements
 
   * New method for generic QR with column pivoting ([#13480]).
 
+  * A new `SparseVector` type allows for one-dimensional sparse arrays. Slicing
+    and reshaping sparse matrices now return vectors when appropriate. The
+    `sparsevec` function returns a one-dimensional sparse vector instead of a
+    one-column sparse matrix.
+
 Deprecated or removed
 ---------------------
 

base/deprecated.jl

@@ -852,3 +852,5 @@ end
 @deprecate cor(X::AbstractMatrix; vardim=1, mean=Base.mean(X, vardim)) corm(X, mean, vardim)
 @deprecate cor(x::AbstractVector, y::AbstractVector; mean=(Base.mean(x), Base.mean(y))) corm(x, mean[1], y, mean[2])
 @deprecate cor(X::AbstractVecOrMat, Y::AbstractVecOrMat; vardim=1, mean=(Base.mean(X, vardim), Base.mean(Y, vardim))) corm(X, mean[1], Y, mean[2], vardim)
+
+@deprecate_binding SparseMatrix SparseArrays

base/docs/helpdb.jl

@@ -2387,17 +2387,6 @@ Right division operator: multiplication of `x` by the inverse of `y` on the righ
 """
 Base.(:(/))
 
-doc"""
-    ldltfact(::Union{SparseMatrixCSC,Symmetric{Float64,SparseMatrixCSC{Flaot64,SuiteSparse_long}},Hermitian{Complex{Float64},SparseMatrixCSC{Complex{Float64},SuiteSparse_long}}}; shift=0, perm=Int[]) -> CHOLMOD.Factor
-
-Compute the `LDLt` factorization of a sparse symmetric or Hermitian matrix. A fill-reducing permutation is used. `F = ldltfact(A)` is most frequently used to solve systems of equations `A*x = b` with `F\b`, but also the methods `diag`, `det`, `logdet` are defined for `F`. You can also extract individual factors from `F`, using `F[:L]`. However, since pivoting is on by default, the factorization is internally represented as `A == P'*L*D*L'*P` with a permutation matrix `P`; using just `L` without accounting for `P` will give incorrect answers. To include the effects of permutation, it's typically preferable to extact "combined" factors like `PtL = F[:PtL]` (the equivalent of `P'*L`) and `LtP = F[:UP]` (the equivalent of `L'*P`). The complete list of supported factors is `:L, :PtL, :D, :UP, :U, :LD, :DU, :PtLD, :DUP`.
-
-Setting optional `shift` keyword argument computes the factorization of `A+shift*I` instead of `A`. If the `perm` argument is nonempty, it should be a permutation of `1:size(A,1)` giving the ordering to use (instead of CHOLMOD's default AMD ordering).
-
-The function calls the C library CHOLMOD and many other functions from the library are wrapped but not exported.
-"""
-ldltfact(A::SparseMatrixCSC; shift=0, perm=Int[])
-
 doc"""
     connect([host],port) -> TcpSocket
 
@@ -7256,13 +7245,6 @@ Tests whether `A` or its elements are of type `T`.
 """
 iseltype
 
-doc"""
-    symperm(A, p)
-
-Return the symmetric permutation of `A`, which is `A[p,p]`. `A` should be symmetric and sparse, where only the upper triangular part of the matrix is stored. This algorithm ignores the lower triangular part of the matrix. Only the upper triangular part of the result is returned as well.
-"""
-symperm
-
 doc"""
     min(x, y, ...)
 
@@ -9993,15 +9975,6 @@ Like permute!, but the inverse of the given permutation is applied.
 """
 ipermute!
 
-doc"""
-```rst
-..  full(S)
-
-Convert a sparse matrix ``S`` into a dense matrix.
-```
-"""
-full(::AbstractSparseMatrix)
-
 doc"""
 ```rst
 ..  full(F)
@@ -10487,13 +10460,6 @@ k]``.)
 """
 eigfact(A,B)
 
-doc"""
-    rowvals(A)
-
-Return a vector of the row indices of `A`, and any modifications to the returned vector will mutate `A` as well. Given the internal storage format of sparse matrices, providing access to how the row indices are stored internally can be useful in conjuction with iterating over structural nonzero values. See `nonzeros(A)` and `nzrange(A, col)`.
-"""
-rowvals
-
 doc"""
     mkdir(path, [mode])
 

base/exports.jl

@@ -20,16 +20,12 @@ export
     BLAS,
     LAPACK,
     Serializer,
-    SparseMatrix,
     Docs,
     Markdown,
 
 # Types
     AbstractChannel,
     AbstractMatrix,
-    AbstractSparseArray,
-    AbstractSparseMatrix,
-    AbstractSparseVector,
     AbstractVector,
     AbstractVecOrMat,
     Array,
@@ -107,7 +103,6 @@ export
     SharedArray,
     SharedMatrix,
     SharedVector,
-    SparseMatrixCSC,
     StatStruct,
     StepRange,
     StridedArray,
@@ -562,7 +557,6 @@ export
     minimum,
     minmax,
     ndims,
-    nnz,
     nonzeros,
     nthperm!,
     nthperm,
@@ -715,21 +709,7 @@ export
     ×,
 
 # sparse
-    etree,
     full,
-    issparse,
-    sparse,
-    sparsevec,
-    spdiagm,
-    speye,
-    spones,
-    sprand,
-    sprandbool,
-    sprandn,
-    spzeros,
-    symperm,
-    rowvals,
-    nzrange,
 
 # bitarrays
     bitpack,
@@ -1434,4 +1414,27 @@ export
     @assert,
     @enum,
     @label,
-    @goto
+    @goto,
+
+# SparseArrays module re-exports
+    SparseArrays,
+    AbstractSparseArray,
+    AbstractSparseMatrix,
+    AbstractSparseVector,
+    SparseMatrixCSC,
+    SparseVector,
+    etree,
+    issparse,
+    sparse,
+    sparsevec,
+    spdiagm,
+    speye,
+    spones,
+    sprand,
+    sprandbool,
+    sprandn,
+    spzeros,
+    symperm,
+    rowvals,
+    nzrange,
+    nnz

base/irrationals.jl

@@ -122,10 +122,9 @@ const golden = φ
 for T in (Irrational, Rational, Integer, Number)
     ^(::Irrational{:e}, x::T) = exp(x)
 end
-for T in (Range, BitArray, SparseMatrixCSC, StridedArray, AbstractArray)
+for T in (Range, BitArray, StridedArray, AbstractArray)
     .^(::Irrational{:e}, x::T) = exp(x)
 end
-^(::Irrational{:e}, x::AbstractMatrix) = expm(x)
 
 log(::Irrational{:e}) = 1 # use 1 to correctly promote expressions like log(x)/log(e)
 log(::Irrational{:e}, x) = log(x)

base/precompile.jl

@@ -294,7 +294,6 @@ precompile(Base.next, (Dict{Symbol,Any},Int))
 precompile(Base.next, (IntSet, Int))
 precompile(Base.next, (UnitRange{Int},Int))
 precompile(Base.nextind, (ASCIIString, Int))
-precompile(Base.nnz, (BitArray{1},))
 precompile(Base.normpath, (ASCIIString, ASCIIString))
 precompile(Base.normpath, (ASCIIString,))
 precompile(Base.normpath, (UTF8String, UTF8String))

base/sparse.jl

@@ -1,28 +1,38 @@
 # This file is a part of Julia. License is MIT: http://julialang.org/license
 
-module SparseMatrix
+module SparseArrays
 
 using Base: Func, AddFun, OrFun, ConjFun, IdFun
 using Base.Sort: Forward
 using Base.LinAlg: AbstractTriangular, PosDefException
 
 import Base: +, -, *, &, |, $, .+, .-, .*, ./, .\, .^, .<, .!=, ==
-import Base: A_mul_B!, Ac_mul_B, Ac_mul_B!, At_mul_B!, A_ldiv_B!
-import Base: @get!, abs, abs2, broadcast, ceil, complex, cond, conj, convert, copy,
-    ctranspose, diagm, exp, expm1, factorize, find, findmax, findmin, findnz, float,
-    full, getindex, hcat, hvcat, imag, indmax, ishermitian, kron, length, log, log1p,
-    max, min, norm, one, promote_eltype, real, reinterpret, reshape, rot180, rotl90,
-    rotr90, round, scale, scale!, setindex!, similar, size, transpose, tril, triu, vcat,
-    vec
+import Base: A_mul_B!, Ac_mul_B, Ac_mul_B!, At_mul_B, At_mul_B!, A_ldiv_B!
+
+import Base: @get!, acos, acosd, acot, acotd, acsch, asech, asin, asind, asinh,
+    atan, atand, atanh, broadcast!, chol, conj!, cos, cosc, cosd, cosh, cospi, cot,
+    cotd, coth, countnz, csc, cscd, csch, ctranspose!, diag, diff, done, dot, eig,
+    exp10, exp2, eye, findn, floor, hash, indmin, inv, issym, istril, istriu, log10,
+    log2, lu, maxabs, minabs, next, sec, secd, sech, show, showarray, sin, sinc,
+    sind, sinh, sinpi, squeeze, start, sum, sumabs, sumabs2, summary, tan, tand,
+    tanh, trace, transpose!, tril!, triu!, trunc, vecnorm, writemime, abs, abs2,
+    broadcast, call, ceil, complex, cond, conj, convert, copy, ctranspose, diagm,
+    exp, expm1, factorize, find, findmax, findmin, findnz, float, full, getindex,
+    hcat, hvcat, imag, indmax, ishermitian, kron, length, log, log1p, max, min,
+    maximum, minimum, norm, one, promote_eltype, real, reinterpret, reshape, rot180,
+    rotl90, rotr90, round, scale, scale!, setindex!, similar, size, transpose, tril,
+    triu, vcat, vec
+
 import Base.Broadcast: eltype_plus, broadcast_shape
 
-export AbstractSparseArray, AbstractSparseMatrix, AbstractSparseVector, SparseMatrixCSC,
-       blkdiag, dense, droptol!, dropzeros!, etree, issparse, nnz, nonzeros, nzrange,
-       rowvals, sparse, sparsevec, spdiagm, speye, spones, sprand, sprandbool, sprandn,
-       spzeros, symperm
+export AbstractSparseArray, AbstractSparseMatrix, AbstractSparseVector,
+    SparseMatrixCSC, SparseVector, blkdiag, dense, droptol!, dropzeros!, etree,
+    issparse, nonzeros, nzrange, rowvals, sparse, sparsevec, spdiagm, speye, spones,
+    sprand, sprandbool, sprandn, spzeros, symperm, nnz
 
 include("sparse/abstractsparse.jl")
 include("sparse/sparsematrix.jl")
+include("sparse/sparsevector.jl")
 include("sparse/csparse.jl")
 
 include("sparse/linalg.jl")
@@ -32,4 +42,4 @@ if Base.USE_GPL_LIBS
     include("sparse/spqr.jl")
 end
 
-end # module SparseMatrix
+end

base/sparse/cholmod.jl

@@ -9,14 +9,14 @@ import Base.LinAlg: (\), A_mul_Bc, A_mul_Bt, Ac_ldiv_B, Ac_mul_B, At_ldiv_B, At_
                  cholfact, det, diag, ishermitian, isposdef,
                  issym, ldltfact, logdet
 
-import Base.SparseMatrix: sparse, nnz
+importall ..SparseArrays
 
 export
     Dense,
     Factor,
     Sparse
 
-using Base.SparseMatrix: AbstractSparseMatrix, SparseMatrixCSC, increment, indtype
+import ..SparseArrays: AbstractSparseMatrix, SparseMatrixCSC, increment, indtype
 
 #########
 # Setup #
@@ -853,6 +853,9 @@ function convert{Tv<:VTypes}(::Type{Sparse}, A::SparseMatrixCSC{Tv,SuiteSparse_l
 
     return o
 end
+
+# convert SparseVectors into CHOLMOD Sparse types through a mx1 CSC matrix
+convert{Tv<:VTypes}(::Type{Sparse}, A::SparseVector{Tv,SuiteSparse_long}) = convert(Sparse, convert(SparseMatrixCSC, A))
 function convert{Tv<:VTypes}(::Type{Sparse}, A::SparseMatrixCSC{Tv,SuiteSparse_long})
     o = Sparse(A, 0)
     # check if array is symmetric and change stype if it is
@@ -994,7 +997,7 @@ function sparse(F::Factor)
         L, d = getLd!(LD)
         A = scale(L, d)*L'
     end
-    SparseMatrix.sortSparseMatrixCSC!(A)
+    SparseArrays.sortSparseMatrixCSC!(A)
     p = get_perm(F)
     if p != [1:s.n;]
         pinv = Array(Int, length(p))
@@ -1216,6 +1219,32 @@ function cholfact(A::Sparse; kws...)
     return F
 end
 
+doc"""
+    ldltfact(::Union{SparseMatrixCSC,Symmetric{Float64,SparseMatrixCSC{Flaot64,SuiteSparse_long}},Hermitian{Complex{Float64},SparseMatrixCSC{Complex{Float64},SuiteSparse_long}}}; shift=0, perm=Int[]) -> CHOLMOD.Factor
+
+Compute the `LDLt` factorization of a sparse symmetric or Hermitian matrix. A
+fill-reducing permutation is used. `F = ldltfact(A)` is most frequently used to
+solve systems of equations `A*x = b` with `F\b`. The returned factorization
+object `F` also supports the methods `diag`, `det`, and `logdet`. You can
+extract individual factors from `F` using `F[:L]`. However, since pivoting is
+on by default, the factorization is internally represented as `A == P'*L*D*L'*P`
+with a permutation matrix `P`; using just `L` without accounting for `P` will
+give incorrect answers. To include the effects of permutation, it's typically
+preferable to extact "combined" factors like `PtL = F[:PtL]` (the equivalent of
+`P'*L`) and `LtP = F[:UP]` (the equivalent of `L'*P`). The complete list of
+supported factors is `:L, :PtL, :D, :UP, :U, :LD, :DU, :PtLD, :DUP`.
+
+Setting optional `shift` keyword argument computes the factorization of
+`A+shift*I` instead of `A`. If the `perm` argument is nonempty, it should be a
+permutation of `1:size(A,1)` giving the ordering to use (instead of CHOLMOD's
+default AMD ordering).
+
+The function calls the C library CHOLMOD and many other functions from the
+library are wrapped but not exported.
+
+"""
+ldltfact(A::SparseMatrixCSC; shift=0, perm=Int[])
+
 function ldltfact(A::Sparse; kws...)
     cm = defaults(common()) # setting the common struct to default values. Should only be done when creating new factorization.
     set_print_level(cm, 0) # no printing from CHOLMOD by default
@@ -1294,13 +1323,15 @@ for (T, f) in ((:Dense, :solve), (:Sparse, :spsolve))
     end
 end
 
+typealias SparseVecOrMat{Tv,Ti} Union{SparseVector{Tv,Ti}, SparseMatrixCSC{Tv,Ti}}
+
 function (\)(L::FactorComponent, b::Vector)
     reshape(convert(Matrix, L\Dense(b)), length(b))
 end
 function (\)(L::FactorComponent, B::Matrix)
     convert(Matrix, L\Dense(B))
 end
-function (\)(L::FactorComponent, B::SparseMatrixCSC)
+function (\)(L::FactorComponent, B::SparseVecOrMat)
     sparse(L\Sparse(B,0))
 end
 
@@ -1311,12 +1342,12 @@ Ac_ldiv_B(L::FactorComponent, B) = ctranspose(L)\B
 (\)(L::Factor, B::Matrix) = convert(Matrix, solve(CHOLMOD_A, L, Dense(B)))
 (\)(L::Factor, B::Sparse) = spsolve(CHOLMOD_A, L, B)
 # When right hand side is sparse, we have to ensure that the rhs is not marked as symmetric.
-(\)(L::Factor, B::SparseMatrixCSC) = sparse(spsolve(CHOLMOD_A, L, Sparse(B, 0)))
+(\)(L::Factor, B::SparseVecOrMat) = sparse(spsolve(CHOLMOD_A, L, Sparse(B, 0)))
 
 Ac_ldiv_B(L::Factor, B::Dense) = solve(CHOLMOD_A, L, B)
 Ac_ldiv_B(L::Factor, B::VecOrMat) = convert(Matrix, solve(CHOLMOD_A, L, Dense(B)))
 Ac_ldiv_B(L::Factor, B::Sparse) = spsolve(CHOLMOD_A, L, B)
-Ac_ldiv_B(L::Factor, B::SparseMatrixCSC) = Ac_ldiv_B(L, Sparse(B))
+Ac_ldiv_B(L::Factor, B::SparseVecOrMat) = Ac_ldiv_B(L, Sparse(B))
 
 ## Other convenience methods
 function diag{Tv}(F::Factor{Tv})
@@ -1398,7 +1429,7 @@ function ishermitian(A::Sparse{Complex{Float64}})
     end
 end
 
-(*){Ti}(A::Symmetric{Float64,SparseMatrixCSC{Float64,Ti}}, B::SparseMatrixCSC{Float64,Ti}) = sparse(Sparse(A)*Sparse(B))
-(*){Ti}(A::Hermitian{Complex{Float64},SparseMatrixCSC{Complex{Float64},Ti}}, B::SparseMatrixCSC{Complex{Float64},Ti}) = sparse(Sparse(A)*Sparse(B))
+(*){Ti}(A::Symmetric{Float64,SparseMatrixCSC{Float64,Ti}}, B::SparseVecOrMat{Float64,Ti}) = sparse(Sparse(A)*Sparse(B))
+(*){Ti}(A::Hermitian{Complex{Float64},SparseMatrixCSC{Complex{Float64},Ti}}, B::SparseVecOrMat{Complex{Float64},Ti}) = sparse(Sparse(A)*Sparse(B))
 
 end #module

base/sparse/csparse.jl

@@ -313,8 +313,17 @@ function csc_permute{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}, pinv::Vector{Ti}, q::Vect
     (C.').' # double transpose to order the columns
 end
 
+
 # based on cs_symperm p. 21, "Direct Methods for Sparse Linear Systems"
 # form A[p,p] for a symmetric A stored in the upper triangle
+doc"""
+    symperm(A, p)
+
+Return the symmetric permutation of `A`, which is `A[p,p]`. `A` should be
+symmetric, sparse, and only contain nonzeros in the upper triangular part of the
+matrix is stored. This algorithm ignores the lower triangular part of the
+matrix. Only the upper triangular part of the result is returned.
+"""
 function symperm{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}, pinv::Vector{Ti})
     m, n = size(A)
     if m != n

base/sparse/linalg.jl

@@ -160,7 +160,7 @@ function spmatmul{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}, B::SparseMatrixCSC{Tv,Ti};
 
     # The Gustavson algorithm does not guarantee the product to have sorted row indices.
     Cunsorted = SparseMatrixCSC(mA, nB, colptrC, rowvalC, nzvalC)
-    C = Base.SparseMatrix.sortSparseMatrixCSC!(Cunsorted, sortindices=sortindices)
+    C = SparseArrays.sortSparseMatrixCSC!(Cunsorted, sortindices=sortindices)
     return C
 end
 
@@ -752,7 +752,7 @@ inv(A::SparseMatrixCSC) = error("The inverse of a sparse matrix can often be den
 
 ## scale methods
 
-# Copy colptr and rowval from one SparseMatrix to another
+# Copy colptr and rowval from one sparse matrix to another
 function copyinds!(C::SparseMatrixCSC, A::SparseMatrixCSC)
     if C.colptr !== A.colptr
         resize!(C.colptr, length(A.colptr))