[WIP] GradedArraysNext

Project.toml

@@ -8,11 +8,14 @@ ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
 BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
 BlockSparseArrays = "2c9a651f-6452-4ace-a6ac-809f4280fbb4"
 Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"
+DiagonalArrays = "74fd4be6-21e2-4f6f-823a-4360d37c7a77"
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 HalfIntegers = "f0d1745a-41c9-11e9-1dd9-e5d34d218721"
+KroneckerArrays = "05d0b138-81bc-4ff7-84be-08becefb1ccc"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MatrixAlgebraKit = "6c742aac-3347-4629-af66-fc926824e5e4"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+SparseArraysBase = "0d5efcca-f356-4864-8770-e1ed8d78f208"
 SplitApplyCombine = "03a91e81-4c3e-53e1-a0a4-9c0c8f19dd66"
 TensorAlgebra = "68bd88dc-f39d-4e12-b2ca-f046b68fcc6a"
 TensorKitSectors = "13a9c161-d5da-41f0-bcbd-e1a08ae0647f"
@@ -30,15 +33,18 @@ ArrayLayouts = "1"
 BlockArrays = "1.6"
 BlockSparseArrays = "0.10"
 Compat = "4.16"
+DiagonalArrays = "0.3.23"
 FillArrays = "1.13"
 HalfIntegers = "1.6"
+KroneckerArrays = "0.3"
 LinearAlgebra = "1.10"
 MatrixAlgebraKit = "0.2, 0.3, 0.4, 0.5"
 Random = "1.10"
 SUNRepresentations = "0.3"
+SparseArraysBase = "0.7.7"
 SplitApplyCombine = "1.2.3"
 TensorAlgebra = "0.3.2, 0.4"
 TensorKitSectors = "0.1, 0.2"
-TensorProducts = "0.1.3"
+TensorProducts = "0.1.8"
 TypeParameterAccessors = "0.4"
 julia = "1.10"

src/GradedArrays.jl

@@ -1,32 +1,37 @@
 module GradedArrays
 
-include("gradedunitrange_interface.jl")
+# exports
+# -------
+export TrivialSector, Z, Z2, U1, O2, SU2, Fib, Ising
+export SectorRange, SectorDelta
+export SectorUnitRange, SectorOneTo, SectorArray, SectorMatrix
+export GradedUnitRange, GradedOneTo, GradedArray
+
+export sector_type
+
+export dag, dual, flip, gradedrange, isdual,
+    sector, sector_multiplicities, sector_multiplicity,
+    sectorrange, sectors, sector_type,
+    space_isequal, ungrade
+
+# imports
+# -------
+using LinearAlgebra: LinearAlgebra, Adjoint
+using KroneckerArrays
+using KroneckerArrays: AbstractKroneckerArray, CartesianProductUnitRange
+using SparseArraysBase: SparseArraysBase, isstored
+using BlockArrays: BlockArrays, Block
+using BlockSparseArrays: AbstractBlockSparseArray, blockrange, @view!
 
 include("sectorrange.jl")
-include("sectorunitrange.jl")
-include("gradedunitrange.jl")
+include("sectorarray.jl")
+include("gradedarray.jl")
 
-include("namedtuple_operations.jl")
 include("sectorproduct.jl")
 
 include("fusion.jl")
-include("gradedarray.jl")
 include("tensoralgebra.jl")
 include("factorizations.jl")
 
-export TrivialSector, Z, Z2, U1, O2, SU2, Fib, Ising
-export dag,
-    dual,
-    flip,
-    gradedrange,
-    isdual,
-    sector,
-    sector_multiplicities,
-    sector_multiplicity,
-    sectorrange,
-    sectors,
-    sector_type,
-    space_isequal,
-    ungrade
 
 end

src/fusion.jl

@@ -1,4 +1,4 @@
-using BlockArrays: Block, blocks
+using BlockArrays: Block, blocks, BlockVector
 using SplitApplyCombine: groupcount
 using TensorProducts: TensorProducts, ⊗, OneToOne, tensor_product
 
@@ -13,28 +13,27 @@ function TensorProducts.tensor_product(
         ::AbelianStyle, sr1::SectorUnitRange, sr2::SectorUnitRange
     )
     s = sector(flip_dual(sr1)) ⊗ sector(flip_dual(sr2))
-    return sectorrange(s, sector_multiplicity(sr1) * sector_multiplicity(sr2))
+    return cartesianrange(s, sector_multiplicity(sr1) * sector_multiplicity(sr2))
 end
 
 function TensorProducts.tensor_product(
         ::NotAbelianStyle, sr1::SectorUnitRange, sr2::SectorUnitRange
     )
-    g0 = sector(flip_dual(sr1)) ⊗ sector(flip_dual(sr2))
-    return gradedrange(
-        sectors(g0) .=>
-            sector_multiplicity(sr1) * sector_multiplicity(sr2) .* sector_multiplicities(g0),
-    )
+    g = sector(flip_dual(sr1)) ⊗ sector(flip_dual(sr2))
+    d₁ = sector_multiplicity(sr1)
+    d₂ = sector_multiplicity(sr2)
+    return blockrange([c × (d₁ * d₂ * d) for (c, d) in zip(sectors(g), sector_multiplicities(g))])
 end
 
 # allow to fuse a Sector with a GradedUnitRange
 function TensorProducts.tensor_product(
-        s::Union{SectorRange, SectorUnitRange}, g::AbstractGradedUnitRange
+        s::Union{SectorRange, SectorUnitRange}, g::GradedUnitRange
     )
     return to_gradedrange(s) ⊗ g
 end
 
 function TensorProducts.tensor_product(
-        g::AbstractGradedUnitRange, s::Union{SectorRange, SectorUnitRange}
+        g::GradedUnitRange, s::Union{SectorRange, SectorUnitRange}
     )
     return g ⊗ to_gradedrange(s)
 end
@@ -62,7 +61,7 @@ end
 # default to tensor_product
 unmerged_tensor_product(a1, a2) = a1 ⊗ a2
 
-function unmerged_tensor_product(a1::AbstractGradedUnitRange, a2::AbstractGradedUnitRange)
+function unmerged_tensor_product(a1::GradedUnitRange, a2::GradedUnitRange)
     new_axes = map(splat(⊗), Iterators.flatten((Iterators.product(blocks(a1), blocks(a2)),)))
     return mortar_axis(new_axes)
 end
@@ -91,22 +90,23 @@ end
 # Used by `TensorAlgebra.unmatricize` in `GradedArraysTensorAlgebraExt`.
 invblockperm(a::Vector{<:Block{1}}) = Block.(invperm(Int.(a)))
 
-function sectormergesort(g::AbstractGradedUnitRange)
+function sectormergesort(g::GradedUnitRange)
     glabels = sectors(g)
     multiplicities = sector_multiplicities(g)
-    new_blocklengths = map(sort(unique(glabels))) do la
-        return la => sum(multiplicities[findall(==(la), glabels)]; init = 0)
+    dict = Dict{eltype(glabels), eltype(multiplicities)}()
+    for (l, m) in zip(glabels, multiplicities)
+        dict[l] = get(dict, l, 0) + m
     end
-    return gradedrange(new_blocklengths)
+
+    total = sort!(collect(pairs(dict)); by = first)
+    return blockrange([c × m for (c, m) in total])
 end
 
 sectormergesort(g::AbstractUnitRange) = g
 
 # tensor_product produces a sorted, non-dual GradedUnitRange
-TensorProducts.tensor_product(g::AbstractGradedUnitRange) = sectormergesort(flip_dual(g))
+TensorProducts.tensor_product(g::GradedUnitRange) = sectormergesort(flip_dual(g))
 
-function TensorProducts.tensor_product(
-        g1::AbstractGradedUnitRange, g2::AbstractGradedUnitRange
-    )
+function TensorProducts.tensor_product(g1::GradedUnitRange, g2::GradedUnitRange)
     return sectormergesort(unmerged_tensor_product(g1, g2))
 end