Add Array function

src/GCPDecompositions.jl

@@ -6,6 +6,7 @@ module GCPDecompositions
 # Imports
 import Base: ndims, size, show, summary
 import Base: getindex
+import Base: AbstractArray, Array
 import LinearAlgebra: norm
 using IntervalSets: Interval
 using Random: default_rng

src/cpd.jl

@@ -87,6 +87,9 @@ function getindex(M::CPD{T,N}, I::Vararg{Int,N}) where {T,N}
 end
 getindex(M::CPD{T,N}, I::CartesianIndex{N}) where {T,N} = getindex(M, Tuple(I)...)
 
+AbstractArray(A::CPD) = reshape(TensorKernels.khatrirao(reverse(A.U)...) * A.λ, size(A))
+Array(A::CPD) = Array(AbstractArray(A))
+
 norm(M::CPD, p::Real = 2) =
     p == 2 ? norm2(M) : norm((M[I] for I in CartesianIndices(size(M))), p)
 function norm2(M::CPD{T,N}) where {T,N}

test/items/cpd.jl

@@ -122,6 +122,20 @@ end
     end
 end
 
+@testitem "Array" begin
+    @testset "N=$N, K=$K" for N in 1:3, K in 1:3
+        T = Float64
+        λfull = T[1, 100, 10000]
+        U1full, U2full, U3full = T[1 2 3; 4 5 6], T[-1 0 1], T[1 2 3; 4 5 6; 7 8 9]
+        λ = λfull[1:K]
+        U = (U1full[:, 1:K], U2full[:, 1:K], U3full[:, 1:K])[1:N]
+        M = CPD(λ, U)
+
+        X = Array(M)
+        @test all(I -> M[I] == X[I], CartesianIndices(X))
+    end
+end
+
 @testitem "norm" begin
     using LinearAlgebra