diff --git a/Project.toml b/Project.toml index 1f524ee..a99bd82 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "GradedArrays" uuid = "bc96ca6e-b7c8-4bb6-888e-c93f838762c2" -version = "0.13.7" +version = "0.13.8" authors = ["ITensor developers and contributors"] [workspace] @@ -44,7 +44,7 @@ SUNRepresentations = "0.3, 0.4" SparseArraysBase = "0.10" SplitApplyCombine = "1.2.3" StridedViews = "0.5" -TensorAlgebra = "0.16" +TensorAlgebra = "0.17" TensorKit = "0.17" TensorKitSectors = "0.3" julia = "1.10" diff --git a/src/abeliangradedarray.jl b/src/abeliangradedarray.jl index 58a6158..d627f32 100644 --- a/src/abeliangradedarray.jl +++ b/src/abeliangradedarray.jl @@ -485,15 +485,26 @@ function LinearAlgebra.isdiag(A::AbelianGradedMatrix) return true end +# Throw unless `sz1` and `sz2` are equal ignoring trailing length-1 axes (an axis beyond one +# size's rank counts as length 1, mirroring `Base.size(A, d)` for `d > ndims(A)`), guarding a +# `reshape` against silently reinterpreting same-length data of a genuinely different shape. +# Vendored from TensorAlgebra rather than reused because it is not part of its public API. +function check_project_size(sz1::Dims, sz2::Dims) + all(i -> get(sz1, i, 1) == get(sz2, i, 1), 1:max(length(sz1), length(sz2))) || throw( + DimensionMismatch("sizes $sz1 and $sz2 differ beyond trailing length-1 axes") + ) + return nothing +end + # Orthogonal projection of a dense source into the symmetry-allowed subspace. # Magnitude-blind: forbidden-block entries of `src` are dropped without inspection. -# Use `TensorAlgebra.checked_projectto!` to verify the discarded weight is small. +# The `TensorAlgebra.project` wrapper verifies the discarded weight is small. function TensorAlgebra.projectto!(dest::AbelianGradedArray, src::AbstractArray) - size(dest) == size(src) || throw( - DimensionMismatch( - "projectto!: dest has size $(size(dest)), src has size $(size(src))" - ) - ) + # Reshape `src` to `size(dest)` (a no-op when the ranks already match), so a lower-rank + # `src` may omit trailing length-1 axes (e.g. an auxiliary flux-canceling leg); a genuine + # shape mismatch errors rather than reinterpreting the data. + check_project_size(size(src), size(dest)) + src = reshape(src, size(dest)) zero!(dest) for b in allowedblocks(axes(dest)) block_ranges = ntuple(d -> axes(dest, d)[Block(Int(Tuple(b)[d]))], ndims(dest)) @@ -502,20 +513,85 @@ function TensorAlgebra.projectto!(dest::AbelianGradedArray, src::AbstractArray) return dest end -# Compare via `Array(dest)` so the generic `isapprox(::AbstractArray, ::AbelianGradedArray)` -# path doesn't fall back to a `src - dest` broadcast that scalar-indexes the block storage. -# `kwargs` (e.g. `atol`/`rtol`) are forwarded to `isapprox`, matching the generic verb. -function TensorAlgebra.checked_projectto!( - dest::AbelianGradedArray, src::AbstractArray; kwargs... +# `allocate_project` with graded axes routes both the codomain-led and the (empty-codomain) +# domain-led cases to `allocate_project_graded`, taking the graded structure from whichever side +# is non-empty, the same two-entry split `similar_map` uses. +function TensorAlgebra.allocate_project( + src::AbstractArray, + codomain_axes::Tuple{GradedOneTo, Vararg{GradedOneTo}}, + domain_axes::Tuple{Vararg{GradedOneTo}} ) - TensorAlgebra.projectto!(dest, src) - isapprox(src, Array(dest); kwargs...) || - throw(InexactError(:checked_projectto!, typeof(dest), src)) - return dest + return allocate_project_graded(src, codomain_axes, domain_axes) +end +function TensorAlgebra.allocate_project( + src::AbstractArray, + codomain_axes::Tuple{}, + domain_axes::Tuple{GradedOneTo, Vararg{GradedOneTo}} + ) + return allocate_project_graded(src, codomain_axes, domain_axes) +end + +# The destination allocation, and where a trailing surplus axis in `src` gets its space derived. +# With no surplus this is plain `similar_map`; a single trailing surplus axis is an auxiliary leg +# appended as the last domain axis, its space derived (see `infer_aux_space`) so the result is +# symmetry-allowed. The trailing position matches how `stack` lays out an operator multiplet and +# the Julia convention that trailing length-1 axes are implicit. For the matricized result `X` the +# Gram matrix `X * X'` contracts the aux away (for a spin multiplet, `X * X' = S·S`, the Casimir). +function allocate_project_graded(src, codomain_axes, domain_axes) + nphys = length(codomain_axes) + length(domain_axes) + ndims(src) <= nphys && + return TensorAlgebra.similar_map(src, codomain_axes, domain_axes) + ndims(src) == nphys + 1 || throw( + ArgumentError( + "`project`: expected at most one trailing auxiliary axis beyond the $nphys \ + given axes, got a rank-$(ndims(src)) input" + ) + ) + aux = infer_aux_space(src, codomain_axes, domain_axes) + return TensorAlgebra.similar_map(src, codomain_axes, (domain_axes..., aux)) +end + +# The space of `src`'s trailing auxiliary axis, derived so the projected result is symmetry- +# allowed. Abelian sectors are one-dimensional, so each length-1 slice along the aux axis carries +# a single charge (`projected_charge`); contiguous equal charges merge into one sector of that +# multiplicity (matching the `TensorMap` backend), while non-contiguous repeats stay separate to +# preserve slice order. A multi-slice aux comes out as a direct sum (e.g. stacking `[S⁺, S⁻]` +# gives `[U1(2), U1(-2)]`, an MPO-virtual-leg structure). +function infer_aux_space(src::AbstractArray, codomain_axes, domain_axes) + aux_dim = length(codomain_axes) + length(domain_axes) + 1 + qs = map(eachslice(src; dims = aux_dim)) do slice + return projected_charge(slice, codomain_axes, domain_axes) + end + ps = Pair{eltype(qs), Int}[] + for q in qs + if isempty(ps) || first(ps[end]) != q + push!(ps, q => 1) + else + ps[end] = q => (last(ps[end]) + 1) + end + end + return gradedrange(ps) +end + +# Net charge of a dense operator, read from its dominant-magnitude entry: find the block +# holding that entry over the stored axes (domain dualized to match `zeros_map`/`similar_map`) +# and fuse that block's per-axis sectors, each with its axis's arrow applied (the same fusion +# `allowedblocks` is built on, so the charge lines up with which blocks `project` keeps). This is +# the abelian per-slice primitive `infer_aux_space` builds the aux space from; the general +# (possibly multi-sector) derivation lives in the `TensorMap` backend. +function projected_charge(src::AbstractArray, codomain_axes, domain_axes) + stored = (codomain_axes..., conj.(domain_axes)...) + src = reshape(src, length.(stored)) + I = Tuple(findmax(abs, src)[2]) + secs = map(stored, I) do ax, i + return eachsectoraxis(ax)[Int(BlockArrays.findblock(ax, i))] + end + return reduce(tensor_product, secs) end -# Materialize the graded array into a dense `Array` for the default -# `checked_projectto!`/`isapprox`-after path. +# Materialize the graded array into a dense `Array`. The checked projection verbs reach this +# through `convert(Array, dest)` to compare elementwise against a dense source, which would +# otherwise fall back to a `src - dest` broadcast that scalar-indexes the block storage. function Base.Array(a::AbelianGradedArray{T, <:Any, N}) where {T, N} dest = zeros(T, size(a)) for bI in eachblockstoredindex(a) @@ -727,17 +803,18 @@ end Construct an `AbelianGradedArray` by projecting the dense data of `a` onto the symmetry-allowed blocks of the graded axes `(ax1, axs...)`, via -`TensorAlgebra.checked_projectto!` (which errors if `a` has weight outside +`TensorAlgebra.project` (which errors if `a` has weight outside the allowed blocks). `a` is reshaped to `length.((ax1, axs...))` first, so a trailing size-1 bond can be supplied implicitly. Each axis carries its own arrow, so index with `dual`/`conj` axes to set duality. """ function Base.getindex(a::AbstractArray, ax1::GradedOneTo, axs::GradedOneTo...) dest_axes = (ax1, axs...) - a_reshaped = reshape(a, length.(dest_axes)) - return TensorAlgebra.checked_projectto!( - similar(a_reshaped, eltype(a), dest_axes), a_reshaped - ) + # Match `a` to the requested axes up to trailing length-1 axes: indexing selects exactly these + # axes, so the surplus-axis derivation branch of `project` must not trigger, and a genuine + # shape mismatch errors rather than reinterpreting the data. + check_project_size(size(a), length.(dest_axes)) + return TensorAlgebra.project(reshape(a, length.(dest_axes)), dest_axes) end # Disambiguate the single-axis case for a concrete `Array`: `Base.getindex(::Array, # ::AbstractUnitRange{<:Integer})` and the projection method above are otherwise equally diff --git a/src/abstractgradedarray.jl b/src/abstractgradedarray.jl index 4a55209..0a4bdf3 100644 --- a/src/abstractgradedarray.jl +++ b/src/abstractgradedarray.jl @@ -117,6 +117,11 @@ function TensorAlgebra.scale!(a::AbstractGradedArray, β::Number) return a end +# The `LinearAlgebra` spelling of blockwise scaling (the generic fallback +# scalar-indexes). +LinearAlgebra.rmul!(a::AbstractGradedArray, β::Number) = TensorAlgebra.scale!(a, β) +LinearAlgebra.lmul!(β::Number, a::AbstractGradedArray) = TensorAlgebra.scale!(a, β) + function TensorAlgebra.zero!(a::AbstractGradedArray) for bI in eachblockstoredindex(a) zero!(view(a, bI)) diff --git a/src/fusedgradedmatrix.jl b/src/fusedgradedmatrix.jl index 176cd31..6eea755 100644 --- a/src/fusedgradedmatrix.jl +++ b/src/fusedgradedmatrix.jl @@ -69,6 +69,27 @@ end # Block-diagonal by construction (one block per sector), so just check each block. LinearAlgebra.isdiag(A::FusedGradedMatrix) = all(LinearAlgebra.isdiag, A.blocks) +# Blockwise copy: the generic `AbstractArray` fallback copies elementwise, which +# scalar-indexes (disallowed for graded arrays). +function Base.copy(A::FusedGradedMatrix) + return FusedGradedMatrix(A.codomain, A.domain, map(copy, A.blocks)) +end + +# Materialize into a dense `Array` blockwise, like the `AbelianGradedArray` method +# (the generic fallback copies elementwise, which scalar-indexes). Sectors present +# in only one of the two axes have no stored block and stay zero. +function Base.Array(a::FusedGradedMatrix{T}) where {T} + dest = zeros(T, size(a)) + rowax, colax = axes(a) + rowsectors, colsectors = collect(keys(a.codomain)), collect(keys(a.domain)) + for (s, b) in pairs(a.blocks) + r = rowax[Block(findfirst(==(s), rowsectors))] + c = colax[Block(findfirst(==(s), colsectors))] + copyto!(view(dest, r, c), b) + end + return dest +end + # Block-diagonal by construction, so any matrix function `f(A) = blkdiag(f(blk_i))` for # each stored block — covers `sqrt`, `exp`, `log`, etc. Routes around the generic # `LinearAlgebra` impls that scalar-index for triangular / Hermitian detection. diff --git a/src/fusedgradedvector.jl b/src/fusedgradedvector.jl index 4713f87..810f282 100644 --- a/src/fusedgradedvector.jl +++ b/src/fusedgradedvector.jl @@ -149,6 +149,18 @@ end # ======================== Accessors ======================== +# Blockwise copyto!: the generic `AbstractArray` fallback copies elementwise, which +# scalar-indexes (disallowed for graded arrays). Also the write path for mutating a +# `MAK.diagview(::FusedGradedMatrix)` (whose blocks alias the matrix diagonals). +function Base.copyto!(dest::FusedGradedVector, src::FusedGradedVector) + keys(dest.blocks) == keys(src.blocks) || + throw(ArgumentError("`copyto!` requires matching sectors")) + for s in keys(src.blocks) + copyto!(dest.blocks[s], src.blocks[s]) + end + return dest +end + BlockArrays.blocklength(v::FusedGradedVector) = length(v.axis) function blocktype(::Type{<:FusedGradedVector{T, S, D}}) where {T, S, D} @@ -276,6 +288,18 @@ function Base.copy(v::FusedGradedVector) return FusedGradedVector(copy(v.axis), map(copy, v.blocks)) end +# Materialize into a dense `Array` blockwise, like the `AbelianGradedArray` method +# (the generic fallback copies elementwise, which scalar-indexes). +function Base.Array(v::FusedGradedVector{T}) where {T} + dest = zeros(T, size(v)) + ax = only(axes(v)) + sectors = collect(keys(v.axis)) + for (s, b) in pairs(v.blocks) + copyto!(view(dest, ax[Block(findfirst(==(s), sectors))]), b) + end + return dest +end + # ======================== LinearAlgebra ====================== function LinearAlgebra.norm(A::FusedGradedVector, p::Real = 2) diff --git a/src/gradedoneto.jl b/src/gradedoneto.jl index 34282e1..f2111dd 100644 --- a/src/gradedoneto.jl +++ b/src/gradedoneto.jl @@ -63,6 +63,9 @@ end trivial(g::GradedOneTo) = trivial(typeof(g)) TensorAlgebra.trivialrange(R::Type{<:GradedOneTo}) = trivial(R) +function TensorAlgebra.trivialrange(::Type{GradedOneTo{S}}, n::Integer) where {S} + return gradedrange([trivial(S) => n]) +end """ gradedrange(xs::AbstractVector{<:Pair}) diff --git a/src/matrixalgebrakit.jl b/src/matrixalgebrakit.jl index d5edd4d..0c4e03d 100644 --- a/src/matrixalgebrakit.jl +++ b/src/matrixalgebrakit.jl @@ -12,6 +12,7 @@ for f in [ :lq_compact, :lq_full, :lq_null, :eig_full, :eig_vals, :eigh_full, :eigh_vals, :left_polar, :right_polar, + :project_hermitian, :project_antihermitian, :project_isometric, ] f! = Symbol(f, :!) @eval function MAK.default_algorithm( @@ -283,6 +284,32 @@ function MAK.initialize_output( return P, Wᴴ end +# Projections +# ----------- +# Same output conventions as the generic implementations: hermitian and +# antihermitian project in place, isometric writes to a fresh output. +function MAK.initialize_output( + ::typeof(MAK.project_hermitian!), + A::FusedGradedMatrix, + alg::GradedBlockAlgorithm + ) + return A +end +function MAK.initialize_output( + ::typeof(MAK.project_antihermitian!), + A::FusedGradedMatrix, + alg::GradedBlockAlgorithm + ) + return A +end +function MAK.initialize_output( + ::typeof(MAK.project_isometric!), + A::FusedGradedMatrix, + alg::GradedBlockAlgorithm + ) + return similar(A) +end + # Truncation support # ------------------ @@ -303,6 +330,26 @@ function MAK.diagonal(v::FusedGradedVector) return FusedGradedMatrix(v.axis, v.axis, diag_blocks) end +# `pow_diag_safe!` for a block-diagonal graded matrix: clamp-power each reduced diagonal +# block. Only the reduced (degeneracy) data is touched, and that is correct even in the +# non-abelian case: a diagonal factor is `Diagonal(λ) ⊗ I` per sector, and `f(A ⊗ I) = +# f(A) ⊗ I`, so the power passes straight to the reduced eigenvalues. This is why the +# diagonal power is well defined here whereas a general element-wise `map!` on a graded +# array is not. +function TensorAlgebra.MatrixAlgebra.pow_diag_safe!( + Dp::FusedGradedMatrix, D::FusedGradedMatrix, p, tol + ) + foreachblock(MAK.diagview(Dp), MAK.diagview(D)) do _, (σp, σ) + return map!(d -> _clamped_pow(d, p, tol), σp, σ) + end + return Dp +end + +# Vendored from `TensorAlgebra.MatrixAlgebra` (not part of its public API): clamp entries +# below `tol` to zero, then raise to `p`; a negative entry above `tol` lets `real(d)^p` +# error for fractional `p`, enforcing the PSD precondition per-power. +_clamped_pow(d, p, tol) = abs(d) < tol ? zero(d) : real(d)^p + # Count how many elements are kept for a given index specification and block size _count_kept(::Colon, n) = n _count_kept(ind::AbstractVector{Bool}, _) = count(ind) diff --git a/test/Project.toml b/test/Project.toml index 1038a3d..d88ce40 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -33,7 +33,7 @@ SUNRepresentations = "0.3, 0.4" SafeTestsets = "0.1" SparseArraysBase = "0.10" Suppressor = "0.2.8" -TensorAlgebra = "0.16" +TensorAlgebra = "0.17" TensorKit = "0.17" TensorKitSectors = "0.3" Test = "1.10" diff --git a/test/test_abelianarray.jl b/test/test_abelianarray.jl index acb521a..8c01b2a 100644 --- a/test/test_abelianarray.jl +++ b/test/test_abelianarray.jl @@ -1,9 +1,9 @@ -using BlockArrays: BlockArrays, Block, blocklength +using BlockArrays: BlockArrays, Block, blocklength, blocklengths using Dictionaries: Dictionary using GradedArrays: GradedArrays, AbelianGradedArray, AbelianSectorArray, - AbstractGradedArray, FusedGradedMatrix, GradedOneTo, SU2, SectorRange, U1, data, - datalengths, dual, eachblockstoredindex, gradedrange, isdual, sectoraxes, sectors, - sectortype + AbstractGradedArray, FusedGradedMatrix, GradedOneTo, SU2, SectorRange, U1, + blockstoredlength, data, datalengths, dual, eachblockstoredindex, gradedrange, isdual, + sectoraxes, sectors, sectortype using LinearAlgebra: LinearAlgebra using Random: Random using TensorAlgebra: TensorAlgebra @@ -527,7 +527,7 @@ end @test !LinearAlgebra.isdiag(b) end -@testset "projectto! / checked_projectto!" begin +@testset "projectto! / project" begin g = gradedrange([U1(0) => 2, U1(1) => 3]) dest = AbelianGradedArray{Float64}(undef, g, dual(g)) @@ -541,17 +541,140 @@ end @test iszero(proj[1:2, 3:5]) @test iszero(proj[3:5, 1:2]) - # `checked_projectto!` accepts a source already in the allowed subspace. + # The checked `project` accepts a source already in the allowed subspace... src_allowed = Array(dest) - dest_ok = AbelianGradedArray{Float64}(undef, g, dual(g)) - @test TensorAlgebra.checked_projectto!(dest_ok, src_allowed) === dest_ok + dest_ok = TensorAlgebra.project(src_allowed, (g,), (g,)) + @test dest_ok isa AbelianGradedArray @test Array(dest_ok) ≈ src_allowed - # ... and rejects one carrying significant forbidden-block weight. + # ... and rejects one carrying significant forbidden-block weight, which the + # unchecked projection drops silently. src_bad = copy(src_allowed) src_bad[1, 5] += 10.0 - dest_bad = AbelianGradedArray{Float64}(undef, g, dual(g)) - @test_throws InexactError TensorAlgebra.checked_projectto!(dest_bad, src_bad) + @test_throws InexactError TensorAlgebra.project(src_bad, (g,), (g,)) + @test Array(TensorAlgebra.unchecked_project(src_bad, (g,), (g,))) ≈ src_allowed + + # A lower-rank `src` may omit trailing length-1 axes; `projectto!` reshapes it. This is + # the auxiliary flux-canceling leg idiom: an extra length-1 axis carries the charge that + # keeps an otherwise-forbidden component (here the `U1(1)` sector) in the allowed subspace. + site = gradedrange([U1(0) => 1, U1(1) => 1]) + aux = gradedrange([U1(0) => 1]) + state = AbelianGradedArray{Float64}(undef, site, aux) + @test TensorAlgebra.projectto!(state, [1.0, 0.0]) === state + @test size(state) == (2, 1) + @test Array(state) == reshape([1.0, 0.0], 2, 1) + @test TensorAlgebra.project([1.0, 0.0], (site, aux)) isa AbelianGradedArray +end + +@testset "project with a derived auxiliary leg" begin + # spin-1/2 site: up = U1(1), down = U1(-1) + g = gradedrange([U1(1) => 1, U1(-1) => 1]) + Splus = [0.0 1.0; 0.0 0.0] # |down> -> |up>, flux +2 + Sminus = [0.0 0.0; 1.0 0.0] # |up> -> |down>, flux -2 + Sz = [0.5 0.0; 0.0 -0.5] # neutral + + # The internal abelian fast path reads the net flux from the dominant entry. + @test GradedArrays.projected_charge(Splus, (g,), (g,)) == U1(2) + @test GradedArrays.projected_charge(Sminus, (g,), (g,)) == U1(-2) + @test GradedArrays.projected_charge(Sz, (g,), (g,)) == U1(0) + + # With all axes given (`src` rank matches the physical axes), `project` is the plain + # projection into the allowed subspace. + @test TensorAlgebra.project(Sz, (g,), (g,)) isa AbelianGradedArray + + # A trailing surplus axis in `src` (here the length-1 last axis of `(2,2,1)`, with only 2 + # physical axes given) is the aux: `project` derives its space into a flux-canceling last + # domain axis, giving a valid zero-total-flux tensor whose squeezed data is the original + # operator. The result's shape matches `src`'s. + t = TensorAlgebra.project(reshape(Splus, 2, 2, 1), (g,), (g,)) + @test t isa AbelianGradedArray + @test size(t) == (2, 2, 1) + @test blockstoredlength(t) == 1 # only the allowed block survives + @test Array(t)[:, :, 1] == Splus + + # A neutral operator still gets an aux, but a trivial one (dummy bond). + tz = TensorAlgebra.project(reshape(Sz, 2, 2, 1), (g,), (g,)) + @test size(tz) == (2, 2, 1) + @test Array(tz)[:, :, 1] == Sz + + # A multi-slice aux derives a direct sum, one charge per slice — the MPO-virtual-leg + # structure: each slice's flux is canceled by its own aux sector, so the whole tensor is + # invariant even though the slices carry different charges. + stack = cat(reshape(Splus, 2, 2, 1), reshape(Sminus, 2, 2, 1); dims = 3) + ts = TensorAlgebra.project(stack, (g,), (g,)) + @test ts isa AbelianGradedArray + @test size(ts) == (2, 2, 2) + @test Array(ts)[:, :, 1] == Splus + @test Array(ts)[:, :, 2] == Sminus + + # ... including a neutral slice mixed with a charged one. + tsz = TensorAlgebra.project( + cat(reshape(Sz, 2, 2, 1), reshape(Splus, 2, 2, 1); dims = 3), (g,), (g,) + ) + @test Array(tsz)[:, :, 1] == Sz + @test Array(tsz)[:, :, 2] == Splus + + # Contiguous equal charges merge into one sector of that multiplicity (matching the + # `TensorMap` backend's `s => m` behavior), and agree with passing the merged aux + # explicitly through the all-given form. + pp = cat(reshape.((Splus, Splus), 2, 2, 1)...; dims = 3) + tpp = TensorAlgebra.project(pp, (g,), (g,)) + @test blocklengths(axes(tpp, 3)) == [2] + @test Array(tpp) == pp + tmerged = TensorAlgebra.project(pp, (g,), (g, gradedrange([U1(2) => 2]))) + @test axes(tmerged) == axes(tpp) + @test Array(tmerged) == pp + + # Non-contiguous repeats stay separate sectors (a `GradedOneTo` permits unmerged + # duplicates), so slice order is always preserved. + repeats = cat(reshape.((Splus, Sz, Splus), 2, 2, 1)...; dims = 3) + trep = TensorAlgebra.project(repeats, (g,), (g,)) + @test size(trep) == (2, 2, 3) + @test blocklengths(axes(trep, 3)) == [1, 1, 1] + @test Array(trep) == repeats + + # Only one trailing surplus axis is supported: more is an error, not a silent flattening. + @test_throws ArgumentError TensorAlgebra.project( + reshape(stack, 2, 2, 2, 1), (g,), (g,) + ) + + # A lower-rank `src` that omits explicitly-given trailing length-1 axes is the trailing-axes + # tolerance, not a surplus axis: it pads, it does not derive an extra leg. The domain aux is + # given codomain-facing (stored dualized), so `[U1(1)]` cancels the charge-1 component. + site = gradedrange([U1(0) => 1, U1(1) => 1]) + saux = gradedrange([U1(1) => 1]) + po = TensorAlgebra.project([0.0, 1.0], (site,), (saux,)) + @test size(po) == (2, 1) + @test Array(po) == reshape([0.0, 1.0], 2, 1) + + # The flat all-codomain (state) form also derives: a stack of basis states with different + # charges gets a multi-sector aux. + ps = TensorAlgebra.project([1.0 0.0; 0.0 1.0], (site,)) + @test size(ps) == (2, 2) + @test Array(ps) == [1.0 0.0; 0.0 1.0] + + # `project` verifies nothing was discarded after the derivation: a slice with weight + # outside its derived (dominant-entry) charge is rejected, where `unchecked_project` + # silently drops it. + junk = copy(reshape(Splus, 2, 2, 1)) + junk[2, 2, 1] = 0.3 + @test_throws InexactError TensorAlgebra.project(junk, (g,), (g,)) + @test Array(TensorAlgebra.unchecked_project(junk, (g,), (g,)))[:, :, 1] == Splus + + # `tryproject` is the nullable sibling of `project`: branch on whether the data is + # invariant in the given axes, falling back to deriving the flux-carrying aux. + v_inv, v_chg = [1.0, 0.0], [0.0, 1.0] + @test TensorAlgebra.tryproject(v_inv, (site,)) isa AbelianGradedArray + @test isnothing(TensorAlgebra.tryproject(v_chg, (site,))) + t_chg = @something TensorAlgebra.tryproject(v_chg, (site,)) TensorAlgebra.project( + reshape(v_chg, 2, 1), (site,) + ) + @test ndims(t_chg) == 2 + @test Array(t_chg) == reshape(v_chg, 2, 1) + t_inv = @something TensorAlgebra.tryproject(v_inv, (site,)) TensorAlgebra.project( + reshape(v_inv, 2, 1), (site,) + ) + @test ndims(t_inv) == 1 end @testset "getindex (project dense onto graded axes)" begin