From 49192f276d0fa9f784f5a9a8902be7dc4ae8c090 Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Wed, 8 Jul 2026 20:01:50 -0400 Subject: [PATCH 01/21] Remove the ITensorsITensorBaseCompat submodule Move its contents to their permanent homes: the operator and named-state system to a TNQS-owned `Ops` submodule, the remaining legacy tensor helpers to `itensors.jl`, and the index and tensor operations ITensorBase now provides (`prime`, `noprime`, `id`, `sim`, `scalartype`) to the upstream versions. Drop the legacy no-ops `disable_warn_order`, `denseblocks`, and `dense`, and name the ITensor trace `itensor_tr` so it stays distinct from `tr` on plain matrices. --- Project.toml | 10 +- docs/src/gates.md | 17 +- .../2dIsing_dynamics_Heisenbergpicture.jl | 10 +- examples/boundarymps.jl | 1 - examples/heavyhexIsing_dynamics.jl | 2 +- .../hexagonal_heisenbergmodel_thermalstate.jl | 4 +- examples/loopcorrections.jl | 1 - src/Apply/apply_gates.jl | 2 +- src/Apply/full_update.jl | 13 +- src/Apply/gate_definitions.jl | 45 ++--- src/Apply/simple_update.jl | 6 +- src/Forms/abstractform.jl | 4 +- src/Forms/quadraticform.jl | 2 +- .../ITensorsITensorBaseCompat.jl | 54 ----- .../abstractbeliefpropagationcache.jl | 4 +- src/MessagePassing/beliefpropagationcache.jl | 12 +- src/MessagePassing/boundarympscache.jl | 4 +- .../ops.jl => Ops.jl} | 44 ++-- src/TensorNetworkQuantumSimulator.jl | 5 +- src/TensorNetworks/abstracttensornetwork.jl | 23 +-- src/TensorNetworks/tensornetwork.jl | 4 +- src/TensorNetworks/tensornetworkstate.jl | 9 +- .../tensornetworkstate_constructors.jl | 16 +- src/contract.jl | 12 +- src/contraction_sequences.jl | 1 - src/entanglement.jl | 2 +- src/expect.jl | 1 - src/imports.jl | 26 +-- src/inner.jl | 18 +- .../itensor.jl => itensors.jl} | 189 +++++++----------- src/rdm.jl | 3 +- src/sampling.jl | 10 +- src/siteinds.jl | 2 +- src/symmetric_gauge.jl | 4 +- src/truncate.jl | 16 +- src/utils.jl | 8 +- test/test_apply.jl | 5 +- test/test_constructors.jl | 12 +- test/test_contraction_sequences.jl | 19 +- 39 files changed, 241 insertions(+), 379 deletions(-) delete mode 100644 src/ITensorsITensorBaseCompat/ITensorsITensorBaseCompat.jl rename src/{ITensorsITensorBaseCompat/ops.jl => Ops.jl} (72%) rename src/{ITensorsITensorBaseCompat/itensor.jl => itensors.jl} (79%) diff --git a/Project.toml b/Project.toml index 9788744..38023f9 100644 --- a/Project.toml +++ b/Project.toml @@ -25,11 +25,19 @@ TensorAlgebra = "68bd88dc-f39d-4e12-b2ca-f046b68fcc6a" TensorOperations = "6aa20fa7-93e2-5fca-9bc0-fbd0db3c71a2" TypeParameterAccessors = "7e5a90cf-f82e-492e-a09b-e3e26432c138" +[sources.ITensorBase] +url = "https://github.com/ITensor/ITensorBase.jl" +rev = "mf/nextgen-followups-round1" + +[sources.TensorAlgebra] +url = "https://github.com/ITensor/TensorAlgebra.jl" +rev = "mf/nextgen-followups-round1" + [compat] Adapt = "4.3.0" Dictionaries = "0.4" Graphs = "1.8.0" -ITensorBase = "0.10.9" +ITensorBase = "0.11" KrylovKit = "0.10.2" LinearAlgebra = "1.11.0" MatrixAlgebraKit = "0.6.8" diff --git a/docs/src/gates.md b/docs/src/gates.md index 1b19529..9d2b5ef 100644 --- a/docs/src/gates.md +++ b/docs/src/gates.md @@ -122,15 +122,16 @@ gate2 = ITensor(my_nn_gate, s[v1], s[v2], s[v1]', s[v2]') For a gate you'll reuse repeatedly, register it once and then use the tuple form like any built-in. Two steps: -1. Define the matrix as an `ITensors.op` method. -2. Call [`register_gate!`](@ref) to tell the dispatcher which keyword arguments your `op` accepts. +1. Define the matrix as an `Ops.op` method. +2. Call [`register_gate!`](@ref) to tell the dispatcher which keyword arguments your `Ops.op` accepts. ```julia -using ITensors: op, OpName, SiteType, Index +using TensorNetworkQuantumSimulator: Ops, OpName, SiteType, @OpName_str, @SiteType_str +using ITensorBase: Index # 1. Define the matrix (the "physics" part) -function ITensors.op(::OpName"FSim", ::SiteType"S=1/2", s1::Index, s2::Index; - θ::Number, ϕ::Number) +function Ops.op(::OpName"FSim", ::SiteType"S=1/2", s1::Index, s2::Index; + θ::Number, ϕ::Number) # ... return a 4-leg ITensor ... end @@ -144,9 +145,9 @@ circuit = [("FSim", [v1, v2], (π/4, π/8))] The keyword arguments map as: -- `paramkeys = (:θ, :ϕ)` — names of the kwargs your `op` expects, in the order they appear in the circuit-tuple parameter. For a single-parameter gate, use a 1-tuple like `(:θ,)`. -- `opname = "FSim"` — defaults to the gate name; override if your circuit-level name should differ from the `OpName` your `op` uses. -- `rescale = identity` — applied to the parameter(s) before forwarding to `op`. Useful when your `op` follows a different convention (e.g., the built-in `"Rxx"` uses `rescale = θ -> θ/2` to bridge our qiskit convention to ITensors'). +- `paramkeys = (:θ, :ϕ)` — names of the kwargs your `Ops.op` expects, in the order they appear in the circuit-tuple parameter. For a single-parameter gate, use a 1-tuple like `(:θ,)`. +- `opname = "FSim"` — defaults to the gate name; override if your circuit-level name should differ from the `OpName` your `Ops.op` uses. +- `rescale = identity` — applied to the parameter(s) before forwarding to `Ops.op`. Useful when your `Ops.op` follows a different convention (e.g., the built-in `"Rxx"` uses `rescale = θ -> θ/2` to bridge our qiskit convention to the half-angle convention of the matrix). To add a qiskit-style alias for your gate (so e.g. `"fsim"` resolves to `"FSim"`), call [`register_alias!`](@ref): diff --git a/examples/2dIsing_dynamics_Heisenbergpicture.jl b/examples/2dIsing_dynamics_Heisenbergpicture.jl index 3f1801d..313b6b9 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 TensorNetworkQuantumSimulator: ITensors, ITensor +using TensorNetworkQuantumSimulator: Ops, ITensor function main() nx, ny = 4, 4 @@ -14,7 +14,7 @@ function main() #Start from the identity operator, then place a single Z on vertex vz ψI = identity_tensornetworkstate(ComplexF64, g, s) ψ0 = copy(ψI) - setindex_preserve!(ψ0, noprime(ψ0[vz] * ITensors.op("Z", s[vz][1])), vz) + setindex_preserve!(ψ0, noprime(ψ0[vz] * Ops.op("Z", s[vz][1])), vz) maxdim, cutoff = 4, 1.0e-14 apply_kwargs = (; maxdim, cutoff, normalize_tensors = false) @@ -34,11 +34,11 @@ function main() ec = edge_color(g, 4) #Ket leg (s[v][1]) gets U† (angle negated), bra leg (s[v][2]) gets U (angle unchanged) so that O -> U'OU - append!(layer, [ITensors.op("Rz", s[v][1]; θ = -h * δt)*ITensors.op("Rz", s[v][2]; θ = h * δt) for v in vertices(g)]) + append!(layer, [Ops.op("Rz", s[v][1]; θ = -h * δt)*Ops.op("Rz", s[v][2]; θ = h * δt) for v in vertices(g)]) for es in ec - append!(layer, [ITensors.op("Rxx", s[src(e)][1], s[dst(e)][1], ϕ = -J * δt)*ITensors.op("Rxx", s[src(e)][2], s[dst(e)][2], ϕ = J * δt) for e in es]) + append!(layer, [Ops.op("Rxx", s[src(e)][1], s[dst(e)][1], ϕ = -J * δt)*Ops.op("Rxx", s[src(e)][2], s[dst(e)][2], ϕ = J * δt) for e in es]) end - append!(layer, [ITensors.op("Rz", s[v][1]; θ = -h * δt)*ITensors.op("Rz", s[v][2]; θ = h * δt) for v in vertices(g)]) + append!(layer, [Ops.op("Rz", s[v][1]; θ = -h * δt)*Ops.op("Rz", s[v][2]; θ = h * δt) for v in vertices(g)]) χinit = maxvirtualdim(ψ) println("Initial bond dimension of the Heisenberg operator is $χinit") diff --git a/examples/boundarymps.jl b/examples/boundarymps.jl index f4ba696..4b45545 100644 --- a/examples/boundarymps.jl +++ b/examples/boundarymps.jl @@ -1,5 +1,4 @@ using TensorNetworkQuantumSimulator -using TensorNetworkQuantumSimulator: ITensors using Random Random.seed!(1634) diff --git a/examples/heavyhexIsing_dynamics.jl b/examples/heavyhexIsing_dynamics.jl index 62b87dc..4f9df1b 100644 --- a/examples/heavyhexIsing_dynamics.jl +++ b/examples/heavyhexIsing_dynamics.jl @@ -1,6 +1,6 @@ using TensorNetworkQuantumSimulator using Statistics -using TensorNetworkQuantumSimulator: ITensor, ITensors +using TensorNetworkQuantumSimulator: ITensor function main() #Define the lattice diff --git a/examples/hexagonal_heisenbergmodel_thermalstate.jl b/examples/hexagonal_heisenbergmodel_thermalstate.jl index 0c9f08e..397a2c0 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 TensorNetworkQuantumSimulator: ITensors, ITensor +using TensorNetworkQuantumSimulator: Ops, ITensor function main() χ = 32 @@ -17,7 +17,7 @@ function main() apply_kwargs= (; maxdim = χ, cutoff = 1e-14, normalize_tensors = false) two_site_gates =ITensor[] for es in ec - append!(two_site_gates, [ITensors.op("Rxxyyzz", s[src(e)][1], s[dst(e)][1], θ = -0.5*J*dβ*im) for e in es]) + append!(two_site_gates, [Ops.op("Rxxyyzz", s[src(e)][1], s[dst(e)][1], θ = -0.5*J*dβ*im) for e in es]) end nsteps = 25 diff --git a/examples/loopcorrections.jl b/examples/loopcorrections.jl index e5f65fd..c30ae5a 100644 --- a/examples/loopcorrections.jl +++ b/examples/loopcorrections.jl @@ -1,5 +1,4 @@ using TensorNetworkQuantumSimulator -using TensorNetworkQuantumSimulator: ITensors using LinearAlgebra: norm diff --git a/src/Apply/apply_gates.jl b/src/Apply/apply_gates.jl index 3c9798b..4f5e045 100644 --- a/src/Apply/apply_gates.jl +++ b/src/Apply/apply_gates.jl @@ -39,7 +39,7 @@ function apply_gates( end function adapt_gate(gate::ITensor, ψ_bpc::BeliefPropagationCache) - gate = scalartype(gate) <: Complex ? ITensors.adapt_scalartype(complex(scalartype(ψ_bpc)), gate) : ITensors.adapt_scalartype(scalartype(ψ_bpc), gate) + gate = scalartype(gate) <: Complex ? adapt_scalartype(complex(scalartype(ψ_bpc)), gate) : adapt_scalartype(scalartype(ψ_bpc), gate) return adapt(unspecify_type_parameters(datatype(ψ_bpc)), gate) end diff --git a/src/Apply/full_update.jl b/src/Apply/full_update.jl index 75b81fc..ed0bf55 100644 --- a/src/Apply/full_update.jl +++ b/src/Apply/full_update.jl @@ -1,5 +1,4 @@ using KrylovKit: linsolve -using .ITensorsITensorBaseCompat: namesetdiff """ full_update(o::ITensor, ψ::TensorNetworkState, v⃗; envs, kwargs...) @@ -82,7 +81,7 @@ function fidelity( envs, ) sequence = contraction_sequence(term1_tns; alg = "optimal") - term1 = ITensors.contract(term1_tns; sequence) + term1 = contract(term1_tns; sequence) term2_tns = vcat( [ @@ -94,10 +93,10 @@ function fidelity( envs, ) sequence = contraction_sequence(term2_tns; alg = "optimal") - term2 = ITensors.contract(term2_tns; sequence) + term2 = contract(term2_tns; sequence) term3_tns = vcat([p_prev, q_prev, prime(dag(p_cur)), prime(dag(q_cur)), gate], envs) sequence = contraction_sequence(term3_tns; alg = "optimal") - term3 = ITensors.contract(term3_tns; sequence) + term3 = contract(term3_tns; sequence) f = scalar(term3) / sqrt(scalar(term1) * scalar(term2)) return f * conj(f) @@ -127,18 +126,18 @@ function optimise_p_q( function b(p::ITensor, q::ITensor, o::ITensor, envs::Vector{ITensor}, r::ITensor) ts = vcat(ITensor[p, q, o, dag(prime(r))], envs) sequence = contraction_sequence(ts; alg = "optimal") - return noprime(ITensors.contract(ts; sequence)) + return noprime(contract(ts; sequence)) end function M_p(envs::Vector{ITensor}, p_q_tensor::ITensor, s_ind, apply_tensor::ITensor) ts = vcat( ITensor[ - p_q_tensor, replaceinds(prime(dag(p_q_tensor)), prime(s_ind), s_ind), apply_tensor, + p_q_tensor, replaceinds(prime(dag(p_q_tensor)), prime.(s_ind), s_ind), apply_tensor, ], envs, ) sequence = contraction_sequence(ts; alg = "optimal") - return noprime(ITensors.contract(ts; sequence)) + return noprime(contract(ts; sequence)) end for i in 1:nfullupdatesweeps b_vec = b(p, q, o, envs, q_cur) diff --git a/src/Apply/gate_definitions.jl b/src/Apply/gate_definitions.jl index c348fcd..13b47f6 100644 --- a/src/Apply/gate_definitions.jl +++ b/src/Apply/gate_definitions.jl @@ -1,15 +1,14 @@ -using .ITensorsITensorBaseCompat: hastags # --- Gate registry ----------------------------------------------------------- # Internal dispatch record for a circuit-tuple gate name. # -# - `opname`: the `OpName` string forwarded to `ITensors.op`. Usually equal to the +# - `opname`: the `OpName` string forwarded to `Ops.op`. Usually equal to the # user-facing key, but kept separate so a registry entry can rename if needed. # - `paramkeys`: keyword names accepted by the underlying `op` definition, e.g. # `(:θ,)`, `(:ϕ,)`, or `(:θ, :β)`. Empty for fixed gates. # - `rescale`: applied to the user-supplied parameter(s) before forwarding. Used -# when our (qiskit) convention differs from the `ITensors.op` convention. For +# when our (qiskit) convention differs from the `Ops.op` convention. For # multi-parameter gates, `rescale` receives and returns a tuple/vector. struct GateSpec opname::String @@ -19,7 +18,7 @@ end GateSpec(opname; paramkeys = (), rescale = identity) = GateSpec(opname, paramkeys, rescale) # Registry of circuit-tuple gates. Adding a new gate is one entry here (plus an -# `ITensors.op` method if upstream doesn't already provide one). +# `Ops.op` method if upstream doesn't already provide one). const GATES = Dict{String, GateSpec}( # Single-qubit fixed "X" => GateSpec("X"), @@ -125,7 +124,7 @@ function toitensor(gate::Tuple, g::NamedGraph, siteinds::Dictionary) # Multi-letter Pauli-string sugar: "XYZ" → X⊗Y⊗Z applied componentwise. # Single-letter "X"/"Y"/"Z" goes through the registry below. if _ispaulistring(name) && length(name) > 1 - t = prod(ITensors.op(string(c), sind) for (c, sind) in zip(name, s_inds)) + t = prod(Ops.op(string(c), sind) for (c, sind) in zip(name, s_inds)) return t, verts end @@ -142,7 +141,7 @@ function toitensor(gate::Tuple, g::NamedGraph, siteinds::Dictionary) end if isempty(spec.paramkeys) - return ITensors.op(spec.opname, s_inds...), verts + return Ops.op(spec.opname, s_inds...), verts end raw = spec.rescale(gate[3]) @@ -151,7 +150,7 @@ function toitensor(gate::Tuple, g::NamedGraph, siteinds::Dictionary) "Gate \"$name\" expects $(length(spec.paramkeys)) parameter(s), got $(length(pvals))." )) kwargs = NamedTuple{spec.paramkeys}(pvals) - return ITensors.op(spec.opname, s_inds...; kwargs...), verts + return Ops.op(spec.opname, s_inds...; kwargs...), verts end # --- Public registration API ------------------------------------------------ @@ -162,7 +161,7 @@ end Register a custom gate `name` so it can be used in circuit-tuple form `(name, vertices, parameter)` with `apply_gates`. -The matrix itself must be defined separately as an `ITensors.op` method whose +The matrix itself must be defined separately as an `Ops.op` method whose `OpName` matches `opname` (defaults to `name`). See "Custom Gates" in the gate docs for a worked example. @@ -172,14 +171,14 @@ in your script's startup, or in a downstream package's `__init__()`. Built-in gates are locked: passing a built-in name throws `ArgumentError`. Choose a different name for your custom gate, or — if you really need a new -matrix under an existing name — define your own `ITensors.op` method directly. +matrix under an existing name — define your own `Ops.op` method directly. Previously user-registered names may be overwritten freely. # Arguments - `name`: name used in circuit tuples. # Keyword Arguments -- `opname`: the `OpName` string forwarded to `ITensors.op`. Defaults to `name`. +- `opname`: the `OpName` string forwarded to `Ops.op`. Defaults to `name`. - `paramkeys`: tuple of keyword names accepted by the underlying `op`, e.g. `(:θ,)` for a single rotation angle, `(:θ, :β)` for a two-parameter gate. Empty (`()`) for non-parametric gates. @@ -197,7 +196,7 @@ function register_gate!( name in BUILTIN_GATES && throw(ArgumentError( "\"$name\" is a built-in gate and cannot be overwritten. " * "Choose a different name for your custom gate, or define your own " * - "`ITensors.op` method directly if you need to override the matrix." + "`Ops.op` method directly if you need to override the matrix." )) GATES[name] = GateSpec(opname, paramkeys, rescale) return name @@ -243,11 +242,11 @@ end # --- In-house gate definitions ---------------------------------------------- """ - ITensors.op(::OpName"xx_plus_yy", ::SiteType"S=1/2"; θ::Number, β::Number) + Ops.op(::OpName"xx_plus_yy", ::SiteType"S=1/2"; θ::Number, β::Number) Gate for rotation by XX+YY at a given angle with Rz rotations either side. Consistent with qiskit. """ -function ITensors.op(::OpName"xx_plus_yy", ::SiteType"S=1/2"; θ::Number, β::Number) +function Ops.op(::OpName"xx_plus_yy", ::SiteType"S=1/2"; θ::Number, β::Number) return [ [1 0 0 0]; [0 cos(θ / 2) -im * sin(θ / 2) * exp(-im * β) 0] @@ -255,15 +254,15 @@ function ITensors.op(::OpName"xx_plus_yy", ::SiteType"S=1/2"; θ::Number, β::Nu [0 0 0 1] ] end -ITensors.op(o::OpName"xx_plus_yy", ::SiteType"Qubit"; θ::Number, β::Number) = - ITensors.op(o, ITensors.SiteType("S=1/2"); θ, β) +Ops.op(o::OpName"xx_plus_yy", ::SiteType"Qubit"; θ::Number, β::Number) = + Ops.op(o, SiteType("S=1/2"); θ, β) """ - ITensors.op(::OpName"Rxxyy", ::SiteType"S=1/2"; θ::Number) + Ops.op(::OpName"Rxxyy", ::SiteType"S=1/2"; θ::Number) Gate for rotation by XXYY at a given angle. """ -function ITensors.op(::OpName"Rxxyy", ::SiteType"S=1/2"; θ = 1) +function Ops.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 @@ -272,20 +271,20 @@ function ITensors.op(::OpName"Rxxyy", ::SiteType"S=1/2"; θ = 1) h = kron(σp, σm) + kron(σm, σp) return exp(-im * θ * h) end -ITensors.op(o::OpName"Rxxyy", ::SiteType"Qubit"; θ::Number) = - ITensors.op(o, ITensors.SiteType("S=1/2"); θ) +Ops.op(o::OpName"Rxxyy", ::SiteType"Qubit"; θ::Number) = + Ops.op(o, SiteType("S=1/2"); θ) """ - ITensors.op(::OpName"Rxxyyzz", ::SiteType"S=1/2"; θ::Number) + Ops.op(::OpName"Rxxyyzz", ::SiteType"S=1/2"; θ::Number) Gate for rotation by XXYYZZ at a given angle. """ -function ITensors.op(::OpName"Rxxyyzz", ::SiteType"S=1/2"; θ = 1) +function Ops.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) = - ITensors.op(o, ITensors.SiteType("S=1/2"); θ) +Ops.op(o::OpName"Rxxyyzz", ::SiteType"Qubit"; θ::Number) = + Ops.op(o, SiteType("S=1/2"); θ) diff --git a/src/Apply/simple_update.jl b/src/Apply/simple_update.jl index f4a3722..ea7dfcc 100644 --- a/src/Apply/simple_update.jl +++ b/src/Apply/simple_update.jl @@ -24,7 +24,7 @@ function simple_update( ) if length(ψ⃗) == 1 - updated_tensors = ITensor[ITensors.apply(o, only(ψ⃗))] + updated_tensors = ITensor[apply(o, only(ψ⃗))] s_values, err = nothing, 0 else # When envs is empty no gauging happens and the cutoff is unused, so fall back to @@ -54,7 +54,7 @@ function simple_update( Qᵥ₂, Rᵥ₂ = qr(ψᵥ₂, uniqueinds(uniqueinds(ψᵥ₂, ψᵥ₁), sᵥ₂)) rᵥ₁ = commoninds(Qᵥ₁, Rᵥ₁) rᵥ₂ = commoninds(Qᵥ₂, Rᵥ₂) - oR = ITensors.apply(o, Rᵥ₁ * Rᵥ₂) + oR = apply(o, Rᵥ₁ * Rᵥ₂) # 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 @@ -81,7 +81,7 @@ function simple_update( if normalize_tensors for ψᵥ in updated_tensors - rmul!(ITensors.data(ψᵥ), inv(norm(ψᵥ))) + rmul!(data(ψᵥ), inv(norm(ψᵥ))) end end diff --git a/src/Forms/abstractform.jl b/src/Forms/abstractform.jl index 70b4c2d..84ae69b 100644 --- a/src/Forms/abstractform.jl +++ b/src/Forms/abstractform.jl @@ -6,8 +6,8 @@ abstract type AbstractForm{V} <: AbstractTensorNetwork{V} end #Forward onto the ket for f in [ :(graph), - :(ITensors.datatype), - :(ITensors.scalartype), + :(datatype), + :(scalartype), :(NamedGraphs.edgeinduced_subgraphs_no_leaves), ] @eval begin diff --git a/src/Forms/quadraticform.jl b/src/Forms/quadraticform.jl index 5c6ab3a..e995c7c 100644 --- a/src/Forms/quadraticform.jl +++ b/src/Forms/quadraticform.jl @@ -16,7 +16,7 @@ function QuadraticForm(ket::TensorNetworkState, f::Function = v -> "I") sinds = siteinds(ket) verts = collect(vertices(ket)) dtype = datatype(ket) - operator_tensors = adapt(dtype).([reduce(prod, ITensor[ITensors.op(f(v), sind) for sind in sinds[v]]) for v in verts]) + operator_tensors = adapt(dtype).([reduce(prod, ITensor[Ops.op(f(v), sind) for sind in sinds[v]]) for v in verts]) operator = TensorNetworkState(Dictionary(verts, operator_tensors)) return QuadraticForm(ket, operator) end diff --git a/src/ITensorsITensorBaseCompat/ITensorsITensorBaseCompat.jl b/src/ITensorsITensorBaseCompat/ITensorsITensorBaseCompat.jl deleted file mode 100644 index c88c5d1..0000000 --- a/src/ITensorsITensorBaseCompat/ITensorsITensorBaseCompat.jl +++ /dev/null @@ -1,54 +0,0 @@ -# Compatibility layer providing the legacy `ITensors.jl` API that TNQS was written -# against, implemented over the next-gen `ITensorBase.jl` backend. Wrapped as a -# submodule so it can later be lifted out of TNQS into a standalone package of the -# same name (the migration aid for the whole ecosystem); the export list below is the -# spec for what that package must cover. -# -# Export vs `public`: the eventual standalone package, used as -# `using ITensorBase; using ITensorsITensorBaseCompat`, would declare the names -# ITensorBase also exports (`apply`, `prime`, `noprime`, `state`) as `public` rather -# than `export`ed, so they don't shadow ITensorBase's versions — a consumer wanting -# the legacy semantics imports them explicitly. Here every name is exported instead: -# TNQS does not `using ITensorBase` broadly, so there is no collision to avoid, and -# `public` would force Julia >= 1.11 while this package still supports 1.10. The -# main module imports `contract` / `inner` / `scalartype` / `datatype` (the generics -# TNQS adds methods to) so it can extend them. -module ITensorsITensorBaseCompat - -# Legacy `ITensors`/`ITensorMPS` exported `truncate`; TNQS extends it. Bind it to -# `Base.truncate` (which is always in scope, so a `using` consumer sees one -# unambiguous binding) and re-export, so `ITensors.truncate` resolves through this -# module and `function ITensors.truncate(...)` extends `Base.truncate`. -import Base: truncate - -include("itensor.jl") -include("ops.jl") - -export - # Index access and set algebra - 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, - # Factorizations - qr, svd, svd_trunc, eigen, factorize, itensor_trunc, - # Storage / element-type accessors - scalartype, datatype, array, data, - # Dense / quantum-number no-ops - denseblocks, dense, hasqns, - # Contraction, inner product, gate application - contract, inner, apply, - # Direct sum and misc legacy helpers - directsum, disable_warn_order, - # Algorithm dispatch tag - Algorithm, @Algorithm_str, - # Tags - hastags, - # Operator / named-state system - state, op, OpName, SiteType, @OpName_str, @SiteType_str, - # Bond truncation (bound to `Base.truncate`) - truncate - -end diff --git a/src/MessagePassing/abstractbeliefpropagationcache.jl b/src/MessagePassing/abstractbeliefpropagationcache.jl index 247b1af..dee7e6a 100644 --- a/src/MessagePassing/abstractbeliefpropagationcache.jl +++ b/src/MessagePassing/abstractbeliefpropagationcache.jl @@ -41,8 +41,8 @@ for f in [ :(bp_factors), :(default_bp_maxiter), :(virtualinds), - :(ITensors.datatype), - :(ITensors.scalartype), + :(datatype), + :(scalartype), :(maxvirtualdim), :(default_message), :(siteinds), diff --git a/src/MessagePassing/beliefpropagationcache.jl b/src/MessagePassing/beliefpropagationcache.jl index 224d993..2934036 100644 --- a/src/MessagePassing/beliefpropagationcache.jl +++ b/src/MessagePassing/beliefpropagationcache.jl @@ -1,8 +1,6 @@ using Dictionaries: Dictionary, set!, delete! using Graphs: AbstractGraph, is_tree, connected_components using NamedGraphs.GraphsExtensions: default_root_vertex, forest_cover, post_order_dfs_edges, forest_cover_edge_sequence, boundary_edges, leaf_vertices, a_star -using .ITensorsITensorBaseCompat: dim, ITensor, delta, Algorithm -using .ITensorsITensorBaseCompat: scalartype using LinearAlgebra: normalize #TODO: Make this show() nicely. @@ -111,7 +109,7 @@ function default_bp_update_kwargs(tn::AbstractTensorNetwork) if is_tree(tn) maxiter, tolerance, verbose = 1, nothing, false else - maxiter, tolerance, verbose = _default_bp_update_maxiter, default_tolerance(ITensors.scalartype(tn)), false + maxiter, tolerance, verbose = _default_bp_update_maxiter, default_tolerance(scalartype(tn)), false end return (; maxiter, tolerance, verbose) end @@ -119,9 +117,9 @@ end default_bp_update_kwargs(bp_cache::BeliefPropagationCache) = default_bp_update_kwargs(network(bp_cache)) function make_hermitian(A::ITensor) - A_inds = ITensors.inds(A) + A_inds = inds(A) @assert length(A_inds) == 2 - return (A + ITensors.swapind(dag(A), first(A_inds), last(A_inds))) / 2 + return (A + swapind(dag(A), first(A_inds), last(A_inds))) / 2 end function rescale_messages!(bp_cache::BeliefPropagationCache, edges::Vector{<:AbstractEdge}) @@ -172,9 +170,9 @@ function loop_correlation(bpc::BeliefPropagationCache, loop::Vector{<:NamedEdge} 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 = matricize(t, Tuple(e_virtualinds) => row_name, Tuple(e_virtualinds_sim) => col_name) t = adapt(Vector{ComplexF64})(t) - t = ITensors.array(t) + t = array(t) λs = reverse(sort(LinearAlgebra.eigvals(t); by = abs)) err = 1 - abs(λs[1]) / sum(abs.(λs)) return err diff --git a/src/MessagePassing/boundarympscache.jl b/src/MessagePassing/boundarympscache.jl index 4e0240b..91c6444 100644 --- a/src/MessagePassing/boundarympscache.jl +++ b/src/MessagePassing/boundarympscache.jl @@ -39,7 +39,7 @@ end default_message_update_alg(bmps_cache::BoundaryMPSCache) = default_bmps_message_update_alg(network(bmps_cache)) default_normalize(alg::Algorithm"fitting") = true -default_tolerance(bmps_cache::BoundaryMPSCache) = default_tolerance(ITensors.scalartype(bmps_cache)) +default_tolerance(bmps_cache::BoundaryMPSCache) = default_tolerance(scalartype(bmps_cache)) _default_boundarymps_update_niters = 50 function set_default_kwargs(alg::Algorithm"fitting", bmps_cache::BoundaryMPSCache) normalize = get(alg.kwargs, :normalize, default_normalize(alg)) @@ -195,7 +195,7 @@ function set_interpartition_messages!( # 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)") + ind = 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. diff --git a/src/ITensorsITensorBaseCompat/ops.jl b/src/Ops.jl similarity index 72% rename from src/ITensorsITensorBaseCompat/ops.jl rename to src/Ops.jl index bed0387..1d2a67c 100644 --- a/src/ITensorsITensorBaseCompat/ops.jl +++ b/src/Ops.jl @@ -1,11 +1,16 @@ -# 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. +# Operator / named-state system, scoped to the `S=1/2` single- and two-qubit gate set and the +# computational/Pauli-basis product states, using the `OpName` / `SiteType` dispatch that gate +# definitions (`Apply/gate_definitions.jl`) and user-registered gates extend. # -# 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. +# `state` is called qualified (`Ops.state`) at call sites because ITensorBase also exports an +# unrelated `state`, and the two do not extend each other. +module Ops -using TensorAlgebra: TensorAlgebra, project, tryproject +using ITensorBase: ITensorBase, Index, id +using TensorAlgebra: project +using ..TensorNetworkQuantumSimulator: project_aux + +export state, op, OpName, SiteType, @OpName_str, @SiteType_str # # Named single-site states. `state(name, i)` returns the named state vector as an @@ -22,32 +27,13 @@ const _STATE_VECTORS = Dict{String, Vector{ComplexF64}}( function state(name::AbstractString, i::Index) v = get(_STATE_VECTORS, name) do return error( - "unknown single-site state \"$name\" (vendored states: $(sort(collect(keys(_STATE_VECTORS)))))" + "unknown single-site state \"$name\" (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 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`. +# over `i`, through the `project_aux` utility. state(v::AbstractVector{<:Number}, i::Index) = project_aux(v, i) # @@ -113,7 +99,7 @@ const _GATE_MATRICES = Dict{String, Matrix{ComplexF64}}( function _gate_matrix(name::AbstractString; kwargs...) name in ("Rx", "Ry", "Rz", "P") || return get(_GATE_MATRICES, name) do return error( - "unknown single-site operator \"$name\" (vendored operators: $(sort(collect(keys(_GATE_MATRICES)))) plus Rx/Ry/Rz/P)" + "unknown single-site operator \"$name\" (operators: $(sort(collect(keys(_GATE_MATRICES)))) plus Rx/Ry/Rz/P)" ) end if name == "Rx" @@ -148,3 +134,5 @@ 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 + +end diff --git a/src/TensorNetworkQuantumSimulator.jl b/src/TensorNetworkQuantumSimulator.jl index 5825b4c..f3621d1 100644 --- a/src/TensorNetworkQuantumSimulator.jl +++ b/src/TensorNetworkQuantumSimulator.jl @@ -1,9 +1,8 @@ module TensorNetworkQuantumSimulator -# Legacy `ITensors` / `ITensorMPS` API, provided over the ITensorBase backend. Included -# before `imports.jl`, which aliases it as `ITensors` and imports the legacy names. -include("ITensorsITensorBaseCompat/ITensorsITensorBaseCompat.jl") +include("itensors.jl") +include("Ops.jl") include("imports.jl") include("siteinds.jl") diff --git a/src/TensorNetworks/abstracttensornetwork.jl b/src/TensorNetworks/abstracttensornetwork.jl index 444f66f..6df43cf 100644 --- a/src/TensorNetworks/abstracttensornetwork.jl +++ b/src/TensorNetworks/abstracttensornetwork.jl @@ -1,5 +1,4 @@ using Graphs: Graphs, has_vertex -import .ITensorsITensorBaseCompat as ITensors using NamedGraphs: NamedGraphs using Adapt @@ -21,7 +20,7 @@ NamedGraphs.edgetype(tn::AbstractTensorNetwork) = NamedGraphs.edgetype(graph(tn) NamedGraphs.vertextype(tn::AbstractTensorNetwork) = NamedGraphs.vertextype(graph(tn)) NamedGraphs.steiner_tree(tn::AbstractTensorNetwork, vs) = NamedGraphs.steiner_tree(graph(tn), vs) -virtualinds(tn::AbstractTensorNetwork, e::NamedEdge) = ITensors.commoninds(tn[src(e)], tn[dst(e)]) +virtualinds(tn::AbstractTensorNetwork, e::NamedEdge) = commoninds(tn[src(e)], tn[dst(e)]) virtualind(tn::AbstractTensorNetwork, e::NamedEdge) = only(virtualinds(tn, e)) function maxvirtualdim(tn::AbstractTensorNetwork) @@ -31,7 +30,7 @@ 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) +function uniqueinds(tn::AbstractTensorNetwork, v) tv_inds = Index[i for i in inds(tn[v])] vns = neighbors(tn, v) isempty(vns) && return tv_inds @@ -50,12 +49,12 @@ function Base.setindex!(tn::AbstractTensorNetwork, value::ITensor, vertex) return tn end -function ITensors.scalartype(tn::AbstractTensorNetwork) +function scalartype(tn::AbstractTensorNetwork) return mapreduce(v -> scalartype(tn[v]), promote_type, vertices(tn)) end -function ITensors.datatype(tn::AbstractTensorNetwork) - return mapreduce(v -> ITensors.datatype(tn[v]), promote_type, vertices(tn)) +function datatype(tn::AbstractTensorNetwork) + return mapreduce(v -> datatype(tn[v]), promote_type, vertices(tn)) end function map_tensors!(f::Function, tn::AbstractTensorNetwork) @@ -77,7 +76,7 @@ end function insert_virtualinds!(tn::AbstractTensorNetwork; bond_dimension::Integer = 1) dtype = datatype(tn) for e in edges(tn) - if isempty(ITensors.commoninds(tn[src(e)], tn[dst(e)])) + if isempty(commoninds(tn[src(e)], tn[dst(e)])) l = Index(bond_dimension) p = adapt(dtype)(onehot(l => 1)) setindex_preserve!(tn, tn[src(e)] * p, src(e)) @@ -94,10 +93,10 @@ end function map_virtualinds!(f::Function, tn::AbstractTensorNetwork) for e in edges(tn) - vinds = ITensors.commoninds(tn[src(e)], tn[dst(e)]) - vinds_sim = f(vinds) - setindex_preserve!(tn, ITensors.replaceinds(tn[src(e)], vinds, vinds_sim), src(e)) - setindex_preserve!(tn, ITensors.replaceinds(tn[dst(e)], vinds, vinds_sim), dst(e)) + vinds = commoninds(tn[src(e)], tn[dst(e)]) + vinds_sim = f.(vinds) + setindex_preserve!(tn, replaceinds(tn[src(e)], vinds, vinds_sim), src(e)) + setindex_preserve!(tn, replaceinds(tn[dst(e)], vinds, vinds_sim), dst(e)) end return tn end @@ -139,7 +138,7 @@ function add(tn1::AbstractTensorNetwork, tn2::AbstractTensorNetwork) tn12v_linkinds = Index[new_edge_indices[e] for e in es_v] setindex_preserve!( - tn12, ITensors.directsum( + tn12, directsum( tn12v_linkinds, tn1[v] => Tuple(tn1v_linkinds), tn2[v] => Tuple(tn2v_linkinds) diff --git a/src/TensorNetworks/tensornetwork.jl b/src/TensorNetworks/tensornetwork.jl index ceff6ef..ee9f4fe 100644 --- a/src/TensorNetworks/tensornetwork.jl +++ b/src/TensorNetworks/tensornetwork.jl @@ -1,7 +1,5 @@ using Dictionaries: Dictionary using Graphs: Graphs -import .ITensorsITensorBaseCompat as ITensors -using .ITensorsITensorBaseCompat: ITensor using NamedGraphs: NamedGraphs, add_edge!, incident_edges using NamedGraphs.GraphsExtensions: rem_edges! using Adapt @@ -68,7 +66,7 @@ function add_tensor!(tn::TensorNetwork, tensor::ITensor, v) end function default_message(tn::TensorNetwork, edge::NamedEdge) - return adapt_like(tn, denseblocks(delta(virtualinds(tn, edge)))) + return adapt_like(tn, delta(virtualinds(tn, edge))) end function bp_factors(tn::TensorNetwork, vertex) diff --git a/src/TensorNetworks/tensornetworkstate.jl b/src/TensorNetworks/tensornetworkstate.jl index 87cfe0a..b76780d 100644 --- a/src/TensorNetworks/tensornetworkstate.jl +++ b/src/TensorNetworks/tensornetworkstate.jl @@ -1,4 +1,3 @@ -using .ITensorsITensorBaseCompat: nameisdisjoint, namesetdiff, random_itensor """ TensorNetworkState{V} <: AbstractTensorNetwork{V} @@ -82,7 +81,7 @@ function norm_factors(tns::TensorNetworkState, verts::Vector; op_strings::Functi tnv_dag = replaceinds(tnv_dag, prime.(sinds), sinds) append!(factors, ITensor[tnv, tnv_dag]) else - op = adapt_like(tnv, ITensors.op(op_strings(v), only(sinds))) + op = adapt_like(tnv, Ops.op(op_strings(v), only(sinds))) append!(factors, ITensor[tnv, tnv_dag, op]) end end @@ -97,7 +96,7 @@ bp_factors(tns::TensorNetworkState, v) = norm_factors(tns, v) # links give a graded message; the legacy `delta` filled a dense diagonal). function default_message(tns::TensorNetworkState, edge::AbstractEdge) linds = virtualinds(tns, edge) - cod, dom = Tuple(linds), Tuple(prime(dag(linds))) + cod, dom = Tuple(linds), Tuple(prime.(dag(linds))) return adapt_like(tns, one(zeros(scalartype(tns), cod..., dom...), cod, dom)) end @@ -171,9 +170,9 @@ function tensornetworkstate(eltype, f::Function, g::AbstractGraph, siteinds::Dic for v in vs tnv = f(v) if tnv isa String - set!(tensors, v, ITensors.adapt_scalartype(eltype)(ITensors.state(f(v), only(siteinds[v])))) + set!(tensors, v, adapt_scalartype(eltype)(Ops.state(f(v), only(siteinds[v])))) elseif tnv isa Vector{<:Number} - set!(tensors, v, ITensors.adapt_scalartype(eltype)(ITensors.state(f(v), only(siteinds[v])))) + set!(tensors, v, adapt_scalartype(eltype)(Ops.state(f(v), only(siteinds[v])))) else error("Unrecognized local state constructor. Currently supported: Strings and Vectors.") end diff --git a/src/TensorNetworks/tensornetworkstate_constructors.jl b/src/TensorNetworks/tensornetworkstate_constructors.jl index c1a5a4a..1d60f3d 100644 --- a/src/TensorNetworks/tensornetworkstate_constructors.jl +++ b/src/TensorNetworks/tensornetworkstate_constructors.jl @@ -19,7 +19,7 @@ first half being the "ket" indices and the second half the "bra" indices; index index `n/2 + i`. """ function identity_tensornetworkstate(eltype, g::NamedGraph, s::Dictionary = siteinds("S=1/2", g; inds_per_site = 2)) - links = Dictionary(edges(g), [ITensors.settags(Index(1), "e$(src(e))_$(dst(e))") for e in edges(g)]) + links = Dictionary(edges(g), [settags(Index(1), "e$(src(e))_$(dst(e))") for e in edges(g)]) links = merge(links, Dictionary(reverse.(edges(g)), [links[e] for e in edges(g)])) ts = Dictionary{vertextype(g), ITensor}() @@ -27,7 +27,7 @@ function identity_tensornetworkstate(eltype, g::NamedGraph, s::Dictionary = site es = incident_edges(g, v; dir = :in) ninds = length(s[v]) ninds % 2 != 0 && error("Odd number of siteinds on vertex $v - don't know how to partition into rows and column") - t = ITensors.delta(eltype, [links[e] for e in es]) + t = delta(eltype, [links[e] for e in es]) if ninds > 0 row_inds, col_inds = s[v][1:(ninds÷2)], s[v][((ninds÷2)+1):ninds] id = identity_tensor(eltype, row_inds, col_inds) @@ -78,11 +78,11 @@ function toriccode_groundstate(n::Int, s::Dictionary = siteinds("S=1/2", named_g west_index = e_dict[NamedEdge(v => (v[1], mod1(v[2]-1, n)))] if iseven(sum(v)) - state = state + (ITensors.onehot(north_index => 1) * ITensors.onehot(east_index => 1) + ITensors.onehot(north_index => 2) * ITensors.onehot(east_index => 2)) * (ITensors.onehot(south_index => 1) * ITensors.onehot(west_index => 1) + ITensors.onehot(south_index => 2) * ITensors.onehot(west_index => 2)) * ITensors.onehot(sv => 1) - state = state + (ITensors.onehot(north_index => 1) * ITensors.onehot(east_index => 1) - ITensors.onehot(north_index => 2) * ITensors.onehot(east_index => 2)) * (ITensors.onehot(south_index => 1) * ITensors.onehot(west_index => 1) - ITensors.onehot(south_index => 2) * ITensors.onehot(west_index => 2)) * ITensors.onehot(sv => 2) + state = state + (onehot(north_index => 1) * onehot(east_index => 1) + onehot(north_index => 2) * onehot(east_index => 2)) * (onehot(south_index => 1) * onehot(west_index => 1) + onehot(south_index => 2) * onehot(west_index => 2)) * onehot(sv => 1) + state = state + (onehot(north_index => 1) * onehot(east_index => 1) - onehot(north_index => 2) * onehot(east_index => 2)) * (onehot(south_index => 1) * onehot(west_index => 1) - onehot(south_index => 2) * onehot(west_index => 2)) * onehot(sv => 2) else - state = state + (ITensors.onehot(north_index => 1) * ITensors.onehot(west_index => 1) + ITensors.onehot(north_index => 2) * ITensors.onehot(west_index => 2)) * (ITensors.onehot(south_index => 1) * ITensors.onehot(east_index => 1) + ITensors.onehot(south_index => 2) * ITensors.onehot(east_index => 2)) * ITensors.onehot(sv => 1) - state = state + (ITensors.onehot(north_index => 1) * ITensors.onehot(west_index => 1) - ITensors.onehot(north_index => 2) * ITensors.onehot(west_index => 2)) * (ITensors.onehot(south_index => 1) * ITensors.onehot(east_index => 1) - ITensors.onehot(south_index => 2) * ITensors.onehot(east_index => 2)) * ITensors.onehot(sv => 2) + state = state + (onehot(north_index => 1) * onehot(west_index => 1) + onehot(north_index => 2) * onehot(west_index => 2)) * (onehot(south_index => 1) * onehot(east_index => 1) + onehot(south_index => 2) * onehot(east_index => 2)) * onehot(sv => 1) + state = state + (onehot(north_index => 1) * onehot(west_index => 1) - onehot(north_index => 2) * onehot(west_index => 2)) * (onehot(south_index => 1) * onehot(east_index => 1) - onehot(south_index => 2) * onehot(east_index => 2)) * onehot(sv => 2) end set!(tensors, v, state) end @@ -111,7 +111,7 @@ Returns a `TensorNetwork` (not a `TensorNetworkState`); contract it to obtain ``Z(β)``. """ function ising_partitionfunction(g::NamedGraph, β::Real; Js::Dictionary = Dictionary(edges(g), [1.0 for e in edges(g)])) - links = Dictionary(edges(g), [ITensors.settags(Index(2), "e$(src(e))_$(dst(e))") for e in edges(g)]) + links = Dictionary(edges(g), [settags(Index(2), "e$(src(e))_$(dst(e))") for e in edges(g)]) links = merge(links, Dictionary(reverse.(edges(g)), [links[e] for e in edges(g)])) # symmetric sqrt of Boltzmann matrix W = exp(β σσ') @@ -133,7 +133,7 @@ function ising_partitionfunction(g::NamedGraph, β::Real; Js::Dictionary = Dicti ts = Dictionary{vertextype(g), ITensor}() for v in vertices(g) es = incident_edges(g, v; dir = :in) - t = ITensors.delta([links[e] for e in es]) + t = delta([links[e] for e in es]) for e in es t = noprime(ITensor(ComplexF64, sqrt_Ws[e], links[e], prime(links[e]))*t) end diff --git a/src/contract.jl b/src/contract.jl index 058bf66..cf72ab1 100644 --- a/src/contract.jl +++ b/src/contract.jl @@ -1,17 +1,17 @@ -function ITensors.contract(alg::Algorithm"exact", tn::AbstractTensorNetwork; contraction_sequence_kwargs = (; alg = "omeinsum", optimizer = GreedyMethod())) +function contract(alg::Algorithm"exact", tn::AbstractTensorNetwork; contraction_sequence_kwargs = (; alg = "omeinsum", optimizer = GreedyMethod())) tn_tensors = [tn[v] for v in vertices(tn)] seq = contraction_sequence(tn_tensors; contraction_sequence_kwargs...) - return scalar(ITensors.contract(tn_tensors; sequence = seq)) + return scalar(contract(tn_tensors; sequence = seq)) end -function ITensors.contract(alg::Algorithm"bp", tn::TensorNetwork; bp_update_kwargs = default_bp_update_kwargs(tn)) +function contract(alg::Algorithm"bp", tn::TensorNetwork; bp_update_kwargs = default_bp_update_kwargs(tn)) return partitionfunction(update(BeliefPropagationCache(tn); bp_update_kwargs...)) end -function ITensors.contract(alg::Algorithm"boundarymps", tn::TensorNetwork; mps_bond_dimension::Integer, bmps_update_kwargs = default_bmps_update_kwargs(tn)) +function contract(alg::Algorithm"boundarymps", tn::TensorNetwork; mps_bond_dimension::Integer, bmps_update_kwargs = default_bmps_update_kwargs(tn)) return partitionfunction(update(BoundaryMPSCache(tn, mps_bond_dimension); bmps_update_kwargs...)) end -function ITensors.contract(tn::AbstractTensorNetwork; alg = "exact", kwargs...) - return ITensors.contract(Algorithm(alg), tn; kwargs...) +function contract(tn::AbstractTensorNetwork; alg = "exact", kwargs...) + return contract(Algorithm(alg), tn; kwargs...) end diff --git a/src/contraction_sequences.jl b/src/contraction_sequences.jl index bcbe3c6..73b572d 100644 --- a/src/contraction_sequences.jl +++ b/src/contraction_sequences.jl @@ -1,4 +1,3 @@ -using .ITensorsITensorBaseCompat: Index, ITensor, @Algorithm_str, inds, noncommoninds, dim using TensorOperations: TensorOperations, optimaltree using OMEinsumContractionOrders: OMEinsumContractionOrders, optimize_code, EinCode, NestedEinsum, TreeSA, GreedyMethod, SABipartite, Treewidth, ExactTreewidth, HyperND diff --git a/src/entanglement.jl b/src/entanglement.jl index 0a93c0c..4276c86 100644 --- a/src/entanglement.jl +++ b/src/entanglement.jl @@ -32,7 +32,7 @@ function matricize(a::ITensor, row_inds = filter(i -> plev(i) ==0, inds(a))) col_inds = prime.(row_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)) + return array(matricize(a, Tuple(row_inds) => row_name, Tuple(col_inds) => col_name)) end """ diff --git a/src/expect.jl b/src/expect.jl index 003f9f7..02fe6f1 100644 --- a/src/expect.jl +++ b/src/expect.jl @@ -4,7 +4,6 @@ function expect( observables::Vector{<:Tuple}; contraction_sequence_kwargs = (; alg = "omeinsum", optimizer = GreedyMethod()) ) - ITensors.disable_warn_order() denom = norm_sqr(alg, ψ; contraction_sequence_kwargs) out = Number[] diff --git a/src/imports.jl b/src/imports.jl index c9c85c6..e404c21 100644 --- a/src/imports.jl +++ b/src/imports.jl @@ -46,26 +46,12 @@ using NamedGraphs.NamedGraphGenerators: named_grid, named_hexagonal_lattice_grap using TensorOperations -# Legacy `ITensors` / `ITensorMPS` API, republished over the ITensorBase backend by the -# `ITensorsITensorBaseCompat` submodule (included before this file). It is aliased as -# `ITensors` so existing `ITensors.foo` call sites and `function ITensors.foo` -# extensions keep working unchanged, and its legacy names are imported for unqualified -# use. Each source file keeps its own `import`/`using` of this module exactly where it -# 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, cat_inds, - sim, dag, prime, noprime, replaceind, replaceinds, dim, swapind, - 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, - directsum, disable_warn_order, - Algorithm, @Algorithm_str, - hastags, - state, op, OpName, SiteType, @OpName_str, @SiteType_str -using ITensorBase: ITensorBase, Index, ITensor, name, plev, tags, unnamed +# TNQS-owned operator / named-state system. `op`/`state` are called qualified (`Ops.op`, +# `Ops.state`) so `state` does not clash with the unrelated `ITensorBase.state`, and gates +# are registered by extending `Ops.op`. The types and string macros are imported for +# unqualified use (gate definitions dispatch on bare `OpName"…"` / `SiteType"…"`). +using .Ops: OpName, SiteType, @OpName_str, @SiteType_str +using ITensorBase: ITensorBase, Index, ITensor, name, noprime, plev, prime, tags, unnamed using TensorAlgebra: trivialrange using TensorAlgebra.MatrixAlgebra: sqrth_invsqrth_safe, sqrth_safe using MatrixAlgebraKit: project_hermitian diff --git a/src/inner.jl b/src/inner.jl index c3c641c..c1d581d 100644 --- a/src/inner.jl +++ b/src/inner.jl @@ -31,22 +31,22 @@ Compute the inner product ⟨ψ|ϕ⟩ between two `TensorNetworkState`s using th ϕ = random_tensornetworkstate(ComplexF32, g, s; bond_dimension = 4) # Exact inner product - ip_exact = ITensors.inner(ψ, ϕ; alg = "exact") + ip_exact = inner(ψ, ϕ; alg = "exact") # Belief propagation inner product - ip_bp = ITensors.inner(ψ, ϕ; alg = "bp") + ip_bp = inner(ψ, ϕ; alg = "bp") # Boundary MPS inner product with bond dimension 10 - ip_bmps = ITensors.inner(ψ, ϕ; alg = "boundarymps", mps_bond_dimension = 10) + ip_bmps = inner(ψ, ϕ; alg = "boundarymps", mps_bond_dimension = 10) ``` """ -function ITensors.inner(ψ::TensorNetworkState, ϕ::TensorNetworkState; alg, kwargs...) +function inner(ψ::TensorNetworkState, ϕ::TensorNetworkState; alg, kwargs...) algorithm_check(ψ, "inner", alg) algorithm_check(ϕ, "inner", alg) return inner(Algorithm(alg), ψ, ϕ; kwargs...) end -function ITensors.inner( +function inner( alg::Algorithm"exact", blf::BilinearForm; contraction_sequence_kwargs = (; alg = "omeinsum", optimizer = GreedyMethod()) ) @@ -55,26 +55,26 @@ function ITensors.inner( return scalar(contract(blf_tensors; sequence = seq)) end -function ITensors.inner(alg::Algorithm, cache::AbstractBeliefPropagationCache; max_configuration_size = nothing) +function inner(alg::Algorithm, cache::AbstractBeliefPropagationCache; max_configuration_size = nothing) tn = network(cache) z = cache_partitionfunction(alg, cache; max_configuration_size) tn isa BilinearForm && return z return state_error("BilinearForm") end -function ITensors.inner(alg::Union{Algorithm"bp", Algorithm"loopcorrections"}, ψ::TensorNetworkState, ϕ::TensorNetworkState; cache_update_kwargs = (;), kwargs...) +function inner(alg::Union{Algorithm"bp", Algorithm"loopcorrections"}, ψ::TensorNetworkState, ϕ::TensorNetworkState; cache_update_kwargs = (;), kwargs...) ψϕ_bpc = BeliefPropagationCache(BilinearForm(ψ, ϕ)) ψϕ_bpc = update(ψϕ_bpc; cache_update_kwargs...) return inner(alg, ψϕ_bpc; kwargs...) end -function ITensors.inner(alg::Algorithm"boundarymps", ψ::TensorNetworkState, ϕ::TensorNetworkState; mps_bond_dimension::Integer, partition_by = "row", cache_update_kwargs = (;), kwargs...) +function inner(alg::Algorithm"boundarymps", ψ::TensorNetworkState, ϕ::TensorNetworkState; mps_bond_dimension::Integer, partition_by = "row", cache_update_kwargs = (;), kwargs...) ψϕ_bmps = BoundaryMPSCache(BilinearForm(ψ, ϕ), mps_bond_dimension; partition_by) cache_update_kwargs = with_default_maxiter(cache_update_kwargs, ψϕ_bmps) ψϕ_bmps = update(ψϕ_bmps; cache_update_kwargs...) return inner(alg, ψϕ_bmps; kwargs...) end -function ITensors.inner(alg::Algorithm"exact", ψ::TensorNetworkState, ϕ::TensorNetworkState) +function inner(alg::Algorithm"exact", ψ::TensorNetworkState, ϕ::TensorNetworkState) return inner(alg, BilinearForm(ψ, ϕ)) end diff --git a/src/ITensorsITensorBaseCompat/itensor.jl b/src/itensors.jl similarity index 79% rename from src/ITensorsITensorBaseCompat/itensor.jl rename to src/itensors.jl index 8c8507f..c3ee4fb 100644 --- a/src/ITensorsITensorBaseCompat/itensor.jl +++ b/src/itensors.jl @@ -1,27 +1,53 @@ -# Compatibility shims bridging the legacy `ITensors.jl` API that TNQS was written -# against to the next-gen `ITensorBase.jl` backend. -# -# Strategy (see ITensorDevelopmentPlans api_migration_map.md): -# - Names below are thin wrappers over ITensorBase / TensorAlgebra / MatrixAlgebraKit. -# - 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. +# TNQS-owned tensor utilities and the legacy `ITensors.jl`-style API TNQS was written +# against, implemented over the next-gen `ITensorBase` / `TensorAlgebra` / `MatrixAlgebraKit` +# stack. ITensorBase keeps most of this API internal (unexported), so we reach for the +# qualified names and define the legacy spellings here as functionality TNQS owns. +import Base: truncate +import ITensorBase: scalartype import MatrixAlgebraKit as MAK +import TensorAlgebra: matricize using Adapt: Adapt -using ITensorBase: ITensorBase, AbstractITensor, ITensor, Index, NamedUnitRange, name, - nameddims, plev, tags, unnamed +using ITensorBase: ITensorBase, AbstractITensor, ITensor, Index, NamedUnitRange, dimnames, + id, name, nameddims, noprime, plev, prime, replacedimnames, sim, tags, unnamed using LinearAlgebra: LinearAlgebra -using TensorAlgebra.MatrixAlgebra: MatrixAlgebra -using TensorAlgebra: TensorAlgebra +using TensorAlgebra: TensorAlgebra, project, tryproject -# 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). +# +# State-vector projection, shared by the `Ops` state constructors and by `onehot`. +# +# Project a raw vector as an `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`). +function project_aux(v::AbstractVector{<:Number}, i::Index) + length(v) == length(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 + +# One-hot vector along `i` at position `p` (legacy `onehot(i => p)`), through `project_aux`, +# 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 project_aux(v, i) +end +onehot(p::Pair{<:Index}) = onehot(Float64, p) + +# Legacy `inds(t; plev, tags)` took index-filtering keywords; `ITensorBase.inds` takes none. +# This forwards to it and applies the legacy filters. function inds(t::AbstractITensor; plev = nothing, tags = nothing) is = ITensorBase.inds(t) isnothing(plev) || (is = filter(i -> ITensorBase.plev(i) == plev, is)) @@ -29,13 +55,12 @@ function inds(t::AbstractITensor; plev = nothing, tags = nothing) return is end -# 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) +# `contract` / `inner` are TNQS operations whose base generics live here and are extended +# for TNQS types (tensor networks, caches) in the respective files. function inner end # 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 +# `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 @@ -105,13 +130,6 @@ noncommoninds(a, b) = namesymdiff(_compat_inds(a), _compat_inds(b)) dim(i::Index) = length(i) dim(is::Union{Tuple, AbstractVector}) = isempty(is) ? 1 : prod(length, is) -# -# Index operations. -# -# 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) - # Conjugate (legacy `dag`): `conj` the tensor, and on bare indices flip the sector arrows # on a graded axis. `conj(::Index)` is id-preserving on the dense backend, so `dag` there # is effectively the identity on indices, matching legacy behavior. @@ -119,21 +137,7 @@ dag(t::AbstractITensor) = conj(t) dag(i::Index) = conj(i) dag(is::Union{Tuple, AbstractVector}) = map(conj, is) -# `prime` / `noprime`: ITensorBase primes an `Index`; TNQS also primes whole ITensors -# (all dimnames). Owned here with a fallback to ITensorBase for the index case. -prime(x) = ITensorBase.prime(x) -function prime(t::AbstractITensor) - return nameddims(unnamed(t), map(ITensorBase.prime, ITensorBase.dimnames(t))) -end -prime(is::Union{Tuple, AbstractVector{<:Index}}) = map(prime, is) -noprime(x) = ITensorBase.noprime(x) -function noprime(t::AbstractITensor) - return nameddims(unnamed(t), map(ITensorBase.noprime, ITensorBase.dimnames(t))) -end -noprime(is::Union{Tuple, AbstractVector{<:Index}}) = map(noprime, is) - -# `replaceind` (singular) maps to a single-pair replacement, forwarding to the -# compat-owned `replaceinds` (below). +# `replaceind` (singular) maps to a single-pair replacement, forwarding to `replaceinds`. replaceind(t, p::Pair) = replaceinds(t, p) replaceind(t, from::Index, to::Index) = replaceinds(t, from => to) @@ -145,8 +149,9 @@ 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` took collection arguments (`replaceinds(t, [i,k], [j,l])`, `... => ...`); -# ITensorBase provides only the pair-splat form, so this compat handles the collection forms. +# ITensorBase provides only the pair-splat form, so this handles the collection forms. # # 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` @@ -156,10 +161,7 @@ cat_inds(xs...) = reduce(vcat, map(_as_index_vec, xs)) # 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( - t, - map(p -> name(first(p)) => name(last(p)), pairs)... - ) + return replacedimnames(t, map(p -> name(first(p)) => name(last(p)), pairs)...) end function replaceinds(t, from::_IndexColl, to::_IndexColl) return replaceinds(t, map(=>, from, to)...) @@ -168,9 +170,6 @@ function replaceinds(t::AbstractITensor, p::Pair{<:_IndexColl, <:_IndexColl}) return replaceinds(t, first(p), last(p)) end -# -# ITensor construction. -# # 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 @@ -198,8 +197,8 @@ random_itensor(is::Union{Tuple, AbstractVector}) = random_itensor(Float64, is) scalar(t::AbstractITensor) = t[] # 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). +# `ITensorNetworkGenerators/delta_network.jl`. A graded/sector-aware `delta` is a stack gap +# (tracked separately). diaglength(a::AbstractArray) = minimum(size(a)) function diagstride(a::AbstractArray) s = 1 @@ -243,15 +242,6 @@ 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. @@ -262,32 +252,22 @@ delta(is::Tuple) = delta(Float64, is) delta(is::Index...) = delta(Float64, is) delta(is::AbstractVector{<:Index}) = delta(Float64, Tuple(is)) -# 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) +# Trace over prime pairs (legacy ITensors `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. Named `itensor_tr` (not a `LinearAlgebra.tr` method, which would be piracy on +# `AbstractITensor`) so it stays distinct from `tr` on plain matrices, which TNQS also calls. +function itensor_tr(t::AbstractITensor) 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)`), 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 project_aux(v, i) -end -onehot(p::Pair{<:Index}) = onehot(Float64, p) - # # Factorizations. These map onto MatrixAlgebraKit's named-tensor methods # (`f(a, codomain, domain)`). The legacy return shapes differ from MAK's, so the # heavy factorization callsites (simple_update / full_update / symmetric_gauge) are -# translated directly rather than fully wrapped here; these aliases cover the simple -# uses. See api_migration_map.md. +# translated directly rather than fully wrapped here; these aliases cover the simple uses. # const svd_trunc = MAK.svd_trunc @@ -324,8 +304,7 @@ end # `eigh_full` returns `D` over two independent fresh indices; we rename its `U`-disjoint # index to `prime(u)` so `D` becomes the legacy prime-paired diagonal. The legacy # 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). +# `trunc=(; ...)` spec; the current callsites pass `cutoff = nothing` (full decomposition). function _eigh( m::AbstractITensor, codomain, @@ -439,24 +418,14 @@ function factorize( return L, R 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 imported from ITensorBase (extended for +# TNQS types elsewhere). `datatype` is the underlying storage array type (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. # -# Storage / element type accessors. -# -# `scalartype` is re-exported above. `datatype` is the underlying storage array type -# (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) = convert(Array, unnamed(T)) data(T::AbstractITensor) = unnamed(T) @@ -479,10 +448,6 @@ end # `swapind`: swap two indices (legacy convenience over `replaceinds`). swapind(T::AbstractITensor, i::Index, j::Index) = replaceinds(T, i => j, j => i) -# Dense no-ops. Legacy QN-storage helpers; on the dense next-gen backend the tensor -# is already dense, so these are identities. (Graded/QN path is a stack gap.) -denseblocks(T::AbstractITensor) = T -dense(T::AbstractITensor) = T # Whether a tensor carries quantum-number (graded) block structure: true when any of its # indices is graded. `loopcorrection` branches on this to pick a contraction-order # algorithm. A graded axis differs from its conjugate (conjugation flips the sector @@ -492,9 +457,6 @@ hasqns(i::Index) = conj(unnamed(i)) != unnamed(i) hasqns(t::AbstractITensor) = any(hasqns, inds(t)) 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 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), ...) @@ -521,14 +483,12 @@ function directsum(out_inds, pairs::Pair...) return out[target...] end -# No index-order warnings in the next-gen stack (legacy `disable_warn_order`). -disable_warn_order(args...) = nothing - # # Algorithm dispatch tag (legacy `Algorithm` / `@Algorithm_str`). `Algorithm"exact"` # is the type `Algorithm{:exact}` (usable in `::Algorithm"exact"` signatures); # `Algorithm("exact"; kwargs...)` constructs an instance carrying keyword options. # Faithful minimal copy of the ITensors/NDTensors helper. +# struct Algorithm{Alg, Kwargs <: NamedTuple} kwargs::Kwargs end @@ -554,6 +514,7 @@ end # stores tags as a `Dict{String, String}`. For the legacy flat-tag usage TNQS has # (site-type labels, link names), we store each token as a keyed tag with empty # value, and `hastags` checks membership. (A fuller tag-compat story is a follow-up.) +# function settags(i::Index, tagstr::AbstractString) for t in split(tagstr, ",") s = String(strip(t)) @@ -574,13 +535,3 @@ function hastags(i::Index, tagstr::AbstractString) haskey(tags(i), String(strip(t))) for t in split(tagstr, ",") if !isempty(strip(t)) ) end - -# TODO (small inline residue — can't be a drop-in shim): -# - `contract` / `inner` / `truncate`: TNQS *extends* these (method definitions), -# so the call sites drop the `ITensors.` qualifier to extend the generics this -# module owns rather than ITensors'. (`truncate` clashes with `Base.truncate`.) -# -# TODO (stack gaps — own roadmap items, not solvable in this file): -# - op / @OpName_str / @SiteType_str (operator/site system) -# - MPS / MPO / boundary-MPS (no next-gen MPS layer yet) -# - hasqns / QN-aware storage (next-gen symmetry via GradedArrays) diff --git a/src/rdm.jl b/src/rdm.jl index 43b8ac1..0ec72bd 100644 --- a/src/rdm.jl +++ b/src/rdm.jl @@ -1,5 +1,5 @@ function normalize_rdm(ρ::ITensor) - return ρ / ITensors.tr(ρ) + return ρ / itensor_tr(ρ) end """ @@ -32,7 +32,6 @@ function reduced_density_matrix( contraction_sequence_kwargs = (; alg = "omeinsum", optimizer = GreedyMethod()), normalize = true ) - ITensors.disable_warn_order() op_string_f = v -> v ∈ verts ? "ρ" : "I" ρ_tensors = norm_factors(ψ, collect(vertices(ψ)); op_strings = op_string_f) seq = contraction_sequence(ρ_tensors; contraction_sequence_kwargs...) diff --git a/src/sampling.jl b/src/sampling.jl index 635d447..9116f0d 100644 --- a/src/sampling.jl +++ b/src/sampling.jl @@ -23,11 +23,11 @@ function sample( ψ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) + ρ = contract(tensors; sequence = seq) - ρ_tr = ITensors.tr(ρ) + ρ_tr = itensor_tr(ρ) ρ *= inv(ρ_tr) - ρ_diag = collect(real.(diag(ITensors.array(ρ)))) + ρ_diag = collect(real.(diag(array(ρ)))) config = StatsBase.sample(1:length(ρ_diag), Weights(ρ_diag)) # config is 1,2,...,d, but we want 0,1...,d-1 for the sample itself set!(bit_string, v, config - 1) @@ -238,10 +238,10 @@ function sample_partition!( ts = [incoming_ms; [ψv, ψvdag]] seq = contraction_sequence(ts; alg = "optimal") ρ = contract(ts; sequence = seq) - ρ_tr = ITensors.tr(ρ) + ρ_tr = itensor_tr(ρ) push!(traces, ρ_tr) ρ *= inv(ρ_tr) - ρ_diag = collect(real.(diag(ITensors.array(ρ)))) + ρ_diag = collect(real.(diag(array(ρ)))) config = StatsBase.sample(1:length(ρ_diag), Weights(ρ_diag)) # config is 1,2,...,d, but we want 0,1...,d-1 for the sample itself set!(bit_string, v, config - 1) diff --git a/src/siteinds.jl b/src/siteinds.jl index 4713612..d08758e 100644 --- a/src/siteinds.jl +++ b/src/siteinds.jl @@ -6,7 +6,7 @@ end function siteinds(sitetype::String, g::AbstractGraph, sitedimension::Integer = site_dimension(sitetype); inds_per_site::Integer = 1) vs = collect(vertices(g)) - return Dictionary{vertextype(g), Vector{<:Index}}(vs, [Index[ITensors.settags(Index(sitedimension), site_tag(sitetype)) for i in 1:inds_per_site] for v in vs]) + return Dictionary{vertextype(g), Vector{<:Index}}(vs, [Index[settags(Index(sitedimension), site_tag(sitetype)) for i in 1:inds_per_site] for v in vs]) end function site_dimension(sitetype::String) diff --git a/src/symmetric_gauge.jl b/src/symmetric_gauge.jl index 7c2bda5..5a9c699 100644 --- a/src/symmetric_gauge.jl +++ b/src/symmetric_gauge.jl @@ -7,7 +7,7 @@ function symmetric_gauge!(bp_cache::BeliefPropagationCache; regularization = 10 edge_ind = commoninds(ψvsrc, ψvdst) edge_ind_p = prime.(edge_ind) - edge_ind_sim = sim(edge_ind) + edge_ind_sim = sim.(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 @@ -48,7 +48,7 @@ function symmetric_gauge!(bp_cache::BeliefPropagationCache; regularization = 10 ψvsrc = ψvsrc * sqrtS ψvdst = ψvdst * sqrtS - new_edge_ind = Index[ITensors.settags(only(v), tags(first(edge_ind)))] + new_edge_ind = Index[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) diff --git a/src/truncate.jl b/src/truncate.jl index 4cec666..c105fa1 100644 --- a/src/truncate.jl +++ b/src/truncate.jl @@ -9,7 +9,7 @@ function truncatable_edge(cache::AbstractBeliefPropagationCache, e::NamedEdge) return true end -function ITensors.truncate(bpc::BeliefPropagationCache; bp_update_kwargs = default_bp_update_kwargs(bpc), maxdim::Integer, cutoff = nothing, edge_color = true, normalize_tensors = true) +function truncate(bpc::BeliefPropagationCache; bp_update_kwargs = default_bp_update_kwargs(bpc), maxdim::Integer, cutoff = nothing, edge_color = true, normalize_tensors = true) bpc = copy(bpc) s = siteinds(network(bpc)) apply_kwargs = (; maxdim, cutoff, normalize_tensors) @@ -21,7 +21,7 @@ function ITensors.truncate(bpc::BeliefPropagationCache; bp_update_kwargs = defau for eg in edge_groups for e in eg if truncatable_edge(bpc, e) - g1, g2 = reduce(*, [ITensors.op("I", sv) for sv in s[src(e)] ]), reduce(*, [ITensors.op("I", sv) for sv in s[dst(e)] ]) + g1, g2 = reduce(*, [Ops.op("I", sv) for sv in s[src(e)] ]), reduce(*, [Ops.op("I", sv) for sv in s[dst(e)] ]) apply_gate!(adapt(dtype)(g1 * g2), bpc; v⃗ = [src(e), dst(e)], apply_kwargs) end end @@ -29,7 +29,7 @@ function ITensors.truncate(bpc::BeliefPropagationCache; bp_update_kwargs = defau end else for e in edges(bpc) - g1, g2 = reduce(*, [ITensors.op("I", sv) for sv in s[src(e)]]), reduce(*, [ITensors.op("I", sv) for sv in s[dst(e)]]) + g1, g2 = reduce(*, [Ops.op("I", sv) for sv in s[src(e)]]), reduce(*, [Ops.op("I", sv) for sv in s[dst(e)]]) apply_gate!(adapt(dtype)(g1 * g2), bpc; v⃗ = [src(e), dst(e)], apply_kwargs) bpc = update(bpc; bp_update_kwargs...) end @@ -37,7 +37,7 @@ function ITensors.truncate(bpc::BeliefPropagationCache; bp_update_kwargs = defau return bpc end -function ITensors.truncate(bmps_cache::BoundaryMPSCache; maxdim::Integer, cutoff = nothing, normalize_tensors = true) +function truncate(bmps_cache::BoundaryMPSCache; maxdim::Integer, cutoff = nothing, normalize_tensors = true) bmps_cache = copy(bmps_cache) s = siteinds(network(bmps_cache)) apply_kwargs = (; maxdim, cutoff) @@ -50,7 +50,7 @@ function ITensors.truncate(bmps_cache::BoundaryMPSCache; maxdim::Integer, cutoff !isempty(seq) && update_partition!(bmps_cache, seq) for e in reverse.(reverse(seq)) if truncatable_edge(bmps_cache, e) - g1, g2 = reduce(*, [ITensors.op("I", sv) for sv in s[src(e)]]), reduce(*, [ITensors.op("I", sv) for sv in s[dst(e)]]) + g1, g2 = reduce(*, [Ops.op("I", sv) for sv in s[src(e)]]), reduce(*, [Ops.op("I", sv) for sv in s[dst(e)]]) envs = incoming_messages(bmps_cache, [src(e), dst(e)]) ρv1, ρv2 = full_update(adapt(dtype)(g1 * g2), network(bmps_cache), [src(e), dst(e)]; envs, apply_kwargs...) if normalize_tensors @@ -71,14 +71,14 @@ function ITensors.truncate(bmps_cache::BoundaryMPSCache; maxdim::Integer, cutoff return bmps_cache end -function ITensors.truncate(alg::Algorithm"bp", tns::TensorNetworkState; kwargs...) +function truncate(alg::Algorithm"bp", tns::TensorNetworkState; kwargs...) bp_cache = BeliefPropagationCache(tns) bp_cache = update(bp_cache) bp_cache = truncate(bp_cache; kwargs...) return network(bp_cache) end -function ITensors.truncate(alg::Algorithm"boundarymps", tns::TensorNetworkState; mps_bond_dimension::Integer, gauge_state = true, kwargs...) +function truncate(alg::Algorithm"boundarymps", tns::TensorNetworkState; mps_bond_dimension::Integer, gauge_state = true, kwargs...) tns = copy(tns) bmps_cache = BoundaryMPSCache(tns, mps_bond_dimension; partition_by = "row", gauge_state) leaves = leaf_vertices(quotient_graph(supergraph(bmps_cache))) @@ -113,7 +113,7 @@ Truncate the virtual indices of tensors in the given `TensorNetworkState` using # Returns - The truncated `TensorNetworkState`. """ -function ITensors.truncate(tns::TensorNetworkState; alg = default_truncate_alg(tns), kwargs...) +function truncate(tns::TensorNetworkState; alg = default_truncate_alg(tns), kwargs...) algorithm_check(tns, "truncate", alg) return truncate(Algorithm(alg), tns; kwargs...) end diff --git a/src/utils.jl b/src/utils.jl index 1d39e0e..790f07b 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -71,20 +71,20 @@ end """ safe_eigen(m::ITensor, args...; kwargs...) - A wrapper around ITensors.eigen that ensures eigen computations are done in Float64/ComplexF64 precision on CPU for better numerical stability. + A wrapper around eigen that ensures eigen computations are done in Float64/ComplexF64 precision on CPU for better numerical stability. """ function safe_eigen(m::ITensor, args...; kwargs...) dtype = datatype(m) e = eltype(m) if e == ComplexF64 || e == Float64 - return ITensors.eigen(m, args...; kwargs...) + return eigen(m, args...; kwargs...) elseif e == Float32 m = adapt(Vector{Float64}, m) - D, U = ITensors.eigen(m, args...; kwargs...) + D, U = eigen(m, args...; kwargs...) return adapt(dtype)(D), adapt(dtype)(U) elseif e == ComplexF32 m = adapt(Vector{ComplexF64}, m) - D, U = ITensors.eigen(m, args...; kwargs...) + D, U = eigen(m, args...; kwargs...) return adapt(dtype)(D), adapt(dtype)(U) end end diff --git a/test/test_apply.jl b/test/test_apply.jl index 120685a..1bf95be 100644 --- a/test/test_apply.jl +++ b/test/test_apply.jl @@ -1,7 +1,6 @@ @eval module $(gensym()) -import TensorNetworkQuantumSimulator as ITensors using ITensorBase: Index -using TensorNetworkQuantumSimulator: datatype, op, @OpName_str, @SiteType_str +using TensorNetworkQuantumSimulator: datatype, Ops, OpName, SiteType, @OpName_str, @SiteType_str using Random using TensorNetworkQuantumSimulator using Test: @testset, @test, @test_throws @@ -59,7 +58,7 @@ end # Define a custom op: a Z-axis rotation under a non-built-in name. # (Same matrix as the built-in "Rz", under a new name, so we can verify # the registered gate dispatches correctly.) - ITensors.op(::ITensors.OpName"MyZRot", ::ITensors.SiteType"S=1/2"; θ::Number) = + Ops.op(::OpName"MyZRot", ::SiteType"S=1/2"; θ::Number) = exp(-im * (θ / 2) * [1 0; 0 -1]) # Register the dispatch info: name "MyZRot" takes a single keyword `θ`. diff --git a/test/test_constructors.jl b/test/test_constructors.jl index 9d09181..e04adc8 100644 --- a/test/test_constructors.jl +++ b/test/test_constructors.jl @@ -3,11 +3,9 @@ using Dictionaries: Dictionary using ITensorBase: Index, inds using Random using TensorNetworkQuantumSimulator -# `random_itensor` / `contract` come from TNQS; alias TNQS as `ITensors` so the -# `ITensors.`-qualified legacy calls resolve. `dag` and `prime` are TNQS-owned compat -# (`prime` extends ITensorBase's to whole tensors / index collections, which -# `map_virtualinds` needs); `inds` is ITensorBase's. -import TensorNetworkQuantumSimulator as ITensors +# `random_itensor` / `contract` are TNQS-owned; reach them via the `TNQS` alias. `dag` is +# TNQS-owned; `prime` and `inds` are ITensorBase's. +const TNQS = TensorNetworkQuantumSimulator using TensorNetworkQuantumSimulator: dag, prime using Test: @testset, @test, @test_throws @@ -17,7 +15,7 @@ using Test: @testset, @test, @test_throws #TensorNetwork construction from tensors i, j, k, l = Index(2), Index(2), Index(2), Index(2) - A, B, C, D = ITensors.random_itensor(i, j), ITensors.random_itensor(j, k), ITensors.random_itensor(k, l), ITensors.random_itensor(l, i) + A, B, C, D = TNQS.random_itensor(i, j), TNQS.random_itensor(j, k), TNQS.random_itensor(k, l), TNQS.random_itensor(l, i) t = TensorNetwork([A, B, C, D]) @test t isa TensorNetwork @test scalartype(t) == eltype(A) @@ -38,7 +36,7 @@ using Test: @testset, @test, @test_throws ψdag = map_virtualinds(prime, map_tensors(dag, ψ)) @test ψdag isa TensorNetwork - @test ITensors.contract(ψdag; alg = "exact") ≈ conj(ITensors.contract(ψ; alg = "exact")) + @test TNQS.contract(ψdag; alg = "exact") ≈ conj(TNQS.contract(ψ; alg = "exact")) v = first(vertices(g)) rem_vertex!(ψ, v) diff --git a/test/test_contraction_sequences.jl b/test/test_contraction_sequences.jl index da0e44b..2de5374 100644 --- a/test/test_contraction_sequences.jl +++ b/test/test_contraction_sequences.jl @@ -4,7 +4,6 @@ using Random using TensorNetworkQuantumSimulator const TNQS = TensorNetworkQuantumSimulator # `random_itensor`/`contract`/`scalar` come from TNQS's compat layer (see test_constructors). -import TensorNetworkQuantumSimulator as ITensors using TensorNetworkQuantumSimulator: scalar using OMEinsumContractionOrders: NestedEinsum, EinCode, getixsv, getiyv using Test: @testset, @test @@ -21,8 +20,8 @@ collect_leaves!(acc, x) = (for y in x; collect_leaves!(acc, y); end; acc) # 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) + A = TNQS.random_itensor(i, j) + B = TNQS.random_itensor(j, k) code, size_dict = TNQS.to_eincode([A, B]) @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) @@ -50,21 +49,21 @@ collect_leaves!(acc, x) = (for y in x; collect_leaves!(acc, y); end; acc) # --- the sequence the backend returns is a *correct* contraction: executing it gives the # same scalar as the independent `optimal` backend. - ref = scalar(ITensors.contract(tensors; sequence = TNQS.contraction_sequence(tensors; alg = "optimal"))) + ref = scalar(TNQS.contract(tensors; sequence = TNQS.contraction_sequence(tensors; alg = "optimal"))) for optimizer in (GreedyMethod(), TreeSA()) seq = TNQS.contraction_sequence(tensors; alg = "omeinsum", optimizer) - @test scalar(ITensors.contract(tensors; sequence = seq)) ≈ ref + @test scalar(TNQS.contract(tensors; sequence = seq)) ≈ ref end # --- open network: result is a tensor with dangling indices (iy non-empty). p, q, r, s, t = Index(2), Index(3), Index(2), Index(3), Index(2) - X = ITensors.random_itensor(p, q) - Y = ITensors.random_itensor(q, r, s) - Z = ITensors.random_itensor(s, t) + X = TNQS.random_itensor(p, q) + Y = TNQS.random_itensor(q, r, s) + Z = TNQS.random_itensor(s, t) open_tensors = [X, Y, Z] # open indices: p, r, t seq_open = TNQS.contraction_sequence(open_tensors; alg = "omeinsum", optimizer = GreedyMethod()) @test sort(collect_leaves!(Int[], seq_open)) == [1, 2, 3] - @test ITensors.contract(open_tensors; sequence = seq_open) ≈ - ITensors.contract(open_tensors; sequence = TNQS.contraction_sequence(open_tensors; alg = "optimal")) + @test TNQS.contract(open_tensors; sequence = seq_open) ≈ + TNQS.contract(open_tensors; sequence = TNQS.contraction_sequence(open_tensors; alg = "optimal")) end end From 1756c27a1966f779a8b661b2593b19d952f9bfbf Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Wed, 8 Jul 2026 20:24:42 -0400 Subject: [PATCH 02/21] Replace the random_itensor shim with randn randn over Index objects already covers what the legacy random_itensor wrapper did, so drop the four-line shim and call randn directly at the construction sites and in the tests. Also updates the states.md example off the removed Index tag-string constructor. Co-Authored-By: Claude Opus 4.8 (1M context) --- docs/src/states.md | 6 +++--- src/TensorNetworks/tensornetwork.jl | 2 +- src/TensorNetworks/tensornetworkstate.jl | 2 +- src/itensors.jl | 6 ------ test/test_constructors.jl | 6 +++--- test/test_contraction_sequences.jl | 12 ++++++------ 6 files changed, 14 insertions(+), 20 deletions(-) diff --git a/docs/src/states.md b/docs/src/states.md index dfd647a..d74c0df 100644 --- a/docs/src/states.md +++ b/docs/src/states.md @@ -8,10 +8,10 @@ A `TensorNetwork` is the simpler of the two types: a collection of ITensors livi ```julia using TensorNetworkQuantumSimulator -using ITensors: Index, random_itensor +using ITensorBase: Index -i, j = Index(2, "i"), Index(2, "j") -t_a, t_b, t_c = random_itensor(i), random_itensor(i, j), random_itensor(j) +i, j = Index(2), Index(2) +t_a, t_b, t_c = randn(i), randn(i, j), randn(j) # Construct from a dictionary of tensors (graph is inferred from shared indices) tn = TensorNetwork(Dictionary(["a", "b", "c"], [t_a, t_b, t_c])) diff --git a/src/TensorNetworks/tensornetwork.jl b/src/TensorNetworks/tensornetwork.jl index ee9f4fe..c05a934 100644 --- a/src/TensorNetworks/tensornetwork.jl +++ b/src/TensorNetworks/tensornetwork.jl @@ -84,7 +84,7 @@ function random_tensornetwork(eltype, g::AbstractGraph; bond_dimension::Integer tensors = Dictionary{vertextype(g), ITensor}() for v in vs is = [l[NamedEdge(v => vn)] for vn in neighbors(g, v)] - set!(tensors, v, random_itensor(eltype, is)) + set!(tensors, v, randn(eltype, is...)) end return TensorNetwork(tensors, g) end diff --git a/src/TensorNetworks/tensornetworkstate.jl b/src/TensorNetworks/tensornetworkstate.jl index b76780d..3b69aa8 100644 --- a/src/TensorNetworks/tensornetworkstate.jl +++ b/src/TensorNetworks/tensornetworkstate.jl @@ -123,7 +123,7 @@ function random_tensornetworkstate(eltype, g::AbstractGraph, siteinds::Dictionar tensors = Dictionary{vertextype(g), ITensor}() for v in vs is = vcat(siteinds[v], [l[NamedEdge(v => vn)] for vn in neighbors(g, v)]) - set!(tensors, v, random_itensor(eltype, is)) + set!(tensors, v, randn(eltype, is...)) end return TensorNetworkState(TensorNetwork(tensors, g), siteinds) end diff --git a/src/itensors.jl b/src/itensors.jl index c3ee4fb..8e2cd06 100644 --- a/src/itensors.jl +++ b/src/itensors.jl @@ -187,12 +187,6 @@ function itensor(array, is...) end itensor(array, is::Union{Tuple, AbstractVector}) = itensor(array, is...) -# Random ITensor over the given indices (legacy `random_itensor`). -random_itensor(eltype::Type, is::Index...) = randn(eltype, is...) -random_itensor(eltype::Type, is::Union{Tuple, AbstractVector}) = randn(eltype, is...) -random_itensor(is::Index...) = random_itensor(Float64, is...) -random_itensor(is::Union{Tuple, AbstractVector}) = random_itensor(Float64, is) - # Rank-0 scalar extraction (legacy `scalar`). scalar(t::AbstractITensor) = t[] diff --git a/test/test_constructors.jl b/test/test_constructors.jl index e04adc8..a48a6da 100644 --- a/test/test_constructors.jl +++ b/test/test_constructors.jl @@ -3,8 +3,8 @@ using Dictionaries: Dictionary using ITensorBase: Index, inds using Random using TensorNetworkQuantumSimulator -# `random_itensor` / `contract` are TNQS-owned; reach them via the `TNQS` alias. `dag` is -# TNQS-owned; `prime` and `inds` are ITensorBase's. +# `contract` is TNQS-owned; reach it via the `TNQS` alias. `dag` is TNQS-owned; `prime`, +# `inds`, and `randn` (over `Index`es) are ITensorBase's. const TNQS = TensorNetworkQuantumSimulator using TensorNetworkQuantumSimulator: dag, prime using Test: @testset, @test, @test_throws @@ -15,7 +15,7 @@ using Test: @testset, @test, @test_throws #TensorNetwork construction from tensors i, j, k, l = Index(2), Index(2), Index(2), Index(2) - A, B, C, D = TNQS.random_itensor(i, j), TNQS.random_itensor(j, k), TNQS.random_itensor(k, l), TNQS.random_itensor(l, i) + A, B, C, D = randn(i, j), randn(j, k), randn(k, l), randn(l, i) t = TensorNetwork([A, B, C, D]) @test t isa TensorNetwork @test scalartype(t) == eltype(A) diff --git a/test/test_contraction_sequences.jl b/test/test_contraction_sequences.jl index 2de5374..2f221da 100644 --- a/test/test_contraction_sequences.jl +++ b/test/test_contraction_sequences.jl @@ -3,7 +3,7 @@ using ITensorBase: Index, name using Random using TensorNetworkQuantumSimulator const TNQS = TensorNetworkQuantumSimulator -# `random_itensor`/`contract`/`scalar` come from TNQS's compat layer (see test_constructors). +# `contract`/`scalar` are TNQS-owned; `randn` (over `Index`es) is ITensorBase's. using TensorNetworkQuantumSimulator: scalar using OMEinsumContractionOrders: NestedEinsum, EinCode, getixsv, getiyv using Test: @testset, @test @@ -20,8 +20,8 @@ collect_leaves!(acc, x) = (for y in x; collect_leaves!(acc, y); end; acc) # 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 = TNQS.random_itensor(i, j) - B = TNQS.random_itensor(j, k) + A = randn(i, j) + B = randn(j, k) code, size_dict = TNQS.to_eincode([A, B]) @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) @@ -57,9 +57,9 @@ collect_leaves!(acc, x) = (for y in x; collect_leaves!(acc, y); end; acc) # --- open network: result is a tensor with dangling indices (iy non-empty). p, q, r, s, t = Index(2), Index(3), Index(2), Index(3), Index(2) - X = TNQS.random_itensor(p, q) - Y = TNQS.random_itensor(q, r, s) - Z = TNQS.random_itensor(s, t) + X = randn(p, q) + Y = randn(q, r, s) + Z = randn(s, t) open_tensors = [X, Y, Z] # open indices: p, r, t seq_open = TNQS.contraction_sequence(open_tensors; alg = "omeinsum", optimizer = GreedyMethod()) @test sort(collect_leaves!(Int[], seq_open)) == [1, 2, 3] From 9dbb8e8088eb080f1cdf1fd2d3918b6234c47c13 Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Wed, 8 Jul 2026 20:24:42 -0400 Subject: [PATCH 03/21] Build the boundary-MPS bond index with trivialrange directly trivialrange on a named range already mints a fresh, backend-preserving range, so wrapping it back in Index(trivialrange(unnamed(...))) is redundant and collapses to trivialrange(...). Co-Authored-By: Claude Opus 4.8 (1M context) --- src/MessagePassing/boundarympscache.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/MessagePassing/boundarympscache.jl b/src/MessagePassing/boundarympscache.jl index 91c6444..2d0e8ef 100644 --- a/src/MessagePassing/boundarympscache.jl +++ b/src/MessagePassing/boundarympscache.jl @@ -195,7 +195,7 @@ function set_interpartition_messages!( # 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 = settags(Index(trivialrange(unnamed(first(inds(m1))), virt_dim)), "m$(i)$(i + 1)") + ind = settags(trivialrange(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. From 76659db38d488b216d9b732dc844ace02c0be018 Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Wed, 8 Jul 2026 21:29:07 -0400 Subject: [PATCH 04/21] Trim the itensors.jl shim toward Base and upstream A mechanical pass replacing legacy shims with the spellings the next-gen stack already provides, shrinking the TNQS-owned shim file: - Remove the dead `itensor` and `algorithm_name` (no callers). - Re-export `TensorAlgebra.scalar` in place of the local `scalar`. - Replace `dim` with `length` at the callsites. - Replace `dag` with `conj`, which covers whole tensors, indices, and (broadcast) index collections. - Drop the `inds(; plev, tags)` filtered wrapper for inline `filter`, taking `inds` from ITensorBase. Co-Authored-By: Claude Opus 4.8 (1M context) --- src/Apply/apply_gates.jl | 2 +- src/Apply/full_update.jl | 20 +++---- src/Apply/simple_update.jl | 4 +- src/Forms/abstractform.jl | 2 +- src/Forms/bilinearform.jl | 2 +- src/Forms/quadraticform.jl | 4 +- src/MessagePassing/beliefpropagationcache.jl | 2 +- src/MessagePassing/boundarympscache.jl | 12 ++-- src/MessagePassing/loopcorrection.jl | 6 +- src/TensorNetworks/abstracttensornetwork.jl | 4 +- src/TensorNetworks/tensornetworkstate.jl | 6 +- src/contraction_sequences.jl | 6 +- src/itensors.jl | 61 ++++---------------- src/sampling.jl | 6 +- src/truncate.jl | 2 +- test/test_constructors.jl | 8 +-- 16 files changed, 53 insertions(+), 94 deletions(-) diff --git a/src/Apply/apply_gates.jl b/src/Apply/apply_gates.jl index 4f5e045..ea0bd84 100644 --- a/src/Apply/apply_gates.jl +++ b/src/Apply/apply_gates.jl @@ -118,7 +118,7 @@ function apply_gate!( # (`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, e, conj(s_values)) setmessage!(ψ_bpc, reverse(e), s_values) end diff --git a/src/Apply/full_update.jl b/src/Apply/full_update.jl index ed0bf55..566783f 100644 --- a/src/Apply/full_update.jl +++ b/src/Apply/full_update.jl @@ -25,7 +25,7 @@ function full_update( ψ[v⃗[2]], uniqueinds(uniqueinds(ψ[v⃗[2]], ψ[v⃗[1]]), uniqueinds(ψ, v⃗[2])) ) - extended_envs = vcat(envs, Qᵥ₁, prime(dag(Qᵥ₁)), Qᵥ₂, prime(dag(Qᵥ₂))) + extended_envs = vcat(envs, Qᵥ₁, prime(conj(Qᵥ₁)), Qᵥ₂, prime(conj(Qᵥ₂))) Rᵥ₁, Rᵥ₂ = optimise_p_q( Rᵥ₁, Rᵥ₂, @@ -69,13 +69,13 @@ function fidelity( p_sind, q_sind = commonind(p_cur, gate), commonind(q_cur, gate) p_sind_sim, q_sind_sim = sim(p_sind), sim(q_sind) gate_sq = - gate * replaceinds(dag(gate), Index[p_sind, q_sind], Index[p_sind_sim, q_sind_sim]) + gate * replaceinds(conj(gate), Index[p_sind, q_sind], Index[p_sind_sim, q_sind_sim]) term1_tns = vcat( [ p_prev, q_prev, - replaceind(prime(dag(p_prev)), prime(p_sind), p_sind_sim), - replaceind(prime(dag(q_prev)), prime(q_sind), q_sind_sim), + replaceind(prime(conj(p_prev)), prime(p_sind), p_sind_sim), + replaceind(prime(conj(q_prev)), prime(q_sind), q_sind_sim), gate_sq, ], envs, @@ -87,14 +87,14 @@ function fidelity( [ p_cur, q_cur, - replaceind(prime(dag(p_cur)), prime(p_sind), p_sind), - replaceind(prime(dag(q_cur)), prime(q_sind), q_sind), + replaceind(prime(conj(p_cur)), prime(p_sind), p_sind), + replaceind(prime(conj(q_cur)), prime(q_sind), q_sind), ], envs, ) sequence = contraction_sequence(term2_tns; alg = "optimal") term2 = contract(term2_tns; sequence) - term3_tns = vcat([p_prev, q_prev, prime(dag(p_cur)), prime(dag(q_cur)), gate], envs) + term3_tns = vcat([p_prev, q_prev, prime(conj(p_cur)), prime(conj(q_cur)), gate], envs) sequence = contraction_sequence(term3_tns; alg = "optimal") term3 = contract(term3_tns; sequence) @@ -103,7 +103,7 @@ function fidelity( end """Do Full Update Sweeping, Optimising the tensors p and q in the presence of the environments envs, -Specifically this functions find the p_cur and q_cur which optimise envs*gate*p*q*dag(prime(p_cur))*dag(prime(q_cur))""" +Specifically this functions find the p_cur and q_cur which optimise envs*gate*p*q*conj(prime(p_cur))*conj(prime(q_cur))""" function optimise_p_q( p::ITensor, q::ITensor, @@ -124,7 +124,7 @@ function optimise_p_q( 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) + ts = vcat(ITensor[p, q, o, conj(prime(r))], envs) sequence = contraction_sequence(ts; alg = "optimal") return noprime(contract(ts; sequence)) end @@ -132,7 +132,7 @@ function optimise_p_q( function M_p(envs::Vector{ITensor}, p_q_tensor::ITensor, s_ind, apply_tensor::ITensor) ts = vcat( ITensor[ - p_q_tensor, replaceinds(prime(dag(p_q_tensor)), prime.(s_ind), s_ind), apply_tensor, + p_q_tensor, replaceinds(prime(conj(p_q_tensor)), prime.(s_ind), s_ind), apply_tensor, ], envs, ) diff --git a/src/Apply/simple_update.jl b/src/Apply/simple_update.jl index ea7dfcc..8518823 100644 --- a/src/Apply/simple_update.jl +++ b/src/Apply/simple_update.jl @@ -71,8 +71,8 @@ function simple_update( 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)]) + Qᵥ₁ = contract([Qᵥ₁; conj.(inv_sqrt_envs_v1)]) + Qᵥ₂ = contract([Qᵥ₂; conj.(inv_sqrt_envs_v2)]) updated_tensors = [Qᵥ₁ * Rᵥ₁, Qᵥ₂ * Rᵥ₂] if normalize_tensors s_values = normalize(s_values) diff --git a/src/Forms/abstractform.jl b/src/Forms/abstractform.jl index 84ae69b..71b0c41 100644 --- a/src/Forms/abstractform.jl +++ b/src/Forms/abstractform.jl @@ -23,7 +23,7 @@ end function default_message(form::AbstractForm, edge::AbstractEdge) cod = virtualinds(ket(form), edge) - dom = dag.(bra_virtualinds(form, edge)) + dom = conj.(bra_virtualinds(form, edge)) return one(similar_map(ket(form)[src(edge)], cod, dom), cod, dom) end diff --git a/src/Forms/bilinearform.jl b/src/Forms/bilinearform.jl index faadbdb..691f3ee 100644 --- a/src/Forms/bilinearform.jl +++ b/src/Forms/bilinearform.jl @@ -19,7 +19,7 @@ function BilinearForm(ket::TensorNetworkState, bra::TensorNetworkState) verts = collect(vertices(ket)) bra = TensorNetworkState(Dictionary(verts, [bra_tensor(bra, v) for v in verts])) operator_tensors = [ - let codomain = dag.(sinds[v]), domain = dag.(prime.(sinds[v])) + let codomain = conj.(sinds[v]), domain = conj.(prime.(sinds[v])) one(similar_map(ket[v], codomain, domain), codomain, domain) end for v in verts ] diff --git a/src/Forms/quadraticform.jl b/src/Forms/quadraticform.jl index e995c7c..72f9342 100644 --- a/src/Forms/quadraticform.jl +++ b/src/Forms/quadraticform.jl @@ -5,9 +5,9 @@ end ket(qf::QuadraticForm) = qf.ket operator(qf::QuadraticForm) = qf.operator -bra(qf::QuadraticForm) = prime(dag(ket(qf))) +bra(qf::QuadraticForm) = prime(conj(ket(qf))) bra_tensor(qf::QuadraticForm, v) = bra_tensor(ket(qf), v) -bra_virtualinds(qf::QuadraticForm, edge::NamedEdge) = dag.(prime.(virtualinds(ket(qf), edge))) +bra_virtualinds(qf::QuadraticForm, edge::NamedEdge) = conj.(prime.(virtualinds(ket(qf), edge))) Base.copy(qf::QuadraticForm) = QuadraticForm(copy(qf.ket), copy(qf.operator)) diff --git a/src/MessagePassing/beliefpropagationcache.jl b/src/MessagePassing/beliefpropagationcache.jl index 2934036..e580c27 100644 --- a/src/MessagePassing/beliefpropagationcache.jl +++ b/src/MessagePassing/beliefpropagationcache.jl @@ -119,7 +119,7 @@ default_bp_update_kwargs(bp_cache::BeliefPropagationCache) = default_bp_update_k function make_hermitian(A::ITensor) A_inds = inds(A) @assert length(A_inds) == 2 - return (A + swapind(dag(A), first(A_inds), last(A_inds))) / 2 + return (A + swapind(conj(A), first(A_inds), last(A_inds))) / 2 end function rescale_messages!(bp_cache::BeliefPropagationCache, edges::Vector{<:AbstractEdge}) diff --git a/src/MessagePassing/boundarympscache.jl b/src/MessagePassing/boundarympscache.jl index 2d0e8ef..36162b0 100644 --- a/src/MessagePassing/boundarympscache.jl +++ b/src/MessagePassing/boundarympscache.jl @@ -133,8 +133,8 @@ function virtual_index_dimension( inds_above = collect(Iterators.flatten(virtualinds.((bmps_cache,), edges_above(bmps_cache, lower_e)))) inds_below = collect(Iterators.flatten(virtualinds.((bmps_cache,), edges_below(bmps_cache, upper_e)))) - x1 = prod(Float64.(dim.(inds_above))) - x2 = prod(Float64.(dim.(inds_below))) + x1 = prod(Float64.(length.(inds_above))) + x2 = prod(Float64.(length.(inds_below))) if network(bmps_cache) isa TensorNetworkState return Int(minimum((x1 * x1, x2 * x2, Float64(mps_bond_dimension(bmps_cache))))) else @@ -200,7 +200,7 @@ function set_interpartition_messages!( # 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 * dag(t)) + setmessage!(bmps_cache, es[i + 1], m2 * conj(t)) end end return bmps_cache @@ -210,8 +210,8 @@ end function switch_message!(bmps_cache::BoundaryMPSCache, e::NamedEdge) ms = messages(bmps_cache) me, mer = message(bmps_cache, e), message(bmps_cache, reverse(e)) - set!(ms, e, dag(mer)) - set!(ms, reverse(e), dag(me)) + set!(ms, e, conj(mer)) + set!(ms, reverse(e), conj(me)) return bmps_cache end @@ -308,7 +308,7 @@ function inserter!( update_e::NamedEdge, m::ITensor ) - setmessage!(bmps_cache, reverse(update_e), dag(m)) + setmessage!(bmps_cache, reverse(update_e), conj(m)) return bmps_cache end diff --git a/src/MessagePassing/loopcorrection.jl b/src/MessagePassing/loopcorrection.jl index dd30010..2ea2013 100644 --- a/src/MessagePassing/loopcorrection.jl +++ b/src/MessagePassing/loopcorrection.jl @@ -38,7 +38,7 @@ function sim_edgeinduced_subgraph(bpc::BeliefPropagationCache, eg) linds_sim = sim.(linds) mer = replaceinds(mer, linds, linds_sim) if network(bpc) isa TensorNetworkState - mer = replaceinds(mer, dag.(prime.(linds)), dag.(prime.(linds_sim))) + mer = replaceinds(mer, conj.(prime.(linds)), conj.(prime.(linds_sim))) end ms = messages(bpc) set!(ms, reverse(e), mer) @@ -60,14 +60,14 @@ function sim_edgeinduced_subgraph(bpc::BeliefPropagationCache, eg) # 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. + # relative to the passed indices, so the columns go in `conj`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 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 - ap = adapt_like(me, identity_tensor(row_inds, dag.(col_inds))) + ap = adapt_like(me, identity_tensor(row_inds, conj.(col_inds))) ap = ap - me * mer push!(antiprojectors, ap) end diff --git a/src/TensorNetworks/abstracttensornetwork.jl b/src/TensorNetworks/abstracttensornetwork.jl index 6df43cf..99c0600 100644 --- a/src/TensorNetworks/abstracttensornetwork.jl +++ b/src/TensorNetworks/abstracttensornetwork.jl @@ -24,7 +24,7 @@ virtualinds(tn::AbstractTensorNetwork, e::NamedEdge) = commoninds(tn[src(e)], tn virtualind(tn::AbstractTensorNetwork, e::NamedEdge) = only(virtualinds(tn, e)) function maxvirtualdim(tn::AbstractTensorNetwork) - return maximum(maximum.([dim.(virtualinds(tn, e)) for e in edges(tn)])) + return maximum(maximum.([length.(virtualinds(tn, e)) for e in edges(tn)])) end # Compare by name, not by `Index` equality: a shared graded link is stored nondual @@ -123,7 +123,7 @@ function add(tn1::AbstractTensorNetwork, tn2::AbstractTensorNetwork) es, [ Index( - dim(only(virtualinds(tn1, e))) + dim(only(virtualinds(tn2, e))), + length(only(virtualinds(tn1, e))) + length(only(virtualinds(tn2, e))), ) for e in es ], ), diff --git a/src/TensorNetworks/tensornetworkstate.jl b/src/TensorNetworks/tensornetworkstate.jl index 3b69aa8..825ec4d 100644 --- a/src/TensorNetworks/tensornetworkstate.jl +++ b/src/TensorNetworks/tensornetworkstate.jl @@ -54,11 +54,11 @@ function Base.setindex!(tns::TensorNetworkState, value::ITensor, v) return tns end -# Bra copy of a tensor: `dag` and prime all legs except `auxinds` — dangling +# Bra copy of a tensor: `conj` 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)) + tdag = conj(prime(t)) return isempty(auxinds) ? tdag : replaceinds(tdag, prime.(auxinds), auxinds) end bra_tensor(tns::TensorNetworkState, v) = bra_tensor(tns[v], auxinds(tns, v)) @@ -96,7 +96,7 @@ bp_factors(tns::TensorNetworkState, v) = norm_factors(tns, v) # links give a graded message; the legacy `delta` filled a dense diagonal). function default_message(tns::TensorNetworkState, edge::AbstractEdge) linds = virtualinds(tns, edge) - cod, dom = Tuple(linds), Tuple(prime.(dag(linds))) + cod, dom = Tuple(linds), Tuple(prime.(conj.(linds))) return adapt_like(tns, one(zeros(scalartype(tns), cod..., dom...), cod, dom)) end diff --git a/src/contraction_sequences.jl b/src/contraction_sequences.jl index 73b572d..3f868db 100644 --- a/src/contraction_sequences.jl +++ b/src/contraction_sequences.jl @@ -7,7 +7,7 @@ using OMEinsumContractionOrders: OMEinsumContractionOrders, optimize_code, EinCo # so the returned sequence still indexes the original tensor list. Only the index sets feed # `optimaltree`; the tensors themselves are never contracted here (the caller contracts the # untouched originals), so no placeholder tensor is needed. -is_trivial_tensor(t::ITensor) = all(i -> dim(i) == 1, inds(t)) +is_trivial_tensor(t::ITensor) = all(i -> length(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 @@ -23,7 +23,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(name(i) => Float64(dim(i)) for t in tensors for i in inds(t)) + inds_to_dims = Dict(name(i) => Float64(length(i)) for t in tensors for i in inds(t)) seq, _ = optimaltree(network, inds_to_dims) seq = typeof(seq) <: Int ? [seq] : seq return seq @@ -44,7 +44,7 @@ function to_eincode(tensors::Vector{<:ITensor}) ixs = map(t -> collect(name.(inds(t))), tensors) LT = eltype(eltype(ixs)) 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)) + size_dict = Dict{LT, Int}(name(i) => length(i) for t in tensors for i in inds(t)) return EinCode(ixs, iy), size_dict end diff --git a/src/itensors.jl b/src/itensors.jl index 8e2cd06..ca61fe1 100644 --- a/src/itensors.jl +++ b/src/itensors.jl @@ -9,9 +9,9 @@ import MatrixAlgebraKit as MAK import TensorAlgebra: matricize using Adapt: Adapt using ITensorBase: ITensorBase, AbstractITensor, ITensor, Index, NamedUnitRange, dimnames, - id, name, nameddims, noprime, plev, prime, replacedimnames, sim, tags, unnamed + id, inds, name, nameddims, noprime, plev, prime, replacedimnames, sim, tags, unnamed using LinearAlgebra: LinearAlgebra -using TensorAlgebra: TensorAlgebra, project, tryproject +using TensorAlgebra: TensorAlgebra, project, scalar, tryproject # # State-vector projection, shared by the `Ops` state constructors and by `onehot`. @@ -46,15 +46,6 @@ function onehot(eltype::Type, (i, p)::Pair{<:Index}) end onehot(p::Pair{<:Index}) = onehot(Float64, p) -# Legacy `inds(t; plev, tags)` took index-filtering keywords; `ITensorBase.inds` takes none. -# This forwards to it and applies the legacy filters. -function inds(t::AbstractITensor; plev = nothing, tags = nothing) - is = ITensorBase.inds(t) - isnothing(plev) || (is = filter(i -> ITensorBase.plev(i) == plev, is)) - isnothing(tags) || (is = filter(i -> hastags(i, tags), is)) - return is -end - # `contract` / `inner` are TNQS operations whose base generics live here and are extended # for TNQS types (tensor networks, caches) in the respective files. function inner end @@ -126,17 +117,6 @@ noncommonind(a, b) = (us = uniqueinds(a, b); isempty(us) ? nothing : first(us)) # Plural: indices not shared by both (symmetric difference). 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) -dim(is::Union{Tuple, AbstractVector}) = isempty(is) ? 1 : prod(length, is) - -# Conjugate (legacy `dag`): `conj` the tensor, and on bare indices flip the sector arrows -# on a graded axis. `conj(::Index)` is id-preserving on the dense backend, so `dag` there -# is effectively the identity on indices, matching legacy behavior. -dag(t::AbstractITensor) = conj(t) -dag(i::Index) = conj(i) -dag(is::Union{Tuple, AbstractVector}) = map(conj, is) - # `replaceind` (singular) maps to a single-pair replacement, forwarding to `replaceinds`. replaceind(t, p::Pair) = replaceinds(t, p) replaceind(t, from::Index, to::Index) = replaceinds(t, from => to) @@ -170,26 +150,6 @@ function replaceinds(t::AbstractITensor, p::Pair{<:_IndexColl, <:_IndexColl}) return replaceinds(t, first(p), last(p)) 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. 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...) - -# Rank-0 scalar extraction (legacy `scalar`). -scalar(t::AbstractITensor) = t[] - # Dense Kronecker delta tensor (legacy `delta`), vendored from ITensorNetworksNext's # `ITensorNetworkGenerators/delta_network.jl`. A graded/sector-aware `delta` is a stack gap # (tracked separately). @@ -247,13 +207,13 @@ delta(is::Index...) = delta(Float64, is) delta(is::AbstractVector{<:Index}) = delta(Float64, Tuple(is)) # Trace over prime pairs (legacy ITensors `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. Named `itensor_tr` (not a `LinearAlgebra.tr` method, which would be piracy on -# `AbstractITensor`) so it stays distinct from `tr` on plain matrices, which TNQS also calls. +# unprimed index with its prime. The domain is built as `conj.(prime.(codomain))`, not +# by filtering `plev == 1` — `plev` filtering does not preserve the pairing order between the +# plev-0 and plev-1 groups. Named `itensor_tr` (not a `LinearAlgebra.tr` method, which would be +# piracy on `AbstractITensor`) so it stays distinct from `tr` on plain matrices, which TNQS also calls. function itensor_tr(t::AbstractITensor) - unprimed = inds(t; plev = 0) - codomain, domain = dag.(unprimed), dag.(prime.(unprimed)) + unprimed = filter(i -> plev(i) == 0, inds(t)) + codomain, domain = conj.(unprimed), conj.(prime.(unprimed)) return scalar(t * one(similar_map(t, codomain, domain), codomain, domain)) end @@ -320,7 +280,7 @@ function _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`) — +# `m = Vt * D * conj(V)` (with `Vt` the relabeling of `V` from `Rinds` to `Linds`) — # no conjugation of `U`. function eigen(m::AbstractITensor, codomain, domain; kwargs...) return _eigh(m, codomain, domain; kwargs...) @@ -329,7 +289,7 @@ end # No-partition form `eigen(m)` matches legacy ITensors' `eigen(A)`, which auto-partitions # `Ris = filterinds(plev = 0)`, `Lis = Ris'`. That orientation is the adjoint view of the # partitioned call, so `eigh_full` returns the conjugate eigenvectors; `conj(U)` recovers -# the convention `symmetric_gauge` uses (`U * D * prime(dag(U)) == m`, a true matrix sqrt). +# the convention `symmetric_gauge` uses (`U * D * prime(conj(U)) == m`, a true matrix sqrt). function eigen(m::AbstractITensor; kwargs...) is = collect(inds(m)) p0 = filter(i -> plev(i) == 0, is) @@ -498,7 +458,6 @@ function Base.getproperty(alg::Algorithm, name::Symbol) getfield(getfield(alg, :kwargs), name) end end -algorithm_name(::Algorithm{Alg}) where {Alg} = Alg macro Algorithm_str(s) return :(Algorithm{$(Expr(:quote, Symbol(s)))}) end diff --git a/src/sampling.jl b/src/sampling.jl index 9116f0d..caec56b 100644 --- a/src/sampling.jl +++ b/src/sampling.jl @@ -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(ρ, dag(onehot(s_ind => config))) + P = adapt_like(ρ, conj(onehot(s_ind => config))) setindex_preserve!(projected_bp_cache, ψv * P, v) if v != last(vertices(ψ)) @@ -246,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(ρ, dag(onehot(s_ind => config))) + P = adapt_like(ρ, conj(onehot(s_ind => config))) q = ρ_diag[config] logq += log(q) Pψv = copy(network(norm_bmps_cache)[v]) * inv(sqrt(q)) * P @@ -274,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], dag(onehot(only(s[v]) => bitstring[v] + 1))) + P = adapt_like(ψproj[v], conj(onehot(only(s[v]) => bitstring[v] + 1))) setindex_preserve!(ψproj, ψproj[v] * P * inv(qv), v) end diff --git a/src/truncate.jl b/src/truncate.jl index c105fa1..830fa53 100644 --- a/src/truncate.jl +++ b/src/truncate.jl @@ -5,7 +5,7 @@ default_truncate_alg(tns::TensorNetworkState) = nothing function truncatable_edge(cache::AbstractBeliefPropagationCache, e::NamedEdge) vinds = virtualinds(cache, e) isempty(vinds) && return false - all([dim(vind) == 1 for vind in vinds]) && return false + all([length(vind) == 1 for vind in vinds]) && return false return true end diff --git a/test/test_constructors.jl b/test/test_constructors.jl index a48a6da..c17da93 100644 --- a/test/test_constructors.jl +++ b/test/test_constructors.jl @@ -3,10 +3,10 @@ using Dictionaries: Dictionary using ITensorBase: Index, inds using Random using TensorNetworkQuantumSimulator -# `contract` is TNQS-owned; reach it via the `TNQS` alias. `dag` is TNQS-owned; `prime`, -# `inds`, and `randn` (over `Index`es) are ITensorBase's. +# `contract` is TNQS-owned; reach it via the `TNQS` alias. `prime`, `inds`, and `randn` +# (over `Index`es) are ITensorBase's; `conj` is `Base`. const TNQS = TensorNetworkQuantumSimulator -using TensorNetworkQuantumSimulator: dag, prime +using TensorNetworkQuantumSimulator: prime using Test: @testset, @test, @test_throws @@ -34,7 +34,7 @@ using Test: @testset, @test, @test_throws @test maxvirtualdim(ψ) == 3 @test all([length(inds(ψ[v])) == degree(g, v) for v in vertices(ψ)]) - ψdag = map_virtualinds(prime, map_tensors(dag, ψ)) + ψdag = map_virtualinds(prime, map_tensors(conj, ψ)) @test ψdag isa TensorNetwork @test TNQS.contract(ψdag; alg = "exact") ≈ conj(TNQS.contract(ψ; alg = "exact")) From ac92b6dfd9b87679af44f8b479577ad6f11d9895 Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Wed, 8 Jul 2026 21:41:25 -0400 Subject: [PATCH 05/21] Strip the explanatory comments from the itensors.jl shim The shim is transitional and shrinking, so the per-symbol rationale is recorded outside the source rather than carried inline. Code is unchanged, and the full test suite still passes. Co-Authored-By: Claude Opus 4.8 (1M context) --- src/itensors.jl | 169 ++---------------------------------------------- 1 file changed, 4 insertions(+), 165 deletions(-) diff --git a/src/itensors.jl b/src/itensors.jl index ca61fe1..9795360 100644 --- a/src/itensors.jl +++ b/src/itensors.jl @@ -1,7 +1,5 @@ -# TNQS-owned tensor utilities and the legacy `ITensors.jl`-style API TNQS was written -# against, implemented over the next-gen `ITensorBase` / `TensorAlgebra` / `MatrixAlgebraKit` -# stack. ITensorBase keeps most of this API internal (unexported), so we reach for the -# qualified names and define the legacy spellings here as functionality TNQS owns. +# TNQS-owned tensor utilities and the legacy `ITensors.jl`-style API, implemented over the +# next-gen `ITensorBase` / `TensorAlgebra` / `MatrixAlgebraKit` stack. import Base: truncate import ITensorBase: scalartype @@ -13,12 +11,6 @@ using ITensorBase: ITensorBase, AbstractITensor, ITensor, Index, NamedUnitRange, using LinearAlgebra: LinearAlgebra using TensorAlgebra: TensorAlgebra, project, scalar, tryproject -# -# State-vector projection, shared by the `Ops` state constructors and by `onehot`. -# -# Project a raw vector as an `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`). function project_aux(v::AbstractVector{<:Number}, i::Index) length(v) == length(i) || error( @@ -26,19 +18,11 @@ function project_aux(v::AbstractVector{<:Number}, i::Index) ) ψ = 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 -# One-hot vector along `i` at position `p` (legacy `onehot(i => p)`), through `project_aux`, -# 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) @@ -46,14 +30,8 @@ function onehot(eltype::Type, (i, p)::Pair{<:Index}) end onehot(p::Pair{<:Index}) = onehot(Float64, p) -# `contract` / `inner` are TNQS operations whose base generics live here and are extended -# for TNQS types (tensor networks, caches) in the respective files. function inner end -# Base contraction of a list of ITensors along a (possibly nested) pairwise `sequence` (legacy -# `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) @@ -61,19 +39,9 @@ end _contract_sequence(tensors, s::Integer) = tensors[s] _contract_sequence(tensors, s) = reduce(*, (_contract_sequence(tensors, x) for x in s)) -# Get the index collection of an `ITensor`, or pass an index collection through -# unchanged. TNQS nests these calls (e.g. `uniqueinds(uniqueinds(...), ...)`), so the -# index-set helpers below accept both tensors and bare index collections. _compat_inds(t::AbstractITensor) = inds(t) _compat_inds(is) = is -# -# 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. -# 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)) @@ -82,13 +50,6 @@ smallisdisjoint(a, b; by = identity) = (kb = Iterators.map(by, b); !any(x -> by( 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) @@ -98,47 +59,22 @@ 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 -# eigenvalue tensor); revisit if a symmetric-difference callsite turns up. 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) = namesymdiff(_compat_inds(a), _compat_inds(b)) -# `replaceind` (singular) maps to a single-pair replacement, forwarding to `replaceinds`. replaceind(t, p::Pair) = replaceinds(t, p) replaceind(t, from::Index, to::Index) = replaceinds(t, from => to) -# Concatenate indices and/or index collections into a single `Vector{Index}`. Legacy -# ITensors code spells this `vcat(i, j)` / `vcat(is, i)`, but next-gen `Index` is a -# `NamedUnitRange` (`<: AbstractVector`), so `vcat` of bare indices concatenates their -# integer ranges instead of collecting the indices. Used where the legacy idiom builds -# an index list from a mix of single indices and collections. _as_index_vec(x::Index) = [x] _as_index_vec(xs) = collect(xs) cat_inds(xs...) = reduce(vcat, map(_as_index_vec, xs)) -# Legacy `replaceinds` took collection arguments (`replaceinds(t, [i,k], [j,l])`, `... => ...`); -# ITensorBase provides only the pair-splat form, so this handles the collection forms. -# -# 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` (`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 replacedimnames(t, map(p -> name(first(p)) => name(last(p)), pairs)...) @@ -150,9 +86,6 @@ function replaceinds(t::AbstractITensor, p::Pair{<:_IndexColl, <:_IndexColl}) return replaceinds(t, first(p), last(p)) end -# Dense Kronecker delta tensor (legacy `delta`), vendored from ITensorNetworksNext's -# `ITensorNetworkGenerators/delta_network.jl`. A graded/sector-aware `delta` is a stack gap -# (tracked separately). diaglength(a::AbstractArray) = minimum(size(a)) function diagstride(a::AbstractArray) s = 1 @@ -182,10 +115,6 @@ 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) @@ -196,9 +125,6 @@ function similar_map(prototype::AbstractITensor, codomain, domain) return similar_map(prototype, scalartype(prototype), codomain, 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)) @@ -206,40 +132,22 @@ delta(is::Tuple) = delta(Float64, is) delta(is::Index...) = delta(Float64, is) delta(is::AbstractVector{<:Index}) = delta(Float64, Tuple(is)) -# Trace over prime pairs (legacy ITensors `tr`): contract with the identity map pairing each -# unprimed index with its prime. The domain is built as `conj.(prime.(codomain))`, not -# by filtering `plev == 1` — `plev` filtering does not preserve the pairing order between the -# plev-0 and plev-1 groups. Named `itensor_tr` (not a `LinearAlgebra.tr` method, which would be -# piracy on `AbstractITensor`) so it stays distinct from `tr` on plain matrices, which TNQS also calls. function itensor_tr(t::AbstractITensor) unprimed = filter(i -> plev(i) == 0, inds(t)) codomain, domain = conj.(unprimed), conj.(prime.(unprimed)) return scalar(t * one(similar_map(t, codomain, domain), codomain, domain)) end -# -# Factorizations. These map onto MatrixAlgebraKit's named-tensor methods -# (`f(a, codomain, domain)`). The legacy return shapes differ from MAK's, so the -# heavy factorization callsites (simple_update / full_update / symmetric_gauge) are -# translated directly rather than fully wrapped here; these aliases cover the simple uses. -# const svd_trunc = MAK.svd_trunc -# Legacy `qr(a, linds)`: `Q` over `(linds..., q)` (isometric), `R` over `(q, rest...)`. -# `linds` may be splatted, a single index, or one collection. function qr(a::AbstractITensor, linds...) left = cat_inds(linds...) right = namesetdiff(inds(a), left) return MAK.qr_compact(a, Tuple(left), Tuple(right)) end -# Legacy `apply(o, ψ)` (gate application): contract `o`'s unprimed legs with `ψ`'s -# matching indices, then drop the prime so `o`'s output legs become `ψ`'s site legs. -# Covers the one- and two-site gates TNQS applies to states without pre-existing primes. apply(o::AbstractITensor, ψ::AbstractITensor) = noprime(o * ψ) -# `svd` / `eigen` accept a single `Index` or a collection as the codomain, and -# `eigen` auto-partitions a 2-index tensor when no split is given (legacy behavior). _astuple(x::Index) = (x,) _astuple(x) = Tuple(x) @@ -250,15 +158,6 @@ function svd(a::AbstractITensor, codomain, domain; kwargs...) return MAK.svd_compact(a, _astuple(codomain), _astuple(domain); kwargs...) end -# Legacy `eigen` here is a hermitian eigendecomposition returning `(D, U)`. It routes -# to `MAK.eigh_full` and then renames indices into the legacy ITensors convention the -# TNQS callsites assume: -# - `U` over `(domain..., u)` where `u` is the eigenvalue index shared with `D`. -# - `D` diagonal over `(prime(u), u)` — a prime pair, like legacy `eigen`. -# `eigh_full` returns `D` over two independent fresh indices; we rename its `U`-disjoint -# index to `prime(u)` so `D` becomes the legacy prime-paired diagonal. The legacy -# truncation kwargs (`cutoff`/`maxdim`/`mindim`) are not yet translated to MAK's -# `trunc=(; ...)` spec; the current callsites pass `cutoff = nothing` (full decomposition). function _eigh( m::AbstractITensor, codomain, @@ -273,23 +172,16 @@ function _eigh( ) cod, dom = _astuple(codomain), _astuple(domain) 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 + u = only(commoninds(D, U)) + t = only(uniqueinds(D, U)) D = replaceinds(D, t => ITensorBase.prime(u)) return D, U end -# Partitioned form `eigen(m, Linds, Rinds)` reproduces legacy ITensors' reconstruction -# `m = Vt * D * conj(V)` (with `Vt` the relabeling of `V` from `Rinds` to `Linds`) — -# no conjugation of `U`. function eigen(m::AbstractITensor, codomain, domain; kwargs...) return _eigh(m, codomain, domain; kwargs...) end -# No-partition form `eigen(m)` matches legacy ITensors' `eigen(A)`, which auto-partitions -# `Ris = filterinds(plev = 0)`, `Lis = Ris'`. That orientation is the adjoint view of the -# partitioned call, so `eigh_full` returns the conjugate eigenvectors; `conj(U)` recovers -# the convention `symmetric_gauge` uses (`U * D * prime(conj(U)) == m`, a true matrix sqrt). function eigen(m::AbstractITensor; kwargs...) is = collect(inds(m)) p0 = filter(i -> plev(i) == 0, is) @@ -301,11 +193,6 @@ function eigen(m::AbstractITensor; kwargs...) return D, conj(U) end -# 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)) @@ -313,11 +200,6 @@ function itensor_trunc(; maxdim = nothing, cutoff = nothing) 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)) @@ -326,12 +208,6 @@ function _bipartition_inds(a::AbstractITensor, linds) return left, right end -# Legacy `factorize(a, linds...; ortho, cutoff, maxdim, tags)` splits `a` into `L * R` -# with the new bond between them. `linds` selects `L`'s indices (besides the bond) and -# may be passed splatted, as a single index, or as one collection. `ortho = "left"` -# makes `L` isometric, `"right"` makes `R` isometric. With no truncation this is a -# plain QR/LQ; with `cutoff`/`maxdim` it is a truncated SVD whose singular values are -# absorbed into the non-isometric factor. `tags`, if given, names the new bond. function factorize( a::AbstractITensor, linds...; @@ -372,48 +248,25 @@ function factorize( return L, R end -# -# Storage / element type accessors. `scalartype` is imported from ITensorBase (extended for -# TNQS types elsewhere). `datatype` is the underlying storage array type (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) = convert(Array, unnamed(T)) data(T::AbstractITensor) = unnamed(T) -# Scalar (element) type conversion as an owned `Adapt` adapter. `ScalarTypeAdaptor` is ours, so -# `adapt(ScalarTypeAdaptor(T), x)` is not piracy, and `adapt_scalartype` mirrors `adapt`'s curried -# and applied forms. Reproduces legacy `adapt(eltype)(t)` scalar conversion (ITensorBase's Adapt -# integration adapts storage/device but leaves the element type alone). struct ScalarTypeAdaptor{T} end ScalarTypeAdaptor(T::Type) = ScalarTypeAdaptor{T}() adapt_scalartype(T::Type) = Adapt.adapt(ScalarTypeAdaptor(T)) adapt_scalartype(T::Type, x) = Adapt.adapt(ScalarTypeAdaptor(T), x) function Adapt.adapt_structure(::ScalarTypeAdaptor{elt}, T::AbstractITensor) where {elt} - # 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 -# `swapind`: swap two indices (legacy convenience over `replaceinds`). swapind(T::AbstractITensor, i::Index, j::Index) = replaceinds(T, i => j, j => i) -# Whether a tensor carries quantum-number (graded) block structure: true when any of its -# indices is graded. `loopcorrection` branches on this to pick a contraction-order -# algorithm. A graded axis differs from its conjugate (conjugation flips the sector -# arrows) while a dense axis is self-conjugate — a dependency-free discriminator that needs -# no GradedArrays import and hardcodes no dense range type. hasqns(i::Index) = conj(unnamed(i)) != unnamed(i) hasqns(t::AbstractITensor) = any(hasqns, inds(t)) hasqns(::Any) = false -# 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) t1, s1 = first(pairs[1]), Tuple(last(pairs[1])) @@ -437,12 +290,6 @@ function directsum(out_inds, pairs::Pair...) return out[target...] end -# -# Algorithm dispatch tag (legacy `Algorithm` / `@Algorithm_str`). `Algorithm"exact"` -# is the type `Algorithm{:exact}` (usable in `::Algorithm"exact"` signatures); -# `Algorithm("exact"; kwargs...)` constructs an instance carrying keyword options. -# Faithful minimal copy of the ITensors/NDTensors helper. -# struct Algorithm{Alg, Kwargs <: NamedTuple} kwargs::Kwargs end @@ -462,12 +309,6 @@ macro Algorithm_str(s) return :(Algorithm{$(Expr(:quote, Symbol(s)))}) end -# -# Tags. Legacy ITensors uses a flat tag set (`Index(dim, "S=1/2,Site")`); ITensorBase -# stores tags as a `Dict{String, String}`. For the legacy flat-tag usage TNQS has -# (site-type labels, link names), we store each token as a keyed tag with empty -# value, and `hastags` checks membership. (A fuller tag-compat story is a follow-up.) -# function settags(i::Index, tagstr::AbstractString) for t in split(tagstr, ",") s = String(strip(t)) @@ -475,8 +316,6 @@ function settags(i::Index, tagstr::AbstractString) end return i end -# Copy a whole tag dictionary onto an index (legacy `settags(i, tags(j))`, used when a -# factorization is asked to give its new bond the same tags as an existing bond). function settags(i::Index, d::AbstractDict) for (k, v) in d i = ITensorBase.settag(i, k, v) From b96afd5f139f584ed24813493d610891a98ae8ba Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Thu, 9 Jul 2026 12:00:18 -0400 Subject: [PATCH 06/21] Remove the local index-set helpers in favor of ITensorBase and Base set-ops The shim carried its own small-collection set operations and the commoninds/uniqueinds wrappers built on them. ITensorBase now provides these as public functions, keyed by index name internally so a graded bond still matches its dual, and NamedUnitRange equality is dual-insensitive, so Base set-ops work directly on index collections. This removes the duplicated helpers: tensor-argument callsites use ITensorBase's commoninds/uniqueinds/unioninds/hascommoninds/noncommoninds, and collection or mixed-argument callsites use Base's setdiff/intersect/union/symdiff/isdisjoint. The thin commonind/noncommonind conveniences stay, returning a single index or nothing. Co-Authored-By: Claude Opus 4.8 (1M context) --- src/Apply/full_update.jl | 8 ++--- src/Apply/simple_update.jl | 6 ++-- src/MessagePassing/loopcorrection.jl | 14 ++++----- src/TensorNetworks/tensornetworkstate.jl | 4 +-- src/contraction_sequences.jl | 2 +- src/itensors.jl | 38 +++++------------------- src/sampling.jl | 2 +- 7 files changed, 26 insertions(+), 48 deletions(-) diff --git a/src/Apply/full_update.jl b/src/Apply/full_update.jl index 566783f..0d90e41 100644 --- a/src/Apply/full_update.jl +++ b/src/Apply/full_update.jl @@ -19,10 +19,10 @@ function full_update( ) Qᵥ₁, Rᵥ₁ = factorize( - ψ[v⃗[1]], uniqueinds(uniqueinds(ψ[v⃗[1]], ψ[v⃗[2]]), uniqueinds(ψ, v⃗[1])) + ψ[v⃗[1]], setdiff(uniqueinds(ψ[v⃗[1]], ψ[v⃗[2]]), uniqueinds(ψ, v⃗[1])) ) Qᵥ₂, Rᵥ₂ = factorize( - ψ[v⃗[2]], uniqueinds(uniqueinds(ψ[v⃗[2]], ψ[v⃗[1]]), uniqueinds(ψ, v⃗[2])) + ψ[v⃗[2]], setdiff(uniqueinds(ψ[v⃗[2]], ψ[v⃗[1]]), uniqueinds(ψ, v⃗[2])) ) extended_envs = vcat(envs, Qᵥ₁, prime(conj(Qᵥ₁)), Qᵥ₂, prime(conj(Qᵥ₂))) @@ -120,8 +120,8 @@ function optimise_p_q( fstart = print_fidelity_loss ? fidelity(envs, p_cur, q_cur, p, q, o) : 0 - 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))))) + 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))))) function b(p::ITensor, q::ITensor, o::ITensor, envs::Vector{ITensor}, r::ITensor) ts = vcat(ITensor[p, q, o, conj(prime(r))], envs) diff --git a/src/Apply/simple_update.jl b/src/Apply/simple_update.jl index 8518823..cbc2c79 100644 --- a/src/Apply/simple_update.jl +++ b/src/Apply/simple_update.jl @@ -50,8 +50,8 @@ function simple_update( ψᵥ₂ = contract([ψ⃗[2]; sqrt_envs_v2]) sᵥ₁ = commoninds(ψ⃗[1], o) sᵥ₂ = commoninds(ψ⃗[2], o) - Qᵥ₁, Rᵥ₁ = qr(ψᵥ₁, uniqueinds(uniqueinds(ψᵥ₁, ψᵥ₂), sᵥ₁)) - Qᵥ₂, Rᵥ₂ = qr(ψᵥ₂, uniqueinds(uniqueinds(ψᵥ₂, ψᵥ₁), sᵥ₂)) + Qᵥ₁, Rᵥ₁ = qr(ψᵥ₁, setdiff(uniqueinds(ψᵥ₁, ψᵥ₂), sᵥ₁)) + Qᵥ₂, Rᵥ₂ = qr(ψᵥ₂, setdiff(uniqueinds(ψᵥ₂, ψᵥ₁), sᵥ₂)) rᵥ₁ = commoninds(Qᵥ₁, Rᵥ₁) rᵥ₂ = commoninds(Qᵥ₂, Rᵥ₂) oR = apply(o, Rᵥ₁ * Rᵥ₂) @@ -59,7 +59,7 @@ function simple_update( # 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, S, V = svd_trunc(oR, union(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) diff --git a/src/MessagePassing/loopcorrection.jl b/src/MessagePassing/loopcorrection.jl index 2ea2013..8a22cfc 100644 --- a/src/MessagePassing/loopcorrection.jl +++ b/src/MessagePassing/loopcorrection.jl @@ -43,9 +43,9 @@ function sim_edgeinduced_subgraph(bpc::BeliefPropagationCache, eg) ms = messages(bpc) set!(ms, reverse(e), mer) t = network(bpc)[src(e)] - # 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) + # On a graded backend the network tensor carries the dual of the message's axis. + # Index equality is dual-insensitive, so `intersect` matches the two. + t_inds = intersect(inds(t), linds) if !isempty(t_inds) t_ind = only(t_inds) t_ind_pos = findfirst(x -> name(x) == name(t_ind), linds) @@ -61,11 +61,11 @@ function sim_edgeinduced_subgraph(bpc::BeliefPropagationCache, eg) # 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 `conj`ed. - row_inds = Index[only(commoninds(me, [l])) for l in linds] - col_inds = Index[only(commoninds(mer, [l])) for l in linds_sim] + row_inds = Index[only(intersect(inds(me), [l])) for l in linds] + col_inds = Index[only(intersect(inds(mer), [l])) for l in linds_sim] if network(bpc) isa TensorNetworkState - 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]) + append!(row_inds, Index[only(intersect(inds(me), [prime(l)])) for l in linds]) + append!(col_inds, Index[only(intersect(inds(mer), [prime(l)])) for l in linds_sim]) end ap = adapt_like(me, identity_tensor(row_inds, conj.(col_inds))) ap = ap - me * mer diff --git a/src/TensorNetworks/tensornetworkstate.jl b/src/TensorNetworks/tensornetworkstate.jl index 825ec4d..dbed9b0 100644 --- a/src/TensorNetworks/tensornetworkstate.jl +++ b/src/TensorNetworks/tensornetworkstate.jl @@ -64,7 +64,7 @@ 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(tns::TensorNetworkState, v) = Index[i for i in setdiff(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, @@ -219,5 +219,5 @@ function tensornetworkstate(f::Function, args...) end function NamedGraphs.vertices(t::ITensor, tns::TensorNetworkState) - return filter(v -> !nameisdisjoint(inds(t), siteinds(tns, v)), collect(vertices(tns))) + return filter(v -> !isdisjoint(inds(t), siteinds(tns, v)), collect(vertices(tns))) end diff --git a/src/contraction_sequences.jl b/src/contraction_sequences.jl index 3f868db..ee32d94 100644 --- a/src/contraction_sequences.jl +++ b/src/contraction_sequences.jl @@ -43,7 +43,7 @@ end function to_eincode(tensors::Vector{<:ITensor}) ixs = map(t -> collect(name.(inds(t))), tensors) LT = eltype(eltype(ixs)) - iy = collect(LT, name.(reduce(noncommoninds, tensors))) + iy = collect(LT, name.(reduce(symdiff, inds.(tensors)))) size_dict = Dict{LT, Int}(name(i) => length(i) for t in tensors for i in inds(t)) return EinCode(ixs, iy), size_dict end diff --git a/src/itensors.jl b/src/itensors.jl index 9795360..74dfbd4 100644 --- a/src/itensors.jl +++ b/src/itensors.jl @@ -2,12 +2,13 @@ # next-gen `ITensorBase` / `TensorAlgebra` / `MatrixAlgebraKit` stack. import Base: truncate -import ITensorBase: scalartype +import ITensorBase: scalartype, uniqueinds import MatrixAlgebraKit as MAK import TensorAlgebra: matricize using Adapt: Adapt -using ITensorBase: ITensorBase, AbstractITensor, ITensor, Index, NamedUnitRange, dimnames, - id, inds, name, nameddims, noprime, plev, prime, replacedimnames, sim, tags, unnamed +using ITensorBase: ITensorBase, AbstractITensor, ITensor, Index, NamedUnitRange, commoninds, + dimnames, hascommoninds, id, inds, name, nameddims, noncommoninds, noprime, plev, prime, + replacedimnames, sim, tags, unioninds, unnamed using LinearAlgebra: LinearAlgebra using TensorAlgebra: TensorAlgebra, project, scalar, tryproject @@ -39,34 +40,11 @@ end _contract_sequence(tensors, s::Integer) = tensors[s] _contract_sequence(tensors, s) = reduce(*, (_contract_sequence(tensors, x) for x in s)) -_compat_inds(t::AbstractITensor) = inds(t) -_compat_inds(is) = is - -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) - -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) - -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)) - +# The single shared/unique index, or `nothing` when there is not exactly one. `commoninds` +# and `uniqueinds` come from ITensorBase, which keys index-set algebra by name so a graded +# bond still matches its dual. commonind(a, b) = (cs = commoninds(a, b); isempty(cs) ? nothing : first(cs)) noncommonind(a, b) = (us = uniqueinds(a, b); isempty(us) ? nothing : first(us)) -noncommoninds(a, b) = namesymdiff(_compat_inds(a), _compat_inds(b)) replaceind(t, p::Pair) = replaceinds(t, p) replaceind(t, from::Index, to::Index) = replaceinds(t, from => to) @@ -142,7 +120,7 @@ const svd_trunc = MAK.svd_trunc function qr(a::AbstractITensor, linds...) left = cat_inds(linds...) - right = namesetdiff(inds(a), left) + right = setdiff(inds(a), left) return MAK.qr_compact(a, Tuple(left), Tuple(right)) end diff --git a/src/sampling.jl b/src/sampling.jl index caec56b..b1cce79 100644 --- a/src/sampling.jl +++ b/src/sampling.jl @@ -203,7 +203,7 @@ function get_one_sample( # 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)...)) + aux = setdiff(inds(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 From 3ca00104de1109c63737d579753e0577ba3de72b Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Thu, 9 Jul 2026 12:29:31 -0400 Subject: [PATCH 07/21] Factor through MatrixAlgebraKit-native decompositions instead of the compat helpers Replaces the hand-rolled factorization compat layer with direct MatrixAlgebraKit-native calls. Adds operator_inds, the codomain/domain prime-level bipartition of an operator, and takes the operator trace through it with tr(a, operator_inds(a)...) rather than contracting against an identity. factorize now wraps left_orth / right_orth, which cover both the exact and truncating cases through trunc, and eigen wraps eigh_full. Drops the qr / svd / svd_trunc / itensor_tr shims and the index-bipartition helpers in favor of MatrixAlgebraKit calls and Base set functions at the callsites. Co-Authored-By: Claude Opus 4.8 (1M context) --- src/Apply/full_update.jl | 5 +- src/Apply/simple_update.jl | 6 +- src/MessagePassing/boundarympscache.jl | 4 +- src/itensors.jl | 95 +++++--------------------- src/rdm.jl | 2 +- src/sampling.jl | 6 +- src/symmetric_gauge.jl | 2 +- 7 files changed, 31 insertions(+), 89 deletions(-) diff --git a/src/Apply/full_update.jl b/src/Apply/full_update.jl index 0d90e41..8bcb8e2 100644 --- a/src/Apply/full_update.jl +++ b/src/Apply/full_update.jl @@ -40,7 +40,7 @@ function full_update( 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, S, V = MAK.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) @@ -114,8 +114,9 @@ function optimise_p_q( envisposdef = true, apply_kwargs..., ) + pq = apply(o, p * q) p_cur, q_cur = factorize( - apply(o, p * q), inds(p); tags = tags(commonind(p, q)), apply_kwargs... + pq, intersect(inds(pq), inds(p)); tags = tags(commonind(p, q)), apply_kwargs... ) fstart = print_fidelity_loss ? fidelity(envs, p_cur, q_cur, p, q, o) : 0 diff --git a/src/Apply/simple_update.jl b/src/Apply/simple_update.jl index cbc2c79..9c24c05 100644 --- a/src/Apply/simple_update.jl +++ b/src/Apply/simple_update.jl @@ -50,8 +50,8 @@ function simple_update( ψᵥ₂ = contract([ψ⃗[2]; sqrt_envs_v2]) sᵥ₁ = commoninds(ψ⃗[1], o) sᵥ₂ = commoninds(ψ⃗[2], o) - Qᵥ₁, Rᵥ₁ = qr(ψᵥ₁, setdiff(uniqueinds(ψᵥ₁, ψᵥ₂), sᵥ₁)) - Qᵥ₂, Rᵥ₂ = qr(ψᵥ₂, setdiff(uniqueinds(ψᵥ₂, ψᵥ₁), sᵥ₂)) + Qᵥ₁, Rᵥ₁ = MAK.qr_compact(ψᵥ₁, setdiff(uniqueinds(ψᵥ₁, ψᵥ₂), sᵥ₁)) + Qᵥ₂, Rᵥ₂ = MAK.qr_compact(ψᵥ₂, setdiff(uniqueinds(ψᵥ₂, ψᵥ₁), sᵥ₂)) rᵥ₁ = commoninds(Qᵥ₁, Rᵥ₁) rᵥ₂ = commoninds(Qᵥ₂, Rᵥ₂) oR = apply(o, Rᵥ₁ * Rᵥ₂) @@ -59,7 +59,7 @@ function simple_update( # 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, union(rᵥ₁, sᵥ₁); trunc = itensor_trunc(; apply_kwargs...)) + U, S, V = MAK.svd_trunc(oR, union(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) diff --git a/src/MessagePassing/boundarympscache.jl b/src/MessagePassing/boundarympscache.jl index 36162b0..1bc6b23 100644 --- a/src/MessagePassing/boundarympscache.jl +++ b/src/MessagePassing/boundarympscache.jl @@ -430,7 +430,7 @@ function generic_apply( end keep = left_link === nothing ? Index[site...] : Index[site..., left_link] - L, R = factorize(T, keep...; ortho = "left", cutoff, maxdim, tags = "Link,l=$i") + L, R = factorize(T, keep; ortho = "left", cutoff, maxdim, tags = "Link,l=$i") push!(out, L) carry = R left_link = only(commoninds(L, R)) @@ -441,7 +441,7 @@ function generic_apply( # Back sweep: right-to-left SVD recompression (optimal truncation of the forward result). for i in length(out):-1:2 bond = only(commoninds(out[i - 1], out[i])) - L, R = factorize(out[i], bond; ortho = "right", cutoff, maxdim, tags = "Link,l=$(i - 1)") + L, R = factorize(out[i], [bond]; ortho = "right", cutoff, maxdim, tags = "Link,l=$(i - 1)") out[i] = R out[i - 1] *= L end diff --git a/src/itensors.jl b/src/itensors.jl index 74dfbd4..65eade2 100644 --- a/src/itensors.jl +++ b/src/itensors.jl @@ -49,10 +49,6 @@ noncommonind(a, b) = (us = uniqueinds(a, b); isempty(us) ? nothing : first(us)) replaceind(t, p::Pair) = replaceinds(t, p) 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)) - const _IndexColl = Union{Tuple{Vararg{Index}}, AbstractVector{<:Index}} function replaceinds(t, pairs::Pair...) return replacedimnames(t, map(p -> name(first(p)) => name(last(p)), pairs)...) @@ -110,64 +106,31 @@ delta(is::Tuple) = delta(Float64, is) delta(is::Index...) = delta(Float64, is) delta(is::AbstractVector{<:Index}) = delta(Float64, Tuple(is)) -function itensor_tr(t::AbstractITensor) - unprimed = filter(i -> plev(i) == 0, inds(t)) - codomain, domain = conj.(unprimed), conj.(prime.(unprimed)) - return scalar(t * one(similar_map(t, codomain, domain), codomain, domain)) -end - -const svd_trunc = MAK.svd_trunc - -function qr(a::AbstractITensor, linds...) - left = cat_inds(linds...) - right = setdiff(inds(a), left) - return MAK.qr_compact(a, Tuple(left), Tuple(right)) +# The codomain/domain bipartition of an operator tensor: each plev-0 index paired with its +# prime. Viewing the operator as this square map is what `tr` and `eigen` factor through. +function operator_inds(a::AbstractITensor) + domain = filter(i -> plev(i) == 0, inds(a)) + return prime.(domain), domain end apply(o::AbstractITensor, ψ::AbstractITensor) = noprime(o * ψ) -_astuple(x::Index) = (x,) -_astuple(x) = Tuple(x) - -function svd(a::AbstractITensor, codomain; kwargs...) - return MAK.svd_compact(a, _astuple(codomain); kwargs...) -end -function svd(a::AbstractITensor, codomain, domain; kwargs...) - return MAK.svd_compact(a, _astuple(codomain), _astuple(domain); kwargs...) -end - -function _eigh( - m::AbstractITensor, - codomain, - domain; - ishermitian = false, - cutoff = nothing - ) +function eigen(m::AbstractITensor, codomain, domain; ishermitian = false, cutoff = nothing) ishermitian || error("the compat `eigen` only supports the hermitian case (ishermitian = true)") isnothing(cutoff) || error( "the compat `eigen` does not yet translate the `cutoff` truncation kwarg to MatrixAlgebraKit's `trunc` spec" ) - cod, dom = _astuple(codomain), _astuple(domain) - D, U = MAK.eigh_full(m, cod, dom) + D, U = MAK.eigh_full(m, codomain, domain) u = only(commoninds(D, U)) t = only(uniqueinds(D, U)) D = replaceinds(D, t => ITensorBase.prime(u)) return D, U end -function eigen(m::AbstractITensor, codomain, domain; kwargs...) - return _eigh(m, codomain, domain; kwargs...) -end - function eigen(m::AbstractITensor; kwargs...) - is = collect(inds(m)) - p0 = filter(i -> plev(i) == 0, is) - p1 = filter(i -> plev(i) != 0, is) - length(p0) == length(p1) || error( - "`eigen` without an index partition expects each plev-0 index to be paired with its prime" - ) - D, U = _eigh(m, Tuple(p1), Tuple(p0); kwargs...) + codomain, domain = operator_inds(m) + D, U = eigen(m, codomain, domain; kwargs...) return D, conj(U) end @@ -178,45 +141,23 @@ function itensor_trunc(; maxdim = nothing, cutoff = nothing) return trunc end -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 - function factorize( a::AbstractITensor, - linds...; + codomain; ortho = "left", cutoff = nothing, maxdim = nothing, tags = nothing ) - 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" - L, R = MAK.lq_compact(a, Tuple(left), Tuple(right)) - else - error( - "compat `factorize` supports ortho = \"left\" / \"right\" (got $(repr(ortho)))" - ) - end + # `left_orth` / `right_orth` take the codomain indices and infer the domain, and `trunc` + # covers both the exact (no cutoff/maxdim) and truncating cases. + trunc = itensor_trunc(; cutoff, maxdim) + if ortho == "left" + L, R = MAK.left_orth(a, codomain; trunc) + elseif ortho == "right" + L, R = MAK.right_orth(a, codomain; trunc) else - 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" - L, R = U * S, Vt - else - error( - "compat `factorize` supports ortho = \"left\" / \"right\" (got $(repr(ortho)))" - ) - end + error("compat `factorize` supports ortho = \"left\" / \"right\" (got $(repr(ortho)))") end if !isnothing(tags) b = only(commoninds(L, R)) diff --git a/src/rdm.jl b/src/rdm.jl index 0ec72bd..2ecb13a 100644 --- a/src/rdm.jl +++ b/src/rdm.jl @@ -1,5 +1,5 @@ function normalize_rdm(ρ::ITensor) - return ρ / itensor_tr(ρ) + return ρ / tr(ρ, operator_inds(ρ)...) end """ diff --git a/src/sampling.jl b/src/sampling.jl index b1cce79..4370615 100644 --- a/src/sampling.jl +++ b/src/sampling.jl @@ -25,7 +25,7 @@ function sample( seq = contraction_sequence(tensors; alg = "optimal") ρ = contract(tensors; sequence = seq) - ρ_tr = itensor_tr(ρ) + ρ_tr = tr(ρ, operator_inds(ρ)...) ρ *= inv(ρ_tr) ρ_diag = collect(real.(diag(array(ρ)))) config = StatsBase.sample(1:length(ρ_diag), Weights(ρ_diag)) @@ -203,7 +203,7 @@ function get_one_sample( # 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 = setdiff(inds(mt), cat_inds(net_inds, (inds(m) for m in outgoing_mps if m !== mt)...)) + aux = setdiff(inds(mt), net_inds, (inds(m) for m in outgoing_mps if m !== mt)...) setmessage!(norm_bmps_cache, e, ITensor[mt, bra_tensor(mt, aux)]) end @@ -238,7 +238,7 @@ function sample_partition!( ts = [incoming_ms; [ψv, ψvdag]] seq = contraction_sequence(ts; alg = "optimal") ρ = contract(ts; sequence = seq) - ρ_tr = itensor_tr(ρ) + ρ_tr = tr(ρ, operator_inds(ρ)...) push!(traces, ρ_tr) ρ *= inv(ρ_tr) ρ_diag = collect(real.(diag(array(ρ)))) diff --git a/src/symmetric_gauge.jl b/src/symmetric_gauge.jl index 5a9c699..225aa01 100644 --- a/src/symmetric_gauge.jl +++ b/src/symmetric_gauge.jl @@ -35,7 +35,7 @@ function symmetric_gauge!(bp_cache::BeliefPropagationCache; regularization = 10 # 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, S, V = MAK.svd_compact(Ce, edge_ind_p; kwargs...) u, v = commoninds(S, U), commoninds(S, V) ψvsrc = replaceinds(ψvsrc, edge_ind, edge_ind_p) * U From c29d08b0d582fc2cecb9b9ba4f014982443bb8a1 Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Thu, 9 Jul 2026 13:13:54 -0400 Subject: [PATCH 08/21] Drop the local index and similar shims for ITensorBase equivalents Removes the lenient `commonind`/`noncommonind` and the `similar_map` wrapper. Call sites now use ITensorBase's `trycommonind`/`trynoncommonind` (the single shared or unique index, or `nothing`) and the map-shaped `similar`, which takes its axes as a vector directly. The boundary-MPS link bonds carry a single `link` tag keyed by the bond index instead of the old two-tag `Link,l=...` string. Co-Authored-By: Claude Opus 4.8 (1M context) --- src/Apply/apply_gates.jl | 4 ++-- src/Apply/full_update.jl | 4 ++-- src/Forms/abstractform.jl | 2 +- src/Forms/bilinearform.jl | 2 +- src/MessagePassing/boundarympscache.jl | 4 ++-- src/itensors.jl | 19 ++----------------- 6 files changed, 10 insertions(+), 25 deletions(-) diff --git a/src/Apply/apply_gates.jl b/src/Apply/apply_gates.jl index ea0bd84..0393add 100644 --- a/src/Apply/apply_gates.jl +++ b/src/Apply/apply_gates.jl @@ -111,8 +111,8 @@ function apply_gate!( if length(v⃗) == 2 v1, v2 = v⃗ e = NamedEdge(v1 => v2) - u = commonind(s_values, first(updated_tensors)) - v = noncommonind(s_values, first(updated_tensors)) + u = trycommonind(s_values, first(updated_tensors)) + v = trynoncommonind(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 diff --git a/src/Apply/full_update.jl b/src/Apply/full_update.jl index 8bcb8e2..8700659 100644 --- a/src/Apply/full_update.jl +++ b/src/Apply/full_update.jl @@ -66,7 +66,7 @@ function fidelity( q_prev::ITensor, gate::ITensor, ) - p_sind, q_sind = commonind(p_cur, gate), commonind(q_cur, gate) + p_sind, q_sind = trycommonind(p_cur, gate), trycommonind(q_cur, gate) p_sind_sim, q_sind_sim = sim(p_sind), sim(q_sind) gate_sq = gate * replaceinds(conj(gate), Index[p_sind, q_sind], Index[p_sind_sim, q_sind_sim]) @@ -116,7 +116,7 @@ function optimise_p_q( ) pq = apply(o, p * q) p_cur, q_cur = factorize( - pq, intersect(inds(pq), inds(p)); tags = tags(commonind(p, q)), apply_kwargs... + pq, intersect(inds(pq), inds(p)); tags = tags(trycommonind(p, q)), apply_kwargs... ) fstart = print_fidelity_loss ? fidelity(envs, p_cur, q_cur, p, q, o) : 0 diff --git a/src/Forms/abstractform.jl b/src/Forms/abstractform.jl index 71b0c41..4a667af 100644 --- a/src/Forms/abstractform.jl +++ b/src/Forms/abstractform.jl @@ -24,7 +24,7 @@ end function default_message(form::AbstractForm, edge::AbstractEdge) cod = virtualinds(ket(form), edge) dom = conj.(bra_virtualinds(form, edge)) - return one(similar_map(ket(form)[src(edge)], cod, dom), cod, dom) + return one(similar(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 691f3ee..dbe7f2e 100644 --- a/src/Forms/bilinearform.jl +++ b/src/Forms/bilinearform.jl @@ -20,7 +20,7 @@ function BilinearForm(ket::TensorNetworkState, bra::TensorNetworkState) bra = TensorNetworkState(Dictionary(verts, [bra_tensor(bra, v) for v in verts])) operator_tensors = [ let codomain = conj.(sinds[v]), domain = conj.(prime.(sinds[v])) - one(similar_map(ket[v], codomain, domain), codomain, domain) + one(similar(ket[v], codomain, domain), codomain, domain) end for v in verts ] operator = TensorNetworkState(Dictionary(verts, operator_tensors)) diff --git a/src/MessagePassing/boundarympscache.jl b/src/MessagePassing/boundarympscache.jl index 1bc6b23..3ee3f9a 100644 --- a/src/MessagePassing/boundarympscache.jl +++ b/src/MessagePassing/boundarympscache.jl @@ -430,7 +430,7 @@ function generic_apply( end keep = left_link === nothing ? Index[site...] : Index[site..., left_link] - L, R = factorize(T, keep; ortho = "left", cutoff, maxdim, tags = "Link,l=$i") + L, R = factorize(T, keep; ortho = "left", cutoff, maxdim, tags = "link" => "$i") push!(out, L) carry = R left_link = only(commoninds(L, R)) @@ -441,7 +441,7 @@ function generic_apply( # Back sweep: right-to-left SVD recompression (optimal truncation of the forward result). for i in length(out):-1:2 bond = only(commoninds(out[i - 1], out[i])) - L, R = factorize(out[i], [bond]; ortho = "right", cutoff, maxdim, tags = "Link,l=$(i - 1)") + L, R = factorize(out[i], [bond]; ortho = "right", cutoff, maxdim, tags = "link" => "$(i - 1)") out[i] = R out[i - 1] *= L end diff --git a/src/itensors.jl b/src/itensors.jl index 65eade2..648abac 100644 --- a/src/itensors.jl +++ b/src/itensors.jl @@ -8,7 +8,7 @@ import TensorAlgebra: matricize using Adapt: Adapt using ITensorBase: ITensorBase, AbstractITensor, ITensor, Index, NamedUnitRange, commoninds, dimnames, hascommoninds, id, inds, name, nameddims, noncommoninds, noprime, plev, prime, - replacedimnames, sim, tags, unioninds, unnamed + replacedimnames, sim, tags, trycommonind, trynoncommonind, unioninds, unnamed using LinearAlgebra: LinearAlgebra using TensorAlgebra: TensorAlgebra, project, scalar, tryproject @@ -40,12 +40,6 @@ end _contract_sequence(tensors, s::Integer) = tensors[s] _contract_sequence(tensors, s) = reduce(*, (_contract_sequence(tensors, x) for x in s)) -# The single shared/unique index, or `nothing` when there is not exactly one. `commoninds` -# and `uniqueinds` come from ITensorBase, which keys index-set algebra by name so a graded -# bond still matches its dual. -commonind(a, b) = (cs = commoninds(a, b); isempty(cs) ? nothing : first(cs)) -noncommonind(a, b) = (us = uniqueinds(a, b); isempty(us) ? nothing : first(us)) - replaceind(t, p::Pair) = replaceinds(t, p) replaceind(t, from::Index, to::Index) = replaceinds(t, from => to) @@ -89,16 +83,6 @@ function diagonaltensor( return nameddims(diagonaltensor(diag, unnamed.(is)), name.(is)) end -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 - 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)) @@ -235,6 +219,7 @@ function settags(i::Index, tagstr::AbstractString) end return i end +settags(i::Index, p::Pair) = ITensorBase.settag(i, first(p), last(p)) function settags(i::Index, d::AbstractDict) for (k, v) in d i = ITensorBase.settag(i, k, v) From 437d2922eeb3c0c8903136cc8afacefc7e4d987a Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Thu, 9 Jul 2026 13:31:54 -0400 Subject: [PATCH 09/21] Label contraction sequences by Index instead of name Now that `Index` equality is name-based, the sequence optimizers can take the `Index`es directly rather than their names. `to_eincode`, the omeinsum path, and the `optimal` dimension dict drop the `name.(...)` conversion, and a shared leg still matches across its two tensors even when stored nondual on one and dual on the other. The test compares against the `Index`es themselves. Co-Authored-By: Claude Opus 4.8 (1M context) --- src/contraction_sequences.jl | 17 ++++++++--------- test/test_contraction_sequences.jl | 13 +++++++------ 2 files changed, 15 insertions(+), 15 deletions(-) diff --git a/src/contraction_sequences.jl b/src/contraction_sequences.jl index ee32d94..bf27c14 100644 --- a/src/contraction_sequences.jl +++ b/src/contraction_sequences.jl @@ -9,13 +9,12 @@ using OMEinsumContractionOrders: OMEinsumContractionOrders, optimize_code, EinCo # untouched originals), so no placeholder tensor is needed. is_trivial_tensor(t::ITensor) = all(i -> length(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). +# The sequence optimizers use indices as opaque labels (plus their dimension). `Index` +# equality is name-based, so a shared leg matches across its two tensors even when it is +# stored nondual on one and dual on the other under a graded backend. function contraction_network(tensors::Vector{<:ITensor}; prune_tensors = false) return map(tensors) do t - is = collect(name.(inds(t))) + is = collect(inds(t)) (prune_tensors && is_trivial_tensor(t)) ? empty(is) : is end end @@ -23,7 +22,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(name(i) => Float64(length(i)) for t in tensors for i in inds(t)) + inds_to_dims = Dict(i => Float64(length(i)) for t in tensors for i in inds(t)) seq, _ = optimaltree(network, inds_to_dims) seq = typeof(seq) <: Int ? [seq] : seq return seq @@ -41,10 +40,10 @@ end #OMEinsumContractionOrders helpers function to_eincode(tensors::Vector{<:ITensor}) - ixs = map(t -> collect(name.(inds(t))), tensors) + ixs = map(t -> collect(inds(t)), tensors) LT = eltype(eltype(ixs)) - iy = collect(LT, name.(reduce(symdiff, inds.(tensors)))) - size_dict = Dict{LT, Int}(name(i) => length(i) for t in tensors for i in inds(t)) + iy = collect(LT, reduce(symdiff, inds.(tensors))) + size_dict = Dict{LT, Int}(i => length(i) for t in tensors for i in inds(t)) return EinCode(ixs, iy), size_dict end diff --git a/test/test_contraction_sequences.jl b/test/test_contraction_sequences.jl index 2f221da..cfa6de7 100644 --- a/test/test_contraction_sequences.jl +++ b/test/test_contraction_sequences.jl @@ -1,5 +1,5 @@ @eval module $(gensym()) -using ITensorBase: Index, name +using ITensorBase: Index using Random using TensorNetworkQuantumSimulator const TNQS = TensorNetworkQuantumSimulator @@ -17,15 +17,16 @@ 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. + # Labels are the `Index`es themselves: `Index` equality is name-based, so a shared + # leg matches across its two tensors even when stored nondual on one and dual on the + # other under a graded backend. i, j, k = Index(2), Index(3), Index(4) A = randn(i, j) B = randn(j, k) code, size_dict = TNQS.to_eincode([A, B]) - @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) + @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) # --- to_contraction_sequence: NestedEinsum -> nested tensor-position tree. Tests our # converter on hand-built trees with known shapes (deterministic, exact). From e03e116c58a5714457ac92dcb3c38f42511c2aee Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Thu, 9 Jul 2026 13:37:51 -0400 Subject: [PATCH 10/21] Compare tensor-network indices by Index equality, not name With `Index` equality name-based, the `uniqueinds(tn, v)` graph-index difference and the loop-correction bond lookup drop their explicit `name(...)` comparisons for `Index` set-ops and equality. A shared graded link, stored nondual on one endpoint and dual on the other, still matches. Co-Authored-By: Claude Opus 4.8 (1M context) --- src/MessagePassing/loopcorrection.jl | 2 +- src/TensorNetworks/abstracttensornetwork.jl | 11 +++++------ 2 files changed, 6 insertions(+), 7 deletions(-) diff --git a/src/MessagePassing/loopcorrection.jl b/src/MessagePassing/loopcorrection.jl index 8a22cfc..6e7acdd 100644 --- a/src/MessagePassing/loopcorrection.jl +++ b/src/MessagePassing/loopcorrection.jl @@ -48,7 +48,7 @@ function sim_edgeinduced_subgraph(bpc::BeliefPropagationCache, eg) t_inds = intersect(inds(t), linds) if !isempty(t_inds) t_ind = only(t_inds) - t_ind_pos = findfirst(x -> name(x) == name(t_ind), linds) + t_ind_pos = findfirst(==(t_ind), linds) t = replaceind(t, t_ind, linds_sim[t_ind_pos]) setindex_preserve!(bpc, t, src(e)) end diff --git a/src/TensorNetworks/abstracttensornetwork.jl b/src/TensorNetworks/abstracttensornetwork.jl index 99c0600..291819a 100644 --- a/src/TensorNetworks/abstracttensornetwork.jl +++ b/src/TensorNetworks/abstracttensornetwork.jl @@ -27,15 +27,14 @@ function maxvirtualdim(tn::AbstractTensorNetwork) return maximum(maximum.([length.(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. +# `setdiff` is by-name (dual-insensitive), so a shared graded link, stored nondual on one +# endpoint and dual on the other, is still excluded here. function uniqueinds(tn::AbstractTensorNetwork, v) - tv_inds = Index[i for i in inds(tn[v])] + tv_inds = collect(inds(tn[v])) vns = neighbors(tn, v) isempty(vns) && return tv_inds - 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) + neighbor_inds = reduce(vcat, [collect(inds(tn[vn])) for vn in vns]) + return setdiff(tv_inds, neighbor_inds) end function setindex_preserve!(tn::AbstractTensorNetwork, value::ITensor, vertex) From 8f8062e5848853070304c5e35ad2a509470751cd Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Thu, 9 Jul 2026 13:40:43 -0400 Subject: [PATCH 11/21] Drop redundant collect around inds `inds` already returns a fresh `Vector`, so the `collect` calls in `to_eincode`, `contraction_network`, and `uniqueinds` were redundant. Use the `inds` result directly. Co-Authored-By: Claude Opus 4.8 (1M context) --- src/TensorNetworks/abstracttensornetwork.jl | 4 ++-- src/contraction_sequences.jl | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/TensorNetworks/abstracttensornetwork.jl b/src/TensorNetworks/abstracttensornetwork.jl index 291819a..a39db4e 100644 --- a/src/TensorNetworks/abstracttensornetwork.jl +++ b/src/TensorNetworks/abstracttensornetwork.jl @@ -30,10 +30,10 @@ end # `setdiff` is by-name (dual-insensitive), so a shared graded link, stored nondual on one # endpoint and dual on the other, is still excluded here. function uniqueinds(tn::AbstractTensorNetwork, v) - tv_inds = collect(inds(tn[v])) + tv_inds = inds(tn[v]) vns = neighbors(tn, v) isempty(vns) && return tv_inds - neighbor_inds = reduce(vcat, [collect(inds(tn[vn])) for vn in vns]) + neighbor_inds = reduce(vcat, [inds(tn[vn]) for vn in vns]) return setdiff(tv_inds, neighbor_inds) end diff --git a/src/contraction_sequences.jl b/src/contraction_sequences.jl index bf27c14..cac3eb3 100644 --- a/src/contraction_sequences.jl +++ b/src/contraction_sequences.jl @@ -14,7 +14,7 @@ is_trivial_tensor(t::ITensor) = all(i -> length(i) == 1, inds(t)) # stored nondual on one and dual on the other under a graded backend. function contraction_network(tensors::Vector{<:ITensor}; prune_tensors = false) return map(tensors) do t - is = collect(inds(t)) + is = inds(t) (prune_tensors && is_trivial_tensor(t)) ? empty(is) : is end end @@ -40,7 +40,7 @@ end #OMEinsumContractionOrders helpers function to_eincode(tensors::Vector{<:ITensor}) - ixs = map(t -> collect(inds(t)), tensors) + ixs = map(inds, tensors) LT = eltype(eltype(ixs)) iy = collect(LT, reduce(symdiff, inds.(tensors))) size_dict = Dict{LT, Int}(i => length(i) for t in tensors for i in inds(t)) From 57593023d4b5ed07d80cd4c8190e045d288568b1 Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Thu, 9 Jul 2026 14:03:31 -0400 Subject: [PATCH 12/21] Retire the compat factorize and eigen wrappers Factorization call sites go straight through MatrixAlgebraKit (`left_orth`/`right_orth`), keeping `itensor_trunc` for the `cutoff`/`maxdim` to `trunc` translation. Where a call named the new bond, the retag is inlined at the call site for now, until tag support lands in the MatrixAlgebraKit/ITensorBase factorizations. The compat `eigen` and the `safe_eigen` precision wrapper around it had no live callers and are deleted. Co-Authored-By: Claude Opus 4.8 (1M context) --- src/Apply/full_update.jl | 11 +++--- src/MessagePassing/boundarympscache.jl | 12 +++++-- src/itensors.jl | 47 +------------------------- src/utils.jl | 20 ----------- 4 files changed, 17 insertions(+), 73 deletions(-) diff --git a/src/Apply/full_update.jl b/src/Apply/full_update.jl index 8700659..5fdd313 100644 --- a/src/Apply/full_update.jl +++ b/src/Apply/full_update.jl @@ -18,10 +18,10 @@ function full_update( apply_kwargs..., ) - Qᵥ₁, Rᵥ₁ = factorize( + Qᵥ₁, Rᵥ₁ = MAK.left_orth( ψ[v⃗[1]], setdiff(uniqueinds(ψ[v⃗[1]], ψ[v⃗[2]]), uniqueinds(ψ, v⃗[1])) ) - Qᵥ₂, Rᵥ₂ = factorize( + Qᵥ₂, Rᵥ₂ = MAK.left_orth( ψ[v⃗[2]], setdiff(uniqueinds(ψ[v⃗[2]], ψ[v⃗[1]]), uniqueinds(ψ, v⃗[2])) ) @@ -115,9 +115,12 @@ function optimise_p_q( apply_kwargs..., ) pq = apply(o, p * q) - p_cur, q_cur = factorize( - pq, intersect(inds(pq), inds(p)); tags = tags(trycommonind(p, q)), apply_kwargs... + p_cur, q_cur = MAK.left_orth( + pq, intersect(inds(pq), inds(p)); trunc = itensor_trunc(; apply_kwargs...) ) + b = only(commoninds(p_cur, q_cur)) + bnew = settags(b, tags(trycommonind(p, q))) + p_cur, q_cur = replaceind(p_cur, b, bnew), replaceind(q_cur, b, bnew) fstart = print_fidelity_loss ? fidelity(envs, p_cur, q_cur, p, q, o) : 0 diff --git a/src/MessagePassing/boundarympscache.jl b/src/MessagePassing/boundarympscache.jl index 3ee3f9a..32a97e3 100644 --- a/src/MessagePassing/boundarympscache.jl +++ b/src/MessagePassing/boundarympscache.jl @@ -282,7 +282,7 @@ function gauge_step!( m1, m2 = message(bmps_cache, e1), message(bmps_cache, e2) @assert !isempty(commoninds(m1, m2)) left_inds = uniqueinds(m1, m2) - m1, Y = factorize(m1, left_inds; ortho = "left", kwargs...) + m1, Y = MAK.left_orth(m1, left_inds; trunc = itensor_trunc(; kwargs...)) m2 = m2 * Y setmessage!(bmps_cache, e1, m1) setmessage!(bmps_cache, e2, m2) @@ -430,7 +430,10 @@ function generic_apply( end keep = left_link === nothing ? Index[site...] : Index[site..., left_link] - L, R = factorize(T, keep; ortho = "left", cutoff, maxdim, tags = "link" => "$i") + L, R = MAK.left_orth(T, keep; trunc = itensor_trunc(; cutoff, maxdim)) + b = only(commoninds(L, R)) + bnew = settags(b, "link" => "$i") + L, R = replaceind(L, b, bnew), replaceind(R, b, bnew) push!(out, L) carry = R left_link = only(commoninds(L, R)) @@ -441,7 +444,10 @@ function generic_apply( # Back sweep: right-to-left SVD recompression (optimal truncation of the forward result). for i in length(out):-1:2 bond = only(commoninds(out[i - 1], out[i])) - L, R = factorize(out[i], [bond]; ortho = "right", cutoff, maxdim, tags = "link" => "$(i - 1)") + L, R = MAK.right_orth(out[i], [bond]; trunc = itensor_trunc(; cutoff, maxdim)) + b = only(commoninds(L, R)) + bnew = settags(b, "link" => "$(i - 1)") + L, R = replaceind(L, b, bnew), replaceind(R, b, bnew) out[i] = R out[i - 1] *= L end diff --git a/src/itensors.jl b/src/itensors.jl index 648abac..73788ee 100644 --- a/src/itensors.jl +++ b/src/itensors.jl @@ -91,7 +91,7 @@ delta(is::Index...) = delta(Float64, is) delta(is::AbstractVector{<:Index}) = delta(Float64, Tuple(is)) # The codomain/domain bipartition of an operator tensor: each plev-0 index paired with its -# prime. Viewing the operator as this square map is what `tr` and `eigen` factor through. +# prime. Viewing the operator as this square map is what `tr` factors through. function operator_inds(a::AbstractITensor) domain = filter(i -> plev(i) == 0, inds(a)) return prime.(domain), domain @@ -99,25 +99,6 @@ end apply(o::AbstractITensor, ψ::AbstractITensor) = noprime(o * ψ) -function eigen(m::AbstractITensor, codomain, domain; ishermitian = false, cutoff = nothing) - ishermitian || - error("the compat `eigen` only supports the hermitian case (ishermitian = true)") - isnothing(cutoff) || error( - "the compat `eigen` does not yet translate the `cutoff` truncation kwarg to MatrixAlgebraKit's `trunc` spec" - ) - D, U = MAK.eigh_full(m, codomain, domain) - u = only(commoninds(D, U)) - t = only(uniqueinds(D, U)) - D = replaceinds(D, t => ITensorBase.prime(u)) - return D, U -end - -function eigen(m::AbstractITensor; kwargs...) - codomain, domain = operator_inds(m) - D, U = eigen(m, codomain, domain; kwargs...) - return D, conj(U) -end - function itensor_trunc(; maxdim = nothing, cutoff = nothing) trunc = MAK.notrunc() isnothing(maxdim) || (trunc &= MAK.truncrank(maxdim)) @@ -125,32 +106,6 @@ function itensor_trunc(; maxdim = nothing, cutoff = nothing) return trunc end -function factorize( - a::AbstractITensor, - codomain; - ortho = "left", - cutoff = nothing, - maxdim = nothing, - tags = nothing - ) - # `left_orth` / `right_orth` take the codomain indices and infer the domain, and `trunc` - # covers both the exact (no cutoff/maxdim) and truncating cases. - trunc = itensor_trunc(; cutoff, maxdim) - if ortho == "left" - L, R = MAK.left_orth(a, codomain; trunc) - elseif ortho == "right" - L, R = MAK.right_orth(a, codomain; trunc) - else - error("compat `factorize` supports ortho = \"left\" / \"right\" (got $(repr(ortho)))") - end - 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 -end - datatype(T::AbstractITensor) = typeof(unnamed(T)) array(T::AbstractITensor) = convert(Array, unnamed(T)) data(T::AbstractITensor) = unnamed(T) diff --git a/src/utils.jl b/src/utils.jl index 790f07b..2fec555 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -69,26 +69,6 @@ function with_default_maxiter(cache_update_kwargs, cache) return (; cache_update_kwargs..., maxiter) end -""" - safe_eigen(m::ITensor, args...; kwargs...) - A wrapper around eigen that ensures eigen computations are done in Float64/ComplexF64 precision on CPU for better numerical stability. -""" -function safe_eigen(m::ITensor, args...; kwargs...) - dtype = datatype(m) - e = eltype(m) - if e == ComplexF64 || e == Float64 - return eigen(m, args...; kwargs...) - elseif e == Float32 - m = adapt(Vector{Float64}, m) - D, U = eigen(m, args...; kwargs...) - return adapt(dtype)(D), adapt(dtype)(U) - elseif e == ComplexF32 - m = adapt(Vector{ComplexF64}, m) - D, U = eigen(m, args...; kwargs...) - return adapt(dtype)(D), adapt(dtype)(U) - end -end - collect_vertices(e::NamedEdge, g::NamedGraph) = collect_vertices([src(e), dst(e)], g) collect_vertices(es::Vector{<:NamedEdge}, g::NamedGraph) = reduce(vcat, [collect_vertices(e, g) for e in es]) From 6a641f042a4507b5298528bb8302ed919444c1e3 Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Thu, 9 Jul 2026 14:07:06 -0400 Subject: [PATCH 13/21] Delete the dead hastags helper Nothing reads tags back (no `hastags`/`gettag` lookups anywhere), so the tag-string membership helper had no callers. Co-Authored-By: Claude Opus 4.8 (1M context) --- src/itensors.jl | 5 ----- 1 file changed, 5 deletions(-) diff --git a/src/itensors.jl b/src/itensors.jl index 73788ee..4876ddf 100644 --- a/src/itensors.jl +++ b/src/itensors.jl @@ -181,8 +181,3 @@ function settags(i::Index, d::AbstractDict) end return i end -function hastags(i::Index, tagstr::AbstractString) - return all( - haskey(tags(i), String(strip(t))) for t in split(tagstr, ",") if !isempty(strip(t)) - ) -end From 32b3252d530ca94b14d33293b1810ef79e6447d5 Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Thu, 9 Jul 2026 14:15:02 -0400 Subject: [PATCH 14/21] Simplify directsum onto Index equality `directsum` matches each tensor's axes to the target order via `findfirst(==(o), inds(t))`, since `Index` equality is name-based, dropping the `name`/`dimnames` detour. The defensive `Tuple`/`collect` wrapping is removed (`inds` returns a fresh `Vector` and the tuple splices work on any collection) and `shared` uses `setdiff`. Equivalence with the previous implementation was verified directly. Co-Authored-By: Claude Opus 4.8 (1M context) --- src/itensors.jl | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/src/itensors.jl b/src/itensors.jl index 4876ddf..e050a08 100644 --- a/src/itensors.jl +++ b/src/itensors.jl @@ -126,17 +126,15 @@ hasqns(t::AbstractITensor) = any(hasqns, inds(t)) hasqns(::Any) = false function directsum(out_inds, pairs::Pair...) - out_inds = Tuple(out_inds) - t1, s1 = first(pairs[1]), Tuple(last(pairs[1])) - shared = Tuple(filter(i -> !(i in s1), collect(inds(t1)))) + t1, s1 = first(pairs[1]), last(pairs[1]) + shared = setdiff(inds(t1), s1) target = (shared..., out_inds...) out = zeros(scalartype(t1), length.(target)) offsets = zeros(Int, length(out_inds)) for p in pairs - t, sinds = first(p), Tuple(last(p)) + t, sinds = first(p), last(p) order = (shared..., sinds...) - cur = collect(ITensorBase.dimnames(t)) - perm = [findfirst(==(name(o)), cur) for o in order] + perm = [findfirst(==(o), inds(t)) for o in order] a = permutedims(unnamed(t), perm) ranges = ( Base.OneTo.(length.(shared))..., From e80f128e01c9f5a12d14121ccb8107ba9f696088 Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Thu, 9 Jul 2026 14:20:28 -0400 Subject: [PATCH 15/21] Add a test for tensor-network addition Covers the `+`/`add` direct-sum path (the only `directsum` caller), which had no test. Builds two states over a shared set of site indices, adds them, and checks the defining property that contracting the sum equals the sum of the contractions. Co-Authored-By: Claude Opus 4.8 (1M context) --- test/runtests.jl | 1 + test/test_add.jl | 26 ++++++++++++++++++++++++++ 2 files changed, 27 insertions(+) create mode 100644 test/test_add.jl diff --git a/test/runtests.jl b/test/runtests.jl index 1cec2bc..163b6c9 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,6 +3,7 @@ using Test @testset "TensorNetworkQuantumSimulator.jl" begin include("test_constructors.jl") + include("test_add.jl") include("test_forms.jl") include("test_expect.jl") include("test_boundarymps.jl") diff --git a/test/test_add.jl b/test/test_add.jl new file mode 100644 index 0000000..9e712dc --- /dev/null +++ b/test/test_add.jl @@ -0,0 +1,26 @@ +@eval module $(gensym()) +using Random +using TensorNetworkQuantumSimulator +const TNQS = TensorNetworkQuantumSimulator +using Test: @testset, @test + +# Fully contract a tensor-network state to a dense tensor over its site indices. +contract_dense(ψ) = TNQS.contract([ψ[v] for v in vertices(ψ)]) + +@testset "Tensor-network addition (direct sum)" begin + Random.seed!(123) + for elt in (Float64, ComplexF64) + g = named_grid((2, 2)) + ψ1 = random_tensornetworkstate(elt, g; bond_dimension = 2) + s = siteinds(ψ1) + # `ψ2` shares `ψ1`'s site indices so the two states can be added. + ψ2 = random_tensornetworkstate(elt, g, s; bond_dimension = 3) + ψ3 = ψ1 + ψ2 + @test graph(ψ3) == g + @test siteinds(ψ3) == s + # Defining property of the direct sum: contracting the sum equals the sum of the + # contractions. + @test contract_dense(ψ3) ≈ contract_dense(ψ1) + contract_dense(ψ2) + end +end +end From 9ca41b7eed8e3fcb423d8648d983ca5377a8664a Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Thu, 9 Jul 2026 14:23:34 -0400 Subject: [PATCH 16/21] Align directsum with unname and use eltype `directsum` gets each tensor's axes via `unname(t, order)`, which aligns to the requested index order and unnames in one step, replacing the explicit `findfirst` permutation over `permutedims(unnamed(t), ...)`. `eltype` replaces `scalartype` for the output element type. Verified against the tensor-network addition test. Co-Authored-By: Claude Opus 4.8 (1M context) --- src/itensors.jl | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/src/itensors.jl b/src/itensors.jl index e050a08..1b7404e 100644 --- a/src/itensors.jl +++ b/src/itensors.jl @@ -129,13 +129,11 @@ function directsum(out_inds, pairs::Pair...) t1, s1 = first(pairs[1]), last(pairs[1]) shared = setdiff(inds(t1), s1) target = (shared..., out_inds...) - out = zeros(scalartype(t1), length.(target)) + out = zeros(eltype(t1), length.(target)) offsets = zeros(Int, length(out_inds)) for p in pairs t, sinds = first(p), last(p) - order = (shared..., sinds...) - perm = [findfirst(==(o), inds(t)) for o in order] - a = permutedims(unnamed(t), perm) + a = ITensorBase.unname(t, (shared..., sinds...)) ranges = ( Base.OneTo.(length.(shared))..., ntuple(k -> (offsets[k] + 1):(offsets[k] + length(sinds[k])), length(sinds))..., From 7c741b846ccde82f138673ec40e875919b541937 Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Thu, 9 Jul 2026 15:17:35 -0400 Subject: [PATCH 17/21] Call ITensorBase and TensorAlgebra APIs directly, drop more itensors.jl wrappers Continues shrinking `src/itensors.jl` by deleting thin compat wrappers and calling the upstream ITensorBase / TensorAlgebra functions at the call sites instead. Index relabeling uses ITensorBase's `replaceinds` directly (in place of the vendored `replaceind`/`replaceinds` shims, with collection relabels spelled `(from .=> to)...`), storage access uses `unnamed`/`Array` and a `datatype` that now delegates to `TensorAlgebra.datatype` (https://github.com/ITensor/ITensorBase.jl/pull/213), and in-place scaling goes through ITensorBase's `rmul!` for `AbstractNamedTensor`. The list-contraction verb `contract` is renamed to `contract_network`, and tag construction passes a `key => value` pair now that the comma-string `settags` form is gone (`site_tag` returns the site-type tag). `hasqns` stays for contraction-optimizer selection. Co-Authored-By: Claude Opus 4.8 (1M context) --- docs/src/states.md | 6 +-- src/Apply/apply_gates.jl | 2 +- src/Apply/full_update.jl | 26 ++++++------- src/Apply/simple_update.jl | 12 +++--- .../abstractbeliefpropagationcache.jl | 4 +- src/MessagePassing/beliefpropagationcache.jl | 8 ++-- src/MessagePassing/boundarympscache.jl | 12 +++--- src/MessagePassing/loopcorrection.jl | 8 ++-- src/TensorNetworkQuantumSimulator.jl | 2 +- src/TensorNetworks/abstracttensornetwork.jl | 4 +- src/TensorNetworks/tensornetworkstate.jl | 4 +- .../tensornetworkstate_constructors.jl | 4 +- src/contract.jl | 12 +++--- src/entanglement.jl | 6 +-- src/expect.jl | 4 +- src/inner.jl | 2 +- src/itensors.jl | 38 ++----------------- src/norm_sqr.jl | 2 +- src/rdm.jl | 4 +- src/sampling.jl | 8 ++-- src/siteinds.jl | 4 +- src/symmetric_gauge.jl | 18 ++++----- test/test_add.jl | 2 +- test/test_beliefpropagation.jl | 4 +- test/test_boundarymps.jl | 4 +- test/test_constructors.jl | 2 +- test/test_contraction_sequences.jl | 8 ++-- 27 files changed, 90 insertions(+), 120 deletions(-) diff --git a/docs/src/states.md b/docs/src/states.md index d74c0df..0aabc1e 100644 --- a/docs/src/states.md +++ b/docs/src/states.md @@ -16,14 +16,14 @@ t_a, t_b, t_c = randn(i), randn(i, j), randn(j) tn = TensorNetwork(Dictionary(["a", "b", "c"], [t_a, t_b, t_c])) # Contracts to a scalar -z = contract(tn; alg = "exact") +z = contract_network(tn; alg = "exact") # Random tensor network with graph-specified connectivity and given bond dimension g = named_grid((3, 3)) tn = random_tensornetwork(Float64, g; bond_dimension = 4) # Should return a scalar -z = contract(tn; alg = "exact") +z = contract_network(tn; alg = "exact") ``` `TensorNetwork` is a useful type for representing objects like classical partition functions or the solutions to counting problems where you don't need the concept of a physical site index. @@ -34,7 +34,7 @@ The library ships with a built-in constructor for the classical Ising-model part g = named_grid((4, 4)) β = 0.4 Z_tn = ising_partitionfunction(g, β) # uniform J = 1 -Z = contract(Z_tn; alg = "bp") # approximate via belief propagation +Z = contract_network(Z_tn; alg = "bp") # approximate via belief propagation # Anisotropic couplings: pass a Dictionary keyed by edges Js = Dictionary(edges(g), [isodd(src(e)[1]) ? 1.0 : 0.5 for e in edges(g)]) diff --git a/src/Apply/apply_gates.jl b/src/Apply/apply_gates.jl index 0393add..5407378 100644 --- a/src/Apply/apply_gates.jl +++ b/src/Apply/apply_gates.jl @@ -117,7 +117,7 @@ function apply_gate!( # 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)) + s_values = replaceinds(s_values, v => prime(u)) setmessage!(ψ_bpc, e, conj(s_values)) setmessage!(ψ_bpc, reverse(e), s_values) end diff --git a/src/Apply/full_update.jl b/src/Apply/full_update.jl index 5fdd313..a431d39 100644 --- a/src/Apply/full_update.jl +++ b/src/Apply/full_update.jl @@ -44,7 +44,7 @@ function full_update( 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 + Rᵥ₁, Rᵥ₂ = U * replaceinds(sqrtS, v => prime(u)), replaceinds(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)) @@ -69,34 +69,34 @@ function fidelity( p_sind, q_sind = trycommonind(p_cur, gate), trycommonind(q_cur, gate) p_sind_sim, q_sind_sim = sim(p_sind), sim(q_sind) gate_sq = - gate * replaceinds(conj(gate), Index[p_sind, q_sind], Index[p_sind_sim, q_sind_sim]) + gate * replaceinds(conj(gate), p_sind => p_sind_sim, q_sind => q_sind_sim) term1_tns = vcat( [ p_prev, q_prev, - replaceind(prime(conj(p_prev)), prime(p_sind), p_sind_sim), - replaceind(prime(conj(q_prev)), prime(q_sind), q_sind_sim), + replaceinds(prime(conj(p_prev)), prime(p_sind) => p_sind_sim), + replaceinds(prime(conj(q_prev)), prime(q_sind) => q_sind_sim), gate_sq, ], envs, ) sequence = contraction_sequence(term1_tns; alg = "optimal") - term1 = contract(term1_tns; sequence) + term1 = contract_network(term1_tns; sequence) term2_tns = vcat( [ p_cur, q_cur, - replaceind(prime(conj(p_cur)), prime(p_sind), p_sind), - replaceind(prime(conj(q_cur)), prime(q_sind), q_sind), + replaceinds(prime(conj(p_cur)), prime(p_sind) => p_sind), + replaceinds(prime(conj(q_cur)), prime(q_sind) => q_sind), ], envs, ) sequence = contraction_sequence(term2_tns; alg = "optimal") - term2 = contract(term2_tns; sequence) + term2 = contract_network(term2_tns; sequence) term3_tns = vcat([p_prev, q_prev, prime(conj(p_cur)), prime(conj(q_cur)), gate], envs) sequence = contraction_sequence(term3_tns; alg = "optimal") - term3 = contract(term3_tns; sequence) + term3 = contract_network(term3_tns; sequence) f = scalar(term3) / sqrt(scalar(term1) * scalar(term2)) return f * conj(f) @@ -120,7 +120,7 @@ function optimise_p_q( ) b = only(commoninds(p_cur, q_cur)) bnew = settags(b, tags(trycommonind(p, q))) - p_cur, q_cur = replaceind(p_cur, b, bnew), replaceind(q_cur, b, bnew) + p_cur, q_cur = replaceinds(p_cur, b => bnew), replaceinds(q_cur, b => bnew) fstart = print_fidelity_loss ? fidelity(envs, p_cur, q_cur, p, q, o) : 0 @@ -130,18 +130,18 @@ function optimise_p_q( function b(p::ITensor, q::ITensor, o::ITensor, envs::Vector{ITensor}, r::ITensor) ts = vcat(ITensor[p, q, o, conj(prime(r))], envs) sequence = contraction_sequence(ts; alg = "optimal") - return noprime(contract(ts; sequence)) + return noprime(contract_network(ts; sequence)) end function M_p(envs::Vector{ITensor}, p_q_tensor::ITensor, s_ind, apply_tensor::ITensor) ts = vcat( ITensor[ - p_q_tensor, replaceinds(prime(conj(p_q_tensor)), prime.(s_ind), s_ind), apply_tensor, + p_q_tensor, replaceinds(prime(conj(p_q_tensor)), (prime.(s_ind) .=> s_ind)...), apply_tensor, ], envs, ) sequence = contraction_sequence(ts; alg = "optimal") - return noprime(contract(ts; sequence)) + return noprime(contract_network(ts; sequence)) end for i in 1:nfullupdatesweeps b_vec = b(p, q, o, envs, q_cur) diff --git a/src/Apply/simple_update.jl b/src/Apply/simple_update.jl index 9c24c05..9aee2b8 100644 --- a/src/Apply/simple_update.jl +++ b/src/Apply/simple_update.jl @@ -46,8 +46,8 @@ function simple_update( 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) - ψᵥ₁ = contract([ψ⃗[1]; sqrt_envs_v1]) - ψᵥ₂ = contract([ψ⃗[2]; sqrt_envs_v2]) + ψᵥ₁ = contract_network([ψ⃗[1]; sqrt_envs_v1]) + ψᵥ₂ = contract_network([ψ⃗[2]; sqrt_envs_v2]) sᵥ₁ = commoninds(ψ⃗[1], o) sᵥ₂ = commoninds(ψ⃗[2], o) Qᵥ₁, Rᵥ₁ = MAK.qr_compact(ψᵥ₁, setdiff(uniqueinds(ψᵥ₁, ψᵥ₂), sᵥ₁)) @@ -63,7 +63,7 @@ function simple_update( 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 + Rᵥ₁, Rᵥ₂ = U * replaceinds(sqrtS, v => prime(u)), replaceinds(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 `ϵ` @@ -71,8 +71,8 @@ function simple_update( total = abs2(norm(oR)) err = iszero(total) ? zero(real(scalartype(oR))) : max(zero(real(scalartype(oR))), 1 - abs2(norm(S)) / total) - Qᵥ₁ = contract([Qᵥ₁; conj.(inv_sqrt_envs_v1)]) - Qᵥ₂ = contract([Qᵥ₂; conj.(inv_sqrt_envs_v2)]) + Qᵥ₁ = contract_network([Qᵥ₁; conj.(inv_sqrt_envs_v1)]) + Qᵥ₂ = contract_network([Qᵥ₂; conj.(inv_sqrt_envs_v2)]) updated_tensors = [Qᵥ₁ * Rᵥ₁, Qᵥ₂ * Rᵥ₂] if normalize_tensors s_values = normalize(s_values) @@ -81,7 +81,7 @@ function simple_update( if normalize_tensors for ψᵥ in updated_tensors - rmul!(data(ψᵥ), inv(norm(ψᵥ))) + rmul!(ψᵥ, inv(norm(ψᵥ))) end end diff --git a/src/MessagePassing/abstractbeliefpropagationcache.jl b/src/MessagePassing/abstractbeliefpropagationcache.jl index dee7e6a..2fdac86 100644 --- a/src/MessagePassing/abstractbeliefpropagationcache.jl +++ b/src/MessagePassing/abstractbeliefpropagationcache.jl @@ -24,7 +24,7 @@ function vertex_scalar(bp_cache::AbstractBeliefPropagationCache, vertex) state = bp_factors(bp_cache, vertex) contract_list = [state; incoming_ms] sequence = contraction_sequence(contract_list; alg = "optimal") - return scalar(contract(contract_list; sequence)) + return scalar(contract_network(contract_list; sequence)) end function edge_scalar( @@ -177,7 +177,7 @@ function updated_message( sequence = contraction_sequence(contract_list; alg = alg.kwargs.sequence_alg) seq_changed = true end - updated_message = contract(contract_list; sequence) + updated_message = contract_network(contract_list; sequence) if alg.kwargs.normalize message_norm = sum(updated_message) diff --git a/src/MessagePassing/beliefpropagationcache.jl b/src/MessagePassing/beliefpropagationcache.jl index e580c27..665c3f4 100644 --- a/src/MessagePassing/beliefpropagationcache.jl +++ b/src/MessagePassing/beliefpropagationcache.jl @@ -119,7 +119,7 @@ default_bp_update_kwargs(bp_cache::BeliefPropagationCache) = default_bp_update_k function make_hermitian(A::ITensor) A_inds = inds(A) @assert length(A_inds) == 2 - return (A + swapind(conj(A), first(A_inds), last(A_inds))) / 2 + return (A + replaceinds(conj(A), first(A_inds) => last(A_inds), last(A_inds) => first(A_inds))) / 2 end function rescale_messages!(bp_cache::BeliefPropagationCache, edges::Vector{<:AbstractEdge}) @@ -159,20 +159,20 @@ function loop_correlation(bpc::BeliefPropagationCache, loop::Vector{<:NamedEdge} if !isempty(t_inds) t_ind = only(t_inds) t_ind_pos = findfirst(x -> x == t_ind, e_virtualinds) - t = replaceind(t, t_ind, e_virtualinds_sim[t_ind_pos]) + t = replaceinds(t, t_ind => e_virtualinds_sim[t_ind_pos]) end push!(local_tensors, t) end tensors = ITensor[local_tensors; reduce(vcat, [bp_factors(bpc, v) for v in setdiff(vs, [src_vertex])]); incoming_messages] seq = contraction_sequence(tensors; alg = "omeinsum", optimizer = GreedyMethod()) - t = contract(tensors; sequence = seq) + t = contract_network(tensors; sequence = seq) row_name = ITensorBase.uniquename(ITensorBase.IndexName) col_name = ITensorBase.uniquename(ITensorBase.IndexName) t = matricize(t, Tuple(e_virtualinds) => row_name, Tuple(e_virtualinds_sim) => col_name) t = adapt(Vector{ComplexF64})(t) - t = array(t) + t = Array(t) λs = reverse(sort(LinearAlgebra.eigvals(t); by = abs)) err = 1 - abs(λs[1]) / sum(abs.(λs)) return err diff --git a/src/MessagePassing/boundarympscache.jl b/src/MessagePassing/boundarympscache.jl index 32a97e3..6e4e5f6 100644 --- a/src/MessagePassing/boundarympscache.jl +++ b/src/MessagePassing/boundarympscache.jl @@ -195,7 +195,7 @@ function set_interpartition_messages!( # 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 = settags(trivialrange(first(inds(m1)), virt_dim), "m$(i)$(i + 1)") + ind = settags(trivialrange(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. @@ -433,7 +433,7 @@ function generic_apply( L, R = MAK.left_orth(T, keep; trunc = itensor_trunc(; cutoff, maxdim)) b = only(commoninds(L, R)) bnew = settags(b, "link" => "$i") - L, R = replaceind(L, b, bnew), replaceind(R, b, bnew) + L, R = replaceinds(L, b => bnew), replaceinds(R, b => bnew) push!(out, L) carry = R left_link = only(commoninds(L, R)) @@ -447,7 +447,7 @@ function generic_apply( L, R = MAK.right_orth(out[i], [bond]; trunc = itensor_trunc(; cutoff, maxdim)) b = only(commoninds(L, R)) bnew = settags(b, "link" => "$(i - 1)") - L, R = replaceind(L, b, bnew), replaceind(R, b, bnew) + L, R = replaceinds(L, b => bnew), replaceinds(R, b => bnew) out[i] = R out[i - 1] *= L end @@ -653,7 +653,7 @@ function path_contract( m != nothing && push!(contract_list, m) sequence = contraction_sequence(contract_list; alg = "optimal") - m = contract(contract_list; sequence) + m = contract_network(contract_list; sequence) prev_edge = e end @@ -662,13 +662,13 @@ function path_contract( append!(contract_list, incoming_ms) push!(contract_list, m) sequence = contraction_sequence(contract_list; alg = "optimal") - numer = contract(contract_list; sequence) + numer = contract_network(contract_list; sequence) else contract_list = norm_factors(network(cache), vs; op_strings = op_string_f) incoming_ms = incoming_messages(cache, only(vs)) append!(contract_list, incoming_ms) sequence = contraction_sequence(contract_list; alg = "optimal") - numer = contract(contract_list; sequence) + numer = contract_network(contract_list; sequence) end return numer, denom diff --git a/src/MessagePassing/loopcorrection.jl b/src/MessagePassing/loopcorrection.jl index 6e7acdd..121dee5 100644 --- a/src/MessagePassing/loopcorrection.jl +++ b/src/MessagePassing/loopcorrection.jl @@ -36,9 +36,9 @@ function sim_edgeinduced_subgraph(bpc::BeliefPropagationCache, eg) mer = message(bpc, reverse(e)) linds = filter(i -> plev(i) == 0, inds(mer)) linds_sim = sim.(linds) - mer = replaceinds(mer, linds, linds_sim) + mer = replaceinds(mer, (linds .=> linds_sim)...) if network(bpc) isa TensorNetworkState - mer = replaceinds(mer, conj.(prime.(linds)), conj.(prime.(linds_sim))) + mer = replaceinds(mer, (conj.(prime.(linds)) .=> conj.(prime.(linds_sim)))...) end ms = messages(bpc) set!(ms, reverse(e), mer) @@ -49,7 +49,7 @@ function sim_edgeinduced_subgraph(bpc::BeliefPropagationCache, eg) if !isempty(t_inds) t_ind = only(t_inds) t_ind_pos = findfirst(==(t_ind), linds) - t = replaceind(t, t_ind, linds_sim[t_ind_pos]) + t = replaceinds(t, t_ind => linds_sim[t_ind_pos]) setindex_preserve!(bpc, t, src(e)) end push!(updated_es, e) @@ -105,7 +105,7 @@ function weight(bpc::BeliefPropagationCache, eg) 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)) + return scalar(contract_network(ts; sequence = seq)) end #Vectorized version of weight diff --git a/src/TensorNetworkQuantumSimulator.jl b/src/TensorNetworkQuantumSimulator.jl index f3621d1..8d7cf87 100644 --- a/src/TensorNetworkQuantumSimulator.jl +++ b/src/TensorNetworkQuantumSimulator.jl @@ -92,7 +92,7 @@ export TensorNetworkState, AbstractTensorNetwork, partitionfunction, - contract, + contract_network, TreeSA, GreedyMethod, SABipartite, diff --git a/src/TensorNetworks/abstracttensornetwork.jl b/src/TensorNetworks/abstracttensornetwork.jl index a39db4e..d39a640 100644 --- a/src/TensorNetworks/abstracttensornetwork.jl +++ b/src/TensorNetworks/abstracttensornetwork.jl @@ -94,8 +94,8 @@ function map_virtualinds!(f::Function, tn::AbstractTensorNetwork) for e in edges(tn) vinds = commoninds(tn[src(e)], tn[dst(e)]) vinds_sim = f.(vinds) - setindex_preserve!(tn, replaceinds(tn[src(e)], vinds, vinds_sim), src(e)) - setindex_preserve!(tn, replaceinds(tn[dst(e)], vinds, vinds_sim), dst(e)) + setindex_preserve!(tn, replaceinds(tn[src(e)], (vinds .=> vinds_sim)...), src(e)) + setindex_preserve!(tn, replaceinds(tn[dst(e)], (vinds .=> vinds_sim)...), dst(e)) end return tn end diff --git a/src/TensorNetworks/tensornetworkstate.jl b/src/TensorNetworks/tensornetworkstate.jl index dbed9b0..f53bc6c 100644 --- a/src/TensorNetworks/tensornetworkstate.jl +++ b/src/TensorNetworks/tensornetworkstate.jl @@ -59,7 +59,7 @@ end # between ket and bra rather than through an operator or a message. function bra_tensor(t::ITensor, auxinds::Vector{<:Index}) tdag = conj(prime(t)) - return isempty(auxinds) ? tdag : replaceinds(tdag, prime.(auxinds), auxinds) + return isempty(auxinds) ? tdag : replaceinds(tdag, (prime.(auxinds) .=> auxinds)...) end bra_tensor(tns::TensorNetworkState, v) = bra_tensor(tns[v], auxinds(tns, v)) @@ -78,7 +78,7 @@ function norm_factors(tns::TensorNetworkState, verts::Vector; op_strings::Functi if op_strings(v) == "ρ" || isempty(sinds) append!(factors, ITensor[tnv, tnv_dag]) elseif op_strings(v) == "I" - tnv_dag = replaceinds(tnv_dag, prime.(sinds), sinds) + tnv_dag = replaceinds(tnv_dag, (prime.(sinds) .=> sinds)...) append!(factors, ITensor[tnv, tnv_dag]) else op = adapt_like(tnv, Ops.op(op_strings(v), only(sinds))) diff --git a/src/TensorNetworks/tensornetworkstate_constructors.jl b/src/TensorNetworks/tensornetworkstate_constructors.jl index 1d60f3d..ae4ff0a 100644 --- a/src/TensorNetworks/tensornetworkstate_constructors.jl +++ b/src/TensorNetworks/tensornetworkstate_constructors.jl @@ -19,7 +19,7 @@ first half being the "ket" indices and the second half the "bra" indices; index index `n/2 + i`. """ function identity_tensornetworkstate(eltype, g::NamedGraph, s::Dictionary = siteinds("S=1/2", g; inds_per_site = 2)) - links = Dictionary(edges(g), [settags(Index(1), "e$(src(e))_$(dst(e))") for e in edges(g)]) + links = Dictionary(edges(g), [settags(Index(1), "e" => "$(src(e))_$(dst(e))") for e in edges(g)]) links = merge(links, Dictionary(reverse.(edges(g)), [links[e] for e in edges(g)])) ts = Dictionary{vertextype(g), ITensor}() @@ -111,7 +111,7 @@ Returns a `TensorNetwork` (not a `TensorNetworkState`); contract it to obtain ``Z(β)``. """ function ising_partitionfunction(g::NamedGraph, β::Real; Js::Dictionary = Dictionary(edges(g), [1.0 for e in edges(g)])) - links = Dictionary(edges(g), [settags(Index(2), "e$(src(e))_$(dst(e))") for e in edges(g)]) + links = Dictionary(edges(g), [settags(Index(2), "e" => "$(src(e))_$(dst(e))") for e in edges(g)]) links = merge(links, Dictionary(reverse.(edges(g)), [links[e] for e in edges(g)])) # symmetric sqrt of Boltzmann matrix W = exp(β σσ') diff --git a/src/contract.jl b/src/contract.jl index cf72ab1..bc452b9 100644 --- a/src/contract.jl +++ b/src/contract.jl @@ -1,17 +1,17 @@ -function contract(alg::Algorithm"exact", tn::AbstractTensorNetwork; contraction_sequence_kwargs = (; alg = "omeinsum", optimizer = GreedyMethod())) +function contract_network(alg::Algorithm"exact", tn::AbstractTensorNetwork; contraction_sequence_kwargs = (; alg = "omeinsum", optimizer = GreedyMethod())) tn_tensors = [tn[v] for v in vertices(tn)] seq = contraction_sequence(tn_tensors; contraction_sequence_kwargs...) - return scalar(contract(tn_tensors; sequence = seq)) + return scalar(contract_network(tn_tensors; sequence = seq)) end -function contract(alg::Algorithm"bp", tn::TensorNetwork; bp_update_kwargs = default_bp_update_kwargs(tn)) +function contract_network(alg::Algorithm"bp", tn::TensorNetwork; bp_update_kwargs = default_bp_update_kwargs(tn)) return partitionfunction(update(BeliefPropagationCache(tn); bp_update_kwargs...)) end -function contract(alg::Algorithm"boundarymps", tn::TensorNetwork; mps_bond_dimension::Integer, bmps_update_kwargs = default_bmps_update_kwargs(tn)) +function contract_network(alg::Algorithm"boundarymps", tn::TensorNetwork; mps_bond_dimension::Integer, bmps_update_kwargs = default_bmps_update_kwargs(tn)) return partitionfunction(update(BoundaryMPSCache(tn, mps_bond_dimension); bmps_update_kwargs...)) end -function contract(tn::AbstractTensorNetwork; alg = "exact", kwargs...) - return contract(Algorithm(alg), tn; kwargs...) +function contract_network(tn::AbstractTensorNetwork; alg = "exact", kwargs...) + return contract_network(Algorithm(alg), tn; kwargs...) end diff --git a/src/entanglement.jl b/src/entanglement.jl index 4276c86..5afa8f7 100644 --- a/src/entanglement.jl +++ b/src/entanglement.jl @@ -32,7 +32,7 @@ function matricize(a::ITensor, row_inds = filter(i -> plev(i) ==0, inds(a))) col_inds = prime.(row_inds) row_name = ITensorBase.uniquename(ITensorBase.IndexName) col_name = ITensorBase.uniquename(ITensorBase.IndexName) - return array(matricize(a, Tuple(row_inds) => row_name, Tuple(col_inds) => col_name)) + return Array(matricize(a, Tuple(row_inds) => row_name, Tuple(col_inds) => col_name)) end """ @@ -87,8 +87,8 @@ function renyi_entropy( ) 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 - ρ = replaceind(ρ, edge_ind_pp, edge_ind_p) + ρ = (m1 * replaceinds(root_m2, edge_ind_p => edge_ind_pp)) * root_m2 + ρ = replaceinds(ρ, edge_ind_pp => edge_ind_p) return renyi_entropy(ρ; α) end diff --git a/src/expect.jl b/src/expect.jl index 02fe6f1..77742c5 100644 --- a/src/expect.jl +++ b/src/expect.jl @@ -16,7 +16,7 @@ function expect( op_string_f = op_string_function(op_strings, vs) ψOψ_tensors = norm_factors(ψ, collect(vertices(ψ)); op_strings = op_string_f) numer_seq = contraction_sequence(ψOψ_tensors; contraction_sequence_kwargs...) - numer = scalar(contract(ψOψ_tensors; sequence = numer_seq)) + numer = scalar(contract_network(ψOψ_tensors; sequence = numer_seq)) push!(out, coeff * (numer / denom)) end return out @@ -71,7 +71,7 @@ function expect( tensors = norm_factors(network(cache), steiner_vs; op_strings = op_string_f) append!(tensors, incoming_ms) seq = contraction_sequence(tensors; alg = "optimal", prune_tensors = true) - return scalar(contract(tensors; sequence = seq)) + return scalar(contract_network(tensors; sequence = seq)) end denom = contract_region(v -> "I") diff --git a/src/inner.jl b/src/inner.jl index c1d581d..ba7404c 100644 --- a/src/inner.jl +++ b/src/inner.jl @@ -52,7 +52,7 @@ function inner( ) blf_tensors = bp_factors(blf, collect(vertices(ket(blf)))) seq = contraction_sequence(blf_tensors; contraction_sequence_kwargs...) - return scalar(contract(blf_tensors; sequence = seq)) + return scalar(contract_network(blf_tensors; sequence = seq)) end function inner(alg::Algorithm, cache::AbstractBeliefPropagationCache; max_configuration_size = nothing) diff --git a/src/itensors.jl b/src/itensors.jl index 1b7404e..6db69ac 100644 --- a/src/itensors.jl +++ b/src/itensors.jl @@ -4,11 +4,11 @@ import Base: truncate import ITensorBase: scalartype, uniqueinds import MatrixAlgebraKit as MAK -import TensorAlgebra: matricize +import TensorAlgebra: datatype, matricize using Adapt: Adapt using ITensorBase: ITensorBase, AbstractITensor, ITensor, Index, NamedUnitRange, commoninds, dimnames, hascommoninds, id, inds, name, nameddims, noncommoninds, noprime, plev, prime, - replacedimnames, sim, tags, trycommonind, trynoncommonind, unioninds, unnamed + replaceinds, sim, tags, trycommonind, trynoncommonind, unioninds, unnamed using LinearAlgebra: LinearAlgebra using TensorAlgebra: TensorAlgebra, project, scalar, tryproject @@ -31,29 +31,13 @@ function onehot(eltype::Type, (i, p)::Pair{<:Index}) end onehot(p::Pair{<:Index}) = onehot(Float64, p) -function inner end - -function contract end -function contract(tensors::AbstractVector; sequence = nothing) +function contract_network end +function contract_network(tensors::AbstractVector; sequence = nothing) return isnothing(sequence) ? reduce(*, tensors) : _contract_sequence(tensors, sequence) end _contract_sequence(tensors, s::Integer) = tensors[s] _contract_sequence(tensors, s) = reduce(*, (_contract_sequence(tensors, x) for x in s)) -replaceind(t, p::Pair) = replaceinds(t, p) -replaceind(t, from::Index, to::Index) = replaceinds(t, from => to) - -const _IndexColl = Union{Tuple{Vararg{Index}}, AbstractVector{<:Index}} -function replaceinds(t, pairs::Pair...) - return replacedimnames(t, map(p -> name(first(p)) => name(last(p)), pairs)...) -end -function replaceinds(t, from::_IndexColl, to::_IndexColl) - return replaceinds(t, map(=>, from, to)...) -end -function replaceinds(t::AbstractITensor, p::Pair{<:_IndexColl, <:_IndexColl}) - return replaceinds(t, first(p), last(p)) -end - diaglength(a::AbstractArray) = minimum(size(a)) function diagstride(a::AbstractArray) s = 1 @@ -106,10 +90,6 @@ function itensor_trunc(; maxdim = nothing, cutoff = nothing) return trunc end -datatype(T::AbstractITensor) = typeof(unnamed(T)) -array(T::AbstractITensor) = convert(Array, unnamed(T)) -data(T::AbstractITensor) = unnamed(T) - struct ScalarTypeAdaptor{T} end ScalarTypeAdaptor(T::Type) = ScalarTypeAdaptor{T}() adapt_scalartype(T::Type) = Adapt.adapt(ScalarTypeAdaptor(T)) @@ -119,11 +99,8 @@ function Adapt.adapt_structure(::ScalarTypeAdaptor{elt}, T::AbstractITensor) whe return nameddims(convert(AbstractArray{elt}, unnamed(T)), ITensorBase.dimnames(T)) end -swapind(T::AbstractITensor, i::Index, j::Index) = replaceinds(T, i => j, j => i) - hasqns(i::Index) = conj(unnamed(i)) != unnamed(i) hasqns(t::AbstractITensor) = any(hasqns, inds(t)) -hasqns(::Any) = false function directsum(out_inds, pairs::Pair...) t1, s1 = first(pairs[1]), last(pairs[1]) @@ -163,13 +140,6 @@ macro Algorithm_str(s) return :(Algorithm{$(Expr(:quote, Symbol(s)))}) end -function settags(i::Index, tagstr::AbstractString) - for t in split(tagstr, ",") - s = String(strip(t)) - isempty(s) || (i = ITensorBase.settag(i, s, "")) - end - return i -end settags(i::Index, p::Pair) = ITensorBase.settag(i, first(p), last(p)) function settags(i::Index, d::AbstractDict) for (k, v) in d diff --git a/src/norm_sqr.jl b/src/norm_sqr.jl index 6ce9824..fe47079 100644 --- a/src/norm_sqr.jl +++ b/src/norm_sqr.jl @@ -66,7 +66,7 @@ function norm_sqr( ) ψIψ_tensors = norm_factors(ψ, collect(vertices(ψ))) denom_seq = contraction_sequence(ψIψ_tensors; contraction_sequence_kwargs...) - return scalar(contract(ψIψ_tensors; sequence = denom_seq)) + return scalar(contract_network(ψIψ_tensors; sequence = denom_seq)) end function norm_sqr(alg::Algorithm, cache::AbstractBeliefPropagationCache; max_configuration_size = nothing) diff --git a/src/rdm.jl b/src/rdm.jl index 2ecb13a..c3c01a4 100644 --- a/src/rdm.jl +++ b/src/rdm.jl @@ -35,7 +35,7 @@ function reduced_density_matrix( op_string_f = v -> v ∈ verts ? "ρ" : "I" ρ_tensors = norm_factors(ψ, collect(vertices(ψ)); op_strings = op_string_f) seq = contraction_sequence(ρ_tensors; contraction_sequence_kwargs...) - ρ = contract(ρ_tensors; sequence = seq) + ρ = contract_network(ρ_tensors; sequence = seq) if normalize ρ = normalize_rdm(ρ) end @@ -58,7 +58,7 @@ function reduced_density_matrix( ρ_tensors = norm_factors(network(cache), steiner_vs; op_strings = op_string_f) append!(ρ_tensors, incoming_ms) seq = contraction_sequence(ρ_tensors; alg = "optimal", prune_tensors = true) - ρ = contract(ρ_tensors; sequence = seq) + ρ = contract_network(ρ_tensors; sequence = seq) if normalize ρ = normalize_rdm(ρ) diff --git a/src/sampling.jl b/src/sampling.jl index 4370615..aa61a5f 100644 --- a/src/sampling.jl +++ b/src/sampling.jl @@ -23,11 +23,11 @@ function sample( ψ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") - ρ = contract(tensors; sequence = seq) + ρ = contract_network(tensors; sequence = seq) ρ_tr = tr(ρ, operator_inds(ρ)...) ρ *= inv(ρ_tr) - ρ_diag = collect(real.(diag(array(ρ)))) + ρ_diag = collect(real.(diag(Array(ρ)))) config = StatsBase.sample(1:length(ρ_diag), Weights(ρ_diag)) # config is 1,2,...,d, but we want 0,1...,d-1 for the sample itself set!(bit_string, v, config - 1) @@ -237,11 +237,11 @@ function sample_partition!( ψvdag = bra_tensor(network(norm_bmps_cache), v) ts = [incoming_ms; [ψv, ψvdag]] seq = contraction_sequence(ts; alg = "optimal") - ρ = contract(ts; sequence = seq) + ρ = contract_network(ts; sequence = seq) ρ_tr = tr(ρ, operator_inds(ρ)...) push!(traces, ρ_tr) ρ *= inv(ρ_tr) - ρ_diag = collect(real.(diag(array(ρ)))) + ρ_diag = collect(real.(diag(Array(ρ)))) config = StatsBase.sample(1:length(ρ_diag), Weights(ρ_diag)) # config is 1,2,...,d, but we want 0,1...,d-1 for the sample itself set!(bit_string, v, config - 1) diff --git a/src/siteinds.jl b/src/siteinds.jl index d08758e..6d614e4 100644 --- a/src/siteinds.jl +++ b/src/siteinds.jl @@ -18,7 +18,7 @@ end function site_tag(sitetype::String) sitetype = replace(lowercase(sitetype), " " => "") - sitetype ∈ ["s=1/2", "qubit", "spin1/2", "spinhalf"] && return "S=1/2" - sitetype ∈ ["qutrit", "s=1", "spin1"] && return "S=1" + sitetype ∈ ["s=1/2", "qubit", "spin1/2", "spinhalf"] && return "SiteType" => "S=1/2" + sitetype ∈ ["qutrit", "s=1", "spin1"] && return "SiteType" => "S=1" error("Don't know how to interpret that site type. Supported: S=1/2, S=1.") end diff --git a/src/symmetric_gauge.jl b/src/symmetric_gauge.jl index 225aa01..1f770b3 100644 --- a/src/symmetric_gauge.jl +++ b/src/symmetric_gauge.jl @@ -27,19 +27,19 @@ function symmetric_gauge!(bp_cache::BeliefPropagationCache; regularization = 10 # 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 + ψ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) + Ce = rootX * replaceinds(rootY, (edge_ind_p .=> edge_ind_sim)...) U, S, V = MAK.svd_compact(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 + ψ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. @@ -49,15 +49,15 @@ function symmetric_gauge!(bp_cache::BeliefPropagationCache; regularization = 10 ψvdst = ψvdst * sqrtS new_edge_ind = Index[settags(only(v), tags(first(edge_ind)))] - ψvsrc = replaceinds(ψvsrc, v, new_edge_ind) - ψvdst = replaceinds(ψvdst, u, new_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) # 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))])) + setmessage!(bp_cache, e, replaceinds(S, only(u) => prime(only(new_edge_ind)), only(v) => only(new_edge_ind))) + setmessage!(bp_cache, reverse(e), replaceinds(S, only(u) => only(new_edge_ind), only(v) => prime(only(new_edge_ind)))) end return bp_cache diff --git a/test/test_add.jl b/test/test_add.jl index 9e712dc..047abb3 100644 --- a/test/test_add.jl +++ b/test/test_add.jl @@ -5,7 +5,7 @@ const TNQS = TensorNetworkQuantumSimulator using Test: @testset, @test # Fully contract a tensor-network state to a dense tensor over its site indices. -contract_dense(ψ) = TNQS.contract([ψ[v] for v in vertices(ψ)]) +contract_dense(ψ) = TNQS.contract_network([ψ[v] for v in vertices(ψ)]) @testset "Tensor-network addition (direct sum)" begin Random.seed!(123) diff --git a/test/test_beliefpropagation.jl b/test/test_beliefpropagation.jl index 160aba8..3f2bdec 100644 --- a/test/test_beliefpropagation.jl +++ b/test/test_beliefpropagation.jl @@ -26,8 +26,8 @@ const TNQS = TensorNetworkQuantumSimulator @test !isempty(messages(ψ_BPC)) @test length(keys(messages(ψ_BPC))) == 2 * length(edges(g)) z_bp = partitionfunction(ψ_BPC) - @test z_bp ≈ contract(ψ; alg = "exact") - @test z_bp ≈ contract(ψ; alg = "bp") + @test z_bp ≈ contract_network(ψ; alg = "exact") + @test z_bp ≈ contract_network(ψ; alg = "bp") end #BP Cache diff --git a/test/test_boundarymps.jl b/test/test_boundarymps.jl index 652aed0..5d4bbb4 100644 --- a/test/test_boundarymps.jl +++ b/test/test_boundarymps.jl @@ -25,8 +25,8 @@ const TNQS = TensorNetworkQuantumSimulator ψ_BMPS = update(ψ_BMPS) @test isempty(TNQS.contraction_sequences(ψ_BMPS)) z_bmps = partitionfunction(ψ_BMPS) - @test z_bmps ≈ contract(ψ; alg = "exact") - @test z_bmps ≈ contract(ψ; alg = "boundarymps", mps_bond_dimension = 4) + @test z_bmps ≈ contract_network(ψ; alg = "exact") + @test z_bmps ≈ contract_network(ψ; alg = "boundarymps", mps_bond_dimension = 4) end #BMPS Cache diff --git a/test/test_constructors.jl b/test/test_constructors.jl index c17da93..0c09d36 100644 --- a/test/test_constructors.jl +++ b/test/test_constructors.jl @@ -36,7 +36,7 @@ using Test: @testset, @test, @test_throws ψdag = map_virtualinds(prime, map_tensors(conj, ψ)) @test ψdag isa TensorNetwork - @test TNQS.contract(ψdag; alg = "exact") ≈ conj(TNQS.contract(ψ; alg = "exact")) + @test TNQS.contract_network(ψdag; alg = "exact") ≈ conj(TNQS.contract_network(ψ; alg = "exact")) v = first(vertices(g)) rem_vertex!(ψ, v) diff --git a/test/test_contraction_sequences.jl b/test/test_contraction_sequences.jl index cfa6de7..a0229d9 100644 --- a/test/test_contraction_sequences.jl +++ b/test/test_contraction_sequences.jl @@ -50,10 +50,10 @@ collect_leaves!(acc, x) = (for y in x; collect_leaves!(acc, y); end; acc) # --- the sequence the backend returns is a *correct* contraction: executing it gives the # same scalar as the independent `optimal` backend. - ref = scalar(TNQS.contract(tensors; sequence = TNQS.contraction_sequence(tensors; alg = "optimal"))) + ref = scalar(TNQS.contract_network(tensors; sequence = TNQS.contraction_sequence(tensors; alg = "optimal"))) for optimizer in (GreedyMethod(), TreeSA()) seq = TNQS.contraction_sequence(tensors; alg = "omeinsum", optimizer) - @test scalar(TNQS.contract(tensors; sequence = seq)) ≈ ref + @test scalar(TNQS.contract_network(tensors; sequence = seq)) ≈ ref end # --- open network: result is a tensor with dangling indices (iy non-empty). @@ -64,7 +64,7 @@ collect_leaves!(acc, x) = (for y in x; collect_leaves!(acc, y); end; acc) open_tensors = [X, Y, Z] # open indices: p, r, t seq_open = TNQS.contraction_sequence(open_tensors; alg = "omeinsum", optimizer = GreedyMethod()) @test sort(collect_leaves!(Int[], seq_open)) == [1, 2, 3] - @test TNQS.contract(open_tensors; sequence = seq_open) ≈ - TNQS.contract(open_tensors; sequence = TNQS.contraction_sequence(open_tensors; alg = "optimal")) + @test TNQS.contract_network(open_tensors; sequence = seq_open) ≈ + TNQS.contract_network(open_tensors; sequence = TNQS.contraction_sequence(open_tensors; alg = "optimal")) end end From 61aa2867996ca94bd4009075b27918c9dfa812d4 Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Thu, 9 Jul 2026 15:20:20 -0400 Subject: [PATCH 18/21] Duck-type settags to accept any iterable of tag pairs Relaxes the `AbstractDict` method to `settags(i, tags)` so it accepts any iterable of `key => value` pairs (a `Vector` or `Tuple` of `Pair`s, not only a `Dict`). The single-`Pair` method stays more specific, so single-tag calls are unaffected. Co-Authored-By: Claude Opus 4.8 (1M context) --- src/itensors.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/itensors.jl b/src/itensors.jl index 6db69ac..d37da78 100644 --- a/src/itensors.jl +++ b/src/itensors.jl @@ -141,8 +141,8 @@ macro Algorithm_str(s) end settags(i::Index, p::Pair) = ITensorBase.settag(i, first(p), last(p)) -function settags(i::Index, d::AbstractDict) - for (k, v) in d +function settags(i::Index, tags) + for (k, v) in tags i = ITensorBase.settag(i, k, v) end return i From 08269c3b04e4868ca96ba28c9582d5042f99515a Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Thu, 9 Jul 2026 15:32:45 -0400 Subject: [PATCH 19/21] Report the exact truncation error from svd_trunc MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `simple_update` and `full_update` now take the truncation error from `svd_trunc`'s `ϵ` (the 2-norm of the discarded singular values, exposed in https://github.com/ITensor/ITensorBase.jl/pull/213) as `(ϵ / ‖M‖)^2`. This is the same relative squared error the code reported before, but computed directly instead of by the `1 - ‖S‖²/‖M‖²` subtraction, which lost all precision when little was discarded. Co-Authored-By: Claude Opus 4.8 (1M context) --- src/Apply/full_update.jl | 12 ++++++------ src/Apply/simple_update.jl | 13 ++++++------- 2 files changed, 12 insertions(+), 13 deletions(-) diff --git a/src/Apply/full_update.jl b/src/Apply/full_update.jl index a431d39..673094b 100644 --- a/src/Apply/full_update.jl +++ b/src/Apply/full_update.jl @@ -40,16 +40,16 @@ function full_update( M = Rᵥ₁ * Rᵥ₂ codomain = inds(Rᵥ₁) # Balanced SVD: split the singular values symmetrically (√S into each factor). - U, S, V = MAK.svd_trunc(M, codomain; trunc = itensor_trunc(; apply_kwargs...)) + U, S, V, ϵ = MAK.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 * replaceinds(sqrtS, v => prime(u)), replaceinds(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) + # Relative squared truncation error, from MatrixAlgebraKit's exact discarded-weight `ϵ` + # (the 2-norm of the discarded singular values) rather than the cancellation-prone + # `1 - ‖S‖²/‖M‖²` norm subtraction. + total = norm(M) + truncation_error = iszero(total) ? zero(real(scalartype(M))) : (ϵ / total)^2 callback(; singular_values = S, truncation_error) end ψᵥ₁ = Qᵥ₁ * Rᵥ₁ diff --git a/src/Apply/simple_update.jl b/src/Apply/simple_update.jl index 9aee2b8..05f15cd 100644 --- a/src/Apply/simple_update.jl +++ b/src/Apply/simple_update.jl @@ -59,18 +59,17 @@ function simple_update( # 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 = MAK.svd_trunc(oR, union(rᵥ₁, sᵥ₁); trunc = itensor_trunc(; apply_kwargs...)) + U, S, V, ϵ = MAK.svd_trunc(oR, union(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 * replaceinds(sqrtS, v => prime(u)), replaceinds(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) + # Relative squared truncation error, from MatrixAlgebraKit's exact discarded-weight `ϵ` + # (the 2-norm of the discarded singular values) rather than the cancellation-prone + # `1 - ‖S‖²/‖oR‖²` norm subtraction. + total = norm(oR) + err = iszero(total) ? zero(real(scalartype(oR))) : (ϵ / total)^2 Qᵥ₁ = contract_network([Qᵥ₁; conj.(inv_sqrt_envs_v1)]) Qᵥ₂ = contract_network([Qᵥ₂; conj.(inv_sqrt_envs_v2)]) updated_tensors = [Qᵥ₁ * Rᵥ₁, Qᵥ₂ * Rᵥ₂] From b3608c7977ec9344d65e04671f77ad8ec887cf95 Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Thu, 9 Jul 2026 17:51:36 -0400 Subject: [PATCH 20/21] Remove matricize(a::ITensor) shim, use upstream auto-naming matricize Deletes the TNQS `matricize(a::ITensor)` helper and calls the new upstream `matricize(a, codomain, domain)` from ITensorBase, which mints the fused output names itself. `renyi_entropy` inlines the call, and the same `uniquename` pattern in `loop_correlation` switches to the upstream form too, so no fresh-name minting leaks into TNQS. `matricize` moves from `import` to `using` since TNQS no longer extends it. Co-Authored-By: Claude Opus 4.8 (1M context) --- src/MessagePassing/beliefpropagationcache.jl | 4 +--- src/entanglement.jl | 9 +-------- src/itensors.jl | 4 ++-- 3 files changed, 4 insertions(+), 13 deletions(-) diff --git a/src/MessagePassing/beliefpropagationcache.jl b/src/MessagePassing/beliefpropagationcache.jl index 665c3f4..ab72bb6 100644 --- a/src/MessagePassing/beliefpropagationcache.jl +++ b/src/MessagePassing/beliefpropagationcache.jl @@ -168,9 +168,7 @@ function loop_correlation(bpc::BeliefPropagationCache, loop::Vector{<:NamedEdge} seq = contraction_sequence(tensors; alg = "omeinsum", optimizer = GreedyMethod()) t = contract_network(tensors; sequence = seq) - row_name = ITensorBase.uniquename(ITensorBase.IndexName) - col_name = ITensorBase.uniquename(ITensorBase.IndexName) - t = matricize(t, Tuple(e_virtualinds) => row_name, Tuple(e_virtualinds_sim) => col_name) + t = matricize(t, e_virtualinds, e_virtualinds_sim) t = adapt(Vector{ComplexF64})(t) t = Array(t) λs = reverse(sort(LinearAlgebra.eigvals(t); by = abs)) diff --git a/src/entanglement.jl b/src/entanglement.jl index 5afa8f7..d841177 100644 --- a/src/entanglement.jl +++ b/src/entanglement.jl @@ -28,13 +28,6 @@ function renyi_entropy(ρ::AbstractMatrix, α::Real; normalize = true) return log(sum(λs .^ α)) / (1 - α) end -function matricize(a::ITensor, row_inds = filter(i -> plev(i) ==0, inds(a))) - col_inds = prime.(row_inds) - row_name = ITensorBase.uniquename(ITensorBase.IndexName) - col_name = ITensorBase.uniquename(ITensorBase.IndexName) - return Array(matricize(a, Tuple(row_inds) => row_name, Tuple(col_inds) => col_name)) -end - """ renyi_entropy(a::ITensor, row_inds = ...; normalize = true, α = 1) @@ -52,7 +45,7 @@ and primed indices are column indices. - `α`: Rényi index (default `1`, i.e. von Neumann entropy). """ function renyi_entropy(a::ITensor, row_inds = filter(i -> plev(i) ==0, inds(a)); normalize = true, α = 1) - return renyi_entropy(matricize(a, row_inds), α) + return renyi_entropy(Array(matricize(a, row_inds, prime.(row_inds))), α) end """ diff --git a/src/itensors.jl b/src/itensors.jl index d37da78..3352b04 100644 --- a/src/itensors.jl +++ b/src/itensors.jl @@ -4,13 +4,13 @@ import Base: truncate import ITensorBase: scalartype, uniqueinds import MatrixAlgebraKit as MAK -import TensorAlgebra: datatype, matricize +import TensorAlgebra: datatype using Adapt: Adapt using ITensorBase: ITensorBase, AbstractITensor, ITensor, Index, NamedUnitRange, commoninds, dimnames, hascommoninds, id, inds, name, nameddims, noncommoninds, noprime, plev, prime, replaceinds, sim, tags, trycommonind, trynoncommonind, unioninds, unnamed using LinearAlgebra: LinearAlgebra -using TensorAlgebra: TensorAlgebra, project, scalar, tryproject +using TensorAlgebra: TensorAlgebra, matricize, project, scalar, tryproject function project_aux(v::AbstractVector{<:Number}, i::Index) length(v) == length(i) || From af6dd7f2beedd803bef11dd7e2f6677316b0f307 Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Thu, 9 Jul 2026 20:20:24 -0400 Subject: [PATCH 21/21] Drop the ITensorBase/TensorAlgebra [sources] pins ITensorBase 0.11.0 and TensorAlgebra 0.17.1 are registered, so drop the branch `[sources]` pins and resolve both from the registry. The existing `[compat]` entries already admit the registered versions. Co-Authored-By: Claude Opus 4.8 (1M context) --- Project.toml | 8 -------- 1 file changed, 8 deletions(-) diff --git a/Project.toml b/Project.toml index 38023f9..294a5a7 100644 --- a/Project.toml +++ b/Project.toml @@ -25,14 +25,6 @@ TensorAlgebra = "68bd88dc-f39d-4e12-b2ca-f046b68fcc6a" TensorOperations = "6aa20fa7-93e2-5fca-9bc0-fbd0db3c71a2" TypeParameterAccessors = "7e5a90cf-f82e-492e-a09b-e3e26432c138" -[sources.ITensorBase] -url = "https://github.com/ITensor/ITensorBase.jl" -rev = "mf/nextgen-followups-round1" - -[sources.TensorAlgebra] -url = "https://github.com/ITensor/TensorAlgebra.jl" -rev = "mf/nextgen-followups-round1" - [compat] Adapt = "4.3.0" Dictionaries = "0.4"