add tests and docs

README.md

@@ -6,7 +6,11 @@ Fast tensor operations using a convenient index notation.
 
 ## What's new
 
-- Fully compatible with Julia v0.6. (v0.5 no longer supported)
+- Addition of a `@tensoropt` macro that will optimize the contraction order of any product of tensors `A[...]*B[...]*C[...]*...` (see below).
+- The `@tensor` macro will reorganize the contraction order of the so-called NCON style of specifying indices is respected, i.e. all contracted indices are labelled by positive integers and all uncontracted indices are specified by negative integers. In that case, tensors will be contracted in such an order that indices with smaller integer label will
+be contracted first
+- Better overall type stability, both in the `@tensor(opt)` environment and with the function based approach. Even the simple function based syntax can now be made type stable by specifing the indices using tuples.
+- Fully compatible with Julia v0.6 and v0.7. (v0.5 no longer supported)
 
 ## Installation
 
@@ -53,6 +57,21 @@ In this example, the labels were specified by arbitrary letters or even longer n
 
 The index pattern is analyzed at compile time and wrapped in appropriate types such that the result of the operation can be computed with a minimal number of temporaries. The use of `@generated` functions further enables to move as much of the label analysis to compile time. You can read more about these topics in the section "Implementation" below.
 
+By default, a contraction of several tensors `A[a,b,c,d,e]*B[b,e,f,g]*C[c,f,i,j]*...` will be evaluted using pairwise contractions from left to right, i.e. as `( (A[a,b,c,d,e] * B[b,e,f,g]) * C[c,f,i,j]) * ...`. However, if one respects the so-called [NCON](https://arxiv.org/abs/1402.0939) style of specifying indices, i.e. positive integers for the contracted indices and negative indices for the open indices, the different factors will be reordered and so that the pairwise tensor contractions contract over indices with smaller integer label first. For example, `D[:] := A[-1,3,1,-2,2]*B[3,2,4,-5]*C[1,4,-4,-3]` will be evaluated as `(A[-1,3,1,-2,2]*C[1,4,-4,-3])*B[3,2,4,-5]`. Furthermore, in that case the indices of the output tensor (`D` in this case) do not need to be specified, and will be chosen as `(-1,-2,-3,-4,-5)`. Any other order is of course still possible by just specifying it.
+
+Furthermore, there is a `@tensoropt` macro which will optimize the contraction order to minimize the total number of multiplications (cost model might change or become choosable in the future). The optimal contraction order will be determined at compile time and will be hard coded in the macro expansion. The cost/size of the different indices can be specified in various ways, and can be integers or some arbitrary polynomial of an abstract variable, e.g. `χ`. In the latter case, the optimization assumes the assymptotic limit of large `χ`.
+
+```julia
+@tensoropt C[a,b,c,d] := A[a,e,c,f]*B[f,d,e,b]
+# cost χ for all indices (a,b,c,d,e,f)
+@tensoropt (a,b,c,e) C[a,b,c,d] := A[a,e,c,f]*B[f,d,e,b]
+# cost χ for indices a,b,c,e, other indices (d,f) have cost 1
+@tensoropt !(a,b,c,e) C[a,b,c,d] := A[a,e,c,f]*B[f,d,e,b]
+# cost 1 for indices a,b,c,e, other indices (d,f) have cost χ
+@tensoropt (a=>χ,b=>χ^2,c=>2*χ,e=>5) C[a,b,c,d] := A[a,e,c,f]*B[f,d,e,b]
+# cost as specified for listed indices, unlisted indices have cost 1 (any symbol for χ can be used)
+```
+
 ### Functions
 
 The elementary tensor operations can also be accessed via functions, mainly for compatibility with older versions of this toolbox. The function-based syntax is also required when the contraction pattern is not known at compile time but is rather determined dynamically.
@@ -178,30 +197,29 @@ The `src` folder of `TensorOperations.jl` contains four subfolders in which the
 
 The folders `functions` and `indexnotation` contain the necessary code to support the function-based syntax and the `@tensor` macro syntax respectively. In particular, the folder `indexnotation` contains one file for the macro, and one file for each of the special types `IndexedObject`, `SumOfIndexedObjects` and `ProductOfIndexedObjects` discussed above. Both the function and macro based syntax are completely general and should work for any multidimensional object, provided a corresponding implementation for `add!`, `trace!` and `contract!` is provided.
 
-Finally, the folder `aux` contains some auxiliary code, such as some metaprogramming tools, the `unique2` function which removes all elements of a list which appear more than once, the definition of the `StridedData` type, the elementary `axpby` operation and the definition of a dedicated `IndexError` type for reporting errors in the index notation. Finally, there is a file `stridedarray.jl` which provides some auxiliary abstraction to interface with the `StridedArray` type of Julia Base. It contains the following function definitions:
+Finally, the folder `auxiliary` contains some auxiliary code, such as some metaprogramming tools, the `unique2` function which removes all elements of a list which appear more than once, the definition of the `StridedData` type, the elementary `axpby` operation and the definition of a dedicated `IndexError` type for reporting errors in the index notation. Finally, there is a file `stridedarray.jl` which provides some auxiliary abstraction to interface with the `StridedArray` type of Julia Base. It contains the following function definitions:
 
 - `numind(A)`
 
   Returns the number of indices of a tensor-like object `A`, i.e. for a multidimensional array (`<:AbstractArray`) we have `numind(A) = ndims(A)`. Also works in type domain.
 
-- `similar_from_indices(T, indices, A, conjA=Val{:N})`
+- `similar_from_indices(T, indices::Tuple{Vararg{Int}}, A, conjA=Val{:N})`
 
   Returns an object similar to `A` which has an `eltype` given by `T` and dimensions/sizes corresponding to a selection of those of `op(A)`, where the selection is specified by `indices` (which contains integer between `1` and `numind(A)`) and `op` is `conj` if `conjA=Val{:C}` or does nothing if `conjA=Val{:N}` (default).
 
-- `similar_from_indices(T, indices, A, B, conjA=Val{:N}, conjB={:N})`
+- `similar_from_indices(T, indices::Tuple{Vararg{Int}}, A, B, conjA=Val{:N}, conjB={:N})`
 
   Returns an object similar to `A` which has an `eltype` given by `T` and dimensions/sizes corresponding to a selection of those of `op(A)` and `op(B)` concatenated, where the selection is specified by `indices` (which contains integers between `1` and `numind(A)+numind(B)` and `op` is `conj` if `conjA` or `conjB` equal `Val{:C}` or does nothing if `conjA` or `conjB` equal `Val{:N}` (default).
 
 - `scalar(C)`
 
   Returns the single element of a tensor-like object with zero dimensions, i.e. if `numind(C)==0`.
 
-In summary, to add support for a user-defined tensor-like type, the functions in the files `aux/stridedarray.jl` and `implementation/stridedarray.jl` should be reimplemented. This should be rather straightforwarded if the type just wraps multidimensional data with a strided storage in memory.
+In summary, to add support for a user-defined tensor-like type, the functions in the files `auxiliary/stridedarray.jl` and `implementation/stridedarray.jl` should be reimplemented. This should be rather straightforwarded if the type just wraps multidimensional data with a strided storage in memory.
 
 ## Planned features
 
-The following features seem like interesting additions to the `TensorOperations.jl` package, and might therefore appear in the future (not necessarily in this order).
+The following features seem like interesting additions to the `TensorOperations.jl` package, and might therefore appear in the future.
 
-- Functionality to contract a large set of tensors, also called a tensor network, including a method to optimize over the contraction order along the lines of [arXiv:1304.6112v6](http://arxiv.org/abs/1304.6112v6).
-- Implementation of a tensor contraction implementation that can immediately call the BLAS microkernels to act on small blocks of the arrays which are permuted into fixed size buffers.
+- Implementation of a tensor contraction implementation that can immediately call the BLAS microkernels to act on small blocks of the arrays which are permuted into fixed size buffers. Or rather, possibility to use some of the recent efforts in this direction as backend (e.g. [TBLIS](https://github.com/devinamatthews/tblis), [HPTT](https://github.com/springer13/hptt) and [TCL](https://github.com/springer13/tcl)).
 - Implementation of a `@rawindex` macro that translates an expression with index notation directly into a raw set of nested loops. While loosing the advantage of cache-friendlyness for large arrays, this can be advantageous for smaller arrays. In addition, this would allow for more general expressions where functions of array entries are computed without having to allocate a temporary array, or where three or more dimensions are simultaneously summed over.

test/methods.jl

@@ -1,51 +1,49 @@
 # test simple methods
 #---------------------
 # test tensorcopy
-A=randn((3,5,4,6))
-p=randperm(4)
-C1=permutedims(A,p)
-C2=tensorcopy(A,(1:4...),(p...))
-@inferred tensorcopy(A,(1:4...),(p...))
+A = randn((3,5,4,6))
+p = randperm(4)
+C1 = permutedims(A,p)
+C2 = @inferred tensorcopy(A,(1:4...),(p...))
 @test vecnorm(C1-C2)<eps()*sqrt(length(C1))*vecnorm(C1+C2)
 @test_throws TensorOperations.IndexError tensorcopy(A,1:3,1:4)
 @test_throws TensorOperations.IndexError tensorcopy(A,[1,2,2,4],1:4)
 
 # test tensoradd
-B=randn((5,6,3,4))
-p=(3,1,4,2)
-C1=tensoradd(A,p,B,(1:4...))
-@inferred tensoradd(A,p,B,(1:4...))
-C2=A+permutedims(B,p)
+B = randn((5,6,3,4))
+p = (3,1,4,2)
+C1 = @inferred tensoradd(A,p,B,(1:4...))
+C2 = A+permutedims(B,p)
 @test vecnorm(C1-C2)<eps()*sqrt(length(C1))*vecnorm(C1+C2)
 @test_throws DimensionMismatch tensoradd(A,1:4,B,1:4)
 
 # test tensortrace
-A=randn(50,100,100)
-C1=tensortrace(A,[:a,:b,:b])
-C2=zeros(50)
-for i=1:50
-    for j=1:100
-        C2[i]+=A[i,j,j]
+A = randn(50,100,100)
+C1 = tensortrace(A,[:a,:b,:b])
+C2 = zeros(50)
+for i = 1:50
+    for j = 1:100
+        C2[i] += A[i,j,j]
     end
 end
 @test vecnorm(C1-C2)<eps()*sqrt(length(C1))*vecnorm(C1+C2)
-A=randn(3,20,5,3,20,4,5)
-C1=tensortrace(A,[:a,:b,:c,:d,:b,:e,:c],[:e,:a,:d])
-C2=zeros(4,3,3)
-for i1=1:4, i2=1:3, i3=1:3
-    for j1=1:20,j2=1:5
-        C2[i1,i2,i3]+=A[i2,j1,j2,i3,j1,i1,j2]
+A = randn(3,20,5,3,20,4,5)
+C1 = @inferred tensortrace(A,(:a,:b,:c,:d,:b,:e,:c),(:e,:a,:d))
+C2 = zeros(4,3,3)
+for i1 = 1:4, i2 = 1:3, i3 = 1:3
+    for j1 = 1:20,j2 = 1:5
+        C2[i1,i2,i3] += A[i2,j1,j2,i3,j1,i1,j2]
     end
 end
 @test vecnorm(C1-C2)<eps()*sqrt(length(C1))*vecnorm(C1+C2)
 
 # test tensorcontract
-A=randn(3,20,5,3,4)
-B=randn(5,6,20,3)
-C1=tensorcontract(A,[:a,:b,:c,:d,:e],B,[:c,:f,:b,:g],[:a,:g,:e,:d,:f];method=:BLAS)
-C2=tensorcontract(A,[:a,:b,:c,:d,:e],B,[:c,:f,:b,:g],[:a,:g,:e,:d,:f];method=:native)
-C3=zeros(3,3,4,3,6)
-for a=1:3, b=1:20, c=1:5, d=1:3, e=1:4, f=1:6, g=1:3
+A = randn(3,20,5,3,4)
+B = randn(5,6,20,3)
+C1 = @inferred tensorcontract(A,(:a,:b,:c,:d,:e),B,(:c,:f,:b,:g),(:a,:g,:e,:d,:f);method = :BLAS)
+C2 = @inferred tensorcontract(A,(:a,:b,:c,:d,:e),B,(:c,:f,:b,:g),(:a,:g,:e,:d,:f);method = :native)
+C3 = zeros(3,3,4,3,6)
+for a = 1:3, b = 1:20, c = 1:5, d = 1:3, e = 1:4, f = 1:6, g = 1:3
     C3[a,g,e,d,f] += A[a,b,c,d,e]*B[c,f,b,g]
 end
 @test vecnorm(C1-C3)<eps()*sqrt(length(C1))*vecnorm(C1+C3)
@@ -54,29 +52,29 @@ end
 @test_throws TensorOperations.IndexError tensorcontract(A,[:a,:b,:c,:a,:e],B,[:c,:f,:b,:g])
 
 # test tensorproduct
-A=randn(5,5,5,5)
-B=rand(Complex128,(5,5,5,5))
-C1=reshape(tensorproduct(A,[1,2,3,4],B,[5,6,7,8],[1,2,5,6,3,4,7,8]),(5*5*5*5,5*5*5*5))
-C2=kron(reshape(B,(25,25)),reshape(A,(25,25)))
+A = randn(5,5,5,5)
+B = rand(Complex128,(5,5,5,5))
+C1 = reshape(@inferred tensorproduct(A,(1,2,3,4),B,(5,6,7,8),(1,2,5,6,3,4,7,8)),(5*5*5*5,5*5*5*5))
+C2 = kron(reshape(B,(25,25)),reshape(A,(25,25)))
 @test vecnorm(C1-C2)<eps()*sqrt(length(C1))*vecnorm(C1+C2)
 @test_throws TensorOperations.IndexError tensorproduct(A,[:a,:b,:c,:d],B,[:d,:e,:f,:g])
 @test_throws TensorOperations.IndexError tensorproduct(A,[:a,:b,:c,:d],B,[:e,:f,:g,:h],[:a,:b,:c,:d,:e,:f,:g,:i])
 
 # test index notation
 #---------------------
-# Da=10
-# Db=15
-# Dc=4
-# Dd=8
-# De=6
-# Df=7
-# Dg=3
-# Dh=2
-# A=rand(Complex128,(Da,Dc,Df,Da,De,Db,Db,Dg))
-# B=rand(Complex128,(Dc,Dh,Dg,De,Dd))
-# C=rand(Complex128,(Dd,Dh,Df))
-# D=rand(Complex128,(Dd,Df,Dh))
-# D[l"d,f,h"]=A[l"a,c,f,a,e,b,b,g"]*B[l"c,h,g,e,d"]+0.5*C[l"d,h,f"]
+# Da = 10
+# Db = 15
+# Dc = 4
+# Dd = 8
+# De = 6
+# Df = 7
+# Dg = 3
+# Dh = 2
+# A = rand(Complex128,(Da,Dc,Df,Da,De,Db,Db,Dg))
+# B = rand(Complex128,(Dc,Dh,Dg,De,Dd))
+# C = rand(Complex128,(Dd,Dh,Df))
+# D = rand(Complex128,(Dd,Df,Dh))
+# D[l"d,f,h"] = A[l"a,c,f,a,e,b,b,g"]*B[l"c,h,g,e,d"]+0.5*C[l"d,h,f"]
 # @test_approx_eq(vecnorm(D),sqrt(abs(scalar(D[l"d,f,h"]*conj(D[l"d,f,h"])))))
 # @test_throws IndexError D[l"a,a,a"]
 # @test_throws IndexError D[l"a,b,c,d"]
@@ -88,13 +86,13 @@ C2=kron(reshape(B,(25,25)),reshape(A,(25,25)))
 # with changing element type and with nontrivial strides
 
 # tensorcopy!
-Abig=randn((30,30,30,30))
-A=view(Abig,1+3*(0:9),2+2*(0:6),5+4*(0:6),4+3*(0:8))
-p=[3,1,4,2]
-Cbig=zeros(Complex128,(50,50,50,50))
-C=view(Cbig,13+(0:6),11+4*(0:9),15+4*(0:8),4+3*(0:6))
-Acopy=tensorcopy(A,1:4,1:4)
-Ccopy=tensorcopy(C,1:4,1:4)
+Abig = randn((30,30,30,30))
+A = view(Abig,1+3*(0:9),2+2*(0:6),5+4*(0:6),4+3*(0:8))
+p = [3,1,4,2]
+Cbig = zeros(Complex128,(50,50,50,50))
+C = view(Cbig,13+(0:6),11+4*(0:9),15+4*(0:8),4+3*(0:6))
+Acopy = tensorcopy(A,1:4,1:4)
+Ccopy = tensorcopy(C,1:4,1:4)
 TensorOperations.tensorcopy!(A,1:4,C,p)
 TensorOperations.tensorcopy!(Acopy,1:4,Ccopy,p)
 @test vecnorm(C-Ccopy)<eps()*sqrt(length(C))*vecnorm(C+Ccopy)
@@ -103,32 +101,32 @@ TensorOperations.tensorcopy!(Acopy,1:4,Ccopy,p)
 @test_throws TensorOperations.IndexError TensorOperations.tensorcopy!(A,1:4,C,[1,1,2,3])
 
 # tensoradd!
-Cbig=zeros(Complex128,(50,50,50,50))
-C=view(Cbig,13+(0:6),11+4*(0:9),15+4*(0:8),4+3*(0:6))
-Acopy=tensorcopy(A,1:4,p)
-Ccopy=tensorcopy(C,1:4,1:4)
-alpha=randn()
-beta=randn()
+Cbig = zeros(Complex128,(50,50,50,50))
+C = view(Cbig,13+(0:6),11+4*(0:9),15+4*(0:8),4+3*(0:6))
+Acopy = tensorcopy(A,1:4,p)
+Ccopy = tensorcopy(C,1:4,1:4)
+alpha = randn()
+beta = randn()
 TensorOperations.tensoradd!(alpha,A,1:4,beta,C,p)
-Ccopy=beta*Ccopy+alpha*Acopy
+Ccopy = beta*Ccopy+alpha*Acopy
 @test vecnorm(C-Ccopy)<eps()*sqrt(length(C))*vecnorm(C+Ccopy)
 @test_throws TensorOperations.IndexError TensorOperations.tensoradd!(1.2,A,1:3,0.5,C,p)
 @test_throws DimensionMismatch TensorOperations.tensoradd!(1.2,A,p,0.5,C,p)
 @test_throws TensorOperations.IndexError TensorOperations.tensoradd!(1.2,A,1:4,0.5,C,[1,1,2,3])
 
 # tensortrace!
-Abig=rand((30,30,30,30))
-A=view(Abig,1+3*(0:8),2+2*(0:14),5+4*(0:6),7+2*(0:8))
-Bbig=rand(Complex128,(50,50))
-B=view(Bbig,13+(0:14),3+5*(0:6))
-Acopy=tensorcopy(A,1:4)
-Bcopy=tensorcopy(B,1:2)
-alpha=randn()
-beta=randn()
+Abig = rand((30,30,30,30))
+A = view(Abig,1+3*(0:8),2+2*(0:14),5+4*(0:6),7+2*(0:8))
+Bbig = rand(Complex128,(50,50))
+B = view(Bbig,13+(0:14),3+5*(0:6))
+Acopy = tensorcopy(A,1:4)
+Bcopy = tensorcopy(B,1:2)
+alpha = randn()
+beta = randn()
 TensorOperations.tensortrace!(alpha,A,[:a,:b,:c,:a],beta,B,[:b,:c])
-Bcopy=beta*Bcopy
-for i=1+(0:8)
-    Bcopy+=alpha*view(A,i,:,:,i)
+Bcopy = beta*Bcopy
+for i = 1+(0:8)
+    Bcopy += alpha*view(A,i,:,:,i)
 end
 @test vecnorm(B-Bcopy)<eps()*vecnorm(B+Bcopy)*sqrt(length(B))
 @test_throws TensorOperations.IndexError TensorOperations.tensortrace!(alpha,A,[:a,:b,:c],beta,B,[:b,:c])
@@ -137,35 +135,35 @@ end
 @test_throws DimensionMismatch TensorOperations.tensortrace!(alpha,A,[:a,:b,:a,:c],beta,B,[:c,:b])
 
 # tensorcontract!
-Abig=rand((30,30,30,30))
-A=view(Abig,1+3*(0:8),2+2*(0:14),5+4*(0:6),7+2*(0:8))
-Bbig=rand(Complex128,(50,50,50))
-B=view(Bbig,3+5*(0:6),7+2*(0:7),13+(0:14))
-Cbig=rand(Complex64,(40,40,40))
-C=view(Cbig,3+2*(0:8),13+(0:8),7+3*(0:7))
-Acopy=tensorcopy(A,1:4)
-Bcopy=tensorcopy(B,1:3)
-Ccopy=tensorcopy(C,1:3)
-alpha=randn()
-beta=randn()
-Ccopy=beta*Ccopy
-for d=1+(0:8),a=1+(0:8),e=1+(0:7)
-    for b=1+(0:14),c=1+(0:6)
-        Ccopy[d,a,e]+=alpha*A[a,b,c,d]*conj(B[c,e,b])
+Abig = rand((30,30,30,30))
+A = view(Abig,1+3*(0:8),2+2*(0:14),5+4*(0:6),7+2*(0:8))
+Bbig = rand(Complex128,(50,50,50))
+B = view(Bbig,3+5*(0:6),7+2*(0:7),13+(0:14))
+Cbig = rand(Complex64,(40,40,40))
+C = view(Cbig,3+2*(0:8),13+(0:8),7+3*(0:7))
+Acopy = tensorcopy(A,1:4)
+Bcopy = tensorcopy(B,1:3)
+Ccopy = tensorcopy(C,1:3)
+alpha = randn()
+beta = randn()
+Ccopy = beta*Ccopy
+for d = 1+(0:8),a = 1+(0:8),e = 1+(0:7)
+    for b = 1+(0:14),c = 1+(0:6)
+        Ccopy[d,a,e] += alpha*A[a,b,c,d]*conj(B[c,e,b])
     end
 end
-TensorOperations.tensorcontract!(alpha,A,[:a,:b,:c,:d],'N',B,[:c,:e,:b],'C',beta,C,[:d,:a,:e];method=:BLAS)
+TensorOperations.tensorcontract!(alpha,A,[:a,:b,:c,:d],'N',B,[:c,:e,:b],'C',beta,C,[:d,:a,:e];method = :BLAS)
 @test vecnorm(C-Ccopy)<eps(Float32)*vecnorm(C+Ccopy)*sqrt(length(C))
-Cbig=rand(Complex64,(40,40,40))
-C=view(Cbig,3+2*(0:8),13+(0:8),7+3*(0:7))
-Ccopy=tensorcopy(C,1:3)
-Ccopy=beta*Ccopy
-for d=1+(0:8),a=1+(0:8),e=1+(0:7)
-    for b=1+(0:14),c=1+(0:6)
-        Ccopy[d,a,e]+=alpha*A[a,b,c,d]*conj(B[c,e,b])
+Cbig = rand(Complex64,(40,40,40))
+C = view(Cbig,3+2*(0:8),13+(0:8),7+3*(0:7))
+Ccopy = tensorcopy(C,1:3)
+Ccopy = beta*Ccopy
+for d = 1+(0:8),a = 1+(0:8),e = 1+(0:7)
+    for b = 1+(0:14),c = 1+(0:6)
+        Ccopy[d,a,e] += alpha*A[a,b,c,d]*conj(B[c,e,b])
     end
 end
-TensorOperations.tensorcontract!(alpha,A,[:a,:b,:c,:d],'N',B,[:c,:e,:b],'C',beta,C,[:d,:a,:e];method=:native)
+TensorOperations.tensorcontract!(alpha,A,[:a,:b,:c,:d],'N',B,[:c,:e,:b],'C',beta,C,[:d,:a,:e];method = :native)
 @test vecnorm(C-Ccopy)<eps(Float32)*vecnorm(C+Ccopy)*sqrt(length(C))
 @test_throws TensorOperations.IndexError TensorOperations.tensorcontract!(alpha,A,[:a,:b,:c,:a],'N',B,[:c,:e,:b],'N',beta,C,[:d,:a,:e])
 @test_throws TensorOperations.IndexError TensorOperations.tensorcontract!(alpha,A,[:a,:b,:c,:d],'N',B,[:c,:b],'N',beta,C,[:d,:a,:e])