Skip to content

Construct BlockSparseArray when slicing with graded unit ranges #36

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 15 commits into from
Feb 4, 2025

Conversation

mtfishman
Copy link
Member

@mtfishman mtfishman commented Feb 3, 2025

This requires ITensor/DerivableInterfaces.jl#14 getting fixed, in progress here: ITensor/DerivableInterfaces.jl#15.

With that fixed, this enables:

using BlockSparseArrays: BlockSparseArrays
using GradedUnitRanges: gradedrange
using SymmetrySectors: U1, dual

r = gradedrange([U1(0) => 2, U1(1) => 3])
a = cat(randn(2, 2), randn(3, 3); dims=(1, 2))
b = a[r, dual(r)]

which gives:

julia> b
2×2-blocked 5×5 BlockSparseArrays.BlockSparseMatrix{}:
  0.5699044959108884  -0.5270234956792514  │    .                    .                    .                
 -0.1657624514089808   0.6835463111834327  │    .                    .                    .                
 ──────────────────────────────────────────┼───────────────────────────────────────────────────────────────
   .                    .                  │  -0.36206253094583096  1.6086749840883592  -1.0391612058066324
   .                    .                  │   0.6750593012545241   1.3531901922413943  -1.1126126506657081
   .                    .                  │   1.2064441481034023   0.5393283008204176  -1.0247393221148995

julia> axes(b)
(GradedOneTo[U(1)[0] => 1:2, U(1)[1] => 3:5], GradedUnitRangeDual[U(1)[0] => 1:2, U(1)[-1] => 3:5])

so slicing a dense array by graded unit ranges creates a block sparse array, and additionally drops zero blocks. Note that it also preserves dual information, which is enabled by ITensor/GradedUnitRanges.jl#10.

This fixes ITensor/GradedUnitRanges.jl#9, and as discussed this will be helpful for constructing abelian symmetric operators and states in https://github.com/ITensor/QuantumOperatorDefinitions.jl.

Copy link

codecov bot commented Feb 4, 2025

Codecov Report

Attention: Patch coverage is 80.23256% with 17 lines in your changes missing coverage. Please review.

Project coverage is 73.83%. Comparing base (37fecf7) to head (82abcf7).
Report is 1 commits behind head on main.

Files with missing lines Patch % Lines
src/abstractblocksparsearray/map.jl 0.00% 6 Missing ⚠️
...tRangesExt/BlockSparseArraysGradedUnitRangesExt.jl 77.77% 4 Missing ⚠️
...ksparsearrayinterface/blocksparsearrayinterface.jl 20.00% 4 Missing ⚠️
src/blocksparsearrayinterface/map.jl 92.30% 3 Missing ⚠️
Additional details and impacted files
@@            Coverage Diff             @@
##             main      #36      +/-   ##
==========================================
- Coverage   74.02%   73.83%   -0.20%     
==========================================
  Files          29       30       +1     
  Lines        1001     1051      +50     
==========================================
+ Hits          741      776      +35     
- Misses        260      275      +15     
Flag Coverage Δ
docs 26.20% <41.86%> (-0.32%) ⬇️

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.

@mtfishman
Copy link
Member Author

mtfishman commented Feb 4, 2025

I believe the only issues remaining are due to a change to DerivableInterfaces.jl's cat function introduced in ITensor/DerivableInterfaces.jl#13. (EDIT: Fixed by ITensor/DerivableInterfaces.jl#17.)

@mtfishman
Copy link
Member Author

I've added tests. There is one final test failure caused by the fact that the changes to the broadcast/map functionality leads to this new behavior:

julia> using BlockSparseArrays: BlockSparseMatrix

julia> using BlockArrays: Block

julia> a = BlockSparseMatrix{Float64}([2, 2], [2, 2])
2×2-blocked 4×4 BlockSparseMatrix{Float64, Matrix{Float64}, SparseArraysBase.SparseMatrixDOK{Matrix{Float64}, BlockSparseArrays.GetUnstoredBlock{Tuple{BlockedOneTo{Int64, Vector{Int64}}, BlockedOneTo{Int64, Vector{Int64}}}}}, Tuple{BlockedOneTo{Int64, Vector{Int64}}, BlockedOneTo{Int64, Vector{Int64}}}}:
 .  .  │  .  .
 .  .  │  .  .
 ──────┼──────
 .  .  │  .  .
 .  .  │  .  .

julia> a[Block(1, 1)] = randn(2, 2)
2×2 Matrix{Float64}:
 -0.528387  -0.36891
  0.451625   1.32933

julia> a .= 0
2×2-blocked 4×4 BlockSparseMatrix{Float64, Matrix{Float64}, SparseArraysBase.SparseMatrixDOK{Matrix{Float64}, BlockSparseArrays.GetUnstoredBlock{Tuple{BlockedOneTo{Int64, Vector{Int64}}, BlockedOneTo{Int64, Vector{Int64}}}}}, Tuple{BlockedOneTo{Int64, Vector{Int64}}, BlockedOneTo{Int64, Vector{Int64}}}}:
 0.0  0.0  │  .  .
 0.0  0.0  │  .  .
 ──────────┼──────
  .    .   │  .  .
  .    .   │  .  .

while ideally it would drop the block, like fill!(a, 0) does:

julia> fill!(a, 0)
2×2-blocked 4×4 BlockSparseMatrix{Float64, Matrix{Float64}, SparseArraysBase.SparseMatrixDOK{Matrix{Float64}, BlockSparseArrays.GetUnstoredBlock{Tuple{BlockedOneTo{Int64, Vector{Int64}}, BlockedOneTo{Int64, Vector{Int64}}}}}, Tuple{BlockedOneTo{Int64, Vector{Int64}}, BlockedOneTo{Int64, Vector{Int64}}}}:
 .  .  │  .  .
 .  .  │  .  .
 ──────┼──────
 .  .  │  .  .
 .  .  │  .  .

I'll fix that and then merge.

@mtfishman mtfishman closed this Feb 4, 2025
@mtfishman mtfishman reopened this Feb 4, 2025
@mtfishman mtfishman merged commit 243e917 into main Feb 4, 2025
21 checks passed
@mtfishman mtfishman deleted the gradedunitrange_slicing branch February 4, 2025 20:49
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

[BUG] Issues slicing dense arrays by graded unit ranges
1 participant