whole lotta change

Author: Jutho

src/fusiontrees/braiding_manipulations.jl

@@ -64,7 +64,7 @@ function artin_braid(f::FusionTree{I, N}, i; inv::Bool = false) where {I, N}
     a = inner_extended[i - 1]
     c = inner_extended[i]
     e = inner_extended[i + 1]
-    c′ = first(a ⊗ d)
+    c′ = only(a ⊗ d)
     coeff = oftype(
         oneT,
         if inv
@@ -122,6 +122,8 @@ function artin_braid(src::FusionTreeBlock{I, N, 0}, i; inv::Bool = false) where
     BraidingStyle(I) isa NoBraiding &&
         throw(SectorMismatch(lazy"Cannot braid sectors $a and $b"))
 
+    T = typeof(oneT)
+    localbraidcache = Dict{NTuple{6, I}, FusionStyle(I) isa MultiplicityFreeFusion ? T : Array{T, 4}}()
     for (col, (f, f₂)) in enumerate(fusiontrees(src))
         inner = f.innerlines
         inner_extended = (uncoupled[1], inner..., coupled′)
@@ -156,10 +158,8 @@ function artin_braid(src::FusionTreeBlock{I, N, 0}, i; inv::Bool = false) where
         e = inner_extended[i + 1]
         if FusionStyle(I) isa MultiplicityFreeFusion
             for c′ in intersect(a ⊗ d, e ⊗ conj(b))
-                coeff = if inv
-                    conj(Rsymbol(d, c, e) * Fsymbol(d, a, b, e, c′, c)) * Rsymbol(d, a, c′)
-                else
-                    Rsymbol(c, d, e) * conj(Fsymbol(d, a, b, e, c′, c) * Rsymbol(a, d, c′))
+                coeff = let k = (a, b, c, d, e, c′)
+                    get!(() -> _artin_braid_local(k, inv), localbraidcache, k)
                 end
                 iszero(coeff) && continue
                 inner′ = TupleTools.setindex(inner, c′, i - 1)
@@ -172,15 +172,14 @@ function artin_braid(src::FusionTreeBlock{I, N, 0}, i; inv::Bool = false) where
                 Rmat1 = inv ? Rsymbol(d, c, e)' : Rsymbol(c, d, e)
                 Rmat2 = inv ? Rsymbol(d, a, c′)' : Rsymbol(a, d, c′)
                 Fmat = Fsymbol(d, a, b, e, c′, c)
+                coeff_tensor = let k = (a, b, c, d, e, c′)
+                    get!(() -> _artin_braid_local(k, inv), localbraidcache, k)
+                end
                 μ = vertices[i - 1]
                 ν = vertices[i]
-                for σ in 1:Nsymbol(a, d, c′)
-                    for λ in 1:Nsymbol(c′, b, e)
-                        coeff = zero(oneT)
-                        for ρ in 1:Nsymbol(d, c, e), κ in 1:Nsymbol(d, a, c′)
-                            coeff += Rmat1[ν, ρ] * conj(Fmat[κ, λ, μ, ρ]) *
-                                conj(Rmat2[σ, κ])
-                        end
+                for λ in 1:size(coeff_tensor, 2)
+                    for σ in 1:size(coeff_tensor, 1)
+                        coeff = coeff_tensor[σ, λ, μ, ν]
                         iszero(coeff) && continue
                         vertices′ = TupleTools.setindex(vertices, σ, i - 1)
                         vertices′ = TupleTools.setindex(vertices′, λ, i)
@@ -193,10 +192,26 @@ function artin_braid(src::FusionTreeBlock{I, N, 0}, i; inv::Bool = false) where
             end
         end
     end
-
     return dst => U
 end
 
+function _artin_braid_local((a, b, c, d, e, c′)::NTuple{6, I}, inv::Bool) where {I}
+    if FusionStyle(I) isa MultiplicityFreeFusion
+        coeff = if inv
+            conj(Rsymbol(d, c, e) * Fsymbol(d, a, b, e, c′, c)) * Rsymbol(d, a, c′)
+        else
+            Rsymbol(c, d, e) * conj(Fsymbol(d, a, b, e, c′, c) * Rsymbol(a, d, c′))
+        end
+        return coeff
+    else
+        Rmat1 = inv ? Rsymbol(d, c, e)' : Rsymbol(c, d, e)
+        Rmat2 = inv ? Rsymbol(d, a, c′)' : Rsymbol(a, d, c′)
+        Fmat = Fsymbol(d, a, b, e, c′, c)
+        @tensor coeff[σ, λ, μ, ν] := Rmat1[ν, ρ] * conj(Fmat[κ, λ, μ, ρ]) * conj(Rmat2[σ, κ])
+        return coeff
+    end
+end
+
 # braid fusion tree
 """
     braid(f::FusionTree{<:Sector, N}, p::NTuple{N, Int}, levels::NTuple{N, Int})
@@ -223,7 +238,7 @@ function braid(f::FusionTree{I, N}, (p, _)::Index2Tuple{N, 0}, (levels, _)::Inde
             for j in 1:(i - 1)
                 if p[j] > p[i]
                     a, b = f.uncoupled[p[j]], f.uncoupled[p[i]]
-                    coeff *= Rsymbol(a, b, first(a ⊗ b))
+                    coeff *= Rsymbol(a, b, only(a ⊗ b))
                 end
             end
         end
@@ -308,19 +323,22 @@ end
         f′, coeff2 = braid(f, p, levels)
         (f₁′, f₂′), coeff3 = repartition((f′, f0), N₁)
         return (f₁′, f₂′) => coeff1 * coeff2 * coeff3
-
     else
         src, (p1, p2), (l1, l2) = key
 
         p = linearizepermutation(p1, p2, numout(src), numin(src))
         levels = (l1..., reverse(l2)...)
 
-        dst, U = repartition(src, numind(src))
+        dst, U′ = repartition(src, numind(src))
+        T = sectorscalartype(I)
+        U = eltype(U′) == T ? U′ : T.(U′) # U′ has fusionscalartype(I) elements
+        Uold = similar(U)
 
         for s in permutation2swaps(p)
+            U, Uold = Uold, U
             inv = levels[s] > levels[s + 1]
             dst, U_tmp = artin_braid(dst, s; inv)
-            U = U_tmp * U
+            U = mul!(U, U_tmp, Uold)
             l = levels[s]
             levels = TupleTools.setindex(levels, levels[s + 1], s)
             levels = TupleTools.setindex(levels, l, s + 1)
@@ -329,8 +347,9 @@ end
         if N₂ == 0
             return dst => U
         else
+            U, Uold = Uold, U
             dst, U_tmp = repartition(dst, N₁)
-            U = U_tmp * U
+            U = mul!(U, U_tmp, Uold)
             return dst => U
         end
     end

src/spaces/productspace.jl

@@ -278,11 +278,7 @@ function insertleftunit(
         u = unitspace(spacetype(P))
     else
         N > 0 || throw(ArgumentError("cannot insert a sensible unit space in the empty product space"))
-        if i == N + 1
-            u = rightunitspace(P[N])
-        else
-            u = leftunitspace(P[i])
-        end
+        u = (i == N + 1) ? rightunitspace(P[N]) : leftunitspace(P[i])
     end
     if dual
         u = TensorKit.dual(u)
@@ -312,7 +308,7 @@ function insertrightunit(
         u = unitspace(spacetype(P))
     else
         N > 0 || throw(ArgumentError("cannot insert a sensible unit space in the empty product space"))
-        u = rightunitspace(P[i])
+        u = (i == 0 ) ? leftunitspace(P[1]) : rightunitspace(P[i])
     end
     if dual
         u = TensorKit.dual(u)

test/factorizations/eig.jl

@@ -17,8 +17,6 @@ for V in spacelist
     @timedtestset "Factorizations with symmetry: $Istr" verbose = true begin
         V1, V2, V3, V4, V5 = V
         W = V1 ⊗ V2
-        @assert !isempty(blocksectors(W))
-        @assert !isempty(intersect(blocksectors(V4), blocksectors(W)))
 
         @testset "Eigenvalue decomposition" begin
             for T in eltypes,

test/factorizations/ortho.jl

@@ -17,13 +17,11 @@ for V in spacelist
     @timedtestset "Factorizations with symmetry: $Istr" verbose = true begin
         V1, V2, V3, V4, V5 = V
         W = V1 ⊗ V2
-        @assert !isempty(blocksectors(W))
-        @assert !isempty(intersect(blocksectors(V4), blocksectors(W)))
 
         @testset "QR decomposition" begin
             for T in eltypes,
                     t in (
-                        rand(T, W, W), rand(T, W, W)', rand(T, W, V4), rand(T, V4, W)',
+                        rand(T, W, W), rand(T, W, W)', rand(T, (V1 ⊗ V2 ⊗ V3), (V4 ⊗ V5)'), rand(T, (V1 ⊗ V2)', (V3 ⊗ V4 ⊗ V5))',
                         DiagonalTensorMap(rand(T, reduceddim(V1)), V1),
                     )
 
@@ -76,7 +74,7 @@ for V in spacelist
         @testset "LQ decomposition" begin
             for T in eltypes,
                     t in (
-                        rand(T, W, W), rand(T, W, W)', rand(T, W, V4), rand(T, V4, W)',
+                        rand(T, W, W), rand(T, W, W)', rand(T, (V1 ⊗ V2), (V3 ⊗ V4 ⊗ V5)'), rand(T, (V1 ⊗ V2 ⊗ V3)', (V4 ⊗ V5))',
                         DiagonalTensorMap(rand(T, reduceddim(V1)), V1),
                     )
 

test/factorizations/projections.jl

@@ -17,46 +17,6 @@ for V in spacelist
     @timedtestset "Factorizations with symmetry: $Istr" verbose = true begin
         V1, V2, V3, V4, V5 = V
         W = V1 ⊗ V2
-        @assert !isempty(blocksectors(W))
-        @assert !isempty(intersect(blocksectors(V4), blocksectors(W)))
-
-        @testset "Condition number and rank" begin
-            for T in eltypes,
-                    t in (
-                        rand(T, W, W), rand(T, W, W)',
-                        rand(T, W, V4), rand(T, V4, W),
-                        rand(T, W, V4)', rand(T, V4, W)',
-                        DiagonalTensorMap(rand(T, reduceddim(V1)), V1),
-                    )
-
-                d1, d2 = dim(codomain(t)), dim(domain(t))
-                r = rank(t)
-                @test r == min(d1, d2)
-                @test typeof(r) == typeof(d1)
-                M = left_null(t)
-                @test @constinferred(rank(M)) + r ≈ d1
-                Mᴴ = right_null(t)
-                @test rank(Mᴴ) + r ≈ d2
-            end
-            for T in eltypes
-                u = unitary(T, V1 ⊗ V2, V1 ⊗ V2)
-                @test @constinferred(cond(u)) ≈ one(real(T))
-                @test @constinferred(rank(u)) == dim(V1 ⊗ V2)
-
-                t = rand(T, zerospace(V1), W)
-                @test rank(t) == 0
-                t2 = rand(T, zerospace(V1) * zerospace(V2), zerospace(V1) * zerospace(V2))
-                @test rank(t2) == 0
-                @test cond(t2) == 0.0
-            end
-            for T in eltypes, t in (rand(T, W, W), rand(T, W, W)')
-                project_hermitian!(t)
-                vals = @constinferred LinearAlgebra.eigvals(t)
-                λmax = maximum(s -> maximum(abs, s), values(vals))
-                λmin = minimum(s -> minimum(abs, s), values(vals))
-                @test cond(t) ≈ λmax / λmin
-            end
-        end
 
         @testset "Hermitian projections" begin
             for T in eltypes,
@@ -92,8 +52,7 @@ for V in spacelist
         @testset "Isometric projections" begin
             for T in eltypes,
                     t in (
-                        randn(T, W, W), randn(T, W, W)',
-                        randn(T, W, V4), randn(T, V4, W)',
+                        rand(T, W, W), rand(T, W, W)', rand(T, (V1 ⊗ V2 ⊗ V3), (V4 ⊗ V5)'), rand(T, (V1 ⊗ V2)', (V3 ⊗ V4 ⊗ V5))',
                     )
                 t2 = project_isometric(t)
                 @test isisometric(t2)

test/factorizations/svd.jl

@@ -17,13 +17,49 @@ for V in spacelist
     @timedtestset "Factorizations with symmetry: $Istr" verbose = true begin
         V1, V2, V3, V4, V5 = V
         W = V1 ⊗ V2
-        @assert !isempty(blocksectors(W))
-        @assert !isempty(intersect(blocksectors(V4), blocksectors(W)))
+
+        @testset "Condition number and rank" begin
+            for T in eltypes,
+                    t in (
+                        rand(T, W, W), rand(T, W, W)',
+                        rand(T, (V1 ⊗ V2 ⊗ V3), (V4 ⊗ V5)'), rand(T, (V1 ⊗ V2)', (V3 ⊗ V4 ⊗ V5))',
+                        rand(T, (V1 ⊗ V2), (V3 ⊗ V4 ⊗ V5)'), rand(T, (V1 ⊗ V2 ⊗ V3)', (V4 ⊗ V5))',
+                        DiagonalTensorMap(rand(T, reduceddim(V1)), V1),
+                    )
+
+                d1, d2 = dim(codomain(t)), dim(domain(t))
+                r = rank(t)
+                @test r == min(d1, d2)
+                @test typeof(r) == typeof(d1)
+                M = left_null(t)
+                @test @constinferred(rank(M)) + r ≈ d1
+                Mᴴ = right_null(t)
+                @test rank(Mᴴ) + r ≈ d2
+            end
+            for T in eltypes
+                u = unitary(T, V1 ⊗ V2, V1 ⊗ V2)
+                @test @constinferred(cond(u)) ≈ one(real(T))
+                @test @constinferred(rank(u)) == dim(V1 ⊗ V2)
+
+                t = rand(T, zerospace(V1), W)
+                @test rank(t) == 0
+                t2 = rand(T, zerospace(V1) * zerospace(V2), zerospace(V1) * zerospace(V2))
+                @test rank(t2) == 0
+                @test cond(t2) == 0.0
+            end
+            for T in eltypes, t in (rand(T, W, W), rand(T, W, W)')
+                project_hermitian!(t)
+                vals = @constinferred LinearAlgebra.eigvals(t)
+                λmax = maximum(s -> maximum(abs, s), values(vals))
+                λmin = minimum(s -> minimum(abs, s), values(vals))
+                @test cond(t) ≈ λmax / λmin
+            end
+        end
 
         @testset "Polar decomposition" begin
             for T in eltypes,
                     t in (
-                        rand(T, W, W), rand(T, W, W)', rand(T, W, V4), rand(T, V4, W)',
+                        rand(T, W, W), rand(T, (V1 ⊗ V2 ⊗ V3), (V4 ⊗ V5)'), rand(T, (V1 ⊗ V2)', (V3 ⊗ V4 ⊗ V5))',
                         DiagonalTensorMap(rand(T, reduceddim(V1)), V1),
                     )
 
@@ -39,7 +75,10 @@ for V in spacelist
             end
 
             for T in eltypes,
-                    t in (rand(T, W, W), rand(T, W, W)', rand(T, V4, W), rand(T, W, V4)')
+                    t in (
+                        rand(T, W, W), rand(T, W, W)', rand(T, (V1 ⊗ V2), (V3 ⊗ V4 ⊗ V5)'), rand(T, (V1 ⊗ V2 ⊗ V3)', (V4 ⊗ V5))',
+                        DiagonalTensorMap(rand(T, reduceddim(V1)), V1),
+                    )
 
                 @assert codomain(t) ≾ domain(t)
                 p, wᴴ = @constinferred right_polar(t)
@@ -57,8 +96,8 @@ for V in spacelist
             for T in eltypes,
                     t in (
                         rand(T, W, W), rand(T, W, W)',
-                        rand(T, W, V4), rand(T, V4, W),
-                        rand(T, W, V4)', rand(T, V4, W)',
+                        rand(T, (V1 ⊗ V2 ⊗ V3), (V4 ⊗ V5)'), rand(T, (V1 ⊗ V2)', (V3 ⊗ V4 ⊗ V5))',
+                        rand(T, (V1 ⊗ V2), (V3 ⊗ V4 ⊗ V5)'), rand(T, (V1 ⊗ V2 ⊗ V3)', (V4 ⊗ V5))',
                         DiagonalTensorMap(rand(T, reduceddim(V1)), V1),
                     )
 
@@ -121,10 +160,10 @@ for V in spacelist
         @testset "truncated SVD" begin
             for T in eltypes,
                     t in (
-                        randn(T, W, W), randn(T, W, W)',
-                        randn(T, W, V4), randn(T, V4, W),
-                        randn(T, W, V4)', randn(T, V4, W)',
-                        DiagonalTensorMap(randn(T, reduceddim(V1)), V1),
+                        rand(T, W, W), rand(T, W, W)',
+                        rand(T, (V1 ⊗ V2 ⊗ V3), (V4 ⊗ V5)'), rand(T, (V1 ⊗ V2)', (V3 ⊗ V4 ⊗ V5))',
+                        rand(T, (V1 ⊗ V2), (V3 ⊗ V4 ⊗ V5)'), rand(T, (V1 ⊗ V2 ⊗ V3)', (V4 ⊗ V5))',
+                        DiagonalTensorMap(rand(T, reduceddim(V1)), V1),
                     )
 
                 @constinferred normalize!(t)

test/mooncake/factorizations.jl

@@ -13,7 +13,6 @@ rng = Random.default_rng()
 spacelist = ad_spacelist(fast_tests)
 eltypes = (Float64, ComplexF64)
 
-
 @timedtestset "Mooncake - Factorizations: $(TensorKit.type_repr(sectortype(eltype(V)))) ($T)" for V in spacelist, T in eltypes
     atol = default_tol(T)
     rtol = default_tol(T)
@@ -31,7 +30,7 @@ eltypes = (Float64, ComplexF64)
         # TODO:
         # Mooncake.TestUtils.test_rule(rng, qr_null, A; atol, rtol, mode, is_primitive = false)
 
-        A = randn(T, V[1] ⊗ V[2] ← V[1])
+        A = randn(T, V[1] ⊗ V[2] ⊗ V[3] ← (V[4] ⊗ V[5])')
 
         Mooncake.TestUtils.test_rule(rng, qr_compact, A; atol, rtol, mode, is_primitive = false)
 
@@ -57,7 +56,7 @@ eltypes = (Float64, ComplexF64)
         # TODO:
         # Mooncake.TestUtils.test_rule(rng, lq_null, A; atol, rtol, mode, is_primitive = false)
 
-        A = randn(T, V[1] ⊗ V[2] ← V[1])
+        A = randn(T, V[1] ⊗ V[2] ← (V[3] ⊗ V[4] ⊗ V[5])')
 
         Mooncake.TestUtils.test_rule(rng, lq_compact, A; atol, rtol, mode, is_primitive = false)
 
@@ -86,7 +85,7 @@ eltypes = (Float64, ComplexF64)
     end
 
     @timedtestset "Singular value decomposition" begin
-        for t in (randn(T, V[1] ← V[1]), randn(T, V[1] ⊗ V[2] ← V[3] ⊗ V[4]))
+        for t in (randn(T, V[1] ← V[1]), randn(T, V[1] ⊗ V[2] ← (V[3] ⊗ V[4] ⊗ V[5])'))
             USVᴴ = svd_compact(t)
             ΔUSVᴴ = Mooncake.randn_tangent(rng, USVᴴ)
             remove_svdgauge_dependence!(ΔUSVᴴ[1], ΔUSVᴴ[3], USVᴴ...)