Incorporate blocks into axis, support offset blocks

.travis.yml

@@ -4,7 +4,6 @@ os:
   - linux
   - osx
 julia:
-  - 1.0
   - 1.1
   - 1.2
   - 1.3
@@ -25,7 +24,7 @@ after_success:
 jobs:
   include:
     - stage: "Documentation"
-      julia: 1.0
+      julia: 1.3
       os: linux
       script:
         - julia --project=docs/ --color=yes -e '

Project.toml

@@ -1,17 +1,19 @@
 name = "BlockArrays"
 uuid = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
-version = "0.10.2"
+version = "0.11"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 
 [compat]
-julia = "1"
+julia = "1.1"
 
 [extras]
 Base64 = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
+FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
+OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Base64", "Test", "SparseArrays"]
+test = ["Base64", "FillArrays", "OffsetArrays", "SparseArrays", "Test"]

docs/src/lib/internals.md

@@ -19,8 +19,7 @@ Pages = ["internals.md"]
 ## Internals
 
 ```@docs
-blockindex2global
-global2blockindex
+BlockedUnitRange
 BlockRange
 BlockIndexRange
 BlockSlice

docs/src/lib/public.md

@@ -30,9 +30,12 @@ AbstractBlockArray
 BlockBoundsError
 Block
 BlockIndex
-nblocks
+blockaxes
+blockisequal
 blocksize
-blocksizes
+blockfirsts
+blocklasts
+blocklengths
 getblock
 getblock!
 setblock!

docs/src/man/abstractblockarrayinterface.md

@@ -1,29 +1,53 @@
-# The `AbstractBlockSizes` interface
+# The block axis interface
 
-In order to follow the `AbstractBlockSizes` the following methods should be implemented:
+A block array's block structure is dictated by its axes, which
+are typically a `BlockedUnitRange` but may also be a `UnitRange`, 
+which is assumed to be a single block, or other type that implements
+the block axis interface.
 
 
 | Methods to implement    | Brief description |
 | :---------------------- | :---------------- |
-| `cumulsizes(A)`      | A Tuple of abstract vectors storing the cumulative block sizes |
-| **Optional methods**    |       
-| `nblocks(A)`            | Tuple of number of blocks in each dimension |
-| `nblocks(A, i)`         | Number of blocks in dimension `i` |
-| `blocksize(A, i)`    | Size of the block at block index `i` |
+| `blockaxes(A)`      | A one-tuple returning a range of blocks specifying the block structure |
+| `getindex(A, K::Block{1})`      | return a unit range of indices in the specified block |
+| `blocklasts(A)`      | Returns the last index of each block |
+| `findblock(A, k)`      | return the block that contains the `k`th entry of `A` 
+
 
 # The `AbstractBlockArray` interface
 
+An arrays block structure is inferred from an axes, and therefore every array
+is in some sense already a block array:
+```julia
+julia> A = randn(5,5)
+5×5 Array{Float64,2}:
+  0.452801   -0.416508   1.17406    1.52575     3.1574  
+  0.413142   -1.34722   -1.28597    0.637721    0.30655 
+  0.34907    -0.887615   0.284972  -0.0212884  -0.225832
+  0.466102   -1.10425    1.49226    0.968436   -2.13637 
+ -0.0971956  -1.7664    -0.592629  -1.48947     1.53418 
+
+julia> A[Block(1,1)]
+5×5 Array{Float64,2}:
+  0.452801   -0.416508   1.17406    1.52575     3.1574  
+  0.413142   -1.34722   -1.28597    0.637721    0.30655 
+  0.34907    -0.887615   0.284972  -0.0212884  -0.225832
+  0.466102   -1.10425    1.49226    0.968436   -2.13637 
+ -0.0971956  -1.7664    -0.592629  -1.48947     1.53418 
+```
+It is possible to override additional functions to improve speed, however.
+
 | Methods to implement    | Brief description |
 | :---------------------- | :---------------- |
-| `blocksizes(A)`         | Return a subtype of `AbstractBlockSizes` |
-| **Optional methods**    |                       
+| **Optional methods**    |           
 | `getblock(A, i...)`     | `X[Block(i...)]`, blocked indexing  |
 | `setblock!(A, v, i...)` | `X[Block(i...)] = v`, blocked index assignment |
 | `getblock!(x, A, i)`    | `X[i]`, blocked index assignment with in place storage in `x` |
 
 For a more thorough description of the methods see the public interface documentation.
 
-With the methods above implemented the following are automatically provided:
+With the methods above implemented the following are automatically provided for arrays
+that are subtypes of `AbstractBlockArray`:
 
 * A pretty printing `show` function that uses unicode lines to split up the blocks:
 ```

docs/src/man/blockarrays.md

@@ -22,7 +22,7 @@ julia> BlockArray(rand(4, 4), [2,2], [1,1,2])
  0.844314  │  0.794279  │  0.0421491  0.683791
 
 julia> block_array_sparse = BlockArray(sprand(4, 5, 0.7), [1,3], [2,3])
-2×2-blocked 4×5 BlockArray{Float64,2,Array{SparseMatrixCSC{Float64,Int64},2},BlockArrays.BlockSizes{2,Array{Int64,1}}}:
+2×2-blocked 4×5 BlockArray{Float64,2,Array{SparseMatrixCSC{Float64,Int64},2},Tuple{BlockedUnitRange{Array{Int64,1}},BlockedUnitRange{Array{Int64,1}}}}:
  0.0341601  0.374187  │  0.0118196  0.299058  0.0     
  ---------------------┼-------------------------------
  0.0945445  0.931115  │  0.0460428  0.0       0.0     
@@ -67,7 +67,7 @@ The `block_type` should be an array type.  It specifies the internal block type,
 
 ```julia
 julia> BlockArray(undef_blocks, SparseVector{Float64, Int}, [1,2])
-2-blocked 3-element BlockArray{Float64,1,Array{SparseVector{Float64,Int64},1},BlockArrays.BlockSizes{1,Tuple{Array{Int64,1}}}}:
+2-blocked 3-element BlockArray{Float64,1,Array{SparseVector{Float64,Int64},1},Tuple{BlockedUnitRange{Array{Int64,1}}}}:
  #undef
  ------
  #undef
@@ -146,7 +146,7 @@ We can also view and modify views of blocks of `BlockArray` using the `view` syn
 julia> A = BlockArray(ones(6), 1:3);
 
 julia> view(A, Block(2))
-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:
+2-element view(::BlockArray{Float64,1,Array{Array{Float64,1},1},Tuple{BlockedUnitRange{Array{Int64,1}}}}, BlockSlice(Block(2),2:3)) with eltype Float64:
  1.0
  1.0
 
@@ -164,7 +164,7 @@ An array can be repacked into a `BlockArray` with `BlockArray(array, block_sizes
 
 ```jl
 julia> block_array_sparse = BlockArray(sprand(4, 5, 0.7), [1,3], [2,3])
-2×2-blocked 4×5 BlockArray{Float64,2,Array{SparseMatrixCSC{Float64,Int64},2},BlockArrays.BlockSizes{2,Tuple{Array{Int64,1}}}}:
+2×2-blocked 4×5 BlockArray{Float64,2,Array{SparseMatrixCSC{Float64,Int64},2},Tuple{BlockedUnitRange{Array{Int64,1}},BlockedUnitRange{Array{Int64,1}}}}:
  0.0341601  0.374187  │  0.0118196  0.299058  0.0     
  ---------------------┼-------------------------------
  0.0945445  0.931115  │  0.0460428  0.0       0.0     

docs/src/man/pseudoblockarrays.md

@@ -83,7 +83,7 @@ We can also view and modify views of blocks of `PseudoBlockArray` using the `vie
 julia> A = PseudoBlockArray(ones(6), 1:3);
 
 julia> view(A, Block(2))
-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:
+2-element view(::PseudoBlockArray{Float64,1,Array{Float64,1},Tuple{BlockedUnitRange{Array{Int64,1}}}}, BlockSlice(Block(2),2:3)) with eltype Float64:
  1.0
  1.0
 

src/BlockArrays.jl

@@ -1,23 +1,24 @@
-__precompile__()
-
 module BlockArrays
 using Base.Cartesian
 using LinearAlgebra
 
 # AbstractBlockArray interface exports
 export AbstractBlockArray, AbstractBlockMatrix, AbstractBlockVector, AbstractBlockVecOrMat
-export Block, getblock, getblock!, setblock!, nblocks, blocksize, blocksizes, blockcheckbounds, BlockBoundsError, BlockIndex
-export BlockRange
+export Block, getblock, getblock!, setblock!
+export blockaxes, blocksize, blocklength, blockcheckbounds, BlockBoundsError, BlockIndex
+export blocklengths, blocklasts, blockfirsts, blockisequal
+export BlockRange, blockedrange, BlockedUnitRange
 
 export BlockArray, BlockMatrix, BlockVector, BlockVecOrMat, mortar
 export PseudoBlockArray, PseudoBlockMatrix, PseudoBlockVector, PseudoBlockVecOrMat
 
-export undef_blocks, undef
+export undef_blocks, undef, findblock, findblockindex
 
 import Base: @propagate_inbounds, Array, to_indices, to_index,
             unsafe_indices, first, last, size, length, unsafe_length,
             unsafe_convert,
             getindex, show,
+            step, 
             broadcast, eltype, convert, similar,
             @_inline_meta, _maybetail, tail, @_propagate_inbounds_meta, reindex,
             RangeIndex, Int, Integer, Number,
@@ -26,22 +27,18 @@ import Base: @propagate_inbounds, Array, to_indices, to_index,
 using Base: ReshapedArray, dataids
 
 
-
 import Base: (:), IteratorSize, iterate, axes1
 import Base.Broadcast: broadcasted, DefaultArrayStyle, AbstractArrayStyle, Broadcasted
 import LinearAlgebra: lmul!, rmul!, AbstractTriangular, HermOrSym, AdjOrTrans,
                         StructuredMatrixStyle
 
 
+include("blockindices.jl")                        
+include("blockaxis.jl")
 include("abstractblockarray.jl")
-
-include("blocksizes.jl")
-include("blockindices.jl")
 include("blockarray.jl")
 include("pseudo_blockarray.jl")
-include("blockrange.jl")
 include("views.jl")
-include("blockindexrange.jl")
 include("show.jl")
 include("blockarrayinterface.jl")
 include("blockbroadcast.jl")

src/abstractblockarray.jl

@@ -22,109 +22,26 @@ const AbstractBlockVector{T} = AbstractBlockArray{T, 1}
 const AbstractBlockVecOrMat{T} = Union{AbstractBlockMatrix{T}, AbstractBlockVector{T}}
 
 block2string(b, s) = string(join(map(string,b), '×'), "-blocked ", Base.dims2string(s))
-Base.summary(a::AbstractBlockArray) = string(block2string(nblocks(a), size(a)), " ", typeof(a))
+_block_summary(a) = string(block2string(blocksize(a), size(a)), " ", typeof(a))
+Base.summary(a::AbstractBlockArray) = _block_summary(a)
 _show_typeof(io, a) = show(io, typeof(a))
-function Base.summary(io::IO, a::AbstractBlockArray)
-    print(io, block2string(nblocks(a), size(a)))
+function _block_summary(io, a)    
+    print(io, block2string(blocksize(a), size(a)))
     print(io, ' ')
     _show_typeof(io, a)
 end
-Base.similar(block_array::AbstractBlockArray{T}) where {T} = similar(block_array, T)
-Base.IndexStyle(::Type{<:AbstractBlockArray}) = IndexCartesian()
-
-"""
-    nblocks(A, [dim...])
-
-Returns a tuple containing the number of blocks in a block array.  Optionally you can specify
-the dimension(s) you want the number of blocks for.
-
-```jldoctest; setup = quote using BlockArrays end
-julia> A =  BlockArray(rand(5,4,6), [1,4], [1,2,1], [1,2,2,1]);
-
-julia> nblocks(A)
-(2, 3, 4)
-
-julia> nblocks(A, 2)
-3
-
-julia> nblocks(A, 3, 2)
-(4, 3)
-```
-"""
-nblocks(block_array::AbstractArray, i::Integer) = nblocks(block_array)[i]
-
-nblocks(block_array::AbstractArray, i::Vararg{Integer, N}) where {N} =
-    nblocks(blocksizes(block_array), i...)
-
-
-"""
-    Block(inds...)
-
-A `Block` is simply a wrapper around a set of indices or enums so that it can be used to dispatch on. By
-indexing a `AbstractBlockArray` with a `Block` the a block at that block index will be returned instead of
-a single element.
-
-```jldoctest; setup = quote using BlockArrays end
-julia> A = BlockArray(ones(2,3), [1, 1], [2, 1])
-2×2-blocked 2×3 BlockArray{Float64,2}:
- 1.0  1.0  │  1.0
- ──────────┼─────
- 1.0  1.0  │  1.0
-
-julia> A[Block(1, 1)]
-1×2 Array{Float64,2}:
- 1.0  1.0
-```
-"""
-struct Block{N, T}
-    n::NTuple{N, T}
-    Block{N, T}(n::NTuple{N, T}) where {N, T} = new{N, T}(n)
-end
+Base.summary(io::IO, a::AbstractBlockArray) = _block_summary(io, a)
 
+# avoid to_shape which complicates axes
+Base.similar(a::AbstractBlockArray{T}) where {T}                             = similar(a, T)
+Base.similar(a::AbstractBlockArray, ::Type{T}) where {T}                     = similar(a, T, axes(a))
+Base.similar(a::AbstractBlockArray{T}, dims::Tuple) where {T}                = similar(a, T, dims)
 
-Block{N, T}(n::Vararg{T, N}) where {N,T} = Block{N, T}(n)
-Block{N}(n::Vararg{T, N}) where {N,T} = Block{N, T}(n)
-Block() = Block{0,Int}()
-Block(n::Vararg{T, N}) where {N,T} = Block{N, T}(n)
-Block{1}(n::Tuple{T}) where {T} = Block{1, T}(n)
-Block{N}(n::NTuple{N, T}) where {N,T} = Block{N, T}(n)
-Block(n::NTuple{N, T}) where {N,T} = Block{N, T}(n)
-
-@inline function Block(blocks::NTuple{N, Block{1, T}}) where {N,T}
-    Block{N, T}(ntuple(i -> blocks[i].n[1], Val(N)))
-end
-
-
-# The following code is taken from CartesianIndex
-@inline (+)(index::Block{N}) where {N} = Block{N}(map(+, index.n))
-@inline (-)(index::Block{N}) where {N} = Block{N}(map(-, index.n))
-
-@inline (+)(index1::Block{N}, index2::Block{N}) where {N} =
-    Block{N}(map(+, index1.n, index2.n))
-@inline (-)(index1::Block{N}, index2::Block{N}) where {N} =
-    Block{N}(map(-, index1.n, index2.n))
-@inline min(index1::Block{N}, index2::Block{N}) where {N} =
-    Block{N}(map(min, index1.n, index2.n))
-@inline max(index1::Block{N}, index2::Block{N}) where {N} =
-    Block{N}(map(max, index1.n, index2.n))
-
-@inline (+)(i::Integer, index::Block) = index+i
-@inline (+)(index::Block{N}, i::Integer) where {N} = Block{N}(map(x->x+i, index.n))
-@inline (-)(index::Block{N}, i::Integer) where {N} = Block{N}(map(x->x-i, index.n))
-@inline (-)(i::Integer, index::Block{N}) where {N} = Block{N}(map(x->i-x, index.n))
-@inline (*)(a::Integer, index::Block{N}) where {N} = Block{N}(map(x->a*x, index.n))
-@inline (*)(index::Block, a::Integer) = *(a,index)
-
-# comparison
-@inline isless(I1::Block{N}, I2::Block{N}) where {N} = Base.IteratorsMD._isless(0, I1.n, I2.n)
-
-# conversions
-convert(::Type{T}, index::Block{1}) where {T<:Number} = convert(T, index.n[1])
-convert(::Type{T}, index::Block) where {T<:Tuple} = convert(T, index.n)
+Base.IndexStyle(::Type{<:AbstractBlockArray}) = IndexCartesian()
 
-Int(index::Block{1}) = Int(index.n[1])
-Integer(index::Block{1}) = index.n[1]
-Number(index::Block{1}) = index.n[1]
+# need to overload axes to return BlockAxis
+@inline size(block_array::AbstractBlockArray) = map(length, axes(block_array))
+@inline axes(block_array::AbstractBlockArray) = throw(error("axes for ", typeof(block_array), " is not implemented"))
 
 """
     getblock(A, inds...)
@@ -222,10 +139,12 @@ struct BlockBoundsError <: Exception
     a::Any
     i::Any
     BlockBoundsError() = new()
-    BlockBoundsError(a::AbstractBlockArray) = new(a)
-    BlockBoundsError(a::AbstractBlockArray, @nospecialize(i)) = new(a,i)
+    BlockBoundsError(a::AbstractArray) = new(a)
+    BlockBoundsError(a::AbstractArray, @nospecialize(i)) = new(a,i)
 end
 
+BlockBoundsError(a::AbstractArray, I::Block) = BlockBoundsError(a, I.n)
+
 function Base.showerror(io::IO, ex::BlockBoundsError)
     print(io, "BlockBoundsError")
     if isdefined(ex, :a)
@@ -250,7 +169,7 @@ specialize this method if they need to provide custom block bounds checking beha
 julia> A = BlockArray(rand(2,3), [1,1], [2,1]);
 
 julia> blockcheckbounds(A, 3, 2)
-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]
+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]
 [...]
 ```
 """
@@ -263,14 +182,12 @@ ERROR: BlockBoundsError: attempt to access 2×2-blocked 2×3 BlockArray{Float64,
 end
 
 @inline function blockcheckbounds(::Type{Bool}, A::AbstractArray{T, N}, i::Vararg{Integer, N}) where {T,N}
-    n = nblocks(A)
+    n = blockaxes(A)
     k = 0
     for idx in 1:N # using enumerate here will allocate
         k += 1
         @inbounds _i = i[idx]
-        if _i <= 0 || _i > n[k]
-            return false
-        end
+        Block(_i) in n[k] || return false
     end
     return true
 end