diff --git a/Project.toml b/Project.toml index 9788744..294a5a7 100644 --- a/Project.toml +++ b/Project.toml @@ -29,7 +29,7 @@ TypeParameterAccessors = "7e5a90cf-f82e-492e-a09b-e3e26432c138" 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/docs/src/states.md b/docs/src/states.md index dfd647a..0aabc1e 100644 --- a/docs/src/states.md +++ b/docs/src/states.md @@ -8,22 +8,22 @@ 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])) # 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/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..5407378 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 @@ -111,14 +111,14 @@ 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 # handling for this message construction is future work. - s_values = replaceind(s_values, v, prime(u)) - setmessage!(ψ_bpc, e, dag(s_values)) + 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 75b81fc..673094b 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...) @@ -19,14 +18,14 @@ function full_update( apply_kwargs..., ) - Qᵥ₁, Rᵥ₁ = factorize( - ψ[v⃗[1]], uniqueinds(uniqueinds(ψ[v⃗[1]], ψ[v⃗[2]]), uniqueinds(ψ, v⃗[1])) + Qᵥ₁, Rᵥ₁ = MAK.left_orth( + ψ[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])) + Qᵥ₂, Rᵥ₂ = MAK.left_orth( + ψ[v⃗[2]], setdiff(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ᵥ₂, @@ -41,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 = 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 * replaceind(sqrtS, v, prime(u)), replaceind(sqrtS, u, prime(u)) * V - # Best-effort truncation error from norms; suffers catastrophic cancellation when little is - # discarded. TODO: expose MatrixAlgebraKit's `ϵ` from `ITensorBase.svd_trunc` and use it here. - total = abs2(norm(M)) - truncation_error = iszero(total) ? zero(real(scalartype(M))) : - max(zero(real(scalartype(M))), 1 - abs2(norm(S)) / total) + Rᵥ₁, Rᵥ₂ = U * replaceinds(sqrtS, v => prime(u)), replaceinds(sqrtS, u => prime(u)) * V + # 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ᵥ₁ @@ -67,44 +66,44 @@ 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(dag(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(dag(p_prev)), prime(p_sind), p_sind_sim), - replaceind(prime(dag(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 = ITensors.contract(term1_tns; sequence) + term1 = contract_network(term1_tns; sequence) term2_tns = vcat( [ p_cur, q_cur, - replaceind(prime(dag(p_cur)), prime(p_sind), p_sind), - replaceind(prime(dag(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 = ITensors.contract(term2_tns; sequence) - term3_tns = vcat([p_prev, q_prev, prime(dag(p_cur)), prime(dag(q_cur)), gate], envs) + 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 = ITensors.contract(term3_tns; sequence) + term3 = contract_network(term3_tns; sequence) f = scalar(term3) / sqrt(scalar(term1) * scalar(term2)) return f * conj(f) 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, @@ -115,30 +114,34 @@ function optimise_p_q( envisposdef = true, apply_kwargs..., ) - p_cur, q_cur = factorize( - apply(o, p * q), inds(p); tags = tags(commonind(p, q)), apply_kwargs... + pq = apply(o, p * q) + 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 = replaceinds(p_cur, b => bnew), replaceinds(q_cur, b => bnew) 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, dag(prime(r))], envs) + ts = vcat(ITensor[p, q, o, conj(prime(r))], envs) sequence = contraction_sequence(ts; alg = "optimal") - return noprime(ITensors.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(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, ) sequence = contraction_sequence(ts; alg = "optimal") - return noprime(ITensors.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/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..05f15cd 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 @@ -46,33 +46,32 @@ 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ᵥ₁ = qr(ψᵥ₁, uniqueinds(uniqueinds(ψᵥ₁, ψᵥ₂), sᵥ₁)) - Qᵥ₂, Rᵥ₂ = qr(ψᵥ₂, uniqueinds(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 = 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 # `(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, ϵ = 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 * 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 `ϵ` - # from `ITensorBase.svd_trunc` and use it here instead. - total = abs2(norm(oR)) - err = iszero(total) ? zero(real(scalartype(oR))) : - max(zero(real(scalartype(oR))), 1 - abs2(norm(S)) / total) - Qᵥ₁ = contract([Qᵥ₁; dag.(inv_sqrt_envs_v1)]) - Qᵥ₂ = contract([Qᵥ₂; dag.(inv_sqrt_envs_v2)]) + # 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ᵥ₂] if normalize_tensors s_values = normalize(s_values) @@ -81,7 +80,7 @@ function simple_update( if normalize_tensors for ψᵥ in updated_tensors - rmul!(ITensors.data(ψᵥ), inv(norm(ψᵥ))) + rmul!(ψᵥ, inv(norm(ψᵥ))) end end diff --git a/src/Forms/abstractform.jl b/src/Forms/abstractform.jl index 70b4c2d..4a667af 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 @@ -23,8 +23,8 @@ end function default_message(form::AbstractForm, edge::AbstractEdge) cod = virtualinds(ket(form), edge) - dom = dag.(bra_virtualinds(form, edge)) - return one(similar_map(ket(form)[src(edge)], cod, dom), cod, dom) + dom = conj.(bra_virtualinds(form, edge)) + 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 faadbdb..dbe7f2e 100644 --- a/src/Forms/bilinearform.jl +++ b/src/Forms/bilinearform.jl @@ -19,8 +19,8 @@ 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])) - one(similar_map(ket[v], codomain, domain), codomain, domain) + let codomain = conj.(sinds[v]), domain = conj.(prime.(sinds[v])) + one(similar(ket[v], codomain, domain), codomain, domain) end for v in verts ] operator = TensorNetworkState(Dictionary(verts, operator_tensors)) diff --git a/src/Forms/quadraticform.jl b/src/Forms/quadraticform.jl index 5c6ab3a..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)) @@ -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/ITensorsITensorBaseCompat/itensor.jl b/src/ITensorsITensorBaseCompat/itensor.jl deleted file mode 100644 index 8c8507f..0000000 --- a/src/ITensorsITensorBaseCompat/itensor.jl +++ /dev/null @@ -1,586 +0,0 @@ -# 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. - -import MatrixAlgebraKit as MAK -using Adapt: Adapt -using ITensorBase: ITensorBase, AbstractITensor, ITensor, Index, NamedUnitRange, name, - nameddims, plev, tags, unnamed -using LinearAlgebra: LinearAlgebra -using TensorAlgebra.MatrixAlgebra: MatrixAlgebra -using TensorAlgebra: TensorAlgebra - -# Legacy `inds(t; plev, tags)` took index-filtering keywords; `ITensorBase.inds` takes none. Compat -# `inds` forwards to it and applies the legacy filters (a compat-owned function, not a pirated method). -function inds(t::AbstractITensor; plev = nothing, tags = nothing) - is = ITensorBase.inds(t) - isnothing(plev) || (is = filter(i -> ITensorBase.plev(i) == plev, is)) - isnothing(tags) || (is = filter(i -> hastags(i, tags), is)) - 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) -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 -# `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) -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)) -smallsymdiff(a, b; by = identity) = vcat(smallsetdiff(a, b; by), smallsetdiff(b, a; by)) -smallisdisjoint(a, b; by = identity) = (kb = Iterators.map(by, b); !any(x -> by(x) ∈ kb, a)) -smallissubset(a, b; by = identity) = (kb = Iterators.map(by, b); all(x -> by(x) ∈ kb, a)) -smallissetequal(a, b; by = identity) = smallissubset(a, b; by) && smallissubset(b, a; by) - -# -# Name-based index-set algebra: the small-collection ops keyed by `IndexName` (`by = name`) rather -# than full `Index` equality. On a graded axis a shared bond appears as an index on one tensor and -# its dual (`conj`) on the other — same name, opposite arrow — so full-`Index` `==` misses it while -# the names match; on the dense backend the two coincide. The `Index`-returning ops give back full -# `Index`es (callers need the space to feed `Index`/`qr`/`replaceinds`). -# -nameisequal(i, j) = name(i) == name(j) -nameintersect(a, b) = smallintersect(a, b; by = name) -namesetdiff(a, b) = smallsetdiff(a, b; by = name) -nameunion(a, b) = smallunion(a, b; by = name) -namesymdiff(a, b) = smallsymdiff(a, b; by = name) -nameisdisjoint(a, b) = smallisdisjoint(a, b; by = name) -nameissubset(a, b) = smallissubset(a, b; by = name) -nameissetequal(a, b) = smallissetequal(a, b; by = name) - -# -# Named-tensor index-set algebra (the `commoninds`/etc. surface). `_compat_inds` coerces a tensor to -# its indices and passes an index collection through unchanged, so these accept either; the `name*` -# ops take it from there (always returning `Vector`s, so a result can be fed back in, as -# `reduce(noncommoninds, tensors)` does). -# -commoninds(a, b) = nameintersect(_compat_inds(a), _compat_inds(b)) -uniqueinds(a, b) = namesetdiff(_compat_inds(a), _compat_inds(b)) -unioninds(a, b) = nameunion(_compat_inds(a), _compat_inds(b)) -hascommoninds(a, b) = !nameisdisjoint(_compat_inds(a), _compat_inds(b)) - -# Singular forms: the one common / one unique index. `noncommonind` is used in TNQS -# as "the index of `a` not shared with `b`" (e.g. the non-contracted leg of an -# 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)) - -# 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) - -# -# 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. -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(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 compat 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 ITensorBase.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 - -# -# 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 -# `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...) - -# 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[] - -# Dense Kronecker delta tensor (legacy `delta`), vendored from ITensorNetworksNext's -# `ITensorNetworkGenerators/delta_network.jl` (TNQS doesn't depend on it for this migration). A -# graded/sector-aware `delta` is a stack gap (tracked separately). -diaglength(a::AbstractArray) = minimum(size(a)) -function diagstride(a::AbstractArray) - s = 1 - p = 1 - for i in 1:(ndims(a) - 1) - p *= size(a, i) - s += p - end - return s -end -function diagindices(a::AbstractArray) - maxdiag = LinearIndices(a)[CartesianIndex(ntuple(Returns(diaglength(a)), ndims(a)))] - return 1:diagstride(a):maxdiag -end -diagview(a::AbstractArray) = @view a[diagindices(a)] - -function diagonaltensor(diag::AbstractVector, ax::Tuple{Vararg{AbstractUnitRange}}) - a = similar(diag, ax) - fill!(a, zero(eltype(a))) - diagview(a) .= diag - return a -end -function diagonaltensor( - diag::AbstractVector, - is::Tuple{NamedUnitRange, Vararg{NamedUnitRange}} - ) - return nameddims(diagonaltensor(diag, unnamed.(is)), name.(is)) -end - -# Allocate an identity-map-shaped tensor over the `codomain`/`domain` indices, following -# `prototype`'s backend/eltype (domain axes dualized in storage). Unlike `Base.one(a, codomain, -# domain)` — which needs `a` to already carry the map's legs — this lets a caller whose prototype -# lacks a leg still build a device-following identity via `one(similar_map(...), codomain, domain)`. -function similar_map(prototype::AbstractITensor, eltype::Type, codomain, domain) - raw = TensorAlgebra.similar_map( - unnamed(prototype), eltype, unnamed.(codomain), unnamed.(domain) - ) - return nameddims(raw, (name.(codomain)..., name.(domain)...)) -end -function similar_map(prototype::AbstractITensor, codomain, domain) - return similar_map(prototype, scalartype(prototype), codomain, domain) -end - -# From-scratch identity map: a dense identity embedded onto the `codomain`/`domain` index -# partition via checked `project`, so the index axes select the backend (dense, graded, -# `TensorMap`). Unlike `one(a, codomain, domain)` it needs no prototype tensor, so it is the -# right primitive when only the indices and an element type are in hand (e.g. `op("I")`). -function id(eltype::Type, codomain, domain) - m = Matrix{eltype}(LinearAlgebra.I, prod(length, codomain), prod(length, domain)) - return TensorAlgebra.project(m, Tuple(codomain), Tuple(domain)) -end - -# Dense Kronecker copy (`delta`) tensor over the index axes. Dense-only: a super-diagonal -# generally cannot be embedded while preserving a nontrivial symmetry, so graded/`TensorMap` -# callers that want an order-2 identity build it via `id`/`one` at the callsite instead. -delta(eltype::Type, is::Tuple) = diagonaltensor(ones(eltype, minimum(length, is)), is) -delta(eltype::Type, is::Index...) = delta(eltype, is) -delta(eltype::Type, is::AbstractVector{<:Index}) = delta(eltype, Tuple(is)) -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) - 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. -# -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) - -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 - -# 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 — see api_migration_map.md; the current callsites pass -# `cutoff = nothing` (full decomposition). -function _eigh( - 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) - u = only(commoninds(D, U)) # eigenvalue index shared with U - t = only(uniqueinds(D, U)) # D's other (independent) index - D = replaceinds(D, t => ITensorBase.prime(u)) - return D, U -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`) — -# 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(dag(U)) == m`, a true matrix sqrt). -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...) - 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)) - isnothing(cutoff) || iszero(cutoff) || (trunc &= MAK.truncerror(; rtol = sqrt(cutoff), p = 2)) - return trunc -end - -# Split `a`'s indices into (left, right) by the requested `linds`, matching by name: -# the caller's `Index` objects may carry the other endpoint's (dual) space on a graded -# backend, so `Index`-equality filtering silently drops them (see the name-matching -# note on the index-set algebra above). Both groups are `a`'s own indices, so the -# spaces handed to the factorization are `a`'s. -function _bipartition_inds(a::AbstractITensor, linds) - lnames = name.(cat_inds(linds...)) - allinds = collect(inds(a)) - left = filter(i -> name(i) ∈ lnames, allinds) - right = filter(i -> name(i) ∉ lnames, allinds) - return left, right -end - -# Legacy `factorize(a, linds...; ortho, cutoff, maxdim, tags)` splits `a` into `L * R` -# 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...; - 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 - 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 - 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 - - - -# -# Index fusion. The legacy `combiner` is retired (no compat shim); call sites fuse -# index groups with the next-gen `matricize(t, row_inds => row_name, col_inds => -# col_name)` (minting each fused name via `uniquename(IndexName)`) or build a fused -# identity with `Base.one`, both of which are graded-capable. `matricize` is the -# `TensorAlgebra` generic, extended by ITensorBase for named tensors. -using TensorAlgebra: matricize - -# -# Storage / element type accessors. -# -# `scalartype` is re-exported above. `datatype` is the underlying storage array type -# (used by `adapt`); `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) - -# 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 -# 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 - -# 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), ...) -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)))) - 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)) - order = (shared..., sinds...) - cur = collect(ITensorBase.dimnames(t)) - perm = [findfirst(==(name(o)), cur) for o in order] - a = permutedims(unnamed(t), perm) - ranges = ( - Base.OneTo.(length.(shared))..., - ntuple(k -> (offsets[k] + 1):(offsets[k] + length(sinds[k])), length(sinds))..., - ) - out[ranges...] .= a - offsets .+= length.(sinds) - end - 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 -Algorithm{Alg}(kwargs::NamedTuple) where {Alg} = Algorithm{Alg, typeof(kwargs)}(kwargs) -Algorithm{Alg}(; kwargs...) where {Alg} = Algorithm{Alg}((; kwargs...)) -Algorithm(alg::Symbol; kwargs...) = Algorithm{alg}(; kwargs...) -Algorithm(alg::AbstractString; kwargs...) = Algorithm(Symbol(alg); kwargs...) -Algorithm(alg::Algorithm) = alg -function Base.getproperty(alg::Algorithm, name::Symbol) - return if name === :kwargs - getfield(alg, :kwargs) - else - 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 - -# -# 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)) - isempty(s) || (i = ITensorBase.settag(i, s, "")) - 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) - 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 - -# 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/MessagePassing/abstractbeliefpropagationcache.jl b/src/MessagePassing/abstractbeliefpropagationcache.jl index 247b1af..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( @@ -41,8 +41,8 @@ for f in [ :(bp_factors), :(default_bp_maxiter), :(virtualinds), - :(ITensors.datatype), - :(ITensors.scalartype), + :(datatype), + :(scalartype), :(maxvirtualdim), :(default_message), :(siteinds), @@ -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 224d993..ab72bb6 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 + 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}) @@ -161,20 +159,18 @@ 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 = ITensors.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 = 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..6e4e5f6 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)) @@ -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 @@ -195,12 +195,12 @@ 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(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. 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 @@ -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) @@ -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 @@ -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,l=$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 = replaceinds(L, b => bnew), replaceinds(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,l=$(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 = replaceinds(L, b => bnew), replaceinds(R, b => bnew) out[i] = R out[i - 1] *= L end @@ -647,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 @@ -656,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 dd30010..121dee5 100644 --- a/src/MessagePassing/loopcorrection.jl +++ b/src/MessagePassing/loopcorrection.jl @@ -36,20 +36,20 @@ 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, 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) 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) - t = replaceind(t, t_ind, linds_sim[t_ind_pos]) + t_ind_pos = findfirst(==(t_ind), linds) + t = replaceinds(t, t_ind => linds_sim[t_ind_pos]) setindex_preserve!(bpc, t, src(e)) end push!(updated_es, e) @@ -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. - row_inds = Index[only(commoninds(me, [l])) for l in linds] - col_inds = Index[only(commoninds(mer, [l])) for l in linds_sim] + # relative to the passed indices, so the columns go in `conj`ed. + 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, dag.(col_inds))) + ap = adapt_like(me, identity_tensor(row_inds, conj.(col_inds))) ap = ap - me * mer push!(antiprojectors, ap) end @@ -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/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..8d7cf87 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") @@ -93,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 444f66f..d39a640 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,22 +20,21 @@ 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) - 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 -# on one endpoint and dual (conjugated) on the other, so the two `Index` objects -# differ even though they name the same bond. -function ITensors.uniqueinds(tn::AbstractTensorNetwork, v) - tv_inds = Index[i for i in inds(tn[v])] +# `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 = 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, [inds(tn[vn]) for vn in vns]) + return setdiff(tv_inds, neighbor_inds) end function setindex_preserve!(tn::AbstractTensorNetwork, value::ITensor, vertex) @@ -50,12 +48,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 +75,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 +92,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 @@ -124,7 +122,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 ], ), @@ -139,7 +137,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..c05a934 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) @@ -86,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 87cfe0a..f53bc6c 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} @@ -55,17 +54,17 @@ 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)) - return isempty(auxinds) ? tdag : replaceinds(tdag, prime.(auxinds), auxinds) + 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)) # 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, @@ -79,10 +78,10 @@ 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, 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.(conj.(linds))) return adapt_like(tns, one(zeros(scalartype(tns), cod..., dom...), cod, dom)) end @@ -124,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 @@ -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 @@ -220,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/TensorNetworks/tensornetworkstate_constructors.jl b/src/TensorNetworks/tensornetworkstate_constructors.jl index c1a5a4a..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), [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..bc452b9 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_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(ITensors.contract(tn_tensors; sequence = seq)) + return scalar(contract_network(tn_tensors; sequence = seq)) end -function ITensors.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 ITensors.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 ITensors.contract(tn::AbstractTensorNetwork; alg = "exact", kwargs...) - return ITensors.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/contraction_sequences.jl b/src/contraction_sequences.jl index bcbe3c6..cac3eb3 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 @@ -8,15 +7,14 @@ 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 -# 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 = inds(t) (prune_tensors && is_trivial_tensor(t)) ? empty(is) : is end end @@ -24,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(dim(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 @@ -42,10 +40,10 @@ end #OMEinsumContractionOrders helpers function to_eincode(tensors::Vector{<:ITensor}) - ixs = map(t -> collect(name.(inds(t))), tensors) + ixs = map(inds, 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)) + 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/src/entanglement.jl b/src/entanglement.jl index 0a93c0c..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 ITensors.array(ITensors.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 """ @@ -87,8 +80,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 003f9f7..77742c5 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[] @@ -17,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 @@ -72,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/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..ba7404c 100644 --- a/src/inner.jl +++ b/src/inner.jl @@ -31,50 +31,50 @@ 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()) ) 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 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/itensors.jl b/src/itensors.jl new file mode 100644 index 0000000..3352b04 --- /dev/null +++ b/src/itensors.jl @@ -0,0 +1,149 @@ +# 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, uniqueinds +import MatrixAlgebraKit as MAK +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, matricize, project, scalar, tryproject + +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 ψ + 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 + +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) + +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)) + +diaglength(a::AbstractArray) = minimum(size(a)) +function diagstride(a::AbstractArray) + s = 1 + p = 1 + for i in 1:(ndims(a) - 1) + p *= size(a, i) + s += p + end + return s +end +function diagindices(a::AbstractArray) + maxdiag = LinearIndices(a)[CartesianIndex(ntuple(Returns(diaglength(a)), ndims(a)))] + return 1:diagstride(a):maxdiag +end +diagview(a::AbstractArray) = @view a[diagindices(a)] + +function diagonaltensor(diag::AbstractVector, ax::Tuple{Vararg{AbstractUnitRange}}) + a = similar(diag, ax) + fill!(a, zero(eltype(a))) + diagview(a) .= diag + return a +end +function diagonaltensor( + diag::AbstractVector, + is::Tuple{NamedUnitRange, Vararg{NamedUnitRange}} + ) + return nameddims(diagonaltensor(diag, unnamed.(is)), name.(is)) +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)) +delta(is::Tuple) = delta(Float64, is) +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` factors 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 * ψ) + +function itensor_trunc(; maxdim = nothing, cutoff = nothing) + trunc = MAK.notrunc() + isnothing(maxdim) || (trunc &= MAK.truncrank(maxdim)) + isnothing(cutoff) || iszero(cutoff) || (trunc &= MAK.truncerror(; rtol = sqrt(cutoff), p = 2)) + return trunc +end + +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} + eltype(T) === elt && return T + return nameddims(convert(AbstractArray{elt}, unnamed(T)), ITensorBase.dimnames(T)) +end + +hasqns(i::Index) = conj(unnamed(i)) != unnamed(i) +hasqns(t::AbstractITensor) = any(hasqns, inds(t)) + +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(eltype(t1), length.(target)) + offsets = zeros(Int, length(out_inds)) + for p in pairs + t, sinds = first(p), last(p) + a = ITensorBase.unname(t, (shared..., sinds...)) + ranges = ( + Base.OneTo.(length.(shared))..., + ntuple(k -> (offsets[k] + 1):(offsets[k] + length(sinds[k])), length(sinds))..., + ) + out[ranges...] .= a + offsets .+= length.(sinds) + end + return out[target...] +end + +struct Algorithm{Alg, Kwargs <: NamedTuple} + kwargs::Kwargs +end +Algorithm{Alg}(kwargs::NamedTuple) where {Alg} = Algorithm{Alg, typeof(kwargs)}(kwargs) +Algorithm{Alg}(; kwargs...) where {Alg} = Algorithm{Alg}((; kwargs...)) +Algorithm(alg::Symbol; kwargs...) = Algorithm{alg}(; kwargs...) +Algorithm(alg::AbstractString; kwargs...) = Algorithm(Symbol(alg); kwargs...) +Algorithm(alg::Algorithm) = alg +function Base.getproperty(alg::Algorithm, name::Symbol) + return if name === :kwargs + getfield(alg, :kwargs) + else + getfield(getfield(alg, :kwargs), name) + end +end +macro Algorithm_str(s) + return :(Algorithm{$(Expr(:quote, Symbol(s)))}) +end + +settags(i::Index, p::Pair) = ITensorBase.settag(i, first(p), last(p)) +function settags(i::Index, tags) + for (k, v) in tags + i = ITensorBase.settag(i, k, v) + end + return i +end 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 43b8ac1..c3c01a4 100644 --- a/src/rdm.jl +++ b/src/rdm.jl @@ -1,5 +1,5 @@ function normalize_rdm(ρ::ITensor) - return ρ / ITensors.tr(ρ) + return ρ / tr(ρ, operator_inds(ρ)...) end """ @@ -32,11 +32,10 @@ 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...) - ρ = contract(ρ_tensors; sequence = seq) + ρ = contract_network(ρ_tensors; sequence = seq) if normalize ρ = normalize_rdm(ρ) end @@ -59,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 635d447..aa61a5f 100644 --- a/src/sampling.jl +++ b/src/sampling.jl @@ -23,16 +23,16 @@ 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_network(tensors; sequence = seq) - ρ_tr = ITensors.tr(ρ) + ρ_tr = tr(ρ, operator_inds(ρ)...) ρ *= 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) 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(ψ)) @@ -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), net_inds, (inds(m) for m in outgoing_mps if m !== mt)...) setmessage!(norm_bmps_cache, e, ITensor[mt, bra_tensor(mt, aux)]) end @@ -237,16 +237,16 @@ 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) - ρ_tr = ITensors.tr(ρ) + ρ = contract_network(ts; sequence = seq) + ρ_tr = tr(ρ, operator_inds(ρ)...) 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) 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/siteinds.jl b/src/siteinds.jl index 4713612..6d614e4 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) @@ -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 7c2bda5..1f770b3 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 @@ -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 = 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 - ψ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. @@ -48,16 +48,16 @@ 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)))] - ψvsrc = replaceinds(ψvsrc, v, new_edge_ind) - ψvdst = replaceinds(ψvdst, u, new_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) 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/src/truncate.jl b/src/truncate.jl index 4cec666..830fa53 100644 --- a/src/truncate.jl +++ b/src/truncate.jl @@ -5,11 +5,11 @@ 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 -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..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 ITensors.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...) - elseif e == Float32 - m = adapt(Vector{Float64}, m) - D, U = ITensors.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...) - 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]) 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..047abb3 --- /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_network([ψ[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 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_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 9d09181..0c09d36 100644 --- a/test/test_constructors.jl +++ b/test/test_constructors.jl @@ -3,12 +3,10 @@ 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 -using TensorNetworkQuantumSimulator: dag, prime +# `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: 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 = 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) @@ -36,9 +34,9 @@ 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 ITensors.contract(ψdag; alg = "exact") ≈ conj(ITensors.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 da0e44b..a0229d9 100644 --- a/test/test_contraction_sequences.jl +++ b/test/test_contraction_sequences.jl @@ -1,10 +1,9 @@ @eval module $(gensym()) -using ITensorBase: Index, name +using ITensorBase: Index using Random using TensorNetworkQuantumSimulator const TNQS = TensorNetworkQuantumSimulator -# `random_itensor`/`contract`/`scalar` come from TNQS's compat layer (see test_constructors). -import TensorNetworkQuantumSimulator as ITensors +# `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 @@ -18,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 = ITensors.random_itensor(i, j) - B = ITensors.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) - @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). @@ -50,21 +50,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_network(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_network(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 = 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] - @test ITensors.contract(open_tensors; sequence = seq_open) ≈ - ITensors.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