From c8fc7f3359b48f14a4734fa6a9c83a0d9356b0c7 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Thu, 7 May 2026 07:52:34 -0400 Subject: [PATCH 01/12] Add periodic getindex/setindex! utility --- src/utility/util.jl | 46 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 46 insertions(+) diff --git a/src/utility/util.jl b/src/utility/util.jl index 02183f33b..d6c11370d 100644 --- a/src/utility/util.jl +++ b/src/utility/util.jl @@ -2,6 +2,52 @@ _next(i, total) = mod1(i + 1, total) _prev(i, total) = mod1(i - 1, total) +@inline _periodic_inds(data, J::NTuple{N, Int}) where {N} = + ntuple(i -> mod1(J[i], size(data, i)), Val(N)) + +""" + periodic_getindex(A, data, I::Tuple) + +Periodic indexing of `A[I...]` backed by `data`. +For `I::Vararg{Int}`, this yields `data[mod1.(I, size(data))...]`. +Mixed index arguments (`Int` / `CartesianIndex`) are supported by flattening with `Base.to_indices`. + +!!! warning + Linear indexing schemes are not supported. +""" +@inline function periodic_getindex(A, data, I::Tuple) + J = to_indices(data, I) + return _periodic_getindex_dispatch(A, data, J) +end + +# Dispatch helper +_periodic_getindex_dispatch(A, data, J::Tuple{Vararg{Int}}) = @inbounds data[_periodic_inds(data, J)...] +_periodic_getindex_dispatch(A, data, J::Tuple) = throw(ArgumentError("Invalid periodic indexing type $(typeof(J))")) + +""" + periodic_setindex!(A, data, v, I::Tuple) + +Periodic indexing assignment `A[I...] = v` backed by `data`, returning `v` so that the +syntax `A[I...] = v` evaluates to `v` in line with Base Julia conventions. +For `I::Vararg{Int}`, this yields `data[mod1.(I, size(data))...] = v`. +Mixed index arguments (`Int` / `CartesianIndex`) are supported by flattening with `Base.to_indices`. + +!!! warning + Linear indexing schemes are not supported. +""" +@inline function periodic_setindex!(A, data, v, I::Tuple) + J = to_indices(data, I) + return _periodic_setindex_dispatch!(A, data, v, J) +end + +# Dispatch helper +function _periodic_setindex_dispatch!(A, data, v, J::Tuple{Vararg{Int}}) + @inbounds data[_periodic_inds(data, J)...] = v + return v +end +_periodic_setindex_dispatch!(_A, data, v, J::Tuple) = + throw(ArgumentError("Invalid periodic indexing type $(typeof(J))")) + # Get next and previous coordinate (direction, row, column), given a direction and going around the environment clockwise function _next_coordinate((dir, row, col), rowsize, colsize) if dir == 1 From a20a3bc7c7950882b1aac3d134f87e4f0b6491c2 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Thu, 7 May 2026 07:57:52 -0400 Subject: [PATCH 02/12] add periodic getters and setters to object types --- src/PEPSKit.jl | 1 + src/environments/bp_environments.jl | 5 +++- src/environments/ctmrg_environments.jl | 40 +++++++++++++++++++++++++ src/environments/suweight.jl | 6 ++-- src/networks/infinitesquarenetwork.jl | 13 ++++---- src/operators/infinitepepo.jl | 24 +++++++-------- src/states/infinitepartitionfunction.jl | 13 ++++---- src/states/infinitepeps.jl | 16 ++++------ 8 files changed, 75 insertions(+), 43 deletions(-) diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl index 56f1a6ce2..d669f02e1 100644 --- a/src/PEPSKit.jl +++ b/src/PEPSKit.jl @@ -133,6 +133,7 @@ using .Defaults: set_scheduler! export set_scheduler! export EighAdjoint, IterEigh, SVDAdjoint, IterSVD, QRAdjoint export CTMRGEnv, SequentialCTMRG, SimultaneousCTMRG +export corner, edge, setcorner!, setedge! export FixedSpaceTruncation, SiteDependentTruncation export HalfInfiniteProjector, FullInfiniteProjector export C4vCTMRG, C4vEighProjector, C4vQRProjector diff --git a/src/environments/bp_environments.jl b/src/environments/bp_environments.jl index ea7ae9301..30d6cdd70 100644 --- a/src/environments/bp_environments.jl +++ b/src/environments/bp_environments.jl @@ -145,7 +145,10 @@ end Base.eltype(::Type{BPEnv{T}}) where {T} = T Base.size(env::BPEnv, args...) = size(env.messages, args...) -Base.getindex(env::BPEnv, args...) = Base.getindex(env.messages, args...) +Base.@propagate_inbounds Base.getindex(env::BPEnv, I...) = + periodic_getindex(env, env.messages, I) +Base.@propagate_inbounds Base.setindex!(env::BPEnv, v, I...) = + periodic_setindex!(env, env.messages, v, I) Base.axes(env::BPEnv, args...) = Base.axes(env.messages, args...) Base.eachindex(env::BPEnv) = eachindex(IndexCartesian(), env.messages) VectorInterface.scalartype(::Type{BPEnv{T}}) where {T} = scalartype(T) diff --git a/src/environments/ctmrg_environments.jl b/src/environments/ctmrg_environments.jl index 8060595ec..17b4594d1 100644 --- a/src/environments/ctmrg_environments.jl +++ b/src/environments/ctmrg_environments.jl @@ -263,6 +263,46 @@ end Base.size(env::CTMRGEnv, args...) = size(env.corners, args...) Base.axes(x::CTMRGEnv, args...) = axes(x.corners, args...) + +""" + corner(env::CTMRGEnv, dir, i, j) + +Return the corner tensor of `env` at direction `dir` and unit-cell coordinates +`(i, j)`. The coordinates are wrapped periodically (`mod1`); the direction is +passed through unchanged. Mixed `Int`/`CartesianIndex` arguments are supported via +`Base.to_indices`, e.g. `corner(env, NORTHWEST, CartesianIndex(r, c))`. +""" +Base.@propagate_inbounds corner(env::CTMRGEnv, I...) = + periodic_getindex(env, env.corners, I) + +""" + edge(env::CTMRGEnv, dir, i, j) + +Return the edge tensor of `env` at direction `dir` and unit-cell coordinates +`(i, j)`. The coordinates are wrapped periodically (`mod1`); the direction is +passed through unchanged. Mixed `Int`/`CartesianIndex` arguments are supported via +`Base.to_indices`. +""" +Base.@propagate_inbounds edge(env::CTMRGEnv, I...) = + periodic_getindex(env, env.edges, I) + +""" + setcorner!(env::CTMRGEnv, value, dir, i, j) -> value + +Store `value` as the corner tensor of `env` at direction `dir` and unit-cell +coordinates `(i, j)`, with periodic wrapping on `(i, j)`. Returns `value`. +""" +Base.@propagate_inbounds setcorner!(env::CTMRGEnv, v, I...) = + periodic_setindex!(env, env.corners, v, I) + +""" + setedge!(env::CTMRGEnv, value, dir, i, j) -> value + +Store `value` as the edge tensor of `env` at direction `dir` and unit-cell +coordinates `(i, j)`, with periodic wrapping on `(i, j)`. Returns `value`. +""" +Base.@propagate_inbounds setedge!(env::CTMRGEnv, v, I...) = + periodic_setindex!(env, env.edges, v, I) function eachcoordinate(x::CTMRGEnv) return collect(Iterators.product(axes(x, 2), axes(x, 3))) end diff --git a/src/environments/suweight.jl b/src/environments/suweight.jl index e052d08d3..1ac771530 100644 --- a/src/environments/suweight.jl +++ b/src/environments/suweight.jl @@ -127,8 +127,10 @@ Base.eltype(W::SUWeight) = eltype(typeof(W)) Base.eltype(::Type{SUWeight{E}}) where {E} = E VI.scalartype(::Type{T}) where {T <: SUWeight} = scalartype(eltype(T)) -Base.getindex(W::SUWeight, args...) = Base.getindex(W.data, args...) -Base.setindex!(W::SUWeight, args...) = (Base.setindex!(W.data, args...); W) +Base.@propagate_inbounds Base.getindex(W::SUWeight, I...) = + periodic_getindex(W, W.data, I) +Base.@propagate_inbounds Base.setindex!(W::SUWeight, v, I...) = + periodic_setindex!(W, W.data, v, I) Base.axes(W::SUWeight, args...) = axes(W.data, args...) Base.iterate(W::SUWeight, args...) = iterate(W.data, args...) diff --git a/src/networks/infinitesquarenetwork.jl b/src/networks/infinitesquarenetwork.jl index 67adaf574..6f09c73ca 100644 --- a/src/networks/infinitesquarenetwork.jl +++ b/src/networks/infinitesquarenetwork.jl @@ -43,10 +43,10 @@ function Base.repeat(n::InfiniteSquareNetwork, counts...) end ## Indexing -Base.getindex(n::InfiniteSquareNetwork, args...) = Base.getindex(unitcell(n), args...) -function Base.setindex!(n::InfiniteSquareNetwork, args...) - return (Base.setindex!(unitcell(n), args...); n) -end +Base.@propagate_inbounds Base.getindex(n::InfiniteSquareNetwork, I...) = + periodic_getindex(n, unitcell(n), I) +Base.@propagate_inbounds Base.setindex!(n::InfiniteSquareNetwork, v, I...) = + periodic_setindex!(n, unitcell(n), v, I) Base.axes(n::InfiniteSquareNetwork, args...) = axes(unitcell(n), args...) eachcoordinate(n::InfiniteSquareNetwork) = collect(Iterators.product(axes(n)...)) function eachcoordinate(n::InfiniteSquareNetwork, dirs) @@ -56,10 +56,7 @@ end ## Spaces TensorKit.spacetype(::Type{T}) where {T <: InfiniteSquareNetwork} = spacetype(eltype(T)) -function virtualspace(n::InfiniteSquareNetwork, r::Int, c::Int, dir) - Nr, Nc = size(n) - return virtualspace(n[mod1(r, Nr), mod1(c, Nc)], dir) -end +virtualspace(n::InfiniteSquareNetwork, r::Int, c::Int, dir) = virtualspace(n[r, c], dir) ## Vector interface diff --git a/src/operators/infinitepepo.jl b/src/operators/infinitepepo.jl index 85270c636..a24af3901 100644 --- a/src/operators/infinitepepo.jl +++ b/src/operators/infinitepepo.jl @@ -138,8 +138,10 @@ function Base.similar(A::InfinitePEPO, T::Type{TorA} = scalartype(A)) where {Tor end Base.repeat(A::InfinitePEPO, counts...) = InfinitePEPO(repeat(unitcell(A), counts...)) -Base.getindex(A::InfinitePEPO, args...) = Base.getindex(unitcell(A), args...) -Base.setindex!(A::InfinitePEPO, args...) = (Base.setindex!(unitcell(A), args...); A) +Base.@propagate_inbounds Base.getindex(A::InfinitePEPO, I...) = + periodic_getindex(A, unitcell(A), I) +Base.@propagate_inbounds Base.setindex!(A::InfinitePEPO, v, I...) = + periodic_setindex!(A, unitcell(A), v, I) Base.axes(A::InfinitePEPO, args...) = axes(unitcell(A), args...) eachcoordinate(A::InfinitePEPO) = collect(Iterators.product(axes(A)...)) function eachcoordinate(A::InfinitePEPO, dirs) @@ -149,18 +151,12 @@ end ## Spaces TensorKit.spacetype(::Type{P}) where {P <: InfinitePEPO} = spacetype(eltype(P)) -function virtualspace(T::InfinitePEPO, r::Int, c::Int, h::Int, dir) - Nr, Nc, Nh = size(T) - return virtualspace(T[mod1(r, Nr), mod1(c, Nc), mod1(h, Nh)], dir) -end -function domain_physicalspace(T::InfinitePEPO, r::Int, c::Int) - Nr, Nc, = size(T) - return domain_physicalspace(T[mod1(r, Nr), mod1(c, Nc), 1]) -end -function codomain_physicalspace(T::InfinitePEPO, r::Int, c::Int) - Nr, Nc, = size(T) - return codomain_physicalspace(T[mod1(r, Nr), mod1(c, Nc), end]) -end +virtualspace(T::InfinitePEPO, r::Int, c::Int, h::Int, dir) = + virtualspace(T[r, c, h], dir) +domain_physicalspace(T::InfinitePEPO, r::Int, c::Int) = + domain_physicalspace(T[r, c, 1]) +codomain_physicalspace(T::InfinitePEPO, r::Int, c::Int) = + codomain_physicalspace(T[r, c, size(T, 3)]) physicalspace(T::InfinitePEPO) = [physicalspace(T, row, col) for row in axes(T, 1), col in axes(T, 2)] function physicalspace(T::InfinitePEPO, r::Int, c::Int) codomain_physicalspace(T, r, c) == domain_physicalspace(T, r, c) || throw( diff --git a/src/states/infinitepartitionfunction.jl b/src/states/infinitepartitionfunction.jl index aab1b5b7c..8adb2517f 100644 --- a/src/states/infinitepartitionfunction.jl +++ b/src/states/infinitepartitionfunction.jl @@ -126,10 +126,10 @@ function Base.repeat(A::InfinitePartitionFunction, counts...) return InfinitePartitionFunction(repeat(unitcell(A), counts...)) end -Base.getindex(A::InfinitePartitionFunction, args...) = Base.getindex(unitcell(A), args...) -function Base.setindex!(A::InfinitePartitionFunction, args...) - return (Base.setindex!(unitcell(A), args...); A) -end +Base.@propagate_inbounds Base.getindex(A::InfinitePartitionFunction, I...) = + periodic_getindex(A, unitcell(A), I) +Base.@propagate_inbounds Base.setindex!(A::InfinitePartitionFunction, v, I...) = + periodic_setindex!(A, unitcell(A), v, I) Base.axes(A::InfinitePartitionFunction, args...) = axes(unitcell(A), args...) eachcoordinate(A::InfinitePartitionFunction) = collect(Iterators.product(axes(A)...)) function eachcoordinate(A::InfinitePartitionFunction, dirs) @@ -139,10 +139,7 @@ end ## Spaces TensorKit.spacetype(::Type{T}) where {T <: InfinitePartitionFunction} = spacetype(eltype(T)) -function virtualspace(n::InfinitePartitionFunction, r::Int, c::Int, dir) - Nr, Nc = size(n) - return virtualspace(n[mod1(r, Nr), mod1(c, Nc)], dir) -end +virtualspace(n::InfinitePartitionFunction, r::Int, c::Int, dir) = virtualspace(n[r, c], dir) ## InfiniteSquareNetwork interface diff --git a/src/states/infinitepeps.jl b/src/states/infinitepeps.jl index 29b7f4242..f4c147549 100644 --- a/src/states/infinitepeps.jl +++ b/src/states/infinitepeps.jl @@ -131,8 +131,10 @@ function Base.similar(A::InfinitePEPS, T::Type{TorA} = scalartype(A)) where {Tor end Base.repeat(A::InfinitePEPS, counts...) = InfinitePEPS(repeat(unitcell(A), counts...)) -Base.getindex(A::InfinitePEPS, args...) = Base.getindex(unitcell(A), args...) -Base.setindex!(A::InfinitePEPS, args...) = (Base.setindex!(unitcell(A), args...); A) +Base.@propagate_inbounds Base.getindex(A::InfinitePEPS, I...) = + periodic_getindex(A, unitcell(A), I) +Base.@propagate_inbounds Base.setindex!(A::InfinitePEPS, v, I...) = + periodic_setindex!(A, unitcell(A), v, I) Base.axes(A::InfinitePEPS, args...) = axes(unitcell(A), args...) eachcoordinate(A::InfinitePEPS) = collect(Iterators.product(axes(A)...)) function eachcoordinate(A::InfinitePEPS, dirs) @@ -143,15 +145,9 @@ end TensorKit.spacetype(::Type{T}) where {T <: InfinitePEPS} = spacetype(eltype(T)) virtualspace(n::InfinitePEPS, dir) = virtualspace.(unitcell(n), dir) -function virtualspace(n::InfinitePEPS, r::Int, c::Int, dir) - Nr, Nc = size(n) - return virtualspace(n[mod1(r, Nr), mod1(c, Nc)], dir) -end +virtualspace(n::InfinitePEPS, r::Int, c::Int, dir) = virtualspace(n[r, c], dir) physicalspace(n::InfinitePEPS) = physicalspace.(unitcell(n)) -function physicalspace(n::InfinitePEPS, r::Int, c::Int) - Nr, Nc = size(n) - return physicalspace(n[mod1(r, Nr), mod1(c, Nc)]) -end +physicalspace(n::InfinitePEPS, r::Int, c::Int) = physicalspace(n[r, c]) ## InfiniteSquareNetwork interface From ed3e3d81386b6db1cedbae2f136c8e41ab796e48 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Thu, 7 May 2026 07:59:42 -0400 Subject: [PATCH 03/12] simplify codebase to avoid modular indexing --- src/algorithms/contractions/absorb_weight.jl | 10 +- .../contractions/bondenv/benv_ctm.jl | 22 ++-- .../contractions/bp_contractions.jl | 104 +++++++++--------- .../contractions/correlator/pepo_1layer.jl | 30 ++--- .../contractions/correlator/peps.jl | 30 ++--- src/algorithms/contractions/ctmrg/gaugefix.jl | 16 +-- .../contractions/ctmrg/renormalize_corner.jl | 8 +- .../contractions/ctmrg/renormalize_edge.jl | 10 +- src/algorithms/contractions/localoperator.jl | 100 ++++++++--------- src/algorithms/correlators.jl | 22 ++-- src/algorithms/ctmrg/sequential.jl | 6 +- src/algorithms/ctmrg/sparse_environments.jl | 28 +++-- src/algorithms/time_evolution/get_cluster.jl | 5 +- src/algorithms/time_evolution/simpleupdate.jl | 21 ++-- .../time_evolution/simpleupdate3site.jl | 5 +- src/algorithms/toolbox.jl | 68 ++++++------ 16 files changed, 232 insertions(+), 253 deletions(-) diff --git a/src/algorithms/contractions/absorb_weight.jl b/src/algorithms/contractions/absorb_weight.jl index 155034c3e..f6b764a5a 100644 --- a/src/algorithms/contractions/absorb_weight.jl +++ b/src/algorithms/contractions/absorb_weight.jl @@ -70,18 +70,16 @@ square root (or inverse square root if `inv = true`). function weight_to_absorb( weights::SUWeight, row::Int, col::Int, ax::Int; inv::Bool = false ) - _, Nr, Nc = size(weights) - r, c = mod1(row, Nr), mod1(col, Nc) pow = inv ? -1 / 2 : 1 / 2 wt = sdiag_pow( if ax == NORTH - weights[2, r, c] + weights[2, row, col] elseif ax == EAST - weights[1, r, c] + weights[1, row, col] elseif ax == SOUTH - weights[2, _next(r, Nr), c] + weights[2, row + 1, col] else # WEST - weights[1, r, _prev(c, Nc)] + weights[1, row, col - 1] end, pow, ) diff --git a/src/algorithms/contractions/bondenv/benv_ctm.jl b/src/algorithms/contractions/bondenv/benv_ctm.jl index 143dce59a..7e34db02e 100644 --- a/src/algorithms/contractions/bondenv/benv_ctm.jl +++ b/src/algorithms/contractions/bondenv/benv_ctm.jl @@ -24,20 +24,14 @@ Axis order: `[DX1 DY1; DX0 DY0]`, as in ``` """ function bondenv_ctm(row::Int, col::Int, X, Y, env::CTMRGEnv) - Nr, Nc = size(env.corners)[[2, 3]] - cm1 = _prev(col, Nc) - cp1 = _next(col, Nc) - cp2 = _next(cp1, Nc) - rm1 = _prev(row, Nr) - rp1 = _next(row, Nr) - c1 = env.corners[1, rm1, cm1] - c2 = env.corners[2, rm1, cp2] - c3 = env.corners[3, rp1, cp2] - c4 = env.corners[4, rp1, cm1] - e1X, e1Y = env.edges[1, rm1, col], env.edges[1, rm1, cp1] - e2 = env.edges[2, row, cp2] - e3X, e3Y = env.edges[3, rp1, col], env.edges[3, rp1, cp1] - e4 = env.edges[4, row, cm1] + c1 = corner(env, 1, row - 1, col - 1) + c2 = corner(env, 2, row - 1, col + 2) + c3 = corner(env, 3, row + 1, col + 2) + c4 = corner(env, 4, row + 1, col - 1) + e1X, e1Y = edge(env, 1, row - 1, col), edge(env, 1, row - 1, col + 1) + e2 = edge(env, 2, row, col + 2) + e3X, e3Y = edge(env, 3, row + 1, col), edge(env, 3, row + 1, col + 1) + e4 = edge(env, 4, row, col - 1) #= index labels C1--χ4--E1X---------χ6---------E1Y--χ8---C2 r-1 diff --git a/src/algorithms/contractions/bp_contractions.jl b/src/algorithms/contractions/bp_contractions.jl index d1de2d3ad..b751cbd75 100644 --- a/src/algorithms/contractions/bp_contractions.jl +++ b/src/algorithms/contractions/bp_contractions.jl @@ -105,14 +105,14 @@ function contract_local_operator1x1( env::BPEnv, ) row, col = Tuple(ind) - M_north = env.messages[NORTH, _prev(row, end), mod1(col, end)] - M_east = env.messages[EAST, mod1(row, end), _next(col, end)] - M_south = env.messages[SOUTH, _next(row, end), mod1(col, end)] - M_west = env.messages[WEST, mod1(row, end), _prev(col, end)] + M_north = env[NORTH, row - 1, col] + M_east = env[EAST, row, col + 1] + M_south = env[SOUTH, row + 1, col] + M_west = env[WEST, row, col - 1] return @autoopt @tensor begin - ket[mod1(row, end), mod1(col, end)][dt; DNt DEt DSt DWt] * - conj(bra[mod1(row, end), mod1(col, end)][db; DNb DEb DSb DWb]) * + ket[row, col][dt; DNt DEt DSt DWt] * + conj(bra[row, col][db; DNb DEb DSb DWb]) * O[db; dt] * M_north[DNt; DNb] * M_east[DEt; DEb] * @@ -125,14 +125,14 @@ function contract_local_norm1x1( ind::CartesianIndex{2}, ket::InfinitePEPS, bra::InfinitePEPS, env::BPEnv ) row, col = Tuple(ind) - M_north = env.messages[NORTH, _prev(row, end), mod1(col, end)] - M_east = env.messages[EAST, mod1(row, end), _next(col, end)] - M_south = env.messages[SOUTH, _next(row, end), mod1(col, end)] - M_west = env.messages[WEST, mod1(row, end), _prev(col, end)] + M_north = env[NORTH, row - 1, col] + M_east = env[EAST, row, col + 1] + M_south = env[SOUTH, row + 1, col] + M_west = env[WEST, row, col - 1] return @autoopt @tensor begin - ket[mod1(row, end), mod1(col, end)][d; DNt DEt DSt DWt] * - conj(bra[mod1(row, end), mod1(col, end)][d; DNb DEb DSb DWb]) * + ket[row, col][d; DNt DEt DSt DWt] * + conj(bra[row, col][d; DNb DEb DSb DWb]) * M_north[DNt; DNb] * M_east[DEt; DEb] * M_south[DSb; DSt] * @@ -148,17 +148,17 @@ function contract_local_operator2x1( env::BPEnv, ) row, col = Tuple(coord) - M_north = env.messages[NORTH, _prev(row, end), mod1(col, end)] - M_northeast = env.messages[EAST, mod1(row, end), _next(col, end)] - M_southeast = env.messages[EAST, _next(row, end), _next(col, end)] - M_south = env.messages[SOUTH, mod1(row + 2, end), mod1(col, end)] - M_southwest = env.messages[WEST, _next(row, end), _prev(col, end)] - M_northwest = env.messages[WEST, mod1(row, end), _prev(col, end)] + M_north = env[NORTH, row - 1, col] + M_northeast = env[EAST, row, col + 1] + M_southeast = env[EAST, row + 1, col + 1] + M_south = env[SOUTH, row + 2, col] + M_southwest = env[WEST, row + 1, col - 1] + M_northwest = env[WEST, row, col - 1] - return @autoopt @tensor ket[mod1(row, end), mod1(col, end)][dNt; DNt DNEt DMt DNWt] * - ket[_next(row, end), mod1(col, end)][dSt; DMt DSEt DSt DSWt] * - conj(bra[mod1(row, end), mod1(col, end)][dNb; DNb DNEb DMb DNWb]) * - conj(bra[_next(row, end), mod1(col, end)][dSb; DMb DSEb DSb DSWb]) * + return @autoopt @tensor ket[row, col][dNt; DNt DNEt DMt DNWt] * + ket[row + 1, col][dSt; DMt DSEt DSt DSWt] * + conj(bra[row, col][dNb; DNb DNEb DMb DNWb]) * + conj(bra[row + 1, col][dSb; DMb DSEb DSb DSWb]) * M_north[DNt; DNb] * M_northeast[DNEt; DNEb] * M_southeast[DSEt; DSEb] * @@ -176,16 +176,16 @@ function contract_local_operator1x2( env::BPEnv, ) row, col = Tuple(coord) - M_west = env.messages[WEST, mod1(row, end), _prev(col, end)] - M_northwest = env.messages[NORTH, _prev(row, end), mod1(col, end)] - M_northeast = env.messages[NORTH, _prev(row, end), _next(col, end)] - M_east = env.messages[EAST, mod1(row, end), mod1(col + 2, end)] - M_southeast = env.messages[SOUTH, _next(row, end), _next(col, end)] - M_southwest = env.messages[SOUTH, _next(row, end), mod1(col, end)] - A_west = ket[mod1(row, end), mod1(col, end)] - Ā_west = bra[mod1(row, end), mod1(col, end)] - A_east = ket[mod1(row, end), _next(col, end)] - Ā_east = bra[mod1(row, end), _next(col, end)] + M_west = env[WEST, row, col - 1] + M_northwest = env[NORTH, row - 1, col] + M_northeast = env[NORTH, row - 1, col + 1] + M_east = env[EAST, row, col + 2] + M_southeast = env[SOUTH, row + 1, col + 1] + M_southwest = env[SOUTH, row + 1, col] + A_west = ket[row, col] + Ā_west = bra[row, col] + A_east = ket[row, col + 1] + Ā_east = bra[row, col + 1] return @autoopt @tensor begin A_west[dWt; DNWt DMt DSWt DWt] * @@ -206,17 +206,17 @@ function contract_local_norm2x1( coord::CartesianIndex{2}, ket::InfinitePEPS, bra::InfinitePEPS, env::BPEnv ) row, col = Tuple(coord) - M_north = env.messages[NORTH, _prev(row, end), mod1(col, end)] - M_northeast = env.messages[EAST, mod1(row, end), _next(col, end)] - M_southeast = env.messages[EAST, _next(row, end), _next(col, end)] - M_south = env.messages[SOUTH, mod1(row + 2, end), mod1(col, end)] - M_southwest = env.messages[WEST, _next(row, end), _prev(col, end)] - M_northwest = env.messages[WEST, mod1(row, end), _prev(col, end)] + M_north = env[NORTH, row - 1, col] + M_northeast = env[EAST, row, col + 1] + M_southeast = env[EAST, row + 1, col + 1] + M_south = env[SOUTH, row + 2, col] + M_southwest = env[WEST, row + 1, col - 1] + M_northwest = env[WEST, row, col - 1] - return @autoopt @tensor ket[mod1(row, end), mod1(col, end)][dN; DNt DNEt DMt DNWt] * - ket[_next(row, end), mod1(col, end)][dS; DMt DSEt DSt DSWt] * - conj(bra[mod1(row, end), mod1(col, end)][dN; DNb DNEb DMb DNWb]) * - conj(bra[_next(row, end), mod1(col, end)][dS; DMb DSEb DSb DSWb]) * + return @autoopt @tensor ket[row, col][dN; DNt DNEt DMt DNWt] * + ket[row + 1, col][dS; DMt DSEt DSt DSWt] * + conj(bra[row, col][dN; DNb DNEb DMb DNWb]) * + conj(bra[row + 1, col][dS; DMb DSEb DSb DSWb]) * M_north[DNt; DNb] * M_northeast[DNEt; DNEb] * M_southeast[DSEt; DSEb] * @@ -230,17 +230,17 @@ function contract_local_norm1x2( ) row, col = Tuple(coord) - M_west = env.messages[WEST, mod1(row, end), _prev(col, end)] - M_northwest = env.messages[NORTH, _prev(row, end), mod1(col, end)] - M_northeast = env.messages[NORTH, _prev(row, end), _next(col, end)] - M_east = env.messages[EAST, mod1(row, end), mod1(col + 2, end)] - M_southeast = env.messages[SOUTH, _next(row, end), _next(col, end)] - M_southwest = env.messages[SOUTH, _next(row, end), mod1(col, end)] + M_west = env[WEST, row, col - 1] + M_northwest = env[NORTH, row - 1, col] + M_northeast = env[NORTH, row - 1, col + 1] + M_east = env[EAST, row, col + 2] + M_southeast = env[SOUTH, row + 1, col + 1] + M_southwest = env[SOUTH, row + 1, col] - A_west = ket[mod1(row, end), mod1(col, end)] - Ā_west = bra[mod1(row, end), mod1(col, end)] - A_east = ket[mod1(row, end), _next(col, end)] - Ā_east = bra[mod1(row, end), _next(col, end)] + A_west = ket[row, col] + Ā_west = bra[row, col] + A_east = ket[row, col + 1] + Ā_east = bra[row, col + 1] return @autoopt @tensor begin A_west[dW; DNWt DMt DSWt DWt] * diff --git a/src/algorithms/contractions/correlator/pepo_1layer.jl b/src/algorithms/contractions/correlator/pepo_1layer.jl index a8e7f0490..580f271d7 100644 --- a/src/algorithms/contractions/correlator/pepo_1layer.jl +++ b/src/algorithms/contractions/correlator/pepo_1layer.jl @@ -5,12 +5,12 @@ function start_correlator( (size(ρ, 3) == 1) || throw(ArgumentError("The input PEPO ρ must have only one layer.")) r, c = Tuple(i) - E_north = env.edges[NORTH, _prev(r, end), mod1(c, end)] - E_south = env.edges[SOUTH, _next(r, end), mod1(c, end)] - E_west = env.edges[WEST, mod1(r, end), _prev(c, end)] - C_northwest = env.corners[NORTHWEST, _prev(r, end), _prev(c, end)] - C_southwest = env.corners[SOUTHWEST, _next(r, end), _prev(c, end)] - t = twistdual(ρ[mod1(r, end), mod1(c, end)], 1:2) + E_north = edge(env, NORTH, r - 1, c) + E_south = edge(env, SOUTH, r + 1, c) + E_west = edge(env, WEST, r, c - 1) + C_northwest = corner(env, NORTHWEST, r - 1, c - 1) + C_southwest = corner(env, SOUTHWEST, r + 1, c - 1) + t = twistdual(ρ[r, c], 1:2) # TODO: part of these contractions is duplicated between the two output tensors, # so could be optimized further @autoopt @tensor Vn[χSE De; χNE] := @@ -32,12 +32,12 @@ function end_correlator_numerator( (size(ρ, 3) == 1) || throw(ArgumentError("The input PEPO ρ must have only one layer.")) r, c = Tuple(j) - E_north = env.edges[NORTH, _prev(r, end), mod1(c, end)] - E_east = env.edges[EAST, mod1(r, end), _next(c, end)] - E_south = env.edges[SOUTH, _next(r, end), mod1(c, end)] - C_northeast = env.corners[NORTHEAST, _prev(r, end), _next(c, end)] - C_southeast = env.corners[SOUTHEAST, _next(r, end), _next(c, end)] - t = twistdual(ρ[mod1(r, end), mod1(c, end)], 1:2) + E_north = edge(env, NORTH, r - 1, c) + E_east = edge(env, EAST, r, c + 1) + E_south = edge(env, SOUTH, r + 1, c) + C_northeast = corner(env, NORTHEAST, r - 1, c + 1) + C_southeast = corner(env, SOUTHEAST, r + 1, c + 1) + t = twistdual(ρ[r, c], 1:2) return @autoopt @tensor V[χSW DW dstring; χNW] * E_south[χSSE DS; χSW] * E_east[χNEE DE; χSEE] * E_north[χNW DN; χNNE] * C_northeast[χNNE; χNEE] * C_southeast[χSEE; χSSE] * @@ -48,9 +48,9 @@ function end_correlator_denominator( j::CartesianIndex{2}, V::CTMRGEdgeTensor{T, S, 2}, env::CTMRGEnv ) where {T, S} r, c = Tuple(j) - C_northeast = env.corners[NORTHEAST, _prev(r, end), _next(c, end)] - E_east = env.edges[EAST, mod1(r, end), _next(c, end)] - C_southeast = env.corners[SOUTHEAST, _next(r, end), _next(c, end)] + C_northeast = corner(env, NORTHEAST, r - 1, c + 1) + E_east = edge(env, EAST, r, c + 1) + C_southeast = corner(env, SOUTHEAST, r + 1, c + 1) return @autoopt @tensor V[χS DE; χN] * C_northeast[χN; χNE] * E_east[χNE DE; χSE] * C_southeast[χSE; χS] end diff --git a/src/algorithms/contractions/correlator/peps.jl b/src/algorithms/contractions/correlator/peps.jl index 697ad2a72..fafdd938a 100644 --- a/src/algorithms/contractions/correlator/peps.jl +++ b/src/algorithms/contractions/correlator/peps.jl @@ -6,12 +6,12 @@ function start_correlator( env::CTMRGEnv, ) r, c = Tuple(i) - E_north = env.edges[NORTH, _prev(r, end), mod1(c, end)] - E_south = env.edges[SOUTH, _next(r, end), mod1(c, end)] - E_west = env.edges[WEST, mod1(r, end), _prev(c, end)] - C_northwest = env.corners[NORTHWEST, _prev(r, end), _prev(c, end)] - C_southwest = env.corners[SOUTHWEST, _next(r, end), _prev(c, end)] - sandwich = (below[mod1(r, end), mod1(c, end)], above[mod1(r, end), mod1(c, end)]) + E_north = edge(env, NORTH, r - 1, c) + E_south = edge(env, SOUTH, r + 1, c) + E_west = edge(env, WEST, r, c - 1) + C_northwest = corner(env, NORTHWEST, r - 1, c - 1) + C_southwest = corner(env, SOUTHWEST, r + 1, c - 1) + sandwich = (below[r, c], above[r, c]) # TODO: part of these contractions is duplicated between the two output tensors, # so could be optimized further @@ -46,12 +46,12 @@ function end_correlator_numerator( env::CTMRGEnv, ) where {T, S} r, c = Tuple(j) - E_north = env.edges[NORTH, _prev(r, end), mod1(c, end)] - E_east = env.edges[EAST, mod1(r, end), _next(c, end)] - E_south = env.edges[SOUTH, _next(r, end), mod1(c, end)] - C_northeast = env.corners[NORTHEAST, _prev(r, end), _next(c, end)] - C_southeast = env.corners[SOUTHEAST, _next(r, end), _next(c, end)] - sandwich = (above[mod1(r, end), mod1(c, end)], below[mod1(r, end), mod1(c, end)]) + E_north = edge(env, NORTH, r - 1, c) + E_east = edge(env, EAST, r, c + 1) + E_south = edge(env, SOUTH, r + 1, c) + C_northeast = corner(env, NORTHEAST, r - 1, c + 1) + C_southeast = corner(env, SOUTHEAST, r + 1, c + 1) + sandwich = (above[r, c], below[r, c]) return @autoopt @tensor V[χSW DWt dstring DWb; χNW] * E_south[χSSE DSt DSb; χSW] * @@ -69,9 +69,9 @@ function end_correlator_denominator( env::CTMRGEnv ) where {T, S} r, c = Tuple(j) - C_northeast = env.corners[NORTHEAST, _prev(r, end), _next(c, end)] - E_east = env.edges[EAST, mod1(r, end), _next(c, end)] - C_southeast = env.corners[SOUTHEAST, _next(r, end), _next(c, end)] + C_northeast = corner(env, NORTHEAST, r - 1, c + 1) + E_east = edge(env, EAST, r, c + 1) + C_southeast = corner(env, SOUTHEAST, r + 1, c + 1) return @autoopt @tensor V[χS DEt DEb; χN] * C_northeast[χN; χNE] * diff --git a/src/algorithms/contractions/ctmrg/gaugefix.jl b/src/algorithms/contractions/ctmrg/gaugefix.jl index 554a8fbb8..b4d37acbb 100644 --- a/src/algorithms/contractions/ctmrg/gaugefix.jl +++ b/src/algorithms/contractions/ctmrg/gaugefix.jl @@ -30,7 +30,7 @@ Apply `fix_gauge_corner` to the northwest corner with appropriate row and column """ function fix_gauge_northwest_corner((row, col), env::CTMRGEnv, signs) return fix_gauge_corner( - env.corners[NORTHWEST, row, col], + corner(env, NORTHWEST, row, col), signs[WEST, row, col], signs[NORTH, row, _next(col, end)], ) @@ -43,7 +43,7 @@ Apply `fix_gauge_corner` to the northeast corner with appropriate row and column """ function fix_gauge_northeast_corner((row, col), env::CTMRGEnv, signs) return fix_gauge_corner( - env.corners[NORTHEAST, row, col], + corner(env, NORTHEAST, row, col), signs[NORTH, row, col], signs[EAST, _next(row, end), col], ) @@ -56,7 +56,7 @@ Apply `fix_gauge_corner` to the southeast corner with appropriate row and column """ function fix_gauge_southeast_corner((row, col), env::CTMRGEnv, signs) return fix_gauge_corner( - env.corners[SOUTHEAST, row, col], + corner(env, SOUTHEAST, row, col), signs[EAST, row, col], signs[SOUTH, row, _prev(col, end)], ) @@ -69,7 +69,7 @@ Apply `fix_gauge_corner` to the southwest corner with appropriate row and column """ function fix_gauge_southwest_corner((row, col), env::CTMRGEnv, signs) return fix_gauge_corner( - env.corners[SOUTHWEST, row, col], + corner(env, SOUTHWEST, row, col), signs[SOUTH, row, col], signs[WEST, _prev(row, end), col], ) @@ -113,7 +113,7 @@ Apply `fix_gauge_edge` to the north edge with appropriate row and column indices """ function fix_gauge_north_edge((row, col), env::CTMRGEnv, signs) return fix_gauge_edge( - env.edges[NORTH, row, col], signs[NORTH, row, col], signs[NORTH, row, _next(col, end)], + edge(env, NORTH, row, col), signs[NORTH, row, col], signs[NORTH, row, _next(col, end)], ) end @@ -124,7 +124,7 @@ Apply `fix_gauge_edge` to the east edge with appropriate row and column indices. """ function fix_gauge_east_edge((row, col), env::CTMRGEnv, signs) return fix_gauge_edge( - env.edges[EAST, row, col], signs[EAST, row, col], signs[EAST, _next(row, end), col] + edge(env, EAST, row, col), signs[EAST, row, col], signs[EAST, _next(row, end), col] ) end @@ -135,7 +135,7 @@ Apply `fix_gauge_edge` to the south edge with appropriate row and column indices """ function fix_gauge_south_edge((row, col), env::CTMRGEnv, signs) return fix_gauge_edge( - env.edges[SOUTH, row, col], signs[SOUTH, row, col], signs[SOUTH, row, _prev(col, end)], + edge(env, SOUTH, row, col), signs[SOUTH, row, col], signs[SOUTH, row, _prev(col, end)], ) end @@ -146,7 +146,7 @@ Apply `fix_gauge_edge` to the west edge with appropriate row and column indices. """ function fix_gauge_west_edge((row, col), env::CTMRGEnv, signs) return fix_gauge_edge( - env.edges[WEST, row, col], signs[WEST, row, col], signs[WEST, _prev(row, end), col] + edge(env, WEST, row, col), signs[WEST, row, col], signs[WEST, _prev(row, end), col] ) end diff --git a/src/algorithms/contractions/ctmrg/renormalize_corner.jl b/src/algorithms/contractions/ctmrg/renormalize_corner.jl index 5522b2849..51e2e49e0 100644 --- a/src/algorithms/contractions/ctmrg/renormalize_corner.jl +++ b/src/algorithms/contractions/ctmrg/renormalize_corner.jl @@ -407,8 +407,8 @@ Apply left projector to southwest corner and south edge. ``` """ function renormalize_southwest_corner((row, col), env::CTMRGEnv, projectors) - C_southwest = env.corners[SOUTHWEST, row, _prev(col, end)] - E_south = env.edges[SOUTH, row, col] + C_southwest = corner(env, SOUTHWEST, row, col - 1) + E_south = edge(env, SOUTH, row, col) P_left = projectors[1][row] return renormalize_southwest_corner(C_southwest, E_south, P_left) end @@ -449,8 +449,8 @@ Apply right projector to northwest corner and north edge. ``` """ function renormalize_northwest_corner((row, col), env::CTMRGEnv, projectors) - C_northwest = env.corners[NORTHWEST, row, _prev(col, end)] - E_north = env.edges[NORTH, row, col] + C_northwest = corner(env, NORTHWEST, row, col - 1) + E_north = edge(env, NORTH, row, col) P_right = projectors[2][_next(row, end)] return renormalize_northwest_corner(C_northwest, E_north, P_right) end diff --git a/src/algorithms/contractions/ctmrg/renormalize_edge.jl b/src/algorithms/contractions/ctmrg/renormalize_edge.jl index 3a4921787..a5f8e2a3f 100644 --- a/src/algorithms/contractions/ctmrg/renormalize_edge.jl +++ b/src/algorithms/contractions/ctmrg/renormalize_edge.jl @@ -19,7 +19,7 @@ function renormalize_north_edge( (row, col), env::CTMRGEnv, P_left, P_right, network::InfiniteSquareNetwork ) return renormalize_north_edge( - env.edges[NORTH, _prev(row, end), col], + edge(env, NORTH, row - 1, col), P_left[NORTH, row, col], P_right[NORTH, row, _prev(col, end)], network[row, col], # so here it's fine @@ -90,7 +90,7 @@ function renormalize_east_edge( (row, col), env::CTMRGEnv, P_left, P_right, network::InfiniteSquareNetwork ) return renormalize_east_edge( - env.edges[EAST, row, _next(col, end)], + edge(env, EAST, row, col + 1), P_left[EAST, row, col], P_right[EAST, _prev(row, end), col], network[row, col], @@ -155,7 +155,7 @@ function renormalize_south_edge( (row, col), env::CTMRGEnv, P_left, P_right, network::InfiniteSquareNetwork ) return renormalize_south_edge( - env.edges[SOUTH, _next(row, end), col], + edge(env, SOUTH, row + 1, col), P_left[SOUTH, row, col], P_right[SOUTH, row, _next(col, end)], network[row, col], @@ -225,7 +225,7 @@ function renormalize_west_edge( # For simultaneous CTMRG scheme (row, col), env::CTMRGEnv, P_left, P_right, network::InfiniteSquareNetwork, ) return renormalize_west_edge( - env.edges[WEST, row, _prev(col, end)], + edge(env, WEST, row, col - 1), P_left[WEST, row, col], P_right[WEST, _next(row, end), col], network[row, col], @@ -235,7 +235,7 @@ function renormalize_west_edge( # For sequential CTMRG scheme (row, col), env::CTMRGEnv, projectors, network::InfiniteSquareNetwork, ) return renormalize_west_edge( - env.edges[WEST, row, _prev(col, end)], + edge(env, WEST, row, col - 1), projectors[1][row], projectors[2][_next(row, end)], network[row, col], diff --git a/src/algorithms/contractions/localoperator.jl b/src/algorithms/contractions/localoperator.jl index ff734f00a..56364cdde 100644 --- a/src/algorithms/contractions/localoperator.jl +++ b/src/algorithms/contractions/localoperator.jl @@ -43,16 +43,16 @@ function _contract_corner_expr(rowrange, colrange) cmin, cmax = extrema(colrange) gridsize = (rmax - rmin + 1, cmax - cmin + 1) - C_NW = :(env.corners[NORTHWEST, mod1($(rmin - 1), end), mod1($(cmin - 1), end)]) + C_NW = :(corner(env, NORTHWEST, $(rmin - 1), $(cmin - 1))) corner_NW = tensorexpr(C_NW, envlabel(WEST, 0), envlabel(NORTH, 0)) - C_NE = :(env.corners[NORTHEAST, mod1($(rmin - 1), end), mod1($(cmax + 1), end)]) + C_NE = :(corner(env, NORTHEAST, $(rmin - 1), $(cmax + 1))) corner_NE = tensorexpr(C_NE, envlabel(NORTH, gridsize[2]), envlabel(EAST, 0)) - C_SE = :(env.corners[SOUTHEAST, mod1($(rmax + 1), end), mod1($(cmax + 1), end)]) + C_SE = :(corner(env, SOUTHEAST, $(rmax + 1), $(cmax + 1))) corner_SE = tensorexpr(C_SE, envlabel(EAST, gridsize[1]), envlabel(SOUTH, gridsize[2])) - C_SW = :(env.corners[SOUTHWEST, mod1($(rmax + 1), end), mod1($(cmin - 1), end)]) + C_SW = :(corner(env, SOUTHWEST, $(rmax + 1), $(cmin - 1))) corner_SW = tensorexpr(C_SW, envlabel(SOUTH, 0), envlabel(WEST, gridsize[1])) return corner_NW, corner_NE, corner_SE, corner_SW @@ -64,7 +64,7 @@ function _contract_edge_expr(rowrange, colrange, height) gridsize = (rmax - rmin + 1, cmax - cmin + 1) edges_N = map(1:gridsize[2]) do i - E_N = :(env.edges[NORTH, mod1($(rmin - 1), end), mod1($(cmin + i - 1), end)]) + E_N = :(edge(env, NORTH, $(rmin - 1), $(cmin + i - 1))) return tensorexpr( E_N, (envlabel(NORTH, i - 1), virtuallabel.(NORTH, ntuple(identity, height), i)...), @@ -73,7 +73,7 @@ function _contract_edge_expr(rowrange, colrange, height) end edges_E = map(1:gridsize[1]) do i - E_E = :(env.edges[EAST, mod1($(rmin + i - 1), end), mod1($(cmax + 1), end)]) + E_E = :(edge(env, EAST, $(rmin + i - 1), $(cmax + 1))) return tensorexpr( E_E, (envlabel(EAST, i - 1), virtuallabel.(EAST, ntuple(identity, height), i)...), @@ -82,7 +82,7 @@ function _contract_edge_expr(rowrange, colrange, height) end edges_S = map(1:gridsize[2]) do i - E_S = :(env.edges[SOUTH, mod1($(rmax + 1), end), mod1($(cmin + i - 1), end)]) + E_S = :(edge(env, SOUTH, $(rmax + 1), $(cmin + i - 1))) return tensorexpr( E_S, (envlabel(SOUTH, i), virtuallabel.(SOUTH, ntuple(identity, height), i)...), @@ -91,7 +91,7 @@ function _contract_edge_expr(rowrange, colrange, height) end edges_W = map(1:gridsize[1]) do i - E_W = :(env.edges[WEST, mod1($(rmin + i - 1), end), mod1($(cmin - 1), end)]) + E_W = :(edge(env, WEST, $(rmin + i - 1), $(cmin - 1))) return tensorexpr( E_W, (envlabel(WEST, i), virtuallabel.(WEST, ntuple(identity, height), i)...), @@ -120,7 +120,7 @@ function _contract_state_expr(rowrange, colrange, height, cartesian_inds = nothi physicallabel(:O, side, inds_id) end return tensorexpr( - :(state[$(side)][mod1($(rmin + i - 1), end), mod1($(cmin + j - 1), end)]), + :(state[$(side)][$(rmin + i - 1), $(cmin + j - 1)]), (physical_label,), ( if i == 1 @@ -172,7 +172,7 @@ function _contract_pepo_state_expr(rowrange, colrange, height, cartesian_inds = else physicallabel(:O, 1, inds_id) end - tensor_name = :(twistdual(state[1][mod1($(rmin + i - 1), end), mod1($(cmin + j - 1), end)], 2)) + tensor_name = :(twistdual(state[1][$(rmin + i - 1), $(cmin + j - 1)], 2)) else physical_label_in = if isnothing(inds_id) physicallabel(:in, i, j) @@ -185,9 +185,9 @@ function _contract_pepo_state_expr(rowrange, colrange, height, cartesian_inds = physicallabel(:O, side, inds_id) end tensor_name = if side == 2 - :(state[2][mod1($(rmin + i - 1), end), mod1($(cmin + j - 1), end)]) + :(state[2][$(rmin + i - 1), $(cmin + j - 1)]) else - :(twistdual(state[1][mod1($(rmin + i - 1), end), mod1($(cmin + j - 1), end)], (1, 2))) + :(twistdual(state[1][$(rmin + i - 1), $(cmin + j - 1)], (1, 2))) end end @@ -396,21 +396,21 @@ function reduced_densitymatrix1x1( row, col = Tuple(inds) # Unpack variables and absorb corners - A = ket[mod1(row, end), mod1(col, end)] - Ā = bra[mod1(row, end), mod1(col, end)] + A = ket[row, col] + Ā = bra[row, col] E_north = - env.edges[NORTH, mod1(row - 1, end), mod1(col, end)] * - twistdual(env.corners[NORTHEAST, mod1(row - 1, end), mod1(col + 1, end)], 1) + edge(env, NORTH, row - 1, col) * + twistdual(corner(env, NORTHEAST, row - 1, col + 1), 1) E_east = - env.edges[EAST, mod1(row, end), mod1(col + 1, end)] * - twistdual(env.corners[SOUTHEAST, mod1(row + 1, end), mod1(col + 1, end)], 1) + edge(env, EAST, row, col + 1) * + twistdual(corner(env, SOUTHEAST, row + 1, col + 1), 1) E_south = - env.edges[SOUTH, mod1(row + 1, end), mod1(col, end)] * - twistdual(env.corners[SOUTHWEST, mod1(row + 1, end), mod1(col - 1, end)], 1) + edge(env, SOUTH, row + 1, col) * + twistdual(corner(env, SOUTHWEST, row + 1, col - 1), 1) E_west = - env.edges[WEST, mod1(row, end), mod1(col - 1, end)] * - twistdual(env.corners[NORTHWEST, mod1(row - 1, end), mod1(col - 1, end)], 1) + edge(env, WEST, row, col - 1) * + twistdual(corner(env, NORTHWEST, row - 1, col - 1), 1) @tensor EE_SW[χSE χNW DSb DWb; DSt DWt] := E_south[χSE DSt DSb; χSW] * E_west[χSW DWt DWb; χNW] @@ -450,25 +450,25 @@ function reduced_densitymatrix2x1( row, col = Tuple(ind) # Unpack variables and absorb corners - A_north = ket[mod1(row, end), mod1(col, end)] - Ā_north = bra[mod1(row, end), mod1(col, end)] - A_south = ket[mod1(row + 1, end), mod1(col, end)] - Ā_south = bra[mod1(row + 1, end), mod1(col, end)] + A_north = ket[row, col] + Ā_north = bra[row, col] + A_south = ket[row + 1, col] + Ā_south = bra[row + 1, col] E_north = - env.edges[NORTH, mod1(row - 1, end), mod1(col, end)] * - twistdual(env.corners[NORTHEAST, mod1(row - 1, end), mod1(col + 1, end)], 1) - E_northeast = env.edges[EAST, mod1(row, end), mod1(col + 1, end)] + edge(env, NORTH, row - 1, col) * + twistdual(corner(env, NORTHEAST, row - 1, col + 1), 1) + E_northeast = edge(env, EAST, row, col + 1) E_southeast = - env.edges[EAST, mod1(row + 1, end), mod1(col + 1, end)] * - twistdual(env.corners[SOUTHEAST, mod1(row + 2, end), mod1(col + 1, end)], 1) + edge(env, EAST, row + 1, col + 1) * + twistdual(corner(env, SOUTHEAST, row + 2, col + 1), 1) E_south = - env.edges[SOUTH, mod1(row + 2, end), mod1(col, end)] * - twistdual(env.corners[SOUTHWEST, mod1(row + 2, end), mod1(col - 1, end)], 1) - E_southwest = env.edges[WEST, mod1(row + 1, end), mod1(col - 1, end)] + edge(env, SOUTH, row + 2, col) * + twistdual(corner(env, SOUTHWEST, row + 2, col - 1), 1) + E_southwest = edge(env, WEST, row + 1, col - 1) E_northwest = - env.edges[WEST, mod1(row, end), mod1(col - 1, end)] * - twistdual(env.corners[NORTHWEST, mod1(row - 1, end), mod1(col - 1, end)], 1) + edge(env, WEST, row, col - 1) * + twistdual(corner(env, NORTHWEST, row - 1, col - 1), 1) @tensor EE_NW[χW χNE DNWt DNt; DNWb DNb] := E_northwest[χW DNWt DNWb; χNW] * E_north[χNW DNt DNb; χNE] @@ -500,25 +500,25 @@ function reduced_densitymatrix1x2( row, col = Tuple(ind) # Unpack variables and absorb corners - A_west = ket[mod1(row, end), mod1(col, end)] - Ā_west = bra[mod1(row, end), mod1(col, end)] - A_east = ket[mod1(row, end), mod1(col + 1, end)] - Ā_east = bra[mod1(row, end), mod1(col + 1, end)] + A_west = ket[row, col] + Ā_west = bra[row, col] + A_east = ket[row, col + 1] + Ā_east = bra[row, col + 1] - E_northwest = env.edges[NORTH, mod1(row - 1, end), mod1(col, end)] + E_northwest = edge(env, NORTH, row - 1, col) E_northeast = - env.edges[NORTH, mod1(row - 1, end), mod1(col + 1, end)] * - twistdual(env.corners[NORTHEAST, mod1(row - 1, end), mod1(col + 2, end)], 1) + edge(env, NORTH, row - 1, col + 1) * + twistdual(corner(env, NORTHEAST, row - 1, col + 2), 1) E_east = - env.edges[EAST, mod1(row, end), mod1(col + 2, end)] * - twistdual(env.corners[SOUTHEAST, mod1(row + 1, end), mod1(col + 2, end)], 1) - E_southeast = env.edges[SOUTH, mod1(row + 1, end), mod1(col + 1, end)] + edge(env, EAST, row, col + 2) * + twistdual(corner(env, SOUTHEAST, row + 1, col + 2), 1) + E_southeast = edge(env, SOUTH, row + 1, col + 1) E_southwest = - env.edges[SOUTH, mod1(row + 1, end), mod1(col, end)] * - twistdual(env.corners[SOUTHWEST, mod1(row + 1, end), mod1(col - 1, end)], 1) + edge(env, SOUTH, row + 1, col) * + twistdual(corner(env, SOUTHWEST, row + 1, col - 1), 1) E_west = - env.edges[WEST, mod1(row, end), mod1(col - 1, end)] * - twistdual(env.corners[NORTHWEST, mod1(row - 1, end), mod1(col - 1, end)], 1) + edge(env, WEST, row, col - 1) * + twistdual(corner(env, NORTHWEST, row - 1, col - 1), 1) @tensor EE_SW[χS χNW DSWt DWt; DSWb DWb] := E_southwest[χS DSWt DSWb; χSW] * E_west[χSW DWt DWb; χNW] diff --git a/src/algorithms/correlators.jl b/src/algorithms/correlators.jl index 674644aab..0c95723b6 100644 --- a/src/algorithms/correlators.jl +++ b/src/algorithms/correlators.jl @@ -30,8 +30,8 @@ function correlator_horizontal( for (k, j) in enumerate(js) # transfer until left of site j while j > i - Etop = env.edges[NORTH, _prev(i[1], end), mod1(i[2], end)] - Ebot = env.edges[SOUTH, _next(i[1], end), mod1(i[2], end)] + Etop = edge(env, NORTH, i[1] - 1, i[2]) + Ebot = edge(env, SOUTH, i[1] + 1, i[2]) sandwich = ( ket[mod1(i[1], end), mod1(i[2], end)], bra[mod1(i[1], end), mod1(i[2], end)], ) @@ -43,11 +43,9 @@ function correlator_horizontal( # compute overlap with operator numerator = end_correlator_numerator(j, Vo, bra, O[2], ket, env) # transfer right of site j - Etop = env.edges[NORTH, _prev(i[1], end), mod1(i[2], end)] - Ebot = env.edges[SOUTH, _next(i[1], end), mod1(i[2], end)] - sandwich = ( - ket[mod1(i[1], end), mod1(i[2], end)], bra[mod1(i[1], end), mod1(i[2], end)], - ) + Etop = edge(env, NORTH, i[1] - 1, i[2]) + Ebot = edge(env, SOUTH, i[1] + 1, i[2]) + sandwich = (ket[i[1], i[2]], bra[i[1], i[2]]) T = edge_transfermatrix(Etop, sandwich, Ebot) if k < length(js) Vo = Vo * T @@ -137,9 +135,9 @@ function correlator_horizontal( for (k, j) in enumerate(js) # transfer until left of site j while j > i - Etop = env.edges[NORTH, _prev(i[1], end), mod1(i[2], end)] - Omid = trace_physicalspaces(ρ[mod1(i[1], end), mod1(i[2], end)]) - Ebot = env.edges[SOUTH, _next(i[1], end), mod1(i[2], end)] + Etop = edge(env, NORTH, i[1] - 1, i[2]) + Omid = trace_physicalspaces(ρ[i[1], i[2]]) + Ebot = edge(env, SOUTH, i[1] + 1, i[2]) T = edge_transfermatrix(Etop, Omid, Ebot) Vo = Vo * T Vn = Vn * T @@ -148,9 +146,9 @@ function correlator_horizontal( # compute overlap with operator numerator = end_correlator_numerator(j, Vo, ρ, O[2], env) # transfer right of site j - Etop = env.edges[NORTH, _prev(i[1], end), mod1(i[2], end)] + Etop = edge(env, NORTH, i[1] - 1, i[2]) Omid = trace_physicalspaces(ρ[mod1(i[1], end), mod1(i[2], end)]) - Ebot = env.edges[SOUTH, _next(i[1], end), mod1(i[2], end)] + Ebot = edge(env, SOUTH, i[1] + 1, i[2]) T = edge_transfermatrix(Etop, Omid, Ebot) if k < length(js) Vo = Vo * T diff --git a/src/algorithms/ctmrg/sequential.jl b/src/algorithms/ctmrg/sequential.jl index 69d45a34a..55ddacfef 100644 --- a/src/algorithms/ctmrg/sequential.jl +++ b/src/algorithms/ctmrg/sequential.jl @@ -94,7 +94,7 @@ function sequential_projectors( ) _, r, c = coordinate r′ = _prev(r, size(env, 2)) - trunc = truncation_strategy(alg, env.edges[WEST, r′, c]) + trunc = truncation_strategy(alg, edge(env, WEST, r′, c)) alg´ = _set_decomposition_truncation(alg, trunc) Q1 = TensorMap(EnlargedCorner(network, env, (SOUTHWEST, r, c))) Q2 = TensorMap(EnlargedCorner(network, env, (NORTHWEST, r′, c))) @@ -129,11 +129,11 @@ function renormalize_sequentially(col::Int, projectors, network, env) for (dir, r, c) in eachcoordinate(network, 1:4) (c == col && dir in [SOUTHWEST, NORTHWEST]) && continue - corners[dir, r, c] = env.corners[dir, r, c] + corners[dir, r, c] = corner(env, dir, r, c) end for (dir, r, c) in eachcoordinate(network, 1:4) (c == col && dir == WEST) && continue - edges[dir, r, c] = env.edges[dir, r, c] + edges[dir, r, c] = edge(env, dir, r, c) end # Apply projectors to renormalize corners and edge diff --git a/src/algorithms/ctmrg/sparse_environments.jl b/src/algorithms/ctmrg/sparse_environments.jl index 5d5526fd9..40990f4a3 100644 --- a/src/algorithms/ctmrg/sparse_environments.jl +++ b/src/algorithms/ctmrg/sparse_environments.jl @@ -26,33 +26,33 @@ function EnlargedCorner(network::InfiniteSquareNetwork, env, coordinates) dir, r, c = coordinates if dir == NORTHWEST return EnlargedCorner( - env.corners[NORTHWEST, _prev(r, end), _prev(c, end)], - env.edges[WEST, r, _prev(c, end)], - env.edges[NORTH, _prev(r, end), c], + corner(env, NORTHWEST, r - 1, c - 1), + edge(env, WEST, r, c - 1), + edge(env, NORTH, r - 1, c), network[r, c], dir, ) elseif dir == NORTHEAST return EnlargedCorner( - env.corners[NORTHEAST, _prev(r, end), _next(c, end)], - env.edges[NORTH, _prev(r, end), c], - env.edges[EAST, r, _next(c, end)], + corner(env, NORTHEAST, r - 1, c + 1), + edge(env, NORTH, r - 1, c), + edge(env, EAST, r, c + 1), network[r, c], dir, ) elseif dir == SOUTHEAST return EnlargedCorner( - env.corners[SOUTHEAST, _next(r, end), _next(c, end)], - env.edges[EAST, r, _next(c, end)], - env.edges[SOUTH, _next(r, end), c], + corner(env, SOUTHEAST, r + 1, c + 1), + edge(env, EAST, r, c + 1), + edge(env, SOUTH, r + 1, c), network[r, c], dir, ) elseif dir == SOUTHWEST return EnlargedCorner( - env.corners[SOUTHWEST, _next(r, end), _prev(c, end)], - env.edges[SOUTH, _next(r, end), c], - env.edges[WEST, r, _prev(c, end)], + corner(env, SOUTHWEST, r + 1, c - 1), + edge(env, SOUTH, r + 1, c), + edge(env, WEST, r, c - 1), network[r, c], dir, ) @@ -424,11 +424,9 @@ struct ColumnEnlargedCorner{TC, TE} end function ColumnEnlargedCorner(env::CTMRGEnv, coordinates) dir, row, col = coordinates - Nc = size(env, 3) if dir == NORTHWEST - cm1 = _prev(col, Nc) return ColumnEnlargedCorner( - env.corners[dir, row, cm1], env.edges[NORTH, row, col], dir + corner(env, dir, row, col - 1), edge(env, NORTH, row, col), dir ) else error("Not implemented.") diff --git a/src/algorithms/time_evolution/get_cluster.jl b/src/algorithms/time_evolution/get_cluster.jl index 81ab6e599..92b1c1480 100644 --- a/src/algorithms/time_evolution/get_cluster.jl +++ b/src/algorithms/time_evolution/get_cluster.jl @@ -99,7 +99,6 @@ Obtain the cluster `Ms` along the (open) path `sites` in `state`. function _get_cluster( state::InfiniteState, sites::Vector{CartesianIndex{2}} ) - Nr, Nc = size(state) # number of sites Ns = length(sites) # number of physical axes @@ -133,9 +132,7 @@ function _get_cluster( end return _mpo_perm(out_ax + Np, in_ax + Np, Nax) end - Ms = map(sites) do site - return state[CartesianIndex(mod1(site[1], Nr), mod1(site[2], Nc))] - end + Ms = map(site -> state[site], sites) return Ms, open_vaxs, perms end diff --git a/src/algorithms/time_evolution/simpleupdate.jl b/src/algorithms/time_evolution/simpleupdate.jl index 9fab38443..c4ab7f6d3 100644 --- a/src/algorithms/time_evolution/simpleupdate.jl +++ b/src/algorithms/time_evolution/simpleupdate.jl @@ -107,9 +107,7 @@ function _su_iter!( B = _bond_rotation(B, dir, rev; inv = true) rev && (s = transpose(s)) # remove environment weights - siteA, siteB = map(sites) do site - return CartesianIndex(mod1(site[1], Nr), mod1(site[2], Nc)) - end + siteA, siteB = sites A = absorb_weight(A, env, siteA[1], siteA[2], open_vaxs[1]; inv = true) B = absorb_weight(B, env, siteB[1], siteB[2], open_vaxs[2]; inv = true) # update tensor dict and weight on current bond @@ -133,8 +131,7 @@ function su_iter( # 1-site gate # TODO: special treatment for bipartite state site = sites[1] - r, c = mod1(site[1], Nr), mod1(site[2], Nc) - state2[r, c] = _apply_sitegate(state2[r, c], gate; alg.purified) + state2[site] = _apply_sitegate(state2[site], gate; alg.purified) elseif length(sites) == 2 (d, r, c), = _nn_bondrev(sites..., (Nr, Nc)) alg.bipartite && r > 1 && continue @@ -142,15 +139,13 @@ function su_iter( ϵ = max(ϵ, ϵ′) (!alg.bipartite) && continue if d == 1 - rp1, cp1 = _next(r, Nr), _next(c, Nc) - state2[rp1, cp1] = copy(state2[r, c]) - state2[rp1, c] = copy(state2[r, cp1]) - env2[1, rp1, cp1] = copy(env2[1, r, c]) + state2[r + 1, c + 1] = copy(state2[r, c]) + state2[r + 1, c] = copy(state2[r, c + 1]) + env2[1, r + 1, c + 1] = copy(env2[1, r, c]) else - rm1, cm1 = _prev(r, Nr), _prev(c, Nc) - state2[rm1, cm1] = copy(state2[r, c]) - state2[r, cm1] = copy(state2[rm1, c]) - env2[2, rm1, cm1] = copy(env2[2, r, c]) + state2[r - 1, c - 1] = copy(state2[r, c]) + state2[r, c - 1] = copy(state2[r - 1, c]) + env2[2, r - 1, c - 1] = copy(env2[2, r, c]) end else # N-site MPO gate (N ≥ 2) diff --git a/src/algorithms/time_evolution/simpleupdate3site.jl b/src/algorithms/time_evolution/simpleupdate3site.jl index 3fe31c98f..28ca01339 100644 --- a/src/algorithms/time_evolution/simpleupdate3site.jl +++ b/src/algorithms/time_evolution/simpleupdate3site.jl @@ -58,10 +58,9 @@ function _su_iter!( end # update state tensors for (M, s, invperm, vaxs) in zip(Ms, sites, invperms, open_vaxs) - s′ = CartesianIndex(mod1(s[1], Nr), mod1(s[2], Nc)) # restore original axes order and remove weights on open axes of the cluster - M′ = absorb_weight(permute(M, invperm), env, s′, vaxs; inv = true) - state[s′] = normalize!(M′, Inf) + M′ = absorb_weight(permute(M, invperm), env, s, vaxs; inv = true) + state[s] = normalize!(M′, Inf) end return ϵ end diff --git a/src/algorithms/toolbox.jl b/src/algorithms/toolbox.jl index d9b7802a3..c532d3e1d 100644 --- a/src/algorithms/toolbox.jl +++ b/src/algorithms/toolbox.jl @@ -96,12 +96,12 @@ Contract around a single site `ind` of a square network using a given CTMRG envi function _contract_site(ind::Tuple{Int, Int}, network, env::CTMRGEnv) r, c = ind return _contract_site( - env.corners[NORTHWEST, _prev(r, end), _prev(c, end)], - env.corners[NORTHEAST, _prev(r, end), _next(c, end)], - env.corners[SOUTHEAST, _next(r, end), _next(c, end)], - env.corners[SOUTHWEST, _next(r, end), _prev(c, end)], - env.edges[NORTH, _prev(r, end), c], env.edges[EAST, r, _next(c, end)], - env.edges[SOUTH, _next(r, end), c], env.edges[WEST, r, _prev(c, end)], + corner(env, NORTHWEST, r - 1, c - 1), + corner(env, NORTHEAST, r - 1, c + 1), + corner(env, SOUTHEAST, r + 1, c + 1), + corner(env, SOUTHWEST, r + 1, c - 1), + edge(env, NORTH, r - 1, c), edge(env, EAST, r, c + 1), + edge(env, SOUTH, r + 1, c), edge(env, WEST, r, c - 1), network[r, c], ) end @@ -147,10 +147,10 @@ environment `env`. """ function _contract_corners(ind::Tuple{Int, Int}, env::CTMRGEnv) r, c = ind - C_NW = env.corners[NORTHWEST, _prev(r, end), _prev(c, end)] - C_NE = env.corners[NORTHEAST, _prev(r, end), c] - C_SE = env.corners[SOUTHEAST, r, c] - C_SW = env.corners[SOUTHWEST, r, _prev(c, end)] + C_NW = corner(env, NORTHWEST, r - 1, c - 1) + C_NE = corner(env, NORTHEAST, r - 1, c) + C_SE = corner(env, SOUTHEAST, r, c) + C_SW = corner(env, SOUTHWEST, r, c - 1) return @tensor C_NW[1; 2] * C_NE[2; 3] * C_SE[3; 4] * C_SW[4; 1] end @@ -163,12 +163,12 @@ CTMRG environment `env`. function _contract_vertical_edges(ind::Tuple{Int, Int}, env::CTMRGEnv) r, c = ind return _contract_vertical_edges( - env.corners[NORTHWEST, _prev(r, end), _prev(c, end)], - env.corners[NORTHEAST, _prev(r, end), c], - env.corners[SOUTHEAST, _next(r, end), c], - env.corners[SOUTHWEST, _next(r, end), _prev(c, end)], - env.edges[EAST, r, c], - env.edges[WEST, r, _prev(c, end)], + corner(env, NORTHWEST, r - 1, c - 1), + corner(env, NORTHEAST, r - 1, c), + corner(env, SOUTHEAST, r + 1, c), + corner(env, SOUTHWEST, r + 1, c - 1), + edge(env, EAST, r, c), + edge(env, WEST, r, c - 1), ) end @generated function _contract_vertical_edges( @@ -206,12 +206,12 @@ CTMRG environment `env`. function _contract_horizontal_edges(ind::Tuple{Int, Int}, env::CTMRGEnv) r, c = ind return _contract_horizontal_edges( - env.corners[NORTHWEST, _prev(r, end), _prev(c, end)], - env.corners[NORTHEAST, _prev(r, end), _next(c, end)], - env.corners[SOUTHEAST, r, _next(c, end)], - env.corners[SOUTHWEST, r, _prev(c, end)], - env.edges[NORTH, _prev(r, end), c], - env.edges[SOUTH, r, c], + corner(env, NORTHWEST, r - 1, c - 1), + corner(env, NORTHEAST, r - 1, c + 1), + corner(env, SOUTHEAST, r, c + 1), + corner(env, SOUTHWEST, r, c - 1), + edge(env, NORTH, r - 1, c), + edge(env, SOUTH, r, c), ) end @generated function _contract_horizontal_edges( @@ -380,12 +380,12 @@ function contract_local_tensor( ) where {C} r, c = inds return _contract_site( - env.corners[NORTHWEST, _prev(r, end), _prev(c, end)], - env.corners[NORTHEAST, _prev(r, end), _next(c, end)], - env.corners[SOUTHEAST, _next(r, end), _next(c, end)], - env.corners[SOUTHWEST, _next(r, end), _prev(c, end)], - env.edges[NORTH, _prev(r, end), c], env.edges[EAST, r, _next(c, end)], - env.edges[SOUTH, _next(r, end), c], env.edges[WEST, r, _prev(c, end)], + corner(env, NORTHWEST, r - 1, c - 1), + corner(env, NORTHEAST, r - 1, c + 1), + corner(env, SOUTHEAST, r + 1, c + 1), + corner(env, SOUTHWEST, r + 1, c - 1), + edge(env, NORTH, r - 1, c), edge(env, EAST, r, c + 1), + edge(env, SOUTH, r + 1, c), edge(env, WEST, r, c - 1), O, ) end @@ -405,12 +405,12 @@ function contract_local_tensor( r, c, h = ind sandwich´ = Base.setindex(network[r, c], O, h + 2) return _contract_site( - env.corners[NORTHWEST, _prev(r, end), _prev(c, end)], - env.corners[NORTHEAST, _prev(r, end), _next(c, end)], - env.corners[SOUTHEAST, _next(r, end), _next(c, end)], - env.corners[SOUTHWEST, _next(r, end), _prev(c, end)], - env.edges[NORTH, _prev(r, end), c], env.edges[EAST, r, _next(c, end)], - env.edges[SOUTH, _next(r, end), c], env.edges[WEST, r, _prev(c, end)], + corner(env, NORTHWEST, r - 1, c - 1), + corner(env, NORTHEAST, r - 1, c + 1), + corner(env, SOUTHEAST, r + 1, c + 1), + corner(env, SOUTHWEST, r + 1, c - 1), + edge(env, NORTH, r - 1, c), edge(env, EAST, r, c + 1), + edge(env, SOUTH, r + 1, c), edge(env, WEST, r, c - 1), sandwich´, ) end From f730ef0f6bf55afd0e50ff8ef789493c3c7cd5cd Mon Sep 17 00:00:00 2001 From: lkdvos Date: Thu, 7 May 2026 09:36:47 -0400 Subject: [PATCH 04/12] bypass non-scalar indexing calls --- src/environments/suweight.jl | 12 ++++++------ src/operators/transfermatrix.jl | 4 ++-- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/src/environments/suweight.jl b/src/environments/suweight.jl index 1ac771530..205ed5992 100644 --- a/src/environments/suweight.jl +++ b/src/environments/suweight.jl @@ -304,18 +304,18 @@ function _rot180_wts_y(wts_y::AbstractMatrix{<:PEPSWeight}) end function Base.rotl90(wts::SUWeight) - wts_y = _rotl90_wts_x(wts[1, :, :]) - wts_x = _rotl90_wts_y(wts[2, :, :]) + wts_y = _rotl90_wts_x(wts.data[1, :, :]) + wts_x = _rotl90_wts_y(wts.data[2, :, :]) return SUWeight(wts_x, wts_y) end function Base.rotr90(wts::SUWeight) - wts_y = _rotr90_wts_x(wts[1, :, :]) - wts_x = _rotr90_wts_y(wts[2, :, :]) + wts_y = _rotr90_wts_x(wts.data[1, :, :]) + wts_x = _rotr90_wts_y(wts.data[2, :, :]) return SUWeight(wts_x, wts_y) end function Base.rot180(wts::SUWeight) - wts_x = _rot180_wts_x(wts[1, :, :]) - wts_y = _rot180_wts_y(wts[2, :, :]) + wts_x = _rot180_wts_x(wts.data[1, :, :]) + wts_y = _rot180_wts_y(wts.data[2, :, :]) return SUWeight(wts_x, wts_y) end diff --git a/src/operators/transfermatrix.jl b/src/operators/transfermatrix.jl index 7520d69a2..a033817a4 100644 --- a/src/operators/transfermatrix.jl +++ b/src/operators/transfermatrix.jl @@ -31,7 +31,7 @@ the direction `dir` faces north, after which its `row`th row from the north is s """ function InfiniteTransferPEPS(T::InfinitePEPS, dir, row) T = rotate_north(T, dir) - return InfiniteTransferPEPS(PeriodicArray(T[row, :])) + return InfiniteTransferPEPS(PeriodicArray(unitcell(T)[mod1(row, end), :])) end """ @@ -90,7 +90,7 @@ north is selected. function InfiniteTransferPEPO(T::InfinitePEPS, O::InfinitePEPO, dir, row) T = rotate_north(T, dir) O = rotate_north(O, dir) - return InfiniteTransferPEPO(PeriodicArray(T[row, :]), PeriodicArray(O[row, :, :])) + return InfiniteTransferPEPO(PeriodicArray(unitcell(T)[mod1(row, end), :]), PeriodicArray(unitcell(O)[mod1(row, end), :, :])) end """ From 0eec4b229700f2aa20ab130105a9ba10c367911f Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Fri, 8 May 2026 11:44:49 +0800 Subject: [PATCH 05/12] Move indexing utils to a separate file --- src/PEPSKit.jl | 1 + src/utility/indexing.jl | 83 ++++++++++++++++++++++++++++++++++++++++ src/utility/util.jl | 84 ----------------------------------------- 3 files changed, 84 insertions(+), 84 deletions(-) create mode 100644 src/utility/indexing.jl diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl index d669f02e1..9fd2b396f 100644 --- a/src/PEPSKit.jl +++ b/src/PEPSKit.jl @@ -41,6 +41,7 @@ using DocStringExtensions include("Defaults.jl") # Include first to allow for docstring interpolation with Defaults values include("utility/util.jl") +include("utility/indexing.jl") include("utility/diffable_threads.jl") include("utility/eigh.jl") include("utility/svd.jl") diff --git a/src/utility/indexing.jl b/src/utility/indexing.jl new file mode 100644 index 000000000..7e1f7a942 --- /dev/null +++ b/src/utility/indexing.jl @@ -0,0 +1,83 @@ +# Get next and previous directional CTMRG environment index, respecting periodicity +_next(i, total) = mod1(i + 1, total) +_prev(i, total) = mod1(i - 1, total) + +@inline _periodic_inds(data, J::NTuple{N, Int}) where {N} = + ntuple(i -> mod1(J[i], size(data, i)), Val(N)) + +""" + periodic_getindex(A, data, I::Tuple) + +Periodic indexing of `A[I...]` backed by `data`. +For `I::Vararg{Int}`, this yields `data[mod1.(I, size(data))...]`. +Mixed index arguments (`Int` / `CartesianIndex`) are supported by flattening with `Base.to_indices`. + +!!! warning + Linear indexing schemes are not supported. +""" +@inline function periodic_getindex(A, data, I::Tuple) + J = to_indices(data, I) + return _periodic_getindex_dispatch(A, data, J) +end + +# Dispatch helper +_periodic_getindex_dispatch(A, data, J::Tuple{Vararg{Int}}) = @inbounds data[_periodic_inds(data, J)...] +_periodic_getindex_dispatch(A, data, J::Tuple) = throw(ArgumentError("Invalid periodic indexing type $(typeof(J))")) + +""" + periodic_setindex!(A, data, v, I::Tuple) + +Periodic indexing assignment `A[I...] = v` backed by `data`, returning `v` so that the +syntax `A[I...] = v` evaluates to `v` in line with Base Julia conventions. +For `I::Vararg{Int}`, this yields `data[mod1.(I, size(data))...] = v`. +Mixed index arguments (`Int` / `CartesianIndex`) are supported by flattening with `Base.to_indices`. + +!!! warning + Linear indexing schemes are not supported. +""" +@inline function periodic_setindex!(A, data, v, I::Tuple) + J = to_indices(data, I) + return _periodic_setindex_dispatch!(A, data, v, J) +end + +# Dispatch helper +function _periodic_setindex_dispatch!(A, data, v, J::Tuple{Vararg{Int}}) + @inbounds data[_periodic_inds(data, J)...] = v + return v +end +_periodic_setindex_dispatch!(_A, data, v, J::Tuple) = + throw(ArgumentError("Invalid periodic indexing type $(typeof(J))")) + +# Get next and previous coordinate (direction, row, column), given a direction and going around the environment clockwise +function _next_coordinate((dir, row, col), rowsize, colsize) + if dir == 1 + return (_next(dir, 4), row, _next(col, colsize)) + elseif dir == 2 + return (_next(dir, 4), _next(row, rowsize), col) + elseif dir == 3 + return (_next(dir, 4), row, _prev(col, colsize)) + elseif dir == 4 + return (_next(dir, 4), _prev(row, rowsize), col) + end +end +function _prev_coordinate((dir, row, col), rowsize, colsize) + if dir == 1 + return (_prev(dir, 4), _next(row, rowsize), col) + elseif dir == 2 + return (_prev(dir, 4), row, _prev(col, colsize)) + elseif dir == 3 + return (_prev(dir, 4), _prev(row, rowsize), col) + elseif dir == 4 + return (_prev(dir, 4), row, _next(col, colsize)) + end +end + +# iterator over each coordinates +""" + eachcoordinate(x, [dirs=1:4]) + +Enumerate all (dir, row, col) pairs. +""" +function eachcoordinate end + +@non_differentiable eachcoordinate(args...) diff --git a/src/utility/util.jl b/src/utility/util.jl index d6c11370d..2c58f8801 100644 --- a/src/utility/util.jl +++ b/src/utility/util.jl @@ -1,87 +1,3 @@ -# Get next and previous directional CTMRG environment index, respecting periodicity -_next(i, total) = mod1(i + 1, total) -_prev(i, total) = mod1(i - 1, total) - -@inline _periodic_inds(data, J::NTuple{N, Int}) where {N} = - ntuple(i -> mod1(J[i], size(data, i)), Val(N)) - -""" - periodic_getindex(A, data, I::Tuple) - -Periodic indexing of `A[I...]` backed by `data`. -For `I::Vararg{Int}`, this yields `data[mod1.(I, size(data))...]`. -Mixed index arguments (`Int` / `CartesianIndex`) are supported by flattening with `Base.to_indices`. - -!!! warning - Linear indexing schemes are not supported. -""" -@inline function periodic_getindex(A, data, I::Tuple) - J = to_indices(data, I) - return _periodic_getindex_dispatch(A, data, J) -end - -# Dispatch helper -_periodic_getindex_dispatch(A, data, J::Tuple{Vararg{Int}}) = @inbounds data[_periodic_inds(data, J)...] -_periodic_getindex_dispatch(A, data, J::Tuple) = throw(ArgumentError("Invalid periodic indexing type $(typeof(J))")) - -""" - periodic_setindex!(A, data, v, I::Tuple) - -Periodic indexing assignment `A[I...] = v` backed by `data`, returning `v` so that the -syntax `A[I...] = v` evaluates to `v` in line with Base Julia conventions. -For `I::Vararg{Int}`, this yields `data[mod1.(I, size(data))...] = v`. -Mixed index arguments (`Int` / `CartesianIndex`) are supported by flattening with `Base.to_indices`. - -!!! warning - Linear indexing schemes are not supported. -""" -@inline function periodic_setindex!(A, data, v, I::Tuple) - J = to_indices(data, I) - return _periodic_setindex_dispatch!(A, data, v, J) -end - -# Dispatch helper -function _periodic_setindex_dispatch!(A, data, v, J::Tuple{Vararg{Int}}) - @inbounds data[_periodic_inds(data, J)...] = v - return v -end -_periodic_setindex_dispatch!(_A, data, v, J::Tuple) = - throw(ArgumentError("Invalid periodic indexing type $(typeof(J))")) - -# Get next and previous coordinate (direction, row, column), given a direction and going around the environment clockwise -function _next_coordinate((dir, row, col), rowsize, colsize) - if dir == 1 - return (_next(dir, 4), row, _next(col, colsize)) - elseif dir == 2 - return (_next(dir, 4), _next(row, rowsize), col) - elseif dir == 3 - return (_next(dir, 4), row, _prev(col, colsize)) - elseif dir == 4 - return (_next(dir, 4), _prev(row, rowsize), col) - end -end -function _prev_coordinate((dir, row, col), rowsize, colsize) - if dir == 1 - return (_prev(dir, 4), _next(row, rowsize), col) - elseif dir == 2 - return (_prev(dir, 4), row, _prev(col, colsize)) - elseif dir == 3 - return (_prev(dir, 4), _prev(row, rowsize), col) - elseif dir == 4 - return (_prev(dir, 4), row, _next(col, colsize)) - end -end - -# iterator over each coordinates -""" - eachcoordinate(x, [dirs=1:4]) - -Enumerate all (dir, row, col) pairs. -""" -function eachcoordinate end - -@non_differentiable eachcoordinate(args...) - # Element-wise multiplication of TensorMaps respecting block structure function _elementwise_mult(a₁::AbstractTensorMap, a₂::AbstractTensorMap) dst = similar(a₁) From e8c3021b7420aace0dde7f40c2d83b7b89eb0491 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Fri, 8 May 2026 11:54:42 +0800 Subject: [PATCH 06/12] Improve SUWeight to CTMRGEnv conversion --- src/environments/suweight.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/environments/suweight.jl b/src/environments/suweight.jl index 205ed5992..440d3ba33 100644 --- a/src/environments/suweight.jl +++ b/src/environments/suweight.jl @@ -332,9 +332,9 @@ function CTMRGEnv(wts::SUWeight) V_env = oneunit(spacetype(wts)) edges = map(Iterators.product(1:4, 1:Nr, 1:Nc)) do (d, r, c) wt_idx = if d == NORTH - CartesianIndex(2, _next(r, Nr), c) + CartesianIndex(2, r + 1, c) elseif d == EAST - CartesianIndex(1, r, _prev(c, Nc)) + CartesianIndex(1, r, c - 1) elseif d == SOUTH CartesianIndex(2, r, c) else # WEST From 7d1bb9d36b5216ef8601d9474ce1b19566081da5 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Fri, 8 May 2026 12:04:44 +0800 Subject: [PATCH 07/12] Periodic indexing for SiteDependentTruncation --- src/algorithms/truncation/truncationschemes.jl | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/src/algorithms/truncation/truncationschemes.jl b/src/algorithms/truncation/truncationschemes.jl index 7cad1dfd6..cdddbf41d 100644 --- a/src/algorithms/truncation/truncationschemes.jl +++ b/src/algorithms/truncation/truncationschemes.jl @@ -32,7 +32,8 @@ struct SiteDependentTruncation{T <: TruncationStrategy} <: TruncationStrategy end end -Base.getindex(trunc::SiteDependentTruncation, args...) = Base.getindex(trunc.truncs, args...) +Base.@propagate_inbounds Base.getindex(trunc::SiteDependentTruncation, I...) = + periodic_getindex(trunc, trunc.truncs, I) # TODO: _is_bipartite(trunc::SiteDependentTruncation) @@ -113,20 +114,20 @@ function _rot180_trunc_y(trunc_y::AbstractMatrix{<:TruncationStrategy}) end function Base.rotl90(trunc::SiteDependentTruncation) - trunc_y = _rotl90_trunc_x(trunc[1, :, :]) - trunc_x = _rotl90_trunc_y(trunc[2, :, :]) + trunc_y = _rotl90_trunc_x(trunc.truncs[1, :, :]) + trunc_x = _rotl90_trunc_y(trunc.truncs[2, :, :]) trunc = stack((trunc_x, trunc_y); dims = 1) return SiteDependentTruncation(trunc) end function Base.rotr90(trunc::SiteDependentTruncation) - trunc_y = _rotr90_trunc_x(trunc[1, :, :]) - trunc_x = _rotr90_trunc_y(trunc[2, :, :]) + trunc_y = _rotr90_trunc_x(trunc.truncs[1, :, :]) + trunc_x = _rotr90_trunc_y(trunc.truncs[2, :, :]) trunc = stack((trunc_x, trunc_y); dims = 1) return SiteDependentTruncation(trunc) end function Base.rot180(trunc::SiteDependentTruncation) - trunc_x = _rot180_trunc_x(trunc[1, :, :]) - trunc_y = _rot180_trunc_y(trunc[2, :, :]) + trunc_x = _rot180_trunc_x(trunc.truncs[1, :, :]) + trunc_y = _rot180_trunc_y(trunc.truncs[2, :, :]) trunc = stack((trunc_x, trunc_y); dims = 1) return SiteDependentTruncation(trunc) end From 638293fd7409fde5aea840483d7e776a65ae27a4 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Fri, 8 May 2026 12:21:00 +0800 Subject: [PATCH 08/12] Improve _is_bipartite --- src/environments/bp_environments.jl | 2 +- src/environments/suweight.jl | 2 +- src/operators/infinitepepo.jl | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/environments/bp_environments.jl b/src/environments/bp_environments.jl index 30d6cdd70..affb71e4b 100644 --- a/src/environments/bp_environments.jl +++ b/src/environments/bp_environments.jl @@ -165,7 +165,7 @@ end function _is_bipartite(env::BPEnv) (size(env, 2) == size(env, 3) == 2) || (return false) for (d, c) in Iterators.product(axes(env, 1), axes(env, 3)) - (env[d, 1, c] == env[d, 2, _next(c, 2)]) || (return false) + (env[d, 1, c] == env[d, 2, c + 1]) || (return false) end return true end diff --git a/src/environments/suweight.jl b/src/environments/suweight.jl index 440d3ba33..dbf7b0059 100644 --- a/src/environments/suweight.jl +++ b/src/environments/suweight.jl @@ -144,7 +144,7 @@ TensorKit.sectortype(::Type{<:SUWeight{T}}) where {T} = sectortype(spacetype(T)) function _is_bipartite(wts::SUWeight) (size(wts, 2) == size(wts, 3) == 2) || (return false) for (d, c) in Iterators.product(1:2, 1:2) - (wts[d, 1, c] == wts[d, 2, _next(c, 2)]) || (return false) + (wts[d, 1, c] == wts[d, 2, c + 1]) || (return false) end return true end diff --git a/src/operators/infinitepepo.jl b/src/operators/infinitepepo.jl index a24af3901..a180a953c 100644 --- a/src/operators/infinitepepo.jl +++ b/src/operators/infinitepepo.jl @@ -202,7 +202,7 @@ end function _is_bipartite(psi::InfiniteState) (size(psi, 1) == size(psi, 2) == 2) || (return false) for (c, h) in Iterators.product(1:2, 1:size(psi, 3)) - (psi[1, c, h] == psi[2, _next(c, 2), h]) || (return false) + (psi[1, c, h] == psi[2, c + 1, h]) || (return false) end return true end From 07cedc3030a1bc1bf967cf651022775e0a24b9b8 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Fri, 8 May 2026 12:21:09 +0800 Subject: [PATCH 09/12] Improve BPEnv indexing --- src/algorithms/bp/beliefpropagation.jl | 8 ++++---- src/algorithms/bp/gaugefix.jl | 2 +- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/algorithms/bp/beliefpropagation.jl b/src/algorithms/bp/beliefpropagation.jl index 0aa900912..5f64ed0b7 100644 --- a/src/algorithms/bp/beliefpropagation.jl +++ b/src/algorithms/bp/beliefpropagation.jl @@ -104,10 +104,10 @@ function update_message(I::CartesianIndex{3}, network::InfiniteSquareNetwork, en (1 <= dir <= 4) || throw(ArgumentError("Invalid direction $dir")) A = network[row, col] - dir == SOUTH || (M_north = env[NORTH, _prev(row, end), col]) - dir == WEST || (M_east = env[EAST, row, _next(col, end)]) - dir == NORTH || (M_south = env[SOUTH, _next(row, end), col]) - dir == EAST || (M_west = env[WEST, row, _prev(col, end)]) + dir == SOUTH || (M_north = env[NORTH, row - 1, col]) + dir == WEST || (M_east = env[EAST, row, col + 1]) + dir == NORTH || (M_south = env[SOUTH, row + 1, col]) + dir == EAST || (M_west = env[WEST, row, col - 1]) return if dir == NORTH contract_north_message(A, M_west, M_north, M_east) diff --git a/src/algorithms/bp/gaugefix.jl b/src/algorithms/bp/gaugefix.jl index 56e5e9f30..74b19409f 100644 --- a/src/algorithms/bp/gaugefix.jl +++ b/src/algorithms/bp/gaugefix.jl @@ -52,7 +52,7 @@ end function _sqrt_bp_messages(I::CartesianIndex{3}, env::BPEnv) dir, row, col = Tuple(I) @assert dir == NORTH || dir == EAST - M12 = env[dir, dir == NORTH ? _prev(row, end) : row, dir == EAST ? _next(col, end) : col] + M12 = env[dir, dir == NORTH ? row - 1 : row, dir == EAST ? col + 1 : col] sqrtM12, isqrtM12 = sqrt_invsqrt(twist(M12, 1)) M21 = env[dir + 2, row, col] sqrtM21, isqrtM21 = sqrt_invsqrt(M21) From 879bb8e1665dd0bc060c2183be5f825d352a0435 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Fri, 8 May 2026 12:37:06 +0800 Subject: [PATCH 10/12] Drop unitcell argument in _nn_bondenv, _get_cluster_trunc --- src/algorithms/ctmrg/sequential.jl | 2 +- src/algorithms/time_evolution/get_cluster.jl | 10 +++++----- src/algorithms/time_evolution/simpleupdate.jl | 11 +++++------ src/algorithms/time_evolution/simpleupdate3site.jl | 10 ++++------ 4 files changed, 15 insertions(+), 18 deletions(-) diff --git a/src/algorithms/ctmrg/sequential.jl b/src/algorithms/ctmrg/sequential.jl index 55ddacfef..775164c7f 100644 --- a/src/algorithms/ctmrg/sequential.jl +++ b/src/algorithms/ctmrg/sequential.jl @@ -93,7 +93,7 @@ function sequential_projectors( coordinate::NTuple{3, Int}, network, env::CTMRGEnv, alg::HalfInfiniteProjector ) _, r, c = coordinate - r′ = _prev(r, size(env, 2)) + r′ = r - 1 trunc = truncation_strategy(alg, edge(env, WEST, r′, c)) alg´ = _set_decomposition_truncation(alg, trunc) Q1 = TensorMap(EnlargedCorner(network, env, (SOUTHWEST, r, c))) diff --git a/src/algorithms/time_evolution/get_cluster.jl b/src/algorithms/time_evolution/get_cluster.jl index 92b1c1480..452c0cd3b 100644 --- a/src/algorithms/time_evolution/get_cluster.jl +++ b/src/algorithms/time_evolution/get_cluster.jl @@ -54,19 +54,19 @@ Given `site1`, `site2` connected by a nearest neighbor bond, return the bond index and whether it is reversed from the standard orientation (`site1` on the west/south of `site2`). """ -function _nn_bondrev(site1::CartesianIndex{2}, site2::CartesianIndex{2}, (Nrow, Ncol)::NTuple{2, Int}) +function _nn_bondrev(site1::CartesianIndex{2}, site2::CartesianIndex{2}) diff = site1 - site2 if diff == CartesianIndex(0, -1) - r, c = mod1(site1[1], Nrow), mod1(site1[2], Ncol) + r, c = site1[1], site1[2] return (1, r, c), false elseif diff == CartesianIndex(0, 1) - r, c = mod1(site2[1], Nrow), mod1(site2[2], Ncol) + r, c = site2[1], site2[2] return (1, r, c), true elseif diff == CartesianIndex(1, 0) - r, c = mod1(site1[1], Nrow), mod1(site1[2], Ncol) + r, c = site1[1], site1[2] return (2, r, c), false elseif diff == CartesianIndex(-1, 0) - r, c = mod1(site2[1], Nrow), mod1(site2[2], Ncol) + r, c = site2[1], site2[2] return (2, r, c), true else error("`site1` and `site2` are not nearest neighbors.") diff --git a/src/algorithms/time_evolution/simpleupdate.jl b/src/algorithms/time_evolution/simpleupdate.jl index c4ab7f6d3..1ae452846 100644 --- a/src/algorithms/time_evolution/simpleupdate.jl +++ b/src/algorithms/time_evolution/simpleupdate.jl @@ -81,13 +81,12 @@ function _su_iter!( state::InfiniteState, gate::NNGate, env::SUWeight, sites::Vector{CartesianIndex{2}}, alg::SimpleUpdate ) - Nr, Nc = size(state) @assert length(sites) == 2 - trunc = only(_get_cluster_trunc(alg.trunc, sites, (Nr, Nc))) + trunc = only(_get_cluster_trunc(alg.trunc, sites)) Ms, open_vaxs, = _get_cluster(state, sites) _absorb_weight!(Ms, sites, open_vaxs, env) # rotate - bond, rev = _nn_bondrev(sites..., (Nr, Nc)) + bond, rev = _nn_bondrev(sites...) dir = first(bond) A, B = _bond_rotation.(Ms, dir, rev; inv = false) # apply gate @@ -124,7 +123,7 @@ function su_iter( state::InfiniteState, circuit::LocalCircuit, alg::SimpleUpdate, env::SUWeight ) - Nr, Nc, = size(state) + Nr = size(state, 2) state2, env2, ϵ = deepcopy(state), deepcopy(env), 0.0 for (sites, gate) in circuit.gates if length(sites) == 1 @@ -133,8 +132,8 @@ function su_iter( site = sites[1] state2[site] = _apply_sitegate(state2[site], gate; alg.purified) elseif length(sites) == 2 - (d, r, c), = _nn_bondrev(sites..., (Nr, Nc)) - alg.bipartite && r > 1 && continue + (d, r, c), = _nn_bondrev(sites...) + alg.bipartite && mod1(r, Nr) > 1 && continue ϵ′ = _su_iter!(state2, gate, env2, sites, alg) ϵ = max(ϵ, ϵ′) (!alg.bipartite) && continue diff --git a/src/algorithms/time_evolution/simpleupdate3site.jl b/src/algorithms/time_evolution/simpleupdate3site.jl index 28ca01339..80947b876 100644 --- a/src/algorithms/time_evolution/simpleupdate3site.jl +++ b/src/algorithms/time_evolution/simpleupdate3site.jl @@ -30,8 +30,7 @@ function _su_iter!( state::InfiniteState, gate::Vector{T}, env::SUWeight, sites::Vector{CartesianIndex{2}}, alg::SimpleUpdate ) where {T <: AbstractTensorMap} - Nr, Nc = size(state) - truncs = _get_cluster_trunc(alg.trunc, sites, (Nr, Nc)) + truncs = _get_cluster_trunc(alg.trunc, sites) Ms, open_vaxs, invperms = _get_cluster_with_weights(state, sites, env) # flip virtual arrows in `Ms` to ← flips = [isdual(space(M, 1)) for M in Iterators.drop(Ms, 1)] @@ -49,7 +48,7 @@ function _su_iter!( _flip_virtuals!(Ms, flips) # update env weights bond_revs = map(sites, Iterators.drop(sites, 1)) do site1, site2 - _nn_bondrev(site1, site2, (Nr, Nc)) + _nn_bondrev(site1, site2) end for (wt, (bond, rev), flip) in zip(wts, bond_revs, flips) wt_new = flip ? _fliptwist_s(wt) : wt @@ -70,11 +69,10 @@ Get the `TruncationStrategy` for each bond in the cluster updated by the Trotter evolution MPO. """ function _get_cluster_trunc( - trunc::TruncationStrategy, sites::Vector{CartesianIndex{2}}, - unitcell::NTuple{2, Int} + trunc::TruncationStrategy, sites::Vector{CartesianIndex{2}} ) return map(sites, Iterators.drop(sites, 1)) do site1, site2 - (d, r, c), rev = _nn_bondrev(site1, site2, unitcell) + (d, r, c), rev = _nn_bondrev(site1, site2) t = truncation_strategy(trunc, d, r, c) if rev && isa(t, TruncationSpace) t = truncspace(flip(t.space)') From 8cc09e39f6d7f9411ff7c903196bafa92a03bdca Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Fri, 8 May 2026 13:25:33 +0800 Subject: [PATCH 11/12] Periodic `physicalspace` for LocalOperator, LocalCircuit --- src/algorithms/time_evolution/trotter_gate.jl | 5 +---- src/operators/localcircuit.jl | 4 +++- src/operators/localoperator.jl | 3 ++- 3 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/algorithms/time_evolution/trotter_gate.jl b/src/algorithms/time_evolution/trotter_gate.jl index 67ae916d3..3e1fa7432 100644 --- a/src/algorithms/time_evolution/trotter_gate.jl +++ b/src/algorithms/time_evolution/trotter_gate.jl @@ -118,10 +118,7 @@ function _trotterize_nnn2site!(gates::Vector, H::LocalOperator, dt::Number) coord = [x1, x3] haskey(H.terms, coord) || continue gate = gate_to_mpo(exp(H.terms[coord] * -dt / 2)) - x2′ = CartesianIndex(mod1.(Tuple(x2), size(H))) - b = TensorKit.BraidingTensor{T}( - physicalspace(H, x2′), left_virtualspace(gate[2]) - ) + b = TensorKit.BraidingTensor{T}(physicalspace(H, x2), left_virtualspace(gate[2])) insert!(gate, 2, TensorMap(b)) push!(gates, [x1, x2, x3] => gate) end diff --git a/src/operators/localcircuit.jl b/src/operators/localcircuit.jl index f60a17ad2..f609aee34 100644 --- a/src/operators/localcircuit.jl +++ b/src/operators/localcircuit.jl @@ -105,7 +105,9 @@ end Return lattice of physical spaces on which the `LocalCircuit` is defined. """ physicalspace(gates::LocalCircuit) = gates.lattice -physicalspace(gates::LocalCircuit, args...) = physicalspace(gates)[args...] +Base.@propagate_inbounds physicalspace(gate::LocalCircuit, I...) = + periodic_getindex(gate, gate.lattice, I) + Base.size(gates::LocalCircuit) = size(physicalspace(gates)) # Equality diff --git a/src/operators/localoperator.jl b/src/operators/localoperator.jl index 251d0c5a5..a2d9baf95 100644 --- a/src/operators/localoperator.jl +++ b/src/operators/localoperator.jl @@ -147,7 +147,8 @@ end Return lattice of physical spaces on which the `LocalOperator` is defined. """ physicalspace(O::LocalOperator) = O.lattice -physicalspace(O::LocalOperator, args...) = physicalspace(O)[args...] +Base.@propagate_inbounds physicalspace(O::LocalOperator, I...) = + periodic_getindex(O, O.lattice, I) Base.size(O::LocalOperator, args...) = size(physicalspace(O), args...) Base.eltype(::Type{LocalOperator{O, S}}) where {O, S} = O From 0ad58a757a6388be4c93057cd490239a5a20a942 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Fri, 8 May 2026 13:36:24 +0800 Subject: [PATCH 12/12] Loose ends --- docs/src/examples/heisenberg_su/index.md | 4 ++-- docs/src/examples/heisenberg_su/main.ipynb | 4 ++-- examples/heisenberg_su/main.jl | 4 ++-- src/algorithms/bp/gaugefix.jl | 10 +++++----- src/algorithms/time_evolution/simpleupdate.jl | 3 +-- test/bondenv/benv_ctm.jl | 4 ++-- test/bp/gaugefix.jl | 2 +- test/ctmrg/suweight.jl | 2 +- test/ctmrg/unitcell.jl | 2 +- test/timeevol/cluster_projectors.jl | 4 ++-- 10 files changed, 19 insertions(+), 20 deletions(-) diff --git a/docs/src/examples/heisenberg_su/index.md b/docs/src/examples/heisenberg_su/index.md index aeba4f389..1225f06c4 100644 --- a/docs/src/examples/heisenberg_su/index.md +++ b/docs/src/examples/heisenberg_su/index.md @@ -67,8 +67,8 @@ else end peps = InfinitePEPS(rand, Float64, physical_space, bond_space; unitcell = (Nr, Nc)); -peps.A[2, 2] = copy(peps.A[1, 1]) ## make initial random state bipartite -peps.A[2, 1] = copy(peps.A[1, 2]) +peps[2, 2] = copy(peps[1, 1]) ## make initial random state bipartite +peps[2, 1] = copy(peps[1, 2]) wts = SUWeight(peps); ```` diff --git a/docs/src/examples/heisenberg_su/main.ipynb b/docs/src/examples/heisenberg_su/main.ipynb index 9296c0a9f..f762d6db9 100644 --- a/docs/src/examples/heisenberg_su/main.ipynb +++ b/docs/src/examples/heisenberg_su/main.ipynb @@ -97,8 +97,8 @@ "end\n", "\n", "peps = InfinitePEPS(rand, Float64, physical_space, bond_space; unitcell = (Nr, Nc));\n", - "peps.A[2, 2] = copy(peps.A[1, 1]) ## make initial random state bipartite\n", - "peps.A[2, 1] = copy(peps.A[1, 2])\n", + "peps[2, 2] = copy(peps[1, 1]) ## make initial random state bipartite\n", + "peps[2, 1] = copy(peps[1, 2])\n", "wts = SUWeight(peps);" ], "metadata": {}, diff --git a/examples/heisenberg_su/main.jl b/examples/heisenberg_su/main.jl index 05891a220..bdc6f2e08 100644 --- a/examples/heisenberg_su/main.jl +++ b/examples/heisenberg_su/main.jl @@ -60,8 +60,8 @@ else end peps = InfinitePEPS(rand, Float64, physical_space, bond_space; unitcell = (Nr, Nc)); -peps.A[2, 2] = copy(peps.A[1, 1]) ## make initial random state bipartite -peps.A[2, 1] = copy(peps.A[1, 2]) +peps[2, 2] = copy(peps[1, 1]) ## make initial random state bipartite +peps[2, 1] = copy(peps[1, 2]) wts = SUWeight(peps); md""" diff --git a/src/algorithms/bp/gaugefix.jl b/src/algorithms/bp/gaugefix.jl index 74b19409f..0f074dada 100644 --- a/src/algorithms/bp/gaugefix.jl +++ b/src/algorithms/bp/gaugefix.jl @@ -25,7 +25,7 @@ function gauge_fix(psi::InfinitePEPS, alg::BPGauge, env::BPEnv) # copy 1st column to 2nd column to eliminate differences # caused by order of applying gauge transformations for r in 1:2 - psi′[_next(r, 2), 2] = copy(psi′[r, 1]) + psi′[r + 1, 2] = copy(psi′[r, 1]) end end return psi′, XXinv @@ -90,10 +90,10 @@ function _bp_gauge_fix!(I::CartesianIndex{3}, psi::InfinitePEPS, env::BPEnv) end if dir == NORTH psi[row, col] = absorb_north_message(psi[row, col], X) - psi[_prev(row, end), col] = absorb_south_message(psi[_prev(row, end), col], invX) + psi[row - 1, col] = absorb_south_message(psi[row - 1, col], invX) elseif dir == EAST psi[row, col] = absorb_east_message(psi[row, col], X) - psi[row, _next(col, end)] = absorb_west_message(psi[row, _next(col, end)], invX) + psi[row, col + 1] = absorb_west_message(psi[row, col + 1], invX) end return psi, X, invX end @@ -122,9 +122,9 @@ to a belief propagation environment. function BPEnv(wts::SUWeight) messages = map(Iterators.product(1:4, axes(wts, 2), axes(wts, 3))) do (d, r, c) wt = if d == NORTH - twist(wts[2, _next(r, end), c], 1) + twist(wts[2, r + 1, c], 1) elseif d == EAST - twist(wts[1, r, _prev(c, end)], 1) + twist(wts[1, r, c - 1], 1) elseif d == SOUTH copy(wts[2, r, c]) else # WEST diff --git a/src/algorithms/time_evolution/simpleupdate.jl b/src/algorithms/time_evolution/simpleupdate.jl index 1ae452846..1fe8b543c 100644 --- a/src/algorithms/time_evolution/simpleupdate.jl +++ b/src/algorithms/time_evolution/simpleupdate.jl @@ -123,7 +123,6 @@ function su_iter( state::InfiniteState, circuit::LocalCircuit, alg::SimpleUpdate, env::SUWeight ) - Nr = size(state, 2) state2, env2, ϵ = deepcopy(state), deepcopy(env), 0.0 for (sites, gate) in circuit.gates if length(sites) == 1 @@ -133,7 +132,7 @@ function su_iter( state2[site] = _apply_sitegate(state2[site], gate; alg.purified) elseif length(sites) == 2 (d, r, c), = _nn_bondrev(sites...) - alg.bipartite && mod1(r, Nr) > 1 && continue + alg.bipartite && iseven(r) && continue ϵ′ = _su_iter!(state2, gate, env2, sites, alg) ϵ = max(ϵ, ϵ′) (!alg.bipartite) && continue diff --git a/test/bondenv/benv_ctm.jl b/test/bondenv/benv_ctm.jl index fcf6ef258..1b23b6f9b 100644 --- a/test/bondenv/benv_ctm.jl +++ b/test/bondenv/benv_ctm.jl @@ -38,8 +38,8 @@ function test_benv_ctm(state::Union{InfinitePEPS, InfinitePEPO}) network = isa(state, InfinitePEPS) ? state : InfinitePEPS(state) env, = leading_boundary(CTMRGEnv(rand, ComplexF64, network, Envspace), network, ctm_alg) for row in 1:Nr, col in 1:Nc - cp1 = PEPSKit._next(col, Nc) - A, B = state.A[row, col], state.A[row, cp1] + cp1 = col + 1 + A, B = state[row, col], state[row, cp1] a, X = PEPSKit.bond_tensor_first(A) b, Y = PEPSKit.bond_tensor_last(B) benv = PEPSKit.bondenv_ctm(row, col, X, Y, env) diff --git a/test/bp/gaugefix.jl b/test/bp/gaugefix.jl index 76d8fdc4e..bad3dc933 100644 --- a/test/bp/gaugefix.jl +++ b/test/bp/gaugefix.jl @@ -46,7 +46,7 @@ using PEPSKit: _next, _is_bipartite peps0 = InfinitePEPS(randn, elt, Pspaces, Nspaces, Espaces) if bipartite for c in 1:2 - peps0[2, c] = copy(peps0[1, _next(c, 2)]) + peps0[2, c] = copy(peps0[1, c + 1]) end end diff --git a/test/ctmrg/suweight.jl b/test/ctmrg/suweight.jl index 12ff64e09..570ea0ee6 100644 --- a/test/ctmrg/suweight.jl +++ b/test/ctmrg/suweight.jl @@ -2,7 +2,7 @@ using Test using Random using TensorKit using PEPSKit -using PEPSKit: str, twistdual, _prev, _next, unitcell +using PEPSKit: str, twistdual, unitcell Vps = Dict( Z2Irrep => Vect[Z2Irrep](0 => 1, 1 => 2), diff --git a/test/ctmrg/unitcell.jl b/test/ctmrg/unitcell.jl index 1b83de2ef..e57c14ec8 100644 --- a/test/ctmrg/unitcell.jl +++ b/test/ctmrg/unitcell.jl @@ -1,7 +1,7 @@ using Test using Random using PEPSKit -using PEPSKit: _prev, _next, ctmrg_iteration, compute_gauge_fix_gauge, ScramblingEnvGauge +using PEPSKit: ctmrg_iteration, compute_gauge_fix_gauge, ScramblingEnvGauge using TensorKit # settings diff --git a/test/timeevol/cluster_projectors.jl b/test/timeevol/cluster_projectors.jl index 6b86f6ee6..6e48ec6b4 100644 --- a/test/timeevol/cluster_projectors.jl +++ b/test/timeevol/cluster_projectors.jl @@ -4,7 +4,7 @@ using PEPSKit using LinearAlgebra using Random import MPSKitModels: hubbard_space -using PEPSKit: sdiag_pow, _cluster_truncate!, _flip_virtuals!, _next +using PEPSKit: sdiag_pow, _cluster_truncate!, _flip_virtuals! using MPSKit: GenericMPSTensor, MPSBondTensor # Utility setup @@ -230,7 +230,7 @@ end peps0 = InfinitePEPS(rand, Float64, Pspace, Vspace, Vspace'; unitcell = (Nr, Nc)) # make initial state bipartite for r in 1:2 - peps0.A[_next(r, 2), 2] = copy(peps0.A[r, 1]) + peps0[r + 1, 2] = copy(peps0[r, 1]) end wts0 = SUWeight(peps0) ham = hubbard_model(Float64, Trivial, U1Irrep, InfiniteSquare(Nr, Nc); t = 1.0, U = 6.0, mu = 3.0)