-
Notifications
You must be signed in to change notification settings - Fork 29
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
[WIP] BlockIndices
#356
base: master
Are you sure you want to change the base?
[WIP] BlockIndices
#356
Conversation
@jishnub let me know what you think of the design of what is implemented, and also what you think the scope of the PR should be, i.e. how much work should go into getting more slicing working vs. iterating on that in future PRs. |
Codecov ReportAttention: Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## master #356 +/- ##
===========================================
+ Coverage 0.00% 94.16% +94.16%
===========================================
Files 16 17 +1
Lines 1480 1559 +79
===========================================
+ Hits 0 1468 +1468
+ Misses 1480 91 -1389 ☔ View full report in Codecov by Sentry. |
In the latest commit, I made it so that more slicing operations return either julia> using BlockArrays
julia> using BlockArrays: BlockIndices
julia> A = BlockArray(randn(4, 4), [2, 2], [2, 2])
2×2-blocked 4×4 BlockMatrix{Float64}:
-1.38012 -1.38908 │ 0.653502 -0.176679
-0.61462 -1.30911 │ -0.797909 0.083623
─────────────────────┼──────────────────────
-0.889299 1.48537 │ -0.83928 0.470424
1.01574 -1.31194 │ -0.282103 0.863026
julia> B = BlockIndices(A)
2×2-blocked 4×4 BlockIndices{2, Tuple{BlockedUnitRange{Vector{Int64}}, BlockedUnitRange{Vector{Int64}}}}:
Block(1, 1)[1, 1] Block(1, 1)[1, 2] │ Block(1, 2)[1, 1] Block(1, 2)[1, 2]
Block(1, 1)[2, 1] Block(1, 1)[2, 2] │ Block(1, 2)[2, 1] Block(1, 2)[2, 2]
──────────────────────────────────────┼──────────────────────────────────────
Block(2, 1)[1, 1] Block(2, 1)[1, 2] │ Block(2, 2)[1, 1] Block(2, 2)[1, 2]
Block(2, 1)[2, 1] Block(2, 1)[2, 2] │ Block(2, 2)[2, 1] Block(2, 2)[2, 2]
julia> B[Block(1, 1)[1:2, 1:2]]
2×2 BlockIndexRange{2, Tuple{UnitRange{Int64}, UnitRange{Int64}}}:
Block(1, 1)[1, 1] Block(1, 1)[1, 2]
Block(1, 1)[2, 1] Block(1, 1)[2, 2]
julia> B[2:3, 2:3]
2×2-blocked 2×2 BlockIndices{2, Tuple{BlockedUnitRange{Vector{Int64}}, BlockedUnitRange{Vector{Int64}}}}:
Block(1, 1)[2, 2] │ Block(1, 2)[2, 1]
───────────────────┼───────────────────
Block(2, 1)[1, 2] │ Block(2, 2)[1, 1]
julia> B[Block(2):Block(2), Block(1):Block(2)]
1×2-blocked 2×4 BlockIndices{2, Tuple{BlockedUnitRange{Vector{Int64}}, BlockedUnitRange{Vector{Int64}}}}:
Block(2, 1)[1, 1] Block(2, 1)[1, 2] │ Block(2, 2)[1, 1] Block(2, 2)[1, 2]
Block(2, 1)[2, 1] Block(2, 1)[2, 2] │ Block(2, 2)[2, 1] Block(2, 2)[2, 2] I'm not sure if it is the best design for that code, but it isn't too complicated. Happy to refactor the code, change names and code style, etc. I also realized that the previous definition of |
I'm not fully sure about retaining block axes of the source in the indexing operation because of the reasons discussed earlier. E.g., this seems confusing: julia> B = BlockIndices(A)
2×2-blocked 4×4 BlockIndices{2, Tuple{BlockedUnitRange{Vector{Int64}}, BlockedUnitRange{Vector{Int64}}}}:
Block(1, 1)[1, 1] Block(1, 1)[1, 2] │ Block(1, 2)[1, 1] Block(1, 2)[1, 2]
Block(1, 1)[2, 1] Block(1, 1)[2, 2] │ Block(1, 2)[2, 1] Block(1, 2)[2, 2]
──────────────────────────────────────┼──────────────────────────────────────
Block(2, 1)[1, 1] Block(2, 1)[1, 2] │ Block(2, 2)[1, 1] Block(2, 2)[1, 2]
Block(2, 1)[2, 1] Block(2, 1)[2, 2] │ Block(2, 2)[2, 1] Block(2, 2)[2, 2]
julia> b = BlockArrays._BlockedUnitRange(1, [1,4])
2-blocked 4-element BlockedUnitRange{Vector{Int64}}:
1
─
2
3
4
julia> B[b, b]
2×2-blocked 4×4 BlockIndices{2, Tuple{BlockedUnitRange{Vector{Int64}}, BlockedUnitRange{Vector{Int64}}}}:
Block(1, 1)[1, 1] Block(1, 1)[1, 2] │ Block(1, 2)[1, 1] Block(1, 2)[1, 2]
Block(1, 1)[2, 1] Block(1, 1)[2, 2] │ Block(1, 2)[2, 1] Block(1, 2)[2, 2]
──────────────────────────────────────┼──────────────────────────────────────
Block(2, 1)[1, 1] Block(2, 1)[1, 2] │ Block(2, 2)[1, 1] Block(2, 2)[1, 2]
Block(2, 1)[2, 1] Block(2, 1)[2, 2] │ Block(2, 2)[2, 1] Block(2, 2)[2, 2] In this, the block axes of julia> A = BlockArray(randn(4, 4), [2, 2], [2, 2])
2×2-blocked 4×4 BlockMatrix{Float64}:
0.374444 0.888294 │ -1.25349 0.0794558
-0.230626 -0.21222 │ -1.50958 -1.98783
──────────────────────┼───────────────────────
0.846922 -0.493969 │ 0.557274 -1.0908
-1.35255 0.300943 │ -1.35591 -0.201546
julia> b = BlockArrays._BlockedUnitRange(1, [1,4])
2-blocked 4-element BlockedUnitRange{Vector{Int64}}:
1
─
2
3
4
julia> A[b, b]
2×2-blocked 4×4 BlockMatrix{Float64}:
0.374444 │ 0.888294 -1.25349 0.0794558
───────────┼──────────────────────────────────
-0.230626 │ -0.21222 -1.50958 -1.98783
0.846922 │ -0.493969 0.557274 -1.0908
-1.35255 │ 0.300943 -1.35591 -0.201546 I think it makes sense to be consistent with the convention here, where slices constructed using julia> B[b, b][Block(1), Block(1)]
2×2 BlockArrays.BlockIndexRange{2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}:
Block(1, 1)[1, 1] Block(1, 1)[1, 2]
Block(1, 1)[2, 1] Block(1, 1)[2, 2]
julia> B[b[Block(1)], b[Block(1)]]
1×1-blocked 1×1 BlockIndices{2, Tuple{BlockedUnitRange{Vector{Int64}}, BlockedUnitRange{Vector{Int64}}}}:
Block(1, 1)[1, 1] where both the two would be equivalent. |
Right, I realize the behavior around slicing in this PR is different from that for other I guess it is a different discussion to decide what slicing with |
In the latest, I made it so the axes of the julia> using BlockArrays
julia> using BlockArrays: BlockIndices
julia> A = BlockArray(randn(4, 4), [2, 2], [2, 2])
2×2-blocked 4×4 BlockMatrix{Float64}:
1.09381 0.613925 │ 0.275561 -0.0245506
0.586465 3.50312 │ 1.51753 -0.123905
─────────────────────┼───────────────────────
0.613102 0.288321 │ -0.121851 -1.01473
-0.387993 0.871633 │ 1.15325 1.0105
julia> B = BlockIndices(A)
2×2-blocked 4×4 BlockIndices{2, Tuple{BlockedUnitRange{Vector{Int64}}, BlockedUnitRange{Vector{Int64}}}, Tuple{BlockedUnitRange{Vector{Int64}}, BlockedUnitRange{Vector{Int64}}}}:
Block(1, 1)[1, 1] Block(1, 1)[1, 2] │ Block(1, 2)[1, 1] Block(1, 2)[1, 2]
Block(1, 1)[2, 1] Block(1, 1)[2, 2] │ Block(1, 2)[2, 1] Block(1, 2)[2, 2]
──────────────────────────────────────┼──────────────────────────────────────
Block(2, 1)[1, 1] Block(2, 1)[1, 2] │ Block(2, 2)[1, 1] Block(2, 2)[1, 2]
Block(2, 1)[2, 1] Block(2, 1)[2, 2] │ Block(2, 2)[2, 1] Block(2, 2)[2, 2]
julia> B[Block(1, 2)[1:2, 1:2]]
2×2 BlockArrays.BlockIndexRange{2, Tuple{UnitRange{Int64}, UnitRange{Int64}}}:
Block(1, 2)[1, 1] Block(1, 2)[1, 2]
Block(1, 2)[2, 1] Block(1, 2)[2, 2]
julia> B[Block.(2:2), Block.(1:2)]
1×2-blocked 2×4 BlockIndices{2, Tuple{BlockedUnitRange{Vector{Int64}}, BlockedUnitRange{Vector{Int64}}}, Tuple{BlockedUnitRange{Vector{Int64}}, BlockedUnitRange{Vector{Int64}}}}:
Block(2, 1)[1, 1] Block(2, 1)[1, 2] │ Block(2, 2)[1, 1] Block(2, 2)[1, 2]
Block(2, 1)[2, 1] Block(2, 1)[2, 2] │ Block(2, 2)[2, 1] Block(2, 2)[2, 2]
julia> B[1:4, 1:4]
1×1-blocked 4×4 BlockIndices{2, Tuple{BlockedUnitRange{Vector{Int64}}, BlockedUnitRange{Vector{Int64}}}, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}:
Block(1, 1)[1, 1] Block(1, 1)[1, 2] Block(1, 2)[1, 1] Block(1, 2)[1, 2]
Block(1, 1)[2, 1] Block(1, 1)[2, 2] Block(1, 2)[2, 1] Block(1, 2)[2, 2]
Block(2, 1)[1, 1] Block(2, 1)[1, 2] Block(2, 2)[1, 1] Block(2, 2)[1, 2]
Block(2, 1)[2, 1] Block(2, 1)[2, 2] Block(2, 2)[2, 1] Block(2, 2)[2, 2]
julia> B[2:3, 2:3]
1×1-blocked 2×2 BlockIndices{2, Tuple{BlockedUnitRange{Vector{Int64}}, BlockedUnitRange{Vector{Int64}}}, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}:
Block(1, 1)[2, 2] Block(1, 2)[2, 1]
Block(2, 1)[1, 2] Block(2, 2)[1, 1]
julia> B[blockedrange([1, 3]), blockedrange([1, 3])]
2×2-blocked 4×4 BlockIndices{2, Tuple{BlockedUnitRange{Vector{Int64}}, BlockedUnitRange{Vector{Int64}}}, Tuple{BlockedUnitRange{Vector{Int64}}, BlockedUnitRange{Vector{Int64}}}}:
Block(1, 1)[1, 1] │ Block(1, 1)[1, 2] Block(1, 2)[1, 1] Block(1, 2)[1, 2]
───────────────────┼─────────────────────────────────────────────────────────
Block(1, 1)[2, 1] │ Block(1, 1)[2, 2] Block(1, 2)[2, 1] Block(1, 2)[2, 2]
Block(2, 1)[1, 1] │ Block(2, 1)[1, 2] Block(2, 2)[1, 1] Block(2, 2)[1, 2]
Block(2, 1)[2, 1] │ Block(2, 1)[2, 2] Block(2, 2)[2, 1] Block(2, 2)[2, 2] This required a slight redesign of the @jishnub let me know what you think. |
One behavior I wasn't sure about was indexing with julia> B[:, :]
2×2-blocked 4×4 PseudoBlockMatrix{BlockIndex{2}}:
Block(1, 1)[1, 1] Block(1, 1)[1, 2] │ Block(1, 2)[1, 1] Block(1, 2)[1, 2]
Block(1, 1)[2, 1] Block(1, 1)[2, 2] │ Block(1, 2)[2, 1] Block(1, 2)[2, 2]
──────────────────────────────────────┼──────────────────────────────────────
Block(2, 1)[1, 1] Block(2, 1)[1, 2] │ Block(2, 2)[1, 1] Block(2, 2)[1, 2]
Block(2, 1)[2, 1] Block(2, 1)[2, 2] │ Block(2, 2)[2, 1] Block(2, 2)[2, 2] (though in principle it should be outputting a julia> A[:, :]
2×2-blocked 4×4 BlockMatrix{Float64}:
1.09381 0.613925 │ 0.275561 -0.0245506
0.586465 3.50312 │ 1.51753 -0.123905
─────────────────────┼───────────────────────
0.613102 0.288321 │ -0.121851 -1.01473
-0.387993 0.871633 │ 1.15325 1.0105 It seems a bit strange to me that |
Regarding the Colon, this is the standard behavior, which differs from indexing using a |
It looks like the following errors at present: julia> A = BlockArray(randn(4, 4), [2, 2], [2, 2]);
julia> B = BlockArrays.BlockIndices(A);
julia> A[B]
ERROR: ArgumentError: unable to check bounds for indices of type BlockIndex{2}
Stacktrace:
[1] checkindex(::Type{Bool}, inds::Base.OneTo{Int64}, i::BlockIndex{2})
@ Base ./abstractarray.jl:759
[2] checkindex
@ ./abstractarray.jl:776 [inlined]
[3] checkbounds
@ ./abstractarray.jl:687 [inlined]
[4] checkbounds
@ ./abstractarray.jl:702 [inlined]
[5] _getindex
@ ./multidimensional.jl:888 [inlined]
[6] getindex(A::BlockMatrix{Float64, Matrix{…}, Tuple{…}}, I::BlockIndices{2, Tuple{…}, Tuple{…}})
@ Base ./abstractarray.jl:1291
[7] top-level scope
@ REPL[14]:1
Some type information was truncated. Use `show(err)` to see complete types. I wonder if this may be made to work? |
I don't think this is handling offset axes correctly: julia> B
2×2-blocked 4×4 BlockIndices{2, Tuple{BlockedUnitRange{Vector{Int64}}, BlockedUnitRange{Vector{Int64}}}, Tuple{BlockedUnitRange{Vector{Int64}}, BlockedUnitRange{Vector{Int64}}}}:
Block(1, 1)[1, 1] Block(1, 1)[1, 2] │ Block(1, 2)[1, 1] Block(1, 2)[1, 2]
Block(1, 1)[2, 1] Block(1, 1)[2, 2] │ Block(1, 2)[2, 1] Block(1, 2)[2, 2]
──────────────────────────────────────┼──────────────────────────────────────
Block(2, 1)[1, 1] Block(2, 1)[1, 2] │ Block(2, 2)[1, 1] Block(2, 2)[1, 2]
Block(2, 1)[2, 1] Block(2, 1)[2, 2] │ Block(2, 2)[2, 1] Block(2, 2)[2, 2]
julia> @view B[Base.IdentityUnitRange(2:4), Base.IdentityUnitRange(2:4)]
1×1-blocked 3×3 BlockIndices{2, Tuple{BlockedUnitRange{Vector{Int64}}, BlockedUnitRange{Vector{Int64}}}, Tuple{Base.IdentityUnitRange{UnitRange{Int64}}, Base.IdentityUnitRange{UnitRange{Int64}}}}:
Block(2, 2)[1, 1] Block(2, 2)[1, 2] #undef
Block(2, 2)[2, 1] Block(2, 2)[2, 2] #undef
#undef #undef #undef |
Yeah, I didn't look into indexing arrays with Good point about offset axes, I also didn't have that use case in mind but hopefully is easy to support. |
In the latest I added support for indexing with julia> using BlockArrays
julia> using BlockArrays: BlockIndices
julia> A = BlockArray(randn(4, 4), [2, 2], [2, 2])
2×2-blocked 4×4 BlockMatrix{Float64}:
-2.27028 1.18245 │ -0.0100978 0.219669
0.534015 -2.69366 │ 0.131658 -0.131192
─────────────────────────┼───────────────────────
0.0449745 -0.288047 │ 1.71181 0.76368
1.27463 -0.00159753 │ 0.0061755 -0.427065
julia> r = BlockArrays._BlockedUnitRange(2, [2, 3])
2-blocked 2-element BlockedUnitRange{Vector{Int64}}:
2
─
3
julia> B = BlockIndices(A)[r, r]
2×2-blocked 2×2 BlockIndices{2, Tuple{BlockedUnitRange{Vector{Int64}}, BlockedUnitRange{Vector{Int64}}}, Tuple{BlockedUnitRange{Vector{Int64}}, BlockedUnitRange{Vector{Int64}}}}:
Block(1, 1)[2, 2] │ Block(1, 2)[2, 1]
───────────────────┼───────────────────
Block(2, 1)[1, 2] │ Block(2, 2)[1, 1]
julia> CartesianIndices(A)[B]
CartesianIndices((2:1:3, 2:1:3))
julia> A[B]
2×2-blocked 2×2 BlockMatrix{Float64}:
-2.69366 │ 0.131658
───────────┼──────────
-0.288047 │ 1.71181 though I think I still need to check that the |
This is an implementation of #346.
TODO:
AbstractArray
withBlockIndices
(i.e.getindex(a::AbstractArray, indices::BlockIndices)
) check that the block indices are compatible, say withblockisequal
.view(a::CartesianIndices, indices::BlockIndices)
.getindex(a::AbstractArray, indices::BlockIndices{1}...)
.BlockIndices
with offset axes ([WIP]BlockIndices
#356 (comment)).Here is a demonstration of what works so far:
This is the current slicing behavior:
Ideally, any slicing operation that can get mapped to slicing with a unit range in each dimension (say slicing with
AbstractUnitRange
,BlockRange
,BlockIndexRange
, etc.) would outputBlockIndices
orBlockIndexRange
, but that isn't implemented yet. I think ideally that would be implemented as something like this:i.e. slicing a
BlockIndices
slices the indices in each dimension, similar to the definition forBase.CartesianIndices
: https://github.com/JuliaLang/julia/blob/64de065a183ac70bb049f7f9e30d790f8845dd2b/base/multidimensional.jl#L384-L391. However, slicing blocked unit ranges doesn't always preserve blocking information right now, see for example #347 and #355.