diff --git a/Changelog.md b/Changelog.md index db8eb1f41..2f702c908 100644 --- a/Changelog.md +++ b/Changelog.md @@ -15,6 +15,24 @@ Features that are planned to be implemented before the release of v1.0.0, in no # Changelog +## v0.14 + +### `DiagonalTensorMap` and `reduceddim` + +This adds a `DiagonalTensorMap` type for representing tensor maps in which all of the +blocks are diagonal. This only makes sense for maps between a single index in the domain and +the codomain (which are furthermore required to be the same), as otherwise the fusion and +splitting trees from the domain and codomain to the blocked sectors would itself be +nondiagonal. This new type will be used to capture the singular values and eigenvalues as +tensor maps in the corresponding decompositions. The number of free parameters in a +`DiagonalTensorMap` instance on vector space `V` is equal to `reduceddim(V)`, a new function +that sums up the degeneracy dimension `dim(V, c)` for each of the sectors `c` in `V`. This +function can be useful by itself and is also exported from the `TensorKit` module. An +instance of `DiagonalTensorMap` can then be created as `DiagonalTensorMap(data, V)` where +`data` is a vector of length `reduceddim(V)`. + +## v0.13 + ### `AbstractTensorMap{E,S,N₁,N₂}` This adds the scalar type as a parameter to the `AbstractTensorMap` type. This is useful in diff --git a/Project.toml b/Project.toml index 0cfc5f922..6aca1592b 100644 --- a/Project.toml +++ b/Project.toml @@ -40,7 +40,7 @@ TensorOperations = "5.1" Test = "1" TestExtras = "0.2,0.3" TupleTools = "1.1" -VectorInterface = "0.4" +VectorInterface = "0.4, 0.5" julia = "1.10" [extras] diff --git a/docs/src/lib/sectors.md b/docs/src/lib/sectors.md index 2e0f763ee..4764094ad 100644 --- a/docs/src/lib/sectors.md +++ b/docs/src/lib/sectors.md @@ -114,8 +114,7 @@ are then used in the index manipulation of `AbstractTensorMap` objects. The foll defined on fusion splitting tree pairs have an associated definition for tensors. ```@docs repartition -transpose(f₁::FusionTree{I}, f₂::FusionTree{I}, - p1::IndexTuple{N₁}, p2::IndexTuple{N₂}) where {I<:Sector,N₁,N₂} -braid(f₁::FusionTree{I}, f₂::FusionTree{I}, levels1::IndexTuple, levels2::IndexTuple, p1::IndexTuple{N₁}, p2::IndexTuple{N₂}) where {I<:Sector,N₁,N₂} -permute(f₁::FusionTree{I}, f₂::FusionTree{I}, p1::IndexTuple{N₁}, p2::IndexTuple{N₂}) where {I<:Sector,N₁,N₂} +transpose(::FusionTree{I}, ::FusionTree{I}, ::IndexTuple{N₁}, ::IndexTuple{N₂}) where {I<:Sector,N₁,N₂} +braid(::FusionTree{I}, ::FusionTree{I}, ::IndexTuple, ::IndexTuple, ::IndexTuple{N₁}, ::IndexTuple{N₂}) where {I<:Sector,N₁,N₂} +permute(::FusionTree{I}, ::FusionTree{I}, ::IndexTuple{N₁}, ::IndexTuple{N₂}) where {I<:Sector,N₁,N₂} ``` diff --git a/docs/src/lib/spaces.md b/docs/src/lib/spaces.md index 8fea48216..2a2b5696d 100644 --- a/docs/src/lib/spaces.md +++ b/docs/src/lib/spaces.md @@ -59,7 +59,7 @@ const CU₁Space = CU1Space const SU₂Space = SU2Space ``` -## Methods +## [Methods](@id s_spacemethods) Methods often apply similar to e.g. spaces and corresponding tensors or tensor maps, e.g.: @@ -70,14 +70,15 @@ sectors hassector dim(::VectorSpace) dim(::ElementarySpace, ::Sector) +reduceddim dim(P::ProductSpace{<:ElementarySpace,N}, sector::NTuple{N,<:Sector}) where {N} dim(::HomSpace) dims -blocksectors(::ProductSpace) +blocksectors(P::ProductSpace{S,N}) where {S,N} blocksectors(::HomSpace) hasblock blockdim -fusiontrees(::ProductSpace, ::Sector) +fusiontrees(P::ProductSpace{S,N}, blocksector::I) where {S,N,I} space ``` @@ -107,10 +108,11 @@ isisomorphic insertunit ``` -There are also specific methods for `HomSpace` instances, that mimic the effect of that -operation on the corresponding tensor maps: +There are also specific methods for `HomSpace` instances, that are used in determining +the resuling `HomSpace` after applying certain tensor operations. ```@docs -permute(::HomSpace, ::Index2Tuple) -compose(::HomSpace{S}, ::HomSpace{S}) where {S} +TensorKit.permute(::HomSpace{S}, ::Index2Tuple{N₁,N₂}) where {S,N₁,N₂} +TensorKit.select(::HomSpace{S}, ::Index2Tuple{N₁,N₂}) where {S,N₁,N₂} +TensorKit.compose(::HomSpace{S}, ::HomSpace{S}) where {S} ``` diff --git a/docs/src/lib/tensors.md b/docs/src/lib/tensors.md index a0cdfe5c6..83ce15214 100644 --- a/docs/src/lib/tensors.md +++ b/docs/src/lib/tensors.md @@ -14,6 +14,7 @@ AbstractTensorMap The following concrete subtypes are provided within the TensorKit library: ```@docs TensorMap +DiagonalTensorMap AdjointTensorMap BraidingTensor ``` @@ -201,6 +202,11 @@ contract! ## `TensorMap` factorizations +The factorisation methods come in two flavors, namely a non-destructive version where you +can specify an additional permutation of the domain and codomain indices before the +factorisation is performed (provided that `sectorstyle(t)` has a symmetric braiding) as +well as a destructive version The non-destructive methods are given first: + ```@docs leftorth rightorth @@ -209,7 +215,13 @@ rightnull tsvd eigh eig +eigen isposdef ``` +The corresponding destructive methods have an exclamation mark at the end of their name, +and only accept the `TensorMap` object as well as the method-specific algorithm and keyword +arguments. + + TODO: document svd truncation types \ No newline at end of file diff --git a/docs/src/man/fusiontrees.md b/docs/src/man/fusiontrees.md deleted file mode 100644 index a75f7a7f8..000000000 --- a/docs/src/man/fusiontrees.md +++ /dev/null @@ -1,434 +0,0 @@ -## [Fusion trees](@id ss_fusiontrees) -The gain in efficiency (both in memory occupation and computation time) obtained from using -symmetric (equivariant) tensor maps is that, by Schur's lemma, they are block diagonal in -the basis of coupled sectors, i.e. they exhibit block sparsity. To exploit this block -diagonal form, it is however essential that we know the basis transform from the individual -(uncoupled) sectors appearing in the tensor product form of the domain and codomain, to the -totally coupled sectors that label the different blocks. We refer to the latter as block -sectors, as we already encountered in the previous section [`blocksectors`](@ref) and -[`blockdim`](@ref) defined on the type [`ProductSpace`](@ref). - -This basis transform consists of a basis of inclusion and projection maps, denoted as -``X^{a_1a_2…a_N}_{c,α}: R_c → R_{a_1} ⊗ R_{a_2} ⊗ … ⊗ R_{a_N}`` and their adjoints -``(X^{a_1a_2…a_N}_{c,α})^†``, such that - -``(X^{a_1a_2…a_N}_{c,α})^† ∘ X^{a_1a_2…a_N}_{c′,α′} = δ_{c,c′} δ_{α,α′} \mathrm{id}_c`` - -and - -``∑_{c,α} X^{a_1a_2…a_N}_{c,α} ∘ (X^{a_1a_2…a_N}_{c,α})^† = \mathrm{id}_{a_1 ⊗ a_2 ⊗ … ⊗ a_N} = \mathrm{id}_{a_1} ⊗ \mathrm{id}_{a_2} ⊗ … ⊗ \mathrm{id}_{a_N} `` - -Fusion trees provide a particular way to construct such a basis. It is useful to know about -the existence of fusion trees and how they are represented, as discussed in the first -subsection. The next two subsections discuss possible manipulations that can be performed -with fusion trees. These are used under the hood when manipulating the indices of tensors, -but a typical user would not need to use these manipulations on fusion trees directly. -Hence, these last two sections can safely be skipped. - -### Canonical representation - -To couple or fuse the different sectors together into a single block sector, we can -sequentially fuse together two sectors into a single coupled sector, which is then fused -with the next uncoupled sector, using the splitting tensors ``X_{a,b}^{c,μ} : R_c → R_a ⊗ -R_b`` and their adjoints. This amounts to the canonical choice of our tensor product, and -for a given tensor mapping from ``(((W_1 ⊗ W_2) ⊗ W_3) ⊗ … )⊗ W_{N_2})`` to ``(((V_1 ⊗ V_2) -⊗ V_3) ⊗ … )⊗ V_{N_1})``, the corresponding fusion and splitting trees take the form - -![double fusion tree](img/tree-simple.svg) - -for the specific case ``N_1=4`` and ``N_2=3``. We can separate this tree into the fusing -part ``(b_1⊗b_2)⊗b_3 → c`` and the splitting part ``c→(((a_1⊗a_2)⊗a_3)⊗a_4)``. Given that -the fusion tree can be considered to be the adjoint of a corresponding splitting tree -``c→(b_1⊗b_2)⊗b_3``, we now first consider splitting trees in isolation. A splitting tree -which goes from one coupled sectors ``c`` to ``N`` uncoupled sectors ``a_1``, ``a_2``, …, -``a_N`` needs ``N-2`` additional internal sector labels ``e_1``, …, ``e_{N-2}``, and, if -`FusionStyle(I) isa GenericFusion`, ``N-1`` additional multiplicity labels ``μ_1``, -…, ``μ_{N-1}``. We henceforth refer to them as vertex labels, as they are associated with -the vertices of the splitting tree. In the case of `FusionStyle(I) isa UniqueFusion`, the -internal sectors ``e_1``, …, ``e_{N-2}`` are completely fixed, for -the generic case they can also take different values. In our abstract -notation of the splitting basis ``X^{a_1a_2…a_N}_{c,α}`` used above, ``α`` can be consided -a collective label, i.e. ``α = (e_1, …, e_{N-2}; μ₁, … ,μ_{N-1})``. Indeed, we can check -the orthogonality condition -``(X^{a_1a_2…a_N}_{c,α})^† ∘ X^{a_1a_2…a_N}_{c′,α′} = δ_{c,c′} δ_{α,α′} \mathrm{id}_c``, -which now forces all internal lines ``e_k`` and vertex labels ``μ_l`` to be the same. - -There is one subtle remark that we have so far ignored. Within the specific subtypes of -`Sector`, we do not explicitly distinguish between ``R_a^*`` (simply denoted as ``a`^*`` -and graphically depicted as an upgoing arrow ``a``) and ``R_{\bar{a}}`` (simply denoted as -``\bar{a}`` and depicted with a downgoing arrow), i.e. between the dual space of ``R_a`` on -which the conjugated irrep acts, or the irrep ``\bar{a}`` to which the complex conjugate of -irrep ``a`` is isomorphic. This distinction is however important, when certain uncoupled -sectors in the fusion tree actually originate from a dual space. We use the isomorphisms -``Z_a:R_a^* → R_{\bar{a}}`` and its adjoint ``Z_a^†:R_{\bar{a}}→R_a^*``, as introduced in -the section on [topological data of a fusion category](@ref ss_topologicalfusion), to build -fusion and splitting trees that take the distinction between irreps and their conjugates -into account. Hence, in the previous example, if e.g. the first and third space in the -codomain and the second space in the domain of the tensor were dual spaces, the actual pair -of splitting and fusion tree would look like - -![extended double fusion tree](img/tree-extended.svg) - -The presence of these isomorphisms will be important when we start to bend lines, to move -uncoupled sectors from the incoming to the outgoing part of the fusion-splitting tree. Note -that we can still represent the fusion tree as the adjoint of a corresponding splitting -tree, because we also use the adjoint of the ``Z`` isomorphisms in the splitting part, and -the ``Z`` isomorphism in the fusion part. Furthermore, the presence of the ``Z`` -isomorphisms does not affect the orthonormality. - -We represent splitting trees and their adjoints using a specific immutable type called -`FusionTree` (which actually represents a splitting tree, but fusion tree is a more common -term), defined as -```julia -struct FusionTree{I<:Sector,N,M,L} - uncoupled::NTuple{N,I} - coupled::I - isdual::NTuple{N,Bool} - innerlines::NTuple{M,I} # fixed to M = N-2 - vertices::NTuple{L,Int} # fixed to L = N-1 -end -``` -Here, the fields are probably self-explanotary. The `isdual` field indicates whether an -isomorphism is present (if the corresponding value is `true`) or not. Note that the field -`uncoupled` contains the sectors coming out of the splitting trees, before the possible -``Z`` isomorphism, i.e. the splitting tree in the above example would have -`sectors = (a₁, a₂, a₃, a₄)`. The `FusionTree` type has a number of basic properties and -capabilities, such as checking for equality with `==` and support for -`hash(f::FusionTree, h::UInt)`, as splitting and fusion trees are used as keys in look-up -tables (i.e. `AbstractDictionary` instances) to look up certain parts of the data of a -tensor. - -`FusionTree` instances are not checked for consistency (i.e. valid fusion rules etc) upon -creation, hence, they are assumed to be created correctly. The most natural way to create -them is by using the `fusiontrees(uncoupled::NTuple{N,I}, coupled::I = one(I))` method, -which returns an iterator over all possible fusion trees from a set of `N` uncoupled -sectors to a given coupled sector, which by default is assumed to be the trivial sector of -that group or fusion category (i.e. the identity object in categorical nomenclature). The -return type of `fusiontrees` is a custom type `FusionTreeIterator` which conforms to the -complete interface of an iterator, and has a custom `length` function that computes the -number of possible fusion trees without iterating over all of them explicitly. This is best -illustrated with some examples - -```@repl sectors -s = Irrep[SU₂](1/2) -collect(fusiontrees((s,s,s,s))) -collect(fusiontrees((s,s,s,s,s), s, (true, false, false, true, false))) -iter = fusiontrees(ntuple(n->s, 16)) -sum(n->1, iter) -length(iter) -@elapsed sum(n->1, iter) -@elapsed length(iter) -s2 = s ⊠ s -collect(fusiontrees((s2,s2,s2,s2))) -``` -Note that `FusionTree` instances are shown (printed) in a way that is valid code to -reproduce them, a property which also holds for both instances of `Sector` and instances of -`VectorSpace`. All of those should be displayed in a way that can be copy pasted as valid -code. Furthermore, we use contact to determine how to print e.g. a sector. In isolation, -`s2` is printed as `(Irrep[SU₂](1/2) ⊠ Irrep[SU₂](1/2))`, however, within the fusion tree, -it is simply printed as `(1/2, 1/2)`, because it will be converted back into a -`ProductSector`, namely `Irrep[SU₂] ⊠ Irrep[SU₂]` by the constructor of -`FusionTree{Irrep[SU₂] ⊠ Irrep[SU₂]}`. - -### Manipulations on a fusion tree - -We now discuss elementary manipulations that we want to perform on or between fusion trees -(where we actually mean splitting trees), which will form the building block for more -general manipulations on a pair of a fusion and splitting tree discussed in the next -subsection, and then for casting a general index manipulation of a tensor map as a linear -operation in the basis of canonically ordered splitting and fusion trees. In this section, -we will ignore the ``Z`` isomorphisms, as they are just trivially reshuffled under the -different operations that we describe. These manipulations are used as low-level methods by -the `TensorMap` methods discussed on the next page. As such, they are not exported by -TensorKit.jl, nor do they overload similarly named methods from Julia Base (see `split` and -`merge` below). - -The first operation we discuss is an elementary braid of two neighbouring sectors -(indices), i.e. a so-called Artin braid or Artin generator of the braid group. Because -these two sectors do not necessarily appear on the same fusion vertex, some recoupling is necessary. -The following represents two different ways to compute the result of such a braid as a -linear combination of new fusion trees in canonical order: - -![artin braid](img/tree-artinbraid.svg) - -While the upper path is the most intuitive, it requires two recouplings or F-moves (one -forward and one reverse). On the other hand, the lower path requires only one (reverse) F- -move, and two R-moves. The latter are less expensive to compute, and so the lower path is -computationally more efficient. However, the end result should be the same, provided the -pentagon and hexagon equations are satisfied. We always assume that these are satisfied for -any new subtype of `Sector`, and it is up to the user to verify that they are when -implementing new custom `Sector` types. This result is implemented in the function -[`artin_braid(f::FusionTree, i; inv = false)`](@ref TensorKit.artin_braid) where `i` -denotes the position of the first sector (i.e. labeled `b` in the above graph) which is then -braided with the sector at position `i+1` in the fusion tree `f`. The keyword argument `inv` -allows to select the inverse braiding operation, which amounts to replacing the R-matrix -with its inverse (or thus, adjoint) in the above steps. The result is returned as a -dictionary with possible output fusion trees as keys and corresponding coefficients as -value. In the case of `FusionStyle(I) isa UniqueFusion`, their is only one resulting fusion -tree, with corresponding coefficient a complex phase (which is one for the bosonic -representation theory of an Abelian group), and the result is a special -`SingletonDict<:AbstractDict`, a `struct` type defined in TensorKit.jl to hold a single key -value pair. - -With the elementary `artin_braid`, we can then compute a more general braid. For this, we -provide an interface - -[`braid(f::FusionTree{I,N}, levels::NTuple{N,Int}, permutation::NTuple{N,Int})`](@ref) - -where the braid is specified as a permutation, such that the new sector at position `i` was -originally at position `permutation[i]`, and where every uncoupled sector is also assigned -a level or depth. The permutation is decomposed into swaps between neighbouring sectors, -and when two sectors are swapped, their respective level will determine whether the left -sector is braided over or under its right neighbor. This interface does not allow to -specify the most general braid, and in particular will never wind one line around another, -but can be used as a more general building block for arbitrary braids than the elementary -Artin generators. A graphical example makes this probably more clear, i.e for -`levels=(1,2,3,4,5)` and `permutation=(5,3,1,4,2)`, the corresponding braid is given by - -![braid interface](img/tree-braidinterface.svg) - -that is, the first sector or space goes to position 3, and crosses over all other lines, -because it has the lowest level (i.e. think of level as depth in the third dimension), and -so forth. We sketch this operation both as a general braid on the left hand side, and as a -particular composition of Artin braids on the right hand side. - -When `BraidingStyle(I) == SymmetricBraiding()`, there is no distinction between applying -the braiding or its inverse (i.e. lines crossing over or under each other in the graphical -notation) and the whole operation simplifies down to a permutation. We then also support -the interface - -[`permute(f::FusionTree{I,N}, permutation::NTuple{N,Int})`](@ref) - -Other manipulations which are sometimes needed are - -* [insertat(f1::FusionTree{I,N₁}, i::Int, f2::FusionTree{I,N₂})](@ref TensorKit.insertat) : - inserts a fusion tree `f2` at the `i`th uncoupled sector of fusion tree `f1` (this - requires that the coupled sector `f2` matches with the `i`th uncoupled sector of `f1`, - and that `!f1.isdual[i]`, i.e. that there is no ``Z``-isomorphism on the `i`th line of - `f1`), and recouple this into a linear combination of trees in canonical order, with - `N₁+N₂-1` uncoupled sectors, i.e. diagrammatically for `i=3` - - ![insertat](img/tree-insertat.svg) - -* [split(f::FusionTree{I,N}, M::Int)](@ref TensorKit.split) : - splits a fusion tree `f` into two trees `f1` and `f2`, such that `f1` has the first `M` - uncoupled sectors of `f`, and `f2` the remaining `N-M`. This function is type stable if `M` is a compile time constant. - - `split(f, M)` is the inverse of `insertat` in the sence that `insertat(f2, 1, f1)` - should return a dictionary with a single key-value pair `f=>1`. Diagrammatically, for - `M=4`, the function `split` returns - - ![split](img/tree-split.svg) - -* [merge(f1::FusionTree{I,N₁}, f2::FusionTree{I,N₂}, c::I, μ=nothing)](@ref TensorKit.merge) : - merges two fusion trees `f1` and `f2` by fusing the coupled sectors of `f1` and `f2` - into a sector `c` (with vertex label `μ` if `FusionStyle(I) == GenericFusion()`), - and reexpressing the result as a linear combination of fusion trees with `N₁+N₂` - uncoupled sectors in canonical order. This is a simple application of `insertat`. - Diagrammatically, this operation is represented as: - - ![merge](img/tree-merge.svg) - -### Manipulations on a splitting - fusion tree pair - -In this subsection we discuss manipulations that act on a splitting and fusion tree pair, -which we will always as two separate trees `f1, f2`, where `f1` is the splitting tree and -`f2` represents the fusion tree, and they should have `f1.coupled == f2.coupled`. - -The most important manipulation on such a pair is to move sectors from one to the other. -Given the canonical order of these trees, we exclusively use the *left duality* (see the -section on [categories](@ref s_categories)), for which the evaluation and coevaluation maps -establish isomorphisms between - -``\mathrm{Hom}((((b_1 ⊗ b_2) ⊗ …) ⊗ b_{N_2}), (((a_1 ⊗ a_2) ⊗ …) ⊗ a_{N_1}))`` - -`` ≂ \mathrm{Hom}((((b_1 ⊗ b_2) ⊗ ...) ⊗ b_{N_2-1}), ((((a_1 ⊗ a_2) ⊗ ...) ⊗ a_{N_1}) ⊗ b_{N_2}^*))`` - -`` ≂ \mathrm{Hom}(1, (((((((a_1 ⊗ a_2) ⊗ ...) ⊗ a_{N_1}) ⊗ b_{N_2}^*) ⊗ …) ⊗ b_2^*) ⊗ b_1^*) )`` - -where the last morphism space is then labeled by the basis of only splitting trees. We can -then use the manipulations from the previous section, and then again use the left duality -to bring this back to a pair of splitting and fusion tree with `N₂′` incoming and `N₁′` -incoming sectors (with `N₁′ + N₂′ == N₁ + N₂`). - -We now discuss how to actually bend lines, and thus, move sectors from the incoming part -(fusion tree) to the outgoing part (splitting tree). Hereby, we exploit the relations -between the (co)evaluation (exact pairing) and the fusion tensors, discussed in -[topological data of a fusion category](@ref ss_topologicalfusion). The main ingredient -that we need is summarized in - -![line bending](img/tree-linebending.svg) - -We will only need the B-symbol and not the A-symbol. Applying the left evaluation on the -second sector of a splitting tensor thus yields a linear combination of fusion tensors -(when `FusionStyle(I) == GenericFusion()`, or just a scalar times the corresponding -fusion tensor otherwise), with corresponding ``Z`` isomorphism. Taking the adjoint of this -relation yields the required relation to transform a fusion tensor into a splitting tensor -with an added ``Z^†`` isomorphism. - -However, we have to be careful if we bend a line on which a ``Z`` isomorphism (or its -adjoint) is already present. Indeed, it is exactly for this operation that we explicitly -need to take the presence of these isomorphisms into account, obtaining the relation - -![dual line bending](img/tree-linebending2.svg) - -Hence, bending an `isdual` sector from the splitting tree to the fusion tree yields an -additional Frobenius-Schur factor, and of course leads to a normal sector (which is no -longer `isdual` and does thus not come with a ``Z``-isomorphism) on the fusion side. We -again use the adjoint of this relation to bend an `isdual` sector from the fusion tree to -the splitting tree. - -The `FusionTree` interface to duality and line bending is given by - -`repartition(f1::FusionTree{I,N₁}, f2::FusionTree{I,N₂}, N::Int)` - -which takes a splitting tree `f1` with `N₁` outgoing sectors, a fusion tree `f2` with `N₂` -incoming sectors, and applies line bending such that the resulting splitting and fusion -trees have `N` outgoing sectors, corresponding to the first `N` sectors out of the list -``(a_1, a_2, …, a_{N_1}, b_{N_2}^*, …, b_{1}^*)`` and `N₁+N₂-N` incoming sectors, -corresponding to the dual of the last `N₁+N₂-N` sectors from the previous list, in reverse. -This return values are correctly inferred if `N` is a compile time constant. - -Graphically, for `N₁ = 4`, `N₂ = 3`, `N = 2` and some particular choice of `isdual` in both -the fusion and splitting tree: - -![repartition](img/tree-repartition.svg) - -The result is returned as a dictionary with keys `(f1′, f2′)` and the corresponding `coeff` -as value. Note that the summation is only over the ``κ_j`` labels, such that, in the case -of `FusionStyle(I) isa MultiplicityFreeFusion`, the linear combination simplifies to -a single term with a scalar coefficient. - -With this basic function, we can now perform arbitrary combinations of braids or -permutations with line bendings, to completely reshuffle where sectors appear. The -interface provided for this is given by - -[`braid(f1::FusionTree{I,N₁}, f2::FusionTree{I,N₂}, levels1::NTuple{N₁,Int}, levels2::NTuple{N₂,Int}, p1::NTuple{N₁′,Int}, p2::NTuple{N₂′,Int})`](@ref) - -where we now have splitting tree `f1` with `N₁` outgoing sectors, a fusion tree `f2` with -`N₂` incoming sectors, `levels1` and `levels2` assign a level or depth to the corresponding -uncoupled sectors in `f1` and `f2`, and we represent the new configuration as a pair `p1` -and `p2`. Together, `(p1..., p2...)` represents a permutation of length `N₁+N₂ = N₁′+N₂′`, -where `p1` indicates which of the original sectors should appear as outgoing sectors in the -new splitting tree and `p2` indicates which appear as incoming sectors in the new fusion -tree. Hereto, we label the uncoupled sectors of `f1` from `1` to `N₁`, followed by the -uncoupled sectors of `f2` from `N₁+1` to `N₁+N₂`. Note that simply repartitioning the -splitting and fusion tree such that e.g. all sectors appear in the new splitting tree (i.e. -are outgoing), amounts to chosing `p1 = (1,..., N₁, N₁+N₂, N₁+N₂-1, ... , N₁+1)` and -`p2=()`, because the duality isomorphism reverses the order of the tensor product. - -This routine is implemented by indeed first making all sectors outgoing using the -`repartition` function discussed above, such that only splitting trees remain, then -braiding those using the routine from the previous subsection such that the new outgoing -sectors appear first, followed by the new incoming sectors (in reverse order), and then -again invoking the `repartition` routine to bring everything in final form. The result is -again returned as a dictionary where the keys are `(f1′,f2′)` and the values the -corresponding coefficients. - -As before, there is a simplified interface for the case where -`BraidingStyle(I) isa SymmetricBraiding` and the levels are not needed. This is simply -given by - -[`permute(f1::FusionTree{I,N₁}, f2::FusionTree{I,N₂}, p1::NTuple{N₁′,Int}, p2::NTuple{N₂′,Int})`](@ref) - -The `braid` and `permute` routines for double fusion trees will be the main access point for -corresponding manipulations on tensors. As a consequence, results from this routine are -memoized, i.e. they are stored in some package-wide 'least-recently used' cache (from -[LRUCache.jl](https://github.com/JuliaCollections/LRUCache.jl)) that can be accessed as -`TensorKit.braidcache`. By default, this cache stores up to `10^5` different `braid` or -`permute` resuls, where one result corresponds to one particular combination of `(f1, f2, -p1, p2, levels1, levels2)`. This should be sufficient for most algorithms. While there are -currently no (official) access methods to change the default settings of this cache (one can -always resort to `resize!(TensorKit.permutecache)` and other methods from LRUCache.jl), this -might change in the future. The use of this cache is however controlled by two constants of -type `RefValue{Bool}`, namely `usebraidcache_abelian` and `usebraidcache_nonabelian`. The -default values are given by `TensorKit.usebraidcache_abelian[] = false` and -`TensorKit.usebraidcache_nonabelian[] = true`, and respectively reflect that the cache is -likely not going to help (or even slow down) fusion trees with -`FusionStyle(f) isa UniqueFusion`, but is probably useful for fusion trees with -`FusionStyle(f) isa MultipleFusion`. One can change these values and test the effect on -their application. - -The existence of `braidcache` also implies that potential inefficiencies in the fusion -tree manipulations (which we nonetheless try to avoid) will not seriously affect -performance of tensor manipulations. - -### Inspecting fusion trees as tensors -For those cases where the fusion and splitting tensors have an explicit representation as -a tensor, i.e. a morphism in the category `Vect` (this essentially coincides with the case -of group representations), this explicit representation can be created, which can be useful -for checking purposes. Hereto, it is necessary that the *splitting tensor* -``X^{ab}_{c,μ}``, i.e. the Clebsch-Gordan coefficients of the group, are encoded via the -routine `fusiontensor(a,b,c [,μ = 1])`, where the last argument is only necessary in -the case of `FusionStyle(I) == GenericFusion()`. We can then convert a -`FusionTree{I,N}` into an `Array`, which will yield a rank `N+1` array where the first `N` -dimensions correspond to the uncoupled sectors, and the last dimension to the coupled -sector. Note that this is mostly useful for the case of `FusionStyle(I) isa MultipleFusion` -groups, as in the case of abelian groups, all irreps are one-dimensional. - -Some examples: -```@repl sectors -s = Irrep[SU₂](1/2) -iter = fusiontrees((s, s, s, s), SU2Irrep(1)) -f = first(iter) -convert(Array, f) - -I ≈ convert(Array, FusionTree((SU₂(1/2),), SU₂(1/2), (false,), ())) -Z = adjoint(convert(Array, FusionTree((SU2Irrep(1/2),), SU2Irrep(1/2), (true,), ()))) -transpose(Z) ≈ frobeniusschur(SU2Irrep(1/2)) * Z - -I ≈ convert(Array, FusionTree((Irrep[SU₂](1),), Irrep[SU₂](1), (false,), ())) -Z = adjoint(convert(Array, FusionTree((Irrep[SU₂](1),), Irrep[SU₂](1), (true,), ()))) -transpose(Z) ≈ frobeniusschur(Irrep[SU₂](1)) * Z - -#check orthogonality -for f1 in iter - for f2 in iter - dotproduct = dot(convert(Array, f1), convert(Array, f2)) - println("< $f1, $f2> = $dotproduct") - end -end -``` -Note that we take the adjoint when computing `Z`, because `convert(Array, f)` assumes `f` -to be splitting tree, which is built using ``Z^†``. Further note that the normalization -(squared) of a fusion tree is given by the dimension of the coupled sector, as we are also -tracing over the ``\mathrm{id}_c`` when checking the orthogonality by computing `dot` of -the corresponding tensors. - -## Fermions - -TODO: Update the documentation for this section. - -Fermionic sectors are represented by the type [`FermionParity`](@ref), which effectively -behaves like a ℤ₂ sector, but with two modifications. Firstly, the exchange of two sectors -with odd fermion parity should yield a minus sign, which is taken care of by virtue of the -R-symbol. This ensures that permuting tensors behave as expected. Secondly, diagrams with -self-crossing lines (aka twists) give rise to a minus sign for odd fermion parity. This is -in essence equivalent to having supertraces, which is what ensures that `@tensor` has a -result that is invariant under permutation of its input tensors. This does however lead to -unwanted minus signs for certain types of diagrams. To avoid this, the `@planar` macro does -not include a supertrace, but requires a manual resolution of all crossings in the diagram. - -## Anyons - -There is currently one example of a `Sector` subtype that has anyonic braiding style, -namely that of the Fibonacci fusion category. It has to (isomorphism classes of) simple -objects, namely the identity `𝟙` and a non-trivial object known as `τ`, with fusion rules -`τ ⊗ τ = 𝟙 ⊕ τ`. Let's summarize the topological data - -```@repl sectors -𝟙 = FibonacciAnyon(:I) -τ = FibonacciAnyon(:τ) -collect(τ ⊗ τ) -FusionStyle(τ) -BraidingStyle(τ) -dim(𝟙) -dim(τ) -F𝟙 = Fsymbol(τ,τ,τ,𝟙,τ,τ) -Fτ = [Fsymbol(τ,τ,τ,τ,𝟙,𝟙) Fsymbol(τ,τ,τ,τ,𝟙,τ); Fsymbol(τ,τ,τ,τ,τ,𝟙) Fsymbol(τ,τ,τ,τ,τ,τ)] -Fτ'*Fτ -polar(x) = rationalize.((abs(x), angle(x)/(2pi))) -Rsymbol(τ,τ,𝟙) |> polar -Rsymbol(τ,τ,τ) |> polar -twist(τ) |> polar -``` diff --git a/docs/src/man/sectors.md b/docs/src/man/sectors.md index 0823b5d56..81dbfd35d 100644 --- a/docs/src/man/sectors.md +++ b/docs/src/man/sectors.md @@ -1028,7 +1028,7 @@ value pair. With the elementary `artin_braid`, we can then compute a more general braid. For this, we provide an interface -[`braid(f::FusionTree{I,N}, levels::NTuple{N,Int}, permutation::NTuple{N,Int})`](@ref) +[`braid(f::FusionTree{I,N}, levels::NTuple{N,Int}, permutation::NTuple{N,Int})`](@ref braid(f::FusionTree{I,N}, levels::NTuple{N,Int}, p::NTuple{N,Int}) where {I<:Sector,N}) where the braid is specified as a permutation, such that the new sector at position `i` was originally at position `permutation[i]`, and where every uncoupled sector is also assigned @@ -1052,7 +1052,7 @@ the braiding or its inverse (i.e. lines crossing over or under each other in the notation) and the whole operation simplifies down to a permutation. We then also support the interface -[`permute(f::FusionTree{I,N}, permutation::NTuple{N,Int})`](@ref) +[`permute(f::FusionTree{I,N}, permutation::NTuple{N,Int})`](@ref permute(f::FusionTree{I,N}, p::NTuple{N,Int}) where {I<:Sector,N}) Other manipulations which are sometimes needed are @@ -1135,7 +1135,7 @@ the splitting tree. The `FusionTree` interface to duality and line bending is given by -`repartition(f1::FusionTree{I,N₁}, f2::FusionTree{I,N₂}, N::Int)` +[`repartition(f1::FusionTree{I,N₁}, f2::FusionTree{I,N₂}, N::Int)`](@ref repartition) which takes a splitting tree `f1` with `N₁` outgoing sectors, a fusion tree `f2` with `N₂` incoming sectors, and applies line bending such that the resulting splitting and fusion @@ -1158,7 +1158,7 @@ With this basic function, we can now perform arbitrary combinations of braids or permutations with line bendings, to completely reshuffle where sectors appear. The interface provided for this is given by -[`braid(f1::FusionTree{I,N₁}, f2::FusionTree{I,N₂}, levels1::NTuple{N₁,Int}, levels2::NTuple{N₂,Int}, p1::NTuple{N₁′,Int}, p2::NTuple{N₂′,Int})`](@ref) +[`braid(f1::FusionTree{I,N₁}, f2::FusionTree{I,N₂}, levels1::NTuple{N₁,Int}, levels2::NTuple{N₂,Int}, p1::NTuple{N₁′,Int}, p2::NTuple{N₂′,Int})`](@ref braid(::FusionTree{I}, ::FusionTree{I}, ::IndexTuple, ::IndexTuple, ::IndexTuple{N₁}, ::IndexTuple{N₂}) where {I<:Sector,N₁,N₂}) where we now have splitting tree `f1` with `N₁` outgoing sectors, a fusion tree `f2` with `N₂` incoming sectors, `levels1` and `levels2` assign a level or depth to the corresponding @@ -1184,7 +1184,7 @@ As before, there is a simplified interface for the case where `BraidingStyle(I) isa SymmetricBraiding` and the levels are not needed. This is simply given by -[`permute(f1::FusionTree{I,N₁}, f2::FusionTree{I,N₂}, p1::NTuple{N₁′,Int}, p2::NTuple{N₂′,Int})`](@ref) +[`permute(f1::FusionTree{I,N₁}, f2::FusionTree{I,N₂}, p1::NTuple{N₁′,Int}, p2::NTuple{N₂′,Int})`](@ref permute(::FusionTree{I}, ::FusionTree{I}, ::IndexTuple{N₁}, ::IndexTuple{N₂}) where {I<:Sector,N₁,N₂}) The `braid` and `permute` routines for double fusion trees will be the main access point for corresponding manipulations on tensors. As a consequence, results from this routine are diff --git a/docs/src/man/spaces.md b/docs/src/man/spaces.md index 8128aa8f4..df59733ee 100644 --- a/docs/src/man/spaces.md +++ b/docs/src/man/spaces.md @@ -172,7 +172,7 @@ InnerProductStyle(ℂ^5) One more important instance of `ElementarySpace` is the `GradedSpace`, which is used to represent a graded complex vector space with Euclidean inner product, where the grading is provided by the irreducible representations of a group, or more generally, the simple -objects of a fusion category. We refer to the subsection on [graded spaces](@ref s_rep) on +objects of a fusion category. We refer to the subsection on [graded spaces](@ref ss_rep) on the [next page](@ref s_sectorsrepfusion) for further information about `GradedSpace`. ## [Composite spaces](@id ss_compositespaces) @@ -301,7 +301,8 @@ isomorphic to `V` but has `isdual(flip(V)) == isdual(V')`, i.e. if `V` is a norm than `flip(V)` is a dual space. `flip(V)` is different from `dual(V)` in the case of [`GradedSpace`](@ref). It is useful to flip a tensor index from a ket to a bra (or vice versa), by contracting that index with a unitary map from `V1` to `flip(V1)`. We refer -to [Index operations](@ref) for further information. Some examples: +to the reference on [vector space methods](@ref s_spacemethods) for further information. +Some examples: ```@repl tensorkit ℝ^3 ≾ ℝ^5 ℂ^3 ≾ (ℂ^5)' diff --git a/docs/src/man/tensors.md b/docs/src/man/tensors.md index 594812201..1d6f14fec 100644 --- a/docs/src/man/tensors.md +++ b/docs/src/man/tensors.md @@ -915,8 +915,8 @@ to either `eig` or `eigh`. Other factorizations that are provided by TensorKit.jl are orthogonal or unitary in nature, and thus always require a `AbstractEuclideanTensorMap`. However, they don't require equal domain and codomain. Let us first discuss the *singular value decomposition*, for which we -define and export the methods [`tsvd`](@ref) and [`tsvd!`](@ref) (where as always, the -latter destroys the input). +define and export the methods [`tsvd`](@ref) and `tsvd!` (where as always, the latter +destroys the input). ```julia U, Σ, Vʰ, ϵ = tsvd(t; trunc = notrunc(), p::Real = 2, @@ -934,7 +934,8 @@ given by `min(fuse(codomain(t)), fuse(domain(t)))`. The singular values are cont and are stored on the diagonal of a (collection of) `DenseMatrix` instance(s), similar to the eigenvalues before. -The keyword argument `trunc` provides a way to control the truncation, and is connected to the keyword argument `p`. The default value `notrunc()` implies no truncation, and thus +The keyword argument `trunc` provides a way to control the truncation, and is connected to +the keyword argument `p`. The default value `notrunc()` implies no truncation, and thus `ϵ = 0`. Other valid options are * `truncerr(η::Real)`: truncates such that the `p`-norm of the truncated singular values diff --git a/src/TensorKit.jl b/src/TensorKit.jl index 569532d37..8cdeaafb6 100644 --- a/src/TensorKit.jl +++ b/src/TensorKit.jl @@ -25,12 +25,13 @@ export Vect, Rep # space constructors export CompositeSpace, ProductSpace # composite spaces export FusionTree export IndexSpace, HomSpace, TensorSpace, TensorMapSpace -export AbstractTensorMap, AbstractTensor, TensorMap, Tensor, BraidingTensor # tensors and tensor properties +export AbstractTensorMap, AbstractTensor, TensorMap, Tensor # tensors and tensor properties +export DiagonalTensorMap, BraidingTensor export TruncationScheme export SpaceMismatch, SectorMismatch, IndexError # error types # general vector space methods -export space, field, dual, dim, dims, fuse, flip, isdual, insertunit, oplus +export space, field, dual, dim, reduceddim, dims, fuse, flip, isdual, insertunit, oplus # partial order for vector spaces export infimum, supremum, isisomorphic, ismonomorphic, isepimorphic @@ -189,6 +190,7 @@ include("tensors/vectorinterface.jl") include("tensors/tensoroperations.jl") include("tensors/treetransformers.jl") include("tensors/indexmanipulations.jl") +include("tensors/diagonal.jl") include("tensors/truncation.jl") include("tensors/factorizations.jl") include("tensors/braidingtensor.jl") diff --git a/src/spaces/gradedspace.jl b/src/spaces/gradedspace.jl index 2d1a92e7a..23c380fda 100644 --- a/src/spaces/gradedspace.jl +++ b/src/spaces/gradedspace.jl @@ -34,7 +34,9 @@ function GradedSpace{I,NTuple{N,Int}}(dims; dual::Bool=false) where {I,N} d = ntuple(n -> 0, N) isset = ntuple(n -> false, N) for (c, dc) in dims - i = findindex(values(I), convert(I, c)) + k = convert(I, c) + i = findindex(values(I), k) + k = dc < 0 && throw(ArgumentError("Sector $k has negative dimension $dc")) isset[i] && throw(ArgumentError("Sector $c appears multiple times")) isset = TupleTools.setindex(isset, true, i) d = TupleTools.setindex(d, dc, i) @@ -50,6 +52,7 @@ function GradedSpace{I,SectorDict{I,Int}}(dims; dual::Bool=false) where {I<:Sect for (c, dc) in dims k = convert(I, c) haskey(d, k) && throw(ArgumentError("Sector $k appears multiple times")) + dc < 0 && throw(ArgumentError("Sector $k has negative dimension $dc")) !iszero(dc) && push!(d, k => dc) end return GradedSpace{I,SectorDict{I,Int}}(d, dual) diff --git a/src/spaces/homspace.jl b/src/spaces/homspace.jl index a2fecea44..608b54041 100644 --- a/src/spaces/homspace.jl +++ b/src/spaces/homspace.jl @@ -129,6 +129,19 @@ Return the `HomSpace` obtained by permuting the indices of the domain and codoma according to the permutation `p₁` and `p₂` respectively. """ function permute(W::HomSpace{S}, (p₁, p₂)::Index2Tuple{N₁,N₂}) where {S,N₁,N₂} + p = (p₁..., p₂...) + TupleTools.isperm(p) && length(p) == numind(W) || + throw(ArgumentError("$((p₁, p₂)) is not a valid permutation for $(W)")) + return select(W, (p₁, p₂)) +end + +""" + select(W::HomSpace, (p₁, p₂)::Index2Tuple{N₁,N₂}) + +Return the `HomSpace` obtained by a selection from the domain and codomain of `W` according +to the indices in `p₁` and `p₂` respectively. +""" +function select(W::HomSpace{S}, (p₁, p₂)::Index2Tuple{N₁,N₂}) where {S,N₁,N₂} cod = ProductSpace{S,N₁}(map(n -> W[n], p₁)) dom = ProductSpace{S,N₂}(map(n -> dual(W[n]), p₂)) return cod ← dom diff --git a/src/spaces/vectorspaces.jl b/src/spaces/vectorspaces.jl index abf6d50ed..1468892e5 100644 --- a/src/spaces/vectorspaces.jl +++ b/src/spaces/vectorspaces.jl @@ -111,6 +111,13 @@ field(V::ElementarySpace) = field(typeof(V)) Return the degeneracy dimension corresponding to the sector `s` of the vector space `V`. """ dim(::ElementarySpace, ::Sector) +@doc """ + reduceddim(V::ElementarySpace) -> Int + +Return the sum of all degeneracy dimensions of the vector space `V`. +""" +reduceddim(V::ElementarySpace) = sum(Base.Fix1(dim, V), sectors(V); init=0) + """ oneunit(V::S) where {S<:ElementarySpace} -> S @@ -256,7 +263,7 @@ field(P::Type{<:CompositeSpace}) = field(spacetype(P)) sectortype(P::Type{<:CompositeSpace}) = sectortype(spacetype(P)) # make ElementarySpace instances behave similar to ProductSpace instances -blocksectors(V::ElementarySpace) = sectors(V) +blocksectors(V::ElementarySpace) = collect(sectors(V)) blockdim(V::ElementarySpace, c::Sector) = dim(V, c) # Specific realizations of ElementarySpace types diff --git a/src/tensors/braidingtensor.jl b/src/tensors/braidingtensor.jl index 921e7f7be..17406b6a2 100644 --- a/src/tensors/braidingtensor.jl +++ b/src/tensors/braidingtensor.jl @@ -178,7 +178,7 @@ blocks(b::BraidingTensor) = blocks(TensorMap(b)) # Index manipulations # ------------------- -has_shared_permute(t::BraidingTensor, args...) = false +has_shared_permute(t::BraidingTensor, ::Index2Tuple) = false function add_transform!(tdst::AbstractTensorMap, tsrc::BraidingTensor, (p₁, p₂)::Index2Tuple, diff --git a/src/tensors/diagonal.jl b/src/tensors/diagonal.jl new file mode 100644 index 000000000..ebae85ccc --- /dev/null +++ b/src/tensors/diagonal.jl @@ -0,0 +1,335 @@ +# DiagonalTensorMap +#==========================================================# +struct DiagonalTensorMap{T,S<:IndexSpace,A<:DenseVector{T}} <: AbstractTensorMap{T,S,1,1} + data::A + domain::S # equals codomain + + # uninitialized constructors + function DiagonalTensorMap{T,S,A}(::UndefInitializer, + dom::S) where {T,S<:IndexSpace,A<:DenseVector{T}} + data = A(undef, reduceddim(dom)) + if !isbitstype(T) + zerovector!(data) + end + return DiagonalTensorMap{T,S,A}(data, dom) + end + # constructors from data + function DiagonalTensorMap{T,S,A}(data::A, + dom::S) where {T,S<:IndexSpace,A<:DenseVector{T}} + T ⊆ field(S) || @warn("scalartype(data) = $T ⊈ $(field(S)))", maxlog = 1) + return new{T,S,A}(data, dom) + end +end + +# Basic methods for characterising a tensor: +#-------------------------------------------- +space(d::DiagonalTensorMap) = d.domain ← d.domain + +storagetype(::Type{<:DiagonalTensorMap{T,S,A}}) where {T,S,A<:DenseVector{T}} = A + +# DiagonalTensorMap constructors +#-------------------------------- +# undef constructors +""" + DiagonalTensorMap{T}(undef, domain::S) where {T,S<:IndexSpace} + # expert mode: select storage type `A` + DiagonalTensorMap{T,S,A}(undef, domain::S) where {T,S<:IndexSpace,A<:DenseVector{T}} + +Construct a `DiagonalTensorMap` with uninitialized data. +""" +function DiagonalTensorMap{T}(::UndefInitializer, V::S) where {T,S<:IndexSpace} + return DiagonalTensorMap{T,S,Vector{T}}(undef, V) +end +DiagonalTensorMap(::UndefInitializer, V::IndexSpace) = DiagonalTensorMap{Float64}(undef, V) + +function DiagonalTensorMap{T}(data::A, V::S) where {T,S<:IndexSpace,A<:DenseVector{T}} + length(data) == reduceddim(V) || + throw(DimensionMismatch("length(data) = $(length(data)) is not compatible with the space $V")) + return DiagonalTensorMap{T,S,A}(data, V) +end + +function DiagonalTensorMap(data::DenseVector{T}, V::IndexSpace) where {T} + return DiagonalTensorMap{T}(data, V) +end + +# TODO: more constructors needed? + +# Special case adjoint: +#----------------------- +Base.adjoint(d::DiagonalTensorMap{<:Real}) = d +Base.adjoint(d::DiagonalTensorMap{<:Complex}) = DiagonalTensorMap(conj(d.data), d.domain) + +# Efficient copy constructors and TensorMap converters +#----------------------------------------------------- +Base.copy(d::DiagonalTensorMap) = typeof(d)(copy(d.data), d.domain) + +function Base.copy!(t::AbstractTensorMap, d::DiagonalTensorMap) + space(t) == space(d) || throw(SpaceMismatch()) + for (c, b) in blocks(d) + copy!(block(t, c), b) + end + return t +end +TensorMap(d::DiagonalTensorMap) = copy!(similar(d), d) +Base.convert(::Type{TensorMap}, d::DiagonalTensorMap) = TensorMap(d) + +function Base.convert(::Type{DiagonalTensorMap{T,S,A}}, + d::DiagonalTensorMap{T,S,A}) where {T,S,A} + return d +end +function Base.convert(D::Type{<:DiagonalTensorMap}, d::DiagonalTensorMap) + return DiagonalTensorMap(convert(storagetype(D), d.data), d.domain) +end + +# Complex, real and imaginary parts +#----------------------------------- +for f in (:real, :imag, :complex) + @eval begin + function Base.$f(d::DiagonalTensorMap) + return DiagonalTensorMap($f(d.data), d.domain) + end + end +end + +# Getting and setting the data at the block level +#------------------------------------------------- +function block(d::DiagonalTensorMap, s::Sector) + sectortype(d) == typeof(s) || throw(SectorMismatch()) + offset = 0 + dom = domain(d)[1] + for c in blocksectors(d) + if c < s + offset += dim(dom, c) + elseif c == s + r = offset .+ (1:dim(dom, c)) + return Diagonal(view(d.data, r)) + else # s not in sectors(t) + break + end + end + return Diagonal(view(d.data, 1:0)) +end + +# TODO: is relying on generic AbstractTensorMap blocks sufficient? + +# Indexing and getting and setting the data at the subblock level +#----------------------------------------------------------------- +@inline function Base.getindex(d::DiagonalTensorMap, + f₁::FusionTree{I,1}, + f₂::FusionTree{I,1}) where {I<:Sector} + s = f₁.uncoupled[1] + s == f₁.coupled == f₂.uncoupled[1] == f₂.coupled || throw(SectorMismatch()) + return block(d, s) + # TODO: do we want a StridedView here? Then we need to allocate a new matrix. +end + +function Base.setindex!(d::DiagonalTensorMap, + v, + f₁::FusionTree{I,1}, + f₂::FusionTree{I,1}) where {I<:Sector} + return copy!(getindex(d, f₁, f₂), v) +end + +function Base.getindex(d::DiagonalTensorMap) + sectortype(d) === Trivial || throw(SectorMismatch()) + return Diagonal(d.data) +end + +# Index manipulations +# ------------------- +function has_shared_permute(d::DiagonalTensorMap, (p₁, p₂)::Index2Tuple) + if p₁ === (1,) && p₂ === (2,) + return true + elseif p₁ === (2,) && p₂ === (1,) # transpose + return sectortype(d) === Trivial + else + return false + end +end + +function permute(d::DiagonalTensorMap, (p₁, p₂)::Index2Tuple{1,1}; + copy::Bool=false) + if p₁ === (1,) && p₂ === (2,) + return copy ? Base.copy(d) : d + elseif p₁ === (2,) && p₂ === (1,) # transpose + if has_shared_permute(d, (p₁, p₂)) # tranpose for bosonic sectors + return DiagonalTensorMap(copy ? Base.copy(d.data) : d.data, dual(d.domain)) + end + d′ = typeof(d)(undef, dual(d.domain)) + for (c, b) in blocks(d) + f = only(fusiontrees(codomain(d), c)) + ((f′, _), coeff) = only(permute(f, f, p₁, p₂)) + c′ = f′.coupled + scale!(block(d′, c′), b, coeff) + end + return d′ + else + throw(ArgumentError("invalid permutation $((p₁, p₂)) for tensor in space $(space(d))")) + end +end + +# VectorInterface +# --------------- +function VectorInterface.zerovector(d::DiagonalTensorMap, ::Type{S}) where {S<:Number} + return DiagonalTensorMap(zerovector(d.data, S), d.domain) +end +function VectorInterface.add(ty::DiagonalTensorMap, tx::DiagonalTensorMap, + α::Number, β::Number) + domain(ty) == domain(tx) || throw(SpaceMismatch("$(space(ty)) ≠ $(space(tx))")) + T = VectorInterface.promote_add(ty, tx, α, β) + return add!(scale!(zerovector(ty, T), ty, β), tx, α) # zerovector instead of similar preserves diagonal structure +end + +# TensorOperations +# ---------------- +function TO.tensoradd_type(TC, A::DiagonalTensorMap, ::Index2Tuple{1,1}, ::Bool) + M = similarstoragetype(A, TC) + return DiagonalTensorMap{TC,spacetype(A),M} +end + +function TO.tensorcontract_type(TC, A::DiagonalTensorMap, ::Index2Tuple{1,1}, ::Bool, + B::DiagonalTensorMap, ::Index2Tuple{1,1}, ::Bool, + ::Index2Tuple{1,1}) + M = similarstoragetype(A, TC) + M == similarstoragetype(B, TC) || + throw(ArgumentError("incompatible storage types:\n$(M) ≠ $(similarstoragetype(B, TC))")) + spacetype(A) == spacetype(B) || throw(SpaceMismatch("incompatible space types")) + return DiagonalTensorMap{TC,spacetype(A),M} +end + +function TO.tensoralloc(::Type{DiagonalTensorMap{T,S,M}}, + structure::TensorMapSpace{S,1,1}, + istemp::Val, + allocator=TO.DefaultAllocator()) where {T,S,M} + domain(structure) == codomain(structure) || throw(ArgumentError("domain ≠ codomain")) + V = only(domain(structure)) + dim = reduceddim(V) + data = TO.tensoralloc(M, dim, istemp, allocator) + return DiagonalTensorMap{T,S,M}(data, V) +end + +# Linear Algebra and factorizations +# --------------------------------- +function one!(d::DiagonalTensorMap) + fill!(d.data, one(eltype(d.data))) + return d +end +function Base.one(d::DiagonalTensorMap) + return DiagonalTensorMap(one.(d.data), d.domain) +end +function Base.zero(d::DiagonalTensorMap) + return DiagonalTensorMap(zero.(d.data), d.domain) +end + +function LinearAlgebra.mul!(dC::DiagonalTensorMap, + dA::DiagonalTensorMap, + dB::DiagonalTensorMap, + α::Number, β::Number) + dC.domain == dA.domain == dB.domain || throw(SpaceMismatch()) + mul!(Diagonal(dC.data), Diagonal(dA.data), Diagonal(dB.data), α, β) + return dC +end + +Base.inv(d::DiagonalTensorMap) = DiagonalTensorMap(inv.(d.data), d.domain) +function Base.:\(d1::DiagonalTensorMap, d2::DiagonalTensorMap) + d1.domain == d2.domain || throw(SpaceMismatch()) + return DiagonalTensorMap(d1.data .\ d2.data, d1.domain) +end +function Base.:/(d1::DiagonalTensorMap, d2::DiagonalTensorMap) + d1.domain == d2.domain || throw(SpaceMismatch()) + return DiagonalTensorMap(d1.data ./ d2.data, d1.domain) +end +function LinearAlgebra.pinv(d::DiagonalTensorMap; kwargs...) + T = eltype(d.data) + atol = get(kwargs, :atol, zero(real(T))) + if iszero(atol) + rtol = get(kwargs, :rtol, zero(real(T))) + else + rtol = sqrt(eps(real(float(oneunit(T))))) * length(d.data) + end + pdata = let tol = max(atol, rtol * maximum(abs, d.data)) + map(x -> abs(x) < tol ? zero(x) : pinv(x), d.data) + end + return DiagonalTensorMap(pdata, d.domain) +end +function LinearAlgebra.isposdef(d::DiagonalTensorMap) + return all(isposdef, d.data) +end + +function eig!(d::DiagonalTensorMap) + return d, one(d) +end +function eigh!(d::DiagonalTensorMap{<:Real}) + return d, one(d) +end +function eigh!(d::DiagonalTensorMap{<:Complex}) + # TODO: should this test for hermiticity? `eigh!(::TensorMap)` also does not do this. + return DiagonalTensorMap(real(d.data), d.domain), one(d) +end + +function leftorth!(d::DiagonalTensorMap; alg=QR(), kwargs...) + @assert alg isa Union{QR,QL} + return one(d), d # TODO: this is only correct for `alg = QR()` or `alg = QL()` +end +function rightorth!(d::DiagonalTensorMap; alg=LQ(), kwargs...) + @assert alg isa Union{LQ,RQ} + return d, one(d) # TODO: this is only correct for `alg = LQ()` or `alg = RQ()` +end +# not much to do here: +leftnull!(d::DiagonalTensorMap; kwargs...) = leftnull!(TensorMap(d); kwargs...) +rightnull!(d::DiagonalTensorMap; kwargs...) = rightnull!(TensorMap(d); kwargs...) + +function tsvd!(d::DiagonalTensorMap; trunc=NoTruncation(), p::Real=2, alg=SDD()) + return _tsvd!(d, alg, trunc, p) +end +# helper function +function _compute_svddata!(d::DiagonalTensorMap, alg::Union{SVD,SDD}) + InnerProductStyle(d) === EuclideanInnerProduct() || throw_invalid_innerproduct(:tsvd!) + I = sectortype(d) + dims = SectorDict{I,Int}() + generator = Base.Iterators.map(blocks(d)) do (c, b) + lb = length(b.diag) + U = zerovector!(similar(b.diag, lb, lb)) + V = zerovector!(similar(b.diag, lb, lb)) + p = sortperm(b.diag; by=abs, rev=true) + for (i, pi) in enumerate(p) + U[pi, i] = MatrixAlgebra.safesign(b.diag[pi]) + V[i, pi] = 1 + end + Σ = abs.(view(b.diag, p)) + dims[c] = lb + return c => (U, Σ, V) + end + SVDdata = SectorDict(generator) + return SVDdata, dims +end + +# matrix functions +for f in + (:exp, :cos, :sin, :tan, :cot, :cosh, :sinh, :tanh, :coth, :atan, :acot, :asinh, :sqrt, + :log, :asin, :acos, :acosh, :atanh, :acoth) + @eval Base.$f(d::DiagonalTensorMap) = DiagonalTensorMap($f.(d.data), d.domain) +end + +# Show +#------ +function Base.summary(io::IO, t::DiagonalTensorMap) + return print(io, "DiagonalTensorMap(", space(t), ")") +end +function Base.show(io::IO, t::DiagonalTensorMap) + summary(io, t) + get(io, :compact, false) && return nothing + println(io, ":") + + if sectortype(t) == Trivial + Base.print_array(io, Diagonal(t.data)) + println(io) + else + for (c, b) in blocks(t) + println(io, "* Data for sector ", c, ":") + Base.print_array(io, b) + println(io) + end + end + return nothing +end diff --git a/src/tensors/indexmanipulations.jl b/src/tensors/indexmanipulations.jl index 5b4cc2ddc..ee7f90fcc 100644 --- a/src/tensors/indexmanipulations.jl +++ b/src/tensors/indexmanipulations.jl @@ -31,12 +31,22 @@ function permute(t::AbstractTensorMap, (p₁, p₂)::Index2Tuple{N₁,N₂}; copy::Bool=false) where {N₁,N₂} space′ = permute(space(t), (p₁, p₂)) # share data if possible + if !copy && p₁ === codomainind(t) && p₂ === domainind(t) + return t + end + # general case + @inbounds begin + return permute!(similar(t, space′), t, (p₁, p₂)) + end +end +function permute(t::TensorMap, (p₁, p₂)::Index2Tuple{N₁,N₂}; copy::Bool=false) where {N₁,N₂} + space′ = permute(space(t), (p₁, p₂)) + # share data if possible if !copy if p₁ === codomainind(t) && p₂ === domainind(t) return t elseif has_shared_permute(t, (p₁, p₂)) - return TensorMap(reshape(t.data, dim(codomain(space′)), dim(domain(space′))), - codomain(space′), domain(space′)) + return TensorMap(t.data, space′) end end # general case @@ -53,6 +63,9 @@ function permute(t::AbstractTensorMap, p::IndexTuple; copy::Bool=false) return permute(t, (p, ()); copy) end +function has_shared_permute(t::AbstractTensorMap, (p₁, p₂)::Index2Tuple) + return (p₁ === codomainind(t) && p₂ === domainind(t)) +end function has_shared_permute(t::TensorMap, (p₁, p₂)::Index2Tuple) if p₁ === codomainind(t) && p₂ === domainind(t) return true diff --git a/src/tensors/linalg.jl b/src/tensors/linalg.jl index 8f543cdab..d60158d44 100644 --- a/src/tensors/linalg.jl +++ b/src/tensors/linalg.jl @@ -345,6 +345,7 @@ end function LinearAlgebra.pinv(t::AbstractTensorMap; kwargs...) T = float(scalartype(t)) tpinv = similar(t, T, domain(t) ← codomain(t)) + # TODO: fix so that `rtol` used total tensor norm instead of per block for (c, b) in blocks(t) copy!(block(tpinv, c), pinv(b; kwargs...)) end diff --git a/src/tensors/tensor.jl b/src/tensors/tensor.jl index 90596086f..7f071f104 100644 --- a/src/tensors/tensor.jl +++ b/src/tensors/tensor.jl @@ -18,6 +18,9 @@ struct TensorMap{T,S<:IndexSpace,N₁,N₂,A<:DenseVector{T}} <: AbstractTensorM A<:DenseVector{T}} d = fusionblockstructure(space).totaldim data = A(undef, d) + if !isbitstype(T) + zerovector!(data) + end return TensorMap{T,S,N₁,N₂,A}(data, space) end diff --git a/src/tensors/tensoroperations.jl b/src/tensors/tensoroperations.jl index 02b7c9021..cf977690f 100644 --- a/src/tensors/tensoroperations.jl +++ b/src/tensors/tensoroperations.jl @@ -59,7 +59,8 @@ end function TO.tensoradd_structure(A::AbstractTensorMap, pA::Index2Tuple{N₁,N₂}, conjA::Bool) where {N₁,N₂} if !conjA - return permute(space(A), pA) + # don't use `permute` as this is also used when indices are traced + return select(space(A), pA) else return TO.tensoradd_structure(adjoint(A), adjointtensorindices(A, pA), false) end @@ -172,7 +173,7 @@ function trace_permute!(tdst::AbstractTensorMap, N₁, N₂ = length(p₁), length(p₂) @boundscheck begin - space(tdst) == permute(space(tsrc), (p₁, p₂)) || + space(tdst) == select(space(tsrc), (p₁, p₂)) || throw(SpaceMismatch("trace: tsrc = $(codomain(tsrc))←$(domain(tsrc)), tdst = $(codomain(tdst))←$(domain(tdst)), p₁ = $(p₁), p₂ = $(p₂)")) all(i -> space(tsrc, q₁[i]) == dual(space(tsrc, q₂[i])), 1:N₃) || diff --git a/src/tensors/vectorinterface.jl b/src/tensors/vectorinterface.jl index c9d83ec82..635eed660 100644 --- a/src/tensors/vectorinterface.jl +++ b/src/tensors/vectorinterface.jl @@ -23,7 +23,7 @@ VectorInterface.zerovector!!(t::AbstractTensorMap) = zerovector!(t) #------------------------- function VectorInterface.scale(t::AbstractTensorMap, α::Number) T = VectorInterface.promote_scale(t, α) - return scale!(similar(t, T), t, α) + return scale!(zerovector(t, T), t, α) end function VectorInterface.scale!(t::AbstractTensorMap, α::Number) for (c, b) in blocks(t) @@ -74,7 +74,7 @@ function VectorInterface.add!(ty::AbstractTensorMap, tx::AbstractTensorMap, α::Number, β::Number) space(ty) == space(tx) || throw(SpaceMismatch("$(space(ty)) ≠ $(space(tx))")) for ((cy, by), (cx, bx)) in zip(blocks(ty), blocks(tx)) - add!(StridedView(by), StridedView(bx), α, β) + add!(by, bx, α, β) end return ty end diff --git a/test/diagonal.jl b/test/diagonal.jl new file mode 100644 index 000000000..944fecb25 --- /dev/null +++ b/test/diagonal.jl @@ -0,0 +1,219 @@ +diagspacelist = ((ℂ^4)', ℂ[Z2Irrep](0 => 2, 1 => 3), + ℂ[FermionNumber](0 => 2, 1 => 2, -1 => 1), + ℂ[SU2Irrep](0 => 2, 1 => 1)', ℂ[FibonacciAnyon](:I => 2, :τ => 2)) + +@testset "DiagonalTensor with domain $V" for V in diagspacelist + @timedtestset "Basic properties and algebra" begin + for T in (Float32, Float64, ComplexF32, ComplexF64, BigFloat) + t = @constinferred DiagonalTensorMap{T}(undef, V) + t = @constinferred DiagonalTensorMap(rand(T, reduceddim(V)), V) + @test @constinferred(hash(t)) == hash(deepcopy(t)) + @test scalartype(t) == T + @test codomain(t) == ProductSpace(V) + @test domain(t) == ProductSpace(V) + @test space(t) == (V ← V) + @test space(t') == (V ← V) + @test dim(t) == dim(space(t)) + @test isa(@constinferred(norm(t)), real(T)) + @test norm(t)^2 ≈ dot(t, t) + α = rand(T) + @test norm(α * t) ≈ abs(α) * norm(t) + @test norm(t + t, 2) ≈ 2 * norm(t, 2) + @test norm(t + t, 1) ≈ 2 * norm(t, 1) + @test norm(t + t, Inf) ≈ 2 * norm(t, Inf) + p = 3 * rand(Float64) + @test norm(t + t, p) ≈ 2 * norm(t, p) + @test norm(t) ≈ norm(t') + + @test t == @constinferred(TensorMap(t)) + @test norm(t + TensorMap(t)) ≈ 2 * norm(t) + + @test norm(zerovector!(t)) == 0 + @test norm(one!(t)) ≈ sqrt(dim(V)) + @test one!(t) == id(V) + @test norm(one!(t) - id(V)) == 0 + + t1 = DiagonalTensorMap(rand(T, reduceddim(V)), V) + t2 = DiagonalTensorMap(rand(T, reduceddim(V)), V) + t3 = DiagonalTensorMap(rand(T, reduceddim(V)), V) + α = rand(T) + β = rand(T) + @test @constinferred(dot(t1, t2)) ≈ conj(dot(t2, t1)) + @test dot(t2, t1) ≈ conj(dot(t2', t1')) + @test dot(t3, α * t1 + β * t2) ≈ α * dot(t3, t1) + β * dot(t3, t2) + end + end + I = sectortype(V) + @timedtestset "Basic linear algebra: test via conversion" begin + for T in (Float32, ComplexF64) + t1 = DiagonalTensorMap(rand(T, reduceddim(V)), V) + t2 = DiagonalTensorMap(rand(T, reduceddim(V)), V) + @test norm(t1, 2) ≈ norm(convert(TensorMap, t1), 2) + @test dot(t2, t1) ≈ dot(convert(TensorMap, t2), convert(TensorMap, t1)) + α = rand(T) + @test convert(TensorMap, α * t1) ≈ α * convert(TensorMap, t1) + @test convert(TensorMap, t1') ≈ convert(TensorMap, t1)' + @test convert(TensorMap, t1 + t2) ≈ + convert(TensorMap, t1) + convert(TensorMap, t2) + end + end + @timedtestset "Real and imaginary parts" begin + for T in (Float64, ComplexF64, ComplexF32) + t = DiagonalTensorMap(rand(T, reduceddim(V)), V) + + tr = @constinferred real(t) + @test scalartype(tr) <: Real + @test real(convert(TensorMap, t)) == convert(TensorMap, tr) + + ti = @constinferred imag(t) + @test scalartype(ti) <: Real + @test imag(convert(TensorMap, t)) == convert(TensorMap, ti) + + tc = @inferred complex(t) + @test scalartype(tc) <: Complex + @test complex(convert(TensorMap, t)) == convert(TensorMap, tc) + + tc2 = @inferred complex(tr, ti) + @test tc2 ≈ tc + end + end + @timedtestset "Tensor conversion" begin + t = @constinferred DiagonalTensorMap(undef, V) + rand!(t.data) + tc = complex(t) + @test convert(typeof(tc), t) == tc + @test typeof(convert(typeof(tc), t)) == typeof(tc) + end + I = sectortype(V) + if BraidingStyle(I) isa SymmetricBraiding + @timedtestset "Permutations" begin + t = DiagonalTensorMap(randn(ComplexF64, reduceddim(V)), V) + t1 = @constinferred permute(t, $(((2,), (1,)))) + if BraidingStyle(sectortype(V)) isa Bosonic + @test t1 ≈ transpose(t) + end + @test convert(TensorMap, t1) == permute(convert(TensorMap, t), (((2,), (1,)))) + t2 = @constinferred permute(t, $(((1, 2), ()))) + @test convert(TensorMap, t2) == permute(convert(TensorMap, t), (((1, 2), ()))) + t3 = @constinferred permute(t, $(((2, 1), ()))) + @test convert(TensorMap, t3) == permute(convert(TensorMap, t), (((2, 1), ()))) + t4 = @constinferred permute(t, $(((), (1, 2)))) + @test convert(TensorMap, t4) == permute(convert(TensorMap, t), (((), (1, 2)))) + t5 = @constinferred permute(t, $(((), (2, 1)))) + @test convert(TensorMap, t5) == permute(convert(TensorMap, t), (((), (2, 1)))) + end + end + @timedtestset "Trace, Multiplication and inverse" begin + t1 = DiagonalTensorMap(rand(Float64, reduceddim(V)), V) + t2 = DiagonalTensorMap(rand(ComplexF64, reduceddim(V)), V) + @test tr(TensorMap(t1)) == @constinferred tr(t1) + @test tr(TensorMap(t2)) == @constinferred tr(t2) + @test TensorMap(@constinferred t1 * t2) ≈ TensorMap(t1) * TensorMap(t2) + @test TensorMap(@constinferred t1 \ t2) ≈ TensorMap(t1) \ TensorMap(t2) + @test TensorMap(@constinferred t1 / t2) ≈ TensorMap(t1) / TensorMap(t2) + @test TensorMap(@constinferred inv(t1)) ≈ inv(TensorMap(t1)) + @test TensorMap(@constinferred pinv(t1)) ≈ pinv(TensorMap(t1)) + @test all(Base.Fix2(isa, DiagonalTensorMap), + (t1 * t2, t1 \ t2, t1 / t2, inv(t1), pinv(t1))) + + u = randn(Float64, V * V' * V, V) + @test u * t1 ≈ u * TensorMap(t1) + @test u / t1 ≈ u / TensorMap(t1) + @test t1 * u' ≈ TensorMap(t1) * u' + @test t1 \ u' ≈ TensorMap(t1) \ u' + end + @timedtestset "Tensor contraction" begin + d = DiagonalTensorMap(rand(ComplexF64, reduceddim(V)), V) + t = TensorMap(d) + A = randn(ComplexF64, V ⊗ V' ⊗ V, V) + B = randn(ComplexF64, V ⊗ V' ⊗ V, V ⊗ V') + if BraidingStyle(I) isa SymmetricBraiding + @tensor C[a b c; d] := A[a b c; e] * d[e, d] + @test C ≈ A * d + @tensor D[a; b] := d[a, c] * d[c, b] + @test D ≈ d * d + @test D isa DiagonalTensorMap + end + @planar E1[-1 -2 -3; -4 -5] := B[-1 -2 -3; 1 -5] * d[1; -4] + @planar E2[-1 -2 -3; -4 -5] := B[-1 -2 -3; 1 -5] * t[1; -4] + @test E1 ≈ E2 + @planar E1[-1 -2 -3; -4 -5] = B[-1 -2 -3; -4 1] * d'[-5; 1] + @planar E2[-1 -2 -3; -4 -5] = B[-1 -2 -3; -4 1] * t'[-5; 1] + @test E1 ≈ E2 + @planar E1[-1 -2 -3; -4 -5] = B[1 -2 -3; -4 -5] * d[-1; 1] + @planar E2[-1 -2 -3; -4 -5] = B[1 -2 -3; -4 -5] * t[-1; 1] + @test E1 ≈ E2 + @planar E1[-1 -2 -3; -4 -5] = B[-1 1 -3; -4 -5] * d[1; -2] + @planar E2[-1 -2 -3; -4 -5] = B[-1 1 -3; -4 -5] * t[1; -2] + @test E1 ≈ E2 + @planar E1[-1 -2 -3; -4 -5] = B[-1 -2 1; -4 -5] * d'[-3; 1] + @planar E2[-1 -2 -3; -4 -5] = B[-1 -2 1; -4 -5] * t'[-3; 1] + @test E1 ≈ E2 + end + @timedtestset "Factorization" begin + for T in (Float32, ComplexF64) + t = DiagonalTensorMap(rand(T, reduceddim(V)), V) + @testset "eig" begin + D, W = @constinferred eig(t) + @test t * W ≈ W * D + t2 = t + t' + D2, V2 = @constinferred eigh(t2) + VdV2 = V2' * V2 + @test VdV2 ≈ one(VdV2) + @test t2 * V2 ≈ V2 * D2 + end + @testset "leftorth with $alg" for alg in (TensorKit.QR(), TensorKit.QL()) + Q, R = @constinferred leftorth(t; alg=alg) + QdQ = Q' * Q + @test QdQ ≈ one(QdQ) + @test Q * R ≈ t + if alg isa Polar + @test isposdef(R) + end + end + @testset "rightorth with $alg" for alg in (TensorKit.RQ(), TensorKit.LQ()) + L, Q = @constinferred rightorth(t; alg=alg) + QQd = Q * Q' + @test QQd ≈ one(QQd) + @test L * Q ≈ t + if alg isa Polar + @test isposdef(L) + end + end + @testset "tsvd with $alg" for alg in (TensorKit.SVD(), TensorKit.SDD()) + U, S, Vᴴ = @constinferred tsvd(t; alg=alg) + UdU = U' * U + @test UdU ≈ one(UdU) + VdV = Vᴴ * Vᴴ' + @test VdV ≈ one(VdV) + @test U * S * Vᴴ ≈ t + end + end + end + @timedtestset "Tensor functions" begin + for T in (Float64, ComplexF64) + d = DiagonalTensorMap(rand(T, reduceddim(V)), V) + # rand is important for positive numbers in the real case, for log and sqrt + t = TensorMap(d) + @test @constinferred exp(d) ≈ exp(t) + @test @constinferred log(d) ≈ log(t) + @test @constinferred sqrt(d) ≈ sqrt(t) + @test @constinferred sin(d) ≈ sin(t) + @test @constinferred cos(d) ≈ cos(t) + @test @constinferred tan(d) ≈ tan(t) + @test @constinferred cot(d) ≈ cot(t) + @test @constinferred sinh(d) ≈ sinh(t) + @test @constinferred cosh(d) ≈ cosh(t) + @test @constinferred tanh(d) ≈ tanh(t) + @test @constinferred coth(d) ≈ coth(t) + @test @constinferred asin(d) ≈ asin(t) + @test @constinferred acos(d) ≈ acos(t) + @test @constinferred atan(d) ≈ atan(t) + @test @constinferred acot(d) ≈ acot(t) + @test @constinferred asinh(d) ≈ asinh(t) + @test @constinferred acosh(one(d) + d) ≈ acosh(one(t) + t) + @test @constinferred atanh(d) ≈ atanh(t) + @test @constinferred acoth(one(t) + d) ≈ acoth(one(d) + t) + end + end +end diff --git a/test/runtests.jl b/test/runtests.jl index 34037d27c..cb491e1ea 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -57,8 +57,11 @@ Ti = time() include("fusiontrees.jl") include("spaces.jl") include("tensors.jl") +include("diagonal.jl") include("planar.jl") -include("ad.jl") +if !(Sys.isapple()) # TODO: remove once we know why this is so slow on macOS + include("ad.jl") +end include("bugfixes.jl") Tf = time() printstyled("Finished all tests in ", diff --git a/test/spaces.jl b/test/spaces.jl index c9d114728..c3890db0c 100644 --- a/test/spaces.jl +++ b/test/spaces.jl @@ -221,6 +221,7 @@ println("------------------------------------") slist = @constinferred sectors(V) @test @constinferred(TensorKit.hassector(V, first(slist))) @test @constinferred(dim(V)) == sum(dim(s) * dim(V, s) for s in slist) + @test @constinferred(reduceddim(V)) == sum(dim(V, s) for s in slist) @constinferred dim(V, first(slist)) if hasfusiontensor(I) @test @constinferred(TensorKit.axes(V)) == Base.OneTo(dim(V))