Remove the compat submodule and shrink the itensors.jl shim onto the upstream API

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"

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):
 

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)])

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")

examples/boundarymps.jl

@@ -1,5 +1,4 @@
 using TensorNetworkQuantumSimulator
-using TensorNetworkQuantumSimulator: ITensors
 
 using Random
 Random.seed!(1634)

examples/heavyhexIsing_dynamics.jl

@@ -1,6 +1,6 @@
 using TensorNetworkQuantumSimulator
 using Statistics
-using TensorNetworkQuantumSimulator: ITensor, ITensors
+using TensorNetworkQuantumSimulator: ITensor
 
 function main()
     #Define the lattice

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

examples/loopcorrections.jl

@@ -1,5 +1,4 @@
 using TensorNetworkQuantumSimulator
-using TensorNetworkQuantumSimulator: ITensors
 
 using LinearAlgebra: norm
 

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
 

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)