add accumulate, accumulate!

- rename cumop! to accumulate! and export it
- add specialized 1d method to accumulate!
- add accumulate
- deprecate cummin, cummax

Author: jw3126

base/abstractarray.jl

@@ -1407,7 +1407,7 @@ function typed_hvcat{T}(::Type{T}, rows::Tuple{Vararg{Int}}, as...)
     T[rs...;]
 end
 
-## Reductions and scans ##
+## Reductions and accumulates ##
 
 function isequal(A::AbstractArray, B::AbstractArray)
     if A === B return true end

base/arraymath.jl

@@ -448,7 +448,7 @@ ctranspose{T}(x::AbstractVector{T}) = T[ ctranspose(v) for i=of_indices(x, OneTo
 
 # see discussion in #18364 ... we try not to widen type of the resulting array
 # from cumsum or cumprod, but in some cases (+, Bool) we may not have a choice.
-rcum_promote_type{T<:Number}(op, ::Type{T}) = promote_op(op, T)
+rcum_promote_type{T<:Number}(op, ::Type{T}) = promote_op(op, T, T)
 rcum_promote_type{T}(op, ::Type{T}) = T
 
 # handle sums of Vector{Bool} and similar.   it would be nice to handle

base/deprecated.jl

@@ -144,6 +144,10 @@ end
 @deprecate chol(A::Number, ::Type{Val{:L}})         ctranspose(chol(A))
 @deprecate chol(A::AbstractMatrix, ::Type{Val{:L}}) ctranspose(chol(A))
 
+@deprecate cummin(A, dim=1) accumulate(min, A, dim=1)
+@deprecate cummax(A, dim=1) accumulate(max, A, dim=1)
+
+
 # Number updates
 
 # rem1 is inconsistent for x==0: The result should both have the same

base/exports.jl

@@ -498,11 +498,15 @@ export
     conj!,
     copy!,
     cummax,
+    cummax!,
     cummin,
+    cummin!,
     cumprod,
     cumprod!,
     cumsum,
     cumsum!,
+    accumulate,
+    accumulate!,
     cumsum_kbn,
     eachindex,
     extrema,

base/multidimensional.jl

@@ -479,51 +479,6 @@ end
     end
 end
 
-for (f, fmod, op) = ((:cummin, :_cummin!, :min), (:cummax, :_cummax!, :max))
-    @eval function ($f)(v::AbstractVector)
-        n = length(v)
-        cur_val = v[1]
-        res = similar(v, n)
-        res[1] = cur_val
-        for i in 2:n
-            cur_val = ($op)(v[i], cur_val)
-            res[i] = cur_val
-        end
-        return res
-    end
-
-    @eval function ($f)(A::AbstractArray, axis::Integer)
-        axis > 0 || throw(ArgumentError("axis must be a positive integer"))
-        res = similar(A)
-        axis > ndims(A) && return copy!(res, A)
-        inds = indices(A)
-        if isempty(inds[axis])
-            return res
-        end
-        R1 = CartesianRange(inds[1:axis-1])
-        R2 = CartesianRange(inds[axis+1:end])
-        ($fmod)(res, A, R1, R2, axis)
-    end
-
-    @eval @noinline function ($fmod)(res, A::AbstractArray, R1::CartesianRange, R2::CartesianRange, axis::Integer)
-        inds = indices(A, axis)
-        i1 = first(inds)
-        for I2 in R2
-            for I1 in R1
-                res[I1, i1, I2] = A[I1, i1, I2]
-            end
-            for i = i1+1:last(inds)
-                for I1 in R1
-                    res[I1, i, I2] = ($op)(A[I1, i, I2], res[I1, i-1, I2])
-                end
-            end
-        end
-        res
-    end
-
-    @eval ($f)(A::AbstractArray) = ($f)(A, 1)
-end
-
 """
     cumsum(A, dim=1)
 
@@ -548,8 +503,9 @@ julia> cumsum(a,2)
  4  9  15
 ```
 """
-cumsum{T}(A::AbstractArray{T}, axis::Integer=1) =  cumsum!(similar(A, Base.rcum_promote_type(+, T)), A, axis)
-cumsum!(B, A::AbstractArray) = cumsum!(B, A, 1)
+cumsum(A::AbstractArray, axis::Integer=1) =  accumulate(+, A, axis)
+cumsum!(B, A, axis::Integer) = accumulate!(+, B, A, axis)
+
 """
     cumprod(A, dim=1)
 
@@ -574,13 +530,64 @@ julia> cumprod(a,2)
  4  20  120
 ```
 """
-cumprod(A::AbstractArray, axis::Integer=1) = cumprod!(similar(A), A, axis)
-cumprod!(B, A) = cumprod!(B, A, 1)
+cumprod(A::AbstractArray, axis::Integer=1) = accumulate(*, A, axis)
+cumprod!(B, A, axis::Integer) = accumulate!(*, B, A, axis)
 
-cumsum!(B, A, axis::Integer) = cumop!(+, B, A, axis)
-cumprod!(B, A, axis::Integer) = cumop!(*, B, A, axis)
+"""
+    accumulate(op, A, dim=1)
+
+Cumulative operation `op` along a dimension `dim` (defaults to 1). See also
+[`accumulate!`](:func:`accumulate!`) 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`](:func:`cumsum`), [`cummin`](:func:`cummin`), [`cummax`](:func:`cummax`), [`cumprod`](:func:`cumprod`)
+
+```jldoctest
+julia> accumulate(+, [1,2,3])
+3-element Array{Int64,1}:
+ 1
+ 3
+ 6
+
+julia> accumulate(*, [1,2,3])
+3-element Array{Int64,1}:
+ 1
+ 2
+ 6
+
+julia> accumulate(min, [1,2,-1])
+3-element Array{Int64,1}:
+  1
+  1
+ -1
+```
+"""
+function accumulate(op, A, axis::Integer=1)
+    out = similar(A, Base.rcum_promote_type(op, eltype(A)))
+    accumulate!(op, out, A, axis)
+end
 
-function cumop!(op, B, A, axis::Integer)
+function accumulate!(op, B, A::AbstractVector, axis::Integer=1)
+    axis > 0 || throw(ArgumentError("axis must be a positive integer"))
+    size(B) == size(A) || throw(DimensionMismatch("shape of B must match A"))
+    isempty(A) && return B
+    axis > 1 && return copy!(B, A)
+    cur_val = A[1]
+    B[1] = cur_val
+    @inbounds for i in 2:length(A)
+        cur_val = op(A[i], cur_val)
+        B[i] = cur_val
+    end
+    return B
+end
+
+"""
+    accumulate!(op, B, A, dim=1)
+
+Cumulative operation `op` on A along a dimension, storing the result in B. The dimension defaults to 1.
+See also [`accumulate`](:func:`accumulate`).
+"""
+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"))
@@ -601,12 +608,12 @@ function cumop!(op, B, A, axis::Integer)
     else
         R1 = CartesianRange(indices(A)[1:axis-1])   # not type-stable
         R2 = CartesianRange(indices(A)[axis+1:end])
-        _cumop!(op, B, A, R1, inds_t[axis], R2) # use function barrier
+        _accumulate!(op, B, A, R1, inds_t[axis], R2) # use function barrier
     end
     return B
 end
 
-@noinline function _cumop!(op, B, A, R1, ind, R2)
+@noinline function _accumulate!(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

doc/stdlib/arrays.rst

@@ -1463,6 +1463,38 @@ Indexing, Assignment, and Concatenation
 Array functions
 ---------------
 
+.. function:: accumulate(op, A, dim=1)
+
+   .. Docstring generated from Julia source
+
+   Cumulative operation ``op`` along a dimension ``dim`` (defaults to 1). See also :func:`accumulate!` 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: :func:`cumsum`\ , :func:`cummin`\ , :func:`cummax`\ , :func:`cumprod`
+
+   .. doctest::
+
+       julia> accumulate(+, [1,2,3])
+       3-element Array{Int64,1}:
+        1
+        3
+        6
+
+       julia> accumulate(*, [1,2,3])
+       3-element Array{Int64,1}:
+        1
+        2
+        6
+
+       julia> accumulate(min, [1,2,-1])
+       3-element Array{Int64,1}:
+         1
+         1
+        -1
+
+.. function:: accumulate!(op, B, A, dim=1)
+
+   .. Docstring generated from Julia source
+
+   Cumulative operation ``op`` on A along a dimension, storing the result in B. The dimension defaults to 1. See also :func:`accumulate`\ .
+
 .. function:: cumprod(A, dim=1)
 
    .. Docstring generated from Julia source

test/arrayops.jl

@@ -1868,3 +1868,42 @@ end
         @test size(a) == size(b)
     end
 end
+
+@testset "accumulate, accumulate!" begin
+
+    @test accumulate(+, [1,2,3]) == [1, 3, 6]
+    @test accumulate(min, [1 2; 3 4], 1) == [1 2; 1 2]
+    @test accumulate(max, [1 2; 3 0], 2) == [1 2; 3 3]
+    @test accumulate(+, Bool[]) == Int[]
+    @test accumulate(*, Bool[]) == Bool[]
+    @test accumulate(+, Float64[]) == Float64[]
+
+    @test accumulate(min, [1, 2, 5, -1, 3, -2]) == [1, 1, 1, -1, -1, -2]
+    @test accumulate(max, [1, 2, 5, -1, 3, -2]) == [1, 2, 5, 5, 5, 5]
+
+    @test accumulate(max, [1 0; 0 1], 1) == [1 0; 1 1]
+    @test accumulate(max, [1 0; 0 1], 2) == [1 1; 0 1]
+    @test accumulate(min, [1 0; 0 1], 1) == [1 0; 0 0]
+    @test accumulate(min, [1 0; 0 1], 2) == [1 0; 0 0]
+
+    N = 5
+    for arr in [rand(Float64, N), rand(Bool, N), rand(-2:2, N)]
+        for (op, cumop) in [(+, cumsum), (*, cumprod)]
+            @inferred accumulate(op, arr)
+            accumulate_arr = accumulate(op, arr)
+            @test accumulate_arr ≈ cumop(arr)
+            @test accumulate_arr[end] ≈ reduce(op, arr)
+            @test accumulate_arr[1] ≈ arr[1]
+            @test accumulate(op, arr, 10) ≈ arr
+
+            if eltype(arr) in [Int, Float64] # eltype of out easy
+                out = similar(arr)
+                @test out ≉ accumulate_arr
+                @test accumulate!(op, out, arr) ≈ accumulate_arr
+                @test out ≈ accumulate_arr
+            end
+
+        end
+    end
+
+end

test/reduce.jl

@@ -268,13 +268,6 @@ let es = sum_kbn(z), es2 = sum_kbn(z[1:10^5])
     @test (es2 - cs[10^5]) < es2 * 1e-13
 end
 
-@test isequal(cummin([1, 2, 5, -1, 3, -2]), [1, 1, 1, -1, -1, -2])
-@test isequal(cummax([1, 2, 5, -1, 3, -2]), [1, 2, 5, 5, 5, 5])
-
-@test isequal(cummax([1 0; 0 1], 1), [1 0; 1 1])
-@test isequal(cummax([1 0; 0 1], 2), [1 1; 0 1])
-@test isequal(cummin([1 0; 0 1], 1), [1 0; 0 0])
-@test isequal(cummin([1 0; 0 1], 2), [1 0; 0 0])
 
 @test sum(collect(map(UInt8,0:255))) == 32640
 @test sum(collect(map(UInt8,254:255))) == 509