Skip to content

Commit

Permalink
Deprecate (cum)sum_kbn
Browse files Browse the repository at this point in the history
This removes `sum_kbn` and `cumsum_kbn` in favor of `sum` and `cumsum`,
respectively, with the `*_kbn` functions moving to a package.

Fixes #24804
  • Loading branch information
ararslan committed Dec 7, 2017
1 parent 037b818 commit 1a4ebd4
Show file tree
Hide file tree
Showing 9 changed files with 20 additions and 150 deletions.
67 changes: 0 additions & 67 deletions base/abstractarraymath.jl
Original file line number Diff line number Diff line change
Expand Up @@ -238,73 +238,6 @@ function circshift(a::AbstractArray, shiftamt)
circshift!(similar(a), a, map(Integer, (shiftamt...,)))
end

# Uses K-B-N summation
function cumsum_kbn(v::AbstractVector{T}) where T<:AbstractFloat
r = similar(v)
if isempty(v); return r; end

inds = indices(v, 1)
i1 = first(inds)
s = r[i1] = v[i1]
c = zero(T)
for i=i1+1:last(inds)
vi = v[i]
t = s + vi
if abs(s) >= abs(vi)
c += ((s-t) + vi)
else
c += ((vi-t) + s)
end
s = t
r[i] = s+c
end
return r
end

# Uses K-B-N summation
# TODO: Needs a separate IndexCartesian method, this is only fast for IndexLinear

"""
cumsum_kbn(A, dim::Integer)
Cumulative sum along a dimension, using the Kahan-Babuska-Neumaier compensated summation
algorithm for additional accuracy.
"""
function cumsum_kbn(A::AbstractArray{T}, axis::Integer) where T<:AbstractFloat
dimsA = size(A)
ndimsA = ndims(A)
axis_size = dimsA[axis]
axis_stride = 1
for i = 1:(axis-1)
axis_stride *= size(A,i)
end

if axis_size <= 1
return A
end

B = similar(A)
C = similar(A)

for i = 1:length(A)
if div(i-1, axis_stride) % axis_size == 0
B[i] = A[i]
C[i] = zero(T)
else
s = B[i-axis_stride]
Ai = A[i]
B[i] = t = s + Ai
if abs(s) >= abs(Ai)
C[i] = C[i-axis_stride] + ((s-t) + Ai)
else
C[i] = C[i-axis_stride] + ((Ai-t) + s)
end
end
end

return B + C
end

## Other array functions ##

"""
Expand Down
5 changes: 4 additions & 1 deletion base/deprecated.jl
Original file line number Diff line number Diff line change
Expand Up @@ -2136,7 +2136,6 @@ end
@deprecate RowVector{T}(n::Tuple{Int,Int}) where {T} RowVector{T}(uninitialized, n)

@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
Expand Down Expand Up @@ -2183,6 +2182,10 @@ end
@deprecate SSHCredentials SSHCredential false
end

# issue #24804
@deprecate_moved sum_kbn "KahanSummation"
@deprecate_moved cumsum_kbn "KahanSummation"

# END 0.7 deprecations

# BEGIN 1.0 deprecations
Expand Down
2 changes: 0 additions & 2 deletions base/exports.jl
Original file line number Diff line number Diff line change
Expand Up @@ -433,7 +433,6 @@ export
cumsum!,
accumulate,
accumulate!,
cumsum_kbn,
eachindex,
extrema,
fill!,
Expand Down Expand Up @@ -527,7 +526,6 @@ export
sub2ind,
sum!,
sum,
sum_kbn,
to_indices,
vcat,
vec,
Expand Down
32 changes: 0 additions & 32 deletions base/reduce.jl
Original file line number Diff line number Diff line change
Expand Up @@ -391,38 +391,6 @@ julia> sum(1:20)
sum(a) = mapreduce(promote_sys_size_add, +, a)
sum(a::AbstractArray{Bool}) = count(a)


# Kahan (compensated) summation: O(1) error growth, at the expense
# of a considerable increase in computational expense.

"""
sum_kbn(A)
Returns the sum of all elements of `A`, using the Kahan-Babuska-Neumaier compensated
summation algorithm for additional accuracy.
"""
function sum_kbn(A)
T = @default_eltype(typeof(A))
c = promote_sys_size_add(zero(T)::T)
i = start(A)
if done(A, i)
return c
end
Ai, i = next(A, i)
s = Ai - c
while !(done(A, i))
Ai, i = next(A, i)
t = s + Ai
if abs(s) >= abs(Ai)
c -= ((s-t) + Ai)
else
c -= ((Ai-t) + s)
end
s = t
end
s - c
end

## prod
"""
prod(f, itr)
Expand Down
2 changes: 1 addition & 1 deletion doc/src/manual/style-guide.md
Original file line number Diff line number Diff line change
Expand Up @@ -159,7 +159,7 @@ uses (e.g. `a[i]::Int`) than to try to pack many alternatives into one type.
* functions are lowercase ([`maximum`](@ref), [`convert`](@ref)) and, when readable, with multiple
words squashed together ([`isequal`](@ref), [`haskey`](@ref)). When necessary, use underscores
as word separators. Underscores are also used to indicate a combination of concepts ([`remotecall_fetch`](@ref)
as a more efficient implementation of `fetch(remotecall(...))`) or as modifiers ([`sum_kbn`](@ref)).
as a more efficient implementation of `fetch(remotecall(...))`) or as modifiers.
* conciseness is valued, but avoid abbreviation ([`indexin`](@ref) rather than `indxin`) as
it becomes difficult to remember whether and how particular words are abbreviated.

Expand Down
2 changes: 0 additions & 2 deletions doc/src/stdlib/arrays.md
Original file line number Diff line number Diff line change
Expand Up @@ -147,7 +147,6 @@ Base.cumprod
Base.cumprod!
Base.cumsum
Base.cumsum!
Base.cumsum_kbn
Base.LinAlg.diff
Base.repeat(::AbstractArray)
Base.rot180
Expand All @@ -156,7 +155,6 @@ Base.rotr90
Base.reducedim
Base.mapreducedim
Base.mapslices
Base.sum_kbn
```

## Combinatorics
Expand Down
23 changes: 0 additions & 23 deletions test/arrayops.jl
Original file line number Diff line number Diff line change
Expand Up @@ -870,29 +870,6 @@ end
@test c[2,2] == A[3,1]+A[3,3]
end

v = [1,1e100,1,-1e100]*1000
v2 = [1,-1e100,1,1e100]*1000

cv = [1,1e100,1e100,2]*1000
cv2 = [1,-1e100,-1e100,2]*1000

@test isequal(cumsum_kbn(v), cv)
@test isequal(cumsum_kbn(v2), cv2)

A = [v reverse(v) v2 reverse(v2)]

c = cumsum_kbn(A, 1)

@test isequal(c[:,1], cv)
@test isequal(c[:,3], cv2)
@test isequal(c[4,:], [2.0, 2.0, 2.0, 2.0]*1000)

c = cumsum_kbn(A, 2)

@test isequal(c[1,:], cv2)
@test isequal(c[3,:], cv)
@test isequal(c[:,4], [2.0,2.0,2.0,2.0]*1000)

@test repeat(BitMatrix(Matrix(I, 2, 2)), inner = (2,1), outer = (1,2)) ==
repeat(Matrix(I, 2, 2), inner = (2,1), outer = (1,2))
end
Expand Down
22 changes: 14 additions & 8 deletions test/offsetarray.jl
Original file line number Diff line number Diff line change
Expand Up @@ -383,14 +383,20 @@ I,J,N = findnz(z)
@test vecnorm(A) vecnorm(parent(A))
@test vecdot(v, v) vecdot(v0, v0)

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)))
# Prior to its removal from Base, cumsum_kbn was used here. To achieve the same level of
# accuracy in the tests, we need to use BigFloats with enlarged precision.
@testset "high-precision array reduction" begin
setprecision(BigFloat, 500) do
v = OffsetArray(BigFloat[1,1e100,1,-1e100], (-3,)) .* 1000
v2 = OffsetArray(BigFloat[1,-1e100,1,1e100], ( 5,)) .* 1000
@test isa(v, OffsetArray)
cv = OffsetArray(BigFloat[1, 1e100, 1e100,2], (-3,)) .* 1000
cv2 = OffsetArray(BigFloat[1,-1e100,-1e100,2], ( 5,)) .* 1000
@test cumsum(v) cv
@test cumsum(v2) cv2
@test sum(v) sum(parent(v))
end
end

io = IOBuffer()
writedlm(io, A)
Expand Down
15 changes: 1 addition & 14 deletions test/reduce.jl
Original file line number Diff line number Diff line change
Expand Up @@ -135,18 +135,6 @@ for f in (sum3, sum4, sum7, sum8)
end
@test typeof(sum(Int8[])) == typeof(sum(Int8[1])) == typeof(sum(Int8[1 7]))

@test sum_kbn([1,1e100,1,-1e100]) === 2.0
@test sum_kbn(Float64[]) === 0.0
@test sum_kbn(i for i=1.0:1.0:10.0) === 55.0
@test sum_kbn(i for i=1:1:10) === 55
@test sum_kbn([1 2 3]) === 6
@test sum_kbn([2+im 3-im]) === 5+0im
@test sum_kbn([1+im 2+3im]) === 3+4im
@test sum_kbn([7 8 9]) === sum_kbn([8 9 7])
@test sum_kbn(i for i=1:1:10) === sum_kbn(i for i=10:-1:1)
@test sum_kbn([-0.0]) === -0.0
@test sum_kbn([-0.0,-0.0]) === -0.0

# check sum(abs, ...) for support of empty collections
@testset "sum(abs, [])" begin
@test @inferred(sum(abs, Float64[])) === 0.0
Expand Down Expand Up @@ -362,14 +350,13 @@ end
## cumsum, cummin, cummax

z = rand(10^6)
let es = sum_kbn(z), es2 = sum_kbn(z[1:10^5])
let es = sum(BigFloat.(z)), es2 = sum(BigFloat.(z[1:10^5]))
@test (es - sum(z)) < es * 1e-13
cs = cumsum(z)
@test (es - cs[end]) < es * 1e-13
@test (es2 - cs[10^5]) < es2 * 1e-13
end


@test sum(collect(map(UInt8,0:255))) == 32640
@test sum(collect(map(UInt8,254:255))) == 509

Expand Down

0 comments on commit 1a4ebd4

Please sign in to comment.