From 360fd408d3590924178e74268180cc22ace440aa Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Tue, 18 Oct 2016 11:30:40 -0500 Subject: [PATCH 1/3] Some vector-indexing fixes for non-1 indices --- base/abstractarray.jl | 10 +++++++++- base/multidimensional.jl | 6 +++--- test/offsetarray.jl | 14 ++++++++++++++ 3 files changed, 26 insertions(+), 4 deletions(-) diff --git a/base/abstractarray.jl b/base/abstractarray.jl index abd6772288d96..6ec8f5cce19d2 100644 --- a/base/abstractarray.jl +++ b/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)) diff --git a/base/multidimensional.jl b/base/multidimensional.jl index 3c7e81728188c..79d444d0cfbb5 100644 --- a/base/multidimensional.jl +++ b/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...)...) diff --git a/test/offsetarray.jl b/test/offsetarray.jl index 458b6fd5ad429..514102f14d234 100644 --- a/test/offsetarray.jl +++ b/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)] From 28c9b2ba92d9a2440b4e6756fa164d8cffd1f036 Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Tue, 18 Oct 2016 11:31:17 -0500 Subject: [PATCH 2/3] Fix region size computation for non-1 indices --- base/deprecated.jl | 22 ++++++++++++ base/reducedim.jl | 72 ++++++++++++++++--------------------- base/sparse/sparsematrix.jl | 4 +-- test/offsetarray.jl | 3 ++ test/reducedim.jl | 4 +-- 5 files changed, 59 insertions(+), 46 deletions(-) diff --git a/base/deprecated.jl b/base/deprecated.jl index e5e1438572983..4165ed3574a3a 100644 --- a/base/deprecated.jl +++ b/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 diff --git a/base/reducedim.jl b/base/reducedim.jl index a6bc026ee0ae5..1d9f82593e54d 100644 --- a/base/reducedim.jl +++ b/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 diff --git a/base/sparse/sparsematrix.jl b/base/sparse/sparsematrix.jl index e1126c3d916e5..46080f5db4a15 100644 --- a/base/sparse/sparsematrix.jl +++ b/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) diff --git a/test/offsetarray.jl b/test/offsetarray.jl index 514102f14d234..2980185182e88 100644 --- a/test/offsetarray.jl +++ b/test/offsetarray.jl @@ -273,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) diff --git a/test/reducedim.jl b/test/reducedim.jl index dee4785878e05..5bbb2bbca6896 100644 --- a/test/reducedim.jl +++ b/test/reducedim.jl @@ -20,7 +20,7 @@ for region in Any[ 1, 2, 3, 4, 5, (1, 2), (1, 3), (1, 4), (2, 3), (2, 4), (3, 4), (1, 2, 3), (1, 3, 4), (2, 3, 4), (1, 2, 3, 4)] # println("region = $region") - r = fill(NaN, Base.reduced_dims(size(Areduc), region)) + r = fill(NaN, map(length, Base.reduced_dims(indices(Areduc), region))) @test sum!(r, Areduc) ≈ safe_sum(Areduc, region) @test prod!(r, Areduc) ≈ safe_prod(Areduc, region) @test maximum!(r, Areduc) ≈ safe_maximum(Areduc, region) @@ -62,7 +62,7 @@ end # Test reduction along first dimension; this is special-cased for # size(A, 1) >= 16 Breduc = rand(64, 3) -r = fill(NaN, Base.reduced_dims(size(Breduc), 1)) +r = fill(NaN, map(length, Base.reduced_dims(indices(Breduc), 1))) @test sum!(r, Breduc) ≈ safe_sum(Breduc, 1) @test sumabs!(r, Breduc) ≈ safe_sumabs(Breduc, 1) @test sumabs2!(r, Breduc) ≈ safe_sumabs2(Breduc, 1) From 750d93bab1cfe80c86cfde8c5bfb376bc6a791c2 Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Sun, 23 Oct 2016 14:50:54 -0500 Subject: [PATCH 3/3] Fix statistics functions for non-1 indices --- base/sort.jl | 15 ++++++++------- base/statistics.jl | 26 ++++++++++++++------------ test/offsetarray.jl | 7 +++++++ 3 files changed, 29 insertions(+), 19 deletions(-) diff --git a/base/sort.jl b/base/sort.jl index 5f1ef98158562..08a6e829bb36a 100644 --- a/base/sort.jl +++ b/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) diff --git a/base/statistics.jl b/base/statistics.jl index 05805f2a2db50..e3766ef5cecaa 100644 --- a/base/statistics.jl +++ b/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]) diff --git a/test/offsetarray.jl b/test/offsetarray.jl index 2980185182e88..88b4a9097d0cd 100644 --- a/test/offsetarray.jl +++ b/test/offsetarray.jl @@ -310,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))