Merge pull request #19015 from JuliaLang/teh/stats_indices

More non-1 indexing fixes

base/abstractarray.jl

@@ -193,13 +193,21 @@ function isassigned(a::AbstractArray, i::Int...)
 end
 
 # used to compute "end" for last index
-function trailingsize(A, n)
+function trailingsize(A::AbstractArray, n)
  s = 1
  for i=n:ndims(A)
  s *= size(A,i)
  end
  return s
 end
+function trailingsize(inds::Indices, n)
+ s = 1
+ for i=n:length(inds)
+ s *= unsafe_length(inds[i])
+ end
+ return s
+end
+# This version is type-stable even if inds is heterogeneous
 function trailingsize(inds::Indices)
  @_inline_meta
  prod(map(unsafe_length, inds))

base/deprecated.jl

@@ -1023,6 +1023,7 @@ end))
 
 @deprecate is (===)
 
+
 @deprecate_binding Filter Iterators.Filter
 @deprecate_binding Zip Iterators.Zip
 @deprecate filter(flt, itr) Iterators.filter(flt, itr)
@@ -1036,4 +1037,25 @@ end))
 # NOTE: Deprecation of Channel{T}() is implemented in channels.jl.
 # To be removed from there when 0.6 deprecations are removed.
 
+# Not exported, but probably better to have deprecations anyway
+function reduced_dims(::Tuple{}, d::Int)
+ d < 1 && throw(ArgumentError("dimension must be ≥ 1, got $d"))
+ ()
+end
+reduced_dims(::Tuple{}, region) = ()
+function reduced_dims(dims::Dims, region)
+ Base.depwarn("`reduced_dims` is deprecated for Dims-tuples; pass `indices` instead", :reduced_dims)
+ map(last, reduced_dims(map(n->OneTo(n), dims), region))
+end
+
+function reduced_dims0(::Tuple{}, d::Int)
+ d < 1 && throw(ArgumentError("dimension must be ≥ 1, got $d"))
+ ()
+end
+reduced_dims0(::Tuple{}, region) = ()
+function reduced_dims0(dims::Dims, region)
+ Base.depwarn("`reduced_dims0` is deprecated for Dims-tuples; pass `indices` instead", :reduced_dims0)
+ map(last, reduced_dims0(map(n->OneTo(n), dims), region))
+end
+
 # End deprecations scheduled for 0.6

base/multidimensional.jl

@@ -221,11 +221,11 @@ index_ndims() = ()
 # trailing indexes -- they behave as though they were never there for the
 # purposes of generalized linear indexing.
 typealias CI0 Union{CartesianIndex{0}, AbstractArray{CartesianIndex{0}}}
-index_lengths(A::AbstractArray, I::Colon) = (length(A),)
+index_lengths(A::AbstractArray, I::Colon) = (_length(A),)
 @inline index_lengths(A::AbstractArray, I...) = index_lengths_dim(A, 1, I...)
 index_lengths_dim(A, dim) = ()
-index_lengths_dim(A, dim, ::Colon) = (trailingsize(A, dim),)
-index_lengths_dim(A, dim, ::Colon, i::CI0, I::CI0...) = (trailingsize(A, dim), index_lengths_dim(A, dim+1, i, I...)...)
+index_lengths_dim(A, dim, ::Colon) = (trailingsize(indices(A), dim),)
+index_lengths_dim(A, dim, ::Colon, i::CI0, I::CI0...) = (trailingsize(indices(A), dim), index_lengths_dim(A, dim+1, i, I...)...)
 @inline index_lengths_dim(A, dim, ::Colon, i, I...) = (_length(indices(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_lengths_dim(A, dim+N, I...)...)

base/reducedim.jl

@@ -3,70 +3,58 @@
 ## Functions to compute the reduced shape
 
 # for reductions that expand 0 dims to 1
-reduced_dims(a::AbstractArray, region) = reduced_dims(size(a), region)
+reduced_dims(a::AbstractArray, region) = reduced_dims(indices(a), region)
 
 # for reductions that keep 0 dims as 0
-reduced_dims0(a::AbstractArray, region) = reduced_dims0(size(a), region)
-
-function reduced_dims{N}(siz::NTuple{N,Int}, d::Int, rd::Int)
- if d < 1
- throw(ArgumentError("dimension must be ≥ 1, got $d"))
- elseif d == 1
- return tuple(rd, siz[d+1:N]...)::typeof(siz)
- elseif 1 < d < N
- return tuple(siz[1:d-1]..., rd, siz[d+1:N]...)::typeof(siz)
- elseif d == N
- return tuple(siz[1:N-1]..., rd)::typeof(siz)
+reduced_dims0(a::AbstractArray, region) = reduced_dims0(indices(a), region)
+
+function reduced_dims{N}(inds::Indices{N}, d::Int, rd::AbstractUnitRange)
+ d < 1 && throw(ArgumentError("dimension must be ≥ 1, got $d"))
+ if d == 1
+ return (oftype(inds[1], rd), tail(inds)...)
+ elseif 1 < d <= N
+ return tuple(inds[1:d-1]..., oftype(inds[d], rd), inds[d+1:N]...)::typeof(inds)
  else
- return siz
+ return inds
  end
 end
-reduced_dims{N}(siz::NTuple{N,Int}, d::Int) = reduced_dims(siz, d, 1)
+reduced_dims{N}(inds::Indices{N}, d::Int) = reduced_dims(inds, d, OneTo(1))
 
-function reduced_dims0{N}(siz::NTuple{N,Int}, d::Int)
- if d < 1
- throw(ArgumentError("dimension must be ≥ 1, got $d"))
- elseif d <= N
- return reduced_dims(siz, d, (siz[d] == 0 ? 0 : 1))
+function reduced_dims0{N}(inds::Indices{N}, d::Int)
+ d < 1 && throw(ArgumentError("dimension must be ≥ 1, got $d"))
+ if d <= N
+ return reduced_dims(inds, d, (inds[d] == OneTo(0) ? OneTo(0) : OneTo(1)))
  else
- return siz
+ return inds
  end
 end
 
-function reduced_dims{N}(siz::NTuple{N,Int}, region)
- rsiz = [siz...]
+function reduced_dims{N}(inds::Indices{N}, region)
+ rinds = [inds...]
  for i in region
  isa(i, Integer) || throw(ArgumentError("reduced dimension(s) must be integers"))
- d = convert(Int, i)::Int
+ d = Int(i)
  if d < 1
  throw(ArgumentError("region dimension(s) must be ≥ 1, got $d"))
  elseif d <= N
- rsiz[d] = 1
+ rinds[d] = oftype(rinds[d], OneTo(1))
  end
  end
- tuple(rsiz...)::typeof(siz)
+ tuple(rinds...)::typeof(inds)
 end
 
-function reduced_dims0{N}(siz::NTuple{N,Int}, region)
- rsiz = [siz...]
+function reduced_dims0{N}(inds::Indices{N}, region)
+ rinds = [inds...]
  for i in region
  isa(i, Integer) || throw(ArgumentError("reduced dimension(s) must be integers"))
- d = convert(Int, i)::Int
+ d = Int(i)
  if d < 1
  throw(ArgumentError("region dimension(s) must be ≥ 1, got $d"))
  elseif d <= N
- rsiz[d] = (rsiz[d] == 0 ? 0 : 1)
+ rinds[d] = oftype(rinds[d], (rinds[d] == OneTo(0) ? OneTo(0) : OneTo(1)))
  end
  end
- tuple(rsiz...)::typeof(siz)
-end
-
-function regionsize(a, region)
- s = 1
- for d in region
- s *= size(a,d)
- end
- s
+ tuple(rinds...)::typeof(inds)
 end
 
 
@@ -388,10 +376,10 @@ For an array input, returns the value and index of the minimum over the given re
 function findmin{T}(A::AbstractArray{T}, region)
  if isempty(A)
  return (similar(A, reduced_dims0(A, region)),
- zeros(Int, reduced_dims0(A, region)))
+ similar(dims->zeros(Int, dims), reduced_dims0(A, region)))
  end
  return findminmax!(<, reducedim_initarray0(A, region, typemax(T)),
- zeros(Int, reduced_dims0(A, region)), A)
+ similar(dims->zeros(Int, dims), reduced_dims0(A, region)), A)
 end
 
 """
@@ -415,10 +403,10 @@ For an array input, returns the value and index of the maximum over the given re
 function findmax{T}(A::AbstractArray{T}, region)
  if isempty(A)
  return (similar(A, reduced_dims0(A,region)),
- zeros(Int, reduced_dims0(A,region)))
+ similar(dims->zeros(Int, dims), reduced_dims0(A,region)))
  end
  return findminmax!(>, reducedim_initarray0(A, region, typemin(T)),
- zeros(Int, reduced_dims0(A, region)), A)
+ similar(dims->zeros(Int, dims), reduced_dims0(A, region)), A)
 end
 
 reducedim1(R, A) = length(indices1(R)) == 1

base/sort.jl

@@ -2,7 +2,7 @@
 
 module Sort
 
-using Base: Order, Checked, copymutable, linearindices, linearindexing, viewindexing, LinearFast
+using Base: Order, Checked, copymutable, linearindices, linearindexing, viewindexing, LinearFast, _length
 
 import
  Base.sort,
@@ -66,7 +66,8 @@ issorted(itr;
  issorted(itr, ord(lt,by,rev,order))
 
 function select!(v::AbstractVector, k::Union{Int,OrdinalRange}, o::Ordering)
- sort!(v, 1, length(v), PartialQuickSort(k), o)
+ inds = indices(v, 1)
+ sort!(v, first(inds), last(inds), PartialQuickSort(k), o)
  v[k]
 end
 select!(v::AbstractVector, k::Union{Int,OrdinalRange};
@@ -187,7 +188,7 @@ searchsorted{T<:Real}(a::Range{T}, x::Real, o::DirectOrdering) =
 
 for s in [:searchsortedfirst, :searchsortedlast, :searchsorted]
  @eval begin
- $s(v::AbstractVector, x, o::Ordering) = $s(v,x,1,length(v),o)
+ $s(v::AbstractVector, x, o::Ordering) = (inds = indices(v, 1); $s(v,x,first(inds),last(inds),o))
  $s(v::AbstractVector, x;
  lt=isless, by=identity, rev::Bool=false, order::Ordering=Forward) =
  $s(v,x,ord(lt,by,rev,order))
@@ -432,7 +433,7 @@ function sort!(v::AbstractVector;
  order::Ordering=Forward)
  ordr = ord(lt,by,rev,order)
  if ordr === Forward && isa(v,Vector) && eltype(v)<:Integer
- n = length(v)
+ n = _length(v)
  if n > 1
  min, max = extrema(v)
  (diff, o1) = sub_with_overflow(max, min)
@@ -478,7 +479,7 @@ sort(v::AbstractVector; kws...) = sort!(copymutable(v); kws...)
 ## selectperm: the permutation to sort the first k elements of an array ##
 
 selectperm(v::AbstractVector, k::Union{Integer,OrdinalRange}; kwargs...) =
- selectperm!(Vector{eltype(k)}(length(v)), v, k; kwargs..., initialized=false)
+ selectperm!(similar(Vector{eltype(k)}, indices(v,1)), v, k; kwargs..., initialized=false)
 
 function selectperm!{I<:Integer}(ix::AbstractVector{I}, v::AbstractVector,
  k::Union{Int, OrdinalRange};
@@ -521,7 +522,7 @@ function sortperm(v::AbstractVector;
  order::Ordering=Forward)
  ordr = ord(lt,by,rev,order)
  if ordr === Forward && isa(v,Vector) && eltype(v)<:Integer
- n = length(v)
+ n = _length(v)
  if n > 1
  min, max = extrema(v)
  (diff, o1) = sub_with_overflow(max, min)
@@ -553,7 +554,7 @@ function sortperm!{I<:Integer}(x::AbstractVector{I}, v::AbstractVector;
  order::Ordering=Forward,
  initialized::Bool=false)
  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))"))
+ throw(ArgumentError("index vector must have the same indices as the source vector, $(indices(x,1)) != $(indices(v,1))"))
  end
  if !initialized
  @inbounds for i = indices(v,1)

base/sparse/sparsematrix.jl

@@ -1864,9 +1864,9 @@ end
 # and computing reductions along columns into SparseMatrixCSC is
 # non-trivial, so use Arrays for output
 Base.reducedim_initarray{R}(A::SparseMatrixCSC, region, v0, ::Type{R}) =
- fill!(Array{R}(Base.reduced_dims(A,region)), v0)
+ fill!(similar(dims->Array{R}(dims), Base.reduced_dims(A,region)), v0)
 Base.reducedim_initarray0{R}(A::SparseMatrixCSC, region, v0, ::Type{R}) =
- fill!(Array{R}(Base.reduced_dims0(A,region)), v0)
+ fill!(similar(dims->Array{R}(dims), Base.reduced_dims0(A,region)), v0)
 
 # General mapreduce
 function _mapreducezeros(f, op, T::Type, nzeros::Int, v0)

base/statistics.jl

@@ -24,12 +24,12 @@ function mean(f::Callable, iterable)
  return total/count
 end
 mean(iterable) = mean(identity, iterable)
-mean(f::Callable, A::AbstractArray) = sum(f, A) / length(A)
-mean(A::AbstractArray) = sum(A) / length(A)
+mean(f::Callable, A::AbstractArray) = sum(f, A) / _length(A)
+mean(A::AbstractArray) = sum(A) / _length(A)
 
 function mean!{T}(R::AbstractArray{T}, A::AbstractArray)
  sum!(R, A; init=true)
- scale!(R, length(R) / length(A))
+ scale!(R, _length(R) / _length(A))
  return R
 end
 
@@ -140,7 +140,7 @@ function centralize_sumabs2!{S,T,N}(R::AbstractArray{S}, A::AbstractArray{T,N},
 end
 
 function varm{T}(A::AbstractArray{T}, m::Number; corrected::Bool=true)
- n = length(A)
+ n = _length(A)
  n == 0 && return convert(real(momenttype(T)), NaN)
  n == 1 && return convert(real(momenttype(T)), abs2(A[1] - m)/(1 - Int(corrected)))
  return centralize_sumabs2(A, m) / (n - Int(corrected))
@@ -150,7 +150,7 @@ function varm!{S}(R::AbstractArray{S}, A::AbstractArray, m::AbstractArray; corre
  if isempty(A)
  fill!(R, convert(S, NaN))
  else
- rn = div(length(A), length(R)) - Int(corrected)
+ rn = div(_length(A), _length(R)) - Int(corrected)
  scale!(centralize_sumabs2!(R, A, m), convert(S, 1/rn))
  end
  return R
@@ -282,7 +282,7 @@ stdm(iterable, m::Number; corrected::Bool=true) =
 _conj{T<:Real}(x::AbstractArray{T}) = x
 _conj(x::AbstractArray) = conj(x)
 
-_getnobs(x::AbstractVector, vardim::Int) = length(x)
+_getnobs(x::AbstractVector, vardim::Int) = _length(x)
 _getnobs(x::AbstractMatrix, vardim::Int) = size(x, vardim)
 
 function _getnobs(x::AbstractVecOrMat, y::AbstractVecOrMat, vardim::Int)
@@ -309,11 +309,11 @@ unscaled_covzm(x::AbstractMatrix, y::AbstractMatrix, vardim::Int) =
 
 # covzm (with centered data)
 
-covzm(x::AbstractVector, corrected::Bool=true) = unscaled_covzm(x) / (length(x) - Int(corrected))
+covzm(x::AbstractVector, corrected::Bool=true) = unscaled_covzm(x) / (_length(x) - Int(corrected))
 covzm(x::AbstractMatrix, vardim::Int=1, corrected::Bool=true) =
  scale!(unscaled_covzm(x, vardim), inv(size(x,vardim) - Int(corrected)))
 covzm(x::AbstractVector, y::AbstractVector, corrected::Bool=true) =
- unscaled_covzm(x, y) / (length(x) - Int(corrected))
+ unscaled_covzm(x, y) / (_length(x) - Int(corrected))
 covzm(x::AbstractVecOrMat, y::AbstractVecOrMat, vardim::Int=1, corrected::Bool=true) =
  scale!(unscaled_covzm(x, y, vardim), inv(_getnobs(x, y, vardim) - Int(corrected)))
 
@@ -568,16 +568,18 @@ function median!{T}(v::AbstractVector{T})
  isnan(x) && return x
  end
  end
- n = length(v)
+ inds = indices(v, 1)
+ n = length(inds)
+ mid = div(first(inds)+last(inds),2)
  if isodd(n)
- return middle(select!(v,div(n+1,2)))
+ return middle(select!(v,mid))
  else
- m = select!(v, div(n,2):div(n,2)+1)
+ m = select!(v, mid:mid+1)
  return middle(m[1], m[2])
  end
 end
 median!{T}(v::AbstractArray{T}) = median!(vec(v))
-median{T}(v::AbstractArray{T}) = median!(copy!(Array{T,1}(length(v)), v))
+median{T}(v::AbstractArray{T}) = median!(copy!(Array{T,1}(_length(v)), v))
 
 """
  median(v[, region])

test/offsetarray.jl

@@ -42,6 +42,20 @@ S = OffsetArray(view(A0, 1:2, 1:2), (-1,2)) # LinearSlow
 @test A[1, [4,3]] == S[1, [4,3]] == [4,2]
 @test A[:, :] == S[:, :] == A
 
+A_3_3 = OffsetArray(Array{Int}(3,3), (-2,-1))
+A_3_3[:, :] = reshape(1:9, 3, 3)
+for i = 1:9 @test A_3_3[i] == i end
+A_3_3[-1:1, 0:2] = reshape(1:9, 3, 3)
+for i = 1:9 @test A_3_3[i] == i end
+A_3_3[:, :] = 1:9
+for i = 1:9 @test A_3_3[i] == i end
+A_3_3[-1:1, 0:2] = 1:9
+for i = 1:9 @test A_3_3[i] == i end
+A_3_3[:] = 1:9
+for i = 1:9 @test A_3_3[i] == i end
+A_3_3[1:9] = 1:9
+for i = 1:9 @test A_3_3[i] == i end
+
 # CartesianIndexing
 @test A[CartesianIndex((0,3))] == S[CartesianIndex((0,3))] == 1
 @test_throws BoundsError A[CartesianIndex(1,1)]
@@ -259,6 +273,9 @@ 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))
+@test maximum(A, 1) == OffsetArray(maximum(parent(A), 1), (0,A.offsets[2]))
+@test maximum(A, 2) == OffsetArray(maximum(parent(A), 2), (A.offsets[1],0))
+@test maximum(A, 1:2) == maximum(parent(A), 1:2)
 C = similar(A)
 cumsum!(C, A, 1)
 @test parent(C) == cumsum(parent(A), 1)
@@ -293,6 +310,13 @@ I,J,N = findnz(z)
 @test find(x->x>0, h) == [-1,1]
 @test find(x->x<0, h) == [-2,0]
 @test find(x->x==0, h) == [2]
+@test mean(A_3_3) == median(A_3_3) == 5
+@test mean(x->2x, A_3_3) == 10
+@test mean(A_3_3, 1) == median(A_3_3, 1) == OffsetArray([2 5 8], (0,A_3_3.offsets[2]))
+@test mean(A_3_3, 2) == median(A_3_3, 2) == OffsetArray([4,5,6]'', (A_3_3.offsets[1],0))
+@test var(A_3_3) == 7.5
+@test std(A_3_3, 1) == OffsetArray([1 1 1], (0,A_3_3.offsets[2]))
+@test std(A_3_3, 2) == OffsetArray([3,3,3]'', (A_3_3.offsets[1],0))
 
 @test_approx_eq vecnorm(v) vecnorm(parent(v))
 @test_approx_eq vecnorm(A) vecnorm(parent(A))