diff --git a/.github/workflows/PR_comment.yml b/.github/workflows/PR_comment.yml new file mode 100644 index 0000000..ea00f15 --- /dev/null +++ b/.github/workflows/PR_comment.yml @@ -0,0 +1,17 @@ +name: Docs preview comment +on: + pull_request: + types: [labeled] + +permissions: + pull-requests: write +jobs: + pr_comment: + runs-on: ubuntu-latest + steps: + - name: Create PR comment + if: github.event_name == 'pull_request' && github.repository == github.event.pull_request.head.repo.full_name && github.event.label.name == 'documentation' + uses: thollander/actions-comment-pull-request@v3 + with: + message: 'After the build completes, the updated documentation will be available [here](https://quantumkithub.github.io/MultiTensorKit.jl/previews/PR${{ github.event.number }}/)' + comment-tag: 'preview-doc' \ No newline at end of file diff --git a/.gitignore b/.gitignore index 10593a9..5f9fc64 100644 --- a/.gitignore +++ b/.gitignore @@ -11,4 +11,3 @@ Manifest.toml benchmark/*.json docs/Manifest.toml docs/build/ -docs/src/index.md diff --git a/Project.toml b/Project.toml index a542085..04d74c0 100644 --- a/Project.toml +++ b/Project.toml @@ -4,25 +4,30 @@ authors = ["Boris De Vos, Laurens Lootens and Lukas Devos"] version = "0.1.0" [deps] -Artifacts = "56f22d72-fd6d-98f1-02f0-08ddc0907c33" -JSON3 = "0f8b85d8-7281-11e9-16c2-39a750bddbf1" +DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab" +Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" +TensorKit = "07d1fe3e-3e46-537d-9eac-e9e13d0d4cec" TensorKitSectors = "13a9c161-d5da-41f0-bcbd-e1a08ae0647f" [compat] Aqua = "0.8.9" -Artifacts = "1.10, 1" -JSON3 = "1.14.1" +DelimitedFiles = "1.9" +Pkg = "1" +Random = "1" SafeTestsets = "0.1" -TensorKitSectors = "0.1.2" +TensorKit = "0.16" +TensorKitSectors = "0.3" Test = "1.10" TestExtras = "0.3" julia = "1.10" [extras] Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" TestExtras = "5ed8adda-3752-4e41-b88a-e8b09835ee3a" [targets] -test = ["Aqua", "SafeTestsets", "Test", "TestExtras"] +test = ["Aqua", "LinearAlgebra", "Random", "SafeTestsets", "Test", "TestExtras"] diff --git a/docs/Project.toml b/docs/Project.toml new file mode 100644 index 0000000..a36ea53 --- /dev/null +++ b/docs/Project.toml @@ -0,0 +1,4 @@ +[deps] +Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244" +MultiTensorKit = "f0555a46-f681-4ef1-bed5-d64870568636" diff --git a/docs/make.jl b/docs/make.jl new file mode 100644 index 0000000..8356949 --- /dev/null +++ b/docs/make.jl @@ -0,0 +1,27 @@ +using Documenter +using DocumenterCitations +# using TensorKitSectors, TensorKit +using MultiTensorKit + +pages = ["Home" => "index.md", + "Manual" => ["man/fusioncats.md", "man/multifusioncats.md", "man/extension.md", "man/implementation.md"], + "Library" => "lib/library.md", + "References" => "references.md"] + +# bibliography +bibpath = joinpath(@__DIR__, "src", "assets", "MTKrefs.bib") +bib = CitationBibliography(bibpath; style=:authoryear) + +mathengine = MathJax3(Dict(:loader => Dict("load" => ["[tex]/physics"]), + :tex => Dict("inlineMath" => [["\$", "\$"], ["\\(", "\\)"]], + "tags" => "ams", + "packages" => ["base", "ams", "autoload", "physics", + "mathtools"]))) + +makedocs(; sitename="MultiTensorKit.jl", modules=[MultiTensorKit], + authors="Boris De Vos, Laurens Lootens and Lukas Devos", + pages=pages, pagesonly=true, plugins=[bib], + format=Documenter.HTML(; prettyurls=true, mathengine=mathengine, + assets=["assets/custom.css"])) + +deploydocs(; repo="https://github.com/QuantumKitHub/MultiTensorKit.jl", push_preview = true) \ No newline at end of file diff --git a/docs/src/assets/MTKrefs.bib b/docs/src/assets/MTKrefs.bib new file mode 100644 index 0000000..7a95ca2 --- /dev/null +++ b/docs/src/assets/MTKrefs.bib @@ -0,0 +1,53 @@ +@book{etingof2016tensor, + title={Tensor Categories}, + author={Etingof, P. and Gelaki, S. and Nikshych, D. and Ostrik, V.}, + isbn={9781470434410}, + lccn={2015006773}, + series={Mathematical Surveys and Monographs}, + url={https://books.google.be/books?id=Z6XLDAAAQBAJ}, + year={2016}, + publisher={American Mathematical Society} +} + +@article{Lootens_2023, + title={Dualities in One-Dimensional Quantum Lattice Models: Symmetric Hamiltonians and Matrix Product Operator Intertwiners}, + volume={4}, + ISSN={2691-3399}, + url={http://dx.doi.org/10.1103/PRXQuantum.4.020357}, + DOI={10.1103/prxquantum.4.020357}, + number={2}, + journal={PRX Quantum}, + publisher={American Physical Society (APS)}, + author={Lootens, Laurens and Delcamp, Clement and Ortiz, Gerardo and Verstraete, Frank}, + year={2023}, + month=jun } + +@misc{etingof2009, + title={Fusion categories and homotopy theory}, + author={Pavel Etingof and Dmitri Nikshych and Victor Ostrik and with an appendix by Ehud Meir}, + year={2009}, + eprint={0909.3140}, + archivePrefix={arXiv}, + primaryClass={math.QA}, + url={https://arxiv.org/abs/0909.3140}, +} + +@misc{henriques2020, + title={Representations of fusion categories and their commutants}, + author={André Henriques and David Penneys}, + year={2020}, + eprint={2004.08271}, + archivePrefix={arXiv}, + primaryClass={math.OA}, + url={https://arxiv.org/abs/2004.08271}, +} + +@misc{Lootens_2024, + title={Entanglement and the density matrix renormalisation group in the generalised Landau paradigm}, + author={Laurens Lootens and Clement Delcamp and Frank Verstraete}, + year={2024}, + eprint={2408.06334}, + archivePrefix={arXiv}, + primaryClass={quant-ph}, + url={https://arxiv.org/abs/2408.06334}, +} \ No newline at end of file diff --git a/docs/src/assets/custom.css b/docs/src/assets/custom.css new file mode 100644 index 0000000..4e38a20 --- /dev/null +++ b/docs/src/assets/custom.css @@ -0,0 +1,47 @@ +.center { + display: block; + margin-left: auto; + margin-right: auto; +} + +/* set fixed non-trivial inversion and hue rotation */ +:root { + --inversionFraction: 100%; + --hueRotation: -180deg; +} + +/* color-invert images */ +.color-invertible { + filter: invert(var(--inversionFraction)) hue-rotate(var(--hueRotation)) !important; + -ms-filter: invert(var(--inversionFraction)) !important; + -webkit-filter: invert(var(--inversionFraction)) hue-rotate(var(--hueRotation)) !important; +} +/* including the logo when we make one +.docs-logo { + filter: invert(var(--inversionFraction)) hue-rotate(var(--hueRotation)) !important; + -ms-filter: invert(var(--inversionFraction)) !important; + -webkit-filter: invert(var(--inversionFraction)) hue-rotate(var(--hueRotation)) !important; +} */ + +/* but disable for the two light themes */ +.theme--documenter-light .color-invertible { + filter: invert(0%) hue-rotate(0deg) !important; + -ms-filter: invert(var(0%)) !important; + -webkit-filter: invert(var(0%)) hue-rotate(0deg) !important; +} +.theme--catppuccin-latte .color-invertible { + filter: invert(0%) hue-rotate(0deg) !important; + -ms-filter: invert(var(0%)) !important; + -webkit-filter: invert(var(0%)) hue-rotate(0deg) !important; +} +/* for the logo as well */ +/* .theme--documenter-light .docs-logo { + filter: invert(0%) hue-rotate(0deg) !important; + -ms-filter: invert(var(0%)) !important; + -webkit-filter: invert(var(0%)) hue-rotate(0deg) !important; +} +.theme--catppuccin-latte .docs-logo { + filter: invert(0%) hue-rotate(0deg) !important; + -ms-filter: invert(var(0%)) !important; + -webkit-filter: invert(var(0%)) hue-rotate(0deg) !important; +} */ \ No newline at end of file diff --git a/docs/src/index.md b/docs/src/index.md new file mode 100644 index 0000000..d112e44 --- /dev/null +++ b/docs/src/index.md @@ -0,0 +1,41 @@ +# MultiTensorKit + +**TensorKit extension to multifusion categories** + +## Package summary +MultiTensorKit.jl provides the user a package to work with multifusion categories, the extension of regular fusion categories where the unit is no longer simple and unique. +Multifusion categories naturally embed the structure of module categories over fusion categories. +Hence, MultiTensorKit.jl allows not only the fusion of objects within the same +fusion category (as TensorKit.jl), but also the fusion with and between module categories over these fusion categories. + +MultiTensorKit.jl is built to be compatible with TensorKit, thus allowing the construction of symmetric tensors with new symmetries due to the module structure. +Through this, +tensor network simulations of quantum many-body systems with aid of [MPSKit.jl](https://github.com/QuantumKitHub/MPSKit.jl) can be performed. + +## Table of contents + +```@contents +Pages = ["man/fusioncats.md", "man/multifusioncats.md", "man/extension.md", "man/implementation.md", "lib/library.md", "references.md"] +Depth = 2 +``` + +## Installation + +MultiTensorKit.jl is currently not registered to the Julia General Registry. +You can install the package as +``` +pkg> add https://github.com/QuantumKitHub/MultiTensorKit.jl.git +``` + +## Usage + +As the name suggests, MultiTensorKit is an extension of [TensorKit.jl](https://github.com/Jutho/TensorKit.jl) and +[TensorKitSectors.jl](https://github.com/QuantumKitHub/TensorKitSectors.jl). +Therefore, we recommend including TensorKit +to your project. +Additionally, MultiTensorKit was made to be functional with [MPSKit.jl](https://github.com/QuantumKitHub/MPSKit.jl) +and [MPSKitModels.jl](https://github.com/QuantumKitHub/MPSKitModels.jl) for Matrix Product State (MPS) calculations, supporting symmetries +which go beyond TensorKit. +In particular, TensorKit v0.16.0 and MPSKit v0.13.9 contain the necessary functionality to deal with the multifusion categorical structure provided by MultiTensorKit. +All these packages are registered in JuliaRegistries and can be added through the package manager. + diff --git a/docs/src/lib/library.md b/docs/src/lib/library.md new file mode 100644 index 0000000..86a8f6c --- /dev/null +++ b/docs/src/lib/library.md @@ -0,0 +1,5 @@ +# Library + +```@autodocs +Modules = [MultiTensorKit] +``` \ No newline at end of file diff --git a/docs/src/man/extension.md b/docs/src/man/extension.md new file mode 100644 index 0000000..5b89c22 --- /dev/null +++ b/docs/src/man/extension.md @@ -0,0 +1,96 @@ +# [MultiTensorKit as an extension to TensorKit](@id extension) + +This section will explain the internal changes to TensorKit which are required to extend the compatibility with fusion categories to multifusion ones. +Users who are unfamiliar with TensorKit are kindly guided towards the [TensorKit tutorial](https://jutho.github.io/TensorKit.jl/stable/man/tutorial/). + +## As a `Sector` +MultiTensorKit is at its core an extension to [TensorKitSectors](https://github.com/QuantumKitHub/TensorKitSectors.jl), as it simply provides a new `Sector`, i.e. simple objects of a multifusion category, named `BimoduleSector`. + +````julia +struct BimoduleSector{Name} <: Sector + i::Int + j::Int + label::Int +end +```` + +`i` and `j` specify which subcategory $\mathcal{C}_{ij}$ we are considering, and `label` selects a particular simple object within that subcategory. +`Name` selects which multifusion category to work with, and is a `Symbol`. +As of now, only `:A4` is available, referring to the multifusion category consisting of $\mathsf{Rep(A_4)}$ as the largest fusion category, and all its Morita dual fusion categories and corresponding bimodule categories. +The fusion of these `BimoduleSector`s is then defined by the fusion rules of the multifusion category. +In particular, +````julia +a = BimoduleSector{:A4}(i, j, label1) +b = BimoduleSector{:A4}(k, l, label2) +a ⊗ b # empty unless j == k +```` + +The fusion rules are read in via the artifact labeled "fusiondata", and are extracted at runtime when calling the fusion of two `BimoduleSector`s (or `Nsymbol`) for the first time. +These data are cached in a hash table for later use. + +A consequence of the multifusion structure is the colorings used in the graphical calculus of fusions. +A natural introduction is the notion of a *left* and *right* unit of some (simple) object in the multifusion category. +These are called with the functions `TensorKitSectors.leftunit` and `TensorKitSectors.rightunit`. +Clearly, for the usual case of just one fusion category, these both coincide with the unique unit object. +For this reason, the `TensorKitSectors.UnitStyle` trait is introduced to differentiate between fusion categories and multifusion categories, returning `SimpleUnit()` and `GenericUnit()`, respectively. +This trait is used to specialise certain functions within TensorKit (see below). +`TensorKitSectors.allunits` is also introduced to return all the simple unit objects of the (multi)fusion category. +This function must be defined for every multifusion category to be compatible with TensorKit. + +Via the fusion rules, the left and right units of all `BimoduleSector`s along with their duals are also extracted and cached. +The fusion behavior that must be imposed for the entire `BimoduleSector` is the most general one that appears in the fusion between any two `BimoduleSector`s. +In particular, if one of the diagonal fusion categories has fusion rules with multiplicities, then the entire `BimoduleSector` must be set to have `TensorKitSectors.FusionStyle(::Type{<:BimoduleSector}) = GenericFusion()`. +This is the case for the $A_4$ `BimoduleSector`, since $\mathsf{Rep(A_4)}$ is one of the diagonal fusion categories and has fusion rules with multiplicities. + +Going beyond the ring structure, the F-symbols are also read in from the artifact, and are stored in a hash table for later use. +The F-symbols are then used to perform F-moves on `BimoduleSector`s, which is required to perform recouplings of fusion trees when doing e.g. contractions of tensors with these categories grading the vector spaces. + +Due to the nature of the multifusion category, it is currently not possible to define a non-trivial braiding on the `BimoduleSector`s, so these are set to be non-braided: `TensorKitSectors.BraidingStyle(::Type{<:BimoduleSector}) = NoBraiding()`. +This is especially important when working with matrix product states, as all algorithms are required to remain planar, since no (half-)braiding is available to perform crossings of fusion trees. + +## As a symmetry in TensorKit +Since `BimoduleSector`s are `Sector`s, they can be used as symmetries in TensorKit. +This way, we can construct symmetric tensors with the symmetries of the multifusion category, which are more general than those of the fusion categories. +In particular, the vector spaces graded by these `BimoduleSector`s are not only graded by the simple objects of the fusion categories, but also by the simple objects of the bimodule categories. +This allows for more general tensor network simulations of quantum many-body systems with symmetries which go beyond those of fusion categories. + +Certain changes within TensorKit were required to make it compatible with the multifusion categorical structure. +In particular, the presence of a simple unit object for every fusion category on the diagonal of the multifusion category, along with the off-diagonal nature of the simple objects of the bimodule categories, required some internal changes to the way unit objects were treated in TensorKit. +Most notably, the unit object is no longer unique, and thus it is of utmost importance that the correct one is considered when performing tensor contractions at the level of the fusion trees. +This is achieved precisely through colorings and the use of `leftunit` and `rightunit`. +For this reason, every fusion tree manipulation which previously involved "the" unit object, now involves the `leftunit` and `rightunit` of some neighboring sector in the manipulation to identify the correct color. +An important example of this is explained in the previous section on [Opposite module categories](@ref opposite_module_categories), namely the mapping of a splitting vertex to a fusion vertex via a B-move. + +When manipulating spaces graded by `BimoduleSector`s, one needs to also be careful of which unit spaces can compose with other graded spaces on which side. +This introduces the functions `TensorKit.leftunitspace` and `TensorKit.rightunitspace`, which check whether the coloring of the (in general) composite object grading the space is uniform, then return the one-dimensional space with the unique left/right unit object consistent with that coloring. +Clearly, `leftunitspace` and `rightunitspace` coincide for fusion categories; this defaults to `TensorKit.unitspace`. +For multifusion categories, the latter function returns the space with the semisimple unit object grading it. +`leftunitspace` and `rightunitspace` are used with the functions `TensorKit.insertleftunit` and `TensorKit.insertrightunit`. +There, based on the index where a unit space is wished to be inserted, a `leftunitspace` or `rightunitspace` is added such that the resulting space remains consistent with the fusion rules of the multifusion category. +Similarly, `TensorKit.removeunit` is used to remove unit spaces, and will explicitly check whether the space contains only unit objects of any color with `TensorKit.isunitspace`. + +## MultiTensorKit compatibility with MPSKit + +This section will briefly explain the changes within MPSKit which are required to make it compatible with MultiTensorKit. +For a more practical explanation, users are kindly guided towards the [Example section](@ref implementation). + +The main change within MPSKit is very similar to the fusion tree manipulations in TensorKit, namely the use of `leftunit` and `rightunit` to identify the correct unit object. +In the case of MPSKit, trivial spaces are used everywhere, from the boundary of a finite matrix product state to the virtual spaces of a Hamiltonian written as an matrix product operator. +Additionally, multiple tensor contractions made use of braiding tensors to perform crossings of legs of the MPSs/MPOs. +Since no (half-)braiding is available for the `BimoduleSector`s, all algorithms had to be made planar, and thus all braiding tensors were removed. +Since all braidings were trivial, this was dealt with by simply removing the braiding tensors and replacing the crossing of legs with a termination and reintroduction of the legs without crossing by aid of `removeunit`, `insertleftunit` and `insertrightunit`. + +At the level of the MPS, the correct unit object can be identified through the use of `leftunit(ψ::AbstractMPS)` and `rightunit(ψ::AbstractMPS)`. +When the virtual space of the MPS is graded by a diagonal `BimoduleSector`, i.ea unitary fusion category, then these all coincide with the unique unit object of that fusion category. +However, when the virtual space is graded by an off-diagonal `BimoduleSector`, i.e. a bimodule category, then the left and right units are different, and thus it is important to identify the correct one when performing MPS algorithms. +For example, Hamiltonians should always have the same coloring as the right unit of the MPS, since they are contracted at the physical level of the MPS. +Similarly, excitations of an MPS are labeled by `BimoduleSector`s with the same coloring as the left unit of the MPS, since the auxiliary charge leg of the excitation is attached on the other side of the MPS to the virtual level. +This is illustrated in the following figure. + +```@raw html + +``` + +Here, $\mathcal{D}$ is the fusion category grading the physical space of the MPS, $\mathcal{R}$ is the right $\mathcal{D}$-module category grading the virtual space of the MPS. +$\mathcal{E}$ is the Morita dual fusion category of $\mathcal{D}$ with respect to $\mathcal{R}$, which is the category labelling the excitations of the MPS. +In this case, the left unit of the MPS is in $\mathcal{E}$, and the right unit in $\mathcal{D}$. \ No newline at end of file diff --git a/docs/src/man/fusioncats.md b/docs/src/man/fusioncats.md new file mode 100644 index 0000000..e38a766 --- /dev/null +++ b/docs/src/man/fusioncats.md @@ -0,0 +1,76 @@ +The manual has been divided into different sections in an attempt to break down the information the user requires to use MultiTensorKit.jl. +We start off with a short summary of fusion category theory. +Users familiar with TensorKit.jl may have read the [Optional introduction to category theory](https://jutho.github.io/TensorKit.jl/stable/man/categories/) in the documentation of TensorKit; this section can then largely be skipped. +Be aware that notation may differ from the literature. + +Afterwards, the extension to multifusion categories is explained, and its relation to (bi)module categories over fusion categories is shown. + +# [Fusion category theory](@id fusion_cat_theory) + +The aim of this section is to explain the bare minimum required to proceed to the next section on multifusion category theory and bimodule categories. +More details can be found in the [TensorKit](https://jutho.github.io/TensorKit.jl/stable/man/categories/) documentation or the book Tensor Categories [etingof2016tensor](@cite). + +Let us start simple and introduce the **fusion ring** $\mathcal{C}$ in a black-box manner. +This ring +* consists of finitely many simple objects $\{ X_1, X_2, ..., X_R \}$, with $R$ the rank of the fusion ring, +* which can be fused with one another: $X_i \otimes X_j \cong \sum_k N_{ij}^k X_k,$ introducing the **N-symbol** $N_{ij}^k \in \mathbb{N}$ in the fusion rules, +* contains a *unique* unit object $1_\mathcal{C}$ which satisfies $1_\mathcal{C} \otimes X_i \cong X_i \otimes 1_\mathcal{C} \cong X_i$ for all objects $X_i \in \mathcal{C}$, +* has a dual object $\overline{X}$ for every object $X$ such that $X \otimes \overline{X} \cong \overline{X} \otimes X \cong 1_\mathcal{C} \oplus ...$, generalising the notion of the inverse element. + +To extend the fusion ring to the **fusion category**, we need to add the following structure: +* Consider only the representatives of isomorphism classes of simple objects $\mathcal{I}_\mathcal{C}$, +* The associator $F^{X_iX_jX_k}: (X_i \otimes X_j) \otimes X_k \xrightarrow{\sim} X_i \otimes (X_j \otimes X_k)$ + which fulfills the famous pentagon equation, +* Morphisms between (simple) objects $\text{Hom}_\mathcal{C}(X_i, X_j)$, which are empty vector spaces unless the objects are isomorphic, the latter then giving $\mathbb{C}$, +* More general morphisms $\text{Hom}_\mathcal{C}(X_i \otimes X_j, X_k) \cong \mathbb{C}^{N_{ij}^k}$. + +This way, we can describe fusion categories by a triple $(\otimes, 1_\mathcal{C}, F)$ of $\mathcal{C}$ defining its monoidal product, unit object and monoidal associator, the latter also commonly called the **F-symbol**. +In particular, the simple objects have their respective quantum dimensions $d_i = \dim(X_i)$ which form their own one-dimensional representation of the fusion algebra: $d_i d_j = \sum_k N_{ij}^k d_k$. +In particular, the unit object always has quantum dimension 1, and all other quantum dimensions are larger or equal to one. +These quantum dimensions are encoded in the F-symbol. +The isomorphisms instead of the equalities are a technical detail, so we drop that notation. + +Vectors in these hom-spaces are graphically denoted as living in the trivalent junction +```@raw html + +``` + + +With the F-symbol, we can perform F-moves: + +```@raw html + +``` + +TensorKit requires the F-symbols to be unitary. +This way, we can interpret the F-symbol $F^{ijk}_l$ as a unitary matrix, and the F-move as a unitary basis transformation. +Unitarity is also useful from a diagrammatic point of view because the category is then equipped with a pivotal and spherical structure. +This essentially means that morphisms can be drawn and moved around freely on a 2-sphere, such that vector spaces can be moved freely from domain (codomain) to codomain (domain). + +## Examples + +### $\mathsf{Vec_G}$ and $\mathsf{Rep(G)}$ +Colloquially speaking, category theory attempts to generalise mathematical structures and their relations in a way that different structures can be treated in an equal manner. +This is noted in particular as fusion category theory encompasses not only finite and compact groups, but also their representations. +We show a table sketching how these are put on equal footing categorically. + +|$\mathsf{Vec_G}$|$\mathsf{Rep(G)}$|Categorical generalisation| +|:---:|:---:|:---:| +|$G$-graded vector spaces $V_1, V_2, ...$|Representations of $G$ $(V_1, \pi_1), (V_2, \pi_2), ...$|Objects| +|$G$-graded preserving linear maps $\phi: V \rightarrow W$|Intertwiners $f: V_1 \rightarrow V_2$, $f \circ \pi_1 = \pi_2 \circ f$|Morphisms $\text{Hom}_\mathcal{C}$| +1d subspaces $\mathbb{C}_{g_1}, \mathbb{C}_{g_2}$: $\text{Hom}_{\mathsf{Vec_G}}(\mathbb{C}_{g_1},\mathbb{C}_{g_2}) = \delta_{g_1g_2}$|Irreps: $\text{Hom}_{\mathsf{Rep(G)}}(\rho_i,\rho_j) = \delta_{ij} \mathbb{C}$ (Schur)|Simple objects: $\text{Hom}_{\mathcal{C}}(a,b) = \delta_{ab}\mathbb{C}$| +$G$-graded tensor product $(V \otimes W)_g = \oplus_{hk=g} V_h \otimes W_k$| $\pi_i \otimes \pi_j \simeq \oplus_i N_{ij}^k\rho_k$ | Direct sum, monoidal product, fusion rules, multiplicity| +$\mathbb{C}_1 \otimes W \simeq W \simeq W \otimes \mathbb{C}_1$ | Trivial rep 1: $1 \otimes \rho = \rho = \rho \otimes 1$ | Monoidal unit $1_\mathcal{C}$ +$\mathbb{C}_g \otimes \mathbb{C}_{g^{-1}} = \mathbb{C}_1 = \mathbb{C}_{g^{-1}} \otimes \mathbb{C}_g$ | $\pi \otimes \overline{\pi} = 1 \oplus ...$ | Dual object +$F:(V \otimes W) \otimes U \xrightarrow{\sim}V \otimes (W \otimes U)$|$F: (\pi_1 \otimes \pi_2) \otimes \pi_3 \xrightarrow{\sim} \pi_1 \otimes (\pi_2 \otimes \pi_3)$|F-symbol| + +### $\mathsf{Fib}$ and $\mathsf{Ising}$ +Arguably the simplest fusion category besides the familiar groups or representations of groups is the Fibonacci fusion category. +This contain 2 simple objects $1$ and $\tau$, with non-trivial fusion rule $\tau \otimes \tau = 1 \oplus \tau$. +This fusion category is in fact **braided** as well, and actually **modular**. + +Another simple fusion category is the Ising category, commonly denoted $\mathsf{Ising}$. +3 simple objects form this category, namely $\{1, \psi, \sigma\}$, where $1$ and $\psi$ behave like the trivial charged representations of $\mathbb{Z}_2$, while $\sigma$ is the $\mathbb{Z}_2$ extension of this. +The fusion rules reflect these: $1 \otimes \psi = \psi = \psi \otimes 1, \sigma \otimes X = \sigma = X \otimes \sigma$ for $X = 1, \psi$, and $\sigma \otimes \sigma = 1 \oplus \psi$. +This fusion category is also **modular**. +Both these modular fusion categories are already implemented in [TensorKitSectors](https://github.com/QuantumKitHub/TensorKitSectors.jl). \ No newline at end of file diff --git a/docs/src/man/img/A4_sym_entanglement_spectrum.svg b/docs/src/man/img/A4_sym_entanglement_spectrum.svg new file mode 100644 index 0000000..7f1de43 --- /dev/null +++ b/docs/src/man/img/A4_sym_entanglement_spectrum.svg @@ -0,0 +1,554 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/docs/src/man/img/Bmove_MF.svg b/docs/src/man/img/Bmove_MF.svg new file mode 100644 index 0000000..d7b2c30 --- /dev/null +++ b/docs/src/man/img/Bmove_MF.svg @@ -0,0 +1,915 @@ + + + + diff --git a/docs/src/man/img/Bmove_lit.svg b/docs/src/man/img/Bmove_lit.svg new file mode 100644 index 0000000..7f2051e --- /dev/null +++ b/docs/src/man/img/Bmove_lit.svg @@ -0,0 +1,731 @@ + + + + diff --git a/docs/src/man/img/Fmove.svg b/docs/src/man/img/Fmove.svg new file mode 100644 index 0000000..455b368 --- /dev/null +++ b/docs/src/man/img/Fmove.svg @@ -0,0 +1,1074 @@ + + + + diff --git a/docs/src/man/img/Fmove_CM.svg b/docs/src/man/img/Fmove_CM.svg new file mode 100644 index 0000000..fc5cc83 --- /dev/null +++ b/docs/src/man/img/Fmove_CM.svg @@ -0,0 +1,663 @@ + + + + diff --git a/docs/src/man/img/Fmove_CMD.svg b/docs/src/man/img/Fmove_CMD.svg new file mode 100644 index 0000000..2db2f12 --- /dev/null +++ b/docs/src/man/img/Fmove_CMD.svg @@ -0,0 +1,688 @@ + + + + diff --git a/docs/src/man/img/Fmove_D.svg b/docs/src/man/img/Fmove_D.svg new file mode 100644 index 0000000..a66cdac --- /dev/null +++ b/docs/src/man/img/Fmove_D.svg @@ -0,0 +1,659 @@ + + + + diff --git a/docs/src/man/img/Fmove_MD.svg b/docs/src/man/img/Fmove_MD.svg new file mode 100644 index 0000000..d2c7733 --- /dev/null +++ b/docs/src/man/img/Fmove_MD.svg @@ -0,0 +1,673 @@ + + + + diff --git a/docs/src/man/img/Fmove_coloring.svg b/docs/src/man/img/Fmove_coloring.svg new file mode 100644 index 0000000..a78e3b3 --- /dev/null +++ b/docs/src/man/img/Fmove_coloring.svg @@ -0,0 +1,2294 @@ + + + + diff --git a/docs/src/man/img/Nsymbol_coloring.svg b/docs/src/man/img/Nsymbol_coloring.svg new file mode 100644 index 0000000..616bfa2 --- /dev/null +++ b/docs/src/man/img/Nsymbol_coloring.svg @@ -0,0 +1,1198 @@ + + + + diff --git a/docs/src/man/img/anyonchain_excitation.svg b/docs/src/man/img/anyonchain_excitation.svg new file mode 100644 index 0000000..241b1be --- /dev/null +++ b/docs/src/man/img/anyonchain_excitation.svg @@ -0,0 +1,214 @@ + + + + diff --git a/docs/src/man/img/homvector.svg b/docs/src/man/img/homvector.svg new file mode 100644 index 0000000..d11f63e --- /dev/null +++ b/docs/src/man/img/homvector.svg @@ -0,0 +1,345 @@ + + + + diff --git a/docs/src/man/img/pentagon_colored.svg b/docs/src/man/img/pentagon_colored.svg new file mode 100644 index 0000000..1d008a9 --- /dev/null +++ b/docs/src/man/img/pentagon_colored.svg @@ -0,0 +1,7977 @@ + + + + diff --git a/docs/src/man/img/qdim_fs_MF.svg b/docs/src/man/img/qdim_fs_MF.svg new file mode 100644 index 0000000..38d13b7 --- /dev/null +++ b/docs/src/man/img/qdim_fs_MF.svg @@ -0,0 +1,382 @@ + + + + diff --git a/docs/src/man/implementation.md b/docs/src/man/implementation.md new file mode 100644 index 0000000..47c1381 --- /dev/null +++ b/docs/src/man/implementation.md @@ -0,0 +1,201 @@ +# [Symmetric tensor networks: $\mathsf{Rep(A_4)}$ as a guiding example](@id implementation) +This tutorial is dedicated to explaining how MultiTensorKit was implemented to be compatible with with TensorKit and MPSKit for matrix product state simulations. +This section should be largely self-contained and is meant to be a practical guide to using MultiTensorKit. +In particular, we will be making a generalised anyonic spin chain. +We will demonstrate how to reproduce the entanglement spectra found in [Lootens_2024](@cite). +The model considered there is a spin-1 Heisenberg model with additional terms to break the usual $\mathsf{U_1}$ symmetry to $\mathsf{Rep(A_4)}$, while having a non-trivial phase diagram and relatively easy Hamiltonian to write down. + +This will be done with the `A4Object = BimoduleSector{:A4}` `Sector`, which is the multifusion category which contains the structure of the module categories over $\mathsf{Rep(A_4)}$. +Since there are 7 module categories, `A4Object` is a $r=7$ multifusion category. +There are 3 fusion categories up to equivalence: +- ``\mathsf{Vec_{A_4}}``: the category of $\mathsf{A_4}$-graded vector spaces. The group $\mathsf{A}_4$ is order $4!/2 = 12$ and has thus 12 simple objects. +- ``\mathsf{Rep(A_4)}``: the irreducible representations of the group $\mathsf{A}_4$, of which there are 4. One is the trivial representation, two are one-dimensional non-trivial and the last is three-dimensional. +- ``\mathsf{Rep(H)}``: the representation category of some Hopf algebra which does not have a name. It has 6 simple objects. + +For this example, we will require the following packages: +````julia +using TensorKit, MultiTensorKit, MPSKit, MPSKitModels, Plots +```` + +## Identifying the simple objects +We first need to select which fusion category we wish to use to grade the physical Hilbert space, and which fusion category to represent e.g. the symmetry category. +In our case, we are interested in selecting $\mathcal{D} = \mathsf{Rep(A_4)}$ for the physical Hilbert space. +We know the module categories over $\mathsf{Rep(G)}$ to be $\mathsf{Rep^\psi(H)}$ for a subgroup $\mathsf{H}$ and 2-cocycle $\psi$. +Thus, the 7 module categories $\mathcal{M}$ one can choose over $\mathsf{Rep(A_4)}$ are + +- ``\mathsf{Rep(A_4)}`` itself as the regular module category, +- ``\mathsf{Vec}``: the category of vector spaces, +- ``\mathsf{Rep(\mathbb{Z}_2)}``, +- ``\mathsf{Rep(\mathbb{Z}_3)}``, +- ``\mathsf{Rep(\mathbb{Z}_2 \times \mathbb{Z}_2)}``, +- ``\mathsf{Rep^\psi(\mathbb{Z}_2 \times \mathbb{Z}_2)}``, +- ``\mathsf{Rep^\psi(A_4)}``. + +When referring to specific fusion and module categories, we will use this non-multifusion notation. + +Now that we have identified the fusion and module categories, we want to select the relevant objects we wish to place in our graded spaces. +Unfortunately, due to the nature of how the N-symbol and F-symbol data are generated, the objects of the fusion subcategories are not ordered such that `label = 1` corresponds to the unit object. +Hence, the simplest way to find the unit object of a fusion subcategory is + +````julia +unit(A4Object(i, i, 1)) +```` + +Left and right units of subcategories are uniquely specified by their fusion rules. +For example, the left unit of a subcategory $\mathcal{C}_{ij}$ is the simple object in $\mathcal{C}_i$ for which + +$$^{}_a \mathbb{1} \times a = a \quad \forall a \in \mathcal{C}_i.$$ + +A similar condition uniquely defines the right unit of a subcategory. +For fusion subcategories, a necessary condition is that the left and right units coincide. + +Identifying the other simple objects of a (not necessarily fusion) category requires more work. +We recommend a combination of the following to uniquely determine any simple object `a`: +- Check the dimension of the simple object: `dim(a)` +- Check the dual of the simple object: `dual(a)` +- Check how this simple object fuses with other (simple) objects: `Nsymbol(a, b, c)` + +The dual object of some simple object $a$ of an arbitrary subcategory $\mathcal{C}_{ij}$ is defined as the unique object $a^* \in \mathcal{C}_{ji}$ satisfying + +$$^{}_a \mathbb{1} \in a \times a^* \quad \text{and} \quad \mathbb{1}_a \in a^* \times a$$ + +with multiplicity 1. +## Constructing the Hamiltonian and matrix product state +TensorKit has been made compatible with the multifusion structure. +With this, we can make `GradedSpace`s whose objects are in `A4Object`: + +````julia +D1 = A4Object(6, 6, 1) # unit in this case +D2 = A4Object(6, 6, 2) # non-trivial 1d irrep +D3 = A4Object(6, 6, 3) # non-trivial 1d irrep +D4 = A4Object(6, 6, 4) # 3d irrep +```` +Since we want to replicate a spin-1 Heisenberg model, it makes sense to use the 3-dimensional irrep to grade the physical space, and thus construct our Hamiltonian. +We don't illustrate here how to derive the considered Hamiltonian in a $\mathsf{Rep(A_4)}$ basis, but simply give it. +For now, we construct a finite-size spin chain; below we will repeat the calculation for an infinite system. + +````julia +P = Vect[A4Object](D4 => 1) # physical space +T = ComplexF64 +# usual Heisenberg part +h1_L = zeros(T, P ⊗ P ← P) +h1_R = zeros(T, P ← P ⊗ P) +block(h1_L, D4) .= [0; 1] +block(h1_R, D4) .= [0 1;] +@planar h1[-1 -2; -3 -4] := h1_L[-1 1; -3] * h1_R[-2; 1 -4] + +# biquadratic term +h2_L = zeros(T, P ⊗ Vect[A4Object](D1 => 1, D2 => 1, D3 => 1) ← P) +h2_R = ones(T, P ← Vect[A4Object](D1 => 1, D2 => 1, D3 => 1) ⊗ P) +block(h2_L, D4) .= [4 / 3; 1 / 3; 1 / 3] +@planar h2[-1 -2; -3 -4] := h2_L[-1 1; -3] * h2_R[-2; 1 -4] + +# anti-commutation term +h3_L = zeros(T, P ⊗ P ← P) +h3_R = zeros(T, P ← P ⊗ P) +block(h3_L, D4) .= [1; 0] +block(h3_R, D4) .= [0 1;] +@planar h3[-1 -2; -3 -4] := h3_L[-1 1; -3] * h3_R[-2; 1 -4] + +L = 60 +J1 = -2.0 # probing the A4 symmetric phase first +J2 = -5.0 +lattice = FiniteChain(L) +H1 = @mpoham sum(-2 * h1{i,j} for (i, j) in nearest_neighbours(lattice)) +H2 = @mpoham sum(h2{i,j} for (i, j) in nearest_neighbours(lattice)) +H3 = @mpoham sum(2im * h3{i,j} for (i, j) in nearest_neighbours(lattice)) + +H = H1 + J1 * H2 + J2 * H3 +```` + +For the matrix product state, we will select $\mathsf{Vec}$ as the module category for now: +````julia +M = A4Object(1, 6, 1) # Vec +```` +and construct the finite MPS: +````julia +D = 40 # bond dimension +V = Vect[A4Object](M => D) +Vb = Vect[A4Object](M => 1) # non-degenerate boundary virtual space +init_mps = FiniteMPS(L, P, V; left = Vb, right = Vb) +```` + +!!! warning "Important" + We must pass on a left and right virtual space to the keyword arguments `left` and `right` of the `FiniteMPS` constructor, since these would by default try to place the unit space of the `Sector`, which for any `BimoduleSector` consists of the semisimple unit. As of now, this will error at the level of the N-symbol call as the colors of the various components do not all match. + +## Two-site DMRG and the entanglement spectrum +We can now look to find the ground state of the Hamiltonian with two-site DMRG. +We use this instead of the "usual" one-site DMRG because the two-site algorithm will smartly fill up the blocks of the local tensor during the sweep, allowing one to initialise as a product state in one block and more likely avoid local minima, a common occurence in symmetric tensor network simulations. +````julia +dmrg2alg = DMRG2(; verbosity = 2, tol = 1e-7, trscheme = trunctol(; atol = 1e-4 / sqrt(dim(M)))) +ψ, _ = find_groundstate(init_mps, H, dmrg2alg) +```` +The truncation scheme keyword argument is mandatory when calling `DMRG2` in MPSKit. +Here, we choose to truncate such that all singular values are larger than $10^{-4}$, while setting the default tolerance for convergence to $10^{-7}$. +The extra factor is required to cancel the normalisation of the singular values by the quantum dimensions. +More information on this can be found in the [MPSKit](https://github.com/QuantumKitHub/MPSKit.jl) documentation. +To run one-site DMRG anyway, use `DMRG` which does not require a truncation scheme. + +Now that we've found the ground state, we can compute the entanglement spectrum in the middle of the chain. +````julia +entanglement_spectrum(ψ, round(Int, L/2)) +```` +This returns a `TensorKit.SectorVector` containing the singular values. +This can be sliced to find the singular values corresponding to a particular simple object. +More information on that can be found in the docstring of `MPSKit.entanglement_spectrum`. +In this case, all singular values correspond to $\mathsf{Vec}$. +We can also immediately return a plot of this data by the following: +````julia +entanglementplot(ψ; site = round(Int, L/2)) +```` +This plot will show the singular values per simple object, as well as include the "effective" bond dimension, which is simply the dimension of the virtual space where we cut the system. +The next section will show this plot along with those when selecting the other module categories. + +## Search for the correct dual model + +Consider a quantum lattice model with its symmetries determining the phase diagram. +For every phase in the phase diagram, the dual model for which the ground state maximally breaks all symmetries spontaneously is the one where the entanglement is minimised and the tensor network is represented most efficiently [Lootens_2024](@cite). +Let us confirm this result, starting with the $A_4$-symmetric phase whose ideal dual model is in the $\mathsf{Rep(A_4)}$ spontaneous symmetry breaking phase. +The code will look exactly the same as above, except the virtual space of the MPS will change to be graded by the other module categories: + +````julia +module_numlabels(i) = MultiTensorKit._numlabels(A4Object, i, 6) +V = Vect[A4Object]((i, 6, label) => D for label in 1:module_numlabels(i)) +Vb = Vect[A4Object](first(sectors(V)) => 1) # not all charges on boundary, play around with what is there +```` + +```@raw html + +``` + +The plot shows the entanglement spectra of the various dual models in the middle of the ground state for the original Hamiltonian probing the $A_4$-symmetric phase, along with the ground state without symmetries. +Every subtitle also mentions the memory required to store the same middle ground state tensor. +This plot should be compared to [Lootens_2024; Figure 2](@cite). + +This can be repeated with other parameter values for $J_1$ and $J_2$ in the Hamiltonian to probe the $\mathbb{Z}_2 \times \mathbb{Z}_2$-symmetric or $A_4$ SPT phase. + +!!! note "Additional functions and keyword arguments" + Certain commonly used functions within MPSKit can pass extra keyword arguments compatible with multifusion MPS simulations. + In particular, the keyword argument `sector` (note the lowercase "s") appears in + - `excitations` with `QuasiparticleAnsatz`: the sector is selected by adding an auxiliary space to the *domain* of each eigenvector of the transfer matrix. Since in a full contraction the domain of the eigenvector lies in the opposite side of the physical space (labeled by objects in $\mathcal{D} = \mathsf{Rep(A_4)}$), the excitations lie in the Morita dual category $\mathcal{E} = \mathcal{D^*_M}$. This will default to the unit of $\mathcal{E}$ if not specified. + - `exact_diagonalization`: the `sector` keyword argument now requires a simple object in $\mathcal{D}$, since this is the fusion category which specifies the bond algebra from which the Hamiltonian is constructed. This is equivalent to adding a charged leg on the leftmost (or rightmost) virtual space of the MPS in conventional MPS cases. This will default to the unit of $\mathcal{D}$ if not specified. + +## Differences with the infinite case +We can repeat the above calcalations also for an infinite system. +The `lattice` variable will change, as well as the MPS constructor and the algorithm: +````julia +lattice = InfiniteChain(1) +init = InfiniteMPS([P], [V]) +inf_alg = VUMPS(; verbosity = 2, tol = 1e-7) +```` + +Besides `VUMPS`, `IDMRG` and `IDMRG2` are as easy to run with the `A4Object` `BimoduleSector`. +It is also clear that boundary terms do not play a role in this case. + +!!! note "More functions for infinite systems" + When dealing with an infinite system, additional information can be retrieved from the matrix product state. + These also require a keyword argument `sector` to be specified, as otherwise the correct unit object is placed by default. + These are + - `transfer_spectrum`: similar to `excitations`, the (partial) transfer matrix spectrum is selected by adding a charged auxiliary space to the transfer matrix eigenvectors. + - `correlation_length`: since this function calls `transfer_spectrum`, the same logic applies. + - `excitations` with `QuasiparticleAnsatz`. \ No newline at end of file diff --git a/docs/src/man/multifusioncats.md b/docs/src/man/multifusioncats.md new file mode 100644 index 0000000..c5d8122 --- /dev/null +++ b/docs/src/man/multifusioncats.md @@ -0,0 +1,279 @@ +# [Extending to multifusion category theory](@id multifusion_cat_theory) + +This section will explain how to go from a fusion category to a multifusion category, as well as why one would want to consider the latter. +Multifusion categories naturally embed the structure of **bimodule categories**. +To explain this, we must start by explaining module categories over fusion categories, following this up with (invertible) bimodule categories, and finishing off with the multifusion structure. + +## Module categories +We start from a fusion category $\mathcal{D}$ defined by the triple $(\otimes_\mathcal{D}, \mathbb{1}_\mathcal{D}, {}^\mathcal{D}\!F)$ with simple objects $\alpha, \beta, ... \in \mathcal{I}_\mathcal{D}$. +We drop the $\mathcal{D}$ subscript when there is no ambiguity concerning the fusion category. +We call its associator + +$${}^\mathcal{D}\!F^{\alpha \beta \gamma}: \alpha \otimes (\beta \otimes \gamma) \rightarrow (\alpha \otimes \beta) \otimes \gamma$$ + + the **monoidal associator**. + An F-move is now graphically portrayed as: + +```@raw html + +``` + +We can consider the **right module category** $\mathcal{M}$ over $\mathcal{D}$, which is a category (not necessarily fusion!) with (isomorphism classes of) simple objects $\mathcal{I}_\mathcal{M} = \{A,B,...\}$, a right action + +$$\triangleleft: \mathcal{M} \times \mathcal{D} \rightarrow \mathcal{M}$$ + +and the **right module associator** + +$${}^\triangleleft\!F^{A\alpha\beta}: A \triangleleft (\alpha \otimes \beta) \rightarrow (A \triangleleft \alpha) \triangleleft \beta.$$ + +An F-move with this module associator can be expressed as: + +```@raw html + +``` + +The module structure of $\mathcal{M}$ is now defined as the triple $(\mathcal{M}, \triangleleft, {}^\triangleleft\!F)$. +The right module associator now satisfies a mixed pentagon equation with ${}^\mathcal{D}\!F$. + +Similarly, we can define a **left module category** $(\mathcal{M}, \triangleright, {}^\triangleright\!F)$ over a fusion category $(\otimes_\mathcal{C}, 1_\mathcal{C}, {}^\mathcal{C}\!F)$. +The functor is now a left action of $\mathcal{C}$ on $\mathcal{M}$ + +$$\triangleright: \mathcal{C} \times \mathcal{M} \rightarrow \mathcal{M},$$ + +while the **left module associator** is a natural isomorphism that acts as + +$${}^\triangleright\!F^{abA}: (a \otimes b) \triangleright A \rightarrow a \triangleright (b \triangleright A)$$ + +for $\mathcal{I}_\mathcal{C} = \{a,b,...\}$. +The left module associator also fulfills a mixed pentagon equation with ${}^\mathcal{C}\!F$. +An F-move with ${}^\mathcal{C}\!F$ takes on the form: + +```@raw html + +``` + +We can combine the concepts of left and right module categories as follows. +Say there are two fusion categories $\mathcal{C}$ and $\mathcal{D}$. +A $(\mathcal{C}, \mathcal{D})$-**bimodule category** is a module category, now defined through a sextuple $(\mathcal{M}, \triangleright, \triangleleft, {}^\triangleright\!F, {}^\triangleleft\!F, {}^{\bowtie}\!F)$ such that $(\mathcal{M}, \triangleright, {}^\triangleright\!F)$ is a left $\mathcal{C}$-module category, $(\mathcal{M}, \triangleleft, {}^\triangleleft\!F)$ is a right $\mathcal{D}$-module category, and with additional structure such that the **bimodule associator** acts as + +$${}^{\bowtie}\!F^{aA\alpha}: (a \triangleright A) \triangleleft \alpha \rightarrow a \triangleright (A \triangleleft \alpha)$$ + +for $a \in \mathcal{I}_\mathcal{C}, \alpha \in \mathcal{I}_\mathcal{D}, A \in \mathcal{I}_\mathcal{M}$. +The bimodule associator fulfills a mixed pentagon equation with the module associators. +An F-move with ${}^{\bowtie}\!F$ is given by: + +```@raw html + +``` +## [Opposite module categories](@id opposite_module_categories) +Consider a fusion category $\mathcal{D}$ and a *right* module category $\mathcal{M}$ over $\mathcal{D}$. +We can define $\mathcal{M}^{\text{op}}$ to be the opposite category of $\mathcal{M}$ [etingof2009](@cite). +Then, $\mathcal{M}^{\text{op}}$ is a *left* module category over $\mathcal{D}$. +A similar statement can be made starting from a left module category and getting an opposite right module category. +In particular, given a $(\mathcal{C}, \mathcal{D})$-bimodule category $\mathcal{M}$ over the fusion categories $\mathcal{C}$ and $\mathcal{D}$, we see immediately that $\mathcal{M}^{\text{op}}$ is a $(\mathcal{D}, \mathcal{C})$-bimodule category. + +Interestingly, due to the opposite actions + +$$\triangleleft^{\text{op}}: \mathcal{D} \times \mathcal{M}^{\text{op}} \rightarrow \mathcal{M}^{\text{op}}$$ + +and + +$$\triangleright^{\text{op}}: \mathcal{M}^{\text{op}} \times \mathcal{C} \rightarrow \mathcal{M}^{\text{op}},$$ + +there is a valid notion of multiplying a module category with its opposite: + +$$\mathcal{M} \times \mathcal{M}^{\text{op}} \rightarrow \mathcal{C}, \quad \mathcal{M}^{\text{op}} \times \mathcal{M} \rightarrow \mathcal{D}.$$ + +## Invertible bimodule categories and Morita equivalence + +The bimodule categories we consider are restricted to be **invertible**. +This means that this bimodule category fulfills the condition $\mathcal{C} \equiv \mathcal{D}^*_\mathcal{M}$, or in other words the two fusion categories $\mathcal{C}$ and $\mathcal{D}$ are **Morita equivalent** (or each other's **Morita dual** with respect to $\mathcal{M}$). +Morita equivalence between $\mathcal{C}$ and $\mathcal{D}$ can be shown to hold if and only if the bimodule category fulfills + +$$\mathcal{M} \boxtimes_\mathcal{D} \mathcal{M}^\text{op} \simeq \mathcal{C}, \quad \text{and} \quad \mathcal{M}^\text{op} \boxtimes_\mathcal{C} \mathcal{M} \simeq \mathcal{D},$$ + +where $\boxtimes_\mathcal{C}$ and $\boxtimes_\mathcal{D}$ denote the Deligne product relative to $\mathcal{C}$ and $\mathcal{D}$ respectively. +This is precisely the invertibility property of the bimodule category. + +## Multifusion categories +A fusion category has the important condition that $\mathsf{End}_\mathcal{C}(1_\mathcal{C}) \cong \mathbb{C}$, i.e., the unit of the fusion category is simple. +If we drop this condition, then we consider a **multifusion category**. +We assume the multifusion category itself to be indecomposable, meaning it is not the direct sum of two non-trivial multifusion categories. +Let us call this multifusion category $\mathcal{C}$. +It will be clear in a moment that this will not be ambiguous. +Its unit can then be decomposed as + +$$1_\mathcal{C} = \bigoplus_{i=1}^r 1_r,$$ + +i.e., it is decomposed into simple unit objects of the subcategories $\mathcal{C}_{ij} \coloneqq 1_i \otimes \mathcal{C} \otimes 1_j$. +With this, we see that the multifusion category itself can be decomposed into its subcategories + +$$\mathcal{C} = \bigoplus_{i,j=1}^r \mathcal{C}_{ij}.$$ + +We call this an $r \times r$ multifusion category. + +We want to consider multifusion categories because **their structure encapsulates that of (bi)module categories**. +Every diagonal category $\mathcal{C}_{ii} \coloneqq \mathcal{C}_i$ (also known as a component category) is a fusion category, and every off-diagonal category $\mathcal{C}_{ij}$ is an invertible $(\mathcal{C}_{i}, \mathcal{C}_{j})$-bimodule category. +That way, as long as we know how the simple objects of the fusion and module categories fuse with one another, and we can determine all the monoidal and module associators, we can treat the multifusion category as one large fusion category with limited fusion rules. +In particular, the tensor product + +$$\otimes_\mathcal{C}: \mathcal{C}_{ij} \times \mathcal{C}_{kl} \rightarrow \delta_{jk}\mathcal{C}_{il}$$ + +takes on the same structure as the product of two matrix units. +This is not a coincidence; there is a deep relation between multifusion categories and matrix algebras [etingof2016tensor; Section 4.3](@cite). + +Given a subcategory $\mathcal{C}_{ij}$, we can define the **left** unit as the unit object of the fusion category $\mathcal{C}_i$, while the **right** unit is the unit object of the fusion category $\mathcal{C}_j$. +In other words, the left unit of $\mathcal{C}_{ij}$ is the unique object of the multifusion category $\mathcal{C}$ for which + +$$1_i \otimes_\mathcal{C} m_{ij} = m_{ij} \quad \forall m_{ij} \in \mathcal{C}_{ij},$$ + +and similarly for the right unit of $\mathcal{C}_{ij}$, + +$$m_{ij} \otimes_\mathcal{C} 1_j = m_{ij} \quad \forall m_{ij} \in \mathcal{C}_{ij}.$$ + +We can also immediately see that for a bimodule subcategory $\mathcal{C}_{ij}$, the opposite (bi)module subcategory $\mathcal{C}_{ij}^{\text{op}} \equiv \mathcal{C}_{ji}$, and as expected, + +$$\mathcal{C}_{ij} \times \mathcal{C}_{ji} \rightarrow \mathcal{C}_i,$$ + +just like what we concluded when considering opposite module categories outside of the multifusion structure. +### 2-category and coloring +Multifusion categories can also be interpreted as 2-categories. +We still interpret the objects of this 2-category the same way. +The 1-morphisms are the subcategories themselves, and the 2-morphisms the morphisms of the multifusion category. +The graphical calculus of monoidal 1-categories can be extended to 2-categories by use of *colorings* [henriques2020](@cite). +We have previously differed between module strands and fusion strands by the color of the strand itself. +However, in 2-categories the strands (1-morphisms) separate regions which are colored based on the fusion categories they are representing. +Since we draw the strands vertically, a single strand results in a left and right region, and the colorings will determine the fusion category which fuses from the left or right with that single strand. +In particular, fusion strands necessarily have the same coloring on the left and right, while module strands have a mismatching coloring. + +The simplest non-trivial fusion diagram is a trivalent junction: + +```@raw html + +``` + +The most general case is the top left figure, where all three regions have a different coloring. +The top middle region having the same coloring from the top left and top right strands follow from the delta function in the tensor product definition. +However, this most general trivalent junction with three colorings will never be needed. +In short, we will always be considering a single bimodule category $\mathcal{C}_{ij}$ at a time, and the only other non-diagonal subcategory which fuses with this is its opposite $\mathcal{C}_{ji}$. +This is displayed in the top middle and right. +Similarly, two colorings are required when considering the fusion between a fusion and module strand, shown in the bottom left and middle figure. +The simplest trivalent junctions boil down to fusions within fusion categories, which is obviously drawn with just one color. +This is shown in the bottom right. + +With this coloring system, we can specify which associator must be called to perform a particular F-move. +Such an F-move would look like + +```@raw html + +``` + +### Why opposite module categories end up being necessary in MultiTensorKit + +One of the common manipulations which can act on a tensor map is the transposition of vector spaces. +We will refer to this as the bending of legs. +One of the elementary bends is the **right bend**, where one of the tensor legs is bent along the right from the codomain to the domain, or vice versa. +At the level of the tensor, a covariant index becomes contravariant, or vice versa. +Similarly, a **left bend** can also be performed, bending the leg along the left. +This guarantees that legs will not cross, preventing braidings which require extra data known as R-symbols. + +Linear algebra tells us that given a (finite-dimensional) vector space $V$ with a basis denoted $\{|e_i\rangle\}$, one can consider the **dual** vector space $V^*$, whose dual basis $\{\langle e_i^*|\}$ satisfies the property $\langle e_j^* | e_i \rangle = \delta_{ij}$. +In the diagrammatic calculus, specifying whether a tensor map leg represents a vector space or its dual is done with an arrow. +Following the TensorKit convention, legs with arrows pointing downwards are vector spaces, and arrows pointing upwards state that we are considering its dual. +In particular, at the level of fusion trees we can also draw arrows on the strands to denote whether we are considering morphisms between objects or dual objects. + +In principle, choosing to bend, e.g., codomain legs to the right and domain legs to left is an arbitrary choice, but would require to distinguish between left and right transposes. +However, TensorKit.jl is implemented in a way that does not differentiate the two. +In particular, we do not have to worry about this when considering categorical symmetries where, in principle, the left and right dual of an object are not equivalent. +This is because this left-right symmetry is guaranteed when considering unitary fusion categories, which is what is done in TensorKit and necessarily in MultiTensorKit. + +For this reason, at the level of the fusion trees the topological move that is performed to bend the legs along the right is called a **B-move**. +Graphically, one can show that this bend boils down to a particular F-move. +The typical equation found in the literature is the following: +```@raw html + +``` + +The reason to only consider B-moves is rooted in the choice of canonical form of fusion trees within TensorKit, where fusions are iterated over from left to right and splittings from right to left. + +Importantly, we identify the dual vector space labeled by a module category with a vector space labeled by the opposite module category. +Consequently, + +$$\mathcal{M}^* \simeq \mathcal{M}^\text{op}.$$ + +In the multifusion setting, this can also be seen graphically. +By keeping track of the colorings and the directions of the arrows of the legs, one can see that we need to slightly modify the expression for the B-move to the following: + +```@raw html + +``` + +where by $\mathbb{1}_k$ we mean the unit of $\mathcal{C}_k$. +Note how we keep the color subscripts for the dual object, but in reality they switch since the arrow direction is reversed. + +### More on the topological data: gauge choices and distilling properties of the subcategories + +An important property of the F-symbols is that they must satisfy the **triangle identities**. +In fusion category theory, this states that isomorphisms between (simple) objects $a$ and the tensor product between $a$ and the unit $\mathbb{1}$ exists, and that + +$$(a \otimes \mathbb{1}) \otimes b \cong a \otimes (\mathbb{1} \otimes b)$$ + +for $b$ in the same fusion category. +This can be straightforwardly generalised to multifusion categories. +This requires a particular partial gauging of these trivalent vertices. + +Besides the triangle identities, the (multi)fusion category must also fulfill the **pentagon equations**. +These encapsulate the two identical manners to evaluating the fusion of four objects in the (multi)fusion category. +Every fusion category's F-symbols must satisfy these individually, but also the (bi)module associators between bimodule and fusion categories. +One can check that, for every pair of fusion categories, their bimodule category and opposite bimodule category, there are 32 pentagon equations to be satisfied. +In the multifusion notation, they can be represented by + +```@raw html + +``` + +The most generic F-move contains 4 colors. +For that reason, MultiTensorKit requires the F-symbol data to be provided as some data file (currently .txt) with 4 + 6 + 4 + 2 = 16 columns. +The first 4 refer to the colors, the following 6 label the simple objects of the corresponding subcategories, the next 4 label multiplicities, and the final 2 provide the real and imaginary value of the F-symbol itself. + +In a similar manner, the N-symbols contain maximally three colors, so these data must provide 3 columns labeling the colors, 3 columns labeling the simple objects and a final column with the dimension of the corresponding vector space. + +Besides the B-move (and closely related A-move, which we do not illustrate), we can also see how the quantum dimension and Frobenius-Schur indicator expressions get modified. +We already know that an F-move of the form $F^{c \bar{c} c}_{c}$ needs to be evaluated for these topological data, which is in fact another gauge fixing. +Graphically, we find + +```@raw html + +``` +In principle, a gauge fixing can be done to set the Frobenius-Schur indicator to $\pm 1$. +However, this assumption is no longer required within TensorKit and can be relaxed to just be a phase. + +The above gauge fixing is a property of the unitary gauge. +By choosing appropriate bases, one can transform the F-symbols of a (multi)fusion category to be unitary matrices. +More details on the importance of unitary topological data can be found in the [TensorKit](https://jutho.github.io/TensorKit.jl/stable/man/categories/#ss_topologicalfusion) documentation. + +### Braiding +Sectors provided by MultiTensorKit do *not* support braiding. +We do this for two reasons. +On the one hand, there is no natural 1-categorical way of defining braidings between the components of the multifusion category. +It is possible that the diagonal fusion categories themselves are braided, but a "componentwise" braiding is unwise to support. +On the other hand, it is entirely possible to write matrix product state manipulations in a planar manner (which has been done in [MPSKit](https://github.com/QuantumKitHub/MPSKit.jl)), thus avoiding the need of a braiding tensor. More information on this can be found in [MultiTensorKit as an extension to TensorKit](@ref extension). + +## Examples of multifusion categories +Without specifying any of the categories, the simplest non-trivial multifusion category is a $2\times 2$ one, and the categories can be organised in a matrix as + +$$\mathcal{C} = \begin{pmatrix} \mathcal{C}_1 & \mathcal{M} \\ \mathcal{M}^{\text{op}} & \mathcal{C}_2\end{pmatrix}.$$ + +We already identified the off-diagonal elements with module categories over the fusion categories on the diagonal. +Accordingly, $\mathcal{M}$ is a $(\mathcal{C}_1, \mathcal{C}_2)$-bimodule category, and $\mathcal{M}^{\text{op}}$ is the opposite module category and a $(\mathcal{C}_2, \mathcal{C}_1)$-bimodule category. + +If we take $\mathcal{C}_1 = \mathcal{C}_2 = \mathsf{Rep}(\mathbb{Z}_2)$ and $\mathcal{M} = \mathsf{Vec}$, then the entire multifusion category is Morita equivalent to $\mathsf{Rep}(\mathbb{Z}_2)$, and we view the $\mathbb{Z}_2$-extension of $\mathsf{Rep}(\mathbb{Z}_2)$ to be precisely the $\mathsf{Ising}$ category [etingof2009; Section 9](@cite) [etingof2016tensor; Example 4.10.5](@cite). +We identify the trivial representation of $\mathsf{Rep}(\mathbb{Z}_2)$ with the unit of $\mathsf{Ising}$, the sign representation with $\psi$ and the unique object of $\mathsf{Vec}$ with the duality object $\sigma$. +One can easily check that the fusion rules of $\mathsf{Ising}$ match with those we expect within $\mathsf{Rep}(\mathbb{Z}_2)$ and with its module category $\mathsf{Vec}$. +Additionally, the fusion between $\mathsf{Vec}$ and $\mathsf{Vec}^\text{op}$ (and vice-versa) giving every object in $\mathcal{C}_1$ ($\mathcal{C}_2$) is consistent with $\sigma \times \sigma^* = 1 + \psi$. +This particular example can be found in [TensorKitSectors](https://github.com/QuantumKitHub/TensorKitSectors.jl). + +This construction can be generalised to $\mathcal{C}_1 = \mathcal{C}_2 = \mathsf{Rep(G)}$ with $\mathsf{G}$ a finite abelian group, such that the entire multifusion category is Morita equivalent to $\mathsf{Rep(G)}$ and can be evaluated as the Tambara-Yamagami category $\mathsf{TY}(\mathsf{G})$ (with positive Frobenius-Schur indicator for our purposes), and $\mathsf{Vec}$ will represent the duality object which squares to all invertible objects of the original group. +To be exact, one of the diagonal fusion categories should be $\mathsf{Vec_G}$ for the correct Morita dual relation, but it is known for abelian groups that this is isomorphic to $\mathsf{Rep(G)}$. + +The example $\mathcal{C}_1 = \mathsf{Rep(A_4)}$ is worked out more in detail in the [guiding example](@ref implementation). \ No newline at end of file diff --git a/docs/src/references.md b/docs/src/references.md new file mode 100644 index 0000000..4b47677 --- /dev/null +++ b/docs/src/references.md @@ -0,0 +1,4 @@ +# References + +```@bibliography +``` \ No newline at end of file diff --git a/src/MultiTensorKit.jl b/src/MultiTensorKit.jl index 34d66bc..635ac69 100644 --- a/src/MultiTensorKit.jl +++ b/src/MultiTensorKit.jl @@ -2,10 +2,14 @@ module MultiTensorKit export BimoduleSector, A4Object -using JSON3 -using Artifacts +using DelimitedFiles +using Pkg +using Pkg.Artifacts using TensorKitSectors +using TensorKit +import TensorKit: hasblock, dim + include("bimodulesector.jl") end diff --git a/src/bimodulesector.jl b/src/bimodulesector.jl index fd52d31..4091914 100644 --- a/src/bimodulesector.jl +++ b/src/bimodulesector.jl @@ -1,80 +1,129 @@ +""" + struct BimoduleSector{Name} <: Sector + BimoduleSector{Name}(i::Int, j::Int, label::Int) + +Represents objects in the component subcategory ``𝒞ᵢⱼ`` of the multifusion category ``𝒞 = ⨁ᵢⱼ 𝒞ᵢⱼ``, +where ``𝒞`` is identified as Name. + +## Fields +- `i::Int`: The row index of the object in the matrix representation of the multifusion category. +- `j::Int`: The column index of the object in the matrix representation of the multifusion category. +- `label::Int`: The label of the object within the component subcategory. +""" +const ImplementedBimoduleSectors = [:A4] + struct BimoduleSector{Name} <: Sector i::Int j::Int label::Int + function BimoduleSector{Name}(i::Int, j::Int, label::Int) where {Name} + Name ∈ ImplementedBimoduleSectors || + throw(ArgumentError("BimoduleSector $Name not implemented")) + i <= size(BimoduleSector{Name}) && j <= size(BimoduleSector{Name}) || + throw(DomainError("object outside the matrix $Name")) + return label <= _numlabels(BimoduleSector{Name}, i, j) ? new{Name}(i, j, label) : + throw(DomainError("label outside category $Name($i, $j)")) + end end -BimoduleSector{Name}(data::NTuple{3,Int}) where {Name} = BimoduleSector{Name}(data...) + +BimoduleSector{Name}(data::NTuple{3, Int}) where {Name} = BimoduleSector{Name}(data...) +BimoduleSectorName(::Type{BimoduleSector{Name}}) where {Name} = Name const A4Object = BimoduleSector{:A4} +Base.convert(::Type{<:BimoduleSector{Name}}, labels::NTuple{3, Int}) where {Name} = BimoduleSector{Name}(labels...) + +function Base.show(io::IO, a::BimoduleSector{Name}) where {Name} + if get(io, :typeinfo, nothing) === typeof(a) + print(io, (a.i, a.j, a.label)) + else + print(io, typeof(a), (a.i, a.j, a.label)) + end + return nothing +end + # Utility implementations # ----------------------- -function Base.isless(a::I, b::I) where {I<:BimoduleSector} +function Base.isless(a::I, b::I) where {I <: BimoduleSector} return isless((a.i, a.j, a.label), (b.i, b.j, b.label)) end Base.hash(a::BimoduleSector, h::UInt) = hash(a.i, hash(a.j, hash(a.label, h))) -function Base.convert(::Type{BimoduleSector{Name}}, d::NTuple{3,Int}) where {Name} +function Base.convert(::Type{BimoduleSector{Name}}, d::NTuple{3, Int}) where {Name} return BimoduleSector{Name}(d...) end -Base.IteratorSize(::Type{SectorValues{<:BimoduleSector}}) = Base.SizeUnknown() +Base.size(::Type{A4Object}) = 7 + +Base.IteratorSize(::Type{<:SectorValues{<:BimoduleSector}}) = Base.SizeUnknown() -# TODO: generalize? -# TODO: numlabels? -function Base.iterate(iter::SectorValues{A4Object}, (I, label)=(1, 1)) - I > 12 * 12 && return nothing - i, j = CartesianIndices((12, 12))[I].I - maxlabel = numlabels(A4Object, i, j) +function Base.iterate(iter::SectorValues{<:BimoduleSector}, (I, label) = (1, 1)) + A = eltype(iter) + s = size(A) + I > s * s && return nothing + i, j = CartesianIndices((s, s))[I].I + maxlabel = _numlabels(A, i, j) return if label > maxlabel iterate(iter, (I + 1, 1)) else - A4Object(i, j, label), (I, label + 1) + A(i, j, label), (I, label + 1) end end +function Base.length(::SectorValues{I}) where {I <: BimoduleSector} + s = size(I) + return sum(_numlabels(I, i, j) for i in 1:s, j in 1:s) +end + TensorKitSectors.FusionStyle(::Type{A4Object}) = GenericFusion() -TensorKitSectors.BraidingStyle(::Type{A4Object}) = NoBraiding() +TensorKitSectors.BraidingStyle(::Type{<:BimoduleSector}) = NoBraiding() +TensorKitSectors.sectorscalartype(::Type{A4Object}) = ComplexF64 -function TensorKitSectors.:⊗(a::A4Object, b::A4Object) +function TensorKitSectors.:⊗(a::I, b::I) where {I <: BimoduleSector} @assert a.j == b.i - Ncache = _get_Ncache(A4Object)[a.i, a.j, b.j] - return A4Object[A4Object(a.i, b.j, c_l) - for (a_l, b_l, c_l) in keys(Ncache) - if (a_l == a.label && b_l == b.label)] + Ncache = _get_Ncache(I)[a.i, a.j, b.j] + return I[ + I(a.i, b.j, c_l) for (a_l, b_l, c_l) in keys(Ncache) + if (a_l == a.label && b_l == b.label) + ] end -function _numlabels(::Type{A4Object}, i, j) - return Ncache = _get_Ncache(A4Object) +function _numlabels(::Type{T}, i, j) where {T <: BimoduleSector} + return length(_get_dual_cache(T)[2][i, j]) end +# User-friendly functions +# ------------------- +#TODO: add functions to identify categories + # Data from files # --------------- const artifact_path = joinpath(artifact"fusiondata", "MultiTensorKit.jl-data-v0.1.5") -function extract_Nsymbol(::Type{A4Object}) - filename = joinpath(artifact_path, "A4", "Nsymbol.json") - isfile(filename) || throw(LoadError(filename, 0, "Nsymbol file not found for $Name")) - json_string = read(filename, String) - Narray = copy(JSON3.read(json_string)) - return map(reshape(Narray, 12, 12, 12)) do x - y = Dict{NTuple{3,Int},Int}() - for (k, v) in x - a, b, c = parse.(Int, split(string(k)[2:(end - 1)], ", ")) - y[(a, b, c)] = v - end - return y +function extract_Nsymbol(::Type{I}) where {I <: BimoduleSector} + name = string(BimoduleSectorName(I)) + filename = joinpath(artifact_path, name, "Nsymbol.txt") + isfile(filename) || throw(LoadError(filename, 0, "Nsymbol file not found for $name")) + Narray = readdlm(filename) # matrix with 7 columns + + data_dict = Dict{NTuple{3, Int}, Dict{NTuple{3, Int}, Int}}() + for row in eachrow(Narray) + i, j, k, a, b, c, N = Int.(@view(row[1:size(I)])) + colordict = get!(data_dict, (i, j, k), Dict{NTuple{3, Int}, Int}()) + push!(colordict, (a, b, c) => N) end + + return data_dict end -const Ncache = IdDict{Type{<:BimoduleSector},Array{Dict{NTuple{3,Int},Int},3}}() +const Ncache = IdDict{Type{<:BimoduleSector}, Dict{NTuple{3, Int}, Dict{NTuple{3, Int}, Int}}}() -function _get_Ncache(::Type{T}) where {T<:BimoduleSector} +function _get_Ncache(::Type{T}) where {T <: BimoduleSector} global Ncache return get!(Ncache, T) do return extract_Nsymbol(T) end end -function TensorKitSectors.Nsymbol(a::I, b::I, c::I) where {I<:A4Object} +function TensorKitSectors.Nsymbol(a::I, b::I, c::I) where {I <: BimoduleSector} # TODO: should this error or return 0? (a.j == b.i && a.i == c.i && b.j == c.j) || throw(ArgumentError("invalid fusion channel")) @@ -82,107 +131,196 @@ function TensorKitSectors.Nsymbol(a::I, b::I, c::I) where {I<:A4Object} return get(_get_Ncache(I)[i, j, k], (a.label, b.label, c.label), 0) end -# TODO: can we define dual for modules? -const Dualcache = IdDict{Type{<:BimoduleSector},Vector{Tuple{Int,Vector{Int}}}}() +const Dualcache = IdDict{Type{<:BimoduleSector}, Tuple{Vector{Int64}, Matrix{Vector{Int64}}}}() -function _get_dual_cache(::Type{T}) where {T<:BimoduleSector} +function _get_dual_cache(::Type{T}) where {T <: BimoduleSector} global Dualcache return get!(Dualcache, T) do return extract_dual(T) end end -function extract_dual(::Type{A4Object}) - N = _get_Ncache(A4Object) - ncats = size(N, 1) - return map(1:ncats) do i +function extract_dual(::Type{I}) where {I <: BimoduleSector} + N = _get_Ncache(I) + ncats = size(I) + Is = zeros(Int, ncats) + + map(1:ncats) do i Niii = N[i, i, i] - nobj = maximum(first, keys(Niii)) - - # find identity object: - # I x I -> I, a x I -> a, I x a -> a - I = findfirst(1:nobj) do u - get(Niii, (u, u, u), 0) == 1 || return false - for j in 1:nobj - get(Niii, (j, u, j), 0) == 1 || return false - get(Niii, (u, j, j), 0) == 1 || return false + nobji = maximum(first, keys(N[i, i, i])) + # want to return a leftunit and rightunit for each entry in multifusion cat + # leftunit/rightunit needs to at least be the unit object within a fusion cat + Is[i] = findfirst(1:nobji) do a + get(Niii, (a, a, a), 0) == 1 || return false # I x I -> I + for othera in 1:nobji + get(Niii, (othera, a, othera), 0) == 1 || return false # a x I -> a + get(Niii, (a, othera, othera), 0) == 1 || return false # I x a -> a + end + + # check leftunit + map(1:ncats) do j + nobjj = maximum(first, keys(N[j, j, j])) + for b in 1:nobjj + get(N[i, j, j], (a, b, b), 0) == 1 || return false # I = leftunit(b) + end + end + + # check rightunit + map(1:ncats) do k + nobjk = maximum(first, keys(N[k, k, k])) + for c in 1:nobjk + get(N[k, i, k], (c, a, c), 0) == 1 || return false # I = rightunit(c) + end end return true end + end + + allduals = Matrix{Vector{Int}}(undef, ncats, ncats) # ncats square matrix of vectors + for i in 1:ncats + nobji = maximum(first, keys(N[i, i, i])) + for j in 1:ncats + allduals[i, j] = Int[] - # find duals - # a x abar -> I - duals = map(1:nobj) do j - return findfirst(1:nobj) do k - return get(Niii, (j, k, I), 0) == 1 + nobjj = maximum(first, keys(N[j, j, j])) + # the nested vectors contain the duals of the objects in 𝒞_ij, which are in C_ji + Niji = N[i, j, i] # 𝒞_ij x 𝒞_ji -> C_ii + Njij = N[j, i, j] # 𝒞_ji x 𝒞_ij -> C_jj + for i_ob in 1:nobji, j_ob in 1:nobjj + get(Niji, (i_ob, j_ob, Is[i]), 0) == 1 || continue # leftunit(c_ij) ∈ c_ij x c_ji + get(Njij, (j_ob, i_ob, Is[j]), 0) == 1 || continue # rightunit(c_ij) ∈ c_ji x c_ij + push!(allduals[i, j], j_ob) end end - - return I, duals end + return Is, allduals end -function Base.one(a::BimoduleSector) - a.i == a.j || error("don't know how to define one for modules") - return A4Object(a.i, a.i, _get_dual_cache(typeof(a))[a.i][1]) -end - -function Base.conj(a::BimoduleSector) - a.i == a.j || error("don't know how to define dual for modules") - return A4Object(a.i, a.i, _get_dual_cache(typeof(a))[a.i][2][a.label]) -end - -function extract_Fsymbol(::Type{A4Object}) - return mapreduce((x, y) -> cat(x, y; dims=1), 1:12) do i - filename = joinpath(artifact_path, "A4", "Fsymbol_$i.json") - @assert isfile(filename) "cannot find $filename" - json_string = read(filename, String) - Farray_part = copy(JSON3.read(json_string)) - return map(enumerate(Farray_part[Symbol(i)])) do (I, x) - j, k, l = Tuple(CartesianIndices((12, 12, 12))[I]) - y = Dict{NTuple{6,Int},Array{ComplexF64,4}}() - for (key, v) in x - a, b, c, d, e, f = parse.(Int, split(string(key)[2:(end - 1)], ", ")) - a_ob, b_ob, c_ob, d_ob, e_ob, f_ob = A4Object.(((i, j, a), (j, k, b), - (k, l, c), (i, l, d), - (i, k, e), (j, l, f))) - result = Array{ComplexF64,4}(undef, - (Nsymbol(a_ob, b_ob, e_ob), - Nsymbol(e_ob, c_ob, d_ob), - Nsymbol(b_ob, c_ob, f_ob), - Nsymbol(a_ob, f_ob, d_ob))) - map!(result, reshape(v, size(result))) do cmplxdict - return complex(cmplxdict[:re], cmplxdict[:im]) - end +function TensorKitSectors.unit(a::BimoduleSector) + a.i == a.j || throw(DomainError("unit of module category ($(a.i), $(a.j)) of $(typeof(a)) is ill-defined")) + return typeof(a)(a.i, a.i, _get_dual_cache(typeof(a))[1][a.i]) +end + +function TensorKitSectors.allunits(::Type{I}) where {I <: BimoduleSector} + s = size(I) + return I[I(i, i, _get_dual_cache(I)[1][i]) for i in 1:s] +end + +function TensorKitSectors.unit(::Type{<:BimoduleSector}) + throw(ArgumentError("unit of Type BimoduleSector doesn't exist")) +end + +function TensorKitSectors.leftunit(a::BimoduleSector) + return typeof(a)(a.i, a.i, _get_dual_cache(typeof(a))[1][a.i]) +end + +function TensorKitSectors.rightunit(a::BimoduleSector) + return typeof(a)(a.j, a.j, _get_dual_cache(typeof(a))[1][a.j]) +end - y[(a, b, c, d, e, f)] = result +function TensorKitSectors.dual(a::BimoduleSector) + return typeof(a)(a.j, a.i, _get_dual_cache(typeof(a))[2][a.i, a.j][a.label]) +end + +function extract_Fsymbol(::Type{I}) where {I <: BimoduleSector} + result = Dict{NTuple{4, Int}, Dict{NTuple{6, Int}, Array{ComplexF64, 4}}}() + name = string(BimoduleSectorName(I)) + filename = joinpath(artifact_path, name, "Fsymbol.txt") + @assert isfile(filename) "cannot find $filename" + + Farray = readdlm(filename) + for ((i, j, k, l), colordict) in convert_Fs(Farray) + result[(i, j, k, l)] = Dict{NTuple{6, Int}, Array{ComplexF64, 4}}() + for ((a, b, c, d, e, f), Fvals) in colordict + a_ob, b_ob, c_ob, d_ob, e_ob, f_ob = I.( + ((i, j, a), (j, k, b), (k, l, c), (i, l, d), (i, k, e), (j, l, f)) + ) + result[(i, j, k, l)][(a, b, c, d, e, f)] = zeros( + ComplexF64, Nsymbol(a_ob, b_ob, e_ob), Nsymbol(e_ob, c_ob, d_ob), + Nsymbol(b_ob, c_ob, f_ob), Nsymbol(a_ob, f_ob, d_ob) + ) + for (K, v) in Fvals + result[(i, j, k, l)][(a, b, c, d, e, f)][K] = v end end end + return result end -const Fcache = IdDict{Type{<:BimoduleSector}, - Array{Dict{NTuple{6,Int},Array{ComplexF64,4}},4}}() +function convert_Fs(Farray_part::Matrix{Float64}) # Farray_part is a matrix with 16 columns + data_dict = Dict{NTuple{4, Int}, Dict{NTuple{6, Int}, Vector{Pair{CartesianIndex{4}, ComplexF64}}}}() + # want to make a Dict with keys (i,j,k,l) and vals + # a Dict with keys (a,b,c,d,e,f) and vals + # a pair of (mu, nu, rho, sigma) and the F value + for row in eachrow(Farray_part) + i, j, k, l, a, b, c, d, e, f, mu, nu, rho, sigma = Int.(@view(row[1:14])) + v = complex(row[15], row[16]) + colordict = get!( + data_dict, (i, j, k, l), Dict{NTuple{6, Int}, Vector{Pair{CartesianIndex{4}, ComplexF64}}}() + ) + Fdict = get!( + colordict, (a, b, c, d, e, f), Vector{Pair{CartesianIndex{4}, ComplexF64}}() + ) + push!(Fdict, CartesianIndex(mu, nu, rho, sigma) => v) + end + return data_dict +end + +const Fcache = IdDict{Type{<:BimoduleSector}, Dict{NTuple{4, Int64}, Dict{NTuple{6, Int64}, Array{ComplexF64, 4}}}}() -function _get_Fcache(::Type{T}) where {T<:BimoduleSector} +function _get_Fcache(::Type{T}) where {T <: BimoduleSector} global Fcache return get!(Fcache, T) do return extract_Fsymbol(T) end end -function TensorKitSectors.Fsymbol(a::I, b::I, c::I, d::I, e::I, - f::I) where {I<:A4Object} - # TODO: should this error or return 0? - (a.j == b.i && b.j == c.i && a.i == d.i && c.j == d.j && - a.i == e.i && b.j == e.j && f.i == a.j && f.j == c.j) || - throw(ArgumentError("invalid fusion channel")) +function TensorKitSectors.Fsymbol(a::I, b::I, c::I, d::I, e::I, f::I) where {I <: BimoduleSector} + # required to keep track of multiplicities where F-move is partially unallowed + # also deals with invalid fusion channels + Nabe = Nsymbol(a, b, e) + Necd = Nsymbol(e, c, d) + Nbcf = Nsymbol(b, c, f) + Nafd = Nsymbol(a, f, d) + + zero_array = zeros(sectorscalartype(I), Nabe, Necd, Nbcf, Nafd) + Nabe > 0 && Necd > 0 && Nbcf > 0 && Nafd > 0 || + return zero_array i, j, k, l = a.i, a.j, b.j, c.j - return get(_get_Fcache(I)[i, j, k, l], - (a.label, b.label, c.label, d.label, e.label, f.label)) do - return zeros(sectorscalartype(A4Object), - (Nsymbol(a, b, e), Nsymbol(e, c, d), Nsymbol(b, c, f), - Nsymbol(a, f, d))) + colordict = _get_Fcache(I)[i, j, k, l] + return get!(colordict, (a.label, b.label, c.label, d.label, e.label, f.label), zero_array) +end + +# interface with TensorKit where necessary +#----------------------------------------- + +# TODO: can remove this once the otimes assert is removed +function TensorKit.fuse(V₁::GradedSpace{I}, V₂::GradedSpace{I}) where {I <: BimoduleSector} + dims = TensorKit.SectorDict{I, Int}() + for a in sectors(V₁), b in sectors(V₂) + a.j == b.i || continue # skip if not compatible + for c in a ⊗ b + dims[c] = get(dims, c, 0) + Nsymbol(a, b, c) * dim(V₁, a) * dim(V₂, b) + end end + return typeof(V₁)(dims) end + +#TODO: these might not be necessary anymore after TensorKit#291 +# check after BlockTensorKit#38 + +# function TensorKit.unitspace(S::SumSpace{<:GradedSpace{<:BimoduleSector}}) +# @assert !isempty(S) "Cannot determine type of empty space" +# return SumSpace(oneunit(first(S.spaces))) # assuming diagonal SumSpace (like in MPSKit) +# end + +# function rightunitspace(S::SumSpace{<:GradedSpace{<:BimoduleSector}}) +# @assert !isempty(S) "Cannot determine type of empty space" +# return SumSpace(rightunitspace(first(S.spaces))) +# end + +# function leftunitspace(S::SumSpace{<:GradedSpace{<:BimoduleSector}}) +# @assert !isempty(S) "Cannot determine type of empty space" +# return SumSpace(leftunitspace(first(S.spaces))) +# end diff --git a/test/runtests.jl b/test/runtests.jl index d5cbe18..6b48268 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -14,6 +14,8 @@ const GROUP = uppercase(if isnothing(arg_id) istestfile(fn) = endswith(fn, ".jl") && startswith(basename(fn), "test_") && !contains(fn, "setup") +include("setup.jl") + @time begin # tests in groups based on folder structure for testgroup in filter(isdir, readdir(@__DIR__)) diff --git a/test/setup.jl b/test/setup.jl new file mode 100644 index 0000000..8c5d2a5 --- /dev/null +++ b/test/setup.jl @@ -0,0 +1,75 @@ +module TestSetup + +export unitarity_test, rand_object, random_fusion, eval_show + +using MultiTensorKit +using TensorKitSectors +using Random + +const MTK = MultiTensorKit + +Random.seed!(1234) + +function unitarity_test(as::V, bs::V, cs::V) where {V <: AbstractVector{<:BimoduleSector}} + @assert all(a.j == b.i for a in as, b in bs) + @assert all(b.j == c.i for b in bs, c in cs) + + for a in as, b in bs, c in cs + for d in ⊗(a, b, c) + es = collect(intersect(⊗(a, b), map(dual, ⊗(c, dual(d))))) + fs = collect(intersect(⊗(b, c), map(dual, ⊗(dual(d), a)))) + Fblocks = Vector{Any}() + for e in es, f in fs + Fs = Fsymbol(a, b, c, d, e, f) + push!(Fblocks, reshape(Fs, (size(Fs, 1) * size(Fs, 2), size(Fs, 3) * size(Fs, 4)))) + end + F = hvcat(length(fs), Fblocks...) + isapprox(F' * F, one(F); atol = 1.0e-12, rtol = 1.0e-12) || return false + end + end + return true +end + +all_objects(::Type{<:BimoduleSector}, i::Int, j::Int) = [I(i, j, k) for k in 1:MTK._numlabels(I, i, j)] + +function rand_object(I::Type{<:BimoduleSector}, i::Int, j::Int) + obs = all_objects(I, i, j) + ob = rand(obs) + while isunit(ob) # unit of any fusion cat avoided + ob = rand(obs) + end + + return ob +end + +function random_fusion(I::Type{<:BimoduleSector}, i::Int, j::Int, ::Val{N}) where {N} # for fusion tree tests + N == 1 && return (rand_object(I, i, j),) + tail = random_fusion(I, i, j, Val(N - 1)) + counter = 0 + + Cs = all_objects(I, i, i) + Ds = all_objects(I, j, j) + Ms = all_objects(I, i, j) + Mops = all_objects(I, j, i) + allobs = vcat(Cs, Ds, Ms, Mops) + s = rand(allobs) + + while isempty(⊗(s, first(tail))) && counter < 40 + counter += 1 + s = (counter < 40) ? rand(allobs) : leftunit(first(tail)) + end + return (s, tail...) +end + +""" + eval_show(x) + +Use `show` to generate a string representation of `x`, then parse and evaluate the resulting expression. +""" +function eval_show(x) + str = sprint(show, x; context = (:module => @__MODULE__)) + ex = Meta.parse(str) + return eval(ex) +end + +end # end of module TestSetup diff --git a/test/test_A4.jl b/test/test_A4.jl index 696ca5b..9b2cada 100644 --- a/test/test_A4.jl +++ b/test/test_A4.jl @@ -1,54 +1,1638 @@ using MultiTensorKit -using TensorKitSectors +using TensorKitSectors, TensorKit using Test, TestExtras +using Random +using LinearAlgebra: LinearAlgebra + +const MTK = MultiTensorKit +const TK = TensorKit + +@isdefined(TestSetup) || include("setup.jl") +using .TestSetup I = A4Object +Istr = TensorKitSectors.type_repr(I) +r = size(I) -@testset "Basic type properties" begin - Istr = TensorKitSectors.type_repr(I) +println("----------------------") +println("| Sector tests |") +println("----------------------") + +@testset "$Istr Basic type properties" verbose = true begin @test eval(Meta.parse(sprint(show, I))) == I @test eval(Meta.parse(TensorKitSectors.type_repr(I))) == I end -@testset "Fusion Category $i" for i in 1:12 - objects = A4Object.(i, i, MultiTensorKit._get_dual_cache(I)[i][2]) - - @testset "Basic properties" begin - s = rand(objects, 3) - @test eval(Meta.parse(sprint(show, s[1]))) == s[1] - @test @constinferred(hash(s[1])) == hash(deepcopy(s[1])) - @test isone(@constinferred(one(s[1]))) - @constinferred dual(s[1]) - @constinferred dim(s[1]) - @constinferred frobeniusschur(s[1]) - @constinferred Bsymbol(s...) - @constinferred Fsymbol(s..., s...) +@testset "$Istr: Value iterator" begin + @test eltype(values(I)) == I + @test_throws ArgumentError unit(I) + sprev = I(1, 1, 1) # first in SectorValues + for (i, s) in enumerate(values(I)) + @test !isless(s, sprev) # confirm compatibility with sort order + @test s == @constinferred (values(I)[i]) + @test findindex(values(I), s) == i + sprev = s + i >= 10 && break + end + @test I(1, 1, 1) == first(values(I)) + @test (@constinferred findindex(values(I), I(1, 1, 1))) == 1 + for s in collect(values(I)) + @test (@constinferred values(I)[findindex(values(I), s)]) == s end +end - @testset "Unitarity of F-move" begin - for a in objects, b in objects, c in objects - for d in ⊗(a, b, c) - es = collect(intersect(⊗(a, b), map(dual, ⊗(c, dual(d))))) - fs = collect(intersect(⊗(b, c), map(dual, ⊗(dual(d), a)))) - Fblocks = Vector{Any}() - for e in es - for f in fs - Fs = Fsymbol(a, b, c, d, e, f) - push!(Fblocks, - reshape(Fs, - (size(Fs, 1) * size(Fs, 2), - size(Fs, 3) * size(Fs, 4)))) - end - end - F = hvcat(length(fs), Fblocks...) - @test isapprox(F' * F, one(F); atol=1e-12, rtol=1e-12) +@testset "$Istr ($i, $j) basic properties" for i in 1:r, j in 1:r + Cii_obs = I.(i, i, MTK._get_dual_cache(I)[2][i, i]) + Cij_obs = I.(i, j, MTK._get_dual_cache(I)[2][i, j]) + Cji_obs = I.(j, i, MTK._get_dual_cache(I)[2][j, i]) + Cjj_obs = I.(j, j, MTK._get_dual_cache(I)[2][j, j]) + c, m, mop, d = rand(Cii_obs), rand(Cij_obs), rand(Cji_obs), rand(Cjj_obs) + + if i == j + @testset "Basic fusion properties" begin + s = rand(Cii_obs, 3) + @test eval(Meta.parse(sprint(show, s[1]))) == s[1] + @test @constinferred(hash(s[1])) == hash(deepcopy(s[1])) + @test isunit(@constinferred(unit(s[1]))) + u = I.(i, i, MTK._get_dual_cache(I)[1][i]) + @test u == @constinferred(leftunit(u)) == @constinferred(rightunit(u)) == + @constinferred(unit(u)) + @test isunit(@constinferred(unit(s[1]))) + @constinferred dual(s[1]) + @test dual(s[1]) == I.(i, i, MTK._get_dual_cache(I)[2][i, i][s[1].label]) + @constinferred dim(s[1]) + @constinferred frobenius_schur_phase(s[1]) + @constinferred frobenius_schur_indicator(s[1]) + @constinferred Nsymbol(s...) + @constinferred Asymbol(s...) + @constinferred Bsymbol(s...) + F = @constinferred Fsymbol(s..., s...) + @test eltype(F) <: @testinferred sectorscalartype(I) + end + else + @testset "Basic module properties" begin + @test eval(Meta.parse(sprint(show, m))) == m + @test @constinferred(hash(m)) == hash(deepcopy(m)) + + @test isunit(m) == false + @test isunit(mop) == false + @test (isunit(@constinferred(leftunit(m))) && isunit(@constinferred(rightunit(m)))) + @test unit(c) == leftunit(m) == rightunit(mop) + @test unit(d) == rightunit(m) == leftunit(mop) + @test_throws DomainError unit(m) + @test_throws DomainError unit(mop) + + @constinferred dual(m) + @test dual(m) == I.(j, i, MTK._get_dual_cache(I)[2][i, j][m.label]) + @test dual(dual(m)) == m + + @constinferred dim(m) + @constinferred frobenius_schur_phase(m) + @constinferred frobenius_schur_indicator(m) + @constinferred Bsymbol(m, mop, c) + @constinferred Fsymbol(mop, m, mop, mop, d, c) + end + + @testset "$Istr Fusion rules" begin + argerr = ArgumentError("invalid fusion channel") + # forbidden fusions + for obs in [(c, d), (d, c), (m, m), (mop, mop), (d, m), (m, c), (mop, d), (c, mop)] + @test_throws AssertionError("a.j == b.i") isempty(⊗(obs...)) + @test_throws argerr Nsymbol(obs..., rand([c, m, mop, d])) end + + # allowed fusions + for obs in [(c, c), (d, d), (m, mop), (mop, m), (c, m), (mop, c), (m, d), (d, mop)] + @test !isempty(⊗(obs...)) + end + + @test Nsymbol(c, unit(c), c) == Nsymbol(d, unit(d), d) == 1 + + @test_throws argerr Nsymbol(m, mop, d) + @test_throws argerr Nsymbol(mop, m, c) + @test_throws argerr Fsymbol(m, mop, m, mop, c, d) end end +end + +println("-----------------------------") +println("| F-symbol data tests |") +println("-----------------------------") + +for i in 1:r, j in 1:r + @testset "Unitarity of $Istr F-move ($i, $j)" begin + if i == j + @testset "Unitarity of fusion F-move ($i, $j)" begin + fusion_objects = I.(i, i, MTK._get_dual_cache(I)[2][i, i]) + @test unitarity_test(fusion_objects, fusion_objects, fusion_objects) + end + end + + i != j || continue # do this part only when off-diagonal + mod_objects = I.(i, j, MTK._get_dual_cache(I)[2][i, j]) + left_fusion_objects = I.(i, i, MTK._get_dual_cache(I)[2][i, i]) + right_fusion_objects = I.(j, j, MTK._get_dual_cache(I)[2][j, j]) + + # C x C x M -> M or D x D x Mop -> Mop + @testset "Unitarity of left module F-move ($i, $j)" begin + @test unitarity_test(left_fusion_objects, left_fusion_objects, mod_objects) + end + + # M x D x D -> M or Mop x C x C -> Mop + @testset "Unitarity of right module F-move ($i, $j)" begin + @test unitarity_test(mod_objects, right_fusion_objects, right_fusion_objects) + end - @testset "Pentagon equation" begin - for a in objects, b in objects, c in objects, d in objects - @test pentagon_equation(a, b, c, d; atol=1e-12, rtol=1e-12) + # C x M x D -> M or D x Mop x C -> Mop + @testset "Unitarity of bimodule F-move ($i, $j)" begin + @test unitarity_test(left_fusion_objects, mod_objects, right_fusion_objects) + end + + @testset "Unitarity of mixed module F-move ($i, $j) and opposite ($j, $i)" begin + modop_objects = I.(j, i, MTK._get_dual_cache(I)[2][j, i]) + + # C x M x Mop -> C or D x Mop x M -> D + @test unitarity_test(left_fusion_objects, mod_objects, modop_objects) + # M x Mop x C -> C or Mop x M x D -> D + @test unitarity_test(mod_objects, modop_objects, left_fusion_objects) + # Mop x C x M -> D or M x D x Mop -> C + @test unitarity_test(modop_objects, left_fusion_objects, mod_objects) + + # M x Mop x M -> M or Mop x M x Mop -> Mop + @test unitarity_test(mod_objects, modop_objects, mod_objects) + end + end +end + +@testset "Triangle equation" begin + objects = collect(values(I)) + for a in objects, b in objects + a.j == b.i || continue # skip if not compatible + @test triangle_equation(a, b; atol = 1.0e-12, rtol = 1.0e-12) + end +end + +@testset "$Istr Pentagon equation" begin + objects = collect(values(I)) + for a in objects + for b in objects + a.j == b.i || continue # skip if not compatible + for c in objects + b.j == c.i || continue # skip if not compatible + for d in objects + c.j == d.i || continue # skip if not compatible + @test pentagon_equation(a, b, c, d; atol = 1.0e-12, rtol = 1.0e-12) + end + end end end end + +### start of TensorKit tests ### + +# println("---------------------------------") +# println("| Multifusion space tests |") +# println("---------------------------------") + +# @timedtestset "Multifusion spaces " verbose = true begin +# @timedtestset "GradedSpace: $(TK.type_repr(Vect[I]))" begin +# gen = (values(I)[k] => (k + 1) for k in 1:length(values(I))) + +# V = GradedSpace(gen) +# @test eval(Meta.parse(type_repr(typeof(V)))) == typeof(V) +# @test eval_show(V) == V +# @test eval_show(V') == V' +# @test V' == GradedSpace(gen; dual = true) +# @test V == @constinferred GradedSpace(gen...) +# @test V' == @constinferred GradedSpace(gen...; dual = true) +# @test V == @constinferred GradedSpace(tuple(gen...)) +# @test V' == @constinferred GradedSpace(tuple(gen...); dual = true) +# @test V == @constinferred GradedSpace(Dict(gen)) +# @test V' == @constinferred GradedSpace(Dict(gen); dual = true) +# @test V == @inferred Vect[I](gen) +# @test V' == @constinferred Vect[I](gen; dual = true) +# @test V == @constinferred Vect[I](gen...) +# @test V' == @constinferred Vect[I](gen...; dual = true) +# @test V == @constinferred Vect[I](Dict(gen)) +# @test V' == @constinferred Vect[I](Dict(gen); dual = true) +# @test @constinferred(hash(V)) == hash(deepcopy(V)) != hash(V') +# @test V == GradedSpace(reverse(collect(gen))...) +# @test eval_show(V) == V +# @test eval_show(typeof(V)) == typeof(V) + +# @test dim(@constinferred(zerospace(V))) == 0 + +# W = @constinferred GradedSpace(unit => 1 for unit in allunits(I)) +# dict = Dict(unit => 1 for unit in allunits(I)) +# @test W == GradedSpace(dict) +# @test W == GradedSpace(push!(dict, randsector(I) => 0)) +# @test @constinferred(zerospace(V)) == GradedSpace(unit => 0 for unit in allunits(I)) +# randunit = rand(collect(allunits(I))) +# @test_throws ArgumentError("Sector $(randunit) appears multiple times") GradedSpace(randunit => 1, randunit => 3) + +# @test isunitspace(W) +# @test @constinferred(unitspace(V)) == W == unitspace(typeof(V)) +# @test_throws ArgumentError leftunitspace(V) +# @test_throws ArgumentError rightunitspace(V) +# @test eval_show(W) == W + +# @test isa(V, VectorSpace) +# @test isa(V, ElementarySpace) +# @test isa(InnerProductStyle(V), HasInnerProduct) +# @test isa(InnerProductStyle(V), EuclideanInnerProduct) +# @test isa(V, GradedSpace) +# @test isa(V, GradedSpace{I}) +# @test @constinferred(dual(V)) == @constinferred(conj(V)) == @constinferred(adjoint(V)) != V +# @test @constinferred(field(V)) == ℂ +# @test @constinferred(sectortype(V)) == I +# slist = @constinferred sectors(V) +# @test @constinferred(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)) + +# @test @constinferred(⊕(V, zerospace(V))) == V +# @test @constinferred(⊕(V, V)) == Vect[I](c => 2dim(V, c) for c in sectors(V)) +# @test @constinferred(⊕(V, V, V, V)) == Vect[I](c => 4dim(V, c) for c in sectors(V)) +# @test @constinferred(⊕(V, unitspace(V))) == Vect[I](c => isunit(c) + dim(V, c) for c in sectors(V)) +# @test @constinferred(fuse(V, unitspace(V))) == V + +# @testset "$Istr ($i, $j) spaces" for i in 1:r, j in 1:r #TODO: look at these tests better +# # space with a single sector +# Wleft = @constinferred Vect[I]((i, i, label) => 1 for label in 1:MTK._numlabels(I, i, i)) +# Wright = @constinferred Vect[I]((j, j, label) => 1 for label in 1:MTK._numlabels(I, j, j)) +# WM = @constinferred Vect[I]((i, j, label) => 1 for label in 1:MTK._numlabels(I, i, j)) +# WMop = @constinferred Vect[I]((j, i, label) => 1 for label in 1:MTK._numlabels(I, j, i)) + +# @test leftunitspace(Wleft) == rightunitspace(Wleft) +# @test leftunitspace(Wright) == rightunitspace(Wright) +# @test @constinferred(leftunitspace(⊕(Wleft, WM))) == leftunitspace(Wleft) +# @test @constinferred(leftunitspace(⊕(Wright, WMop))) == leftunitspace(Wright) +# @test @constinferred(rightunitspace(⊕(Wright, WM))) == rightunitspace(Wright) +# @test @constinferred(rightunitspace(⊕(Wleft, WMop))) == rightunitspace(Wleft) + +# if i != j # some tests specialised for modules + +# # sensible direct sums and fuses +# ul, ur = unit(I(i, i, 1)), unit(I(j, j, 1)) +# @test @constinferred(⊕(Wleft, WM)) == +# Vect[I](c => 1 for c in sectors(V) if leftunit(c) == ul == rightunit(c) || (c.i == i && c.j == j)) +# @test @constinferred(⊕(Wright, WMop)) == +# Vect[I](c => 1 for c in sectors(V) if leftunit(c) == ur == rightunit(c) || (c.i == j && c.j == i)) +# @test @constinferred(⊕(Wright, WM)) == +# Vect[I](c => 1 for c in sectors(V) if rightunit(c) == ur == leftunit(c) || (c.i == i && c.j == j)) +# @test @constinferred(⊕(Wleft, WMop)) == +# Vect[I](c => 1 for c in sectors(V) if rightunit(c) == ul == leftunit(c) || (c.i == j && c.j == i)) +# # round needed below because of numerical F-symbols not being integer when they should be +# # although this test might be stupid, because I'm assuming integer qdims bc everything's a group or irrep on the diagonal +# @test @constinferred(fuse(Wleft, WM)) == Vect[I](c => round(Int, dim(Wleft)) for c in sectors(WM)) # this might be wrong +# @test @constinferred(fuse(Wright, WMop)) == Vect[I](c => round(Int, dim(Wright)) for c in sectors(WMop)) # same + +# # less sensible fuse +# @test @constinferred(fuse(Wleft, WMop)) == fuse(Wright, WM) == +# Vect[I](c => 0 for c in sectors(V)) + +# for W in [WM, WMop, Wright] +# @test infimum(Wleft, W) == Vect[I](c => 0 for c in sectors(V)) +# end +# else +# @test @constinferred(⊕(Wleft, Wright)) == +# Vect[I](c => 2 for c in sectors(V) if c.i == c.j == i) +# @test @constinferred(fuse(Wleft, WMop)) == fuse(Wright, WM) +# end + +# for W in [Wleft, Wright] +# @test @constinferred(⊕(W, rightunitspace(W))) == +# Vect[I](c => isunit(c) + dim(W, c) for c in sectors(W)) +# @test @constinferred(fuse(W, rightunitspace(W))) == W +# end +# end + +# d = Dict{I, Int}() +# for a in sectors(V), b in sectors(V) +# a.j == b.i || continue # skip if not compatible +# for c in a ⊗ b +# d[c] = get(d, c, 0) + dim(V, a) * dim(V, b) * Nsymbol(a, b, c) +# end +# end +# @test @constinferred(fuse(V, V)) == GradedSpace(d) +# @test @constinferred(flip(V)) == +# Vect[I](conj(c) => dim(V, c) for c in sectors(V))' +# @test flip(V) ≅ V +# @test flip(V) ≾ V +# @test flip(V) ≿ V +# @test @constinferred(⊕(V, V)) == @constinferred supremum(V, ⊕(V, V)) +# @test V == @constinferred infimum(V, ⊕(V, V)) +# @test V ≺ ⊕(V, V) +# @test !(V ≻ ⊕(V, V)) + +# randlen = rand(1:length(values(I))) +# s = rand(collect(values(I))[randlen:end]) # such that dim(V, s) > randlen +# @test infimum(V, GradedSpace(s => randlen)) == GradedSpace(s => randlen) +# @test_throws SpaceMismatch (⊕(V, V')) +# end + +# @timedtestset "HomSpace with $(TK.type_repr(Vect[I])) involving ($i, $j)" for i in 1:r, j in 1:r +# V1, V2, V3, V4, V5 = ( +# Vect[I]((i, i, label) => 1 for label in 1:MTK._numlabels(I, i, i)), +# Vect[I]((i, j, label) => 1 for label in 1:MTK._numlabels(I, i, j)), +# Vect[I]((i, i, label) => 1 for label in 1:MTK._numlabels(I, i, i)), # same as V1 +# Vect[I]((i, j, 1) => 3), +# Vect[I]((j, j, label) => 1 for label in 1:MTK._numlabels(I, j, j)), +# ) +# W = HomSpace(V1 ⊗ V2, V3 ⊗ V4 ⊗ V5) + +# @test W == (V3 ⊗ V4 ⊗ V5 → V1 ⊗ V2) +# @test W == (V1 ⊗ V2 ← V3 ⊗ V4 ⊗ V5) +# @test W' == (V1 ⊗ V2 → V3 ⊗ V4 ⊗ V5) +# @test eval(Meta.parse(sprint(show, W))) == W +# @test eval(Meta.parse(sprint(show, typeof(W)))) == typeof(W) +# @test spacetype(W) == typeof(V1) +# @test sectortype(W) == sectortype(V1) +# @test W[1] == V1 +# @test W[2] == V2 +# @test W[3] == V3' +# @test W[4] == V4' +# @test W[5] == V5' + +# @test @constinferred(hash(W)) == hash(deepcopy(W)) != hash(W') +# @test W == deepcopy(W) +# @test W == @constinferred permute(W, ((1, 2), (3, 4, 5))) +# @test permute(W, ((2, 4, 5), (3, 1))) == (V2 ⊗ V4' ⊗ V5' ← V3 ⊗ V1') +# @test (V1 ⊗ V2 ← V1 ⊗ V2) == @constinferred TK.compose(W, W') + +# @test (V1 ⊗ V2 ← V3 ⊗ V4 ⊗ V5 ⊗ rightunitspace(V5)) == +# @constinferred(insertleftunit(W)) == +# @constinferred(insertrightunit(W)) +# @test @constinferred(removeunit(insertleftunit(W), $(numind(W) + 1))) == W +# @test_throws BoundsError insertrightunit(W, 6) +# @test_throws BoundsError insertleftunit(W, 0) + +# @test (V1 ⊗ V2 ⊗ rightunitspace(V2) ← V3 ⊗ V4 ⊗ V5) == +# @constinferred(insertrightunit(W, 2)) +# @test (V1 ⊗ V2 ← leftunitspace(V3) ⊗ V3 ⊗ V4 ⊗ V5) == +# @constinferred(insertleftunit(W, 3)) +# @test @constinferred(removeunit(insertleftunit(W, 3), 3)) == W +# @test_throws ArgumentError @constinferred(insertrightunit(one(V1) ← V1, 0)) # should I specify it's the other error? +# @test_throws ArgumentError insertleftunit(one(V1) ← V1, 0) +# end +# end + +# println("---------------------------------------") +# println("| Multifusion fusion tree tests |") +# println("---------------------------------------") + +# @timedtestset "Fusion trees for $(TK.type_repr(I)) involving ($i, $j)" verbose = true for i in 1:r, j in 1:r +# N = 6 +# out = random_fusion(I, i, j, Val(N)) +# isdual = ntuple(n -> rand(Bool), N) +# in = rand(collect(⊗(out...))) # will be in 𝒞ⱼⱼ with this choice of out + +# numtrees = length(fusiontrees(out, in, isdual)) # will be 1 for i != j +# @test numtrees == count(n -> true, fusiontrees(out, in, isdual)) + +# it = @constinferred fusiontrees(out, in, isdual) +# @constinferred Nothing iterate(it) +# f, s = iterate(it) +# @constinferred Nothing iterate(it, s) +# @test f == @constinferred first(it) +# @testset "Fusion tree $Istr: printing" begin +# @test eval(Meta.parse(sprint(show, f))) == f +# end + +# C0, D0 = unit(I(i, i, 1)), unit(I(j, j, 1)) +# @testset "Fusion tree $Istr: constructor properties" for u in (C0, D0) +# @constinferred FusionTree((), u, (), (), ()) +# @constinferred FusionTree((u,), u, (false,), (), ()) +# @constinferred FusionTree((u, u), u, (false, false), (), (1,)) +# @constinferred FusionTree((u, u, u), u, (false, false, false), (u,), (1, 1)) +# @constinferred FusionTree( +# (u, u, u, u), u, (false, false, false, false), (u, u), (1, 1, 1) +# ) +# @test_throws MethodError FusionTree((u, u, u), u, (false, false), (u,), (1, 1)) +# @test_throws MethodError FusionTree( +# (u, u, u), u, (false, false, false), (u, u), (1, 1) +# ) +# @test_throws MethodError FusionTree( +# (u, u, u), u, (false, false, false), (u,), (1, 1, 1) +# ) +# @test_throws MethodError FusionTree((u, u, u), u, (false, false, false), (), (1,)) + +# f = FusionTree((u, u, u), u, (false, false, false), (u,), (1, 1)) +# @test sectortype(f) == I +# @test length(f) == 3 +# @test FusionStyle(f) == FusionStyle(I) +# @test BraidingStyle(f) == BraidingStyle(I) + +# if FusionStyle(I) isa UniqueFusion +# @constinferred FusionTree((), u, ()) +# @constinferred FusionTree((u,), u, (false,)) +# @constinferred FusionTree((u, u), u, (false, false)) +# @constinferred FusionTree((u, u, u), u) +# if UnitStyle(I) isa SimpleUnit +# @constinferred FusionTree((u, u, u, u)) +# else +# @test_throws ArgumentError FusionTree((u, u, u, u)) +# end +# @test_throws MethodError FusionTree((u, u), u, (false, false, false)) +# else +# @test_throws ArgumentError FusionTree((), u, ()) +# @test_throws ArgumentError FusionTree((u,), u, (false,)) +# @test_throws ArgumentError FusionTree((u, u), u, (false, false)) +# @test_throws ArgumentError FusionTree((u, u, u), u) +# if I <: ProductSector && UnitStyle(I) isa GenericUnit +# @test_throws DomainError FusionTree((u, u, u, u)) +# else +# @test_throws ArgumentError FusionTree((u, u, u, u)) +# end +# end +# end + +# @testset "Fusion tree $Istr: insertat" begin +# N = 4 +# out2 = random_fusion(I, i, j, Val(N)) +# in2 = rand(collect(⊗(out2...))) +# isdual2 = ntuple(n -> rand(Bool), N) +# f2 = rand(collect(fusiontrees(out2, in2, isdual2))) +# for k in 1:N +# out1 = random_fusion(I, i, j, Val(N)) # guaranteed good fusion +# out1 = Base.setindex(out1, in2, i) # can lead to poor fusion +# while isempty(⊗(out1...)) # TODO: better way to do this? +# out1 = random_fusion(I, i, j, Val(N)) +# out1 = Base.setindex(out1, in2, i) +# end +# in1 = rand(collect(⊗(out1...))) +# isdual1 = ntuple(n -> rand(Bool), N) +# isdual1 = Base.setindex(isdual1, false, k) +# f1 = rand(collect(fusiontrees(out1, in1, isdual1))) + +# trees = @constinferred TK.insertat(f1, k, f2) +# @test norm(values(trees)) ≈ 1 + +# f1a, f1b = @constinferred TK.split(f1, $k) +# @test length(TK.insertat(f1b, 1, f1a)) == 1 +# @test first(TK.insertat(f1b, 1, f1a)) == (f1 => 1) + +# # no braid tests for non-hardcoded example +# end +# end +# # no planar trace tests + +# @testset "Fusion tree $Istr: elementary artin braid" begin +# N = length(out) +# isdual = ntuple(n -> rand(Bool), N) +# # no general artin braid test + +# # not sure how useful this test is, it does the trivial braiding (choice of out) +# f = rand(collect(it)) # in this case the 1 tree +# d1 = TK.artin_braid(f, 2) # takes unit C0 with current out +# d2 = empty(d1) +# for (f1, coeff1) in d1 +# for (f2, coeff2) in TK.artin_braid(f1, 3) +# d2[f2] = get(d2, f2, zero(coeff1)) + coeff2 * coeff1 +# end +# end +# d1 = d2 +# d2 = empty(d1) +# for (f1, coeff1) in d1 +# for (f2, coeff2) in TK.artin_braid(f1, 3; inv = true) +# d2[f2] = get(d2, f2, zero(coeff1)) + coeff2 * coeff1 +# end +# end +# d1 = d2 +# d2 = empty(d1) +# for (f1, coeff1) in d1 +# for (f2, coeff2) in TK.artin_braid(f1, 2; inv = true) +# d2[f2] = get(d2, f2, zero(coeff1)) + coeff2 * coeff1 +# end +# end +# d1 = d2 +# for (f1, coeff1) in d1 +# if f1 == f +# @test coeff1 ≈ 1 +# else +# @test isapprox(coeff1, 0; atol = 1.0e-12, rtol = 1.0e-12) +# end +# end +# end + +# # no braiding and permuting test +# @testset "Fusion tree $Istr: merging" begin +# N = 3 +# out1 = random_fusion(I, i, j, N) +# out2 = random_fusion(I, i, j, N) +# in1 = rand(collect(⊗(out1...))) +# in2 = rand(collect(⊗(out2...))) +# tp = ⊗(in1, in2) # messy solution but it works +# while isempty(tp) +# out1 = random_fusion(I, i, j, Val(N)) +# out2 = random_fusion(I, i, j, Val(N)) +# in1 = rand(collect(⊗(out1...))) +# in2 = rand(collect(⊗(out2...))) +# tp = ⊗(in1, in2) +# end + +# f1 = rand(collect(fusiontrees(out1, in1))) +# f2 = rand(collect(fusiontrees(out2, in2))) + + +# @test dim(in1) * dim(in2) ≈ sum( +# abs2(coeff) * dim(c) for c in in1 ⊗ in2 +# for μ in 1:Nsymbol(in1, in2, c) +# for (f, coeff) in TK.merge(f1, f2, c, μ) +# ) +# # no merge and braid interplay tests +# end + +# # double fusion tree tests +# N = 4 +# out = random_fusion(I, i, j, Val(N)) +# out2 = random_fusion(I, i, j, Val(N)) +# tp = ⊗(out...) +# tp2 = ⊗(out2...) +# while isempty(intersect(tp, tp2)) # guarantee fusion to same coloring +# out2 = random_fusion(I, i, j, Val(N)) +# tp2 = ⊗(out2...) +# end +# @test_throws ArgumentError fusiontrees((out..., map(dual, out)...)) +# incoming = rand(collect(intersect(tp, tp2))) +# f1 = rand(collect(fusiontrees(out, incoming, ntuple(n -> rand(Bool), N)))) +# f2 = rand(collect(fusiontrees(out2, incoming, ntuple(n -> rand(Bool), N)))) # no permuting + +# @testset "Double fusion tree $Istr: repartitioning" begin +# for n in 0:(2 * N) +# d = @constinferred TK.repartition(f1, f2, $n) +# @test dim(incoming) ≈ +# sum(abs2(coef) * dim(f1.coupled) for ((f1, f2), coef) in d) +# d2 = Dict{typeof((f1, f2)), valtype(d)}() +# for ((f1′, f2′), coeff) in d +# for ((f1′′, f2′′), coeff2) in TK.repartition(f1′, f2′, N) +# d2[(f1′′, f2′′)] = get(d2, (f1′′, f2′′), zero(coeff)) + coeff2 * coeff +# end +# end +# for ((f1′, f2′), coeff2) in d2 +# if f1 == f1′ && f2 == f2′ +# @test coeff2 ≈ 1 +# else +# @test isapprox(coeff2, 0; atol = 1.0e-12, rtol = 1.0e-12) +# end +# end +# end +# end + +# # no double fusion tree permutation tests + +# # very slow for (1, 6), (3, 4), (3, 5), (3, 6), (5, 6), (6, 1), (6, 5), (7, 1), (7, 4), (7, 6) +# @testset "Double fusion tree $Istr: transposition" begin +# for n in 0:(2N) +# i0 = rand(1:(2N)) +# p = mod1.(i0 .+ (1:(2N)), 2N) +# ip = mod1.(-i0 .+ (1:(2N)), 2N) +# p′ = tuple(getindex.(Ref(vcat(1:N, (2N):-1:(N + 1))), p)...) +# p1, p2 = p′[1:n], p′[(2N):-1:(n + 1)] +# ip′ = tuple(getindex.(Ref(vcat(1:n, (2N):-1:(n + 1))), ip)...) +# ip1, ip2 = ip′[1:N], ip′[(2N):-1:(N + 1)] + +# d = @constinferred transpose(f1, f2, p1, p2) +# @test dim(incoming) ≈ +# sum(abs2(coef) * dim(f1.coupled) for ((f1, f2), coef) in d) +# d2 = Dict{typeof((f1, f2)), valtype(d)}() +# for ((f1′, f2′), coeff) in d +# d′ = transpose(f1′, f2′, ip1, ip2) +# for ((f1′′, f2′′), coeff2) in d′ +# d2[(f1′′, f2′′)] = get(d2, (f1′′, f2′′), zero(coeff)) + coeff2 * coeff +# end +# end +# for ((f1′, f2′), coeff2) in d2 +# if f1 == f1′ && f2 == f2′ +# @test coeff2 ≈ 1 +# else +# @test abs(coeff2) < 1.0e-12 +# end +# end +# end +# end + +# @testset "Double fusion tree $Istr: planar trace" begin +# d1 = transpose(f1, f1, (N + 1, 1:N..., ((2N):-1:(N + 3))...), (N + 2,)) +# f1front, = TK.split(f1, N - 1) +# T = sectorscalartype(I) +# d2 = Dict{typeof((f1front, f1front)), T}() +# for ((f1′, f2′), coeff′) in d1 +# for ((f1′′, f2′′), coeff′′) in +# TK.planar_trace( +# f1′, f2′, (2:N...,), (1, ((2N):-1:(N + 3))...), (N + 1,), +# (N + 2,) +# ) +# coeff = coeff′ * coeff′′ +# d2[(f1′′, f2′′)] = get(d2, (f1′′, f2′′), zero(coeff)) + coeff +# end +# end +# for ((f1_, f2_), coeff) in d2 +# if (f1_, f2_) == (f1front, f1front) +# @test coeff ≈ dim(f1.coupled) / dim(f1front.coupled) +# else +# @test abs(coeff) < 1.0e-12 +# end +# end +# end +# end + +# println("-------------------------------------------") +# println("| Multifusion diagonal tensor tests |") +# println("-------------------------------------------") + +# V = Vect[I](values(I)[k] => 1 for k in 1:length(values(I))) + +# @timedtestset "DiagonalTensor" begin +# @timedtestset "Basic properties and algebra" begin +# for T in (Float32, Float64, ComplexF32, ComplexF64, BigFloat) +# # constructors +# t = @constinferred DiagonalTensorMap{T}(undef, V) +# t = @constinferred DiagonalTensorMap(rand(T, reduceddim(V)), V) +# t2 = @constinferred DiagonalTensorMap{T}(undef, space(t)) +# @test space(t2) == space(t) +# @test_throws ArgumentError DiagonalTensorMap{T}(undef, V^2 ← V) +# t2 = @constinferred DiagonalTensorMap{T}(undef, domain(t)) +# @test space(t2) == space(t) +# @test_throws ArgumentError DiagonalTensorMap{T}(undef, V^2) +# # properties +# @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)) +# # blocks +# bs = @constinferred blocks(t) +# (c, b1), state = @constinferred Nothing iterate(bs) +# @test c == first(blocksectors(V ← V)) +# next = @constinferred Nothing iterate(bs, state) +# b2 = @constinferred block(t, first(blocksectors(t))) +# @test b1 == b2 +# @test eltype(bs) === Pair{typeof(c), typeof(b1)} +# @test typeof(b1) === TensorKit.blocktype(t) +# # basic linear algebra +# @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) +# if T != BigFloat # seems broken for now +# @test norm(one!(t) - id(V)) == 0 +# end + +# 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 + +# @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) +# # element type conversion +# tc = complex(t) +# @test convert(typeof(tc), t) == tc +# @test typeof(convert(typeof(tc), t)) == typeof(tc) +# # to and from generic TensorMap +# td = DiagonalTensorMap(TensorMap(t)) +# @test t == td +# @test typeof(td) == typeof(t) +# 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)) +# ) +# # no V * V' * V ← V or V^2 ← V tests due to Nsymbol erroring where fusion is forbidden +# end +# @timedtestset "Tensor contraction " for i in 1:r +# W = Vect[I]((i, i, label) => 2 for label in 1:MTK._numlabels(I, i, i)) + +# d = DiagonalTensorMap(rand(ComplexF64, reduceddim(W)), W) +# t = TensorMap(d) +# A = randn(ComplexF64, W ⊗ W' ⊗ W, W) +# B = randn(ComplexF64, W ⊗ W' ⊗ W, W ⊗ W') # empty for modules so untested + +# @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 "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 + +# # no conversion tests because no fusion tensor +# # no permute tests: NoBraiding() + +# println("---------------------------------------") +# println("Tensors with symmetry: $Istr") +# println("---------------------------------------") + +# @timedtestset "Tensors with symmetry involving $Istr ($i, $j)" verbose = true for i in 1:r, j in 1:r +# isdiag = i == j + +# VC = ( +# Vect[I]((i, i, label) => 1 for label in 1:MTK._numlabels(I, i, i)), +# Vect[I](unit(I(i, i, 1)) => 1, rand_object(I, i, i) => 1), # avoids OOMs? +# Vect[I](unit(I(i, i, 1)) => 2, rand_object(I, i, i) => 1), +# Vect[I]((i, i, label) => 1 for label in 1:MTK._numlabels(I, i, i)), +# Vect[I](unit(I(i, i, 1)) => 2, rand_object(I, i, i) => 3), +# ) + +# VM = Vect[I]((i, j, label) => 1 for label in 1:MTK._numlabels(I, i, j)) # all module objects + +# VM1 = ( +# Vect[I](unit(I(i, i, 1)) => 1, rand_object(I, i, i) => 1), # written so V1 ⊗ V2 ← V3 ⊗ V4 ⊗ V5 works +# Vect[I](rand_object(I, i, j) => 2), # generally less blocksectors +# Vect[I](unit(I(i, i, 1)) => 1, rand_object(I, i, i) => 1), +# VM, # important that V4 is module-graded +# Vect[I](unit(I(j, j, 1)) => 2, rand_object(I, j, j) => 1), +# ) + +# VM2 = ( +# Vect[I](rand_object(I, i, j) => 2), # second set where module is V1 here +# Vect[I](unit(I(j, j, 1)) => 1, rand_object(I, j, j) => 1), +# Vect[I](unit(I(i, i, 1)) => 1, rand_object(I, i, i) => 1), +# VM, +# Vect[I](unit(I(j, j, 1)) => 2, rand_object(I, j, j) => 1), +# ) + +# Vcol = isdiag ? (VC,) : (VM1, VM2) # avoid duplicate runs + +# for V in Vcol # TODO: add enumerate to keep track of potential erroring space +# V1, V2, V3, V4, V5 = V +# @timedtestset "Basic tensor properties" begin +# W = isdiag ? V1 ⊗ V2 ⊗ V3 ⊗ V4 ⊗ V5 : V3 ⊗ V4 ⊗ V5 # fusion matters +# for T in (Int, Float32, Float64, ComplexF32, ComplexF64, BigFloat) +# t = @constinferred zeros(T, W) # empty for i != j b/c blocks are module-graded +# @test @constinferred(hash(t)) == hash(deepcopy(t)) +# @test scalartype(t) == T +# @test norm(t) == 0 +# @test codomain(t) == W +# @test space(t) == (W ← one(W)) +# @test domain(t) == one(W) +# @test typeof(t) == TensorMap{T, spacetype(t), length(W), 0, Vector{T}} +# # blocks +# bs = @constinferred blocks(t) +# if !isempty(bs) +# (c, b1), state = @constinferred Nothing iterate(bs) # errors if fusion gives empty data +# # @test c == first(blocksectors(W)) # unit doesn't have label 1 +# next = @constinferred Nothing iterate(bs, state) +# b2 = @constinferred block(t, first(blocksectors(t))) +# @test b1 == b2 +# @test eltype(bs) === Pair{typeof(c), typeof(b1)} +# @test typeof(b1) === TK.blocktype(t) +# @test typeof(c) === sectortype(t) +# end +# end +# end + +# @timedtestset "Tensor Dict conversion" begin +# W = V1 ⊗ V2 ← V3 ⊗ V4 ⊗ V5 # rewritten to be compatible with module fusion +# for T in (Int, Float32, ComplexF64) +# t = @constinferred rand(T, W) +# d = convert(Dict, t) +# @test t == convert(TensorMap, d) +# end +# end + +# @timedtestset "Basic linear algebra" begin +# W = V1 ⊗ V2 ← V3 ⊗ V4 ⊗ V5 +# for T in (Float32, ComplexF64) +# t = @constinferred rand(T, W) # fusion matters here +# @test scalartype(t) == T +# @test space(t) == W +# @test space(t') == W' +# @test dim(t) == dim(space(t)) +# @test codomain(t) == codomain(W) +# @test domain(t) == domain(W) +# # blocks for adjoint +# bs = @constinferred blocks(t') +# (c, b1), state = @constinferred Nothing iterate(bs) +# @test c == first(blocksectors(W')) +# next = @constinferred Nothing iterate(bs, state) +# b2 = @constinferred block(t', first(blocksectors(t'))) +# @test b1 == b2 +# @test eltype(bs) === Pair{typeof(c), typeof(b1)} +# @test typeof(b1) === TensorKit.blocktype(t') +# @test typeof(c) === sectortype(t) +# # linear algebra +# @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') + +# t2 = @constinferred rand!(similar(t)) +# β = rand(T) +# @test @constinferred(dot(β * t2, α * t)) ≈ conj(β) * α * conj(dot(t, t2)) +# @test dot(t2, t) ≈ conj(dot(t, t2)) +# @test dot(t2, t) ≈ conj(dot(t2', t')) +# @test dot(t2, t) ≈ dot(t', t2') + +# if !isempty(blocksectors(V2 ⊗ V1)) +# i1 = @constinferred(isomorphism(T, V1 ⊗ V2, V2 ⊗ V1)) # can't reverse fusion here when modules are involved +# i2 = @constinferred(isomorphism(Vector{T}, V2 ⊗ V1, V1 ⊗ V2)) +# @test i1 * i2 == @constinferred(id(T, V1 ⊗ V2)) +# @test i2 * i1 == @constinferred(id(Vector{T}, V2 ⊗ V1)) +# end + +# w = @constinferred isometry(T, V1 ⊗ (rightunitspace(V1) ⊕ rightunitspace(V1)), V1) +# @test dim(w) == 2 * dim(V1 ← V1) +# @test w' * w == id(Vector{T}, V1) +# @test w * w' == (w * w')^2 +# end +# end + +# @timedtestset "Trivial space insertion and removal" begin +# W = V1 ⊗ V2 ← V3 ⊗ V4 ⊗ V5 +# for T in (Float32, ComplexF64) +# t = @constinferred rand(T, W) # fusion matters here +# t2 = @constinferred insertleftunit(t) +# @test t2 == @constinferred insertrightunit(t) +# @test space(t2) == insertleftunit(space(t)) +# @test @constinferred(removeunit(t2, $(numind(t2)))) == t +# t3 = @constinferred insertleftunit(t; copy = true) +# @test t3 == @constinferred insertrightunit(t; copy = true) +# @test @constinferred(removeunit(t3, $(numind(t3)))) == t + +# @test numind(t2) == numind(t) + 1 +# @test scalartype(t2) === T +# @test t.data === t2.data + +# @test t.data !== t3.data +# for (c, b) in blocks(t) +# @test b == block(t3, c) +# end + +# t4 = @constinferred insertrightunit(t, 3; dual = true) +# @test numin(t4) == numin(t) + 1 && numout(t4) == numout(t) +# for (c, b) in blocks(t) +# @test b == block(t4, c) +# end +# @test @constinferred(removeunit(t4, 4)) == t + +# t5 = @constinferred insertleftunit(t, 4; dual = true) +# @test numin(t5) == numin(t) + 1 && numout(t5) == numout(t) +# for (c, b) in blocks(t) +# @test b == block(t5, c) +# end +# @test @constinferred(removeunit(t5, 4)) == t +# end +# end + +# @timedtestset "Tensor conversion" begin +# W = V1 ⊗ V2 +# t = @constinferred randn(W ← W) # fusion matters here +# @test typeof(convert(TensorMap, t')) == typeof(t) +# tc = complex(t) +# @test convert(typeof(tc), t) == tc +# @test typeof(convert(typeof(tc), t)) == typeof(tc) +# @test typeof(convert(typeof(tc), t')) == typeof(tc) +# @test Base.promote_typeof(t, tc) == typeof(tc) +# @test Base.promote_typeof(tc, t) == typeof(tc + t) +# end + +# @timedtestset "Full trace: test self-consistency" begin +# t = rand(ComplexF64, V1 ⊗ V2 ← V1 ⊗ V2) # avoid permutes +# ss = @constinferred tr(t) +# @test conj(ss) ≈ tr(t') +# @planar s2 = t[a b; a b] +# @planar t3[a; b] := t[a c; b c] +# @planar s3 = t3[a; a] + +# @test ss ≈ s2 +# @test ss ≈ s3 +# end + +# @timedtestset "Partial trace: test self-consistency" begin +# t = rand(ComplexF64, V3 ⊗ V4 ⊗ V5 ← V3 ⊗ V4 ⊗ V5) # compatible with module fusion +# @planar t2[a; b] := t[c a d; c b d] +# @planar t4[a b; c d] := t[e a b; e c d] +# @planar t5[a; b] := t4[a c; b c] +# @test t2 ≈ t5 +# end + +# @timedtestset "Trace and contraction" begin #TODO: find some version of this that works for off-diagonal case +# t1 = rand(ComplexF64, V3 ⊗ V4 ⊗ V5) +# t2 = rand(ComplexF64, V3 ⊗ V4 ⊗ V5) +# t3 = t1 ⊗ t2' +# # if all(a.i != a.j for a in blocksectors(t3)) +# # replace!(x -> rand(ComplexF64), t3.data) # otherwise full of zeros in off-diagonal case +# # end +# if all(a.i == a.j for a in blocksectors(t3)) +# @planar ta[b; a] := conj(t2[x, a, y]) * t1[x, b, y] # works for diagonal case +# @planar tb[a; b] := t3[x a y; x b y] +# @test ta ≈ tb +# end +# end + +# @timedtestset "Multiplication of isometries: test properties" begin +# W2 = V4 ⊗ V5 +# W1 = W2 ⊗ (rightunitspace(V5) ⊕ rightunitspace(V5)) +# for T in (Float64, ComplexF64) +# t1 = randisometry(T, W1, W2) +# t2 = randisometry(T, W2 ← W2) +# @test isisometric(t1) +# @test isunitary(t2) +# P = t1 * t1' +# @test P * P ≈ P +# end +# end + +# @timedtestset "Multiplication and inverse: test compatibility" begin +# W1 = V1 ⊗ V2 +# W2 = V3 ⊗ V4 ⊗ V5 +# for T in (Float64, ComplexF64) +# t1 = rand(T, W1, W1) +# t2 = rand(T, W2 ← W2) +# t = rand(T, W1, W2) +# @test t1 * (t1 \ t) ≈ t +# @test (t / t2) * t2 ≈ t +# @test t1 \ one(t1) ≈ inv(t1) +# @test one(t1) / t1 ≈ pinv(t1) +# @test_throws SpaceMismatch inv(t) +# @test_throws SpaceMismatch t2 \ t +# @test_throws SpaceMismatch t / t1 +# tp = pinv(t) * t +# @test tp ≈ tp * tp +# end +# end + +# @timedtestset "diag/diagm" begin +# W = V1 ⊗ V2 ← V3 ⊗ V4 ⊗ V5 +# t = randn(ComplexF64, W) +# d = LinearAlgebra.diag(t) +# D = LinearAlgebra.diagm(codomain(t), domain(t), d) +# @test LinearAlgebra.isdiag(D) +# @test LinearAlgebra.diag(D) == d +# end + +# @timedtestset "Sylvester equation" begin +# for T in (Float32, ComplexF64) +# tA = rand(T, V1 ⊗ V2, V1 ⊗ V2) # rewritten for modules +# tB = rand(T, V4 ⊗ V5, V4 ⊗ V5) +# tA = 3 // 2 * leftorth(tA; alg = TK.Polar())[1] +# tB = 1 // 5 * leftorth(tB; alg = TK.Polar())[1] +# tC = rand(T, V1 ⊗ V2, V4 ⊗ V5) +# t = @constinferred sylvester(tA, tB, tC) +# @test codomain(t) == V1 ⊗ V2 +# @test domain(t) == V4 ⊗ V5 +# @test norm(tA * t + t * tB + tC) < +# (norm(tA) + norm(tB) + norm(tC)) * eps(real(T))^(2 / 3) +# end +# end + +# @timedtestset "Tensor product: test via norm preservation" begin # OOMs over here with full spaces +# for T in (Float32, ComplexF64) +# if !isempty(blocksectors(V2 ⊗ V1)) +# t1 = rand(T, V2 ⊗ V3 ⊗ V1, V1 ⊗ V2) +# t2 = rand(T, V2 ⊗ V1 ⊗ V3, V1 ⊗ V1) +# else +# t1 = rand(T, V3 ⊗ V4 ⊗ V5, V1 ⊗ V2) +# t2 = rand(T, V5' ⊗ V4' ⊗ V3', V2' ⊗ V1') +# end +# t = @constinferred (t1 ⊗ t2) +# @test norm(t) ≈ norm(t1) * norm(t2) +# end +# end + +# # TODO: should this test exist? +# @timedtestset "Tensor product: test via tensor contraction" begin +# # W = V3 ⊗ V4 ⊗ V5 ← V1 ⊗ V2 +# W = V4 ← V1 ⊗ V2 # less costly +# for T in (Float32, ComplexF64) +# if !isdiag +# t1 = rand(T, W) +# t2 = rand(T, V4' ← V2' ⊗ V1') +# # t2 = rand(T, V5' ⊗ V4' ⊗ V3', V2' ⊗ V1') # same as previous test +# # @planar t′[1 2 3 6 7 8; 4 5 9 10] := t1[1 2 3; 4 5] * t2[6 7 8; 9 10] +# @planar t′[1 4; 2 3 5 6] := t1[1; 2 3] * t2[4; 5 6] +# else +# t1 = rand(T, V2 ⊗ V3, V1) +# t2 = rand(T, V2, V1 ⊗ V3) +# @planar t′[1 2 4; 3 5 6] := t1[1 2; 3] * t2[4; 5 6] +# end +# t = @constinferred (t1 ⊗ t2) +# @test t ≈ t′ +# end +# end +# end + +# @timedtestset "Tensor absorption" begin +# # absorbing small into large +# if !isempty(blocksectors(V2 ⊗ V3)) +# t1 = zeros(V1 ⊕ V1, V2 ⊗ V3) +# t2 = rand(V1, V2 ⊗ V3) +# else +# t1 = zeros(V1 ⊕ V2, V3 ⊗ V4 ⊗ V5) +# t2 = rand(V1, V3 ⊗ V4 ⊗ V5) +# end +# t3 = @constinferred absorb(t1, t2) +# @test norm(t3) ≈ norm(t2) +# @test norm(t1) == 0 +# t4 = @constinferred absorb!(t1, t2) +# @test t1 === t4 +# @test t3 ≈ t4 + +# # absorbing large into small +# if !isempty(blocksectors(V2 ⊗ V3)) +# t1 = rand(V1 ⊕ V1, V2 ⊗ V3) +# t2 = zeros(V1, V2 ⊗ V3) +# else +# t1 = rand(V1 ⊕ V2, V3 ⊗ V4 ⊗ V5) +# t2 = zeros(V1, V3 ⊗ V4 ⊗ V5) +# end +# t3 = @constinferred absorb(t2, t1) +# @test norm(t3) < norm(t1) +# @test norm(t2) == 0 +# t4 = @constinferred absorb!(t2, t1) +# @test t2 === t4 +# @test t3 ≈ t4 +# end +# end + +# println("---------------------------------------") +# println("Factorizations with symmetry: $Istr") +# println("---------------------------------------") + +# @timedtestset "Factorizations with symmetry involving $Istr ($i, $j)" verbose = true for i in 1:r, j in 1:r +# isdiag = i == j + +# VC = ( +# Vect[I]((i, i, label) => 1 for label in 1:MTK._numlabels(I, i, i)), +# Vect[I](unit(I(i, i, 1)) => 1, rand_object(I, i, i) => 1), # avoids OOMs? +# Vect[I](unit(I(i, i, 1)) => 2, rand_object(I, i, i) => 1), +# Vect[I]((i, i, label) => 1 for label in 1:MTK._numlabels(I, i, i)), +# Vect[I](unit(I(i, i, 1)) => 2, rand_object(I, i, i) => 3), +# ) + +# VM = Vect[I]((i, j, label) => 1 for label in 1:MTK._numlabels(I, i, j)) # all module objects + +# VM1 = ( +# Vect[I](unit(I(i, i, 1)) => 1, rand_object(I, i, i) => 1), # written so V1 ⊗ V2 ← V3 ⊗ V4 ⊗ V5 works +# Vect[I](rand_object(I, i, j) => 2), # generally less blocksectors +# Vect[I](unit(I(i, i, 1)) => 1, rand_object(I, i, i) => 1), +# VM, # important that V4 is module-graded +# Vect[I](unit(I(j, j, 1)) => 2, rand_object(I, j, j) => 1), +# ) + +# VM2 = ( +# Vect[I](rand_object(I, i, j) => 2), # second set where module is V1 here +# Vect[I](unit(I(j, j, 1)) => 1, rand_object(I, j, j) => 1), +# Vect[I](unit(I(i, i, 1)) => 1, rand_object(I, i, i) => 1), +# VM, +# Vect[I](unit(I(j, j, 1)) => 2, rand_object(I, j, j) => 1), +# ) + +# Vs = isdiag ? (VC,) : (VM1, VM2) # avoid duplicate runs + +# # some fail for (2, 2), (3, 3), (6, 6) +# # rightorth RQ(pos) and Polar (fail) for 2nd space +# # leftorth with QL(pos) and Polar for 1st space +# # leftnull QR for 1st space +# # cond and rank leftnull for 1st space + +# # factorization tests require equal objects in blocksectors of domain and codomain, so just put them all +# # FIXME: not sure if still needed +# # VC_all = fill(Vect[I]((i, i, label) => 1 for label in 1:MTK._numlabels(I, i, i)), 5) + +# # VM1_all = (Vect[I](unit(I(i, i, 1)) => 1, rand_object(I, i, i) => 1), +# # VM, +# # Vect[I](unit(I(i, i, 1)) => 1, rand_object(I, i, i) => 1), +# # VM, +# # Vect[I](unit(I(j, j, 1)) => 1, rand_object(I, j, j) => 2) +# # ) + +# # VM2_all = (VM, +# # Vect[I](unit(I(j, j, 1)) => 1, rand_object(I, j, j) => 1), +# # Vect[I](unit(I(i, i, 1)) => 1, rand_object(I, i, i) => 1), +# # VM, +# # Vect[I](unit(I(j, j, 1)) => 2, rand_object(I, j, j) => 2) +# # ) + +# # fact_Vs = (i != j) ? (VM1_all, VM2_all) : (VC_all,) + +# @timedtestset "Factorization" for V in Vs +# V1, V2, V3, V4, V5 = V +# W = V1 ⊗ V2 +# @assert !isempty(blocksectors(W)) +# @assert !isempty(intersect(blocksectors(V4), blocksectors(W))) + +# @testset "QR decomposition" begin +# for T in eltypes, +# t in ( +# rand(T, W, W), rand(T, W, W)', rand(T, W, V4), rand(T, V4, W)', +# DiagonalTensorMap(rand(T, reduceddim(V1)), V1), +# ) + +# Q, R = @constinferred qr_full(t) +# @test Q * R ≈ t +# @test isunitary(Q) + +# Q, R = @constinferred qr_compact(t) +# @test Q * R ≈ t +# @test isisometric(Q) + +# Q, R = @constinferred left_orth(t) +# @test Q * R ≈ t +# @test isisometric(Q) + +# N = @constinferred qr_null(t) +# @test isisometric(N) +# @test norm(N' * t) ≈ 0 atol = 100 * eps(norm(t)) + +# N = @constinferred left_null(t) +# @test isisometric(N) +# @test norm(N' * t) ≈ 0 atol = 100 * eps(norm(t)) +# end + +# # empty tensor +# for T in eltypes +# t = rand(T, V1 ⊗ V2, zerospace(V1)) + +# Q, R = @constinferred qr_full(t) +# @test Q * R ≈ t +# @test isunitary(Q) +# @test dim(R) == dim(t) == 0 + +# Q, R = @constinferred qr_compact(t) +# @test Q * R ≈ t +# @test isisometric(Q) +# @test dim(Q) == dim(R) == dim(t) + +# Q, R = @constinferred left_orth(t) +# @test Q * R ≈ t +# @test isisometric(Q) +# @test dim(Q) == dim(R) == dim(t) + +# N = @constinferred qr_null(t) +# @test isunitary(N) +# @test norm(N' * t) ≈ 0 atol = 100 * eps(norm(t)) +# end +# end + +# @testset "LQ decomposition" begin +# for T in eltypes, +# t in ( +# rand(T, W, W), rand(T, W, W)', rand(T, W, V4), rand(T, V4, W)', +# DiagonalTensorMap(rand(T, reduceddim(V1)), V1), +# ) + +# L, Q = @constinferred lq_full(t) +# @test L * Q ≈ t +# @test isunitary(Q) + +# L, Q = @constinferred lq_compact(t) +# @test L * Q ≈ t +# @test isisometric(Q; side = :right) + +# L, Q = @constinferred right_orth(t) +# @test L * Q ≈ t +# @test isisometric(Q; side = :right) + +# Nᴴ = @constinferred lq_null(t) +# @test isisometric(Nᴴ; side = :right) +# @test norm(t * Nᴴ') ≈ 0 atol = 100 * eps(norm(t)) +# end + +# for T in eltypes +# # empty tensor +# t = rand(T, zerospace(V1), V1 ⊗ V2) + +# L, Q = @constinferred lq_full(t) +# @test L * Q ≈ t +# @test isunitary(Q) +# @test dim(L) == dim(t) == 0 + +# L, Q = @constinferred lq_compact(t) +# @test L * Q ≈ t +# @test isisometric(Q; side = :right) +# @test dim(Q) == dim(L) == dim(t) + +# L, Q = @constinferred right_orth(t) +# @test L * Q ≈ t +# @test isisometric(Q; side = :right) +# @test dim(Q) == dim(L) == dim(t) + +# Nᴴ = @constinferred lq_null(t) +# @test isunitary(Nᴴ) +# @test norm(t * Nᴴ') ≈ 0 atol = 100 * eps(norm(t)) +# end +# end + +# @testset "Polar decomposition" begin +# for T in eltypes, +# t in ( +# rand(T, W, W), rand(T, W, W)', rand(T, W, V4), rand(T, V4, W)', +# DiagonalTensorMap(rand(T, reduceddim(V1)), V1), +# ) + +# @assert domain(t) ≾ codomain(t) +# w, p = @constinferred left_polar(t) +# @test w * p ≈ t +# @test isisometric(w) +# @test isposdef(p) + +# w, p = @constinferred left_orth(t; alg = :polar) +# @test w * p ≈ t +# @test isisometric(w) +# end + +# for T in eltypes, +# t in (rand(T, W, W), rand(T, W, W)', rand(T, V4, W), rand(T, W, V4)') + +# @assert codomain(t) ≾ domain(t) +# p, wᴴ = @constinferred right_polar(t) +# @test p * wᴴ ≈ t +# @test isisometric(wᴴ; side = :right) +# @test isposdef(p) + +# p, wᴴ = @constinferred right_orth(t; alg = :polar) +# @test p * wᴴ ≈ t +# @test isisometric(wᴴ; side = :right) +# end +# end + +# @testset "SVD" begin +# for T in eltypes, +# t in ( +# rand(T, W, W), rand(T, W, W)', +# rand(T, W, V4), rand(T, V4, W), +# rand(T, W, V4)', rand(T, V4, W)', +# DiagonalTensorMap(rand(T, reduceddim(V1)), V1), +# ) + +# u, s, vᴴ = @constinferred svd_full(t) +# @test u * s * vᴴ ≈ t +# @test isunitary(u) +# @test isunitary(vᴴ) + +# u, s, vᴴ = @constinferred svd_compact(t) +# @test u * s * vᴴ ≈ t +# @test isisometric(u) +# @test isposdef(s) +# @test isisometric(vᴴ; side = :right) + +# s′ = @constinferred svd_vals(t) +# @test s′ ≈ diagview(s) +# @test s′ isa TensorKit.SectorVector + +# v, c = @constinferred left_orth(t; alg = :svd) +# @test v * c ≈ t +# @test isisometric(v) + +# c, vᴴ = @constinferred right_orth(t; alg = :svd) +# @test c * vᴴ ≈ t +# @test isisometric(vᴴ; side = :right) + +# N = @constinferred left_null(t; alg = :svd) +# @test isisometric(N) +# @test norm(N' * t) ≈ 0 atol = 100 * eps(norm(t)) + +# N = @constinferred left_null(t; trunc = (; atol = 100 * eps(norm(t)))) +# @test isisometric(N) +# @test norm(N' * t) ≈ 0 atol = 100 * eps(norm(t)) + +# Nᴴ = @constinferred right_null(t; alg = :svd) +# @test isisometric(Nᴴ; side = :right) +# @test norm(t * Nᴴ') ≈ 0 atol = 100 * eps(norm(t)) + +# Nᴴ = @constinferred right_null(t; trunc = (; atol = 100 * eps(norm(t)))) +# @test isisometric(Nᴴ; side = :right) +# @test norm(t * Nᴴ') ≈ 0 atol = 100 * eps(norm(t)) +# end + +# # empty tensor +# for T in eltypes, t in (rand(T, W, zerospace(V1)), rand(T, zerospace(V1), W)) +# U, S, Vᴴ = @constinferred svd_full(t) +# @test U * S * Vᴴ ≈ t +# @test isunitary(U) +# @test isunitary(Vᴴ) + +# U, S, Vᴴ = @constinferred svd_compact(t) +# @test U * S * Vᴴ ≈ t +# @test dim(U) == dim(S) == dim(Vᴴ) == dim(t) == 0 +# end +# end + +# @testset "truncated SVD" begin +# for T in eltypes, +# t in ( +# randn(T, W, W), randn(T, W, W)', +# randn(T, W, V4), randn(T, V4, W), +# randn(T, W, V4)', randn(T, V4, W)', +# DiagonalTensorMap(randn(T, reduceddim(V1)), V1), +# ) + +# @constinferred normalize!(t) + +# U, S, Vᴴ, ϵ = @constinferred svd_trunc(t; trunc = notrunc()) +# @test U * S * Vᴴ ≈ t +# @test ϵ ≈ 0 +# @test isisometric(U) +# @test isisometric(Vᴴ; side = :right) + +# # dimension of S is a float for IsingBimodule +# nvals = round(Int, dim(domain(S)) / 2) +# trunc = truncrank(nvals) +# U1, S1, Vᴴ1, ϵ1 = @constinferred svd_trunc(t; trunc) +# @test t * Vᴴ1' ≈ U1 * S1 +# @test isisometric(U1) +# @test isisometric(Vᴴ1; side = :right) +# @test norm(t - U1 * S1 * Vᴴ1) ≈ ϵ1 atol = eps(real(T))^(4 / 5) +# @test dim(domain(S1)) <= nvals + +# λ = minimum(diagview(S1)) +# trunc = trunctol(; atol = λ - 10eps(λ)) +# U2, S2, Vᴴ2, ϵ2 = @constinferred svd_trunc(t; trunc) +# @test t * Vᴴ2' ≈ U2 * S2 +# @test isisometric(U2) +# @test isisometric(Vᴴ2; side = :right) +# @test norm(t - U2 * S2 * Vᴴ2) ≈ ϵ2 atol = eps(real(T))^(4 / 5) +# @test minimum(diagview(S1)) >= λ +# @test U2 ≈ U1 +# @test S2 ≈ S1 +# @test Vᴴ2 ≈ Vᴴ1 +# @test ϵ1 ≈ ϵ2 + +# trunc = truncspace(space(S2, 1)) +# U3, S3, Vᴴ3, ϵ3 = @constinferred svd_trunc(t; trunc) +# @test t * Vᴴ3' ≈ U3 * S3 +# @test isisometric(U3) +# @test isisometric(Vᴴ3; side = :right) +# @test norm(t - U3 * S3 * Vᴴ3) ≈ ϵ3 atol = eps(real(T))^(4 / 5) +# @test space(S3, 1) ≾ space(S2, 1) + +# for trunc in (truncerror(; atol = ϵ2), truncerror(; rtol = ϵ2 / norm(t))) +# U4, S4, Vᴴ4, ϵ4 = @constinferred svd_trunc(t; trunc) +# @test t * Vᴴ4' ≈ U4 * S4 +# @test isisometric(U4) +# @test isisometric(Vᴴ4; side = :right) +# @test norm(t - U4 * S4 * Vᴴ4) ≈ ϵ4 atol = eps(real(T))^(4 / 5) +# @test ϵ4 ≤ ϵ2 +# end + +# trunc = truncrank(nvals) & trunctol(; atol = λ - 10eps(λ)) +# U5, S5, Vᴴ5, ϵ5 = @constinferred svd_trunc(t; trunc) +# @test t * Vᴴ5' ≈ U5 * S5 +# @test isisometric(U5) +# @test isisometric(Vᴴ5; side = :right) +# @test norm(t - U5 * S5 * Vᴴ5) ≈ ϵ5 atol = eps(real(T))^(4 / 5) +# @test minimum(diagview(S5)) >= λ +# @test dim(domain(S5)) ≤ nvals +# end +# end + +# @testset "Eigenvalue decomposition" begin +# for T in eltypes, +# t in ( +# rand(T, V1, V1), rand(T, W, W), rand(T, W, W)', +# DiagonalTensorMap(rand(T, reduceddim(V1)), V1), +# ) + +# d, v = @constinferred eig_full(t) +# @test t * v ≈ v * d + +# d′ = @constinferred eig_vals(t) +# @test d′ ≈ diagview(d) +# @test d′ isa TensorKit.SectorVector + +# vdv = project_hermitian!(v' * v) +# @test @constinferred isposdef(vdv) +# t isa DiagonalTensorMap || @test !isposdef(t) # unlikely for non-hermitian map + +# nvals = round(Int, dim(domain(t)) / 2) +# d, v = @constinferred eig_trunc(t; trunc = truncrank(nvals)) +# @test t * v ≈ v * d +# @test dim(domain(d)) ≤ nvals + +# t2 = @constinferred project_hermitian(t) +# D, V = eigen(t2) +# @test isisometric(V) +# D̃, Ṽ = @constinferred eigh_full(t2) +# @test D ≈ D̃ +# @test V ≈ Ṽ +# λ = minimum(real, diagview(D)) +# @test cond(Ṽ) ≈ one(real(T)) +# @test isposdef(t2) == isposdef(λ) +# @test isposdef(t2 - λ * one(t2) + 0.1 * one(t2)) +# @test !isposdef(t2 - λ * one(t2) - 0.1 * one(t2)) + +# d, v = @constinferred eigh_full(t2) +# @test t2 * v ≈ v * d +# @test isunitary(v) + +# d′ = @constinferred eigh_vals(t2) +# @test d′ ≈ diagview(d) +# @test d′ isa TensorKit.SectorVector + +# λ = minimum(real, diagview(d)) +# @test cond(v) ≈ one(real(T)) +# @test isposdef(t2) == isposdef(λ) +# @test isposdef(t2 - λ * one(t) + 0.1 * one(t2)) +# @test !isposdef(t2 - λ * one(t) - 0.1 * one(t2)) + +# d, v = @constinferred eigh_trunc(t2; trunc = truncrank(nvals)) +# @test t2 * v ≈ v * d +# @test dim(domain(d)) ≤ nvals +# end +# end + +# @testset "Condition number and rank" begin +# for T in eltypes, +# t in ( +# rand(T, W, W), rand(T, W, W)', +# rand(T, W, V4), rand(T, V4, W), +# rand(T, W, V4)', rand(T, V4, W)', +# DiagonalTensorMap(rand(T, reduceddim(V1)), V1), +# ) + +# d1, d2 = dim(codomain(t)), dim(domain(t)) +# r = rank(t) +# @test r == min(d1, d2) +# @test typeof(r) == typeof(d1) +# M = left_null(t) +# @test @constinferred(rank(M)) + r ≈ d1 +# Mᴴ = right_null(t) +# @test rank(Mᴴ) + r ≈ d2 +# end +# for T in eltypes +# u = unitary(T, V1 ⊗ V2, V1 ⊗ V2) +# @test @constinferred(cond(u)) ≈ one(real(T)) +# @test @constinferred(rank(u)) == dim(V1 ⊗ V2) + +# t = rand(T, zerospace(V1), W) +# @test rank(t) == 0 +# t2 = rand(T, zerospace(V1) * zerospace(V2), zerospace(V1) * zerospace(V2)) +# @test rank(t2) == 0 +# @test cond(t2) == 0.0 +# end +# for T in eltypes, t in (rand(T, W, W), rand(T, W, W)') +# project_hermitian!(t) +# vals = @constinferred LinearAlgebra.eigvals(t) +# λmax = maximum(s -> maximum(abs, s), values(vals)) +# λmin = minimum(s -> minimum(abs, s), values(vals)) +# @test cond(t) ≈ λmax / λmin +# end +# end + +# @testset "Hermitian projections" begin +# for T in eltypes, +# t in ( +# rand(T, V1, V1), rand(T, W, W), rand(T, W, W)', +# DiagonalTensorMap(rand(T, reduceddim(V1)), V1), +# ) +# normalize!(t) +# noisefactor = eps(real(T))^(3 / 4) + +# th = (t + t') / 2 +# ta = (t - t') / 2 +# tc = copy(t) + +# th′ = @constinferred project_hermitian(t) +# @test ishermitian(th′) +# @test th′ ≈ th +# @test t == tc +# th_approx = th + noisefactor * ta +# @test !ishermitian(th_approx) || (T <: Real && t isa DiagonalTensorMap) +# @test ishermitian(th_approx; atol = 10 * noisefactor) + +# ta′ = project_antihermitian(t) +# @test isantihermitian(ta′) +# @test ta′ ≈ ta +# @test t == tc +# ta_approx = ta + noisefactor * th +# @test !isantihermitian(ta_approx) +# @test isantihermitian(ta_approx; atol = 10 * noisefactor) || (T <: Real && t isa DiagonalTensorMap) +# end +# end + +# @testset "Isometric projections" begin +# for T in eltypes, +# t in ( +# randn(T, W, W), randn(T, W, W)', +# randn(T, W, V4), randn(T, V4, W)', +# ) +# t2 = project_isometric(t) +# @test isisometric(t2) +# t3 = project_isometric(t2) +# @test t3 ≈ t2 # stability of the projection +# @test t2 * (t2' * t) ≈ t + +# tc = similar(t) +# t3 = @constinferred project_isometric!(copy!(tc, t), t2) +# @test t3 === t2 +# @test isisometric(t2) + +# # test that t2 is closer to A then any other isometry +# for k in 1:10 +# δt = randn!(similar(t)) +# t3 = project_isometric(t + δt / 100) +# @test norm(t - t3) > norm(t - t2) +# end +# end +# end +# end +# end diff --git a/test/test_aqua.jl b/test/test_aqua.jl index 14e5b7c..f57c8ba 100644 --- a/test/test_aqua.jl +++ b/test/test_aqua.jl @@ -3,5 +3,5 @@ using Aqua: Aqua using Test: @testset @testset "Code quality (Aqua.jl)" begin - Aqua.test_all(MultiTensorKit) + Aqua.test_all(MultiTensorKit) end