Implement rotation for InfiniteWeightPEPS

src/algorithms/time_evolution/simpleupdate.jl

@@ -78,12 +78,6 @@ function su_iter(
             "iPEPS unit cell size for simple update should be no smaller than (2, 2)."
         ),
     )
-    # TODO: make algorithm independent on the choice of dual in the network
-    for (r, c) in Iterators.product(1:Nr, 1:Nc)
-        @assert [isdual(space(peps.vertices[r, c], ax)) for ax in 1:5] == [0, 1, 1, 0, 0]
-        @assert [isdual(space(peps.weights[1, r, c], ax)) for ax in 1:2] == [0, 1]
-        @assert [isdual(space(peps.weights[2, r, c], ax)) for ax in 1:2] == [0, 1]
-    end
     peps2 = deepcopy(peps)
     gate_mirrored = mirror_antidiag(gate)
     for direction in 1:2

src/networks/tensors.jl

@@ -86,6 +86,7 @@ function PEPSTensor(
     return f(T, ℂ^Pspace ← ℂ^Nspace ⊗ ℂ^Espace ⊗ (ℂ^Sspace)' ⊗ (ℂ^Wspace)')
 end
 
+mirror_antidiag(t::PEPSTensor) = permute(t, ((1,), (3, 2, 5, 4)))
 Base.rotl90(t::PEPSTensor) = permute(t, ((1,), (3, 4, 5, 2)))
 Base.rotr90(t::PEPSTensor) = permute(t, ((1,), (5, 2, 3, 4)))
 Base.rot180(t::PEPSTensor) = permute(t, ((1,), (4, 5, 2, 3)))

src/states/infiniteweightpeps.jl

@@ -10,7 +10,10 @@
 """
     struct SUWeight{E<:PEPSWeight}
 
-Schmidt bond weights used in simple/cluster update. Weight elements are always real.
+Schmidt bond weights used in simple/cluster update. 
+Weight elements are always real and non-negative.
+The domain and codomain of each weight matrix 
+must be an un-dualed `ElementarySapce`.
 
 ## Fields
 
@@ -22,11 +25,21 @@
 """
 struct SUWeight{E<:PEPSWeight}
     data::Array{E,3}
+    SUWeight{E}(data::Array{E,3}) where {E} = new{E}(data)
+end
 
-    function SUWeight(data::Array{E,3}) where {E<:PEPSWeight}
-        @assert eltype(data[1]) <: Real
-        return new{E}(data)
+function SUWeight(data::Array{E,3}) where {E<:PEPSWeight}
+    scalartype(data) <: Real || error("Weight elements must be real numbers.")
+    for wt in data
+        isa(wt, DiagonalTensorMap) ||
+            error("Each weight matrix should be a DiagonalTensorMap")
+        domain(wt, 1) == codomain(wt, 1) ||
+            error("Domain and codomain of each weight matrix must be the same.")
+        !isdual(codomain(wt, 1)) ||
+            error("Domain and codomain of each weight matrix cannot be a dual space.")
+        all(wt.data .>= 0) || error("Weight elements must be non-negative.")
     end
+    return SUWeight{E}(data)
 end
 
 function SUWeight(wts_mats::AbstractMatrix{E}...) where {E<:PEPSWeight}
@@ -50,6 +63,18 @@
 Base.getindex(W::SUWeight, args...) = Base.getindex(W.data, args...)
 Base.setindex!(W::SUWeight, args...) = (Base.setindex!(W.data, args...); W)
 Base.axes(W::SUWeight, args...) = axes(W.data, args...)
+Base.iterate(W::SUWeight, args...) = iterate(W.data, args...)
+
+## (Approximate) equality
+function Base.:(==)(wts1::SUWeight, wts2::SUWeight)
+    return wts1.data == wts2.data
+end
+function Base.isapprox(wts1::SUWeight, wts2::SUWeight; kwargs...)
+    for (wt1, wt2) in zip(wts1, wts2)
+        !isapprox(wt1, wt2; kwargs...) && return false
+    end
+    return true
+end
 
 function compare_weights(wts1::SUWeight, wts2::SUWeight)
     @assert size(wts1) == size(wts2)
@@ -72,6 +97,20 @@
 
 Represents an infinite projected entangled-pair state on a 2D square lattice
 consisting of vertex tensors and bond weights.
+The vertex tensor, x-weight and y-weight at row `i`, column `j` 
+are defined as (the numbers show the axis order)
+```
+        2
+        ↓
+        yᵢⱼ
+        ↓
+        1
+        2
+        ↓
+    5←-Tᵢⱼ←-3  1←-xᵢⱼ←-2
+        ↓ ↘
+        4   1
+```
 
 ## Fields
 
@@ -92,6 +131,7 @@
     ) where {T<:PEPSTensor,E<:PEPSWeight}
         @assert size(vertices) == size(weights)[2:end]
         Nr, Nc = size(vertices)
+        # check space matching between vertex tensors and weight matrices
         for (r, c) in Iterators.product(1:Nr, 1:Nc)
             space(weights[2, r, c], 1)' == space(vertices[r, c], 2) || throw(
                 SpaceMismatch("South space of bond weight y$((r, c)) does not match.")
@@ -136,6 +176,9 @@
 function InfiniteWeightPEPS(
     f, T, Pspaces::M, Nspaces::M, Espaces::M=Nspaces
 ) where {M<:AbstractMatrix{<:Union{Int,ElementarySpace}}}
+    @assert all(!isdual, Pspaces)
+    @assert all(!isdual, Nspaces)
+    @assert all(!isdual, Espaces)
     vertices = InfinitePEPS(f, T, Pspaces, Nspaces, Espaces).A
     Nr, Nc = size(vertices)
     weights = map(Iterators.product(1:2, 1:Nr, 1:Nc)) do (d, r, c)
@@ -296,17 +339,201 @@
     return InfiniteWeightPEPS(peps.A, SUWeight(weights))
 end
 
+## (Approximate) equality (gauge freedom is not allowed)
+function Base.:(==)(peps1::InfiniteWeightPEPS, peps2::InfiniteWeightPEPS)
+    return peps1.vertices == peps2.vertices && peps1.weights == peps2.weights
+end
+function Base.isapprox(peps1::InfiniteWeightPEPS, peps2::InfiniteWeightPEPS; kwargs...)
+    for (v1, v2) in zip(peps1.vertices, peps2.vertices)
+        !isapprox(v1, v2; kwargs...) && return false
+    end
+    !isapprox(peps1.weights, peps2.weights; kwargs...) && return false
+    return true
+end
+
+# Mirroring and rotation
+#= Example: 3 x 3 network
+
+- Original
+    ```
+            |         |         |
+            y₁₁       y₁₂       y₁₃
+            |         |         |
+    ..x₁₃...┼---x₁₁---┼---x₁₂---┼---x₁₃---
+            |         |         |           2
+            y₂₁       y₂₂       y₂₃         ↓
+            |         |         |           y
+    ..x₂₃...┼---x₂₁---┼---x₂₂---┼---x₂₃---  ↓
+            |         |         |           1
+            y₃₁       y₃₂       y₃₃
+            |         |         |           1 ←- x ←- 2
+    ..x₃₃...┼---x₃₁---┼---x₃₂---┼---x₃₃---
+            :         :         :
+            y₁₁       y₁₂       y₁₃
+            :         :         :
+    ```
+
+- After `mirror_antidiag`, x/y-weights are exchanged. 
+    ```
+            |         |         |
+            x₃₃       x₂₃       x₁₃
+            |         |         |
+    ..y₁₃...┼---y₃₃---┼---y₂₃---┼---y₁₃---
+            |         |         |           2
+            x₃₂       x₂₂       x₁₂         ↓
+            |         |         |           x
+    ..y₁₂...┼---y₃₂---┼---y₂₂---┼---y₁₂---  ↓
+            |         |         |           1
+            x₃₁       x₂₁       x₁₁
+            |         |         |           1 ←- y ←- 2
+    ..y₁₁...┼---y₃₁---┼---y₂₁---┼---y₁₁---
+            :         :         :
+            x₃₃       x₂₃       x₁₃
+            :         :         :
+    ```
+    No further operations are needed. 
+
+- After `rotl90`, x/y-weights are exchanged. 
+    ```
+            |         |         |
+            x₁₃       x₂₃       x₃₃
+            |         |         |
+    --y₁₃---┼---y₂₃---┼---y₃₃---┼...y₁₃...
+            |         |         |                     2
+            x₁₂       x₂₂       x₃₂                   ↓
+            |         |         |                     x
+    --y₁₂---┼---y₂₂---┼---y₃₂---┼...y₁₂...            ↓
+            |         |         |                     1
+            x₁₁       x₂₁       x₃₁
+            |         |         |           2 -→ y -→ 1
+    --y₁₁---┼---y₂₁---┼---y₃₁---┼...y₁₁...
+            :         :         :
+            x₁₃       x₂₃       x₃₃
+            :         :         :
+    ```
+    We need to further:
+    - Move 1st column of x-weights to the last column.
+    - Permute axes of x-weights.
+    - Flip x-arrows from → to ←. 
+
+- After `rotr90`, x/y-weights are exchanged. 
+    ```
+            :         :         :
+            x₃₃       x₂₃       x₁₃
+            :         :         :
+    ..y₁₁...┼---y₃₁---┼---y₂₁---┼---y₁₁---
+            |         |         |           1 ←- y ←- 2
+            x₃₁       x₂₁       x₁₁
+            |         |         |           1
+    ..y₁₂...┼---y₃₂---┼---y₂₂---┼---y₁₂---  ↑
+            |         |         |           x
+            x₃₂       x₂₂       x₁₂         ↑
+            |         |         |           2
+    ..y₁₃...┼---y₃₃---┼---y₂₃---┼---y₁₃---
+            |         |         |
+            x₃₃       x₂₃       x₁₃
+            |         |         |
+    ```
+    We need to further:
+    - Move last row of y-weights to the 1st row. 
+    - Permute axes of y-weights. 
+    - Flip y-arrows from ↑ to ↓. 
+
+After `rot180`, x/y-weights are not exchanged. 
+    ```
+            :         :         :
+            y₁₃       y₁₂       y₁₁
+            :         :         :
+    --x₃₃---┼---x₃₂---┼---x₃₁---┼...x₃₃...
+            |         |         |           2 -→ x -→ 1
+            y₃₃       y₃₂       y₃₁
+            |         |         |                     1
+    --x₂₃---┼---x₂₂---┼---x₂₁---┼...x₂₃...            ↑
+            |         |         |                     y
+            y₂₃       y₂₂       y₂₁                   ↑
+            |         |         |                     2
+    --x₁₃---┼---x₁₂---┼---x₁₁---┼...x₁₃...
+            |         |         |
+            y₁₃       y₁₂       y₁₁
+            |         |         |
+    ```
+    We need to further:
+    - Move 1st column of x-weights to the last column.
+    - Move last row of y-weights to the 1st row.
+    - Permute axes of all weights and twist their axis 1. 
+    - Flip x-arrows from → to ←, and y-arrows from ↑ to ↓. 
+=#
+
+"""
+Mirror an `SUWeight` by its anti-diagonal line.
+"""
+function mirror_antidiag(wts::SUWeight)
+    weights2_x = mirror_antidiag(wts[2, :, :])
+    weights2_y = mirror_antidiag(wts[1, :, :])
+    return SUWeight(weights2_x, weights2_y)
+end
+function Base.rotl90(wts::SUWeight)
+    wts_x = circshift(rotl90(wts[2, :, :]), (0, -1))
+    for (i, wt) in enumerate(wts_x)
+        wts_x[i] = DiagonalTensorMap(flip(permute(wt, ((2,), (1,))), (1, 2)))
+    end
+    wts_y = rotl90(wts[1, :, :])
+    return SUWeight(wts_x, wts_y)
+end
+function Base.rotr90(wts::SUWeight)
+    wts_x = rotr90(wts[2, :, :])
+    wts_y = circshift(rotr90(wts[1, :, :]), (1, 0))
+    for (i, wt) in enumerate(wts_y)
+        wts_y[i] = DiagonalTensorMap(flip(permute(wt, ((2,), (1,))), (1, 2)))
+    end
+    return SUWeight(wts_x, wts_y)
+end
+function Base.rot180(wts::SUWeight)
+    wts_x = circshift(rot180(wts[1, :, :]), (0, -1))
+    wts_y = circshift(rot180(wts[2, :, :]), (1, 0))
+    for (i, wt) in enumerate(wts_x)
+        wts_x[i] = DiagonalTensorMap(flip(permute(wt, ((2,), (1,))), (1, 2)))
+    end
+    for (i, wt) in enumerate(wts_y)
+        wts_y[i] = DiagonalTensorMap(flip(permute(wt, ((2,), (1,))), (1, 2)))
+    end
+    return SUWeight(wts_x, wts_y)
+end
+
 """
     mirror_antidiag(peps::InfiniteWeightPEPS)
 
-Mirror the unit cell of an iPEPS with weights by its anti-diagonal line.
+Mirror an `InfiniteWeightPEPS` by its anti-diagonal line.
 """
 function mirror_antidiag(peps::InfiniteWeightPEPS)
     vertices2 = mirror_antidiag(peps.vertices)
     for (i, t) in enumerate(vertices2)
-        vertices2[i] = permute(t, ((1,), (3, 2, 5, 4)))
+        vertices2[i] = mirror_antidiag(t)
+    end
+    weights2 = mirror_antidiag(peps.weights)
+    return InfiniteWeightPEPS(vertices2, weights2)
+end
+function Base.rotl90(peps::InfiniteWeightPEPS)
+    vertices2 = rotl90(peps.vertices)
+    for (i, t) in enumerate(vertices2)
+        vertices2[i] = flip(rotl90(t), (3, 5))
+    end
+    weights2 = rotl90(peps.weights)
+    return InfiniteWeightPEPS(vertices2, weights2)
+end
+function Base.rotr90(peps::InfiniteWeightPEPS)
+    vertices2 = rotr90(peps.vertices)
+    for (i, t) in enumerate(vertices2)
+        vertices2[i] = flip(rotr90(t), (2, 4))
+    end
+    weights2 = rotr90(peps.weights)
+    return InfiniteWeightPEPS(vertices2, weights2)
+end
+function Base.rot180(peps::InfiniteWeightPEPS)
+    vertices2 = rot180(peps.vertices)
+    for (i, t) in enumerate(vertices2)
+        vertices2[i] = flip(rot180(t), Tuple(2:5))
     end
-    weights2_x = mirror_antidiag(peps.weights[2, :, :])
-    weights2_y = mirror_antidiag(peps.weights[1, :, :])
-    return InfiniteWeightPEPS(vertices2, weights2_x, weights2_y)
+    weights2 = rot180(peps.weights)
+    return InfiniteWeightPEPS(vertices2, weights2)
 end

src/utility/util.jl

@@ -102,6 +102,19 @@
     return U * sqrt_S, sqrt_S * V
 end
 
+"""
+    flip_svd(u::AbstractTensorMap, s::DiagonalTensorMap, vh::AbstractTensorMap)
+
+Given `tsvd` result `u ← s ← vh`, flip the arrow between the three tensors 
+to `u2 → s2 → vh2` such that
+```
+    u * s * vh = (@tensor t2[-1; -2] := u2[-1; 1] * s2[1; 2] * vh2[2; -2])
+```
+"""
+function flip_svd(u::AbstractTensorMap, s::DiagonalTensorMap, vh::AbstractTensorMap)
+    return flip(u, 2), DiagonalTensorMap(flip(s, (1, 2))), flip(vh, 1)
+end
+
 """
     twistdual(t::AbstractTensorMap, i)
     twistdual!(t::AbstractTensorMap, i)

test/runtests.jl

@@ -61,6 +61,9 @@ end
         @time @safetestset "Norm-preserving retractions" begin
             include("utility/retractions.jl")
         end
+        @time @safetestset "Rotation of InfiniteWeightPEPS" begin
+            include("utility/iwpeps_rotation.jl")
+        end
     end
     if GROUP == "ALL" || GROUP == "EXAMPLES"
         @time @safetestset "Transverse Field Ising model" begin