From a0f30e6601f35d5edde89fc7a5d7aecf175eb9e8 Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Mon, 5 Sep 2016 14:30:25 -0400 Subject: [PATCH 01/10] cumsum fixes (fixes #18363 and #18336) --- base/arraymath.jl | 15 ++++++++------- base/multidimensional.jl | 2 +- test/arrayops.jl | 9 +++++++++ 3 files changed, 18 insertions(+), 8 deletions(-) diff --git a/base/arraymath.jl b/base/arraymath.jl index 4d9b9a8f09ed3..4333a0e5050e6 100644 --- a/base/arraymath.jl +++ b/base/arraymath.jl @@ -446,9 +446,6 @@ ctranspose{T<:Real}(A::AbstractVecOrMat{T}) = transpose(A) transpose(x::AbstractVector) = [ transpose(v) for i=of_indices(x, OneTo(1)), v in x ] ctranspose{T}(x::AbstractVector{T}) = T[ ctranspose(v) for i=of_indices(x, OneTo(1)), v in x ] -_cumsum_type{T<:Number}(v::AbstractArray{T}) = typeof(+zero(T)) -_cumsum_type(v) = typeof(v[1]+v[1]) - for (f, f!, fp, op) = ((:cumsum, :cumsum!, :cumsum_pairwise!, :+), (:cumprod, :cumprod!, :cumprod_pairwise!, :*) ) # in-place cumsum of c = s+v[range(i1,n)], using pairwise summation @@ -472,12 +469,16 @@ for (f, f!, fp, op) = ((:cumsum, :cumsum!, :cumsum_pairwise!, :+), @eval function ($f!)(result::AbstractVector, v::AbstractVector) n = length(v) if n == 0; return result; end - ($fp)(v, result, $(op==:+ ? :(zero(first(v))) : :(one(first(v)))), first(indices(v,1)), n) + li = linearindices(v) + li != linearindices(result) && throw(BoundsError()) + i1 = first(li) + @inbounds result[i1] = v1 = v[i1] + n == 1 && return result + ($fp)(v, result, v1, i1+1, n-1) return result end - @eval function ($f)(v::AbstractVector) - c = $(op===:+ ? (:(similar(v,_cumsum_type(v)))) : (:(similar(v)))) - return ($f!)(c, v) + @eval function ($f){T}(v::AbstractVector{T}) + return ($f!)(similar(v), v) end end diff --git a/base/multidimensional.jl b/base/multidimensional.jl index 10a495213f639..e8ae5312fd8b8 100644 --- a/base/multidimensional.jl +++ b/base/multidimensional.jl @@ -527,7 +527,7 @@ julia> cumsum(a,2) 4 9 15 ``` """ -cumsum(A::AbstractArray, axis::Integer=1) = cumsum!(similar(A, Base._cumsum_type(A)), A, axis) +cumsum(A::AbstractArray, axis::Integer=1) = cumsum!(similar(A), A, axis) cumsum!(B, A::AbstractArray) = cumsum!(B, A, 1) """ cumprod(A, dim=1) diff --git a/test/arrayops.jl b/test/arrayops.jl index 05f23a06a76d5..f953d5da90fac 100644 --- a/test/arrayops.jl +++ b/test/arrayops.jl @@ -1689,6 +1689,15 @@ end @test cumsum([1 2; 3 4], 2) == [1 3; 3 7] @test cumsum([1 2; 3 4], 3) == [1 2; 3 4] +# issue #18363 +@test_throws BoundsError cumsum!([0,0], 1:4) +@test cumsum(Any[]) == Any[] && isa(cumsum(Any[]), Vector{Any}) +@test cumsum(Any[1, 2.3]) == [1, 3.3] + +#issue #18336 +@test cumsum([-0.0, -0.0])[1] === cumsum([-0.0, -0.0])[2] === -0.0 +@test cumprod(-0.0im + (0:0))[1] === Complex(0.0, -0.0) + module TestNLoops15895 using Base.Cartesian From 4263ffb6718bc1d71467f5b5f13f91b7e753093a Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Mon, 5 Sep 2016 15:29:46 -0400 Subject: [PATCH 02/10] use r_promote_type, similar to r_promote, in cumsum --- base/arraymath.jl | 2 +- base/multidimensional.jl | 2 +- base/reduce.jl | 6 ++++++ test/arrayops.jl | 2 ++ 4 files changed, 10 insertions(+), 2 deletions(-) diff --git a/base/arraymath.jl b/base/arraymath.jl index 4333a0e5050e6..f797882cdb8e5 100644 --- a/base/arraymath.jl +++ b/base/arraymath.jl @@ -479,6 +479,6 @@ for (f, f!, fp, op) = ((:cumsum, :cumsum!, :cumsum_pairwise!, :+), end @eval function ($f){T}(v::AbstractVector{T}) - return ($f!)(similar(v), v) + return ($f!)(similar(v, r_promote_type($op, T)), v) end end diff --git a/base/multidimensional.jl b/base/multidimensional.jl index e8ae5312fd8b8..3197b412e6740 100644 --- a/base/multidimensional.jl +++ b/base/multidimensional.jl @@ -527,7 +527,7 @@ julia> cumsum(a,2) 4 9 15 ``` """ -cumsum(A::AbstractArray, axis::Integer=1) = cumsum!(similar(A), A, axis) +cumsum{T}(A::AbstractArray{T}, axis::Integer=1) = cumsum!(similar(A, Base.r_promote_type(+, T)), A, axis) cumsum!(B, A::AbstractArray) = cumsum!(B, A, 1) """ cumprod(A, dim=1) diff --git a/base/reduce.jl b/base/reduce.jl index f437f38fe80c7..d23966c6476ae 100644 --- a/base/reduce.jl +++ b/base/reduce.jl @@ -33,6 +33,12 @@ r_promote(::typeof(min), x::WidenReduceResult) = r_promote(scalarmin, x) r_promote(::typeof(max), x) = r_promote(scalarmax, x) r_promote(::typeof(min), x) = r_promote(scalarmin, x) +# like r_promote but acting on types T rather than on values x +r_promote_type(op, T::Type) = T +r_promote_type{T<:WidenReduceResult}(::typeof(+), ::Type{T}) = widen(T) +r_promote_type{T<:WidenReduceResult}(::typeof(*), ::Type{T}) = widen(T) +r_promote_type{T<:Number}(::typeof(+), ::Type{T}) = typeof(zero(T)+zero(T)) +r_promote_type{T<:Number}(::typeof(*), ::Type{T}) = typeof(one(T)*one(T)) ## foldl && mapfoldl diff --git a/test/arrayops.jl b/test/arrayops.jl index f953d5da90fac..ea6c9b8d1be88 100644 --- a/test/arrayops.jl +++ b/test/arrayops.jl @@ -1693,6 +1693,8 @@ end @test_throws BoundsError cumsum!([0,0], 1:4) @test cumsum(Any[]) == Any[] && isa(cumsum(Any[]), Vector{Any}) @test cumsum(Any[1, 2.3]) == [1, 3.3] +@test cumsum([true,true,true]) == [1,2,3] +@test cumsum(0x00:0xff)[end] === 0x00007f80 #issue #18336 @test cumsum([-0.0, -0.0])[1] === cumsum([-0.0, -0.0])[2] === -0.0 From 4fca7c1e1faadabbf4634c9c7e4931efb6a1e79c Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Mon, 5 Sep 2016 15:52:29 -0400 Subject: [PATCH 03/10] combine r_promote and r_promote_type --- base/reduce.jl | 31 +++++++++++-------------------- 1 file changed, 11 insertions(+), 20 deletions(-) diff --git a/base/reduce.jl b/base/reduce.jl index d23966c6476ae..7e8d97d059f2a 100644 --- a/base/reduce.jl +++ b/base/reduce.jl @@ -15,30 +15,21 @@ end typealias CommonReduceResult Union{UInt64,UInt128,Int64,Int128,Float32,Float64} typealias WidenReduceResult Union{SmallSigned, SmallUnsigned, Float16} -# r_promote: promote x to the type of reduce(op, [x]) -r_promote(op, x::WidenReduceResult) = widen(x) -r_promote(op, x) = x -r_promote(::typeof(+), x::WidenReduceResult) = widen(x) -r_promote(::typeof(*), x::WidenReduceResult) = widen(x) -r_promote(::typeof(+), x::Number) = oftype(x + zero(x), x) -r_promote(::typeof(*), x::Number) = oftype(x * one(x), x) -r_promote(::typeof(+), x) = x -r_promote(::typeof(*), x) = x -r_promote(::typeof(scalarmax), x::WidenReduceResult) = x -r_promote(::typeof(scalarmin), x::WidenReduceResult) = x -r_promote(::typeof(scalarmax), x) = x -r_promote(::typeof(scalarmin), x) = x -r_promote(::typeof(max), x::WidenReduceResult) = r_promote(scalarmax, x) -r_promote(::typeof(min), x::WidenReduceResult) = r_promote(scalarmin, x) -r_promote(::typeof(max), x) = r_promote(scalarmax, x) -r_promote(::typeof(min), x) = r_promote(scalarmin, x) - -# like r_promote but acting on types T rather than on values x -r_promote_type(op, T::Type) = T +# r_promote_type: promote T to the type of reduce(op, ::Array{T}) +# (some "extra" methods are required here to avoid ambiguity warnings) +r_promote_type{T}(op, ::Type{T}) = T +r_promote_type{T<:WidenReduceResult}(op, ::Type{T}) = widen(T) r_promote_type{T<:WidenReduceResult}(::typeof(+), ::Type{T}) = widen(T) r_promote_type{T<:WidenReduceResult}(::typeof(*), ::Type{T}) = widen(T) r_promote_type{T<:Number}(::typeof(+), ::Type{T}) = typeof(zero(T)+zero(T)) r_promote_type{T<:Number}(::typeof(*), ::Type{T}) = typeof(one(T)*one(T)) +r_promote_type{T<:WidenReduceResult}(::typeof(scalarmax), ::Type{T}) = T +r_promote_type{T<:WidenReduceResult}(::typeof(scalarmin), ::Type{T}) = T +r_promote_type{T<:WidenReduceResult}(::typeof(max), ::Type{T}) = T +r_promote_type{T<:WidenReduceResult}(::typeof(min), ::Type{T}) = T + +# r_promote: promote x to the type of reduce(op, [x]) +r_promote{T}(op, x::T) = convert(r_promote_type(op, T), x) ## foldl && mapfoldl From deef68d3e2df44ff778bd4a41551da6535ad96e0 Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Mon, 5 Sep 2016 16:17:48 -0400 Subject: [PATCH 04/10] NEWS for breaking change --- NEWS.md | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/NEWS.md b/NEWS.md index 0cd0ee6f1f31e..6a45ae9628d88 100644 --- a/NEWS.md +++ b/NEWS.md @@ -24,6 +24,10 @@ This section lists changes that do not have deprecation warnings. * Keyword arguments are processed left-to-right: if the same keyword is specified more than once, the rightmost occurrence takes precedence ([#17785]). + * The `cumsum` and `cumprod` functions now use the same type promotion as `sum` and `prod`, + and in particular their numeric result types are now widened to at least 32 bits; + you can use `cumsum!` or `cumprod!` to specify a different type ([#18364]). + Library improvements -------------------- @@ -648,3 +652,4 @@ Language tooling improvements [#17546]: https://github.com/JuliaLang/julia/issues/17546 [#17668]: https://github.com/JuliaLang/julia/issues/17668 [#17785]: https://github.com/JuliaLang/julia/issues/17785 +[#18364]: https://github.com/JuliaLang/julia/issues/18364 From 1980212ece58bb5198a81c65ff18b6d6bad53f1a Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Mon, 5 Sep 2016 17:24:07 -0400 Subject: [PATCH 05/10] use length(linearindices), and also throw error on empty case to be consistent with non-empty case for size mismatches --- base/arraymath.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/base/arraymath.jl b/base/arraymath.jl index f797882cdb8e5..a77744e9f95b7 100644 --- a/base/arraymath.jl +++ b/base/arraymath.jl @@ -467,10 +467,10 @@ for (f, f!, fp, op) = ((:cumsum, :cumsum!, :cumsum_pairwise!, :+), end @eval function ($f!)(result::AbstractVector, v::AbstractVector) - n = length(v) - if n == 0; return result; end li = linearindices(v) li != linearindices(result) && throw(BoundsError()) + n = length(li) + if n == 0; return result; end i1 = first(li) @inbounds result[i1] = v1 = v[i1] n == 1 && return result From 6821c022d04faedd8e7418b9bae86143d315eb2f Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Mon, 5 Sep 2016 20:40:57 -0400 Subject: [PATCH 06/10] specify cumsum return type for sparse matrix constructor --- base/sparse/sparsematrix.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/base/sparse/sparsematrix.jl b/base/sparse/sparsematrix.jl index d150815d7a478..6d4f366f287e7 100644 --- a/base/sparse/sparsematrix.jl +++ b/base/sparse/sparsematrix.jl @@ -353,7 +353,7 @@ function sparse_IJ_sorted!{Ti<:Integer}(I::AbstractVector{Ti}, J::AbstractVector end end - colptr = cumsum(cols) + colptr = cumsum!(similar(cols), cols) # Allow up to 20% slack if ndups > 0.2*L From f3b745f3564d89e38b6be2685723cc8493757873 Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Tue, 6 Sep 2016 15:39:08 -0400 Subject: [PATCH 07/10] make rcum_promote_type less eager to widen than r_promote_type --- NEWS.md | 4 ---- base/arraymath.jl | 11 ++++++++++- base/multidimensional.jl | 2 +- test/arrayops.jl | 7 ++++--- 4 files changed, 15 insertions(+), 9 deletions(-) diff --git a/NEWS.md b/NEWS.md index 6a45ae9628d88..3b3608c6af525 100644 --- a/NEWS.md +++ b/NEWS.md @@ -24,10 +24,6 @@ This section lists changes that do not have deprecation warnings. * Keyword arguments are processed left-to-right: if the same keyword is specified more than once, the rightmost occurrence takes precedence ([#17785]). - * The `cumsum` and `cumprod` functions now use the same type promotion as `sum` and `prod`, - and in particular their numeric result types are now widened to at least 32 bits; - you can use `cumsum!` or `cumprod!` to specify a different type ([#18364]). - Library improvements -------------------- diff --git a/base/arraymath.jl b/base/arraymath.jl index a77744e9f95b7..4983ca6edb2db 100644 --- a/base/arraymath.jl +++ b/base/arraymath.jl @@ -446,6 +446,15 @@ ctranspose{T<:Real}(A::AbstractVecOrMat{T}) = transpose(A) transpose(x::AbstractVector) = [ transpose(v) for i=of_indices(x, OneTo(1)), v in x ] ctranspose{T}(x::AbstractVector{T}) = T[ ctranspose(v) for i=of_indices(x, OneTo(1)), v in x ] +# 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}(op, ::Type{T}) = T + +# handle sums of Vector{Bool} and similar. it would be nice to handle +# any AbstractArray here, but it's not clear how that would be possible +rcum_promote_type{T,N}(op, ::Type{Array{T,N}}) = Array{rcum_promote_type(op,T), N} + for (f, f!, fp, op) = ((:cumsum, :cumsum!, :cumsum_pairwise!, :+), (:cumprod, :cumprod!, :cumprod_pairwise!, :*) ) # in-place cumsum of c = s+v[range(i1,n)], using pairwise summation @@ -479,6 +488,6 @@ for (f, f!, fp, op) = ((:cumsum, :cumsum!, :cumsum_pairwise!, :+), end @eval function ($f){T}(v::AbstractVector{T}) - return ($f!)(similar(v, r_promote_type($op, T)), v) + return ($f!)(similar(v, rcum_promote_type($op, T)), v) end end diff --git a/base/multidimensional.jl b/base/multidimensional.jl index 3197b412e6740..591e0c46c89a9 100644 --- a/base/multidimensional.jl +++ b/base/multidimensional.jl @@ -527,7 +527,7 @@ julia> cumsum(a,2) 4 9 15 ``` """ -cumsum{T}(A::AbstractArray{T}, axis::Integer=1) = cumsum!(similar(A, Base.r_promote_type(+, T)), A, axis) +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) """ cumprod(A, dim=1) diff --git a/test/arrayops.jl b/test/arrayops.jl index ea6c9b8d1be88..0fb70e267b94d 100644 --- a/test/arrayops.jl +++ b/test/arrayops.jl @@ -1691,10 +1691,11 @@ end # issue #18363 @test_throws BoundsError cumsum!([0,0], 1:4) -@test cumsum(Any[]) == Any[] && isa(cumsum(Any[]), Vector{Any}) -@test cumsum(Any[1, 2.3]) == [1, 3.3] +@test cumsum(Any[])::Vector{Any} == Any[] +@test cumsum(Any[1, 2.3])::Vector{Any} == [1, 3.3] == cumsum(Real[1, 2.3])::Vector{Real} @test cumsum([true,true,true]) == [1,2,3] -@test cumsum(0x00:0xff)[end] === 0x00007f80 +@test cumsum(0x00:0xff)[end] === 0x80 # overflow +@test cumsum([[true], [true], [false]])::Vector{Vector{Int}} == [[1], [2], [2]] #issue #18336 @test cumsum([-0.0, -0.0])[1] === cumsum([-0.0, -0.0])[2] === -0.0 From e4a09882b2b3753e169e257dd1a47a004a773362 Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Tue, 6 Sep 2016 15:42:54 -0400 Subject: [PATCH 08/10] rm obsolete NEWS link --- NEWS.md | 1 - 1 file changed, 1 deletion(-) diff --git a/NEWS.md b/NEWS.md index 3b3608c6af525..0cd0ee6f1f31e 100644 --- a/NEWS.md +++ b/NEWS.md @@ -648,4 +648,3 @@ Language tooling improvements [#17546]: https://github.com/JuliaLang/julia/issues/17546 [#17668]: https://github.com/JuliaLang/julia/issues/17668 [#17785]: https://github.com/JuliaLang/julia/issues/17785 -[#18364]: https://github.com/JuliaLang/julia/issues/18364 From fc6643ece48824e955b5676dc4ee09af4d0306fe Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Wed, 7 Sep 2016 19:26:11 -0400 Subject: [PATCH 09/10] change BoundsError to DimensionMismatch --- base/arraymath.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/base/arraymath.jl b/base/arraymath.jl index 4983ca6edb2db..5fbe638d37691 100644 --- a/base/arraymath.jl +++ b/base/arraymath.jl @@ -477,7 +477,7 @@ for (f, f!, fp, op) = ((:cumsum, :cumsum!, :cumsum_pairwise!, :+), @eval function ($f!)(result::AbstractVector, v::AbstractVector) li = linearindices(v) - li != linearindices(result) && throw(BoundsError()) + li != linearindices(result) && throw(DimensionMismatch("input and output array sizes and indices must match")) n = length(li) if n == 0; return result; end i1 = first(li) From 9ed4912d90c16c283a1edbd00d68cf3606180e55 Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Thu, 8 Sep 2016 10:02:41 -0400 Subject: [PATCH 10/10] fix test to match DimensionMismatch exception type --- test/arrayops.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/arrayops.jl b/test/arrayops.jl index 0fb70e267b94d..74014aeb4e3ff 100644 --- a/test/arrayops.jl +++ b/test/arrayops.jl @@ -1690,7 +1690,7 @@ end @test cumsum([1 2; 3 4], 3) == [1 2; 3 4] # issue #18363 -@test_throws BoundsError cumsum!([0,0], 1:4) +@test_throws DimensionMismatch cumsum!([0,0], 1:4) @test cumsum(Any[])::Vector{Any} == Any[] @test cumsum(Any[1, 2.3])::Vector{Any} == [1, 3.3] == cumsum(Real[1, 2.3])::Vector{Real} @test cumsum([true,true,true]) == [1,2,3]