more fixes and first tests

Author: Jutho

Project.toml

@@ -40,7 +40,7 @@ TensorOperations = "5.1"
 Test = "1"
 TestExtras = "0.2,0.3"
 TupleTools = "1.1"
-VectorInterface = "0.4"
+VectorInterface = "0.4,0.5"
 julia = "1.10"
 
 [extras]

src/TensorKit.jl

@@ -31,7 +31,7 @@ export TruncationScheme
 export SpaceMismatch, SectorMismatch, IndexError # error types
 
 # general vector space methods
-export space, field, dual, dim, dims, fuse, flip, isdual, insertunit, oplus
+export space, field, dual, dim, reduceddim, dims, fuse, flip, isdual, insertunit, oplus
 
 # partial order for vector spaces
 export infimum, supremum, isisomorphic, ismonomorphic, isepimorphic

src/spaces/gradedspace.jl

@@ -34,7 +34,9 @@ function GradedSpace{I,NTuple{N,Int}}(dims; dual::Bool=false) where {I,N}
     d = ntuple(n -> 0, N)
     isset = ntuple(n -> false, N)
     for (c, dc) in dims
-        i = findindex(values(I), convert(I, c))
+        k = convert(I, c)
+        i = findindex(values(I), k)
+        k = dc < 0 && throw(ArgumentError("Sector $k has negative dimension $dc"))
         isset[i] && throw(ArgumentError("Sector $c appears multiple times"))
         isset = TupleTools.setindex(isset, true, i)
         d = TupleTools.setindex(d, dc, i)
@@ -50,6 +52,7 @@ function GradedSpace{I,SectorDict{I,Int}}(dims; dual::Bool=false) where {I<:Sect
     for (c, dc) in dims
         k = convert(I, c)
         haskey(d, k) && throw(ArgumentError("Sector $k appears multiple times"))
+        dc < 0 && throw(ArgumentError("Sector $k has negative dimension $dc"))
         !iszero(dc) && push!(d, k => dc)
     end
     return GradedSpace{I,SectorDict{I,Int}}(d, dual)

src/spaces/vectorspaces.jl

@@ -111,6 +111,13 @@ field(V::ElementarySpace) = field(typeof(V))
 Return the degeneracy dimension corresponding to the sector `s` of the vector space `V`.
 """ dim(::ElementarySpace, ::Sector)
 
+@doc """
+    reduceddim(V::ElementarySpace) -> Int
+
+Return the sum of all degeneracy dimensions of the vector space `V`.
+"""
+reduceddim(V::ElementarySpace) = sum(c -> dim(V, c), sectors(V); init=0)
+
 """
     oneunit(V::S) where {S<:ElementarySpace} -> S
 

src/tensors/diagtensor.jl

@@ -8,6 +8,9 @@ struct DiagonalTensorMap{T,S<:IndexSpace,A<:DenseVector{T}} <: AbstractTensorMap
     function DiagonalTensorMap{T,S,A}(::UndefInitializer,
                                       dom::S) where {T,S<:IndexSpace,A<:DenseVector{T}}
         data = A(undef, reduceddim(dom))
+        if !isbitstype(T)
+            zerovector!(data)
+        end
         return DiagonalTensorMap{T,S,A}(data, dom)
     end
     # constructors from data
@@ -17,7 +20,6 @@ struct DiagonalTensorMap{T,S<:IndexSpace,A<:DenseVector{T}} <: AbstractTensorMap
         return new{T,S,A}(data, dom)
     end
 end
-reduceddim(V::IndexSpace) = sum(c -> dim(V, c), sectors(V); init=0)
 
 # Basic methods for characterising a tensor:
 #--------------------------------------------
@@ -43,6 +45,7 @@ Construct a `DiagonalTensorMap` with uninitialized data.
 function DiagonalTensorMap{T}(::UndefInitializer, V::S) where {T,S<:IndexSpace}
     return DiagonalTensorMap{T,S,Vector{T}}(undef, V)
 end
+DiagonalTensorMap(::UndefInitializer, V::IndexSpace) = DiagonalTensorMap{Float64}(undef, V)
 
 function DiagonalTensorMap{T}(data::A, V::S) where {T,S<:IndexSpace,A<:DenseVector{T}}
     length(data) == reduceddim(V) ||
@@ -75,6 +78,14 @@ end
 TensorMap(d::DiagonalTensorMap) = copy!(similar(d), d)
 Base.convert(::Type{TensorMap}, d::DiagonalTensorMap) = TensorMap(d)
 
+function Base.convert(::Type{DiagonalTensorMap{T,S,A}},
+                      d::DiagonalTensorMap{T,S,A}) where {T,S,A}
+    return d
+end
+function Base.convert(D::Type{<:DiagonalTensorMap}, d::DiagonalTensorMap)
+    return DiagonalTensorMap(convert(storagetype(D), d.data), d.domain)
+end
+
 # Complex, real and imaginary parts
 #-----------------------------------
 for f in (:real, :imag, :complex)
@@ -87,22 +98,21 @@ end
 
 # Getting and setting the data at the block level
 #-------------------------------------------------
-blocksectors(d::DiagonalTensorMap) = blocksectors(d.domain)
-
 function block(d::DiagonalTensorMap, s::Sector)
     sectortype(d) == typeof(s) || throw(SectorMismatch())
     offset = 0
     dom = domain(d)[1]
-    for c in sectors(dom)
+    for c in blocksectors(d)
         if c < s
             offset += dim(dom, c)
         elseif c == s
             r = offset .+ (1:dim(dom, c))
             return Diagonal(view(d.data, r))
         else # s not in sectors(t)
-            return Diagonal(view(d.data, 1:0))
+            break
         end
     end
+    return Diagonal(view(d.data, 1:0))
 end
 
 # TODO: is relying on generic AbstractTensorMap blocks sufficient?
@@ -133,29 +143,33 @@ end
 # Index manipulations
 # -------------------
 function has_shared_permute(d::DiagonalTensorMap, (p₁, p₂)::Index2Tuple)
-    if p₁ === codomainind(d) && p₂ === domainind(d)
+    if p₁ === (1,) && p₂ === (2,)
         return true
-    elseif BraidingStyle(sectortype(d)) isa Bosonic
-        # TODO: is this always correct? transpose has no effect for bosonic sectors?
-        return p₁ === domainind(d) && p₂ === codomainind(d)
+    elseif p₁ === (2,) && p₂ === (1,) # transpose
+        return sectortype(d) === Trivial
     else
         return false
     end
 end
 
-function permute(d::DiagonalTensorMap, (p₁, p₂)::Index2Tuple{N₁,N₂};
-                 copy::Bool=false) where {N₁,N₂}
-    if !copy
-        if p₁ === codomainind(d) && p₂ === domainind(d)
-            return d
-        elseif has_shared_permute(d, (p₁, p₂)) # tranpose for bosonic sectors
-            return DiagonalTensorMap(d.data, dual(d.domain))
+function permute(d::DiagonalTensorMap, (p₁, p₂)::Index2Tuple{1,1};
+                 copy::Bool=false)
+    if p₁ === (1,) && p₂ === (2,)
+        return copy ? Base.copy(d) : d
+    elseif p₁ === (2,) && p₂ === (1,) # transpose
+        if has_shared_permute(d, (p₁, p₂)) # tranpose for bosonic sectors
+            return DiagonalTensorMap(copy ? Base.copy(d.data) : d.data, dual(d.domain))
         end
-    end
-    # general case
-    space′ = permute(space(d), (p₁, p₂))
-    @inbounds begin
-        return permute!(similar(d, space′), d, (p₁, p₂))
+        d′ = typeof(d)(undef, dual(d.domain))
+        for (c, b) in blocks(d)
+            f = only(fusiontrees(codomain(d), c))
+            ((f′, _), coeff) = only(permute(f, f, p₁, p₂))
+            c′ = f′.coupled
+            scale!(block(d′, c′), b, coeff)
+        end
+        return d′
+    else
+        throw(ArgumentError("invalid permutation $((p₁, p₂)) for tensor in space $(space(d))"))
     end
 end
 

src/tensors/indexmanipulations.jl

@@ -40,7 +40,6 @@ function permute(t::AbstractTensorMap, (p₁, p₂)::Index2Tuple{N₁,N₂};
     end
 end
 function permute(t::TensorMap, (p₁, p₂)::Index2Tuple{N₁,N₂}; copy::Bool=false) where {N₁,N₂}
-    @show (p₁, p₂)
     space′ = permute(space(t), (p₁, p₂))
     # share data if possible
     if !copy

src/tensors/tensor.jl

@@ -18,6 +18,9 @@ struct TensorMap{T,S<:IndexSpace,N₁,N₂,A<:DenseVector{T}} <: AbstractTensorM
                                                                            A<:DenseVector{T}}
         d = fusionblockstructure(space).totaldim
         data = A(undef, d)
+        if !isbitstype(T)
+            zerovector!(data)
+        end
         return TensorMap{T,S,N₁,N₂,A}(data, space)
     end
 

test/diagtensors.jl

@@ -0,0 +1,106 @@
+diagspacelist = ((ℂ^4)', ℂ[Z2Irrep](0 => 2, 1 => 3),
+                 ℂ[FermionNumber](0 => 2, 1 => 2, -1 => 1)',
+                 ℂ[SU2Irrep](0 => 2, 1 => 1), ℂ[FibonacciAnyon](:I => 2, :τ => 2))
+
+@testset "DiagonalTensor with domain $V" for V in diagspacelist
+    @timedtestset "Basic properties and algebra" begin
+        for T in (Float32, Float64, ComplexF32, ComplexF64, BigFloat)
+            t = @constinferred DiagonalTensorMap{T}(undef, V)
+            t = @constinferred DiagonalTensorMap(rand(T, reduceddim(V)), V)
+            @test @constinferred(hash(t)) == hash(deepcopy(t))
+            @test scalartype(t) == T
+            @test codomain(t) == ProductSpace(V)
+            @test domain(t) == ProductSpace(V)
+            @test space(t) == (V ← V)
+            @test space(t') == (V ← V)
+            @test dim(t) == dim(space(t))
+            @test isa(@constinferred(norm(t)), real(T))
+            @test norm(t)^2 ≈ dot(t, t)
+            α = rand(T)
+            @test norm(α * t) ≈ abs(α) * norm(t)
+            @test norm(t + t, 2) ≈ 2 * norm(t, 2)
+            @test norm(t + t, 1) ≈ 2 * norm(t, 1)
+            @test norm(t + t, Inf) ≈ 2 * norm(t, Inf)
+            p = 3 * rand(Float64)
+            @test norm(t + t, p) ≈ 2 * norm(t, p)
+            @test norm(t) ≈ norm(t')
+
+            @test t == @constinferred(TensorMap(t))
+            @test norm(t + TensorMap(t)) ≈ 2 * norm(t)
+
+            @test norm(zerovector!(t)) == 0
+            @test norm(one!(t)) ≈ sqrt(dim(V))
+            @test one!(t) == id(V)
+            @test norm(one!(t) - id(V)) == 0
+
+            t1 = DiagonalTensorMap(rand(T, reduceddim(V)), V)
+            t2 = DiagonalTensorMap(rand(T, reduceddim(V)), V)
+            t3 = DiagonalTensorMap(rand(T, reduceddim(V)), V)
+            α = rand(T)
+            β = rand(T)
+            @test @constinferred(dot(t1, t2)) ≈ conj(dot(t2, t1))
+            @test dot(t2, t1) ≈ conj(dot(t2', t1'))
+            @test dot(t3, α * t1 + β * t2) ≈ α * dot(t3, t1) + β * dot(t3, t2)
+        end
+    end
+    I = sectortype(V)
+    @timedtestset "Basic linear algebra: test via conversion" begin
+        for T in (Float32, ComplexF64)
+            t1 = DiagonalTensorMap(rand(T, reduceddim(V)), V)
+            t2 = DiagonalTensorMap(rand(T, reduceddim(V)), V)
+            @test norm(t1, 2) ≈ norm(convert(TensorMap, t1), 2)
+            @test dot(t2, t1) ≈ dot(convert(TensorMap, t2), convert(TensorMap, t1))
+            α = rand(T)
+            @test convert(TensorMap, α * t1) ≈ α * convert(TensorMap, t1)
+            @test convert(TensorMap, t1') ≈ convert(TensorMap, t1)'
+            @test convert(TensorMap, t1 + t2) ≈
+                  convert(TensorMap, t1) + convert(TensorMap, t2)
+        end
+    end
+    @timedtestset "Real and imaginary parts" begin
+        for T in (Float64, ComplexF64, ComplexF32)
+            t = DiagonalTensorMap(rand(T, reduceddim(V)), V)
+
+            tr = @constinferred real(t)
+            @test scalartype(tr) <: Real
+            @test real(convert(TensorMap, t)) == convert(TensorMap, tr)
+
+            ti = @constinferred imag(t)
+            @test scalartype(ti) <: Real
+            @test imag(convert(TensorMap, t)) == convert(TensorMap, ti)
+
+            tc = @inferred complex(t)
+            @test scalartype(tc) <: Complex
+            @test complex(convert(TensorMap, t)) == convert(TensorMap, tc)
+
+            tc2 = @inferred complex(tr, ti)
+            @test tc2 ≈ tc
+        end
+    end
+    @timedtestset "Tensor conversion" begin
+        t = @constinferred DiagonalTensorMap(undef, V)
+        rand!(t.data)
+        tc = complex(t)
+        @test convert(typeof(tc), t) == tc
+        @test typeof(convert(typeof(tc), t)) == typeof(tc)
+    end
+    I = sectortype(V)
+    if BraidingStyle(I) isa SymmetricBraiding
+        @timedtestset "Permutations" begin
+            t = DiagonalTensorMap(randn(ComplexF64, reduceddim(V)), V)
+            t1 = @constinferred permute(t, $(((2,), (1,))))
+            if BraidingStyle(sectortype(V)) isa Bosonic
+                @test t1 == transpose(t)
+            end
+            @test convert(TensorMap, t1) == permute(convert(TensorMap, t), (((2,), (1,))))
+            t2 = @constinferred permute(t, $(((1, 2), ())))
+            @test convert(TensorMap, t2) == permute(convert(TensorMap, t), (((1, 2), ())))
+            t3 = @constinferred permute(t, $(((2, 1), ())))
+            @test convert(TensorMap, t3) == permute(convert(TensorMap, t), (((2, 1), ())))
+            t4 = @constinferred permute(t, $(((), (1, 2))))
+            @test convert(TensorMap, t4) == permute(convert(TensorMap, t), (((), (1, 2))))
+            t5 = @constinferred permute(t, $(((), (2, 1))))
+            @test convert(TensorMap, t5) == permute(convert(TensorMap, t), (((), (2, 1))))
+        end
+    end
+end

test/runtests.jl

@@ -57,6 +57,7 @@ Ti = time()
 include("fusiontrees.jl")
 include("spaces.jl")
 include("tensors.jl")
+include("diagtensors.jl")
 include("planar.jl")
 include("ad.jl")
 include("bugfixes.jl")

test/spaces.jl

@@ -221,6 +221,7 @@ println("------------------------------------")
         slist = @constinferred sectors(V)
         @test @constinferred(TensorKit.hassector(V, first(slist)))
         @test @constinferred(dim(V)) == sum(dim(s) * dim(V, s) for s in slist)
+        @test @constinferred(reduceddim(V)) == sum(dim(V, s) for s in slist)
         @constinferred dim(V, first(slist))
         if hasfusiontensor(I)
             @test @constinferred(TensorKit.axes(V)) == Base.OneTo(dim(V))