Jutho / TensorKit.jl

A Julia package for large-scale tensor computations, with a hint of category theory
MIT License
218 stars 38 forks source link

Constructing a `TensorMap` with a different shape to input data fails when index space is symmetric #143

Closed jack-dunham closed 1 month ago

jack-dunham commented 1 month ago

Since v0.12.5, attempting to construct TensorMap with data of a different (but compatible) shape to that specified by the provided codomain and domain will result in an error when the space is non-trivial.

In [1]: using TensorKit; s = U1Space(-1=>1,1=>1); TensorMap(zeros(4,4),s*s,s*s)
ERROR: BoundsError: attempt to access 4×4 StridedViews.StridedView{Float64, 2, Matrix{Float64}, typeof(identity)} at index [2:2, 1:1, 2:2, 1:1]
Stacktrace:
 [1] throw_boundserror(A::StridedViews.StridedView{Float64, 2, Matrix{Float64}, typeof(identity)}, I::NTuple{4, UnitRange{Int64}})
   @ Base ./abstractarray.jl:734
 [2] checkbounds
   @ ./abstractarray.jl:699 [inlined]
 [3] _getindex
   @ ./multidimensional.jl:888 [inlined]
 [4] getindex
   @ ./abstractarray.jl:1288 [inlined]
 [5] project_symmetric!(t::TensorMap{GradedSpace{…}, 2, 2, U1Irrep, TensorKit.SortedVectorDict{…}, FusionTree{…}, FusionTree{…}}, data::Matrix{Float64})
   @ TensorKit ~/.julia/packages/TensorKit/Fz2UW/src/tensors/tensor.jl:387
 [6] TensorMap(data::Matrix{Float64}, codom::ProductSpace{GradedSpace{U1Irrep, TensorKit.SortedVectorDict{U1Irrep, Int64}}, 2}, dom::ProductSpace{GradedSpace{U1Irrep, TensorKit.SortedVectorDict{U1Irrep, Int64}}, 2}; tol::Float64)
   @ TensorKit ~/.julia/packages/TensorKit/Fz2UW/src/tensors/tensor.jl:361
 [7] TensorMap(data::Matrix{Float64}, codom::ProductSpace{GradedSpace{U1Irrep, TensorKit.SortedVectorDict{U1Irrep, Int64}}, 2}, dom::ProductSpace{GradedSpace{U1Irrep, TensorKit.SortedVectorDict{U1Irrep, Int64}}, 2})
   @ TensorKit ~/.julia/packages/TensorKit/Fz2UW/src/tensors/tensor.jl:346
 [8] top-level scope
   @ REPL[19]:1
 [9] top-level scope
   @ ~/.julia/juliaup/julia-1.10.0+0.aarch64.apple.darwin14/share/julia/stdlib/v1.10/REPL/src/REPL.jl:1428
Some type information was truncated. Use `show(err)` to see complete types.

This works fine for ComplexSpace; the data automatically gets reshaped to the shape of destination TensorMap:

In [2]: s = ComplexSpace(2); TensorMap(zeros(4,4),s*s,s*s)
TensorMap((ℂ^2 ⊗ ℂ^2) ← (ℂ^2 ⊗ ℂ^2)):
[:, :, 1, 1] =
 0.0  0.0
 0.0  0.0

[:, :, 2, 1] =
 0.0  0.0
 0.0  0.0

[:, :, 1, 2] =
 0.0  0.0
 0.0  0.0

[:, :, 2, 2] =
 0.0  0.0
 0.0  0.0

As this is an inconsistency between types of IndexSpace and a breaking change since v0.12.4, would it be convenient to reshape the data within the constructor before calling project_symmetric!?

Thanks for all the work on this great package, Jack

lkdvos commented 1 month ago

My apologies for the slow reply.

It indeed looks like I overlooked this use-case when I changed this. I will try and fix this in the master, and then attempt to backport it to 0.12.5 as a hotfix.

lkdvos commented 1 month ago

This should now be fixed as of version v0.12.6, tagged as a patch release such that you should be able to automatically update. I hope this resolves your issues, feel free to ping me if anything is still missing.