Skip to content

Show improvements + subblocks iterator #304

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 21 commits into
base: main
Choose a base branch
from
+389 −184
Open

Show improvements + subblocks iterator #304

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
21 commits
Select commit Hold shift + click to select a range
8ca49b8
rework gradedspace show
lkdvos Oct 14, 2025
a6e091a
rework tensor show
lkdvos Oct 14, 2025
f48e980
add convenience functions
lkdvos Oct 17, 2025
a2b595f
rework tensor show again
lkdvos Oct 17, 2025
9710451
rework adjoint tensor show
lkdvos Oct 17, 2025
f121193
update blockshow
lkdvos Oct 17, 2025
070dab2
add `SubblockIterator`
lkdvos Oct 17, 2025
53f0b61
rework diagonal tensor show
lkdvos Oct 17, 2025
fe20385
some indexing reworking
lkdvos Oct 17, 2025
1586be7
adjoint indexing
lkdvos Oct 17, 2025
5fcd778
diagonal indexing
lkdvos Oct 17, 2025
a7a6af4
remove double export
lkdvos Oct 17, 2025
cdd6649
small doc changes
lkdvos Oct 17, 2025
077de13
remove vertical gradedspace printing
lkdvos Oct 17, 2025
4bbdb93
retry gradedspace printing
lkdvos Oct 17, 2025
1c156fe
improve errors
lkdvos Oct 17, 2025
53a5dc9
remove ambiguous function
lkdvos Oct 17, 2025
e1ef363
simplify dual
lkdvos Oct 23, 2025
9266ed9
use reduced dim
lkdvos Oct 23, 2025
0285d36
address pr comments
lkdvos Oct 24, 2025
4c77cb3
code review
lkdvos Oct 25, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions Project.toml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
MatrixAlgebraKit = "6c742aac-3347-4629-af66-fc926824e5e4"
OhMyThreads = "67456a42-1dca-4109-a031-0a68de7e3ad5"
PackageExtensionCompat = "65ce6f38-6b18-4e1d-a461-8949797d7930"
Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
ScopedValues = "7e506255-f358-4e82-b7e4-beb19740aa63"
Strided = "5e0ebb24-38b0-5f93-81fe-25c709ecae67"
Expand Down Expand Up @@ -37,6 +38,7 @@ LinearAlgebra = "1"
MatrixAlgebraKit = "0.5.0"
OhMyThreads = "0.8.0"
PackageExtensionCompat = "1"
Printf = "1"
Random = "1"
SafeTestsets = "0.1"
ScopedValues = "1.3.0"
Expand Down
6 changes: 3 additions & 3 deletions docs/src/lib/tensors.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -118,15 +118,15 @@ blocks

To access the data associated with a specific fusion tree pair, you can use:
```@docs
Base.getindex(::TensorMap{T,S,N₁,N₂}, ::FusionTree{I,N₁}, ::FusionTree{I,N₂}) where {T,S,N₁,N₂,I<:Sector}
Base.setindex!(::TensorMap{T,S,N₁,N₂}, ::Any, ::FusionTree{I,N₁}, ::FusionTree{I,N₂}) where {T,S,N₁,N₂,I<:Sector}
Base.getindex(::AbstractTensorMap, ::FusionTree, ::FusionTree)
Base.setindex!(::AbstractTensorMap, ::Any, ::FusionTree, ::FusionTree)
```

For a tensor `t` with `FusionType(sectortype(t)) isa UniqueFusion`, fusion trees are
completely determined by the outcoming sectors, and the data can be accessed in a more
straightforward way:
```@docs
Base.getindex(::TensorMap, ::Tuple{I,Vararg{I}}) where {I<:Sector}
Base.getindex(::AbstractTensorMap, ::Tuple{I,Vararg{I}}) where {I<:Sector}
```

For tensor `t` with `sectortype(t) == Trivial`, the data can be accessed and manipulated
Expand Down
3 changes: 2 additions & 1 deletion src/TensorKit.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -57,7 +57,7 @@ export ℤ₂Space, ℤ₃Space, ℤ₄Space, U₁Space, CU₁Space, SU₂Space
# Export tensor map methods
export domain, codomain, numind, numout, numin, domainind, codomainind, allind
export spacetype, storagetype, scalartype, tensormaptype
export blocksectors, blockdim, block, blocks
export blocksectors, blockdim, block, blocks, subblocks, subblock

# random methods for constructor
export randisometry, randisometry!, rand, rand!, randn, randn!
Expand Down Expand Up @@ -127,6 +127,7 @@ using Base: @boundscheck, @propagate_inbounds, @constprop,
tuple_type_head, tuple_type_tail, tuple_type_cons,
SizeUnknown, HasLength, HasShape, IsInfinite, EltypeUnknown, HasEltype
using Base.Iterators: product, filter
using Printf: @sprintf

using LinearAlgebra: LinearAlgebra, BlasFloat
using LinearAlgebra: norm, dot, normalize, normalize!, tr,
Expand Down
82 changes: 68 additions & 14 deletions src/spaces/gradedspace.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -197,21 +197,75 @@ function supremum(V₁::GradedSpace{I}, V₂::GradedSpace{I}) where {I <: Sector
)
end

function Base.show(io::IO, V::GradedSpace{I}) where {I <: Sector}
print(io, type_repr(typeof(V)), "(")
separator = ""
comma = ", "
io2 = IOContext(io, :typeinfo => I)
for c in sectors(V)
if isdual(V)
print(io2, separator, dual(c), "=>", dim(V, c))
else
print(io2, separator, c, "=>", dim(V, c))
end
separator = comma
Base.summary(io::IO, V::GradedSpace) = print(io, type_repr(typeof(V)))

function Base.show(io::IO, V::GradedSpace)
pre = (get(io, :typeinfo, Any)::DataType == typeof(V) ? "" : type_repr(typeof(V)))

io = IOContext(io, :typeinfo => Pair{sectortype(V), Int})

pre *= "("
if isdual(V)
post = ")'"
V = dual(V)
else
post = ")"
end
hdots = " \u2026 "
sep = ", "
sepsize = length(sep)

limited = get(io, :limit, false)::Bool
screenwidth = limited ? displaysize(io)[2] : typemax(Int)
screenwidth -= length(pre) + length(post)

cs = reshape(collect([c => dim(V, c) for c in sectors(V)]), 1, :)
ncols = length(cs)

maxpossiblecols = screenwidth ÷ (1 + sepsize)
if ncols > maxpossiblecols
cols = [1:(maxpossiblecols - 1); (ncols - maxpossiblecols + 1):ncols]
else
cols = collect(1:ncols)
end

A = Base.alignment(io, cs, [1], cols, screenwidth, screenwidth, sepsize, ncols)

if ncols <= length(A) # fits on screen
print(io, pre)
Base.print_matrix_row(io, cs, A, 1, cols, sep, ncols)
print(io, post)
else # doesn't fit on screen
half = (screenwidth - length(hdots) + 1) ÷ 2 + 1
Ralign = reverse(Base.alignment(io, cs, [1], reverse(cols), half, half, sepsize, ncols))
half = screenwidth - sum(map(sum, Ralign)) - (length(Ralign) - 1) * sepsize - length(hdots)
Lalign = Base.alignment(io, cs, [1], cols, half, half, sepsize, ncols)
print(io, pre)
Base.print_matrix_row(io, cs, Lalign, 1, cols[1:length(Lalign)], sep, ncols)
print(io, hdots)
Base.print_matrix_row(io, cs, Ralign, 1, (length(cs) - length(Ralign)) .+ cols, sep, length(cs))
print(io, post)
end
print(io, ")")
V.dual && print(io, "'")

return nothing
end

function Base.show(io::IO, ::MIME"text/plain", V::GradedSpace)
# print small summary, e.g.
# d-element Vect[I] or d-element dual(Vect[I])
d = reduceddim(V)
print(io, d, "-element ")
isdual(V) && print(io, "dual(")
print(io, type_repr(typeof(V)))
isdual(V) && print(io, ")")
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I still think this d-element printing style is not really appropriate. When is this method called; is this in notebooks?

How about some printing style for an object like V = SU2Space(0=>3, 1=>2, 2=>1):

"Rep[SU₂] object with dimension 14 (reduced dimension 6)"

or for V'

"Rep[SU₂] dual object with dimension 14 (reduced dimension 6)"

where also the extra information in parentheses can be omitted if dim(V) == reduceddim(V).

Just a suggestion, let me know what you think.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Doing another deep dive into the implementations and documentation, this is what I can find:

The different show methods are mostly used like this:

  • show(io, MIME"text/plain"(), x) is what gets called mostly in interactive usage, i.e. it ends up being the implementation when you display something to an output that can support text. (I.e. REPL, or notebook output, ...). This is supposed to be for human consumption
  • Container types generally implement 3-argument show by calling show(io, MIME"text/plain"(), x) for elements x, with :compact => true set. In this case you should avoid \n to not break the printing of the container.
  • :limit => true can be set to print … in place of most elements in a container.
  • show(io, x) is the default implementation/fallback, which is typically used for computer consumption (i.e. parseable)

For example Array handles this for the 3-arg version by writing a summary first, and then, on new line(s) prints the data. If compact => true, the summary is omitted and the [x; y;; ...] single-line notation is used to show the data.

So, what I would like is:

  • The 3-arg show for a space gives me a quick summary with the most important data, which is, in my opinion, the sectortype, some measure of size and the duality. Then, the sectors and dimensions can be used afterwards to fully specify the space
  • The 3-arg show with :compact => true should not contain any \n. We can decide to either show Vect[I](x => d, ...) or show the summary, depending on what we think is the most useful information to show at a glance in a container of spaces. This is the thing that would get called e.g. when showing the spaces of an MPS, which you could think of as a container space

I don't actually feel to strongly about how this information is worded exactly, but I do think showing the dimension somewhere is a strict improvement, for example in MPS or PEPS methods, the relevant spaces to inspect are the virtual spaces, which tend to become completely unreadable when there are many sectors. In this case, simply truncating the sectors and going for the … reduces the screen clutter, but then there is no more information on the screen.

I have no real preference about the dim vs reduced dim, I argue that the information is mostly in comparing different values anyways (I want to see if my dimension increased or decreased during a DMRG run), but given that we are typically not using reduceddim in many places I feel like to not make this number confusing we would have to be a little explicit about it.

"Rep[SU₂] dual object with dimension 14 (reduced dimension 6)" seems a bit long, I think one of the dimensions is sufficient? I would also prefer to avoid the categoric "object" language when this isn't necessary, so maybe to suggest either alternative:

"(dual) Rep[SU₂] of dim d"
"(dual) Rep[SU₂] of reduceddim d"

This gives all information, would be Float dim-friendly since the floating point number is at the end, avoids (more) use of category language and is reasonably short. dim and reduceddim are idiomatic since we are already using that as a function name, so it cannot be confusing what this number is referring to.

Copy link
Member

@Jutho Jutho Oct 25, 2025 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The use of 'object' was supposed to be such that only category-minded people (like you) think about categorical objects 😄 , and other people think of a computer-science object of type Rep[SU₂]. Though the possible occurance of dual then spoils that a little bit. I just think the one guideline should be not to take Array as an example, as really spaces are supposed to be like small bits types (even if GradedSpace isn't really. Furthermore, Array still prints out all elements even with :compact => true, so by that guideline we can still print all the sectors. But I agree that this is not what you want to see when an MPS is printed. However, I do think it is strange if ℂ^4 is printed just as this, and Rep[SU₂](...) is printed as dual Rep[SU₂] of dim d. I kind of dislike the implicit implication that it is an instance and not a type.

How about the following compromise: With :compact => true, a V::GradedSpace is printed as

type_repr(typeof(V)) * "(…)" * (isdual(V) ? "'" : "") * " with dim $(dim(V))"

Depending on the length of the information that you want to show, you can add additional information such as " with $(length(sectors(V)) sectors and dim $(dim(V))" but this is really up to you.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am still confused by the show methods. If you just have an interactive REPL session, it goes via display(V) and ends up calling show(stdout, MIME"text/plain", V). Is that than also the version that should check for :limit? So we might have :compact=>false but :limit=>true, so we would have a second line printing the sectors? In fact, for show(stdout, MIME"text/plain", V), I wouldn't mind if the sectors are printed underneath each other as in your original design.

The version without MIME type, show(io, V), is only called if you call explicitly show(V). That one, even for Array, prints out all information and is supposed to be parseable. So can that just remain the current version that prints out all the sectors? For Array, the only difference between :compact true or false in that case is that it truncates floating point numbers to less accuracy, but that doesn't really make sense for sector dimensions.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Oh wait, show(io, V) does also listen to :limit. Pff, this is annoying.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think (this is my interpretation, it's all very badly defined/documented) that :compact means no linebreaks and optionally more dense representations (less floating point digits, or x=>y instead of x => y), and :limit means respecting the displaysize.
The version without MIME type is only called if you explicitly show(V) or print(io, V), where I think it is supposed to be parseable but can still be :limitted if you like


compact = get(io, :compact, false)::Bool
(iszero(d) || compact) && return nothing

# print detailed sector information
print(io, ":\n ")
ioc = IOContext(io, :typeinfo => typeof(V))
show(ioc, V)
return nothing
end

Expand Down
184 changes: 184 additions & 0 deletions src/tensors/abstracttensor.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -250,6 +250,12 @@ Return an iterator over all splitting - fusion tree pairs of a tensor.
"""
fusiontrees(t::AbstractTensorMap) = fusionblockstructure(t).fusiontreelist

fusiontreetype(t::AbstractTensorMap) = fusiontreetype(typeof(t))
function fusiontreetype(::Type{T}) where {T <: AbstractTensorMap}
I = sectortype(T)
return Tuple{fusiontreetype(I, numout(T)), fusiontreetype(I, numin(T))}
end

# auxiliary function
@inline function trivial_fusiontree(t::AbstractTensorMap)
sectortype(t) === Trivial ||
Expand Down Expand Up @@ -295,6 +301,145 @@ function blocktype(::Type{T}) where {T <: AbstractTensorMap}
return Core.Compiler.return_type(block, Tuple{T, sectortype(T)})
end

# tensor data: subblock access
# ----------------------------
@doc """
subblocks(t::AbstractTensorMap)

Return an iterator over all subblocks of a tensor, i.e. all fusiontrees and their
corresponding tensor subblocks.

See also [`subblock`](@ref), [`fusiontrees`](@ref), and [`hassubblock`](@ref).
"""
subblocks(t::AbstractTensorMap) = SubblockIterator(t, fusiontrees(t))

const _doc_subblock = """
Return a view into the data of `t` corresponding to the splitting - fusion tree pair
`(f₁, f₂)`. In particular, this is an `AbstractArray{T}` with `T = scalartype(t)`, of size
`(dims(codomain(t), f₁.uncoupled)..., dims(codomain(t), f₂.uncoupled)...)`.

Whenever `FusionStyle(sectortype(t)) isa UniqueFusion` , it is also possible to provide only
the external `sectors`, in which case the fusion tree pair will be constructed automatically.
"""

@doc """
subblock(t::AbstractTensorMap, (f₁, f₂)::Tuple{FusionTree,FusionTree})
subblock(t::AbstractTensorMap, sectors::Tuple{Vararg{Sector}})

$_doc_subblock

In general, new tensor types should provide an implementation of this function for the
fusion tree signature.

See also [`subblocks`](@ref) and [`fusiontrees`](@ref).
""" subblock

Base.@propagate_inbounds function subblock(t::AbstractTensorMap, sectors::Tuple{I, Vararg{I}}) where {I <: Sector}
# input checking
I === sectortype(t) || throw(SectorMismatch("Not a valid sectortype for this tensor."))
FusionStyle(I) isa UniqueFusion ||
throw(SectorMismatch("Indexing with sectors is only possible for unique fusion styles."))
length(sectors) == numind(t) || throw(ArgumentError("invalid number of sectors"))

# convert to fusiontrees
s₁ = TupleTools.getindices(sectors, codomainind(t))
s₂ = map(dual, TupleTools.getindices(sectors, domainind(t)))
c1 = length(s₁) == 0 ? unit(I) : (length(s₁) == 1 ? s₁[1] : first(⊗(s₁...)))
@boundscheck begin
hassector(codomain(t), s₁) && hassector(domain(t), s₂) || throw(BoundsError(t, sectors))
c2 = length(s₂) == 0 ? unit(I) : (length(s₂) == 1 ? s₂[1] : first(⊗(s₂...)))
c2 == c1 || throw(SectorMismatch("Not a valid fusion channel for this tensor"))
end
f₁ = FusionTree(s₁, c1, map(isdual, tuple(codomain(t)...)))
f₂ = FusionTree(s₂, c1, map(isdual, tuple(domain(t)...)))
return @inbounds subblock(t, (f₁, f₂))
end
Base.@propagate_inbounds function subblock(t::AbstractTensorMap, sectors::Tuple)
return subblock(t, map(Base.Fix1(convert, sectortype(t)), sectors))
end
# attempt to provide better error messages
function subblock(t::AbstractTensorMap, (f₁, f₂)::Tuple{FusionTree, FusionTree})
(sectortype(t)) == sectortype(f₁) == sectortype(f₂) ||
throw(SectorMismatch("Not a valid sectortype for this tensor."))
numout(t) == length(f₁) && numin(t) == length(f₂) ||
throw(DimensionMismatch("Invalid number of fusiontree legs for this tensor."))
throw(MethodError(subblock, (t, (f₁, f₂))))
end

@doc """
subblocktype(t)
subblocktype(::Type{T})

Return the type of the tensor subblocks of a tensor.
""" subblocktype

function subblocktype(::Type{T}) where {T <: AbstractTensorMap}
return Core.Compiler.return_type(subblock, Tuple{T, fusiontreetype(T)})
end
subblocktype(t) = subblocktype(typeof(t))
subblocktype(T::Type) = throw(MethodError(subblocktype, (T,)))

# Indexing behavior
# -----------------
@doc """
Base.view(t::AbstractTensorMap, sectors::Tuple{Vararg{Sector}})
Base.view(t::AbstractTensorMap, f₁::FusionTree, f₂::FusionTree)

$_doc_subblock

!!! note
Contrary to Julia's array types, the default indexing behavior is to return a view
into the tensor data. As a result, `getindex` and `view` have the same behavior.

See also [`subblock`](@ref), [`subblocks`](@ref) and [`fusiontrees`](@ref).
""" Base.view(::AbstractTensorMap, ::Tuple{I, Vararg{I}}) where {I <: Sector},
Base.view(::AbstractTensorMap, ::FusionTree, ::FusionTree)

@inline Base.view(t::AbstractTensorMap, sectors::Tuple{I, Vararg{I}}) where {I <: Sector} =
subblock(t, sectors)
@inline Base.view(t::AbstractTensorMap, f₁::FusionTree, f₂::FusionTree) =
subblock(t, (f₁, f₂))
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do you think it is useful to have view methods? I kind of consider this to be AbstractArray specific. I am a bit hesitant to provide to many alternative interfaces.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think if I would redesign everything right now, I would have gone with view instead of getindex to begin with, since we are technically performing "method-piracy" on Base.getindex, which is supposed to yield a copy.
I don't want to change this now (that would be really painful), but I do find it easier to say "we have made the decision to make Base.getindex(::AbstractTensorMap) default to behave like Base.view(::AbstractTensorMap)", since these two concepts (as words) already exist and have meaning to me.
This is more or less the only reason why I added this, I don't really see many people using it, but I also don't think it hurts to have this here.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

First an aside, your remark about the "method piracy" on Base.getindex is probably valid for StridedView, though in that case I do think it is well documented that slicing it will produce a view instead of a new copy.

For AbstractTensorMap, I do not agree, because I do not consider Base.getindex as a method that is specific to AbstractArray, and therefore I do not think it should follow the same rules. Many types have a getindex that have nothing to do with AbstractArray. If you would have some kind of dictionary that stores blocks as values with fusion tree pairs as keys, it would be very natural that indexing with the key would return the original block, without making a copy. I think this is the analogy I have in my head when I think of (Abstract)TensorMap, so I really don't think of them as AbstractArrays. It seems ITensors have the same behavior: T[Block((1,3))] returns an AbstractArray subtype that is actually providing a view into the data.

In fact, I believe returning a view is necessary if we want some syntax like T[(f1, f2)] .= 1 to work. Here, it becomes messy with Array, since A[1:2, 1:3] .= value will also overwrite the data, but I really have no idea what this is lowered to. I assume something like setindex!.(Ref(A), value, ...) where I don't know what ... needs to be so that it is not considered as just a one-dimensional broadcast. But this anyways cannot work for our case.

So while I don't think our TensorMap getindex is problematic, I agree that using a separate method like subblock makes the distinction even more clear.

However, also recycling view for this is again hitting the fundamental disagreement we seem to have. I don't like to think of (Abstract)TensorMap as an AbstractArray subtype, and so I had opted not to use any methods that seem to be specific to AbstractArray in Julia base (ndims, size, axes, view, ...) for AbstractTensorMap. However, AbstractTensorMap definitely is a container type, and so the basic methods for general containers (getindex and setindex!) are available to us, with the key types and behavior that seems appropriate for our use case.

It seems ITensors are a bit in a mixed category: they listen to size and ndims, and you can index them with integer indices, but not with ranges and you cannot apply view. Do you have any insight into the user experience and typical user requests in that regard?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think most of the complaints surrounding ITensors seem to land in the same category: the lack of support for slicing is the most-often quoted nuisance, and the opaque manner in which the data is stored when introducing symmetries is not great.
Overall, my feeling is that people tend to expect these things to "just work", and they care less about the performance than they do about the convenience of not having to alter their habits of how to code or interact with arrays. I would say that people are receptive to the "this works, but if you want a fast/better implementation you have to work a little for it"


# by default getindex returns views
@doc """
Base.getindex(t::AbstractTensorMap, sectors::Tuple{Vararg{Sector}})
t[sectors]
Base.getindex(t::AbstractTensorMap, f₁::FusionTree, f₂::FusionTree)
t[f₁, f₂]

$_doc_subblock

!!! warning
Contrary to Julia's array types, the default behavior is to return a view into the tensor data.
As a result, modifying the view will modify the data in the tensor.

See also [`subblock`](@ref), [`subblocks`](@ref) and [`fusiontrees`](@ref).
""" Base.getindex(::AbstractTensorMap, ::Tuple{I, Vararg{I}}) where {I <: Sector},
Base.getindex(::AbstractTensorMap, ::FusionTree, ::FusionTree)

@inline Base.getindex(t::AbstractTensorMap, sectors::Tuple{I, Vararg{I}}) where {I <: Sector} =
view(t, sectors)
@inline Base.getindex(t::AbstractTensorMap, f₁::FusionTree, f₂::FusionTree) =
view(t, f₁, f₂)

@doc """
Base.setindex!(t::AbstractTensorMap, v, sectors::Tuple{Vararg{Sector}})
t[sectors] = v
Base.setindex!(t::AbstractTensorMap, v, f₁::FusionTree, f₂::FusionTree)
t[f₁, f₂] = v

Copies `v` into the data slice of `t` corresponding to the splitting - fusion tree pair `(f₁, f₂)`.
By default, `v` can be any object that can be copied into the view associated with `t[f₁, f₂]`.

See also [`subblock`](@ref), [`subblocks`](@ref) and [`fusiontrees`](@ref).
""" Base.setindex!(::AbstractTensorMap, ::Any, ::Tuple{I, Vararg{I}}) where {I <: Sector},
Base.setindex!(::AbstractTensorMap, ::Any, ::FusionTree, ::FusionTree)

@inline Base.setindex!(t::AbstractTensorMap, v, sectors::Tuple{I, Vararg{I}}) where {I <: Sector} =
copy!(view(t, sectors), v)
@inline Base.setindex!(t::AbstractTensorMap, v, f₁::FusionTree, f₂::FusionTree) =
copy!(view(t, (f₁, f₂)), v)

# Derived indexing behavior for tensors with trivial symmetry
#-------------------------------------------------------------
using TensorKit.Strided: SliceIndex
Expand Down Expand Up @@ -480,6 +625,10 @@ end

# Conversion to Array:
#----------------------
Base.ndims(t::AbstractTensorMap) = numind(t)
Base.size(t::AbstractTensorMap) = ntuple(Base.Fix1(size, t), numind(t))
Base.size(t::AbstractTensorMap, i) = dim(space(t, i))
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am still not convinced about the usefulness of this. This is mostly aimed at making it more friendly for newcomers? But shouldn't we stimulate them to just learn from the start that TensorMaps are not just AbstractArrays and that they should just check the documentation to learn about the interface. I am happy to for example have a summary of the interface in the README or whathever can help.

One reason why I think this might not be a good idea is that this forces length(t) = prod(size(t)) = dim(codomain) * dim(domain). But there might be an argument to be made for length(t) = dim(codomain ← domain) = dim(t). So once we establish this is part of our public interface, it will be very hard to revert this if we ever come up with a better use case for these definitions.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I definitely agree with the sentiment to steer people towards our interface, but I don't think we currently have one for this exact concept, and (dims(codomain(t))..., dims(domain(t))...) is a bit less convenient in my opinion. Since I think it is quite obvious what size would return, it is a bit weird to me to get a MethodError. I don't think we as a library should be deciding if it is allowed to use size or not, and because of type piracy reasons this is basically the only place this method can be defined, so we might as well add it?

I have to say that I am not sure I agree that your reason is really related to this. I don't think length(t) = prod(size(t)) is expected, and I don't see why we should hold off on defining size just because there might be other functions that could have multiple meanings.

Ultimately, it's mostly just about being pragmatic. It's useful functionality, it's a short, easily available function name that is easy to remember and it's clear what it does, and I see very little harm in overloading it

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Would using dims(::HomSpace) be an adequate solution, or you really think users will want to use size? I have to say that this has never come up as an open issue on Github, but maybe users complain to you directly?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I have had a couple people mention that they found it surprising that it doesn't work (technically, I quote: "this is why no one is using TensorKit")
I'm okay with dims(::HomSpace), (can we also have dims(::AbstractTensorMap)?), but still have a hard time motivating "we didn't implement size because or tensor type would look to much like an array which we want to avoid as it might cause confusion". Maybe I'm just reading too much into it, but this seems like a very easy thing to fix and just provide?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sure, if dims(::HomSpace) work, than dims(::AbstractTensorMap) must return the same. But I guess I will have to give up my fight against size then. Maybe we need an allowsize like CUDA.allowscalar ? 😄 .

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I especially look forward to the warning on the first call:

julia >size(::TensorMap)
WARNING: a TensorMap is not an AbstractArray and has a different interface; expect non-integer dimensions and other unexpected behavior. RTFM!

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Jokes aside, I would actually be quite happy with a TensorKit.allowscalar to allow scalar getindex/setindex for abelian and getindex for non-abelian tensors too 😁


# probably not optimized for speed, only for checking purposes
function Base.convert(::Type{Array}, t::AbstractTensorMap)
I = sectortype(t)
Expand All @@ -499,3 +648,38 @@ function Base.convert(::Type{Array}, t::AbstractTensorMap)
return A
end
end

# Show and friends
# ----------------

function Base.dims2string(V::HomSpace)
str_cod = numout(V) == 0 ? "()" : join(dim.(codomain(V)), '×')
str_dom = numin(V) == 0 ? "()" : join(dim.(domain(V)), '×')
return str_cod * "←" * str_dom
end

function Base.summary(io::IO, t::AbstractTensorMap)
V = space(t)
print(io, Base.dims2string(V), " ")
Base.showarg(io, t, true)
return nothing
end

# Human-readable:
function Base.show(io::IO, ::MIME"text/plain", t::AbstractTensorMap)
# 1) show summary: typically d₁×d₂×… ← d₃×d₄×… $(typeof(t)):
summary(io, t)
println(io, ":")

# 2) show spaces
# println(io, " space(t):")
println(io, " codomain: ", codomain(t))
println(io, " domain: ", domain(t))

# 3) [optional]: show data
get(io, :compact, true) && return nothing
ioc = IOContext(io, :typeinfo => sectortype(t))
println(io, "\n\n blocks(t):")
show_blocks(io, MIME"text/plain"(), blocks(t))
return nothing
end
Loading
Loading