From ef5aef44864fd65223443140993d45ef4f329218 Mon Sep 17 00:00:00 2001 From: Jeffrey Sarnoff Date: Mon, 11 Dec 2017 15:48:56 -0500 Subject: [PATCH 01/15] vec -> avec (vec is a function) --- README.md | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/README.md b/README.md index bbfca5d..07fb045 100644 --- a/README.md +++ b/README.md @@ -14,14 +14,14 @@ These functions were formerly part of Julia's Base library. ```julia julia> using KahanSummation -julia> vec = [1.0, 1.0e16, -1.0e16, -0.5]; +julia> avec = [1.0, 1.0e16, -1.0e16, -0.5]; -julia> sum(vec), sum_kbn(vec) +julia> sum(avec), sum_kbn(avec) (-0.5, 0.5) -julia> vec = [1.0, 1.0e16, 1.0, -1.0e16]; +julia> avec = [1.0, 1.0e16, 1.0, -1.0e16]; -julia> cumsum_kbn(vec) .- cumsum(vec) +julia> cumsum_kbn(avec) .- cumsum(avec) 4-element Array{Float64,1}: 0.0 0.0 From 4876eb7bb193107e82cb276ecbe53eccfd31c903 Mon Sep 17 00:00:00 2001 From: Jeffrey Sarnoff Date: Mon, 11 Dec 2017 16:14:24 -0500 Subject: [PATCH 02/15] https://github.com/JuliaMath/KahanSummation.jl/issues/3 (in part) --- src/KahanSummation.jl | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/KahanSummation.jl b/src/KahanSummation.jl index fde6700..e898aff 100644 --- a/src/KahanSummation.jl +++ b/src/KahanSummation.jl @@ -86,10 +86,11 @@ function cumsum_kbn(v::AbstractVector{T}) where T<:AbstractFloat end """ - sum_kbn(A) + sum_kbn([f,] A) Return the sum of all elements of `A`, using the Kahan-Babuska-Neumaier compensated -summation algorithm for additional accuracy. +summation algorithm for additional accuracy. When a function `f` is supplied, the +result of calling f on each element of `A` is summed. """ function sum_kbn(A) T = @default_eltype(typeof(A)) @@ -113,4 +114,6 @@ function sum_kbn(A) s - c end +sum_kbn(f, A) = f.(A) + end # module From 32b1e76a3f47af6f561fbc6d29a824109f3d623d Mon Sep 17 00:00:00 2001 From: Jeffrey Sarnoff Date: Mon, 11 Dec 2017 16:57:58 -0500 Subject: [PATCH 03/15] Update KahanSummation.jl --- src/KahanSummation.jl | 22 +++++++++++++++++++++- 1 file changed, 21 insertions(+), 1 deletion(-) diff --git a/src/KahanSummation.jl b/src/KahanSummation.jl index e898aff..1ae4982 100644 --- a/src/KahanSummation.jl +++ b/src/KahanSummation.jl @@ -114,6 +114,26 @@ function sum_kbn(A) s - c end -sum_kbn(f, A) = f.(A) +function sum_kbn(f,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 = f(Ai) - c + while !(done(A, i)) + Ai, i = next(A, i) + fAi = f(Ai); t = s + fAi + if abs(s) >= abs(fAi) + c -= ((s-t) + fAi) + else + c -= ((fAi-t) + s) + end + s = t + end + s - c +end end # module From 863d46c265c0a231fb94436492200e7fa2ef3470 Mon Sep 17 00:00:00 2001 From: Jeffrey Sarnoff Date: Mon, 11 Dec 2017 17:27:27 -0500 Subject: [PATCH 04/15] test sum_kbn(f, A) --- test/runtests.jl | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/test/runtests.jl b/test/runtests.jl index 133e2e5..9e2ef0d 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -46,3 +46,13 @@ end @test sum_kbn([-0.0]) === -0.0 @test sum_kbn([-0.0,-0.0]) === -0.0 end + +@testset "sum_kbn(f,A)" begin + @test sum_kbn(x->x*x, [1.0]) == 1.0 + @test sum_kbn(x->x*x, [2.0]) == 4.0 + @test sum_kbn(x->x+1, Float64[]) == 0.0 + @test sum_kbn(sqrt, i for i=1.0:1.0:10.0) === 22.4682781862041 + @test sum_kbn(x->x*im,[1+im 2+3im]) == -4+3im + @test sum_kbn(sqrt, [i for i in 1.0:1.0:10.0]) == sum_kbn(sqrt, [i for i in 10.0:-1.0:1.0]) + @test sum_kbn([1.0, 1.0e16, 1.0, -1.0e16]) == sum_kbn(x->x+1, [1.0, 1.0e16, 1.0, -1.0e16] - 1.0) +end From d3451641f749d91f4c7e2a3805e189f7a7d0d16f Mon Sep 17 00:00:00 2001 From: Jeffrey Sarnoff Date: Mon, 11 Dec 2017 17:30:33 -0500 Subject: [PATCH 05/15] respond to comment --- src/KahanSummation.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/KahanSummation.jl b/src/KahanSummation.jl index 1ae4982..4254735 100644 --- a/src/KahanSummation.jl +++ b/src/KahanSummation.jl @@ -136,4 +136,6 @@ function sum_kbn(f,A) s - c end +sum_kbn(A) = sum_kbn(identity, A) + end # module From ee5abb3ff91f4bdd45b9e791b4181f98782842af Mon Sep 17 00:00:00 2001 From: Jeffrey Sarnoff Date: Mon, 11 Dec 2017 17:37:46 -0500 Subject: [PATCH 06/15] fixup --- src/KahanSummation.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/KahanSummation.jl b/src/KahanSummation.jl index 4254735..1fb0caa 100644 --- a/src/KahanSummation.jl +++ b/src/KahanSummation.jl @@ -136,6 +136,6 @@ function sum_kbn(f,A) s - c end -sum_kbn(A) = sum_kbn(identity, A) +sum_kbn(identity, A) = sum_kbn(A) end # module From 423929c62a44707f83dc028c570b5a3356acbbe5 Mon Sep 17 00:00:00 2001 From: Jeffrey Sarnoff Date: Mon, 11 Dec 2017 17:48:34 -0500 Subject: [PATCH 07/15] correct for Method definition sum_kbn(Any, Any) overwrite --- src/KahanSummation.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/KahanSummation.jl b/src/KahanSummation.jl index 1fb0caa..f7757c1 100644 --- a/src/KahanSummation.jl +++ b/src/KahanSummation.jl @@ -114,7 +114,7 @@ function sum_kbn(A) s - c end -function sum_kbn(f,A) +function sum_kbn(f::Function, A) T = @default_eltype(typeof(A)) c = promote_sys_size_add(zero(T)::T) i = start(A) From 3d930138030e32d157a01c730a4414f3b36170c2 Mon Sep 17 00:00:00 2001 From: Jeffrey Sarnoff Date: Mon, 11 Dec 2017 17:50:32 -0500 Subject: [PATCH 08/15] change name of testset --- test/runtests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index 9e2ef0d..b0917ba 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -47,7 +47,7 @@ end @test sum_kbn([-0.0,-0.0]) === -0.0 end -@testset "sum_kbn(f,A)" begin +@testset "sum_kbn_fn" begin @test sum_kbn(x->x*x, [1.0]) == 1.0 @test sum_kbn(x->x*x, [2.0]) == 4.0 @test sum_kbn(x->x+1, Float64[]) == 0.0 From 733f0eafabc7355f9eda58d8efed37b947170c73 Mon Sep 17 00:00:00 2001 From: Jeffrey Sarnoff Date: Mon, 11 Dec 2017 17:52:08 -0500 Subject: [PATCH 09/15] combine test sets --- test/runtests.jl | 3 --- 1 file changed, 3 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index b0917ba..655238e 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -44,10 +44,7 @@ end @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 -end -@testset "sum_kbn_fn" begin @test sum_kbn(x->x*x, [1.0]) == 1.0 @test sum_kbn(x->x*x, [2.0]) == 4.0 @test sum_kbn(x->x+1, Float64[]) == 0.0 From f8915eb7b625a2c53d4f43ae3afbd1542b5dcc9e Mon Sep 17 00:00:00 2001 From: Jeffrey Sarnoff Date: Mon, 11 Dec 2017 18:02:54 -0500 Subject: [PATCH 10/15] test identity --- test/runtests.jl | 16 ++++++++++------ 1 file changed, 10 insertions(+), 6 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index 655238e..efd91b2 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -45,11 +45,15 @@ end @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(x->x*x, [1.0]) == 1.0 - @test sum_kbn(x->x*x, [2.0]) == 4.0 - @test sum_kbn(x->x+1, Float64[]) == 0.0 + @test sum_kbn(identity, [1,1e100,1,-1e100]) === 2.0 + + @test sum_kbn([-0.0]) === -0.0 + @test sum_kbn(x->x*x, [1.0]) === 1.0 + @test sum_kbn(x->x*x, [2.0]) === 4.0 + @test sum_kbn(x->x+1, Float64[]) === 0.0 + @test sum_kbn(x->x+1, [6 7 8]) === sum_kbn([8 9 7]) @test sum_kbn(sqrt, i for i=1.0:1.0:10.0) === 22.4682781862041 - @test sum_kbn(x->x*im,[1+im 2+3im]) == -4+3im - @test sum_kbn(sqrt, [i for i in 1.0:1.0:10.0]) == sum_kbn(sqrt, [i for i in 10.0:-1.0:1.0]) - @test sum_kbn([1.0, 1.0e16, 1.0, -1.0e16]) == sum_kbn(x->x+1, [1.0, 1.0e16, 1.0, -1.0e16] - 1.0) + @test sum_kbn(x->x*im,[1+im 2+3im]) === -4+3im + @test sum_kbn(sqrt, [i for i in 1.0:1.0:10.0]) === sum_kbn(sqrt, [i for i in 10.0:-1.0:1.0]) + @test sum_kbn([1.0, 1.0e16, 1.0, -1.0e16]) === sum_kbn(x->x+1, [1.0, 1.0e16, 1.0, -1.0e16] - 1.0) end From 39bc2743a93912a0dc4b0143c2e8dd3dc266202b Mon Sep 17 00:00:00 2001 From: Jeffrey Sarnoff Date: Mon, 11 Dec 2017 18:09:50 -0500 Subject: [PATCH 11/15] try to make coveralls happy --- test/runtests.jl | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index efd91b2..c5e71cd 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -46,7 +46,16 @@ end @test sum_kbn([-0.0]) === -0.0 @test sum_kbn(identity, [1,1e100,1,-1e100]) === 2.0 - + @test sum_kbn(identity, Float64[]) === 0.0 + @test sum_kbn(identity, i for i=1.0:1.0:10.0) === 55.0 + @test sum_kbn(identity, i for i=1:1:10) === 55 + @test sum_kbn(identity, [1 2 3]) === 6 + @test sum_kbn(identity, [2+im 3-im]) === 5+0im + @test sum_kbn(identity, [1+im 2+3im]) === 3+4im + @test sum_kbn(identity, [7 8 9]) === sum_kbn([8 9 7]) + @test sum_kbn(identity, i for i=1:1:10) === sum_kbn(i for i=10:-1:1) + @test sum_kbn(identity, [-0.0]) === -0.0 + @test sum_kbn([-0.0]) === -0.0 @test sum_kbn(x->x*x, [1.0]) === 1.0 @test sum_kbn(x->x*x, [2.0]) === 4.0 From ad48228cb798a9c02d10abf9c90b22d3b07c5fa8 Mon Sep 17 00:00:00 2001 From: Jeffrey Sarnoff Date: Mon, 11 Dec 2017 18:21:08 -0500 Subject: [PATCH 12/15] fold special case for identity into sum_kbn --- src/KahanSummation.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/KahanSummation.jl b/src/KahanSummation.jl index f7757c1..05ebf28 100644 --- a/src/KahanSummation.jl +++ b/src/KahanSummation.jl @@ -115,6 +115,7 @@ function sum_kbn(A) end function sum_kbn(f::Function, A) + f === identity && return sum_kbn(A) T = @default_eltype(typeof(A)) c = promote_sys_size_add(zero(T)::T) i = start(A) @@ -136,6 +137,4 @@ function sum_kbn(f::Function, A) s - c end -sum_kbn(identity, A) = sum_kbn(A) - end # module From 03c5bc2e79b76875b54bc994c6fe1acd7c9774a8 Mon Sep 17 00:00:00 2001 From: Jeffrey Sarnoff Date: Tue, 12 Dec 2017 14:15:07 -0500 Subject: [PATCH 13/15] Update KahanSummation.jl --- src/KahanSummation.jl | 37 +++++++++++-------------------------- 1 file changed, 11 insertions(+), 26 deletions(-) diff --git a/src/KahanSummation.jl b/src/KahanSummation.jl index 05ebf28..304a099 100644 --- a/src/KahanSummation.jl +++ b/src/KahanSummation.jl @@ -65,13 +65,19 @@ function cumsum_kbn(A::AbstractArray{T}, axis::Integer) where T<:AbstractFloat end function cumsum_kbn(v::AbstractVector{T}) where T<:AbstractFloat - r = similar(v) + r = similar(v)ou can just use BenchmarkTools’ built-in setup argument which is always run before the code to be timed: + +@btime sort!(x) setup=(x = rand(100)) evals=1 +Note that evals=1 is important to ensure that the sort!(x) call is run exactly once after each setup call. isempty(v) && return r inds = indices(v, 1) i1 = first(inds) s = r[i1] = v[i1] c = zero(T) - for i = i1+1:last(inds) + for i = i1+1:last(inds)ou can just use BenchmarkTools’ built-in setup argument which is always run before the code to be timed: + +@btime sort!(x) setup=(x = rand(100)) evals=1 +Note that evals=1 is important to ensure that the sort!(x) call is run exactly once after each setup call. vi = v[i] t = s + vi if abs(s) >= abs(vi) @@ -92,31 +98,8 @@ Return the sum of all elements of `A`, using the Kahan-Babuska-Neumaier compensa summation algorithm for additional accuracy. When a function `f` is supplied, the result of calling f on each element of `A` is summed. """ -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 - function sum_kbn(f::Function, A) - f === identity && return sum_kbn(A) - T = @default_eltype(typeof(A)) + T = Base.return_types(f, (@default_eltype(typeof(A)),))[1] c = promote_sys_size_add(zero(T)::T) i = start(A) if done(A, i) @@ -137,4 +120,6 @@ function sum_kbn(f::Function, A) s - c end +sum_kbn(A) = sum_kbn(identity, A) + end # module From cbbcbccaa2a7a5c16b6c668afe16891b7af5d9ce Mon Sep 17 00:00:00 2001 From: Jeffrey Sarnoff Date: Tue, 12 Dec 2017 14:18:59 -0500 Subject: [PATCH 14/15] maybe ok --- src/KahanSummation.jl | 68 ++++++++++++++++++++++++++++++++++++++----- 1 file changed, 60 insertions(+), 8 deletions(-) diff --git a/src/KahanSummation.jl b/src/KahanSummation.jl index 304a099..bb7118e 100644 --- a/src/KahanSummation.jl +++ b/src/KahanSummation.jl @@ -64,20 +64,72 @@ function cumsum_kbn(A::AbstractArray{T}, axis::Integer) where T<:AbstractFloat return B + C end -function cumsum_kbn(v::AbstractVector{T}) where T<:AbstractFloat - r = similar(v)ou can just use BenchmarkTools’ built-in setup argument which is always run before the code to be timed: +if VERSION >= v"0.7.0-DEV.3000" # TODO: More specific bound + if isdefined(Base, :sum_kbn) # Deprecated + import Base: sum_kbn, cumsum_kbn + else + export sum_kbn, cumsum_kbn + end +end + +if isdefined(Base, Symbol("@default_eltype")) + using Base: @default_eltype +else + macro default_eltype(itr) + quote + Core.Inference.return_type(first, Tuple{$(esc(itr))}) + end + end +end + +if isdefined(Base, :promote_sys_size_add) + using Base: promote_sys_size_add +else + promote_sys_size_add(x::T) where {T} = Base.r_promote(+, zero(T)::T) +end + +""" + 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 + axis_size <= 1 && return A + 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 -@btime sort!(x) setup=(x = rand(100)) evals=1 -Note that evals=1 is important to ensure that the sort!(x) call is run exactly once after each setup call. +function cumsum_kbn(v::AbstractVector{T}) where T<:AbstractFloat + r = similar(v) isempty(v) && return r inds = indices(v, 1) i1 = first(inds) s = r[i1] = v[i1] c = zero(T) - for i = i1+1:last(inds)ou can just use BenchmarkTools’ built-in setup argument which is always run before the code to be timed: - -@btime sort!(x) setup=(x = rand(100)) evals=1 -Note that evals=1 is important to ensure that the sort!(x) call is run exactly once after each setup call. + for i = i1+1:last(inds) vi = v[i] t = s + vi if abs(s) >= abs(vi) From 5cbd33a3970fcfb0c0a2f5d7c6a0cf4035a1abd6 Mon Sep 17 00:00:00 2001 From: Jeffrey Sarnoff Date: Tue, 12 Dec 2017 14:39:12 -0500 Subject: [PATCH 15/15] add more cumsum_kbr(A, dim) tests --- test/runtests.jl | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/test/runtests.jl b/test/runtests.jl index c5e71cd..6d33231 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -31,6 +31,12 @@ import KahanSummation: sum_kbn, cumsum_kbn @test isequal(c[1,:], cv2) @test isequal(c[3,:], cv) @test isequal(c[:,4], [2.0,2.0,2.0,2.0]*1000) + + amat = reshape(collect(1.0:16.0),(4,4)) + + @test cumsum_kbn(amat, 2)[3,:] == [3.0, 10.0, 21.0, 36.0] + @test diag(cumsum_kbn(amat, 1)) == [1.0, 11.0, 30.0, 58.0] + end @testset "sum_kbn" begin