Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Deprecate (cum)sum_kbn #24869

Merged
merged 1 commit into from
Dec 8, 2017
Merged

Deprecate (cum)sum_kbn #24869

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
67 changes: 0 additions & 67 deletions base/abstractarraymath.jl
View file
Open in desktop
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
View file
Open in desktop
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
View file
Open in desktop
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
View file
Open in desktop
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
View file
Open in desktop
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
View file
Open in desktop
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
View file
Open in desktop
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
View file
Open in desktop
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)))
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Perhaps modify rather than eradicate these tests? I conjecture that testing OffsetArray is these tests' purpose rather than testing [cum]sum_kbn specifically :).

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Interestingly, even if I convert these to BigFloats, using cumsum rather than cumsum_kbn seems to give a different answer:

julia> v = OffsetArray(BigFloat[1,1e100,1,-1e100], (-3,))*1000
OffsetArray{BigFloat,1,Array{BigFloat,1}} with indices -2:1:
  1.0e+03
  1.000000000000000015902891109759918046836080856394528138978132755774783877217038e+103
  1.0e+03
 -1.000000000000000015902891109759918046836080856394528138978132755774783877217038e+103

julia> cumsum(v)
OffsetArray{BigFloat,1,Array{BigFloat,1}} with indices -2:1:
 1.0e+03
 1.000000000000000015902891109759918046836080856394528138978132755774783877217038e+103
 1.000000000000000015902891109759918046836080856394528138978132755774783877217038e+103
 1.0e+03

julia> cumsum_kbn(v)
OffsetArray{BigFloat,1,Array{BigFloat,1}} with indices -2:1:
 1.0e+03
 1.000000000000000015902891109759918046836080856394528138978132755774783877217038e+103
 1.000000000000000015902891109759918046836080856394528138978132755774783877217038e+103
 2.0e+03

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Perhaps need higher precision, and/or cumsum doesn't preserve precision in the accumulation?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'd just use isapprox here. The fact that the result was exact before must have been due to the way the KBN algorithm worked in that particular case, but we don't really care here since the goal is just to test that it works on OffsetArray.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The problem with isapprox in this case is that 1000 and 2000 aren't approximately equal. 😛

Copy link
Member

@nalimilan nalimilan Dec 4, 2017 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

OK, I misread. That's just because BigFloat doesn't have enough precision by default. It works with a higher precision:

julia> setprecision(BigFloat, 500)
500

julia> cumsum(BigFloat[1,1e100,1,-1e100]*1000)
4-element Array{BigFloat,1}:
 1.0e+03                                                                                                    
 1.0000000000000000159028911097599180468360808563945281389781327557747838772170381060813469985856815105e+103
 1.0000000000000000159028911097599180468360808563945281389781327557747838772170381060813469985856815106e+103
 2.0e+03                                                                                                    

julia> cumsum_kbn(BigFloat[1,1e100,1,-1e100]*1000)
4-element Array{BigFloat,1}:
 1.0e+03                                                                                                    
 1.0000000000000000159028911097599180468360808563945281389781327557747838772170381060813469985856815105e+103
 1.0000000000000000159028911097599180468360808563945281389781327557747838772170381060813469985856815106e+103
 2.0e+03

But it's really weird to mix testing OffsetArray with testing precision with very high values. Maybe split them into two separate tests, with a comment saying what is being tested?

# 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
View file
Open in desktop
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