diff --git a/.gitignore b/.gitignore index 8c960ec..3f02ca7 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,4 @@ *.jl.cov *.jl.*.cov *.jl.mem +Manifest.toml diff --git a/.travis.yml b/.travis.yml index 96805b8..a15f4c7 100644 --- a/.travis.yml +++ b/.travis.yml @@ -1,36 +1,6 @@ -## Documentation: http://docs.travis-ci.com/user/languages/julia/ language: julia -os: - - linux - - osx julia: - - 1.0 - - 1.1 - - nightly -notifications: - email: false -git: - depth: 99999999 + - 1 -## uncomment the following lines to allow failures on nightly julia -## (tests will run but not make your overall status red) -matrix: - allow_failures: - - julia: nightly - -## uncomment and modify the following lines to manually install system packages -#addons: -# apt: # apt-get for linux -# packages: -# - gfortran -#before_script: # homebrew for mac -# - if [ $TRAVIS_OS_NAME = osx ]; then brew install gcc; fi - -## uncomment the following lines to override the default test script -#script: -# - julia -e 'Pkg.clone(pwd()); Pkg.build("TensorToolbox"); Pkg.test("TensorToolbox"; coverage=true)' -after_success: - # push coverage results to Coveralls - - julia -e 'cd(Pkg.dir("TensorToolbox")); Pkg.add("Coverage"); using Coverage; Coveralls.submit(Coveralls.process_folder())' - # push coverage results to Codecov - - julia -e 'cd(Pkg.dir("TensorToolbox")); Pkg.add("Coverage"); using Coverage; Codecov.submit(Codecov.process_folder())' +env: + - JULIA_NUM_THREADS=2 diff --git a/Project.toml b/Project.toml index 5487a8e..2487e1c 100644 --- a/Project.toml +++ b/Project.toml @@ -1,10 +1,12 @@ name = "TensorToolbox" uuid = "9c690861-8ade-587a-897e-15364bc6f718" -version = "1.0.1" author = ["Lana Periša lana.perisa@fesb.hr"] +version = "1.0.1" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] diff --git a/src/TensorToolbox.jl b/src/TensorToolbox.jl index 3ffd1a4..a784301 100644 --- a/src/TensorToolbox.jl +++ b/src/TensorToolbox.jl @@ -13,5 +13,6 @@ include("ktensor.jl") include("dimtree.jl") include("htensor.jl") include("TTtensor.jl") +include("sptensor.jl") end #module diff --git a/src/ktensor.jl b/src/ktensor.jl index 71f338b..95ea55e 100644 --- a/src/ktensor.jl +++ b/src/ktensor.jl @@ -258,7 +258,7 @@ function innerprod(X1::ktensor,X2::ktensor) sum(inpr[:]) end function innerprod(X1::ktensor,X2) - @assert(isa(X2,ttensor) || isa(X2,Array), "Inner product not available for type $(typeof(X2)).") + #@assert(isa(X2,ttensor) || isa(X2,Array), "Inner product not available for type $(typeof(X2)).") @assert(size(X1)==size(X2),"Dimension mismatch.") inpr=0 for r=1:ncomponents(X1) diff --git a/src/sptensor.jl b/src/sptensor.jl new file mode 100644 index 0000000..f47359f --- /dev/null +++ b/src/sptensor.jl @@ -0,0 +1,202 @@ +# Inspired by: +# +# BW Bader and TG Kolda. Efficient MATLAB Computations with Sparse +# and Factored Tensors, SIAM J Scientific Computing 30:205-231, 2007. +# DOI: 10.1137/060676489. [BibTeX] + +export SparseTensor, ttt, sptenrand + +import SparseArrays + +mutable struct SparseTensor{T,N} <: AbstractArray{T,N} + dict::Dict{Dims{N},T} + # dict::Dict{Tuple{Int,Vararg{Int}},T} # Julia doesn't like this + dims::Dims{N} + function SparseTensor(v::Array{E,1},subs::Array{Int,2},dims="auto") where E <: Number + newsubs = unique(subs,dims=1) + + # This feels like it should be slow + dict = Dict(s => zero(eltype(v)) for s in Tuple.(eachrow(newsubs))) + for i in 1:length(v) + dict[Tuple(subs[i,:])] += v[i] + end + + #dims = CartesianIndex(maximum(newsubs, dims=1)...) + SparseTensor(dict,dims == "auto" ? maximum(subs,dims=1) |> Tuple : dims) + end + SparseTensor(d::Dict,dims) = new{d |> values |> eltype, length(dims)}(d,dims) + SparseTensor(d::Dict) = SparseTensor(d,maximum(_indarray(d),dims=1) |> Tuple) + SparseTensor(a::AbstractArray) = SparseTensor(Dict(Tuple(ind) => a[ind] for ind in findall(x->x!==zero(x),a)),size(a)) +end + +function _indarray(d::Dict) + reduce(hcat,d |> keys .|> collect) |> permutedims +end +_indarray(t::SparseTensor) = _indarray(t.dict) + +Base.size(t::SparseTensor) = t.dims + +# TODO: add error checking - out of bounds etc. +function Base.getindex(A::SparseTensor{T,N}, inds::Vararg{Int,N}) where {N,T} + get(A.dict,inds,zero(eltype(A))) +end +# Bug: Julia asks for (1,1) indices for rank-1 tensors - this is a workaround but not a pretty one +# and this doesn't seem to do anything +# Base.IndexStyle(::SparseTensor{T,1}) where T = IndexLinear() +function Base.getindex(A::SparseTensor{T,1}, inds::Vararg{Int,N}) where {N,T} + get(A.dict,(inds[1],),zero(eltype(A))) +end + +function Base.setindex!(A::SparseTensor, value, inds::Vararg{Int,N}) where N + A.dict[inds] = value +end + +# If you need exactly n non-zeroes, provide subrng: n,dims -> Array{Int,2} +# E.g. +# ``` +# subrng=(n,dims)->begin +# reduce(hcat, StatsBase.sample( +# Iterators.product([range(1;length=e) for e in dims]...) |> collect, +# n; +# replace=false +# ) .|> collect) |> permutedims +# end +# ``` (NB: very slow for large numbers of dimensions - ideally we want a sampler that can operate on an iterator) +function sptenrand(dims::Array{Int,1},n::Int;rng=rand,subrng=_subrng) + subs = subrng(n,dims) + vals = rng(size(subs)[1]) + SparseTensor(vals,subs,dims|>Tuple) +end + +# Generates n indices into dims and discards repeats +_subrng(n,dims::Array{Int,1}) = unique(round.((rand(Float64,n,length(dims))) .* (dims' .-1)) .+ 1 .|> Int, dims=1) + +sptenrand(dims,d::AbstractFloat;kwargs...) = sptenrand(dims,d*prod(dims)|>round|>Int;kwargs...) + +# Broadcasting - inspired by https://github.com/andyferris/Dictionaries.jl/blob/9b22a254683260354fda8e32291727dfad040106/src/broadcast.jl#L94 +# See also: https://docs.julialang.org/en/v1/manual/interfaces/index.html +Base.Broadcast.broadcastable(t::SparseTensor) = t + +struct SparseTensorBroadcastStyle{N} <: Base.Broadcast.AbstractArrayStyle{N} +end + +Base.Broadcast.BroadcastStyle(::Type{<:SparseTensor{T,N}}) where {T,N} = SparseTensorBroadcastStyle{N}() + +# Fall back to dense array if one of the broadcast arguments is a dense array +Base.Broadcast.BroadcastStyle(::SparseTensorBroadcastStyle{D}, a::Base.Broadcast.DefaultArrayStyle{N}) where {D,N} = Base.Broadcast.DefaultArrayStyle{max(N,D)}() + +# TODO: +# - sparse .* dense = sparse +# - sparse .+ dense = dense, but it can be sped up by being s .+ d = begin c = copy(d); for (k,v) in s.dict; c[k...] .+ v; end +# (need to make own SparseDense broadcast style to replace the DefaultArrayStyle above and then make the function below specialise on f={Base.:*,Base.:+} + +function Base.Broadcast.broadcasted(::SparseTensorBroadcastStyle, f, args...) + # max() currently fails for (2,2,2),(3,1,3) - should return (3,2,3) or maybe just an error + SparseTensor(merge_default(f,getfield.(args,:dict)...),max(getfield.(args,:dims)...)) +end + +" Merge two dictionaries such that ans[k] = f((get(l,k,0),get(r,k,0)) for all k in keys(l) ∪ keys(r)" +function merge_default(f,l,r) + merged = promote_type(typeof(l),typeof(r))() + for k in keys(l) + merged[k] = f(l[k],get(r,k,zero(valtype(r)))) + end + for k in keys(r) + haskey(merged, k) && continue + merged[k] = f(zero(valtype(l)),r[k]) + end + merged +end + +# I'm not sure what these do but they're not good +# Base.Broadcast.materialize(t::SparseTensor) = copy(t) +# Base.Broadcast.materialize!(out::SparseTensor, t::SparseTensor) = copyto!(out, t) + + +# Base.:* - see ttt.m +# Outer product: (l⊗r)_{i_1,...,i_n,j_1,...,j_n} = l_{i_1,...,i_n} r_{j_1,...,j_n} (i.e. not the Kronecker product) +function ttt(l::SparseTensor, r::SparseTensor) + SparseTensor(Dict((a..., b...) => l[a...] * r[b...] for (a, b) in Iterators.product(l.dict |> keys,r.dict |> keys)),(l.dims...,r.dims...)) +end + +function Base.Array(t::SparseTensor) + a = zeros(size(t)) + for (k,v) in t.dict + a[k...] = v + end + a +end + +# To implement for ttv: +# - permutedims +# - reshape +# - sparse-vector multiplication +# function ttv2(X::AbstractArray{<:Number,N},V::VectorCell,modes::AbstractVector{<:Integer}) where N +# remmodes=setdiff(1:N,modes)' +# if N > 1 +# X=permutedims(X,[remmodes modes']) +# @show [remmodes modes'] +# @show X +# end +# sz=size(X) +# @show sz +# if length(modes) < length(V) +# V=V[modes] +# end +# M=N +# for n=length(modes):-1:1 +# X=reshape(X,prod(sz[1:M-1]),sz[M]) +# @show X +# X=X*V[n] +# @show X +# M-=1 +# end +# if M>0 +# X=reshape(X,sz[1:M]) +# end +# X +# end + +# TODO: +# - consider automatically backing off to CSC for all SparseTensor{T,2} +function SparseArrays.sparse(t::SparseTensor{T,2}) where T + i = _indarray(t) + sparse(i[:,1],i[:,2],t.dict|>values|>collect,size(t)...) +end + +Base.:*(t::SparseTensor{T,2},v::AbstractArray{T2,1}) where {T,T2} = begin + SparseTensor(SparseArrays.sparse(t)*v) +end + +Base.promote_rule(::Type{SparseTensor{T1,D}},::Type{Array{T2,D}}) where {T1,T2,D} = SparseTensor{promote_type(T1,T2),D} + +Base.permutedims(t::SparseTensor,dims) = SparseTensor(Dict(k[vec(dims)] => v for (k,v) in t.dict), t.dims) + +# heart of ttv for s(2,2,2) and v(2) is thus: +# julia> reshape(reshape(permutedims(s, [2 3 1]),2*2,2)*v,2,2) +# +# which uses dense arrays. Reshape has a special type itself. +# +# TS = promote_op(matprod, T, S) # where matprod = x*y + x*y +# julia> Base.promote_op(+,SparseTensor{Int,2},Array{Int,2}) +# Array{Int64,2} +# +# Even though we have the promotion rule above. It uses julia> Core.Compiler.return_type(+,Tuple{T1,T2}) internally. +# Ideally we should be able to avoid it altoghter. + +mutable struct ReshapedSparseTensor{T,D} # <: AbstractSparseTensor + # Will need custom get/set index + underlying::SparseTensor{T,D} + dims::Dims{D} +end + +function Base.reshape(t::SparseTensor,dims::Dims) +end + +# TODO: Support slices / Colon() + +# TODO: stop implementing these myself. I should just write a TensorToolbox.jl compatible type and then optimise when I can be bothered. _headdesk_ +# (bonus: it would let me check my maths) diff --git a/test/runtests.jl b/test/runtests.jl index 4d37ca7..1714191 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -10,5 +10,6 @@ include("ktensorTest.jl") include("dimtreeTest.jl") include("htensorTest.jl") include("TTtensorTest.jl") +include("sptensorTest.jl") println("\n\n") diff --git a/test/sptensorTest.jl b/test/sptensorTest.jl new file mode 100644 index 0000000..734a7ce --- /dev/null +++ b/test/sptensorTest.jl @@ -0,0 +1,16 @@ +import Random +Random.seed!(1337) + +println("\n") + +sp = sptenrand([10, 10, 10], 0.01) +dense = Array(sp) + +println("... Testing sparse and dense tensor product equivalence.") +@test ttt(sp,sp) == ttt(dense,dense) + +println("... Testing sparse tensor decomposition") +@test cp_als(reduce(ttt,SparseTensor.([[0, 0, 1],[0, 1, 0],[1, 0, 0]])),1) == ktensor([[0 0 1],[0 1 0],[1 0 0]].|>permutedims .|> x->x .|> Float64) + +println("... Testing sparse tensor broadcast") +@test SparseTensor([0 0 1]) .* SparseTensor([1 0 0]) == [0 0 0]