Testsuite @planar (#76)

* Add Trivial anyon

* SpaceMismatches print more information

* Add `checkcontractible` for TensorMaps

* Add planar tests

* Add `tensorcost`

src/sectors/anyons.jl

@@ -1,3 +1,35 @@
+# PlanarTrivial
+"""
+    struct PlanarTrivial <: Sector
+    PlanarTrivial()
+
+Represents a trivial anyon sector, i.e. a trivial sector without braiding. This is mostly
+useful for testing.
+"""
+struct PlanarTrivial <: Sector end
+
+Base.IteratorSize(::Type{SectorValues{PlanarTrivial}}) = HasLength()
+Base.length(::SectorValues{PlanarTrivial}) = 1
+Base.iterate(::SectorValues{PlanarTrivial}, i = 0) = i == 0 ? (PlanarTrivial(), 1) : nothing
+function Base.getindex(::SectorValues{PlanarTrivial}, i::Int)
+    return i == 1 ? PlanarTrivial() :
+           throw(BoundsError(values(PlanarTrivial), i))
+end
+findindex(::SectorValues{PlanarTrivial}, c::PlanarTrivial) = 1
+Base.isless(::PlanarTrivial, ::PlanarTrivial) = false
+
+Base.one(::Type{PlanarTrivial}) = PlanarTrivial()
+Base.conj(::PlanarTrivial) = PlanarTrivial()
+
+FusionStyle(::Type{PlanarTrivial}) = UniqueFusion()
+BraidingStyle(::Type{PlanarTrivial}) = NoBraiding()
+Base.isreal(::Type{PlanarTrivial}) = true
+
+Nsymbol(::Vararg{PlanarTrivial,3}) = 1
+Fsymbol(::Vararg{PlanarTrivial,6}) = 1
+
+⊗(::PlanarTrivial, ::PlanarTrivial) = (PlanarTrivial(),)
+
 # FibonacciAnyons
 """
     struct FibonacciAnyon <: Sector

src/spaces/planarspace.jl

@@ -0,0 +1,11 @@
+struct PlanarNumbers end # this is probably not a field?
+const ℙ = PlanarNumbers()
+Base.show(io::IO, ::PlanarNumbers) = print(io, "ℙ")
+
+# convenience constructor
+Base.:^(::PlanarNumbers, d::Int) = Vect[PlanarTrivial](PlanarTrivial() => d)
+
+# convenience show
+function Base.show(io::IO, V::GradedSpace{PlanarTrivial})
+    return print(io, isdual(V) ? "(ℙ^$(dim(V)))'" : "ℙ^$(dim(V))")
+end

src/spaces/vectorspaces.jl

@@ -255,6 +255,7 @@ include("generalspace.jl")
 # space with internal structure corresponding to the irreducible representations of
 # a group, or more generally, the simple objects of a fusion category.
 include("gradedspace.jl")
+include("planarspace.jl")
 
 # Specific realizations of CompositeSpace types
 #-----------------------------------------------

src/tensors/linalg.jl

@@ -233,7 +233,7 @@ function LinearAlgebra.mul!(tC::AbstractTensorMap,
                             tA::AbstractTensorMap,
                             tB::AbstractTensorMap, α = true, β = false)
     if !(codomain(tC) == codomain(tA) && domain(tC) == domain(tB) && domain(tA) == codomain(tB))
-        throw(SpaceMismatch())
+        throw(SpaceMismatch("$(space(tC)) ≠ $(space(tA)) * $(space(tB))"))
     end
     for c in blocksectors(tC)
         if hasblock(tA, c) # then also tB should have such a block

src/tensors/tensoroperations.jl

@@ -136,6 +136,17 @@ function TO.tensorcontract_structure(pC::Index2Tuple{N₁,N₂},
     return dom → cod
 end
 
+function TO.checkcontractible(tA::AbstractTensorMap{S}, iA::Int, conjA::Symbol,
+                              tB::AbstractTensorMap{S}, iB::Int, conjB::Symbol, label) where {S}
+    sA = TO.tensorstructure(tA, iA, conjA)'
+    sB = TO.tensorstructure(tB, iB, conjB)
+    sA == sB ||
+        throw(SpaceMismatch("incompatible spaces for $label: $sA ≠ $sB"))
+    return nothing
+end
+
+TO.tensorcost(t::AbstractTensorMap, i::Int) = dim(space(t, i))
+
 #----------------
 # IMPLEMENTATONS
 #----------------

src/tensors/vectorinterface.jl

@@ -62,7 +62,7 @@ function VectorInterface.add(ty::AbstractTensorMap, tx::AbstractTensorMap, α::O
     return VectorInterface.add!(scale!(similar(ty, T), ty, β), tx, α)
 end
 function VectorInterface.add!(ty::AbstractTensorMap, tx::AbstractTensorMap, α::ONumber=_one, β::ONumber=_one)
-    space(ty) == space(tx) || throw(SpaceMismatch())
+    space(ty) == space(tx) || throw(SpaceMismatch("$(space(ty)) ≠ $(space(tx))"))
     for c in blocksectors(tx)
         VectorInterface.add!(block(ty, c), block(tx, c), α, β)
     end

test/planar.jl

@@ -0,0 +1,133 @@
+using TensorKit, TensorOperations, Test
+using TensorKit: planaradd!, planartrace!, planarcontract!
+using TensorKit: PlanarTrivial, ℙ
+
+"""
+    force_planar(obj)
+
+Replace an object with a planar equivalent -- i.e. one that disallows braiding.
+"""
+force_planar(V::ComplexSpace) = isdual(V) ? (ℙ^dim(V))' : ℙ^dim(V)
+function force_planar(V::GradedSpace)
+    return GradedSpace((c ⊠ PlanarTrivial() => dim(V, c) for c in sectors(V))..., isdual(V))
+end
+force_planar(V::ProductSpace) = mapreduce(force_planar, ⊗, V)
+function force_planar(tsrc::TensorMap{ComplexSpace})
+    tdst = TensorMap(undef, eltype(tsrc),
+                     force_planar(codomain(tsrc)) ← force_planar(domain(tsrc)))
+    copyto!(blocks(tdst)[PlanarTrivial()], blocks(tsrc)[Trivial()])
+    return tdst
+end
+function force_planar(tsrc::TensorMap{<:GradedSpace})
+    tdst = TensorMap(undef, eltype(tsrc),
+                     force_planar(codomain(tsrc)) ← force_planar(domain(tsrc)))
+    for (c, b) in blocks(tsrc)
+        copyto!(blocks(tdst)[c ⊠ PlanarTrivial()], b)
+    end
+    return tdst
+end
+
+@testset "planar methods" verbose = true begin
+    @testset "planaradd" begin
+        A = TensorMap(randn, ℂ^2 ⊗ ℂ^3 ← ℂ^6 ⊗ ℂ^5 ⊗ ℂ^4)
+        C = TensorMap(randn, (ℂ^5)' ⊗ (ℂ^6)' ← ℂ^4 ⊗ (ℂ^3)' ⊗ (ℂ^2)')
+        A′ = force_planar(A)
+        C′ = force_planar(C)
+        pC = ((4, 3), (5, 2, 1))
+
+        @test force_planar(tensoradd!(C, pC, A, :N, true, true)) ≈
+              planaradd!(C′, pC, A′, true, true)
+    end
+
+    @testset "planartrace" begin
+        A = TensorMap(randn, ℂ^2 ⊗ ℂ^3 ← ℂ^2 ⊗ ℂ^5 ⊗ ℂ^4)
+        C = TensorMap(randn, (ℂ^5)' ⊗ ℂ^3 ← ℂ^4)
+        A′ = force_planar(A)
+        C′ = force_planar(C)
+        pA = ((1,), (3,))
+        pC = ((4, 2), (5,))
+
+        @test force_planar(tensortrace!(C, pC, A, pA, :N, true, true)) ≈
+              planartrace!(C′, pC, A′, pA, true, true)
+    end
+
+    @testset "planarcontract" begin
+        A = TensorMap(randn, ℂ^2 ⊗ ℂ^3 ← ℂ^2 ⊗ ℂ^5 ⊗ ℂ^4)
+        B = TensorMap(randn, ℂ^2 ⊗ ℂ^4 ← ℂ^4 ⊗ ℂ^3)
+        C = TensorMap(randn, (ℂ^5)' ⊗ (ℂ^2)' ⊗ ℂ^2 ← (ℂ^2)' ⊗ ℂ^4)
+
+        A′ = force_planar(A)
+        B′ = force_planar(B)
+        C′ = force_planar(C)
+
+        pA = ((1, 3, 4), (5, 2))
+        pB = ((2, 4), (1, 3))
+        pC = ((3, 2, 1), (4, 5))
+
+        @test force_planar(tensorcontract!(C, pC, A, pA, :N, B, pB, :N, true, true)) ≈
+            planarcontract!(C′, pC, A′, pA, B′, pB, true, true)
+    end
+end
+
+@testset "@planar" verbose=true begin
+    T = ComplexF64
+    @testset "MPS networks" begin
+        P = ℂ^2
+        Vmps = ℂ^12
+        Vmpo = ℂ^4
+
+        # ∂AC
+        # -------
+        x = TensorMap(randn, T, Vmps ⊗ P ← Vmps)
+        O = TensorMap(randn, T, Vmpo ⊗ P ← P ⊗ Vmpo)
+        GL = TensorMap(randn, T, Vmps ⊗ Vmpo' ← Vmps)
+        GR = TensorMap(randn, T, Vmps ⊗ Vmpo ← Vmps)
+
+        x′ = force_planar(x)
+        O′ = force_planar(O)
+        GL′ = force_planar(GL)
+        GR′ = force_planar(GR)
+
+        @tensor y[-1 -2; -3] := GL[-1 2; 1] * x[1 3; 4] * O[2 -2; 3 5] * GR[4 5; -3]
+        @planar y′[-1 -2; -3] := GL′[-1 2; 1] * x′[1 3; 4] * O′[2 -2; 3 5] * GR′[4 5; -3]
+        @test force_planar(y) ≈ y′
+
+        # ∂AC2
+        # -------
+        x2 = TensorMap(randn, T, Vmps ⊗ P ← Vmps ⊗ P')
+        x2′ = force_planar(x2)
+        @tensor contractcheck=true y2[-1 -2; -3 -4] := GL[-1 7; 6] * x2[6 5; 1 3] * O[7 -2; 5 4] * O[4 -4; 3 2] * GR[1 2; -3]
+        @planar y2′[-1 -2; -3 -4] := GL′[-1 7; 6] * x2′[6 5; 1 3] * O′[7 -2; 5 4] * O′[4 -4; 3 2] * GR′[1 2; -3]
+        @test force_planar(y2) ≈ y2′
+    end
+
+    @testset "MERA networks" begin
+        Vmera = ℂ^2
+
+        u = TensorMap(randn, T, Vmera ⊗ Vmera ← Vmera ⊗ Vmera)
+        w = TensorMap(randn, T, Vmera ⊗ Vmera ← Vmera)
+        ρ = TensorMap(randn, T, Vmera ⊗ Vmera ⊗ Vmera ← Vmera ⊗ Vmera ⊗ Vmera)
+        h = TensorMap(randn, T, Vmera ⊗ Vmera ⊗ Vmera ← Vmera ⊗ Vmera ⊗ Vmera)
+
+        u′ = force_planar(u)
+        w′ = force_planar(w)
+        ρ′ = force_planar(ρ)
+        h′ = force_planar(h)
+
+        @tensor begin
+            C = (((((((h[9 3 4; 5 1 2] * u[1 2; 7 12]) * conj(u[3 4; 11 13])) *
+                     (u[8 5; 15 6] * w[6 7; 19])) *
+                    (conj(u[8 9; 17 10]) * conj(w[10 11; 22]))) *
+                   ((w[12 14; 20] * conj(w[13 14; 23])) * ρ[18 19 20; 21 22 23])) *
+                  w[16 15; 18]) * conj(w[16 17; 21]))
+        end
+        @planar begin
+            C′ = (((((((h′[9 3 4; 5 1 2] * u′[1 2; 7 12]) * conj(u′[3 4; 11 13])) *
+                     (u′[8 5; 15 6] * w′[6 7; 19])) *
+                    (conj(u′[8 9; 17 10]) * conj(w′[10 11; 22]))) *
+                   ((w′[12 14; 20] * conj(w′[13 14; 23])) * ρ′[18 19 20; 21 22 23])) *
+                  w′[16 15; 18]) * conj(w′[16 17; 21]))
+        end
+        @test C ≈ C′
+    end
+end

test/runtests.jl

@@ -61,6 +61,7 @@ include("sectors.jl")
 include("fusiontrees.jl")
 include("spaces.jl")
 include("tensors.jl")
+include("planar.jl")
 Tf = time()
 printstyled("Finished all tests in ",
             string(round((Tf - Ti) / 60; sigdigits=3)),