diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl index 827a16689..c6436fb09 100644 --- a/src/PEPSKit.jl +++ b/src/PEPSKit.jl @@ -28,12 +28,15 @@ using LoggingExtras import TupleTools using MPSKit -using MPSKit: MPSTensor, MPOTensor, GenericMPSTensor, MPSBondTensor, ProductTransferMatrix +using MPSKit: + MPSTensor, MPOTensor, GenericMPSTensor, MPSBondTensor, + ProductTransferMatrix, TransferMatrix using MPSKit: InfiniteEnvironments using MPSKit: DynamicTol, updatetol import MPSKit.DynamicTols: _updatetol import MPSKit: tensorexpr, leading_boundary, loginit!, logiter!, logfinish!, logcancel!, physicalspace import MPSKit: infinite_temperature_density_matrix +import MPSKit: fuser using TensorKitTensors: fuse_charge import TensorKitTensors.SpinOperators as SO @@ -79,6 +82,7 @@ include("operators/infinitepepo.jl") include("operators/transfermatrix.jl") include("operators/localoperator.jl") include("operators/localcircuit.jl") + include("operators/lattices/squarelattice.jl") include("operators/models.jl") @@ -120,6 +124,12 @@ include("algorithms/contractions/correlator/peps.jl") include("algorithms/contractions/correlator/pepo_purified.jl") include("algorithms/contractions/correlator/pepo_1layer.jl") +include("algorithms/contractions/mpo_path/routing.jl") +include("algorithms/contractions/mpo_path/pepo_1layer.jl") +include("algorithms/contractions/patch/tools.jl") +include("algorithms/contractions/patch/pepo_1layer.jl") +include("algorithms/contractions/patch/twosite/pepo_1layer.jl") + include("algorithms/ctmrg/sparse_environments.jl") include("algorithms/ctmrg/ctmrg.jl") include("algorithms/ctmrg/projectors/projectors.jl") @@ -160,6 +170,8 @@ include("algorithms/expectation_value/network_value.jl") include("algorithms/expectation_value/correlator_adapters.jl") include("algorithms/expectation_value/correlators.jl") include("algorithms/toolbox.jl") +include("algorithms/expval_approx.jl") +include("algorithms/correlator_approx.jl") include("algorithms/optimization/implicit_differentiation.jl") include("algorithms/optimization/preconditioning.jl") @@ -178,9 +190,10 @@ export FixedSpaceTruncation, SiteDependentTruncation export HalfInfiniteProjector, FullInfiniteProjector export C4vCTMRG, C4vEighProjector, C4vQRProjector export initialize_random_c4v_env, initialize_singlet_c4v_env -export LocalOperator, physicalspace +export LocalOperator, MPOTerm, physicalspace export product_peps -export reduced_densitymatrix, expectation_value, network_value, cost_function +export reduced_densitymatrix, expectation_value_approx, correlator_approx +export expectation_value, network_value, cost_function export correlator, correlation_length export leading_boundary export PEPSOptimize, FixedPointGradient, GeomSum, ManualIter, ImplicitGradient diff --git a/src/algorithms/contractions/mpo_path/pepo_1layer.jl b/src/algorithms/contractions/mpo_path/pepo_1layer.jl new file mode 100644 index 000000000..e93e8dbad --- /dev/null +++ b/src/algorithms/contractions/mpo_path/pepo_1layer.jl @@ -0,0 +1,247 @@ +""" +Check that the physical legs of a first OBC-MPO tensor match the PEPO site's physical space. +""" +function _check_pepo_first_physicalspace(A, op) + physicalspace(A) == space(op, 1) == space(op, 2)' || + throw(SpaceMismatch("first MPO tensor physical space does not match PEPO site")) + return nothing +end + +""" +Check that the physical legs of a last OBC-MPO tensor match the PEPO site's physical space. +""" +function _check_pepo_last_physicalspace(A, op) + physicalspace(A) == space(op, 2) == space(op, 3)' || + throw(SpaceMismatch("last MPO tensor physical space does not match PEPO site")) + return nothing +end + +""" +Check that the physical legs of a middle MPO tensor match the PEPO site's physical space. +""" +function _check_pepo_middle_physicalspace(A, op) + physicalspace(A) == space(op, 2) == space(op, 3)' || + throw(SpaceMismatch("middle MPO tensor physical space does not match PEPO site")) + return nothing +end + +""" +Convert a symbolic cardinal path direction to the corresponding PEPO virtual-leg index. +""" +function _mpo_path_direction(direction::Symbol) + direction === :north && return NORTH + direction === :east && return EAST + direction === :south && return SOUTH + direction === :west && return WEST + throw(ArgumentError("invalid MPO path direction: $direction")) +end + +""" +Return the tensor-expression label for the PEPO virtual leg in a cardinal direction. +""" +function _mpo_path_virtual_label(direction::Symbol) + direction === :north && return :N + direction === :east && return :E + direction === :south && return :S + direction === :west && return :W + throw(ArgumentError("invalid MPO path direction: $direction")) +end + +""" +Canonicalize an incoming MPO fuser so its fused PEPO leg has standard dualness. +""" +function _mpo_path_incoming_fuser(F, direction::Int) + direction in (NORTH, EAST) && return F + direction in (SOUTH, WEST) && return twist(flip(F, 1), 1) + throw(ArgumentError("invalid MPO path direction index: $direction")) +end + +""" +Canonicalize an outgoing MPO fuser so its fused PEPO leg has standard dualness and braiding. +""" +function _mpo_path_outgoing_fuser(F, direction::Int) + direction in (NORTH, EAST) && return twist(flip(F, 1), 1) + direction in (SOUTH, WEST) && return twist(F, 3) + throw(ArgumentError("invalid MPO path direction index: $direction")) +end + +""" +Build the `@tensor` labels from `(direction, suffix)` pairs +used to fuse MPO virtual strings. + +- Each direction is `:north`, `:east`, `:south`, or `:west`. +- Suffix `:l` marks an incoming MPO bond and `:r` an outgoing one. + +Examples: + +- `((:east, :r),)` produces `[W S; N Er]`. +- `((:west, :l), (:north, :r))` produces `[Wl S; Nr E]`. +""" +function _mpo_path_result_expr(directions) + labels = [:N, :E, :S, :W] + for (direction, suffix) in directions + index = _mpo_path_direction(direction) + labels[index] = Symbol(labels[index], suffix) + end + return tensorexpr(:t, (labels[WEST], labels[SOUTH]), (labels[NORTH], labels[EAST])) +end + +""" +Act the first tensor `op` of an OBC-MPO on PEPO tensor `A` and fuse the outgoing MPO string with the virtual space of `A` along `direction`. +""" +function mpo_path_first(A::PEPOTensor, op, direction::Val{D}) where {D} + _check_pepo_first_physicalspace(A, op) + direction_index = _mpo_path_direction(D) + A′ = twistdual(A, 2) + F = _mpo_path_outgoing_fuser( + fuser(storagetype(A), domain(A, direction_index)', space(op, 3)), + direction_index, + ) + return _mpo_path_first(A′, op, F, direction) +end + +""" +Contract a prepared PEPO tensor with the first MPO tensor and its outgoing fuser. +""" +@generated function _mpo_path_first(A, op, F, ::Val{direction}) where {direction} + virtual_label = _mpo_path_virtual_label(direction) + fused_label = Symbol(virtual_label, :r) + + result_e = _mpo_path_result_expr(((direction, :r),)) + op_e = tensorexpr(:op, :dout, (:din, :r)) + A_e = tensorexpr(:A, (:din, :dout), (:N, :E, :S, :W)) + F_e = tensorexpr(:F, fused_label, (virtual_label, :r)) + rhs = Expr(:call, :*, op_e, A_e, F_e) + return macroexpand( + @__MODULE__, :(return @tensoropt $result_e := $rhs) + ) +end + +""" +Act the last tensor `op` of an OBC-MPO on PEPO tensor `A` and fuse the incoming MPO string with the virtual space of `A` along `direction`. +""" +function mpo_path_last(A::PEPOTensor, op, direction::Val{D}) where {D} + _check_pepo_last_physicalspace(A, op) + direction_index = _mpo_path_direction(D) + A′ = twistdual(A, 2) + F = _mpo_path_incoming_fuser( + fuser(storagetype(A), domain(A, direction_index), space(op, 1)'), + direction_index, + ) + return _mpo_path_last(A′, op, F, direction) +end + +""" +Contract a prepared PEPO tensor with the last MPO tensor and its incoming fuser. +""" +@generated function _mpo_path_last(A, op, F, ::Val{direction}) where {direction} + virtual_label = _mpo_path_virtual_label(direction) + fused_label = Symbol(virtual_label, :l) + + result_e = _mpo_path_result_expr(((direction, :l),)) + F_e = Expr(:call, :conj, tensorexpr(:F, fused_label, (virtual_label, :l))) + op_e = tensorexpr(:op, (:l, :dout), :din) + A_e = tensorexpr(:A, (:din, :dout), (:N, :E, :S, :W)) + rhs = Expr(:call, :*, F_e, op_e, A_e) + return macroexpand( + @__MODULE__, :(return @tensoropt $result_e := $rhs) + ) +end + +""" +Act the middle tensor `op` of an MPO on PEPO tensor `A` and fuse the incoming and outgoing MPO strings with the virtual space of `A` along `directions = (incoming, outgoing)`. +""" +function mpo_path_middle(A::PEPOTensor, op, directions::Val{D}) where {D} + incoming, outgoing = D + incoming == outgoing && + throw(ArgumentError("MPO path should enter and exit in different directions")) + _check_pepo_middle_physicalspace(A, op) + incoming_index = _mpo_path_direction(incoming) + outgoing_index = _mpo_path_direction(outgoing) + A′ = twistdual(A, 2) + Fin = _mpo_path_incoming_fuser( + fuser(storagetype(A), domain(A, incoming_index), space(op, 1)'), + incoming_index, + ) + Fout = _mpo_path_outgoing_fuser( + fuser(storagetype(A), domain(A, outgoing_index)', space(op, 4)), + outgoing_index, + ) + return _mpo_path_middle(A′, op, Fin, Fout, directions) +end + +""" +Contract a prepared PEPO tensor with a middle MPO tensor and its incoming and outgoing fusers. +""" +@generated function _mpo_path_middle( + A, op, Fin, Fout, ::Val{directions} + ) where {directions} + incoming, outgoing = directions + incoming_label = _mpo_path_virtual_label(incoming) + outgoing_label = _mpo_path_virtual_label(outgoing) + fused_incoming_label = Symbol(incoming_label, :l) + fused_outgoing_label = Symbol(outgoing_label, :r) + + result_e = _mpo_path_result_expr(((incoming, :l), (outgoing, :r))) + Fin_e = Expr( + :call, :conj, + tensorexpr(:Fin, fused_incoming_label, (incoming_label, :l)), + ) + op_e = tensorexpr(:op, (:l, :dout), (:din, :r)) + A_e = tensorexpr(:A, (:din, :dout), (:N, :E, :S, :W)) + Fout_e = tensorexpr( + :Fout, fused_outgoing_label, (outgoing_label, :r) + ) + rhs = Expr(:call, :*, Fin_e, op_e, A_e, Fout_e) + return macroexpand( + @__MODULE__, :(return @tensoropt $result_e := $rhs) + ) +end + +""" +Route an MPO virtual string with `stringspace` through a PEPO tensor `A` +along `directions = (incoming, outgoing)`. +""" +@generated function mpo_path_string( + A::PEPOTensor, stringspace::ElementarySpace, ::Val{directions} + ) where {directions} + incoming, outgoing = directions + incoming == outgoing && + throw(ArgumentError("MPO path should enter and exit in different directions")) + + incoming_index = _mpo_path_direction(incoming) + outgoing_index = _mpo_path_direction(outgoing) + incoming_label = _mpo_path_virtual_label(incoming) + outgoing_label = _mpo_path_virtual_label(outgoing) + fused_incoming_label = Symbol(incoming_label, :l) + fused_outgoing_label = Symbol(outgoing_label, :r) + + result_e = _mpo_path_result_expr(((incoming, :l), (outgoing, :r))) + Fin_e = Expr( + :call, :conj, + tensorexpr(:Fin, fused_incoming_label, (incoming_label, :l)), + ) + O_e = tensorexpr(:O, (:W, :S), (:N, :E)) + I_e = tensorexpr(:I, :l, :r) + Fout_e = tensorexpr( + :Fout, fused_outgoing_label, (outgoing_label, :r) + ) + rhs = Expr(:call, :*, Fin_e, O_e, I_e, Fout_e) + contraction = macroexpand( + @__MODULE__, :(return @tensoropt $result_e := $rhs) + ) + + return quote + O = trace_physicalspaces(A) + I = id(storagetype(A), stringspace) + Fin = _mpo_path_incoming_fuser( + fuser(storagetype(A), domain(A, $incoming_index), stringspace'), + $incoming_index, + ) + Fout = _mpo_path_outgoing_fuser( + fuser(storagetype(A), domain(A, $outgoing_index)', stringspace'), + $outgoing_index, + ) + $contraction + end +end diff --git a/src/algorithms/contractions/mpo_path/routing.jl b/src/algorithms/contractions/mpo_path/routing.jl new file mode 100644 index 000000000..81995118a --- /dev/null +++ b/src/algorithms/contractions/mpo_path/routing.jl @@ -0,0 +1,206 @@ +""" +A `path => mpo` pair with one MPO tensor per site of a non-self-intersecting nearest-neighbor lattice path. +Intermediate sites carry inserted braiding tensors that propagate the MPO bond. +These properties are checked by routing validation, not enforced by this alias. +""" +const RoutedMPOTerm = Pair{<:LatticePath, <:MPOTerm} + +""" +Convert a dense term to an exact MPO in column-snake order and expand its nearest-neighbor path. + +Visit occupied columns from left to right, alternating between increasing and decreasing row indices within each column. +Permute the operator's physical legs into this order before decomposing it into an MPO. +For example, the six operator sites below are visited as `A → B → C → D → E → F`, with row indices increasing downward: + +```text + col 1 col 2 col 3 +row 1 A + --→-- E + ↓ ↑ ↓ +row 2 + D + + ↓ ↑ ↓ +row 3 B + F + ↓ ↑ +row 4 + --→-- C +``` + +The path is constructed successively between these ordered sites: + +1. `A → B`: descend from `(1, 1)` through `(2, 1)` to `(3, 1)`. +2. `B → C`: connect the column bottoms at row `max(3, 4) = 4`, passing through `(4, 1)` before reaching `(4, 2)`. +3. `C → D`: ascend through `(3, 2)` to `(2, 2)`. +4. `D → E`: connect the column tops at row `min(2, 1) = 1`, passing through `(1, 2)` before reaching `(1, 3)`. +5. `E → F`: descend through `(2, 3)` to `(3, 3)`. + +The `+` sites carry inserted braiding tensors that propagate the MPO bond without an additional physical operator. +Connecting successive columns at their outermost endpoint row avoids retracing the path: a horizontal-first connection from `B` to `C` would visit `(3, 2)`, which the later `C → D` segment needs. +""" +function _route_mpo_term( + sites::LatticePath, op::AbstractTensorMap, + lattice::Matrix{<:ElementarySpace}, + ) + isempty(sites) && throw(ArgumentError("an operator term requires at least one site")) + allunique(sites) || throw(ArgumentError("operator sites should be unique")) + length(sites) == numout(op) == numin(op) || + throw(ArgumentError("number of operator legs should match the number of sites")) + Nr, Nc = size(lattice) + for (k, site) in enumerate(sites) + V = lattice[mod1(site[1], Nr), mod1(site[2], Nc)] + V == codomain(op)[k] == domain(op)[k] || + throw(SpaceMismatch("operator physical space does not match lattice site $site")) + end + length(sites) == 1 && return copy(sites) => [op] + + columns = sort!(unique(site[2] for site in sites)) + permutation = Int[] + for (k, col) in enumerate(columns) + indices = findall(site -> site[2] == col, sites) + sort!(indices; by = i -> sites[i][1], rev = iseven(k)) + append!(permutation, indices) + end + ordered_sites = sites[permutation] + N = length(sites) + ordered_op = permute(op, (Tuple(permutation), Tuple(permutation .+ N))) + mpo = gate_to_mpo(ordered_op; trunc = notrunc()) + path = _dense_mpo_path(ordered_sites) + return _expand_mpo_path(ordered_sites, mpo, path, lattice) +end + +""" +Preserve an explicit MPO's factor order while connecting its sites with simple nearest-neighbor paths. +""" +function _route_mpo_term( + sites::LatticePath, mpo::MPOTerm, + lattice::Matrix{<:ElementarySpace}, + ) + _validate_mpo_term(sites, mpo, lattice) + path = _ordered_mpo_path(sites) + return _expand_mpo_path(sites, mpo, path, lattice) +end + +""" +Connect snake-ordered columns beyond their bottom or top endpoints to avoid retracing sparse columns. +""" +function _dense_mpo_path(sites::LatticePath) + path = CartesianIndex{2}[first(sites)] + column_index = 1 + for (from, to) in zip(sites, Iterators.drop(sites, 1)) + if from[2] == to[2] + append!(path, Iterators.drop(_l_path(from, to), 1)) + else + row = isodd(column_index) ? max(from[1], to[1]) : min(from[1], to[1]) + corners = (CartesianIndex(row, from[2]), CartesianIndex(row, to[2]), to) + for corner in corners + last(path) == corner && continue + append!(path, Iterators.drop(_l_path(last(path), corner), 1)) + end + column_index += 1 + end + end + return path +end + +""" +Route ordered sites with horizontal-first L paths, falling back to vertical-first paths without backtracking. +""" +function _ordered_mpo_path(sites::LatticePath) + isempty(sites) && throw(ArgumentError("an MPO term requires at least one site")) + allunique(sites) || throw(ArgumentError("operator sites should be unique")) + path = CartesianIndex{2}[first(sites)] + measured = Set(sites) + visited = Set(path) + for (from, to) in zip(sites, Iterators.drop(sites, 1)) + segment = _l_path(from, to) + if !_mpo_segment_is_free(segment, measured, visited) + segment = _l_path(from, to; horizontal_first = false) + _mpo_segment_is_free(segment, measured, visited) || throw( + ArgumentError( + "routing this MPO ordering is not implemented: no nonintersecting L path from $from to $to", + ), + ) + end + append!(path, Iterators.drop(segment, 1)) + union!(visited, segment) + end + return path +end + +""" +Check that a connecting segment neither revisits a site nor crosses another operator site. +""" +function _mpo_segment_is_free(segment, measured::Set, visited::Set) + any(site -> site in visited, @view segment[2:end]) && return false + return all(site -> !(site in measured), @view segment[2:(end - 1)]) +end + +""" +Insert braiding tensors at unmeasured path sites using periodic physical spaces and the neighboring MPO bond. +""" +function _expand_mpo_path( + sites::LatticePath, mpo::MPOTerm, + path::LatticePath, lattice::Matrix{<:ElementarySpace}, + ) + _validate_mpo_term(sites, mpo, lattice) + allunique(path) || throw(ArgumentError("the MPO path should not intersect itself")) + positions = indexin(sites, path) + all(!isnothing, positions) && issorted(positions) || + throw(ArgumentError("the MPO path should contain operator sites in MPO order")) + first(positions) == 1 && last(positions) == length(path) || + throw(ArgumentError("the MPO path should start and end at operator sites")) + for (from, to) in zip(path, Iterators.drop(path, 1)) + _step_direction(from, to) + end + + Nr, Nc = size(lattice) + k = 1 + expanded_mpo = AbstractTensorMap[] + for site in path + if site == sites[k] + op = mpo[k] + k += 1 + push!(expanded_mpo, op) + continue + end + previous = mpo[k - 1] + V = lattice[mod1(site[1], Nr), mod1(site[2], Nc)] + bond = _mpo_right_stringspace(previous)' + braid = TensorKit.BraidingTensor{scalartype(previous)}(V, bond) + push!(expanded_mpo, copy!(similar(previous, scalartype(previous), space(braid)), braid)) + end + return copy(path) => expanded_mpo +end + +""" +Return the right MPO string space in the local tensor's stored leg orientation. +""" +_mpo_right_stringspace(op::AbstractTensorMap) = space(op, numind(op)) + +""" +Build a shortest nearest-neighbor L path with the chosen horizontal or vertical first step. +""" +function _l_path( + start::CartesianIndex{2}, stop::CartesianIndex{2}; horizontal_first::Bool = true + ) + start == stop && throw(ArgumentError("MPO path sites should be unique")) + path = CartesianIndex{2}[start] + axes = horizontal_first ? (2, 1) : (1, 2) + for axis in axes + while last(path)[axis] != stop[axis] + step = sign(stop[axis] - last(path)[axis]) + delta = axis == 1 ? CartesianIndex(step, 0) : CartesianIndex(0, step) + push!(path, last(path) + delta) + end + end + return path +end + +""" +Return the cardinal direction of a nearest-neighbor step, rejecting longer or stationary steps. +""" +function _step_direction(from::CartesianIndex{2}, to::CartesianIndex{2}) + delta = to - from + delta == CartesianIndex(0, 1) && return :east + delta == CartesianIndex(0, -1) && return :west + delta == CartesianIndex(1, 0) && return :south + delta == CartesianIndex(-1, 0) && return :north + throw(ArgumentError("MPO path should use nearest-neighbor steps")) +end diff --git a/src/algorithms/contractions/patch/pepo_1layer.jl b/src/algorithms/contractions/patch/pepo_1layer.jl new file mode 100644 index 000000000..a366a7775 --- /dev/null +++ b/src/algorithms/contractions/patch/pepo_1layer.jl @@ -0,0 +1,179 @@ +# Approximate finite-patch contractions for single-layer PEPO networks. + +""" +Validate that the network is a single-layer PEPO and that the sweep direction is supported. +""" +function _check_patch_inputs(ρ::InfinitePEPO, direction::Symbol) + size(ρ, 3) == 1 || throw(DimensionMismatch("only single-layer PEPO contractions are supported")) + direction in (:auto, :rows, :columns) || + throw(ArgumentError("invalid sweep direction: $direction")) + return nothing +end + +""" +Return a PEPO and CTMRG environment with standard virtual-space dualness without mutating the inputs. +""" +function standardize_dualness(ρ::InfinitePEPO, env::CTMRGEnv) + isdual_easts, isdual_norths = _check_virtual_dualness(ρ) + all(isdual_easts) && all(isdual_norths) && return ρ, env + + nrows, ncols = size(ρ, 1), size(ρ, 2) + tensors = map(CartesianIndices(unitcell(ρ))) do site + row, col, layer = Tuple(site) + directions = Int[] + !isdual_norths[row, col, layer] && push!(directions, NORTH) + !isdual_easts[row, col, layer] && push!(directions, EAST) + !isdual_norths[_next(row, nrows), col, layer] && push!(directions, SOUTH) + !isdual_easts[row, _prev(col, ncols), layer] && push!(directions, WEST) + A = unitcell(ρ)[site] + return isempty(directions) ? A : flip_virtualspace(A, directions) + end + ρ′ = InfinitePEPO(tensors) + + edges = map(CartesianIndices(env.edges)) do index + direction, row, col = Tuple(index) + should_flip = if direction == NORTH + !isdual_norths[_next(row, nrows), col, 1] + elseif direction == EAST + !isdual_easts[row, _prev(col, ncols), 1] + elseif direction == SOUTH + !isdual_norths[row, col, 1] + else + !isdual_easts[row, col, 1] + end + E = env.edges[index] + return should_flip ? flip(E, 2) : E + end + env′ = CTMRGEnv(copy(env.corners), edges) + return ρ′, env′ +end + +""" +Contract a routed MPO term in its enclosing patch, rotating column sweeps into row sweeps. +""" +function _expectation_value_approx( + ρ::InfinitePEPO, routed::RoutedMPOTerm, env::CTMRGEnv, + alg::PatchApprox, direction::Symbol, + ) + _check_patch_inputs(ρ, direction) + rowrange, colrange = _patch_ranges(first(routed)) + sweep = if direction === :auto + length(colrange) >= length(rowrange) ? :rows : :columns + else + direction + end + if sweep === :rows + return _expectation_value_approx_rows( + ρ, routed, env, rowrange, colrange, alg + ) + else + unitcell = size(ρ)[1:2] + path, mpo = routed + rotated = siterotl90.(path, Ref(unitcell)) => mpo + rotated_rowrange, rotated_colrange = _patch_ranges(first(rotated)) + return _expectation_value_approx_rows( + rotl90(ρ), rotated, rotl90(env), + rotated_rowrange, rotated_colrange, alg + ) + end +end + +""" +Build a local tensor for row MPOs without observables by tracing the PEPO physical legs. +""" +function _patch_site_tensor( + ρ::InfinitePEPO, ::Nothing, row::Int, col::Int, + ) + return trace_physicalspaces(ρ[row, col, 1]) +end + +""" +Insert the MPO factor at a path site, or trace the PEPO physical legs at an off-path site. +""" +function _patch_site_tensor( + ρ::InfinitePEPO, routed::RoutedMPOTerm, + row::Int, col::Int, + ) + A = ρ[row, col, 1] + path, mpo = routed + k = findfirst(==(CartesianIndex(row, col)), path) + isnothing(k) && return trace_physicalspaces(A) + + if k == 1 + direction = _step_direction(path[1], path[2]) + return mpo_path_first(A, mpo[k], Val(direction)) + elseif k == length(path) + direction = _step_direction(path[end], path[end - 1]) + return mpo_path_last(A, mpo[k], Val(direction)) + else + incoming = _step_direction(path[k], path[k - 1]) + outgoing = _step_direction(path[k], path[k + 1]) + return mpo_path_middle(A, mpo[k], Val((incoming, outgoing))) + end +end + +""" +Contract and normalize a routed MPO term using row-oriented patch boundary contractions. +""" +function _expectation_value_approx_rows( + ρ::InfinitePEPO, routed::RoutedMPOTerm, env::CTMRGEnv, + rowrange::UnitRange{Int}, colrange::UnitRange{Int}, alg::PatchApprox, + ) + ρ, env = standardize_dualness(ρ, env) + numerator = _contract_patch_rows(ρ, routed, env, rowrange, colrange, alg) + norm = _contract_patch_rows(ρ, nothing, env, rowrange, colrange, alg) + return numerator / norm +end + +""" +Contract a complete PEPO patch row by row from north to south, optionally inserting a routed MPO term. +""" +function _contract_patch_rows( + ρ::InfinitePEPO, routed::Union{Nothing, RoutedMPOTerm}, + env::CTMRGEnv, rowrange::UnitRange{Int}, colrange::UnitRange{Int}, + alg::PatchApprox, + ) + ψ = _north_boundary_mps(env, first(rowrange), colrange) + for row in rowrange + W = _row_mpo(ρ, routed, env, row, colrange) + ψ = _approximate(W, ψ, alg) + end + south = _south_boundary_mps(env, last(rowrange), colrange) + return dot_noconj(south, ψ) +end + +""" +Build one finite row MPO from west/east CTMRG edges and the PEPO tensors inside the patch. + +Convention of west, east CTM edges and the PF tensors: +``` + [1 2; 3] [1 2; 3 4] [1 2; 3] + 3 3 1 + ↓ ↓ ↑ + E₄-←-2 1-←-O-←-4 2-←-C₂ + ↓ ↓ ↑ + 1 2 3 +``` +Legs 1, 3 need to be flipped to match standard MPS convention +""" +function _row_mpo( + ρ::InfinitePEPO, routed::Union{Nothing, RoutedMPOTerm}, + env::CTMRGEnv, row::Int, colrange::UnitRange{Int}, + ) + cmin, cmax = first(colrange), last(colrange) + W = repartition(edge(env, WEST, row, cmin - 1), 1, 2) + tensors = [insertleftunit(W, 1)] + append!( + tensors, + ( + _patch_site_tensor(ρ, routed, row, col) + for col in colrange + ), + ) + E = permute( + flip(edge(env, EAST, row, cmax + 1), (1, 3)), + ((2, 3), (1,)) + ) + push!(tensors, insertrightunit(E, 3)) + return FiniteMPO(tensors) +end diff --git a/src/algorithms/contractions/patch/tools.jl b/src/algorithms/contractions/patch/tools.jl new file mode 100644 index 000000000..1a9b0e586 --- /dev/null +++ b/src/algorithms/contractions/patch/tools.jl @@ -0,0 +1,80 @@ +""" +Bundle the zip-up contraction and optional DMRG refinement +algorithms used after each finite patch MPO-MPS contraction. +""" +struct PatchApprox{Z, D} + zipup::Z + dmrg::D +end + +# TODO: generalize the following to multi-layer networks + +""" +Apply a finite MPO to a finite MPS with zip-up truncation and optional DMRG refinement. +""" +function _approximate(W::FiniteMPO, ψ::FiniteMPS, alg::PatchApprox) + ψ′, = approximate((W, ψ), alg.zipup) + isnothing(alg.dmrg) && return ψ′ + ψ′, = approximate(ψ′, (W, ψ), alg.dmrg) + return ψ′ +end + +""" +Build the finite MPS representing the north CTMRG boundary of a patch. + +Convention of CTM tensors on the north boundary is +``` + [1; 2] [1 2; 3] [1; 2] + C₁-←-2 1-←-E₁-←-3 1-←-C₂ + ↓ ↓ ↑ + 1 2 2 +``` +Leg 2 of C₂ needs to be flipped to have a non-dual physical space. +""" +function _north_boundary_mps( + env::CTMRGEnv, row::Int, colrange::UnitRange{Int}, + ) + r = row - 1 + cmin, cmax = first(colrange), last(colrange) + Cwest = insertleftunit(corner(env, NORTHWEST, r, cmin - 1), 1) + tensors = [Cwest] + append!(tensors, (edge(env, NORTH, r, col) for col in colrange)) + Ceast = repartition( + flip(corner(env, NORTHEAST, r, cmax + 1), 2), 2, 0 + ) + push!(tensors, insertleftunit(Ceast, 3)) + return FiniteMPS(tensors) +end + +""" +Build the finite MPS representing the south CTMRG boundary of a patch, +but with dual physical legs, and sites ordered from east to west. + +Convention of CTM tensors on the south boundary is +``` + [1; 2] [1 2; 3] [1; 2] + 2 2 1 + ↓ ↓ ↑ + C₄-→-1 3-→-E₃-→-1 2-→-C₃ +``` +Leg 1 of C₃ needs to be flipped to have a dual physical space. + +North site `k` pairs with south site `N + 1 - k`. +``` + west east + north: 1 ← … ← N - 1 ← N + south: N → … → 2 → 1 +``` +Viewed after a 180° rotation of the entire network, this is a north boundary ordered from west to east as usual, only with dual physical legs. +Because of this reversed site order, `AL` tensors lie to the right (east) of the canonical center in the patch, while `AR` tensors lie to its left (west). +""" +function _south_boundary_mps(env::CTMRGEnv, row::Int, colrange::UnitRange{Int}) + r = row + 1 + cmin, cmax = first(colrange), last(colrange) + Ceast = insertleftunit(flip(corner(env, SOUTHEAST, r, cmax + 1), 1), 1) + tensors = [Ceast] + append!(tensors, (edge(env, SOUTH, r, col) for col in reverse(colrange))) + Cwest = repartition(corner(env, SOUTHWEST, r, cmin - 1), 2, 0) + push!(tensors, insertleftunit(Cwest, 3)) + return FiniteMPS(tensors) +end diff --git a/src/algorithms/contractions/patch/twosite/pepo_1layer.jl b/src/algorithms/contractions/patch/twosite/pepo_1layer.jl new file mode 100644 index 000000000..b747045bd --- /dev/null +++ b/src/algorithms/contractions/patch/twosite/pepo_1layer.jl @@ -0,0 +1,173 @@ +""" +Validate operator spaces and rotate column sweeps into the row-oriented contraction. +""" +function _correlator_approx( + ρ::InfinitePEPO, op::AbstractTensorMap, + source::CartesianIndex{2}, targets::Vector{CartesianIndex{2}}, + env::CTMRGEnv, alg::PatchApprox, direction::Symbol, + ) + _check_patch_inputs(ρ, direction) + numout(op) == numin(op) == 2 || + throw(ArgumentError("correlator_approx requires a two-site operator")) + for (leg, sites) in enumerate(((source,), targets)), site in sites + V = physicalspace(ρ, Tuple(site)...) + V == codomain(op)[leg] == domain(op)[leg] || + throw(SpaceMismatch("operator physical space does not match PEPO site $site")) + end + if direction === :columns + unitcell = size(ρ)[1:2] + source = siterotl90(source, unitcell) + targets = siterotl90.(targets, Ref(unitcell)) + ρ, env = rotl90(ρ), rotl90(env) + end + return _correlator_approx_rows(ρ, op, source, targets, env, alg) +end + +""" +Measure targets in patches of fixed width ending at each target row, propagating observable-free and open-string north states together. +""" +function _correlator_approx_rows( + ρ::InfinitePEPO, op::AbstractTensorMap, + source::CartesianIndex{2}, targets::Vector{CartesianIndex{2}}, + env::CTMRGEnv, alg::PatchApprox, + ) + ρ, env = standardize_dualness(ρ, env) + rowrange, colrange = _patch_ranges([source; targets]) + targets_by_row = _twosite_targets_by_row(targets) + mpo = gate_to_mpo(op; trunc = notrunc()) + stringspace = space(mpo[2], 1) + values = zeros(promote_type(scalartype(op), scalartype(ρ), scalartype(env)), length(targets)) + north = _north_boundary_mps(env, source[1], colrange) + plain_north = north + for row in rowrange + W = _row_mpo(ρ, nothing, env, row, colrange) + if haskey(targets_by_row, row) + south = _south_boundary_mps(env, row, colrange) + norm = dot_noconj(south, W, plain_north) + _contract_twosite_target_row!( + values, ρ, mpo, source, targets_by_row[row], north, south, W, colrange + ) + for k in Base.values(targets_by_row[row]) + values[k] /= norm + end + end + row == last(rowrange) && break + A = ρ[row, source[2], 1] + tensor = row == source[1] ? mpo_path_first(A, mpo[1], Val(:south)) : + mpo_path_string(A, stringspace, Val((:north, :south))) + plain_north = _approximate(W, plain_north, alg) + parent(W)[_patch_mps_site(source[2], colrange)] = tensor + north = _approximate(W, north, alg) + end + return values +end + +""" +Contract all targets in one row `W` with a shared `north` and `south` boundary MPS, writing results into `numerators`. +""" +function _contract_twosite_target_row!( + numerators::Vector{<:Number}, ρ::InfinitePEPO, + mpo::AbstractVector{<:AbstractTensorMap}, + source::CartesianIndex{2}, targets::Dict{CartesianIndex{2}, Int}, + north::FiniteMPS, south::FiniteMPS, W::FiniteMPO, colrange::UnitRange{Int}, + ) + N = length(south) + row = first(keys(targets))[1] + source_site = _patch_mps_site(source[2], colrange) + envs = _patch_edge_environments(south, W, north, source_site) + stringspace = space(mpo[2], 1) + + # close the target right at the column of the incoming string + same_col = get(targets, CartesianIndex(row, source[2]), nothing) + if !isnothing(same_col) + target_tensor = mpo_path_last(ρ[row, source[2], 1], mpo[2], Val(:north)) + value = _contract_patch_site(envs, north, south, source_site, target_tensor) + numerators[same_col] = value + end + + # Close targets on the right of the incoming string from left to right + right_targets = [target for target in keys(targets) if target[2] > source[2]] + if !isempty(right_targets) + sort!(right_targets; by = x -> x[2]) + A = ρ[row, source[2], 1] + source_tensor = if row == source[1] + mpo_path_first(A, mpo[1], Val(:east)) + else + mpo_path_string(A, stringspace, Val((:north, :east))) + end + left = envs.lefts[source_site] * + edge_transfermatrix(north.AC[source_site], source_tensor, south.AC[south_site(source_site, N)]) + previous_col = source[2] + for target in right_targets + target_col = target[2] + for col in (previous_col + 1):(target_col - 1) + site = _patch_mps_site(col, colrange) + string_tensor = mpo_path_string(ρ[row, col, 1], stringspace, Val((:west, :east))) + left = left * edge_transfermatrix(north.AR[site], string_tensor, south.AL[south_site(site, N)]) + end + + target_site = _patch_mps_site(target_col, colrange) + target_tensor = mpo_path_last(ρ[row, target_col, 1], mpo[2], Val(:west)) + target_left = left * edge_transfermatrix(north.AR[target_site], target_tensor, south.AL[south_site(target_site, N)]) + value = _contract_transfer_boundaries(target_left, envs.rights[target_site - source_site + 1]) + numerators[targets[target]] = value + + string_tensor = mpo_path_string(ρ[row, target_col, 1], stringspace, Val((:west, :east))) + left = left * edge_transfermatrix(north.AR[target_site], string_tensor, south.AL[south_site(target_site, N)]) + previous_col = target_col + end + end + + # Close targets on the left of the incoming string from right to left + left_targets = [target for target in keys(targets) if target[2] < source[2]] + if !isempty(left_targets) + sort!(left_targets; by = x -> x[2], rev = true) + A = ρ[row, source[2], 1] + source_tensor = if row == source[1] + mpo_path_first(A, mpo[1], Val(:west)) + else + mpo_path_string(A, stringspace, Val((:north, :west))) + end + right = edge_transfermatrix( + north.AC[source_site], source_tensor, south.AC[south_site(source_site, N)] + ) * first(envs.rights) + previous_col = source[2] + for target in left_targets + target_col = target[2] + for col in (previous_col - 1):-1:(target_col + 1) + site = _patch_mps_site(col, colrange) + string_tensor = mpo_path_string(ρ[row, col, 1], stringspace, Val((:east, :west))) + right = edge_transfermatrix(north.AL[site], string_tensor, south.AR[south_site(site, N)]) * right + end + + target_site = _patch_mps_site(target_col, colrange) + target_tensor = mpo_path_last(ρ[row, target_col, 1], mpo[2], Val(:east)) + target_right = edge_transfermatrix(north.AL[target_site], target_tensor, south.AR[south_site(target_site, N)]) * right + value = _contract_transfer_boundaries(envs.lefts[target_site], target_right) + numerators[targets[target]] = value + + string_tensor = mpo_path_string(ρ[row, target_col, 1], stringspace, Val((:east, :west))) + right = edge_transfermatrix(north.AL[target_site], string_tensor, south.AR[south_site(target_site, N)]) * right + previous_col = target_col + end + end + return numerators +end + +""" +Map a PEPO column to its finite-MPS site, accounting for the additional west CTM edge. +""" +_patch_mps_site(col::Int, colrange::UnitRange{Int}) = col - first(colrange) + 2 + +""" +Contract one modified row site between precomputed left and right MPS environments. +""" +function _contract_patch_site( + envs::NamedTuple, north::FiniteMPS, south::FiniteMPS, + site::Int, tensor::MPOTensor, + ) + N = length(south) + left = envs.lefts[site] * + edge_transfermatrix(north.AC[site], tensor, south.AC[south_site(site, N)]) + return _contract_transfer_boundaries(left, envs.rights[site - length(envs.lefts) + 1]) +end diff --git a/src/algorithms/contractions/transfer.jl b/src/algorithms/contractions/transfer.jl index 19ac6db38..b5440f4c7 100644 --- a/src/algorithms/contractions/transfer.jl +++ b/src/algorithms/contractions/transfer.jl @@ -247,3 +247,99 @@ function edge_transfer_right( Ebot[χ_SE D_S; χ_SW] * O[D_W D_S; D_N D_E] end + +""" +Map north-boundary site `k` to the corresponding site in the east-to-west south boundary. +`N` is the finite patch length, including the two boundary sites. +""" +south_site(k::Int, N::Int) = N + 1 - k + +""" +Construct the endpoint identities for a finite north-W-south sandwich. +``` + ┌-←-- north --←-┐ + | | | + L-←---- W ----←-R + | | | + └-→-- south --→-┘ +``` +""" +function _patch_edge_boundaries(south::FiniteMPS, W::FiniteMPO, north::FiniteMPS) + N = length(north) + length(south) == length(W) == N || throw(DimensionMismatch("row and boundary lengths must match")) + left = isomorphism( + storagetype(north.AL[1]), domain(south.AR[N])[1] ⊗ space(W[1], 1)', space(north.AL[1], 1) + ) + right = isomorphism( + storagetype(north.AR[N]), domain(north.AR[N])[1] ⊗ domain(W[N])[2], space(south.AL[1], 1) + ) + return left, right +end + +""" +Build left and right environments outside `site` in the north-W-south sandwich. +`lefts[k]` contains columns before `k`; `rights[k - site + 1]` contains columns after `k`. +Only valid entries are stored, with lengths `site` and `length(north) - site + 1`, respectively. + +Note that sites in `south` are ordered from right (east) to left (west). +""" +function _patch_edge_environments(south::FiniteMPS, W::FiniteMPO, north::FiniteMPS, site::Int) + left, right = _patch_edge_boundaries(south, W, north) + N = length(north) + lefts, rights = [left], [right] + for k in 1:(site - 1) + push!(lefts, last(lefts) * edge_transfermatrix(north.AL[k], W[k], south.AR[south_site(k, N)])) + end + for k in N:-1:(site + 1) + push!(rights, edge_transfermatrix(north.AR[k], W[k], south.AL[south_site(k, N)]) * last(rights)) + end + reverse!(rights) + return (; lefts, rights) +end + +""" +Contract opposite-oriented boundary MPSs without conjugation. +""" +function dot_noconj(south::FiniteMPS, north::FiniteMPS) + N = length(north) + length(south) == N || throw(DimensionMismatch("boundary lengths must match")) + right = isomorphism(storagetype(north.AR[N]), domain(north.AR[N]), space(south.AL[1], 1)) + for k in N:-1:1 + top = k == 1 ? north.AC[k] : north.AR[k] + bottom = k == 1 ? south.AC[N] : south.AL[south_site(k, N)] + right = edge_transfermatrix(top, bottom) * right + end + left = isomorphism(storagetype(right), domain(south.AC[N]), space(north.AC[1], 1)) + return tr(left * right) +end + +""" +Contract a row MPO between opposite-oriented boundary MPSs without conjugation. +""" +function dot_noconj(south::FiniteMPS, W::FiniteMPO, north::FiniteMPS) + left, right = _patch_edge_boundaries(south, W, north) + N = length(north) + for k in N:-1:1 + top = k == 1 ? north.AC[k] : north.AR[k] + bottom = k == 1 ? south.AC[N] : south.AL[south_site(k, N)] + right = edge_transfermatrix(top, W[k], bottom) * right + end + return _contract_transfer_boundaries(left, right) +end + +""" +Contract the left and right transfer-matrix environments to a scalar. +``` + (north) + ┌-←-- 3 --←-┐ + | | + L-←-- 2 --←-R + | | + └-→-- 1 --→-┘ + (south) +``` +""" +function _contract_transfer_boundaries(left::MPSTensor, right::MPSTensor) + # The three bonds close around the patch without crossing + return @tensor left[1 2; 3] * right[3 2; 1] +end diff --git a/src/algorithms/correlator_approx.jl b/src/algorithms/correlator_approx.jl new file mode 100644 index 000000000..25b59cb7c --- /dev/null +++ b/src/algorithms/correlator_approx.jl @@ -0,0 +1,59 @@ +# Approximate finite-patch two-site correlators +# ---------------------------------------------- + +""" +$(SIGNATURES) + +Approximately measure a dense two-site operator between a fixed first site `i` and second sites `js` in a single-layer PEPO. + +- Operator leg 1 acts at `i`, and leg 2 acts at each second site. +- The contraction direction is chosen automatically. + Row-by-row contraction proceeds southward and requires every `j[1] ≥ i[1]`; column-by-column contraction proceeds westward and requires every `j[2] ≤ i[2]`. + If both are possible, consider the smallest rectangle containing `i` and all sites in `js`: use rows if it is square or wider than tall, and columns otherwise. + If neither direction is possible, throw an `ArgumentError`. +- The collection `js` must be nonempty, contain no duplicates, and exclude `i`. + A collection returns a vector in `vec(js)` order; a single second site returns a scalar. +- For row-by-row contraction, all patches have the same width, covering the columns of `i` and all sites in `js`, but each patch ends at the row of the measured site `j`. + For column-by-column contraction, all patches cover the same rows, but each patch ends at the column of `j`. + Each patch has its own normalization, calculated without the operator using the same contraction arrangement. + The last row or column is contracted without further truncation against the CTMRG boundary immediately beyond it. +- `trunc` controls boundary-MPS truncation; by default, it limits the rank to the largest CTMRG boundary dimension. + After each zipup step, `maxiter` DMRG sweeps refine the result (default 1; use 0 to disable refinement). +""" +function correlator_approx( + ρ::InfinitePEPO, op::AbstractTensorMap, + i::CartesianIndex{2}, j::CartesianIndex{2}, env::CTMRGEnv; + trunc = _approx_trunc(env), maxiter::Int = 1, + ) + return only(correlator_approx(ρ, op, i, j:j, env; trunc, maxiter)) +end + +function correlator_approx( + ρ::InfinitePEPO, op::AbstractTensorMap, + i::CartesianIndex{2}, js::CoordCollection{2}, env::CTMRGEnv; + trunc = _approx_trunc(env), maxiter::Int = 1, + ) + isempty(js) && throw(ArgumentError("correlator_approx requires at least one second site")) + allunique(js) || throw(ArgumentError("second sites should be unique")) + i in js && throw(ArgumentError("second sites should be distinct from the first site")) + rowrange, colrange = _patch_ranges([i; vec(js)]) + rows, columns = first(rowrange) == i[1], last(colrange) == i[2] + rows || columns || throw(ArgumentError("no valid sweep: second sites must all be at or south of the first row, or all at or west of the first column")) + direction = rows && (!columns || length(colrange) >= length(rowrange)) ? :rows : :columns + return _correlator_approx( + ρ, op, i, collect(vec(js)), env, + PatchApprox(Zipup(; trunc), _approx_dmrg(maxiter)), direction + ) +end + +""" +Group targets by row, retaining each target's original result position. +""" +function _twosite_targets_by_row(targets::Vector{CartesianIndex{2}}) + targets_by_row = Dict{Int, Dict{CartesianIndex{2}, Int}}() + for (position, target) in enumerate(targets) + row_targets = get!(Dict{CartesianIndex{2}, Int}, targets_by_row, target[1]) + row_targets[target] = position + end + return targets_by_row +end diff --git a/src/algorithms/expval_approx.jl b/src/algorithms/expval_approx.jl new file mode 100644 index 000000000..813cc3201 --- /dev/null +++ b/src/algorithms/expval_approx.jl @@ -0,0 +1,77 @@ +# Approximate finite-patch expectation values +# -------------------------------------------- + +""" +$(SIGNATURES) + +Approximately measure a `LocalOperator` in a single-layer PEPO using finite boundary MPS zipup sweeps. +Each term is normalized in its own enclosing patch, and the resulting contributions are summed. +Dense terms are reordered into a column-wise snake and decomposed into MPOs without truncation. +Explicit [`MPOTerm`](@ref) factors retain their given order and are connected by non-self-intersecting nearest-neighbor paths. +Routing tries horizontal-first and then vertical-first shortest connections; more complicated explicit MPO routes raise a not-implemented error. +Single-site terms are evaluated exactly, and an empty operator returns zero. + +- By default, `direction = :auto` selects north-to-south sweeps for wide and square patches, and east-to-west for tall patches, separately for each term. + Specify `direction = :rows` or `:columns` to override this choice. +- The zipup truncation is controlled by `trunc`, which defaults to `truncrank(χ)` with `χ` the largest CTMRG boundary dimension. + This keyword does not truncate the operator decomposition. +- After each zipup step, the result is refined by a single-site DMRG approximation step with `maxiter` sweeps. + Set `maxiter = 0` to disable this refinement. +""" +function expectation_value_approx( + ρ::InfinitePEPO, O::LocalOperator, env::CTMRGEnv; + trunc = _approx_trunc(env), maxiter::Int = 1, direction::Symbol = :auto, + ) + _check_patch_inputs(ρ, direction) + checklattice(ρ, O) + alg = PatchApprox(Zipup(; trunc), _approx_dmrg(maxiter)) + isempty(O.terms) && return zero(promote_type(scalartype(ρ), scalartype(env))) + term_vals = map(collect(O.terms)) do (sites, term) + return _local_expectation_value_approx(sites, ρ, term, env, alg, direction) + end + return sum(term_vals) +end + +""" +Evaluate a dense or MPO term, using exact single-site contractions and routing larger terms. +""" +function _local_expectation_value_approx( + sites::Vector{CartesianIndex{2}}, ρ::InfinitePEPO, + term::Union{AbstractTensorMap, MPOTerm}, env::CTMRGEnv, + alg::PatchApprox, direction::Symbol, + ) + if _local_term_iszero(term) + return zero(promote_type(scalartype(ρ), scalartype(env), _local_term_scalartype(term))) + end + if length(sites) == 1 + op = term isa AbstractTensorMap ? term : only(term) + return local_expectation_value(sites, ρ, op, env) + end + routed = _route_mpo_term(sites, term, physicalspace(ρ)) + return _expectation_value_approx(ρ, routed, env, alg, direction) +end + + +""" +Return the largest CTMRG boundary-space dimension appearing in the corner tensors. +""" +function _ctmrg_boundary_chi(env::CTMRGEnv) + χ = 0 + for C in env.corners + χ = max(χ, dim(space(C, 1)), dim(space(C, 2))) + end + return χ +end + +""" +Construct the default rank truncation from the largest CTMRG boundary dimension. +""" +_approx_trunc(env::CTMRGEnv) = truncrank(_ctmrg_boundary_chi(env)) + +""" +Construct the optional one-site DMRG refinement, or disable refinement for zero iterations. +""" +function _approx_dmrg(maxiter::Int) + maxiter >= 0 || throw(ArgumentError("maxiter should be nonnegative")) + return iszero(maxiter) ? nothing : DMRG(; maxiter, verbosity = 0) +end diff --git a/src/operators/localoperator.jl b/src/operators/localoperator.jl index c35e5ddb9..35e632637 100644 --- a/src/operators/localoperator.jl +++ b/src/operators/localoperator.jl @@ -1,10 +1,19 @@ # Hamiltonian consisting of local terms # ------------------------------------- +""" +An open-boundary MPO represented by an ordered vector of tensor maps. +The endpoint partitions are `(1, 2)` and `(2, 1)`, with `(2, 2)` tensors in between; a one-site MPO has partition `(1, 1)`. +Dense tensors can be converted with `gate_to_mpo`. +MPO terms are supported by `expectation_value_approx`; exact expectation values and vectors of tensor-product factors are not implemented. +""" +const MPOTerm{T} = AbstractVector{T} where {T <: AbstractTensorMap} + """ $(TYPEDEF) A sum of local operators acting on a lattice. The lattice is stored as a matrix of vector spaces, and the terms are stored as a `Dict` of indices mapping to operators. +Terms can be dense tensor maps or `MPOTerm`s, whose site and factor ordering is preserved. ## Fields @@ -50,10 +59,23 @@ end # Default to Any for eltype: needs to be abstract anyways so not that much to gain LocalOperator(lattice, terms) = LocalOperator{Any}(lattice, terms) LocalOperator(lattice, terms::Pair...) = LocalOperator(lattice, terms) -# TODO: add terms beyond AbstractTensorMap -# e.g. tensor product of 1-site operators, MPOs -add_term!(operator::LocalOperator, inds::Tuple, term::AbstractTensorMap) = add_term!(operator, collect(inds), term) -add_term!(operator::LocalOperator, inds::Vector, term::AbstractTensorMap) = add_term!(operator, map(CartesianIndex{2}, inds), term) + +""" +Sort operator sites using the default `CartesianIndex` ordering. +When the order changes, permute the corresponding output and input physical legs by the same ordering. +""" +function _sort_op_sites(sites::Vector{CartesianIndex{2}}, op::AbstractTensorMap) + issorted(sites) && return sites, op + order = sortperm(sites) + sites′ = sites[order] + op′ = permute(op, (Tuple(order), Tuple(order) .+ numout(op))) + return sites′, op′ +end + +add_term!(operator::LocalOperator, inds::Tuple, term::Union{AbstractTensorMap, MPOTerm}; kwargs...) = + add_term!(operator, collect(inds), term; kwargs...) +add_term!(operator::LocalOperator, inds::AbstractVector, term::Union{AbstractTensorMap, MPOTerm}; kwargs...) = + add_term!(operator, CartesianIndex{2}[CartesianIndex{2}(ind) for ind in inds], term; kwargs...) function add_term!( operator::LocalOperator, inds::Vector{CartesianIndex{2}}, term::AbstractTensorMap; atol = zero(real(scalartype(term))), @@ -68,17 +90,14 @@ function add_term!( end norm(term) <= atol && return operator # skip adding negligible terms - # permute input - if !issorted(inds) - I = sortperm(inds) - inds = inds[I] - term = permute(term, (Tuple(I), Tuple(I) .+ numout(term))) - end + inds, term = _sort_op_sites(copy(inds), term) # translate coordinates _shift_into_unitcell!(inds, size(operator)) if haskey(operator.terms, inds) + operator.terms[inds] isa MPOTerm && + throw(ArgumentError("Accumulating terms with the same sites is not implemented for MPO terms.")) operator.terms[inds] = VI.add!!(operator.terms[inds], term) else operator.terms[inds] = term @@ -87,6 +106,51 @@ function add_term!( return operator end +""" +Validate an ordered MPO's sites, tensor partitions, physical spaces, and adjacent bond spaces. +""" +function _validate_mpo_term(sites, term::MPOTerm, lattice::AbstractMatrix{<:ElementarySpace}) + isempty(term) && throw(ArgumentError("An MPO term should contain at least one tensor.")) + length(sites) == length(term) || + throw(ArgumentError("The MPO should contain one tensor for every operator site.")) + allunique(sites) || throw(ArgumentError("`inds` should not contain repeated coordinates.")) + + N = length(term) + for (k, op) in enumerate(term) + expected = N == 1 ? (1, 1) : k == 1 ? (1, 2) : k == N ? (2, 1) : (2, 2) + (numout(op), numin(op)) == expected || + throw(ArgumentError("MPO tensor $k should have partition $expected.")) + site = sites[k] + physical = lattice[mod1(site[1], size(lattice, 1)), mod1(site[2], size(lattice, 2))] + physical == codomain(op)[numout(op)] == domain(op)[1] || + throw(SpaceMismatch("MPO physical space does not match lattice site $site.")) + end + for k in 1:(N - 1) + domain(term[k])[numin(term[k])] == codomain(term[k + 1])[1] || + throw(SpaceMismatch("Incompatible MPO bond spaces between tensors $k and $(k + 1).")) + end + return nothing +end + +""" +Insert an ordered MPO term, copying its coordinate and factor containers before translation. +""" +function add_term!(operator::LocalOperator, inds::Vector{CartesianIndex{2}}, term::MPOTerm) + _validate_mpo_term(inds, term, physicalspace(operator)) + _local_term_iszero(term) && return operator + inds = _shift_into_unitcell!(copy(inds), size(operator)) + haskey(operator.terms, inds) && + throw(ArgumentError("Accumulating terms with the same sites is not implemented for MPO terms.")) + operator.terms[inds] = collect(term) + return operator +end + +""" +Detect an identically zero dense tensor or an MPO containing a zero factor. +""" +_local_term_iszero(term::AbstractTensorMap) = iszero(norm(term)) +_local_term_iszero(term::MPOTerm) = any(_local_term_iszero, term) + """ checklattice(Bool, args...) @@ -154,16 +218,27 @@ Base.eltype(::Type{LocalOperator{O, S}}) where {O, S} = O # Real and imaginary part # ----------------------- function Base.real(O::LocalOperator) + any(term -> term isa MPOTerm, values(O.terms)) && + throw(ArgumentError("Taking the real part is not implemented for MPO terms.")) return LocalOperator(O.lattice, (sites => real(op) for (sites, op) in O.terms)...) end function Base.imag(O::LocalOperator) + any(term -> term isa MPOTerm, values(O.terms)) && + throw(ArgumentError("Taking the imaginary part is not implemented for MPO terms.")) return LocalOperator(O.lattice, (sites => imag(op) for (sites, op) in O.terms)...) end # Linear Algebra # -------------- +""" +Scale a local term, applying an MPO's scalar to its first factor only. +""" +_scale_local_term(α::Number, term::AbstractTensorMap) = α * term +_scale_local_term(α::Number, term::MPOTerm) = + AbstractTensorMap[k == 1 ? α * tensor : tensor for (k, tensor) in enumerate(term)] + Base.:*(α::Number, O::LocalOperator) = - LocalOperator(physicalspace(O), inds => α * operator for (inds, operator) in O.terms) + LocalOperator(physicalspace(O), inds => _scale_local_term(α, term) for (inds, term) in O.terms) Base.:*(O::LocalOperator, α::Number) = α * O Base.:/(O::LocalOperator, α::Number) = O * inv(α) @@ -171,7 +246,16 @@ Base.:\(α::Number, O::LocalOperator) = inv(α) * O function Base.:+(O1::LocalOperator, O2::LocalOperator) checklattice(O1, O2) - return LocalOperator(physicalspace(O1), mergewith(VI.add, O1.terms, O2.terms)) + return LocalOperator(physicalspace(O1), mergewith(_add_local_terms, O1.terms, O2.terms)) +end + +""" +Accumulate dense terms while rejecting addition involving an MPO at the same sites. +""" +function _add_local_terms(term1, term2) + (term1 isa MPOTerm || term2 isa MPOTerm) && + throw(ArgumentError("Accumulating terms with the same sites is not implemented for MPO terms.")) + return VI.add(term1, term2) end Base.:-(O::LocalOperator) = -1 * O @@ -182,9 +266,14 @@ Base.:-(O1::LocalOperator, O2::LocalOperator) = O1 + (-O2) # Since we allow abstract types in T, value and type domain might not match function VI.scalartype(operator::LocalOperator) - return promote_type((scalartype(term[2]) for term in operator.terms)...) + return promote_type((_local_term_scalartype(term) for term in values(operator.terms))...) end +""" +Return the promoted scalar type of a dense tensor or the factors of an MPO. +""" +_local_term_scalartype(term::AbstractTensorMap) = scalartype(term) +_local_term_scalartype(term::MPOTerm) = promote_type((scalartype(tensor) for tensor in term)...) # Equivalence # ----------- diff --git a/src/utility/indexing.jl b/src/utility/indexing.jl index cc5fa1249..de0b166b8 100644 --- a/src/utility/indexing.jl +++ b/src/utility/indexing.jl @@ -1,3 +1,9 @@ +""" +An ordered sequence of sites on a two-dimensional lattice. +Nearest-neighbor connectivity and nonintersection are checked by routing functions, not enforced by this alias. +""" +const LatticePath = AbstractVector{CartesianIndex{2}} + # Get next and previous directional CTMRG environment index, respecting periodicity _next(i, total) = mod1(i + 1, total) _prev(i, total) = mod1(i - 1, total) diff --git a/test/testsuite/TestSuite.jl b/test/testsuite/TestSuite.jl index 018a68a82..29967b264 100644 --- a/test/testsuite/TestSuite.jl +++ b/test/testsuite/TestSuite.jl @@ -239,6 +239,30 @@ module ToolboxDensityMatrices end using .ToolboxDensityMatrices +module ToolboxExpvalApprox + include("toolbox/expval_approx.jl") + export toolbox_expval_approx, toolbox_expval_approx_localoperator +end +using .ToolboxExpvalApprox + +module ToolboxCorrelatorApprox + include("toolbox/correlator_approx.jl") + export toolbox_correlator_approx +end +using .ToolboxCorrelatorApprox + +module ToolboxCorrelatorApproxPhys + include("toolbox/correlator_approx_phys.jl") + export toolbox_correlator_approx_phys +end +using .ToolboxCorrelatorApproxPhys + +module ToolboxMPORouting + include("toolbox/mpo_routing.jl") + export toolbox_mpo_routing_fusers, toolbox_mpo_routing_identities, toolbox_mpo_routing_paths +end +using .ToolboxMPORouting + # Utility # ------- module UtilityCorrelator diff --git a/test/testsuite/toolbox/correlator_approx.jl b/test/testsuite/toolbox/correlator_approx.jl new file mode 100644 index 000000000..fad3f9c02 --- /dev/null +++ b/test/testsuite/toolbox/correlator_approx.jl @@ -0,0 +1,104 @@ +using TensorKit +using PEPSKit +using MPSKit +using Test +using Adapt +using Random + +const CI = CartesianIndex + +""" +Contract `⟨op⟩` between `i` and each site in `js` independently without caching, using fixed width and a patch ending at each target row. +""" +function _adaptive_patch_reference(op::AbstractTensorMap, i::CI{2}, js, ρ, env) + lattice = physicalspace(ρ) + observables = [PEPSKit._route_mpo_term([i, j], op, lattice) for j in js] + _, colrange = PEPSKit._patch_ranges([i; js]) + alg = PEPSKit.PatchApprox(Zipup(; trunc = notrunc()), nothing) + return map(observables, js) do observable, j + rowrange = i[1]:j[1] + norm = PEPSKit._contract_patch_rows(ρ, nothing, env, rowrange, colrange, alg) + numerator = PEPSKit._contract_patch_rows(ρ, observable, env, rowrange, colrange, alg) + return numerator / norm + end +end + +# Exercise both U(1) symmetry and fermionic signs with small physical, virtual, and boundary spaces. +spaces = Dict( + U1Irrep => ( + U1Space(1 => 2, -1 => 1), + U1Space(0 => 1, 1 => 1), + U1Space(0 => 1, 1 => 1), + ), + FermionParity => ( + Vect[FermionParity](0 => 1, 1 => 1), + Vect[FermionParity](0 => 1, 1 => 1), + Vect[FermionParity](0 => 1, 1 => 1), + ) +) + +""" +Check approximate correlators, sweep selection, and adaptive patch normalization. +""" +function toolbox_correlator_approx(AT) + return @testset "Single-layer PEPO ($S) ($AT)" for S in keys(spaces) + # Use a reproducible random PEPO on a rectangular unit cell and disable boundary truncation. + Random.seed!(1234) + d, D, χ = spaces[S] + ρ = adapt(AT, InfinitePEPO(d, D; unitcell = (2, 3, 1))) + env = CTMRGEnv(InfinitePartitionFunction(ρ), χ) + trunc = notrunc() + i = CI(-1, 0) + # Unsorted targets include both same-row neighbors of i and targets to the east, west, and directly south on row 1. + # Row 0 has no measurement targets but is still contracted, testing MPO-string propagation through a row without targets. + js = [CI(1, 1), CI(-1, -1), CI(1, 0), CI(-1, 1), CI(1, -1)] + + # The normalized expectation of the identity must be one at every target. + id² = adapt(AT, isomorphism(d, d) ⊗ isomorphism(d, d)) + @test correlator_approx(ρ, id², i, js, env; trunc, maxiter = 0) ≈ ones(length(js)) + + # This target distribution only permits row sweeps. + # Compare against independent contractions in target order. + O² = adapt(AT, rand(ComplexF64, d^2, d^2)) + vals_ref = _adaptive_patch_reference(O², i, js, ρ, env) + vals_rows = correlator_approx(ρ, O², i, js, env; trunc, maxiter = 0) + @test vals_rows ≈ vals_ref + + # CartesianIndices are flattened exactly as in correlator. + grid = CartesianIndices((1:1, -1:1)) + @test correlator_approx(ρ, O², i, grid, env; trunc, maxiter = 0) ≈ + vals_ref[[5, 3, 1]] + + # The rotated distribution only permits column sweeps. + # Check rotation and automatic selection together. + cell = size(ρ)[1:2] + ic = PEPSKit.siterotr90(i, cell) + jcs = PEPSKit.siterotr90.(js, Ref(cell)) + ρc, envc = rotr90(ρ), rotr90(env) + @test correlator_approx(ρc, O², ic, jcs, envc; trunc, maxiter = 0) ≈ vals_rows + + # Identity normalization must survive truncation for every target row and both sweep orientations. + for maxiter in (0, 1) + @test correlator_approx(ρ, id², i, js, env; trunc = truncrank(2), maxiter) ≈ ones(length(js)) + @test correlator_approx(ρc, id², ic, jcs, envc; trunc = truncrank(2), maxiter) ≈ ones(length(js)) + end + + # Test auto sweep direction choice in for southwest targets. + selection_trunc = truncrank(2) + selection_alg = PEPSKit.PatchApprox(Zipup(; trunc = selection_trunc), nothing) + for (offset, direction) in ((CI(1, -2), :rows), (CI(2, -1), :columns), (CI(1, -1), :rows)) + j = i + offset + expected = only(PEPSKit._correlator_approx(ρ, O², i, [j], env, selection_alg, direction)) + @test correlator_approx(ρ, O², i, j, env; trunc = selection_trunc, maxiter = 0) == expected + end + + # Extending the depth must preserve earlier values when the width and selected direction stay fixed. + extended_js = [js; CI(2, 0)] + baseline = correlator_approx(ρ, O², i, js, env; trunc = selection_trunc, maxiter = 0) + extended = correlator_approx(ρ, O², i, extended_js, env; trunc = selection_trunc, maxiter = 0) + @test extended[1:length(js)] ≈ baseline + extended_jcs = PEPSKit.siterotr90.(extended_js, Ref(cell)) + extended_columns = correlator_approx(ρc, O², ic, extended_jcs, envc; trunc = selection_trunc, maxiter = 0) + @test extended_columns ≈ extended + end +end diff --git a/test/testsuite/toolbox/correlator_approx_phys.jl b/test/testsuite/toolbox/correlator_approx_phys.jl new file mode 100644 index 000000000..d54a56400 --- /dev/null +++ b/test/testsuite/toolbox/correlator_approx_phys.jl @@ -0,0 +1,42 @@ +using TensorKit +using PEPSKit +using Test +using Adapt +using TensorKitTensors.SpinOperators: S_exchange + +const CI = CartesianIndex + +""" +Compare approximate and exact correlators in an evolved physical state. +""" +function toolbox_correlator_approx_phys(AT) + return @testset "correlator_approx for physical state ($AT)" begin + Nr, Nc = 2, 2 + lattice = InfiniteSquare(Nr, Nc) + sym = Trivial + ham = adapt(AT, j1_j2_model(Float64, sym, lattice; J1 = 1.0, J2 = 0.0, sublattice = false)) + op = adapt(AT, S_exchange(Float64, sym)) + lattice = physicalspace(ham) + + ρ = PEPSKit.infinite_temperature_density_matrix(ham) + state_trunc = truncrank(4) & truncerror(; atol = 1.0e-12) + su_alg = SimpleUpdate(; trunc = state_trunc, purified = false) + ρ, = time_evolve(ρ, ham, 1.0e-2, 50, su_alg, SUWeight(ρ)) + + network = InfinitePartitionFunction(ρ) + env = initialize_ctmrg_environment(network, ProductStateInitialization()) + env_trunc = truncrank(8) & truncerror(; atol = 1.0e-12) + env, = leading_boundary(env, network; alg = :SequentialCTMRG, trunc = env_trunc) + + i = CI(1, 1) + js = [CI(1, 2), CI(2, 0), CI(2, 1), CI(2, 2), CI(3, 2)] + cor_exact = map(js) do j + O = LocalOperator(lattice, (i, j) => op) + return expectation_value(ρ, O, env) + end + cor_trunc = correlator_approx(ρ, op, i, js, env; trunc = env_trunc, maxiter = 1) + @info "Exact:" cor_exact + @info "Approx:" cor_trunc + @test cor_trunc ≈ cor_exact rtol = 1.0e-3 + end +end diff --git a/test/testsuite/toolbox/expval_approx.jl b/test/testsuite/toolbox/expval_approx.jl new file mode 100644 index 000000000..699eca1cd --- /dev/null +++ b/test/testsuite/toolbox/expval_approx.jl @@ -0,0 +1,120 @@ +using TensorKit +using PEPSKit +using MPSKit +using Test +using Adapt +using Random + +const CI = CartesianIndex + +spaces = Dict( + U1Irrep => ( + U1Space(1 => 2, -1 => 1), + U1Space(1 => 1, 0 => 1, -1 => 2), + U1Space(1 => 1, 0 => 1, -1 => 2), + ), + FermionParity => ( + Vect[FermionParity](0 => 1, 1 => 1), + Vect[FermionParity](0 => 1, 1 => 2), + Vect[FermionParity](0 => 2, 1 => 2), + ), +) + +sites_list = ( + [CI(1, 1), CI(1, 2)], # horizontal + [CI(1, 1), CI(2, 1)], # vertical + [CI(1, 1), CI(2, 2)], # turned + [CI(2, 2), CI(1, 1)], # reversed turned + [CI(2, 1), CI(1, 1), CI(1, 2), CI(2, 2)], # U-shaped +) + +""" +Check approximate expectation values against exact single-layer PEPO contractions. +""" +function toolbox_expval_approx(AT) + return @testset "Single-layer PEPO ($S) ($AT)" for S in keys(spaces) + Random.seed!(1234) + + d, D, χ = spaces[S] + ρ = adapt(AT, InfinitePEPO(d, D; unitcell = (2, 2, 1))) + env = CTMRGEnv(InfinitePartitionFunction(ρ), χ) + lattice = physicalspace(ρ) + trunc = notrunc() + + # Dense terms may change MPO order when forming a snake; explicit MPOs must keep theirs. + # Exact contractions expose permutation or braiding errors in both paths and sweep orientations. + for sites in sites_list + n = length(sites) + op = adapt(AT, randn(ComplexF64, d^n → d^n)) + mpo = PEPSKit.gate_to_mpo(op; trunc) + dense = LocalOperator(lattice, sites => op) + factorized = LocalOperator(lattice, sites => mpo) + exact = expectation_value(ρ, dense, env) + for observable in (dense, factorized), direction in (:rows, :columns) + @test expectation_value_approx( + ρ, observable, env; trunc, maxiter = 0, direction + ) ≈ exact + end + # Square patches must default to row sweeps. + if sites == sites_list[3] + auto = expectation_value_approx(ρ, dense, env; trunc = truncrank(2), maxiter = 0) + rows = expectation_value_approx(ρ, dense, env; trunc = truncrank(2), maxiter = 0, direction = :rows) + @test auto == rows + end + end + end +end + +""" +Check single-site dispatch, sparse snakes, and sums of independently normalized terms. +""" +function toolbox_expval_approx_localoperator(AT) + return @testset "LocalOperator interface ($AT)" begin + Random.seed!(4321) + d = Vect[FermionParity](0 => 1, 1 => 1) + ρ = adapt(AT, InfinitePEPO(d, d; unitcell = (2, 2, 1))) + env = CTMRGEnv(InfinitePartitionFunction(ρ), d) + lattice = physicalspace(ρ) + trunc = notrunc() + + # Single-site terms bypass MPO decomposition, which requires at least two sites. + op1 = adapt(AT, randn(ComplexF64, d → d)) + one_site = LocalOperator(lattice, [CI(1, 1)] => op1) + one_factor = LocalOperator(lattice, [CI(1, 1)] => [op1]) + for observable in (one_site, one_factor) + @test expectation_value_approx(ρ, observable, env; trunc, maxiter = 0) ≈ + expectation_value(ρ, one_site, env) + end + + # A greedy horizontal-first route revisits a site for this support. + # The dense snake must connect outside the column's endpoints and preserve the operator. + op3 = adapt(AT, randn(ComplexF64, d^3 → d^3)) + snake = LocalOperator(lattice, [CI(1, 0), CI(0, 1), CI(2, 1)] => op3) + @test expectation_value_approx(ρ, snake, env; trunc, maxiter = 0) ≈ + expectation_value(ρ, snake, env) + + # These terms use different normalization patches; the gapped MPO also inserts a string. + # Non-unit-cell coordinates and both sweeps exercise translation and rotation of the expanded path. + op2 = adapt(AT, randn(ComplexF64, d^2 → d^2)) + sites2 = [CI(-1, 0), CI(-1, 2)] + mixed = LocalOperator(lattice, [CI(1, 1)] => op1, sites2 => PEPSKit.gate_to_mpo(op2; trunc)) + dense_sum = LocalOperator(lattice, [CI(1, 1)] => op1, sites2 => op2) + for direction in (:rows, :columns) + @test expectation_value_approx(ρ, mixed, env; trunc, maxiter = 0, direction) ≈ + expectation_value(ρ, dense_sum, env) + end + + # Complex scaling leaves real factors in the MPO; expansion must retain a valid mixed-scalar container. + real_op = real(op2) + real_mpo = LocalOperator(lattice, sites2 => PEPSKit.gate_to_mpo(real_op; trunc)) + coefficient = 2 + 3im + @test expectation_value_approx(ρ, coefficient * real_mpo, env; trunc, maxiter = 0) ≈ + coefficient * expectation_value(ρ, LocalOperator(lattice, sites2 => real_op), env) + + # Constructors discard zero terms, so an empty sum must produce a scalar zero without contracting. + @test iszero(expectation_value_approx(ρ, LocalOperator(lattice), env)) + # Reject a mismatched unit cell even when the individual operator's physical space matches. + wrong_lattice = LocalOperator(fill(d, 1, 1), [CI(1, 1)] => op1) + @test_throws ArgumentError expectation_value_approx(ρ, wrong_lattice, env) + end +end diff --git a/test/testsuite/toolbox/mpo_routing.jl b/test/testsuite/toolbox/mpo_routing.jl new file mode 100644 index 000000000..96cbf851d --- /dev/null +++ b/test/testsuite/toolbox/mpo_routing.jl @@ -0,0 +1,153 @@ +using PEPSKit +using TensorKit +using Test +using Adapt +using Random + +const directions = (:north, :east, :south, :west) + +spaces = Dict( + U1Irrep => ( + Rep[U₁](0 => 1, 1 => 1), + Rep[U₁](0 => 1, 1 => 1, -1 => 1), + Rep[U₁](1 => 1), + ), + FermionParity => ( + Vect[FermionParity](0 => 1, 1 => 1), + Vect[FermionParity](0 => 1, 1 => 1), + Vect[FermionParity](1 => 1), + ), +) + +""" +Check fermionic fuser flips and twists for each MPO routing direction. +""" +function toolbox_mpo_routing_fusers(AT) + return @testset "Fuser flips and twists ($AT)" begin + d = Vect[FermionParity](0 => 1, 1 => 1) + D = Vect[FermionParity](0 => 1, 1 => 1) + ρ = adapt(AT, InfinitePEPO(d, D; unitcell = (2, 2, 1))) + A = ρ[1, 1, 1] + Bh = ρ[1, 2, 1] + Bv = ρ[2, 1, 1] + A′ = PEPSKit.twistdual(A, 2) + Bh′ = PEPSKit.twistdual(Bh, 2) + Bv′ = PEPSKit.twistdual(Bv, 2) + + op = adapt(AT, randn(ComplexF64, d^2 → d^2)) + mpo = PEPSKit.gate_to_mpo(op; trunc = notrunc()) + + first_tensor = PEPSKit.mpo_path_first(A, first(mpo), Val(:east)) + last_tensor = PEPSKit.mpo_path_last(Bh, last(mpo), Val(:west)) + @tensor exact[W1 S1 S2; N1 N2 E2] := op[po1 po2; pi1 pi2] * + A′[pi1 po1; N1 x S1 W1] * Bh′[pi2 po2; N2 E2 S2 x] + @tensor routed[W1 S1 S2; N1 N2 E2] := + first_tensor[W1 S1; N1 x] * last_tensor[x S2; N2 E2] + @test routed ≈ exact + + first_tensor = PEPSKit.mpo_path_first(Bh, first(mpo), Val(:west)) + last_tensor = PEPSKit.mpo_path_last(A, last(mpo), Val(:east)) + @tensor exact[W2 S1 S2; N1 E1 N2] := op[po1 po2; pi1 pi2] * + Bh′[pi1 po1; N1 E1 S1 x] * A′[pi2 po2; N2 x S2 W2] + @tensor routed[W2 S1 S2; N1 E1 N2] := + first_tensor[x S1; N1 E1] * last_tensor[W2 S2; N2 x] + @test routed ≈ exact + + first_tensor = PEPSKit.mpo_path_first(A, first(mpo), Val(:south)) + last_tensor = PEPSKit.mpo_path_last(Bv, last(mpo), Val(:north)) + @tensor exact[W1 W2 S2; N1 E1 E2] := op[po1 po2; pi1 pi2] * + A′[pi1 po1; N1 E1 x W1] * Bv′[pi2 po2; x E2 S2 W2] + @tensor routed[W1 W2 S2; N1 E1 E2] := + first_tensor[W1 x; N1 E1] * last_tensor[W2 S2; x E2] + @test routed ≈ exact + + first_tensor = PEPSKit.mpo_path_first(Bv, first(mpo), Val(:north)) + last_tensor = PEPSKit.mpo_path_last(A, last(mpo), Val(:south)) + @tensor exact[W1 S1 W2; E1 N2 E2] := op[po1 po2; pi1 pi2] * + Bv′[pi1 po1; x E1 S1 W1] * A′[pi2 po2; N2 E2 x W2] + @tensor routed[W1 S1 W2; E1 N2 E2] := + first_tensor[W1 S1; x E1] * last_tensor[W2 x; N2 E2] + @test routed ≈ exact + end +end + +""" +Check MPO endpoint shapes and string-routing identities for each symmetry. +""" +function toolbox_mpo_routing_identities(AT) + return @testset "Routing identities ($S) ($AT)" for S in keys(spaces) + Random.seed!(1234) + d, D, stringspace = spaces[S] + ρ = adapt(AT, InfinitePEPO(d, D; unitcell = (1, 1, 1))) + op = adapt(AT, rand(ComplexF64, d^2 → d^2)) + mpo = PEPSKit.gate_to_mpo(op; trunc = notrunc()) + + A = ρ[1, 1, 1] + for direction in directions + first_tensor = PEPSKit.mpo_path_first(A, first(mpo), Val(direction)) + last_tensor = PEPSKit.mpo_path_last(A, last(mpo), Val(direction)) + @test (numout(first_tensor), numin(first_tensor)) == (2, 2) + @test (numout(last_tensor), numin(last_tensor)) == (2, 2) + end + + middle = adapt(AT, TensorMap(TensorKit.BraidingTensor{ComplexF64}(d, stringspace))) + for incoming in directions, outgoing in directions + incoming == outgoing && continue + tensor = PEPSKit.mpo_path_middle(A, middle, Val((incoming, outgoing))) + string_tensor = PEPSKit.mpo_path_string( + A, stringspace, Val((incoming, outgoing)) + ) + @test (numout(tensor), numin(tensor)) == (2, 2) + @test string_tensor ≈ tensor + end + end +end + +""" +Check routing choices, dense leg permutations, and periodic string insertion. +""" +function toolbox_mpo_routing_paths(AT) + return @testset "MPO paths ($AT)" begin + Random.seed!(1234) + CI = CartesianIndex + + # The first horizontal L crosses a later operator site, so routing must try the vertical L. + # Conversely, a collinear out-of-order MPO cannot be connected by either supported L path. + fallback_sites = CI.([(1, 1), (3, 3), (1, 2)]) + @test PEPSKit._ordered_mpo_path(fallback_sites) == + CI.([(1, 1), (2, 1), (3, 1), (3, 2), (3, 3), (2, 3), (1, 3), (1, 2)]) + @test_throws r"routing this MPO ordering is not implemented" PEPSKit._ordered_mpo_path( + CI.([(1, 1), (1, 3), (1, 2)]) + ) + + # Reconstructing on unequal physical spaces catches mismatches between the snake order and either group of permuted physical legs. + lattice = [ℂ^2 ℂ^3; ℂ^1 ℂ^2] + sites = CI.([(2, 2), (1, 1), (1, 2), (2, 1)]) + physical = foldl(⊗, lattice[sites]) + op = adapt(AT, randn(ComplexF64, physical ← physical)) + path, expanded = PEPSKit._route_mpo_term(sites, op, lattice) + @test path == sites[[2, 4, 1, 3]] + @tensor reconstructed[p1 p2 p3 p4; q1 q2 q3 q4] := + expanded[1][p1; q1 a] * expanded[2][a p2; q2 b] * + expanded[3][b p3; q3 c] * expanded[4][c p4; q4] + @test reconstructed ≈ permute(op, ((2, 4, 1, 3), (6, 8, 5, 7))) + + # Intermediate sites have different physical spaces from the endpoints and lie outside the unit cell. + # Each inserted braid must therefore use periodic lookup and the MPO's storage type. + @testset "Periodic strings ($S)" for S in keys(spaces) + physical, _, _ = spaces[S] + intermediate = physical ⊕ physical + lattice = [physical intermediate; intermediate physical] + op = adapt(AT, randn(ComplexF64, physical^2 ← physical^2)) + mpo = PEPSKit.gate_to_mpo(op; trunc = notrunc()) + path, expanded = PEPSKit._route_mpo_term(CI.([(0, 0), (2, 2)]), mpo, lattice) + for k in 2:(length(path) - 1) + site = path[k] + V = lattice[mod1(site[1], 2), mod1(site[2], 2)] + braid = adapt(AT, TensorMap(TensorKit.BraidingTensor{ComplexF64}(V, space(mpo[2], 1)))) + @test expanded[k] ≈ braid + end + @test all(t -> storagetype(t) == storagetype(mpo[1]), expanded) + end + end +end diff --git a/test/toolbox/correlator_approx.jl b/test/toolbox/correlator_approx.jl new file mode 100644 index 000000000..9ad275d5d --- /dev/null +++ b/test/toolbox/correlator_approx.jl @@ -0,0 +1,11 @@ +using Test +using PEPSKit + +@isdefined(TestSuite) || include("../testsuite/TestSuite.jl") +using .TestSuite + +is_buildkite = get(ENV, "BUILDKITE", "false") == "true" + +if !is_buildkite + TestSuite.toolbox_correlator_approx(Vector) +end diff --git a/test/toolbox/correlator_approx_phys.jl b/test/toolbox/correlator_approx_phys.jl new file mode 100644 index 000000000..135b330b0 --- /dev/null +++ b/test/toolbox/correlator_approx_phys.jl @@ -0,0 +1,11 @@ +using Test +using PEPSKit + +@isdefined(TestSuite) || include("../testsuite/TestSuite.jl") +using .TestSuite + +is_buildkite = get(ENV, "BUILDKITE", "false") == "true" + +if !is_buildkite + TestSuite.toolbox_correlator_approx_phys(Vector) +end diff --git a/test/toolbox/expval_approx.jl b/test/toolbox/expval_approx.jl new file mode 100644 index 000000000..ca8fb80ec --- /dev/null +++ b/test/toolbox/expval_approx.jl @@ -0,0 +1,12 @@ +using Test +using PEPSKit + +@isdefined(TestSuite) || include("../testsuite/TestSuite.jl") +using .TestSuite + +is_buildkite = get(ENV, "BUILDKITE", "false") == "true" + +if !is_buildkite + TestSuite.toolbox_expval_approx(Vector) + TestSuite.toolbox_expval_approx_localoperator(Vector) +end diff --git a/test/toolbox/mpo_routing.jl b/test/toolbox/mpo_routing.jl new file mode 100644 index 000000000..49c08c02b --- /dev/null +++ b/test/toolbox/mpo_routing.jl @@ -0,0 +1,13 @@ +using Test +using PEPSKit + +@isdefined(TestSuite) || include("../testsuite/TestSuite.jl") +using .TestSuite + +is_buildkite = get(ENV, "BUILDKITE", "false") == "true" + +if !is_buildkite + TestSuite.toolbox_mpo_routing_fusers(Vector) + TestSuite.toolbox_mpo_routing_identities(Vector) + TestSuite.toolbox_mpo_routing_paths(Vector) +end diff --git a/test/types/localoperator.jl b/test/types/localoperator.jl index 412d9b9a2..4a356fbbd 100644 --- a/test/types/localoperator.jl +++ b/test/types/localoperator.jl @@ -52,6 +52,89 @@ if !is_buildkite @test length(op5.terms) == 2 end + @testset "LocalOperator MPO bookkeeping" begin + d = ℂ^2 + lattice = fill(d, 2, 2) + dense = randn(Float64, d ⊗ d ← d ⊗ d) + mpo = PEPSKit.gate_to_mpo(dense; trunc = notrunc()) + sites = CartesianIndex.([(3, 4), (3, 3)]) + original_sites, original_mpo = copy(sites), copy(mpo) + operator = LocalOperator(lattice, sites => mpo) + stored_sites, stored_mpo = only(operator.terms) + + # Construction translates a copy of the coordinates and must retain the caller's MPO order. + # Mutating either input container afterwards must not alter the stored term. + @test sites == original_sites + reverse!(sites) + reverse!(mpo) + @test stored_sites == CartesianIndex.([(1, 2), (1, 1)]) + @test stored_mpo == original_mpo + + # Scaling every factor would multiply the operator by α^N; only one factor may change. + # A complex coefficient must also widen the container without mutating the original tensors. + α = 2 + 3im + snapshot = deepcopy(stored_mpo) + scaled_mpo = last(only((α * operator).terms)) + @test first(scaled_mpo) ≈ α * first(stored_mpo) + @test last(scaled_mpo) === last(stored_mpo) + @test stored_mpo == snapshot + onsite = randn(ComplexF64, d ← d) + mixed = operator + LocalOperator(lattice, ((2, 1),) => onsite) + @test scalartype(mixed) == ComplexF64 + + # Taking real/imaginary parts factorwise does not give the real/imaginary part of an MPO. + @test_throws ArgumentError real(mixed) + @test_throws ArgumentError imag(mixed) + + # Insertion and operator addition use different accumulation paths; both must reject MPO sums. + # Dense terms sharing the same sites must still accumulate normally. + inds = CartesianIndex.([(1, 1), (1, 2)]) + dense_operator = LocalOperator(lattice, inds => dense) + mpo_operator = LocalOperator(lattice, inds => original_mpo) + @test_throws ArgumentError PEPSKit.add_term!(mpo_operator, inds, dense) + @test_throws ArgumentError PEPSKit.add_term!(dense_operator, inds, original_mpo) + @test_throws ArgumentError mpo_operator + mpo_operator + @test_throws ArgumentError mpo_operator + dense_operator + @test_throws ArgumentError dense_operator + mpo_operator + @test last(only((dense_operator + dense_operator).terms)) ≈ 2 * dense + + # Drop zero factors before they reach boundary-MPS normalization, where they could cause 0/0. + @test isempty(LocalOperator(lattice, inds => [zero(first(original_mpo)), last(original_mpo)]).terms) + + # Unequal physical spaces expose accidental factor reordering during coordinate transformations. + lattice = [ℂ^2 ℂ^3 ℂ^4; ℂ^5 ℂ^6 ℂ^7] + dense = randn(Float64, ℂ^6 ⊗ ℂ^2 ← ℂ^6 ⊗ ℂ^2) + mpo = PEPSKit.gate_to_mpo(dense; trunc = notrunc()) + operator = LocalOperator(lattice, ((2, 2), (1, 1)) => mpo) + @test rotr90(rotl90(operator)) == operator + @test last(only(rotr90(operator).terms)) == mpo + sites, term = only(operator.terms) + @test repeat(operator, 2, 1).terms == Dict(sites => term, (sites .+ CartesianIndex(2, 0)) => term) + end + + @testset "MPO validation" begin + d, b1, b2 = ℂ^2, ℂ^3, ℂ^4 + lattice = fill(d, 2, 2) + first_tensor = randn(Float64, d ← d ⊗ b1) + last_tensor = randn(Float64, b1 ⊗ d ← d) + mpo = [first_tensor, last_tensor] + sites = ((1, 1), (1, 2)) + + # Reject malformed chains at construction, before any routing or contraction is attempted. + @test_throws ArgumentError LocalOperator(lattice, CartesianIndex{2}[] => AbstractTensorMap[]) + @test_throws ArgumentError LocalOperator(lattice, ((1, 1),) => mpo) + @test_throws ArgumentError LocalOperator(lattice, ((1, 1), (1, 1)) => mpo) + @test_throws ArgumentError LocalOperator(lattice, sites => reverse(mpo)) + + # Bond matching and the two physical legs are independent constraints. + wrong_bond = randn(Float64, b2 ⊗ d ← d) + wrong_input = randn(Float64, b1 ⊗ d ← ℂ^3) + wrong_output = randn(Float64, b1 ⊗ ℂ^3 ← d) + @test_throws SpaceMismatch LocalOperator(lattice, sites => [first_tensor, wrong_bond]) + @test_throws SpaceMismatch LocalOperator(lattice, sites => [first_tensor, wrong_input]) + @test_throws SpaceMismatch LocalOperator(lattice, sites => [first_tensor, wrong_output]) + end + @testset "Charge shifting" begin lattice = InfiniteSquare(1, 1) elt = ComplexF64 diff --git a/test/types/standardize_dualness.jl b/test/types/standardize_dualness.jl new file mode 100644 index 000000000..5b394e39f --- /dev/null +++ b/test/types/standardize_dualness.jl @@ -0,0 +1,39 @@ +using TensorKit +using PEPSKit +using Test + +@testset "Standardize virtual-space dualness" begin + d = U1Space(0 => 1) + D = U1Space(0 => 1) + χ = U1Space(0 => 1) + Pspaces = fill(d, 2, 2, 1) + + patterns = ( + (fill(D', 2, 2, 1), fill(D', 2, 2, 1)), + ( + reshape([D, D', D', D], 2, 2, 1), + reshape([D', D, D, D'], 2, 2, 1), + ), + ) + for (Nspaces, Espaces) in patterns + ρ = InfinitePEPO(randn, ComplexF64, Pspaces, Nspaces, Espaces) + env = CTMRGEnv(randn, ComplexF64, InfinitePartitionFunction(ρ), χ) + standardized_ρ, standardized_env = PEPSKit.standardize_dualness(ρ, env) + + for row in axes(ρ, 1), col in axes(ρ, 2) + A = standardized_ρ[row, col, 1] + edge_coordinates = ( + PEPSKit.NORTH => (row - 1, col), + PEPSKit.EAST => (row, col + 1), + PEPSKit.SOUTH => (row + 1, col), + PEPSKit.WEST => (row, col - 1), + ) + for (direction, coordinates) in edge_coordinates + V = PEPSKit.virtualspace(A, direction) + E = PEPSKit.edge(standardized_env, direction, coordinates...) + @test isdual(V) == (direction in (PEPSKit.NORTH, PEPSKit.EAST)) + @test V == space(E, 2)' + end + end + end +end