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

Jump to bottom

cumsum fixes (fixes #18363 and #18336) #18364

Merged
merged 10 commits into from
Sep 8, 2016
Merged

cumsum fixes (fixes #18363 and #18336) #18364

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
a0f30e6
cumsum fixes (fixes #18363 and #18336)
stevengj Sep 5, 2016
4263ffb
use r_promote_type, similar to r_promote, in cumsum
stevengj Sep 5, 2016
4fca7c1
combine r_promote and r_promote_type
stevengj Sep 5, 2016
deef68d
NEWS for breaking change
stevengj Sep 5, 2016
1980212
use length(linearindices), and also throw error on empty case to be c…
stevengj Sep 5, 2016
6821c02
specify cumsum return type for sparse matrix constructor
stevengj Sep 6, 2016
f3b745f
make rcum_promote_type less eager to widen than r_promote_type
stevengj Sep 6, 2016
e4a0988
rm obsolete NEWS link
stevengj Sep 6, 2016
fc6643e
change BoundsError to DimensionMismatch
stevengj Sep 7, 2016
9ed4912
fix test to match DimensionMismatch exception type
stevengj Sep 8, 2016
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
24 changes: 17 additions & 7 deletions base/arraymath.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -446,8 +446,14 @@ 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])
# 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}
Copy link
Contributor

@pabloferz pabloferz Sep 8, 2016 •
edited
Loading

Choose a reason for hiding this comment

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

I think you can do

rcum_promote_type{T}(op, ::Type{T}) = promote_eltype_op(op, T)
rcum_promote_type{T}(op, ::Type{Array{T,N}) = Array{rcum_promote_type(op, T), N}

without specializing for <:Number

Copy link
Member Author

@stevengj stevengj Sep 8, 2016 •
edited
Loading

Choose a reason for hiding this comment

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

promote_eltype_op gives the wrong answers for non-Number arguments too. e.g. it gives Base.promote_eltype_op(+, Range{Int}) --> Int when I want Range{Int}.

Copy link
Contributor

@pabloferz pabloferz Sep 8, 2016 •
edited
Loading

Choose a reason for hiding this comment

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

Yeah, that's why I left the second definition. Maybe some day, if we get triangular dispatch or something like it, that will be easier.

EDIT: I see that for Range that won't work unless we also have something like triangular dispatch.


for (f, f!, fp, op) = ((:cumsum, :cumsum!, :cumsum_pairwise!, :+),
(:cumprod, :cumprod!, :cumprod_pairwise!, :*) )
Expand All @@ -470,14 +476,18 @@ for (f, f!, fp, op) = ((:cumsum, :cumsum!, :cumsum_pairwise!, :+),
end

@eval function ($f!)(result::AbstractVector, v::AbstractVector)
n = length(v)
li = linearindices(v)
Copy link
Sponsor Member

Choose a reason for hiding this comment

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

Best would be to have n = length(li) since that's guaranteed to work even if v has non-1 indices.

Copy link
Member Author

Choose a reason for hiding this comment

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

Doesn't length(v) return the number of elements in the array regardless? Or is length(v) equivalent to last(linearindices(v))? The latter would seem odd to me.

Copy link
Member Author

Choose a reason for hiding this comment

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

Ah, I see that the length documentation says it returns "For ordered, indexable collections, the maximum index i for which getindex(collection, i) is valid."

Okay, changed it to length(li) as you suggest. It still seems strange to me.

(Note also that this cumsum method is only for AbstractVector, and the length documentation says that it "Returns the number of elements in A." for any AbstractArray.)

Copy link
Sponsor Member

Choose a reason for hiding this comment

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

http://docs.julialang.org/en/latest/devdocs/offset-arrays/#background

Though if this isn't to be backported, then perhaps you could just stick with the original, since the situation with length is just for 0.5. However, packages probably haven't started preparing for 0.6 yet, so they may not support length yet for unconventional arrays.

Copy link
Sponsor Member

Choose a reason for hiding this comment

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

I think the definition/documentation for length should be revisited. It makes more sense to me for it to always return an element count, and have nothing to do with indices. The docs should say "returns the number of elements generated by an iterator".

li != linearindices(result) && throw(DimensionMismatch("input and output array sizes and indices must match"))
n = length(li)
if n == 0; return result; end
($fp)(v, result, $(op==:+ ? :(zero(first(v))) : :(one(first(v)))), first(indices(v,1)), n)
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, rcum_promote_type($op, T)), v)
end
end
2 changes: 1 addition & 1 deletion base/multidimensional.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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{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)
Expand Down
31 changes: 14 additions & 17 deletions base/reduce.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,24 +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)
# 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

Expand Down
2 changes: 1 addition & 1 deletion base/sparse/sparsematrix.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
12 changes: 12 additions & 0 deletions test/arrayops.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1689,6 +1689,18 @@ 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 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]
@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
@test cumprod(-0.0im + (0:0))[1] === Complex(0.0, -0.0)

module TestNLoops15895

using Base.Cartesian
Expand Down