Skip to content

Commit 3addf8f

Browse files
dlfivefiftytkf
andauthored
Incorporate blocks into axis, support offset blocks (#95)
* Create BlockAxis * use searchsortedfirst for simple block lookup * Update blockaxis.jl * BlockAxis working * start deleting BlockSizes * BlockIndex with BlockAxis * basic BlockArray works * update PseudoBlockArray * tests starting to pass * work on views * improvements to view * blockaxes returns a BlockRange * tests pass!(?) * broadcast now partially working * Fix broadcast blocksize * test_blockarrays now passes * fix blockaxis BlockRange indexing * minor fixes * Redesign as CumsumBlockRange * Block-views now working * Most tests pass * Add reshape * tests pass! * Add show for CumsumBlockRange * Fixes from BlockBandedMatrices * Improve similar * CumsumBlockRange constructor renamed blockedrange to avoid ambiguity with conversion * use sortedunion with special overrides * Facilitate changing BroadcastStyle * fix tests * Drop Julia v1.0 * Update .travis.yml * update docs * add show * update docs * Update blockaxis.jl * Update internals.md * Update src/blockbroadcast.jl Co-Authored-By: Takafumi Arakaki <[email protected]> * _block_cumsum -> blocklasts * Update test_blockarrayinterface.jl * add docs, test blockfirst and convert * add blockisequal * increase coverage Co-authored-by: Takafumi Arakaki <[email protected]>
1 parent b4e9812 commit 3addf8f

26 files changed

+1309
-1127
lines changed

.travis.yml

Lines changed: 1 addition & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -4,7 +4,6 @@ os:
44
- linux
55
- osx
66
julia:
7-
- 1.0
87
- 1.1
98
- 1.2
109
- 1.3
@@ -25,7 +24,7 @@ after_success:
2524
jobs:
2625
include:
2726
- stage: "Documentation"
28-
julia: 1.0
27+
julia: 1.3
2928
os: linux
3029
script:
3130
- julia --project=docs/ --color=yes -e '

Project.toml

Lines changed: 5 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1,17 +1,19 @@
11
name = "BlockArrays"
22
uuid = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
3-
version = "0.10.2"
3+
version = "0.11"
44

55
[deps]
66
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
77

88
[compat]
9-
julia = "1"
9+
julia = "1.1"
1010

1111
[extras]
1212
Base64 = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
13+
FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
14+
OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
1315
SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
1416
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
1517

1618
[targets]
17-
test = ["Base64", "Test", "SparseArrays"]
19+
test = ["Base64", "FillArrays", "OffsetArrays", "SparseArrays", "Test"]

docs/src/lib/internals.md

Lines changed: 1 addition & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -19,8 +19,7 @@ Pages = ["internals.md"]
1919
## Internals
2020

2121
```@docs
22-
blockindex2global
23-
global2blockindex
22+
BlockedUnitRange
2423
BlockRange
2524
BlockIndexRange
2625
BlockSlice

docs/src/lib/public.md

Lines changed: 5 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -30,9 +30,12 @@ AbstractBlockArray
3030
BlockBoundsError
3131
Block
3232
BlockIndex
33-
nblocks
33+
blockaxes
34+
blockisequal
3435
blocksize
35-
blocksizes
36+
blockfirsts
37+
blocklasts
38+
blocklengths
3639
getblock
3740
getblock!
3841
setblock!

docs/src/man/abstractblockarrayinterface.md

Lines changed: 34 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -1,29 +1,53 @@
1-
# The `AbstractBlockSizes` interface
1+
# The block axis interface
22

3-
In order to follow the `AbstractBlockSizes` the following methods should be implemented:
3+
A block array's block structure is dictated by its axes, which
4+
are typically a `BlockedUnitRange` but may also be a `UnitRange`,
5+
which is assumed to be a single block, or other type that implements
6+
the block axis interface.
47

58

69
| Methods to implement | Brief description |
710
| :---------------------- | :---------------- |
8-
| `cumulsizes(A)` | A Tuple of abstract vectors storing the cumulative block sizes |
9-
| **Optional methods** |
10-
| `nblocks(A)` | Tuple of number of blocks in each dimension |
11-
| `nblocks(A, i)` | Number of blocks in dimension `i` |
12-
| `blocksize(A, i)` | Size of the block at block index `i` |
11+
| `blockaxes(A)` | A one-tuple returning a range of blocks specifying the block structure |
12+
| `getindex(A, K::Block{1})` | return a unit range of indices in the specified block |
13+
| `blocklasts(A)` | Returns the last index of each block |
14+
| `findblock(A, k)` | return the block that contains the `k`th entry of `A`
15+
1316

1417
# The `AbstractBlockArray` interface
1518

19+
An arrays block structure is inferred from an axes, and therefore every array
20+
is in some sense already a block array:
21+
```julia
22+
julia> A = randn(5,5)
23+
5×5 Array{Float64,2}:
24+
0.452801 -0.416508 1.17406 1.52575 3.1574
25+
0.413142 -1.34722 -1.28597 0.637721 0.30655
26+
0.34907 -0.887615 0.284972 -0.0212884 -0.225832
27+
0.466102 -1.10425 1.49226 0.968436 -2.13637
28+
-0.0971956 -1.7664 -0.592629 -1.48947 1.53418
29+
30+
julia> A[Block(1,1)]
31+
5×5 Array{Float64,2}:
32+
0.452801 -0.416508 1.17406 1.52575 3.1574
33+
0.413142 -1.34722 -1.28597 0.637721 0.30655
34+
0.34907 -0.887615 0.284972 -0.0212884 -0.225832
35+
0.466102 -1.10425 1.49226 0.968436 -2.13637
36+
-0.0971956 -1.7664 -0.592629 -1.48947 1.53418
37+
```
38+
It is possible to override additional functions to improve speed, however.
39+
1640
| Methods to implement | Brief description |
1741
| :---------------------- | :---------------- |
18-
| `blocksizes(A)` | Return a subtype of `AbstractBlockSizes` |
19-
| **Optional methods** |
42+
| **Optional methods** |
2043
| `getblock(A, i...)` | `X[Block(i...)]`, blocked indexing |
2144
| `setblock!(A, v, i...)` | `X[Block(i...)] = v`, blocked index assignment |
2245
| `getblock!(x, A, i)` | `X[i]`, blocked index assignment with in place storage in `x` |
2346

2447
For a more thorough description of the methods see the public interface documentation.
2548

26-
With the methods above implemented the following are automatically provided:
49+
With the methods above implemented the following are automatically provided for arrays
50+
that are subtypes of `AbstractBlockArray`:
2751

2852
* A pretty printing `show` function that uses unicode lines to split up the blocks:
2953
```

docs/src/man/blockarrays.md

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -22,7 +22,7 @@ julia> BlockArray(rand(4, 4), [2,2], [1,1,2])
2222
0.8443140.7942790.0421491 0.683791
2323

2424
julia> block_array_sparse = BlockArray(sprand(4, 5, 0.7), [1,3], [2,3])
25-
2×2-blocked 4×5 BlockArray{Float64,2,Array{SparseMatrixCSC{Float64,Int64},2},BlockArrays.BlockSizes{2,Array{Int64,1}}}:
25+
2×2-blocked 4×5 BlockArray{Float64,2,Array{SparseMatrixCSC{Float64,Int64},2},Tuple{BlockedUnitRange{Array{Int64,1}},BlockedUnitRange{Array{Int64,1}}}}:
2626
0.0341601 0.3741870.0118196 0.299058 0.0
2727
----------------------------------------------------
2828
0.0945445 0.9311150.0460428 0.0 0.0
@@ -67,7 +67,7 @@ The `block_type` should be an array type. It specifies the internal block type,
6767

6868
```julia
6969
julia> BlockArray(undef_blocks, SparseVector{Float64, Int}, [1,2])
70-
2-blocked 3-element BlockArray{Float64,1,Array{SparseVector{Float64,Int64},1},BlockArrays.BlockSizes{1,Tuple{Array{Int64,1}}}}:
70+
2-blocked 3-element BlockArray{Float64,1,Array{SparseVector{Float64,Int64},1},Tuple{BlockedUnitRange{Array{Int64,1}}}}:
7171
#undef
7272
------
7373
#undef
@@ -146,7 +146,7 @@ We can also view and modify views of blocks of `BlockArray` using the `view` syn
146146
julia> A = BlockArray(ones(6), 1:3);
147147
148148
julia> view(A, Block(2))
149-
2-element view(::BlockArray{Float64,1,Array{Array{Float64,1},1},BlockArrays.BlockSizes{1,Tuple{Array{Int64,1}}}}, BlockSlice(Block{1,Int64}((2,)),2:3)) with eltype Float64:
149+
2-element view(::BlockArray{Float64,1,Array{Array{Float64,1},1},Tuple{BlockedUnitRange{Array{Int64,1}}}}, BlockSlice(Block(2),2:3)) with eltype Float64:
150150
1.0
151151
1.0
152152
@@ -164,7 +164,7 @@ An array can be repacked into a `BlockArray` with `BlockArray(array, block_sizes
164164

165165
```jl
166166
julia> block_array_sparse = BlockArray(sprand(4, 5, 0.7), [1,3], [2,3])
167-
2×2-blocked 4×5 BlockArray{Float64,2,Array{SparseMatrixCSC{Float64,Int64},2},BlockArrays.BlockSizes{2,Tuple{Array{Int64,1}}}}:
167+
2×2-blocked 4×5 BlockArray{Float64,2,Array{SparseMatrixCSC{Float64,Int64},2},Tuple{BlockedUnitRange{Array{Int64,1}},BlockedUnitRange{Array{Int64,1}}}}:
168168
0.0341601 0.3741870.0118196 0.299058 0.0
169169
----------------------------------------------------
170170
0.0945445 0.9311150.0460428 0.0 0.0

docs/src/man/pseudoblockarrays.md

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -83,7 +83,7 @@ We can also view and modify views of blocks of `PseudoBlockArray` using the `vie
8383
julia> A = PseudoBlockArray(ones(6), 1:3);
8484
8585
julia> view(A, Block(2))
86-
2-element view(::PseudoBlockArray{Float64,1,Array{Float64,1},BlockArrays.BlockSizes{1,Tuple{Array{Int64,1}}}}, BlockSlice(Block{1,Int64}((2,)),2:3)) with eltype Float64:
86+
2-element view(::PseudoBlockArray{Float64,1,Array{Float64,1},Tuple{BlockedUnitRange{Array{Int64,1}}}}, BlockSlice(Block(2),2:3)) with eltype Float64:
8787
1.0
8888
1.0
8989

src/BlockArrays.jl

Lines changed: 8 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -1,23 +1,24 @@
1-
__precompile__()
2-
31
module BlockArrays
42
using Base.Cartesian
53
using LinearAlgebra
64

75
# AbstractBlockArray interface exports
86
export AbstractBlockArray, AbstractBlockMatrix, AbstractBlockVector, AbstractBlockVecOrMat
9-
export Block, getblock, getblock!, setblock!, nblocks, blocksize, blocksizes, blockcheckbounds, BlockBoundsError, BlockIndex
10-
export BlockRange
7+
export Block, getblock, getblock!, setblock!
8+
export blockaxes, blocksize, blocklength, blockcheckbounds, BlockBoundsError, BlockIndex
9+
export blocklengths, blocklasts, blockfirsts, blockisequal
10+
export BlockRange, blockedrange, BlockedUnitRange
1111

1212
export BlockArray, BlockMatrix, BlockVector, BlockVecOrMat, mortar
1313
export PseudoBlockArray, PseudoBlockMatrix, PseudoBlockVector, PseudoBlockVecOrMat
1414

15-
export undef_blocks, undef
15+
export undef_blocks, undef, findblock, findblockindex
1616

1717
import Base: @propagate_inbounds, Array, to_indices, to_index,
1818
unsafe_indices, first, last, size, length, unsafe_length,
1919
unsafe_convert,
2020
getindex, show,
21+
step,
2122
broadcast, eltype, convert, similar,
2223
@_inline_meta, _maybetail, tail, @_propagate_inbounds_meta, reindex,
2324
RangeIndex, Int, Integer, Number,
@@ -26,22 +27,18 @@ import Base: @propagate_inbounds, Array, to_indices, to_index,
2627
using Base: ReshapedArray, dataids
2728

2829

29-
3030
import Base: (:), IteratorSize, iterate, axes1
3131
import Base.Broadcast: broadcasted, DefaultArrayStyle, AbstractArrayStyle, Broadcasted
3232
import LinearAlgebra: lmul!, rmul!, AbstractTriangular, HermOrSym, AdjOrTrans,
3333
StructuredMatrixStyle
3434

3535

36+
include("blockindices.jl")
37+
include("blockaxis.jl")
3638
include("abstractblockarray.jl")
37-
38-
include("blocksizes.jl")
39-
include("blockindices.jl")
4039
include("blockarray.jl")
4140
include("pseudo_blockarray.jl")
42-
include("blockrange.jl")
4341
include("views.jl")
44-
include("blockindexrange.jl")
4542
include("show.jl")
4643
include("blockarrayinterface.jl")
4744
include("blockbroadcast.jl")

src/abstractblockarray.jl

Lines changed: 20 additions & 103 deletions
Original file line numberDiff line numberDiff line change
@@ -22,109 +22,26 @@ const AbstractBlockVector{T} = AbstractBlockArray{T, 1}
2222
const AbstractBlockVecOrMat{T} = Union{AbstractBlockMatrix{T}, AbstractBlockVector{T}}
2323

2424
block2string(b, s) = string(join(map(string,b), '×'), "-blocked ", Base.dims2string(s))
25-
Base.summary(a::AbstractBlockArray) = string(block2string(nblocks(a), size(a)), " ", typeof(a))
25+
_block_summary(a) = string(block2string(blocksize(a), size(a)), " ", typeof(a))
26+
Base.summary(a::AbstractBlockArray) = _block_summary(a)
2627
_show_typeof(io, a) = show(io, typeof(a))
27-
function Base.summary(io::IO, a::AbstractBlockArray)
28-
print(io, block2string(nblocks(a), size(a)))
28+
function _block_summary(io, a)
29+
print(io, block2string(blocksize(a), size(a)))
2930
print(io, ' ')
3031
_show_typeof(io, a)
3132
end
32-
Base.similar(block_array::AbstractBlockArray{T}) where {T} = similar(block_array, T)
33-
Base.IndexStyle(::Type{<:AbstractBlockArray}) = IndexCartesian()
34-
35-
"""
36-
nblocks(A, [dim...])
37-
38-
Returns a tuple containing the number of blocks in a block array. Optionally you can specify
39-
the dimension(s) you want the number of blocks for.
40-
41-
```jldoctest; setup = quote using BlockArrays end
42-
julia> A = BlockArray(rand(5,4,6), [1,4], [1,2,1], [1,2,2,1]);
43-
44-
julia> nblocks(A)
45-
(2, 3, 4)
46-
47-
julia> nblocks(A, 2)
48-
3
49-
50-
julia> nblocks(A, 3, 2)
51-
(4, 3)
52-
```
53-
"""
54-
nblocks(block_array::AbstractArray, i::Integer) = nblocks(block_array)[i]
55-
56-
nblocks(block_array::AbstractArray, i::Vararg{Integer, N}) where {N} =
57-
nblocks(blocksizes(block_array), i...)
58-
59-
60-
"""
61-
Block(inds...)
62-
63-
A `Block` is simply a wrapper around a set of indices or enums so that it can be used to dispatch on. By
64-
indexing a `AbstractBlockArray` with a `Block` the a block at that block index will be returned instead of
65-
a single element.
66-
67-
```jldoctest; setup = quote using BlockArrays end
68-
julia> A = BlockArray(ones(2,3), [1, 1], [2, 1])
69-
2×2-blocked 2×3 BlockArray{Float64,2}:
70-
1.0 1.0 │ 1.0
71-
──────────┼─────
72-
1.0 1.0 │ 1.0
73-
74-
julia> A[Block(1, 1)]
75-
1×2 Array{Float64,2}:
76-
1.0 1.0
77-
```
78-
"""
79-
struct Block{N, T}
80-
n::NTuple{N, T}
81-
Block{N, T}(n::NTuple{N, T}) where {N, T} = new{N, T}(n)
82-
end
33+
Base.summary(io::IO, a::AbstractBlockArray) = _block_summary(io, a)
8334

35+
# avoid to_shape which complicates axes
36+
Base.similar(a::AbstractBlockArray{T}) where {T} = similar(a, T)
37+
Base.similar(a::AbstractBlockArray, ::Type{T}) where {T} = similar(a, T, axes(a))
38+
Base.similar(a::AbstractBlockArray{T}, dims::Tuple) where {T} = similar(a, T, dims)
8439

85-
Block{N, T}(n::Vararg{T, N}) where {N,T} = Block{N, T}(n)
86-
Block{N}(n::Vararg{T, N}) where {N,T} = Block{N, T}(n)
87-
Block() = Block{0,Int}()
88-
Block(n::Vararg{T, N}) where {N,T} = Block{N, T}(n)
89-
Block{1}(n::Tuple{T}) where {T} = Block{1, T}(n)
90-
Block{N}(n::NTuple{N, T}) where {N,T} = Block{N, T}(n)
91-
Block(n::NTuple{N, T}) where {N,T} = Block{N, T}(n)
92-
93-
@inline function Block(blocks::NTuple{N, Block{1, T}}) where {N,T}
94-
Block{N, T}(ntuple(i -> blocks[i].n[1], Val(N)))
95-
end
96-
97-
98-
# The following code is taken from CartesianIndex
99-
@inline (+)(index::Block{N}) where {N} = Block{N}(map(+, index.n))
100-
@inline (-)(index::Block{N}) where {N} = Block{N}(map(-, index.n))
101-
102-
@inline (+)(index1::Block{N}, index2::Block{N}) where {N} =
103-
Block{N}(map(+, index1.n, index2.n))
104-
@inline (-)(index1::Block{N}, index2::Block{N}) where {N} =
105-
Block{N}(map(-, index1.n, index2.n))
106-
@inline min(index1::Block{N}, index2::Block{N}) where {N} =
107-
Block{N}(map(min, index1.n, index2.n))
108-
@inline max(index1::Block{N}, index2::Block{N}) where {N} =
109-
Block{N}(map(max, index1.n, index2.n))
110-
111-
@inline (+)(i::Integer, index::Block) = index+i
112-
@inline (+)(index::Block{N}, i::Integer) where {N} = Block{N}(map(x->x+i, index.n))
113-
@inline (-)(index::Block{N}, i::Integer) where {N} = Block{N}(map(x->x-i, index.n))
114-
@inline (-)(i::Integer, index::Block{N}) where {N} = Block{N}(map(x->i-x, index.n))
115-
@inline (*)(a::Integer, index::Block{N}) where {N} = Block{N}(map(x->a*x, index.n))
116-
@inline (*)(index::Block, a::Integer) = *(a,index)
117-
118-
# comparison
119-
@inline isless(I1::Block{N}, I2::Block{N}) where {N} = Base.IteratorsMD._isless(0, I1.n, I2.n)
120-
121-
# conversions
122-
convert(::Type{T}, index::Block{1}) where {T<:Number} = convert(T, index.n[1])
123-
convert(::Type{T}, index::Block) where {T<:Tuple} = convert(T, index.n)
40+
Base.IndexStyle(::Type{<:AbstractBlockArray}) = IndexCartesian()
12441

125-
Int(index::Block{1}) = Int(index.n[1])
126-
Integer(index::Block{1}) = index.n[1]
127-
Number(index::Block{1}) = index.n[1]
42+
# need to overload axes to return BlockAxis
43+
@inline size(block_array::AbstractBlockArray) = map(length, axes(block_array))
44+
@inline axes(block_array::AbstractBlockArray) = throw(error("axes for ", typeof(block_array), " is not implemented"))
12845

12946
"""
13047
getblock(A, inds...)
@@ -222,10 +139,12 @@ struct BlockBoundsError <: Exception
222139
a::Any
223140
i::Any
224141
BlockBoundsError() = new()
225-
BlockBoundsError(a::AbstractBlockArray) = new(a)
226-
BlockBoundsError(a::AbstractBlockArray, @nospecialize(i)) = new(a,i)
142+
BlockBoundsError(a::AbstractArray) = new(a)
143+
BlockBoundsError(a::AbstractArray, @nospecialize(i)) = new(a,i)
227144
end
228145

146+
BlockBoundsError(a::AbstractArray, I::Block) = BlockBoundsError(a, I.n)
147+
229148
function Base.showerror(io::IO, ex::BlockBoundsError)
230149
print(io, "BlockBoundsError")
231150
if isdefined(ex, :a)
@@ -250,7 +169,7 @@ specialize this method if they need to provide custom block bounds checking beha
250169
julia> A = BlockArray(rand(2,3), [1,1], [2,1]);
251170
252171
julia> blockcheckbounds(A, 3, 2)
253-
ERROR: BlockBoundsError: attempt to access 2×2-blocked 2×3 BlockArray{Float64,2,Array{Array{Float64,2},2},BlockArrays.BlockSizes{2,Tuple{Array{Int64,1},Array{Int64,1}}}} at block index [3,2]
172+
ERROR: BlockBoundsError: attempt to access 2×2-blocked 2×3 BlockArray{Float64,2,Array{Array{Float64,2},2},Tuple{BlockedUnitRange{Array{Int64,1}},BlockedUnitRange{Array{Int64,1}}}} at block index [3,2]
254173
[...]
255174
```
256175
"""
@@ -263,14 +182,12 @@ ERROR: BlockBoundsError: attempt to access 2×2-blocked 2×3 BlockArray{Float64,
263182
end
264183

265184
@inline function blockcheckbounds(::Type{Bool}, A::AbstractArray{T, N}, i::Vararg{Integer, N}) where {T,N}
266-
n = nblocks(A)
185+
n = blockaxes(A)
267186
k = 0
268187
for idx in 1:N # using enumerate here will allocate
269188
k += 1
270189
@inbounds _i = i[idx]
271-
if _i <= 0 || _i > n[k]
272-
return false
273-
end
190+
Block(_i) in n[k] || return false
274191
end
275192
return true
276193
end

0 commit comments

Comments
 (0)