diff --git a/NEWS.md b/NEWS.md index b13d5b6831d08..7af6632eb04aa 100644 --- a/NEWS.md +++ b/NEWS.md @@ -27,6 +27,11 @@ New language features * Function return type syntax `function f()::T` has been added ([#1090]). Values returned from a function with such a declaration will be converted to the specified type `T`. + * Experimental support for arrays with indexing starting at values + different from 1. The array types are expected to be defined in + packages, but now Julia provides an API for writing generic + algorithms for arbitrary indexing schemes ([#16260]). + Language changes ---------------- @@ -130,7 +135,7 @@ Library improvements * `cov` and `cor` don't use keyword arguments anymore and are therefore now type stable ([#13465]). - * Linear algebra: + * Arrays and linear algebra: * All dimensions indexed by scalars are now dropped, whereas previously only trailing scalar dimensions would be omitted from the result ([#13612]). @@ -268,3 +273,4 @@ Deprecated or removed [#16403]: https://github.com/JuliaLang/julia/issues/16403 [#16481]: https://github.com/JuliaLang/julia/issues/16481 [#16731]: https://github.com/JuliaLang/julia/issues/16731 +[#16280]: https://github.com/JuliaLang/julia/issues/16260 diff --git a/base/abstractarray.jl b/base/abstractarray.jl index ad3ea1605be48..c75ad4bb41d3f 100644 --- a/base/abstractarray.jl +++ b/base/abstractarray.jl @@ -6,6 +6,8 @@ typealias AbstractVector{T} AbstractArray{T,1} typealias AbstractMatrix{T} AbstractArray{T,2} typealias AbstractVecOrMat{T} Union{AbstractVector{T}, AbstractMatrix{T}} typealias RangeIndex Union{Int, Range{Int}, UnitRange{Int}, Colon} +typealias UnitRangeInteger{T<:Integer} UnitRange{T} +typealias Indices{N} NTuple{N,UnitRangeInteger} ## Basic functions ## @@ -21,6 +23,46 @@ end size{T,N}(t::AbstractArray{T,N}, d) = d <= N ? size(t)[d] : 1 size{N}(x, d1::Integer, d2::Integer, dx::Vararg{Integer, N}) = (size(x, d1), size(x, d2), ntuple(k->size(x, dx[k]), Val{N})...) +""" + indices(A, d) + +Returns the valid range of indices for array `A` along dimension `d`. +""" +function indices(A::AbstractArray, d) + @_inline_meta + 1:size(A,d) +end +""" + indices(A) + +Returns the tuple of valid indices for array `A`. +""" +indices{T,N}(A::AbstractArray{T,N}) = _indices((), A) +_indices{T,N}(out::NTuple{N}, A::AbstractArray{T,N}) = out +function _indices(out, A::AbstractArray) + @_inline_meta + _indices((out..., indices(A, length(out)+1)), A) +end +# This simpler implementation suffers from #16327 +# function indices{T,N}(A::AbstractArray{T,N}) +# @_inline_meta +# ntuple(d->indices(A, d), Val{N}) +# end +indices1(A) = (@_inline_meta; indices(A, 1)) +""" + linearindices(A) + +Returns a `UnitRange` specifying the valid range of indices for `A[i]` +where `i` is an `Int`. For arrays with conventional indexing (indices +start at 1), or any multidimensional array, this is `1:length(A)`; +however, for one-dimensional arrays with unconventional indices, this +is `indices(A, 1)`. + +Calling this function is the "safe" way to write algorithms that +exploit linear indexing. +""" +linearindices(A) = 1:length(A) +linearindices(A::AbstractVector) = indices1(A) eltype{T}(::Type{AbstractArray{T}}) = T eltype{T,N}(::Type{AbstractArray{T,N}}) = T elsize{T}(::AbstractArray{T}) = sizeof(T) @@ -86,8 +128,48 @@ linearindexing(A::AbstractArray, B::AbstractArray...) = linearindexing(linearind linearindexing(::LinearFast, ::LinearFast) = LinearFast() linearindexing(::LinearIndexing, ::LinearIndexing) = LinearSlow() +abstract IndicesBehavior +immutable IndicesStartAt1 <: IndicesBehavior end # indices 1:size(A,d) +immutable IndicesUnitRange <: IndicesBehavior end # arb UnitRange indices +immutable IndicesList <: IndicesBehavior end # indices like (:cat, :dog, :mouse) + +indicesbehavior(A::AbstractArray) = indicesbehavior(typeof(A)) +indicesbehavior{T<:AbstractArray}(::Type{T}) = IndicesStartAt1() +indicesbehavior(::Number) = IndicesStartAt1() + +abstract IndicesPerformance +immutable IndicesFast1D <: IndicesPerformance end # indices(A, d) is fast +immutable IndicesSlow1D <: IndicesPerformance end # indices(A) is better than indices(A,d) + +indicesperformance(A::AbstractArray) = indicesperformance(typeof(A)) +indicesperformance{T<:AbstractArray}(::Type{T}) = IndicesFast1D() + +""" + shape(A) + +Returns a tuple specifying the "shape" of array `A`. For arrays with +conventional indexing (indices start at 1), this is equivalent to +`size(A)`; otherwise it is equivalent to `incides(A)`. +""" +shape(a) = shape(indicesbehavior(a), a) +""" + shape(A, d) + +Specifies the "shape" of the array `A` along dimension `d`. For arrays +with conventional indexing (starting at 1), this is equivalent to +`size(A, d)`; for arrays with unconventional indexing (indexing may +start at something different from 1), it is equivalent to `indices(A, +d)`. +""" +shape(a, d) = shape(indicesbehavior(a), a, d) +shape(::IndicesStartAt1, a) = size(a) +shape(::IndicesStartAt1, a, d) = size(a, d) +shape(::IndicesBehavior, a) = indices(a) +shape(::IndicesBehavior, a, d) = indices(a, d) + ## Bounds checking ## @generated function trailingsize{T,N,n}(A::AbstractArray{T,N}, ::Type{Val{n}}) + (isa(n, Int) && isa(N, Int)) || error("Must have concrete type") n > N && return 1 ex = :(size(A, $n)) for m = n+1:N @@ -97,116 +179,276 @@ linearindexing(::LinearIndexing, ::LinearIndexing) = LinearSlow() end # check along a single dimension -checkbounds(::Type{Bool}, sz::Integer, i) = throw(ArgumentError("unable to check bounds for indices of type $(typeof(i))")) -checkbounds(::Type{Bool}, sz::Integer, i::Real) = 1 <= i <= sz -checkbounds(::Type{Bool}, sz::Integer, ::Colon) = true -function checkbounds(::Type{Bool}, sz::Integer, r::Range) +""" + checkindex(Bool, inds::UnitRange, index) + +Return `true` if the given `index` is within the bounds of +`inds`. Custom types that would like to behave as indices for all +arrays can extend this method in order to provide a specialized bounds +checking implementation. +""" +checkindex(::Type{Bool}, inds::UnitRange, i) = throw(ArgumentError("unable to check bounds for indices of type $(typeof(i))")) +checkindex(::Type{Bool}, inds::UnitRange, i::Real) = (first(inds) <= i) & (i <= last(inds)) +checkindex(::Type{Bool}, inds::UnitRange, ::Colon) = true +function checkindex(::Type{Bool}, inds::UnitRange, r::Range) @_propagate_inbounds_meta - isempty(r) | (checkbounds(Bool, sz, first(r)) & checkbounds(Bool, sz, last(r))) + isempty(r) | (checkindex(Bool, inds, first(r)) & checkindex(Bool, inds, last(r))) end -checkbounds{N}(::Type{Bool}, sz::Integer, I::AbstractArray{Bool,N}) = N == 1 && length(I) == sz -function checkbounds(::Type{Bool}, sz::Integer, I::AbstractArray) +checkindex{N}(::Type{Bool}, indx::UnitRange, I::AbstractArray{Bool,N}) = N == 1 && indx == indices1(I) +function checkindex(::Type{Bool}, inds::UnitRange, I::AbstractArray) @_inline_meta b = true for i in I - b &= checkbounds(Bool, sz, i) + b &= checkindex(Bool, inds, i) end b end -# check all dimensions -function checkbounds{N,T}(::Type{Bool}, sz::NTuple{N,Integer}, I1::T, I...) +# check all indices/dimensions +# To make extension easier, avoid specializing checkbounds on index types +""" + checkbounds(Bool, array, indexes...) + +Return `true` if the specified `indexes` are in bounds for the given `array`. Subtypes of +`AbstractArray` should specialize this method if they need to provide custom bounds checking +behaviors. +""" +function checkbounds(::Type{Bool}, A::AbstractArray, I...) + @_inline_meta + _checkbounds(Bool, indicesperformance(A), A, I...) +end +function _checkbounds(::Type{Bool}, ::IndicesSlow1D, A::AbstractArray, I...) @_inline_meta - checkbounds(Bool, sz[1], I1) & checkbounds(Bool, tail(sz), I...) + checkbounds_indices(indices(A), I) end -checkbounds{T<:Integer}(::Type{Bool}, sz::Tuple{T}, I1) = (@_inline_meta; checkbounds(Bool, sz[1], I1)) -checkbounds{N}(::Type{Bool}, sz::NTuple{N,Integer}, I1) = (@_inline_meta; checkbounds(Bool, prod(sz), I1)) -checkbounds{N}(::Type{Bool}, sz::NTuple{N,Integer}) = (@_inline_meta; checkbounds(Bool, sz, 1)) # for a[] +_checkbounds(::Type{Bool}, ::IndicesSlow1D, A::AbstractArray, I::AbstractArray{Bool}) = _chkbnds(A, (true,), I) + +# Bounds-checking for arrays for which indices(A) is faster than indices(A, d) +checkbounds_indices(::Tuple{}, ::Tuple{}) = true +checkbounds_indices(::Tuple{}, I::Tuple{Any}) = (@_inline_meta; checkindex(Bool, 1:1, I[1])) +checkbounds_indices(::Tuple{}, I::Tuple) = (@_inline_meta; checkindex(Bool, 1:1, I[1]) & checkbounds_indices((), tail(I))) +checkbounds_indices(inds::Tuple{Any}, I::Tuple{Any}) = (@_inline_meta; checkindex(Bool, inds[1], I[1])) +checkbounds_indices(inds::Tuple, I::Tuple{Any}) = (@_inline_meta; checkindex(Bool, 1:prod(map(dimlength, inds)), I[1])) +checkbounds_indices(inds::Tuple, I::Tuple) = (@_inline_meta; checkindex(Bool, inds[1], I[1]) & checkbounds_indices(tail(inds), tail(I))) -checkbounds(::Type{Bool}, sz::Tuple{}, i) = (@_inline_meta; checkbounds(Bool, 1, i)) -function checkbounds(::Type{Bool}, sz::Tuple{}, i, I...) +# Bounds-checking for arrays for which indices(A, d) is fast +function _checkbounds(::Type{Bool}, ::IndicesFast1D, A::AbstractArray, I...) @_inline_meta - checkbounds(Bool, 1, i) & checkbounds(Bool, (), I...) + # checked::NTuple{M} means we have checked dimensions 1:M-1, now + # need to check dimension M. checked[M] indicates whether all the + # previous ones are in-bounds. + # By growing checked, it allows us to test whether we've processed + # the same number of dimensions as the array, even while + # supporting CartesianIndex + _chkbnds(A, (true,), I...) +end +checkbounds(::Type{Bool}, A::AbstractArray) = checkbounds(Bool, A, 1) # 0-d case +# Single logical array indexing: +_chkbnds(A::AbstractArray, ::NTuple{1,Bool}, I::AbstractArray{Bool}) = indices(A) == indices(I) +_chkbnds(A::AbstractArray, ::NTuple{1,Bool}, I::AbstractVector{Bool}) = length(A) == length(I) +_chkbnds(A::AbstractVector, ::NTuple{1,Bool}, I::AbstractArray{Bool}) = length(A) == length(I) +_chkbnds(A::AbstractVector, ::NTuple{1,Bool}, I::AbstractVector{Bool}) = indices(A) == indices(I) +# Linear indexing: +function _chkbnds(A::AbstractVector, ::NTuple{1,Bool}, I) + @_inline_meta + checkindex(Bool, indices1(A), I) +end +function _chkbnds(A::AbstractArray, ::NTuple{1,Bool}, I) + @_inline_meta + checkindex(Bool, 1:length(A), I) +end +# When all indices have been checked: +_chkbnds{M}(A, checked::NTuple{M,Bool}) = checked[M] +# When the number of indices matches the array dimensionality: +function _chkbnds{T,N}(A::AbstractArray{T,N}, checked::NTuple{N,Bool}, I1) + @_inline_meta + checked[N] & checkindex(Bool, indices(A, N), I1) +end +# When the last checked dimension is not equal to the array dimensionality: +# TODO: for #14770 (deprecating partial linear indexing), change to 1:trailingsize(...) to 1:1 +function _chkbnds{T,N,M}(A::AbstractArray{T,N}, checked::NTuple{M,Bool}, I1) + @_inline_meta + checked[M] & checkindex(Bool, 1:trailingsize(A, Val{M}), I1) +end +# Checking an interior dimension: +function _chkbnds{T,N,M}(A::AbstractArray{T,N}, checked::NTuple{M,Bool}, I1, I...) + @_inline_meta + # grow checked by one + newchecked = (checked..., checked[M] & checkindex(Bool, indices(A, M), I1)) + _chkbnds(A, newchecked, I...) end -# Prevent allocation of a GC frame by hiding the BoundsError in a noinline function + throw_boundserror(A, I) = (@_noinline_meta; throw(BoundsError(A, I))) -# Don't define index types on checkbounds to make extending easier -checkbounds(A::AbstractArray, I...) = (@_inline_meta; _internal_checkbounds(A, I...)) -# The internal function is named _internal_checkbounds since there had been a -# _checkbounds previously that meant something different. -_internal_checkbounds(A::AbstractArray) = _internal_checkbounds(A,1) -_internal_checkbounds(A::AbstractArray, I::AbstractArray{Bool}) = size(A) == size(I) || throw_boundserror(A, I) -_internal_checkbounds(A::AbstractArray, I::AbstractVector{Bool}) = length(A) == length(I) || throw_boundserror(A, I) -function _internal_checkbounds(A::AbstractArray, I1, I...) - # having I1 seems important for good codegen +""" + checkbounds(array, indexes...) + +Throw an error if the specified `indexes` are not in bounds for the given `array`. +""" +function checkbounds(A::AbstractArray, I...) @_inline_meta - checkbounds(Bool, size(A), I1, I...) || throw_boundserror(A, (I1, I...)) + checkbounds(Bool, A, I...) || throw_boundserror(A, I) end +checkbounds(A::AbstractArray) = checkbounds(A, 1) # 0-d case # See also specializations in multidimensional ## Constructors ## # default arguments to similar() +typealias SimIdx Union{Integer,UnitRangeInteger} +""" + similar(array, [element_type=eltype(array)], [dims=size(array)]) + +Create an uninitialized mutable array with the given element type and size, based upon the +given source array. The second and third arguments are both optional, defaulting to the +given array's `eltype` and `size`. The dimensions may be specified either as a single tuple +argument or as a series of integer arguments. + +Custom AbstractArray subtypes may choose which specific array type is best-suited to return +for the given element type and dimensionality. If they do not specialize this method, the +default is an `Array{element_type}(dims...)`. + +For example, `similar(1:10, 1, 4)` returns an uninitialized `Array{Int,2}` since ranges are +neither mutable nor support 2 dimensions: + + julia> similar(1:10, 1, 4) + 1×4 Array{Int64,2}: + 4419743872 4374413872 4419743888 0 + +Conversely, `similar(trues(10,10), 2)` returns an uninitialized `BitVector` with two +elements since `BitArray`s are both mutable and can support 1-dimensional arrays: + + julia> similar(trues(10,10), 2) + 2-element BitArray{1}: + false + false + +Since `BitArray`s can only store elements of type `Bool`, however, if you request a +different element type it will create a regular `Array` instead: + + julia> similar(falses(10), Float64, 2, 4) + 2×4 Array{Float64,2}: + 2.18425e-314 2.18425e-314 2.18425e-314 2.18425e-314 + 2.18425e-314 2.18425e-314 2.18425e-314 2.18425e-314 + +See also `allocate_for`. +""" similar{T}(a::AbstractArray{T}) = similar(a, T) -similar( a::AbstractArray, T::Type) = similar(a, T, size(a)) -similar{T}(a::AbstractArray{T}, dims::DimsInteger) = similar(a, T, dims) -similar{T}(a::AbstractArray{T}, dims::Integer...) = similar(a, T, dims) -similar( a::AbstractArray, T::Type, dims::Integer...) = similar(a, T, dims) -# similar creates an Array by default +similar( a::AbstractArray, T::Type) = _similar(indicesbehavior(a), a, T) +similar{T}(a::AbstractArray{T}, dims::Tuple) = similar(a, T, dims) +similar{T}(a::AbstractArray{T}, dims::SimIdx...) = similar(a, T, dims) +similar( a::AbstractArray, T::Type, dims::SimIdx...) = similar(a, T, dims) similar( a::AbstractArray, T::Type, dims::DimsInteger) = similar(a, T, convert(Dims, dims)) -similar( a::AbstractArray, T::Type, dims::Dims) = Array{T}(dims) +# similar creates an Array by default +similar( a::AbstractArray, T::Type, dims::Dims) = Array(T, dims) + +_similar(::IndicesStartAt1, a::AbstractArray, T::Type) = similar(a, T, size(a)) +_similar(::IndicesBehavior, a::AbstractArray, T::Type) = similar(a, T, indices(a)) + +""" + allocate_for(storagetype, referencearray, [shape]) + +Create an uninitialized mutable array analogous to that specified by +`storagetype`, but with type and shape specified by the final two +arguments. The main purpose of this function is to support allocation +of arrays that may have unconventional indexing (starting at other +than 1), as determined by `referencearray` and the optional `shape` +information. + + **Examples**: + + allocate_for(Array{Int}, A) + +creates an array that "acts like" an `Array{Int}` (and might indeed be +backed by one), but which is indexed identically to `A`. If `A` has +conventional indexing, this will likely just call +`Array{Int}(size(A))`, but if `A` has unconventional indexing then the +indices of the result will match `A`. + + allocate_for(BitArray, A, (shape(A, 2),)) + +would create a 1-dimensional logical array whose indices match those +of the columns of `A`. + +The main purpose of the `referencearray` argument is to select a +particular array type supporting unconventional indexing (as it is +possible that several different ones will be simultaneously in use). + +See also `similar`. +""" +allocate_for(f, a, shape::Union{SimIdx,Tuple{Vararg{SimIdx}}}) = f(shape) +allocate_for(f, a) = allocate_for(f, a, shape(a)) +# allocate_for when passed multiple arrays. Necessary for broadcast, etc. +function allocate_for(f, as::Tuple, shape::Union{SimIdx,Tuple{Vararg{SimIdx}}}) + @_inline_meta + a = promote_indices(as...) + allocate_for(f, a, shape) +end + +promote_indices(a) = a +function promote_indices(a, b, c...) + @_inline_meta + promote_indices(promote_indices(a, b), c...) +end +# overload this to return true for your type, e.g., +# promote_indices(a::OffsetArray, b::OffsetArray) = a +promote_indices(a::AbstractArray, b::AbstractArray) = _promote_indices(indicesbehavior(a), indicesbehavior(b), a, b) +_promote_indices(::IndicesStartAt1, ::IndicesStartAt1, a, b) = a +_promote_indices(::IndicesBehavior, ::IndicesBehavior, a, b) = throw(ArgumentError("types $(typeof(a)) and $(typeof(b)) do not have promote_indices defined")) +promote_indices(a::Number, b::AbstractArray) = b +promote_indices(a::AbstractArray, b::Number) = a + +# Strip off the index-changing container---this assumes that `parent` +# performs such an operation. TODO: since few things in Base need this, it +# would be great to find a way to eliminate this function. +normalize_indices(A) = normalize_indices(indicesbehavior(A), A) +normalize_indices(::IndicesStartAt1, A) = A +normalize_indices(::IndicesBehavior, A) = parent(A) ## from general iterable to any array function copy!(dest::AbstractArray, src) - i = 1 + destiter = eachindex(dest) + state = start(destiter) for x in src + i, state = next(destiter, state) dest[i] = x - i += 1 end return dest end -# if src is not an AbstractArray, moving to the offset might be O(n) -function copy!(dest::AbstractArray, doffs::Integer, src) - doffs < 1 && throw(BoundsError(dest, doffs)) - st = start(src) - i, dmax = doffs, length(dest) - while !done(src, st) - i > dmax && throw(BoundsError(dest, i)) - val, st = next(src, st) - @inbounds dest[i] = val +function copy!(dest::AbstractArray, dstart::Integer, src) + i = Int(dstart) + for x in src + dest[i] = x i += 1 end return dest end # copy from an some iterable object into an AbstractArray -function copy!(dest::AbstractArray, doffs::Integer, src, soffs::Integer) - if (doffs < 1) | (soffs < 1) - doffs < 1 && throw(BoundsError(dest, doffs)) - throw(ArgumentError(string("source start offset (",soffs,") is < 1"))) +function copy!(dest::AbstractArray, dstart::Integer, src, sstart::Integer) + if (sstart < 1) + throw(ArgumentError(string("source start offset (",sstart,") is < 1"))) end st = start(src) - for j = 1:(soffs-1) + for j = 1:(sstart-1) if done(src, st) throw(ArgumentError(string("source has fewer elements than required, ", - "expected at least ",soffs,", got ",j-1))) + "expected at least ",sstart,", got ",j-1))) end _, st = next(src, st) end dn = done(src, st) if dn throw(ArgumentError(string("source has fewer elements than required, ", - "expected at least ",soffs,", got ",soffs-1))) + "expected at least ",sstart,", got ",sstart-1))) end - i, dmax = doffs, length(dest) + i = Int(dstart) while !dn - i > dmax && throw(BoundsError(dest, i)) val, st = next(src, st) - @inbounds dest[i] = val + dest[i] = val i += 1 dn = done(src, st) end @@ -214,24 +456,24 @@ function copy!(dest::AbstractArray, doffs::Integer, src, soffs::Integer) end # this method must be separate from the above since src might not have a length -function copy!(dest::AbstractArray, doffs::Integer, src, soffs::Integer, n::Integer) - n < 0 && throw(BoundsError(dest, n)) +function copy!(dest::AbstractArray, dstart::Integer, src, sstart::Integer, n::Integer) + n < 0 && throw(ArgumentError(string("tried to copy n=", n, " elements, but n should be nonnegative"))) n == 0 && return dest - dmax = doffs + n - 1 - if (dmax > length(dest)) | (doffs < 1) | (soffs < 1) - doffs < 1 && throw(BoundsError(dest, doffs)) - soffs < 1 && throw(ArgumentError(string("source start offset (",soffs,") is < 1"))) - throw(BoundsError(dest, dmax)) + dmax = dstart + n - 1 + inds = linearindices(dest) + if (dstart ∉ inds || dmax ∉ inds) | (sstart < 1) + sstart < 1 && throw(ArgumentError(string("source start offset (",sstart,") is < 1"))) + throw(BoundsError(dest, dstart:dmax)) end st = start(src) - for j = 1:(soffs-1) + for j = 1:(sstart-1) if done(src, st) throw(ArgumentError(string("source has fewer elements than required, ", - "expected at least ",soffs,", got ",j-1))) + "expected at least ",sstart,", got ",j-1))) end _, st = next(src, st) end - i = doffs + i = Int(dstart) while i <= dmax && !done(src, st) val, st = next(src, st) @inbounds dest[i] = val @@ -248,17 +490,17 @@ copy!(dest::AbstractArray, src::AbstractArray) = copy!(linearindexing(dest), dest, linearindexing(src), src) function copy!(::LinearIndexing, dest::AbstractArray, ::LinearIndexing, src::AbstractArray) - n = length(src) - n > length(dest) && throw(BoundsError(dest, n)) - @inbounds for i = 1:n + destinds, srcinds = linearindices(dest), linearindices(src) + isempty(srcinds) || (first(srcinds) ∈ destinds && last(srcinds) ∈ destinds) || throw(BoundsError(dest, srcinds)) + @inbounds for i in srcinds dest[i] = src[i] end return dest end function copy!(::LinearIndexing, dest::AbstractArray, ::LinearSlow, src::AbstractArray) - n = length(src) - n > length(dest) && throw(BoundsError(dest, n)) + destinds, srcinds = linearindices(dest), linearindices(src) + isempty(srcinds) || (first(srcinds) ∈ destinds && last(srcinds) ∈ destinds) || throw(BoundsError(dest, srcinds)) i = 0 @inbounds for a in src dest[i+=1] = a @@ -266,31 +508,34 @@ function copy!(::LinearIndexing, dest::AbstractArray, ::LinearSlow, src::Abstrac return dest end -function copy!(dest::AbstractArray, doffs::Integer, src::AbstractArray) - copy!(dest, doffs, src, 1, length(src)) +function copy!(dest::AbstractArray, dstart::Integer, src::AbstractArray) + copy!(dest, dstart, src, first(linearindices(src)), length(src)) end -function copy!(dest::AbstractArray, doffs::Integer, src::AbstractArray, soffs::Integer) - soffs > length(src) && throw(BoundsError(src, soffs)) - copy!(dest, doffs, src, soffs, length(src)-soffs+1) +function copy!(dest::AbstractArray, dstart::Integer, src::AbstractArray, sstart::Integer) + srcinds = linearindices(src) + sstart ∈ srcinds || throw(BoundsError(src, sstart)) + copy!(dest, dstart, src, sstart, last(srcinds)-sstart+1) end -function copy!(dest::AbstractArray, doffs::Integer, - src::AbstractArray, soffs::Integer, +function copy!(dest::AbstractArray, dstart::Integer, + src::AbstractArray, sstart::Integer, n::Integer) n == 0 && return dest - n < 0 && throw(BoundsError(src, n)) - soffs+n-1 > length(src) && throw(BoundsError(src, soffs+n-1)) - doffs+n-1 > length(dest) && throw(BoundsError(dest, doffs+n-1)) - doffs < 1 && throw(BoundsError(dest, doffs)) - soffs < 1 && throw(BoundsError(src, soffs)) - @inbounds for i = 0:(n-1) #Fixme iter - dest[doffs+i] = src[soffs+i] + n < 0 && throw(ArgumentError(string("tried to copy n=", n, " elements, but n should be nonnegative"))) + destinds, srcinds = linearindices(dest), linearindices(src) + (dstart ∈ destinds && dstart+n-1 ∈ destinds) || throw(BoundsError(dest, dstart:dstart+n-1)) + (sstart ∈ srcinds && sstart+n-1 ∈ srcinds) || throw(BoundsError(src, sstart:sstart+n-1)) + @inbounds for i = 0:(n-1) + dest[dstart+i] = src[sstart+i] end return dest end -copy(a::AbstractArray) = copymutable(a) +function copy(a::AbstractArray) + @_propagate_inbounds_meta + copymutable(a) +end function copy!{R,S}(B::AbstractVecOrMat{R}, ir_dest::Range{Int}, jr_dest::Range{Int}, A::AbstractVecOrMat{S}, ir_src::Range{Int}, jr_src::Range{Int}) @@ -340,7 +585,10 @@ function copy_transpose!{R,S}(B::AbstractVecOrMat{R}, ir_dest::Range{Int}, jr_de return B end -copymutable(a::AbstractArray) = copy!(similar(a), a) +function copymutable(a::AbstractArray) + @_propagate_inbounds_meta + copy!(similar(a), a) +end copymutable(itr) = collect(itr) """ copymutable(a) @@ -360,9 +608,9 @@ zero{T}(x::AbstractArray{T}) = fill!(similar(x), zero(T)) # While the definitions for LinearFast are all simple enough to inline on their # own, LinearSlow's CartesianRange is more complicated and requires explicit # inlining. -start(A::AbstractArray) = (@_inline_meta(); itr = eachindex(A); (itr, start(itr))) -next(A::AbstractArray,i) = (@_inline_meta(); (idx, s) = next(i[1], i[2]); (A[idx], (i[1], s))) -done(A::AbstractArray,i) = done(i[1], i[2]) +start(A::AbstractArray) = (@_inline_meta; itr = eachindex(A); (itr, start(itr))) +next(A::AbstractArray,i) = (@_propagate_inbounds_meta; (idx, s) = next(i[1], i[2]); (A[idx], (i[1], s))) +done(A::AbstractArray,i) = (@_propagate_inbounds_meta; done(i[1], i[2])) # eachindex iterates over all indices. LinearSlow definitions are later. eachindex(A::AbstractArray) = (@_inline_meta(); eachindex(linearindexing(A), A)) @@ -375,7 +623,7 @@ function eachindex(A::AbstractArray, B::AbstractArray...) @_inline_meta eachindex(linearindexing(A,B...), A, B...) end -eachindex(::LinearFast, A::AbstractArray) = 1:length(A) +eachindex(::LinearFast, A::AbstractArray) = linearindices(A) function eachindex(::LinearFast, A::AbstractArray, B::AbstractArray...) @_inline_meta 1:_maxlength(A, B...) @@ -419,7 +667,7 @@ end # a data Ref. they just map the array element type to the pointer type for # convenience in cases that work. pointer{T}(x::AbstractArray{T}) = unsafe_convert(Ptr{T}, x) -pointer{T}(x::AbstractArray{T}, i::Integer) = (@_inline_meta; unsafe_convert(Ptr{T},x) + (i-1)*elsize(x)) +pointer{T}(x::AbstractArray{T}, i::Integer) = (@_inline_meta; unsafe_convert(Ptr{T},x) + (i-first(linearindices(x)))*elsize(x)) ## Approach: @@ -445,21 +693,30 @@ end _getindex(::LinearIndexing, A::AbstractArray, I...) = error("indexing $(typeof(A)) with types $(typeof(I)) is not supported") ## LinearFast Scalar indexing: canonical method is one Int -_getindex(::LinearFast, A::AbstractArray, ::Int) = error("indexing not defined for ", typeof(A)) +_getindex(::LinearFast, A::AbstractVector, ::Int) = error("indexing not defined for ", typeof(A)) +_getindex(::LinearFast, A::AbstractArray, ::Int) = error("indexing not defined for ", typeof(A)) +_getindex{T}(::LinearFast, A::AbstractArray{T,0}) = A[1] _getindex(::LinearFast, A::AbstractArray, i::Real) = (@_propagate_inbounds_meta; getindex(A, to_index(i))) function _getindex{T,N}(::LinearFast, A::AbstractArray{T,N}, I::Vararg{Real,N}) # We must check bounds for sub2ind; so we can then use @inbounds @_inline_meta J = to_indexes(I...) @boundscheck checkbounds(A, J...) - @inbounds r = getindex(A, sub2ind(size(A), J...)) + @inbounds r = getindex(A, sub2ind(A, J...)) + r +end +function _getindex(::LinearFast, A::AbstractVector, I1::Real, I::Real...) + @_inline_meta + J = to_indexes(I1, I...) + @boundscheck checkbounds(A, J...) + @inbounds r = getindex(A, J[1]) r end function _getindex(::LinearFast, A::AbstractArray, I::Real...) # TODO: DEPRECATE FOR #14770 @_inline_meta J = to_indexes(I...) @boundscheck checkbounds(A, J...) - @inbounds r = getindex(A, sub2ind(size(A), J...)) + @inbounds r = getindex(A, sub2ind(A, J...)) r end @@ -471,7 +728,7 @@ function _getindex(::LinearSlow, A::AbstractArray, i::Real) # ind2sub requires all dimensions to be > 0; may as well just check bounds @_inline_meta @boundscheck checkbounds(A, i) - @inbounds r = getindex(A, ind2sub(size(A), to_index(i))...) + @inbounds r = getindex(A, ind2sub(A, to_index(i))...) r end @generated function _getindex{T,AN}(::LinearSlow, A::AbstractArray{T,AN}, I::Real...) # TODO: DEPRECATE FOR #14770 @@ -517,20 +774,28 @@ _setindex!(::LinearIndexing, A::AbstractArray, v, I...) = error("indexing $(type ## LinearFast Scalar indexing _setindex!(::LinearFast, A::AbstractArray, v, ::Int) = error("indexed assignment not defined for ", typeof(A)) +_setindex!{T}(::LinearFast, A::AbstractArray{T,0}, v) = (@_propagate_inbounds_meta; setindex!(A, v, 1)) _setindex!(::LinearFast, A::AbstractArray, v, i::Real) = (@_propagate_inbounds_meta; setindex!(A, v, to_index(i))) function _setindex!{T,N}(::LinearFast, A::AbstractArray{T,N}, v, I::Vararg{Real,N}) # We must check bounds for sub2ind; so we can then use @inbounds @_inline_meta J = to_indexes(I...) @boundscheck checkbounds(A, J...) - @inbounds r = setindex!(A, v, sub2ind(size(A), J...)) + @inbounds r = setindex!(A, v, sub2ind(A, J...)) + r +end +function _setindex!(::LinearFast, A::AbstractArray, v, I1::Real, I::Real...) + @_inline_meta + J = to_indexes(I1, I...) + @boundscheck checkbounds(A, J...) + @inbounds r = setindex!(A, v, J[1]) r end function _setindex!(::LinearFast, A::AbstractArray, v, I::Real...) # TODO: DEPRECATE FOR #14770 @_inline_meta J = to_indexes(I...) @boundscheck checkbounds(A, J...) - @inbounds r = setindex!(A, v, sub2ind(size(A), J...)) + @inbounds r = setindex!(A, v, sub2ind(A, J...)) r end @@ -541,7 +806,7 @@ function _setindex!(::LinearSlow, A::AbstractArray, v, i::Real) # ind2sub requires all dimensions to be > 0; may as well just check bounds @_inline_meta @boundscheck checkbounds(A, i) - @inbounds r = setindex!(A, v, ind2sub(size(A), to_index(i))...) + @inbounds r = setindex!(A, v, ind2sub(A, to_index(i))...) r end @generated function _setindex!{T,AN}(::LinearSlow, A::AbstractArray{T,AN}, v, I::Real...) # TODO: DEPRECATE FOR #14770 @@ -576,11 +841,21 @@ end typealias RangeVecIntList{A<:AbstractVector{Int}} Union{Tuple{Vararg{Union{Range, AbstractVector{Int}}}}, AbstractVector{UnitRange{Int}}, AbstractVector{Range{Int}}, AbstractVector{A}} -get(A::AbstractArray, i::Integer, default) = checkbounds(Bool, length(A), i) ? A[i] : default +get(A::AbstractArray, i::Integer, default) = checkbounds(Bool, A, i) ? A[i] : default get(A::AbstractArray, I::Tuple{}, default) = similar(A, typeof(default), 0) -get(A::AbstractArray, I::Dims, default) = checkbounds(Bool, size(A), I...) ? A[I...] : default +get(A::AbstractArray, I::Dims, default) = checkbounds(Bool, A, I...) ? A[I...] : default +function get!{T}(X::AbstractVector{T}, A::AbstractVector, I::Union{Range, AbstractVector{Int}}, default::T) + # 1d is not linear indexing + ind = findin(I, indices1(A)) + X[ind] = A[I[ind]] + Xind = indices1(X) + X[first(Xind):first(ind)-1] = default + X[last(ind)+1:last(Xind)] = default + X +end function get!{T}(X::AbstractArray{T}, A::AbstractArray, I::Union{Range, AbstractVector{Int}}, default::T) + # Linear indexing ind = findin(I, 1:length(A)) X[ind] = A[I[ind]] X[1:first(ind)-1] = default @@ -588,8 +863,9 @@ function get!{T}(X::AbstractArray{T}, A::AbstractArray, I::Union{Range, Abstract X end -get(A::AbstractArray, I::Range, default) = get!(similar(A, typeof(default), length(I)), A, I, default) +get(A::AbstractArray, I::Range, default) = get!(similar(A, typeof(default), index_shape(A, I)), A, I, default) +# TODO: DEPRECATE FOR #14770 (just the partial linear indexing part) function get!{T}(X::AbstractArray{T}, A::AbstractArray, I::RangeVecIntList, default::T) fill!(X, default) dst, src = indcopy(size(A), I) @@ -597,7 +873,7 @@ function get!{T}(X::AbstractArray{T}, A::AbstractArray, I::RangeVecIntList, defa X end -get(A::AbstractArray, I::RangeVecIntList, default) = get!(similar(A, typeof(default), map(length, I)...), A, I, default) +get(A::AbstractArray, I::RangeVecIntList, default) = get!(similar(A, typeof(default), index_shape(A, I...)), A, I, default) ## structured matrix methods ## replace_in_print_matrix(A::AbstractMatrix,i::Integer,j::Integer,s::AbstractString) = s @@ -922,7 +1198,7 @@ end function isequal(A::AbstractArray, B::AbstractArray) if A === B return true end - if size(A) != size(B) + if indices(A) != indices(B) return false end if isa(A,Range) != isa(B,Range) @@ -945,7 +1221,7 @@ function lexcmp(A::AbstractArray, B::AbstractArray) end function (==)(A::AbstractArray, B::AbstractArray) - if size(A) != size(B) + if indices(A) != indices(B) return false end if isa(A,Range) != isa(B,Range) @@ -959,69 +1235,121 @@ function (==)(A::AbstractArray, B::AbstractArray) return true end -sub2ind(dims::Tuple{Vararg{Integer}}) = 1 -sub2ind(dims::Tuple{Vararg{Integer}}, I::Integer...) = _sub2ind(dims,I) -@generated function _sub2ind{N,M}(dims::NTuple{N,Integer}, I::NTuple{M,Integer}) - meta = Expr(:meta,:inline) - ex = :(I[$M] - 1) - for i = M-1:-1:1 - if i > N - ex = :(I[$i] - 1 + $ex) - else - ex = :(I[$i] - 1 + dims[$i]*$ex) - end - end - Expr(:block, meta,:($ex + 1)) +# sub2ind and ind2sub +# fallbacks +function sub2ind(A::AbstractArray, I...) + @_inline_meta + sub2ind(indicesbehavior(A), A, I...) +end +function sub2ind(::IndicesStartAt1, A::AbstractArray, I...) + @_inline_meta + sub2ind(size(A), I...) +end +function sub2ind(::IndicesBehavior, A::AbstractArray, I...) + @_inline_meta + sub2ind(indices(A), I...) +end +function ind2sub(A::AbstractArray, ind) + @_inline_meta + ind2sub(indicesbehavior(A), A, ind) +end +function ind2sub(::IndicesStartAt1, A::AbstractArray, ind) + @_inline_meta + ind2sub(size(A), ind) +end +function ind2sub(::IndicesBehavior, A::AbstractArray, ind) + @_inline_meta + ind2sub(indices(A), ind) end -@generated function ind2sub{N}(dims::NTuple{N,Integer}, ind::Integer) - meta = Expr(:meta,:inline) - N==0 && return :($meta; ind==1 ? () : throw(BoundsError())) - exprs = Expr[:(ind = ind-1)] - for i = 1:N-1 - push!(exprs,:(ind2 = div(ind,dims[$i]))) - push!(exprs,Expr(:(=),Symbol(:s,i),:(ind-dims[$i]*ind2+1))) - push!(exprs,:(ind=ind2)) - end - push!(exprs,Expr(:(=),Symbol(:s,N),:(ind+1))) - Expr(:block,meta,exprs...,Expr(:tuple,[Symbol(:s,i) for i=1:N]...)) +sub2ind(::Tuple{}) = 1 +sub2ind(::DimsInteger) = 1 +sub2ind(::Indices) = 1 +sub2ind(::Tuple{}, I::Integer...) = (@_inline_meta; _sub2ind((), 1, 1, I...)) +sub2ind(dims::DimsInteger, I::Integer...) = (@_inline_meta; _sub2ind(dims, 1, 1, I...)) +sub2ind(inds::Indices, I::Integer...) = (@_inline_meta; _sub2ind(inds, 1, 1, I...)) +# In 1d, there's a question of whether we're doing cartesian indexing or linear indexing. Support only the former. +sub2ind(inds::Indices{1}, I::Integer...) = throw(ArgumentError("Linear indexing is not defined for one-dimensional arrays")) + +_sub2ind(::Any, L, ind) = ind +function _sub2ind(::Tuple{}, L, ind, i::Integer, I::Integer...) + @_inline_meta + _sub2ind((), L, ind+(i-1)*L, I...) +end +function _sub2ind(inds, L, ind, i::Integer, I::Integer...) + @_inline_meta + r1 = inds[1] + _sub2ind(tail(inds), nextL(L, r1), ind+offsetin(i, r1)*L, I...) end -ind2sub(a::AbstractArray, ind::Integer) = ind2sub(size(a), ind) -sub2ind(a::AbstractArray, I::Integer...) = sub2ind(size(a), I...) +nextL(L, l::Integer) = L*l +nextL(L, r::UnitRange) = L*unsafe_length(r) +offsetin(i, l::Integer) = i-1 +offsetin(i, r::UnitRange) = i-first(r) +unsafe_length(r::UnitRange) = r.stop-r.start+1 -function sub2ind{T<:Integer}(dims::Tuple{Vararg{Integer}}, I::AbstractVector{T}...) - N = length(dims) - M = length(I[1]) - indices = Array{T}(length(I[1])) - copy!(indices,I[1]) +ind2sub(::Tuple{}, ind::Integer) = (@_inline_meta; ind == 1 ? () : throw(BoundsError())) +ind2sub(dims::DimsInteger, ind::Integer) = (@_inline_meta; _ind2sub((), dims, ind-1)) +ind2sub(inds::Indices, ind::Integer) = (@_inline_meta; _ind2sub((), inds, ind-1)) +ind2sub(inds::Indices{1}, ind::Integer) = throw(ArgumentError("Linear indexing is not defined for one-dimensional arrays")) - s = dims[1] - for j=2:length(I) - Ij = I[j] - for (i, k) in zip(eachindex(indices), eachindex(Ij)) - indices[i] += s*(Ij[k]-1) - end - s *= (j <= N ? dims[j] : 1) +_ind2sub(::Tuple{}, ::Tuple{}, ind) = (ind+1,) +function _ind2sub(out, indslast::NTuple{1}, ind) + @_inline_meta + (out..., _lookup(ind, indslast[1])) +end +function _ind2sub(out, inds, ind) + @_inline_meta + r1 = inds[1] + indnext, f, l = _div(ind, r1) + _ind2sub((out..., ind-l*indnext+f), tail(inds), indnext) +end + +_lookup(ind, d::Integer) = ind+1 +_lookup(ind, r::UnitRange) = ind+first(r) +_div(ind, d::Integer) = div(ind, d), 1, d +_div(ind, r::UnitRange) = (d = unsafe_length(r); (div(ind, d), first(r), d)) + +smart_ind2sub(shape::NTuple{1}, ind) = (ind,) +smart_ind2sub(shape, ind) = ind2sub(shape, ind) +smart_sub2ind(shape::NTuple{1}, i) = (i,) +smart_sub2ind(shape, I...) = (@_inline_meta; sub2ind(shape, I...)) + +# Vectorized forms +function sub2ind{N,T<:Integer}(inds::Union{Dims{N},Indices{N}}, I::AbstractVector{T}...) + I1 = I[1] + Iinds = indices1(I1) + for j = 2:length(I) + indices1(I[j]) == Iinds || throw(DimensionMismatch("indices of I[1] ($(Iinds)) does not match indices of I[$j] ($(indices1(I[j])))")) + end + Iout = similar(I1) + _sub2ind!(Iout, inds, Iinds, I) + Iout +end + +function _sub2ind!(Iout, inds, Iinds, I) + @_noinline_meta + for i in Iinds + # Iout[i] = sub2ind(inds, map(Ij->Ij[i], I)...) + Iout[i] = sub2ind_vec(inds, i, I) end - return indices + Iout end -function ind2sub{N,T<:Integer}(dims::NTuple{N,Integer}, ind::AbstractVector{T}) +sub2ind_vec(inds, i, I) = (@_inline_meta; _sub2ind_vec(inds, (), i, I...)) +_sub2ind_vec(inds, out, i, I1, I...) = (@_inline_meta; _sub2ind_vec(inds, (out..., I1[i]), i, I...)) +_sub2ind_vec(inds, out, i) = (@_inline_meta; sub2ind(inds, out...)) + +function ind2sub{N,T<:Integer}(inds::Union{Dims{N},Indices{N}}, ind::AbstractVector{T}) M = length(ind) - t = NTuple{N,Vector{T}}(ntuple(n->Array{T}(M),N)) - copy!(t[1],ind) - for j = 1:N-1 - d = dims[j] - tj = t[j] - tj2 = t[j+1] - for i = 1:M - ind2 = div(tj[i]-1, d) - tj[i] -= d*ind2 - tj2[i] = ind2+1 + t = ntuple(n->similar(ind),Val{N}) + for (i,idx) in enumerate(ind) # FIXME: change to eachindexvalue + sub = ind2sub(inds, idx) + for j = 1:N + t[j][i] = sub[j] end end - return t + t end function ind2sub!{T<:Integer}(sub::Array{T}, dims::Tuple{Vararg{T}}, ind::T) @@ -1060,21 +1388,19 @@ function mapslices(f, A::AbstractArray, dims::AbstractVector) return map(f,A) end - dimsA = [size(A)...] + dimsA = [shape(A)...] ndimsA = ndims(A) alldims = [1:ndimsA;] otherdims = setdiff(alldims, dims) - idx = Array{Any}(ndimsA) - fill!(idx, 1) - Asliceshape = tuple(dimsA[dims]...) + idx = Any[first(ind) for ind in indices(A)] itershape = tuple(dimsA[otherdims]...) for d in dims - idx[d] = 1:size(A,d) + idx[d] = Colon() end - r1 = f(reshape(A[idx...], Asliceshape)) + r1 = f(slice(A, idx...)) # determine result size and allocate Rsize = copy(dimsA) @@ -1082,27 +1408,31 @@ function mapslices(f, A::AbstractArray, dims::AbstractVector) if !isa(r1, AbstractArray) || ndims(r1) == 0 r1 = [r1] end - Rsize[dims] = [size(r1)...; ones(Int,max(0,length(dims)-ndims(r1)))] - R = similar(r1, tuple(Rsize...)) + nextra = max(0,length(dims)-ndims(r1)) + if eltype(Rsize) == Int + Rsize[dims] = [size(r1)..., ntuple(d->1, nextra)...] + else + Rsize[dims] = [indices(r1)..., ntuple(d->1:1, nextra)...] + end + R = similar(r1, tuple(Rsize...,)) - ridx = Array{Any}(ndims(R)) - fill!(ridx, 1) + ridx = Any[map(first, indices(R))...] for d in dims - ridx[d] = 1:size(R,d) + ridx[d] = indices(R,d) end R[ridx...] = r1 - first = true + isfirst = true nidx = length(otherdims) for I in CartesianRange(itershape) - if first - first = false + if isfirst + isfirst = false else for i in 1:nidx idx[otherdims[i]] = ridx[otherdims[i]] = I.I[i] end - R[ridx...] = f(reshape(A[idx...], Asliceshape)) + R[ridx...] = f(slice(A, idx...)) end end @@ -1147,8 +1477,7 @@ ith_all(i, ::Tuple{}) = () ith_all(i, as) = (as[1][i], ith_all(i, tail(as))...) function map_n!{F}(f::F, dest::AbstractArray, As) - n = length(As[1]) - for i = 1:n #Fixme iter, one might make a @generated function here + for i = linearindices(As[1]) dest[i] = f(ith_all(i, As)...) end return dest diff --git a/base/abstractarraymath.jl b/base/abstractarraymath.jl index 53b70ea611180..d4b2e48e885d8 100644 --- a/base/abstractarraymath.jl +++ b/base/abstractarraymath.jl @@ -61,7 +61,7 @@ imag{T<:Real}(x::AbstractArray{T}) = zero(x) # index A[:,:,...,i,:,:,...] where "i" is in dimension "d" # TODO: more optimized special cases slicedim(A::AbstractArray, d::Integer, i) = - A[[ n==d ? i : (1:size(A,n)) for n in 1:ndims(A) ]...] + A[[ n==d ? i : (indices(A,n)) for n in 1:ndims(A) ]...] function flipdim(A::AbstractVector, d::Integer) d > 0 || throw(ArgumentError("dimension to flip must be positive")) @@ -71,8 +71,7 @@ end function flipdim(A::AbstractArray, d::Integer) nd = ndims(A) - sd = d > nd ? 1 : size(A, d) - if sd == 1 || isempty(A) + if d > nd || isempty(A) return copy(A) end B = similar(A) @@ -80,16 +79,18 @@ function flipdim(A::AbstractArray, d::Integer) for i = 1:nd nnd += Int(size(A,i)==1 || i==d) end + inds = indices(A, d) + sd = first(inds)+last(inds) if nnd==nd # flip along the only non-singleton dimension - for i = 1:sd - B[i] = A[sd+1-i] + for i in inds + B[i] = A[sd-i] end return B end - alli = [ 1:size(B,n) for n in 1:nd ] - for i = 1:sd - B[[ n==d ? sd+1-i : alli[n] for n in 1:nd ]...] = slicedim(A, d, i) + alli = [ indices(B,n) for n in 1:nd ] + for i in inds + B[[ n==d ? sd-i : alli[n] for n in 1:nd ]...] = slicedim(A, d, i) end return B end @@ -107,13 +108,14 @@ end # Uses K-B-N summation function cumsum_kbn{T<:AbstractFloat}(v::AbstractVector{T}) - n = length(v) - r = similar(v, n) - if n == 0; return r; end + r = similar(v) + if isempty(v); return r; end - s = r[1] = v[1] + inds = indices(v, 1) + i1 = first(inds) + s = r[i1] = v[i1] c = zero(T) - for i=2:n #Fixme iter + for i=i1+1:last(inds) vi = v[i] t = s + vi if abs(s) >= abs(vi) diff --git a/base/array.jl b/base/array.jl index bdd145a21b816..0dea6e44e9bbe 100644 --- a/base/array.jl +++ b/base/array.jl @@ -56,7 +56,8 @@ end function copy!{T}(dest::Array{T}, doffs::Integer, src::Array{T}, soffs::Integer, n::Integer) n == 0 && return dest - if n < 0 || soffs < 1 || doffs < 1 || soffs+n-1 > length(src) || doffs+n-1 > length(dest) + n > 0 || throw(ArgumentError(string("tried to copy n=", n, " elements, but n should be nonnegative"))) + if soffs < 1 || doffs < 1 || soffs+n-1 > length(src) || doffs+n-1 > length(dest) throw(BoundsError()) end unsafe_copy!(dest, doffs, src, soffs, n) @@ -299,8 +300,8 @@ end ## Iteration ## start(A::Array) = 1 -next(a::Array,i) = (a[i],i+1) -done(a::Array,i) = i == length(a)+1 +next(a::Array,i) = (@_propagate_inbounds_meta; (a[i],i+1)) +done(a::Array,i) = (@_inline_meta; i == length(a)+1) ## Indexing: getindex ## @@ -917,6 +918,7 @@ function findin(a, b) end # Copying subregions +# TODO: DEPRECATE FOR #14770 function indcopy(sz::Dims, I::Vector) n = length(I) s = sz[n] diff --git a/base/arraymath.jl b/base/arraymath.jl index 235bffe676296..cb3cd4625e26e 100644 --- a/base/arraymath.jl +++ b/base/arraymath.jl @@ -59,21 +59,21 @@ for f in (:+, :-, :div, :mod, :&, :|, :$) return F end function ($f){S,T}(A::AbstractArray{S}, B::Range{T}) - F = similar(A, promote_op($f,S,T), promote_shape(size(A),size(B))) + F = similar(A, promote_op($f,S,T), promote_shape(A,B)) for (iF, iA, iB) in zip(eachindex(F), eachindex(A), eachindex(B)) @inbounds F[iF] = ($f)(A[iA], B[iB]) end return F end function ($f){S,T}(A::Range{S}, B::AbstractArray{T}) - F = similar(B, promote_op($f,S,T), promote_shape(size(A),size(B))) + F = similar(B, promote_op($f,S,T), promote_shape(A,B)) for (iF, iA, iB) in zip(eachindex(F), eachindex(A), eachindex(B)) @inbounds F[iF] = ($f)(A[iA], B[iB]) end return F end function ($f){S,T}(A::AbstractArray{S}, B::AbstractArray{T}) - F = similar(A, promote_op($f,S,T), promote_shape(size(A),size(B))) + F = similar(A, promote_op($f,S,T), promote_shape(A,B)) for (iF, iA, iB) in zip(eachindex(F), eachindex(A), eachindex(B)) @inbounds F[iF] = ($f)(A[iA], B[iB]) end @@ -211,26 +211,29 @@ function flipdim{T}(A::Array{T}, d::Integer) end function rotl90(A::AbstractMatrix) - m,n = size(A) - B = similar(A,(n,m)) - for i=1:m, j=1:n #Fixme iter - B[n-j+1,i] = A[i,j] + B = similar_transpose(A) + ind2 = indices(A,2) + n = first(ind2)+last(ind2) + for i=indices(A,1), j=ind2 + B[n-j,i] = A[i,j] end return B end function rotr90(A::AbstractMatrix) - m,n = size(A) - B = similar(A,(n,m)) - for i=1:m, j=1:n #Fixme iter - B[j,m-i+1] = A[i,j] + B = similar_transpose(A) + ind1 = indices(A,1) + m = first(ind1)+last(ind1) + for i=ind1, j=indices(A,2) + B[j,m-i] = A[i,j] end return B end function rot180(A::AbstractMatrix) - m,n = size(A) B = similar(A) - for i=1:m, j=1:n #Fixme iter - B[m-i+1,n-j+1] = A[i,j] + ind1, ind2 = indices(A,1), indices(A,2) + m, n = first(ind1)+last(ind1), first(ind2)+last(ind2) + for j=ind2, i=ind1 + B[m-i,n-j] = A[i,j] end return B end @@ -243,96 +246,65 @@ end rotr90(A::AbstractMatrix, k::Integer) = rotl90(A,-k) rot180(A::AbstractMatrix, k::Integer) = mod(k, 2) == 1 ? rot180(A) : copy(A) +similar_transpose(A::AbstractMatrix) = similar_transpose(indicesbehavior(A), A) +similar_transpose(::IndicesStartAt1, A::AbstractMatrix) = similar(A, (size(A,2), size(A,1))) +similar_transpose(::IndicesBehavior, A::AbstractMatrix) = similar(A, (indices(A,2), indices(A,1))) ## Transpose ## -const transposebaselength=64 -function transpose!(B::AbstractMatrix,A::AbstractMatrix) - m, n = size(A) - size(B,1) == n && size(B,2) == m || throw(DimensionMismatch("transpose")) - - if m*n<=4*transposebaselength - @inbounds begin - for j = 1:n #Fixme iter - for i = 1:m #Fixme iter - B[j,i] = transpose(A[i,j]) - end - end - end - else - transposeblock!(B,A,m,n,0,0) - end - return B -end +transpose!(B::AbstractMatrix, A::AbstractMatrix) = transpose_f!(transpose, B, A) +ctranspose!(B::AbstractMatrix, A::AbstractMatrix) = transpose_f!(ctranspose, B, A) function transpose!(B::AbstractVector, A::AbstractMatrix) - length(B) == length(A) && size(A,1) == 1 || throw(DimensionMismatch("transpose")) + indices(B,1) == indices(A,2) && indices(A,1) == 1:1 || throw(DimensionMismatch("transpose")) copy!(B, A) end function transpose!(B::AbstractMatrix, A::AbstractVector) - length(B) == length(A) && size(B,1) == 1 || throw(DimensionMismatch("transpose")) + indices(B,2) == indices(A,1) && indices(B,1) == 1:1 || throw(DimensionMismatch("transpose")) copy!(B, A) end -function transposeblock!(B::AbstractMatrix,A::AbstractMatrix,m::Int,n::Int,offseti::Int,offsetj::Int) - if m*n<=transposebaselength - @inbounds begin - for j = offsetj+(1:n) #Fixme iter - for i = offseti+(1:m) #Fixme iter - B[j,i] = transpose(A[i,j]) - end - end - end - elseif m>n - newm=m>>1 - transposeblock!(B,A,newm,n,offseti,offsetj) - transposeblock!(B,A,m-newm,n,offseti+newm,offsetj) - else - newn=n>>1 - transposeblock!(B,A,m,newn,offseti,offsetj) - transposeblock!(B,A,m,n-newn,offseti,offsetj+newn) - end - return B +function ctranspose!(B::AbstractVector, A::AbstractMatrix) + indices(B,1) == indices(A,2) && indices(A,1) == 1:1 || throw(DimensionMismatch("transpose")) + ccopy!(B, A) end -function ctranspose!(B::AbstractMatrix,A::AbstractMatrix) - m, n = size(A) - size(B,1) == n && size(B,2) == m || throw(DimensionMismatch("transpose")) +function ctranspose!(B::AbstractMatrix, A::AbstractVector) + indices(B,2) == indices(A,1) && indices(B,1) == 1:1 || throw(DimensionMismatch("transpose")) + ccopy!(B, A) +end + +const transposebaselength=64 +function transpose_f!(f,B::AbstractMatrix,A::AbstractMatrix) + indices(B,1) == indices(A,2) && indices(B,2) == indices(A,1) || throw(DimensionMismatch(string(f))) + m, n = size(A) if m*n<=4*transposebaselength @inbounds begin - for j = 1:n #Fixme iter - for i = 1:m #Fixme iter - B[j,i] = ctranspose(A[i,j]) + for j = indices(A,2) + for i = indices(A,1) + B[j,i] = f(A[i,j]) end end end else - ctransposeblock!(B,A,m,n,0,0) + transposeblock!(f,B,A,m,n,first(indices(A,1))-1,first(indices(A,2))-1) end return B end -function ctranspose!(B::AbstractVector, A::AbstractMatrix) - length(B) == length(A) && size(A,1) == 1 || throw(DimensionMismatch("transpose")) - ccopy!(B, A) -end -function ctranspose!(B::AbstractMatrix, A::AbstractVector) - length(B) == length(A) && size(B,1) == 1 || throw(DimensionMismatch("transpose")) - ccopy!(B, A) -end -function ctransposeblock!(B::AbstractMatrix,A::AbstractMatrix,m::Int,n::Int,offseti::Int,offsetj::Int) +function transposeblock!(f,B::AbstractMatrix,A::AbstractMatrix,m::Int,n::Int,offseti::Int,offsetj::Int) if m*n<=transposebaselength @inbounds begin - for j = offsetj+(1:n) #Fixme iter - for i = offseti+(1:m) #Fixme iter - B[j,i] = ctranspose(A[i,j]) + for j = offsetj+(1:n) + for i = offseti+(1:m) + B[j,i] = f(A[i,j]) end end end elseif m>n newm=m>>1 - ctransposeblock!(B,A,newm,n,offseti,offsetj) - ctransposeblock!(B,A,m-newm,n,offseti+newm,offsetj) + transposeblock!(f,B,A,newm,n,offseti,offsetj) + transposeblock!(f,B,A,m-newm,n,offseti+newm,offsetj) else newn=n>>1 - ctransposeblock!(B,A,m,newn,offseti,offsetj) - ctransposeblock!(B,A,m,n-newn,offseti,offsetj+newn) + transposeblock!(f,B,A,m,newn,offseti,offsetj) + transposeblock!(f,B,A,m,n-newn,offseti,offsetj+newn) end return B end @@ -343,11 +315,11 @@ function ccopy!(B, A) end function transpose(A::AbstractMatrix) - B = similar(A, size(A, 2), size(A, 1)) + B = similar_transpose(A) transpose!(B, A) end function ctranspose(A::AbstractMatrix) - B = similar(A, size(A, 2), size(A, 1)) + B = similar_transpose(A) ctranspose!(B, A) end ctranspose{T<:Real}(A::AbstractVecOrMat{T}) = transpose(A) @@ -366,7 +338,7 @@ for (f, f!, fp, op) = ((:cumsum, :cumsum!, :cumsum_pairwise!, :+), if n < 128 @inbounds s_ = v[i1] @inbounds c[i1] = ($op)(s, s_) - for i = i1+1:i1+n-1 #Fixme iter + for i = i1+1:i1+n-1 @inbounds s_ = $(op)(s_, v[i]) @inbounds c[i] = $(op)(s, s_) end @@ -381,7 +353,7 @@ for (f, f!, fp, op) = ((:cumsum, :cumsum!, :cumsum_pairwise!, :+), @eval function ($f!)(result::AbstractVector, v::AbstractVector) n = length(v) if n == 0; return result; end - ($fp)(v, result, $(op==:+ ? :(zero(v[1])) : :(one(v[1]))), 1, n) + ($fp)(v, result, $(op==:+ ? :(zero(first(v))) : :(one(first(v)))), first(indices(v,1)), n) return result end diff --git a/base/broadcast.jl b/base/broadcast.jl index 8b96d8fd4e499..0c0165fac300e 100644 --- a/base/broadcast.jl +++ b/base/broadcast.jl @@ -2,163 +2,129 @@ module Broadcast -using ..Cartesian -using Base: promote_op, promote_eltype, promote_eltype_op, @get!, _msk_end, unsafe_bitgetindex +using Base.Cartesian +using Base: promote_op, promote_eltype, promote_eltype_op, @get!, _msk_end, unsafe_bitgetindex, shape, linearindices, allocate_for, tail, dimlength import Base: .+, .-, .*, ./, .\, .//, .==, .<, .!=, .<=, .÷, .%, .<<, .>>, .^ export broadcast, broadcast!, bitbroadcast export broadcast_getindex, broadcast_setindex! ## Broadcasting utilities ## -droparg1(a, args...) = args - -longer_tuple(x::Tuple{}, retx::Tuple, y::Tuple{}, rety::Tuple) = retx -longer_tuple(x::Tuple{}, retx::Tuple, y::Tuple, rety::Tuple) = rety -longer_tuple(x::Tuple, retx::Tuple, y::Tuple{}, rety::Tuple) = retx -longer_tuple(x::Tuple, retx::Tuple, y::Tuple, rety::Tuple) = - longer_tuple(droparg1(x...), retx, droparg1(y...), rety) -longer_tuple(x::Tuple, y::Tuple) = longer_tuple(x, x, y, y) - -longer_size(x::Union{AbstractArray,Number}) = size(x) -longer_size(x::Union{AbstractArray,Number}, y::Union{AbstractArray,Number}...) = - longer_tuple(size(x), longer_size(y...)) - -# Calculate the broadcast shape of the arguments, or error if incompatible +## Calculate the broadcast shape of the arguments, or error if incompatible +# array inputs broadcast_shape() = () -function broadcast_shape(As::Union{AbstractArray,Number}...) - sz = longer_size(As...) - nd = length(sz) - bshape = ones(Int, nd) - for A in As - for d = 1:ndims(A) - n = size(A, d) - if n != 1 - if bshape[d] == 1 - bshape[d] = n - elseif bshape[d] != n - throw(DimensionMismatch("arrays could not be broadcast to a common size")) - end - end - end - end - return tuple(bshape...)::typeof(sz) +broadcast_shape(A) = shape(A) +@inline broadcast_shape(A, B...) = broadcast_shape((), shape(A), map(shape, B)...) +# shape inputs +broadcast_shape(shape::Tuple) = shape +@inline broadcast_shape(shape::Tuple, shape1::Tuple, shapes::Tuple...) = broadcast_shape(_bcs((), shape, shape1), shapes...) +# _bcs consolidates two shapes into a single output shape +_bcs(out, ::Tuple{}, ::Tuple{}) = out +@inline _bcs(out, ::Tuple{}, newshape) = _bcs((out..., newshape[1]), (), tail(newshape)) +@inline _bcs(out, shape, ::Tuple{}) = _bcs((out..., shape[1]), tail(shape), ()) +@inline function _bcs(out, shape, newshape) + newout = _bcs1(shape[1], newshape[1]) + _bcs((out..., newout), tail(shape), tail(newshape)) end - -# Check that all arguments are broadcast compatible with shape -function check_broadcast_shape(shape::Dims, As::Union{AbstractArray,Number}...) - for A in As - if ndims(A) > length(shape) - throw(DimensionMismatch("cannot broadcast array to have fewer dimensions")) - end - for k in 1:ndims(A) - n, nA = shape[k], size(A, k) - if n != nA != 1 - throw(DimensionMismatch("array could not be broadcast to match destination")) - end - end - end +# _bcs1 handles the logic for a single dimension +_bcs1(a::Integer, b::Integer) = a == 1 ? b : (b == 1 ? a : (a == b ? a : throw(DimensionMismatch("arrays could not be broadcast to a common size")))) +_bcs1(a::Integer, b) = a == 1 ? b : (first(b) == 1 && last(b) == a ? b : throw(DimensionMismatch("arrays could not be broadcast to a common size"))) +_bcs1(a, b::Integer) = _bcs1(b, a) +_bcs1(a, b) = _bcsm(b, a) ? b : (_bcsm(a, b) ? a : throw(DimensionMismatch("arrays could not be broadcast to a common size"))) +# _bcsm tests whether the second index is consistent with the first +_bcsm(a, b) = a == b || (dimlength(b) == 1 && first(b) == first(a)) +_bcsm(a, b::Number) = b == 1 +_bcsm(a::Number, b::Number) = a == b || b == 1 + +## Check that all arguments are broadcast compatible with shape +## Check that all arguments are broadcast compatible with shape +# comparing one input against a shape +check_broadcast_shape(shp) = nothing +check_broadcast_shape(shp, A) = check_broadcast_shape(shp, shape(A)) +check_broadcast_shape(::Tuple{}, ::Tuple{}) = nothing +check_broadcast_shape(shp, ::Tuple{}) = nothing +check_broadcast_shape(::Tuple{}, Ashp::Tuple) = throw(DimensionMismatch("cannot broadcast array to have fewer dimensions")) +function check_broadcast_shape(shp, Ashp::Tuple) + _bcsm(shp[1], Ashp[1]) || throw(DimensionMismatch("array could not be broadcast to match destination")) + check_broadcast_shape(tail(shp), tail(Ashp)) end - -## Broadcasting core -# Generate the body for a broadcasting function f_broadcast!(B, A1, A2, ..., A$narrays), -# using function f, output B, and inputs As... -# B must have already been set to the appropriate size. - -# version using cartesian indexing -function gen_broadcast_body_cartesian(nd::Int, narrays::Int, f) - F = Expr(:quote, f) - quote - @assert ndims(B) == $nd - @ncall $narrays check_broadcast_shape size(B) k->A_k - @nloops($nd, i, B, - d->(@nexprs $narrays k->(j_d_k = size(A_k, d) == 1 ? 1 : i_d)), # pre - begin # body - @nexprs $narrays k->(@inbounds v_k = @nref $nd A_k d->j_d_k) - @inbounds (@nref $nd B i) = (@ncall $narrays $F v) - end) - end -end - -# version using start/next for iterating over the arguments -function gen_broadcast_body_iter(nd::Int, narrays::Int, f) - F = Expr(:quote, f) - quote - @assert ndims(B) == $nd - @ncall $narrays check_broadcast_shape size(B) k->A_k - @nexprs 1 d->(@nexprs $narrays k->(state_k_0 = state_k_{$nd} = start(A_k))) - @nexprs $nd d->(@nexprs $narrays k->(skip_k_d = size(A_k, d) == 1)) - @nloops($nd, i, B, - d->(@nexprs $narrays k->(state_k_{d-1} = state_k_d)), # pre - d->(@nexprs $narrays k->(skip_k_d || (state_k_d = state_k_0))), # post - begin # body - @nexprs $narrays k->(@inbounds (v_k, state_k_0) = next(A_k, state_k_0)) - @inbounds (@nref $nd B i) = (@ncall $narrays $F v) - end) - end +# comparing many inputs +@inline function check_broadcast_shape(shp, A, As...) + check_broadcast_shape(shp, A) + check_broadcast_shape(shp, As...) end -## Broadcasting cores specialized for returning a BitArray - +## Indexing manipulations +# newindex(I, rule) replaces a CartesianIndex with something that is +# appropriate for a particular array/scalar. `rule` is a tuple that +# describes the manipulations that should be made. +@inline newindex(I::CartesianIndex, ::Tuple{}) = 1 # for scalars +@inline newindex(I::CartesianIndex, indexmap) = CartesianIndex(_newindex((), I.I, indexmap...)) +@inline _newindex(out, I) = out # can truncate if indexmap is shorter than I +@inline _newindex(out, I, keep::Bool, indexmap...) = _newindex((out..., ifelse(keep, I[1], 1)), tail(I), indexmap...) + +newindexer(sz, x::Number) = () +@inline newindexer(sz, A) = _newindexer(sz, size(A)) +@inline _newindexer(sz, szA::Tuple{}) = () +@inline _newindexer(sz, szA) = (sz[1] == szA[1], _newindexer(tail(sz), tail(szA))...) + +# map(x->newindexer(sz, x), As), but see #15276 +map_newindexer(sz, ::Tuple{}) = () +@inline map_newindexer(sz, As) = (newindexer(sz, As[1]), map_newindexer(sz, tail(As))...) + +# For output BitArrays const bitcache_chunks = 64 # this can be changed const bitcache_size = 64 * bitcache_chunks # do not change this dumpbitcache(Bc::Vector{UInt64}, bind::Int, C::Vector{Bool}) = Base.copy_to_bitarray_chunks!(Bc, ((bind - 1) << 6) + 1, C, 1, min(bitcache_size, (length(Bc)-bind+1) << 6)) -# using cartesian indexing -function gen_broadcast_body_cartesian_tobitarray(nd::Int, narrays::Int, f) - F = Expr(:quote, f) +## Broadcasting core +# nargs encodes the number of As arguments (which matches the number +# of indexmaps). The first two type parameters are to ensure specialization. +@generated function _broadcast!{M,AT,nargs}(f, B::AbstractArray, indexmaps::M, As::AT, ::Type{Val{nargs}}) quote - @assert ndims(B) == $nd - @ncall $narrays check_broadcast_shape size(B) k->A_k - C = Array{Bool}(bitcache_size) - Bc = B.chunks - ind = 1 - cind = 1 - @nloops($nd, i, B, - d->(@nexprs $narrays k->(j_d_k = size(A_k, d) == 1 ? 1 : i_d)), # pre - begin # body - @nexprs $narrays k->(@inbounds v_k = @nref $nd A_k d->j_d_k) - @inbounds C[ind] = (@ncall $narrays $F v) - ind += 1 - if ind > bitcache_size - dumpbitcache(Bc, cind, C) - cind += bitcache_chunks - ind = 1 - end - end) - if ind > 1 - @inbounds C[ind:bitcache_size] = false - dumpbitcache(Bc, cind, C) + $(Expr(:meta, :noinline)) + # destructure the indexmaps and As tuples + @nexprs $nargs i->(A_i = As[i]) + @nexprs $nargs i->(imap_i = indexmaps[i]) + @simd for I in CartesianRange(indices(B)) + # reverse-broadcast the indices + @nexprs $nargs i->(I_i = newindex(I, imap_i)) + # extract array values + @nexprs $nargs i->(@inbounds val_i = A_i[I_i]) + # call the function and store the result + @inbounds B[I] = @ncall $nargs f val end end end -# using start/next -function gen_broadcast_body_iter_tobitarray(nd::Int, narrays::Int, f) - F = Expr(:quote, f) +# For BitArray outputs, we cache the result in a "small" Vector{Bool}, +# and then copy in chunks into the output +@generated function _broadcast!{M,AT,nargs}(f, B::BitArray, indexmaps::M, As::AT, ::Type{Val{nargs}}) quote - @assert ndims(B) == $nd - @ncall $narrays check_broadcast_shape size(B) k->A_k - C = Array{Bool}(bitcache_size) + $(Expr(:meta, :noinline)) + # destructure the indexmaps and As tuples + @nexprs $nargs i->(A_i = As[i]) + @nexprs $nargs i->(imap_i = indexmaps[i]) + C = Vector{Bool}(bitcache_size) Bc = B.chunks ind = 1 cind = 1 - @nexprs 1 d->(@nexprs $narrays k->(state_k_0 = state_k_{$nd} = start(A_k))) - @nexprs $nd d->(@nexprs $narrays k->(skip_k_d = size(A_k, d) == 1)) - @nloops($nd, i, B, - d->(@nexprs $narrays k->(state_k_{d-1} = state_k_d)), # pre - d->(@nexprs $narrays k->(skip_k_d || (state_k_d = state_k_0))), # post - begin # body - @nexprs $narrays k->(@inbounds (v_k, state_k_0) = next(A_k, state_k_0)) - @inbounds C[ind] = (@ncall $narrays $F v) - ind += 1 - if ind > bitcache_size - dumpbitcache(Bc, cind, C) - cind += bitcache_chunks - ind = 1 - end - end) + @simd for I in CartesianRange(indices(B)) + # reverse-broadcast the indices + @nexprs $nargs i->(I_i = newindex(I, imap_i)) + # extract array values + @nexprs $nargs i->(@inbounds val_i = A_i[I_i]) + # call the function and store the result + @inbounds C[ind] = @ncall $nargs f val + ind += 1 + if ind > bitcache_size + dumpbitcache(Bc, cind, C) + cind += bitcache_chunks + ind = 1 + end + end if ind > 1 @inbounds C[ind:bitcache_size] = false dumpbitcache(Bc, cind, C) @@ -166,60 +132,17 @@ function gen_broadcast_body_iter_tobitarray(nd::Int, narrays::Int, f) end end -function gen_broadcast_function(genbody::Function, nd::Int, narrays::Int, f) - As = [Symbol("A_",i) for i = 1:narrays] - body = genbody(nd, narrays, f) - @eval let - local _F_ - function _F_(B, $(As...)) - $body - end - _F_ - end -end - -function gen_broadcast_function_tobitarray(genbody::Function, nd::Int, narrays::Int, f) - As = [Symbol("A_",i) for i = 1:narrays] - body = genbody(nd, narrays, f) - @eval let - local _F_ - function _F_(B::BitArray, $(As...)) - $body - end - _F_ - end -end - -for (Bsig, Asig, gbf, gbb) in - ((BitArray , Union{Array,BitArray,Number} , - :gen_broadcast_function_tobitarray, :gen_broadcast_body_iter_tobitarray ), - (Any , Union{Array,BitArray,Number} , - :gen_broadcast_function , :gen_broadcast_body_iter ), - (BitArray , Any , - :gen_broadcast_function_tobitarray, :gen_broadcast_body_cartesian_tobitarray), - (Any , Any , - :gen_broadcast_function , :gen_broadcast_body_cartesian )) - - @eval let cache = Dict{Any,Dict{Int,Dict{Int,Any}}}() - global broadcast! - function broadcast!(f, B::$Bsig, As::$Asig...) - nd = ndims(B) - narrays = length(As) - - cache_f = @get! cache f Dict{Int,Dict{Int,Any}}() - cache_f_na = @get! cache_f narrays Dict{Int,Any}() - func = @get! cache_f_na nd $gbf($gbb, nd, narrays, f) - - func(B, As...) - B - end - end # let broadcast_cache +@inline function broadcast!{nargs}(f, B::AbstractArray, As::Vararg{Any,nargs}) + check_broadcast_shape(shape(B), As...) + sz = size(B) + mapindex = map(x->newindexer(sz, x), As) + _broadcast!(f, B, mapindex, As, Val{nargs}) + B end +@inline broadcast(f, As...) = broadcast!(f, allocate_for(Array{promote_eltype_op(f, As...)}, As, broadcast_shape(As...)), As...) -broadcast(f, As...) = broadcast!(f, Array{promote_eltype_op(f, As...)}(broadcast_shape(As...)), As...) - -bitbroadcast(f, As...) = broadcast!(f, BitArray(broadcast_shape(As...)), As...) +@inline bitbroadcast(f, As...) = broadcast!(f, allocate_for(BitArray, As, broadcast_shape(As...)), As...) broadcast_getindex(src::AbstractArray, I::AbstractArray...) = broadcast_getindex!(Array{eltype(src)}(broadcast_shape(I...)), src, I...) @generated function broadcast_getindex!(dest::AbstractArray, src::AbstractArray, I::AbstractArray...) @@ -307,77 +230,83 @@ function .^(A::AbstractArray, B::AbstractArray) broadcast!(^, Array{promote_op(^, eltype(A), eltype(B))}(broadcast_shape(A, B)), A, B) end -## element-wise comparison operators returning BitArray ## - -for (f, scalarf, bitf, bitfbody) in ((:.==, :(==), :biteq , :(~a $ b)), - (:.< , :< , :bitlt , :(~a & b)), - (:.!=, :!= , :bitneq, :( a $ b)), - (:.<=, :<= , :bitle , :(~a | b))) - @eval begin - ($f)(A::AbstractArray, B::AbstractArray) = bitbroadcast($scalarf, A, B) - ($bitf)(a::UInt64, b::UInt64) = $bitfbody - function ($f)(A::AbstractArray{Bool}, B::AbstractArray{Bool}) - local shape - try - shape = promote_shape(size(A), size(B)) - catch - return bitbroadcast($scalarf, A, B) - end - F = BitArray(shape) - Fc = F.chunks - Ac = BitArray(A).chunks - Bc = BitArray(B).chunks - if !isempty(Ac) && !isempty(Bc) - for i = 1:length(Fc) - 1 - Fc[i] = ($bitf)(Ac[i], Bc[i]) - end - Fc[end] = ($bitf)(Ac[end], Bc[end]) & _msk_end(F) - end - return F +# ## element-wise comparison operators returning BitArray ## + +.==(A::AbstractArray, B::AbstractArray) = bitbroadcast(==, A, B) + .<(A::AbstractArray, B::AbstractArray) = bitbroadcast(<, A, B) +.!=(A::AbstractArray, B::AbstractArray) = bitbroadcast(!=, A, B) +.<=(A::AbstractArray, B::AbstractArray) = bitbroadcast(<=, A, B) + +function broadcast_bitarrays(scalarf, bitf, A::AbstractArray{Bool}, B::AbstractArray{Bool}) + local shape + try + shape = promote_shape(size(A), size(B)) + catch + return bitbroadcast(scalarf, A, B) + end + F = BitArray(shape) + Fc = F.chunks + Ac = BitArray(A).chunks + Bc = BitArray(B).chunks + if !isempty(Ac) && !isempty(Bc) + for i = 1:length(Fc) - 1 + Fc[i] = (bitf)(Ac[i], Bc[i]) + end + Fc[end] = (bitf)(Ac[end], Bc[end]) & _msk_end(F) + end + return F +end + +biteq(a::UInt64, b::UInt64) = ~a $ b +bitlt(a::UInt64, b::UInt64) = ~a & b +bitneq(a::UInt64, b::UInt64) = a $ b +bitle(a::UInt64, b::UInt64) = ~a | b + +.==(A::AbstractArray{Bool}, B::AbstractArray{Bool}) = broadcast_bitarrays(==, biteq, A, B) + .<(A::AbstractArray{Bool}, B::AbstractArray{Bool}) = broadcast_bitarrays(<, bitlt, A, B) +.!=(A::AbstractArray{Bool}, B::AbstractArray{Bool}) = broadcast_bitarrays(!=, bitneq, A, B) +.<=(A::AbstractArray{Bool}, B::AbstractArray{Bool}) = broadcast_bitarrays(<=, bitle, A, B) + +function bitcache(op, A, B, refA, refB, l::Int, ind::Int, C::Vector{Bool}) + left = l - ind + 1 + @inbounds begin + for j = 1:min(bitcache_size, left) + C[j] = (op)(refA(A, ind), refB(B, ind)) + ind += 1 end + C[left+1:bitcache_size] = false end + return ind end # note: the following are not broadcasting, but need to be defined here to avoid # ambiguity warnings -for (f, cachef, scalarf) in ((:.==, :bitcache_eq , :(==)), - (:.< , :bitcache_lt , :< ), - (:.!=, :bitcache_neq, :!= ), - (:.<=, :bitcache_le , :<= )) - @eval ($cachef)(A::AbstractArray, B::AbstractArray, l::Int, ind::Int, C::Vector{Bool}) = 0 - for (sigA, sigB, expA, expB, shape) in ((:Any, :AbstractArray, - :A, :(B[ind]), - :(size(B))), - (:AbstractArray, :Any, - :(A[ind]), :B, - :(size(A)))) +for (f, scalarf) in ((:.==, :(==)), + (:.< , :< ), + (:.!=, :!= ), + (:.<=, :<= )) + for (sigA, sigB, active, refA, refB) in ((:Any, :AbstractArray, :B, + :((A,ind)->A), :((B,ind)->B[ind])), + (:AbstractArray, :Any, :A, + :((A,ind)->A[ind]), :((B,ind)->B))) + shape = :(shape($active)) @eval begin - function ($cachef)(A::$sigA, B::$sigB, l::Int, ind::Int, C::Vector{Bool}) - left = l - ind + 1 - @inbounds begin - for j = 1:min(bitcache_size, left) - C[j] = ($scalarf)($expA, $expB) - ind += 1 - end - C[left+1:bitcache_size] = false - end - return ind - end function ($f)(A::$sigA, B::$sigB) - F = BitArray($shape) + P = allocate_for(BitArray, $active, $shape) + F = parent(P) l = length(F) l == 0 && return F Fc = F.chunks C = Array{Bool}(bitcache_size) - ind = 1 + ind = first(linearindices($active)) cind = 1 for i = 1:div(l + bitcache_size - 1, bitcache_size) - ind = ($cachef)(A, B, l, ind, C) + ind = bitcache($scalarf, A, B, $refA, $refB, l, ind, C) dumpbitcache(Fc, cind, C) cind += bitcache_chunks end - return F + return P end end end diff --git a/base/cartesian.jl b/base/cartesian.jl index a387ef1bb61e7..16d20243c574d 100644 --- a/base/cartesian.jl +++ b/base/cartesian.jl @@ -41,7 +41,7 @@ end function _nloops(N::Int, itersym::Symbol, arraysym::Symbol, args::Expr...) @gensym d - _nloops(N, itersym, :($d->1:size($arraysym, $d)), args...) + _nloops(N, itersym, :($d->indices($arraysym, $d)), args...) end function _nloops(N::Int, itersym::Symbol, rangeexpr::Expr, args::Expr...) diff --git a/base/combinatorics.jl b/base/combinatorics.jl index 575d5b777072a..dde0d857f6f38 100644 --- a/base/combinatorics.jl +++ b/base/combinatorics.jl @@ -147,6 +147,7 @@ function invperm(a::AbstractVector) end b end +invperm(a::Tuple) = (invperm([a...])...,) #XXX This function should be moved to Combinatorics.jl but is currently used by Base.DSP. """ diff --git a/base/datafmt.jl b/base/datafmt.jl index ee25ed408a2e8..a9bd3708560da 100644 --- a/base/datafmt.jl +++ b/base/datafmt.jl @@ -563,12 +563,11 @@ function writedlm(io::IO, a::AbstractMatrix, dlm; opts...) optsd = val_opts(opts) quotes = get(optsd, :quotes, true) pb = PipeBuffer() - nr = size(a, 1) - nc = size(a, 2) - for i = 1:nr # fixme (iter): improve if timholy/ArrayIteration.jl is merged into Base - for j = 1:nc + lastc = last(indices(a, 2)) + for i = indices(a, 1) + for j = indices(a, 2) writedlm_cell(pb, a[i, j], dlm, quotes) - j == nc ? write(pb,'\n') : print(pb,dlm) + j == lastc ? write(pb,'\n') : print(pb,dlm) end (nb_available(pb) > (16*1024)) && write(io, takebuf_array(pb)) end diff --git a/base/deprecated.jl b/base/deprecated.jl index a4e9b92e4d79f..ded568ac7785b 100644 --- a/base/deprecated.jl +++ b/base/deprecated.jl @@ -683,7 +683,7 @@ function hist2d!{HT}(H::AbstractArray{HT,2}, v::AbstractMatrix, if init fill!(H, zero(HT)) end - for i = 1:size(v,1) # fixme (iter): update when #15459 is done + for i = indices(v,1) x = searchsortedfirst(edg1, v[i,1]) - 1 y = searchsortedfirst(edg2, v[i,2]) - 1 if 1 <= x <= n && 1 <= y <= m @@ -713,6 +713,24 @@ hist2d(v::AbstractMatrix) = hist2d(v, sturges(size(v,1))) @deprecate(pointer_to_string(p::Ptr{UInt8}, own::Bool=false), unsafe_wrap(String, p, own)) +function checkbounds(::Type{Bool}, sz::Integer, i) + depwarn("checkbounds(Bool, size(A, d), i) is deprecated, use checkindex(Bool, indices(A, d), i).", :checkbounds) + checkbounds(Bool, 1:sz, i) +end +immutable FakeArray{T,N} <: AbstractArray{T,N} + dims::NTuple{N,Int} +end +size(A::FakeArray) = A.dims +function checkbounds{N,T}(::Type{Bool}, sz::NTuple{N,Integer}, I1::T, I...) + depwarn("checkbounds(Bool, size(A), I...) is deprecated, use checkbounds(Bool, A, I...).", :checkbounds) + checkbounds(Bool, FakeArray(sz), I1, I...) +end + +function first(::Colon) + depwarn("first(:) is no longer unambiguous, call Base._first(:, A, dim)", :first) + 1 +end + # During the 0.5 development cycle, do not add any deprecations below this line # To be deprecated in 0.6 diff --git a/base/docs/helpdb/Base.jl b/base/docs/helpdb/Base.jl index 0c77ff7f6602b..d92fa6d3a9ca7 100644 --- a/base/docs/helpdb/Base.jl +++ b/base/docs/helpdb/Base.jl @@ -4064,24 +4064,6 @@ Element-wise less-than-or-equals comparison operator. """ Base.:(.<=) -""" - checkbounds(array, indexes...) - -Throw an error if the specified indexes are not in bounds for the given array. Subtypes of -`AbstractArray` should specialize this method if they need to provide custom bounds checking -behaviors. -""" -checkbounds(array, indexes...) - -""" - checkbounds(::Type{Bool}, dimlength::Integer, index) - -Return a `Bool` describing if the given index is within the bounds of the given dimension -length. Custom types that would like to behave as indices for all arrays can extend this -method in order to provide a specialized bounds checking implementation. -""" -checkbounds(::Type{Bool}, ::Integer, index) - """ asec(x) @@ -6091,43 +6073,6 @@ i.e. a character whose category code begins with 'N'. """ isnumber -""" - similar(array, [element_type=eltype(array)], [dims=size(array)]) - -Create an uninitialized mutable array with the given element type and size, based upon the -given source array. The second and third arguments are both optional, defaulting to the -given array's `eltype` and `size`. The dimensions may be specified either as a single tuple -argument or as a series of integer arguments. - -Custom AbstractArray subtypes may choose which specific array type is best-suited to return -for the given element type and dimensionality. If they do not specialize this method, the -default is an `Array{element_type}(dims...)`. - -For example, `similar(1:10, 1, 4)` returns an uninitialized `Array{Int,2}` since ranges are -neither mutable nor support 2 dimensions: - - julia> similar(1:10, 1, 4) - 1×4 Array{Int64,2}: - 4419743872 4374413872 4419743888 0 - -Conversely, `similar(trues(10,10), 2)` returns an uninitialized `BitVector` with two -elements since `BitArray`s are both mutable and can support 1-dimensional arrays: - - julia> similar(trues(10,10), 2) - 2-element BitArray{1}: - false - false - -Since `BitArray`s can only store elements of type `Bool`, however, if you request a -different element type it will create a regular `Array` instead: - - julia> similar(falses(10), Float64, 2, 4) - 2×4 Array{Float64,2}: - 2.18425e-314 2.18425e-314 2.18425e-314 2.18425e-314 - 2.18425e-314 2.18425e-314 2.18425e-314 2.18425e-314 -""" -similar - """ copy(x) diff --git a/base/essentials.jl b/base/essentials.jl index f64418883a0e5..1d1ba215b1865 100644 --- a/base/essentials.jl +++ b/base/essentials.jl @@ -172,6 +172,8 @@ start(v::SimpleVector) = 1 next(v::SimpleVector,i) = (v[i],i+1) done(v::SimpleVector,i) = (i > v.length) isempty(v::SimpleVector) = (v.length == 0) +indices(v::SimpleVector, d) = d == 1 ? (1:length(v)) : (1:1) +linearindices(v::SimpleVector) = indices(v, 1) function ==(v1::SimpleVector, v2::SimpleVector) length(v1)==length(v2) || return false diff --git a/base/exports.jl b/base/exports.jl index 6e9dbd92a106e..68621d887f607 100644 --- a/base/exports.jl +++ b/base/exports.jl @@ -484,6 +484,7 @@ export zeta, # arrays + allocate_for, bitbroadcast, broadcast!, broadcast, @@ -491,6 +492,7 @@ export broadcast_setindex!, cat, checkbounds, + checkindex, circshift, clamp!, colon, @@ -526,6 +528,7 @@ export hvcat, ind2sub, indexin, + indices, indmax, indmin, invperm, @@ -535,6 +538,7 @@ export isperm, issorted, last, + linearindices, linspace, logspace, mapslices, @@ -579,6 +583,7 @@ export searchsortedlast, select!, select, + shape, shuffle, shuffle!, size, diff --git a/base/markdown/GitHub/table.jl b/base/markdown/GitHub/table.jl index b02d2ecaded51..e4e74f152a0f4 100644 --- a/base/markdown/GitHub/table.jl +++ b/base/markdown/GitHub/table.jl @@ -88,7 +88,7 @@ padding(width, twidth, a) = function padcells!(rows, align; len = length, min = 0) widths = colwidths(rows, len = len, min = min) - for i = 1:length(rows), j = 1:length(rows[1]) # fixme (iter): can we make indexing more general here? + for i = 1:length(rows), j = indices(rows[1],1) cell = rows[i][j] lpad, rpad = padding(len(cell), widths[j], align[j]) rows[i][j] = " "^lpad * cell * " "^rpad @@ -105,13 +105,13 @@ _dash(width, align) = function plain(io::IO, md::Table) cells = mapmap(plaininline, md.rows) padcells!(cells, md.align, len = length, min = 3) - for i = 1:length(cells) # fixme (iter): can we make indexing more general here? + for i = indices(cells,1) print(io, "| ") join(io, cells[i], " | ") println(io, " |") if i == 1 print(io, "|") - join(io, [_dash(length(cells[i][j]), md.align[j]) for j = 1:length(cells[1])], "|") + join(io, [_dash(length(cells[i][j]), md.align[j]) for j = indices(cells[1],1)], "|") println(io, "|") end end diff --git a/base/math.jl b/base/math.jl index 426f436c179eb..c7c0516a5658e 100644 --- a/base/math.jl +++ b/base/math.jl @@ -42,7 +42,7 @@ clamp{X,L,H}(x::X, lo::L, hi::H) = clamp{T}(x::AbstractArray{T,1}, lo, hi) = [clamp(xx, lo, hi) for xx in x] clamp{T}(x::AbstractArray{T,2}, lo, hi) = - [clamp(x[i,j], lo, hi) for i in 1:size(x,1), j in 1:size(x,2)] # fixme (iter): change to `eachindex` when #15459 is ready + [clamp(x[i,j], lo, hi) for i in indices(x,1), j in indices(x,2)] clamp{T}(x::AbstractArray{T}, lo, hi) = reshape([clamp(xx, lo, hi) for xx in x], size(x)) diff --git a/base/multidimensional.jl b/base/multidimensional.jl index ec6d31dcae18a..52b078a89c5bd 100644 --- a/base/multidimensional.jl +++ b/base/multidimensional.jl @@ -147,44 +147,17 @@ end # IteratorsMD using .IteratorsMD -# Bounds-checking specialization -# Specializing for a fixed number of arguments provides a ~25% -# improvement over the general definitions in abstractarray.jl - -# This is annoying, but we must first define logical indexing to avoid ambiguities -_internal_checkbounds(A::AbstractVector, I::AbstractVector{Bool}) = length(A) == length(I) || throw_boundserror(A, I) -_internal_checkbounds(A::AbstractVector, I::AbstractArray{Bool}) = size(A) == size(I) || throw_boundserror(A, I) - -for N = 1:5 - args = [:($(Symbol(:I, d))) for d = 1:N] - targs = [:($(Symbol(:I, d))::Union{Colon,Number,AbstractArray}) for d = 1:N] # prevent co-opting the CartesianIndex version - exs = [:(checkbounds(Bool, size(A, $d), $(args[d]))) for d = 1:N] - cbexpr = exs[1] - for d = 2:N - cbexpr = :($(exs[d]) & $cbexpr) - end - @eval begin - function checkbounds(A::AbstractArray, $(args...)) - @_inline_meta - _internal_checkbounds(A, $(args...)) - end - function _internal_checkbounds{T}(A::AbstractArray{T,$N}, $(targs...)) - @_inline_meta - ($cbexpr) || throw_boundserror(A, ($(args...),)) - end - end -end - -# Bounds-checking with CartesianIndex -@inline function checkbounds(::Type{Bool}, ::Tuple{}, I1::CartesianIndex) - checkbounds(Bool, (), I1.I...) -end -@inline function checkbounds(::Type{Bool}, sz::Tuple{}, I1::CartesianIndex, I...) - checkbounds(Bool, (), I1.I..., I...) -end -@inline function checkbounds(::Type{Bool}, sz::Dims, I1::CartesianIndex, I...) - checkbounds(Bool, sz, I1.I..., I...) -end +## Bounds-checking with CartesianIndex +# Ambiguity with linear indexing: +@inline _chkbnds(A::AbstractVector, checked::NTuple{1,Bool}, I::CartesianIndex) = _chkbnds(A, checked, I.I...) +@inline _chkbnds(A::AbstractArray, checked::NTuple{1,Bool}, I::CartesianIndex) = _chkbnds(A, checked, I.I...) +# Generic bounds checking +@inline _chkbnds{T,N}(A::AbstractArray{T,N}, checked::NTuple{N,Bool}, I1::CartesianIndex, I...) = _chkbnds(A, checked, I1.I..., I...) +@inline _chkbnds{T,N,M}(A::AbstractArray{T,N}, checked::NTuple{M,Bool}, I1::CartesianIndex, I...) = _chkbnds(A, checked, I1.I..., I...) + +@inline checkbounds_indices(::Tuple{}, I::Tuple{CartesianIndex,Vararg{Any}}) = checkbounds_indices((), (I[1].I..., tail(I)...)) +@inline checkbounds_indices(inds::Tuple{Any}, I::Tuple{CartesianIndex,Vararg{Any}}) = checkbounds_indices(inds, (I[1].I..., tail(I)...)) +@inline checkbounds_indices(inds::Tuple, I::Tuple{CartesianIndex,Vararg{Any}}) = checkbounds_indices(inds, (I[1].I..., tail(I)...)) # Recursively compute the lengths of a list of indices, without dropping scalars # These need to be inlined for more than 3 indexes @@ -194,22 +167,38 @@ index_lengths_dim(A, dim) = () index_lengths_dim(A, dim, ::Colon) = (trailingsize(A, dim),) @inline index_lengths_dim(A, dim, ::Colon, i, I...) = (size(A, dim), index_lengths_dim(A, dim+1, i, I...)...) @inline index_lengths_dim(A, dim, ::Real, I...) = (1, index_lengths_dim(A, dim+1, I...)...) -@inline index_lengths_dim{N}(A, dim, ::CartesianIndex{N}, I...) = (1, index_shape_dim(A, dim+N, I...)...) +@inline index_lengths_dim{N}(A, dim, ::CartesianIndex{N}, I...) = (1, index_lengths_dim(A, dim+N, I...)...) @inline index_lengths_dim(A, dim, i::AbstractArray, I...) = (length(i), index_lengths_dim(A, dim+1, I...)...) @inline index_lengths_dim(A, dim, i::AbstractArray{Bool}, I...) = (sum(i), index_lengths_dim(A, dim+1, I...)...) @inline index_lengths_dim{N}(A, dim, i::AbstractArray{CartesianIndex{N}}, I...) = (length(i), index_lengths_dim(A, dim+N, I...)...) # shape of array to create for getindex() with indexes I, dropping scalars -index_shape(A::AbstractArray, I::Colon) = (length(A),) -@inline index_shape(A::AbstractArray, I...) = index_shape_dim(A, 1, I...) -index_shape_dim(A, dim, I::Real...) = () -index_shape_dim(A, dim, ::Colon) = (trailingsize(A, dim),) -@inline index_shape_dim(A, dim, ::Colon, i, I...) = (size(A, dim), index_shape_dim(A, dim+1, i, I...)...) -@inline index_shape_dim(A, dim, ::Real, I...) = (index_shape_dim(A, dim+1, I...)...) -@inline index_shape_dim{N}(A, dim, ::CartesianIndex{N}, I...) = (index_shape_dim(A, dim+N, I...)...) -@inline index_shape_dim(A, dim, i::AbstractArray, I...) = (size(i)..., index_shape_dim(A, dim+1, I...)...) -@inline index_shape_dim(A, dim, i::AbstractArray{Bool}, I...) = (sum(i), index_shape_dim(A, dim+1, I...)...) -@inline index_shape_dim{N}(A, dim, i::AbstractArray{CartesianIndex{N}}, I...) = (size(i)..., index_shape_dim(A, dim+N, I...)...) +# Rather than use an Integer dim, we grow a tuple (true, true, ...) +# whose length is equal to the dimension we're to process next. This +# allows us to dispatch, which is important for the type-stability of +# the lines involving Colon as the final index. +index_shape(A::AbstractVector, I::Colon) = shape(A) +index_shape(A::AbstractArray, I::Colon) = (length(A),) +@inline index_shape(A::AbstractArray, I...) = index_shape_dim(A, (true,), I...) +@inline index_shape_dim(A, dim, ::Colon) = (trailingsize(A, length(dim)),) +@inline index_shape_dim{T,N}(A::AbstractArray{T,N}, dim::NTuple{N}, ::Colon) = (shape(A, N),) +@inline index_shape_dim(A, dim, I::Real...) = () +@inline index_shape_dim(A, dim, ::Colon, i, I...) = (shape(A, length(dim)), index_shape_dim(A, (dim...,true), i, I...)...) +@inline index_shape_dim(A, dim, ::Real, I...) = (index_shape_dim(A, (dim...,true), I...)...) +@inline index_shape_dim{N}(A, dim, ::CartesianIndex{N}, I...) = (index_shape_dim(A, (dim...,ntuple(d->true,Val{N})...), I...)...) +@inline index_shape_dim(A, dim, i::AbstractArray, I...) = (shape(i)..., index_shape_dim(A, (dim...,true), I...)...) +@inline index_shape_dim(A, dim, i::AbstractArray{Bool}, I...) = (sum(i), index_shape_dim(A, (dim...,true), I...)...) +@inline index_shape_dim{N}(A, dim, i::AbstractArray{CartesianIndex{N}}, I...) = (shape(i)..., index_shape_dim(A, (dim...,ntuple(d->true,Val{N})...), I...)...) + +@inline decolon(A::AbstractVector, ::Colon) = (indices(A,1),) +@inline decolon(A::AbstractArray, ::Colon) = (1:length(A),) +@inline decolon(A::AbstractArray, I...) = decolon_dim(A, (true,), I...) +@inline decolon_dim(A::AbstractArray, dim) = () +@inline decolon_dim{T,N}(A::AbstractArray{T,N}, dim::NTuple{N}, ::Colon) = (indices(A, N),) +@inline decolon_dim(A, dim, ::Colon) = (1:trailingsize(A, length(dim)),) +@inline decolon_dim(A::AbstractArray, dim, i1, I...) = (i1, decolon_dim(A, (dim...,true), I...)...) +@inline decolon_dim{N}(A::AbstractArray, dim, i1::AbstractArray{CartesianIndex{N}}, I...) = (i1, decolon_dim(A, (dim...,ntuple(d->true,Val{N})...), I...)...) +@inline decolon_dim(A::AbstractArray, dim, ::Colon, I...) = (indices(A, length(dim)), decolon_dim(A, (dim...,true), I...)...) ### From abstractarray.jl: Internal multidimensional indexing definitions ### # These are not defined on directly on getindex to avoid @@ -242,7 +231,7 @@ end @nexprs $N d->(I_d = to_index(I[d])) shape = @ncall $N index_shape A I dest = similar(A, shape) - size(dest) == shape || throw_checksize_error(dest, shape) + size(dest) == map(dimlength, shape) || throw_checksize_error(dest, shape) @ncall $N _unsafe_getindex! dest A I end end @@ -251,14 +240,14 @@ end function _unsafe_getindex(::LinearIndexing, src::AbstractArray, I::AbstractArray{Bool}) shape = index_shape(src, I) dest = similar(src, shape) - size(dest) == shape || throw_checksize_error(dest, shape) + size(dest) == map(dimlength, shape) || throw_checksize_error(dest, shape) D = eachindex(dest) Ds = start(D) for (b, s) in zip(I, eachindex(src)) - if b + @inbounds if b d, Ds = next(D, Ds) - @inbounds dest[d] = src[s] + dest[d] = src[s] end end dest @@ -288,17 +277,21 @@ end N = length(I) quote $(Expr(:meta, :inline)) + @nexprs $N d->(I_d = I[d]) + J = @ncall $N decolon src I + @nexprs $N d->(J_d = J[d]) D = eachindex(dest) Ds = start(D) - idxlens = index_lengths(src, I...) - @nloops $N i d->(1:idxlens[d]) d->(@inbounds j_d = getindex(I[d], i_d)) begin + @inbounds @nloops $N j d->J_d begin d, Ds = next(D, Ds) - @inbounds dest[d] = @ncall $N getindex src j + dest[d] = @ncall $N getindex src j end dest end end +dimlength(r::Range) = length(r) +dimlength(i::Integer) = i @noinline throw_checksize_error(A, sz) = throw(DimensionMismatch("output array is the wrong size; expected $sz, got $(size(A))")) ## setindex! ## @@ -331,12 +324,12 @@ function _unsafe_setindex!(::LinearIndexing, A::AbstractArray, x, I::AbstractArr X = _iterable(x) Xs = start(X) c = 0 - for (iA, i) in zip(eachindex(A), eachindex(I)) - @inbounds Ii = I[i] + @inbounds for (iA, i) in zip(eachindex(A), eachindex(I)) + Ii = I[i] if Ii done(X, Xs) && throw_setindex_mismatch(x, c+1) (v, Xs) = next(X, Xs) - @inbounds A[iA] = v + A[iA] = v c += 1 end end @@ -369,10 +362,12 @@ end @nexprs $N d->(I_d = I[d]) idxlens = @ncall $N index_lengths A I @ncall $N setindex_shape_check X (d->idxlens[d]) + J = @ncall $N decolon A I + @nexprs $N d->(J_d = J[d]) Xs = start(X) - @nloops $N i d->(1:idxlens[d]) d->(@inbounds j_d = I_d[i_d]) begin + @inbounds @nloops $N j d->J_d begin v, Xs = next(X, Xs) - @inbounds @ncall $N setindex! A v j + @ncall $N setindex! A v j end A end @@ -447,42 +442,48 @@ cumsum!(B, A::AbstractArray) = cumsum!(B, A, 1) cumprod(A::AbstractArray, axis::Integer=1) = cumprod!(similar(A), A, axis) cumprod!(B, A) = cumprod!(B, A, 1) -for (f, op) in ((:cumsum!, :+), - (:cumprod!, :*)) - @eval begin - @generated function ($f){T,N}(B, A::AbstractArray{T,N}, axis::Integer) - quote - if size(B, axis) < 1 - return B - end - size(B) == size(A) || throw(DimensionMismatch("Size of B must match A")) - if axis > N - copy!(B, A) - return B - end - if axis == 1 - # We can accumulate to a temporary variable, which allows register usage and will be slightly faster - @inbounds @nloops $N i d->(d > 1 ? (1:size(A,d)) : (1:1)) begin - tmp = convert(eltype(B), @nref($N, A, i)) - @nref($N, B, i) = tmp - for i_1 = 2:size(A,1) - tmp = ($($op))(tmp, @nref($N, A, i)) - @nref($N, B, i) = tmp - end - end - else - @nexprs $N d->(isaxis_d = axis == d) - # Copy the initial element in each 1d vector along dimension `axis` - @inbounds @nloops $N i d->(d == axis ? (1:1) : (1:size(A,d))) @nref($N, B, i) = @nref($N, A, i) - # Accumulate - @inbounds @nloops $N i d->((1+isaxis_d):size(A, d)) d->(j_d = i_d - isaxis_d) begin - @nref($N, B, i) = ($($op))(@nref($N, B, j), @nref($N, A, i)) - end - end - B +cumsum!(B, A, axis::Integer) = cumop!(+, B, A, axis) +cumprod!(B, A, axis::Integer) = cumop!(*, B, A, axis) + +function cumop!(op, B, A, axis::Integer) + if size(B, axis) < 1 + return B + end + indices(B) == indices(A) || throw(DimensionMismatch("Shape of B must match A")) + if axis > ndims(A) + copy!(B, A) + return B + end + if axis == 1 + # We can accumulate to a temporary variable, which allows register usage and will be slightly faster + ind1 = indices(A,1) + @inbounds for I in CartesianRange(tail(indices(A))) + tmp = convert(eltype(B), A[first(ind1), I]) + B[first(ind1), I] = tmp + for i_1 = first(ind1)+1:last(ind1) + tmp = op(tmp, A[i_1, I]) + B[i_1, I] = tmp end end + else + R1 = CartesianRange(indices(A)[1:axis-1]) # not type-stable + R2 = CartesianRange(indices(A)[axis+1:end]) + _cumop!(op, B, A, R1, indices(A, axis), R2) # use function barrier end + return B +end + +@noinline function _cumop!(op, B, A, R1, ind, R2) + # Copy the initial element in each 1d vector along dimension `axis` + i = first(ind) + @inbounds for J in R2, I in R1 + B[I, i, J] = A[I, i, J] + end + # Accumulate + @inbounds for J in R2, i in first(ind)+1:last(ind), I in R1 + B[I, i, J] = op(B[I, i-1, J], A[I, i, J]) + end + B end ### from abstractarray.jl @@ -496,13 +497,14 @@ function fill!{T}(A::AbstractArray{T}, x) end function copy!{T,N}(dest::AbstractArray{T,N}, src::AbstractArray{T,N}) - length(dest) >= length(src) || throw(BoundsError()) - for (Isrc, Idest) in zip(eachindex(src), eachindex(dest)) - @inbounds dest[Idest] = src[Isrc] + @boundscheck checkbounds(dest, indices(src)...) + for I in eachindex(linearindexing(src,dest), src) + @inbounds dest[I] = src[I] end dest end + ### BitArrays ## getindex @@ -510,7 +512,7 @@ end # contiguous multidimensional indexing: if the first dimension is a range, # we can get some performance from using copy_chunks! @inline function _unsafe_getindex!(X::BitArray, B::BitArray, I0::Union{UnitRange{Int},Colon}) - copy_chunks!(X.chunks, 1, B.chunks, first(I0), index_lengths(B, I0)[1]) + copy_chunks!(X.chunks, 1, B.chunks, _first(I0, B, :), index_lengths(B, I0)[1]) return X end @@ -521,7 +523,7 @@ end $(Expr(:meta, :inline)) @nexprs $N d->(I_d = I[d]) - f0 = first(I0) + f0 = _first(I0, B, 1) l0 = size(X, 1) gap_lst_1 = 0 @@ -531,7 +533,7 @@ end @nexprs $N d->begin stride *= size(B, d) stride_lst_d = stride - ind += stride * (first(I_d) - 1) + ind += stride * (_first(I_d, B, d) - 1) gap_lst_{d+1} *= stride end @@ -580,7 +582,7 @@ end l0 = index_lengths(B, I0)[1] setindex_shape_check(X, l0) l0 == 0 && return B - f0 = first(I0) + f0 = _first(I0, B, :) copy_to_bitarray_chunks!(B.chunks, f0, X, 1, l0) return B end @@ -590,7 +592,7 @@ end y = Bool(x) l0 = index_lengths(B, I0)[1] l0 == 0 && return B - f0 = first(I0) + f0 = _first(I0, B, :) fill_chunks!(B.chunks, y, f0, l0) return B end @@ -606,7 +608,7 @@ end idxlens = @ncall $N index_lengths B I0 d->I[d] @ncall $N setindex_shape_check X idxlens[1] d->idxlens[d+1] isempty(X) && return B - f0 = first(I0) + f0 = _first(I0, B, 1) l0 = idxlens[1] gap_lst_1 = 0 @@ -616,7 +618,7 @@ end @nexprs $N d->begin stride *= size(B, d) stride_lst_d = stride - ind += stride * (first(I[d]) - 1) + ind += stride * (_first(I[d], B, d) - 1) gap_lst_{d+1} *= stride end @@ -645,7 +647,7 @@ end y = Bool(x) idxlens = @ncall $N index_lengths B I0 d->I[d] - f0 = first(I0) + f0 = _first(I0, B, 1) l0 = idxlens[1] l0 == 0 && return B @nexprs $N d->(isempty(I[d]) && return B) @@ -657,7 +659,7 @@ end @nexprs $N d->begin stride *= size(B, d) stride_lst_d = stride - ind += stride * (first(I[d]) - 1) + ind += stride * (_first(I[d], B, d) - 1) gap_lst_{d+1} *= stride end @@ -722,12 +724,12 @@ function permutedims(B::StridedArray, perm) end function checkdims_perm{TP,TB,N}(P::AbstractArray{TP,N}, B::AbstractArray{TB,N}, perm) - dimsB = size(B) + indsB = indices(B) length(perm) == N || throw(ArgumentError("expected permutation of size $N, but length(perm)=$(length(perm))")) isperm(perm) || throw(ArgumentError("input is not a permutation")) - dimsP = size(P) + indsP = indices(P) for i = 1:length(perm) - dimsP[i] == dimsB[perm[i]] || throw(DimensionMismatch("destination tensor of incorrect size")) + indsP[i] == indsB[perm[i]] || throw(DimensionMismatch("destination tensor of incorrect size")) end nothing end @@ -784,7 +786,7 @@ If `dim` is specified, returns unique regions of the array `itr` along `dim`. @generated function unique{T,N}(A::AbstractArray{T,N}, dim::Int) quote 1 <= dim <= $N || return copy(A) - hashes = zeros(UInt, size(A, dim)) + hashes = allocate_for(inds->zeros(UInt, inds), A, shape(A, dim)) # Compute hash for each row k = 0 @@ -793,15 +795,15 @@ If `dim` is specified, returns unique regions of the array `itr` along `dim`. end # Collect index of first row for each hash - uniquerow = Array{Int}(size(A, dim)) + uniquerow = allocate_for(Array{Int}, A, shape(A, dim)) firstrow = Dict{Prehashed,Int}() - for k = 1:size(A, dim) # fixme (iter): use `eachindex(A, dim)` after #15459 is implemented + for k = indices(A, dim) uniquerow[k] = get!(firstrow, Prehashed(hashes[k]), k) end uniquerows = collect(values(firstrow)) # Check for collisions - collided = falses(size(A, dim)) + collided = allocate_for(falses, A, shape(A, dim)) @inbounds begin @nloops $N i A d->(if d == dim k = i_d @@ -816,11 +818,11 @@ If `dim` is specified, returns unique regions of the array `itr` along `dim`. end if any(collided) - nowcollided = BitArray(size(A, dim)) + nowcollided = allocate_for(BitArray, A, shape(A, dim)) while any(collided) # Collect index of first row for each collided hash empty!(firstrow) - for j = 1:size(A, dim) # fixme (iter): use `eachindex(A, dim)` after #15459 is implemented + for j = indices(A, dim) collided[j] || continue uniquerow[j] = get!(firstrow, Prehashed(hashes[j]), j) end @@ -847,6 +849,6 @@ If `dim` is specified, returns unique regions of the array `itr` along `dim`. end end - @nref $N A d->d == dim ? sort!(uniquerows) : (1:size(A, d)) + @nref $N A d->d == dim ? sort!(uniquerows) : (indices(A, d)) end end diff --git a/base/number.jl b/base/number.jl index ecfa2a749bc93..ab33caf0ebd25 100644 --- a/base/number.jl +++ b/base/number.jl @@ -6,6 +6,8 @@ isinteger(x::Integer) = true size(x::Number) = () size(x::Number,d) = convert(Int,d)<1 ? throw(BoundsError()) : 1 +indices(x::Number) = () +indices(x::Number,d) = convert(Int,d)<1 ? throw(BoundsError()) : (1:1) eltype{T<:Number}(::Type{T}) = T ndims(x::Number) = 0 ndims{T<:Number}(::Type{T}) = 0 diff --git a/base/operators.jl b/base/operators.jl index b488e31a2052b..6570a14e3e612 100644 --- a/base/operators.jl +++ b/base/operators.jl @@ -373,6 +373,23 @@ function promote_shape(a::Dims, b::Dims) return a end +function promote_shape(a::AbstractArray, b::AbstractArray) + if ndims(a) < ndims(b) + return promote_shape(b, a) + end + for i=1:ndims(b) + if indices(a, i) != indices(b, i) + throw(DimensionMismatch("dimensions must match")) + end + end + for i=ndims(b)+1:ndims(a) + if indices(a, i) != 1:1 + throw(DimensionMismatch("dimensions must match")) + end + end + return shape(a) +end + function throw_setindex_mismatch(X, I) if length(I) == 1 throw(DimensionMismatch("tried to assign $(length(X)) elements to $(I[1]) destinations")) diff --git a/base/permuteddimsarray.jl b/base/permuteddimsarray.jl index be951c87a6e8f..887c10898cf27 100644 --- a/base/permuteddimsarray.jl +++ b/base/permuteddimsarray.jl @@ -2,76 +2,98 @@ module PermutedDimsArrays +using Base: IndicesStartAt1, IndicesBehavior, indicesbehavior + export permutedims +# Some day we will want storage-order-aware iteration, so put perm in the parameters immutable PermutedDimsArray{T,N,AA<:AbstractArray,perm} <: AbstractArray{T,N} parent::AA + iperm::NTuple{N,Int} dims::NTuple{N,Int} - function PermutedDimsArray(data::AA, p::AbstractVector) + function PermutedDimsArray(data::AA) + # TODO optimize isperm & invperm for low dimensions? isa(perm, NTuple{N,Int}) || error("perm must be an NTuple{$N,Int}") (length(perm) == N && isperm(perm)) || throw(ArgumentError(string(perm, " is not a valid permutation of dimensions 1:", N))) - for i = 1:N - perm[p[i]] == i || throw(ArgumentError("size permutation must be the inverse of perm")) - end - new(data, ([size(data)...][p]...,)) + iperm = invperm(perm) + new(data, iperm, genperm(size(data), perm)) end end -function PermutedDimsArray{T,N}(data::AbstractArray{T,N}, p::AbstractVector) - length(p) == N || throw(ArgumentError(string(p, " is not a valid permutation of dimensions 1:", N))) - PermutedDimsArray{T,N,typeof(data),(invperm(p)...,)}(data, p) +function PermutedDimsArray{T,N}(data::AbstractArray{T,N}, perm) + length(perm) == N || throw(ArgumentError(string(p, " is not a valid permutation of dimensions 1:", N))) + PermutedDimsArray{T,N,typeof(data),(perm...,)}(data) end +Base.parent(A::PermutedDimsArray) = A.parent Base.size(A::PermutedDimsArray) = A.dims +Base.indices{T,N,AA,perm}(A::PermutedDimsArray{T,N,AA,perm}, d) = indices(parent(A), perm[d]) -@inline function Base.getindex{T,N,AA,perm}(A::PermutedDimsArray{T,N,AA,perm}, I::Int...) +@inline function Base.getindex{T,N}(A::PermutedDimsArray{T,N}, I::Vararg{Int,N}) @boundscheck checkbounds(A, I...) - @inbounds val = getindex(A.parent, _pda_reindex(perm, (), I)...) + @inbounds val = getindex(A.parent, genperm(I, A.iperm)...) val end -@inline function Base.setindex!{T,N,AA,perm}(A::PermutedDimsArray{T,N,AA,perm}, val, I::Int...) +@inline function Base.setindex!{T,N}(A::PermutedDimsArray{T,N}, val, I::Vararg{Int,N}) @boundscheck checkbounds(A, I...) - @inbounds setindex!(A.parent, val, _pda_reindex(perm, (), I)...) + @inbounds setindex!(A.parent, val, genperm(I, A.iperm)...) val end -@inline _pda_reindex(::Tuple{}, out, I) = out -@inline _pda_reindex(perm, out, I) = _pda_reindex(Base.tail(perm), (out..., I[perm[1]]), I) +# Could use ntuple(d->I[perm[d]], Val{N}) once #15276 is solved +@inline genperm{N}(I::NTuple{N}, perm::Dims{N}) = _genperm((), I, perm...) +_genperm(out, I) = out +@inline _genperm(out, I, p, perm...) = _genperm((out..., I[p]), I, perm...) +@inline genperm(I, perm::AbstractVector{Int}) = genperm(I, (perm...,)) function Base.permutedims{T,N}(A::AbstractArray{T,N}, perm) - sz::NTuple{N,Int} = size(A)[perm] - dest = similar(A, sz) + dest = similar_permute(A, perm) permutedims!(dest, A, perm) end +similar_permute(A::AbstractArray, perm) = similar_permute(indicesbehavior(A), A, perm) +similar_permute{T,N}(::IndicesStartAt1, A::AbstractArray{T,N}, perm) = similar(A, genperm(size(A), perm)) +similar_permute{T,N}(::IndicesBehavior, A::AbstractArray{T,N}, perm) = similar(A, genperm(indices(A), perm)) + function Base.permutedims!(dest, src::AbstractArray, perm) Base.checkdims_perm(dest, src, perm) P = PermutedDimsArray(dest, invperm(perm)) + _copy!(P, src) + return dest +end + +function Base.copy!{T,N}(dest::PermutedDimsArray{T,N}, src::AbstractArray{T,N}) + checkbounds(dest, indices(src)...) + _copy!(dest, src) +end +Base.copy!(dest::PermutedDimsArray, src::AbstractArray) = _copy!(dest, src) + +function _copy!{T,N,AA,perm}(P::PermutedDimsArray{T,N,AA,perm}, src) # If dest/src are "close to dense," then it pays to be cache-friendly. # Determine the first permuted dimension d = 0 # d+1 will hold the first permuted dimension of src while d < ndims(src) && perm[d+1] == d+1 d += 1 end - sz1 = size(src)[1:d] - if prod(sz1) > 1 - copy!(P, src) + if d == ndims(src) + copy!(parent(P), src) # it's not permuted else - R1 = CartesianRange(sz1) + R1 = CartesianRange(indices(src)[1:d]) d1 = findfirst(perm, d+1) # first permuted dim of dest - R2 = CartesianRange(size(src)[d+2:d1-1]) - R3 = CartesianRange(size(src)[d1+1:end]) + R2 = CartesianRange(indices(src)[d+2:d1-1]) + R3 = CartesianRange(indices(src)[d1+1:end]) _permutedims!(P, src, R1, R2, R3, d+1, d1) end - dest + return P end @noinline function _permutedims!(P::PermutedDimsArray, src, R1::CartesianRange{CartesianIndex{0}}, R2, R3, ds, dp) - for jo in 1:8:size(src, dp), io in 1:8:size(src, ds) + ip, is = indices(src, dp), indices(src, ds) + for jo in first(ip):8:last(ip), io in first(is):8:last(is) for I3 in R3, I2 in R2 - for j in jo:min(jo+7, size(src, dp)) - for i in io:min(io+7, size(src, ds)) + for j in jo:min(jo+7, last(ip)) + for i in io:min(io+7, last(is)) @inbounds P[i, I2, j, I3] = src[i, I2, j, I3] end end @@ -81,10 +103,11 @@ end end @noinline function _permutedims!(P::PermutedDimsArray, src, R1, R2, R3, ds, dp) - for jo in 1:8:size(src, dp), io in 1:8:size(src, ds) + ip, is = indices(src, dp), indices(src, ds) + for jo in first(ip):8:last(ip), io in first(is):8:last(is) for I3 in R3, I2 in R2 - for j in jo:min(jo+7, size(src, dp)) - for i in io:min(io+7, size(src, ds)) + for j in jo:min(jo+7, last(ip)) + for i in io:min(io+7, last(is)) for I1 in R1 @inbounds P[I1, i, I2, j, I3] = src[I1, i, I2, j, I3] end diff --git a/base/range.jl b/base/range.jl index 0135ebb720c8f..9dfe9a394060d 100644 --- a/base/range.jl +++ b/base/range.jl @@ -2,8 +2,8 @@ ## 1-dimensional ranges ## -typealias Dims Tuple{Vararg{Int}} -typealias DimsInteger Tuple{Vararg{Integer}} +typealias Dims{N} NTuple{N,Int} +typealias DimsInteger{N} NTuple{N,Integer} abstract Range{T} <: AbstractArray{T,1} @@ -820,5 +820,6 @@ function in(x, r::Range) n >= 1 && n <= length(r) && r[n] == x end +in{T<:Integer}(x::Integer, r::UnitRange{T}) = (first(r) <= x) & (x <= last(r)) in{T<:Integer}(x, r::Range{T}) = isinteger(x) && !isempty(r) && x>=minimum(r) && x<=maximum(r) && (mod(convert(T,x),step(r))-mod(first(r),step(r)) == 0) in(x::Char, r::Range{Char}) = !isempty(r) && x >= minimum(r) && x <= maximum(r) && (mod(Int(x) - Int(first(r)), step(r)) == 0) diff --git a/base/reduce.jl b/base/reduce.jl index 9ba63426ac9d4..679fa4b9bb7d8 100644 --- a/base/reduce.jl +++ b/base/reduce.jl @@ -231,13 +231,13 @@ sumabs2(a) = mapreduce(abs2, +, a) # Kahan (compensated) summation: O(1) error growth, at the expense # of a considerable increase in computational expense. function sum_kbn{T<:AbstractFloat}(A::AbstractArray{T}) - n = length(A) c = r_promote(+, zero(T)::T) - if n == 0 + if isempty(A) return c end - s = A[1] + c - for i in 2:n + inds = linearindices(A) + s = A[first(inds)] + c + for i in first(inds)+1:last(inds) @inbounds Ai = A[i] t = s + Ai if abs(s) >= abs(Ai) diff --git a/base/reducedim.jl b/base/reducedim.jl index 026db61f393bf..7f82a176cac98 100644 --- a/base/reducedim.jl +++ b/base/reducedim.jl @@ -166,6 +166,8 @@ function check_reducedims(R, A) sRi = size(R, i) sAi = size(A, i) if sRi == 1 + first(indices(R, i)) == first(indices(A, i)) || + throw(DimensionMismatch("Reduction along dimension $i must use lower indices")) if sAi > 1 if had_nonreduc lsiz = 0 # to reduce along i, but some previous dimensions were non-reducing @@ -174,8 +176,8 @@ function check_reducedims(R, A) end end else - sRi == sAi || - throw(DimensionMismatch("Reduction on array of size $(size(A)) with output of size $(size(R))")) + indices(R, i) == indices(A, i) || + throw(DimensionMismatch("Reduction on array with indices $(indices(A)) with output with indices $(indices(R))")) had_nonreduc = true end end @@ -190,29 +192,29 @@ function _mapreducedim!{T,N}(f, op, R::AbstractArray, A::AbstractArray{T,N}) if has_fast_linear_indexing(A) && lsiz > 16 # use mapreduce_impl, which is probably better tuned to achieve higher performance nslices = div(length(A), lsiz) - ibase = 0 + ibase = first(linearindices(A))-1 for i = 1:nslices @inbounds R[i] = op(R[i], mapreduce_impl(f, op, A, ibase+1, ibase+lsiz)) ibase += lsiz end - elseif size(R, 1) == 1 && sizA1 > 1 + return R + end + IRmax = dims_tail(map(last, indices(R)), A) + if size(R, 1) == 1 && sizA1 > 1 # keep the accumulator as a local variable when reducing along the first dimension - sizeR1 = size_skip1(size(R), A) - sizeA1 = size_skip1(size(A), A) - @inbounds for IA in CartesianRange(sizeA1) - IR = min(sizeR1, IA) - r = R[1,IR] - @simd for i = 1:size(A, 1) # fixme (iter): update when #15459 is implemented (and if it does't affect @simd) + i1 = first(indices(A, 1)) + @inbounds for IA in CartesianRange(tail(indices(A))) + IR = min(IA, IRmax) + r = R[i1,IR] + @simd for i in indices(A, 1) r = op(r, f(A[i, IA])) end - R[1,IR] = r + R[i1,IR] = r end else - sizeR1 = Base.size_skip1(size(R), A) - sizeA1 = Base.size_skip1(size(A), A) - @inbounds for IA in CartesianRange(sizeA1) - IR = min(IA, sizeR1) - @simd for i = 1:size(A, 1) # fixme (iter): update when #15459 is implemented (and if it does't affect @simd) + @inbounds for IA in CartesianRange(tail(indices(A))) + IR = min(IA, IRmax) + @simd for i in indices(A, 1) R[i,IR] = op(R[i,IR], f(A[i,IA])) end end @@ -271,22 +273,21 @@ end function findminmax!{T,N}(f, Rval, Rind, A::AbstractArray{T,N}) (isempty(Rval) || isempty(A)) && return Rval, Rind - (ndims(Rval) <= N && ndims(Rind) <= N) || throw(DimensionMismatch("Cannot find-reduce $(ndims(A))-dimensional array to $(ndims(Rval)),$(ndims(Rind)) dimensions")) + check_reducedims(Rval, A) for i = 1:N - (size(Rval, i) == size(A, i) || size(Rval, i) == 1) || throw(DimensionMismatch("Find-reduction on array of size $(size(A)) with output of size $(size(Rval))")) - size(Rval, i) == size(Rind, i) || throw(DimensionMismatch("Find-reduction: outputs must be of the same size")) + indices(Rval, i) == indices(Rind, i) || throw(DimensionMismatch("Find-reduction: outputs must have the same indices")) end # If we're reducing along dimension 1, for efficiency we can make use of a temporary. # Otherwise, keep the result in Rval/Rind so that we traverse A in storage order. - sizeR1 = size_skip1(size(Rval), A) - sizeA1 = size_skip1(size(A), A) + IRmax = dims_tail(map(last, indices(Rval)), A) k = 0 if size(Rval, 1) < size(A, 1) - @inbounds for IA in CartesianRange(sizeA1) - IR = min(sizeR1, IA) - tmpRv = Rval[1,IR] - tmpRi = Rind[1,IR] - for i = 1:size(A,1) # fixme (iter): update when #15459 is implemented (and if it does't affect @simd) + i1 = first(indices(A, 1)) + @inbounds for IA in CartesianRange(tail(indices(A))) + IR = min(IRmax, IA) + tmpRv = Rval[i1,IR] + tmpRi = Rind[i1,IR] + for i in indices(A,1) k += 1 tmpAv = A[i,IA] if f(tmpAv, tmpRv) @@ -294,13 +295,13 @@ function findminmax!{T,N}(f, Rval, Rind, A::AbstractArray{T,N}) tmpRi = k end end - Rval[1,IR] = tmpRv - Rind[1,IR] = tmpRi + Rval[i1,IR] = tmpRv + Rind[i1,IR] = tmpRi end else - @inbounds for IA in CartesianRange(sizeA1) - IR = min(sizeR1, IA) - for i = 1:size(A, 1) # fixme (iter): update when #15459 is implemented (and if it does't affect @simd) + @inbounds for IA in CartesianRange(tail(indices(A))) + IR = min(IRmax, IA) + for i in indices(A, 1) k += 1 tmpAv = A[i,IA] if f(tmpAv, Rval[i,IR]) @@ -358,7 +359,6 @@ function findmax{T}(A::AbstractArray{T}, region) zeros(Int, reduced_dims0(A, region)), A) end -size_skip1{T}(dims::Tuple{}, Aref::AbstractArray{T,0}) = CartesianIndex(()) -size_skip1{T,N}(dims::NTuple{N,Int}, Aref::AbstractArray{T,N}) = CartesianIndex(skip1(dims...)) -@inline size_skip1{T,M,N}(dims::NTuple{M,Int}, Aref::AbstractArray{T,N}) = size_skip1(tuple(dims..., 1), Aref) -skip1(x, t...) = t +dims_tail{T}(dims::Tuple{}, Aref::AbstractArray{T,0}) = CartesianIndex(()) +dims_tail{T,N}(dims::NTuple{N,Int}, Aref::AbstractArray{T,N}) = CartesianIndex(tail(dims)) +@inline dims_tail{T,M,N}(dims::NTuple{M,Int}, Aref::AbstractArray{T,N}) = dims_tail(tuple(dims..., 1), Aref) diff --git a/base/reshapedarray.jl b/base/reshapedarray.jl index b3746a143bc07..a28a7fb2d4721 100644 --- a/base/reshapedarray.jl +++ b/base/reshapedarray.jl @@ -36,6 +36,10 @@ start(R::ReshapedArrayIterator) = start(R.iter) end length(R::ReshapedArrayIterator) = length(R.iter) +reshape(parent::AbstractArray, ref::AbstractArray) = reshape(indicesbehavior(ref), parent, ref) +reshape(::IndicesStartAt1, parent::AbstractArray, ref::AbstractArray) = reshape(parent, size(ref)) +reshape(::IndicesBehavior, parent::AbstractArray, ref::AbstractArray) = reshape(parent, indices(ref)) + reshape(parent::AbstractArray, dims::Dims) = _reshape(parent, dims) reshape(parent::AbstractArray, len::Integer) = reshape(parent, (Int(len),)) reshape(parent::AbstractArray, dims::Int...) = reshape(parent, dims) diff --git a/base/serialize.jl b/base/serialize.jl index 65ed2d5654710..07281d1a0aaee 100644 --- a/base/serialize.jl +++ b/base/serialize.jl @@ -220,7 +220,10 @@ end trimmedsize(V) = index_lengths(V.parent, V.indexes...) -_trimmedsubarray{T,N,P,I,LD}(A, V::SubArray{T,N,P,I,LD}, newindexes) = SubArray{T,N,P,I,LD}(A, newindexes, size(V), 1, 1) +function _trimmedsubarray{T,N,P,I,LD}(A, V::SubArray{T,N,P,I,LD}, newindexes) + LD && return SubArray{T,N,P,I,LD}(A, newindexes, size(V), Base.compute_offset1(A, 1, newindexes), 1) + SubArray{T,N,P,I,LD}(A, newindexes, size(V), 0, 0) +end _trimmedsubarray(A, V, newindexes, index::ViewIndex, indexes...) = _trimmedsubarray(A, V, (newindexes..., trimmedindex(V.parent, length(newindexes)+1, index)), indexes...) trimmedindex(P, d, i::Real) = oftype(i, 1) diff --git a/base/show.jl b/base/show.jl index 85a450aa1baf0..deacf6191ec02 100644 --- a/base/show.jl +++ b/base/show.jl @@ -290,8 +290,8 @@ function show(io::IO, l::LambdaInfo) show(lambda_io, body) end -function show_delim_array(io::IO, itr::Union{AbstractArray,SimpleVector}, op, delim, cl, delim_one, - i1=1, l=length(itr)) +function show_delim_array(io::IO, itr::Union{AbstractArray,SimpleVector}, op, delim, cl, + delim_one, i1=first(linearindices(itr)), l=last(linearindices(itr))) print(io, op) if !show_circular(io, itr) recur_io = IOContext(io, SHOWN_SET=itr, multiline=false) @@ -301,7 +301,7 @@ function show_delim_array(io::IO, itr::Union{AbstractArray,SimpleVector}, op, de newline = true first = true i = i1 - if l > 0 + if l >= i1 while true if !isassigned(itr, i) print(io, undef_ref_str) @@ -313,7 +313,7 @@ function show_delim_array(io::IO, itr::Union{AbstractArray,SimpleVector}, op, de show(recur_io, x) end i += 1 - if i > i1+l-1 + if i > l delim_one && first && print(io, delim) break end @@ -1234,6 +1234,7 @@ is specified as string sep. function print_matrix_row(io::IO, X::AbstractVecOrMat, A::Vector, i::Integer, cols::AbstractVector, sep::AbstractString) + isempty(A) || first(indices(cols,1)) == 1 || throw(DimensionMismatch("indices of cols ($(indices(cols,1))) must start at 1")) for k = 1:length(A) j = cols[k] if isassigned(X,Int(i),Int(j)) # isassigned accepts only `Int` indices @@ -1303,25 +1304,30 @@ function print_matrix(io::IO, X::AbstractVecOrMat, @assert strwidth(hdots) == strwidth(ddots) sepsize = length(sep) m, n = size(X,1), size(X,2) + rowsA, colsA = collect(indices(X,1)), collect(indices(X,2)) # To figure out alignments, only need to look at as many rows as could # fit down screen. If screen has at least as many rows as A, look at A. # If not, then we only need to look at the first and last chunks of A, # each half a screen height in size. halfheight = div(screenheight,2) - rowsA = m <= screenheight ? (1:m) : [1:halfheight; m-div(screenheight-1,2)+1:m] + if m > screenheight + rowsA = [rowsA[1:halfheight]; rowsA[m-div(screenheight-1,2)+1:m]] + end # Similarly for columns, only necessary to get alignments for as many # columns as could conceivably fit across the screen maxpossiblecols = div(screenwidth, 1+sepsize) - colsA = n <= maxpossiblecols ? (1:n) : [1:maxpossiblecols; (n-maxpossiblecols+1):n] + if n > maxpossiblecols + colsA = [colsA[1:maxpossiblecols]; colsA[(n-maxpossiblecols+1):n]] + end A = alignment(io, X, rowsA, colsA, screenwidth, screenwidth, sepsize) # Nine-slicing is accomplished using print_matrix_row repeatedly if m <= screenheight # rows fit vertically on screen if n <= length(A) # rows and cols fit so just print whole matrix in one piece for i in rowsA - print(io, i == 1 ? pre : presp) + print(io, i == first(rowsA) ? pre : presp) print_matrix_row(io, X,A,i,colsA,sep) - print(io, i == m ? post : postsp) - if i != m; println(io, ); end + print(io, i == last(rowsA) ? post : postsp) + if i != last(rowsA); println(io, ); end end else # rows fit down screen but cols don't, so need horizontal ellipsis c = div(screenwidth-length(hdots)+1,2)+1 # what goes to right of ellipsis @@ -1329,25 +1335,25 @@ function print_matrix(io::IO, X::AbstractVecOrMat, c = screenwidth - sum(map(sum,Ralign)) - (length(Ralign)-1)*sepsize - length(hdots) Lalign = alignment(io, X, rowsA, colsA, c, c, sepsize) # alignments for left of ellipsis for i in rowsA - print(io, i == 1 ? pre : presp) - print_matrix_row(io, X,Lalign,i,1:length(Lalign),sep) - print(io, i % hmod == 1 ? hdots : repeat(" ", length(hdots))) + print(io, i == first(rowsA) ? pre : presp) + print_matrix_row(io, X,Lalign,i,colsA[1:length(Lalign)],sep) + print(io, (i - first(rowsA)) % hmod == 0 ? hdots : repeat(" ", length(hdots))) print_matrix_row(io, X,Ralign,i,n-length(Ralign)+colsA,sep) - print(io, i == m ? post : postsp) - if i != m; println(io, ); end + print(io, i == last(rowsA) ? post : postsp) + if i != last(rowsA); println(io, ); end end end else # rows don't fit so will need vertical ellipsis if n <= length(A) # rows don't fit, cols do, so only vertical ellipsis for i in rowsA - print(io, i == 1 ? pre : presp) + print(io, i == first(rowsA) ? pre : presp) print_matrix_row(io, X,A,i,colsA,sep) - print(io, i == m ? post : postsp) + print(io, i == last(rowsA) ? post : postsp) if i != rowsA[end]; println(io, ); end - if i == halfheight - print(io, i == 1 ? pre : presp) + if i == rowsA[halfheight] + print(io, i == first(rowsA) ? pre : presp) print_matrix_vdots(io, vdots,A,sep,vmod,1) - println(io, i == m ? post : postsp) + println(io, i == last(rowsA) ? post : postsp) end end else # neither rows nor cols fit, so use all 3 kinds of dots @@ -1357,18 +1363,18 @@ function print_matrix(io::IO, X::AbstractVecOrMat, Lalign = alignment(io, X, rowsA, colsA, c, c, sepsize) r = mod((length(Ralign)-n+1),vmod) # where to put dots on right half for i in rowsA - print(io, i == 1 ? pre : presp) - print_matrix_row(io, X,Lalign,i,1:length(Lalign),sep) - print(io, i % hmod == 1 ? hdots : repeat(" ", length(hdots))) + print(io, i == first(rowsA) ? pre : presp) + print_matrix_row(io, X,Lalign,i,colsA[1:length(Lalign)],sep) + print(io, (i - first(rowsA)) % hmod == 0 ? hdots : repeat(" ", length(hdots))) print_matrix_row(io, X,Ralign,i,n-length(Ralign)+colsA,sep) - print(io, i == m ? post : postsp) + print(io, i == last(rowsA) ? post : postsp) if i != rowsA[end]; println(io, ); end - if i == halfheight - print(io, i == 1 ? pre : presp) + if i == rowsA[halfheight] + print(io, i == first(rowsA) ? pre : presp) print_matrix_vdots(io, vdots,Lalign,sep,vmod,1) print(io, ddots) print_matrix_vdots(io, vdots,Ralign,sep,vmod,r) - println(io, i == m ? post : postsp) + println(io, i == last(rowsA) ? post : postsp) end end end @@ -1401,18 +1407,20 @@ function show_nd(io::IO, a::AbstractArray, print_matrix, label_slices) if isempty(a) return end - tail = size(a)[3:end] + tail = indices(a)[3:end] nd = ndims(a)-2 for I in CartesianRange(tail) idxs = I.I if limit for i = 1:nd ii = idxs[i] - if size(a,i+2) > 10 - if ii == 4 && all(x->x==1,idxs[1:i-1]) + ind = tail[i] + if length(ind) > 10 + if ii == ind[4] && all(d->idxs[d]==first(tail[d]),1:i-1) for j=i+1:nd szj = size(a,j+2) - if szj>10 && 3 < idxs[j] <= szj-3 + indj = tail[j] + if szj>10 && first(indj)+2 < idxs[j] <= last(indj)-3 @goto skip end end @@ -1420,7 +1428,7 @@ function show_nd(io::IO, a::AbstractArray, print_matrix, label_slices) print(io, "...\n\n") @goto skip end - if 3 < ii <= size(a,i+2)-3 + if ind[3] < ii <= ind[end-3] @goto skip end end @@ -1431,9 +1439,9 @@ function show_nd(io::IO, a::AbstractArray, print_matrix, label_slices) for i = 1:(nd-1); print(io, "$(idxs[i]), "); end println(io, idxs[end], "] =") end - slice = sub(a, 1:size(a,1), 1:size(a,2), idxs...) + slice = sub(a, indices(a,1), indices(a,2), idxs...) print_matrix(io, slice) - print(io, idxs == tail ? "" : "\n\n") + print(io, idxs == map(last,tail) ? "" : "\n\n") @label skip end end @@ -1449,15 +1457,15 @@ function print_matrix_repr(io, X::AbstractArray) end nr, nc = size(X,1), size(X,2) rdots, cdots = false, false - rr1, rr2 = 1:nr, 1:0 - cr1, cr2 = 1:nc, 1:0 + rr1, rr2 = indices(X,1), 1:0 + cr1, cr2 = indices(X,2), 1:0 if limit if nr > 4 - rr1, rr2 = 1:2, nr-1:nr + rr1, rr2 = rr1[1:2], rr1[nr-1:nr] rdots = true end if nc > 4 - cr1, cr2 = 1:2, nc-1:nc + cr1, cr2 = cr1[1:2], cr1[nc-1:nc] cdots = true end end @@ -1474,8 +1482,8 @@ function print_matrix_repr(io, X::AbstractArray) show(io, el) end end - if last(cr) == nc - i < nr && print(io, "; ") + if last(cr) == last(indices(X,2)) + i < last(indices(X,1)) && print(io, "; ") elseif cdots print(io, " \u2026 ") end @@ -1577,14 +1585,16 @@ function show_vector(io::IO, v, opn, cls) end print(io, prefix) if limited && length(v) > 20 - show_delim_array(io, v, opn, ",", "", false, 1, 10) + inds = _indices1(v) + show_delim_array(io, v, opn, ",", "", false, inds[1], inds[1]+9) print(io, " \u2026 ") - n = length(v) - show_delim_array(io, v, "", ",", cls, false, n-9, 10) + show_delim_array(io, v, "", ",", cls, false, inds[end-9], inds[end]) else show_delim_array(io, v, opn, ",", cls, false) end end +_indices1(v::AbstractArray) = indices(v,1) +_indices1(iter) = 1:length(iter) # printing BitArrays diff --git a/base/sort.jl b/base/sort.jl index 62ffe0c9e8ee1..1abe23f6fdbac 100644 --- a/base/sort.jl +++ b/base/sort.jl @@ -2,7 +2,7 @@ module Sort -using Base.Order, Base.copymutable +using Base: Order, copymutable, linearindices, allocate_for, shape, linearindexing, viewindexing, LinearFast import Base.sort, @@ -399,7 +399,10 @@ end defalg(v::AbstractArray) = DEFAULT_STABLE defalg{T<:Number}(v::AbstractArray{T}) = DEFAULT_UNSTABLE -sort!(v::AbstractVector, alg::Algorithm, order::Ordering) = sort!(v,1,length(v),alg,order) +function sort!(v::AbstractVector, alg::Algorithm, order::Ordering) + inds = indices(v,1) + sort!(v,first(inds),last(inds),alg,order) +end function sort!(v::AbstractVector; alg::Algorithm=defalg(v), @@ -426,7 +429,7 @@ function selectperm!{I<:Integer}(ix::AbstractVector{I}, v::AbstractVector, order::Ordering=Forward, initialized::Bool=false) if !initialized - @inbounds for i = 1:length(ix) + @inbounds for i = indices(ix,1) ix[i] = i end end @@ -444,7 +447,11 @@ function sortperm(v::AbstractVector; by=identity, rev::Bool=false, order::Ordering=Forward) - sort!(collect(1:length(v)), alg, Perm(ord(lt,by,rev,order),v)) + p = Base.allocate_for(Vector{Int}, v, shape(v, 1)) + for (i,ind) in zip(eachindex(p), indices(v, 1)) + p[i] = ind + end + sort!(p, alg, Perm(ord(lt,by,rev,order),v)) end function sortperm!{I<:Integer}(x::AbstractVector{I}, v::AbstractVector; @@ -454,13 +461,11 @@ function sortperm!{I<:Integer}(x::AbstractVector{I}, v::AbstractVector; rev::Bool=false, order::Ordering=Forward, initialized::Bool=false) - lx = length(x) - lv = length(v) - if lx != lv - throw(ArgumentError("index vector must be the same length as the source vector, $lx != $lv")) + if indices(x,1) != indices(v,1) + throw(ArgumentError("index vector must be the same length as the source vector, $(indices(x,1)) != $(indices(v,1))")) end if !initialized - @inbounds for i = 1:length(v) + @inbounds for i = indices(v,1) x[i] = i end end @@ -487,31 +492,48 @@ function sort(A::AbstractArray, dim::Integer; else Av = A[:] sort_chunks!(Av, size(A,1), alg, order) - reshape(Av, size(A)) + reshape(Av, A) end end @noinline function sort_chunks!(Av, n, alg, order) - for s = 1:n:length(Av) + inds = linearindices(Av) + for s = first(inds):n:last(inds) sort!(Av, s, s+n-1, alg, order) end Av end function sortrows(A::AbstractMatrix; kws...) - c = 1:size(A,2) - rows = [ sub(A,i,c) for i=1:size(A,1) ] # fixme (iter): update when #15459 is done + inds = indices(A,1) + T = slicetypeof(A, inds, :) + rows = allocate_for(Vector{T}, A, shape(A, 1)) + for i in inds + rows[i] = slice(A, i, :) + end p = sortperm(rows; kws..., order=Lexicographic) A[p,:] end function sortcols(A::AbstractMatrix; kws...) - r = 1:size(A,1) - cols = [ sub(A,r,i) for i=1:size(A,2) ] # fixme (iter): update when #15459 is done + inds = indices(A,2) + T = slicetypeof(A, :, inds) + cols = allocate_for(Vector{T}, A, shape(A, 2)) + for i in inds + cols[i] = slice(A, :, i) + end p = sortperm(cols; kws..., order=Lexicographic) A[:,p] end +function slicetypeof{T,N}(A::AbstractArray{T,N}, i1, i2) + I = (slice_dummy(i1),slice_dummy(i2)) + fast = isa(linearindexing(viewindexing(I), linearindexing(A)), LinearFast) + SubArray{T,1,typeof(A),typeof(I),fast} +end +slice_dummy(::Colon) = Colon() +slice_dummy{T}(::UnitRange{T}) = one(T) + ## fast clever sorting for floats ## module Float @@ -539,7 +561,7 @@ lt{T<:Floats}(::Right, x::T, y::T) = slt_int(unbox(T,x),unbox(T,y)) isnan(o::DirectOrdering, x::Floats) = (x!=x) isnan(o::Perm, i::Int) = isnan(o.order,o.data[i]) -function nans2left!(v::AbstractVector, o::Ordering, lo::Int=1, hi::Int=length(v)) +function nans2left!(v::AbstractVector, o::Ordering, lo::Int=first(indices(v,1)), hi::Int=last(indices(v,1))) i = lo @inbounds while i <= hi && isnan(o,v[i]) i += 1 @@ -554,7 +576,7 @@ function nans2left!(v::AbstractVector, o::Ordering, lo::Int=1, hi::Int=length(v) end return i, hi end -function nans2right!(v::AbstractVector, o::Ordering, lo::Int=1, hi::Int=length(v)) +function nans2right!(v::AbstractVector, o::Ordering, lo::Int=first(indices(v,1)), hi::Int=last(indices(v,1))) i = hi @inbounds while lo <= i && isnan(o,v[i]) i -= 1 @@ -595,7 +617,7 @@ end fpsort!(v::AbstractVector, a::Sort.PartialQuickSort, o::Ordering) = - sort!(v, 1, length(v), a, o) + sort!(v, first(indices(v,1)), last(indices(v,1)), a, o) sort!{T<:Floats}(v::AbstractVector{T}, a::Algorithm, o::DirectOrdering) = fpsort!(v,a,o) sort!{O<:DirectOrdering,T<:Floats}(v::Vector{Int}, a::Algorithm, o::Perm{O,Vector{T}}) = fpsort!(v,a,o) diff --git a/base/sparse/cholmod.jl b/base/sparse/cholmod.jl index 44934977b7145..ca90b2aa521b5 100644 --- a/base/sparse/cholmod.jl +++ b/base/sparse/cholmod.jl @@ -959,7 +959,11 @@ function convert{T}(::Type{Matrix{T}}, D::Dense{T}) a = Array{T}(s.nrow, s.ncol) copy!(a, D) end -function Base.copy!(dest::AbstractArray, D::Dense) + +Base.copy!(dest::Base.PermutedDimsArrays.PermutedDimsArray, src::Dense) = _copy!(dest, src) # ambig +Base.copy!(dest::AbstractArray, D::Dense) = _copy!(dest, D) + +function _copy!(dest::AbstractArray, D::Dense) s = unsafe_load(D.p) n = s.nrow*s.ncol n <= length(dest) || throw(BoundsError(dest, n)) diff --git a/base/stacktraces.jl b/base/stacktraces.jl index 341cd4f0f8d20..d7fce2c5365f7 100644 --- a/base/stacktraces.jl +++ b/base/stacktraces.jl @@ -182,7 +182,7 @@ function show_spec_linfo(io::IO, frame::StackFrame) end first = true print(io, '(') - for i = 2:length(params) # fixme (iter): `eachindex` with offset? + for i = 2:length(params) first || print(io, ", ") first = false print(io, "::", params[i]) diff --git a/base/statistics.jl b/base/statistics.jl index 197cb795e3796..b474af31fed46 100644 --- a/base/statistics.jl +++ b/base/statistics.jl @@ -91,41 +91,43 @@ centralize_sumabs2(A::AbstractArray, m::Number) = centralize_sumabs2(A::AbstractArray, m::Number, ifirst::Int, ilast::Int) = mapreduce_impl(centralizedabs2fun(m), +, A, ifirst, ilast) -@generated function centralize_sumabs2!{S,T,N}(R::AbstractArray{S}, A::AbstractArray{T,N}, means::AbstractArray) - quote - # following the implementation of _mapreducedim! at base/reducedim.jl - lsiz = check_reducedims(R,A) - isempty(R) || fill!(R, zero(S)) - isempty(A) && return R - @nextract $N sizeR d->size(R,d) - sizA1 = size(A, 1) - - if has_fast_linear_indexing(A) && lsiz > 16 - # use centralize_sumabs2, which is probably better tuned to achieve higher performance - nslices = div(length(A), lsiz) - ibase = 0 - for i = 1:nslices - @inbounds R[i] = centralize_sumabs2(A, means[i], ibase+1, ibase+lsiz) - ibase += lsiz - end - elseif size(R, 1) == 1 && sizA1 > 1 - # keep the accumulator as a local variable when reducing along the first dimension - @nloops $N i d->(d>1? (1:size(A,d)) : (1:1)) d->(j_d = sizeR_d==1 ? 1 : i_d) begin - @inbounds r = (@nref $N R j) - @inbounds m = (@nref $N means j) - for i_1 = 1:sizA1 # fixme (iter): change when #15459 is done - @inbounds r += abs2((@nref $N A i) - m) - end - @inbounds (@nref $N R j) = r +function centralize_sumabs2!{S,T,N}(R::AbstractArray{S}, A::AbstractArray{T,N}, means::AbstractArray) + # following the implementation of _mapreducedim! at base/reducedim.jl + lsiz = check_reducedims(R,A) + isempty(R) || fill!(R, zero(S)) + isempty(A) && return R + sizA1 = size(A, 1) + + if has_fast_linear_indexing(A) && lsiz > 16 + nslices = div(length(A), lsiz) + ibase = first(linearindices(A))-1 + for i = 1:nslices + @inbounds R[i] = centralize_sumabs2(A, means[i], ibase+1, ibase+lsiz) + ibase += lsiz + end + return R + end + IRmax = dims_tail(map(last, indices(R)), A) + if size(R, 1) == 1 && sizA1 > 1 + i1 = first(indices(A, 1)) + @inbounds for IA in CartesianRange(tail(indices(A))) + IR = min(IA, IRmax) + r = R[i1,IR] + m = means[i1,IR] + @simd for i in indices(A, 1) + r += abs2(A[i,IA] - m) end - else - # general implementation - @nloops $N i A d->(j_d = sizeR_d==1 ? 1 : i_d) begin - @inbounds (@nref $N R j) += abs2((@nref $N A i) - (@nref $N means j)) + R[i1,IR] = r + end + else + @inbounds for IA in CartesianRange(tail(indices(A))) + IR = min(IA, IRmax) + @simd for i in indices(A, 1) + R[i,IR] += abs2(A[i,IA] - means[i,IR]) end end - return R end + return R end function varm{T}(A::AbstractArray{T}, m::Number; corrected::Bool=true) diff --git a/base/subarray.jl b/base/subarray.jl index cd41df6ff50ec..aea6dcf80bfdd 100644 --- a/base/subarray.jl +++ b/base/subarray.jl @@ -9,11 +9,11 @@ immutable SubArray{T,N,P,I,L} <: AbstractArray{T,N} parent::P indexes::I dims::NTuple{N,Int} - first_index::Int # for linear indexing and pointer, only valid when L==true + offset1::Int # for linear indexing and pointer, only valid when L==true stride1::Int # used only for linear indexing - function SubArray(parent, indexes, dims, first_index, stride1) + function SubArray(parent, indexes, dims, offset1, stride1) check_parent_index_match(parent, indexes) - new(parent, indexes, dims, first_index, stride1) + new(parent, indexes, dims, offset1, stride1) end end # Compute the linear indexability of the indices, and combine it with the linear indexing of the parent @@ -24,8 +24,9 @@ function SubArray{P, I, N}(::LinearSlow, parent::P, indexes::I, dims::NTuple{N, SubArray{eltype(P), N, P, I, false}(parent, indexes, dims, 0, 0) end function SubArray{P, I, N}(::LinearFast, parent::P, indexes::I, dims::NTuple{N, Int}) - # Compute the first index and stride - SubArray{eltype(P), N, P, I, true}(parent, indexes, dims, compute_first_index(parent, indexes), compute_stride1(parent, indexes)) + # Compute the stride and offset + stride1 = compute_stride1(parent, indexes) + SubArray{eltype(P), N, P, I, true}(parent, indexes, dims, compute_offset1(parent, stride1, indexes), stride1) end check_parent_index_match{T,N}(parent::AbstractArray{T,N}, indexes::NTuple{N}) = nothing @@ -102,7 +103,7 @@ end function unsafe_slice{T,N}(A::AbstractArray{T,N}, I::Vararg{ViewIndex,N}) @_inline_meta J = to_indexes(I...) - SubArray(A, J, index_shape(A, J...)) + SubArray(A, J, map(dimlength, index_shape(A, J...))) end keep_leading_scalars(T::Tuple{}) = T @@ -179,20 +180,12 @@ function getindex{T,N}(V::SlowSubArray{T,N}, I::Vararg{Real,N}) @inbounds r = V.parent[reindex(V, V.indexes, to_indexes(I...))...] r end -# Explicitly define scalar linear indexing -- this is needed so the nonscalar -# indexing methods don't take precedence here -function getindex(V::SlowSubArray, i::Real) - @_inline_meta - @boundscheck checkbounds(V, i) - @inbounds r = V.parent[reindex(V, V.indexes, ind2sub(size(V), to_index(i)))...] - r -end typealias FastSubArray{T,N,P,I} SubArray{T,N,P,I,true} function getindex(V::FastSubArray, i::Real) @_inline_meta @boundscheck checkbounds(V, i) - @inbounds r = V.parent[V.first_index + V.stride1*(to_index(i)-1)] + @inbounds r = V.parent[V.offset1 + V.stride1*to_index(i)] r end # We can avoid a multiplication if the first parent index is a Colon or UnitRange @@ -200,24 +193,10 @@ typealias FastContiguousSubArray{T,N,P,I<:Tuple{Union{Colon, UnitRange}, Vararg{ function getindex(V::FastContiguousSubArray, i::Real) @_inline_meta @boundscheck checkbounds(V, i) - @inbounds r = V.parent[V.first_index + to_index(i)-1] - r -end -# Just like the slow case, explicitly define scalar indexing at N dims, too -function getindex{T,N}(V::FastSubArray{T,N}, I::Vararg{Real,N}) - @_inline_meta - @boundscheck checkbounds(V, I...) - @inbounds r = getindex(V, sub2ind(size(V), to_indexes(I...)...)) + @inbounds r = V.parent[V.offset1 + to_index(i)] r end -# Nonscalar indexing just copies a view -getindex{T,N}(V::SubArray{T,N}, i::ViewIndex) = (@_propagate_inbounds_meta; copy(slice(V, i))) -getindex{T,N}(V::SubArray{T,N}, I::Vararg{ViewIndex,N}) = (@_propagate_inbounds_meta; copy(slice(V, I...))) - -# Setindex is similar, but since we don't specially define non-scalar methods -# we only need to define the canonical methods. We don't need to worry about -# e.g., linear indexing for SlowSubArray since the fallbacks can do their thing function setindex!{T,N}(V::SlowSubArray{T,N}, x, I::Vararg{Real,N}) @_inline_meta @boundscheck checkbounds(V, I...) @@ -227,16 +206,15 @@ end function setindex!(V::FastSubArray, x, i::Real) @_inline_meta @boundscheck checkbounds(V, i) - @inbounds V.parent[V.first_index + V.stride1*(to_index(i)-1)] = x + @inbounds V.parent[V.offset1 + V.stride1*to_index(i)] = x V end function setindex!(V::FastContiguousSubArray, x, i::Real) @_inline_meta @boundscheck checkbounds(V, i) - @inbounds V.parent[V.first_index + to_index(i)-1] = x + @inbounds V.parent[V.offset1 + to_index(i)] = x V end -# Nonscalar setindex! falls back to the defaults function unsafe_slice{T,N}(V::SubArray{T,N}, I::Vararg{ViewIndex,N}) @_inline_meta @@ -257,7 +235,6 @@ getindex(::Colon, i) = to_index(i) unsafe_getindex(::Colon, i) = to_index(i) step(::Colon) = 1 -first(::Colon) = 1 isempty(::Colon) = false in(::Integer, ::Colon) = true @@ -289,42 +266,59 @@ iscontiguous{S<:SubArray}(::Type{S}) = false iscontiguous{F<:FastContiguousSubArray}(::Type{F}) = true # Fast linear SubArrays have their first index cached -first_index(V::FastSubArray) = V.first_index +first_index(V::FastSubArray) = V.offset1 + V.stride1 first_index(V::SubArray) = first_index(V.parent, V.indexes) function first_index(P::AbstractArray, indexes::Tuple) - f = 1 + f = first(linearindices(P)) s = 1 for i = 1:length(indexes) - f += (first(indexes[i])-1)*s + f += (_first(indexes[i], P, i)-first(indices(P, i)))*s s *= size(P, i) end f end +_first(::Colon, P, ::Colon) = first(linearindices(P)) +_first(i, P, ::Colon) = first(i) +_first(::Colon, P, d) = first(indices(P, d)) +_first(i, P, d) = first(i) # Computing the first index simply steps through the indices, accumulating the # sum of index each multiplied by the parent's stride. # The running sum is `f`; the cumulative stride product is `s`. -compute_first_index(parent, I::Tuple) = compute_first_index(1, 1, parent, 1, I) -compute_first_index(f, s, parent, dim, I::Tuple{Real, Vararg{Any}}) = - (@_inline_meta; compute_first_index(f + (I[1]-1)*s, s*size(parent, dim), parent, dim+1, tail(I))) -compute_first_index(f, s, parent, dim, I::Tuple{NoSlice, Vararg{Any}}) = - (@_inline_meta; compute_first_index(f + (I[1].i-1)*s, s*size(parent, dim), parent, dim+1, tail(I))) +# If the result is one-dimensional and it's a Colon, then linear +# indexing uses the indices along the given dimension. Otherwise +# linear indexing always starts with 1. +compute_offset1(parent, stride1::Integer, I::Tuple) = (@_inline_meta; compute_offset1(parent, stride1, find_extended_dims(I)..., I)) +compute_offset1(parent, stride1::Integer, dims::Tuple{Int}, inds::Tuple{Colon}, I::Tuple) = compute_linindex(parent, I) - stride1*first(indices(parent, dims[1])) # index-preserving case +compute_offset1(parent, stride1::Integer, dims, inds, I::Tuple) = compute_linindex(parent, I) - stride1 # linear indexing starts with 1 + +compute_linindex(parent, I) = compute_linindex(1, 1, parent, 1, I) +compute_linindex(f, s, parent, dim, I::Tuple{Real, Vararg{Any}}) = + (@_inline_meta; compute_linindex(f + (I[1]-first(indices(parent,dim)))*s, s*size(parent, dim), parent, dim+1, tail(I))) +compute_linindex(f, s, parent, dim, I::Tuple{NoSlice, Vararg{Any}}) = + (@_inline_meta; compute_linindex(f + (I[1].i-first(indices(parent,dim)))*s, s*size(parent, dim), parent, dim+1, tail(I))) # Just splat out the cartesian indices and continue -compute_first_index(f, s, parent, dim, I::Tuple{AbstractCartesianIndex, Vararg{Any}}) = - (@_inline_meta; compute_first_index(f, s, parent, dim, (I[1].I..., tail(I)...))) -compute_first_index(f, s, parent, dim, I::Tuple{Colon, Vararg{Any}}) = - (@_inline_meta; compute_first_index(f, s*size(parent, dim), parent, dim+1, tail(I))) -compute_first_index(f, s, parent, dim, I::Tuple{Any, Vararg{Any}}) = - (@_inline_meta; compute_first_index(f + (first(I[1])-1)*s, s*size(parent, dim), parent, dim+1, tail(I))) -compute_first_index(f, s, parent, dim, I::Tuple{}) = f - +compute_linindex(f, s, parent, dim, I::Tuple{AbstractCartesianIndex, Vararg{Any}}) = + (@_inline_meta; compute_linindex(f, s, parent, dim, (I[1].I..., tail(I)...))) +compute_linindex(f, s, parent, dim, I::Tuple{Colon, Vararg{Any}}) = + (@_inline_meta; compute_linindex(f, s*size(parent, dim), parent, dim+1, tail(I))) +compute_linindex(f, s, parent, dim, I::Tuple{Any, Vararg{Any}}) = + (@_inline_meta; compute_linindex(f + (first(I[1])-first(indices(parent,dim)))*s, s*size(parent, dim), parent, dim+1, tail(I))) +compute_linindex(f, s, parent, dim, I::Tuple{}) = f + +find_extended_dims(I) = (@_inline_meta; _find_extended_dims((), (), 1, I...)) +_find_extended_dims(dims, inds, dim) = dims, inds +_find_extended_dims(dims, inds, dim, ::Real, I...) = _find_extended_dims(dims, inds, dim+1, I...) +_find_extended_dims(dims, inds, dim, ::NoSlice, I...) = _find_extended_dims(dims, inds, dim+1, I...) +_find_extended_dims(dims, inds, dim, i1::AbstractCartesianIndex, I...) = _find_extended_dims(dims, inds, dim, i1.I..., I...) +_find_extended_dims(dims, inds, dim, i1, I...) = _find_extended_dims((dims..., dim), (inds..., i1), dim+1, I...) unsafe_convert{T,N,P,I<:Tuple{Vararg{Union{RangeIndex, NoSlice}}}}(::Type{Ptr{T}}, V::SubArray{T,N,P,I}) = unsafe_convert(Ptr{T}, V.parent) + (first_index(V)-1)*sizeof(T) -pointer(V::FastSubArray, i::Int) = pointer(V.parent, V.first_index + V.stride1*(i-1)) -pointer(V::FastContiguousSubArray, i::Int) = pointer(V.parent, V.first_index + i-1) -pointer(V::SubArray, i::Int) = pointer(V, ind2sub(size(V), i)) +pointer(V::FastSubArray, i::Int) = pointer(V.parent, V.offset1 + V.stride1*i) +pointer(V::FastContiguousSubArray, i::Int) = pointer(V.parent, V.offset1 + i) +pointer(V::SubArray, i::Int) = pointer(V, smart_ind2sub(shape(V), i)) function pointer{T,N,P<:Array,I<:Tuple{Vararg{Union{RangeIndex, NoSlice}}}}(V::SubArray{T,N,P,I}, is::Tuple{Vararg{Int}}) index = first_index(V) @@ -335,6 +329,38 @@ function pointer{T,N,P<:Array,I<:Tuple{Vararg{Union{RangeIndex, NoSlice}}}}(V::S return pointer(V.parent, index) end +# indices of the parent are preserved for ::Colon indices, otherwise +# they are taken from the range/vector +# Since bounds-checking is performance-critical and uses +# indices, it's worth optimizing these implementations thoroughly +indices(S::SubArray, d::Integer) = 1 <= d <= ndims(S) ? indices(S)[d] : (d > ndims(S) ? (1:1) : error("dimension $d out of range")) +indices(S::SubArray) = (@_inline_meta; _indices(indicesbehavior(parent(S)), S)) +_indices(::IndicesStartAt1, S::SubArray) = (@_inline_meta; map(s->1:s, size(S))) +_indices(::IndicesBehavior, S::SubArray) = (@_inline_meta; _indices((), 1, S, S.indexes...)) +_indices(out::Tuple, dim, S::SubArray) = out +_indices(out::Tuple, dim, S::SubArray, i1, I...) = (@_inline_meta; _indices((out..., 1:length(i1)), dim+1, S, I...)) +_indices(out::Tuple, dim, S::SubArray, ::Real, I...) = (@_inline_meta; _indices(out, dim+1, S, I...)) +_indices(out::Tuple, dim, S::SubArray, ::Colon, I...) = (@_inline_meta; _indices((out..., indices(parent(S), dim)), dim+1, S, I...)) +indices1{T}(S::SubArray{T,0}) = 1:1 +indices1(S::SubArray) = (@_inline_meta; _indices1(indicesbehavior(parent(S)), S)) +_indices1(::IndicesStartAt1, S::SubArray) = 1:S.dims[1] +_indices1(::IndicesBehavior, S::SubArray) = (@_inline_meta; _indices1(S, 1, S.indexes...)) +_indices1(S::SubArray, dim, i1, I...) = (@_inline_meta; 1:length(i1)) +_indices1(S::SubArray, dim, i1::Real, I...) = (@_inline_meta; _indices1(S, dim+1, I...)) +_indices1(S::SubArray, dim, i1::Colon, I...) = (@_inline_meta; indices(parent(S), dim)) + +# Moreover, incides(S) is fast but indices(S, d) is slower +indicesperformance{T<:SubArray}(::Type{T}) = IndicesSlow1D() + +indicesbehavior(S::SubArray) = _indicesbehavior(indicesbehavior(parent(S)), S) +_indicesbehavior(::IndicesStartAt1, S::SubArray) = IndicesStartAt1() +_indicesbehavior(::IndicesBehavior, S::SubArray) = (@_inline_meta; _indicesbehavior(keepcolon((), S.indexes...), S)) +keepcolon(out) = out +keepcolon(out, ::Colon, I...) = (@_inline_meta; (Colon(),)) +keepcolon(out, i1, I...) = (@_inline_meta; keepcolon(out, I...)) +_indicesbehavior(::Tuple{}, S::SubArray) = IndicesStartAt1() +_indicesbehavior(::Tuple{Colon}, S::SubArray) = indicesbehavior(parent(S)) + ## Compatability # deprecate? function parentdims(s::SubArray) diff --git a/doc/manual/arrays.rst b/doc/manual/arrays.rst index 8dea6274d3a6f..67b95d890c14c 100644 --- a/doc/manual/arrays.rst +++ b/doc/manual/arrays.rst @@ -50,7 +50,9 @@ Function Description :func:`length(A) ` the number of elements in ``A`` :func:`ndims(A) ` the number of dimensions of ``A`` :func:`size(A) ` a tuple containing the dimensions of ``A`` -:func:`size(A,n) ` the size of ``A`` in a particular dimension +:func:`size(A,n) ` the size of ``A`` along a particular dimension +:func:`indices(A) ` a tuple containing the valid indices of ``A`` +:func:`indices(A,n) ` a range expressing the valid indices along dimension ``n`` :func:`eachindex(A) ` an efficient iterator for visiting each position in ``A`` :func:`stride(A,k) ` the stride (linear index distance between adjacent elements) along dimension ``k`` :func:`strides(A) ` a tuple of the strides in each dimension diff --git a/doc/manual/interfaces.rst b/doc/manual/interfaces.rst index a12e1e2ccb3f2..5ede56eb24f8c 100644 --- a/doc/manual/interfaces.rst +++ b/doc/manual/interfaces.rst @@ -147,25 +147,29 @@ While this is starting to support more of the :ref:`indexing operations supporte Abstract Arrays --------------- -============================================================================ ============================================ ======================================================================================= -Methods to implement Brief description -============================================================================ ============================================ ======================================================================================= -:func:`size(A) ` Returns a tuple containing the dimensions of A -:func:`Base.linearindexing{T\<:YourType}(::Type{T}) ` Returns either ``Base.LinearFast()`` or ``Base.LinearSlow()``. See the description below. -:func:`getindex(A, i::Int) ` (if ``LinearFast``) Linear scalar indexing -:func:`getindex(A, I::Vararg{Int, N}) ` (if ``LinearSlow``, where ``N = ndims(A)``) N-dimensional scalar indexing -:func:`setindex!(A, v, i::Int) ` (if ``LinearFast``) Scalar indexed assignment -:func:`setindex!(A, v, I::Vararg{Int, N}) ` (if ``LinearSlow``, where ``N = ndims(A)``) N-dimensional scalar indexed assignment -**Optional methods** **Default definition** **Brief description** -:func:`getindex(A, I...) ` defined in terms of scalar :func:`getindex` :ref:`Multidimensional and nonscalar indexing ` -:func:`setindex!(A, I...) ` defined in terms of scalar :func:`setindex!` :ref:`Multidimensional and nonscalar indexed assignment ` -:func:`start`/:func:`next`/:func:`done` defined in terms of scalar :func:`getindex` Iteration -:func:`length(A) ` ``prod(size(A))`` Number of elements -:func:`similar(A) ` ``similar(A, eltype(A), size(A))`` Return a mutable array with the same shape and element type -:func:`similar(A, ::Type{S}) ` ``similar(A, S, size(A))`` Return a mutable array with the same shape and the specified element type -:func:`similar(A, dims::NTuple{Int}) ` ``similar(A, eltype(A), dims)`` Return a mutable array with the same element type and the specified dimensions -:func:`similar(A, ::Type{S}, dims::NTuple{Int}) ` ``Array{S}(dims)`` Return a mutable array with the specified element type and dimensions -============================================================================ ============================================ ======================================================================================= +===================================================================== ============================================ ======================================================================================= +Methods to implement Brief description +===================================================================== ============================================ ======================================================================================= +:func:`size(A) ` Returns a tuple containing the dimensions of ``A`` +:func:`getindex(A, i::Int) ` (if ``LinearFast``) Linear scalar indexing +:func:`getindex(A, I::Vararg{Int, N}) ` (if ``LinearSlow``, where ``N = ndims(A)``) N-dimensional scalar indexing +:func:`setindex!(A, v, i::Int) ` (if ``LinearFast``) Scalar indexed assignment +:func:`setindex!(A, v, I::Vararg{Int, N}) ` (if ``LinearSlow``, where ``N = ndims(A)``) N-dimensional scalar indexed assignment +**Optional methods** **Default definition** **Brief description** +:func:`Base.linearindexing(::Type) ` ``Base.LinearSlow()`` Returns either ``Base.LinearFast()`` or ``Base.LinearSlow()``. See the description below. +:func:`getindex(A, I...) ` defined in terms of scalar :func:`getindex` :ref:`Multidimensional and nonscalar indexing ` +:func:`setindex!(A, I...) ` defined in terms of scalar :func:`setindex!` :ref:`Multidimensional and nonscalar indexed assignment ` +:func:`start`/:func:`next`/:func:`done` defined in terms of scalar :func:`getindex` Iteration +:func:`length(A) ` ``prod(size(A))`` Number of elements +:func:`similar(A) ` ``similar(A, eltype(A), size(A))`` Return a mutable array with the same shape and element type +:func:`similar(A, ::Type{S}) ` ``similar(A, S, size(A))`` Return a mutable array with the same shape and the specified element type +:func:`similar(A, dims::NTuple{Int}) ` ``similar(A, eltype(A), dims)`` Return a mutable array with the same element type and size `dims` +:func:`similar(A, ::Type{S}, dims::NTuple{Int}) ` ``Array{S}(dims)`` Return a mutable array with the specified element type and size +**Non-traditional indices** **Default definition** **Brief description** +:func:`Base.indicesbehavior(::Type) ` ``Base.IndicesStartAt1()`` Trait with values ``IndicesStartAt1()``, ``IndicesUnitRange()``, ``IndicesList()`` +:func:`indices(A, d) ` ``1:size(A, d)`` Return the range of valid indices along dimension ``d`` +:func:`Base.similar(A, ::Type{S}, inds::NTuple{Ind}) ` ``similar(A, S, map(Base.dimlength, inds))`` Return a mutable array with the specified indices ``inds`` (see below for discussion of ``Ind``) +===================================================================== ============================================ ======================================================================================= If a type is defined as a subtype of ``AbstractArray``, it inherits a very large set of rich behaviors including iteration and multidimensional indexing built on top of single-element access. See the :ref:`arrays manual page ` and :ref:`standard library section ` for more supported methods. @@ -261,7 +265,7 @@ The result of indexing an ``AbstractArray`` can itself be an array (for instance 1.0 4.0 7.0 2.0 5.0 8.0 -In this example it is accomplished by defining ``Base.similar{T}(A::SparseArray, ::Type{T}, dims::Dims)`` to create the appropriate wrapped array. For this to work it's important that ``SparseArray`` is mutable (supports ``setindex!``). :func:`similar` is also used to allocate result arrays for arithmetic on ``AbstractArrays``, for instance: +In this example it is accomplished by defining ``Base.similar{T}(A::SparseArray, ::Type{T}, dims::Dims)`` to create the appropriate wrapped array. (Note that while ``similar`` supports 1- and 2-argument forms, in most case you only need to specialize the 3-argument form.) For this to work it's important that ``SparseArray`` is mutable (supports ``setindex!``). :func:`similar` is also used to allocate result arrays for arithmetic on ``AbstractArrays``, for instance: .. doctest:: @@ -283,3 +287,12 @@ In addition to all the iterable and indexable methods from above, these types ca julia> dot(A[:,1],A[:,2]) 32.0 + +If you are defining an array type that allows non-traditional indexing +(indices that start at something other than 1), you should specialize +``indices`` and ``indicesbehavior``. You should also specialize +``similar`` so that the ``dims`` argument (ordinarily a ``Dims`` +size-tuple) can be a mixture of ``Integer`` and ``UnitRange`` objects; +the ``Integer`` entries imply that the indexing starts from 1, whereas +the dimensions encoded with ``UnitRange`` may have arbitrary starting +index. diff --git a/doc/stdlib/arrays.rst b/doc/stdlib/arrays.rst index e931fccf07134..a17da81e6ba61 100644 --- a/doc/stdlib/arrays.rst +++ b/doc/stdlib/arrays.rst @@ -31,6 +31,30 @@ Basic functions julia> size(A,3,2) (4,3) +.. function:: indices(A) + + .. Docstring generated from Julia source + + Returns the tuple of valid indices for array ``A``\ . + +.. function:: indices(A, d) + + .. Docstring generated from Julia source + + Returns the valid range of indices for array ``A`` along dimension ``d``\ . + +.. function:: shape(A) + + .. Docstring generated from Julia source + + Returns a tuple specifying the "shape" of array ``A``\ . For arrays with conventional indexing (indices start at 1), this is equivalent to ``size(A)``\ ; otherwise it is equivalent to ``incides(A)``\ . + +.. function:: shape(A, d) + + .. Docstring generated from Julia source + + Specifies the "shape" of the array ``A`` along dimension ``d``\ . For arrays with conventional indexing (starting at 1), this is equivalent to ``size(A, d)``\ ; for arrays with unconventional indexing (indexing may start at something different from 1), it is equivalent to ``indices(A, d)``\ . + .. function:: length(A) -> Integer .. Docstring generated from Julia source @@ -72,6 +96,14 @@ Basic functions If you supply more than one ``AbstractArray`` argument, ``eachindex`` will create an iterable object that is fast for all arguments (a ``UnitRange`` if all inputs have fast linear indexing, a CartesianRange otherwise). If the arrays have different sizes and/or dimensionalities, ``eachindex`` returns an iterable that spans the largest range along each dimension. +.. function:: linearindices(A) + + .. Docstring generated from Julia source + + Returns a ``UnitRange`` specifying the valid range of indices for ``A[i]`` where ``i`` is an ``Int``\ . For arrays with conventional indexing (indices start at 1), or any multidimensional array, this is ``1:length(A)``\ ; however, for one-dimensional arrays with unconventional indices, this is ``indices(A, 1)``\ . + + Calling this function is the "safe" way to write algorithms that exploit linear indexing. + .. function:: Base.linearindexing(A) .. Docstring generated from Julia source @@ -251,6 +283,32 @@ Constructors 2.18425e-314 2.18425e-314 2.18425e-314 2.18425e-314 2.18425e-314 2.18425e-314 2.18425e-314 2.18425e-314 + See also ``allocate_for``\ . + +.. function:: allocate_for(storagetype, referencearray, [shape]) + + .. Docstring generated from Julia source + + Create an uninitialized mutable array analogous to that specified by ``storagetype``\ , but with type and shape specified by the final two arguments. The main purpose of this function is to support allocation of arrays that may have unconventional indexing (starting at other than 1), as determined by ``referencearray`` and the optional ``shape`` information. + + .. code-block:: julia + + **Examples**: + + allocate_for(Array{Int}, A) + + creates an array that "acts like" an ``Array{Int}`` (and might indeed be backed by one), but which is indexed identically to ``A``\ . If ``A`` has conventional indexing, this will likely just call ``Array{Int}(size(A))``\ , but if ``A`` has unconventional indexing then the indices of the result will match ``A``\ . + + .. code-block:: julia + + allocate_for(BitArray, A, (shape(A, 2),)) + + would create a 1-dimensional logical array whose indices match those of the columns of ``A``\ . + + The main purpose of the ``referencearray`` argument is to select a particular array type supporting unconventional indexing (as it is possible that several different ones will be simultaneously in use). + + See also ``similar``\ . + .. function:: reinterpret(type, A) .. Docstring generated from Julia source @@ -568,13 +626,19 @@ Indexing, Assignment, and Concatenation .. Docstring generated from Julia source - Throw an error if the specified indexes are not in bounds for the given array. Subtypes of ``AbstractArray`` should specialize this method if they need to provide custom bounds checking behaviors. + Throw an error if the specified ``indexes`` are not in bounds for the given ``array``\ . + +.. function:: checkbounds(Bool, array, indexes...) + + .. Docstring generated from Julia source + + Return ``true`` if the specified ``indexes`` are in bounds for the given ``array``\ . Subtypes of ``AbstractArray`` should specialize this method if they need to provide custom bounds checking behaviors. -.. function:: checkbounds(::Type{Bool}, dimlength::Integer, index) +.. function:: checkindex(Bool, inds::UnitRange, index) .. Docstring generated from Julia source - Return a ``Bool`` describing if the given index is within the bounds of the given dimension length. Custom types that would like to behave as indices for all arrays can extend this method in order to provide a specialized bounds checking implementation. + Return ``true`` if the given ``index`` is within the bounds of ``inds``\ . Custom types that would like to behave as indices for all arrays can extend this method in order to provide a specialized bounds checking implementation. .. function:: randsubseq(A, p) -> Vector diff --git a/test/abstractarray.jl b/test/abstractarray.jl index 5ec49feec9fc8..9a60572374a6e 100644 --- a/test/abstractarray.jl +++ b/test/abstractarray.jl @@ -1,5 +1,84 @@ # This file is a part of Julia. License is MIT: http://julialang.org/license +# Bounds checking +A = rand(3,3,3) +@test checkbounds(Bool, A, 1, 1, 1) == true +@test checkbounds(Bool, A, 1, 3, 3) == true +@test checkbounds(Bool, A, 1, 4, 3) == false +@test checkbounds(Bool, A, 1, -1, 3) == false +@test checkbounds(Bool, A, 1, 9) == true # partial linear indexing +@test checkbounds(Bool, A, 1, 10) == false # partial linear indexing +@test checkbounds(Bool, A, 1) == true +@test checkbounds(Bool, A, 27) == true +@test checkbounds(Bool, A, 28) == false +@test checkbounds(Bool, A, 2, 2, 2, 1) == true +@test checkbounds(Bool, A, 2, 2, 2, 2) == false + +@test Base.checkbounds_indices((1:5, 1:5), (2,2)) +@test !Base.checkbounds_indices((1:5, 1:5), (7,2)) +@test !Base.checkbounds_indices((1:5, 1:5), (2,0)) +@test Base.checkbounds_indices((1:5, 1:5), (13,)) +@test !Base.checkbounds_indices((1:5, 1:5), (26,)) +@test Base.checkbounds_indices((1:5, 1:5), (2,2,1)) +@test !Base.checkbounds_indices((1:5, 1:5), (2,2,2)) + +# sub2ind & ind2sub +# 0-dimensional +for i = 1:4 + @test sub2ind((), i) == i +end +@test sub2ind((), 2, 2) == 3 +@test ind2sub((), 1) == () +@test_throws BoundsError ind2sub((), 2) +# 1-dimensional +for i = 1:4 + @test sub2ind((3,), i) == i + @test ind2sub((3,), i) == (i,) +end +@test sub2ind((3,), 2, 2) == 5 +@test_throws MethodError ind2sub((3,), 2, 2) +# ambiguity btw cartesian indexing and linear indexing in 1d when +# indices may be nontraditional +@test_throws ArgumentError sub2ind((1:3,), 2) +@test_throws ArgumentError ind2sub((1:3,), 2) +# 2-dimensional +k = 0 +for j = 1:3, i = 1:4 + @test sub2ind((4,3), i, j) == (k+=1) + @test ind2sub((4,3), k) == (i,j) + @test sub2ind((1:4,1:3), i, j) == k + @test ind2sub((1:4,1:3), k) == (i,j) + @test sub2ind((0:3,3:5), i-1, j+2) == k + @test ind2sub((0:3,3:5), k) == (i-1, j+2) +end +# Delete when partial linear indexing is deprecated (#14770) +@test sub2ind((4,3), 7) == 7 +@test sub2ind((1:4,1:3), 7) == 7 +@test sub2ind((0:3,3:5), 7) == 8 +# 3-dimensional +l = 0 +for k = 1:2, j = 1:3, i = 1:4 + @test sub2ind((4,3,2), i, j, k) == (l+=1) + @test ind2sub((4,3,2), l) == (i,j,k) + @test sub2ind((1:4,1:3,1:2), i, j, k) == l + @test ind2sub((1:4,1:3,1:2), l) == (i,j,k) + @test sub2ind((0:3,3:5,-101:-100), i-1, j+2, k-102) == l + @test ind2sub((0:3,3:5,-101:-100), l) == (i-1, j+2, k-102) +end + +A = reshape(collect(1:9), (3,3)) +@test ind2sub(size(A), 6) == (3,2) +@test sub2ind(size(A), 3, 2) == 6 +@test ind2sub(A, 6) == (3,2) +@test sub2ind(A, 3, 2) == 6 + +# PR #9256 +function pr9256() + m = [1 2 3; 4 5 6; 7 8 9] + ind2sub(m, 6) +end +@test pr9256() == (3,2) + # token type on which to dispatch testing methods in order to avoid potential # name conflicts elsewhere in the base test suite type TestAbstractArray end @@ -254,12 +333,13 @@ end function test_in_bounds(::Type{TestAbstractArray}) n = rand(2:5) - dims = tuple(rand(2:5, n)...) - len = prod(dims) + sz = rand(2:5, n) + len = prod(sz) + A = zeros(sz...) for i in 1:len - @test checkbounds(Bool, dims, i) == true + @test checkbounds(Bool, A, i) == true end - @test checkbounds(Bool, dims, len + 1) == false + @test checkbounds(Bool, A, len + 1) == false end type UnimplementedFastArray{T, N} <: AbstractArray{T, N} end @@ -525,3 +605,9 @@ A = TSlowNIndexes(rand(2,2)) #16381 @inferred size(rand(3,2,1), 2, 1) @inferred size(rand(3,2,1), 2, 1, 3) + +@test @inferred(indices(rand(3,2))) == (1:3,1:2) +@test @inferred(indices(rand(3,2,1))) == (1:3,1:2,1:1) +@test @inferred(indices(rand(3,2), 1)) == 1:3 +@test @inferred(indices(rand(3,2), 2)) == 1:2 +@test @inferred(indices(rand(3,2), 3)) == 1:1 diff --git a/test/arrayops.jl b/test/arrayops.jl index 712d591bb7cb0..a556260ad942a 100644 --- a/test/arrayops.jl +++ b/test/arrayops.jl @@ -1104,18 +1104,6 @@ function i7197() ind2sub(size(S), 5) end @test i7197() == (2,2) -A = reshape(collect(1:9), (3,3)) -@test ind2sub(size(A), 6) == (3,2) -@test sub2ind(size(A), 3, 2) == 6 -@test ind2sub(A, 6) == (3,2) -@test sub2ind(A, 3, 2) == 6 - -# PR #9256 -function pr9256() - m = [1 2 3; 4 5 6; 7 8 9] - ind2sub(m, 6) -end -@test pr9256() == (3,2) # PR #8622 and general indexin test function pr8622() diff --git a/test/broadcast.jl b/test/broadcast.jl index f20ce3cb977c4..357711fbd585d 100644 --- a/test/broadcast.jl +++ b/test/broadcast.jl @@ -1,5 +1,58 @@ # This file is a part of Julia. License is MIT: http://julialang.org/license +module TestBroadcastInternals + +using Base.Broadcast: broadcast_shape, check_broadcast_shape, newindex, _bcs, _bcsm +using Base.Test + +@test @inferred(_bcs((), (3,5), (3,5))) == (3,5) +@test @inferred(_bcs((), (3,1), (3,5))) == (3,5) +@test @inferred(_bcs((), (3,), (3,5))) == (3,5) +@test @inferred(_bcs((), (3,5), (3,))) == (3,5) +@test_throws DimensionMismatch _bcs((), (3,5), (4,5)) +@test_throws DimensionMismatch _bcs((), (3,5), (3,4)) +@test @inferred(_bcs((), (-1:1, 2:5), (-1:1, 2:5))) == (-1:1, 2:5) +@test @inferred(_bcs((), (-1:1, 2:5), (1, 2:5))) == (-1:1, 2:5) +@test @inferred(_bcs((), (-1:1, 1), (1, 2:5))) == (-1:1, 2:5) +@test @inferred(_bcs((), (-1:1,), (-1:1, 2:5))) == (-1:1, 2:5) +@test_throws DimensionMismatch _bcs((), (-1:1, 2:6), (-1:1, 2:5)) +@test_throws DimensionMismatch _bcs((), (-1:1, 2:5), (2, 2:5)) + +@test @inferred(broadcast_shape(zeros(3,4), zeros(3,4))) == (3,4) +@test @inferred(broadcast_shape(zeros(3,4), zeros(3))) == (3,4) +@test @inferred(broadcast_shape(zeros(3), zeros(3,4))) == (3,4) +@test @inferred(broadcast_shape(zeros(3), zeros(1,4), zeros(1))) == (3,4) + +check_broadcast_shape((3,5), zeros(3,5)) +check_broadcast_shape((3,5), zeros(3,1)) +check_broadcast_shape((3,5), zeros(3)) +check_broadcast_shape((3,5), zeros(3,5), zeros(3)) +check_broadcast_shape((3,5), zeros(3,5), 1) +check_broadcast_shape((3,5), 5, 2) +@test_throws DimensionMismatch check_broadcast_shape((3,5), zeros(2,5)) +@test_throws DimensionMismatch check_broadcast_shape((3,5), zeros(3,4)) +@test_throws DimensionMismatch check_broadcast_shape((3,5), zeros(3,4,2)) +@test_throws DimensionMismatch check_broadcast_shape((3,5), zeros(3,5), zeros(2)) + +check_broadcast_shape((-1:1, 6:9), (-1:1, 6:9)) +check_broadcast_shape((-1:1, 6:9), (-1:1, 1)) +check_broadcast_shape((-1:1, 6:9), (1, 6:9)) +@test_throws DimensionMismatch check_broadcast_shape((-1:1, 6:9), (-1, 6:9)) +@test_throws DimensionMismatch check_broadcast_shape((-1:1, 6:9), (-1:1, 6)) +check_broadcast_shape((-1:1, 6:9), 1) +check_broadcast_shape((-1:1, 6:9), zeros(1,1)) + +ci(x) = CartesianIndex(x) +@test @inferred(newindex(ci((2,2)), (true, true))) == ci((2,2)) +@test @inferred(newindex(ci((2,2)), (true, false))) == ci((2,1)) +@test @inferred(newindex(ci((2,2)), (false, true))) == ci((1,2)) +@test @inferred(newindex(ci((2,2)), (false, false))) == ci((1,1)) +@test @inferred(newindex(ci((2,2)), (true,))) == ci((2,)) +@test @inferred(newindex(ci((2,2)), (false,))) == ci((1,)) +@test @inferred(newindex(ci((2,2)), ())) == 1 + +end + function as_sub(x::AbstractVector) y = similar(x, eltype(x), tuple(([size(x)...]*2)...)) y = sub(y, 2:2:length(y)) diff --git a/test/choosetests.jl b/test/choosetests.jl index fc9258e881cde..6a9fb13463df7 100644 --- a/test/choosetests.jl +++ b/test/choosetests.jl @@ -33,7 +33,7 @@ function choosetests(choices = []) "markdown", "base64", "serialize", "misc", "threads", "enums", "cmdlineargs", "i18n", "workspace", "libdl", "int", "checked", "intset", "floatfuncs", "compile", "parallel", "inline", - "boundscheck", "error", "ambiguous" + "boundscheck", "error", "ambiguous", "offsetarray" ] if Base.USE_GPL_LIBS diff --git a/test/copy.jl b/test/copy.jl index 0d6b2c0fecf2a..7be5a705d5f57 100644 --- a/test/copy.jl +++ b/test/copy.jl @@ -37,7 +37,7 @@ for (dest, src, bigsrc, emptysrc, res) in [ @test_throws exc copy!(dest, 1, x, idx, 1) end - @test_throws BoundsError copy!(dest, 1, src(), 1, -1) + @test_throws ArgumentError copy!(dest, 1, src(), 1, -1) @test_throws BoundsError copy!(dest, bigsrc()) diff --git a/test/offsetarray.jl b/test/offsetarray.jl new file mode 100644 index 0000000000000..d7bee902c22ac --- /dev/null +++ b/test/offsetarray.jl @@ -0,0 +1,359 @@ +# This file is a part of Julia. License is MIT: http://julialang.org/license + +# OffsetArrays (arrays with indexing that doesn't start at 1) + +# This test file is designed to exercise support for generic indexing, +# even though offset arrays aren't implemented in Base. + +module OAs + +using Base: SimIdx, Indices, LinearSlow, LinearFast + +export OffsetArray + +immutable OffsetArray{T,N,AA<:AbstractArray} <: AbstractArray{T,N} + parent::AA + offsets::NTuple{N,Int} +end +typealias OffsetVector{T,AA<:AbstractArray} OffsetArray{T,1,AA} + +OffsetArray{T,N}(A::AbstractArray{T,N}, offsets::NTuple{N,Int}) = OffsetArray{T,N,typeof(A)}(A, offsets) +OffsetArray{T,N}(A::AbstractArray{T,N}, offsets::Vararg{Int,N}) = OffsetArray(A, offsets) + +(::Type{OffsetArray{T,N}}){T,N}(inds::Indices{N}) = OffsetArray{T,N,Array{T,N}}(Array{T,N}(map(Base.dimlength, inds)), map(indsoffset, inds)) +(::Type{OffsetArray{T}}){T,N}(inds::Indices{N}) = OffsetArray{T,N}(inds) + +Base.linearindexing{T<:OffsetArray}(::Type{T}) = Base.linearindexing(parenttype(T)) +Base.indicesbehavior{T<:OffsetArray}(::Type{T}) = Base.IndicesUnitRange() +parenttype{T,N,AA}(::Type{OffsetArray{T,N,AA}}) = AA +parenttype(A::OffsetArray) = parenttype(typeof(A)) + +Base.parent(A::OffsetArray) = A.parent +Base.size(A::OffsetArray) = size(parent(A)) +Base.eachindex(::LinearSlow, A::OffsetArray) = CartesianRange(indices(A)) +Base.eachindex(::LinearFast, A::OffsetVector) = indices(A, 1) +Base.summary(A::OffsetArray) = string(typeof(A))*" with indices "*string(indices(A)) + +# Implementations of indices and indices1. Since bounds-checking is +# performance-critical and relies on indices, these are usually worth +# optimizing thoroughly. +@inline Base.indices(A::OffsetArray, d) = 1 <= d <= length(A.offsets) ? (o = A.offsets[d]; (1+o:size(parent(A),d)+o)) : (1:1) +@inline Base.indices(A::OffsetArray) = _indices((), A) # would rather use ntuple, but see #15276 +@inline _indices{T,N}(out::NTuple{N}, A::OffsetArray{T,N}) = out +@inline _indices(out, A::OffsetArray) = (d = length(out)+1; o = A.offsets[d]; _indices((out..., (1+o:size(parent(A),d)+o)), A)) +# By optimizing indices1 we can avoid a branch on the dim-check +Base.indices1{T}(A::OffsetArray{T,0}) = 1:1 +@inline Base.indices1(A::OffsetArray) = (o = A.offsets[1]; 1+o:size(parent(A),1)+o) + +function Base.similar(A::OffsetArray, T::Type, dims::Dims) + B = similar(parent(A), T, dims) +end +function Base.similar(A::AbstractArray, T::Type, inds::Tuple{Vararg{SimIdx}}) + B = similar(A, T, map(Base.dimlength, inds)) + OffsetArray(B, map(indsoffset, inds)) +end + +Base.allocate_for(f, A::OffsetArray, shape::SimIdx) = OffsetArray(f(Base.dimlength(shape)), (indsoffset(shape),)) +Base.allocate_for(f, A::OffsetArray, shape::Tuple{Vararg{SimIdx}}) = OffsetArray(f(map(Base.dimlength, shape)), map(indsoffset, shape)) +Base.promote_indices(a::OffsetArray, b::OffsetArray) = a + +Base.reshape(A::AbstractArray, inds::Indices) = OffsetArray(reshape(A, map(Base.dimlength, inds)), map(indsoffset, inds)) + +@inline function Base.getindex{T,N}(A::OffsetArray{T,N}, I::Vararg{Int,N}) + @boundscheck checkbounds(A, I...) + @inbounds ret = parent(A)[offset(A.offsets, I)...] + ret +end +@inline function Base._getindex(::LinearFast, A::OffsetVector, i::Int) + @boundscheck checkbounds(A, i) + @inbounds ret = parent(A)[offset(A.offsets, (i,))[1]] + ret +end +@inline function Base._getindex(::LinearFast, A::OffsetArray, i::Int) + @boundscheck checkbounds(A, i) + @inbounds ret = parent(A)[i] + ret +end +@inline function Base.setindex!{T,N}(A::OffsetArray{T,N}, val, I::Vararg{Int,N}) + @boundscheck checkbounds(A, I...) + @inbounds parent(A)[offset(A.offsets, I)...] = val + val +end +@inline function Base._setindex!(::LinearFast, A::OffsetVector, val, i::Int) + @boundscheck checkbounds(A, i) + @inbounds parent(A)[offset(A.offsets, (i,))[1]] = val + val +end +@inline function Base._setindex!(::LinearFast, A::OffsetArray, val, i::Int) + @boundscheck checkbounds(A, i) + @inbounds parent(A)[i] = val + val +end + +# Computing a shifted index (subtracting the offset) +offset{N}(offsets::NTuple{N,Int}, inds::NTuple{N,Int}) = _offset((), offsets, inds) +_offset(out, ::Tuple{}, ::Tuple{}) = out +@inline _offset(out, offsets, inds) = _offset((out..., inds[1]-offsets[1]), Base.tail(offsets), Base.tail(inds)) + +indsoffset(r::Range) = first(r) - 1 +indsoffset(i::Integer) = 0 + +end + +using OAs + +# Basics +A0 = [1 3; 2 4] +A = OffsetArray(A0, (-1,2)) # LinearFast +S = OffsetArray(slice(A0, 1:2, 1:2), (-1,2)) # LinearSlow +@test indices(A) == indices(S) == (0:1, 3:4) +@test A[0,3] == A[1] == S[0,3] == S[1] == 1 +@test A[1,3] == A[2] == S[1,3] == S[2] == 2 +@test A[0,4] == A[3] == S[0,4] == S[3] == 3 +@test A[1,4] == A[4] == S[1,4] == S[4] == 4 +@test_throws BoundsError A[1,1] +@test_throws BoundsError S[1,1] + +# Vector indexing +@test A[:, 3] == S[:, 3] == OffsetArray([1,2], (A.offsets[1],)) +@test A[:, 4] == S[:, 4] == OffsetArray([3,4], (A.offsets[1],)) +@test_throws BoundsError A[:, 1] +@test_throws BoundsError S[:, 1] +@test A[0, :] == S[0, :] == OffsetArray([1,3], (A.offsets[2],)) +@test A[1, :] == S[1, :] == OffsetArray([2,4], (A.offsets[2],)) +@test_throws BoundsError A[2, :] +@test_throws BoundsError S[2, :] +@test A[0:1, 3] == S[0:1, 3] == [1,2] +@test A[[1,0], 3] == S[[1,0], 3] == [2,1] +@test A[0, 3:4] == S[0, 3:4] == [1,3] +@test A[1, [4,3]] == S[1, [4,3]] == [4,2] +@test A[:, :] == S[:, :] == A + +# CartesianIndexing +@test A[CartesianIndex((0,3))] == S[CartesianIndex((0,3))] == 1 +@test_throws BoundsError A[CartesianIndex(1,1)] +@test_throws BoundsError S[CartesianIndex(1,1)] +@test eachindex(A) == 1:4 +@test eachindex(S) == CartesianRange((0:1,3:4)) + +# slice +S = slice(A, :, 3) +@test S == OffsetArray([1,2], (A.offsets[1],)) +@test S[0] == 1 +@test S[1] == 2 +@test_throws BoundsError S[2] +S = slice(A, 0, :) +@test S == OffsetArray([1,3], (A.offsets[2],)) +@test S[3] == 1 +@test S[4] == 3 +@test_throws BoundsError S[1] +S = slice(A, 0:0, 4) +@test S == [3] +@test S[1] == 3 +@test_throws BoundsError S[0] +S = slice(A, 1, 3:4) +@test S == [2,4] +@test S[1] == 2 +@test S[2] == 4 +@test_throws BoundsError S[3] +S = slice(A, :, :) +@test S == A +@test S[0,3] == S[1] == 1 +@test S[1,3] == S[2] == 2 +@test S[0,4] == S[3] == 3 +@test S[1,4] == S[4] == 4 +@test_throws BoundsError S[1,1] + +# iteration +for (a,d) in zip(A, A0) + @test a == d +end + + +# show +io = IOBuffer() +show(io, A) +str = takebuf_string(io) +@test str == "[1 3; 2 4]" +show(io, S) +str = takebuf_string(io) +@test str == "[1 3; 2 4]" +show(IOContext(io, multiline=true), A) +strs = split(strip(takebuf_string(io)), '\n') +@test strs[2] == " 1 3" +@test strs[3] == " 2 4" +v = OffsetArray(rand(3), (-2,)) +show(io, v) +str = takebuf_string(io) +show(io, parent(v)) +@test str == takebuf_string(io) +function cmp_showf(printfunc, io, A) + ioc = IOContext(io, limit=true, multiline=true, compact=true) + printfunc(ioc, A) + str1 = takebuf_string(io) + printfunc(ioc, parent(A)) + str2 = takebuf_string(io) + @test str1 == str2 +end +cmp_showf(Base.print_matrix, io, OffsetArray(rand(5,5), (10,-9))) # rows&cols fit +cmp_showf(Base.print_matrix, io, OffsetArray(rand(10^3,5), (10,-9))) # columns fit +cmp_showf(Base.print_matrix, io, OffsetArray(rand(5,10^3), (10,-9))) # rows fit +cmp_showf(Base.print_matrix, io, OffsetArray(rand(10^3,10^3), (10,-9))) # neither fits + +# Similar +B = similar(A, Float32) +@test isa(B, OffsetArray{Float32,2}) +@test size(B) == size(A) +@test indices(B) == indices(A) +B = similar(A, (3,4)) +@test isa(B, Array{Int,2}) +@test size(B) == (3,4) +@test indices(B) == (1:3, 1:4) +B = similar(A, (-3:3,4)) +@test isa(B, OffsetArray{Int,2}) +@test indices(B) == (-3:3, 1:4) +B = similar(parent(A), (-3:3,4)) +@test isa(B, OffsetArray{Int,2}) +@test indices(B) == (-3:3, 1:4) + +# Indexing with OffsetArray indices +i1 = OffsetArray([2,1], (-5,)) +i1 = OffsetArray([2,1], -5) +b = A0[i1, 1] +@test indices(b) == (-4:-3,) +@test b[-4] == 2 +@test b[-3] == 1 +b = A0[1,i1] +@test indices(b) == (-4:-3,) +@test b[-4] == 3 +@test b[-3] == 1 + +# logical indexing +@test A[A .> 2] == [3,4] + +# copy! +a = OffsetArray{Int}((-3:-1,)) +fill!(a, -1) +copy!(a, (1,2)) # non-array iterables +@test a[-3] == 1 +@test a[-2] == 2 +@test a[-1] == -1 +fill!(a, -1) +copy!(a, -2, (1,2)) +@test a[-3] == -1 +@test a[-2] == 1 +@test a[-1] == 2 +@test_throws BoundsError copy!(a, 1, (1,2)) +fill!(a, -1) +copy!(a, -2, (1,2,3), 2) +@test a[-3] == -1 +@test a[-2] == 2 +@test a[-1] == 3 +@test_throws BoundsError copy!(a, -2, (1,2,3), 1) +fill!(a, -1) +copy!(a, -2, (1,2,3), 1, 2) +@test a[-3] == -1 +@test a[-2] == 1 +@test a[-1] == 2 + +b = 1:2 # copy between AbstractArrays +bo = OffsetArray(1:2, (-3,)) +@test_throws BoundsError copy!(a, b) +fill!(a, -1) +copy!(a, bo) +@test a[-3] == -1 +@test a[-2] == 1 +@test a[-1] == 2 +fill!(a, -1) +copy!(a, -2, bo) +@test a[-3] == -1 +@test a[-2] == 1 +@test a[-1] == 2 +@test_throws BoundsError copy!(a, -4, bo) +@test_throws BoundsError copy!(a, -1, bo) +fill!(a, -1) +copy!(a, -3, b, 2) +@test a[-3] == 2 +@test a[-2] == a[-1] == -1 +@test_throws BoundsError copy!(a, -3, b, 1, 4) +am = OffsetArray{Int}((1:1, 7:9)) # for testing linear indexing +fill!(am, -1) +copy!(am, b) +@test am[1] == 1 +@test am[2] == 2 +@test am[3] == -1 +@test am[1,7] == 1 +@test am[1,8] == 2 +@test am[1,9] == -1 + +dest = similar(am) +map!(+, dest, am, am) +@test dest[1,7] == 2 +@test dest[1,8] == 4 +@test dest[1,9] == -2 + +A = OffsetArray(rand(4,4), (-3,5)) +@test maximum(A) == maximum(parent(A)) +@test minimum(A) == minimum(parent(A)) +@test extrema(A) == extrema(parent(A)) +C = similar(A) +cumsum!(C, A, 1) +@test parent(C) == cumsum(parent(A), 1) +@test parent(cumsum(A, 1)) == cumsum(parent(A), 1) +cumsum!(C, A, 2) +@test parent(C) == cumsum(parent(A), 2) +R = similar(A, (-2:-2, 6:9)) +maximum!(R, A) +@test parent(R) == maximum(parent(A), 1) +R = similar(A, (-2:1, 6:6)) +maximum!(R, A) +@test parent(R) == maximum(parent(A), 2) +amin, iamin = findmin(A) +pmin, ipmin = findmin(parent(A)) +@test amin == pmin +@test A[iamin] == amin +@test amin == parent(A)[ipmin] +amax, iamax = findmax(A) +pmax, ipmax = findmax(parent(A)) +@test amax == pmax +@test A[iamax] == amax +@test amax == parent(A)[ipmax] + +v = OffsetArray([1,1e100,1,-1e100], (-3,))*1000 +v2 = OffsetArray([1,-1e100,1,1e100], (5,))*1000 +@test isa(v, OffsetArray) +cv = OffsetArray([1,1e100,1e100,2], (-3,))*1000 +cv2 = OffsetArray([1,-1e100,-1e100,2], (5,))*1000 +@test isequal(cumsum_kbn(v), cv) +@test isequal(cumsum_kbn(v2), cv2) +@test isequal(sum_kbn(v), sum_kbn(parent(v))) + +io = IOBuffer() +writedlm(io, A) +seek(io, 0) +@test readdlm(io, eltype(A)) == parent(A) + +amin, amax = extrema(parent(A)) +@test clamp(A, (amax+amin)/2, amax) == clamp(parent(A), (amax+amin)/2, amax) + +@test unique(A, 1) == parent(A) +@test unique(A, 2) == parent(A) +v = OffsetArray(rand(8), (-2,)) +@test sort(v) == OffsetArray(sort(parent(v)), v.offsets) +@test sortrows(A) == OffsetArray(sortrows(parent(A)), A.offsets) +@test sortcols(A) == OffsetArray(sortcols(parent(A)), A.offsets) +@test sort(A, 1) == OffsetArray(sort(parent(A), 1), A.offsets) +@test sort(A, 2) == OffsetArray(sort(parent(A), 2), A.offsets) + +@test mapslices(v->sort(v), A, 1) == OffsetArray(mapslices(v->sort(v), parent(A), 1), A.offsets) +@test mapslices(v->sort(v), A, 2) == OffsetArray(mapslices(v->sort(v), parent(A), 2), A.offsets) + +@test rotl90(A) == OffsetArray(rotl90(parent(A)), A.offsets[[2,1]]) +@test rotr90(A) == OffsetArray(rotr90(parent(A)), A.offsets[[2,1]]) +@test flipdim(A, 1) == OffsetArray(flipdim(parent(A), 1), A.offsets) +@test flipdim(A, 2) == OffsetArray(flipdim(parent(A), 2), A.offsets) + +@test A+1 == OffsetArray(parent(A)+1, A.offsets) +@test 2*A == OffsetArray(2*parent(A), A.offsets) +@test A+A == OffsetArray(parent(A)+parent(A), A.offsets) +@test A.*A == OffsetArray(parent(A).*parent(A), A.offsets) diff --git a/test/serialize.jl b/test/serialize.jl index d0befa587d20a..d7a1fd996e3b8 100644 --- a/test/serialize.jl +++ b/test/serialize.jl @@ -243,6 +243,7 @@ module ArrayWrappers immutable ArrayWrapper{T,N,A<:AbstractArray} <: AbstractArray{T,N} data::A end +ArrayWrapper{T,N}(data::AbstractArray{T,N}) = ArrayWrapper{T,N,typeof(data)}(data) Base.size(A::ArrayWrapper) = size(A.data) Base.size(A::ArrayWrapper, d) = size(A.data, d) Base.getindex(A::ArrayWrapper, i::Real...) = getindex(A.data, i...) @@ -251,7 +252,7 @@ end let A = rand(3,4) for B in (sub(A, :, 2:4), slice(A, 2, 1:3)) - C = ArrayWrappers.ArrayWrapper{Float64,2,typeof(B)}(B) + C = ArrayWrappers.ArrayWrapper(B) io = IOBuffer() serialize(io, C) seek(io, 0)