diff --git a/src/ITensors.jl b/src/ITensors.jl index 408da78..e25b668 100644 --- a/src/ITensors.jl +++ b/src/ITensors.jl @@ -14,9 +14,9 @@ module ITensors # re-exports `NDTensors`, so `using ITensorsITensorBaseCompat.ITensors: scalartype` must # also work — and bring in the ones this submodule's own methods build on # (`scalartype` / `datatype`). +import ..NDTensors: map_diag, map_diag! import Base: truncate -using ..NDTensors: @Algorithm_str, Algorithm, data, datatype, dense, denseblocks, map_diag, - map_diag!, scalartype +using ..NDTensors: @Algorithm_str, Algorithm, data, datatype, dense, denseblocks, scalartype include("itensor.jl") @@ -38,6 +38,8 @@ export qr, svd, eigen, factorize, factorize_svd, # Diagonal manipulation map_diag, map_diag!, + # Operator exponential + exp, # Storage / element-type accessors scalartype, datatype, array, data, # Dense / quantum-number no-ops diff --git a/src/NDTensors.jl b/src/NDTensors.jl index 3b32b2b..1e468c1 100644 --- a/src/NDTensors.jl +++ b/src/NDTensors.jl @@ -5,7 +5,7 @@ module NDTensors using BackendSelection: @Algorithm_str, Algorithm -using ITensorBase: ITensorBase, AbstractITensor, inds, unnamed +using ITensorBase: ITensorBase, AbstractITensor, unnamed # # Storage / element type accessors. `scalartype` is the scalar (element) type; @@ -21,35 +21,12 @@ denseblocks(T::AbstractITensor) = T dense(T::AbstractITensor) = T # -# Diagonal manipulation. Legacy `map_diag(f, T)` applies `f` to the diagonal of a -# (diagonal-like) tensor; used on factorization spectra (singular values / -# eigenvalues). Reproduced via the same diagonal machinery as `delta`. -_diagcartesian(arr, k) = CartesianIndex(ntuple(Returns(k), ndims(arr))) -function map_diag(f, T::AbstractITensor) - arr = copy(unnamed(T)) - for k in 1:minimum(size(arr)) - idx = _diagcartesian(arr, k) - arr[idx] = f(arr[idx]) - end - return arr[inds(T)...] -end -function map_diag!(f, T::AbstractITensor) - arr = unnamed(T) - for k in 1:minimum(size(arr)) - idx = _diagcartesian(arr, k) - arr[idx] = f(arr[idx]) - end - return T -end -# Out-of-place-into-`dest` form `map_diag!(f, dest, src)`: write `f` of `src`'s diagonal -# onto `dest`'s diagonal (TNQS calls it with `dest === src` for an in-place diagonal map). -function map_diag!(f, dest::AbstractITensor, src::AbstractITensor) - d, s = unnamed(dest), unnamed(src) - for k in 1:minimum(size(s)) - d[_diagcartesian(d, k)] = f(s[_diagcartesian(s, k)]) - end - return dest -end +# Diagonal manipulation. Legacy `map_diag(f, T)` / `map_diag!` apply `f` to a tensor's +# diagonal (used on factorization spectra). These names belong to `NDTensors` in the real +# ecosystem, so the generic functions live here; the `AbstractITensor` methods are defined +# in the sibling `ITensors` submodule, the layer that owns the tensor type. +function map_diag end +function map_diag! end # The algorithm dispatch tag (legacy `Algorithm` / `@Algorithm_str`) comes from # BackendSelection.jl, which NDTensors re-exports in the real ecosystem; imported above diff --git a/src/SiteTypes.jl b/src/SiteTypes.jl index 0de9aec..be4ef25 100644 --- a/src/SiteTypes.jl +++ b/src/SiteTypes.jl @@ -14,7 +14,6 @@ # (`src/lib/SiteTypes/`), re-exported by `ITensors`; reproduced as a submodule here. module SiteTypes -using ..ITensors: combinedind, combiner, dag, inds using ITensorBase: ITensorBase, AbstractITensor, Index using LinearAlgebra: LinearAlgebra @@ -140,27 +139,6 @@ function op(::OpName"CPHASE", ::SiteType"S=1/2"; ϕ) return ComplexF64[1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 cis(ϕ)] end -# TYPE PIRACY (temporary): extends `Base.exp` for an operator `ITensor` (matricize over the -# index/prime pairs, exponentiate, rebuild). Legacy ITensors provided `exp(::ITensor)`; gates -# defined as `exp` of a Hamiltonian operator rely on it (e.g. a user -# `op(::OpName"MyZRot", ...) = exp(-im θ/2 * op("Z", s))`). To de-pirate: make this compat-owned -# (an `exp` in this module, not a `Base.exp` method), inferring the prime-pair codomain/domain and -# forwarding to ITensorBase's matricization `exp(a, dimnames_codomain, dimnames_domain)` — which -# already exists and is graded-capable, so this dense combiner-based version goes away. Not an -# upstream candidate (the upstream matricization `exp` is the target, not a `Base.exp(::ITensor)`). -function Base.exp(t::AbstractITensor) - p0 = filter(i -> ITensorBase.plev(i) == 0, collect(inds(t))) - isempty(p0) && error("exp(::ITensor) expects indices paired as (i, prime(i))") - p1 = map(ITensorBase.prime, p0) - cr, cc = combiner(Tuple(p1)), combiner(Tuple(p0)) - tc = (t * cr) * cc - rci, cci = combinedind(cr), combinedind(cc) - d = length(rci) - M = ComplexF64[tc[rci => a, cci => b] for a in 1:d, b in 1:d] - E = exp(M)[rci, cci] - return (E * dag(cr)) * dag(cc) -end - export @OpName_str, @SiteType_str, OpName, SiteType, op, state end diff --git a/src/itensor.jl b/src/itensor.jl index c022f09..7885adc 100644 --- a/src/itensor.jl +++ b/src/itensor.jl @@ -3,9 +3,9 @@ # # Strategy (see ITensorDevelopmentPlans api_migration_map.md): # - Names below are thin wrappers over ITensorBase / TensorAlgebra / MatrixAlgebraKit. -# - `combiner`, the factorization return shapes, `map_diag`, the operator/SiteType -# system, and the boundary-MPS (ITensorMPS) paths are NOT wrapped here; they need -# callsite translation or upstream stack work and are tracked separately. +# - The factorization return shapes, the operator/SiteType system, and the boundary-MPS +# (ITensorMPS) paths are NOT wrapped here; they need callsite translation or upstream +# stack work and are tracked separately. # # ITensorBase keeps most of this API internal (unexported), so we reach for the # qualified names and re-publish the legacy spellings into this namespace. @@ -162,13 +162,6 @@ end itensor(array, is) = array[is...] itensor(array, is...) = array[is...] -# TYPE PIRACY (temporary, compat-owned — NOT an upstream candidate): adds a rank-0 -# `ITensor(x::Number)` constructor, which ITensorBase deliberately omits and does not plan to -# support. Legacy ITensors uses it as a multiplicative identity to seed a product accumulator -# (`out = ITensor(1); out *= t; ...`). Kept here for now; the accumulator call sites get rewritten -# to a different pattern later, retiring this method rather than upstreaming it. -ITensorBase.ITensor(x::Number) = nameddims(fill(x), ()) - # 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...) @@ -487,21 +480,6 @@ end # the `NDTensors` submodule and are imported into this one.) array(T::AbstractITensor) = unnamed(T) -# TYPE PIRACY (temporary, compat-owned — NOT an upstream candidate): extends -# `Adapt.adapt_structure` for `AbstractITensor` with an eltype target. Using -# `Adapt.adapt_structure` for eltype *conversion* is an abuse of Adapt.jl (Adapt is for -# storage/device adaptation, not changing the scalar type), so this does not belong upstream. -# Kept here for now; the eltype-conversion call sites get rewritten with a different pattern -# later, retiring this shim rather than upstreaming it. -# -# Legacy `adapt(eltype)(t)` converts an ITensor's scalar (element) type. ITensorBase's -# Adapt integration adapts the storage array/device type but leaves the element type -# alone, so reproduce the eltype conversion for a `Number` target (used by -# `adapt(eltype)(state(...))` to build typed product states). -function Adapt.adapt_structure(::Type{elt}, T::AbstractITensor) where {elt <: Number} - 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) @@ -568,29 +546,57 @@ function settags(i::Index, d::AbstractDict) end return i end -# TYPE PIRACY (temporary, compat-owned — NOT an upstream candidate): the two methods below add -# `ITensorBase.Index` constructors taking a tag string / tag dict, a legacy positional-tag form -# ITensorBase does not plan to support (the next-gen spelling passes tags via the `tags` keyword -# argument). Kept here for now; the call sites get modernized to the `tags` kwarg later, retiring -# these methods rather than upstreaming them. -# -# Legacy positional tagged-index constructor `Index(dim, "tag")`. -ITensorBase.Index(dim::Integer, tagstr::AbstractString) = settags(Index(dim), tagstr) -# Build a fresh index carrying a tag dictionary (legacy `Index(dim, tags(i))`, where -# the next-gen `tags` returns a `Dict{String, String}`). -function ITensorBase.Index(dim::Integer, tags::AbstractDict) - i = Index(dim) - for (k, v) in tags - 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 +# +# Diagonal manipulation. The `map_diag` / `map_diag!` generics belong to `NDTensors` +# (imported into this module); their `AbstractITensor` methods live here at the ITensor +# layer. +_diagcartesian(arr, k) = CartesianIndex(ntuple(Returns(k), ndims(arr))) +function map_diag(f, T::AbstractITensor) + arr = copy(unnamed(T)) + for k in 1:minimum(size(arr)) + idx = _diagcartesian(arr, k) + arr[idx] = f(arr[idx]) + end + return arr[inds(T)...] +end +function map_diag!(f, T::AbstractITensor) + arr = unnamed(T) + for k in 1:minimum(size(arr)) + idx = _diagcartesian(arr, k) + arr[idx] = f(arr[idx]) + end + return T +end +# Out-of-place-into-`dest` form `map_diag!(f, dest, src)`: write `f` of `src`'s diagonal +# onto `dest`'s diagonal (TNQS calls it with `dest === src` for an in-place diagonal map). +function map_diag!(f, dest::AbstractITensor, src::AbstractITensor) + d, s = unnamed(dest), unnamed(src) + for k in 1:minimum(size(s)) + d[_diagcartesian(d, k)] = f(s[_diagcartesian(s, k)]) + end + return dest +end + +# +# Owned `exp` (avoids type piracy: `AbstractITensor` is not ours, so we do not add a +# `Base.exp` method for it). Legacy `ITensors.exp(::ITensor)` exponentiates an operator +# ITensor over its `(i, prime(i))` index pairs by forwarding to ITensorBase's +# graded-capable matricization `Base.exp(a, codomain, domain)`. Other arguments forward +# to `Base.exp`. +exp(x) = Base.exp(x) +function exp(t::AbstractITensor) + p0 = filter(i -> ITensorBase.plev(i) == 0, collect(inds(t))) + isempty(p0) && error("exp(::ITensor) expects indices paired as (i, prime(i))") + p1 = map(ITensorBase.prime, p0) + return Base.exp(t, Tuple(p1), Tuple(p0)) +end + # TODO (small inline residue — can't be a drop-in shim): # - `contract` / `inner` / `truncate`: consumers *extend* these (method definitions), # so the call sites drop the `ITensors.` qualifier to extend the generics this