diff --git a/.JuliaFormatter.toml b/.JuliaFormatter.toml deleted file mode 100644 index 4c49a86..0000000 --- a/.JuliaFormatter.toml +++ /dev/null @@ -1,3 +0,0 @@ -# See https://domluna.github.io/JuliaFormatter.jl/stable/ for a list of options -style = "blue" -indent = 2 diff --git a/.github/workflows/CompatHelper.yml b/.github/workflows/CompatHelper.yml index 456fa05..0614de9 100644 --- a/.github/workflows/CompatHelper.yml +++ b/.github/workflows/CompatHelper.yml @@ -2,7 +2,7 @@ name: "CompatHelper" on: schedule: - - cron: 0 0 * * * + - cron: '0 0 * * *' workflow_dispatch: permissions: contents: write diff --git a/.github/workflows/FormatCheck.yml b/.github/workflows/FormatCheck.yml index 3f78afc..1525861 100644 --- a/.github/workflows/FormatCheck.yml +++ b/.github/workflows/FormatCheck.yml @@ -1,11 +1,14 @@ name: "Format Check" on: - push: - branches: - - 'main' - tags: '*' - pull_request: + pull_request_target: + paths: ['**/*.jl'] + types: [opened, synchronize, reopened, ready_for_review] + +permissions: + contents: read + actions: write + pull-requests: write jobs: format-check: diff --git a/.gitignore b/.gitignore index 10593a9..7085ca8 100644 --- a/.gitignore +++ b/.gitignore @@ -9,6 +9,10 @@ .vscode/ Manifest.toml benchmark/*.json +dev/ +docs/LocalPreferences.toml docs/Manifest.toml docs/build/ docs/src/index.md +examples/LocalPreferences.toml +test/LocalPreferences.toml diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 88bc8b4..3fc4743 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -1,5 +1,5 @@ ci: - skip: [julia-formatter] + skip: [runic] repos: - repo: https://github.com/pre-commit/pre-commit-hooks @@ -11,7 +11,7 @@ repos: - id: end-of-file-fixer exclude_types: [markdown] # incompatible with Literate.jl -- repo: "https://github.com/domluna/JuliaFormatter.jl" - rev: v2.1.6 +- repo: https://github.com/fredrikekre/runic-pre-commit + rev: v2.0.1 hooks: - - id: "julia-formatter" + - id: runic diff --git a/Project.toml b/Project.toml index 46a935f..3e3d339 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "FusionTensors" uuid = "e16ca583-1f51-4df0-8e12-57d32947d33e" authors = ["ITensor developers and contributors"] -version = "0.5.11" +version = "0.5.12" [deps] Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697" diff --git a/docs/make.jl b/docs/make.jl index 5aaf7e8..643c700 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,22 +1,22 @@ using FusionTensors: FusionTensors using Documenter: Documenter, DocMeta, deploydocs, makedocs -DocMeta.setdocmeta!(FusionTensors, :DocTestSetup, :(using FusionTensors); recursive=true) +DocMeta.setdocmeta!(FusionTensors, :DocTestSetup, :(using FusionTensors); recursive = true) include("make_index.jl") makedocs(; - modules=[FusionTensors], - authors="ITensor developers and contributors", - sitename="FusionTensors.jl", - format=Documenter.HTML(; - canonical="https://itensor.github.io/FusionTensors.jl", - edit_link="main", - assets=["assets/favicon.ico", "assets/extras.css"], - ), - pages=["Home" => "index.md", "Reference" => "reference.md"], + modules = [FusionTensors], + authors = "ITensor developers and contributors", + sitename = "FusionTensors.jl", + format = Documenter.HTML(; + canonical = "https://itensor.github.io/FusionTensors.jl", + edit_link = "main", + assets = ["assets/favicon.ico", "assets/extras.css"], + ), + pages = ["Home" => "index.md", "Reference" => "reference.md"], ) deploydocs(; - repo="github.com/ITensor/FusionTensors.jl", devbranch="main", push_preview=true + repo = "github.com/ITensor/FusionTensors.jl", devbranch = "main", push_preview = true ) diff --git a/docs/make_index.jl b/docs/make_index.jl index 2ae72dc..7ae18d3 100644 --- a/docs/make_index.jl +++ b/docs/make_index.jl @@ -2,20 +2,20 @@ using Literate: Literate using FusionTensors: FusionTensors function ccq_logo(content) - include_ccq_logo = """ + include_ccq_logo = """ ```@raw html Flatiron Center for Computational Quantum Physics logo. Flatiron Center for Computational Quantum Physics logo. ``` """ - content = replace(content, "{CCQ_LOGO}" => include_ccq_logo) - return content + content = replace(content, "{CCQ_LOGO}" => include_ccq_logo) + return content end Literate.markdown( - joinpath(pkgdir(FusionTensors), "examples", "README.jl"), - joinpath(pkgdir(FusionTensors), "docs", "src"); - flavor=Literate.DocumenterFlavor(), - name="index", - postprocess=ccq_logo, + joinpath(pkgdir(FusionTensors), "examples", "README.jl"), + joinpath(pkgdir(FusionTensors), "docs", "src"); + flavor = Literate.DocumenterFlavor(), + name = "index", + postprocess = ccq_logo, ) diff --git a/docs/make_readme.jl b/docs/make_readme.jl index bb48eba..ce51321 100644 --- a/docs/make_readme.jl +++ b/docs/make_readme.jl @@ -2,20 +2,20 @@ using Literate: Literate using FusionTensors: FusionTensors function ccq_logo(content) - include_ccq_logo = """ + include_ccq_logo = """ Flatiron Center for Computational Quantum Physics logo. """ - content = replace(content, "{CCQ_LOGO}" => include_ccq_logo) - return content + content = replace(content, "{CCQ_LOGO}" => include_ccq_logo) + return content end Literate.markdown( - joinpath(pkgdir(FusionTensors), "examples", "README.jl"), - joinpath(pkgdir(FusionTensors)); - flavor=Literate.CommonMarkFlavor(), - name="README", - postprocess=ccq_logo, + joinpath(pkgdir(FusionTensors), "examples", "README.jl"), + joinpath(pkgdir(FusionTensors)); + flavor = Literate.CommonMarkFlavor(), + name = "README", + postprocess = ccq_logo, ) diff --git a/examples/README.jl b/examples/README.jl index 42cb0ae..7a57397 100644 --- a/examples/README.jl +++ b/examples/README.jl @@ -1,5 +1,5 @@ # # FusionTensors.jl -# +# # [![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://itensor.github.io/FusionTensors.jl/stable/) # [![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://itensor.github.io/FusionTensors.jl/dev/) # [![Build Status](https://github.com/ITensor/FusionTensors.jl/actions/workflows/Tests.yml/badge.svg?branch=main)](https://github.com/ITensor/FusionTensors.jl/actions/workflows/Tests.yml?query=branch%3Amain) diff --git a/src/fusion_trees/clebsch_gordan_tensors.jl b/src/fusion_trees/clebsch_gordan_tensors.jl index 37575d4..3815564 100644 --- a/src/fusion_trees/clebsch_gordan_tensors.jl +++ b/src/fusion_trees/clebsch_gordan_tensors.jl @@ -6,132 +6,132 @@ using HalfIntegers: half using WignerSymbols: clebschgordan using GradedArrays: - AbelianStyle, - AbstractSector, - NotAbelianStyle, - O2, - SU, - SymmetryStyle, - dual, - istrivial, - quantum_dimension, - sector_label, - sectors, - trivial, - zero_odd + AbelianStyle, + AbstractSector, + NotAbelianStyle, + O2, + SU, + SymmetryStyle, + dual, + istrivial, + quantum_dimension, + sector_label, + sectors, + trivial, + zero_odd using TensorAlgebra: contract using TensorProducts: ⊗ function symbol_1j(s::AbstractSector) - cgt = clebsch_gordan_tensor(s, dual(s), trivial(s), 1) - return sqrt(quantum_dimension(s)) * cgt[:, :, 1] + cgt = clebsch_gordan_tensor(s, dual(s), trivial(s), 1) + return sqrt(quantum_dimension(s)) * cgt[:, :, 1] end function clebsch_gordan_tensor( - s1::AbstractSector, - s2::AbstractSector, - s3::AbstractSector, - arrow1::Bool, - arrow2::Bool, - inner_mult_index::Int, -) - cgt = clebsch_gordan_tensor(s1, s2, s3, inner_mult_index) - if arrow1 - flip1 = symbol_1j(s1) - cgt = contract((1, 2, 3), flip1, (4, 1), cgt, (4, 2, 3)) - end - if arrow2 - flip2 = symbol_1j(s2) - cgt = contract((1, 2, 3), flip2, (4, 2), cgt, (1, 4, 3)) - end - return cgt + s1::AbstractSector, + s2::AbstractSector, + s3::AbstractSector, + arrow1::Bool, + arrow2::Bool, + inner_mult_index::Int, + ) + cgt = clebsch_gordan_tensor(s1, s2, s3, inner_mult_index) + if arrow1 + flip1 = symbol_1j(s1) + cgt = contract((1, 2, 3), flip1, (4, 1), cgt, (4, 2, 3)) + end + if arrow2 + flip2 = symbol_1j(s2) + cgt = contract((1, 2, 3), flip2, (4, 2), cgt, (1, 4, 3)) + end + return cgt end function clebsch_gordan_tensor(s1::S, s2::S, s3::S, outer_mult_index::Int) where {S} - return clebsch_gordan_tensor(SymmetryStyle(S), s1, s2, s3, outer_mult_index) + return clebsch_gordan_tensor(SymmetryStyle(S), s1, s2, s3, outer_mult_index) end function clebsch_gordan_tensor( - ::AbelianStyle, s1::S, s2::S, s3::S, outer_mult_index::Int -) where {S} - @assert outer_mult_index == 1 - return s1 ⊗ s2 == s3 ? ones((1, 1, 1)) : zeros((1, 1, 1)) + ::AbelianStyle, s1::S, s2::S, s3::S, outer_mult_index::Int + ) where {S} + @assert outer_mult_index == 1 + return s1 ⊗ s2 == s3 ? ones((1, 1, 1)) : zeros((1, 1, 1)) end function clebsch_gordan_tensor(::NotAbelianStyle, s1::O2, s2::O2, s3::O2, ::Int) - return clebsch_gordan_tensor(s1, s2, s3) # no outer multiplicity + return clebsch_gordan_tensor(s1, s2, s3) # no outer multiplicity end function clebsch_gordan_tensor(s1::O2, s2::O2, s3::O2) - d1 = quantum_dimension(s1) - d2 = quantum_dimension(s2) - d3 = quantum_dimension(s3) - cgt = zeros((d1, d2, d3)) - s3 ∉ sectors(s1 ⊗ s2) && return cgt + d1 = quantum_dimension(s1) + d2 = quantum_dimension(s2) + d3 = quantum_dimension(s3) + cgt = zeros((d1, d2, d3)) + s3 ∉ sectors(s1 ⊗ s2) && return cgt - # adapted from TensorKit - l1 = sector_label(s1) - l2 = sector_label(s2) - l3 = sector_label(s3) - if l3 <= 0 # 0even or 0odd - if l1 <= 0 && l2 <= 0 - cgt[1, 1, 1, 1] = 1.0 - else - if istrivial(s3) - cgt[1, 2, 1, 1] = 1.0 / sqrt(2) - cgt[2, 1, 1, 1] = 1.0 / sqrt(2) - else - cgt[1, 2, 1, 1] = 1.0 / sqrt(2) - cgt[2, 1, 1, 1] = -1.0 / sqrt(2) - end + # adapted from TensorKit + l1 = sector_label(s1) + l2 = sector_label(s2) + l3 = sector_label(s3) + if l3 <= 0 # 0even or 0odd + if l1 <= 0 && l2 <= 0 + cgt[1, 1, 1, 1] = 1.0 + else + if istrivial(s3) + cgt[1, 2, 1, 1] = 1.0 / sqrt(2) + cgt[2, 1, 1, 1] = 1.0 / sqrt(2) + else + cgt[1, 2, 1, 1] = 1.0 / sqrt(2) + cgt[2, 1, 1, 1] = -1.0 / sqrt(2) + end + end + elseif l1 <= 0 # 0even or 0odd + cgt[1, 1, 1, 1] = 1.0 + cgt[1, 2, 2, 1] = s1 == zero_odd(O2) ? -1.0 : 1.0 + elseif l2 == 0 + cgt[1, 1, 1, 1] = 1.0 + cgt[2, 1, 2, 1] = s2 == zero_odd(O2) ? -1.0 : 1.0 + elseif l3 == l1 + l2 + cgt[1, 1, 1, 1] = 1.0 + cgt[2, 2, 2, 1] = 1.0 + elseif l3 == l1 - l2 + cgt[1, 2, 1, 1] = 1.0 + cgt[2, 1, 2, 1] = 1.0 + elseif l3 == l2 - l1 + cgt[2, 1, 1, 1] = 1.0 + cgt[1, 2, 2, 1] = 1.0 end - elseif l1 <= 0 # 0even or 0odd - cgt[1, 1, 1, 1] = 1.0 - cgt[1, 2, 2, 1] = s1 == zero_odd(O2) ? -1.0 : 1.0 - elseif l2 == 0 - cgt[1, 1, 1, 1] = 1.0 - cgt[2, 1, 2, 1] = s2 == zero_odd(O2) ? -1.0 : 1.0 - elseif l3 == l1 + l2 - cgt[1, 1, 1, 1] = 1.0 - cgt[2, 2, 2, 1] = 1.0 - elseif l3 == l1 - l2 - cgt[1, 2, 1, 1] = 1.0 - cgt[2, 1, 2, 1] = 1.0 - elseif l3 == l2 - l1 - cgt[2, 1, 1, 1] = 1.0 - cgt[1, 2, 2, 1] = 1.0 - end - return cgt + return cgt end function clebsch_gordan_tensor(::NotAbelianStyle, s1::SU{2}, s2::SU{2}, s3::SU{2}, ::Int) - return clebsch_gordan_tensor(s1, s2, s3) # no outer multiplicity + return clebsch_gordan_tensor(s1, s2, s3) # no outer multiplicity end function clebsch_gordan_tensor(s1::SU{2}, s2::SU{2}, s3::SU{2}) - d1 = quantum_dimension(s1) - d2 = quantum_dimension(s2) - d3 = quantum_dimension(s3) - j1 = half(d1 - 1) - j2 = half(d2 - 1) - j3 = half(d3 - 1) - cgtensor = Array{Float64,3}(undef, (d1, d2, d3)) - for (i, j, k) in Iterators.product(1:d1, 1:d2, 1:d3) - m1 = j1 - i + 1 - m2 = j2 - j + 1 - m3 = j3 - k + 1 - cgtensor[i, j, k] = clebschgordan(j1, m1, j2, m2, j3, m3) - end - return cgtensor + d1 = quantum_dimension(s1) + d2 = quantum_dimension(s2) + d3 = quantum_dimension(s3) + j1 = half(d1 - 1) + j2 = half(d2 - 1) + j3 = half(d3 - 1) + cgtensor = Array{Float64, 3}(undef, (d1, d2, d3)) + for (i, j, k) in Iterators.product(1:d1, 1:d2, 1:d3) + m1 = j1 - i + 1 + m2 = j2 - j + 1 + m3 = j3 - k + 1 + cgtensor[i, j, k] = clebschgordan(j1, m1, j2, m2, j3, m3) + end + return cgtensor end function clebsch_gordan_tensor( - ::NotAbelianStyle, s1::SU{3}, s2::SU{3}, s3::SU{3}, outer_mult_index::Int -) - d1 = quantum_dimension(s1) - d2 = quantum_dimension(s2) - d3 = quantum_dimension(s3) - cgtensor = zeros(d1, d2, d3) - # dummy - return cgtensor + ::NotAbelianStyle, s1::SU{3}, s2::SU{3}, s3::SU{3}, outer_mult_index::Int + ) + d1 = quantum_dimension(s1) + d2 = quantum_dimension(s2) + d3 = quantum_dimension(s3) + cgtensor = zeros(d1, d2, d3) + # dummy + return cgtensor end diff --git a/src/fusion_trees/fusiontree.jl b/src/fusion_trees/fusiontree.jl index 461c8e2..b448d6f 100644 --- a/src/fusion_trees/fusiontree.jl +++ b/src/fusion_trees/fusiontree.jl @@ -4,21 +4,21 @@ # compatibility with TensorKit conventions? using GradedArrays: - GradedArrays, - ×, - AbstractGradedUnitRange, - AbstractSector, - SectorProduct, - arguments, - flip, - flip_dual, - isdual, - nsymbol, - sector_multiplicities, - sector_type, - sectors, - to_gradedrange, - trivial + GradedArrays, + ×, + AbstractGradedUnitRange, + AbstractSector, + SectorProduct, + arguments, + flip, + flip_dual, + isdual, + nsymbol, + sector_multiplicities, + sector_type, + sectors, + to_gradedrange, + trivial using TensorAlgebra: flatten_tuples using TensorProducts: ⊗ @@ -49,27 +49,27 @@ using TensorProducts: ⊗ # The interface uses AbstractGradedArrays as input for interface simplicity # however only sectors are used and blocklengths are never read. -struct SectorFusionTree{S,N,M} - leaves::NTuple{N,S} # TBD rename outer_sectors or leave_sectors? - arrows::NTuple{N,Bool} - root_sector::S - branch_sectors::NTuple{M,S} # M = N-1 - outer_multiplicity_indices::NTuple{M,Int} # M = N-1 - - # TBD could have branch_sectors with length N-2 - # currently first(branch_sectors) == first(leaves) - # redundant but allows for simpler, generic grow_tree code - - function SectorFusionTree( - leaves, arrows, root_sector, branch_sectors, outer_multiplicity_indices - ) - N = length(leaves) - @assert length(branch_sectors) == max(0, N - 1) - @assert length(outer_multiplicity_indices) == max(0, N - 1) - return new{typeof(root_sector),length(leaves),length(branch_sectors)}( - leaves, arrows, root_sector, branch_sectors, outer_multiplicity_indices - ) - end +struct SectorFusionTree{S, N, M} + leaves::NTuple{N, S} # TBD rename outer_sectors or leave_sectors? + arrows::NTuple{N, Bool} + root_sector::S + branch_sectors::NTuple{M, S} # M = N-1 + outer_multiplicity_indices::NTuple{M, Int} # M = N-1 + + # TBD could have branch_sectors with length N-2 + # currently first(branch_sectors) == first(leaves) + # redundant but allows for simpler, generic grow_tree code + + function SectorFusionTree( + leaves, arrows, root_sector, branch_sectors, outer_multiplicity_indices + ) + N = length(leaves) + @assert length(branch_sectors) == max(0, N - 1) + @assert length(outer_multiplicity_indices) == max(0, N - 1) + return new{typeof(root_sector), length(leaves), length(branch_sectors)}( + leaves, arrows, root_sector, branch_sectors, outer_multiplicity_indices + ) + end end # getters @@ -83,83 +83,83 @@ outer_multiplicity_indices(f::SectorFusionTree) = f.outer_multiplicity_indices Base.convert(T::Type{<:Array}, f::SectorFusionTree) = convert(T, to_array(f)) function Base.isless(f1::SectorFusionTree, f2::SectorFusionTree) - return isless(leaves(f1), leaves(f2)) || - isless(arrows(f1), arrows(f2)) || - isless(root_sector(f1), root_sector(f2)) || - isless(branch_sectors(f1), branch_sectors(f2)) || - isless(outer_multiplicity_indices(f1), outer_multiplicity_indices(f2)) + return isless(leaves(f1), leaves(f2)) || + isless(arrows(f1), arrows(f2)) || + isless(root_sector(f1), root_sector(f2)) || + isless(branch_sectors(f1), branch_sectors(f2)) || + isless(outer_multiplicity_indices(f1), outer_multiplicity_indices(f2)) end -Base.length(::SectorFusionTree{<:Any,N}) where {N} = N +Base.length(::SectorFusionTree{<:Any, N}) where {N} = N # GradedArrays interface GradedArrays.sector_type(::Type{<:SectorFusionTree{S}}) where {S} = S function GradedArrays.flip(f::SectorFusionTree) - return SectorFusionTree( - dual.(leaves(f)), - .!arrows(f), - dual(root_sector(f)), - dual.(branch_sectors(f)), - outer_multiplicity_indices(f), - ) + return SectorFusionTree( + dual.(leaves(f)), + .!arrows(f), + dual(root_sector(f)), + dual.(branch_sectors(f)), + outer_multiplicity_indices(f), + ) end function GradedArrays.:×(f1::SectorFusionTree, f2::SectorFusionTree) - @assert arrows(f1) == arrows(f2) - product_leaves = .×(leaves(f1), leaves(f2)) - product_root_sector = root_sector(f1) × root_sector(f2) - product_branch_sectors = .×(branch_sectors(f1), branch_sectors(f2)) - product_outer_multiplicity_indices = outer_multiplicity_kron.( - Base.tail(leaves(f1)), - branch_sectors(f1), - (Base.tail(branch_sectors(f1))..., root_sector(f1)), - outer_multiplicity_indices(f1), - outer_multiplicity_indices(f2), - ) - return SectorFusionTree( - product_leaves, - arrows(f1), - product_root_sector, - product_branch_sectors, - product_outer_multiplicity_indices, - ) + @assert arrows(f1) == arrows(f2) + product_leaves = .×(leaves(f1), leaves(f2)) + product_root_sector = root_sector(f1) × root_sector(f2) + product_branch_sectors = .×(branch_sectors(f1), branch_sectors(f2)) + product_outer_multiplicity_indices = outer_multiplicity_kron.( + Base.tail(leaves(f1)), + branch_sectors(f1), + (Base.tail(branch_sectors(f1))..., root_sector(f1)), + outer_multiplicity_indices(f1), + outer_multiplicity_indices(f2), + ) + return SectorFusionTree( + product_leaves, + arrows(f1), + product_root_sector, + product_branch_sectors, + product_outer_multiplicity_indices, + ) end function GradedArrays.arguments(f::SectorFusionTree{<:SectorProduct}) - transposed_indices = outer_multiplicity_split.( - Base.tail(leaves(f)), - branch_sectors(f), - (Base.tail(branch_sectors(f))..., root_sector(f)), - outer_multiplicity_indices(f), - ) - arguments_root = arguments(root_sector(f)) - arguments_leaves = arguments.(leaves(f)) - arguments_branch_sectors = arguments.(branch_sectors(f)) - # TODO way to avoid explicit ntuple? - # works fine for Tuple and NamedTuple SectorProduct - return ntuple( - i -> SectorFusionTree( - getindex.(arguments_leaves, i), - arrows(f), - arguments_root[i], - getindex.(arguments_branch_sectors, i), - getindex.(transposed_indices, i), - ), - length(arguments_root), - ) + transposed_indices = outer_multiplicity_split.( + Base.tail(leaves(f)), + branch_sectors(f), + (Base.tail(branch_sectors(f))..., root_sector(f)), + outer_multiplicity_indices(f), + ) + arguments_root = arguments(root_sector(f)) + arguments_leaves = arguments.(leaves(f)) + arguments_branch_sectors = arguments.(branch_sectors(f)) + # TODO way to avoid explicit ntuple? + # works fine for Tuple and NamedTuple SectorProduct + return ntuple( + i -> SectorFusionTree( + getindex.(arguments_leaves, i), + arrows(f), + arguments_root[i], + getindex.(arguments_branch_sectors, i), + getindex.(transposed_indices, i), + ), + length(arguments_root), + ) end -function GradedArrays.arguments(f::SectorFusionTree{<:SectorProduct,0}) - return map(arg -> SectorFusionTree((), (), arg, (), ()), arguments(root_sector(f))) +function GradedArrays.arguments(f::SectorFusionTree{<:SectorProduct, 0}) + return map(arg -> SectorFusionTree((), (), arg, (), ()), arguments(root_sector(f))) end -function GradedArrays.arguments(f::SectorFusionTree{<:SectorProduct,1}) - arguments_root = arguments(root_sector(f)) - arguments_leave = arguments(only(leaves(f))) - # use map(keys) to stay agnostic with respect to SectorProduct implementation - return map(keys(arguments_root)) do k - return SectorFusionTree((arguments_leave[k],), arrows(f), arguments_root[k], (), ()) - end +function GradedArrays.arguments(f::SectorFusionTree{<:SectorProduct, 1}) + arguments_root = arguments(root_sector(f)) + arguments_leave = arguments(only(leaves(f))) + # use map(keys) to stay agnostic with respect to SectorProduct implementation + return map(keys(arguments_root)) do k + return SectorFusionTree((arguments_leave[k],), arrows(f), arguments_root[k], (), ()) + end end # TBD change type depending on AbelianStyle? @@ -167,19 +167,19 @@ fusiontree_eltype(::Type{<:AbstractSector}) = Float64 # constructors function build_trees(legs::Vararg{AbstractGradedUnitRange}) - # construct all authorized trees for each outer block in legs - tree_arrows = isdual.(legs) - return mapreduce(vcat, Iterators.product(sectors.(flip_dual.(legs))...)) do it - return build_trees(it, tree_arrows) - end + # construct all authorized trees for each outer block in legs + tree_arrows = isdual.(legs) + return mapreduce(vcat, Iterators.product(sectors.(flip_dual.(legs))...)) do it + return build_trees(it, tree_arrows) + end end function build_trees( - sectors_to_fuse::NTuple{N,<:AbstractSector}, arrows_to_fuse::NTuple{N,Bool} -) where {N} - # construct all authorized trees with fixed outer sectors - trees = [SectorFusionTree(first(sectors_to_fuse), first(arrows_to_fuse))] - return recursive_build_trees(trees, Base.tail(sectors_to_fuse), Base.tail(arrows_to_fuse)) + sectors_to_fuse::NTuple{N, <:AbstractSector}, arrows_to_fuse::NTuple{N, Bool} + ) where {N} + # construct all authorized trees with fixed outer sectors + trees = [SectorFusionTree(first(sectors_to_fuse), first(arrows_to_fuse))] + return recursive_build_trees(trees, Base.tail(sectors_to_fuse), Base.tail(arrows_to_fuse)) end # @@ -188,155 +188,155 @@ end # --------------- SectorProduct helper functions --------------- function outer_multiplicity_kron( - sec1, sec2, fused, outer_multiplicity1, outer_multiplicity2 -) - n = nsymbol(sec1, sec2, fused) - linear_inds = LinearIndices((n, outer_multiplicity2)) - return linear_inds[outer_multiplicity1, outer_multiplicity2] + sec1, sec2, fused, outer_multiplicity1, outer_multiplicity2 + ) + n = nsymbol(sec1, sec2, fused) + linear_inds = LinearIndices((n, outer_multiplicity2)) + return linear_inds[outer_multiplicity1, outer_multiplicity2] end function outer_multiplicity_split( - sec1::S, sec2::S, fused::S, outer_mult_index::Integer -) where {S<:SectorProduct} - args1 = arguments(sec1) - args2 = arguments(sec2) - args12 = arguments(fused) - nsymbols = Tuple(map(nsymbol, args1, args2, args12)) # CartesianIndices requires explicit Tuple - return Tuple(CartesianIndices(nsymbols)[outer_mult_index]) + sec1::S, sec2::S, fused::S, outer_mult_index::Integer + ) where {S <: SectorProduct} + args1 = arguments(sec1) + args2 = arguments(sec2) + args12 = arguments(fused) + nsymbols = Tuple(map(nsymbol, args1, args2, args12)) # CartesianIndices requires explicit Tuple + return Tuple(CartesianIndices(nsymbols)[outer_mult_index]) end # --------------- Build trees --------------- # zero leg: need S to get sector type information -function SectorFusionTree{S}() where {S<:AbstractSector} - return SectorFusionTree((), (), trivial(S), (), ()) +function SectorFusionTree{S}() where {S <: AbstractSector} + return SectorFusionTree((), (), trivial(S), (), ()) end -function SectorFusionTree{S}(::Tuple{}, ::Tuple{}) where {S<:AbstractSector} - return SectorFusionTree((), (), trivial(S), (), ()) +function SectorFusionTree{S}(::Tuple{}, ::Tuple{}) where {S <: AbstractSector} + return SectorFusionTree((), (), trivial(S), (), ()) end # one leg function SectorFusionTree(sect::AbstractSector, arrow::Bool) - return SectorFusionTree((sect,), (arrow,), sect, (), ()) + return SectorFusionTree((sect,), (arrow,), sect, (), ()) end -function braid_tuples(t1::Tuple{Vararg{Any,N}}, t2::Tuple{Vararg{Any,N}}) where {N} - t12 = (t1, t2) - nested = ntuple(i -> getindex.(t12, i), N) - return flatten_tuples(nested) +function braid_tuples(t1::Tuple{Vararg{Any, N}}, t2::Tuple{Vararg{Any, N}}) where {N} + t12 = (t1, t2) + nested = ntuple(i -> getindex.(t12, i), N) + return flatten_tuples(nested) end function append_tree_leave( - parent_tree::SectorFusionTree, - branch_sector::AbstractSector, - level_arrow::Bool, - child_root_sector, - outer_mult, -) - child_leaves = (leaves(parent_tree)..., branch_sector) - child_arrows = (arrows(parent_tree)..., level_arrow) - child_branch_sectors = (branch_sectors(parent_tree)..., root_sector(parent_tree)) - child_outer_mul = (outer_multiplicity_indices(parent_tree)..., outer_mult) - return SectorFusionTree( - child_leaves, child_arrows, child_root_sector, child_branch_sectors, child_outer_mul - ) + parent_tree::SectorFusionTree, + branch_sector::AbstractSector, + level_arrow::Bool, + child_root_sector, + outer_mult, + ) + child_leaves = (leaves(parent_tree)..., branch_sector) + child_arrows = (arrows(parent_tree)..., level_arrow) + child_branch_sectors = (branch_sectors(parent_tree)..., root_sector(parent_tree)) + child_outer_mul = (outer_multiplicity_indices(parent_tree)..., outer_mult) + return SectorFusionTree( + child_leaves, child_arrows, child_root_sector, child_branch_sectors, child_outer_mul + ) end function fuse_next_sector( - parent_tree::SectorFusionTree, branch_sector::AbstractSector, level_arrow::Bool -) - new_space = to_gradedrange(root_sector(parent_tree) ⊗ branch_sector) - return mapreduce( - vcat, zip(sectors(new_space), sector_multiplicities(new_space)) - ) do (la, n) - return [ - append_tree_leave(parent_tree, branch_sector, level_arrow, la, outer_mult) for - outer_mult in 1:n - ] - end + parent_tree::SectorFusionTree, branch_sector::AbstractSector, level_arrow::Bool + ) + new_space = to_gradedrange(root_sector(parent_tree) ⊗ branch_sector) + return mapreduce( + vcat, zip(sectors(new_space), sector_multiplicities(new_space)) + ) do (la, n) + return [ + append_tree_leave(parent_tree, branch_sector, level_arrow, la, outer_mult) for + outer_mult in 1:n + ] + end end function recursive_build_trees( - old_trees::Vector, sectors_to_fuse::Tuple, arrows_to_fuse::Tuple -) - next_level_trees = mapreduce(vcat, old_trees) do tree - return fuse_next_sector(tree, first(sectors_to_fuse), first(arrows_to_fuse)) - end - return recursive_build_trees( - next_level_trees, Base.tail(sectors_to_fuse), Base.tail(arrows_to_fuse) - ) + old_trees::Vector, sectors_to_fuse::Tuple, arrows_to_fuse::Tuple + ) + next_level_trees = mapreduce(vcat, old_trees) do tree + return fuse_next_sector(tree, first(sectors_to_fuse), first(arrows_to_fuse)) + end + return recursive_build_trees( + next_level_trees, Base.tail(sectors_to_fuse), Base.tail(arrows_to_fuse) + ) end function recursive_build_trees(trees::Vector, ::Tuple{}, ::Tuple{}) - return trees + return trees end # --------------- convert to Array --------------- -to_array(::SectorFusionTree{<:Any,0}) = ones(1) +to_array(::SectorFusionTree{<:Any, 0}) = ones(1) function to_array(f::SectorFusionTree) - # init with dummy trivial leg to get arrow correct and deal with size-1 case - cgt1 = clebsch_gordan_tensor( - trivial(sector_type(f)), first(leaves(f)), first(leaves(f)), false, first(arrows(f)), 1 - ) - tree_tensor = cgt1[1, :, :] - return grow_tensor_tree(tree_tensor, f) + # init with dummy trivial leg to get arrow correct and deal with size-1 case + cgt1 = clebsch_gordan_tensor( + trivial(sector_type(f)), first(leaves(f)), first(leaves(f)), false, first(arrows(f)), 1 + ) + tree_tensor = cgt1[1, :, :] + return grow_tensor_tree(tree_tensor, f) end function to_array(f::SectorFusionTree{<:SectorProduct}) - args = convert.(Array, arguments(f)) - return reduce(_tensor_kron, args) + args = convert.(Array, arguments(f)) + return reduce(_tensor_kron, args) end # LinearAlgebra.kron does not allow input for ndims>2 -function _tensor_kron(a::AbstractArray{<:Any,N}, b::AbstractArray{<:Any,N}) where {N} - t1 = ntuple(_ -> 1, N) - sha = braid_tuples(size(a), t1) - shb = braid_tuples(t1, size(b)) - c = reshape(a, sha) .* reshape(b, shb) - return reshape(c, size(a) .* size(b)) +function _tensor_kron(a::AbstractArray{<:Any, N}, b::AbstractArray{<:Any, N}) where {N} + t1 = ntuple(_ -> 1, N) + sha = braid_tuples(size(a), t1) + shb = braid_tuples(t1, size(b)) + c = reshape(a, sha) .* reshape(b, shb) + return reshape(c, size(a) .* size(b)) end function contract_clebsch_gordan(tree_tensor::AbstractArray, cgt::AbstractArray) - N = ndims(tree_tensor) - return contract( - (ntuple(identity, N - 1)..., N + 1, N + 2), - tree_tensor, - ntuple(identity, N), - cgt, - (N, N + 1, N + 2), - ) + N = ndims(tree_tensor) + return contract( + (ntuple(identity, N - 1)..., N + 1, N + 2), + tree_tensor, + ntuple(identity, N), + cgt, + (N, N + 1, N + 2), + ) end # specialized code when branch_sector is empty -function grow_tensor_tree(tree_tensor::AbstractArray{<:Real,2}, ::SectorFusionTree{<:Any,1}) - return tree_tensor +function grow_tensor_tree(tree_tensor::AbstractArray{<:Real, 2}, ::SectorFusionTree{<:Any, 1}) + return tree_tensor end function grow_tensor_tree( - tree_tensor::AbstractArray{<:Real,N}, f::SectorFusionTree -) where {N} - cgt = clebsch_gordan_tensor( - branch_sectors(f)[N - 1], - leaves(f)[N], - branch_sectors(f)[N], - false, - arrows(f)[N], - outer_multiplicity_indices(f)[N - 1], - ) - next_level_tree = contract_clebsch_gordan(tree_tensor, cgt) - return grow_tensor_tree(next_level_tree, f) + tree_tensor::AbstractArray{<:Real, N}, f::SectorFusionTree + ) where {N} + cgt = clebsch_gordan_tensor( + branch_sectors(f)[N - 1], + leaves(f)[N], + branch_sectors(f)[N], + false, + arrows(f)[N], + outer_multiplicity_indices(f)[N - 1], + ) + next_level_tree = contract_clebsch_gordan(tree_tensor, cgt) + return grow_tensor_tree(next_level_tree, f) end function grow_tensor_tree( - tree_tensor::AbstractArray{<:Real,N}, f::SectorFusionTree{<:Any,N} -) where {N} - cgt = clebsch_gordan_tensor( - last(branch_sectors(f)), - last(leaves(f)), - root_sector(f), - false, - last(arrows(f)), - last(outer_multiplicity_indices(f)), - ) - return contract_clebsch_gordan(tree_tensor, cgt) + tree_tensor::AbstractArray{<:Real, N}, f::SectorFusionTree{<:Any, N} + ) where {N} + cgt = clebsch_gordan_tensor( + last(branch_sectors(f)), + last(leaves(f)), + root_sector(f), + false, + last(arrows(f)), + last(outer_multiplicity_indices(f)), + ) + return contract_clebsch_gordan(tree_tensor, cgt) end diff --git a/src/fusiontensor/array_cast.jl b/src/fusiontensor/array_cast.jl index e505f19..201838e 100644 --- a/src/fusiontensor/array_cast.jl +++ b/src/fusiontensor/array_cast.jl @@ -12,20 +12,20 @@ using TensorAlgebra: contract #### cast from symmetric to array function BlockSparseArrays.BlockSparseArray(ft::FusionTensor) - return to_array(ft) + return to_array(ft) end # ================================= Low level interface ================================== function to_fusiontensor( - array::AbstractArray, - codomain_legs::Tuple{Vararg{AbstractGradedUnitRange}}, - domain_legs::Tuple{Vararg{AbstractGradedUnitRange}}; - atol::Real=0, - rtol::Real=rtoldefault(array), -) - bounds = blocklengths.((codomain_legs..., domain_legs...)) - blockarray = BlockedArray(array, bounds...) - return to_fusiontensor(blockarray, codomain_legs, domain_legs; atol, rtol) + array::AbstractArray, + codomain_legs::Tuple{Vararg{AbstractGradedUnitRange}}, + domain_legs::Tuple{Vararg{AbstractGradedUnitRange}}; + atol::Real = 0, + rtol::Real = rtoldefault(array), + ) + bounds = blocklengths.((codomain_legs..., domain_legs...)) + blockarray = BlockedArray(array, bounds...) + return to_fusiontensor(blockarray, codomain_legs, domain_legs; atol, rtol) end rtoldefault(a::AbstractArray) = rtoldefault(eltype(a)) @@ -33,53 +33,53 @@ rtoldefault(arrayt::Type{<:AbstractArray}) = rtoldefault(eltype(arrayt)) rtoldefault(elt::Type{<:Number}) = 10 * eps(real(float(elt))) function checknorm(ft::FusionTensor, a::AbstractArray, atol::Real, rtol::Real) - return isapprox(norm(ft), norm(a); atol, rtol) || - throw(InexactError(:FusionTensor, typeof(ft), a)) + return isapprox(norm(ft), norm(a); atol, rtol) || + throw(InexactError(:FusionTensor, typeof(ft), a)) end function to_fusiontensor( - blockarray::AbstractBlockArray, - codomain_legs::Tuple{Vararg{AbstractGradedUnitRange}}, - domain_legs::Tuple{Vararg{AbstractGradedUnitRange}}; - atol::Real=0, - rtol::Real=rtoldefault(blockarray), -) - # input validation - if length(codomain_legs) + length(domain_legs) != ndims(blockarray) # compile time - throw(DomainError("legs are incompatible with array ndims")) - end - if quantum_dimension.((codomain_legs..., domain_legs...)) != size(blockarray) - throw(DomainError("legs dimensions are incompatible with array")) - end - - ft = to_fusiontensor_no_checknorm(blockarray, codomain_legs, domain_legs) - - # if blockarray is not G-invariant, norm(ft) < norm(blockarray) - checknorm(ft, blockarray, atol, rtol) - return ft + blockarray::AbstractBlockArray, + codomain_legs::Tuple{Vararg{AbstractGradedUnitRange}}, + domain_legs::Tuple{Vararg{AbstractGradedUnitRange}}; + atol::Real = 0, + rtol::Real = rtoldefault(blockarray), + ) + # input validation + if length(codomain_legs) + length(domain_legs) != ndims(blockarray) # compile time + throw(DomainError("legs are incompatible with array ndims")) + end + if quantum_dimension.((codomain_legs..., domain_legs...)) != size(blockarray) + throw(DomainError("legs dimensions are incompatible with array")) + end + + ft = to_fusiontensor_no_checknorm(blockarray, codomain_legs, domain_legs) + + # if blockarray is not G-invariant, norm(ft) < norm(blockarray) + checknorm(ft, blockarray, atol, rtol) + return ft end function to_fusiontensor_no_checknorm( - blockarray::AbstractBlockArray, - codomain_legs::Tuple{Vararg{AbstractGradedUnitRange}}, - domain_legs::Tuple{Vararg{AbstractGradedUnitRange}}, -) - ft = FusionTensor{eltype(blockarray)}(undef, codomain_legs, domain_legs) - for (f1, f2) in keys(trees_block_mapping(ft)) - b = findblock(ft, f1, f2) - ft[f1, f2] = contract_fusion_trees(blockarray[b], f1, f2) - end - return ft + blockarray::AbstractBlockArray, + codomain_legs::Tuple{Vararg{AbstractGradedUnitRange}}, + domain_legs::Tuple{Vararg{AbstractGradedUnitRange}}, + ) + ft = FusionTensor{eltype(blockarray)}(undef, codomain_legs, domain_legs) + for (f1, f2) in keys(trees_block_mapping(ft)) + b = findblock(ft, f1, f2) + ft[f1, f2] = contract_fusion_trees(blockarray[b], f1, f2) + end + return ft end function to_array(ft::FusionTensor) - bounds = blocklengths.((codomain_axes(ft)..., domain_axes(ft)...)) - bsa = BlockSparseArray{eltype(ft)}(undef, blockedrange.(bounds)) - for (f1, f2) in keys(trees_block_mapping(ft)) - b = findblock(ft, f1, f2) - bsa[b] = bsa[b] + contract_fusion_trees(ft, f1, f2) # init block when needed - end - return bsa + bounds = blocklengths.((codomain_axes(ft)..., domain_axes(ft)...)) + bsa = BlockSparseArray{eltype(ft)}(undef, blockedrange.(bounds)) + for (f1, f2) in keys(trees_block_mapping(ft)) + b = findblock(ft, f1, f2) + bsa[b] = bsa[b] + contract_fusion_trees(ft, f1, f2) # init block when needed + end + return bsa end # ===================================== Internals ======================================== @@ -87,107 +87,107 @@ end #----------------------------------- misc tools ---------------------------------------- function degen_dims_shape(dense_shape::Tuple{Vararg{Int}}, f::SectorFusionTree) - dims = quantum_dimension.(leaves(f)) - mults = dense_shape .÷ dims - return braid_tuples(mults, dims) + dims = quantum_dimension.(leaves(f)) + mults = dense_shape .÷ dims + return braid_tuples(mults, dims) end function split_degen_dims( - array_block::AbstractArray, f1::SectorFusionTree, f2::SectorFusionTree -) - array_block_split_shape = ( - degen_dims_shape(size(array_block)[begin:length(f1)], f1)..., - degen_dims_shape(size(array_block)[(length(f1) + 1):end], f2)..., - ) - return reshape(array_block, array_block_split_shape) + array_block::AbstractArray, f1::SectorFusionTree, f2::SectorFusionTree + ) + array_block_split_shape = ( + degen_dims_shape(size(array_block)[begin:length(f1)], f1)..., + degen_dims_shape(size(array_block)[(length(f1) + 1):end], f2)..., + ) + return reshape(array_block, array_block_split_shape) end function merge_degen_dims(split_array_block::AbstractArray) - s0 = size(split_array_block) - array_shape = - ntuple(i -> s0[2 * i - 1], length(s0) ÷ 2) .* ntuple(i -> s0[2 * i], length(s0) ÷ 2) - array_block = reshape(split_array_block, array_shape) - return array_block + s0 = size(split_array_block) + array_shape = + ntuple(i -> s0[2 * i - 1], length(s0) ÷ 2) .* ntuple(i -> s0[2 * i], length(s0) ÷ 2) + array_block = reshape(split_array_block, array_shape) + return array_block end #---------------------------------- cast from array --------------------------------------- function contract_fusion_trees( - array_block::AbstractArray, f1::SectorFusionTree, f2::SectorFusionTree -) - # start from an array outer block with e.g. N=6 axes divided into N_DO=3 ndims_codomain - # and N_CO=3 ndims_domain. Each leg k can be decomposed as a product of external an - # multiplicity extk and a quantum dimension dimk - # - # ------------------------------array_block------------------------------- - # | | | | | | - # ext1*dim1 ext2*dim2 ext3*dim3 ext4*dim4 ext5*dim5 ext6*dim6 - # - - # each leg of this this array outer block can now be opened to form a 2N-dims tensor. - # note that this 2N-dims form is only defined at the level of the outer block, - # not for a larger block. - # - # ------------------------------split_array_block------------------------- - # | | | | | | - # / \ / \ / \ / \ / \ / \ - # / \ / \ / \ / \ / \ / \ - # ext1 dim1 ext2 dim2 ext3 dim3 ext4 dim4 ext5 dim5 ext6 dim6 - # - N = ndims(array_block) - - split_array_block = split_degen_dims(array_block, f1, f2) - dim_sec = quantum_dimension(root_sector(f1)) - p = contract_singlet_projector(f1, f2) - - # ------------------------------split_array_block------------------------- - # | | | | | | - # / \ / \ / \ / \ / \ / \ - # / \ / \ / \ / \ / \ / \ - # ext1 dim1 ext2 dim2 ext3 dim3 ext4 dim4 ext5 dim5 ext6 dim6 - # \ | | \__ | | - # \________ | | \__ | | - # \| | \ _| | - # \___________| \___________| - # \ | - # \----------------dim_sec---------------- / - res = contract( - ntuple(i -> 2 * i - 1, N), - split_array_block, - ntuple(identity, 2 * N), - p, - ntuple(i -> 2 * i, N), - ) - # normalization factor - res ./= dim_sec - return res + array_block::AbstractArray, f1::SectorFusionTree, f2::SectorFusionTree + ) + # start from an array outer block with e.g. N=6 axes divided into N_DO=3 ndims_codomain + # and N_CO=3 ndims_domain. Each leg k can be decomposed as a product of external an + # multiplicity extk and a quantum dimension dimk + # + # ------------------------------array_block------------------------------- + # | | | | | | + # ext1*dim1 ext2*dim2 ext3*dim3 ext4*dim4 ext5*dim5 ext6*dim6 + # + + # each leg of this this array outer block can now be opened to form a 2N-dims tensor. + # note that this 2N-dims form is only defined at the level of the outer block, + # not for a larger block. + # + # ------------------------------split_array_block------------------------- + # | | | | | | + # / \ / \ / \ / \ / \ / \ + # / \ / \ / \ / \ / \ / \ + # ext1 dim1 ext2 dim2 ext3 dim3 ext4 dim4 ext5 dim5 ext6 dim6 + # + N = ndims(array_block) + + split_array_block = split_degen_dims(array_block, f1, f2) + dim_sec = quantum_dimension(root_sector(f1)) + p = contract_singlet_projector(f1, f2) + + # ------------------------------split_array_block------------------------- + # | | | | | | + # / \ / \ / \ / \ / \ / \ + # / \ / \ / \ / \ / \ / \ + # ext1 dim1 ext2 dim2 ext3 dim3 ext4 dim4 ext5 dim5 ext6 dim6 + # \ | | \__ | | + # \________ | | \__ | | + # \| | \ _| | + # \___________| \___________| + # \ | + # \----------------dim_sec---------------- / + res = contract( + ntuple(i -> 2 * i - 1, N), + split_array_block, + ntuple(identity, 2 * N), + p, + ntuple(i -> 2 * i, N), + ) + # normalization factor + res ./= dim_sec + return res end function contract_singlet_projector(f1::SectorFusionTree, f2::SectorFusionTree) - f1_array = convert(Array, f1) - f2_array = convert(Array, flip(f2)) - N_CO = length(f1) - N_DO = length(f2) - return contract( - ntuple(identity, N_CO + N_DO), - f1_array, - (ntuple(identity, N_CO)..., N_CO + N_DO + 1), - f2_array, - (ntuple(i -> i + N_CO, N_DO)..., N_CO + N_DO + 1), - ) + f1_array = convert(Array, f1) + f2_array = convert(Array, flip(f2)) + N_CO = length(f1) + N_DO = length(f2) + return contract( + ntuple(identity, N_CO + N_DO), + f1_array, + (ntuple(identity, N_CO)..., N_CO + N_DO + 1), + f2_array, + (ntuple(i -> i + N_CO, N_DO)..., N_CO + N_DO + 1), + ) end #----------------------------------- cast to array ---------------------------------------- function contract_fusion_trees(ft::FusionTensor, f1::SectorFusionTree, f2::SectorFusionTree) - N = ndims(ft) - charge_block = view(ft, f1, f2) - p = contract_singlet_projector(f1, f2) + N = ndims(ft) + charge_block = view(ft, f1, f2) + p = contract_singlet_projector(f1, f2) - # outer product of charge block and fusion tree - labels1 = ntuple(identity, N) - labels2 = ntuple(i -> i + N, N) - perm = braid_tuples(ntuple(identity, N), ntuple(i -> i + N, N)) - split_array_block = contract(perm, charge_block, labels1, p, labels2) + # outer product of charge block and fusion tree + labels1 = ntuple(identity, N) + labels2 = ntuple(i -> i + N, N) + perm = braid_tuples(ntuple(identity, N), ntuple(i -> i + N, N)) + split_array_block = contract(perm, charge_block, labels1, p, labels2) - return merge_degen_dims(split_array_block) + return merge_degen_dims(split_array_block) end diff --git a/src/fusiontensor/base_interface.jl b/src/fusiontensor/base_interface.jl index 36b7c89..7615daf 100644 --- a/src/fusiontensor/base_interface.jl +++ b/src/fusiontensor/base_interface.jl @@ -12,22 +12,22 @@ Base.:*(ft::FusionTensor, x::Number) = x * ft # tensor contraction is a block data_matrix product. function Base.:*(left::FusionTensor, right::FusionTensor) - checkspaces_dual(domain_axes(left), codomain_axes(right)) - new_data_matrix = data_matrix(left) * data_matrix(right) - return FusionTensor(new_data_matrix, codomain_axes(left), domain_axes(right)) + checkspaces_dual(domain_axes(left), codomain_axes(right)) + new_data_matrix = data_matrix(left) * data_matrix(right) + return FusionTensor(new_data_matrix, codomain_axes(left), domain_axes(right)) end # tensor addition is a block data_matrix add. function Base.:+(left::FusionTensor, right::FusionTensor) - checkspaces(axes(left), axes(right)) - return set_data_matrix(left, data_matrix(left) + data_matrix(right)) + checkspaces(axes(left), axes(right)) + return set_data_matrix(left, data_matrix(left) + data_matrix(right)) end Base.:-(ft::FusionTensor) = set_data_matrix(ft, -data_matrix(ft)) function Base.:-(left::FusionTensor, right::FusionTensor) - checkspaces(axes(left), axes(right)) - return set_data_matrix(left, data_matrix(left) - data_matrix(right)) + checkspaces(axes(left), axes(right)) + return set_data_matrix(left, data_matrix(left) - data_matrix(right)) end Base.:/(ft::FusionTensor, x::Number) = set_data_matrix(ft, data_matrix(ft) / x) @@ -37,17 +37,17 @@ Base.Array(ft::FusionTensor) = Array(to_array(ft)) # adjoint is costless: dual axes, swap codomain and domain, take data_matrix adjoint. # data_matrix coeff are not modified (beyond complex conjugation) function transpose_mapping(d::Dict) - return Dict([reverse(flip.(k)) => transpose_mapping(v) for (k, v) in d]) + return Dict([reverse(flip.(k)) => transpose_mapping(v) for (k, v) in d]) end function transpose_mapping(b::BlockIndexRange{2}) - new_block = Block(reverse(Tuple(Block(b)))) - return new_block[reverse(b.indices)...] + new_block = Block(reverse(Tuple(Block(b)))) + return new_block[reverse(b.indices)...] end function Base.adjoint(ft::FusionTensor) - new_axes = FusionTensorAxes(dual.(domain_axes(ft)), dual.(codomain_axes(ft))) - return FusionTensor( - adjoint(data_matrix(ft)), new_axes, transpose_mapping(trees_block_mapping(ft)) - ) + new_axes = FusionTensorAxes(dual.(domain_axes(ft)), dual.(codomain_axes(ft))) + return FusionTensor( + adjoint(data_matrix(ft)), new_axes, transpose_mapping(trees_block_mapping(ft)) + ) end Base.axes(ft::FusionTensor) = ft.axes @@ -57,20 +57,20 @@ Base.conj(ft::FusionTensor{<:Real}) = ft # same object for real element type Base.conj(ft::FusionTensor) = set_data_matrix(ft, conj(data_matrix(ft))) function Base.copy(ft::FusionTensor) - return FusionTensor(copy(data_matrix(ft)), copy(axes(ft)), copy(trees_block_mapping(ft))) + return FusionTensor(copy(data_matrix(ft)), copy(axes(ft)), copy(trees_block_mapping(ft))) end function Base.deepcopy(ft::FusionTensor) - return FusionTensor( - deepcopy(data_matrix(ft)), deepcopy(axes(ft)), deepcopy(trees_block_mapping(ft)) - ) + return FusionTensor( + deepcopy(data_matrix(ft)), deepcopy(axes(ft)), deepcopy(trees_block_mapping(ft)) + ) end Base.eachindex(::FusionTensor) = throw(MethodError(eachindex, (FusionTensor,))) function Base.getindex(ft::FusionTensor, f1::SectorFusionTree, f2::SectorFusionTree) - charge_matrix = data_matrix(ft)[trees_block_mapping(ft)[f1, f2]] - return reshape(charge_matrix, charge_block_size(ft, f1, f2)) + charge_matrix = data_matrix(ft)[trees_block_mapping(ft)[f1, f2]] + return reshape(charge_matrix, charge_block_size(ft, f1, f2)) end Base.imag(ft::FusionTensor) = set_data_matrix(ft, imag(data_matrix(ft))) @@ -78,66 +78,66 @@ Base.imag(ft::FusionTensor) = set_data_matrix(ft, imag(data_matrix(ft))) Base.permutedims(ft::FusionTensor, args...) = fusiontensor_permutedims(ft, args...) function Base.permutedims!(ftdst::FusionTensor, ftsrc::FusionTensor, args...) - return fusiontensor_permutedims!(ftdst, ftsrc, args...) + return fusiontensor_permutedims!(ftdst, ftsrc, args...) end Base.real(ft::FusionTensor{<:Real}) = ft # same object Base.real(ft::FusionTensor) = set_data_matrix(ft, real(data_matrix(ft))) function Base.setindex!( - ft::FusionTensor, a::AbstractArray, f1::SectorFusionTree, f2::SectorFusionTree -) - return view(ft, f1, f2) .= a + ft::FusionTensor, a::AbstractArray, f1::SectorFusionTree, f2::SectorFusionTree + ) + return view(ft, f1, f2) .= a end function Base.similar(ft::FusionTensor, T::Type) - # reuse trees_block_mapping + # reuse trees_block_mapping - # some fusion trees have Float64 eltype: need compatible type - @assert promote_type(T, fusiontree_eltype(sector_type(ft))) === T - mat = initialize_data_matrix(T, codomain_axis(ft), domain_axis(ft)) - return FusionTensor(mat, axes(ft), trees_block_mapping(ft)) + # some fusion trees have Float64 eltype: need compatible type + @assert promote_type(T, fusiontree_eltype(sector_type(ft))) === T + mat = initialize_data_matrix(T, codomain_axis(ft), domain_axis(ft)) + return FusionTensor(mat, axes(ft), trees_block_mapping(ft)) end function Base.similar( - ::FusionTensor, ::Type, t::Tuple{AbstractGradedUnitRange,Vararg{AbstractGradedUnitRange}} -) - throw(MethodError(similar, (typeof(t),))) + ::FusionTensor, ::Type, t::Tuple{AbstractGradedUnitRange, Vararg{AbstractGradedUnitRange}} + ) + throw(MethodError(similar, (typeof(t),))) end function Base.similar(::FusionTensor, ::Type, ::Tuple{}) - throw(MethodError(similar, (Tuple{},))) + throw(MethodError(similar, (Tuple{},))) end function Base.similar( - ft::FusionTensor, ::Type{T}, new_axes::Tuple{<:Tuple,<:Tuple} -) where {T} - return similar(ft, T, tuplemortar(new_axes)) + ft::FusionTensor, ::Type{T}, new_axes::Tuple{<:Tuple, <:Tuple} + ) where {T} + return similar(ft, T, tuplemortar(new_axes)) end function Base.similar(ft::FusionTensor, ::Type{T}, new_axes::BlockedTuple{2}) where {T} - return similar(ft, T, FusionTensorAxes(new_axes)) + return similar(ft, T, FusionTensorAxes(new_axes)) end function Base.similar(ft::FusionTensor, new_axes::FusionTensorAxes) - return similar(ft, eltype(ft), new_axes) + return similar(ft, eltype(ft), new_axes) end function Base.similar(::FusionTensor, ::Type{T}, new_axes::FusionTensorAxes) where {T} - return FusionTensor{T}(undef, new_axes) + return FusionTensor{T}(undef, new_axes) end function Base.show(io::IO, ft::FusionTensor) - return print(io, "$(ndims(ft))-dim FusionTensor with size $(size(ft))") + return print(io, "$(ndims(ft))-dim FusionTensor with size $(size(ft))") end function Base.show(io::IO, ::MIME"text/plain", ft::FusionTensor) - print(io, ft) - print(" and axes:") - for ax in axes(ft) - print(io, "\n", ax) - end - return nothing + print(io, ft) + print(" and axes:") + for ax in axes(ft) + print(io, "\n", ax) + end + return nothing end Base.size(ft::FusionTensor) = quantum_dimension.(axes(ft)) function Base.view(ft::FusionTensor, f1::SectorFusionTree, f2::SectorFusionTree) - charge_matrix = @view! data_matrix(ft)[trees_block_mapping(ft)[f1, f2]] - return reshape(charge_matrix, charge_block_size(ft, f1, f2)) + charge_matrix = @view! data_matrix(ft)[trees_block_mapping(ft)[f1, f2]] + return reshape(charge_matrix, charge_block_size(ft, f1, f2)) end diff --git a/src/fusiontensor/fusiontensor.jl b/src/fusiontensor/fusiontensor.jl index 6f9d93d..fcc232e 100644 --- a/src/fusiontensor/fusiontensor.jl +++ b/src/fusiontensor/fusiontensor.jl @@ -3,23 +3,23 @@ using BlockArrays: AbstractBlockMatrix, BlockArrays, BlockIndexRange, blocklength, findblock using BlockSparseArrays: - AbstractBlockSparseMatrix, BlockSparseArray, eachblockstoredindex, to_block_indices + AbstractBlockSparseMatrix, BlockSparseArray, eachblockstoredindex, to_block_indices using GradedArrays: - AbstractGradedUnitRange, - SectorProduct, - TrivialSector, - dual, - findfirstblock, - flip, - flip_dual, - gradedrange, - isdual, - map_sectors, - sector_multiplicity, - sector_type, - sectormergesort, - sectors, - space_isequal + AbstractGradedUnitRange, + SectorProduct, + TrivialSector, + dual, + findfirstblock, + flip, + flip_dual, + gradedrange, + isdual, + map_sectors, + sector_multiplicity, + sector_type, + sectormergesort, + sectors, + space_isequal using LinearAlgebra: UniformScaling using Random: Random, AbstractRNG, randn! using TensorAlgebra: BlockedTuple, trivial_axis, tuplemortar, length_codomain, length_domain @@ -28,115 +28,115 @@ using TypeParameterAccessors: type_parameters # ======================================= Misc =========================================== function flip_domain(nonflipped_col_axis, nonflipped_trees_to_ranges) - col_axis = dual(nonflipped_col_axis) - domain_trees_to_ranges_mapping = Dict( - map(((tree, v),) -> flip(tree) => v, collect(nonflipped_trees_to_ranges)) - ) - return col_axis, domain_trees_to_ranges_mapping + col_axis = dual(nonflipped_col_axis) + domain_trees_to_ranges_mapping = Dict( + map(((tree, v),) -> flip(tree) => v, collect(nonflipped_trees_to_ranges)) + ) + return col_axis, domain_trees_to_ranges_mapping end -function fuse_axes(::Type{S}, ::Tuple{}) where {S<:AbstractSector} - fused_axis = trivial_axis(S) - trees_to_ranges_mapping = Dict([SectorFusionTree{S}() => Block(1)[1:1]]) - return fused_axis, trees_to_ranges_mapping +function fuse_axes(::Type{S}, ::Tuple{}) where {S <: AbstractSector} + fused_axis = trivial_axis(S) + trees_to_ranges_mapping = Dict([SectorFusionTree{S}() => Block(1)[1:1]]) + return fused_axis, trees_to_ranges_mapping end function fuse_axes(::Type, outer_legs::Tuple{Vararg{AbstractGradedUnitRange}}) - fusion_trees_mult = fusion_trees_external_multiplicities(outer_legs) - fused_leg, trees_to_ranges_mapping = compute_inner_ranges(fusion_trees_mult) - return fused_leg, trees_to_ranges_mapping + fusion_trees_mult = fusion_trees_external_multiplicities(outer_legs) + fused_leg, trees_to_ranges_mapping = compute_inner_ranges(fusion_trees_mult) + return fused_leg, trees_to_ranges_mapping end function fusion_trees_external_multiplicities( - outer_legs::Tuple{Vararg{AbstractGradedUnitRange}} -) - return Iterators.flatten( - block_fusion_trees_external_multiplicities.(Iterators.product(blocks.(outer_legs)...)) - ) + outer_legs::Tuple{Vararg{AbstractGradedUnitRange}} + ) + return Iterators.flatten( + block_fusion_trees_external_multiplicities.(Iterators.product(blocks.(outer_legs)...)) + ) end function block_fusion_trees_external_multiplicities(it::Tuple{Vararg{AbstractUnitRange}}) - block_sectors = only.(sectors.(flip_dual.(it))) - block_mult = prod(sector_multiplicity.(it)) - return build_trees(block_sectors, isdual.(it)) .=> block_mult + block_sectors = only.(sectors.(flip_dual.(it))) + block_mult = prod(sector_multiplicity.(it)) + return build_trees(block_sectors, isdual.(it)) .=> block_mult end function compute_inner_ranges(fusion_trees_mult) - fused_leg = sectormergesort( - gradedrange(root_sector.(first.(fusion_trees_mult)) .=> last.(fusion_trees_mult)) - ) - range_mapping = Dict{type_parameters(eltype(fusion_trees_mult), 1),typeof(Block(1)[1:1])}() - fused_sectors = sectors(fused_leg) - shifts = ones(Int, blocklength(fused_leg)) - for (f, m) in fusion_trees_mult - s = root_sector(f) - i = findfirst(==(s), fused_sectors) - range_mapping[f] = Block(i)[shifts[i]:(shifts[i] + m - 1)] - shifts[i] += m - end - return fused_leg, range_mapping + fused_leg = sectormergesort( + gradedrange(root_sector.(first.(fusion_trees_mult)) .=> last.(fusion_trees_mult)) + ) + range_mapping = Dict{type_parameters(eltype(fusion_trees_mult), 1), typeof(Block(1)[1:1])}() + fused_sectors = sectors(fused_leg) + shifts = ones(Int, blocklength(fused_leg)) + for (f, m) in fusion_trees_mult + s = root_sector(f) + i = findfirst(==(s), fused_sectors) + range_mapping[f] = Block(i)[shifts[i]:(shifts[i] + m - 1)] + shifts[i] += m + end + return fused_leg, range_mapping end function intersect_codomain_domain( - codomain_trees_to_ranges_mapping::Dict{<:SectorFusionTree,<:BlockIndexRange{1}}, - domain_trees_to_ranges_mapping::Dict{<:SectorFusionTree,<:BlockIndexRange{1}}, -) - return Dict( - map( - Iterators.filter( - t -> root_sector(first(first(t))) == dual(root_sector(first(t[2]))), - Iterators.product(codomain_trees_to_ranges_mapping, domain_trees_to_ranges_mapping), - ), - ) do t - return first.(t) => BlockIndexRange(last.(t)) - end, - ) + codomain_trees_to_ranges_mapping::Dict{<:SectorFusionTree, <:BlockIndexRange{1}}, + domain_trees_to_ranges_mapping::Dict{<:SectorFusionTree, <:BlockIndexRange{1}}, + ) + return Dict( + map( + Iterators.filter( + t -> root_sector(first(first(t))) == dual(root_sector(first(t[2]))), + Iterators.product(codomain_trees_to_ranges_mapping, domain_trees_to_ranges_mapping), + ), + ) do t + return first.(t) => BlockIndexRange(last.(t)) + end, + ) end function initialize_data_matrix( - elt::Type{<:Number}, - codomain_axis::AbstractGradedUnitRange, - domain_axis::AbstractGradedUnitRange, -) - @assert sector_type(codomain_axis) == sector_type(domain_axis) - # non-abelian fusion trees have float eltype: need compatible type - promoted = promote_type(elt, fusiontree_eltype(sector_type(domain_axis))) - mat = BlockSparseArray{promoted}( - undef, - blockedrange(sector_multiplicities(codomain_axis)), - blockedrange(sector_multiplicities(domain_axis)), - ) - row_sectors = sectors(codomain_axis) - col_sectors = sectors(domain_axis) - row_block_indices = findall(in(col_sectors), row_sectors) - col_block_indices = findall(in(row_sectors), col_sectors) - for rc in zip(row_block_indices, col_block_indices) - mat[Block(rc)] = mat[Block(rc)] - end - return mat + elt::Type{<:Number}, + codomain_axis::AbstractGradedUnitRange, + domain_axis::AbstractGradedUnitRange, + ) + @assert sector_type(codomain_axis) == sector_type(domain_axis) + # non-abelian fusion trees have float eltype: need compatible type + promoted = promote_type(elt, fusiontree_eltype(sector_type(domain_axis))) + mat = BlockSparseArray{promoted}( + undef, + blockedrange(sector_multiplicities(codomain_axis)), + blockedrange(sector_multiplicities(domain_axis)), + ) + row_sectors = sectors(codomain_axis) + col_sectors = sectors(domain_axis) + row_block_indices = findall(in(col_sectors), row_sectors) + col_block_indices = findall(in(row_sectors), col_sectors) + for rc in zip(row_block_indices, col_block_indices) + mat[Block(rc)] = mat[Block(rc)] + end + return mat end # ==================================== Definitions ======================================= -struct FusionTensor{T,N,Axes<:FusionTensorAxes,Mat<:AbstractMatrix{T},Mapping} <: - AbstractArray{T,N} - data_matrix::Mat - axes::Axes - trees_block_mapping::Mapping - - # inner constructor to impose constraints on types - function FusionTensor{T,N,Axes,Mat,Mapping}( - mat, legs, trees_block_mapping - ) where {T,N,Axes,Mat,Mapping} - S = sector_type(legs) - @assert keytype(trees_block_mapping) <: - Tuple{<:SectorFusionTree{S},<:SectorFusionTree{S}} - return new{T,N,Axes,Mat,Mapping}(mat, legs, trees_block_mapping) - end +struct FusionTensor{T, N, Axes <: FusionTensorAxes, Mat <: AbstractMatrix{T}, Mapping} <: + AbstractArray{T, N} + data_matrix::Mat + axes::Axes + trees_block_mapping::Mapping + + # inner constructor to impose constraints on types + function FusionTensor{T, N, Axes, Mat, Mapping}( + mat, legs, trees_block_mapping + ) where {T, N, Axes, Mat, Mapping} + S = sector_type(legs) + @assert keytype(trees_block_mapping) <: + Tuple{<:SectorFusionTree{S}, <:SectorFusionTree{S}} + return new{T, N, Axes, Mat, Mapping}(mat, legs, trees_block_mapping) + end end -const FusionMatrix{T,Axes,Mat,Mapping} = FusionTensor{ - T,2,Axes,Mapping -} where {BT<:BlockedTuple{2,(1, 1)},Axes<:FusionTensorAxes{BT}} +const FusionMatrix{T, Axes, Mat, Mapping} = FusionTensor{ + T, 2, Axes, Mapping, +} where {BT <: BlockedTuple{2, (1, 1)}, Axes <: FusionTensorAxes{BT}} # ===================================== Accessors ======================================== @@ -146,114 +146,114 @@ trees_block_mapping(ft::FusionTensor) = ft.trees_block_mapping # ==================================== Constructors ====================================== function FusionTensor( - mat::AbstractMatrix, - legs::FusionTensorAxes, - trees_block_mapping::Dict{<:Tuple{<:SectorFusionTree,<:SectorFusionTree}}, -) - return FusionTensor{ - eltype(mat),length(legs),typeof(legs),typeof(mat),typeof(trees_block_mapping) - }( - mat, legs, trees_block_mapping - ) + mat::AbstractMatrix, + legs::FusionTensorAxes, + trees_block_mapping::Dict{<:Tuple{<:SectorFusionTree, <:SectorFusionTree}}, + ) + return FusionTensor{ + eltype(mat), length(legs), typeof(legs), typeof(mat), typeof(trees_block_mapping), + }( + mat, legs, trees_block_mapping + ) end # empty matrix function FusionTensor{T}(::UndefInitializer, legs::FusionTensorAxes) where {T} - S = sector_type(legs) - row_axis, codomain_trees_to_ranges = fuse_axes(S, codomain(legs)) - col_axis, domain_trees_to_ranges = flip_domain(fuse_axes(S, dual.(domain(legs)))...) - - mat = initialize_data_matrix(T, row_axis, col_axis) - tree_to_block_mapping = intersect_codomain_domain( - codomain_trees_to_ranges, domain_trees_to_ranges - ) - return FusionTensor(mat, legs, tree_to_block_mapping) + S = sector_type(legs) + row_axis, codomain_trees_to_ranges = fuse_axes(S, codomain(legs)) + col_axis, domain_trees_to_ranges = flip_domain(fuse_axes(S, dual.(domain(legs)))...) + + mat = initialize_data_matrix(T, row_axis, col_axis) + tree_to_block_mapping = intersect_codomain_domain( + codomain_trees_to_ranges, domain_trees_to_ranges + ) + return FusionTensor(mat, legs, tree_to_block_mapping) end #constructor from precomputed data_matrix function FusionTensor(mat::AbstractMatrix, legs::FusionTensorAxes) - # init with empty data_matrix to construct trees_block_mapping - ft = FusionTensor{eltype(mat)}(undef, legs) - for b in eachblockstoredindex(mat) - b in eachblockstoredindex(data_matrix(ft)) || - throw(ArgumentError("matrix block $b is not allowed")) - data_matrix(ft)[b] = mat[b] - end - return ft + # init with empty data_matrix to construct trees_block_mapping + ft = FusionTensor{eltype(mat)}(undef, legs) + for b in eachblockstoredindex(mat) + b in eachblockstoredindex(data_matrix(ft)) || + throw(ArgumentError("matrix block $b is not allowed")) + data_matrix(ft)[b] = mat[b] + end + return ft end FusionTensor(x, legs::BlockedTuple{2}) = FusionTensor(x, FusionTensorAxes(legs)) function FusionTensor{T}(x, legs::BlockedTuple{2}) where {T} - return FusionTensor{T}(x, FusionTensorAxes(legs)) + return FusionTensor{T}(x, FusionTensorAxes(legs)) end # constructor from split axes function FusionTensor( - x, - codomain_legs::Tuple{Vararg{AbstractGradedUnitRange}}, - domain_legs::Tuple{Vararg{AbstractGradedUnitRange}}, -) - return FusionTensor(x, tuplemortar((codomain_legs, domain_legs))) + x, + codomain_legs::Tuple{Vararg{AbstractGradedUnitRange}}, + domain_legs::Tuple{Vararg{AbstractGradedUnitRange}}, + ) + return FusionTensor(x, tuplemortar((codomain_legs, domain_legs))) end function FusionTensor{T}( - x, - codomain_legs::Tuple{Vararg{AbstractGradedUnitRange}}, - domain_legs::Tuple{Vararg{AbstractGradedUnitRange}}, -) where {T} - return FusionTensor{T}(x, tuplemortar((codomain_legs, domain_legs))) + x, + codomain_legs::Tuple{Vararg{AbstractGradedUnitRange}}, + domain_legs::Tuple{Vararg{AbstractGradedUnitRange}}, + ) where {T} + return FusionTensor{T}(x, tuplemortar((codomain_legs, domain_legs))) end # specific constructors function Base.zeros(::Type{T}, fta::FusionTensorAxes) where {T} - ft = FusionTensor{T}(undef, fta) - foreach(m -> fill!(m, zero(T)), eachstoredblock(data_matrix(ft))) - return ft + ft = FusionTensor{T}(undef, fta) + foreach(m -> fill!(m, zero(T)), eachstoredblock(data_matrix(ft))) + return ft end Base.zeros(fta::FusionTensorAxes) = zeros(Float64, fta) function Base.randn(rng::AbstractRNG, ::Type{T}, fta::FusionTensorAxes) where {T} - ft = FusionTensor{T}(undef, fta) - foreach(m -> randn!(rng, m), eachstoredblock(data_matrix(ft))) - return ft + ft = FusionTensor{T}(undef, fta) + foreach(m -> randn!(rng, m), eachstoredblock(data_matrix(ft))) + return ft end Base.randn(rng::AbstractRNG, fta::FusionTensorAxes) = randn(rng, Float64, fta) Base.randn(::Type{T}, fta::FusionTensorAxes) where {T} = randn(Random.default_rng(), T, fta) Base.randn(fta::FusionTensorAxes) = randn(Float64, fta) function FusionTensor{T}( - s::UniformScaling, codomain_legs::Tuple{Vararg{AbstractGradedUnitRange}} -) where {T} - fta = FusionTensorAxes(codomain_legs, dual.(codomain_legs)) - ft = FusionTensor{T}(undef, fta) - for m in eachstoredblock(data_matrix(ft)) - m .= s(size(m, 1)) - end - return ft + s::UniformScaling, codomain_legs::Tuple{Vararg{AbstractGradedUnitRange}} + ) where {T} + fta = FusionTensorAxes(codomain_legs, dual.(codomain_legs)) + ft = FusionTensor{T}(undef, fta) + for m in eachstoredblock(data_matrix(ft)) + m .= s(size(m, 1)) + end + return ft end function FusionTensor( - s::UniformScaling, codomain_legs::Tuple{Vararg{AbstractGradedUnitRange}} -) - return FusionTensor{Float64}(s, codomain_legs) + s::UniformScaling, codomain_legs::Tuple{Vararg{AbstractGradedUnitRange}} + ) + return FusionTensor{Float64}(s, codomain_legs) end # ================================ BlockArrays interface ================================= function BlockArrays.findblock(ft::FusionTensor, f1::SectorFusionTree, f2::SectorFusionTree) - # find outer block corresponding to fusion trees - @assert typeof((f1, f2)) === keytype(trees_block_mapping(ft)) - b1 = findfirstblock.(flip_dual.(codomain_axes(ft)), leaves(f1)) - b2 = findfirstblock.(flip_dual.(domain_axes(ft)), leaves(f2)) - return Block(Int.(b1)..., Int.(b2)...) + # find outer block corresponding to fusion trees + @assert typeof((f1, f2)) === keytype(trees_block_mapping(ft)) + b1 = findfirstblock.(flip_dual.(codomain_axes(ft)), leaves(f1)) + b2 = findfirstblock.(flip_dual.(domain_axes(ft)), leaves(f2)) + return Block(Int.(b1)..., Int.(b2)...) end # ============================== GradedArrays interface ================================== -function GradedArrays.sector_type(::Type{FT}) where {FT<:FusionTensor} - return sector_type(type_parameters(FT, 3)) +function GradedArrays.sector_type(::Type{FT}) where {FT <: FusionTensor} + return sector_type(type_parameters(FT, 3)) end -SymmetryStyle(::Type{FT}) where {FT<:FusionTensor} = SymmetryStyle(sector_type(FT)) +SymmetryStyle(::Type{FT}) where {FT <: FusionTensor} = SymmetryStyle(sector_type(FT)) # ============================== FusionTensor interface ================================== @@ -271,6 +271,6 @@ ndims_codomain(ft::FusionTensor) = length_codomain(axes(ft)) ndims_domain(ft::FusionTensor) = length_domain(axes(ft)) function charge_block_size(ft::FusionTensor, f1::SectorFusionTree, f2::SectorFusionTree) - b = Tuple(findblock(ft, f1, f2)) - return ntuple(i -> sector_multiplicity(axes(ft, i)[b[i]]), ndims(ft)) + b = Tuple(findblock(ft, f1, f2)) + return ntuple(i -> sector_multiplicity(axes(ft, i)[b[i]]), ndims(ft)) end diff --git a/src/fusiontensor/fusiontensoraxes.jl b/src/fusiontensor/fusiontensoraxes.jl index 46d2478..0e999f5 100644 --- a/src/fusiontensor/fusiontensoraxes.jl +++ b/src/fusiontensor/fusiontensoraxes.jl @@ -1,20 +1,20 @@ using BlockArrays: BlockArrays using GradedArrays: - GradedArrays, - AbstractGradedUnitRange, - AbstractSector, - SymmetryStyle, - TrivialSector, - dual, - sector_type, - trivial + GradedArrays, + AbstractGradedUnitRange, + AbstractSector, + SymmetryStyle, + TrivialSector, + dual, + sector_type, + trivial using TensorAlgebra: - TensorAlgebra, - AbstractBlockPermutation, - AbstractBlockTuple, - BlockedTuple, - length_codomain, - length_domain + TensorAlgebra, + AbstractBlockPermutation, + AbstractBlockTuple, + BlockedTuple, + length_codomain, + length_domain using TensorProducts: ⊗ using TypeParameterAccessors: type_parameters @@ -23,46 +23,46 @@ using TypeParameterAccessors: type_parameters promote_sector_type(legs::Tuple) = promote_sector_type(legs...) function promote_sector_type(legs...) - # fuse trivial sectors to produce unified type - # avoid depending on GradedArrays internals - return sector_type(⊗(trivial.(legs)...)) + # fuse trivial sectors to produce unified type + # avoid depending on GradedArrays internals + return sector_type(⊗(trivial.(legs)...)) end -promote_sectors(legs::NTuple{<:Any,<:AbstractGradedUnitRange}) = legs # nothing to do +promote_sectors(legs::NTuple{<:Any, <:AbstractGradedUnitRange}) = legs # nothing to do function promote_sectors(legs) - T = promote_sector_type(legs) - # fuse with trivial to insert all missing arguments inside each GradedAxis - # avoid depending on GradedArrays internals - s0 = trivial(T) - return map_sectors.(s -> only(sectors(to_gradedrange(tensor_product(s0, s)))), legs) + T = promote_sector_type(legs) + # fuse with trivial to insert all missing arguments inside each GradedAxis + # avoid depending on GradedArrays internals + s0 = trivial(T) + return map_sectors.(s -> only(sectors(to_gradedrange(tensor_product(s0, s)))), legs) end function promote_sectors(bt::BlockedTuple) - return BlockedTuple{blocklength(bt),blocklengths(bt)}(promote_sectors(Tuple(bt))) + return BlockedTuple{blocklength(bt), blocklengths(bt)}(promote_sectors(Tuple(bt))) end # ==================================== Definitions ======================================= # TBD explicit axis type as type parameters? -struct FusionTensorAxes{BT<:BlockedTuple{2}} - outer_axes::BT +struct FusionTensorAxes{BT <: BlockedTuple{2}} + outer_axes::BT - function FusionTensorAxes{BT}(bt) where {BT} - @assert BT === typeof(promote_sectors(bt)) - return new{BT}(bt) - end + function FusionTensorAxes{BT}(bt) where {BT} + @assert BT === typeof(promote_sectors(bt)) + return new{BT}(bt) + end end # ==================================== Constructors ====================================== function FusionTensorAxes(bt::BlockedTuple{2}) - promoted = promote_sectors(bt) - return FusionTensorAxes{typeof(promoted)}(promoted) + promoted = promote_sectors(bt) + return FusionTensorAxes{typeof(promoted)}(promoted) end function FusionTensorAxes(codomain_legs, domain_legs) - return FusionTensorAxes(tuplemortar((codomain_legs, domain_legs))) + return FusionTensorAxes(tuplemortar((codomain_legs, domain_legs))) end # ============================== TensorAlgebra interface ================================= @@ -76,17 +76,17 @@ TensorAlgebra.length_domain(fta::FusionTensorAxes) = length(domain(fta)) # ================================== Base interface ====================================== for f in [ - :(broadcastable), :(Tuple), :(axes), :(firstindex), :(lastindex), :(iterate), :(length) -] - @eval Base.$f(fta::FusionTensorAxes) = Base.$f(BlockedTuple(fta)) + :(broadcastable), :(Tuple), :(axes), :(firstindex), :(lastindex), :(iterate), :(length), + ] + @eval Base.$f(fta::FusionTensorAxes) = Base.$f(BlockedTuple(fta)) end for f in [:(getindex), :(iterate)] - @eval Base.$f(fta::FusionTensorAxes, i) = $f(BlockedTuple(fta), i) + @eval Base.$f(fta::FusionTensorAxes, i) = $f(BlockedTuple(fta), i) end function Base.getindex(fta::FusionTensorAxes, bp::AbstractBlockPermutation) - return FusionTensorAxes(BlockedTuple(fta)[bp]) + return FusionTensorAxes(BlockedTuple(fta)[bp]) end Base.copy(fta::FusionTensorAxes) = FusionTensorAxes(copy.(BlockedTuple(fta))) @@ -94,39 +94,39 @@ Base.copy(fta::FusionTensorAxes) = FusionTensorAxes(copy.(BlockedTuple(fta))) Base.deepcopy(fta::FusionTensorAxes) = FusionTensorAxes(deepcopy.(BlockedTuple(fta))) function Base.:(==)(a::FusionTensorAxes, b::FusionTensorAxes) - blocklengths(a) != blocklengths(b) && return false - for i in 1:length(a) - !space_isequal(a[i], b[i]) && return false - end - return true + blocklengths(a) != blocklengths(b) && return false + for i in 1:length(a) + !space_isequal(a[i], b[i]) && return false + end + return true end # ================================ BlockArrays interface ================================= for f in [:(blocklength), :(blocklengths), :(blocks)] - @eval BlockArrays.$f(fta::FusionTensorAxes) = $f(BlockedTuple(fta)) + @eval BlockArrays.$f(fta::FusionTensorAxes) = $f(BlockedTuple(fta)) end # ============================== GradedArrays interface ================================== function GradedArrays.sector_type( - ::Type{FTA} -) where {BT<:BlockedTuple{2,(0, 0)},FTA<:FusionTensorAxes{BT}} - return TrivialSector + ::Type{FTA} + ) where {BT <: BlockedTuple{2, (0, 0)}, FTA <: FusionTensorAxes{BT}} + return TrivialSector end -function GradedArrays.sector_type(::Type{FTA}) where {BT,FTA<:FusionTensorAxes{BT}} - return sector_type(type_parameters(type_parameters(BT, 3), 1)) +function GradedArrays.sector_type(::Type{FTA}) where {BT, FTA <: FusionTensorAxes{BT}} + return sector_type(type_parameters(type_parameters(BT, 3), 1)) end -function GradedArrays.SymmetryStyle(::Type{FTA}) where {FTA<:FusionTensorAxes} - return SymmetryStyle(sector_type(FTA)) +function GradedArrays.SymmetryStyle(::Type{FTA}) where {FTA <: FusionTensorAxes} + return SymmetryStyle(sector_type(FTA)) end function GradedArrays.checkspaces( - ::Type{Bool}, left::FusionTensorAxes, right::FusionTensorAxes -) - return left == right + ::Type{Bool}, left::FusionTensorAxes, right::FusionTensorAxes + ) + return left == right end # ============================== FusionTensor interface ================================== @@ -136,15 +136,15 @@ codomain(fta::FusionTensorAxes) = fta[Block(1)] domain(fta::FusionTensorAxes) = fta[Block(2)] function fused_codomain(fta::FusionTensorAxes) - if length_codomain(fta) == 0 - return trivial_axis(fta) - end - return ⊗(codomain(fta)...) + if length_codomain(fta) == 0 + return trivial_axis(fta) + end + return ⊗(codomain(fta)...) end function fused_domain(fta::FusionTensorAxes) - if length_domain(fta) == 0 - return dual(trivial_axis(fta)) - end - return dual(⊗(dual.(domain(fta))...)) + if length_domain(fta) == 0 + return dual(trivial_axis(fta)) + end + return dual(⊗(dual.(domain(fta))...)) end diff --git a/src/fusiontensor/linear_algebra_interface.jl b/src/fusiontensor/linear_algebra_interface.jl index 2d2a7ca..cd76156 100644 --- a/src/fusiontensor/linear_algebra_interface.jl +++ b/src/fusiontensor/linear_algebra_interface.jl @@ -11,48 +11,48 @@ using GradedArrays: checkspaces, checkspaces_dual, quantum_dimension, sectors # impose matching type and number of axes at compile time # impose matching axes at run time function LinearAlgebra.mul!( - C::FusionMatrix, A::FusionMatrix, B::FusionMatrix, α::Number, β::Number -) - - # compile time checks - if ndims_domain(A) != ndims_codomain(B) - throw(codomainError("Incompatible tensor structures for A and B")) - end - if ndims_codomain(A) != ndims_codomain(C) - throw(codomainError("Incompatible tensor structures for A and C")) - end - if ndims_domain(B) != ndims_domain(C) - throw(codomainError("Incompatible tensor structures for B and C")) - end - - # input validation - checkspaces_dual(domain_axes(A), codomain_axes(B)) - checkspaces(codomain_axes(C), codomain_axes(A)) - checkspaces(domain_axes(C), domain_axes(B)) - mul!(data_matrix(C), data_matrix(A), data_matrix(B), α, β) - return C + C::FusionMatrix, A::FusionMatrix, B::FusionMatrix, α::Number, β::Number + ) + + # compile time checks + if ndims_domain(A) != ndims_codomain(B) + throw(codomainError("Incompatible tensor structures for A and B")) + end + if ndims_codomain(A) != ndims_codomain(C) + throw(codomainError("Incompatible tensor structures for A and C")) + end + if ndims_domain(B) != ndims_domain(C) + throw(codomainError("Incompatible tensor structures for B and C")) + end + + # input validation + checkspaces_dual(domain_axes(A), codomain_axes(B)) + checkspaces(codomain_axes(C), codomain_axes(A)) + checkspaces(domain_axes(C), domain_axes(B)) + mul!(data_matrix(C), data_matrix(A), data_matrix(B), α, β) + return C end -function LinearAlgebra.norm(ft::FusionTensor, p::Real=2) - m = data_matrix(ft) - row_sectors = sectors(codomain_axis(ft)) - np = sum(eachblockstoredindex(m); init=zero(real(eltype(ft)))) do b - return quantum_dimension(row_sectors[Int(first(Tuple(b)))]) * norm(m[b], p)^p - end - return np^(1 / p) +function LinearAlgebra.norm(ft::FusionTensor, p::Real = 2) + m = data_matrix(ft) + row_sectors = sectors(codomain_axis(ft)) + np = sum(eachblockstoredindex(m); init = zero(real(eltype(ft)))) do b + return quantum_dimension(row_sectors[Int(first(Tuple(b)))]) * norm(m[b], p)^p + end + return np^(1 / p) end -LinearAlgebra.normalize(ft::FusionTensor, p::Real=2) = ft / norm(ft, p) +LinearAlgebra.normalize(ft::FusionTensor, p::Real = 2) = ft / norm(ft, p) -function LinearAlgebra.normalize!(ft::FusionTensor, p::Real=2) - data_matrix(ft) ./= norm(ft, p) - return ft +function LinearAlgebra.normalize!(ft::FusionTensor, p::Real = 2) + data_matrix(ft) ./= norm(ft, p) + return ft end function LinearAlgebra.tr(ft::FusionTensor) - m = data_matrix(ft) - row_sectors = sectors(codomain_axis(ft)) - return sum(eachblockstoredindex(m); init=zero(eltype(ft))) do b - return quantum_dimension(row_sectors[Int(first(Tuple(b)))]) * tr(m[b]) - end + m = data_matrix(ft) + row_sectors = sectors(codomain_axis(ft)) + return sum(eachblockstoredindex(m); init = zero(eltype(ft))) do b + return quantum_dimension(row_sectors[Int(first(Tuple(b)))]) * tr(m[b]) + end end diff --git a/src/fusiontensor/tensor_algebra_interface.jl b/src/fusiontensor/tensor_algebra_interface.jl index 1c29973..c5f3fd3 100644 --- a/src/fusiontensor/tensor_algebra_interface.jl +++ b/src/fusiontensor/tensor_algebra_interface.jl @@ -6,67 +6,69 @@ using BlockArrays: Block using GradedArrays: space_isequal using TensorAlgebra: - TensorAlgebra, - AbstractBlockPermutation, - BlockedTrivialPermutation, - BlockedTuple, - FusionStyle, - Matricize, - blockedperm, - genperm, - matricize, - unmatricize + TensorAlgebra, + AbstractBlockPermutation, + BlockedTrivialPermutation, + BlockedTuple, + FusionStyle, + Matricize, + blockedperm, + genperm, + matricize, + unmatricize const MATRIX_FUNCTIONS = [ - :exp, - :cis, - :log, - :sqrt, - :cbrt, - :cos, - :sin, - :tan, - :csc, - :sec, - :cot, - :cosh, - :sinh, - :tanh, - :csch, - :sech, - :coth, - :acos, - :asin, - :atan, - :acsc, - :asec, - :acot, - :acosh, - :asinh, - :atanh, - :acsch, - :asech, - :acoth, + :exp, + :cis, + :log, + :sqrt, + :cbrt, + :cos, + :sin, + :tan, + :csc, + :sec, + :cot, + :cosh, + :sinh, + :tanh, + :csch, + :sech, + :coth, + :acos, + :asin, + :atan, + :acsc, + :asec, + :acot, + :acosh, + :asinh, + :atanh, + :acsch, + :asech, + :acoth, ] function TensorAlgebra.output_axes( - ::typeof(contract), - biperm_dest::AbstractBlockPermutation{2}, - a1::FusionTensor, - biperm1::AbstractBlockPermutation{2}, - a2::FusionTensor, - biperm2::AbstractBlockPermutation{2}, -) - axes_codomain, axes_contracted = blocks(axes(a1)[biperm1]) - axes_contracted2, axes_domain = blocks(axes(a2)[biperm2]) - @assert all(space_isequal.(dual.(axes_contracted), axes_contracted2)) - flat_axes = genperm((axes_codomain..., axes_domain...), Tuple(biperm_dest)) - return FusionTensorAxes( - tuplemortar(( - flat_axes[begin:length_codomain(biperm_dest)], - flat_axes[(length_codomain(biperm_dest) + 1):end], - )), - ) + ::typeof(contract), + biperm_dest::AbstractBlockPermutation{2}, + a1::FusionTensor, + biperm1::AbstractBlockPermutation{2}, + a2::FusionTensor, + biperm2::AbstractBlockPermutation{2}, + ) + axes_codomain, axes_contracted = blocks(axes(a1)[biperm1]) + axes_contracted2, axes_domain = blocks(axes(a2)[biperm2]) + @assert all(space_isequal.(dual.(axes_contracted), axes_contracted2)) + flat_axes = genperm((axes_codomain..., axes_domain...), Tuple(biperm_dest)) + return FusionTensorAxes( + tuplemortar( + ( + flat_axes[begin:length_codomain(biperm_dest)], + flat_axes[(length_codomain(biperm_dest) + 1):end], + ) + ), + ) end struct FusionTensorFusionStyle <: FusionStyle end @@ -74,46 +76,46 @@ struct FusionTensorFusionStyle <: FusionStyle end TensorAlgebra.FusionStyle(::Type{<:FusionTensor}) = FusionTensorFusionStyle() function TensorAlgebra.matricize( - ::FusionTensorFusionStyle, ft::AbstractArray, biperm::BlockedTrivialPermutation{2} -) - blocklengths(biperm) == blocklengths(axes(ft)) || - throw(ArgumentError("Invalid trivial biperm")) - return FusionTensor(data_matrix(ft), (codomain_axis(ft),), (domain_axis(ft),)) + ::FusionTensorFusionStyle, ft::AbstractArray, biperm::BlockedTrivialPermutation{2} + ) + blocklengths(biperm) == blocklengths(axes(ft)) || + throw(ArgumentError("Invalid trivial biperm")) + return FusionTensor(data_matrix(ft), (codomain_axis(ft),), (domain_axis(ft),)) end function TensorAlgebra.unmatricize(::FusionTensorFusionStyle, m, blocked_axes) - return FusionTensor(data_matrix(m), blocked_axes) + return FusionTensor(data_matrix(m), blocked_axes) end function TensorAlgebra.permuteblockeddims( - ft::FusionTensor, biperm::AbstractBlockPermutation -) - return permutedims(ft, biperm) + ft::FusionTensor, biperm::AbstractBlockPermutation + ) + return permutedims(ft, biperm) end function TensorAlgebra.permuteblockeddims!( - a::FusionTensor, b::FusionTensor, biperm::AbstractBlockPermutation -) - return permutedims!(a, b, biperm) + a::FusionTensor, b::FusionTensor, biperm::AbstractBlockPermutation + ) + return permutedims!(a, b, biperm) end # TODO define custom broadcast rules function TensorAlgebra.unmatricizeadd!(a_dest::FusionTensor, a_dest_mat, invbiperm, α, β) - a12 = unmatricize(a_dest_mat, axes(a_dest), invbiperm) - data_matrix(a_dest) .= α .* data_matrix(a12) .+ β .* data_matrix(a_dest) - return a_dest + a12 = unmatricize(a_dest_mat, axes(a_dest), invbiperm) + data_matrix(a_dest) .= α .* data_matrix(a12) .+ β .* data_matrix(a_dest) + return a_dest end for f in MATRIX_FUNCTIONS - @eval begin - function TensorAlgebra.$f( - a::FusionTensor, biperm::AbstractBlockPermutation{2}; kwargs... - ) - a_mat = matricize(a, biperm) - permuted_axes = axes(a)[biperm] - checkspaces_dual(codomain(permuted_axes), domain(permuted_axes)) - fa_mat = set_data_matrix(a_mat, Base.$f(data_matrix(a_mat); kwargs...)) - return unmatricize(fa_mat, permuted_axes) + @eval begin + function TensorAlgebra.$f( + a::FusionTensor, biperm::AbstractBlockPermutation{2}; kwargs... + ) + a_mat = matricize(a, biperm) + permuted_axes = axes(a)[biperm] + checkspaces_dual(codomain(permuted_axes), domain(permuted_axes)) + fa_mat = set_data_matrix(a_mat, Base.$f(data_matrix(a_mat); kwargs...)) + return unmatricize(fa_mat, permuted_axes) + end end - end end diff --git a/src/permutedims/permutedims.jl b/src/permutedims/permutedims.jl index 1208df4..f24cae3 100644 --- a/src/permutedims/permutedims.jl +++ b/src/permutedims/permutedims.jl @@ -7,71 +7,71 @@ using GradedArrays: AbelianStyle, NotAbelianStyle, SymmetryStyle, checkspaces using TensorAlgebra: AbstractBlockPermutation, permmortar # permutedims with 1 tuple of 2 separate tuples -function fusiontensor_permutedims(ft, new_leg_dims::Tuple{Tuple,Tuple}) - return fusiontensor_permutedims(ft, new_leg_dims...) +function fusiontensor_permutedims(ft, new_leg_dims::Tuple{Tuple, Tuple}) + return fusiontensor_permutedims(ft, new_leg_dims...) end -function fusiontensor_permutedims!(ftdst, ftsrc, new_leg_dims::Tuple{Tuple,Tuple}) - return fusiontensor_permutedims!(ftdst, ftsrc, new_leg_dims...) +function fusiontensor_permutedims!(ftdst, ftsrc, new_leg_dims::Tuple{Tuple, Tuple}) + return fusiontensor_permutedims!(ftdst, ftsrc, new_leg_dims...) end # permutedims with 2 separate tuples function fusiontensor_permutedims(ft, new_codomain_dims::Tuple, new_domain_dims::Tuple) - biperm = permmortar((new_codomain_dims, new_domain_dims)) - return fusiontensor_permutedims(ft, biperm) + biperm = permmortar((new_codomain_dims, new_domain_dims)) + return fusiontensor_permutedims(ft, biperm) end function fusiontensor_permutedims!( - ftdst, ftsrc, new_codomain_dims::Tuple, new_domain_dims::Tuple -) - biperm = permmortar((new_codomain_dims, new_domain_dims)) - return fusiontensor_permutedims!(ftdst, ftsrc, biperm) + ftdst, ftsrc, new_codomain_dims::Tuple, new_domain_dims::Tuple + ) + biperm = permmortar((new_codomain_dims, new_domain_dims)) + return fusiontensor_permutedims!(ftdst, ftsrc, biperm) end # permutedims with BlockedPermutation function fusiontensor_permutedims(ft, biperm::AbstractBlockPermutation{2}) - ndims(ft) == length(biperm) || throw(ArgumentError("Invalid permutation length")) - ftdst = similar(ft, axes(ft)[biperm]) - fusiontensor_permutedims!(ftdst, ft, biperm) - return ftdst + ndims(ft) == length(biperm) || throw(ArgumentError("Invalid permutation length")) + ftdst = similar(ft, axes(ft)[biperm]) + fusiontensor_permutedims!(ftdst, ft, biperm) + return ftdst end function fusiontensor_permutedims!(ftdst, ftsrc, biperm::AbstractBlockPermutation{2}) - ndims(ftsrc) == length(biperm) || throw(ArgumentError("Invalid permutation length")) - blocklengths(axes(ftdst)) == blocklengths(biperm) || - throw(ArgumentError("Destination tensor does not match bipermutation")) - checkspaces(axes(ftdst), axes(ftsrc)[biperm]) + ndims(ftsrc) == length(biperm) || throw(ArgumentError("Invalid permutation length")) + blocklengths(axes(ftdst)) == blocklengths(biperm) || + throw(ArgumentError("Destination tensor does not match bipermutation")) + checkspaces(axes(ftdst), axes(ftsrc)[biperm]) - # early return for identity operation. Also handle tricky 0-dim case. - if ndims_codomain(ftdst) == ndims_codomain(ftsrc) # compile time - if Tuple(biperm) == ntuple(identity, ndims(ftdst)) - copy!(data_matrix(ftdst), data_matrix(ftsrc)) - return ftdst + # early return for identity operation. Also handle tricky 0-dim case. + if ndims_codomain(ftdst) == ndims_codomain(ftsrc) # compile time + if Tuple(biperm) == ntuple(identity, ndims(ftdst)) + copy!(data_matrix(ftdst), data_matrix(ftsrc)) + return ftdst + end end - end - return _fusiontensor_permutedims!(SymmetryStyle(ftdst), ftdst, ftsrc, Tuple(biperm)) + return _fusiontensor_permutedims!(SymmetryStyle(ftdst), ftdst, ftsrc, Tuple(biperm)) end # =============================== Internal ============================================= function _fusiontensor_permutedims!(::AbelianStyle, ftdst, ftsrc, flatperm) - # abelian case: all unitary blocks are 1x1 identity matrices - # compute_unitary is only called to get block positions - unitary = compute_unitary(ftdst, ftsrc, flatperm) - for ((old_trees, new_trees), _) in unitary - new_block = view(ftdst, new_trees...) - old_block = view(ftsrc, old_trees...) - @strided new_block .= permutedims(old_block, flatperm) - end - return ftdst + # abelian case: all unitary blocks are 1x1 identity matrices + # compute_unitary is only called to get block positions + unitary = compute_unitary(ftdst, ftsrc, flatperm) + for ((old_trees, new_trees), _) in unitary + new_block = view(ftdst, new_trees...) + old_block = view(ftsrc, old_trees...) + @strided new_block .= permutedims(old_block, flatperm) + end + return ftdst end function _fusiontensor_permutedims!(::NotAbelianStyle, ftdst, ftsrc, flatperm) - foreach(m -> fill!(m, zero(eltype(ftdst))), eachstoredblock(data_matrix(ftdst))) - unitary = compute_unitary(ftdst, ftsrc, flatperm) - for ((old_trees, new_trees), coeff) in unitary - new_block = view(ftdst, new_trees...) - old_block = view(ftsrc, old_trees...) - @strided new_block .+= coeff .* permutedims(old_block, flatperm) - end - return ftdst + foreach(m -> fill!(m, zero(eltype(ftdst))), eachstoredblock(data_matrix(ftdst))) + unitary = compute_unitary(ftdst, ftsrc, flatperm) + for ((old_trees, new_trees), coeff) in unitary + new_block = view(ftdst, new_trees...) + old_block = view(ftsrc, old_trees...) + @strided new_block .+= coeff .* permutedims(old_block, flatperm) + end + return ftdst end diff --git a/src/permutedims/unitaries.jl b/src/permutedims/unitaries.jl index 6fbc6c7..d18f8af 100644 --- a/src/permutedims/unitaries.jl +++ b/src/permutedims/unitaries.jl @@ -5,52 +5,52 @@ using LRUCache: LRU using GradedArrays: quantum_dimension -const unitary_cache = LRU{Any,Float64}(; maxsize=10000) # TBD size +const unitary_cache = LRU{Any, Float64}(; maxsize = 10000) # TBD size # ====================================== Interface ======================================= function compute_unitary( - new_ft::FusionTensor{T,N}, old_ft::FusionTensor{T,N}, flatperm::NTuple{N,Int} -) where {T,N} - return compute_unitary_clebsch_gordan(new_ft, old_ft, flatperm) + new_ft::FusionTensor{T, N}, old_ft::FusionTensor{T, N}, flatperm::NTuple{N, Int} + ) where {T, N} + return compute_unitary_clebsch_gordan(new_ft, old_ft, flatperm) end # =========================== Constructor from Clebsch-Gordan ============================ function overlap_fusion_trees( - old_trees::Tuple{SectorFusionTree{S},SectorFusionTree{S}}, - new_trees::Tuple{SectorFusionTree{S},SectorFusionTree{S}}, - flatperm::Tuple{Vararg{Integer}}, -) where {S} - old_proj = contract_singlet_projector(old_trees...) - new_proj = contract_singlet_projector(new_trees...) - a = contract((), new_proj, flatperm, old_proj, ntuple(identity, ndims(new_proj))) - return a[] / quantum_dimension(root_sector(first(new_trees))) + old_trees::Tuple{SectorFusionTree{S}, SectorFusionTree{S}}, + new_trees::Tuple{SectorFusionTree{S}, SectorFusionTree{S}}, + flatperm::Tuple{Vararg{Integer}}, + ) where {S} + old_proj = contract_singlet_projector(old_trees...) + new_proj = contract_singlet_projector(new_trees...) + a = contract((), new_proj, flatperm, old_proj, ntuple(identity, ndims(new_proj))) + return a[] / quantum_dimension(root_sector(first(new_trees))) end function cached_unitary_coeff( - old_trees::Tuple{SectorFusionTree{S},SectorFusionTree{S}}, - new_trees::Tuple{SectorFusionTree{S},SectorFusionTree{S}}, - flatperm::Tuple{Vararg{Integer}}, -) where {S} - get!(unitary_cache, (old_trees..., new_trees..., flatperm)) do - overlap_fusion_trees(old_trees, new_trees, flatperm) - end + old_trees::Tuple{SectorFusionTree{S}, SectorFusionTree{S}}, + new_trees::Tuple{SectorFusionTree{S}, SectorFusionTree{S}}, + flatperm::Tuple{Vararg{Integer}}, + ) where {S} + return get!(unitary_cache, (old_trees..., new_trees..., flatperm)) do + overlap_fusion_trees(old_trees, new_trees, flatperm) + end end function compute_unitary_clebsch_gordan( - new_ft::FusionTensor{T,N}, old_ft::FusionTensor{T,N}, flatperm::NTuple{N,Int} -) where {T,N} - unitary = Dict{ - Tuple{keytype(trees_block_mapping(old_ft)),keytype(trees_block_mapping(new_ft))},Float64 - }() - for old_trees in keys(trees_block_mapping(old_ft)) - old_outer = Tuple(findblock(old_ft, old_trees...)) - swapped_old_block = Block(getindex.(Ref(Tuple(old_outer)), flatperm)) - for new_trees in keys(trees_block_mapping(new_ft)) - swapped_old_block != findblock(new_ft, new_trees...) && continue - unitary[old_trees, new_trees] = cached_unitary_coeff(old_trees, new_trees, flatperm) + new_ft::FusionTensor{T, N}, old_ft::FusionTensor{T, N}, flatperm::NTuple{N, Int} + ) where {T, N} + unitary = Dict{ + Tuple{keytype(trees_block_mapping(old_ft)), keytype(trees_block_mapping(new_ft))}, Float64, + }() + for old_trees in keys(trees_block_mapping(old_ft)) + old_outer = Tuple(findblock(old_ft, old_trees...)) + swapped_old_block = Block(getindex.(Ref(Tuple(old_outer)), flatperm)) + for new_trees in keys(trees_block_mapping(new_ft)) + swapped_old_block != findblock(new_ft, new_trees...) && continue + unitary[old_trees, new_trees] = cached_unitary_coeff(old_trees, new_trees, flatperm) + end end - end - return unitary + return unitary end # ================================= Constructor from 6j ================================== diff --git a/test/runtests.jl b/test/runtests.jl index 98b2d2b..0008050 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -6,60 +6,62 @@ using Suppressor: Suppressor const pat = r"(?:--group=)(\w+)" arg_id = findfirst(contains(pat), ARGS) const GROUP = uppercase( - if isnothing(arg_id) - get(ENV, "GROUP", "ALL") - else - only(match(pat, ARGS[arg_id]).captures) - end, + if isnothing(arg_id) + get(ENV, "GROUP", "ALL") + else + only(match(pat, ARGS[arg_id]).captures) + end, ) "match files of the form `test_*.jl`, but exclude `*setup*.jl`" function istestfile(fn) - return endswith(fn, ".jl") && startswith(basename(fn), "test_") && !contains(fn, "setup") + return endswith(fn, ".jl") && startswith(basename(fn), "test_") && !contains(fn, "setup") end "match files of the form `*.jl`, but exclude `*_notest.jl` and `*setup*.jl`" function isexamplefile(fn) - return endswith(fn, ".jl") && !endswith(fn, "_notest.jl") && !contains(fn, "setup") + return endswith(fn, ".jl") && !endswith(fn, "_notest.jl") && !contains(fn, "setup") end @time begin - # tests in groups based on folder structure - for testgroup in filter(isdir, readdir(@__DIR__)) - if GROUP == "ALL" || GROUP == uppercase(testgroup) - groupdir = joinpath(@__DIR__, testgroup) - for file in filter(istestfile, readdir(groupdir)) - filename = joinpath(groupdir, file) - @eval @safetestset $file begin - include($filename) + # tests in groups based on folder structure + for testgroup in filter(isdir, readdir(@__DIR__)) + if GROUP == "ALL" || GROUP == uppercase(testgroup) + groupdir = joinpath(@__DIR__, testgroup) + for file in filter(istestfile, readdir(groupdir)) + filename = joinpath(groupdir, file) + @eval @safetestset $file begin + include($filename) + end + end end - end end - end - # single files in top folder - for file in filter(istestfile, readdir(@__DIR__)) - (file == basename(@__FILE__)) && continue # exclude this file to avoid infinite recursion - @eval @safetestset $file begin - include($file) + # single files in top folder + for file in filter(istestfile, readdir(@__DIR__)) + (file == basename(@__FILE__)) && continue # exclude this file to avoid infinite recursion + @eval @safetestset $file begin + include($file) + end end - end - # test examples - examplepath = joinpath(@__DIR__, "..", "examples") - for (root, _, files) in walkdir(examplepath) - contains(chopprefix(root, @__DIR__), "setup") && continue - for file in filter(isexamplefile, files) - filename = joinpath(root, file) - @eval begin - @safetestset $file begin - $(Expr( - :macrocall, - GlobalRef(Suppressor, Symbol("@suppress")), - LineNumberNode(@__LINE__, @__FILE__), - :(include($filename)), - )) + # test examples + examplepath = joinpath(@__DIR__, "..", "examples") + for (root, _, files) in walkdir(examplepath) + contains(chopprefix(root, @__DIR__), "setup") && continue + for file in filter(isexamplefile, files) + filename = joinpath(root, file) + @eval begin + @safetestset $file begin + $( + Expr( + :macrocall, + GlobalRef(Suppressor, Symbol("@suppress")), + LineNumberNode(@__LINE__, @__FILE__), + :(include($filename)), + ) + ) + end + end end - end end - end end diff --git a/test/setup.jl b/test/setup.jl index b36de0e..f4b0063 100644 --- a/test/setup.jl +++ b/test/setup.jl @@ -2,41 +2,41 @@ using LinearAlgebra: Adjoint using BlockArrays: blocklengths using BlockSparseArrays: BlockSparseMatrix, eachblockstoredindex using FusionTensors: - FusionTensor, - codomain_axes, - data_matrix, - domain_axes, - domain_axis, - codomain_axis, - ndims_codomain, - ndims_domain + FusionTensor, + codomain_axes, + data_matrix, + domain_axes, + domain_axis, + codomain_axis, + ndims_codomain, + ndims_domain using GradedArrays: - checkspaces, checkspaces_dual, dual, sectors, sector_multiplicities, space_isequal + checkspaces, checkspaces_dual, dual, sectors, sector_multiplicities, space_isequal function check_sanity(ft::FusionTensor) - nca = ndims_domain(ft) - @assert nca == length(domain_axes(ft)) "ndims_domain does not match domain_axes" - @assert nca <= ndims(ft) "invalid ndims_domain" + nca = ndims_domain(ft) + @assert nca == length(domain_axes(ft)) "ndims_domain does not match domain_axes" + @assert nca <= ndims(ft) "invalid ndims_domain" - nda = ndims_codomain(ft) - @assert nda == length(codomain_axes(ft)) "ndims_codomain does not match codomain_axes" - @assert nda <= ndims(ft) "invalid ndims_codomain" - @assert nda + nca == ndims(ft) "invalid ndims" + nda = ndims_codomain(ft) + @assert nda == length(codomain_axes(ft)) "ndims_codomain does not match codomain_axes" + @assert nda <= ndims(ft) "invalid ndims_codomain" + @assert nda + nca == ndims(ft) "invalid ndims" - @assert length(axes(ft)) == ndims(ft) "ndims does not match axes" - checkspaces(axes(ft)[begin:nda], codomain_axes(ft)) - checkspaces(axes(ft)[(nda + 1):end], domain_axes(ft)) + @assert length(axes(ft)) == ndims(ft) "ndims does not match axes" + checkspaces(axes(ft)[begin:nda], codomain_axes(ft)) + checkspaces(axes(ft)[(nda + 1):end], domain_axes(ft)) - m = data_matrix(ft) - @assert ndims(m) == 2 "invalid data_matrix ndims" - row_axis = codomain_axis(ft) - column_axis = domain_axis(ft) - @assert sector_multiplicities(row_axis) == blocklengths(axes(m, 1)) "invalid row_axis" - @assert sector_multiplicities(column_axis) == blocklengths(axes(m, 2)) "invalid column_axis" + m = data_matrix(ft) + @assert ndims(m) == 2 "invalid data_matrix ndims" + row_axis = codomain_axis(ft) + column_axis = domain_axis(ft) + @assert sector_multiplicities(row_axis) == blocklengths(axes(m, 1)) "invalid row_axis" + @assert sector_multiplicities(column_axis) == blocklengths(axes(m, 2)) "invalid column_axis" - for b in eachblockstoredindex(m) - ir, ic = Int.(Tuple(b)) - @assert sectors(row_axis)[ir] == sectors(dual(column_axis))[ic] "forbidden block" - end - return nothing + for b in eachblockstoredindex(m) + ir, ic = Int.(Tuple(b)) + @assert sectors(row_axis)[ir] == sectors(dual(column_axis))[ic] "forbidden block" + end + return nothing end diff --git a/test/test_aqua.jl b/test/test_aqua.jl index ba30be2..0ff7f63 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(FusionTensors) + Aqua.test_all(FusionTensors) end diff --git a/test/test_array_cast.jl b/test/test_array_cast.jl index 049b0d7..b3fddbc 100644 --- a/test/test_array_cast.jl +++ b/test/test_array_cast.jl @@ -10,402 +10,404 @@ using TensorProducts: tensor_product include("setup.jl") @testset "Trivial FusionTensor" begin - @testset "trivial matrix" begin - g = gradedrange([TrivialSector() => 1]) - gb = dual(g) - m = ones((1, 1)) - ft = to_fusiontensor(m, (g,), (gb,)) - @test size(data_matrix(ft)) == (1, 1) - @test blocksize(data_matrix(ft)) == (1, 1) - @test data_matrix(ft)[1, 1] ≈ 1.0 - @test isnothing(check_sanity(ft)) - @test Array(ft) ≈ m - @test Array(adjoint(ft)) ≈ m - - for elt in (Int, UInt32, Float32) - m = ones(elt, (1, 1)) - ft = to_fusiontensor(m, (g,), (gb,)) - @test eltype(ft) === Float64 - @test Array(ft) ≈ m + @testset "trivial matrix" begin + g = gradedrange([TrivialSector() => 1]) + gb = dual(g) + m = ones((1, 1)) + ft = to_fusiontensor(m, (g,), (gb,)) + @test size(data_matrix(ft)) == (1, 1) + @test blocksize(data_matrix(ft)) == (1, 1) + @test data_matrix(ft)[1, 1] ≈ 1.0 + @test isnothing(check_sanity(ft)) + @test Array(ft) ≈ m + @test Array(adjoint(ft)) ≈ m + + for elt in (Int, UInt32, Float32) + m = ones(elt, (1, 1)) + ft = to_fusiontensor(m, (g,), (gb,)) + @test eltype(ft) === Float64 + @test Array(ft) ≈ m + end + + for elt in (ComplexF32, ComplexF64) + m = ones(elt, (1, 1)) + ft = to_fusiontensor(m, (g,), (gb,)) + @test eltype(ft) === ComplexF64 + @test Array(ft) ≈ m + end end - for elt in (ComplexF32, ComplexF64) - m = ones(elt, (1, 1)) - ft = to_fusiontensor(m, (g,), (gb,)) - @test eltype(ft) === ComplexF64 - @test Array(ft) ≈ m + @testset "several axes, one block" begin + g1 = gradedrange([TrivialSector() => 2]) + g2 = gradedrange([TrivialSector() => 3]) + g3 = gradedrange([TrivialSector() => 4]) + g4 = gradedrange([TrivialSector() => 2]) + codomain_legs = (g1, g2) + domain_legs = dual.((g3, g4)) + t = convert.(Float64, reshape(collect(1:48), (2, 3, 4, 2))) + ft = to_fusiontensor(t, codomain_legs, domain_legs) + @test size(data_matrix(ft)) == (6, 8) + @test blocksize(data_matrix(ft)) == (1, 1) + @test data_matrix(ft)[Block(1, 1)] ≈ reshape(t, (6, 8)) + @test isnothing(check_sanity(ft)) + @test Array(ft) ≈ t + @test Array(adjoint(ft)) ≈ permutedims(t, (3, 4, 1, 2)) end - end - - @testset "several axes, one block" begin - g1 = gradedrange([TrivialSector() => 2]) - g2 = gradedrange([TrivialSector() => 3]) - g3 = gradedrange([TrivialSector() => 4]) - g4 = gradedrange([TrivialSector() => 2]) - codomain_legs = (g1, g2) - domain_legs = dual.((g3, g4)) - t = convert.(Float64, reshape(collect(1:48), (2, 3, 4, 2))) - ft = to_fusiontensor(t, codomain_legs, domain_legs) - @test size(data_matrix(ft)) == (6, 8) - @test blocksize(data_matrix(ft)) == (1, 1) - @test data_matrix(ft)[Block(1, 1)] ≈ reshape(t, (6, 8)) - @test isnothing(check_sanity(ft)) - @test Array(ft) ≈ t - @test Array(adjoint(ft)) ≈ permutedims(t, (3, 4, 1, 2)) - end end @testset "Abelian FusionTensor" begin - @testset "trivial matrix" begin - g = gradedrange([U1(0) => 1]) - gb = dual(g) - m = ones((1, 1)) - ft = to_fusiontensor(m, (g,), (gb,)) - @test size(data_matrix(ft)) == (1, 1) - @test blocksize(data_matrix(ft)) == (1, 1) - @test data_matrix(ft)[1, 1] ≈ 1.0 - @test isnothing(check_sanity(ft)) - @test Array(ft) ≈ m - @test Array(adjoint(ft)) ≈ m - end - - @testset "non self-conjugate matrix" begin - g = gradedrange([U1(1) => 2]) - gb = dual(g) - m = ones((2, 2)) - ft = to_fusiontensor(m, (g,), (gb,)) - @test size(data_matrix(ft)) == (2, 2) - @test blocksize(data_matrix(ft)) == (1, 1) - @test data_matrix(ft)[Block(1, 1)] ≈ m - @test isnothing(check_sanity(ft)) - @test Array(ft) ≈ m - @test Array(adjoint(ft)) ≈ m - end - - @testset "2-block identity" begin - g = gradedrange([U1(1) => 1, U1(2) => 2]) - codomain_legs = (g,) - domain_legs = (dual(g),) - dense = Array{Float64}(LinearAlgebra.I(3)) - ft = to_fusiontensor(dense, codomain_legs, domain_legs) - @test size(data_matrix(ft)) == (3, 3) - @test blocksize(data_matrix(ft)) == (2, 2) - @test data_matrix(ft)[Block(1, 1)] ≈ ones((1, 1)) - @test data_matrix(ft)[Block(2, 2)] ≈ LinearAlgebra.I(2) - @test data_matrix(ft) ≈ dense - @test isnothing(check_sanity(ft)) - @test Array(ft) ≈ dense - @test Array(adjoint(ft)) ≈ adjoint(dense) - - @test_throws BoundsError to_fusiontensor( - dense, (gradedrange([U1(1) => 1, U1(2) => 3]),), domain_legs - ) - @test_throws MethodError to_fusiontensor(dense, (g, g), domain_legs) - - ba = BlockedArray(dense, [1, 2], [1, 2]) - @test_throws DomainError to_fusiontensor( - ba, (gradedrange([U1(1) => 1, U1(2) => 3]),), domain_legs - ) - @test_throws DomainError to_fusiontensor(ba, (g, g), domain_legs) - dense[1, 2] = 1 # forbidden - @test_throws InexactError to_fusiontensor(dense, codomain_legs, domain_legs) - end - - @testset "several axes, one block" begin - g1 = gradedrange([U1(1) => 2]) - g2 = gradedrange([U1(2) => 3]) - g3 = gradedrange([U1(3) => 4]) - g4 = gradedrange([U1(0) => 2]) - codomain_legs = (g1, g2) - domain_legs = dual.((g3, g4)) - t = convert.(Float64, reshape(collect(1:48), (2, 3, 4, 2))) - ft = to_fusiontensor(t, codomain_legs, domain_legs) - @test size(data_matrix(ft)) == (6, 8) - @test blocksize(data_matrix(ft)) == (1, 1) - @test data_matrix(ft)[Block(1, 1)] ≈ reshape(t, (6, 8)) - @test isnothing(check_sanity(ft)) - @test Array(ft) ≈ t - @test Array(adjoint(ft)) ≈ permutedims(t, (3, 4, 1, 2)) - end - - @testset "several axes, several blocks" begin - g1 = gradedrange([U1(1) => 2, U1(2) => 2]) - g2 = gradedrange([U1(2) => 3, U1(3) => 2]) - g3 = gradedrange([U1(3) => 4, U1(4) => 1]) - g4 = gradedrange([U1(0) => 2, U1(2) => 1]) - codomain_legs = (g1, g2) - domain_legs = dual.((g3, g4)) - dense = zeros((4, 5, 5, 3)) - dense[1:2, 1:3, 1:4, 1:2] .= 1.0 - dense[3:4, 1:3, 5:5, 1:2] .= 2.0 - dense[1:2, 4:5, 5:5, 1:2] .= 3.0 - dense[3:4, 4:5, 1:4, 3:3] .= 4.0 - ft = to_fusiontensor(dense, codomain_legs, domain_legs) - @test size(data_matrix(ft)) == (20, 15) - @test blocksize(data_matrix(ft)) == (3, 4) - @test norm(ft) ≈ norm(dense) - @test isnothing(check_sanity(ft)) - @test Array(ft) ≈ dense - @test Array(adjoint(ft)) ≈ permutedims(dense, (3, 4, 1, 2)) - end - - @testset "mixing dual and nondual" begin - g1 = gradedrange([U1(-1) => 1, U1(0) => 1, U1(1) => 1]) - g2 = gradedrange([U1(0) => 1, U1(1) => 2, U1(2) => 3]) - g3 = gradedrange([U1(0) => 2, U1(1) => 2, U1(3) => 1]) - g4 = gradedrange([U1(-1) => 1, U1(0) => 2, U1(1) => 1]) - codomain_legs = (g1,) - domain_legs = (dual(g2), dual(g3), g4) - dense = zeros(ComplexF64, (3, 6, 5, 4)) - dense[2:2, 1:1, 1:2, 2:3] .= 1.0im - ft = to_fusiontensor(dense, codomain_legs, domain_legs) - @test size(data_matrix(ft)) == (3, 120) - @test blocksize(data_matrix(ft)) == (3, 8) - @test norm(ft) ≈ norm(dense) - @test isnothing(check_sanity(ft)) - @test Array(ft) ≈ dense - @test Array(adjoint(ft)) ≈ conj(permutedims(dense, (2, 3, 4, 1))) - end - - @testset "Less than 2 axes" begin - g = gradedrange([U1(0) => 1, U1(1) => 2, U1(2) => 3]) - v = zeros((6,)) - v[1] = 1.0 - - ft1 = to_fusiontensor(v, (g,), ()) - @test isnothing(check_sanity(ft1)) - @test ndims(ft1) == 1 - @test vec(Array(data_matrix(ft1))) ≈ v - @test Array(ft1) ≈ v - @test Array(adjoint(ft1)) ≈ v - - ft2 = to_fusiontensor(v, (), (dual(g),)) - @test isnothing(check_sanity(ft2)) - @test ndims(ft2) == 1 - @test vec(Array(data_matrix(ft2))) ≈ v - @test Array(ft2) ≈ v - @test Array(adjoint(ft2)) ≈ v - - ft3 = to_fusiontensor(v, (dual(g),), ()) - @test isnothing(check_sanity(ft3)) - @test Array(ft3) ≈ v - @test Array(adjoint(ft3)) ≈ v - - ft4 = to_fusiontensor(v, (), (g,)) - @test isnothing(check_sanity(ft4)) - @test Array(ft4) ≈ v - @test Array(adjoint(ft4)) ≈ v - - if VERSION >= v"1.11" # https://github.com/JuliaLang/julia/issues/52615 - zerodim = ones(()) - ft = to_fusiontensor(zerodim, (), ()) - @test ft isa FusionTensor - @test ndims(ft) == 0 - @test isnothing(check_sanity(ft)) - @test size(data_matrix(ft)) == (1, 1) - @test data_matrix(ft)[1, 1] ≈ 1.0 - @test Array(ft) ≈ zerodim - @test ndims(Array(ft)) == 0 + @testset "trivial matrix" begin + g = gradedrange([U1(0) => 1]) + gb = dual(g) + m = ones((1, 1)) + ft = to_fusiontensor(m, (g,), (gb,)) + @test size(data_matrix(ft)) == (1, 1) + @test blocksize(data_matrix(ft)) == (1, 1) + @test data_matrix(ft)[1, 1] ≈ 1.0 + @test isnothing(check_sanity(ft)) + @test Array(ft) ≈ m + @test Array(adjoint(ft)) ≈ m + end + + @testset "non self-conjugate matrix" begin + g = gradedrange([U1(1) => 2]) + gb = dual(g) + m = ones((2, 2)) + ft = to_fusiontensor(m, (g,), (gb,)) + @test size(data_matrix(ft)) == (2, 2) + @test blocksize(data_matrix(ft)) == (1, 1) + @test data_matrix(ft)[Block(1, 1)] ≈ m + @test isnothing(check_sanity(ft)) + @test Array(ft) ≈ m + @test Array(adjoint(ft)) ≈ m + end + + @testset "2-block identity" begin + g = gradedrange([U1(1) => 1, U1(2) => 2]) + codomain_legs = (g,) + domain_legs = (dual(g),) + dense = Array{Float64}(LinearAlgebra.I(3)) + ft = to_fusiontensor(dense, codomain_legs, domain_legs) + @test size(data_matrix(ft)) == (3, 3) + @test blocksize(data_matrix(ft)) == (2, 2) + @test data_matrix(ft)[Block(1, 1)] ≈ ones((1, 1)) + @test data_matrix(ft)[Block(2, 2)] ≈ LinearAlgebra.I(2) + @test data_matrix(ft) ≈ dense + @test isnothing(check_sanity(ft)) + @test Array(ft) ≈ dense + @test Array(adjoint(ft)) ≈ adjoint(dense) + + @test_throws BoundsError to_fusiontensor( + dense, (gradedrange([U1(1) => 1, U1(2) => 3]),), domain_legs + ) + @test_throws MethodError to_fusiontensor(dense, (g, g), domain_legs) + + ba = BlockedArray(dense, [1, 2], [1, 2]) + @test_throws DomainError to_fusiontensor( + ba, (gradedrange([U1(1) => 1, U1(2) => 3]),), domain_legs + ) + @test_throws DomainError to_fusiontensor(ba, (g, g), domain_legs) + dense[1, 2] = 1 # forbidden + @test_throws InexactError to_fusiontensor(dense, codomain_legs, domain_legs) + end + + @testset "several axes, one block" begin + g1 = gradedrange([U1(1) => 2]) + g2 = gradedrange([U1(2) => 3]) + g3 = gradedrange([U1(3) => 4]) + g4 = gradedrange([U1(0) => 2]) + codomain_legs = (g1, g2) + domain_legs = dual.((g3, g4)) + t = convert.(Float64, reshape(collect(1:48), (2, 3, 4, 2))) + ft = to_fusiontensor(t, codomain_legs, domain_legs) + @test size(data_matrix(ft)) == (6, 8) + @test blocksize(data_matrix(ft)) == (1, 1) + @test data_matrix(ft)[Block(1, 1)] ≈ reshape(t, (6, 8)) + @test isnothing(check_sanity(ft)) + @test Array(ft) ≈ t + @test Array(adjoint(ft)) ≈ permutedims(t, (3, 4, 1, 2)) + end + + @testset "several axes, several blocks" begin + g1 = gradedrange([U1(1) => 2, U1(2) => 2]) + g2 = gradedrange([U1(2) => 3, U1(3) => 2]) + g3 = gradedrange([U1(3) => 4, U1(4) => 1]) + g4 = gradedrange([U1(0) => 2, U1(2) => 1]) + codomain_legs = (g1, g2) + domain_legs = dual.((g3, g4)) + dense = zeros((4, 5, 5, 3)) + dense[1:2, 1:3, 1:4, 1:2] .= 1.0 + dense[3:4, 1:3, 5:5, 1:2] .= 2.0 + dense[1:2, 4:5, 5:5, 1:2] .= 3.0 + dense[3:4, 4:5, 1:4, 3:3] .= 4.0 + ft = to_fusiontensor(dense, codomain_legs, domain_legs) + @test size(data_matrix(ft)) == (20, 15) + @test blocksize(data_matrix(ft)) == (3, 4) + @test norm(ft) ≈ norm(dense) + @test isnothing(check_sanity(ft)) + @test Array(ft) ≈ dense + @test Array(adjoint(ft)) ≈ permutedims(dense, (3, 4, 1, 2)) + end + + @testset "mixing dual and nondual" begin + g1 = gradedrange([U1(-1) => 1, U1(0) => 1, U1(1) => 1]) + g2 = gradedrange([U1(0) => 1, U1(1) => 2, U1(2) => 3]) + g3 = gradedrange([U1(0) => 2, U1(1) => 2, U1(3) => 1]) + g4 = gradedrange([U1(-1) => 1, U1(0) => 2, U1(1) => 1]) + codomain_legs = (g1,) + domain_legs = (dual(g2), dual(g3), g4) + dense = zeros(ComplexF64, (3, 6, 5, 4)) + dense[2:2, 1:1, 1:2, 2:3] .= 1.0im + ft = to_fusiontensor(dense, codomain_legs, domain_legs) + @test size(data_matrix(ft)) == (3, 120) + @test blocksize(data_matrix(ft)) == (3, 8) + @test norm(ft) ≈ norm(dense) + @test isnothing(check_sanity(ft)) + @test Array(ft) ≈ dense + @test Array(adjoint(ft)) ≈ conj(permutedims(dense, (2, 3, 4, 1))) + end + + @testset "Less than 2 axes" begin + g = gradedrange([U1(0) => 1, U1(1) => 2, U1(2) => 3]) + v = zeros((6,)) + v[1] = 1.0 + + ft1 = to_fusiontensor(v, (g,), ()) + @test isnothing(check_sanity(ft1)) + @test ndims(ft1) == 1 + @test vec(Array(data_matrix(ft1))) ≈ v + @test Array(ft1) ≈ v + @test Array(adjoint(ft1)) ≈ v + + ft2 = to_fusiontensor(v, (), (dual(g),)) + @test isnothing(check_sanity(ft2)) + @test ndims(ft2) == 1 + @test vec(Array(data_matrix(ft2))) ≈ v + @test Array(ft2) ≈ v + @test Array(adjoint(ft2)) ≈ v + + ft3 = to_fusiontensor(v, (dual(g),), ()) + @test isnothing(check_sanity(ft3)) + @test Array(ft3) ≈ v + @test Array(adjoint(ft3)) ≈ v + + ft4 = to_fusiontensor(v, (), (g,)) + @test isnothing(check_sanity(ft4)) + @test Array(ft4) ≈ v + @test Array(adjoint(ft4)) ≈ v + + if VERSION >= v"1.11" # https://github.com/JuliaLang/julia/issues/52615 + zerodim = ones(()) + ft = to_fusiontensor(zerodim, (), ()) + @test ft isa FusionTensor + @test ndims(ft) == 0 + @test isnothing(check_sanity(ft)) + @test size(data_matrix(ft)) == (1, 1) + @test data_matrix(ft)[1, 1] ≈ 1.0 + @test Array(ft) ≈ zerodim + @test ndims(Array(ft)) == 0 + end end - end end @testset "O(2) FusionTensor" begin - @testset "trivial matrix" begin - g = gradedrange([O2(0) => 1]) - gb = dual(g) - m = ones((1, 1)) - ft = to_fusiontensor(m, (gb,), (g,)) - @test size(data_matrix(ft)) == (1, 1) - @test blocksize(data_matrix(ft)) == (1, 1) - @test data_matrix(ft)[1, 1] ≈ 1.0 - @test isnothing(check_sanity(ft)) - @test Array(ft) ≈ m - @test Array(adjoint(ft)) ≈ m - end - - @testset "spin 1/2 S.S" begin - g2 = gradedrange([O2(1//2) => 1]) - g2b = dual(g2) - - # identity - id2 = LinearAlgebra.I((2)) - ft = to_fusiontensor(id2, (g2,), (g2b,)) - @test norm(ft) ≈ √2 - @test isnothing(check_sanity(ft)) - @test Array(ft) ≈ id2 - @test Array(adjoint(ft)) ≈ id2 - - # S⋅S - sds22 = reshape( - [ - [0.25, 0.0, 0.0, 0.0] - [0.0, -0.25, 0.5, 0.0] - [0.0, 0.5, -0.25, 0.0] - [0.0, 0.0, 0.0, 0.25] - ], - (2, 2, 2, 2), - ) - dense, codomain_legs, domain_legs = sds22, (g2, g2), (g2b, g2b) - ft = to_fusiontensor(dense, codomain_legs, domain_legs) - @test norm(ft) ≈ √3 / 2 - @test isnothing(check_sanity(ft)) - @test Array(ft) ≈ sds22 - @test Array(adjoint(ft)) ≈ sds22 - - # dual over one spin. This changes the dense coefficients but not the FusionTensor ones - sds22b = reshape( - [ - [-0.25, 0.0, 0.0, -0.5] - [0.0, 0.25, 0.0, 0.0] - [0.0, 0.0, 0.25, 0.0] - [-0.5, 0.0, 0.0, -0.25] - ], - (2, 2, 2, 2), - ) - sds22b_codomain_legs = (g2, g2b) - dense, codomain_legs, domain_legs = sds22b, (g2, g2b), (g2b, g2) - ftb = to_fusiontensor(dense, codomain_legs, domain_legs) - @test norm(ftb) ≈ √3 / 2 - @test isnothing(check_sanity(ft)) - @test Array(ftb) ≈ sds22b - @test Array(adjoint(ftb)) ≈ sds22b - - # no domain axis - dense, codomain_legs, domain_legs = sds22, (g2, g2, g2b, g2b), () - ft = to_fusiontensor(dense, codomain_legs, domain_legs) - @test isnothing(check_sanity(ft)) - @test Array(ft) ≈ sds22 - @test Array(adjoint(ft)) ≈ sds22 - - # no codomain axis - dense, codomain_legs, domain_legs = sds22, (), (g2, g2, g2b, g2b) - ft = to_fusiontensor(dense, codomain_legs, domain_legs) - @test isnothing(check_sanity(ft)) - @test Array(ft) ≈ sds22 - @test Array(adjoint(ft)) ≈ sds22 - end + @testset "trivial matrix" begin + g = gradedrange([O2(0) => 1]) + gb = dual(g) + m = ones((1, 1)) + ft = to_fusiontensor(m, (gb,), (g,)) + @test size(data_matrix(ft)) == (1, 1) + @test blocksize(data_matrix(ft)) == (1, 1) + @test data_matrix(ft)[1, 1] ≈ 1.0 + @test isnothing(check_sanity(ft)) + @test Array(ft) ≈ m + @test Array(adjoint(ft)) ≈ m + end + + @testset "spin 1/2 S.S" begin + g2 = gradedrange([O2(1 // 2) => 1]) + g2b = dual(g2) + + # identity + id2 = LinearAlgebra.I((2)) + ft = to_fusiontensor(id2, (g2,), (g2b,)) + @test norm(ft) ≈ √2 + @test isnothing(check_sanity(ft)) + @test Array(ft) ≈ id2 + @test Array(adjoint(ft)) ≈ id2 + + # S⋅S + sds22 = reshape( + [ + [0.25, 0.0, 0.0, 0.0] + [0.0, -0.25, 0.5, 0.0] + [0.0, 0.5, -0.25, 0.0] + [0.0, 0.0, 0.0, 0.25] + ], + (2, 2, 2, 2), + ) + dense, codomain_legs, domain_legs = sds22, (g2, g2), (g2b, g2b) + ft = to_fusiontensor(dense, codomain_legs, domain_legs) + @test norm(ft) ≈ √3 / 2 + @test isnothing(check_sanity(ft)) + @test Array(ft) ≈ sds22 + @test Array(adjoint(ft)) ≈ sds22 + + # dual over one spin. This changes the dense coefficients but not the FusionTensor ones + sds22b = reshape( + [ + [-0.25, 0.0, 0.0, -0.5] + [0.0, 0.25, 0.0, 0.0] + [0.0, 0.0, 0.25, 0.0] + [-0.5, 0.0, 0.0, -0.25] + ], + (2, 2, 2, 2), + ) + sds22b_codomain_legs = (g2, g2b) + dense, codomain_legs, domain_legs = sds22b, (g2, g2b), (g2b, g2) + ftb = to_fusiontensor(dense, codomain_legs, domain_legs) + @test norm(ftb) ≈ √3 / 2 + @test isnothing(check_sanity(ft)) + @test Array(ftb) ≈ sds22b + @test Array(adjoint(ftb)) ≈ sds22b + + # no domain axis + dense, codomain_legs, domain_legs = sds22, (g2, g2, g2b, g2b), () + ft = to_fusiontensor(dense, codomain_legs, domain_legs) + @test isnothing(check_sanity(ft)) + @test Array(ft) ≈ sds22 + @test Array(adjoint(ft)) ≈ sds22 + + # no codomain axis + dense, codomain_legs, domain_legs = sds22, (), (g2, g2, g2b, g2b) + ft = to_fusiontensor(dense, codomain_legs, domain_legs) + @test isnothing(check_sanity(ft)) + @test Array(ft) ≈ sds22 + @test Array(adjoint(ft)) ≈ sds22 + end end @testset "SU(2) FusionTensor" begin - @testset "trivial matrix" begin - g = gradedrange([SU2(0) => 1]) - gb = dual(g) - m = ones((1, 1)) - ft = to_fusiontensor(m, (gb,), (g,)) - @test size(data_matrix(ft)) == (1, 1) - @test blocksize(data_matrix(ft)) == (1, 1) - @test data_matrix(ft)[1, 1] ≈ 1.0 - @test isnothing(check_sanity(ft)) - @test Array(ft) ≈ m - @test Array(adjoint(ft)) ≈ m - end - - @testset "spin 1/2 S.S" begin - g2 = gradedrange([SU2(1 / 2) => 1]) - g2b = dual(g2) - - # identity - id2 = LinearAlgebra.I((2)) - ft = to_fusiontensor(id2, (g2,), (g2b,)) - @test norm(ft) ≈ √2 - @test isnothing(check_sanity(ft)) - @test Array(ft) ≈ id2 - @test Array(adjoint(ft)) ≈ id2 - - # S⋅S - sds22 = reshape( - [ - [0.25, 0.0, 0.0, 0.0] - [0.0, -0.25, 0.5, 0.0] - [0.0, 0.5, -0.25, 0.0] - [0.0, 0.0, 0.0, 0.25] - ], - (2, 2, 2, 2), - ) - dense, codomain_legs, domain_legs = sds22, (g2, g2), (g2b, g2b) - ft = to_fusiontensor(dense, codomain_legs, domain_legs) - @test norm(ft) ≈ √3 / 2 - @test isnothing(check_sanity(ft)) - @test Array(ft) ≈ sds22 - @test Array(adjoint(ft)) ≈ sds22 - - # dual over one spin. This changes the dense coefficients but not the FusionTensor ones - sds22b = reshape( - [ - [-0.25, 0.0, 0.0, -0.5] - [0.0, 0.25, 0.0, 0.0] - [0.0, 0.0, 0.25, 0.0] - [-0.5, 0.0, 0.0, -0.25] - ], - (2, 2, 2, 2), - ) - sds22b_codomain_legs = (g2, g2b) - dense, codomain_legs, domain_legs = sds22b, (g2, g2b), (g2b, g2) - ftb = to_fusiontensor(dense, codomain_legs, domain_legs) - @test norm(ftb) ≈ √3 / 2 - @test isnothing(check_sanity(ft)) - @test Array(ftb) ≈ sds22b - @test Array(adjoint(ftb)) ≈ sds22b - - # no domain axis - dense, codomain_legs, domain_legs = sds22, (g2b, g2b, g2, g2), () - ft = to_fusiontensor(dense, codomain_legs, domain_legs) - @test isnothing(check_sanity(ft)) - @test Array(ft) ≈ sds22 - @test Array(adjoint(ft)) ≈ sds22 - - # no codomain axis - dense, codomain_legs, domain_legs = sds22, (), (g2b, g2b, g2, g2) - ft = to_fusiontensor(dense, codomain_legs, domain_legs) - @test isnothing(check_sanity(ft)) - @test Array(ft) ≈ sds22 - @test Array(adjoint(ft)) ≈ sds22 - end - - @testset "large identity" begin - g = reduce(tensor_product, (SU2(1 / 2), SU2(1 / 2), SU2(1 / 2))) - N = 3 - codomain_legs = ntuple(_ -> g, N) - domain_legs = dual.(codomain_legs) - d = 8 - dense = reshape(LinearAlgebra.I(d^N), ntuple(_ -> d, 2 * N)) - ft = to_fusiontensor(dense, codomain_legs, domain_legs) - @test isnothing(check_sanity(ft)) - @test Array(ft) ≈ dense - @test Array(adjoint(ft)) ≈ dense - end + @testset "trivial matrix" begin + g = gradedrange([SU2(0) => 1]) + gb = dual(g) + m = ones((1, 1)) + ft = to_fusiontensor(m, (gb,), (g,)) + @test size(data_matrix(ft)) == (1, 1) + @test blocksize(data_matrix(ft)) == (1, 1) + @test data_matrix(ft)[1, 1] ≈ 1.0 + @test isnothing(check_sanity(ft)) + @test Array(ft) ≈ m + @test Array(adjoint(ft)) ≈ m + end + + @testset "spin 1/2 S.S" begin + g2 = gradedrange([SU2(1 / 2) => 1]) + g2b = dual(g2) + + # identity + id2 = LinearAlgebra.I((2)) + ft = to_fusiontensor(id2, (g2,), (g2b,)) + @test norm(ft) ≈ √2 + @test isnothing(check_sanity(ft)) + @test Array(ft) ≈ id2 + @test Array(adjoint(ft)) ≈ id2 + + # S⋅S + sds22 = reshape( + [ + [0.25, 0.0, 0.0, 0.0] + [0.0, -0.25, 0.5, 0.0] + [0.0, 0.5, -0.25, 0.0] + [0.0, 0.0, 0.0, 0.25] + ], + (2, 2, 2, 2), + ) + dense, codomain_legs, domain_legs = sds22, (g2, g2), (g2b, g2b) + ft = to_fusiontensor(dense, codomain_legs, domain_legs) + @test norm(ft) ≈ √3 / 2 + @test isnothing(check_sanity(ft)) + @test Array(ft) ≈ sds22 + @test Array(adjoint(ft)) ≈ sds22 + + # dual over one spin. This changes the dense coefficients but not the FusionTensor ones + sds22b = reshape( + [ + [-0.25, 0.0, 0.0, -0.5] + [0.0, 0.25, 0.0, 0.0] + [0.0, 0.0, 0.25, 0.0] + [-0.5, 0.0, 0.0, -0.25] + ], + (2, 2, 2, 2), + ) + sds22b_codomain_legs = (g2, g2b) + dense, codomain_legs, domain_legs = sds22b, (g2, g2b), (g2b, g2) + ftb = to_fusiontensor(dense, codomain_legs, domain_legs) + @test norm(ftb) ≈ √3 / 2 + @test isnothing(check_sanity(ft)) + @test Array(ftb) ≈ sds22b + @test Array(adjoint(ftb)) ≈ sds22b + + # no domain axis + dense, codomain_legs, domain_legs = sds22, (g2b, g2b, g2, g2), () + ft = to_fusiontensor(dense, codomain_legs, domain_legs) + @test isnothing(check_sanity(ft)) + @test Array(ft) ≈ sds22 + @test Array(adjoint(ft)) ≈ sds22 + + # no codomain axis + dense, codomain_legs, domain_legs = sds22, (), (g2b, g2b, g2, g2) + ft = to_fusiontensor(dense, codomain_legs, domain_legs) + @test isnothing(check_sanity(ft)) + @test Array(ft) ≈ sds22 + @test Array(adjoint(ft)) ≈ sds22 + end + + @testset "large identity" begin + g = reduce(tensor_product, (SU2(1 / 2), SU2(1 / 2), SU2(1 / 2))) + N = 3 + codomain_legs = ntuple(_ -> g, N) + domain_legs = dual.(codomain_legs) + d = 8 + dense = reshape(LinearAlgebra.I(d^N), ntuple(_ -> d, 2 * N)) + ft = to_fusiontensor(dense, codomain_legs, domain_legs) + @test isnothing(check_sanity(ft)) + @test Array(ft) ≈ dense + @test Array(adjoint(ft)) ≈ dense + end end @testset "U(1)×SU(2) FusionTensor" begin - for d in 1:6 # any spin dimension - s = SU2((d - 1)//2) # d = 2s+1 - D = d + 1 - tRVB = zeros((d, D, D, D, D)) # tensor RVB SU(2) for spin s - for i in 1:d - tRVB[i, i + 1, 1, 1, 1] = 1.0 - tRVB[i, 1, i + 1, 1, 1] = 1.0 - tRVB[i, 1, 1, i + 1, 1] = 1.0 - tRVB[i, 1, 1, 1, i + 1] = 1.0 + for d in 1:6 # any spin dimension + s = SU2((d - 1) // 2) # d = 2s+1 + D = d + 1 + tRVB = zeros((d, D, D, D, D)) # tensor RVB SU(2) for spin s + for i in 1:d + tRVB[i, i + 1, 1, 1, 1] = 1.0 + tRVB[i, 1, i + 1, 1, 1] = 1.0 + tRVB[i, 1, 1, i + 1, 1] = 1.0 + tRVB[i, 1, 1, 1, i + 1] = 1.0 + end + + gd = gradedrange([SectorProduct(s, U1(3)) => 1]) + codomain_legs = (dual(gd),) + gD = gradedrange([SectorProduct(SU2(0), U1(1)) => 1, SectorProduct(s, U1(0)) => 1]) + domain_legs = (gD, gD, gD, gD) + ft = to_fusiontensor(tRVB, codomain_legs, domain_legs) + @test isnothing(check_sanity(ft)) + @test Array(ft) ≈ tRVB + + # same with NamedTuples + gd_nt = gradedrange([SectorProduct(; S = s, N = U1(3)) => 1]) + codomain_legs_nt = (dual(gd_nt),) + gD_nt = gradedrange( + [ + SectorProduct(; S = SU2(0), N = U1(1)) => 1, SectorProduct(; S = s, N = U1(0)) => 1, + ] + ) + domain_legs_nt = (gD_nt, gD_nt, gD_nt, gD_nt) + ft_nt = to_fusiontensor(tRVB, codomain_legs_nt, domain_legs_nt) + @test isnothing(check_sanity(ft_nt)) + @test Array(ft_nt) ≈ tRVB end - - gd = gradedrange([SectorProduct(s, U1(3)) => 1]) - codomain_legs = (dual(gd),) - gD = gradedrange([SectorProduct(SU2(0), U1(1)) => 1, SectorProduct(s, U1(0)) => 1]) - domain_legs = (gD, gD, gD, gD) - ft = to_fusiontensor(tRVB, codomain_legs, domain_legs) - @test isnothing(check_sanity(ft)) - @test Array(ft) ≈ tRVB - - # same with NamedTuples - gd_nt = gradedrange([SectorProduct(; S=s, N=U1(3)) => 1]) - codomain_legs_nt = (dual(gd_nt),) - gD_nt = gradedrange([ - SectorProduct(; S=SU2(0), N=U1(1)) => 1, SectorProduct(; S=s, N=U1(0)) => 1 - ]) - domain_legs_nt = (gD_nt, gD_nt, gD_nt, gD_nt) - ft_nt = to_fusiontensor(tRVB, codomain_legs_nt, domain_legs_nt) - @test isnothing(check_sanity(ft_nt)) - @test Array(ft_nt) ≈ tRVB - end end diff --git a/test/test_basics.jl b/test/test_basics.jl index 44476ef..be017ce 100644 --- a/test/test_basics.jl +++ b/test/test_basics.jl @@ -3,30 +3,30 @@ using Test: @test, @test_throws, @testset using BlockArrays: Block using BlockSparseArrays: BlockSparseArray, eachblockstoredindex using FusionTensors: - FusionMatrix, - FusionTensor, - FusionTensorAxes, - codomain_axes, - data_matrix, - domain_axes, - FusionTensor, - codomain_axis, - domain_axis, - ndims_domain, - ndims_codomain + FusionMatrix, + FusionTensor, + FusionTensorAxes, + codomain_axes, + data_matrix, + domain_axes, + FusionTensor, + codomain_axis, + domain_axis, + ndims_domain, + ndims_codomain using GradedArrays: - U1, - SU2, - SectorProduct, - TrivialSector, - Z, - checkspaces, - checkspaces_dual, - dual, - flip, - gradedrange, - sector_type, - space_isequal + U1, + SU2, + SectorProduct, + TrivialSector, + Z, + checkspaces, + checkspaces_dual, + dual, + flip, + gradedrange, + sector_type, + space_isequal using TensorAlgebra: tuplemortar using TensorProducts: tensor_product using LinearAlgebra: LinearAlgebra, norm @@ -35,342 +35,344 @@ using Random: Random include("setup.jl") @testset "Fusion matrix" begin - g1 = gradedrange([U1(0) => 1, U1(1) => 2, U1(2) => 3]) - g2 = dual(gradedrange([U1(0) => 2, U1(1) => 2, U1(3) => 1])) - - fta = FusionTensorAxes((g1,), (g2,)) - ft0 = FusionTensor{Float64}(undef, fta) - @test ft0 isa FusionTensor - @test ft0 isa FusionMatrix - @test space_isequal(codomain_axis(ft0), g1) - @test space_isequal(domain_axis(ft0), g2) - - # check dual convention when initializing data_matrix - ft0 = FusionTensor{Float64}(undef, (g1,), (g2,)) - @test ft0 isa FusionTensor - @test space_isequal(codomain_axis(ft0), g1) - @test space_isequal(domain_axis(ft0), g2) - - m = BlockSparseArray{Float64}(undef, g1, g2) - ft1 = FusionTensor(m, (g1,), (g2,)) - - # getters - @test data_matrix(ft1) == m - @test axes(ft1) == FusionTensorAxes((g1,), (g2,)) - - # misc - @test checkspaces(codomain_axes(ft1), (g1,)) - @test checkspaces(domain_axes(ft1), (g2,)) - @test ndims_codomain(ft1) == 1 - @test ndims_domain(ft1) == 1 - @test size(data_matrix(ft1)) == (6, 5) - @test space_isequal(codomain_axis(ft1), g1) - @test space_isequal(domain_axis(ft1), g2) - @test isnothing(check_sanity(ft0)) - @test isnothing(check_sanity(ft1)) - @test sector_type(ft1) === U1{Int} - @test sector_type(typeof(ft1)) === U1{Int} - - m1 = BlockSparseArray{Float64}(undef, g1, g2) - m1[Block(2, 1)] = ones(2, 2) # forbidden - @test_throws ArgumentError FusionTensor(m1, (g1,), (g2,)) - - # Base methods - @test eltype(ft1) === Float64 - @test length(ft1) == 30 - @test ndims(ft1) == 2 - @test size(ft1) == tuplemortar(((6,), (5,))) - @test_throws MethodError eachindex(ft1) - - # copy - ft2 = copy(ft1) - @test isnothing(check_sanity(ft2)) - @test ft2 !== ft1 - @test data_matrix(ft2) == data_matrix(ft1) - @test data_matrix(ft2) !== data_matrix(ft1) - @test checkspaces(codomain_axes(ft2), codomain_axes(ft1)) - @test checkspaces(domain_axes(ft2), domain_axes(ft1)) - - ft2 = deepcopy(ft1) - @test ft2 !== ft1 - @test data_matrix(ft2) == data_matrix(ft1) - @test data_matrix(ft2) !== data_matrix(ft1) - @test checkspaces(codomain_axes(ft2), codomain_axes(ft1)) - @test checkspaces(domain_axes(ft2), domain_axes(ft1)) - - # similar - ft2 = similar(ft1) - @test isnothing(check_sanity(ft2)) - @test eltype(ft2) == Float64 - @test checkspaces(codomain_axes(ft2), codomain_axes(ft1)) - @test checkspaces(domain_axes(ft2), domain_axes(ft1)) - - ft3 = similar(ft1, ComplexF64) - @test isnothing(check_sanity(ft3)) - @test eltype(ft3) == ComplexF64 - @test checkspaces(codomain_axes(ft3), codomain_axes(ft1)) - @test checkspaces(domain_axes(ft3), domain_axes(ft1)) - - @test_throws AssertionError similar(ft1, Int) - - ft5 = similar(ft1, ComplexF32, ((g1, g1), (g2,))) - @test isnothing(check_sanity(ft5)) - @test eltype(ft5) == ComplexF64 - @test checkspaces(codomain_axes(ft5), (g1, g1)) - @test checkspaces(domain_axes(ft5), (g2,)) - - ft5 = similar(ft1, ComplexF32, tuplemortar(((g1, g1), (g2,)))) - @test isnothing(check_sanity(ft5)) - @test eltype(ft5) == ComplexF64 - @test checkspaces(codomain_axes(ft5), (g1, g1)) - @test checkspaces(domain_axes(ft5), (g2,)) + g1 = gradedrange([U1(0) => 1, U1(1) => 2, U1(2) => 3]) + g2 = dual(gradedrange([U1(0) => 2, U1(1) => 2, U1(3) => 1])) + + fta = FusionTensorAxes((g1,), (g2,)) + ft0 = FusionTensor{Float64}(undef, fta) + @test ft0 isa FusionTensor + @test ft0 isa FusionMatrix + @test space_isequal(codomain_axis(ft0), g1) + @test space_isequal(domain_axis(ft0), g2) + + # check dual convention when initializing data_matrix + ft0 = FusionTensor{Float64}(undef, (g1,), (g2,)) + @test ft0 isa FusionTensor + @test space_isequal(codomain_axis(ft0), g1) + @test space_isequal(domain_axis(ft0), g2) + + m = BlockSparseArray{Float64}(undef, g1, g2) + ft1 = FusionTensor(m, (g1,), (g2,)) + + # getters + @test data_matrix(ft1) == m + @test axes(ft1) == FusionTensorAxes((g1,), (g2,)) + + # misc + @test checkspaces(codomain_axes(ft1), (g1,)) + @test checkspaces(domain_axes(ft1), (g2,)) + @test ndims_codomain(ft1) == 1 + @test ndims_domain(ft1) == 1 + @test size(data_matrix(ft1)) == (6, 5) + @test space_isequal(codomain_axis(ft1), g1) + @test space_isequal(domain_axis(ft1), g2) + @test isnothing(check_sanity(ft0)) + @test isnothing(check_sanity(ft1)) + @test sector_type(ft1) === U1{Int} + @test sector_type(typeof(ft1)) === U1{Int} + + m1 = BlockSparseArray{Float64}(undef, g1, g2) + m1[Block(2, 1)] = ones(2, 2) # forbidden + @test_throws ArgumentError FusionTensor(m1, (g1,), (g2,)) + + # Base methods + @test eltype(ft1) === Float64 + @test length(ft1) == 30 + @test ndims(ft1) == 2 + @test size(ft1) == tuplemortar(((6,), (5,))) + @test_throws MethodError eachindex(ft1) + + # copy + ft2 = copy(ft1) + @test isnothing(check_sanity(ft2)) + @test ft2 !== ft1 + @test data_matrix(ft2) == data_matrix(ft1) + @test data_matrix(ft2) !== data_matrix(ft1) + @test checkspaces(codomain_axes(ft2), codomain_axes(ft1)) + @test checkspaces(domain_axes(ft2), domain_axes(ft1)) + + ft2 = deepcopy(ft1) + @test ft2 !== ft1 + @test data_matrix(ft2) == data_matrix(ft1) + @test data_matrix(ft2) !== data_matrix(ft1) + @test checkspaces(codomain_axes(ft2), codomain_axes(ft1)) + @test checkspaces(domain_axes(ft2), domain_axes(ft1)) + + # similar + ft2 = similar(ft1) + @test isnothing(check_sanity(ft2)) + @test eltype(ft2) == Float64 + @test checkspaces(codomain_axes(ft2), codomain_axes(ft1)) + @test checkspaces(domain_axes(ft2), domain_axes(ft1)) + + ft3 = similar(ft1, ComplexF64) + @test isnothing(check_sanity(ft3)) + @test eltype(ft3) == ComplexF64 + @test checkspaces(codomain_axes(ft3), codomain_axes(ft1)) + @test checkspaces(domain_axes(ft3), domain_axes(ft1)) + + @test_throws AssertionError similar(ft1, Int) + + ft5 = similar(ft1, ComplexF32, ((g1, g1), (g2,))) + @test isnothing(check_sanity(ft5)) + @test eltype(ft5) == ComplexF64 + @test checkspaces(codomain_axes(ft5), (g1, g1)) + @test checkspaces(domain_axes(ft5), (g2,)) + + ft5 = similar(ft1, ComplexF32, tuplemortar(((g1, g1), (g2,)))) + @test isnothing(check_sanity(ft5)) + @test eltype(ft5) == ComplexF64 + @test checkspaces(codomain_axes(ft5), (g1, g1)) + @test checkspaces(domain_axes(ft5), (g2,)) end @testset "More than 2 axes" begin - g1 = gradedrange([U1(0) => 1, U1(1) => 2, U1(2) => 3]) - g2 = gradedrange([U1(0) => 2, U1(1) => 2, U1(3) => 1]) - g3 = dual(gradedrange([U1(-1) => 1, U1(0) => 2, U1(1) => 1])) - g4 = dual(gradedrange([U1(-1) => 1, U1(0) => 1, U1(1) => 1])) - gr = tensor_product(g1, g2) - gc = dual(tensor_product(dual(g3), dual(g4))) - m2 = BlockSparseArray{Float64}(undef, gr, gc) - ft = FusionTensor(m2, (g1, g2), (g3, g4)) - - @test ft isa FusionTensor - @test !(ft isa FusionMatrix) - @test data_matrix(ft) == m2 - @test checkspaces(codomain_axes(ft), (g1, g2)) - @test checkspaces(domain_axes(ft), (g3, g4)) - - @test axes(ft) == FusionTensorAxes(tuplemortar(((g1, g2), (g3, g4)))) - @test ndims_codomain(ft) == 2 - @test ndims_domain(ft) == 2 - @test size(data_matrix(ft)) == (30, 12) - @test space_isequal(codomain_axis(ft), gr) - @test space_isequal(domain_axis(ft), gc) - @test isnothing(check_sanity(ft)) - - @test ndims(ft) == 4 - @test size(ft) == tuplemortar(((6, 5), (4, 3))) + g1 = gradedrange([U1(0) => 1, U1(1) => 2, U1(2) => 3]) + g2 = gradedrange([U1(0) => 2, U1(1) => 2, U1(3) => 1]) + g3 = dual(gradedrange([U1(-1) => 1, U1(0) => 2, U1(1) => 1])) + g4 = dual(gradedrange([U1(-1) => 1, U1(0) => 1, U1(1) => 1])) + gr = tensor_product(g1, g2) + gc = dual(tensor_product(dual(g3), dual(g4))) + m2 = BlockSparseArray{Float64}(undef, gr, gc) + ft = FusionTensor(m2, (g1, g2), (g3, g4)) + + @test ft isa FusionTensor + @test !(ft isa FusionMatrix) + @test data_matrix(ft) == m2 + @test checkspaces(codomain_axes(ft), (g1, g2)) + @test checkspaces(domain_axes(ft), (g3, g4)) + + @test axes(ft) == FusionTensorAxes(tuplemortar(((g1, g2), (g3, g4)))) + @test ndims_codomain(ft) == 2 + @test ndims_domain(ft) == 2 + @test size(data_matrix(ft)) == (30, 12) + @test space_isequal(codomain_axis(ft), gr) + @test space_isequal(domain_axis(ft), gc) + @test isnothing(check_sanity(ft)) + + @test ndims(ft) == 4 + @test size(ft) == tuplemortar(((6, 5), (4, 3))) end @testset "Less than 2 axes" begin - g1 = gradedrange([U1(0) => 1, U1(1) => 2, U1(2) => 3]) - - # one row axis - ft1 = FusionTensor{Float64}(undef, (g1,), ()) - @test ft1 isa FusionTensor - @test !(ft1 isa FusionMatrix) - @test ndims_codomain(ft1) == 1 - @test ndims_domain(ft1) == 0 - @test ndims(ft1) == 1 - @test size(ft1) == tuplemortar(((6,), ())) - @test size(data_matrix(ft1)) == (6, 1) - @test isnothing(check_sanity(ft1)) - @test sector_type(ft1) === sector_type(g1) - - # one column axis - ft2 = FusionTensor{Float64}(undef, (), (g1,)) - @test ft2 isa FusionTensor - @test !(ft2 isa FusionMatrix) - @test ndims_codomain(ft2) == 0 - @test ndims_domain(ft2) == 1 - @test ndims(ft2) == 1 - @test size(ft2) == tuplemortar(((), (6,))) - @test size(data_matrix(ft2)) == (1, 6) - @test isnothing(check_sanity(ft2)) - @test sector_type(ft2) === sector_type(g1) - - # zero axis - ft3 = FusionTensor{Float64}(undef, (), ()) - @test ft3 isa FusionTensor - @test !(ft3 isa FusionMatrix) - @test ndims_codomain(ft3) == 0 - @test ndims_domain(ft3) == 0 - @test ndims(ft3) == 0 - @test size(ft3) == tuplemortar(((), ())) - @test size(data_matrix(ft3)) == (1, 1) - @test isnothing(check_sanity(ft3)) - @test sector_type(ft3) === TrivialSector - - ft4 = FusionTensor{Float64}(undef, (g1, g1), ()) - @test ft4 isa FusionTensor - @test !(ft4 isa FusionMatrix) - ft5 = FusionTensor{Float64}(undef, (), (g1, g1)) - @test ft5 isa FusionTensor - @test !(ft5 isa FusionMatrix) + g1 = gradedrange([U1(0) => 1, U1(1) => 2, U1(2) => 3]) + + # one row axis + ft1 = FusionTensor{Float64}(undef, (g1,), ()) + @test ft1 isa FusionTensor + @test !(ft1 isa FusionMatrix) + @test ndims_codomain(ft1) == 1 + @test ndims_domain(ft1) == 0 + @test ndims(ft1) == 1 + @test size(ft1) == tuplemortar(((6,), ())) + @test size(data_matrix(ft1)) == (6, 1) + @test isnothing(check_sanity(ft1)) + @test sector_type(ft1) === sector_type(g1) + + # one column axis + ft2 = FusionTensor{Float64}(undef, (), (g1,)) + @test ft2 isa FusionTensor + @test !(ft2 isa FusionMatrix) + @test ndims_codomain(ft2) == 0 + @test ndims_domain(ft2) == 1 + @test ndims(ft2) == 1 + @test size(ft2) == tuplemortar(((), (6,))) + @test size(data_matrix(ft2)) == (1, 6) + @test isnothing(check_sanity(ft2)) + @test sector_type(ft2) === sector_type(g1) + + # zero axis + ft3 = FusionTensor{Float64}(undef, (), ()) + @test ft3 isa FusionTensor + @test !(ft3 isa FusionMatrix) + @test ndims_codomain(ft3) == 0 + @test ndims_domain(ft3) == 0 + @test ndims(ft3) == 0 + @test size(ft3) == tuplemortar(((), ())) + @test size(data_matrix(ft3)) == (1, 1) + @test isnothing(check_sanity(ft3)) + @test sector_type(ft3) === TrivialSector + + ft4 = FusionTensor{Float64}(undef, (g1, g1), ()) + @test ft4 isa FusionTensor + @test !(ft4 isa FusionMatrix) + ft5 = FusionTensor{Float64}(undef, (), (g1, g1)) + @test ft5 isa FusionTensor + @test !(ft5 isa FusionMatrix) end @testset "specific constructors" begin - g1 = gradedrange([U1(0) => 1, U1(1) => 2, U1(2) => 3]) - g2 = gradedrange([U1(0) => 2, U1(1) => 2, U1(3) => 1]) - g3 = gradedrange([U1(-1) => 1, U1(0) => 2, U1(1) => 1]) - g4 = gradedrange([U1(-1) => 1, U1(0) => 1, U1(1) => 1]) - - fta = FusionTensorAxes((g1,), (g2, g3)) - @test zeros(fta) isa FusionTensor{Float64,3} - @test zeros(ComplexF64, fta) isa FusionTensor{ComplexF64,3} - - rng = Random.default_rng() - ft1 = randn(rng, ComplexF64, fta) - @test ft1 isa FusionTensor{ComplexF64,3} - @test all(!=(0), data_matrix(ft1)[Block(1, 5)]) - @test randn(rng, fta) isa FusionTensor{Float64,3} - @test randn(ComplexF64, fta) isa FusionTensor{ComplexF64,3} - @test randn(fta) isa FusionTensor{Float64,3} - - ft2 = FusionTensor(LinearAlgebra.I, (g1, g2)) - @test ft2 isa FusionTensor{Float64,4} - @test axes(ft2) == FusionTensorAxes((g1, g2), dual.((g1, g2))) - @test collect(eachblockstoredindex(data_matrix(ft2))) == map(i -> Block(i, i), 1:6) - for i in 1:6 - m = data_matrix(ft2)[Block(i, i)] - @test m == LinearAlgebra.I(size(m, 1)) - end - - ft2 = FusionTensor(3 * LinearAlgebra.I, (g1, g2)) - @test ft2 isa FusionTensor{Float64,4} - @test axes(ft2) == FusionTensorAxes((g1, g2), dual.((g1, g2))) - @test collect(eachblockstoredindex(data_matrix(ft2))) == map(i -> Block(i, i), 1:6) - for i in 1:6 - m = data_matrix(ft2)[Block(i, i)] - @test m == 3 * LinearAlgebra.I(size(m, 1)) - end - - @test FusionTensor{ComplexF64}(LinearAlgebra.I, (g1, g2)) isa FusionTensor{ComplexF64,4} + g1 = gradedrange([U1(0) => 1, U1(1) => 2, U1(2) => 3]) + g2 = gradedrange([U1(0) => 2, U1(1) => 2, U1(3) => 1]) + g3 = gradedrange([U1(-1) => 1, U1(0) => 2, U1(1) => 1]) + g4 = gradedrange([U1(-1) => 1, U1(0) => 1, U1(1) => 1]) + + fta = FusionTensorAxes((g1,), (g2, g3)) + @test zeros(fta) isa FusionTensor{Float64, 3} + @test zeros(ComplexF64, fta) isa FusionTensor{ComplexF64, 3} + + rng = Random.default_rng() + ft1 = randn(rng, ComplexF64, fta) + @test ft1 isa FusionTensor{ComplexF64, 3} + @test all(!=(0), data_matrix(ft1)[Block(1, 5)]) + @test randn(rng, fta) isa FusionTensor{Float64, 3} + @test randn(ComplexF64, fta) isa FusionTensor{ComplexF64, 3} + @test randn(fta) isa FusionTensor{Float64, 3} + + ft2 = FusionTensor(LinearAlgebra.I, (g1, g2)) + @test ft2 isa FusionTensor{Float64, 4} + @test axes(ft2) == FusionTensorAxes((g1, g2), dual.((g1, g2))) + @test collect(eachblockstoredindex(data_matrix(ft2))) == map(i -> Block(i, i), 1:6) + for i in 1:6 + m = data_matrix(ft2)[Block(i, i)] + @test m == LinearAlgebra.I(size(m, 1)) + end + + ft2 = FusionTensor(3 * LinearAlgebra.I, (g1, g2)) + @test ft2 isa FusionTensor{Float64, 4} + @test axes(ft2) == FusionTensorAxes((g1, g2), dual.((g1, g2))) + @test collect(eachblockstoredindex(data_matrix(ft2))) == map(i -> Block(i, i), 1:6) + for i in 1:6 + m = data_matrix(ft2)[Block(i, i)] + @test m == 3 * LinearAlgebra.I(size(m, 1)) + end + + @test FusionTensor{ComplexF64}(LinearAlgebra.I, (g1, g2)) isa FusionTensor{ComplexF64, 4} end @testset "Base operations" begin - g1 = gradedrange([U1(0) => 1, U1(1) => 2, U1(2) => 3]) - g2 = gradedrange([U1(0) => 2, U1(1) => 2, U1(3) => 1]) - g3 = gradedrange([U1(-1) => 1, U1(0) => 2, U1(1) => 1]) - g4 = gradedrange([U1(-1) => 1, U1(0) => 1, U1(1) => 1]) - ft3 = randn(FusionTensorAxes((g1, g2), (g3, g4))) - @test ft3 isa FusionTensor{Float64,4} - @test norm(ft3) ≉ 0 - @test isnothing(check_sanity(ft3)) - - ft4 = +ft3 - @test ft4 === ft3 # same object - - ft4 = -ft3 - @test isnothing(check_sanity(ft4)) - @test codomain_axes(ft4) === codomain_axes(ft3) - @test domain_axes(ft4) === domain_axes(ft3) - @test norm(ft4) ≈ norm(ft3) - @test norm(ft4 + ft3) ≈ 0.0 - - ft4 = ft3 + ft3 - @test codomain_axes(ft4) === codomain_axes(ft3) - @test domain_axes(ft4) === domain_axes(ft3) - @test space_isequal(codomain_axis(ft4), codomain_axis(ft3)) - @test space_isequal(domain_axis(ft4), domain_axis(ft3)) - @test isnothing(check_sanity(ft4)) - @test norm(ft4) ≈ 2norm(ft3) - - ft4 = ft3 - ft3 - @test codomain_axes(ft4) === codomain_axes(ft3) - @test domain_axes(ft4) === domain_axes(ft3) - @test space_isequal(codomain_axis(ft4), codomain_axis(ft3)) - @test space_isequal(domain_axis(ft4), domain_axis(ft3)) - @test isnothing(check_sanity(ft4)) - @test norm(ft4) ≈ 0.0 - - ft4 = 2 * ft3 - @test codomain_axes(ft4) === codomain_axes(ft3) - @test domain_axes(ft4) === domain_axes(ft3) - @test space_isequal(codomain_axis(ft4), codomain_axis(ft3)) - @test space_isequal(domain_axis(ft4), domain_axis(ft3)) - @test isnothing(check_sanity(ft4)) - @test eltype(ft4) == Float64 - @test norm(ft4) ≈ 2norm(ft3) - - ft4 = 2.0 * ft3 - @test codomain_axes(ft4) === codomain_axes(ft3) - @test domain_axes(ft4) === domain_axes(ft3) - @test space_isequal(codomain_axis(ft4), codomain_axis(ft3)) - @test space_isequal(domain_axis(ft4), domain_axis(ft3)) - @test isnothing(check_sanity(ft4)) - @test eltype(ft4) == Float64 - - ft4 = ft3 / 2.0 - @test codomain_axes(ft4) === codomain_axes(ft3) - @test domain_axes(ft4) === domain_axes(ft3) - @test space_isequal(codomain_axis(ft4), codomain_axis(ft3)) - @test space_isequal(domain_axis(ft4), domain_axis(ft3)) - @test isnothing(check_sanity(ft4)) - @test eltype(ft4) == Float64 - @test norm(ft4) ≈ norm(ft3) / 2.0 - - ft5 = (1.0 + 2.0im) * ft3 - @test codomain_axes(ft5) === codomain_axes(ft3) - @test domain_axes(ft5) === domain_axes(ft3) - @test space_isequal(codomain_axis(ft5), codomain_axis(ft3)) - @test space_isequal(domain_axis(ft5), domain_axis(ft3)) - @test isnothing(check_sanity(ft4)) - @test eltype(ft5) == ComplexF64 - @test norm(ft5) ≈ √5 * norm(ft3) - - @test conj(ft3) === ft3 # same object - @test real(ft3) === ft3 - @test norm(imag(ft3)) == 0 - - @test conj(ft5) isa FusionTensor{ComplexF64,4} - @test real(ft5) isa FusionTensor{Float64,4} - @test imag(ft3) isa FusionTensor{Float64,4} - @test conj(ft5) ≈ (1.0 - 2.0im) * ft3 - @test real(ft5) ≈ ft3 - @test imag(ft5) ≈ 2ft3 - - ft6 = conj(ft5) - @test ft6 !== ft5 # different object - @test isnothing(check_sanity(ft6)) - @test codomain_axes(ft6) === codomain_axes(ft5) - @test domain_axes(ft6) === domain_axes(ft5) - @test space_isequal(codomain_axis(ft6), codomain_axis(ft5)) - @test space_isequal(domain_axis(ft6), domain_axis(ft5)) - @test eltype(ft6) == ComplexF64 - @test ft6 + ft5 ≈ 2 * real(ft5) - - ad = adjoint(ft3) - @test ad isa FusionTensor - @test ndims_codomain(ad) == 2 - @test ndims_domain(ad) == 2 - @test space_isequal(dual(g1), domain_axes(ad)[1]) - @test space_isequal(dual(g2), domain_axes(ad)[2]) - @test space_isequal(dual(g3), codomain_axes(ad)[1]) - @test space_isequal(dual(g4), codomain_axes(ad)[2]) - @test isnothing(check_sanity(ad)) - - ft7 = FusionTensor{Float64}(undef, (g1,), (g2, g3, g4)) - @test_throws ArgumentError ft7 + ft3 - @test_throws ArgumentError ft7 - ft3 - @test_throws ArgumentError ft7 * ft3 + g1 = gradedrange([U1(0) => 1, U1(1) => 2, U1(2) => 3]) + g2 = gradedrange([U1(0) => 2, U1(1) => 2, U1(3) => 1]) + g3 = gradedrange([U1(-1) => 1, U1(0) => 2, U1(1) => 1]) + g4 = gradedrange([U1(-1) => 1, U1(0) => 1, U1(1) => 1]) + ft3 = randn(FusionTensorAxes((g1, g2), (g3, g4))) + @test ft3 isa FusionTensor{Float64, 4} + @test norm(ft3) ≉ 0 + @test isnothing(check_sanity(ft3)) + + ft4 = +ft3 + @test ft4 === ft3 # same object + + ft4 = -ft3 + @test isnothing(check_sanity(ft4)) + @test codomain_axes(ft4) === codomain_axes(ft3) + @test domain_axes(ft4) === domain_axes(ft3) + @test norm(ft4) ≈ norm(ft3) + @test norm(ft4 + ft3) ≈ 0.0 + + ft4 = ft3 + ft3 + @test codomain_axes(ft4) === codomain_axes(ft3) + @test domain_axes(ft4) === domain_axes(ft3) + @test space_isequal(codomain_axis(ft4), codomain_axis(ft3)) + @test space_isequal(domain_axis(ft4), domain_axis(ft3)) + @test isnothing(check_sanity(ft4)) + @test norm(ft4) ≈ 2norm(ft3) + + ft4 = ft3 - ft3 + @test codomain_axes(ft4) === codomain_axes(ft3) + @test domain_axes(ft4) === domain_axes(ft3) + @test space_isequal(codomain_axis(ft4), codomain_axis(ft3)) + @test space_isequal(domain_axis(ft4), domain_axis(ft3)) + @test isnothing(check_sanity(ft4)) + @test norm(ft4) ≈ 0.0 + + ft4 = 2 * ft3 + @test codomain_axes(ft4) === codomain_axes(ft3) + @test domain_axes(ft4) === domain_axes(ft3) + @test space_isequal(codomain_axis(ft4), codomain_axis(ft3)) + @test space_isequal(domain_axis(ft4), domain_axis(ft3)) + @test isnothing(check_sanity(ft4)) + @test eltype(ft4) == Float64 + @test norm(ft4) ≈ 2norm(ft3) + + ft4 = 2.0 * ft3 + @test codomain_axes(ft4) === codomain_axes(ft3) + @test domain_axes(ft4) === domain_axes(ft3) + @test space_isequal(codomain_axis(ft4), codomain_axis(ft3)) + @test space_isequal(domain_axis(ft4), domain_axis(ft3)) + @test isnothing(check_sanity(ft4)) + @test eltype(ft4) == Float64 + + ft4 = ft3 / 2.0 + @test codomain_axes(ft4) === codomain_axes(ft3) + @test domain_axes(ft4) === domain_axes(ft3) + @test space_isequal(codomain_axis(ft4), codomain_axis(ft3)) + @test space_isequal(domain_axis(ft4), domain_axis(ft3)) + @test isnothing(check_sanity(ft4)) + @test eltype(ft4) == Float64 + @test norm(ft4) ≈ norm(ft3) / 2.0 + + ft5 = (1.0 + 2.0im) * ft3 + @test codomain_axes(ft5) === codomain_axes(ft3) + @test domain_axes(ft5) === domain_axes(ft3) + @test space_isequal(codomain_axis(ft5), codomain_axis(ft3)) + @test space_isequal(domain_axis(ft5), domain_axis(ft3)) + @test isnothing(check_sanity(ft4)) + @test eltype(ft5) == ComplexF64 + @test norm(ft5) ≈ √5 * norm(ft3) + + @test conj(ft3) === ft3 # same object + @test real(ft3) === ft3 + @test norm(imag(ft3)) == 0 + + @test conj(ft5) isa FusionTensor{ComplexF64, 4} + @test real(ft5) isa FusionTensor{Float64, 4} + @test imag(ft3) isa FusionTensor{Float64, 4} + @test conj(ft5) ≈ (1.0 - 2.0im) * ft3 + @test real(ft5) ≈ ft3 + @test imag(ft5) ≈ 2ft3 + + ft6 = conj(ft5) + @test ft6 !== ft5 # different object + @test isnothing(check_sanity(ft6)) + @test codomain_axes(ft6) === codomain_axes(ft5) + @test domain_axes(ft6) === domain_axes(ft5) + @test space_isequal(codomain_axis(ft6), codomain_axis(ft5)) + @test space_isequal(domain_axis(ft6), domain_axis(ft5)) + @test eltype(ft6) == ComplexF64 + @test ft6 + ft5 ≈ 2 * real(ft5) + + ad = adjoint(ft3) + @test ad isa FusionTensor + @test ndims_codomain(ad) == 2 + @test ndims_domain(ad) == 2 + @test space_isequal(dual(g1), domain_axes(ad)[1]) + @test space_isequal(dual(g2), domain_axes(ad)[2]) + @test space_isequal(dual(g3), codomain_axes(ad)[1]) + @test space_isequal(dual(g4), codomain_axes(ad)[2]) + @test isnothing(check_sanity(ad)) + + ft7 = FusionTensor{Float64}(undef, (g1,), (g2, g3, g4)) + @test_throws ArgumentError ft7 + ft3 + @test_throws ArgumentError ft7 - ft3 + @test_throws ArgumentError ft7 * ft3 end @testset "missing SectorProduct" begin - g1 = gradedrange([SectorProduct(U1(1)) => 1]) - g2 = gradedrange([SectorProduct(U1(1), SU2(1//2)) => 1]) - g3 = gradedrange([SectorProduct(U1(1), SU2(1//2), Z{2}(1)) => 1]) - S = sector_type(g3) - - ft = FusionTensor{Float64}(undef, (g1,), (dual(g2), dual(g3))) - @test sector_type(ft) === S - gr = gradedrange([SectorProduct(U1(1), SU2(0), Z{2}(0)) => 1]) - @test space_isequal(codomain_axis(ft), gr) - gc = gradedrange([ - SectorProduct(U1(2), SU2(0), Z{2}(1)) => 1, SectorProduct(U1(2), SU2(1), Z{2}(1)) => 1 - ]) - @test space_isequal(domain_axis(ft), dual(gc)) - - gA = gradedrange([SectorProduct(; A=U1(1)) => 1]) - gB = gradedrange([SectorProduct(; B=SU2(1//2)) => 1]) - gC = gradedrange([SectorProduct(; C=Z{2}(0)) => 1]) - gABC = tensor_product(gA, gB, gC) - S = sector_type(gABC) - - ft = FusionTensor{Float64}(undef, (gA, gB), (dual(gA), dual(gB), gC)) - @test sector_type(ft) === S - @test space_isequal(codomain_axis(ft), gABC) - @test space_isequal(domain_axis(ft), dual(gABC)) + g1 = gradedrange([SectorProduct(U1(1)) => 1]) + g2 = gradedrange([SectorProduct(U1(1), SU2(1 // 2)) => 1]) + g3 = gradedrange([SectorProduct(U1(1), SU2(1 // 2), Z{2}(1)) => 1]) + S = sector_type(g3) + + ft = FusionTensor{Float64}(undef, (g1,), (dual(g2), dual(g3))) + @test sector_type(ft) === S + gr = gradedrange([SectorProduct(U1(1), SU2(0), Z{2}(0)) => 1]) + @test space_isequal(codomain_axis(ft), gr) + gc = gradedrange( + [ + SectorProduct(U1(2), SU2(0), Z{2}(1)) => 1, SectorProduct(U1(2), SU2(1), Z{2}(1)) => 1, + ] + ) + @test space_isequal(domain_axis(ft), dual(gc)) + + gA = gradedrange([SectorProduct(; A = U1(1)) => 1]) + gB = gradedrange([SectorProduct(; B = SU2(1 // 2)) => 1]) + gC = gradedrange([SectorProduct(; C = Z{2}(0)) => 1]) + gABC = tensor_product(gA, gB, gC) + S = sector_type(gABC) + + ft = FusionTensor{Float64}(undef, (gA, gB), (dual(gA), dual(gB), gC)) + @test sector_type(ft) === S + @test space_isequal(codomain_axis(ft), gABC) + @test space_isequal(domain_axis(ft), dual(gABC)) end diff --git a/test/test_contraction.jl b/test/test_contraction.jl index c07d2e9..b80d14d 100644 --- a/test/test_contraction.jl +++ b/test/test_contraction.jl @@ -3,146 +3,146 @@ using Test: @test, @testset, @test_broken, @test_throws using BlockSparseArrays: BlockSparseArray using FusionTensors: - FusionTensors, FusionTensor, FusionTensorAxes, domain_axes, codomain_axes, to_fusiontensor + FusionTensors, FusionTensor, FusionTensorAxes, domain_axes, codomain_axes, to_fusiontensor using GradedArrays: SU2, U1, dual, gradedrange using TensorAlgebra: - TensorAlgebra, contract, matricize, permmortar, tuplemortar, unmatricize, unmatricize! + TensorAlgebra, contract, matricize, permmortar, tuplemortar, unmatricize, unmatricize! include("setup.jl") @testset "abelian matricize" begin - g1 = gradedrange([U1(0) => 1, U1(1) => 2, U1(2) => 3]) - g2 = gradedrange([U1(0) => 2, U1(1) => 2, U1(3) => 1]) - g3 = gradedrange([U1(-1) => 1, U1(0) => 2, U1(1) => 1]) - g4 = gradedrange([U1(-1) => 1, U1(0) => 1, U1(1) => 1]) - - ft1 = randn(FusionTensorAxes((g1, g2), (dual(g3), dual(g4)))) - m = matricize(ft1, (1, 2), (3, 4)) - ft2 = unmatricize(m, axes(ft1)) - @test ft1 ≈ ft2 - - biperm = permmortar(((3,), (1, 2, 4))) - m2 = matricize(ft1, biperm) - ft_dest = similar(ft1, axes(ft1)[biperm]) - unmatricize!(ft_dest, m2, permmortar(((1,), (2, 3, 4)))) - @test ft_dest ≈ permutedims(ft1, biperm) - - ft2 = similar(ft1) - unmatricize!(ft2, m2, biperm) - @test ft1 ≈ ft2 + g1 = gradedrange([U1(0) => 1, U1(1) => 2, U1(2) => 3]) + g2 = gradedrange([U1(0) => 2, U1(1) => 2, U1(3) => 1]) + g3 = gradedrange([U1(-1) => 1, U1(0) => 2, U1(1) => 1]) + g4 = gradedrange([U1(-1) => 1, U1(0) => 1, U1(1) => 1]) + + ft1 = randn(FusionTensorAxes((g1, g2), (dual(g3), dual(g4)))) + m = matricize(ft1, (1, 2), (3, 4)) + ft2 = unmatricize(m, axes(ft1)) + @test ft1 ≈ ft2 + + biperm = permmortar(((3,), (1, 2, 4))) + m2 = matricize(ft1, biperm) + ft_dest = similar(ft1, axes(ft1)[biperm]) + unmatricize!(ft_dest, m2, permmortar(((1,), (2, 3, 4)))) + @test ft_dest ≈ permutedims(ft1, biperm) + + ft2 = similar(ft1) + unmatricize!(ft2, m2, biperm) + @test ft1 ≈ ft2 end @testset "non-abelian matricize" begin - g1 = gradedrange([SU2(0) => 1, SU2(1//2) => 2, SU2(1) => 3]) - g2 = gradedrange([SU2(0) => 2, SU2(1//2) => 2, SU2(3//2) => 1]) - g3 = gradedrange([SU2(1//2) => 1, SU2(1) => 2, SU2(2) => 1]) - g4 = gradedrange([SU2(0) => 1, SU2(1) => 1, SU2(3//2) => 1]) - - ft1 = randn(FusionTensorAxes((g1, g2), (dual(g3), dual(g4)))) - m = matricize(ft1, (1, 2), (3, 4)) - ft2 = unmatricize(m, axes(ft1)) - @test ft1 ≈ ft2 - - biperm = permmortar(((3,), (1, 2, 4))) - m2 = matricize(ft1, biperm) - ft_dest = similar(ft1, axes(ft1)[biperm]) - unmatricize!(ft_dest, m2, permmortar(((1,), (2, 3, 4)))) - @test ft_dest ≈ permutedims(ft1, biperm) - - ft2 = similar(ft1) - unmatricize!(ft2, m2, biperm) - @test ft1 ≈ ft2 + g1 = gradedrange([SU2(0) => 1, SU2(1 // 2) => 2, SU2(1) => 3]) + g2 = gradedrange([SU2(0) => 2, SU2(1 // 2) => 2, SU2(3 // 2) => 1]) + g3 = gradedrange([SU2(1 // 2) => 1, SU2(1) => 2, SU2(2) => 1]) + g4 = gradedrange([SU2(0) => 1, SU2(1) => 1, SU2(3 // 2) => 1]) + + ft1 = randn(FusionTensorAxes((g1, g2), (dual(g3), dual(g4)))) + m = matricize(ft1, (1, 2), (3, 4)) + ft2 = unmatricize(m, axes(ft1)) + @test ft1 ≈ ft2 + + biperm = permmortar(((3,), (1, 2, 4))) + m2 = matricize(ft1, biperm) + ft_dest = similar(ft1, axes(ft1)[biperm]) + unmatricize!(ft_dest, m2, permmortar(((1,), (2, 3, 4)))) + @test ft_dest ≈ permutedims(ft1, biperm) + + ft2 = similar(ft1) + unmatricize!(ft2, m2, biperm) + @test ft1 ≈ ft2 end @testset "Matrix functions" begin - sds22 = [ - 0.25 0.0 0.0 0.0 - 0.0 -0.25 0.5 0.0 - 0.0 0.5 -0.25 0.0 - 0.0 0.0 0.0 0.25 - ] - t = reshape(sds22, (2, 2, 2, 2)) - g2 = gradedrange([SU2(1//2) => 1]) - ft = to_fusiontensor(t, (g2, g2), (dual(g2), dual(g2))) - for f in setdiff(FusionTensors.MATRIX_FUNCTIONS, [:acoth, :cbrt]) - t2 = reshape((@eval Base.$f)(sds22), (2, 2, 2, 2)) - ft2 = to_fusiontensor(t2, (g2, g2), (dual(g2), dual(g2))) - @test (@eval TensorAlgebra.$f)(ft, (1, 2, 3, 4), (1, 2), (3, 4)) ≈ ft2 - end - @test_throws ArgumentError TensorAlgebra.exp(ft, (1, 2, 3, 4), (1, 2, 3), (4,)) - @test_throws ArgumentError TensorAlgebra.exp(ft, (1, 2, 3, 4), (1, 3), (2, 4)) + sds22 = [ + 0.25 0.0 0.0 0.0 + 0.0 -0.25 0.5 0.0 + 0.0 0.5 -0.25 0.0 + 0.0 0.0 0.0 0.25 + ] + t = reshape(sds22, (2, 2, 2, 2)) + g2 = gradedrange([SU2(1 // 2) => 1]) + ft = to_fusiontensor(t, (g2, g2), (dual(g2), dual(g2))) + for f in setdiff(FusionTensors.MATRIX_FUNCTIONS, [:acoth, :cbrt]) + t2 = reshape((@eval Base.$f)(sds22), (2, 2, 2, 2)) + ft2 = to_fusiontensor(t2, (g2, g2), (dual(g2), dual(g2))) + @test (@eval TensorAlgebra.$f)(ft, (1, 2, 3, 4), (1, 2), (3, 4)) ≈ ft2 + end + @test_throws ArgumentError TensorAlgebra.exp(ft, (1, 2, 3, 4), (1, 2, 3), (4,)) + @test_throws ArgumentError TensorAlgebra.exp(ft, (1, 2, 3, 4), (1, 3), (2, 4)) end @testset "contraction" begin - g1 = gradedrange([U1(0) => 1, U1(1) => 2, U1(2) => 3]) - g2 = gradedrange([U1(0) => 2, U1(1) => 2, U1(3) => 1]) - g3 = gradedrange([U1(-1) => 1, U1(0) => 2, U1(1) => 1]) - g4 = gradedrange([U1(-1) => 1, U1(0) => 1, U1(1) => 1]) + g1 = gradedrange([U1(0) => 1, U1(1) => 2, U1(2) => 3]) + g2 = gradedrange([U1(0) => 2, U1(1) => 2, U1(3) => 1]) + g3 = gradedrange([U1(-1) => 1, U1(0) => 2, U1(1) => 1]) + g4 = gradedrange([U1(-1) => 1, U1(0) => 1, U1(1) => 1]) - ft1 = FusionTensor{Float64}(undef, (g1, g2), (g3, g4)) - @test isnothing(check_sanity(ft1)) + ft1 = FusionTensor{Float64}(undef, (g1, g2), (g3, g4)) + @test isnothing(check_sanity(ft1)) - ft2 = FusionTensor{Float64}(undef, dual.((g3, g4)), (g1,)) - @test isnothing(check_sanity(ft2)) + ft2 = FusionTensor{Float64}(undef, dual.((g3, g4)), (g1,)) + @test isnothing(check_sanity(ft2)) - ft3 = ft1 * ft2 # tensor contraction - @test isnothing(check_sanity(ft3)) - @test domain_axes(ft3) === domain_axes(ft2) - @test codomain_axes(ft3) === codomain_axes(ft1) + ft3 = ft1 * ft2 # tensor contraction + @test isnothing(check_sanity(ft3)) + @test domain_axes(ft3) === domain_axes(ft2) + @test codomain_axes(ft3) === codomain_axes(ft1) - # test LinearAlgebra.mul! with in-place matrix product - m1 = randn(FusionTensorAxes((g1,), (g2,))) - m2 = randn(FusionTensorAxes((dual(g2),), (g3,))) - m3 = FusionTensor{Float64}(undef, (g1,), (g3,)) + # test LinearAlgebra.mul! with in-place matrix product + m1 = randn(FusionTensorAxes((g1,), (g2,))) + m2 = randn(FusionTensorAxes((dual(g2),), (g3,))) + m3 = FusionTensor{Float64}(undef, (g1,), (g3,)) - mul!(m3, m1, m2, 2.0, 0.0) - @test m3 ≈ 2m1 * m2 + mul!(m3, m1, m2, 2.0, 0.0) + @test m3 ≈ 2m1 * m2 end @testset "TensorAlgebra.contract interface" begin - g1 = gradedrange([U1(0) => 1, U1(1) => 2, U1(2) => 3]) - g2 = gradedrange([U1(0) => 2, U1(1) => 2, U1(3) => 1]) - g3 = gradedrange([U1(-1) => 1, U1(0) => 2, U1(1) => 1]) - g4 = gradedrange([U1(-1) => 1, U1(0) => 1, U1(1) => 1]) - - ft1 = randn(FusionTensorAxes((g1, g2), (g3, g4))) - ft2 = randn(FusionTensorAxes(dual.((g3, g4)), (dual(g1),))) - ft3 = randn(FusionTensorAxes(dual.((g3, g4)), dual.((g1, g2)))) - - ft4, legs = contract(ft1, (1, 2, 3, 4), ft2, (3, 4, 5)) - @test legs == tuplemortar(((1, 2), (5,))) - @test isnothing(check_sanity(ft4)) - @test domain_axes(ft4) === domain_axes(ft2) - @test codomain_axes(ft4) === codomain_axes(ft1) - @test ft4 ≈ ft1 * ft2 - - ft5 = contract((1, 2, 5), ft1, (1, 2, 3, 4), ft2, (3, 4, 5)) - @test isnothing(check_sanity(ft5)) - @test ndims_codomain(ft5) === 3 - @test ndims_domain(ft5) === 0 - @test permutedims(ft5, (1, 2), (3,)) ≈ ft4 - - ft6 = contract(tuplemortar(((1, 2), (5,))), ft1, (1, 2, 3, 4), ft2, (3, 4, 5)) - @test isnothing(check_sanity(ft6)) - @test ft4 ≈ ft6 - - @test permutedims(ft1, (), (1, 2, 3, 4)) * permutedims(ft3, (3, 4, 1, 2), ()) isa - FusionTensor{Float64,0} - ft7, legs = contract(ft1, (1, 2, 3, 4), ft3, (3, 4, 1, 2)) - @test legs == tuplemortar(((), ())) - @test ft7 isa FusionTensor{Float64,0} - - # include permutations - ft6 = contract(tuplemortar(((5, 1), (2,))), ft1, (1, 2, 3, 4), ft2, (3, 4, 5)) - @test isnothing(check_sanity(ft6)) - @test permutedims(ft6, (2, 3), (1,)) ≈ ft4 - - ft8 = contract( - tuplemortar(((-3,), (-1, -2, -4))), ft1, (-1, 1, -2, 2), ft3, (-3, 2, -4, 1) - ) - left = permutedims(ft1, (1, 3), (2, 4)) - right = permutedims(ft3, (4, 2), (1, 3)) - lrprod = left * right - newft = permutedims(lrprod, (3,), (1, 2, 4)) - @test newft ≈ ft8 + g1 = gradedrange([U1(0) => 1, U1(1) => 2, U1(2) => 3]) + g2 = gradedrange([U1(0) => 2, U1(1) => 2, U1(3) => 1]) + g3 = gradedrange([U1(-1) => 1, U1(0) => 2, U1(1) => 1]) + g4 = gradedrange([U1(-1) => 1, U1(0) => 1, U1(1) => 1]) + + ft1 = randn(FusionTensorAxes((g1, g2), (g3, g4))) + ft2 = randn(FusionTensorAxes(dual.((g3, g4)), (dual(g1),))) + ft3 = randn(FusionTensorAxes(dual.((g3, g4)), dual.((g1, g2)))) + + ft4, legs = contract(ft1, (1, 2, 3, 4), ft2, (3, 4, 5)) + @test legs == tuplemortar(((1, 2), (5,))) + @test isnothing(check_sanity(ft4)) + @test domain_axes(ft4) === domain_axes(ft2) + @test codomain_axes(ft4) === codomain_axes(ft1) + @test ft4 ≈ ft1 * ft2 + + ft5 = contract((1, 2, 5), ft1, (1, 2, 3, 4), ft2, (3, 4, 5)) + @test isnothing(check_sanity(ft5)) + @test ndims_codomain(ft5) === 3 + @test ndims_domain(ft5) === 0 + @test permutedims(ft5, (1, 2), (3,)) ≈ ft4 + + ft6 = contract(tuplemortar(((1, 2), (5,))), ft1, (1, 2, 3, 4), ft2, (3, 4, 5)) + @test isnothing(check_sanity(ft6)) + @test ft4 ≈ ft6 + + @test permutedims(ft1, (), (1, 2, 3, 4)) * permutedims(ft3, (3, 4, 1, 2), ()) isa + FusionTensor{Float64, 0} + ft7, legs = contract(ft1, (1, 2, 3, 4), ft3, (3, 4, 1, 2)) + @test legs == tuplemortar(((), ())) + @test ft7 isa FusionTensor{Float64, 0} + + # include permutations + ft6 = contract(tuplemortar(((5, 1), (2,))), ft1, (1, 2, 3, 4), ft2, (3, 4, 5)) + @test isnothing(check_sanity(ft6)) + @test permutedims(ft6, (2, 3), (1,)) ≈ ft4 + + ft8 = contract( + tuplemortar(((-3,), (-1, -2, -4))), ft1, (-1, 1, -2, 2), ft3, (-3, 2, -4, 1) + ) + left = permutedims(ft1, (1, 3), (2, 4)) + right = permutedims(ft3, (4, 2), (1, 3)) + lrprod = left * right + newft = permutedims(lrprod, (3,), (1, 2, 4)) + @test newft ≈ ft8 end diff --git a/test/test_examples.jl b/test/test_examples.jl index 380a628..5cf46f4 100644 --- a/test/test_examples.jl +++ b/test/test_examples.jl @@ -3,6 +3,6 @@ using FusionTensors: FusionTensors using Test: @test, @testset @testset "examples" begin - include(joinpath(pkgdir(FusionTensors), "examples", "README.jl")) + include(joinpath(pkgdir(FusionTensors), "examples", "README.jl")) end end diff --git a/test/test_fusion_trees.jl b/test/test_fusion_trees.jl index 565997b..7e2aa90 100644 --- a/test/test_fusion_trees.jl +++ b/test/test_fusion_trees.jl @@ -4,161 +4,161 @@ using LinearAlgebra: I using BlockArrays: BlockArrays using FusionTensors: - SectorFusionTree, - arrows, - branch_sectors, - build_trees, - fusiontree_eltype, - leaves, - outer_multiplicity_indices, - root_sector + SectorFusionTree, + arrows, + branch_sectors, + build_trees, + fusiontree_eltype, + leaves, + outer_multiplicity_indices, + root_sector using GradedArrays: - ×, SectorProduct, SU, SU2, TrivialSector, arguments, dual, flip, gradedrange, sector_type + ×, SectorProduct, SU, SU2, TrivialSector, arguments, dual, flip, gradedrange, sector_type @testset "Trivial fusion trees" begin - q = TrivialSector() - f = SectorFusionTree{TrivialSector}() - @test arrows(f) == () - @test leaves(f) == () - @test root_sector(f) == q - @test branch_sectors(f) == () - @test sector_type(f) == TrivialSector - @test outer_multiplicity_indices(f) == () - @test convert(Array, f) ≈ ones((1,)) - - f = only(build_trees((q,), (true,))) - @test arrows(f) == (true,) - @test leaves(f) == (q,) - @test root_sector(f) == q - @test branch_sectors(f) == () - @test outer_multiplicity_indices(f) == () - @test convert(Array, f) ≈ ones((1, 1)) - - f = only(build_trees((q, q), (true, false))) - @test arrows(f) == (true, false) - @test leaves(f) == (q, q) - @test root_sector(f) == q - @test branch_sectors(f) == (q,) - @test outer_multiplicity_indices(f) == (1,) - @test convert(Array, f) ≈ ones((1, 1, 1)) - - @test fusiontree_eltype(sector_type(f)) === eltype(convert(Array, f)) + q = TrivialSector() + f = SectorFusionTree{TrivialSector}() + @test arrows(f) == () + @test leaves(f) == () + @test root_sector(f) == q + @test branch_sectors(f) == () + @test sector_type(f) == TrivialSector + @test outer_multiplicity_indices(f) == () + @test convert(Array, f) ≈ ones((1,)) + + f = only(build_trees((q,), (true,))) + @test arrows(f) == (true,) + @test leaves(f) == (q,) + @test root_sector(f) == q + @test branch_sectors(f) == () + @test outer_multiplicity_indices(f) == () + @test convert(Array, f) ≈ ones((1, 1)) + + f = only(build_trees((q, q), (true, false))) + @test arrows(f) == (true, false) + @test leaves(f) == (q, q) + @test root_sector(f) == q + @test branch_sectors(f) == (q,) + @test outer_multiplicity_indices(f) == (1,) + @test convert(Array, f) ≈ ones((1, 1, 1)) + + @test fusiontree_eltype(sector_type(f)) === eltype(convert(Array, f)) end @testset "SU(2) SectorFusionTree" begin - j2 = SU2(1//2) - - f = only(build_trees((j2,), (false,))) - @test arrows(f) == (false,) - @test leaves(f) == (j2,) - @test root_sector(f) == j2 - @test branch_sectors(f) == () - @test outer_multiplicity_indices(f) == () - @test sector_type(f) == typeof(j2) - @test convert(Array, f) ≈ I(2) - @test fusiontree_eltype(sector_type(f)) === eltype(convert(Array, f)) - - f = only(build_trees((j2,), (true,))) - @test arrows(f) == (true,) - @test convert(Array, f) ≈ [0 -1; 1 0] - - trees = build_trees((j2, j2), (false, false)) - @test length(trees) == 2 - f1 = first(trees) - @test root_sector(f1) == SU2(0) - @test branch_sectors(f1) == (SU2(1//2),) - @test outer_multiplicity_indices(f1) == (1,) - @test convert(Array, f1) ≈ [0 1/sqrt(2); -1/sqrt(2) 0] - - f3 = last(trees) - @test root_sector(f3) == SU2(1) - @test branch_sectors(f3) == (SU2(1//2),) - @test outer_multiplicity_indices(f3) == (1,) - t = zeros((2, 2, 3)) - t[1, 1, 1] = 1 - t[1, 2, 2] = 1 / sqrt(2) - t[2, 1, 2] = 1 / sqrt(2) - t[2, 2, 3] = 1 - @test convert(Array, f3) ≈ t - - trees = build_trees((j2, j2, j2), (false, false, false)) - @test length(trees) == 3 - f12, f32, f34 = trees - @test f12 < f32 < f34 - @test root_sector(f12) == SU2(1//2) - @test root_sector(f32) == SU2(1//2) - @test root_sector(f34) == SU2(3//2) - - @test branch_sectors(f12) == (SU2(1//2), SU2(0)) - @test branch_sectors(f32) == (SU2(1//2), SU2(1)) - @test branch_sectors(f34) == (SU2(1//2), SU2(1)) + j2 = SU2(1 // 2) + + f = only(build_trees((j2,), (false,))) + @test arrows(f) == (false,) + @test leaves(f) == (j2,) + @test root_sector(f) == j2 + @test branch_sectors(f) == () + @test outer_multiplicity_indices(f) == () + @test sector_type(f) == typeof(j2) + @test convert(Array, f) ≈ I(2) + @test fusiontree_eltype(sector_type(f)) === eltype(convert(Array, f)) + + f = only(build_trees((j2,), (true,))) + @test arrows(f) == (true,) + @test convert(Array, f) ≈ [0 -1; 1 0] + + trees = build_trees((j2, j2), (false, false)) + @test length(trees) == 2 + f1 = first(trees) + @test root_sector(f1) == SU2(0) + @test branch_sectors(f1) == (SU2(1 // 2),) + @test outer_multiplicity_indices(f1) == (1,) + @test convert(Array, f1) ≈ [0 1 / sqrt(2); -1 / sqrt(2) 0] + + f3 = last(trees) + @test root_sector(f3) == SU2(1) + @test branch_sectors(f3) == (SU2(1 // 2),) + @test outer_multiplicity_indices(f3) == (1,) + t = zeros((2, 2, 3)) + t[1, 1, 1] = 1 + t[1, 2, 2] = 1 / sqrt(2) + t[2, 1, 2] = 1 / sqrt(2) + t[2, 2, 3] = 1 + @test convert(Array, f3) ≈ t + + trees = build_trees((j2, j2, j2), (false, false, false)) + @test length(trees) == 3 + f12, f32, f34 = trees + @test f12 < f32 < f34 + @test root_sector(f12) == SU2(1 // 2) + @test root_sector(f32) == SU2(1 // 2) + @test root_sector(f34) == SU2(3 // 2) + + @test branch_sectors(f12) == (SU2(1 // 2), SU2(0)) + @test branch_sectors(f32) == (SU2(1 // 2), SU2(1)) + @test branch_sectors(f34) == (SU2(1 // 2), SU2(1)) end @testset "SU(3) SectorFusionTree" begin - # convention: irreps are already dualed if needed, arrows do not affect them. They only - # affect the basis on which the tree projects for self-dual irreps. - f3 = SU{3}((1, 0)) - c3 = dual(f3) - - trees = build_trees((f3, f3), (false, false)) - @test root_sector.(trees) == [SU{3}((1, 1)), SU{3}((2, 0))] - - trees = build_trees((f3, f3), (true, false)) - @test root_sector.(trees) == [SU{3}((1, 1)), SU{3}((2, 0))] - - trees = build_trees((f3, f3), (false, true)) - @test root_sector.(trees) == [SU{3}((1, 1)), SU{3}((2, 0))] - - trees = build_trees((f3, f3), (true, true)) - @test root_sector.(trees) == [SU{3}((1, 1)), SU{3}((2, 0))] - - trees = build_trees((f3, c3), (false, false)) - @test root_sector.(trees) == [SU{3}((0, 0)), SU{3}((2, 1))] - - trees = build_trees((c3, c3), (false, false)) - @test root_sector.(trees) == [SU{3}((1, 0)), SU{3}((2, 2))] - - # test GradedUnitRange interface - g = gradedrange([SU{3}((1, 0)) => 1]) - trees = build_trees(g, g) - @test root_sector.(trees) == [SU{3}((1, 1)), SU{3}((2, 0))] - trees = build_trees(g, flip(g)) - @test root_sector.(trees) == [SU{3}((1, 1)), SU{3}((2, 0))] - trees = build_trees(g, flip(dual(g))) - @test root_sector.(trees) == [SU{3}((0, 0)), SU{3}((2, 1))] - trees = build_trees(g, dual(g)) - @test root_sector.(trees) == [SU{3}((0, 0)), SU{3}((2, 1))] - - # test outer outer_multiplicity > 1 - a8 = SU{3}((2, 1)) - trees = build_trees((a8, a8), (false, false)) - @test length(trees) == 6 - f = first(trees) - @test root_sector(f) == SU{3}((0, 0)) - @test sector_type(f) == typeof(a8) - - f8a = trees[2] - f8b = trees[3] - @test root_sector(f8a) == a8 - @test root_sector(f8b) == a8 - @test branch_sectors(f8a) == (a8,) - @test branch_sectors(f8b) == (a8,) - @test outer_multiplicity_indices(f8a) == (1,) - @test outer_multiplicity_indices(f8b) == (2,) + # convention: irreps are already dualed if needed, arrows do not affect them. They only + # affect the basis on which the tree projects for self-dual irreps. + f3 = SU{3}((1, 0)) + c3 = dual(f3) + + trees = build_trees((f3, f3), (false, false)) + @test root_sector.(trees) == [SU{3}((1, 1)), SU{3}((2, 0))] + + trees = build_trees((f3, f3), (true, false)) + @test root_sector.(trees) == [SU{3}((1, 1)), SU{3}((2, 0))] + + trees = build_trees((f3, f3), (false, true)) + @test root_sector.(trees) == [SU{3}((1, 1)), SU{3}((2, 0))] + + trees = build_trees((f3, f3), (true, true)) + @test root_sector.(trees) == [SU{3}((1, 1)), SU{3}((2, 0))] + + trees = build_trees((f3, c3), (false, false)) + @test root_sector.(trees) == [SU{3}((0, 0)), SU{3}((2, 1))] + + trees = build_trees((c3, c3), (false, false)) + @test root_sector.(trees) == [SU{3}((1, 0)), SU{3}((2, 2))] + + # test GradedUnitRange interface + g = gradedrange([SU{3}((1, 0)) => 1]) + trees = build_trees(g, g) + @test root_sector.(trees) == [SU{3}((1, 1)), SU{3}((2, 0))] + trees = build_trees(g, flip(g)) + @test root_sector.(trees) == [SU{3}((1, 1)), SU{3}((2, 0))] + trees = build_trees(g, flip(dual(g))) + @test root_sector.(trees) == [SU{3}((0, 0)), SU{3}((2, 1))] + trees = build_trees(g, dual(g)) + @test root_sector.(trees) == [SU{3}((0, 0)), SU{3}((2, 1))] + + # test outer outer_multiplicity > 1 + a8 = SU{3}((2, 1)) + trees = build_trees((a8, a8), (false, false)) + @test length(trees) == 6 + f = first(trees) + @test root_sector(f) == SU{3}((0, 0)) + @test sector_type(f) == typeof(a8) + + f8a = trees[2] + f8b = trees[3] + @test root_sector(f8a) == a8 + @test root_sector(f8b) == a8 + @test branch_sectors(f8a) == (a8,) + @test branch_sectors(f8b) == (a8,) + @test outer_multiplicity_indices(f8a) == (1,) + @test outer_multiplicity_indices(f8b) == (2,) end @testset "SU(2)×SU(3) SectorFusionTree" begin - j2 = SU2(1//2) - a8 = SU{3}((2, 1)) - s = j2 × a8 - trees = build_trees((s, s), (false, false)) - @test length(trees) == 12 - f = first(trees) - @test sector_type(f) == typeof(s) - argument_trees = arguments(f) - @test length(argument_trees) == 2 - f2, f3 = argument_trees - @test sector_type(f2) == typeof(j2) - @test sector_type(f3) == typeof(a8) - @test f == f2 × f3 + j2 = SU2(1 // 2) + a8 = SU{3}((2, 1)) + s = j2 × a8 + trees = build_trees((s, s), (false, false)) + @test length(trees) == 12 + f = first(trees) + @test sector_type(f) == typeof(s) + argument_trees = arguments(f) + @test length(argument_trees) == 2 + f2, f3 = argument_trees + @test sector_type(f2) == typeof(j2) + @test sector_type(f3) == typeof(a8) + @test f == f2 × f3 end diff --git a/test/test_fusiontensoraxes.jl b/test/test_fusiontensoraxes.jl index f40df4e..d0c66b7 100644 --- a/test/test_fusiontensoraxes.jl +++ b/test/test_fusiontensoraxes.jl @@ -5,113 +5,113 @@ using BlockArrays: Block, blockedrange, blocklength, blocklengths, blocks using TensorAlgebra: BlockedTuple, length_codomain, trivial_axis, tuplemortar using FusionTensors: - FusionTensorAxes, - ndims_domain, - ndims_codomain, - codomain, - domain, - fused_codomain, - fused_domain, - promote_sector_type, - promote_sectors + FusionTensorAxes, + ndims_domain, + ndims_codomain, + codomain, + domain, + fused_codomain, + fused_domain, + promote_sector_type, + promote_sectors using GradedArrays: - ×, - U1, - SectorProduct, - TrivialSector, - SU2, - checkspaces, - dual, - gradedrange, - sector_type, - space_isequal + ×, + U1, + SectorProduct, + TrivialSector, + SU2, + checkspaces, + dual, + gradedrange, + sector_type, + space_isequal @testset "misc FusionTensors.jl" begin - g1 = gradedrange([U1(0) => 1]) - @test promote_sector_type(U1(1), U1(1)) === typeof(U1(1)) - @test promote_sector_type(g1, U1(1)) === typeof(U1(1)) - @test promote_sector_type(g1, g1) === typeof(U1(1)) - @test promote_sector_type((g1, g1)) === typeof(U1(1)) - - sNS = SectorProduct(; N=U1(1), S=SU2(1 / 2)) - gN = gradedrange([(; N=U1(1)) => 1]) - gS = gradedrange([(; S=SU2(1 / 2)) => 1]) - @test promote_sector_type(gN, gS) == typeof(sNS) - - @test promote_sectors((gN, gS)) isa NTuple{2,typeof(gradedrange([sNS => 1]))} + g1 = gradedrange([U1(0) => 1]) + @test promote_sector_type(U1(1), U1(1)) === typeof(U1(1)) + @test promote_sector_type(g1, U1(1)) === typeof(U1(1)) + @test promote_sector_type(g1, g1) === typeof(U1(1)) + @test promote_sector_type((g1, g1)) === typeof(U1(1)) + + sNS = SectorProduct(; N = U1(1), S = SU2(1 / 2)) + gN = gradedrange([(; N = U1(1)) => 1]) + gS = gradedrange([(; S = SU2(1 / 2)) => 1]) + @test promote_sector_type(gN, gS) == typeof(sNS) + + @test promote_sectors((gN, gS)) isa NTuple{2, typeof(gradedrange([sNS => 1]))} end @testset "FusionTensorAxes" begin - s2 = SU2(1//2) - g2 = gradedrange([s2 => 1]) - g2b = dual(g2) - - bt = tuplemortar(((g2, g2), (g2b, g2b))) - fta = FusionTensorAxes(bt) - - @test fta isa FusionTensorAxes - @test BlockedTuple(fta) == bt - - @test Tuple(fta) == Tuple(bt) - @test space_isequal(only(axes(fta)), blockedrange([2, 2])) - @test iterate(fta) == (g2, 2) - @test iterate(fta, 1) == (g2, 2) - @test iterate(fta, 2) == (g2, 3) - @test iterate(fta, 3) == (g2b, 4) - @test iterate(fta, 4) == (g2b, 5) - @test isnothing(iterate(fta, 5)) - - @test length(fta) == 4 - @test space_isequal(fta[1], g2) - @test space_isequal(fta[2], g2) - @test space_isequal(fta[3], g2b) - @test space_isequal(fta[4], g2b) - @test length(fta[Block(1)]) == 2 - @test all(map(r -> space_isequal(r, g2), fta[Block(1)])) - @test length(fta[Block(2)]) == 2 - @test all(map(r -> space_isequal(r, g2b), fta[Block(2)])) - @test length.(fta) == tuplemortar(((2, 2), (2, 2))) - - @test blocklength(fta) == 2 - @test blocklengths(fta) == (2, 2) - @test blocks(fta) == blocks(bt) - @test length_codomain(fta) == 2 - - @test sector_type(fta) === sector_type(g2) - @test length(codomain(fta)) == 2 - @test all(map(r -> space_isequal(r, g2), codomain(fta))) - @test length(domain(fta)) == 2 - @test all(map(r -> space_isequal(r, g2b), domain(fta))) - @test space_isequal(fused_codomain(fta), g2 ⊗ g2) - @test space_isequal(fused_domain(fta), dual(g2 ⊗ g2)) - @test space_isequal(trivial_axis(fta), trivial_axis(typeof(s2))) - - @test fta == fta - @test copy(fta) == fta - @test deepcopy(fta) == fta - @test fta != FusionTensorAxes(tuplemortar(((g2, g2), (g2b, g2)))) - @test fta != FusionTensorAxes(tuplemortar(((g2, g2, g2b), (g2b,)))) - - @test fta == FusionTensorAxes((g2, g2), (g2b, g2b)) - @test checkspaces(fta, fta) - @test_throws ArgumentError checkspaces( - fta, FusionTensorAxes(tuplemortar(((g2, g2), (g2b, g2)))) - ) + s2 = SU2(1 // 2) + g2 = gradedrange([s2 => 1]) + g2b = dual(g2) + + bt = tuplemortar(((g2, g2), (g2b, g2b))) + fta = FusionTensorAxes(bt) + + @test fta isa FusionTensorAxes + @test BlockedTuple(fta) == bt + + @test Tuple(fta) == Tuple(bt) + @test space_isequal(only(axes(fta)), blockedrange([2, 2])) + @test iterate(fta) == (g2, 2) + @test iterate(fta, 1) == (g2, 2) + @test iterate(fta, 2) == (g2, 3) + @test iterate(fta, 3) == (g2b, 4) + @test iterate(fta, 4) == (g2b, 5) + @test isnothing(iterate(fta, 5)) + + @test length(fta) == 4 + @test space_isequal(fta[1], g2) + @test space_isequal(fta[2], g2) + @test space_isequal(fta[3], g2b) + @test space_isequal(fta[4], g2b) + @test length(fta[Block(1)]) == 2 + @test all(map(r -> space_isequal(r, g2), fta[Block(1)])) + @test length(fta[Block(2)]) == 2 + @test all(map(r -> space_isequal(r, g2b), fta[Block(2)])) + @test length.(fta) == tuplemortar(((2, 2), (2, 2))) + + @test blocklength(fta) == 2 + @test blocklengths(fta) == (2, 2) + @test blocks(fta) == blocks(bt) + @test length_codomain(fta) == 2 + + @test sector_type(fta) === sector_type(g2) + @test length(codomain(fta)) == 2 + @test all(map(r -> space_isequal(r, g2), codomain(fta))) + @test length(domain(fta)) == 2 + @test all(map(r -> space_isequal(r, g2b), domain(fta))) + @test space_isequal(fused_codomain(fta), g2 ⊗ g2) + @test space_isequal(fused_domain(fta), dual(g2 ⊗ g2)) + @test space_isequal(trivial_axis(fta), trivial_axis(typeof(s2))) + + @test fta == fta + @test copy(fta) == fta + @test deepcopy(fta) == fta + @test fta != FusionTensorAxes(tuplemortar(((g2, g2), (g2b, g2)))) + @test fta != FusionTensorAxes(tuplemortar(((g2, g2, g2b), (g2b,)))) + + @test fta == FusionTensorAxes((g2, g2), (g2b, g2b)) + @test checkspaces(fta, fta) + @test_throws ArgumentError checkspaces( + fta, FusionTensorAxes(tuplemortar(((g2, g2), (g2b, g2)))) + ) end @testset "Empty FusionTensorAxes" begin - fta = FusionTensorAxes(tuplemortar(((), ()))) - @test fta isa FusionTensorAxes - - @test length(fta) == 0 - @test isempty(fta) - @test blocklength(fta) == 2 - @test blocklengths(fta) == (0, 0) - @test sector_type(fta) == TrivialSector - @test length_codomain(fta) == 0 - - @test codomain(fta) == () - @test space_isequal(fused_codomain(fta), trivial_axis(TrivialSector)) - @test domain(fta) == () - @test space_isequal(fused_domain(fta), dual(trivial_axis(TrivialSector))) + fta = FusionTensorAxes(tuplemortar(((), ()))) + @test fta isa FusionTensorAxes + + @test length(fta) == 0 + @test isempty(fta) + @test blocklength(fta) == 2 + @test blocklengths(fta) == (0, 0) + @test sector_type(fta) == TrivialSector + @test length_codomain(fta) == 0 + + @test codomain(fta) == () + @test space_isequal(fused_codomain(fta), trivial_axis(TrivialSector)) + @test domain(fta) == () + @test space_isequal(fused_domain(fta), dual(trivial_axis(TrivialSector))) end diff --git a/test/test_linear_algebra.jl b/test/test_linear_algebra.jl index 9b9f66c..2944eab 100644 --- a/test/test_linear_algebra.jl +++ b/test/test_linear_algebra.jl @@ -10,41 +10,41 @@ using GradedArrays: O2, SU2, TrivialSector, U1, dual, gradedrange include("setup.jl") @testset "LinearAlgebra interface" begin - sds22 = [ - 0.25 0.0 0.0 0.0 - 0.0 -0.25 0.5 0.0 - 0.0 0.5 -0.25 0.0 - 0.0 0.0 0.0 0.25 - ] - sdst = reshape(sds22, (2, 2, 2, 2)) - - g0 = gradedrange([TrivialSector() => 2]) - gu1 = gradedrange([U1(1) => 1, U1(-1) => 1]) - go2 = gradedrange([O2(1 / 2) => 1]) - gsu2 = gradedrange([SU2(1 / 2) => 1]) - - for g in [g0, gu1, go2, gsu2] - ft = to_fusiontensor(sdst, (g, g), (dual(g), dual(g))) - @test isnothing(check_sanity(ft)) - @test norm(ft) ≈ √3 / 2 - @test norm(ft, 2) ≈ √3 / 2 - @test norm(ft, 2.0) ≈ √3 / 2 - @test isapprox(tr(ft), 0; atol=eps(Float64)) - - ft2 = normalize(ft) - @test norm(ft2) ≈ 1.0 - @test norm(ft) ≈ √3 / 2 # unaffected by normalize - @test ft ≈ √3 / 2 * ft2 - normalize!(ft) - @test norm(ft) ≈ 1.0 - end - - for g in [g0, gu1] - ft = to_fusiontensor(sdst, (g, g), (dual(g), dual(g))) - @test norm(ft, 1) ≈ 2.0 - end - for g in [go2, gsu2] - ft = to_fusiontensor(sdst, (g, g), (dual(g), dual(g))) - @test norm(ft, 1) ≈ 1.5 - end + sds22 = [ + 0.25 0.0 0.0 0.0 + 0.0 -0.25 0.5 0.0 + 0.0 0.5 -0.25 0.0 + 0.0 0.0 0.0 0.25 + ] + sdst = reshape(sds22, (2, 2, 2, 2)) + + g0 = gradedrange([TrivialSector() => 2]) + gu1 = gradedrange([U1(1) => 1, U1(-1) => 1]) + go2 = gradedrange([O2(1 / 2) => 1]) + gsu2 = gradedrange([SU2(1 / 2) => 1]) + + for g in [g0, gu1, go2, gsu2] + ft = to_fusiontensor(sdst, (g, g), (dual(g), dual(g))) + @test isnothing(check_sanity(ft)) + @test norm(ft) ≈ √3 / 2 + @test norm(ft, 2) ≈ √3 / 2 + @test norm(ft, 2.0) ≈ √3 / 2 + @test isapprox(tr(ft), 0; atol = eps(Float64)) + + ft2 = normalize(ft) + @test norm(ft2) ≈ 1.0 + @test norm(ft) ≈ √3 / 2 # unaffected by normalize + @test ft ≈ √3 / 2 * ft2 + normalize!(ft) + @test norm(ft) ≈ 1.0 + end + + for g in [g0, gu1] + ft = to_fusiontensor(sdst, (g, g), (dual(g), dual(g))) + @test norm(ft, 1) ≈ 2.0 + end + for g in [go2, gsu2] + ft = to_fusiontensor(sdst, (g, g), (dual(g), dual(g))) + @test norm(ft, 1) ≈ 1.5 + end end diff --git a/test/test_permutedims.jl b/test/test_permutedims.jl index fc2886b..91bf15a 100644 --- a/test/test_permutedims.jl +++ b/test/test_permutedims.jl @@ -2,226 +2,226 @@ using Test: @test, @testset, @test_broken, @test_throws using BlockArrays: blocks using FusionTensors: - FusionTensor, - FusionTensorAxes, - data_matrix, - codomain_axis, - domain_axis, - ndims_domain, - ndims_codomain, - to_fusiontensor + FusionTensor, + FusionTensorAxes, + data_matrix, + codomain_axis, + domain_axis, + ndims_domain, + ndims_codomain, + to_fusiontensor using GradedArrays: ×, O2, U1, SectorProduct, SU2, dual, gradedrange, space_isequal using TensorAlgebra: permmortar, tuplemortar include("setup.jl") function naive_permutedims(ft, biperm) - @assert ndims(ft) == length(biperm) + @assert ndims(ft) == length(biperm) - # naive permute: cast to dense, permutedims, cast to FusionTensor - arr = Array(ft) - permuted_arr = permutedims(arr, Tuple(biperm)) - permuted = to_fusiontensor(permuted_arr, blocks(axes(ft)[biperm])...) - return permuted + # naive permute: cast to dense, permutedims, cast to FusionTensor + arr = Array(ft) + permuted_arr = permutedims(arr, Tuple(biperm)) + permuted = to_fusiontensor(permuted_arr, blocks(axes(ft)[biperm])...) + return permuted end @testset "Abelian permutedims" begin - @testset "dummy" begin - g1 = gradedrange([U1(0) => 1, U1(1) => 2, U1(2) => 3]) - g2 = gradedrange([U1(0) => 2, U1(1) => 2, U1(3) => 1]) - g3 = gradedrange([U1(-1) => 1, U1(0) => 2, U1(1) => 1]) - g4 = gradedrange([U1(-1) => 1, U1(0) => 1, U1(1) => 1]) - ftaxes1 = FusionTensorAxes((g1, g2), (dual(g3), dual(g4))) - - for elt in (Float64, ComplexF64) - ft1 = randn(elt, ftaxes1) - @test isnothing(check_sanity(ft1)) - - # test permutedims interface - ft2 = permutedims(ft1, (1, 2), (3, 4)) # trivial with 2 tuples - @test ft2 ≈ ft1 - @test ft2 !== ft1 - @test data_matrix(ft2) !== data_matrix(ft1) # check copy - @test data_matrix(ft2) == data_matrix(ft1) # check copy - - ft2 = permutedims(ft1, ((1, 2), (3, 4))) # trivial with tuple of 2 tuples - @test ft2 ≈ ft1 - @test ft2 !== ft1 - @test data_matrix(ft2) !== data_matrix(ft1) # check copy - @test data_matrix(ft2) == data_matrix(ft1) # check copy - - biperm = permmortar(((1, 2), (3, 4))) - ft2 = permutedims(ft1, biperm) # trivial with biperm - @test ft2 ≈ ft1 - @test ft2 !== ft1 - @test data_matrix(ft2) !== data_matrix(ft1) # check copy - @test data_matrix(ft2) == data_matrix(ft1) # check copy - - ft3 = permutedims(ft1, (4,), (1, 2, 3)) - @test ft3 !== ft1 - @test ft3 isa FusionTensor{elt,4} - @test axes(ft3) == FusionTensorAxes((dual(g4),), (g1, g2, dual(g3))) - @test isnothing(check_sanity(ft3)) - - ft4 = permutedims(ft3, (2, 3), (4, 1)) - @test axes(ft1) == axes(ft4) - @test space_isequal(codomain_axis(ft1), codomain_axis(ft4)) - @test space_isequal(domain_axis(ft1), domain_axis(ft4)) - @test ft4 ≈ ft1 - - # test permutedims! interface - ft2 = randn(elt, axes(ft1)) - permutedims!(ft2, ft1, (1, 2), (3, 4)) - @test ft2 ≈ ft1 - @test data_matrix(ft2) !== data_matrix(ft1) # check copy - @test data_matrix(ft2) == data_matrix(ft1) # check copy - - ft2 = randn(elt, axes(ft1)) - permutedims!(ft2, ft1, ((1, 2), (3, 4))) - @test ft2 ≈ ft1 - @test data_matrix(ft2) !== data_matrix(ft1) # check copy - @test data_matrix(ft2) == data_matrix(ft1) # check copy - - ft2 = randn(elt, axes(ft1)) - permutedims!(ft2, ft1, biperm) - @test ft2 ≈ ft1 - @test data_matrix(ft2) !== data_matrix(ft1) # check copy - @test data_matrix(ft2) == data_matrix(ft1) # check copy - - # test clean errors - ft2 = randn(elt, axes(ft1)) - @test_throws MethodError permutedims(ft1, (2, 3, 4, 1)) - @test_throws ArgumentError permutedims(ft1, (2, 3), (5, 4, 1)) - @test_throws MethodError permutedims!(ft2, ft1, (2, 3, 4, 1)) - @test_throws ArgumentError permutedims!(ft2, ft1, (2, 3), (5, 4, 1)) - @test_throws ArgumentError permutedims!(ft2, ft1, (1, 2, 3), (4,)) - @test_throws ArgumentError permutedims!(ft2, ft1, (1, 2), (4, 3)) + @testset "dummy" begin + g1 = gradedrange([U1(0) => 1, U1(1) => 2, U1(2) => 3]) + g2 = gradedrange([U1(0) => 2, U1(1) => 2, U1(3) => 1]) + g3 = gradedrange([U1(-1) => 1, U1(0) => 2, U1(1) => 1]) + g4 = gradedrange([U1(-1) => 1, U1(0) => 1, U1(1) => 1]) + ftaxes1 = FusionTensorAxes((g1, g2), (dual(g3), dual(g4))) + + for elt in (Float64, ComplexF64) + ft1 = randn(elt, ftaxes1) + @test isnothing(check_sanity(ft1)) + + # test permutedims interface + ft2 = permutedims(ft1, (1, 2), (3, 4)) # trivial with 2 tuples + @test ft2 ≈ ft1 + @test ft2 !== ft1 + @test data_matrix(ft2) !== data_matrix(ft1) # check copy + @test data_matrix(ft2) == data_matrix(ft1) # check copy + + ft2 = permutedims(ft1, ((1, 2), (3, 4))) # trivial with tuple of 2 tuples + @test ft2 ≈ ft1 + @test ft2 !== ft1 + @test data_matrix(ft2) !== data_matrix(ft1) # check copy + @test data_matrix(ft2) == data_matrix(ft1) # check copy + + biperm = permmortar(((1, 2), (3, 4))) + ft2 = permutedims(ft1, biperm) # trivial with biperm + @test ft2 ≈ ft1 + @test ft2 !== ft1 + @test data_matrix(ft2) !== data_matrix(ft1) # check copy + @test data_matrix(ft2) == data_matrix(ft1) # check copy + + ft3 = permutedims(ft1, (4,), (1, 2, 3)) + @test ft3 !== ft1 + @test ft3 isa FusionTensor{elt, 4} + @test axes(ft3) == FusionTensorAxes((dual(g4),), (g1, g2, dual(g3))) + @test isnothing(check_sanity(ft3)) + + ft4 = permutedims(ft3, (2, 3), (4, 1)) + @test axes(ft1) == axes(ft4) + @test space_isequal(codomain_axis(ft1), codomain_axis(ft4)) + @test space_isequal(domain_axis(ft1), domain_axis(ft4)) + @test ft4 ≈ ft1 + + # test permutedims! interface + ft2 = randn(elt, axes(ft1)) + permutedims!(ft2, ft1, (1, 2), (3, 4)) + @test ft2 ≈ ft1 + @test data_matrix(ft2) !== data_matrix(ft1) # check copy + @test data_matrix(ft2) == data_matrix(ft1) # check copy + + ft2 = randn(elt, axes(ft1)) + permutedims!(ft2, ft1, ((1, 2), (3, 4))) + @test ft2 ≈ ft1 + @test data_matrix(ft2) !== data_matrix(ft1) # check copy + @test data_matrix(ft2) == data_matrix(ft1) # check copy + + ft2 = randn(elt, axes(ft1)) + permutedims!(ft2, ft1, biperm) + @test ft2 ≈ ft1 + @test data_matrix(ft2) !== data_matrix(ft1) # check copy + @test data_matrix(ft2) == data_matrix(ft1) # check copy + + # test clean errors + ft2 = randn(elt, axes(ft1)) + @test_throws MethodError permutedims(ft1, (2, 3, 4, 1)) + @test_throws ArgumentError permutedims(ft1, (2, 3), (5, 4, 1)) + @test_throws MethodError permutedims!(ft2, ft1, (2, 3, 4, 1)) + @test_throws ArgumentError permutedims!(ft2, ft1, (2, 3), (5, 4, 1)) + @test_throws ArgumentError permutedims!(ft2, ft1, (1, 2, 3), (4,)) + @test_throws ArgumentError permutedims!(ft2, ft1, (1, 2), (4, 3)) + end end - end - - @testset "Many axes" begin - g1 = gradedrange([U1(1) => 2, U1(2) => 2]) - g2 = gradedrange([U1(2) => 3, U1(3) => 2]) - g3 = gradedrange([U1(3) => 4, U1(4) => 1]) - g4 = gradedrange([U1(0) => 2, U1(2) => 1]) - codomain_legs = (g1, g2) - domain_legs = dual.((g3, g4)) - arr = zeros(ComplexF64, (4, 5, 5, 3)) - arr[1:2, 1:3, 1:4, 1:2] .= 1.0im - arr[3:4, 1:3, 5:5, 1:2] .= 2.0 - arr[1:2, 4:5, 5:5, 1:2] .= 3.0 - arr[3:4, 4:5, 1:4, 3:3] .= 4.0 - ft = to_fusiontensor(arr, codomain_legs, domain_legs) - biperm = permmortar(((3,), (2, 4, 1))) - - ftp = permutedims(ft, biperm) - @test ftp ≈ naive_permutedims(ft, biperm) - ftpp = permutedims(ftp, (4, 2), (1, 3)) - @test ftpp ≈ ft - - ft2 = adjoint(ft) - ftp2 = permutedims(ft2, biperm) - @test ftp2 ≈ naive_permutedims(ft2, biperm) - ftpp2 = permutedims(ftp2, (4, 2), (1, 3)) - @test ftpp2 ≈ ft2 - @test adjoint(ftpp2) ≈ ft - end - - @testset "Less than two axes" begin - if VERSION >= v"1.11" - ft0 = to_fusiontensor(ones(()), (), ()) - ft0p = permutedims(ft0, (), ()) - @test ft0p isa FusionTensor{Float64,0} - @test data_matrix(ft0p) ≈ data_matrix(ft0) - @test ft0p ≈ ft0 - - @test permutedims(ft0, ((), ())) isa FusionTensor{Float64,0} - @test permutedims(ft0, permmortar(((), ()))) isa FusionTensor{Float64,0} + + @testset "Many axes" begin + g1 = gradedrange([U1(1) => 2, U1(2) => 2]) + g2 = gradedrange([U1(2) => 3, U1(3) => 2]) + g3 = gradedrange([U1(3) => 4, U1(4) => 1]) + g4 = gradedrange([U1(0) => 2, U1(2) => 1]) + codomain_legs = (g1, g2) + domain_legs = dual.((g3, g4)) + arr = zeros(ComplexF64, (4, 5, 5, 3)) + arr[1:2, 1:3, 1:4, 1:2] .= 1.0im + arr[3:4, 1:3, 5:5, 1:2] .= 2.0 + arr[1:2, 4:5, 5:5, 1:2] .= 3.0 + arr[3:4, 4:5, 1:4, 3:3] .= 4.0 + ft = to_fusiontensor(arr, codomain_legs, domain_legs) + biperm = permmortar(((3,), (2, 4, 1))) + + ftp = permutedims(ft, biperm) + @test ftp ≈ naive_permutedims(ft, biperm) + ftpp = permutedims(ftp, (4, 2), (1, 3)) + @test ftpp ≈ ft + + ft2 = adjoint(ft) + ftp2 = permutedims(ft2, biperm) + @test ftp2 ≈ naive_permutedims(ft2, biperm) + ftpp2 = permutedims(ftp2, (4, 2), (1, 3)) + @test ftpp2 ≈ ft2 + @test adjoint(ftpp2) ≈ ft end - g = gradedrange([U1(0) => 1, U1(1) => 2, U1(2) => 3]) - v = zeros((6,)) - v[1] = 1.0 - biperm = permmortar(((), (1,))) - ft1 = to_fusiontensor(v, (g,), ()) - ft2 = permutedims(ft1, biperm) - @test isnothing(check_sanity(ft2)) - @test ft2 ≈ naive_permutedims(ft1, biperm) - ft3 = permutedims(ft2, (1,), ()) - @test ft1 ≈ ft3 - end + @testset "Less than two axes" begin + if VERSION >= v"1.11" + ft0 = to_fusiontensor(ones(()), (), ()) + ft0p = permutedims(ft0, (), ()) + @test ft0p isa FusionTensor{Float64, 0} + @test data_matrix(ft0p) ≈ data_matrix(ft0) + @test ft0p ≈ ft0 + + @test permutedims(ft0, ((), ())) isa FusionTensor{Float64, 0} + @test permutedims(ft0, permmortar(((), ()))) isa FusionTensor{Float64, 0} + end + + g = gradedrange([U1(0) => 1, U1(1) => 2, U1(2) => 3]) + v = zeros((6,)) + v[1] = 1.0 + biperm = permmortar(((), (1,))) + ft1 = to_fusiontensor(v, (g,), ()) + ft2 = permutedims(ft1, biperm) + @test isnothing(check_sanity(ft2)) + @test ft2 ≈ naive_permutedims(ft1, biperm) + ft3 = permutedims(ft2, (1,), ()) + @test ft1 ≈ ft3 + end end @testset "Non-abelian permutedims" begin - sds22 = reshape( - [ - [0.25, 0.0, 0.0, 0.0] - [0.0, -0.25, 0.5, 0.0] - [0.0, 0.5, -0.25, 0.0] - [0.0, 0.0, 0.0, 0.25] - ], - (2, 2, 2, 2), - ) - - sds22b = reshape( - [ - [-0.25, 0.0, 0.0, -0.5] - [0.0, 0.25, 0.0, 0.0] - [0.0, 0.0, 0.25, 0.0] - [-0.5, 0.0, 0.0, -0.25] - ], - (2, 2, 2, 2), - ) - - for g2 in ( - gradedrange([O2(1//2) => 1]), - dual(gradedrange([O2(1//2) => 1])), - gradedrange([SU2(1//2) => 1]), - dual(gradedrange([SU2(1//2) => 1])), - ) - g2b = dual(g2) - for biperm in [ - permmortar(((2, 1), (3, 4))), - permmortar(((3, 1), (2, 4))), - permmortar(((3, 1, 4), (2,))), - ] - ft = to_fusiontensor(sds22, (g2, g2), (g2b, g2b)) - @test permutedims(ft, biperm) ≈ naive_permutedims(ft, biperm) - @test permutedims(adjoint(ft), biperm) ≈ naive_permutedims(adjoint(ft), biperm) - - ft = to_fusiontensor(sds22b, (g2, g2b), (g2b, g2)) - @test permutedims(ft, biperm) ≈ naive_permutedims(ft, biperm) - @test permutedims(adjoint(ft), biperm) ≈ naive_permutedims(adjoint(ft), biperm) - end - for biperm in [permmortar(((1, 2, 3, 4), ())), permmortar(((), (3, 1, 2, 4)))] - ft = to_fusiontensor(sds22, (g2, g2), (g2b, g2b)) - @test permutedims(ft, biperm) ≈ naive_permutedims(ft, biperm) + sds22 = reshape( + [ + [0.25, 0.0, 0.0, 0.0] + [0.0, -0.25, 0.5, 0.0] + [0.0, 0.5, -0.25, 0.0] + [0.0, 0.0, 0.0, 0.25] + ], + (2, 2, 2, 2), + ) + + sds22b = reshape( + [ + [-0.25, 0.0, 0.0, -0.5] + [0.0, 0.25, 0.0, 0.0] + [0.0, 0.0, 0.25, 0.0] + [-0.5, 0.0, 0.0, -0.25] + ], + (2, 2, 2, 2), + ) + + for g2 in ( + gradedrange([O2(1 // 2) => 1]), + dual(gradedrange([O2(1 // 2) => 1])), + gradedrange([SU2(1 // 2) => 1]), + dual(gradedrange([SU2(1 // 2) => 1])), + ) + g2b = dual(g2) + for biperm in [ + permmortar(((2, 1), (3, 4))), + permmortar(((3, 1), (2, 4))), + permmortar(((3, 1, 4), (2,))), + ] + ft = to_fusiontensor(sds22, (g2, g2), (g2b, g2b)) + @test permutedims(ft, biperm) ≈ naive_permutedims(ft, biperm) + @test permutedims(adjoint(ft), biperm) ≈ naive_permutedims(adjoint(ft), biperm) + + ft = to_fusiontensor(sds22b, (g2, g2b), (g2b, g2)) + @test permutedims(ft, biperm) ≈ naive_permutedims(ft, biperm) + @test permutedims(adjoint(ft), biperm) ≈ naive_permutedims(adjoint(ft), biperm) + end + for biperm in [permmortar(((1, 2, 3, 4), ())), permmortar(((), (3, 1, 2, 4)))] + ft = to_fusiontensor(sds22, (g2, g2), (g2b, g2b)) + @test permutedims(ft, biperm) ≈ naive_permutedims(ft, biperm) + end end - end end @testset "SectorProduct permutedims" begin - d = 2 - D = 3 - tRVB = zeros((d, D, D, D, D)) # tensor RVB SU(2) for spin s - for i in 1:d - tRVB[i, i + 1, 1, 1, 1] = 1.0 - tRVB[i, 1, i + 1, 1, 1] = 1.0 - tRVB[i, 1, 1, i + 1, 1] = 1.0 - tRVB[i, 1, 1, 1, i + 1] = 1.0 - end - - gd = gradedrange([SU2(1//2) × U1(3) => 1]) - gD = dual(gradedrange([SU2(0) × U1(1) => 1, SU2(1//2) × U1(0) => 1])) - ft = to_fusiontensor(tRVB, (gd,), (gD, gD, gD, gD)) - @test Array(ft) ≈ tRVB - for biperm in [ - permmortar(((1,), (2, 3, 4, 5))), - permmortar(((1, 2, 3), (4, 5))), - permmortar(((3, 1, 4), (2, 5))), - permmortar(((), (2, 4, 1, 5, 3))), - permmortar(((2, 4, 1, 5, 3), ())), - ] - @test permutedims(ft, biperm) ≈ naive_permutedims(ft, biperm) - end + d = 2 + D = 3 + tRVB = zeros((d, D, D, D, D)) # tensor RVB SU(2) for spin s + for i in 1:d + tRVB[i, i + 1, 1, 1, 1] = 1.0 + tRVB[i, 1, i + 1, 1, 1] = 1.0 + tRVB[i, 1, 1, i + 1, 1] = 1.0 + tRVB[i, 1, 1, 1, i + 1] = 1.0 + end + + gd = gradedrange([SU2(1 // 2) × U1(3) => 1]) + gD = dual(gradedrange([SU2(0) × U1(1) => 1, SU2(1 // 2) × U1(0) => 1])) + ft = to_fusiontensor(tRVB, (gd,), (gD, gD, gD, gD)) + @test Array(ft) ≈ tRVB + for biperm in [ + permmortar(((1,), (2, 3, 4, 5))), + permmortar(((1, 2, 3), (4, 5))), + permmortar(((3, 1, 4), (2, 5))), + permmortar(((), (2, 4, 1, 5, 3))), + permmortar(((2, 4, 1, 5, 3), ())), + ] + @test permutedims(ft, biperm) ≈ naive_permutedims(ft, biperm) + end end