diff --git a/docs/src/lib/spaces.md b/docs/src/lib/spaces.md index 554e40c4..244a85f5 100644 --- a/docs/src/lib/spaces.md +++ b/docs/src/lib/spaces.md @@ -106,7 +106,6 @@ fuse ismonomorphic isepimorphic isisomorphic -insertunit ``` There are also specific methods for `HomSpace` instances, that are used in determining @@ -116,4 +115,6 @@ the resuling `HomSpace` after applying certain tensor operations. TensorKit.permute(::HomSpace{S}, ::Index2Tuple{N₁,N₂}) where {S,N₁,N₂} TensorKit.select(::HomSpace{S}, ::Index2Tuple{N₁,N₂}) where {S,N₁,N₂} TensorKit.compose(::HomSpace{S}, ::HomSpace{S}) where {S} +insertleftunit(::HomSpace, ::Int) +insertrightunit(::HomSpace, ::Int) ``` diff --git a/docs/src/lib/tensors.md b/docs/src/lib/tensors.md index 83ce1521..0ad8cf5a 100644 --- a/docs/src/lib/tensors.md +++ b/docs/src/lib/tensors.md @@ -175,6 +175,8 @@ braid(::AbstractTensorMap, ::Index2Tuple, ::IndexTuple; ::Bool) transpose(::AbstractTensorMap, ::Index2Tuple; ::Bool) repartition(::AbstractTensorMap, ::Int, ::Int; ::Bool) twist(::AbstractTensorMap, ::Int; ::Bool) +insertleftunit(::AbstractTensorMap, ::Int) +insertrightunit(::AbstractTensorMap, ::Int) ``` ```@docs @@ -224,4 +226,4 @@ and only accept the `TensorMap` object as well as the method-specific algorithm arguments. -TODO: document svd truncation types \ No newline at end of file +TODO: document svd truncation types diff --git a/src/TensorKit.jl b/src/TensorKit.jl index 8cdeaafb..3b103ea1 100644 --- a/src/TensorKit.jl +++ b/src/TensorKit.jl @@ -31,7 +31,8 @@ export TruncationScheme export SpaceMismatch, SectorMismatch, IndexError # error types # general vector space methods -export space, field, dual, dim, reduceddim, dims, fuse, flip, isdual, insertunit, oplus +export space, field, dual, dim, reduceddim, dims, fuse, flip, isdual, oplus, + insertleftunit, insertrightunit, removeunit # partial order for vector spaces export infimum, supremum, isisomorphic, ismonomorphic, isepimorphic diff --git a/src/auxiliary/auxiliary.jl b/src/auxiliary/auxiliary.jl index e4eb90fb..297fd3ab 100644 --- a/src/auxiliary/auxiliary.jl +++ b/src/auxiliary/auxiliary.jl @@ -41,6 +41,9 @@ function _kron(A, B) return C end +@noinline _boundserror(P, i) = throw(BoundsError(P, i)) +@noinline _nontrivialspaceerror(P, i) = throw(ArgumentError(lazy"Attempting to remove a non-trivial space $(P[i])")) + # Compat implementation: @static if VERSION < v"1.7" macro constprop(setting, ex) diff --git a/src/auxiliary/deprecate.jl b/src/auxiliary/deprecate.jl index b6a15899..fa7667b2 100644 --- a/src/auxiliary/deprecate.jl +++ b/src/auxiliary/deprecate.jl @@ -56,4 +56,6 @@ end Base.@deprecate EuclideanProduct() EuclideanInnerProduct() +Base.@deprecate insertunit(P::ProductSpace, args...; kwargs...) insertleftunit(args...; kwargs...) + #! format: on diff --git a/src/spaces/homspace.jl b/src/spaces/homspace.jl index c6569597..b67ba51b 100644 --- a/src/spaces/homspace.jl +++ b/src/spaces/homspace.jl @@ -174,6 +174,56 @@ function compose(W::HomSpace{S}, V::HomSpace{S}) where {S} return HomSpace(codomain(W), domain(V)) end +""" + insertleftunit(W::HomSpace, i::Int=numind(W) + 1; conj=false, dual=false) + +Insert a trivial vector space, isomorphic to the underlying field, at position `i`. +More specifically, adds a left monoidal unit or its dual. + +See also [`insertrightunit`](@ref), [`removeunit`](@ref). +""" +@constprop :aggressive function insertleftunit(W::HomSpace, i::Int=numind(W) + 1; + conj::Bool=false, dual::Bool=false) + if i ≤ numout(W) + return insertleftunit(codomain(W), i; conj, dual) ← domain(W) + else + return codomain(W) ← insertleftunit(domain(W), i - numout(W); conj, dual) + end +end + +""" + insertrightunit(W::HomSpace, i::Int=numind(W); conj=false, dual=false) + +Insert a trivial vector space, isomorphic to the underlying field, after position `i`. +More specifically, adds a right monoidal unit or its dual. + +See also [`insertleftunit`](@ref), [`removeunit`](@ref). +""" +@constprop :aggressive function insertrightunit(W::HomSpace, i::Int=numind(W); + conj::Bool=false, dual::Bool=false) + if i ≤ numout(W) + return insertrightunit(codomain(W), i; conj, dual) ← domain(W) + else + return codomain(W) ← insertrightunit(domain(W), i - numout(W); conj, dual) + end +end + +""" + removeunit(P::HomSpace, i::Int) + +This removes a trivial tensor product factor at position `1 ≤ i ≤ N`. +For this to work, that factor has to be isomorphic to the field of scalars. + +This operation undoes the work of [`insertleftunit`](@ref) or [`insertrightunit`](@ref). +""" +@constprop :aggressive function removeunit(P::HomSpace, i::Int) + if i ≤ numout(P) + return removeunit(codomain(P), i) ← domain(P) + else + return codomain(P) ← removeunit(domain(P), i - numout(P)) + end +end + # Block and fusion tree ranges: structure information for building tensors #-------------------------------------------------------------------------- struct FusionBlockStructure{I,N,F₁,F₂} diff --git a/src/spaces/productspace.jl b/src/spaces/productspace.jl index 11fceebf..b0736791 100644 --- a/src/spaces/productspace.jl +++ b/src/spaces/productspace.jl @@ -246,15 +246,15 @@ fuse(P::ProductSpace{S,0}) where {S<:ElementarySpace} = oneunit(S) fuse(P::ProductSpace{S}) where {S<:ElementarySpace} = fuse(P.spaces...) """ - insertunit(P::ProductSpace, i::Int = length(P)+1; dual = false, conj = false) + insertleftunit(P::ProductSpace, i::Int=length(P) + 1; conj=false, dual=false) -For `P::ProductSpace{S,N}`, this adds an extra tensor product factor at position -`1 <= i <= N+1` (last position by default) which is just the `S`-equivalent of the -underlying field of scalars, i.e. `oneunit(S)`. With the keyword arguments, one can choose -to insert the conjugated or dual space instead, which are all isomorphic to the field of -scalars. +Insert a trivial vector space, isomorphic to the underlying field, at position `i`. +More specifically, adds a left monoidal unit or its dual. + +See also [`insertrightunit`](@ref), [`removeunit`](@ref). """ -function insertunit(P::ProductSpace, i::Int=length(P) + 1; dual=false, conj=false) +function insertleftunit(P::ProductSpace, i::Int=length(P) + 1; + conj::Bool=false, dual::Bool=false) u = oneunit(spacetype(P)) if dual u = TensorKit.dual(u) @@ -265,6 +265,40 @@ function insertunit(P::ProductSpace, i::Int=length(P) + 1; dual=false, conj=fals return ProductSpace(TupleTools.insertafter(P.spaces, i - 1, (u,))) end +""" + insertrightunit(P::ProductSpace, i::Int=lenght(P); conj=false, dual=false) + +Insert a trivial vector space, isomorphic to the underlying field, after position `i`. +More specifically, adds a right monoidal unit or its dual. + +See also [`insertleftunit`](@ref), [`removeunit`](@ref). +""" +function insertrightunit(P::ProductSpace, i::Int=length(P); + conj::Bool=false, dual::Bool=false) + u = oneunit(spacetype(P)) + if dual + u = TensorKit.dual(u) + end + if conj + u = TensorKit.conj(u) + end + return ProductSpace(TupleTools.insertafter(P.spaces, i, (u,))) +end + +""" + removeunit(P::ProductSpace, i::Int) + +This removes a trivial tensor product factor at position `1 ≤ i ≤ N`. +For this to work, that factor has to be isomorphic to the field of scalars. + +This operation undoes the work of [`insertunit`](@ref). +""" +function removeunit(P::ProductSpace, i::Int) + 1 ≤ i ≤ length(P) || _boundserror(P, i) + isisomorphic(P[i], oneunit(P[i])) || _nontrivialspaceerror(P, i) + return ProductSpace{spacetype(P)}(TupleTools.deleteat(P.spaces, i)) +end + # Functionality for extracting and iterating over spaces #-------------------------------------------------------- Base.length(P::ProductSpace) = length(P.spaces) diff --git a/src/tensors/indexmanipulations.jl b/src/tensors/indexmanipulations.jl index a8bdb30a..55cbd64c 100644 --- a/src/tensors/indexmanipulations.jl +++ b/src/tensors/indexmanipulations.jl @@ -292,6 +292,84 @@ See [`twist!`](@ref) for storing the result in place. """ twist(t::AbstractTensorMap, i; inv::Bool=false) = twist!(copy(t), i; inv) +""" + insertleftunit(tsrc::AbstractTensorMap, i::Int=numind(t) + 1; + conj=false, dual=false, copy=false) -> tdst + +Insert a trivial vector space, isomorphic to the underlying field, at position `i`. +More specifically, adds a left monoidal unit or its dual. + +If `copy=false`, `tdst` might share data with `tsrc` whenever possible. Otherwise, a copy is always made. + +See also [`insertrightunit`](@ref) and [`removeunit`](@ref). +""" +@constprop :aggressive function insertleftunit(t::AbstractTensorMap, + i::Int=numind(t) + 1; copy::Bool=true, + conj::Bool=false, dual::Bool=false) + W = insertleftunit(space(t), i; conj, dual) + tdst = similar(t, W) + for (c, b) in blocks(t) + copy!(block(tdst, c), b) + end + return tdst +end +@constprop :aggressive function insertleftunit(t::TensorMap, i::Int=numind(t) + 1; + copy::Bool=false, + conj::Bool=false, dual::Bool=false) + W = insertleftunit(space(t), i; conj, dual) + return TensorMap{scalartype(t)}(copy ? Base.copy(t.data) : t.data, W) +end + +""" + insertrightunit(tsrc::AbstractTensorMap, i::Int=numind(t); + conj=false, dual=false, copy=false) -> tdst + +Insert a trivial vector space, isomorphic to the underlying field, after position `i`. +More specifically, adds a right monoidal unit or its dual. + +If `copy=false`, `tdst` might share data with `tsrc` whenever possible. Otherwise, a copy is always made. + +See also [`insertleftunit`](@ref) and [`removeunit`](@ref). +""" +@constprop :aggressive function insertrightunit(t::AbstractTensorMap, i::Int=numind(t); + copy::Bool=true, kwargs...) + W = insertrightunit(space(t), i; kwargs...) + tdst = similar(t, W) + for (c, b) in blocks(t) + copy!(block(tdst, c), b) + end + return tdst +end +@constprop :aggressive function insertrightunit(t::TensorMap, i::Int=numind(t); + copy::Bool=false, kwargs...) + W = insertrightunit(space(t), i; kwargs...) + return TensorMap{scalartype(t)}(copy ? Base.copy(t.data) : t.data, W) +end + +""" + removeunit(tsrc::AbstractTensorMap, i::Int; copy=false) -> tdst + +This removes a trivial tensor product factor at position `1 ≤ i ≤ N`. +For this to work, that factor has to be isomorphic to the field of scalars. + +If `copy=false`, `tdst` might share data with `tsrc` whenever possible. Otherwise, a copy is always made. + +This operation undoes the work of [`insertunit`](@ref). +""" +@constprop :aggressive function removeunit(t::TensorMap, i::Int; copy::Bool=false) + W = removeunit(space(t), i) + return TensorMap{scalartype(t)}(copy ? Base.copy(t.data) : t.data, W) +end +@constprop :aggressive function removeunit(t::AbstractTensorMap, i::Int; + copy::Bool=true) + W = removeunit(space(t), i) + tdst = similar(t, W) + for (c, b) in blocks(t) + copy!(block(tdst, c), b) + end + return tdst +end + # Fusing and splitting # TODO: add functionality for easy fusing and splitting of tensor indices diff --git a/test/spaces.jl b/test/spaces.jl index 57592ca3..e7a2bf55 100644 --- a/test/spaces.jl +++ b/test/spaces.jl @@ -278,7 +278,10 @@ println("------------------------------------") @test @constinferred(⊗(V1 ⊗ V2, V3 ⊗ V4)) == P @test @constinferred(⊗(V1, V2, V3 ⊗ V4)) == P @test @constinferred(⊗(V1, V2 ⊗ V3, V4)) == P - @test @constinferred(insertunit(P, 3)) == V1 * V2 * oneunit(V1) * V3 * V4 + @test V1 * V2 * oneunit(V1) * V3 * V4 == + @constinferred(insertleftunit(P, 3)) == + @constinferred(insertrightunit(P, 2)) + @test @constinferred(removeunit(V1 * V2 * oneunit(V1)' * V3 * V4, 3)) == P @test fuse(V1, V2', V3) ≅ V1 ⊗ V2' ⊗ V3 @test fuse(V1, V2', V3) ≾ V1 ⊗ V2' ⊗ V3 @test fuse(V1, V2', V3) ≿ V1 ⊗ V2' ⊗ V3 @@ -338,7 +341,10 @@ println("------------------------------------") @test @constinferred(*(V1, V2, V3)) == P @test @constinferred(⊗(V1, V2, V3)) == P @test @constinferred(adjoint(P)) == dual(P) == V3' ⊗ V2' ⊗ V1' - @test @constinferred(insertunit(P, 3; conj=true)) == V1 * V2 * oneunit(V1)' * V3 + @test V1 * V2 * oneunit(V1)' * V3 == + @constinferred(insertleftunit(P, 3; conj=true)) == + @constinferred(insertrightunit(P, 2; conj=true)) + @test P == @constinferred(removeunit(insertleftunit(P, 3), 3)) @test fuse(V1, V2', V3) ≅ V1 ⊗ V2' ⊗ V3 @test fuse(V1, V2', V3) ≾ V1 ⊗ V2' ⊗ V3 ≾ fuse(V1 ⊗ V2' ⊗ V3) @test fuse(V1, V2') ⊗ V3 ≾ V1 ⊗ V2' ⊗ V3 @@ -419,5 +425,21 @@ println("------------------------------------") @test W == @constinferred permute(W, ((1, 2), (3, 4, 5))) @test permute(W, ((2, 4, 5), (3, 1))) == (V2 ⊗ V4' ⊗ V5' ← V3 ⊗ V1') @test (V1 ⊗ V2 ← V1 ⊗ V2) == @constinferred TensorKit.compose(W, W') + @test (V1 ⊗ V2 ← V3 ⊗ V4 ⊗ V5 ⊗ oneunit(V5)) == + @constinferred(insertleftunit(W)) == + @constinferred(insertrightunit(W)) + @test @constinferred(removeunit(insertleftunit(W), $(numind(W) + 1))) == W + @test (V1 ⊗ V2 ← V3 ⊗ V4 ⊗ V5 ⊗ oneunit(V5)') == + @constinferred(insertleftunit(W; conj=true)) == + @constinferred(insertrightunit(W; conj=true)) + @test (oneunit(V1) ⊗ V1 ⊗ V2 ← V3 ⊗ V4 ⊗ V5) == + @constinferred(insertleftunit(W, 1)) == + @constinferred(insertrightunit(W, 0)) + @test (V1 ⊗ V2 ⊗ oneunit(V1) ← V3 ⊗ V4 ⊗ V5) == + @constinferred(insertrightunit(W, 2)) + @test (V1 ⊗ V2 ← oneunit(V1) ⊗ V3 ⊗ V4 ⊗ V5) == @constinferred(insertleftunit(W, 3)) + @test @constinferred(removeunit(insertleftunit(W, 3), 3)) == W + @test @constinferred(insertrightunit(one(V1) ← V1, 0)) == (oneunit(V1) ← V1) + @test_throws BoundsError insertleftunit(one(V1) ← V1, 0) end end diff --git a/test/tensors.jl b/test/tensors.jl index 4c3328e6..67565d34 100644 --- a/test/tensors.jl +++ b/test/tensors.jl @@ -173,6 +173,38 @@ for V in spacelist @test w * w' == (w * w')^2 end end + @timedtestset "Trivial spaces" begin + W = V1 ⊗ V2 ⊗ V3 ← V4 ⊗ V5 + for T in (Float32, ComplexF64) + t = @constinferred rand(T, W) + t2 = @constinferred insertleftunit(t) + @test t2 == @constinferred insertrightunit(t) + @test numind(t2) == numind(t) + 1 + @test space(t2) == insertleftunit(space(t)) + @test scalartype(t2) === T + @test t.data === t2.data + @test @constinferred(removeunit(t2, $(numind(t2)))) == t + t3 = @constinferred insertleftunit(t; copy=true) + @test t3 == @constinferred insertrightunit(t; copy=true) + @test t.data !== t3.data + for (c, b) in blocks(t) + @test b == block(t3, c) + end + @test @constinferred(removeunit(t3, $(numind(t3)))) == t + t4 = @constinferred insertrightunit(t, 3; dual=true) + @test numin(t4) == numin(t) && numout(t4) == numout(t) + 1 + for (c, b) in blocks(t) + @test b == block(t4, c) + end + @test @constinferred(removeunit(t4, 4)) == t + t5 = @constinferred insertleftunit(t, 4; dual=true) + @test numin(t5) == numin(t) + 1 && numout(t5) == numout(t) + for (c, b) in blocks(t) + @test b == block(t5, c) + end + @test @constinferred(removeunit(t5, 4)) == t + end + end if hasfusiontensor(I) @timedtestset "Basic linear algebra: test via conversion" begin W = V1 ⊗ V2 ⊗ V3 ← V4 ⊗ V5