diff --git a/NEWS.md b/NEWS.md index d3a5986905af1..1cee9e536ba0f 100644 --- a/NEWS.md +++ b/NEWS.md @@ -663,6 +663,9 @@ Deprecated or removed * `a:b` is deprecated for constructing a `StepRange` when `a` and `b` have physical units (Dates and Times). Use `a:s:b`, where `s = Dates.Day(1)` or `s = Dates.Second(1)`. + * `cumsum`, `cumprod`, `accumulate`, and their mutating versions now require a `dim` + argument instead of defaulting to using the first dimension ([#24684]). + Command-line option changes --------------------------- diff --git a/base/abstractarraymath.jl b/base/abstractarraymath.jl index 5d8cfc5607ad8..352d438525a26 100644 --- a/base/abstractarraymath.jl +++ b/base/abstractarraymath.jl @@ -265,12 +265,12 @@ end # TODO: Needs a separate IndexCartesian method, this is only fast for IndexLinear """ - cumsum_kbn(A, [dim::Integer=1]) + cumsum_kbn(A, dim::Integer) Cumulative sum along a dimension, using the Kahan-Babuska-Neumaier compensated summation -algorithm for additional accuracy. The dimension defaults to 1. +algorithm for additional accuracy. """ -function cumsum_kbn(A::AbstractArray{T}, axis::Integer=1) where T<:AbstractFloat +function cumsum_kbn(A::AbstractArray{T}, axis::Integer) where T<:AbstractFloat dimsA = size(A) ndimsA = ndims(A) axis_size = dimsA[axis] diff --git a/base/deprecated.jl b/base/deprecated.jl index 88c604a38d923..18e9ee6d93d41 100644 --- a/base/deprecated.jl +++ b/base/deprecated.jl @@ -2105,6 +2105,9 @@ end @deprecate chol!(x::Number, uplo) chol(x) false end +@deprecate cumsum(A::AbstractArray) cumsum(A, 1) +@deprecate cumsum_kbn(A::AbstractArray) cumsum_kbn(A, 1) +@deprecate cumprod(A::AbstractArray) cumprod(A, 1) # issue #16307 @deprecate finalizer(o, f::Function) finalizer(f, o) diff --git a/base/multidimensional.jl b/base/multidimensional.jl index 8df892a1e3fe4..df0f9c70db4d2 100644 --- a/base/multidimensional.jl +++ b/base/multidimensional.jl @@ -675,25 +675,25 @@ function accumulate_pairwise(op, v::AbstractVector{T}) where T return accumulate_pairwise!(op, out, v) end -function cumsum!(out, v::AbstractVector, axis::Integer=1) +function cumsum!(out, v::AbstractVector, dim::Integer) # we dispatch on the possibility of numerical stability issues - _cumsum!(out, v, axis, TypeArithmetic(eltype(out))) + _cumsum!(out, v, dim, TypeArithmetic(eltype(out))) end -function _cumsum!(out, v, axis, ::ArithmeticRounds) - axis == 1 ? accumulate_pairwise!(+, out, v) : copy!(out, v) +function _cumsum!(out, v, dim, ::ArithmeticRounds) + dim == 1 ? accumulate_pairwise!(+, out, v) : copy!(out, v) end -function _cumsum!(out, v, axis, ::ArithmeticUnknown) - _cumsum!(out, v, axis, ArithmeticRounds()) +function _cumsum!(out, v, dim, ::ArithmeticUnknown) + _cumsum!(out, v, dim, ArithmeticRounds()) end -function _cumsum!(out, v, axis, ::TypeArithmetic) - axis == 1 ? accumulate!(+, out, v) : copy!(out, v) +function _cumsum!(out, v, dim, ::TypeArithmetic) + dim == 1 ? accumulate!(+, out, v) : copy!(out, v) end """ - cumsum(A, dim=1) + cumsum(A, dim::Integer) -Cumulative sum along a dimension `dim`. See also [`cumsum!`](@ref) +Cumulative sum along the dimension `dim`. See also [`cumsum!`](@ref) to use a preallocated output array, both for performance and to control the precision of the output (e.g. to avoid overflow). @@ -714,22 +714,52 @@ julia> cumsum(a,2) 4 9 15 ``` """ -function cumsum(A::AbstractArray{T}, axis::Integer=1) where T +function cumsum(A::AbstractArray{T}, dim::Integer) where T out = similar(A, rcum_promote_type(+, T)) - cumsum!(out, A, axis) + cumsum!(out, A, dim) end """ - cumsum!(B, A, dim::Integer=1) + cumsum(x::AbstractVector) -Cumulative sum of `A` along a dimension, storing the result in `B`. See also [`cumsum`](@ref). +Cumulative sum a vector. See also [`cumsum!`](@ref) +to use a preallocated output array, both for performance and to control the precision of the +output (e.g. to avoid overflow). + +```jldoctest +julia> cumsum([1, 1, 1]) +3-element Array{Int64,1}: + 1 + 2 + 3 + +julia> cumsum([fill(1, 2) for i in 1:3]) +3-element Array{Array{Int64,1},1}: + [1, 1] + [2, 2] + [3, 3] +``` +""" +cumsum(x::AbstractVector) = cumsum(x, 1) + +""" + cumsum!(B, A, dim::Integer) + +Cumulative sum of `A` along the dimension `dim`, storing the result in `B`. See also [`cumsum`](@ref). +""" +cumsum!(B, A, dim::Integer) = accumulate!(+, B, A, dim) + +""" + cumsum!(y::AbstractVector, x::AbstractVector) + +Cumulative sum of a vector `x`, storing the result in `y`. See also [`cumsum`](@ref). """ -cumsum!(B, A, axis::Integer=1) = accumulate!(+, B, A, axis) +cumsum!(y::AbstractVector, x::AbstractVector) = cumsum!(y, x, 1) """ - cumprod(A, dim=1) + cumprod(A, dim::Integer) -Cumulative product along a dimension `dim`. See also +Cumulative product along the dimension `dim`. See also [`cumprod!`](@ref) to use a preallocated output array, both for performance and to control the precision of the output (e.g. to avoid overflow). @@ -750,20 +780,79 @@ julia> cumprod(a,2) 4 20 120 ``` """ -cumprod(A::AbstractArray, axis::Integer=1) = accumulate(*, A, axis) +cumprod(A::AbstractArray, dim::Integer) = accumulate(*, A, dim) + +""" + cumprod(x::AbstractVector) + +Cumulative product of a vector. See also +[`cumprod!`](@ref) to use a preallocated output array, both for performance and +to control the precision of the output (e.g. to avoid overflow). + +```jldoctest +julia> cumprod(fill(1//2, 3)) +3-element Array{Rational{Int64},1}: + 1//2 + 1//4 + 1//8 + +julia> cumprod([fill(1//3, 2, 2) for i in 1:3]) +3-element Array{Array{Rational{Int64},2},1}: + Rational{Int64}[1//3 1//3; 1//3 1//3] + Rational{Int64}[2//9 2//9; 2//9 2//9] + Rational{Int64}[4//27 4//27; 4//27 4//27] +``` +""" +cumprod(x::AbstractVector) = cumprod(x, 1) """ - cumprod!(B, A, dim::Integer=1) + cumprod!(B, A, dim::Integer) -Cumulative product of `A` along a dimension, storing the result in `B`. +Cumulative product of `A` along the dimension `dim`, storing the result in `B`. See also [`cumprod`](@ref). """ -cumprod!(B, A, axis::Integer=1) = accumulate!(*, B, A, axis) +cumprod!(B, A, dim::Integer) = accumulate!(*, B, A, dim) """ - accumulate(op, A, dim=1) + cumprod!(y::AbstractVector, x::AbstractVector) -Cumulative operation `op` along a dimension `dim`. See also +Cumulative product of a vector `x`, storing the result in `y`. +See also [`cumprod`](@ref). +""" +cumprod!(y::AbstractVector, x::AbstractVector) = cumprod!(y, x, 1) + +""" + accumulate(op, A, dim::Integer) + +Cumulative operation `op` along the dimension `dim`. See also +[`accumulate!`](@ref) to use a preallocated output array, both for performance and +to control the precision of the output (e.g. to avoid overflow). For common operations +there are specialized variants of `accumulate`, see: +[`cumsum`](@ref), [`cumprod`](@ref) + +```jldoctest +julia> accumulate(+, fill(1, 3, 3), 1) +3×3 Array{Int64,2}: + 1 1 1 + 2 2 2 + 3 3 3 + +julia> accumulate(+, fill(1, 3, 3), 2) +3×3 Array{Int64,2}: + 1 2 3 + 1 2 3 + 1 2 3 +``` +""" +function accumulate(op, A, dim::Integer) + out = similar(A, rcum_promote_type(op, eltype(A))) + accumulate!(op, out, A, dim) +end + +""" + accumulate(op, x::AbstractVector) + +Cumulative operation `op` on a vector. See also [`accumulate!`](@ref) to use a preallocated output array, both for performance and to control the precision of the output (e.g. to avoid overflow). For common operations there are specialized variants of `accumulate`, see: @@ -783,14 +872,67 @@ julia> accumulate(*, [1,2,3]) 6 ``` """ -function accumulate(op, A, axis::Integer=1) - out = similar(A, rcum_promote_type(op, eltype(A))) - accumulate!(op, out, A, axis) +accumulate(op, x::AbstractVector) = accumulate(op, x, 1) + +""" + accumulate!(op, B, A, dim::Integer) + +Cumulative operation `op` on `A` along the dimension `dim`, storing the result in `B`. +See also [`accumulate`](@ref). +""" +function accumulate!(op, B, A, dim::Integer) + dim > 0 || throw(ArgumentError("dim must be a positive integer")) + inds_t = indices(A) + indices(B) == inds_t || throw(DimensionMismatch("shape of B must match A")) + dim > ndims(A) && return copy!(B, A) + isempty(inds_t[dim]) && return B + if dim == 1 + # We can accumulate to a temporary variable, which allows + # register usage and will be slightly faster + ind1 = inds_t[1] + @inbounds for I in CartesianRange(tail(inds_t)) + 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:dim-1]) # not type-stable + R2 = CartesianRange(indices(A)[dim+1:end]) + _accumulate!(op, B, A, R1, inds_t[dim], R2) # use function barrier + end + return B end +""" + accumulate!(op, y, x::AbstractVector) +Cumulative operation `op` on a vector `x`, storing the result in `y`. +See also [`accumulate`](@ref). """ - accumulate(op, v0, A) +function accumulate!(op::Op, y, x::AbstractVector) where Op + isempty(x) && return y + v1 = first(x) + _accumulate1!(op, y, v1, x, 1) +end + +@noinline function _accumulate!(op, B, A, R1, ind, R2) + # Copy the initial element in each 1d vector along dimension `dim` + ii = first(ind) + @inbounds for J in R2, I in R1 + B[I, ii, J] = A[I, ii, 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 + +""" + accumulate(op, v0, x::AbstractVector) Like `accumulate`, but using a starting element `v0`. The first entry of the result will be `op(v0, first(A))`. @@ -810,30 +952,23 @@ julia> accumulate(min, 0, [1,2,-1]) -1 ``` """ -function accumulate(op, v0, A, axis::Integer=1) - T = rcum_promote_type(op, typeof(v0), eltype(A)) - out = similar(A, T) - accumulate!(op, out, v0, A, 1) -end - -function accumulate!(op::Op, B, A::AbstractVector, axis::Integer=1) where Op - isempty(A) && return B - v1 = first(A) - _accumulate1!(op, B, v1, A, axis) +function accumulate(op, v0, x::AbstractVector) + T = rcum_promote_type(op, typeof(v0), eltype(x)) + out = similar(x, T) + accumulate!(op, out, v0, x) end -function accumulate!(op, B, v0, A::AbstractVector, axis::Integer=1) - isempty(A) && return B - v1 = op(v0, first(A)) - _accumulate1!(op, B, v1, A, axis) +function accumulate!(op, y, v0, x::AbstractVector) + isempty(x) && return y + v1 = op(v0, first(x)) + _accumulate1!(op, y, v1, x, 1) end - -function _accumulate1!(op, B, v1, A::AbstractVector, axis::Integer=1) - axis > 0 || throw(ArgumentError("axis must be a positive integer")) +function _accumulate1!(op, B, v1, A::AbstractVector, dim::Integer) + dim > 0 || throw(ArgumentError("dim must be a positive integer")) inds = linearindices(A) inds == linearindices(B) || throw(DimensionMismatch("linearindices of A and B don't match")) - axis > 1 && return copy!(B, A) + dim > 1 && return copy!(B, A) i1 = inds[1] cur_val = v1 B[i1] = cur_val @@ -844,51 +979,6 @@ function _accumulate1!(op, B, v1, A::AbstractVector, axis::Integer=1) return B end -""" - accumulate!(op, B, A, dim=1) - -Cumulative operation `op` on `A` along a dimension, storing the result in `B`. -See also [`accumulate`](@ref). -""" -function accumulate!(op, B, A, axis::Integer=1) - axis > 0 || throw(ArgumentError("axis must be a positive integer")) - inds_t = indices(A) - indices(B) == inds_t || throw(DimensionMismatch("shape of B must match A")) - axis > ndims(A) && return copy!(B, A) - isempty(inds_t[axis]) && return B - if axis == 1 - # We can accumulate to a temporary variable, which allows - # register usage and will be slightly faster - ind1 = inds_t[1] - @inbounds for I in CartesianRange(tail(inds_t)) - 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]) - _accumulate!(op, B, A, R1, inds_t[axis], R2) # use function barrier - end - return B -end - -@noinline function _accumulate!(op, B, A, R1, ind, R2) - # Copy the initial element in each 1d vector along dimension `axis` - ii = first(ind) - @inbounds for J in R2, I in R1 - B[I, ii, J] = A[I, ii, 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 """ diff --git a/test/arrayops.jl b/test/arrayops.jl index f81d51103335e..e37ed04117946 100644 --- a/test/arrayops.jl +++ b/test/arrayops.jl @@ -632,7 +632,6 @@ end A1 = reshape(repmat([1,2],1,12),2,3,4) A2 = reshape(repmat([1 2 3],2,4),2,3,4) A3 = reshape(repmat([1 2 3 4],6,1),2,3,4) - @test isequal(cumsum(A),A1) @test isequal(cumsum(A,1),A1) @test isequal(cumsum(A,2),A2) @test isequal(cumsum(A,3),A3) @@ -931,7 +930,7 @@ end end # issue #2342 -@test isequal(cumsum([1 2 3]), [1 2 3]) +@test isequal(cumsum([1 2 3], 1), [1 2 3]) @testset "set-like operations" begin @test isequal(union([1,2,3], [4,3,4]), [1,2,3,4])