diff --git a/Project.toml b/Project.toml index 672c0ad..5349d2c 100644 --- a/Project.toml +++ b/Project.toml @@ -21,6 +21,7 @@ SimpleGraphConverter = "205b04f2-f585-4877-a239-566270b3f673" SplitApplyCombine = "03a91e81-4c3e-53e1-a0a4-9c0c8f19dd66" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" +TensorAlgebra = "68bd88dc-f39d-4e12-b2ca-f046b68fcc6a" TensorOperations = "6aa20fa7-93e2-5fca-9bc0-fbd0db3c71a2" TypeParameterAccessors = "7e5a90cf-f82e-492e-a09b-e3e26432c138" @@ -28,7 +29,7 @@ TypeParameterAccessors = "7e5a90cf-f82e-492e-a09b-e3e26432c138" Adapt = "4.3.0" Dictionaries = "0.4" Graphs = "1.8.0" -ITensorBase = "0.10.2" +ITensorBase = "0.10.9" KrylovKit = "0.10.2" LinearAlgebra = "1.11.0" MatrixAlgebraKit = "0.6.8" @@ -40,6 +41,7 @@ SimpleGraphConverter = "0.1.0" SplitApplyCombine = "1.2.3" Statistics = "1.11.1" StatsBase = "0.34.4" +TensorAlgebra = "0.17" TensorOperations = "5.2" TypeParameterAccessors = "0.3.10, 0.4" julia = "1.10" diff --git a/examples/2dIsing_dynamics_Heisenbergpicture.jl b/examples/2dIsing_dynamics_Heisenbergpicture.jl index f9e6f51..3f1801d 100644 --- a/examples/2dIsing_dynamics_Heisenbergpicture.jl +++ b/examples/2dIsing_dynamics_Heisenbergpicture.jl @@ -1,7 +1,7 @@ using TensorNetworkQuantumSimulator using Graphs: center using TensorNetworkQuantumSimulator: setindex_preserve!, noprime -using ITensors: ITensors, ITensor +using TensorNetworkQuantumSimulator: ITensors, ITensor function main() nx, ny = 4, 4 diff --git a/examples/boundarymps.jl b/examples/boundarymps.jl index 4b45545..f4ba696 100644 --- a/examples/boundarymps.jl +++ b/examples/boundarymps.jl @@ -1,4 +1,5 @@ using TensorNetworkQuantumSimulator +using TensorNetworkQuantumSimulator: ITensors using Random Random.seed!(1634) diff --git a/examples/heavyhexIsing_dynamics.jl b/examples/heavyhexIsing_dynamics.jl index e6a02c7..62b87dc 100644 --- a/examples/heavyhexIsing_dynamics.jl +++ b/examples/heavyhexIsing_dynamics.jl @@ -1,5 +1,6 @@ using TensorNetworkQuantumSimulator using Statistics +using TensorNetworkQuantumSimulator: ITensor, ITensors function main() #Define the lattice diff --git a/examples/hexagonal_heisenbergmodel_thermalstate.jl b/examples/hexagonal_heisenbergmodel_thermalstate.jl index 3a1c76d..0c9f08e 100644 --- a/examples/hexagonal_heisenbergmodel_thermalstate.jl +++ b/examples/hexagonal_heisenbergmodel_thermalstate.jl @@ -1,6 +1,6 @@ using TensorNetworkQuantumSimulator using TensorNetworkQuantumSimulator: scalar_factors_quotient, TensorNetworkQuantumSimulator, freenergy -using ITensors: ITensors, ITensor +using TensorNetworkQuantumSimulator: ITensors, ITensor function main() χ = 32 diff --git a/examples/loopcorrections.jl b/examples/loopcorrections.jl index e908f0a..e5f65fd 100644 --- a/examples/loopcorrections.jl +++ b/examples/loopcorrections.jl @@ -1,4 +1,6 @@ using TensorNetworkQuantumSimulator +using TensorNetworkQuantumSimulator: ITensors + using LinearAlgebra: norm using Random diff --git a/src/Apply/apply_gates.jl b/src/Apply/apply_gates.jl index 2225e71..0f0d1f8 100644 --- a/src/Apply/apply_gates.jl +++ b/src/Apply/apply_gates.jl @@ -111,11 +111,13 @@ function apply_gate!( if length(v⃗) == 2 v1, v2 = v⃗ e = NamedEdge(v1 => v2) - ind2 = commonind(s_values, first(updated_tensors)) - δuv = dag(copy(s_values)) - δuv = replaceind(δuv, ind2, prime(ind2)) - map_diag!(sign, δuv, δuv) - s_values = denseblocks(s_values) * denseblocks(δuv) + u = commonind(s_values, first(updated_tensors)) + v = noncommonind(s_values, first(updated_tensors)) + # The new messages are the singular values over the bond-and-prime pair. + # MatrixAlgebraKit singular values are nonnegative, so the legacy sign fix + # (`s * sign(s)` via `map_diag!`) was a no-op and is dropped; fermionic sign + # handling for this message construction is future work. + s_values = replaceind(s_values, v, prime(u)) setmessage!(ψ_bpc, e, dag(s_values)) setmessage!(ψ_bpc, reverse(e), s_values) end diff --git a/src/Apply/full_update.jl b/src/Apply/full_update.jl index e2f626c..75b81fc 100644 --- a/src/Apply/full_update.jl +++ b/src/Apply/full_update.jl @@ -1,4 +1,5 @@ using KrylovKit: linsolve +using .ITensorsITensorBaseCompat: namesetdiff """ full_update(o::ITensor, ψ::TensorNetworkState, v⃗; envs, kwargs...) @@ -37,15 +38,20 @@ function full_update( apply_kwargs..., ) if symmetrize - singular_values! = Ref{ITensor}() - Rᵥ₁, Rᵥ₂, spec = factorize_svd( - Rᵥ₁ * Rᵥ₂, - inds(Rᵥ₁); - ortho = "none", - singular_values!, - apply_kwargs..., - ) - callback(; singular_values = singular_values![], truncation_error = spec.truncerr) + M = Rᵥ₁ * Rᵥ₂ + codomain = inds(Rᵥ₁) + # Balanced SVD: split the singular values symmetrically (√S into each factor). + U, S, V = svd_trunc(M, codomain; trunc = itensor_trunc(; apply_kwargs...)) + u = only(commoninds(U, S)) + v = only(commoninds(S, V)) + sqrtS = sqrth_safe(S, (u,), (v,); atol = 0, rtol = 0) + Rᵥ₁, Rᵥ₂ = U * replaceind(sqrtS, v, prime(u)), replaceind(sqrtS, u, prime(u)) * V + # Best-effort truncation error from norms; suffers catastrophic cancellation when little is + # discarded. TODO: expose MatrixAlgebraKit's `ϵ` from `ITensorBase.svd_trunc` and use it here. + total = abs2(norm(M)) + truncation_error = iszero(total) ? zero(real(scalartype(M))) : + max(zero(real(scalartype(M))), 1 - abs2(norm(S)) / total) + callback(; singular_values = S, truncation_error) end ψᵥ₁ = Qᵥ₁ * Rᵥ₁ ψᵥ₂ = Qᵥ₂ * Rᵥ₂ @@ -115,8 +121,8 @@ function optimise_p_q( fstart = print_fidelity_loss ? fidelity(envs, p_cur, q_cur, p, q, o) : 0 - qs_ind = setdiff(inds(q_cur), collect(Iterators.flatten(inds.(vcat(envs, p_cur))))) - ps_ind = setdiff(inds(p_cur), collect(Iterators.flatten(inds.(vcat(envs, q_cur))))) + qs_ind = namesetdiff(inds(q_cur), collect(Iterators.flatten(inds.(vcat(envs, p_cur))))) + ps_ind = namesetdiff(inds(p_cur), collect(Iterators.flatten(inds.(vcat(envs, q_cur))))) function b(p::ITensor, q::ITensor, o::ITensor, envs::Vector{ITensor}, r::ITensor) ts = vcat(ITensor[p, q, o, dag(prime(r))], envs) diff --git a/src/Apply/gate_definitions.jl b/src/Apply/gate_definitions.jl index 7762de4..c348fcd 100644 --- a/src/Apply/gate_definitions.jl +++ b/src/Apply/gate_definitions.jl @@ -263,8 +263,13 @@ ITensors.op(o::OpName"xx_plus_yy", ::SiteType"Qubit"; θ::Number, β::Number) = Gate for rotation by XXYY at a given angle. """ -function ITensors.op(::OpName"Rxxyy", ::SiteType"S=1/2", s1::Index, s2::Index; θ = 1) - h = 0.5 * (op("X", s1) * op("X", s2) + op("Y", s1) * op("Y", s2)) +function ITensors.op(::OpName"Rxxyy", ::SiteType"S=1/2"; θ = 1) + # Built as one two-site matrix in the manifestly charge-conserving σ± form, + # ½(XX + YY) = σ⁺σ⁻ + σ⁻σ⁺, rather than from single-site `op("X", s)` factors: + # the gate conserves U(1) charge, but a standalone `X` does not, so the factored + # construction has no symmetric representation even though the sum does. + σp, σm = [0.0 1.0; 0.0 0.0], [0.0 0.0; 1.0 0.0] + h = kron(σp, σm) + kron(σm, σp) return exp(-im * θ * h) end ITensors.op(o::OpName"Rxxyy", ::SiteType"Qubit"; θ::Number) = @@ -275,8 +280,11 @@ ITensors.op(o::OpName"Rxxyy", ::SiteType"Qubit"; θ::Number) = Gate for rotation by XXYYZZ at a given angle. """ -function ITensors.op(::OpName"Rxxyyzz", ::SiteType"S=1/2", s1::Index, s2::Index; θ = 1) - h = 0.5 * (op("X", s1) * op("X", s2) + op("Y", s1) * op("Y", s2) + op("Z", s1) * op("Z", s2)) +function ITensors.op(::OpName"Rxxyyzz", ::SiteType"S=1/2"; θ = 1) + # One two-site matrix in the σ± form, for the same reason as `Rxxyy` above: + # ½(XX + YY + ZZ) = σ⁺σ⁻ + σ⁻σ⁺ + ½ ZZ. + σp, σm, σz = [0.0 1.0; 0.0 0.0], [0.0 0.0; 1.0 0.0], [1.0 0.0; 0.0 -1.0] + h = kron(σp, σm) + kron(σm, σp) + 0.5 * kron(σz, σz) return exp(-im * θ * h) end ITensors.op(o::OpName"Rxxyyzz", ::SiteType"Qubit"; θ::Number) = diff --git a/src/Apply/simple_update.jl b/src/Apply/simple_update.jl index 40b7644..f4a3722 100644 --- a/src/Apply/simple_update.jl +++ b/src/Apply/simple_update.jl @@ -35,8 +35,14 @@ function simple_update( envs_v2 = filter(env -> hascommoninds(env, ψ⃗[2]), envs) @assert all(ndims(env) == 2 for env in vcat(envs_v1, envs_v2)) - sqrt_inv_sqrt_envs_v1 = pseudo_sqrt_inv_sqrt.(envs_v1; cutoff = sqrt_cutoff) - sqrt_inv_sqrt_envs_v2 = pseudo_sqrt_inv_sqrt.(envs_v2; cutoff = sqrt_cutoff) + # The environments are hermitian only up to numerical noise, so project before + # the square roots (which require hermitian input). + sqrt_invsqrt = env -> sqrth_invsqrth_safe( + project_hermitian(env, (inds(env)[1],), (inds(env)[2],)), + (inds(env)[1],), (inds(env)[2],); atol = sqrt_cutoff, rtol = 0 + ) + sqrt_inv_sqrt_envs_v1 = map(sqrt_invsqrt, envs_v1) + sqrt_inv_sqrt_envs_v2 = map(sqrt_invsqrt, envs_v2) sqrt_envs_v1, inv_sqrt_envs_v1 = first.(sqrt_inv_sqrt_envs_v1), last.(sqrt_inv_sqrt_envs_v1) sqrt_envs_v2, inv_sqrt_envs_v2 = first.(sqrt_inv_sqrt_envs_v2), last.(sqrt_inv_sqrt_envs_v2) @@ -49,16 +55,22 @@ function simple_update( rᵥ₁ = commoninds(Qᵥ₁, Rᵥ₁) rᵥ₂ = commoninds(Qᵥ₂, Rᵥ₂) oR = ITensors.apply(o, Rᵥ₁ * Rᵥ₂) - singular_values! = Ref{ITensor}() - Rᵥ₁, Rᵥ₂, spec = factorize_svd( - oR, - unioninds(rᵥ₁, sᵥ₁); - ortho = "none", - singular_values!, - apply_kwargs..., - ) - err = spec.truncerr - s_values = singular_values![] + # Balanced SVD: split the singular values symmetrically (√S into each factor) so neither + # side is isometric. The bond stays on `prime(u)` (keeping `u`'s name), so once this + # function `noprime`s its result the bond becomes `u`, which the returned `s_values` (over + # `(u, v)`) still shares for `apply_gate!`'s bond-message construction. + U, S, V = svd_trunc(oR, unioninds(rᵥ₁, sᵥ₁); trunc = itensor_trunc(; apply_kwargs...)) + u = only(commoninds(U, S)) + v = only(commoninds(S, V)) + sqrtS = sqrth_safe(S, (u,), (v,); atol = 0, rtol = 0) + Rᵥ₁, Rᵥ₂ = U * replaceind(sqrtS, v, prime(u)), replaceind(sqrtS, u, prime(u)) * V + s_values = S + # Best-effort truncation error from norms (SVD preserves the Frobenius norm); suffers + # catastrophic cancellation when little is discarded. TODO: expose MatrixAlgebraKit's `ϵ` + # from `ITensorBase.svd_trunc` and use it here instead. + total = abs2(norm(oR)) + err = iszero(total) ? zero(real(scalartype(oR))) : + max(zero(real(scalartype(oR))), 1 - abs2(norm(S)) / total) Qᵥ₁ = contract([Qᵥ₁; dag.(inv_sqrt_envs_v1)]) Qᵥ₂ = contract([Qᵥ₂; dag.(inv_sqrt_envs_v2)]) updated_tensors = [Qᵥ₁ * Rᵥ₁, Qᵥ₂ * Rᵥ₂] diff --git a/src/Forms/abstractform.jl b/src/Forms/abstractform.jl index e1c9ebd..70b4c2d 100644 --- a/src/Forms/abstractform.jl +++ b/src/Forms/abstractform.jl @@ -22,7 +22,9 @@ function virtualinds(form::AbstractForm, edge::NamedEdge) end function default_message(form::AbstractForm, edge::AbstractEdge) - return adapt_like(form, denseblocks(delta(virtualinds(form, edge)))) + cod = virtualinds(ket(form), edge) + dom = dag.(bra_virtualinds(form, edge)) + return one(similar_map(ket(form)[src(edge)], cod, dom), cod, dom) end function bp_factors(form::AbstractForm, verts::Vector) diff --git a/src/Forms/bilinearform.jl b/src/Forms/bilinearform.jl index 9dfbdc1..faadbdb 100644 --- a/src/Forms/bilinearform.jl +++ b/src/Forms/bilinearform.jl @@ -14,12 +14,15 @@ Base.copy(blf::BilinearForm) = BilinearForm(copy(blf.ket), copy(blf.operator), c #Constructor, bra is taken to be in the vector space of ket so the dual is taken function BilinearForm(ket::TensorNetworkState, bra::TensorNetworkState) - dtype = datatype(ket) @assert graph(ket) == graph(bra) - bra = map_tensors(t -> dag(prime(t)), bra) sinds = siteinds(ket) verts = collect(vertices(ket)) - operator_tensors = [adapt(dtype)(reduce(*, ITensor[denseblocks(delta(sind, prime(dag(sind)))) for sind in sinds[v]])) for v in verts] + bra = TensorNetworkState(Dictionary(verts, [bra_tensor(bra, v) for v in verts])) + operator_tensors = [ + let codomain = dag.(sinds[v]), domain = dag.(prime.(sinds[v])) + one(similar_map(ket[v], codomain, domain), codomain, domain) + end for v in verts + ] operator = TensorNetworkState(Dictionary(verts, operator_tensors)) return BilinearForm(ket, operator, bra) end diff --git a/src/Forms/quadraticform.jl b/src/Forms/quadraticform.jl index 12f4b0e..5c6ab3a 100644 --- a/src/Forms/quadraticform.jl +++ b/src/Forms/quadraticform.jl @@ -6,7 +6,7 @@ end ket(qf::QuadraticForm) = qf.ket operator(qf::QuadraticForm) = qf.operator bra(qf::QuadraticForm) = prime(dag(ket(qf))) -bra_tensor(qf::QuadraticForm, v) = dag(prime(ket(qf)[v])) +bra_tensor(qf::QuadraticForm, v) = bra_tensor(ket(qf), v) bra_virtualinds(qf::QuadraticForm, edge::NamedEdge) = dag.(prime.(virtualinds(ket(qf), edge))) Base.copy(qf::QuadraticForm) = QuadraticForm(copy(qf.ket), copy(qf.operator)) diff --git a/src/ITensorsITensorBaseCompat/ITensorsITensorBaseCompat.jl b/src/ITensorsITensorBaseCompat/ITensorsITensorBaseCompat.jl index da3b3b5..f38a706 100644 --- a/src/ITensorsITensorBaseCompat/ITensorsITensorBaseCompat.jl +++ b/src/ITensorsITensorBaseCompat/ITensorsITensorBaseCompat.jl @@ -26,15 +26,14 @@ include("ops.jl") export # Index access and set algebra - inds, commoninds, commonind, uniqueinds, noncommonind, noncommoninds, unioninds, hascommoninds, + inds, commoninds, commonind, uniqueinds, noncommonind, noncommoninds, unioninds, + hascommoninds, # Index operations sim, dag, prime, noprime, replaceind, replaceinds, dim, swapind, # ITensor construction - itensor, random_itensor, scalar, delta, onehot, combiner, combinedind, + itensor, random_itensor, scalar, delta, onehot, # Factorizations - qr, svd, eigen, factorize, factorize_svd, - # Diagonal manipulation - map_diag, map_diag!, + qr, svd, svd_trunc, eigen, factorize, itensor_trunc, # Storage / element-type accessors scalartype, datatype, array, data, # Dense / quantum-number no-ops diff --git a/src/ITensorsITensorBaseCompat/itensors.jl b/src/ITensorsITensorBaseCompat/itensors.jl index d9200b0..e2f7aa0 100644 --- a/src/ITensorsITensorBaseCompat/itensors.jl +++ b/src/ITensorsITensorBaseCompat/itensors.jl @@ -3,9 +3,11 @@ # # Strategy (see ITensorDevelopmentPlans api_migration_map.md): # - Names below are thin wrappers over ITensorBase / TensorAlgebra / MatrixAlgebraKit. -# - `combiner`, the factorization return shapes, `map_diag`, the operator/SiteType -# system, and the boundary-MPS (ITensorMPS) paths are NOT wrapped here; they need -# callsite translation or upstream stack work and are tracked separately. +# - The factorization return shapes, the operator/SiteType system, and +# the boundary-MPS (ITensorMPS) paths are NOT wrapped here; they need callsite +# translation or upstream stack work and are tracked separately. The legacy +# `combiner` is retired rather than wrapped: call sites use the next-gen fusion +# primitives (`matricize` below, the fused identity `Base.one`) directly. # # ITensorBase keeps most of this API internal (unexported), so we reach for the # qualified names and re-publish the legacy spellings into the TNQS namespace. @@ -15,11 +17,11 @@ using Adapt: Adapt using ITensorBase: ITensorBase, AbstractITensor, ITensor, Index, NamedUnitRange, name, nameddims, plev, tags, unnamed using LinearAlgebra: LinearAlgebra +using TensorAlgebra.MatrixAlgebra: MatrixAlgebra +using TensorAlgebra: TensorAlgebra -# Legacy `inds(t; plev, tags)` accepted index-filtering keywords; the next-gen -# `ITensorBase.inds` takes none. Own a compat `inds` that forwards to `ITensorBase.inds` -# and applies the legacy filters — a compat-owned function (like `prime`/`dag`), not a -# pirated method on `ITensorBase.inds`. +# Legacy `inds(t; plev, tags)` took index-filtering keywords; `ITensorBase.inds` takes none. Compat +# `inds` forwards to it and applies the legacy filters (a compat-owned function, not a pirated method). function inds(t::AbstractITensor; plev = nothing, tags = nothing) is = ITensorBase.inds(t) isnothing(plev) || (is = filter(i -> ITensorBase.plev(i) == plev, is)) @@ -27,20 +29,15 @@ function inds(t::AbstractITensor; plev = nothing, tags = nothing) return is end -# Functions TNQS adds methods to for its own types (`scalartype`, `contract`, -# `inner`, `datatype`). Rather than extend the upstream generics, TNQS owns these: -# `scalartype` falls back to ITensorBase for tensors/arrays; `contract`/`inner` are -# TNQS operations with their own base methods defined in the library. The per-file -# method definitions extend these module-owned functions. +# Module-owned functions TNQS extends for its own types: `scalartype` falls back to ITensorBase for +# tensors/arrays; `contract`/`inner` are TNQS operations whose base methods live in the library. scalartype(x) = ITensorBase.scalartype(x) function inner end -# Base contraction of a list of ITensors following a (possibly nested) pairwise -# contraction `sequence` (legacy `ITensors.contract(tensors; sequence)`); leaves are -# integer indices into `tensors`. TNQS adds the network-level `contract` methods. -# The list is typed `AbstractVector` (not `AbstractVector{<:AbstractITensor}`) because -# callers build it by concatenation that can widen to `Vector{Any}` (e.g. splicing in -# an empty environment list). +# Base contraction of a list of ITensors along a (possibly nested) pairwise `sequence` (legacy +# `ITensors.contract(tensors; sequence)`); leaves are integer indices into `tensors`. Typed +# `AbstractVector` (not `{<:AbstractITensor}`) because callers concatenate in ways that widen to +# `Vector{Any}` (e.g. splicing in an empty environment list). function contract end function contract(tensors::AbstractVector; sequence = nothing) return isnothing(sequence) ? reduce(*, tensors) : _contract_sequence(tensors, sequence) @@ -55,22 +52,46 @@ _compat_inds(t::AbstractITensor) = inds(t) _compat_inds(is) = is # -# Index-set algebra. Legacy ITensors provides these as set operations over the -# indices of its arguments; in the next-gen stack the same falls out of Base set ops -# on `inds` (index identity/equality is by id). +# Small-collection set operations keyed by a transform `by` (elements compare equal when +# `by(x) == by(y)`). These scan linearly rather than building `Set`s: Base's `Set`-based ops hash, +# and hashing an `Index` over a TensorKit `GradedSpace` can fall back to iterating the space, which +# throws. The intersect/setdiff/union/symdiff forms return elements of the first argument as +# `Vector`s. # -# Matching is by index *name* (`IndexName`), not by full `Index` equality: two indices -# name the same leg when their names match. On a graded axis a shared bond appears as an -# index on one tensor and its dual (`conj`) on the other — same name, opposite arrow — so -# `==` on the full `Index` misses it while the names still match. On the dense backend -# name-matching and `==` coincide, so this is a strict generalization. The full `Index` is -# returned, not the bare name: callers take `dim`/`tags` of the result and feed it to -# `Index`/`qr`/`replaceinds`, which need the space. Results are always `Vector`s. -_indvec(x) = collect(_compat_inds(x)) -commoninds(a, b) = (nb = name.(_indvec(b)); filter(i -> name(i) ∈ nb, _indvec(a))) -uniqueinds(a, b) = (nb = name.(_indvec(b)); filter(i -> name(i) ∉ nb, _indvec(a))) -unioninds(a, b) = vcat(_indvec(a), uniqueinds(b, a)) -hascommoninds(a, b) = !isempty(commoninds(a, b)) +smallintersect(a, b; by = identity) = (kb = Iterators.map(by, b); [x for x in a if by(x) ∈ kb]) +smallsetdiff(a, b; by = identity) = (kb = Iterators.map(by, b); [x for x in a if by(x) ∉ kb]) +smallunion(a, b; by = identity) = vcat(collect(a), smallsetdiff(b, a; by)) +smallsymdiff(a, b; by = identity) = vcat(smallsetdiff(a, b; by), smallsetdiff(b, a; by)) +smallisdisjoint(a, b; by = identity) = (kb = Iterators.map(by, b); !any(x -> by(x) ∈ kb, a)) +smallissubset(a, b; by = identity) = (kb = Iterators.map(by, b); all(x -> by(x) ∈ kb, a)) +smallissetequal(a, b; by = identity) = smallissubset(a, b; by) && smallissubset(b, a; by) + +# +# Name-based index-set algebra: the small-collection ops keyed by `IndexName` (`by = name`) rather +# than full `Index` equality. On a graded axis a shared bond appears as an index on one tensor and +# its dual (`conj`) on the other — same name, opposite arrow — so full-`Index` `==` misses it while +# the names match; on the dense backend the two coincide. The `Index`-returning ops give back full +# `Index`es (callers need the space to feed `Index`/`qr`/`replaceinds`). +# +nameisequal(i, j) = name(i) == name(j) +nameintersect(a, b) = smallintersect(a, b; by = name) +namesetdiff(a, b) = smallsetdiff(a, b; by = name) +nameunion(a, b) = smallunion(a, b; by = name) +namesymdiff(a, b) = smallsymdiff(a, b; by = name) +nameisdisjoint(a, b) = smallisdisjoint(a, b; by = name) +nameissubset(a, b) = smallissubset(a, b; by = name) +nameissetequal(a, b) = smallissetequal(a, b; by = name) + +# +# Named-tensor index-set algebra (the `commoninds`/etc. surface). `_compat_inds` coerces a tensor to +# its indices and passes an index collection through unchanged, so these accept either; the `name*` +# ops take it from there (always returning `Vector`s, so a result can be fed back in, as +# `reduce(noncommoninds, tensors)` does). +# +commoninds(a, b) = nameintersect(_compat_inds(a), _compat_inds(b)) +uniqueinds(a, b) = namesetdiff(_compat_inds(a), _compat_inds(b)) +unioninds(a, b) = nameunion(_compat_inds(a), _compat_inds(b)) +hascommoninds(a, b) = !nameisdisjoint(_compat_inds(a), _compat_inds(b)) # Singular forms: the one common / one unique index. `noncommonind` is used in TNQS # as "the index of `a` not shared with `b`" (e.g. the non-contracted leg of an @@ -78,7 +99,7 @@ hascommoninds(a, b) = !isempty(commoninds(a, b)) commonind(a, b) = (cs = commoninds(a, b); isempty(cs) ? nothing : first(cs)) noncommonind(a, b) = (us = uniqueinds(a, b); isempty(us) ? nothing : first(us)) # Plural: indices not shared by both (symmetric difference). -noncommoninds(a, b) = vcat(uniqueinds(a, b), uniqueinds(b, a)) +noncommoninds(a, b) = namesymdiff(_compat_inds(a), _compat_inds(b)) # Index dimension (legacy `dim`). `dim(i)` is the length; `dim(is)` the product. dim(i::Index) = length(i) @@ -87,7 +108,7 @@ dim(is::Union{Tuple, AbstractVector}) = isempty(is) ? 1 : prod(length, is) # # Index operations. # -# A fresh index with the same length (legacy `sim`). +# A fresh index with the same length, tags, and prime level (legacy `sim`). sim(i::Index) = ITensorBase.uniquename(i) sim(is::Union{Tuple, AbstractVector{<:Index}}) = map(sim, is) @@ -124,20 +145,15 @@ replaceind(t, from::Index, to::Index) = replaceinds(t, from => to) _as_index_vec(x::Index) = [x] _as_index_vec(xs) = collect(xs) cat_inds(xs...) = reduce(vcat, map(_as_index_vec, xs)) -# Legacy `replaceinds` accepted collection arguments (`replaceinds(t, [i, k], [j, l])` -# and `replaceinds(t, [i, k] => [j, l])`); ITensorBase provides only the pair-splat form -# (`replaceinds(t, i => j, k => l)`). Own a compat `replaceinds` that handles the -# collection forms and forwards the pair-splat form. +# Legacy `replaceinds` took collection arguments (`replaceinds(t, [i,k], [j,l])`, `... => ...`); +# ITensorBase provides only the pair-splat form, so this compat handles the collection forms. # -# The base case relabels *names* via `replacedimnames`, not `ITensorBase.replaceinds`. -# The latter replaces the index *space* (routing through `axes`/`getindex`), which -# scalar-indexes a graded axis and errors; a pure name relabel is space-agnostic and works -# on every backend. Keys are stripped to their `IndexName` because `replacedimnames` -# matches on `dimnames` and silently no-ops on a raw `Index` key. +# The base case relabels *names* via `replacedimnames`, not `ITensorBase.replaceinds` (which replaces +# the index *space*, scalar-indexing a graded axis and erroring). Keys are stripped to `IndexName` +# because `replacedimnames` matches on `dimnames` and silently no-ops on a raw `Index` key. # -# NB: an `Index` is itself an `AbstractVector` (it's a `NamedUnitRange`), so the -# collection types are constrained to `AbstractVector{<:Index}` to avoid capturing a -# bare `Index` (an `AbstractVector{Int}`) and iterating over its integer range. +# NB: an `Index` is itself an `AbstractVector` (`NamedUnitRange`), so the collection types are +# constrained to `AbstractVector{<:Index}` to avoid capturing a bare `Index` and iterating its range. const _IndexColl = Union{Tuple{Vararg{Index}}, AbstractVector{<:Index}} function replaceinds(t, pairs::Pair...) return ITensorBase.replacedimnames( @@ -158,15 +174,23 @@ end # Legacy `itensor(array, inds)`: inherit the index spaces. NB: `ITensor(array, inds)` # with raw `Index` objects is intentionally NOT supported by ITensorBase (the space # is underdefined); use the indexing form, which inherits the indices' spaces, or -# `ITensor(array, name.(inds))` to take the space from the array. -itensor(array, is) = array[is...] -itensor(array, is...) = array[is...] - -# TYPE PIRACY (temporary, compat-owned — NOT an upstream candidate): adds a rank-0 -# `ITensor(x::Number)` constructor, which ITensorBase deliberately omits and does not plan to -# support. Legacy ITensors uses it as a multiplicative identity to seed a product accumulator -# (`out = ITensor(1); out *= t; ...`). Kept here for now; the accumulator call sites get rewritten -# to a different pattern later, retiring this method rather than upstreaming it. +# `ITensor(array, name.(inds))` to take the space from the array. Like the legacy +# constructor, a matching total length is accepted and reshaped to the index +# dimensions (e.g. a `d^2 × d^2` two-site gate matrix over four site legs). +function itensor(array, is...) + length(array) == prod(length, is) || + throw( + DimensionMismatch( + "array with $(length(array)) elements cannot fill indices of dimensions $(length.(is))" + ) + ) + return reshape(array, map(length, is))[is...] +end +itensor(array, is::Union{Tuple, AbstractVector}) = itensor(array, is...) + +# TYPE PIRACY (temporary, compat-owned — NOT upstream): rank-0 `ITensor(x::Number)`, which +# ITensorBase deliberately omits. Legacy ITensors uses it to seed a product accumulator +# (`out = ITensor(1); out *= t`). Retired once those call sites are rewritten. ITensorBase.ITensor(x::Number) = nameddims(fill(x), ()) # Random ITensor over the given indices (legacy `random_itensor`). @@ -178,12 +202,9 @@ random_itensor(is::Union{Tuple, AbstractVector}) = random_itensor(Float64, is) # Rank-0 scalar extraction (legacy `scalar`). scalar(t::AbstractITensor) = t[] -# Dense Kronecker delta tensor. ITensorBase deliberately omits the `delta` tensor -# type that legacy ITensors had; this dense version is vendored from -# ITensorNetworksNext's `ITensorNetworkGenerators/delta_network.jl` (a dense delta -# defined on a graph generator there), copied because TNQS doesn't depend on -# ITensorNetworksNext for this migration. A graded/sector-aware `delta` is a stack -# gap (tracked separately). +# Dense Kronecker delta tensor (legacy `delta`), vendored from ITensorNetworksNext's +# `ITensorNetworkGenerators/delta_network.jl` (TNQS doesn't depend on it for this migration). A +# graded/sector-aware `delta` is a stack gap (tracked separately). diaglength(a::AbstractArray) = minimum(size(a)) function diagstride(a::AbstractArray) s = 1 @@ -213,6 +234,32 @@ function diagonaltensor( return nameddims(diagonaltensor(diag, unnamed.(is)), name.(is)) end +# Allocate an identity-map-shaped tensor over the `codomain`/`domain` indices, following +# `prototype`'s backend/eltype (domain axes dualized in storage). Unlike `Base.one(a, codomain, +# domain)` — which needs `a` to already carry the map's legs — this lets a caller whose prototype +# lacks a leg still build a device-following identity via `one(similar_map(...), codomain, domain)`. +function similar_map(prototype::AbstractITensor, eltype::Type, codomain, domain) + raw = TensorAlgebra.similar_map( + unnamed(prototype), eltype, unnamed.(codomain), unnamed.(domain) + ) + return nameddims(raw, (name.(codomain)..., name.(domain)...)) +end +function similar_map(prototype::AbstractITensor, codomain, domain) + return similar_map(prototype, scalartype(prototype), codomain, domain) +end + +# From-scratch identity map: a dense identity embedded onto the `codomain`/`domain` index +# partition via checked `project`, so the index axes select the backend (dense, graded, +# `TensorMap`). Unlike `one(a, codomain, domain)` it needs no prototype tensor, so it is the +# right primitive when only the indices and an element type are in hand (e.g. `op("I")`). +function id(eltype::Type, codomain, domain) + m = Matrix{eltype}(LinearAlgebra.I, prod(length, codomain), prod(length, domain)) + return TensorAlgebra.project(m, Tuple(codomain), Tuple(domain)) +end + +# Dense Kronecker copy (`delta`) tensor over the index axes. Dense-only: a super-diagonal +# generally cannot be embedded while preserving a nontrivial symmetry, so graded/`TensorMap` +# callers that want an order-2 identity build it via `id`/`one` at the callsite instead. delta(eltype::Type, is::Tuple) = diagonaltensor(ones(eltype, minimum(length, is)), is) delta(eltype::Type, is::Index...) = delta(eltype, is) delta(eltype::Type, is::AbstractVector{<:Index}) = delta(eltype, Tuple(is)) @@ -220,25 +267,23 @@ delta(is::Tuple) = delta(Float64, is) delta(is::Index...) = delta(Float64, is) delta(is::AbstractVector{<:Index}) = delta(Float64, Tuple(is)) -# Trace of an ITensor over its prime pairs (legacy `tr`): contract each unprimed index -# with its prime via a `delta`, the same construction `normalize_rdm` uses. This is the -# index-paired definition — independent of storage order — rather than `sum` of the -# underlying array's storage diagonal, which is only the trace for a rank-2 tensor whose -# legs happen to be ordered to align. Accessed qualified (`ITensors.tr`) so it doesn't -# shadow `LinearAlgebra.tr`, which TNQS still calls on plain matrices. +# Trace over prime pairs (legacy `tr`): contract with the identity map pairing each unprimed index +# with its prime. The domain is built as `dag.(prime.(codomain))`, not `inds(t; plev=1)` — `plev` +# filtering does not preserve the pairing order between the plev-0 and plev-1 groups. Accessed +# qualified (`ITensors.tr`) so it doesn't shadow `LinearAlgebra.tr`, which TNQS calls on matrices. function tr(t::AbstractITensor) - out = copy(t) - for i in inds(t; plev = 0) - out *= Adapt.adapt(datatype(t), delta(i, prime(i))) - end - return scalar(out) + unprimed = inds(t; plev = 0) + codomain, domain = dag.(unprimed), dag.(prime.(unprimed)) + return scalar(t * one(similar_map(t, codomain, domain), codomain, domain)) end -# One-hot vector along `i` at position `p` (legacy `onehot(i => p)`). +# One-hot vector along `i` at position `p` (legacy `onehot(i => p)`), through `project_aux` +# (in `ops.jl`), which follows the index's backend and carries the basis vector's charge on +# a derived auxiliary leg. function onehot(eltype::Type, (i, p)::Pair{<:Index}) v = zeros(eltype, length(i)) v[p] = one(eltype) - return v[i] + return project_aux(v, i) end onehot(p::Pair{<:Index}) = onehot(Float64, p) @@ -255,7 +300,7 @@ const svd_trunc = MAK.svd_trunc # `linds` may be splatted, a single index, or one collection. function qr(a::AbstractITensor, linds...) left = cat_inds(linds...) - right = setdiff(collect(inds(a)), left) + right = namesetdiff(inds(a), left) return MAK.qr_compact(a, Tuple(left), Tuple(right)) end @@ -286,7 +331,7 @@ end # truncation kwargs (`cutoff`/`maxdim`/`mindim`) are not yet translated to MAK's # `trunc=(; ...)` spec — see api_migration_map.md; the current callsites pass # `cutoff = nothing` (full decomposition). -function _hermitian_eigh( +function _eigh( m::AbstractITensor, codomain, domain; @@ -299,11 +344,6 @@ function _hermitian_eigh( "the compat `eigen` does not yet translate the `cutoff` truncation kwarg to MatrixAlgebraKit's `trunc` spec" ) cod, dom = _astuple(codomain), _astuple(domain) - # `MAK.eigh_full` rejects a matrix that is hermitian only up to numerical noise, but - # the caller asserts `ishermitian = true`. Project onto the hermitian part first - # (the conjugate transpose swaps codomain ↔ domain), matching legacy `eigen`'s - # treat-as-hermitian behavior. - m = (m + replaceinds(conj(m), (cod..., dom...), (dom..., cod...))) / 2 D, U = MAK.eigh_full(m, cod, dom) u = only(commoninds(D, U)) # eigenvalue index shared with U t = only(uniqueinds(D, U)) # D's other (independent) index @@ -312,10 +352,10 @@ function _hermitian_eigh( end # Partitioned form `eigen(m, Linds, Rinds)` reproduces legacy ITensors' reconstruction -# `m = Vt * D * dag(V)` (with `Vt` the relabeling of `V` from `Rinds` to `Linds`), which -# is what `eigendecomp`/`pseudo_sqrt_inv_sqrt` rely on — no conjugation of `U`. +# `m = Vt * D * dag(V)` (with `Vt` the relabeling of `V` from `Rinds` to `Linds`) — +# no conjugation of `U`. function eigen(m::AbstractITensor, codomain, domain; kwargs...) - return _hermitian_eigh(m, codomain, domain; kwargs...) + return _eigh(m, codomain, domain; kwargs...) end # No-partition form `eigen(m)` matches legacy ITensors' `eigen(A)`, which auto-partitions @@ -329,24 +369,33 @@ function eigen(m::AbstractITensor; kwargs...) length(p0) == length(p1) || error( "`eigen` without an index partition expects each plev-0 index to be paired with its prime" ) - D, U = _hermitian_eigh(m, Tuple(p1), Tuple(p0); kwargs...) + D, U = _eigh(m, Tuple(p1), Tuple(p0); kwargs...) return D, conj(U) end -# Translate the legacy `cutoff`/`maxdim` truncation kwargs to a MatrixAlgebraKit -# `trunc` strategy. ITensors `cutoff` discards the smallest singular values whose -# summed squares are a fraction ≤ cutoff of the total, i.e. a relative 2-norm -# truncation error of `sqrt(cutoff)` on the singular-value vector — MAK's -# `truncerror(; rtol = sqrt(cutoff), p = 2)`. `maxdim` caps the kept rank -# (`truncrank`). When both apply, the intersection keeps the more aggressive of the -# two, matching ITensors. Returns `nothing` when neither truncates. -function _trunc_spec(cutoff, maxdim) - specs = [] - isnothing(maxdim) || maxdim ≥ typemax(Int) || push!(specs, MAK.truncrank(maxdim)) - isnothing(cutoff) || iszero(cutoff) || - push!(specs, MAK.truncerror(; rtol = sqrt(cutoff), p = 2)) - isempty(specs) && return nothing - return reduce(&, specs) +# ITensors-style truncation as a MatrixAlgebraKit strategy. `maxdim` caps the kept rank +# (`truncrank`); `cutoff` discards the smallest singular values whose summed squares are a fraction +# ≤ cutoff of the total, i.e. a relative 2-norm error of `sqrt(cutoff)` on the singular-value vector +# (`truncerror(; rtol = sqrt(cutoff), p = 2)`). Either kwarg may be `nothing` (or trivial) to drop +# that constraint; with neither, the result is `notrunc()`. +function itensor_trunc(; maxdim = nothing, cutoff = nothing) + trunc = MAK.notrunc() + isnothing(maxdim) || (trunc &= MAK.truncrank(maxdim)) + isnothing(cutoff) || iszero(cutoff) || (trunc &= MAK.truncerror(; rtol = sqrt(cutoff), p = 2)) + return trunc +end + +# Split `a`'s indices into (left, right) by the requested `linds`, matching by name: +# the caller's `Index` objects may carry the other endpoint's (dual) space on a graded +# backend, so `Index`-equality filtering silently drops them (see the name-matching +# note on the index-set algebra above). Both groups are `a`'s own indices, so the +# spaces handed to the factorization are `a`'s. +function _bipartition_inds(a::AbstractITensor, linds) + lnames = name.(cat_inds(linds...)) + allinds = collect(inds(a)) + left = filter(i -> name(i) ∈ lnames, allinds) + right = filter(i -> name(i) ∉ lnames, allinds) + return left, right end # Legacy `factorize(a, linds...; ortho, cutoff, maxdim, tags)` splits `a` into `L * R` @@ -363,11 +412,9 @@ function factorize( maxdim = nothing, tags = nothing ) - allinds = collect(inds(a)) - left = filter(∈(allinds), cat_inds(linds...)) - right = setdiff(allinds, left) - trunc = _trunc_spec(cutoff, maxdim) - if isnothing(trunc) + left, right = _bipartition_inds(a, linds) + notruncation = (isnothing(cutoff) || iszero(cutoff)) && isnothing(maxdim) + if notruncation if ortho == "left" L, R = MAK.qr_compact(a, Tuple(left), Tuple(right)) elseif ortho == "right" @@ -378,7 +425,7 @@ function factorize( ) end else - U, S, Vt = MAK.svd_trunc(a, Tuple(left), Tuple(right); trunc) + U, S, Vt = MAK.svd_trunc(a, Tuple(left), Tuple(right); trunc = itensor_trunc(; cutoff, maxdim)) if ortho == "left" L, R = U, S * Vt elseif ortho == "right" @@ -397,141 +444,38 @@ function factorize( return L, R end -# Absorb an SVD `(U, S, Vt)` into two factors per `ortho`: `"left"` makes `U` isometric, -# `"right"` makes `Vt` isometric, `"none"` splits `S = √S · √S` so the weight is shared. -# For `"none"` the bond lands on `prime(u)` (where `u`/`v` are the SVD's left/right -# indices) and the separately-returned `singular_values!` stays over the unprimed -# `(u, v)` — matching legacy ITensors so a caller that `noprime`s the factors (e.g. -# `simple_update`) ends up with the bond on `u`, which `S` still shares. -function _absorb_svd(U, S, Vt, ortho) - if ortho == "left" - return U, S * Vt - elseif ortho == "right" - return U * S, Vt - elseif ortho == "none" - u = only(commoninds(U, S)) - v = only(commoninds(S, Vt)) - up = prime(u) - sqrtσ = sqrt.(abs.(diagview(unnamed(S)))) - return U * diagonaltensor(sqrtσ, (u, up)), Vt * diagonaltensor(sqrtσ, (v, up)) - end - return error( - "compat `factorize_svd` supports ortho = \"left\" / \"right\" / \"none\" (got $(repr(ortho)))" - ) -end -# Legacy `factorize_svd(a, linds...; ortho, singular_values!, cutoff, maxdim, tags)`: -# an always-SVD factorization returning `(L, R, spec)`. `spec.truncerr` is the fraction -# of squared spectral weight discarded (SVD preserves the Frobenius norm, so the total -# weight is `norm(a)^2`). If `singular_values!` is a `Ref`, it is filled with `S`. -function factorize_svd( - a::AbstractITensor, - linds...; - ortho = "left", - singular_values! = nothing, - cutoff = nothing, - maxdim = nothing, - tags = nothing - ) - allinds = collect(inds(a)) - left = filter(∈(allinds), cat_inds(linds...)) - right = setdiff(allinds, left) - trunc = _trunc_spec(cutoff, maxdim) - U, S, Vt = if isnothing(trunc) - MAK.svd_compact(a, Tuple(left), Tuple(right)) - else - MAK.svd_trunc(a, Tuple(left), Tuple(right); trunc) - end - total = abs2(LinearAlgebra.norm(a)) - kept = sum(abs2, diagview(unnamed(S))) - truncerr = if iszero(total) - zero(real(scalartype(a))) - else - max(zero(kept / total), 1 - kept / total) - end - isnothing(singular_values!) || (singular_values![] = S) - L, R = _absorb_svd(U, S, Vt, ortho) - if !isnothing(tags) - b = only(commoninds(L, R)) - bnew = settags(b, tags) - L, R = replaceind(L, b, bnew), replaceind(R, b, bnew) - end - return L, R, (; truncerr) -end - -# -# Combiner / index fusion. Legacy `combiner(is...)` returns a tensor that, when -# contracted, fuses `is` into one index; `dag(C)` splits it back, and -# `combinedind(C)` recovers that fused index. Reproduced densely as a reshaped -# identity (`C[is..., c] = δ(c, flatten(is))`), which the next-gen fusion -# primitives back via ordinary contraction. The fused index carries a sentinel tag -# so `combinedind` can find it. Dense only; a graded combiner is a stack gap. -const _COMBINER_TAG = "combiner" -function combiner(is::Tuple; eltype::Type = Float64) - d = prod(length, is) - c = ITensorBase.settag(Index(d), _COMBINER_TAG, "combined") - return reshape(Matrix{eltype}(LinearAlgebra.I, d, d), (length.(is)..., d))[is..., c] -end -combiner(is::Index...; kwargs...) = combiner(is; kwargs...) -combiner(is::AbstractVector{<:Index}; kwargs...) = combiner(Tuple(is); kwargs...) -function combinedind(c::AbstractITensor) - return only( - filter(i -> get(ITensorBase.tags(i), _COMBINER_TAG, "") == "combined", inds(c)) - ) -end # -# Diagonal manipulation. Legacy `map_diag(f, T)` applies `f` to the diagonal of a -# (diagonal-like) tensor; used on factorization spectra (singular values / -# eigenvalues). Reproduced via the same diagonal machinery as `delta`. -_diagcartesian(arr, k) = CartesianIndex(ntuple(Returns(k), ndims(arr))) -function map_diag(f, T::AbstractITensor) - arr = copy(unnamed(T)) - for k in 1:minimum(size(arr)) - idx = _diagcartesian(arr, k) - arr[idx] = f(arr[idx]) - end - return arr[inds(T)...] -end -function map_diag!(f, T::AbstractITensor) - arr = unnamed(T) - for k in 1:minimum(size(arr)) - idx = _diagcartesian(arr, k) - arr[idx] = f(arr[idx]) - end - return T -end -# Out-of-place-into-`dest` form `map_diag!(f, dest, src)`: write `f` of `src`'s diagonal -# onto `dest`'s diagonal (TNQS calls it with `dest === src` for an in-place diagonal map). -function map_diag!(f, dest::AbstractITensor, src::AbstractITensor) - d, s = unnamed(dest), unnamed(src) - for k in 1:minimum(size(s)) - d[_diagcartesian(d, k)] = f(s[_diagcartesian(s, k)]) - end - return dest -end +# Index fusion. The legacy `combiner` is retired (no compat shim); call sites fuse +# index groups with the next-gen `matricize(t, row_inds => row_name, col_inds => +# col_name)` (minting each fused name via `uniquename(IndexName)`) or build a fused +# identity with `Base.one`, both of which are graded-capable. `matricize` is the +# `TensorAlgebra` generic, extended by ITensorBase for named tensors. +using TensorAlgebra: matricize # # Storage / element type accessors. # # `scalartype` is re-exported above. `datatype` is the underlying storage array type -# (used by `adapt`); `array` / `data` expose the plain unnamed array. +# (used by `adapt`); `data` exposes the plain unnamed array. `array` densifies (legacy +# `array` materialized a dense array from any storage): a no-op on a dense backend, +# while graded / `TensorMap` storage converts through its canonical flat basis, so +# positions agree with `onehot` / `project` on the same axes. datatype(T::AbstractITensor) = typeof(unnamed(T)) -array(T::AbstractITensor) = unnamed(T) +array(T::AbstractITensor) = convert(Array, unnamed(T)) data(T::AbstractITensor) = unnamed(T) -# TYPE PIRACY (temporary, compat-owned — NOT an upstream candidate): extends -# `Adapt.adapt_structure` for `AbstractITensor` with an eltype target. Using -# `Adapt.adapt_structure` for eltype *conversion* is an abuse of Adapt.jl (Adapt is for -# storage/device adaptation, not changing the scalar type), so this does not belong upstream. -# Kept here for now; the eltype-conversion call sites get rewritten with a different pattern -# later, retiring this shim rather than upstreaming it. -# -# Legacy `adapt(eltype)(t)` converts an ITensor's scalar (element) type. ITensorBase's -# Adapt integration adapts the storage array/device type but leaves the element type -# alone, so reproduce the eltype conversion for a `Number` target (TNQS uses -# `adapt(eltype)(state(...))` to build typed product states). +# TYPE PIRACY (temporary, compat-owned — NOT upstream): extends `Adapt.adapt_structure` for +# `AbstractITensor` with a `Number` eltype target, reproducing legacy `adapt(eltype)(t)` scalar-type +# conversion (ITensorBase's Adapt integration adapts storage/device but leaves the eltype alone). +# Using Adapt for eltype conversion is an abuse, so it does not belong upstream; retired once the +# call sites move to a different pattern. function Adapt.adapt_structure(::Type{elt}, T::AbstractITensor) where {elt <: Number} + # Short-circuit the no-op case: a non-`AbstractArray` backend (e.g. a `TensorMap`) + # has no `convert(AbstractArray{elt}, ...)` method, but needs none when the element + # type already matches. + eltype(T) === elt && return T return nameddims(convert(AbstractArray{elt}, unnamed(T)), ITensorBase.dimnames(T)) end @@ -554,10 +498,8 @@ hasqns(::Any) = false # The operator / named-state system (`op` / `state`) is vendored separately in # `ops.jl`, included right after this file by the module file. -# Direct sum (legacy `directsum`): block-diagonal placement of several tensors along -# specified axes, with the non-summed (shared) axes preserved. Vendored densely -# because the next-gen stack has no `directsum` yet — a real missing ITensorBase -# feature (tracked); drop this once it lands upstream. +# Direct sum (legacy `directsum`): block-diagonal placement of tensors along the summed axes, shared +# axes preserved. Vendored densely — the next-gen stack has no `directsum` yet (tracked upstream). # directsum(out_inds, (t1 => summed_inds1), (t2 => summed_inds2), ...) function directsum(out_inds, pairs::Pair...) out_inds = Tuple(out_inds) @@ -630,13 +572,11 @@ function settags(i::Index, d::AbstractDict) end return i end -# TYPE PIRACY (temporary, compat-owned — NOT an upstream candidate): the two methods below add -# `ITensorBase.Index` constructors taking a tag string / tag dict, a legacy positional-tag form -# ITensorBase does not plan to support (the next-gen spelling passes tags via the `tags` keyword -# argument). Kept here for now; the call sites get modernized to the `tags` kwarg later, retiring -# these methods rather than upstreaming them. -# -# Legacy positional tagged-index constructor `Index(dim, "tag")`. +# TYPE PIRACY (temporary, compat-owned — NOT upstream): the two methods below add `Index` +# constructors taking a tag string / tag dict (legacy positional-tag form; next-gen passes tags via +# the `tags` kwarg). Retired once the call sites move to the kwarg form. Only the `dim` first-argument +# form is pirated: a range/space first argument collides with ITensorBase's own two-argument `Index` +# constructors, so tagging an axis-minted index spells the two steps (`settags(Index(r), "tag")`). ITensorBase.Index(dim::Integer, tagstr::AbstractString) = settags(Index(dim), tagstr) # Build a fresh index carrying a tag dictionary (legacy `Index(dim, tags(i))`, where # the next-gen `tags` returns a `Dict{String, String}`). diff --git a/src/ITensorsITensorBaseCompat/ops.jl b/src/ITensorsITensorBaseCompat/ops.jl index 079564c..0a17db5 100644 --- a/src/ITensorsITensorBaseCompat/ops.jl +++ b/src/ITensorsITensorBaseCompat/ops.jl @@ -1,14 +1,11 @@ -# Vendored minimal operator / named-state system (legacy `ITensors.op` / `ITensors.state`). +# Vendored minimal operator / named-state system (legacy `ITensors.op` / `ITensors.state`), scoped to +# the `S=1/2` single- and two-qubit gate set and the computational/Pauli-basis product states. # -# Round-1 scope is the single-qubit (`S=1/2`) gate set and the computational/Pauli-basis -# product states — enough for `siteinds` / `tensornetworkstate` state construction and -# single-site `expect`. The full two-qubit gate set, the parametric and in-house gates, -# and the Heisenberg/PTM (PauliPropagation) path belong to the gate-application round -# (`Apply/`), which is still excluded; they are not vendored here. -# -# Legacy ITensors exposes these through the `OpName` / `SiteType` dispatch system. TNQS -# only ever calls `op(name, sites...; kwargs...)` and `state(name, site)`, so the vendor -# is a plain name-keyed lookup rather than a reimplementation of that dispatch machinery. +# Legacy ITensors exposes these through the `OpName` / `SiteType` dispatch system. TNQS only ever +# calls `op(name, sites...; kwargs...)` and `state(name, site)`, so this vendors a name-keyed lookup +# over that machinery rather than reimplementing it. + +using TensorAlgebra: TensorAlgebra, project, tryproject # # Named single-site states. `state(name, i)` returns the named state vector as an @@ -20,16 +17,38 @@ const _STATE_VECTORS = Dict{String, Vector{ComplexF64}}( "X+" => [1, 1] / sqrt(2), "+" => [1, 1] / sqrt(2), "X-" => [1, -1] / sqrt(2), "-" => [1, -1] / sqrt(2), "Y+" => [1, im] / sqrt(2), - "Y-" => [1, -im] / sqrt(2), + "Y-" => [1, -im] / sqrt(2) ) function state(name::AbstractString, i::Index) v = get(_STATE_VECTORS, name) do - error("unknown single-site state \"$name\" (vendored states: $(sort(collect(keys(_STATE_VECTORS)))))") + return error( + "unknown single-site state \"$name\" (vendored states: $(sort(collect(keys(_STATE_VECTORS)))))" + ) end + return state(v, i) +end +# Project a raw vector as a state `ITensor` over `i`, adding an auxiliary leg only when the +# vector cannot live in the flux-zero space over `i` alone. The index axis selects the +# backend (dense, graded, `TensorMap`). Shared by the state constructors and `onehot`. +function project_aux(v::AbstractVector{<:Number}, i::Index) length(v) == length(i) || - error("state \"$name\" has dimension $(length(v)) but the site index has dimension $(length(i))") - return collect(v)[i] + error( + "state vector has dimension $(length(v)) but the site index has dimension $(length(i))" + ) + ψ = tryproject(v, (i,)) + isnothing(ψ) || return ψ + # The vector carries a charge under `i`'s grading (e.g. "Dn" on a U(1) site, or a + # one-hot on a graded link), so it can't live in the flux-zero space over `i` alone. + # Carry the charge on an explicit length-1 auxiliary leg: project with a trailing axis + # so the backend derives the leg's sector from the vector, then wrap the derived axis + # in a freshly named `Index`. + raw = project(reshape(v, (length(v), 1)), (unnamed(i),), ()) + aux = Index(TensorAlgebra.axes(raw, 2)) + return nameddims(raw, (ITensorBase.name(i), ITensorBase.name(aux))) end +# Vector form (legacy `ITensor(v, i)` for a state vector): the state vector as an `ITensor` +# over `i`. +state(v::AbstractVector{<:Number}, i::Index) = project_aux(v, i) # # Operators. Legacy ITensors exposes operators through the `OpName` / `SiteType` @@ -53,13 +72,15 @@ macro SiteType_str(name) end # Embed a `d^n × d^n` operator matrix (computational basis, first site most -# significant) into an `ITensor` over `(prime.(sites)..., sites...)`. +# significant) into an `ITensor` with codomain `prime.(sites)` (outputs) and +# domain `sites` (inputs). `project` is checked, so on graded sites an operator +# that is not symmetric under the site index's grading throws an `InexactError`. function _op_matrix_to_itensor(M::AbstractMatrix, sites::Tuple) ds = length.(sites) n = length(sites) A = reshape(Matrix{ComplexF64}(M), (reverse(ds)..., reverse(ds)...)) A = permutedims(A, (reverse(1:n)..., reverse((n + 1):(2n))...)) - return A[ITensorBase.prime.(sites)..., sites...] + return project(A, ITensorBase.prime.(sites), sites) end # Top-level `op(name, sites...; kwargs...)`. The identity is dimension-general (used @@ -68,7 +89,7 @@ end function op(name::AbstractString, sites::Index...; kwargs...) if name in ("I", "Id") && length(sites) == 1 s = only(sites) - return Matrix{ComplexF64}(LinearAlgebra.I, length(s), length(s))[ITensorBase.prime(s), s] + return id(ComplexF64, (ITensorBase.prime(s),), (s,)) end return op(OpName(name), SiteType("S=1/2"), sites...; kwargs...) end @@ -87,15 +108,17 @@ const _GATE_MATRICES = Dict{String, Matrix{ComplexF64}}( "Z" => [1 0; 0 -1], "H" => [1 1; 1 -1] / sqrt(2), "S" => [1 0; 0 im], - "T" => [1 0; 0 cis(π / 4)], + "T" => [1 0; 0 cis(π / 4)] ) function _gate_matrix(name::AbstractString; kwargs...) name in ("Rx", "Ry", "Rz", "P") || return get(_GATE_MATRICES, name) do - error("unknown single-site operator \"$name\" (vendored operators: $(sort(collect(keys(_GATE_MATRICES)))) plus Rx/Ry/Rz/P)") + return error( + "unknown single-site operator \"$name\" (vendored operators: $(sort(collect(keys(_GATE_MATRICES)))) plus Rx/Ry/Rz/P)" + ) end if name == "Rx" θ = kwargs[:θ] - return ComplexF64[cos(θ / 2) (-im*sin(θ / 2)); (-im*sin(θ / 2)) cos(θ / 2)] + return ComplexF64[cos(θ / 2) (-im * sin(θ / 2)); (-im * sin(θ / 2)) cos(θ / 2)] elseif name == "Ry" θ = kwargs[:θ] return ComplexF64[cos(θ / 2) (-sin(θ / 2)); sin(θ / 2) cos(θ / 2)] @@ -108,7 +131,9 @@ function _gate_matrix(name::AbstractString; kwargs...) end end # Matrix method for any single-qubit name handled by `_gate_matrix`. -op(::OpName{Name}, ::SiteType"S=1/2"; kwargs...) where {Name} = _gate_matrix(String(Name); kwargs...) +function op(::OpName{Name}, ::SiteType"S=1/2"; kwargs...) where {Name} + return _gate_matrix(String(Name); kwargs...) +end # Two-qubit gate matrices (computational basis, first site most significant). const _σx = ComplexF64[0 1; 1 0] @@ -117,25 +142,20 @@ const _σz = ComplexF64[1 0; 0 -1] op(::OpName"Rxx", ::SiteType"S=1/2"; ϕ) = exp(-im * ϕ * kron(_σx, _σx)) op(::OpName"Ryy", ::SiteType"S=1/2"; ϕ) = exp(-im * ϕ * kron(_σy, _σy)) op(::OpName"Rzz", ::SiteType"S=1/2"; ϕ) = exp(-im * ϕ * kron(_σz, _σz)) -op(::OpName"CPHASE", ::SiteType"S=1/2"; ϕ) = ComplexF64[1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 cis(ϕ)] +function op(::OpName"CPHASE", ::SiteType"S=1/2"; ϕ) + return ComplexF64[1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 cis(ϕ)] +end +function op(::OpName"SWAP", ::SiteType"S=1/2") + return Float64[1 0 0 0; 0 0 1 0; 0 1 0 0; 0 0 0 1] +end -# TYPE PIRACY (temporary): extends `Base.exp` for an operator `ITensor` (matricize over the -# index/prime pairs, exponentiate, rebuild). Legacy ITensors provided `exp(::ITensor)`; gates -# defined as `exp` of a Hamiltonian operator rely on it (e.g. a user -# `op(::OpName"MyZRot", ...) = exp(-im θ/2 * op("Z", s))`). To de-pirate: make this compat-owned -# (an `exp` in this module, not a `Base.exp` method), inferring the prime-pair codomain/domain and -# forwarding to ITensorBase's matricization `exp(a, dimnames_codomain, dimnames_domain)` — which -# already exists and is graded-capable, so this dense combiner-based version goes away. Not an -# upstream candidate (the upstream matricization `exp` is the target, not a `Base.exp(::ITensor)`). +# TYPE PIRACY (temporary): extends `Base.exp` for an operator `ITensor`, inferring the prime-pair +# codomain/domain and forwarding to ITensorBase's matricization `exp(a, codomain, domain)` +# (graded-capable). Gates defined as `exp` of a Hamiltonian operator rely on it (e.g. a user +# `op(::OpName"MyZRot", ...) = exp(-im θ/2 * op("Z", s))`). De-pirate by making it compat-owned. function Base.exp(t::AbstractITensor) p0 = filter(i -> ITensorBase.plev(i) == 0, collect(inds(t))) isempty(p0) && error("exp(::ITensor) expects indices paired as (i, prime(i))") p1 = map(ITensorBase.prime, p0) - cr, cc = combiner(Tuple(p1)), combiner(Tuple(p0)) - tc = (t * cr) * cc - rci, cci = combinedind(cr), combinedind(cc) - d = length(rci) - M = ComplexF64[tc[rci => a, cci => b] for a in 1:d, b in 1:d] - E = exp(M)[rci, cci] - return (E * dag(cr)) * dag(cc) + return exp(t, Tuple(p1), Tuple(p0)) end diff --git a/src/MessagePassing/beliefpropagationcache.jl b/src/MessagePassing/beliefpropagationcache.jl index f864504..224d993 100644 --- a/src/MessagePassing/beliefpropagationcache.jl +++ b/src/MessagePassing/beliefpropagationcache.jl @@ -170,8 +170,9 @@ function loop_correlation(bpc::BeliefPropagationCache, loop::Vector{<:NamedEdge} seq = contraction_sequence(tensors; alg = "omeinsum", optimizer = GreedyMethod()) t = contract(tensors; sequence = seq) - row_combiner, col_combiner = ITensors.combiner(e_virtualinds), ITensors.combiner(e_virtualinds_sim) - t = t * row_combiner * col_combiner + row_name = ITensorBase.uniquename(ITensorBase.IndexName) + col_name = ITensorBase.uniquename(ITensorBase.IndexName) + t = ITensors.matricize(t, Tuple(e_virtualinds) => row_name, Tuple(e_virtualinds_sim) => col_name) t = adapt(Vector{ComplexF64})(t) t = ITensors.array(t) λs = reverse(sort(LinearAlgebra.eigvals(t); by = abs)) diff --git a/src/MessagePassing/boundarympscache.jl b/src/MessagePassing/boundarympscache.jl index 47e09c8..fcef5d3 100644 --- a/src/MessagePassing/boundarympscache.jl +++ b/src/MessagePassing/boundarympscache.jl @@ -191,11 +191,16 @@ function set_interpartition_messages!( end for i in 1:(length(es) - 1) virt_dim = virtual_index_dimension(bmps_cache, es[i], es[i + 1]) - ind = Index(virt_dim, "m$(i)$(i + 1)") m1, m2 = message(bmps_cache, es[i]), message(bmps_cache, es[i + 1]) - t = adapt_like(m1, dense(delta(ind))) + # The stitching leg is minted trivial (all weight in the charge-0 sector) + # so it follows the messages' backend; the all-ones filling matches the + # legacy dense `delta(ind)`. + ind = ITensors.settags(Index(trivialrange(unnamed(first(inds(m1))), virt_dim)), "m$(i)$(i + 1)") + t = fill!(similar(m1, (ind,)), true) + # The two copies of the stitching leg contract against each other along the + # message MPS, so one side takes the conjugate. setmessage!(bmps_cache, es[i], m1 * t) - setmessage!(bmps_cache, es[i + 1], m2 * t) + setmessage!(bmps_cache, es[i + 1], m2 * dag(t)) end end return bmps_cache @@ -403,7 +408,7 @@ function generic_apply( mps::Dictionary{Int, <:ITensor}, right_inds::Vector{<:Vector{<:Index}}; cutoff = 0.0, - maxdim = typemax(Int), + maxdim = nothing, normalize = true, ) b = length(mpo) diff --git a/src/MessagePassing/loopcorrection.jl b/src/MessagePassing/loopcorrection.jl index cc5298c..dd30010 100644 --- a/src/MessagePassing/loopcorrection.jl +++ b/src/MessagePassing/loopcorrection.jl @@ -19,6 +19,14 @@ end function sim_edgeinduced_subgraph(bpc::BeliefPropagationCache, eg) bpc = copy(bpc) vs = collect(vertices(eg)) + # Aux (dangling non-physical) legs recorded before the bond relabeling below: the + # relabeled bond dangles in the modified network, and the live classification in + # `norm_factors` would misread it as a charge leg and unprime it on the bra. + aux = if network(bpc) isa TensorNetworkState + Dictionary(vs, [auxinds(network(bpc), v) for v in vs]) + else + nothing + end es = unique(collect(Iterators.flatten(boundary_edges(bpc, [v]; dir = :out) for v in vs))) updated_es = NamedEdge[] @@ -35,31 +43,37 @@ function sim_edgeinduced_subgraph(bpc::BeliefPropagationCache, eg) ms = messages(bpc) set!(ms, reverse(e), mer) t = network(bpc)[src(e)] - t_inds = filter(i -> i ∈ linds, inds(t)) + # Match by name: on a graded backend the network tensor carries the dual of + # the message's axis, so a full `Index` comparison never matches. + t_inds = commoninds(t, linds) if !isempty(t_inds) t_ind = only(t_inds) - t_ind_pos = findfirst(x -> x == t_ind, linds) + t_ind_pos = findfirst(x -> name(x) == name(t_ind), linds) t = replaceind(t, t_ind, linds_sim[t_ind_pos]) setindex_preserve!(bpc, t, src(e)) end push!(updated_es, e) if e ∈ edges(eg) || reverse(e) ∈ edges(eg) - row_inds, col_inds = linds, linds_sim + me = message(bpc, e) + # The identity legs take their axes from the actual message legs (`me` + # for the rows, the relabeled `mer` for the columns) so `ap - me * mer` + # lines up on every backend: the two messages carry mutually dual copies + # of the bond axes. The domain of the fused identity comes out dualized + # relative to the passed indices, so the columns go in `dag`ed. + row_inds = Index[only(commoninds(me, [l])) for l in linds] + col_inds = Index[only(commoninds(mer, [l])) for l in linds_sim] if network(bpc) isa TensorNetworkState - row_inds = vcat(row_inds, dag.(prime.(row_inds))) - col_inds = vcat(col_inds, dag.(prime.(col_inds))) + append!(row_inds, Index[only(commoninds(me, [prime(l)])) for l in linds]) + append!(col_inds, Index[only(commoninds(mer, [prime(l)])) for l in linds_sim]) end - row_combiner, col_combiner = combiner(row_inds), combiner(col_inds) - ap = - adapt_like(message(bpc, e), denseblocks(delta(combinedind(col_combiner), dag(combinedind(row_combiner))))) - ap = ap * row_combiner * dag(col_combiner) - ap = ap - message(bpc, e) * mer + ap = adapt_like(me, identity_tensor(row_inds, dag.(col_inds))) + ap = ap - me * mer push!(antiprojectors, ap) end end end - return bpc, antiprojectors + return bpc, antiprojectors, aux end #Get the all edges incident to the region specified by the vector of edges passed @@ -81,10 +95,14 @@ end function weight(bpc::BeliefPropagationCache, eg) vs = collect(vertices(eg)) es = collect(edges(eg)) - bpc, antiprojectors = sim_edgeinduced_subgraph(bpc, eg) + bpc, antiprojectors, aux = sim_edgeinduced_subgraph(bpc, eg) incoming_ms = ITensor[message(bpc, e) for e in boundary_edges(bpc, es)] - local_tensors = collect(Iterators.flatten(bp_factors(bpc, v) for v in vs)) + local_tensors = if isnothing(aux) + collect(Iterators.flatten(bp_factors(bpc, v) for v in vs)) + else + collect(Iterators.flatten(norm_factors(network(bpc), [v]; auxinds_f = u -> aux[u]) for v in vs)) + end ts = [incoming_ms; local_tensors; antiprojectors] seq = any(hasqns.(ts)) ? contraction_sequence(ts; alg = "optimal") : contraction_sequence(ts; alg = "omeinsum", optimizer = GreedyMethod()) return scalar(contract(ts; sequence = seq)) diff --git a/src/TensorNetworks/abstracttensornetwork.jl b/src/TensorNetworks/abstracttensornetwork.jl index 0203abb..444f66f 100644 --- a/src/TensorNetworks/abstracttensornetwork.jl +++ b/src/TensorNetworks/abstracttensornetwork.jl @@ -28,13 +28,15 @@ function maxvirtualdim(tn::AbstractTensorNetwork) return maximum(maximum.([dim.(virtualinds(tn, e)) for e in edges(tn)])) end +# Compare by name, not by `Index` equality: a shared graded link is stored nondual +# on one endpoint and dual (conjugated) on the other, so the two `Index` objects +# differ even though they name the same bond. function ITensors.uniqueinds(tn::AbstractTensorNetwork, v) tv_inds = Index[i for i in inds(tn[v])] vns = neighbors(tn, v) isempty(vns) && return tv_inds - neighbor_inds = reduce(vcat, [Index[i for i in inds(tn[vn])] for vn in vns]) - is = setdiff(tv_inds, neighbor_inds) - return is + neighbor_names = reduce(vcat, [[name(i) for i in inds(tn[vn])] for vn in vns]) + return filter(i -> name(i) ∉ neighbor_names, tv_inds) end function setindex_preserve!(tn::AbstractTensorNetwork, value::ITensor, vertex) @@ -105,24 +107,6 @@ function map_virtualinds(f::Function, tn::AbstractTensorNetwork) return map_virtualinds!(f, tn) end -function combine_virtualinds!(tn::AbstractTensorNetwork) - dtype = datatype(tn) - for e in edges(tn) - vinds = ITensors.commoninds(tn[src(e)], tn[dst(e)]) - if length(vinds) > 1 - C = adapt(dtype)(ITensors.combiner(vinds)) - setindex_preserve!(tn, tn[src(e)] * C, src(e)) - setindex_preserve!(tn, tn[dst(e)] * C, dst(e)) - end - end - return tn -end - -function combine_virtualinds(tn::AbstractTensorNetwork) - tn = copy(tn) - return combine_virtualinds!(tn) -end - """Add two tensornetworks together. The network structures need to be have the same graph structure""" function add(tn1::AbstractTensorNetwork, tn2::AbstractTensorNetwork) @assert graph(tn1) == graph(tn2) diff --git a/src/TensorNetworks/tensornetworkstate.jl b/src/TensorNetworks/tensornetworkstate.jl index b4cbc44..01ee66b 100644 --- a/src/TensorNetworks/tensornetworkstate.jl +++ b/src/TensorNetworks/tensornetworkstate.jl @@ -1,4 +1,4 @@ -using .ITensorsITensorBaseCompat: random_itensor +using .ITensorsITensorBaseCompat: nameisdisjoint, namesetdiff, random_itensor """ TensorNetworkState{V} <: AbstractTensorNetwork{V} @@ -22,6 +22,14 @@ tensors(tns::TensorNetworkState) = tensors(tensornetwork(tns)) Base.copy(tns::TensorNetworkState) = TensorNetworkState(copy(tensornetwork(tns)), copy(siteinds(tns))) TensorNetworkState(tn::TensorNetwork) = TensorNetworkState(tn, siteinds(tn)) +# Widen a concretely-typed site-index dictionary to the field type. +function TensorNetworkState(tn::TensorNetwork{V}, sinds::Dictionary) where {V} + s = Dictionary{V, Vector{<:Index}}() + for v in keys(sinds) + set!(s, v, sinds[v]) + end + return TensorNetworkState(tn, s) +end TensorNetworkState(tensors::Dictionary, g::NamedGraph) = TensorNetworkState(TensorNetwork(tensors, g)) TensorNetworkState(tensors::Union{Dictionary, Vector{<:ITensor}}) = TensorNetworkState(TensorNetwork(tensors)) @@ -47,12 +55,27 @@ function Base.setindex!(tns::TensorNetworkState, value::ITensor, v) return tns end -function norm_factors(tns::TensorNetworkState, verts::Vector; op_strings::Function = v -> "I") +# Bra copy of a tensor: `dag` and prime all legs except `auxinds` — dangling +# non-physical legs (e.g. a charged state's charge leg), which always pair directly +# between ket and bra rather than through an operator or a message. +function bra_tensor(t::ITensor, auxinds::Vector{<:Index}) + tdag = dag(prime(t)) + return isempty(auxinds) ? tdag : replaceinds(tdag, prime.(auxinds), auxinds) +end +bra_tensor(tns::TensorNetworkState, v) = bra_tensor(tns[v], auxinds(tns, v)) + +# The dangling non-physical legs of a vertex tensor: dangling legs that are not site indices. +auxinds(tns::TensorNetworkState, v) = Index[i for i in namesetdiff(uniqueinds(tns, v), siteinds(tns, v))] + +# `auxinds_f` overrides the live aux-leg classification. The loop-correction weights +# need this: they deliberately relabel a bond so it dangles in the modified network, +# and the live classification would misread it as a charge leg. +function norm_factors(tns::TensorNetworkState, verts::Vector; op_strings::Function = v -> "I", auxinds_f::Function = v -> auxinds(tns, v)) factors = ITensor[] for v in verts sinds = siteinds(tns, v) tnv = tns[v] - tnv_dag = dag(prime(tnv)) + tnv_dag = bra_tensor(tnv, auxinds_f(v)) if op_strings(v) == "ρ" || isempty(sinds) append!(factors, ITensor[tnv, tnv_dag]) elseif op_strings(v) == "I" @@ -69,9 +92,13 @@ end norm_factors(tns::TensorNetworkState, v; kwargs...) = norm_factors(tns, [v]; kwargs...) bp_factors(tns::TensorNetworkState, v) = norm_factors(tns, v) +# The flat starting message is the identity between the ket links and their bra +# copies, built as an identity operator so it follows the links' backend (graded +# links give a graded message; the legacy `delta` filled a dense diagonal). function default_message(tns::TensorNetworkState, edge::AbstractEdge) linds = virtualinds(tns, edge) - return adapt_like(tns, denseblocks(delta(vcat(linds, prime(dag(linds)))))) + cod, dom = Tuple(linds), Tuple(prime(dag(linds))) + return adapt_like(tns, one(zeros(scalartype(tns), cod..., dom...), cod, dom)) end """ @@ -146,18 +173,22 @@ function tensornetworkstate(eltype, f::Function, g::AbstractGraph, siteinds::Dic if tnv isa String set!(tensors, v, adapt(eltype)(ITensors.state(f(v), only(siteinds[v])))) elseif tnv isa Vector{<:Number} - set!(tensors, v, adapt(eltype)(ITensors.ITensor(f(v), only(siteinds[v])))) + set!(tensors, v, adapt(eltype)(ITensors.state(f(v), only(siteinds[v])))) else error("Unrecognized local state constructor. Currently supported: Strings and Vectors.") end end - l = Dict(e => Index(1) for e in edges(g)) + # Trivial links, minted to match the site indices' backend (`trivialrange` of a + # graded range is the length-1 trivial-sector range) so graded product states + # stay graded. The dst side takes the conjugate so the pair contracts to 1. for e in edges(g) - tensors[src(e)] *= onehot(eltype, l[e] => 1) - tensors[dst(e)] *= onehot(eltype, l[e] => 1) + l = Index(trivialrange(unnamed(only(siteinds[src(e)])))) + x = fill!(similar(tensors[src(e)], (l,)), true) + tensors[src(e)] *= x + tensors[dst(e)] *= conj(x) end - return TensorNetworkState(tensors, g) + return TensorNetworkState(TensorNetwork(tensors, g), siteinds) end """ @@ -189,6 +220,5 @@ function tensornetworkstate(f::Function, args...) end function NamedGraphs.vertices(t::ITensor, tns::TensorNetworkState) - t_inds = inds(t) - return filter(v -> !isempty(intersect(t_inds, siteinds(tns, v))), collect(vertices(tns))) + return filter(v -> !nameisdisjoint(inds(t), siteinds(tns, v)), collect(vertices(tns))) end diff --git a/src/contraction_sequences.jl b/src/contraction_sequences.jl index f3f5732..bcbe3c6 100644 --- a/src/contraction_sequences.jl +++ b/src/contraction_sequences.jl @@ -10,9 +10,13 @@ using OMEinsumContractionOrders: OMEinsumContractionOrders, optimize_code, EinCo # untouched originals), so no placeholder tensor is needed. is_trivial_tensor(t::ITensor) = all(i -> dim(i) == 1, inds(t)) +# The sequence optimizers only use indices as opaque labels (plus their dimension), so +# hand them the index names: a shared leg is stored nondual on one tensor and dual on the +# other, and comparing whole `Index` objects would see two different labels (or, for a +# space-backed index, fail to compare at all). function contraction_network(tensors::Vector{<:ITensor}; prune_tensors = false) return map(tensors) do t - is = collect(inds(t)) + is = collect(name.(inds(t))) (prune_tensors && is_trivial_tensor(t)) ? empty(is) : is end end @@ -20,7 +24,7 @@ end function contraction_sequence(::Algorithm"optimal", tensors::Vector{<:ITensor}; prune_tensors = false) network = contraction_network(tensors; prune_tensors) #Converting dims to Float64 to minimize overflow issues - inds_to_dims = Dict(i => Float64(dim(i)) for i in unique(Iterators.flatten(network))) + inds_to_dims = Dict(name(i) => Float64(dim(i)) for t in tensors for i in inds(t)) seq, _ = optimaltree(network, inds_to_dims) seq = typeof(seq) <: Int ? [seq] : seq return seq @@ -38,10 +42,10 @@ end #OMEinsumContractionOrders helpers function to_eincode(tensors::Vector{<:ITensor}) - ixs = map(t -> collect(inds(t)), tensors) + ixs = map(t -> collect(name.(inds(t))), tensors) LT = eltype(eltype(ixs)) - iy = collect(LT, reduce(noncommoninds, tensors)) - size_dict = Dict{LT, Int}(i => dim(i) for ix in ixs for i in ix) + iy = collect(LT, name.(reduce(noncommoninds, tensors))) + size_dict = Dict{LT, Int}(name(i) => dim(i) for t in tensors for i in inds(t)) return EinCode(ixs, iy), size_dict end diff --git a/src/entanglement.jl b/src/entanglement.jl index 33c0f34..0a93c0c 100644 --- a/src/entanglement.jl +++ b/src/entanglement.jl @@ -30,7 +30,9 @@ end function matricize(a::ITensor, row_inds = filter(i -> plev(i) ==0, inds(a))) col_inds = prime.(row_inds) - return ITensors.array(a * ITensors.combiner(row_inds) * ITensors.combiner(col_inds)) + row_name = ITensorBase.uniquename(ITensorBase.IndexName) + col_name = ITensorBase.uniquename(ITensorBase.IndexName) + return ITensors.array(ITensors.matricize(a, Tuple(row_inds) => row_name, Tuple(col_inds) => col_name)) end """ @@ -78,7 +80,11 @@ function renyi_entropy( ee = 0 m1, m2 = message(bp_cache, e), message(bp_cache, reverse(e)) edge_ind = only(virtualinds(bp_cache, e)) - root_m2 = first(pseudo_sqrt_inv_sqrt(m2)) + root_m2 = sqrth_safe( + project_hermitian(m2, (inds(m2)[1],), (inds(m2)[2],)), + (inds(m2)[1],), (inds(m2)[2],); + atol = 10 * eps(real(scalartype(m2))), rtol = 0 + ) edge_ind_p, edge_ind_pp = prime(edge_ind), prime(prime(edge_ind)) ρ = (m1 * replaceind(root_m2, edge_ind_p, edge_ind_pp)) * root_m2 diff --git a/src/imports.jl b/src/imports.jl index a804b80..c9c85c6 100644 --- a/src/imports.jl +++ b/src/imports.jl @@ -54,11 +54,10 @@ using TensorOperations # referenced `ITensors`, so the per-file imports mirror the original ITensors-based code. import .ITensorsITensorBaseCompat as ITensors using .ITensorsITensorBaseCompat: - inds, commoninds, commonind, uniqueinds, noncommonind, noncommoninds, unioninds, hascommoninds, + inds, commoninds, commonind, uniqueinds, noncommonind, noncommoninds, unioninds, hascommoninds, cat_inds, sim, dag, prime, noprime, replaceind, replaceinds, dim, swapind, - itensor, random_itensor, scalar, delta, onehot, combiner, combinedind, - qr, svd, eigen, factorize, factorize_svd, - map_diag, map_diag!, + itensor, random_itensor, scalar, delta, similar_map, onehot, + qr, svd, svd_trunc, eigen, factorize, itensor_trunc, scalartype, datatype, array, data, denseblocks, dense, hasqns, contract, inner, apply, @@ -66,7 +65,10 @@ using .ITensorsITensorBaseCompat: Algorithm, @Algorithm_str, hastags, state, op, OpName, SiteType, @OpName_str, @SiteType_str -using ITensorBase: ITensorBase, Index, ITensor, name, plev, tags +using ITensorBase: ITensorBase, Index, ITensor, name, plev, tags, unnamed +using TensorAlgebra: trivialrange +using TensorAlgebra.MatrixAlgebra: sqrth_invsqrth_safe, sqrth_safe +using MatrixAlgebraKit: project_hermitian using Adapt: adapt diff --git a/src/rdm.jl b/src/rdm.jl index bf942e0..43b8ac1 100644 --- a/src/rdm.jl +++ b/src/rdm.jl @@ -1,10 +1,5 @@ function normalize_rdm(ρ::ITensor) - dtype = datatype(ρ) - tr_ρ = copy(ρ) - for i in inds(ρ; plev = 0) - tr_ρ *= adapt(dtype)(delta(i, prime(i))) - end - return ρ / scalar(tr_ρ) + return ρ / ITensors.tr(ρ) end """ diff --git a/src/sampling.jl b/src/sampling.jl index 1051fea..635d447 100644 --- a/src/sampling.jl +++ b/src/sampling.jl @@ -20,7 +20,7 @@ function sample( bit_string = Dictionary{keytype(vertices(ψ)), Int}() for v in vertices(ψ) tensors = incoming_messages(projected_bp_cache, v) - ψv, ψv_dag = network(projected_bp_cache)[v], dag(prime(network(projected_bp_cache)[v])) + ψv, ψv_dag = network(projected_bp_cache)[v], bra_tensor(network(projected_bp_cache), v) push!(tensors, ψv, ψv_dag) seq = contraction_sequence(tensors; alg = "optimal") ρ = ITensors.contract(tensors; sequence = seq) @@ -32,7 +32,7 @@ function sample( # config is 1,2,...,d, but we want 0,1...,d-1 for the sample itself set!(bit_string, v, config - 1) s_ind = inds(ρ)[findfirst(i -> plev(i) == 0, inds(ρ))] - P = adapt_like(ρ, onehot(s_ind => config)) + P = adapt_like(ρ, dag(onehot(s_ind => config))) setindex_preserve!(projected_bp_cache, ψv * P, v) if v != last(vertices(ψ)) @@ -198,7 +198,13 @@ function get_one_sample( es = sorted_edges(norm_bmps_cache, pe) for (i, e) in enumerate(es) - setmessage!(norm_bmps_cache, e, ITensor[outgoing_mps[i], prime(dag(outgoing_mps[i]))]) + mt = outgoing_mps[i] + # Legs unique to `mt` (not facing the network, not MPS chain bonds) are + # dangling non-physical legs the applied row carried into the message + # (e.g. charge legs of projected sites); the bra copy pairs them. + net_inds = virtualinds(network(norm_bmps_cache), e) + aux = uniqueinds(mt, cat_inds(net_inds, (inds(m) for m in outgoing_mps if m !== mt)...)) + setmessage!(norm_bmps_cache, e, ITensor[mt, bra_tensor(mt, aux)]) end incoming_mps = outgoing_mps @@ -228,7 +234,7 @@ function sample_partition!( !isnothing(prev_v) && update_partition!(norm_bmps_cache, [NamedEdge(prev_v => v)]) incoming_ms = incoming_messages(norm_bmps_cache, [v]) ψv = network(norm_bmps_cache)[v] - ψvdag = dag(prime(ψv)) + ψvdag = bra_tensor(network(norm_bmps_cache), v) ts = [incoming_ms; [ψv, ψvdag]] seq = contraction_sequence(ts; alg = "optimal") ρ = contract(ts; sequence = seq) @@ -240,7 +246,7 @@ function sample_partition!( # config is 1,2,...,d, but we want 0,1...,d-1 for the sample itself set!(bit_string, v, config - 1) s_ind = inds(ρ)[findfirst(i -> plev(i) == 0, inds(ρ))] - P = adapt_like(ρ, onehot(s_ind => config)) + P = adapt_like(ρ, dag(onehot(s_ind => config))) q = ρ_diag[config] logq += log(q) Pψv = copy(network(norm_bmps_cache)[v]) * inv(sqrt(q)) * P @@ -268,7 +274,7 @@ function certify_sample( s = siteinds(ψ) qv = sqrt(exp(inv(oftype(logq, length(vertices(ψ)))) * logq)) for v in vertices(ψ) - P = adapt_like(ψproj[v], onehot(only(s[v]) => bitstring[v] + 1)) + P = adapt_like(ψproj[v], dag(onehot(only(s[v]) => bitstring[v] + 1))) setindex_preserve!(ψproj, ψproj[v] * P * inv(qv), v) end diff --git a/src/symmetric_gauge.jl b/src/symmetric_gauge.jl index f5aa129..7c2bda5 100644 --- a/src/symmetric_gauge.jl +++ b/src/symmetric_gauge.jl @@ -1,5 +1,4 @@ function symmetric_gauge!(bp_cache::BeliefPropagationCache; regularization = 10 * eps(real(scalartype(bp_cache))), kwargs...) - dtype = datatype(bp_cache) tn = network(bp_cache) !(tn isa TensorNetworkState) && error("Can only transform TensorNetworkStates to the symmetric gauge") for e in edges(tn) @@ -7,49 +6,58 @@ function symmetric_gauge!(bp_cache::BeliefPropagationCache; regularization = 10 ψvsrc, ψvdst = tn[vsrc], tn[vdst] edge_ind = commoninds(ψvsrc, ψvdst) + edge_ind_p = prime.(edge_ind) edge_ind_sim = sim(edge_ind) - X_D, X_U = safe_eigen(message(bp_cache, e); ishermitian = true, cutoff = nothing) - Y_D, Y_U = safe_eigen(message(bp_cache, reverse(e)); ishermitian = true, cutoff = nothing) - X_D, Y_D = ITensors.map_diag(x -> x + regularization, X_D), - ITensors.map_diag(x -> x + regularization, Y_D) - - rootX_D, rootY_D = ITensors.map_diag(x -> sqrt(x), X_D), ITensors.map_diag(x -> sqrt(x), Y_D) - inv_rootX_D, inv_rootY_D = ITensors.map_diag(x -> inv(sqrt(x)), X_D), - ITensors.map_diag(x -> inv(sqrt(x)), Y_D) - rootX = X_U * rootX_D * prime(dag(X_U)) - rootY = Y_U * rootY_D * prime(dag(Y_U)) - inv_rootX = X_U * inv_rootX_D * prime(dag(X_U)) - inv_rootY = Y_U * inv_rootY_D * prime(dag(Y_U)) - - ψvsrc, ψvdst = noprime(ψvsrc * inv_rootX), noprime(ψvdst * inv_rootY) - - Ce = rootX - Ce = Ce * replaceinds(rootY, edge_ind, edge_ind_sim) - - U, S, V = svd(Ce, edge_ind; kwargs...) - - new_edge_ind = Index[Index(dim(commoninds(S, U)), tags(first(edge_ind)))] - - ψvsrc = replaceinds(ψvsrc * U, commoninds(S, U), new_edge_ind) - ψvdst = replaceinds(ψvdst, edge_ind, edge_ind_sim) - ψvdst = replaceinds(ψvdst * V, commoninds(V, S), new_edge_ind) - - - S = replaceinds( - S, - [commoninds(S, U)..., commoninds(S, V)...] => - [new_edge_ind..., prime(new_edge_ind)...], + # Hermitian square roots (and pseudo-inverse roots) of the two messages, as + # operators with codomain on the unprimed and domain on the primed copies of the + # bond indices. The messages are hermitian only up to numerical noise, so project + # first. Eigenvalues below `regularization` are clamped to zero. + me = project_hermitian(message(bp_cache, e), Tuple(edge_ind), Tuple(edge_ind_p)) + mer = project_hermitian(message(bp_cache, reverse(e)), Tuple(edge_ind), Tuple(edge_ind_p)) + rootX, inv_rootX = sqrth_invsqrth_safe( + me, Tuple(edge_ind), Tuple(edge_ind_p); + atol = regularization, rtol = 0, + ) + rootY, inv_rootY = sqrth_invsqrth_safe( + mer, Tuple(edge_ind), Tuple(edge_ind_p); + atol = regularization, rtol = 0, ) - sqrtS = ITensors.map_diag(sqrt, S) - ψvsrc = noprime(ψvsrc * sqrtS) - ψvdst = noprime(ψvdst * sqrtS) + # Absorb the inverse roots through their domain legs: move the state's bond + # indices onto the primed names so each pairs with its dual copy, and the + # result comes back out on the unprimed (codomain) names. + ψvsrc = replaceinds(ψvsrc, edge_ind, edge_ind_p) * inv_rootX + ψvdst = replaceinds(ψvdst, edge_ind, edge_ind_p) * inv_rootY + + # Bond matrix: contract the codomain legs of the two roots (one from each side + # of the edge, mutually dual). Its open legs are then dual to the state + # tensors' bond legs, so the SVD factors below absorb without any flips. + Ce = rootX * replaceinds(rootY, edge_ind_p, edge_ind_sim) + + U, S, V = svd(Ce, edge_ind_p; kwargs...) + u, v = commoninds(S, U), commoninds(S, V) + + ψvsrc = replaceinds(ψvsrc, edge_ind, edge_ind_p) * U + ψvdst = replaceinds(ψvdst, edge_ind, edge_ind_sim) * V + + # Split the singular values symmetrically into both endpoints. `S`'s legs are + # dual to the tensors' new bond legs on both sides, so these contract directly. + # `S` is diagonal, so this hits the diagonal fast path on every backend. + sqrtS = sqrth_safe(S, Tuple(u), Tuple(v); atol = 0, rtol = 0) + ψvsrc = ψvsrc * sqrtS + ψvdst = ψvdst * sqrtS + + new_edge_ind = Index[ITensors.settags(only(v), tags(first(edge_ind)))] + ψvsrc = replaceinds(ψvsrc, v, new_edge_ind) + ψvdst = replaceinds(ψvdst, u, new_edge_ind) setindex_preserve!(bp_cache, ψvsrc, vsrc) setindex_preserve!(bp_cache, ψvdst, vdst) - setmessage!(bp_cache, e, S) - setmessage!(bp_cache, reverse(e), dag(S)) + # The gauged network's messages are the singular values, with the unprimed leg + # carrying the producing side's bond copy (the message convention). + setmessage!(bp_cache, e, replaceinds(S, [only(u), only(v)], [prime(only(new_edge_ind)), only(new_edge_ind)])) + setmessage!(bp_cache, reverse(e), replaceinds(S, [only(u), only(v)], [only(new_edge_ind), prime(only(new_edge_ind))])) end return bp_cache diff --git a/src/utils.jl b/src/utils.jl index 3e61306..1d39e0e 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -15,32 +15,14 @@ function is_ring_graph(g::AbstractGraph) return is_line_graph(g_mod) end -function pseudo_sqrt_inv_sqrt(M::ITensor; cutoff = 10 * eps(real(scalartype(M)))) - @assert length(inds(M)) == 2 - Q, D, Qdag = eigendecomp(M, inds(M)[1], inds(M)[2]; ishermitian = true) - D_sqrt = ITensors.map_diag(x -> iszero(x) || abs(x) < cutoff ? 0 : sqrt(x), D) - D_inv_sqrt = ITensors.map_diag(x -> iszero(x) || abs(x) < cutoff ? 0 : inv(sqrt(x)), D) - M_sqrt = Q * D_sqrt * Qdag - M_inv_sqrt = Q * D_inv_sqrt * Qdag - return M_sqrt, M_inv_sqrt -end - -#TODO: Make this work for non-hermitian A -function eigendecomp(A::ITensor, linds, rinds; ishermitian = false, kwargs...) - @assert ishermitian - D, U = safe_eigen(A, linds, rinds; ishermitian, kwargs...) - ul, ur = noncommonind(D, U), commonind(D, U) - Ul = replaceinds(U, vcat([rinds], [ur]), vcat([linds], [ul])) - return Ul, D, dag(U) -end - # Adapt `t` to the storage datatype (eltype + device) of `ref`. adapt_like(ref, t) = adapt(datatype(ref))(t) +# The fused identity between the two index groups, via `Base.one` (graded-capable); +# `one` takes a prototype tensor for the names/eltype/spaces, which `zeros` supplies. function identity_tensor(eltype, row_inds::Vector{<:Index}, col_inds::Vector{<:Index}) - c_row, c_col = ITensors.combiner(row_inds),ITensors.combiner(col_inds) - t= ITensors.denseblocks(ITensors.delta(eltype, ITensors.combinedind(c_row), ITensors.combinedind(c_col))) - return (t * c_row)*c_col + row_is, col_is = Tuple(row_inds), Tuple(col_inds) + return one(zeros(eltype, row_is, col_is), row_is, col_is) end identity_tensor(row_inds::Vector{<:Index}, col_inds::Vector{<:Index}) = identity_tensor(Float64, row_inds, col_inds) diff --git a/test/runtests.jl b/test/runtests.jl index 641e25a..1cec2bc 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -11,4 +11,5 @@ using Test include("test_sampling.jl") include("test_truncate.jl") include("test_contraction_sequences.jl") + include("test_examples.jl") end diff --git a/test/test_contraction_sequences.jl b/test/test_contraction_sequences.jl index 0ff76f6..da0e44b 100644 --- a/test/test_contraction_sequences.jl +++ b/test/test_contraction_sequences.jl @@ -1,5 +1,5 @@ @eval module $(gensym()) -using ITensorBase: Index +using ITensorBase: Index, name using Random using TensorNetworkQuantumSimulator const TNQS = TensorNetworkQuantumSimulator @@ -18,13 +18,15 @@ collect_leaves!(acc, x) = (for y in x; collect_leaves!(acc, y); end; acc) # --- to_eincode: ITensors -> (EinCode, size_dict). Tests the omeinsum-specific # input conversion directly, so a silent fallback to another backend can't pass it. + # Labels are index names: a shared leg's `Index` differs between its two tensors + # under a graded backend (nondual vs dual), so names are the backend-stable label. i, j, k = Index(2), Index(3), Index(4) A = ITensors.random_itensor(i, j) B = ITensors.random_itensor(j, k) code, size_dict = TNQS.to_eincode([A, B]) - @test Set(Set.(getixsv(code))) == Set([Set([i, j]), Set([j, k])]) # per-tensor index sets - @test Set(getiyv(code)) == Set([i, k]) # open indices (j is contracted) - @test size_dict == Dict(i => 2, j => 3, k => 4) + @test Set(Set.(getixsv(code))) == Set([Set(name.([i, j])), Set(name.([j, k]))]) # per-tensor index sets + @test Set(getiyv(code)) == Set(name.([i, k])) # open indices (j is contracted) + @test size_dict == Dict(name(i) => 2, name(j) => 3, name(k) => 4) # --- to_contraction_sequence: NestedEinsum -> nested tensor-position tree. Tests our # converter on hand-built trees with known shapes (deterministic, exact).