diff --git a/src/BlockArraysExtensions/BlockArraysExtensions.jl b/src/BlockArraysExtensions/BlockArraysExtensions.jl index 46d015a..46250ea 100644 --- a/src/BlockArraysExtensions/BlockArraysExtensions.jl +++ b/src/BlockArraysExtensions/BlockArraysExtensions.jl @@ -56,104 +56,6 @@ BlockArrays.blocks(a::NonBlockedArray) = SingleBlockView(parent(a)) const NonBlockedVector{T,Array} = NonBlockedArray{T,1,Array} NonBlockedVector(array::AbstractVector) = NonBlockedArray(array) -# BlockIndices works around an issue that the indices of BlockSlice -# are restricted to AbstractUnitRange{Int}. -struct BlockIndices{B,T<:Integer,I<:AbstractVector{T}} <: AbstractVector{T} - blocks::B - indices::I -end -for f in (:axes, :unsafe_indices, :axes1, :first, :last, :size, :length, :unsafe_length) - @eval Base.$f(S::BlockIndices) = Base.$f(S.indices) -end -Base.getindex(S::BlockIndices, i::Integer) = getindex(S.indices, i) -function Base.getindex(S::BlockIndices, i::BlockSlice{<:Block{1}}) - # TODO: Check that `i.indices` is consistent with `S.indices`. - # It seems like this isn't handling the case where `i` is a - # subslice of a block correctly (i.e. it ignores `i.indices`). - @assert length(S.indices[Block(i)]) == length(i.indices) - return BlockSlice(S.blocks[Int(Block(i))], S.indices[Block(i)]) -end - -# This is used in slicing like: -# a = BlockSparseArray{Float64}([2, 2, 2, 2], [2, 2, 2, 2]) -# I = BlockedVector([Block(4), Block(3), Block(2), Block(1)], [2, 2]) -# a[I, I] -function Base.getindex( - S::BlockIndices{<:AbstractBlockVector{<:Block{1}}}, i::BlockSlice{<:Block{1}} -) - # TODO: Check for conistency of indices. - # Wrapping the indices in `NonBlockedVector` reinterprets the blocked indices - # as a single block, since the result shouldn't be blocked. - return NonBlockedVector(BlockIndices(S.blocks[Block(i)], S.indices[Block(i)])) -end -function Base.getindex( - S::BlockIndices{<:BlockedVector{<:Block{1},<:BlockRange{1}}}, i::BlockSlice{<:Block{1}} -) - return i -end -# Views of `BlockIndices` are eager. -# This fixes an issue in Julia 1.11 where reindexing defaults to using views. -Base.view(S::BlockIndices, i) = S[i] - -# Used in indexing such as: -# ```julia -# a = BlockSparseArray{Float64}([2, 2, 2, 2], [2, 2, 2, 2]) -# I = BlockedVector([Block(4), Block(3), Block(2), Block(1)], [2, 2]) -# b = @view a[I, I] -# @view b[Block(1, 1)[1:2, 2:2]] -# ``` -# This is similar to the definition: -# @interface BlockSparseArrayInterface() to_indices(a, inds, I::Tuple{UnitRange{<:Integer},Vararg{Any}}) -function Base.getindex( - a::NonBlockedVector{<:Integer,<:BlockIndices}, I::UnitRange{<:Integer} -) - ax = only(axes(parent(a).indices)) - brs = to_blockindices(ax, I) - inds = blockedunitrange_getindices(ax, I) - return NonBlockedVector(parent(a)[BlockSlice(brs, inds)]) -end - -function Base.getindex(S::BlockIndices, i::BlockSlice{<:BlockRange{1}}) - # TODO: Check that `i.indices` is consistent with `S.indices`. - # TODO: Turn this into a `blockedunitrange_getindices` definition. - subblocks = S.blocks[Int.(i.block)] - subindices = mortar( - map(1:length(i.block)) do I - r = blocks(i.indices)[I] - return S.indices[first(r)]:S.indices[last(r)] - end, - ) - return BlockIndices(subblocks, subindices) -end - -# Used when performing slices like: -# @views a[[Block(2), Block(1)]][2:4, 2:4] -function Base.getindex(S::BlockIndices, i::BlockSlice{<:BlockVector{<:BlockIndex{1}}}) - subblocks = mortar( - map(blocks(i.block)) do br - return S.blocks[Int(Block(br))][only(br.indices)] - end, - ) - subindices = mortar( - map(blocks(i.block)) do br - S.indices[br] - end, - ) - return BlockIndices(subblocks, subindices) -end - -# Similar to the definition of `BlockArrays.BlockSlices`: -# ```julia -# const BlockSlices = Union{Base.Slice,BlockSlice{<:BlockRange{1}}} -# ``` -# but includes `BlockIndices`, where the blocks aren't contiguous. -const BlockSliceCollection = Union{ - Base.Slice,BlockSlice{<:BlockRange{1}},BlockIndices{<:Vector{<:Block{1}}} -} -const SubBlockSliceCollection = BlockIndices{ - <:BlockVector{<:BlockIndex{1},<:Vector{<:BlockIndexRange{1}}} -} - # TODO: This is type piracy. This is used in `reindex` when making # views of blocks of sliced block arrays, for example: # ```julia @@ -326,10 +228,6 @@ function blockrange(axis::AbstractUnitRange, r::BlockSlice) return blockrange(axis, r.block) end -function blockrange(a::AbstractUnitRange, r::BlockIndices) - return blockrange(a, r.blocks) -end - function blockrange(axis::AbstractUnitRange, r::Block{1}) return r:r end @@ -360,6 +258,22 @@ function blockrange(axis::AbstractUnitRange, r::NonBlockedVector) return Block(1):Block(1) end +# Find the blocks associated with the block slice. +# This results from slices like `[Block(1), Block(2)]` +# or `[Block(1)[1:2], Block(2)[2:3]]`. +# This would be easier once something like: +# https://github.com/JuliaArrays/BlockArrays.jl/pull/459 +# is implemented. +function blockrange(axis::AbstractUnitRange, r::AbstractBlockVector) + inds = map(blocks(r)) do b_r + findfirst(blocks(axis)) do b_axis + return b_r ⊆ b_axis + end + end + any(isnothing, inds) && error("Block not found.") + return Block.(inds) +end + function blockrange(axis::AbstractUnitRange, r) return error("Slicing not implemented for range of type `$(typeof(r))`.") end @@ -397,10 +311,6 @@ function blockindices( return Base.OneTo(length(axis[block])) end -function blockindices(a::AbstractUnitRange, b::Block, r::BlockIndices) - return blockindices(a, b, r.blocks) -end - function blockindices( a::AbstractUnitRange, b::Block, diff --git a/src/abstractblocksparsearray/map.jl b/src/abstractblocksparsearray/map.jl index 724b0ff..d1d5a1c 100644 --- a/src/abstractblocksparsearray/map.jl +++ b/src/abstractblocksparsearray/map.jl @@ -22,19 +22,6 @@ function reblock( return @view parent(a)[map(I -> I.array, parentindices(a))...] end -# TODO: Make this more general, independent of `AbstractBlockSparseArray`. -function reblock( - a::SubArray{ - <:Any, - <:Any, - <:AbstractBlockSparseArray, - <:Tuple{Vararg{BlockIndices{<:AbstractBlockVector{<:Block{1}}}}}, - }, -) - # Remove the blocking. - return @view parent(a)[map(I -> Vector(I.blocks), parentindices(a))...] -end - function Base.map!(f, a_dest::AbstractArray, a_srcs::AnyAbstractBlockSparseArray...) @interface interface(a_dest, a_srcs...) map!(f, a_dest, a_srcs...) return a_dest diff --git a/src/abstractblocksparsearray/views.jl b/src/abstractblocksparsearray/views.jl index f245198..d847536 100644 --- a/src/abstractblocksparsearray/views.jl +++ b/src/abstractblocksparsearray/views.jl @@ -119,18 +119,6 @@ function BlockArrays.viewblock( return Base.invoke(view, Tuple{AbstractArray,Vararg{Any}}, a, I...) end -function Base.view( - a::SubArray{ - T, - N, - <:AbstractBlockSparseArray{T,N}, - <:Tuple{Vararg{Union{BlockSliceCollection,SubBlockSliceCollection},N}}, - }, - block::Union{Block{N},BlockIndexRange{N}}, -) where {T,N} - return viewblock(a, block) -end - # Fix ambiguity error with BlockArrays.jl for slices like # `a = BlockSparseArray{Float64}(undef, [2, 2], [2, 2]); @view a[:, :]`. function Base.view( @@ -145,29 +133,6 @@ function Base.view( return viewblock(a, block) end -function Base.view( - a::SubArray{ - T, - N, - <:AbstractBlockSparseArray{T,N}, - <:Tuple{Vararg{Union{BlockSliceCollection,SubBlockSliceCollection},N}}, - }, - block::Vararg{Union{Block{1},BlockIndexRange{1}},N}, -) where {T,N} - return viewblock(a, block...) -end -function BlockArrays.viewblock( - a::SubArray{ - T, - N, - <:AbstractBlockSparseArray{T,N}, - <:Tuple{Vararg{Union{BlockSliceCollection,SubBlockSliceCollection},N}}, - }, - block::Union{Block{N},BlockIndexRange{N}}, -) where {T,N} - return viewblock(a, to_tuple(block)...) -end - # Fixes ambiguity error with `AnyAbstractBlockSparseArray` definition. function Base.view( a::SubArray{ @@ -187,31 +152,6 @@ function Base.view( return viewblock(a, block...) end -# XXX: TODO: Distinguish if a sub-view of the block needs to be taken! -# Define a new `SubBlockSlice` which is used in: -# `@interface BlockSparseArrayInterface() to_indices(a, inds, I::Tuple{UnitRange{<:Integer},Vararg{Any}})` -# in `blocksparsearrayinterface/blocksparsearrayinterface.jl`. -# TODO: Define `@interface BlockSparseArrayInterface() viewblock`. -function BlockArrays.viewblock( - a::SubArray{T,N,<:AbstractBlockSparseArray{T,N},<:Tuple{Vararg{BlockSliceCollection,N}}}, - block::Vararg{Block{1},N}, -) where {T,N} - I = CartesianIndex(Int.(block)) - # TODO: Use `eachblockstoredindex`. - if I ∈ eachstoredindex(blocks(a)) - return blocks(a)[I] - end - return BlockView(parent(a), Block.(Base.reindex(parentindices(blocks(a)), Tuple(I)))) -end - -function to_blockindexrange( - a::BlockIndices{<:BlockVector{<:BlockIndex{1},<:Vector{<:BlockIndexRange{1}}}}, - I::Block{1}, -) - # TODO: Ideally we would just use `a.blocks[I]` but that doesn't - # work right now. - return blocks(a.blocks)[Int(I)] -end function to_blockindexrange( a::Base.Slice{<:AbstractBlockedUnitRange{<:Integer}}, I::Block{1} ) @@ -219,32 +159,6 @@ function to_blockindexrange( return I end -function BlockArrays.viewblock( - a::SubArray{ - T, - N, - <:AbstractBlockSparseArray{T,N}, - <:Tuple{Vararg{Union{BlockSliceCollection,SubBlockSliceCollection},N}}, - }, - block::Vararg{Block{1},N}, -) where {T,N} - brs = ntuple(dim -> to_blockindexrange(parentindices(a)[dim], block[dim]), ndims(a)) - return @view parent(a)[brs...] -end - -# TODO: Define `@interface BlockSparseArrayInterface() viewblock`. -function BlockArrays.viewblock( - a::SubArray{ - T, - N, - <:AbstractBlockSparseArray{T,N}, - <:Tuple{Vararg{Union{BlockSliceCollection,SubBlockSliceCollection},N}}, - }, - block::Vararg{BlockIndexRange{1},N}, -) where {T,N} - return view(viewblock(a, Block.(block)...), map(b -> only(b.indices), block)...) -end - # Block slice of the result of slicing `@view a[2:5, 2:5]`. # TODO: Move this to `BlockArraysExtensions`. const BlockedSlice = BlockSlice{ diff --git a/src/abstractblocksparsearray/wrappedabstractblocksparsearray.jl b/src/abstractblocksparsearray/wrappedabstractblocksparsearray.jl index 0403b78..647484d 100644 --- a/src/abstractblocksparsearray/wrappedabstractblocksparsearray.jl +++ b/src/abstractblocksparsearray/wrappedabstractblocksparsearray.jl @@ -40,13 +40,6 @@ function Base.to_indices( return @interface BlockSparseArrayInterface() to_indices(a, inds, I) end -# a[[Block(2), Block(1)], [Block(2), Block(1)]] -function Base.to_indices( - a::AnyAbstractBlockSparseArray, inds, I::Tuple{Vector{<:Block{1}},Vararg{Any}} -) - return @interface BlockSparseArrayInterface() to_indices(a, inds, I) -end - # a[BlockVector([Block(2), Block(1)], [2]), BlockVector([Block(2), Block(1)], [2])] # a[BlockedVector([Block(2), Block(1)], [2]), BlockedVector([Block(2), Block(1)], [2])] function Base.to_indices( @@ -66,13 +59,6 @@ function Base.to_indices( return @interface BlockSparseArrayInterface() to_indices(a, inds, I) end -# a[[Block(1)[1:2], Block(2)[1:2]], [Block(1)[1:2], Block(2)[1:2]]] -function Base.to_indices( - a::AnyAbstractBlockSparseArray, inds, I::Tuple{Vector{<:BlockIndexRange{1}},Vararg{Any}} -) - return to_indices(a, inds, (mortar(I[1]), Base.tail(I)...)) -end - # BlockArrays `AbstractBlockArray` interface function BlockArrays.blocks(a::AnyAbstractBlockSparseArray) @interface BlockSparseArrayInterface() blocks(a) diff --git a/src/blocksparsearrayinterface/blocksparsearrayinterface.jl b/src/blocksparsearrayinterface/blocksparsearrayinterface.jl index 1b1af5f..608a0fb 100644 --- a/src/blocksparsearrayinterface/blocksparsearrayinterface.jl +++ b/src/blocksparsearrayinterface/blocksparsearrayinterface.jl @@ -145,33 +145,6 @@ end return (inds[1][I[1]], to_indices(a, Base.tail(inds), Base.tail(I))...) end -# a[[Block(2), Block(1)], [Block(2), Block(1)]] -@interface ::AbstractBlockSparseArrayInterface function Base.to_indices( - a, inds, I::Tuple{Vector{<:Block{1}},Vararg{Any}} -) - I1 = BlockIndices(I[1], blockedunitrange_getindices(inds[1], I[1])) - return (I1, to_indices(a, Base.tail(inds), Base.tail(I))...) -end - -# a[mortar([Block(1)[1:2], Block(2)[1:3]]), mortar([Block(1)[1:2], Block(2)[1:3]])] -# a[[Block(1)[1:2], Block(2)[1:3]], [Block(1)[1:2], Block(2)[1:3]]] -@interface ::AbstractBlockSparseArrayInterface function Base.to_indices( - a, inds, I::Tuple{BlockVector{<:BlockIndex{1},<:Vector{<:BlockIndexRange{1}}},Vararg{Any}} -) - I1 = BlockIndices(I[1], blockedunitrange_getindices(inds[1], I[1])) - return (I1, to_indices(a, Base.tail(inds), Base.tail(I))...) -end - -# a[BlockVector([Block(2), Block(1)], [2]), BlockVector([Block(2), Block(1)], [2])] -# Permute and merge blocks. -# TODO: This isn't merging blocks yet, that needs to be implemented that. -@interface ::AbstractBlockSparseArrayInterface function Base.to_indices( - a, inds, I::Tuple{AbstractBlockVector{<:Block{1}},Vararg{Any}} -) - I1 = BlockIndices(I[1], blockedunitrange_getindices(inds[1], I[1])) - return (I1, to_indices(a, Base.tail(inds), Base.tail(I))...) -end - # TODO: Need to implement this! function block_merge end @@ -424,14 +397,7 @@ function SparseArraysBase.getunstoredindex( end to_blocks_indices(I::BlockSlice{<:BlockRange{1}}) = Int.(I.block) -to_blocks_indices(I::BlockIndices{<:Vector{<:Block{1}}}) = Int.(I.blocks) to_blocks_indices(I::Base.Slice{<:BlockedOneTo}) = Base.OneTo(blocklength(I.indices)) -@interface ::AbstractBlockSparseArrayInterface function BlockArrays.blocks( - a::SubArray{<:Any,<:Any,<:Any,<:Tuple{Vararg{BlockSliceCollection}}} -) - return @view blocks(parent(a))[map(to_blocks_indices, parentindices(a))...] -end - using BlockArrays: BlocksView SparseArraysBase.storedlength(a::BlocksView) = length(a)