Move to ITensorBase 0.14 and TensorAlgebra 0.20

Project.toml

@@ -27,7 +27,7 @@ VectorInterface = "409d34a3-91d5-4945-b6ec-7529ddf182d8"
 Adapt = "4.3.0"
 Dictionaries = "0.4"
 Graphs = "1.8.0"
-ITensorBase = "0.13.7"
+ITensorBase = "0.14.0"
 KrylovKit = "0.10.2"
 LinearAlgebra = "1.10"
 MatrixAlgebraKit = "0.6.8"
@@ -37,7 +37,7 @@ SimpleGraphAlgorithms = "0.6.0"
 SimpleGraphConverter = "0.1.0"
 SplitApplyCombine = "1.2.3"
 StatsBase = "0.34.4"
-TensorAlgebra = "0.17.6"
+TensorAlgebra = "0.20.0"
 TypeParameterAccessors = "0.3.10, 0.4"
 VectorInterface = "0.5"
 julia = "1.10"

src/Ops.jl

@@ -8,7 +8,7 @@ module Ops
 
 using ITensorBase: ITensorBase, Index, id
 using TensorAlgebra: project
-using ..TensorNetworkQuantumSimulator: project_aux
+using ..TensorNetworkQuantumSimulator: maybe_project_aux
 
 export state, op, OpName, SiteType, @OpName_str, @SiteType_str
 
@@ -33,8 +33,8 @@ function state(name::AbstractString, i::Index)
     return state(v, i)
 end
 # Vector form (legacy `ITensor(v, i)` for a state vector): the state vector as an `ITensor`
-# over `i`, through the `project_aux` utility.
-state(v::AbstractVector{<:Number}, i::Index) = project_aux(v, i)
+# over `i`, through the `maybe_project_aux` utility.
+state(v::AbstractVector{<:Number}, i::Index) = maybe_project_aux(v, i)
 
 #
 # Operators. Legacy ITensors exposes operators through the `OpName` / `SiteType`
@@ -59,7 +59,7 @@ end
 
 # Embed a `d^n × d^n` operator matrix (computational basis, first site most
 # significant) into an `ITensor` with codomain `prime.(sites)` (outputs) and
-# domain `sites` (inputs). Routed through `project_aux` (like `state`): a graded
+# domain `sites` (inputs). Routed through `maybe_project_aux` (like `state`): a graded
 # operator that is odd under the site grading (e.g. `X`/`Y`, bare `c`/`c†`) has no
 # parity-block-diagonal `(prime.(sites), sites)` map, so it gets a trailing auxiliary
 # charge leg instead of throwing an `InexactError`. Even operators and dense (ungraded)
@@ -69,7 +69,7 @@ function _op_matrix_to_itensor(M::AbstractMatrix, sites::Tuple)
     n = length(sites)
     A = reshape(Matrix{ComplexF64}(M), (reverse(ds)..., reverse(ds)...))
     A = permutedims(A, (reverse(1:n)..., reverse((n + 1):(2n))...))
-    return project_aux(A, ITensorBase.prime.(sites), sites)
+    return maybe_project_aux(A, ITensorBase.prime.(sites), sites)
 end
 
 # Top-level `op(name, sites...; kwargs...)`. The identity is dimension-general (used

src/itensors.jl

@@ -4,28 +4,27 @@
 import MatrixAlgebraKit as MAK
 using Adapt: Adapt
 using ITensorBase: ITensorBase, ITensor, Index, inds, noprime, plev, prime, space, unnamed
-using TensorAlgebra: project, tryproject
+using TensorAlgebra: TensorAlgebra, project, tryproject
 
 # Project `a` onto the symmetry-restricted space given by `codomain`/`domain`, mirroring
 # `TensorAlgebra.project`/`tryproject` (three-argument operator form; `a` is indexed positionally
-# as `(codomain..., domain...)`). When the plain projection is charge-forbidden (a parity-odd state
-# or operator), reshaping to a trailing length-1 axis lets `project` derive a named auxiliary index
-# that absorbs the residual charge, so the returned ITensor carries that extra aux leg last instead
-# of throwing.
-function project_aux(a::AbstractArray, codomain, domain)
-    return @something tryproject(a, codomain, domain) project(
-        reshape(a, (size(a)..., 1)), codomain, domain
-    )
+# as `(codomain..., domain...)`). The `maybe_` is the whole point and is what distinguishes this from
+# `TensorAlgebra.project_aux`, which appends the auxiliary axis unconditionally: here the plain
+# projection is tried first, so an allowed state or operator keeps its natural rank, and only a
+# charge-forbidden one (a parity-odd state or operator) falls back to deriving a named auxiliary
+# index that absorbs the residual charge, carried last, instead of throwing.
+function maybe_project_aux(a::AbstractArray, codomain, domain)
+    return @something tryproject(a, codomain, domain) TensorAlgebra.project_aux(a, codomain, domain)
 end
-project_aux(a::AbstractArray, codomain) = project_aux(a, codomain, ())
-project_aux(v::AbstractVector{<:Number}, i::Index) = project_aux(v, (i,))
+maybe_project_aux(a::AbstractArray, codomain) = maybe_project_aux(a, codomain, ())
+maybe_project_aux(v::AbstractVector{<:Number}, i::Index) = maybe_project_aux(v, (i,))
 
 # Build two tensors that share one contractible auxiliary leg, so `t1 * t2` gives the fermion
-# string of e.g. `c†ᵢcⱼ` (a `c†` and a `c`) with no `flip`. Projecting `a1` with a trailing dummy
-# axis derives a named aux `Index` that absorbs `a1`'s residual charge, and feeding that same
-# `Index` into `a2`'s domain makes both carry it, so `t1 * t2` contracts them by name.
+# string of e.g. `c†ᵢcⱼ` (a `c†` and a `c`) with no `flip`. `TensorAlgebra.project_aux` on `a1`
+# derives a named aux `Index` that absorbs `a1`'s residual charge, and feeding that same `Index` into
+# `a2`'s domain makes both carry it, so `t1 * t2` contracts them by name.
 function project_pair(a1::AbstractArray, codomain1, domain1, a2::AbstractArray, codomain2, domain2)
-    p1 = project(reshape(a1, (size(a1)..., 1)), codomain1, domain1)
+    p1 = TensorAlgebra.project_aux(a1, codomain1, domain1)
     aux = last(inds(p1))
     p2 = project(reshape(a2, (size(a2)..., 1)), codomain2, (domain2..., aux))
     return p1, p2
@@ -34,7 +33,7 @@ end
 function onehot(eltype::Type, (i, p)::Pair{<:Index})
     v = zeros(eltype, length(i))
     v[p] = one(eltype)
-    return project_aux(v, i)
+    return maybe_project_aux(v, i)
 end
 onehot(p::Pair{<:Index}) = onehot(Float64, p)