First draft of FusionTensor #5

ogauthe · 2024-12-19T14:12:01Z

This PR provides a first draft for FusionTensors. It defines a generic interface to implement non-abelian symmetries. It currently supports O(2), SU(2), any abelian group defined in SymmetrySectors. It allows to construct, permute and contract tensors.

Work to be done in future PR:

remove FusedAxes, currently only used in FusionTensor initialization
handle dual in a more explicit way when fusing codomain and domain fusion trees
store unitary as a sparse matrix
improve caching of structural data
maybe improve detection of a charge block size
maybe moving Clebsch-Gordan tensor definition to SymmetrySectors

Work depending on other packages

implement linear algebra decomposition such as SVD (depends on [Enhancement] Add singular value decomposition BlockSparseArrays.jl#16)
rewrite TensorAlgebra interface once contract preserves domain and codomain information
move nsymbol to SymmetrySectors (define nsymbol SymmetrySectors.jl#4)

ogauthe · 2024-12-19T14:17:59Z

I have issues with building test both on local and on GitHub actions. Locally, I get an error

  LoadError: ArgumentError: Package LinearAlgebra not found in current path.

although the package is included in the Project.toml. I guess there was an issue with SafeTestsets configuration when I merged the new skeleton. @lkdvos do you have a clue?

src/fusiontensor/array_cast.jl

mtfishman · 2024-12-19T14:28:59Z

src/fusiontensor/base_interface.jl

+  return FusionTensor(
+    x * data_matrix(ft), codomain_axes(ft), domain_axes(ft), trees_block_mapping(ft)
+  )


This and other functions can be simplified with the follow "field setting" code pattern:

using Accessors: @set 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))

In general I like to define field getter/setter functions like that for the fields of types I define since they come in handy and can simplify a lot of code, and make for better generic code.

mtfishman · 2024-12-19T14:32:15Z

src/fusiontensor/base_interface.jl

+
+# tensor addition is a block data_matrix add.
+function Base.:+(left::FusionTensor, right::FusionTensor)
+  @assert matching_axes(axes(left), axes(right))


Probably would be good to wrap this into a function, and throw a more specific error:

function check_matching_axes(ax1, ax2) matching_axes(ax1, ax2) || throw(DimensionMismatch("axes don't match.")) return nothing end

If you want to go the Base way, similar to checkbounds, you could define

checkaxes(Bool, ax1, ax2) -> Bool checkaxes(ax1, ax2) = checkaxes(Bool, ax1, ax2) || throw(DimensionMismatch(lazy"$ax1 does not match $ax2 ))

I forgot that is the pattern used in Base, that seems reasonable. I always assume checkf(...) means the function will throw an error or not but the pattern of passing in the desired output type is a good one to keep in mind.

mtfishman · 2024-12-19T14:32:52Z

src/fusiontensor/array_cast.jl

+  return ft
+end
+
+function cast_to_array(ft::FusionTensor)


I think we can just call this to_array(...).

If there's no ambiguity, could you also just write Base.Array(::FusionTensor)?

I think the issue with that is that it actually outputs a BlockSparseArray. We could call it BlockSparseArrays.BlockSparseArray. But I also like the idea of having a generic converter from FusionTensor to some kind of non-fused array, which is agnostic from the user perspective about what the non-fused format will be. So maybe:

function BlockSparseArrays.BlockSparseArray(ft::FusionTensor) # This function implementation we are discussing end # Maybe come up with a better name unfuse(ft::FusionTensor) = BlockSparseArray(ft)

Also, it seems like we should support:

function BlockSparseArrays.BlockSparseArray(ft::FusionTensor{<:Any,N}, axes::Tuple{Vararg{AbstractUnitRange,N}}) where {N} # ... end

so that you can attach graded axes to the output, say to specify a set of abelian symmetries. Right now I see it just makes it blocked but doesn't have any sectors on it which makes sense since as we've discussed it isn't necessarily well defined what abelian group the original symmetry group will reduce to (though if it is abelian to start with we could just use the same abelian group by default).

Is there a way to choose a reasonable guess for an abelian group given a non-abelian group? We could even just define a default abelian group for certain non-abelian groups, and if the defaults aren't defined the function doesn't specify sectors on the axes or errors.

mtfishman · 2024-12-19T14:34:01Z

src/fusiontensor/base_interface.jl

+# eachindex is automatically defined for AbstractArray. We do not want it.
+Base.eachindex(::FusionTensor) = error("eachindex not defined for FusionTensor")


It seems ok to me to define it, but then restrict scalar indexing by default, i.e. there is no harm in defining this but still leaving Base.getindex(::FusionTensor{<:Any,N}, I::Vararg{Int,N}) where {N} and Base.setindex!(::FusionTensor{<:Any,N}, value, I::Vararg{Int,N}) where {N} undefined for now.

Though I do think it would be nice to define scalar getindex/setindex! eventually for the sake of testing, even if they are inefficient, and we can disable them by default and provide a macro @allowscalar to enable them in a scope, like how the GPU array libraries do it (see https://cuda.juliagpu.org/stable/usage/workflow/#UsageWorkflowScalar).

Would the correct definition be:

function Base.eachindex(::FusionTensor) ax = blockedrange.(block_dimensions.((codomain_axes(ft)..., domain_axes(ft)...))) return CartesianIndices(ax) end

? But also shouldn't we define:

function Base.axes(::FusionTensor) return blockedrange.(block_dimensions.((codomain_axes(ft)..., domain_axes(ft)...))) end

? I.e. I would think these should be defined so that axes(::FusionTensor), size(::FusionTensor), eachindex(::FusionTensor) match what you would get for those functions if you cast the FusionTensor to an array, i.e. axes(ft) == axes(to_array(ft)).

Something that is confusing about the current set of definitions is that size(ft) != length.(axes(ft)). So I think we may just need different kinds of size/axes definitions, i.e.:

fusiontensor_axes(ft::FusionTensor) = (codomain_axes(ft)..., domain_axes(ft)...) fusiontensor_size(ft::FusionTensor) = length.(fusiontensor_axes(ft)) Base.axes(ft::FusionTensor) = blockedrange.(block_dimensions.(fusiontensor_axes(ft))) Base.size(ft::FusionTensor) = length.(axes(ft))

Probably fusiontensor_axes/fusiontensor_size isn't the best naming convention but you get the idea.

It seems like we need some terminology for this, like reduced_axes/reduced_size, fused_axes/fused_size (though that might not be good since it could be mistaken for the axes of the data_matrix), etc.

An issue with the definition of axes I suggest above is that you don't get the property that similar(a, axes(a)) gets you back an object of type a, i.e. it can't recreate that the object should be a FusionTensor from the axes.

Maybe there should be a special type:

struct FusionTensorAxes{N,Ncodomain,Ndomain,CodomainAxes<:Tuple{Vararg{AbstractUnitRange,Ncodomain}},DomainAxes<:Tuple{Vararg{AbstractUnitRange,Ndomain}}} codomain_axes::CodomainAxes domain_axes::DomainAxes end

that is the output of axes(::FusionTensor). Note it can't just be a Tuple{Tuple,Tuple} since it is more like a blocked tuple since we want a flat structure of the axes, and in fact we should probably make a BlockedTuple type to be the output of axes(::FusionTensor), related to the TensorAlgebra.BlockedPermutation type. That doesn't quite get us the property size(ft::FusionTensor) == length.(axes(ft)) when there are non-abelian symmetries, somehow we want graded axes with non-abelian symmetries to have two kinds of lengths. I know that is the purpose of block_dimensions to some extent but I'm still hung up on what the default behavior should be, it seems like we want a way to wrap a GradedUnitRange in a "view" that makes it act like its abelian counterpart in terms of the lengths, block length, etc. since that is what is expected by the Julia Base API.

Also, it slightly breaks the abstraction we have with the data matrix of the fusion tensor, since the size of the axes won't be the same as the sizes of the actual data we are storing. In principle that's ok, but I'm not clear on how we will handle that in practice so we'll have to think that through.

Maybe we can define a special AbstractBlockSparseMatrix subtype that conceptually has block lengths that correspond to the vector space dimensions, but in practice only stores block matrices corresponding to the number of independent parameters, i.e.:

struct GroupSymmetricBlockSparseMatrix{T,M<:AbstractBlockSparseMatrix{T},Axes<:Tuple{AbstractGradedUnitRange,AbstractGradedUnitRange}} <: AbstractBlockSparseMatrix{T} data_matrix::M axes::Axes end

(Probably there is a better name but this is just for the sake of demonstration.) Alternatively, we could just not attach the sector information to the data_matrix of the FusionTensor and leave the sectors implicit, which I see is basically the strategy TensorKit.TensorMap takes, but I like the idea of keeping it explicit if possible, and maybe that GroupSymmetricBlockSparseMatrix abstraction can be helpful in practice.

I think there is no real way to get around the size of the data matrix differing from the actual size, and this is really what you would want.
I do think that maybe the abstraction is that the datamatrix really is Dict{SymmetrySector,Matrix} (or maybe a vector of pairs, depending on if you want to allow duplications or not). The actual matrix really is the direct sum of kron(identity(quantumdim(sector)), matrix), so the combination of just the data blocks by themselves is a bit of an ill-defined object. For example, the eigenvalues and singular values have to be duplicated, and correspondingly the norm is not simply the norm of the separate blocks. In that sense, I would kind of say that thinking of the "summed size" of the stored data matrices might hurt more then it helps, since you never really want to handle the object that is just the direct sum of these blocks instead of with the kronecker products in there.

I see your point, but also the abstraction of putting the minimal data in a block sparse matrix does work for a variety of operations, for example tensor/matrix multiplication, SVD, etc. can be implemented as those operations mapped over the blocks of the block sparse/diagonal data matrix, though as you say for some of them like SVD when you are truncating you have to take into account that the singular values are repeated. So my motivations/goals would be:

Use the BlockSparseMatrix data structure for the data matrix so we can leverage block sparse SVD, matrix multiplication, etc. that we are already defining in BlockSparseArrays.jl anyway so we can minimize code duplication, take advantage of optimizations we do in BlockSparseArrays.jl, take advantage of the interfaces we are developing for handling block sparse arrays and block arrays such as iterating and slicing blocks, etc., and

Make the abstraction as explicit as possible through the wrapper type GroupSymmetricBlockSparseMatrix, since in the end this is just another sparse data structure with rules about what elements you store and how those map to the elements of the actual matrix in the full space.

So with that kind of GroupSymmetricBlockSparseMatrix design, the way I would see it is that we could achieve what you are describing in a more explicit way, by defining things like:

function Base.:*(a::GroupSymmetricBlockSparseMatrix, b::GroupSymmetricBlockSparseMatrix) return GroupSymmetricBlockSparseMatrix(parent(a) * parent(b), axes(a, 1), axes(b, 2)) end function LinearAlgebra.norm(a::GroupSymmetricBlockSparseMatrix) # Norm of the blocks of `a` scaled by the quantum dimension, say using a mapreduce over stored blocks end function LinearAlgebra.svd(a::GroupSymmetricBlockSparseMatrix) U, S, V = svd(parent(a)) # Attach sector information from `axes(a)` to `U`, `S`, and `V`, wrapping them in `GroupSymmetricBlockSparseMatrix`/`GroupSymmetricBlockSparseVector` end function truncated_svd(a::GroupSymmetricBlockSparseMatrix) # Perform `truncated_svd` taking into account the additional degeneracy of the singular values end

etc. Then in that way, those corresponding operations on the FusionTensor can just map (after permutation) to the same operation on the GroupSymmetricBlockSparseMatrix stored inside the FusionTensor.

EDIT: To be clear, basically GroupSymmetricBlockSparseMatrix (or whatever you want to call it, maybe InvariantBlockMatrix is a better name) would store the minimal data as a BlockDiagonalMatrix, but implicitly act like the block sparse matrix where the blocks are repeated as you wrote above, i.e. for each sector sector and block matrix matrix, it implicitly has blocks kron(identity(quantumdim(sector)), matrix).

I like that wrapper idea. So with that change, and with further features completed for the BlockSparseMatrix type would a lot of code in this package become duplicate and get removed? (Just to make sure I'm following.)

I'm glad to see this get caught about the difference between length.(axes(ft)) and size(ft) since this new type represents one of the biggest changes of any storage type we've had so far, relative to the other ones. I just mean in the sense that it's still not 100% clear to me from an ITensor user level what dim will return on tensor indices and how maxdim will be interpreted in the SVD etc. I was looking forward to discussing what the best choices are there.

I like that wrapper idea. So with that change, and with further features completed for the BlockSparseMatrix type would a lot of code in this package become duplicate and get removed? (Just to make sure I'm following.)

That would be the goal, see a prototype I shared in the "symmetric_tensors" Slack channel. The current design already tries to ensure as much functionality as possible gets implemented in the BlockSparseArrays.jl library (or more specifically whatever can be mapped to what one would reasonably consider to be generic block sparse array logic and functionality), this design proposal is meant to ensure even more of it can be, while taking the size/axes seriously.

I'm glad to see this get caught about the difference between length.(axes(ft)) and size(ft) since this new type represents one of the biggest changes of any storage type we've had so far, relative to the other ones. I just mean in the sense that it's still not 100% clear to me from an ITensor user level what dim will return on tensor indices and how maxdim will be interpreted in the SVD etc. I was looking forward to discussing what the best choices are there.

I think basically from the user perspective they should be defined such that they map to the equivalent definition for the abelian counterpart of the same tensor, from a generic code perspective a user should be able to switch between U1 and SU2 symmetries for the same system and keep the truncation parameters the same, be able to interpret the tensor sizes/bond dimensions in the same way, etc. (obviously they won't get the exact same result since the variational space changes, but things like tolerances in the SVD, sizes of tensors, etc. should have the same interpretation). So, that is another reason that we should define things such that length.(axes(ft)) == size(ft) == length.(axes(BlockSparseArray(ft))) == size(BlockSparseArray(ft)).

mtfishman · 2024-12-19T14:34:37Z

src/fusiontensor/base_interface.jl

+# 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, f1f2::Tuple{<:SectorFusionTree,<:SectorFusionTree})


Is there a reason to keep both the Tuple and splatted version? I would prefer just the splatted version.

mtfishman · 2024-12-19T14:35:22Z

src/fusiontensor/base_interface.jl

+  return view(ft, f1, f2) .= a
+end
+
+Base.ndims(::FusionTensor{T,N}) where {T,N} = N


Isn't this already defined since FusionTensor{T,N} <: AbstractArray{T,N}?

mtfishman · 2024-12-19T14:36:29Z

src/fusiontensor/base_interface.jl

+function Base.similar(ft::FusionTensor)
+  mat = similar(data_matrix(ft))
+  initialize_allowed_sectors!(mat)
+  return FusionTensor(mat, codomain_axes(ft), domain_axes(ft), trees_block_mapping(ft))
+end


I think you can just define this as:

Base.similar(ft::FusionTensor) = similar(ft, eltype(ft))

If it's an AbstractArray, the default will always map it to similar(ft, eltype(ft), axes(ft)), so you actually would only have to define that.

Good point, though maybe in this case we need both similar(::FusionTensor, ::Type) and similar(::FusionTensor, ::Type, axes), since axes(ft) doesn't output bipartitioned axes (that also connects to my question about how we should define axes(::FusionTensor) in the first place...).

So, related to the discussion elsewhere about the axes definition that was initiated by the eachindex definition, we should aim to output bipartitioned axes from axes(::FusionTensor) (maybe as a BlockedTuple type, related to the TensorAlgebra.BlockedPermutation type or something like that), and then only define similar(::FusionTensor, ::Type, axes) as suggested by @lkdvos, and then rely on the Base definitions of similar(::FusionTensor) and similar(::FusionTensor, ::Type) that handle the forwarding.

mtfishman · 2024-12-19T14:40:44Z

I have issues with building test both on local and on GitHub actions. Locally, I get an error
  LoadError: ArgumentError: Package LinearAlgebra not found in current path.
although the package is included in the Project.toml. I guess there was an issue with SafeTestsets configuration when I merged the new skeleton. @lkdvos do you have a clue?

Maybe it needs to be added to test/Project.toml as well. Dependencies of the package aren't automatically seen as dependencies when testing (obviously they are implicitly used when installing the package, but that doesn't count when you do using LinearAlgebra in the tests since test/Project.toml is the environment that is used in the tests).

mtfishman · 2024-12-19T14:43:25Z

src/fusiontensor/fusiontensor.jl

+
+# empty matrix
+function FusionTensor(
+  data_type::Type,


I would call this elt, i.e. it is the element type (I would interpret data_type as typeof(data_matrix(ft))).

mtfishman · 2024-12-19T14:50:09Z

Looks like a great start, thanks @ogauthe. That's a good list of issues to address in future PRs. Can you raise separate issues for those so we can track them?

In this PR I'll focus on things you didn't list in the first post, and also mostly focus on "shallow" issues like names of functions and things like that. Also as I go along with the review I might raise broader issues to discuss as they come to mind, but we can split some of those off into issues to track.

mtfishman · 2024-12-19T14:58:34Z

src/fusion_trees/fusiontree.jl

+  return map(
+    k -> SectorFusionTree((arguments_leave[k],), arrows(f), arguments_root[k], (), ()),
+    keys(arguments_root),
+  )


Suggested change

return map(

k -> SectorFusionTree((arguments_leave[k],), arrows(f), arguments_root[k], (), ()),

keys(arguments_root),

)

return map(keys(arguments_root)) do k

return SectorFusionTree((arguments_leave[k],), arrows(f), arguments_root[k], (), ()),

end

You can also remove quite a few dict lookups by mapping over the pairs and not indexing arguments_root

mtfishman · 2024-12-19T15:01:36Z

src/fusion_trees/fusiontree.jl

+  return flatten_tuples(nested)
+end
+
+function grow_tree(


Could you write a comment in the code about what this function is doing?

mtfishman · 2024-12-19T15:02:00Z

src/fusion_trees/fusiontree.jl

+  end
+end
+
+function build_trees(old_trees::Vector, sectors_to_fuse::Tuple, arrows_to_fuse::Tuple)


Same thing here, please write a comment in the code about what this function is doing.

I guess this one is more obvious since I see that it is expanding thick fusion trees to thin/sector fusion trees, but maybe we can also have a more explicit name like expand_to_sector_fusion_trees.

mtfishman · 2024-12-19T15:04:50Z

src/fusion_trees/fusiontree.jl

+end
+
+# --------------- convert to Array  ---------------
+to_tensor(::SectorFusionTree{<:Any,0}) = ones(1)


I would call it to_array just to stick to Julia notation, I assume it is always outputting some kind of AbstractArray.

Array(::SectorFusionTree)? AbstractArray(::SectorFusionTree)?

Both of those make sense, we could define both, i.e.:

Base.Array(t::SectorFusionTree) = ... Base.AbstractArray(t::SectorFusionTree) = Array(t)

mtfishman · 2024-12-19T15:08:39Z

src/fusiontensor/fusiontensor.jl

+end
+
+# initialize with already computed data_matrix
+function fusiontensor(


As discussed, I think we need to decide on a naming convention for these constructors. For now my suggestion would be to call any constructor that takes in a data matrix FusionTensor (since those are more literally constructors, i.e. they are doing minimal manipulation of the input data), and any function that is converting from a dense array to_fusiontensor (since those are more like converters rather than constructors).

So with that convention, this one should be called FusionTensor.

mtfishman · 2024-12-19T15:12:10Z

src/fusiontensor/array_cast.jl

+
+#### cast from symmetric to array
+function BlockSparseArrays.BlockSparseArray(ft::FusionTensor)
+  return cast_to_array(ft)


Suggested change

return cast_to_array(ft)

return to_array(ft)

mtfishman · 2024-12-19T15:19:37Z

Also @lkdvos @emstoudenmire please look this over and leave comments as well.

mtfishman · 2024-12-19T15:21:39Z

src/fusiontensor/array_cast.jl

+# =================================  High level interface  =================================
+
+#### cast from array to symmetric
+function FusionTensor(


As discussed above, for now let's call this to_fusiontensor to distinguish it from the constructors that take a data matrix as the first argument.

Project.toml

codecov · 2024-12-19T18:17:20Z

Codecov Report

Attention: Patch coverage is 0% with 565 lines in your changes missing coverage. Please review.

Project coverage is 0.00%. Comparing base (b3c4f5a) to head (423b275).

Files with missing lines	Patch %	Lines
src/fusion_trees/fusiontree.jl	0.00%	113 Missing ⚠️
src/fusiontensor/fusiontensor.jl	0.00%	90 Missing ⚠️
src/fusiontensor/array_cast.jl	0.00%	76 Missing ⚠️
src/fusion_trees/clebsch_gordan_tensors.jl	0.00%	75 Missing ⚠️
src/fusiontensor/base_interface.jl	0.00%	74 Missing ⚠️
src/fusiontensor/fusedaxes.jl	0.00%	41 Missing ⚠️
src/fusiontensor/linear_algebra_interface.jl	0.00%	37 Missing ⚠️
src/permutedims/permutedims.jl	0.00%	29 Missing ⚠️
src/permutedims/unitaries.jl	0.00%	21 Missing ⚠️
src/fusiontensor/tensor_algebra_interface.jl	0.00%	9 Missing ⚠️

Additional details and impacted files

@@          Coverage Diff          @@
##           main      #5    +/-   ##
=====================================
  Coverage      0   0.00%            
=====================================
  Files         0      10    +10     
  Lines         0     565   +565     
=====================================
- Misses        0     565   +565

Flag	Coverage Δ
docs	`0.00% <0.00%> (?)`

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

ogauthe · 2024-12-19T19:00:55Z

So I added test dependencies and I am now able to pass tests locally with (FusionTensors) pkg> test FusionTensors. However they still fail on GitHub actions. Am I missing something? I cannot find any clear error message. I wonder if I need to specify somewhere to get packages from ITensorRegistry?

mtfishman · 2024-12-19T19:07:31Z

So I added test dependencies and I am now able to pass tests locally with (FusionTensors) pkg> test FusionTensors. However they still fail on GitHub actions. Am I missing something? I cannot find any clear error message. I wonder if I need to specify somewhere to get packages from ITensorRegistry?

I fixed it, the URL for BlockSparseArrays.jl was wrong in the ITensorRegistry. Now the only failure is an Aqua failure in Julia 1.10, we could just disable that for now.

mtfishman · 2024-12-19T19:53:14Z

src/fusiontensor/array_cast.jl

+get_tol(a::AbstractArray) = get_tol(real(eltype(a)))
+get_tol(T::Type{<:Integer}) = get_tol(Float64)
+get_tol(T::Type{<:Real}) = 10 * eps(T)


A naming convention I have been using is defining functions that provide default values for a keyword argument k as default_k(...), so here it would be default_tol(...).

Also, here's a suggestion for making it a bit more compact and general:

default_tol(a::AbstractArray) = default_tol(eltype(a)) default_tol(arrayt::Type{<:AbstractArray}) = default_tol(eltype(arrayt)) default_tol(elt::Type{<:Number}) = 10 * eps(real(float(elt)))

It might be nice to specify atol or rtol too, as this can become relevant for SU(N), or even larger irreps of SU(2), where the CGCs can become tiny

lkdvos · 2024-12-19T20:12:45Z

src/fusion_trees/clebsch_gordan_tensors.jl

+  cgt = zeros((d1, d2, d3))
+  s3 ∉ blocklabels(fusion_product(s1, s2)) && return cgt
+
+  # adapted from TensorKit


Given that this is adapted from TensorKit anyways, how would you feel about:

function clebsch_gordan_tensor(s1, s2, s3) s1_tk, s2_tk, s3_tk = convert.(TensorKitSectors.Sector, (s1, s2, s3)) cgt = TensorKitSectors.fusiontensor(s1_tk, s2_tk, s3_tk) return reshape(cgt, size(cgt, 1), size(cgt, 2), size(cgt, 3)) # discard trailing singleton dimension end

I don't think there is a lot of value in having the implementation copied here, it's definitely more maintenance and error-prone, so it would be cool if we could immediately share this.
(This is a registered light-weight package that only contains the symmetry data, so the dependency is quite light, and this automatically inherits many of the symmetries)
https://github.com/QuantumKitHub/TensorKitSectors.jl

That sounds reasonable to me.

lkdvos · 2024-12-19T20:14:18Z

src/fusion_trees/clebsch_gordan_tensors.jl

+  return sqrt(quantum_dimension(s)) * cgt[:, :, 1]
+end
+
+function clebsch_gordan_tensor(


If you use symbol_1j, it might be more internally consistent to then also use the symbol_3j and symbol_6j notations instead of clebsch_gordan_tensor? In general, clebsch gordan seems to be very often related specifically to SU(N) anyways.

lkdvos · 2024-12-19T20:19:00Z

src/fusion_trees/fusiontree.jl

+# The interface uses AbstractGradedUnitRanges as input for interface simplicity
+# however only blocklabels are used and blocklengths are never read.
+
+struct SectorFusionTree{S,N,M}


I think it would be nice to be consistent with the fieldnames: either add _sector(s) everywhere or nowhere? Given that sector is already in the typename, it might not give that much additional information to add it in the fields as well, so it could maybe be leaves, arrows, root, branches and vertices/multiplicities/... I usually get a bit confused by outer_multiplicity_indices because these don't actually live on the "outer" charges, but only inside of the tree. Ultimately I don't feel too strongly about this though, just sharing my opinion.

Those changes all sound reasonable to me. I agree the name outer_multiplicity_indices is a bit confusing.

lkdvos · 2024-12-19T20:19:48Z

src/fusion_trees/fusiontree.jl

+  # TBD could have branch_sectors with length N-2
+  # currently first(branch_sectors) == first(leaves)


If deciding to keep this, definitely should make sure that this is properly documented :)

lkdvos · 2024-12-19T20:21:39Z

src/fusion_trees/fusiontree.jl

+Base.length(::SectorFusionTree{<:Any,N}) where {N} = N
+
+# GradedUnitRanges interface
+GradedUnitRanges.sector_type(::SectorFusionTree{S}) where {S} = S


These kinds of functions are typically implemented in the type-domain instead:

Suggested change

GradedUnitRanges.sector_type(::SectorFusionTree{S}) where {S} = S

GradedUnitRanges.sector_type(ft::SectorFusionTree) = sectortype(typeof(ft))

GradedUnitRanges.sector_type(::Type{<:SectorFusionTree{S}}) where {S} = S

lkdvos · 2024-12-19T20:50:06Z

src/fusiontensor/fusiontensor.jl

+  data_type::Type{<:Number}, codomain_fused_axes::FusedAxes, domain_fused_axes::FusedAxes
+)
+  # fusion trees have Float64 eltype: need compatible type
+  promoted = promote_type(data_type, Float64)


probably don't want to hard-code that here, since this is not always the case (could also be Complex, and for abelian symmetries it is even convenient to use Bool). eltype(fusiontensortype) or something similar might work?

lkdvos · 2024-12-19T20:52:47Z

src/fusiontensor/linear_algebra_interface.jl

+    b -> quantum_dimension(row_sectors[Int(first(Tuple(b)))]) * norm(m[b])^2,
+    +,
+    eachblockstoredindex(m);
+    init=0.0,


Suggested change

init=0.0,

init=zero(real(scalartype(ft))),

lkdvos · 2024-12-19T20:56:16Z

src/fusiontensor/linear_algebra_interface.jl

+  n2 = mapreduce(
+    b -> quantum_dimension(row_sectors[Int(first(Tuple(b)))]) * norm(m[b])^2,
+    +,
+    eachblockstoredindex(m);
+    init=0.0,


I tend to find sum a little easier to read, but that's just a personal preference

Suggested change

n2 = mapreduce(

b -> quantum_dimension(row_sectors[Int(first(Tuple(b)))]) * norm(m[b])^2,

+,

eachblockstoredindex(m);

init=0.0,

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

lkdvos · 2024-12-19T20:56:42Z

src/fusiontensor/linear_algebra_interface.jl

+    b -> quantum_dimension(row_sectors[Int(first(Tuple(b)))]) * tr(m[b]),
+    +,
+    eachblockstoredindex(m);
+    init=eltype(ft)(0),


Suggested change

init=eltype(ft)(0),

init=zero(eltype(ft)),

lkdvos · 2024-12-19T20:58:05Z

src/fusiontensor/linear_algebra_interface.jl

+function LinearAlgebra.qr(ft::FusionTensor)
+  qmat, rmat = block_qr(data_matrix(ft))
+  qtens = FusionTensor(qmat, codomain_axes(ft), (axes(qmat)[1],))
+  rtens = FusionTensor(rmat, (axes(rmat)[0],), domain_axes(ft))


Why is this indexed with 0? also, it can probably be axes(rmat, 0) ?

@lkdvos that looks right but this looks like it is just prototype/stand-in code anyway since block_qr isn't defined anywhere.

mtfishman · 2024-12-19T22:00:59Z

src/fusiontensor/base_interface.jl

+  return nothing
+end
+
+Base.size(ft::FusionTensor) = quantum_dimension.(axes(ft))


Related to the discussion elsewhere, probably we should redefine length(::AbstractGradedUnitRange) to output the number of basis vectors in the vector space as opposed to the number of independent parameters, so then this can be defined as: Base.size(ft::FusionTensor) = length.(axes(ft)) (though I think we won't need to define this overload since Base defines size automatically for AbstractArrays as long as you define axes).

ogauthe added 28 commits December 4, 2024 20:10

split project

3dd5ff2

write intersect_sectors

11a5d88

define tree_to_block_mapping

afed27a

define getindex

5ed5b8e

array cast with FusionTree

d66540c

clean Base

a9dc29f

reuse tree_block_mapping in Base operations

d7edb71

working permutedims

4a8f9ba

fix contract

e75d783

put imports on top of files

705b7aa

move check_sanity to tests

deaf6b1

rewrite test FusionTree

8ab786e

cleaning

563bae1

cache unitary coeff

b35757b

fix tests

0299e4b

simplify cast_to_array

4647030

update to BlockSparseArrays 0.2.0

e53bffd

use map_blocklabels

2794d60

fix signatures for fusiontensor_permutedims

a23939c

define unsafe_cast_from_array

f9b0daa

use @strided

d816193

rename FusionTree SectorFusionTree

899bed8

use @view!

3ef076d

use fusiontensor when data_matrix is provided

cdd3a2f

Merge branch 'main' into dev

9ed1f1a

move tests

4f746b8

check matrix axes in fusiontensor

6528b63

remove unused variable

cc9c7dc

mtfishman reviewed Dec 19, 2024

View reviewed changes

src/fusiontensor/array_cast.jl Outdated Show resolved Hide resolved

mtfishman reviewed Dec 19, 2024

View reviewed changes

src/fusiontensor/array_cast.jl Outdated Show resolved Hide resolved

mtfishman reviewed Dec 19, 2024

View reviewed changes

fix tests

067403b

mtfishman reviewed Dec 19, 2024

View reviewed changes

Project.toml Outdated Show resolved Hide resolved

Loosen LinearAlgebra compat

ee06e3f

fix names in array_cast

837e151

mtfishman reviewed Dec 19, 2024

View reviewed changes

lkdvos reviewed Dec 19, 2024

View reviewed changes

mtfishman reviewed Dec 19, 2024

View reviewed changes

Bump TypeParameterAccessors to v0.2

423b275

		# eachindex is automatically defined for AbstractArray. We do not want it.
		Base.eachindex(::FusionTensor) = error("eachindex not defined for FusionTensor")

		# TBD could have branch_sectors with length N-2
		# currently first(branch_sectors) == first(leaves)

	GradedUnitRanges.sector_type(::SectorFusionTree{S}) where {S} = S
	GradedUnitRanges.sector_type(ft::SectorFusionTree) = sectortype(typeof(ft))
	GradedUnitRanges.sector_type(::Type{<:SectorFusionTree{S}}) where {S} = S

First draft of FusionTensor #5

Are you sure you want to change the base?

First draft of FusionTensor #5

Conversation

ogauthe commented Dec 19, 2024

ogauthe commented Dec 19, 2024

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

mtfishman Dec 19, 2024 • edited Loading

Choose a reason for hiding this comment

mtfishman Dec 19, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

mtfishman Dec 19, 2024 • edited Loading

Choose a reason for hiding this comment

mtfishman Dec 19, 2024 • edited Loading

Choose a reason for hiding this comment

mtfishman Dec 19, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

mtfishman Dec 19, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

mtfishman Dec 20, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

mtfishman Dec 19, 2024 • edited Loading

Choose a reason for hiding this comment

mtfishman commented Dec 19, 2024

Choose a reason for hiding this comment

mtfishman commented Dec 19, 2024

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

mtfishman commented Dec 19, 2024

Choose a reason for hiding this comment

codecov bot commented Dec 19, 2024 • edited Loading

Codecov Report

ogauthe commented Dec 19, 2024

mtfishman commented Dec 19, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

mtfishman Dec 19, 2024 •

edited

Loading

mtfishman Dec 19, 2024 •

edited

Loading

mtfishman Dec 19, 2024 •

edited

Loading

mtfishman Dec 19, 2024 •

edited

Loading

mtfishman Dec 19, 2024 •

edited

Loading

mtfishman Dec 19, 2024 •

edited

Loading

mtfishman Dec 20, 2024 •

edited

Loading

mtfishman Dec 19, 2024 •

edited

Loading

codecov bot commented Dec 19, 2024 •

edited

Loading

mtfishman commented Dec 19, 2024 •

edited

Loading