Support GradedArrays 0.17 and show an index's space

Project.toml

@@ -1,6 +1,6 @@
 name = "ITensorBase"
 uuid = "4795dd04-0d67-49bb-8f44-b89c448a1dc7"
-version = "0.15.3"
+version = "0.15.4"
 authors = ["ITensor developers <support@itensor.org> and contributors"]
 
 [workspace]
@@ -44,7 +44,7 @@ Adapt = "4.1.1"
 ArrayLayouts = "1.11"
 Combinatorics = "1"
 ConstructionBase = "1.6"
-GradedArrays = "0.16.6"
+GradedArrays = "0.17"
 LinearAlgebra = "1.10"
 MatrixAlgebraKit = "0.2, 0.3, 0.4, 0.5, 0.6"
 Mooncake = "0.4.202, 0.5"

README.md

@@ -7,6 +7,11 @@
 [![Code Style](https://img.shields.io/badge/code_style-ITensor-purple)](https://github.com/ITensor/ITensorFormatter.jl)
 [![Aqua](https://raw.githubusercontent.com/JuliaTesting/Aqua.jl/master/badge.svg)](https://github.com/JuliaTesting/Aqua.jl)
 
+A next-generation rewrite of [ITensors.jl](https://github.com/ITensor/ITensors.jl), allowing
+arbitrary array backends and a wider range of symmetries. Built on top of
+[TensorAlgebra.jl](https://github.com/ITensor/TensorAlgebra.jl), with group symmetric tensors
+powered by [GradedArrays.jl](https://github.com/ITensor/GradedArrays.jl).
+
 ## Support
 
 <picture>

docs/Project.toml

@@ -1,5 +1,7 @@
 [deps]
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+DocumenterInterLinks = "d12716ef-a0f6-4df4-a9f1-a5a34e75c656"
+GradedArrays = "bc96ca6e-b7c8-4bb6-888e-c93f838762c2"
 ITensorBase = "4795dd04-0d67-49bb-8f44-b89c448a1dc7"
 ITensorFormatter = "b6bf39f1-c9d3-4bad-aad8-593d802f65fd"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
@@ -12,6 +14,8 @@ path = ".."
 
 [compat]
 Documenter = "1"
+DocumenterInterLinks = "1"
+GradedArrays = "0.17"
 ITensorBase = "0.15"
 ITensorFormatter = "0.2.27"
 Literate = "2"

docs/make.jl

@@ -1,4 +1,5 @@
 using Documenter: Documenter, DocMeta, deploydocs, makedocs
+using DocumenterInterLinks: InterLinks
 using ITensorBase
 using ITensorFormatter: ITensorFormatter
 
@@ -10,6 +11,8 @@ DocMeta.setdocmeta!(ITensorBase, :DocTestSetup, :(using ITensorBase); recursive
 
 ITensorFormatter.make_index!(pkgdir(ITensorBase))
 
+links = InterLinks("GradedArrays" => "https://itensor.github.io/GradedArrays.jl/dev/")
+
 makedocs(;
     modules = [ITensorBase],
     authors = "ITensor developers <support@itensor.org> and contributors",
@@ -22,9 +25,11 @@ makedocs(;
     pages = [
         "Home" => "index.md",
         "User Interface" => "user_interface.md",
+        "Symmetric tensors" => "symmetries.md",
         "Developer Interface" => "dev_interface.md",
         "Reference" => "reference.md",
-    ]
+    ],
+    plugins = [links]
 )
 
 deploydocs(;

docs/src/symmetries.md

@@ -0,0 +1,67 @@
+# Symmetric tensors
+
+```@meta
+CurrentModule = ITensorBase
+```
+
+ITensorBase supports tensors that are symmetric under group actions by wrapping ITensors around
+[GradedArrays.jl](https://github.com/ITensor/GradedArrays.jl). To get started, build
+[`Index`](@ref) objects out of `sector => multiplicity` pairs and pass them to the standard Julia
+array constructors (`randn`, `zeros`, and so on):
+
+```@example symmetries
+using ITensorBase: Index, inds
+using GradedArrays: U1, dual, isdual
+
+i = Index([U1(0) => 1, U1(1) => 2])
+j = Index([U1(0) => 2, U1(1) => 1])
+k = Index([U1(0) => 1, U1(1) => 1])
+
+a = randn(i, dual(j))
+```
+
+A `GradedArray` only stores the symmetry-allowed blocks.
+
+These tensors support contraction, multiplication by a scalar, and addition.
+
+```@example symmetries
+b = randn(j, dual(k))
+a * b
+```
+
+```@example symmetries
+2 * a
+```
+
+```@example symmetries
+c = randn(i, dual(j))
+a + c
+```
+
+## Duality
+
+`dual` flips the duality of an index, and `isdual` returns whether an index is dual. Indices
+can only contract with ones that have opposite duality, for example the `j` Index of `b`
+contracts with the `dual(j)` Index of `a`.
+
+```@example symmetries
+isdual.(inds(a))
+```
+
+Note that indices of a `GradedArray` are partitioned into a codomain and a domain, and the
+`GradedArray` stores the block diagonal matrix corresponding to the bipartitioning of the
+indices. When printing, by convention domain indices are implicitly dual (the format and
+conventions are compatible with those from
+[TensorKit.jl](https://github.com/QuantumKitHub/TensorKit.jl)). For more information see the
+documentation on [graded arrays](@extref GradedArrays :doc:`user_interface/graded_arrays`).
+
+## Available symmetries
+
+Some standard symmetries are available such as `Z2`, `fU1` (fermionic `U(1)`), and `SU2`. See
+[symmetry sectors](@extref GradedArrays Symmetry-sectors) for the complete list and more details.
+
+You can use named sectors to conserve a product of symmetries.
+
+```@example symmetries
+Index([(; charge = U1(0), spin = U1(1)) => 1, (; charge = U1(1), spin = U1(0)) => 2])
+```

examples/README.jl

@@ -7,6 +7,11 @@
 # [![Code Style](https://img.shields.io/badge/code_style-ITensor-purple)](https://github.com/ITensor/ITensorFormatter.jl)
 # [![Aqua](https://raw.githubusercontent.com/JuliaTesting/Aqua.jl/master/badge.svg)](https://github.com/JuliaTesting/Aqua.jl)
 
+# A next-generation rewrite of [ITensors.jl](https://github.com/ITensor/ITensors.jl), allowing
+# arbitrary array backends and a wider range of symmetries. Built on top of
+# [TensorAlgebra.jl](https://github.com/ITensor/TensorAlgebra.jl), with group symmetric tensors
+# powered by [GradedArrays.jl](https://github.com/ITensor/GradedArrays.jl).
+
 # ## Support
 #
 # {CCQ_LOGO}

ext/ITensorBaseGradedArraysExt.jl

@@ -1,13 +1,19 @@
 module ITensorBaseGradedArraysExt
 
-using GradedArrays: FusedGradedDiagonal, FusedGradedMatrix, SectorRange
+using GradedArrays: GradedArrays as GA, FusedGradedDiagonal, FusedGradedMatrix
 using ITensorBase: ITensorBase, NamedTensor, name, uniquename, unnamed
 using Random: AbstractRNG, default_rng
 using TensorAlgebra: TensorAlgebra as TA
-using TensorKitSectors: Sector
+using TensorKitSectors: TensorKitSectors as TKS
 
 const NamedUnitRange = ITensorBase.NamedUnitRange
 
+# The `sector => multiplicity` pairs `gradedrange` takes, so an `Index` over a graded space
+# prints what it was written with rather than the `gradedrange(...)` call that built the range.
+function ITensorBase.from_range(g::GA.AbstractGradedOneTo)
+    return [s => m for (s, m) in zip(GA.sectors(g), GA.datalengths(g))]
+end
+
 # GradedArrays defines `unmatricize` for its fused matrices with untyped axes, which ties with
 # ITensorBase's methods on named axes. Restate those for each fused matrix type so the named
 # unmatricize of a graded matrix has a unique most-specific method.
@@ -39,8 +45,8 @@ end
 
 # Flux-canceling constructors at the `Index` level: delegate to the GradedArrays flux backend on
 # the unnamed axes, then reattach names, so the flux convention lives only in the backend. The
-# sector may be a bare `TensorKitSectors.Sector` or a `SectorRange`; this is an extension because
-# ITensorBase does not depend on the sector types.
+# sector may be a `TKS.Sector` or a `GA.Sector`, the same pair GradedArrays' own flux constructors
+# dispatch on. This is an extension because ITensorBase does not depend on the sector types.
 
 # Name the delegated result: the physical-leg names followed by a fresh name for the dangling aux
 # leg, minted of the legs' name type (not hardcoded to `IndexName`).
@@ -53,7 +59,7 @@ end
 # Three signature groups, each carrying a named physical axis so overloading `Base` is not piracy:
 # nonempty codomain with a (possibly empty) domain, the codomain-only form, and empty codomain with
 # a nonempty domain. The all-empty flux-only case has no named leg and is left to the backend.
-for S in (Sector, SectorRange)
+for S in (TKS.Sector, GA.Sector)
     # Nonempty codomain, domain given (possibly empty).
     for f in (:rand, :randn)
         @eval begin

src/abstractnamedtensor.jl

@@ -1496,8 +1496,9 @@ end
 # Copy of `Base.dims2string` defined in `show.jl`.
 function dims_to_string(d)
     isempty(d) && return "0-dimensional"
-    length(d) == 1 && return "$(d[1])-element"
-    return join(map(string, d), '×')
+    strs = map(x -> sprint(show, x; context = :compact => true), d)
+    length(d) == 1 && return "$(only(strs))-element"
+    return join(strs, '×')
 end
 
 function concretetype_to_string_truncated(

src/index.jl

@@ -424,12 +424,28 @@ function primestring(plev)
     end
 end
 
+# The space, so a graded index shows its sectors and its arrow rather than just a total length.
+# A compact context asks for the length instead, which is what a tensor's summary line uses to
+# stay readable with one entry per leg.
 function Base.show(io::IO, i::Index)
-    lenstr = "length=$(length(i))"
+    sp = space(i)
+    # A dual index prints as `dual` of the non-dual one, which is the call that makes it, rather
+    # than as a `dual` around the space, which is not a call at all once the space is a vector of
+    # `sector => multiplicity` pairs.
+    nondual = TA.isdual(sp) ? TA.dual(sp) : sp
+    spacestr = if get(io, :compact, false)
+        "length=$(length(i))"
+    else
+        # The space the `Index` was written with rather than the range it stores. `:typeinfo`
+        # drops the element-type prefix a vector of pairs would otherwise carry.
+        spec = from_range(nondual)
+        sprint(show, spec; context = IOContext(io, :typeinfo => typeof(spec)))
+    end
     idstr = "|id=$(shortid(uuid(i)))"
     tagsstr = !isempty(tags_stored(i)) ? "|$(tagsstring(tags_stored(i)))" : ""
     primestr = primestring(plev(i))
-    str = "Index($(lenstr)$(idstr)$(tagsstr))$(primestr)"
+    str = "Index($(spacestr)$(idstr)$(tagsstr))$(primestr)"
+    TA.isdual(sp) && (str = "dual($(str))")
     print(io, str)
     return nothing
 end

src/namedunitrange.jl

@@ -38,6 +38,14 @@ dropped. Equal to [`unnamed`](@ref) for a `NamedUnitRange`.
 """
 space(i::NamedUnitRange) = unnamed(i)
 
+# The space `r` was built from, the inverse of `to_range`: what a range is written as rather
+# than what it stores, so a `Base.OneTo` gives back its length and, with GradedArrays loaded, a
+# graded space gives back its `sector => multiplicity` pairs. Used for printing, and overloaded
+# by the extension of whichever package defines the space. Belongs beside `to_range` in
+# TensorAlgebra, where every package that takes a space specification could reach it.
+from_range(r) = r
+from_range(r::Base.OneTo) = length(r)
+
 # Construct from a space, minting a fresh name of the requested flavor. The space is
 # anything `to_range` accepts (an `Integer`, an existing range, or a sector-pair vector
 # when GradedArrays is loaded), so `Index(2)`, `Index(1:3)`, and

test/Project.toml

@@ -32,7 +32,7 @@ AbstractTrees = "0.4.5"
 Adapt = "4"
 Aqua = "0.8.9"
 Combinatorics = "1"
-GradedArrays = "0.16.6"
+GradedArrays = "0.17"
 ITensorBase = "0.15"
 ITensorPkgSkeleton = "0.3.42"
 JLArrays = "0.2, 0.3"

test/test_basics.jl

@@ -196,11 +196,11 @@ using UUIDs: UUID
     @testset "show" begin
         i = Index(2)
         @test sprint(show, "text/plain", i) ==
-            "Index(length=2|id=$(first(string(uuid(i)), 8)))"
+            "Index(2|id=$(first(string(uuid(i)), 8)))"
 
         i = settag(Index(2), "X", "Y")
         @test sprint(show, "text/plain", i) ==
-            "Index(length=2|id=$(first(string(uuid(i)), 8))|X=>Y)"
+            "Index(2|id=$(first(string(uuid(i)), 8))|X=>Y)"
     end
     @testset "whole-tensor index manipulation" begin
         elt = Float64

test/test_gradedarraysext.jl

@@ -1,5 +1,5 @@
-using GradedArrays: U1, sectors
-using ITensorBase: ITensorBase, ITensor, Index, align, inds, names, prime, space, unnamed
+using GradedArrays: U1, fU1, sectors
+using ITensorBase: ITensorBase, ITensor, Index, align, inds, prime, space, unnamed, uuid
 using StableRNGs: StableRNG
 using TensorAlgebra: TensorAlgebra, dual, isdual, matricize, project, project_aux,
     tryproject, tryproject_aux, unchecked_project, unchecked_project_aux, unmatricize
@@ -40,16 +40,14 @@ using Test: @test, @test_throws, @testset
     @test length(inds(randn(rng, U1(1), (i, j)))) == 3
     @test length(inds(randn(rng, U1(1), (i,), (j,)))) == 3
 
-    # A bare `TensorKitSectors.Sector` (fermionic) works as the flux.
-    s = [
-        Index([FermionNumber(0) => 2, FermionNumber(1) => 2]; tags = "s" => "$n") for
-            n in 1:4
-    ]
+    # The flux may be a bare `TensorKitSectors.Sector` even where the axes are graded by the
+    # GradedArrays sector, and the aux leg comes back carrying the GradedArrays one.
+    s = [Index([fU1(0) => 2, fU1(1) => 2]; tags = "s" => "$n") for n in 1:4]
     t = randn(rng, elt, FermionNumber(2), (s[1], s[2], s[3], s[4]))
     @test length(inds(t)) == 5
     auxt = only(setdiff(collect(inds(t)), s))
     @test isdual(auxt) && length(auxt) == 1 &&
-        only(sectors(space(auxt))) == FermionNumber(2)
+        only(sectors(space(auxt))) == fU1(2)
 
     # `zeros`/`ones`/`fill` mirror `randn` (`fill` takes the value first). Each carries the
     # flux on an aux leg the same way.
@@ -123,7 +121,7 @@ end
     @test m isa AbstractMatrix{elt}
     @test size(m) == (length(i) * length(j), length(k))
     rt = unmatricize(m, (i, j), (k,))
-    @test names(rt) == names(a)
+    @test ITensorBase.names(rt) == ITensorBase.names(a)
     @test isdual(inds(rt)[3])
     @test unnamed(rt) ≈ unnamed(a)
 end
@@ -187,3 +185,14 @@ end
     # Naming the dimensions flat claims no split, so it stays available.
     @test ITensor(m, (i, dual(j))) isa ITensor
 end
+
+# An `Index` prints the space it was written with, so a graded one shows its
+# `sector => multiplicity` pairs rather than the `gradedrange(...)` call that built the range,
+# and a dual one prints as `dual` of the index rather than of the pairs.
+@testset "GradedArraysExt Index show" begin
+    i = Index([U1(0) => 1, U1(1) => 2])
+    @test sprint(show, "text/plain", i) ==
+        "Index([U1(0) => 1, U1(1) => 2]|id=$(first(string(uuid(i)), 8)))"
+    @test sprint(show, "text/plain", dual(i)) ==
+        "dual(Index([U1(0) => 1, U1(1) => 2]|id=$(first(string(uuid(i)), 8))))"
+end

test/test_linearalgebra.jl

@@ -1,5 +1,5 @@
 import LinearAlgebra as LA
-using ITensorBase: Named, names, unname, unnamed
+using ITensorBase: ITensorBase, Named, unname, unnamed
 using Test: @test, @testset
 
 @testset "LinearAlgebra (eltype=$(elt))" for elt in
@@ -14,5 +14,5 @@ using Test: @test, @testset
     @test unnamed(LA.lmul!(2, copy(a))) ≈ 2 * unnamed(a)
     @test unnamed(LA.rdiv!(copy(a), 2)) ≈ unnamed(a) / 2
     @test unnamed(LA.ldiv!(2, copy(a))) ≈ 2 \ unnamed(a)
-    @test LA.dot(a, b) ≈ LA.dot(unnamed(a), unname(b, names(a)))
+    @test LA.dot(a, b) ≈ LA.dot(unnamed(a), unname(b, ITensorBase.names(a)))
 end

test/test_mooncakeext.jl

@@ -1,4 +1,4 @@
-using ITensorBase: Name, NamedTensor, NamedUnitRange, inds, name, nameperm, names,
+using ITensorBase: ITensorBase, Name, NamedTensor, NamedUnitRange, inds, name, nameperm,
     names_setdiff, to_inds, uniquename
 using LinearAlgebra: mul!
 using Mooncake: Mooncake
@@ -23,8 +23,8 @@ using Test: @test, @testset
         Mooncake.TestUtils.test_rule(
             rng, nameperm, a1, (i,), (j,); mode, is_primitive
         )
-        Mooncake.TestUtils.test_rule(rng, names, a1; mode, is_primitive)
-        Mooncake.TestUtils.test_rule(rng, names, a1, 1; mode, is_primitive)
+        Mooncake.TestUtils.test_rule(rng, ITensorBase.names, a1; mode, is_primitive)
+        Mooncake.TestUtils.test_rule(rng, ITensorBase.names, a1, 1; mode, is_primitive)
         Mooncake.TestUtils.test_rule(rng, inds, a1; mode, is_primitive)
         Mooncake.TestUtils.test_rule(rng, inds, a1, 1; mode, is_primitive)
         Mooncake.TestUtils.test_rule(