diff --git a/Project.toml b/Project.toml index 1d82813..ae16ecc 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,39 @@ name = "FusionTensors" uuid = "e16ca583-1f51-4df0-8e12-57d32947d33e" authors = ["ITensor developers and contributors"] -version = "0.1.0" +version = "0.2.0" + +[deps] +Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697" +BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e" +BlockSparseArrays = "2c9a651f-6452-4ace-a6ac-809f4280fbb4" +BroadcastMapConversion = "4a4adec5-520f-4750-bb37-d5e66b4ddeb2" +GradedUnitRanges = "e2de450a-8a67-46c7-b59c-01d5a3d041c5" +HalfIntegers = "f0d1745a-41c9-11e9-1dd9-e5d34d218721" +LRUCache = "8ac3fa9e-de4c-5943-b1dc-09c6b5f20637" +LabelledNumbers = "f856a3a6-4152-4ec4-b2a7-02c1a55d7993" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +SparseArraysBase = "0d5efcca-f356-4864-8770-e1ed8d78f208" +Strided = "5e0ebb24-38b0-5f93-81fe-25c709ecae67" +SymmetrySectors = "f8a8ad64-adbc-4fce-92f7-ffe2bb36a86e" +TensorAlgebra = "68bd88dc-f39d-4e12-b2ca-f046b68fcc6a" +TypeParameterAccessors = "7e5a90cf-f82e-492e-a09b-e3e26432c138" +WignerSymbols = "9f57e263-0b3d-5e2e-b1be-24f2bb48858b" [compat] +Accessors = "0.1.39" +BlockArrays = "1.2.0" +BlockSparseArrays = "0.2.0" +BroadcastMapConversion = "0.1.0" +GradedUnitRanges = "0.1.1" +HalfIntegers = "1.6.0" +LRUCache = "1.6.1" +LabelledNumbers = "0.1.0" +LinearAlgebra = "1.10.0" +SparseArraysBase = "0.2.0" +Strided = "2.2.0" +SymmetrySectors = "0.1.1" +TensorAlgebra = "0.1.0" +TypeParameterAccessors = "0.2.0" +WignerSymbols = "2.0.0" julia = "1.10" diff --git a/src/FusionTensors.jl b/src/FusionTensors.jl index 336ea97..d17d579 100644 --- a/src/FusionTensors.jl +++ b/src/FusionTensors.jl @@ -1,5 +1,14 @@ module FusionTensors -# Write your package code here. +include("fusion_trees/fusiontree.jl") +include("fusion_trees/clebsch_gordan_tensors.jl") +include("fusiontensor/fusedaxes.jl") +include("fusiontensor/fusiontensor.jl") +include("fusiontensor/base_interface.jl") +include("fusiontensor/array_cast.jl") +include("fusiontensor/linear_algebra_interface.jl") +include("fusiontensor/tensor_algebra_interface.jl") +include("permutedims/unitaries.jl") +include("permutedims/permutedims.jl") end diff --git a/src/fusion_trees/clebsch_gordan_tensors.jl b/src/fusion_trees/clebsch_gordan_tensors.jl new file mode 100644 index 0000000..fd426e6 --- /dev/null +++ b/src/fusion_trees/clebsch_gordan_tensors.jl @@ -0,0 +1,136 @@ +# This file defines Clebsch-Gordan tensors +# one tensor is defined from 3 simple objects s1, s2 and s3 +# and contains the coefficients fusing s1 ⊗ s2 -> s3 + +using HalfIntegers: half +using WignerSymbols: clebschgordan + +using GradedUnitRanges: dual, fusion_product +using SymmetrySectors: + ⊗, + AbelianStyle, + AbstractSector, + NotAbelianStyle, + SymmetryStyle, + O2, + SU, + istrivial, + quantum_dimension, + sector_label, + trivial, + zero_odd +using TensorAlgebra: contract + +function symbol_1j(s::AbstractSector) + 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 +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) +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)) +end + +function clebsch_gordan_tensor(::NotAbelianStyle, s1::O2, s2::O2, s3::O2, ::Int) + 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 ∉ blocklabels(fusion_product(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 + 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 +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 +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 +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 +end diff --git a/src/fusion_trees/fusiontree.jl b/src/fusion_trees/fusiontree.jl new file mode 100644 index 0000000..6d2294d --- /dev/null +++ b/src/fusion_trees/fusiontree.jl @@ -0,0 +1,324 @@ +# This file defines fusion trees for any abelian or non-abelian group + +# TBD +# compatibility with TensorKit conventions? + +using GradedUnitRanges: + AbstractGradedUnitRange, GradedUnitRanges, fusion_product, isdual, sector_type +using SymmetrySectors: ×, AbstractSector, SectorProduct, SymmetrySectors, arguments, trivial +using TensorAlgebra: flatten_tuples + +# +# A fusion tree fuses N sectors sec1, secN onto one sector fused_sec. A given set of +# sectors and arrow directions (as defined by a given outer block) contains several fusion +# trees that typically fuse to several sectors (in the abelian group case, there is only one) +# irrep in the fusion ring and each of them corresponds to a single "thin" fusion tree with +# +# +# +# / +# sec123 +# /\ +# / \ +# sec12 \ +# /\ \ +# / \ \ +# sec1 sec2 sec3 +# +# +# +# +# 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. +# +# +# The interface uses AbstractGradedUnitRanges as input for interface simplicity +# however only blocklabels 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 +end + +# getters +arrows(f::SectorFusionTree) = f.arrows +leaves(f::SectorFusionTree) = f.leaves +root_sector(f::SectorFusionTree) = f.root_sector +branch_sectors(f::SectorFusionTree) = f.branch_sectors +outer_multiplicity_indices(f::SectorFusionTree) = f.outer_multiplicity_indices + +# Base interface +Base.convert(T::Type{<:Array}, f::SectorFusionTree) = convert(T, to_array(f)) +Base.isless(f1::SectorFusionTree, f2::SectorFusionTree) = isless(to_tuple(f1), to_tuple(f2)) +Base.length(::SectorFusionTree{<:Any,N}) where {N} = N + +# GradedUnitRanges interface +GradedUnitRanges.sector_type(::SectorFusionTree{S}) where {S} = S + +function build_trees(legs::Vararg{AbstractGradedUnitRange}) + tree_arrows = isdual.(legs) + sectors = blocklabels.(legs) + return mapreduce(vcat, CartesianIndices(blocklength.(legs))) do it + block_sectors = getindex.(sectors, Tuple(it)) # why not type stable? + return build_trees(block_sectors, tree_arrows) + end +end + +# SymmetrySectors interface +function SymmetrySectors.:×(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, + ) +end + +function SymmetrySectors.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), + ) +end + +function SymmetrySectors.arguments(f::SectorFusionTree{<:SectorProduct,0}) + return map(arg -> SectorFusionTree((), (), arg, (), ()), arguments(root_sector(f))) +end + +function SymmetrySectors.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 + +# +# ===================================== Internals ======================================== +# +# --------------- misc --------------- +function to_tuple(f::SectorFusionTree) + return ( + leaves(f)..., + arrows(f)..., + root_sector(f), + branch_sectors(f)..., + outer_multiplicity_indices(f)..., + ) +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] +end + +# TODO move to GradedUnitRanges +function nsymbol(s1::AbstractSector, s2::AbstractSector, s3::AbstractSector) + full_space = fusion_product(s1, s2) + x = findfirst(==(s3), blocklabels(full_space)) + isnothing(x) && return 0 # OR labelled(0, s3)? + return Int(blocklengths(full_space)[x]) +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]) +end + +# --------------- Build trees --------------- +# zero leg: need S to get sector type information +function SectorFusionTree{S}() where {S<:AbstractSector} + return SectorFusionTree((), (), trivial(S), (), ()) +end +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, (), ()) +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) +end + +function grow_tree( + 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 grow_tree( + parent_tree::SectorFusionTree, branch_sector::AbstractSector, level_arrow::Bool +) + new_space = fusion_product(root_sector(parent_tree), branch_sector) + return mapreduce(vcat, zip(blocklabels(new_space), blocklengths(new_space))) do (la, n) + return [ + grow_tree(parent_tree, branch_sector, level_arrow, la, outer_mult) for + outer_mult in 1:n + ] + end +end + +function build_trees(old_trees::Vector, sectors_to_fuse::Tuple, arrows_to_fuse::Tuple) + next_level_trees = mapreduce(vcat, old_trees) do tree + return grow_tree(tree, first(sectors_to_fuse), first(arrows_to_fuse)) + end + return build_trees( + next_level_trees, Base.tail(sectors_to_fuse), Base.tail(arrows_to_fuse) + ) +end + +function build_trees(trees::Vector, ::Tuple{}, ::Tuple{}) + return trees +end + +function build_trees( + sectors_to_fuse::NTuple{N,<:AbstractSector}, arrows_to_fuse::NTuple{N,Bool} +) where {N} + trees = [SectorFusionTree(first(sectors_to_fuse), first(arrows_to_fuse))] + return build_trees(trees, Base.tail(sectors_to_fuse), Base.tail(arrows_to_fuse)) +end + +# --------------- convert to Array --------------- +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) +end + +function to_array(f::SectorFusionTree{<:SectorProduct}) + 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)) +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), + ) +end + +# specialized code when branch_sector is empty +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) +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) +end diff --git a/src/fusiontensor/array_cast.jl b/src/fusiontensor/array_cast.jl new file mode 100644 index 0000000..afd31da --- /dev/null +++ b/src/fusiontensor/array_cast.jl @@ -0,0 +1,200 @@ +# This file defines interface to cast from and to generic array + +using BlockArrays: AbstractBlockArray, BlockedArray, blockedrange, blocklengths, findblock + +using BlockSparseArrays: BlockSparseArrays, BlockSparseArray +using GradedUnitRanges: AbstractGradedUnitRange, blocklabels +using SymmetrySectors: block_dimensions, quantum_dimension +using TensorAlgebra: contract + +# ================================= High level interface ================================= + +#### cast from array to symmetric +function FusionTensor( + array::AbstractArray, + codomain_legs::Tuple{Vararg{AbstractGradedUnitRange}}, + domain_legs::Tuple{Vararg{AbstractGradedUnitRange}}, +) + return to_fusiontensor(array, codomain_legs, domain_legs) +end + +#### cast from symmetric to array +function BlockSparseArrays.BlockSparseArray(ft::FusionTensor) + return to_array(ft) +end + +# ================================= Low level interface ================================== +function to_fusiontensor( + array::AbstractArray, + codomain_legs::Tuple{Vararg{AbstractGradedUnitRange}}, + domain_legs::Tuple{Vararg{AbstractGradedUnitRange}}, +) + bounds = block_dimensions.((codomain_legs..., domain_legs...)) + blockarray = BlockedArray(array, bounds...) + return to_fusiontensor(blockarray, codomain_legs, domain_legs) +end + +rtoldefault(a::AbstractArray) = rtoldefault(eltype(a)) +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=atol, rtol=rtol) || throw( + InexactError( + :FusionTensor, typeof(a), typeof(codomain_axes(ft)), typeof(domain_axes(ft)) + ), + ) +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 +end + +function to_fusiontensor_no_checknorm( + blockarray::AbstractBlockArray, + codomain_legs::Tuple{Vararg{AbstractGradedUnitRange}}, + domain_legs::Tuple{Vararg{AbstractGradedUnitRange}}, +) + ft = FusionTensor(eltype(blockarray), 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 = block_dimensions.((codomain_axes(ft)..., domain_axes(ft)...)) + bsa = BlockSparseArray{eltype(ft)}(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 ======================================== + +#----------------------------------- 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) +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) +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 +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---------------- / + return contract( + ntuple(i -> 2 * i - 1, N), + split_array_block, + ntuple(identity, 2 * N), + p, + ntuple(i -> 2 * i, N), + 1 / dim_sec, # normalization factor + ) +end + +function contract_singlet_projector(f1::SectorFusionTree, f2::SectorFusionTree) + f1_array = convert(Array, f1) + f2_array = convert(Array, 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) + + # TODO use contract once it supports outer product + swapped = reshape(charge_block, :, 1) * reshape(p, 1, :) + block_shape = (size(charge_block)..., size(p)...) + perm = braid_tuples(ntuple(identity, N), ntuple(i -> i + N, N)) + split_array_block = permutedims(reshape(swapped, block_shape), perm) + + return merge_degen_dims(split_array_block) +end diff --git a/src/fusiontensor/base_interface.jl b/src/fusiontensor/base_interface.jl new file mode 100644 index 0000000..7df5a4d --- /dev/null +++ b/src/fusiontensor/base_interface.jl @@ -0,0 +1,121 @@ +# This files defines Base functions for FusionTensor + +using Accessors: @set + +using BlockSparseArrays: @view! + +set_data_matrix(ft::FusionTensor, data_matrix) = @set ft.data_matrix = data_matrix + +Base.:*(x::Number, ft::FusionTensor) = set_data_matrix(ft, x * data_matrix(ft)) +Base.:*(ft::FusionTensor, x::Number) = set_data_matrix(ft, x * data_matrix(ft)) + +# tensor contraction is a block data_matrix product. +function Base.:*(left::FusionTensor, right::FusionTensor) + checkaxes_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 + +Base.:+(ft::FusionTensor) = ft + +# tensor addition is a block data_matrix add. +function Base.:+(left::FusionTensor, right::FusionTensor) + checkaxes(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) + checkaxes(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) + +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) +transpose_mapping(d::Dict) = Dict([reverse(k) => transpose_mapping(v) for (k, v) in d]) +function transpose_mapping(b::BlockIndexRange{2}) + new_block = Block(reverse(Tuple(Block(b)))) + return new_block[reverse(b.indices)...] +end +function Base.adjoint(ft::FusionTensor) + return FusionTensor( + adjoint(data_matrix(ft)), + dual.(domain_axes(ft)), + dual.(codomain_axes(ft)), + transpose_mapping(trees_block_mapping(ft)), + ) +end + +Base.axes(ft::FusionTensor) = (codomain_axes(ft)..., domain_axes(ft)...) + +# conj is defined as coefficient wise complex conjugation, without axis dual +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.(codomain_axes(ft)), + copy.(domain_axes(ft)), + copy(trees_block_mapping(ft)), + ) +end + +function Base.deepcopy(ft::FusionTensor) + return FusionTensor( + deepcopy(data_matrix(ft)), + deepcopy.(codomain_axes(ft)), + deepcopy.(domain_axes(ft)), + deepcopy(trees_block_mapping(ft)), + ) +end + +# eachindex is automatically defined for AbstractArray. We do not want it. +Base.eachindex(::FusionTensor) = error("eachindex not defined for 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)) +end + +function Base.setindex!( + ft::FusionTensor, a::AbstractArray, f1::SectorFusionTree, f2::SectorFusionTree +) + return view(ft, f1, f2) .= a +end + +Base.permutedims(ft::FusionTensor, args...) = fusiontensor_permutedims(ft, args...) + +function Base.similar(ft::FusionTensor, ::Type{T}) where {T} + # fusion trees have Float64 eltype: need compatible type + @assert promote_type(T, Float64) === T + mat = similar(data_matrix(ft), T) + initialize_allowed_sectors!(mat) + return FusionTensor(mat, codomain_axes(ft), domain_axes(ft), trees_block_mapping(ft)) +end + +function Base.similar(::FusionTensor, ::Type{T}, new_axes::Tuple{<:Tuple,<:Tuple}) where {T} + return FusionTensor(T, new_axes[1], new_axes[2]) +end + +Base.show(io::IO, ft::FusionTensor) = print(io, "$(ndims(ft))-dim FusionTensor") + +function Base.show(io::IO, ::MIME"text/plain", ft::FusionTensor) + println(io, "$(ndims(ft))-dim FusionTensor with axes:") + for ax in axes(ft) + println(io, 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)) +end diff --git a/src/fusiontensor/fusedaxes.jl b/src/fusiontensor/fusedaxes.jl new file mode 100644 index 0000000..f3fa152 --- /dev/null +++ b/src/fusiontensor/fusedaxes.jl @@ -0,0 +1,98 @@ +# This file defines helper functions to access FusionTensor internal structures + +using BlockArrays: Block, BlockIndexRange, blocklength, blocklengths + +using BlockSparseArrays: to_block_indices +using GradedUnitRanges: + AbstractGradedUnitRange, GradedUnitRanges, blocklabels, blockmergesort, gradedrange +using SymmetrySectors: AbstractSector, trivial + +struct FusedAxes{A,B,C} + outer_axes::A + fused_axis::B + trees_to_ranges_mapping::C + + function FusedAxes( + outer_legs::NTuple{N,AbstractGradedUnitRange{LA}}, + fused_axis::AbstractGradedUnitRange{LA}, + trees_to_ranges_mapping::Dict{<:SectorFusionTree{<:AbstractSector,N}}, + ) where {N,LA} + return new{typeof(outer_legs),typeof(fused_axis),typeof(trees_to_ranges_mapping)}( + outer_legs, fused_axis, trees_to_ranges_mapping + ) + end +end + +# getters +fused_axis(fa::FusedAxes) = fa.fused_axis +fusion_trees(fa::FusedAxes) = keys(trees_to_ranges_mapping(fa)) +trees_to_ranges_mapping(fa::FusedAxes) = fa.trees_to_ranges_mapping + +# Base interface +Base.axes(fa::FusedAxes) = fa.outer_axes +Base.ndims(fa::FusedAxes) = length(axes(fa)) + +# GradedUnitRanges interface +GradedUnitRanges.blocklabels(fa::FusedAxes) = blocklabels(fused_axis(fa)) + +# constructors +function FusedAxes{S}(::Tuple{}) where {S<:AbstractSector} + fused_axis = gradedrange([trivial(S) => 1]) + trees_to_ranges_mapping = Dict([SectorFusionTree{S}() => Block(1)[1:1]]) + return FusedAxes((), fused_axis, trees_to_ranges_mapping) +end + +function FusedAxes{S}( + outer_legs::Tuple{Vararg{AbstractGradedUnitRange}} +) where {S<:AbstractSector} + fusion_trees_mult = fusion_trees_external_multiplicities(outer_legs) + + fused_leg, range_mapping = compute_inner_ranges(fusion_trees_mult) + return FusedAxes(outer_legs, fused_leg, range_mapping) +end + +function fusion_trees_external_multiplicities( + outer_legs::Tuple{Vararg{AbstractGradedUnitRange}} +) + tree_arrows = isdual.(outer_legs) + return mapreduce(vcat, CartesianIndices(blocklength.(outer_legs))) do it + block_sectors = map((g, i) -> blocklabels(g)[i], outer_legs, Tuple(it)) + block_mult = mapreduce((g, i) -> blocklengths(g)[i], *, outer_legs, Tuple(it); init=1) + return build_trees(block_sectors, tree_arrows) .=> block_mult + end +end + +function compute_inner_ranges( + fusion_trees_mult::AbstractVector{<:Pair{<:SectorFusionTree,<:Integer}} +) + fused_leg = blockmergesort( + gradedrange(root_sector.(first.(fusion_trees_mult)) .=> last.(fusion_trees_mult)) + ) + range_mapping = Dict{fieldtype(eltype(fusion_trees_mult), 1),typeof(Block(1)[1:1])}() + fused_sectors = blocklabels(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 to_blockindexrange(b1::BlockIndexRange{1}, b2::BlockIndexRange{1}) + t = (b1, b2) + return Block(Block.(t))[to_block_indices.(t)...] +end + +function intersect_sectors(left::FusedAxes, right::FusedAxes) + return Dict( + map( + t -> first.(t) => to_blockindexrange(last.(t)...), + Iterators.filter( + t -> root_sector(first(t[1])) == root_sector(first(t[2])), + Iterators.product(trees_to_ranges_mapping(left), trees_to_ranges_mapping(right)), + ), + ), + ) +end diff --git a/src/fusiontensor/fusiontensor.jl b/src/fusiontensor/fusiontensor.jl new file mode 100644 index 0000000..8132820 --- /dev/null +++ b/src/fusiontensor/fusiontensor.jl @@ -0,0 +1,192 @@ +# This file defines struct FusionTensor and constructors + +using BlockArrays: AbstractBlockMatrix, BlockArrays, blocklength, findblock + +using BlockSparseArrays: AbstractBlockSparseMatrix, BlockSparseArray, eachblockstoredindex +using GradedUnitRanges: + AbstractGradedUnitRange, + blocklabels, + dual, + isdual, + map_blocklabels, + sector_type, + space_isequal +using SymmetrySectors: SectorProduct, TrivialSector + +struct FusionTensor{T,N,CoDomainAxes,DomainAxes,Mat,Mapping} <: AbstractArray{T,N} + data_matrix::Mat + codomain_axes::CoDomainAxes + domain_axes::DomainAxes + trees_block_mapping::Mapping + + # inner constructor to impose constraints on types + # TBD replace codomain_legs with FusedAxes(codomain_legs)? + function FusionTensor( + mat::AbstractMatrix, + codomain_legs::Tuple{Vararg{AbstractGradedUnitRange}}, + domain_legs::Tuple{Vararg{AbstractGradedUnitRange}}, + trees_block_mapping::Dict, + ) + S = sector_type(axes(mat, 1)) + @assert sector_type(axes(mat, 2)) === S + @assert keytype(trees_block_mapping) <: + Tuple{<:SectorFusionTree{S},<:SectorFusionTree{S}} + @assert all(sector_type.(codomain_legs) .=== S) + @assert all(sector_type.(domain_legs) .=== S) + return new{ + eltype(mat), + length(codomain_legs) + length(domain_legs), + typeof(codomain_legs), + typeof(domain_legs), + typeof(mat), + typeof(trees_block_mapping), + }( + mat, codomain_legs, domain_legs, trees_block_mapping + ) + end +end + +# getters +data_matrix(ft::FusionTensor) = ft.data_matrix +codomain_axes(ft::FusionTensor) = ft.codomain_axes +domain_axes(ft::FusionTensor) = ft.domain_axes +trees_block_mapping(ft::FusionTensor) = ft.trees_block_mapping + +# misc access +ndims_codomain(ft::FusionTensor) = length(codomain_axes(ft)) +ndims_domain(ft::FusionTensor) = length(domain_axes(ft)) + +matrix_size(ft::FusionTensor) = quantum_dimension.(axes(data_matrix(ft))) +matrix_row_axis(ft::FusionTensor) = first(axes(data_matrix(ft))) +matrix_column_axis(ft::FusionTensor) = last(axes(data_matrix(ft))) +function charge_block_size(ft::FusionTensor, f1::SectorFusionTree, f2::SectorFusionTree) + b = Tuple(findblock(ft, f1, f2)) + return ntuple(i -> Int(length(axes(ft)[i][b[i]])), ndims(ft)) +end + +# GradedUnitRanges interface +GradedUnitRanges.sector_type(ft::FusionTensor) = sector_type(matrix_row_axis(ft)) + +# BlockArrays interface +function BlockArrays.findblock(ft::FusionTensor, f1::SectorFusionTree, f2::SectorFusionTree) + # find outer block corresponding to fusion trees + @assert ndims_codomain(ft) == length(f1) + @assert ndims_domain(ft) == length(f2) + @assert sector_type(ft) == sector_type(f1) + @assert sector_type(ft) == sector_type(f2) + b1 = ntuple( + i -> findfirst(==(leaves(f1)[i]), blocklabels(codomain_axes(ft)[i])), ndims_codomain(ft) + ) + b2 = ntuple( + i -> findfirst(==(leaves(f2)[i]), blocklabels(dual(domain_axes(ft)[i]))), + ndims_domain(ft), + ) + return Block(b1..., b2...) +end + +sanitize_axes(::Tuple{}) = () +function sanitize_axes(raw_legs::Tuple{Vararg{AbstractGradedUnitRange}}) + legs = promote_sectors(typeof(first(raw_legs)), raw_legs) + @assert all(check_unique_blocklabels.(legs)) + return legs +end + +function check_unique_blocklabels(g::AbstractGradedUnitRange) + return length(unique(blocklabels(g))) == blocklength(g) +end + +function promote_sectors( + ::Type{<:AbstractGradedUnitRange{LA}}, legs::Tuple{Vararg{AbstractGradedUnitRange{LA}}} +) where {LA} # nothing to do + return legs +end + +function promote_sectors( + ::Type{<:AbstractGradedUnitRange}, legs::Tuple{Vararg{AbstractGradedUnitRange}} +) + T = promote_sector_type(legs) + # fuse with trivial to insert all missing arguments inside each GradedAxis + # avoid depending on SymmetrySectors internals + s0 = trivial(T) + unified_legs = map_blocklabels.(s -> only(blocklabels(fusion_product(s0, s))), legs) + return unified_legs +end + +function promote_sector_type(legs) + # fuse trivial sectors to produce unified type + # avoid depending on SymmetrySectors internals + return sector_type(fusion_product(trivial.(legs)...)) +end + +# initialize with already computed data_matrix +function fusiontensor( + mat::AbstractMatrix, + codomain_legs::Tuple{Vararg{AbstractGradedUnitRange}}, + domain_legs::Tuple{Vararg{AbstractGradedUnitRange}}, +) + # init with empty data_matrix to construct trees_block_mapping + ft = FusionTensor(eltype(mat), codomain_legs, domain_legs) + @assert space_isequal(matrix_row_axis(ft), axes(mat, 1)) + @assert space_isequal(matrix_column_axis(ft), axes(mat, 2)) + for b in eachblockstoredindex(mat) + @assert b in eachblockstoredindex(data_matrix(ft)) # check matrix block is allowed + data_matrix(ft)[b] = mat[b] + end + return ft +end + +# empty matrix +function FusionTensor( + elt::Type, + codomain_legs_raw::Tuple{Vararg{AbstractGradedUnitRange}}, + domain_legs_raw::Tuple{Vararg{AbstractGradedUnitRange}}, +) + legs = sanitize_axes((codomain_legs_raw..., domain_legs_raw...)) + S = sector_type(first(legs)) + codomain_legs = legs[begin:length(codomain_legs_raw)] + domain_legs = legs[(length(codomain_legs_raw) + 1):end] + codomain_fused_axes = FusedAxes{S}(codomain_legs) + domain_fused_axes = FusedAxes{S}(dual.(domain_legs)) + mat = initialize_data_matrix(elt, codomain_fused_axes, domain_fused_axes) + tree_to_block_mapping = intersect_sectors(codomain_fused_axes, domain_fused_axes) + return FusionTensor(mat, codomain_legs, domain_legs, tree_to_block_mapping) +end + +function FusionTensor(elt::Type, ::Tuple{}, ::Tuple{}) + codomain_fused_axes = FusedAxes{TrivialSector}(()) + domain_fused_axes = FusedAxes{TrivialSector}(()) + mat = initialize_data_matrix(elt, codomain_fused_axes, domain_fused_axes) + tree_to_block_mapping = intersect_sectors(codomain_fused_axes, domain_fused_axes) + return FusionTensor(mat, (), (), tree_to_block_mapping) +end + +function initialize_data_matrix( + elt::Type{<:Number}, codomain_fused_axes::FusedAxes, domain_fused_axes::FusedAxes +) + # fusion trees have Float64 eltype: need compatible type + promoted = promote_type(elt, Float64) + mat_row_axis = fused_axis(codomain_fused_axes) + mat_col_axis = dual(fused_axis(domain_fused_axes)) + mat = BlockSparseArray{promoted}(mat_row_axis, mat_col_axis) + initialize_allowed_sectors!(mat) + return mat +end + +function initialize_allowed_sectors!(mat::AbstractMatrix) + row_sectors = blocklabels(axes(mat, 1)) + col_sectors = blocklabels(dual(axes(mat, 2))) + 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 +end + +checkaxes_dual(axes1, axes2) = checkaxes(axes1, dual.(axes2)) +function checkaxes(ax1, ax2) + return checkaxes(Bool, ax1, ax2) || + throw(DimensionMismatch(lazy"$ax1 does not match $ax2")) +end +function checkaxes(::Type{Bool}, axes1, axes2) + return length(axes1) == length(axes2) && all(space_isequal.(axes1, axes2)) +end diff --git a/src/fusiontensor/linear_algebra_interface.jl b/src/fusiontensor/linear_algebra_interface.jl new file mode 100644 index 0000000..7edee18 --- /dev/null +++ b/src/fusiontensor/linear_algebra_interface.jl @@ -0,0 +1,68 @@ +# This file defines linalg for FusionTensor + +using LinearAlgebra: LinearAlgebra, mul!, norm, tr + +using BlockArrays: Block, blocks + +using BlockSparseArrays: eachblockstoredindex +using GradedUnitRanges: blocklabels +using SymmetrySectors: quantum_dimension + +# allow to contract with different eltype and let BlockSparseArray ensure compatibility +# impose matching type and number of axes at compile time +# impose matching axes at run time +# TODO remove this once TensorAlgebra.contract can be used? +function LinearAlgebra.mul!( + C::FusionTensor, A::FusionTensor, B::FusionTensor, α::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 + checkaxes_dual(domain_axes(A), codomain_axes(B)) + checkaxes(codomain_axes(C), codomain_axes(A)) + checkaxes(domain_axes(C), domain_axes(B)) + mul!(data_matrix(C), data_matrix(A), data_matrix(B), α, β) + return C +end + +function LinearAlgebra.norm(ft::FusionTensor) + m = data_matrix(ft) + row_sectors = blocklabels(matrix_row_axis(ft)) + n2 = sum(eachblockstoredindex(m); init=zero(real(eltype(ft)))) do b + return quantum_dimension(row_sectors[Int(first(Tuple(b)))]) * norm(m[b])^2 + end + return sqrt(n2) +end + +function LinearAlgebra.tr(ft::FusionTensor) + m = data_matrix(ft) + row_sectors = blocklabels(matrix_row_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 + +function LinearAlgebra.qr(ft::FusionTensor) + qmat, rmat = block_qr(data_matrix(ft)) + qtens = FusionTensor(qmat, codomain_axes(ft), (axes(qmat, 2),)) + rtens = FusionTensor(rmat, (axes(rmat, 1),), domain_axes(ft)) + return qtens, rtens +end + +function LinearAlgebra.svd(ft::FusionTensor) + umat, s, vmat = block_svd(data_matrix(ft)) + utens = FusionTensor(umat, codomain_axes(ft), (axes(umat, 2),)) + stens = FusionTensor(s, (axes(umat, 1),), (axes(vmat, 2),)) + vtens = FusionTensor(vmat, (axes(vmat, 1),), domain_axes(ft)) + return utens, stens, vtens +end diff --git a/src/fusiontensor/tensor_algebra_interface.jl b/src/fusiontensor/tensor_algebra_interface.jl new file mode 100644 index 0000000..cba822c --- /dev/null +++ b/src/fusiontensor/tensor_algebra_interface.jl @@ -0,0 +1,48 @@ +# This file defines TensorAlgebra interface for a FusionTensor + +using LinearAlgebra: mul! + +using BlockArrays: Block + +using TensorAlgebra: BlockedPermutation, Matricize, TensorAlgebra + +# TBD how to deal with inner contraction = no ouput axis? +function TensorAlgebra.allocate_output( + ::typeof(contract), + biperm_dest::BlockedPermutation{2}, + a1::FusionTensor{T1,N}, + biperm1::BlockedPermutation{2,N}, + a2::FusionTensor{T2,M}, + biperm2::BlockedPermutation{2,M}, + α::Number=true, +) where {T1,T2,N,M} + axes_dest = ( + (i -> axes(a1)[i]).(biperm1[Block(1)]), (i -> axes(a2)[i]).(biperm2[Block(2)]) + ) + return similar(a1, promote_type(eltype(a1), eltype(a2), typeof(α)), axes_dest) +end + +# TBD do really I need to defined these as I cannot use them in contract! and has to redefine it? +# TensorAlgebra.fusedims(ft::FusionTensor, perm::BlockedPermutation) = permutedims(ft, perm) +# function TensorAlgebra.splitdims(ft1::FusionTensor, ft2::FusionTensor, blockedperm::BlockedPermutation) +# function TensorAlgebra.splitdims!(ft1::FusionTensor, ft2::FusionTensor, blockedperm::BlockedPermutation) + +# I cannot use contract! from TensorAlgebra/src/contract/contract_matricize/contract.jl +# as it calls _mul!, which I should not overload. +# TBD I can also overload higher up and do not allow use of different algorithms +function TensorAlgebra.contract!( + alg::Matricize, + a_dest::FusionTensor, + biperm_dest::BlockedPermutation, + a1::FusionTensor, + biperm1::BlockedPermutation, + a2::FusionTensor, + biperm2::BlockedPermutation, + α::Number, + β::Number, +) + a1_perm = permutedims(a1, biperm1) + a2_perm = permutedims(a2, biperm2) + mul!(a_dest, a1_perm, a2_perm, α, β) + return a_dest +end diff --git a/src/permutedims/permutedims.jl b/src/permutedims/permutedims.jl new file mode 100644 index 0000000..2bc7031 --- /dev/null +++ b/src/permutedims/permutedims.jl @@ -0,0 +1,58 @@ +# This file defines permutedims for a FusionTensor + +using BlockArrays: blocklengths +using Strided: Strided, @strided + +using TensorAlgebra: BlockedPermutation, blockedperm, blockpermute + +function naive_permutedims(ft, biperm::BlockedPermutation{2}) + @assert ndims(ft) == length(biperm) + new_codomain_legs, new_domain_legs = blockpermute(axes(ft), biperm) + + # naive permute: cast to dense, permutedims, cast to FusionTensor + arr = Array(ft) + permuted_arr = permutedims(arr, Tuple(biperm)) + permuted = FusionTensor(permuted_arr, new_codomain_legs, new_domain_legs) + return permuted +end + +# permutedims with 1 tuple of 2 separate tuples +function fusiontensor_permutedims(ft, new_leg_indices::Tuple{Tuple,Tuple}) + return fusiontensor_permutedims(ft, new_leg_indices...) +end + +# permutedims with 2 separate tuples +function fusiontensor_permutedims( + ft, new_codomain_indices::Tuple, new_domain_indices::Tuple +) + biperm = blockedperm(new_codomain_indices, new_domain_indices) + return fusiontensor_permutedims(ft, biperm) +end + +function fusiontensor_permutedims(ft, biperm::BlockedPermutation{2}) + @assert ndims(ft) == length(biperm) + + # early return for identity operation. Do not copy. Also handle tricky 0-dim case. + if ndims_codomain(ft) == first(blocklengths(biperm)) # compile time + if Tuple(biperm) == ntuple(identity, ndims(ft)) + return ft + end + end + + new_codomain_legs, new_domain_legs = blockpermute(axes(ft), biperm) + new_ft = FusionTensor(eltype(ft), new_codomain_legs, new_domain_legs) + fusiontensor_permutedims!(new_ft, ft, Tuple(biperm)) + return new_ft +end + +function fusiontensor_permutedims!( + new_ft::FusionTensor{T,N}, old_ft::FusionTensor{T,N}, flatperm::NTuple{N,Integer} +) where {T,N} + unitary = compute_unitary(new_ft, old_ft, flatperm) + for p in unitary + old_trees, new_trees = first(p) + new_block = view(new_ft, new_trees...) + old_block = view(old_ft, old_trees...) + @strided new_block .+= last(p) .* permutedims(old_block, flatperm) + end +end diff --git a/src/permutedims/unitaries.jl b/src/permutedims/unitaries.jl new file mode 100644 index 0000000..5e9dbf0 --- /dev/null +++ b/src/permutedims/unitaries.jl @@ -0,0 +1,64 @@ +# This file defines unitaries to be used in permutedims + +using BlockArrays: Block, findblock +using LRUCache: LRU + +using SymmetrySectors: quantum_dimension + +const unitary_cache = LRU{ + Tuple{ + SectorFusionTree,SectorFusionTree,SectorFusionTree,SectorFusionTree,Tuple{Vararg{Int}} + }, + 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) +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))) +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 +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) + end + end + return unitary +end + +# ================================= Constructor from 6j ================================== +# dummy diff --git a/test/Project.toml b/test/Project.toml index e9c291e..110a927 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,7 +1,14 @@ [deps] Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" +BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e" +BlockSparseArrays = "2c9a651f-6452-4ace-a6ac-809f4280fbb4" +FusionTensors = "e16ca583-1f51-4df0-8e12-57d32947d33e" +GradedUnitRanges = "e2de450a-8a67-46c7-b59c-01d5a3d041c5" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb" +SymmetrySectors = "f8a8ad64-adbc-4fce-92f7-ffe2bb36a86e" +TensorAlgebra = "68bd88dc-f39d-4e12-b2ca-f046b68fcc6a" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] diff --git a/test/basics/setup.jl b/test/basics/setup.jl new file mode 100644 index 0000000..70a0025 --- /dev/null +++ b/test/basics/setup.jl @@ -0,0 +1,56 @@ +using LinearAlgebra: Adjoint + +using BlockSparseArrays: BlockSparseMatrix, eachblockstoredindex +using FusionTensors: + FusionTensor, + codomain_axes, + data_matrix, + domain_axes, + checkaxes, + checkaxes_dual, + matrix_column_axis, + matrix_row_axis, + ndims_codomain, + ndims_domain +using GradedUnitRanges: blocklabels, dual, space_isequal + +function check_data_matrix_axes( + mat::BlockSparseMatrix, codomain_legs::Tuple, domain_legs::Tuple +) + ft0 = FusionTensor(Float64, codomain_legs, domain_legs) + @assert space_isequal(matrix_row_axis(ft0), axes(mat, 1)) + @assert space_isequal(matrix_column_axis(ft0), axes(mat, 2)) +end + +function check_data_matrix_axes(mat::Adjoint, codomain_legs::Tuple, domain_legs::Tuple) + return check_data_matrix_axes(adjoint(mat), dual.(domain_legs), dual.(codomain_legs)) +end + +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" + + 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" + checkaxes(axes(ft)[begin:nda], codomain_axes(ft)) + checkaxes(axes(ft)[(nda + 1):end], domain_axes(ft)) + + m = data_matrix(ft) + @assert ndims(m) == 2 "invalid data_matrix ndims" + row_axis = matrix_row_axis(ft) + column_axis = matrix_column_axis(ft) + @assert row_axis === axes(m, 1) "invalid row_axis" + @assert column_axis === axes(m, 2) "invalid column_axis" + check_data_matrix_axes(m, codomain_axes(ft), domain_axes(ft)) + + for b in eachblockstoredindex(m) + ir, ic = Int.(Tuple(b)) + @assert blocklabels(row_axis)[ir] == blocklabels(dual(column_axis))[ic] "forbidden block" + end + return nothing +end diff --git a/test/basics/test_array_cast.jl b/test/basics/test_array_cast.jl new file mode 100644 index 0000000..ea34b3e --- /dev/null +++ b/test/basics/test_array_cast.jl @@ -0,0 +1,415 @@ +@eval module $(gensym()) +using LinearAlgebra: LinearAlgebra, norm +using Test: @test, @test_broken, @test_throws, @testset + +using BlockArrays: Block, BlockedArray, blocksize + +using FusionTensors: FusionTensor, data_matrix +using GradedUnitRanges: dual, fusion_product, gradedrange +using SymmetrySectors: O2, SectorProduct, SU2, TrivialSector, U1 + +include("setup.jl") + +@testset "Trivial FusionTensor" begin + @testset "trivial matrix" begin + g = gradedrange([TrivialSector() => 1]) + gb = dual(g) + m = ones((1, 1)) + ft = 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 = FusionTensor(m, (g,), (gb,)) + @test eltype(ft) === Float64 + @test Array(ft) ≈ m + end + + for elt in (ComplexF32, ComplexF64) + m = ones(elt, (1, 1)) + ft = FusionTensor(m, (g,), (gb,)) + @test eltype(ft) === ComplexF64 + @test Array(ft) ≈ m + 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 = 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 = 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 = 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 = 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 FusionTensor( + dense, (gradedrange([U1(1) => 1, U1(2) => 3]),), domain_legs + ) + @test_throws MethodError FusionTensor(dense, (g, g), domain_legs) + + ba = BlockedArray(dense, [1, 2], [1, 2]) + @test_throws DomainError FusionTensor( + ba, (gradedrange([U1(1) => 1, U1(2) => 3]),), domain_legs + ) + @test_throws DomainError FusionTensor(ba, (g, g), domain_legs) + dense[1, 2] = 1 # forbidden + @test_throws InexactError 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 = 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 = 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 = 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 = 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 = 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 = FusionTensor(v, (dual(g),), ()) + @test isnothing(check_sanity(ft3)) + @test Array(ft3) ≈ v + @test Array(adjoint(ft3)) ≈ v + + ft4 = FusionTensor(v, (), (g,)) + @test isnothing(check_sanity(ft4)) + @test Array(ft4) ≈ v + @test Array(adjoint(ft4)) ≈ v + + zerodim = ones(()) + if VERSION < v"1.11" + @test_broken FusionTensor(zerodim, (), ()) isa FusionTensor # https://github.com/JuliaLang/julia/issues/52615 + else + # TODO fix: add specialized method, maybe fix TensorAlgebra + @test_broken 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_broken Array(ft) ≈ zerodim # cannot create zerodim BlockSparseArray + 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 = 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 = 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 = 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 = 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 = 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 = 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 = 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 = 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 = 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 = 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 = 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 = 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(fusion_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 = 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 + 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 = 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 = FusionTensor(tRVB, codomain_legs_nt, domain_legs_nt) + @test isnothing(check_sanity(ft_nt)) + @test Array(ft_nt) ≈ tRVB + end +end +end diff --git a/test/basics/test_basics.jl b/test/basics/test_basics.jl index a6d9b13..f644f9c 100644 --- a/test/basics/test_basics.jl +++ b/test/basics/test_basics.jl @@ -1,6 +1,267 @@ -using FusionTensors: FusionTensors -using Test: @test, @testset +using Test: @test, @test_throws, @testset -@testset "FusionTensors" begin - # Tests go here. +using BlockSparseArrays: BlockSparseArray +using FusionTensors: + FusionTensor, + codomain_axes, + data_matrix, + domain_axes, + fusiontensor, + checkaxes, + checkaxes_dual, + matrix_column_axis, + matrix_row_axis, + matrix_size, + ndims_domain, + ndims_codomain +using GradedUnitRanges: + blockmergesort, dual, flip, fusion_product, gradedrange, sector_type, space_isequal +using SymmetrySectors: U1, SU2, SectorProduct, Z + +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])) + + # check dual convention when initializing data_matrix + ft0 = FusionTensor(Float64, (g1,), (g2,)) + @test space_isequal(matrix_row_axis(ft0), g1) + @test space_isequal(matrix_column_axis(ft0), g2) + + m = BlockSparseArray{Float64}(g1, g2) + ft1 = fusiontensor(m, (g1,), (g2,)) + + # getters + @test data_matrix(ft1) == m + @test checkaxes(codomain_axes(ft1), (g1,)) + @test checkaxes(domain_axes(ft1), (g2,)) + + # misc + @test checkaxes(axes(ft1), (g1, g2)) + @test ndims_codomain(ft1) == 1 + @test ndims_domain(ft1) == 1 + @test matrix_size(ft1) == (6, 5) + @test space_isequal(matrix_row_axis(ft1), g1) + @test space_isequal(matrix_column_axis(ft1), g2) + @test isnothing(check_sanity(ft0)) + @test isnothing(check_sanity(ft1)) + @test sector_type(ft1) === U1{Int} + + # Base methods + @test eltype(ft1) === Float64 + @test length(ft1) == 30 + @test ndims(ft1) == 2 + @test size(ft1) == (6, 5) + + # 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 checkaxes(codomain_axes(ft2), codomain_axes(ft1)) + @test checkaxes(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 checkaxes(codomain_axes(ft2), codomain_axes(ft1)) + @test checkaxes(domain_axes(ft2), domain_axes(ft1)) + + # similar + ft2 = similar(ft1) + @test isnothing(check_sanity(ft2)) + @test eltype(ft2) == Float64 + @test checkaxes(codomain_axes(ft2), codomain_axes(ft1)) + @test checkaxes(domain_axes(ft2), domain_axes(ft1)) + + ft3 = similar(ft1, ComplexF64) + @test isnothing(check_sanity(ft3)) + @test eltype(ft3) == ComplexF64 + @test checkaxes(codomain_axes(ft3), codomain_axes(ft1)) + @test checkaxes(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 checkaxes(codomain_axes(ft5), (g1, g1)) + @test checkaxes(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 = fusion_product(g1, g2) + gc = dual(fusion_product(dual(g3), dual(g4))) + m2 = BlockSparseArray{Float64}(gr, gc) + ft = fusiontensor(m2, (g1, g2), (g3, g4)) + + @test data_matrix(ft) == m2 + @test checkaxes(codomain_axes(ft), (g1, g2)) + @test checkaxes(domain_axes(ft), (g3, g4)) + + @test axes(ft) == (g1, g2, g3, g4) + @test ndims_codomain(ft) == 2 + @test ndims_domain(ft) == 2 + @test matrix_size(ft) == (30, 12) + @test space_isequal(matrix_row_axis(ft), gr) + @test space_isequal(matrix_column_axis(ft), gc) + @test isnothing(check_sanity(ft)) + + @test ndims(ft) == 4 + @test size(ft) == (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, (g1,), ()) + @test ndims_codomain(ft1) == 1 + @test ndims_domain(ft1) == 0 + @test ndims(ft1) == 1 + @test size(ft1) == (6,) + @test size(data_matrix(ft1)) == (6, 1) + @test isnothing(check_sanity(ft1)) + + # one column axis + ft2 = FusionTensor(Float64, (), (g1,)) + @test ndims_codomain(ft2) == 0 + @test ndims_domain(ft2) == 1 + @test ndims(ft2) == 1 + @test size(ft2) == (6,) + @test size(data_matrix(ft2)) == (1, 6) + @test isnothing(check_sanity(ft2)) + + # zero axis + ft3 = FusionTensor(Float64, (), ()) + @test ndims_codomain(ft3) == 0 + @test ndims_domain(ft3) == 0 + @test ndims(ft3) == 0 + @test size(ft3) == () + @test size(data_matrix(ft3)) == (1, 1) + @test isnothing(check_sanity(ft3)) +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 = FusionTensor(Float64, (g1, g2), (g3, g4)) + @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) + + ft4 = ft3 + ft3 + @test codomain_axes(ft4) === codomain_axes(ft3) + @test domain_axes(ft4) === domain_axes(ft3) + @test space_isequal(matrix_row_axis(ft4), matrix_row_axis(ft3)) + @test space_isequal(matrix_column_axis(ft4), matrix_column_axis(ft3)) + @test isnothing(check_sanity(ft4)) + + ft4 = ft3 - ft3 + @test codomain_axes(ft4) === codomain_axes(ft3) + @test domain_axes(ft4) === domain_axes(ft3) + @test space_isequal(matrix_row_axis(ft4), matrix_row_axis(ft3)) + @test space_isequal(matrix_column_axis(ft4), matrix_column_axis(ft3)) + @test isnothing(check_sanity(ft4)) + + ft4 = 2 * ft3 + @test codomain_axes(ft4) === codomain_axes(ft3) + @test domain_axes(ft4) === domain_axes(ft3) + @test space_isequal(matrix_row_axis(ft4), matrix_row_axis(ft3)) + @test space_isequal(matrix_column_axis(ft4), matrix_column_axis(ft3)) + @test isnothing(check_sanity(ft4)) + @test eltype(ft4) == Float64 + + ft4 = 2.0 * ft3 + @test codomain_axes(ft4) === codomain_axes(ft3) + @test domain_axes(ft4) === domain_axes(ft3) + @test space_isequal(matrix_row_axis(ft4), matrix_row_axis(ft3)) + @test space_isequal(matrix_column_axis(ft4), matrix_column_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(matrix_row_axis(ft4), matrix_row_axis(ft3)) + @test space_isequal(matrix_column_axis(ft4), matrix_column_axis(ft3)) + @test isnothing(check_sanity(ft4)) + @test eltype(ft4) == Float64 + + ft5 = 2.0im * ft3 + @test codomain_axes(ft5) === codomain_axes(ft3) + @test domain_axes(ft5) === domain_axes(ft3) + @test space_isequal(matrix_row_axis(ft5), matrix_row_axis(ft3)) + @test space_isequal(matrix_column_axis(ft5), matrix_column_axis(ft3)) + @test isnothing(check_sanity(ft4)) + @test eltype(ft5) == ComplexF64 + + ft4 = conj(ft3) + @test ft4 === ft3 # same object + + 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(matrix_row_axis(ft6), matrix_row_axis(ft5)) + @test space_isequal(matrix_column_axis(ft6), matrix_column_axis(ft5)) + @test eltype(ft6) == ComplexF64 + + 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, (g1,), (g2, g3, g4)) + @test_throws DimensionMismatch ft7 + ft3 + @test_throws DimensionMismatch ft7 - ft3 + @test_throws DimensionMismatch ft7 * ft3 +end + +@testset "mising 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, (g1,), (dual(g2), dual(g3))) + @test sector_type(ft) === S + gr = gradedrange([SectorProduct(U1(1), SU2(0), Z{2}(0)) => 1]) + @test space_isequal(matrix_row_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(matrix_column_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 = fusion_product(fusion_product(gA, gB), gC) + S = sector_type(gABC) + + ft = FusionTensor(Float64, (gA, gB), (dual(gA), dual(gB), gC)) + @test sector_type(ft) === S + @test space_isequal(matrix_row_axis(ft), gABC) + @test space_isequal(matrix_column_axis(ft), dual(gABC)) end diff --git a/test/basics/test_contraction.jl b/test/basics/test_contraction.jl new file mode 100644 index 0000000..3342a8d --- /dev/null +++ b/test/basics/test_contraction.jl @@ -0,0 +1,68 @@ +@eval module $(gensym()) +using LinearAlgebra: mul! +using Test: @test, @testset, @test_broken + +using BlockSparseArrays: BlockSparseArray +using FusionTensors: FusionTensor, domain_axes, codomain_axes +using GradedUnitRanges: dual, gradedrange +using SymmetrySectors: U1 +using TensorAlgebra: contract + +include("setup.jl") + +@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]) + + ft1 = FusionTensor(Float64, (g1, g2), (g3, g4)) + @test isnothing(check_sanity(ft1)) + + ft2 = FusionTensor(Float64, 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) + + # test LinearAlgebra.mul! with in-place matrix product + mul!(ft3, ft1, ft2) + @test isnothing(check_sanity(ft3)) + @test domain_axes(ft3) === domain_axes(ft2) + @test codomain_axes(ft3) === codomain_axes(ft1) + + mul!(ft3, ft1, ft2, 1.0, 1.0) + @test isnothing(check_sanity(ft2)) + @test domain_axes(ft3) === domain_axes(ft2) + @test codomain_axes(ft3) === codomain_axes(ft1) +end + +@testset "TensorAlgebra 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 = FusionTensor(Float64, (g1, g2), (g3, g4)) + ft2 = FusionTensor(Float64, dual.((g3, g4)), (dual(g1),)) + ft3 = FusionTensor(Float64, dual.((g3, g4)), dual.((g1, g2))) + + ft4, legs = contract(ft1, (1, 2, 3, 4), ft2, (3, 4, 5)) + @test legs == (1, 2, 5) + @test isnothing(check_sanity(ft4)) + @test domain_axes(ft4) === domain_axes(ft2) + @test codomain_axes(ft4) === codomain_axes(ft1) + + ft5 = contract((1, 2, 5), ft1, (1, 2, 3, 4), ft2, (3, 4, 5)) + @test isnothing(check_sanity(ft5)) + + # biperm is not allowed + @test_broken contract(((1, 2), (5,)), ft1, (1, 2, 3, 4), ft2, (3, 4, 5)) + + @test permutedims(ft1, (), (1, 2, 3, 4)) * permutedims(ft3, (3, 4, 1, 2), ()) isa + FusionTensor{Float64,0} + @test_broken contract(ft1, (1, 2, 3, 4), ft3, (3, 4, 1, 2)) +end +end diff --git a/test/basics/test_fusion_trees.jl b/test/basics/test_fusion_trees.jl new file mode 100644 index 0000000..9043e0e --- /dev/null +++ b/test/basics/test_fusion_trees.jl @@ -0,0 +1,127 @@ +@eval module $(gensym()) +using Test: @test, @testset +using LinearAlgebra: I + +using BlockArrays: BlockArrays + +using FusionTensors: + SectorFusionTree, + arrows, + branch_sectors, + build_trees, + leaves, + outer_multiplicity_indices, + root_sector +using GradedUnitRanges: blocklabels, fusion_product, sector_type +using SymmetrySectors: ×, SectorProduct, SU, SU2, TrivialSector, arguments + +@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)) +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) + + 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 + 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 +end +end diff --git a/test/basics/test_linear_algebra.jl b/test/basics/test_linear_algebra.jl new file mode 100644 index 0000000..6fc33b0 --- /dev/null +++ b/test/basics/test_linear_algebra.jl @@ -0,0 +1,39 @@ +@eval module $(gensym()) +using LinearAlgebra: norm, tr +using Test: @test, @testset + +using BlockArrays: BlockArrays + +using BlockSparseArrays: BlockSparseArrays +using FusionTensors: FusionTensor +using GradedUnitRanges: dual, gradedrange +using SymmetrySectors: U1, SU2, TrivialSector + +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]) + gsu2 = gradedrange([SU2(1 / 2) => 1]) + + for g in [g0, gu1, gsu2] + ft0 = FusionTensor(Float64, (g, g), (dual(g), dual(g))) + @test isnothing(check_sanity(ft0)) + @test norm(ft0) == 0 + @test tr(ft0) == 0 + + ft = FusionTensor(sdst, (g, g), (dual(g), dual(g))) + @test isnothing(check_sanity(ft)) + @test norm(ft) ≈ √3 / 2 + @test isapprox(tr(ft), 0; atol=eps(Float64)) + end +end +end diff --git a/test/basics/test_permutedims.jl b/test/basics/test_permutedims.jl new file mode 100644 index 0000000..12c7695 --- /dev/null +++ b/test/basics/test_permutedims.jl @@ -0,0 +1,159 @@ +@eval module $(gensym()) +using Test: @test, @testset, @test_broken + +using FusionTensors: + FusionTensor, + data_matrix, + checkaxes, + matrix_column_axis, + matrix_row_axis, + naive_permutedims, + ndims_domain, + ndims_codomain +using GradedUnitRanges: dual, gradedrange, space_isequal +using SymmetrySectors: O2, U1, SectorProduct, SU2 +using TensorAlgebra: blockedperm + +include("setup.jl") + +@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]) + + for elt in (Float64, ComplexF64) + ft1 = FusionTensor(elt, (g1, g2), dual.((g3, g4))) + @test isnothing(check_sanity(ft1)) + + # test permutedims interface + ft2 = permutedims(ft1, (1, 2), (3, 4)) # trivial with 2 tuples + @test ft2 === ft1 # same object + + ft2 = permutedims(ft1, ((1, 2), (3, 4))) # trivial with tuple of 2 tuples + @test ft2 === ft1 # same object + + biperm = blockedperm((1, 2), (3, 4)) + ft2 = permutedims(ft1, biperm) # trivial with biperm + @test ft2 === ft1 # same object + + ft3 = permutedims(ft1, (4,), (1, 2, 3)) + @test ft3 !== ft1 + @test ft3 isa FusionTensor{elt,4} + @test checkaxes(axes(ft3), (dual(g4), g1, g2, dual(g3))) + @test ndims_domain(ft3) == 3 + @test ndims_codomain(ft3) == 1 + @test ndims(ft3) == 4 + @test isnothing(check_sanity(ft3)) + + ft4 = permutedims(ft3, (2, 3), (4, 1)) + @test checkaxes(axes(ft1), axes(ft4)) + @test space_isequal(matrix_column_axis(ft1), matrix_column_axis(ft4)) + @test space_isequal(matrix_row_axis(ft1), matrix_row_axis(ft4)) + @test ft4 ≈ ft1 + 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 = FusionTensor(arr, codomain_legs, domain_legs) + biperm = blockedperm((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" + @test_broken FusionTensor(ones(()), (), ()) isa FusionTensor + #= 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, blockedperm((), ())) isa FusionTensor{Float64,0} + =# + end + + g = gradedrange([U1(0) => 1, U1(1) => 2, U1(2) => 3]) + v = zeros((6,)) + v[1] = 1.0 + biperm = blockedperm((), (1,)) + ft1 = 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])), + gradedrange([SectorProduct(SU2(1//2), U1(0)) => 1]), + gradedrange([SectorProduct(SU2(1//2), SU2(0)) => 1]), + ) + g2b = dual(g2) + for biperm in [ + blockedperm((2, 1), (3, 4)), blockedperm((3, 1), (2, 4)), blockedperm((3, 1, 4), (2,)) + ] + ft = 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 = 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 [blockedperm((1, 2, 3, 4), ()), blockedperm((), (3, 1, 2, 4))] + ft = FusionTensor(sds22, (g2, g2), (g2b, g2b)) + @test permutedims(ft, biperm) ≈ naive_permutedims(ft, biperm) + end + end +end +end