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

stevengj · 2016-09-05T18:33:32Z

This fixes #18336 and #18363 — there were an embarrassing number of problems with cumsum and cumprod for corner cases. (My fault, I think, since I wrote this routine back in #4039.)

The breaking change is that cumsum(v) now uses similar(v, rcum_promote_type(+, eltype(v))) rather than the broken _cumsum_type(v) function, where rcum_promote_type is a new function that calls promote_op for numbers but which leaves most other types alone. I don't think this should change any working cases for cumsum of numeric types, but it might break some unusual user-defined types in which + produces a different type.

(The user can still control the result type manually by passing an array to cumsum!, however.)

stevengj · 2016-09-05T18:45:17Z

Also fixes #13244.

stevengj · 2016-09-05T19:02:24Z

(I thought about using r_promote or similar to widen the cumsum result. However, widening makes more sense for a reduce operation, where you are only widening a single scalar result, than in a function like this where you are allocating a whole array. Also, the user can always widen manually if desired by passing the desired array type to cumsum!, unlike reduce.)

stevengj · 2016-09-05T19:06:01Z

As @andreasnoack pointed out, this does not produce the expected result for Bool arrays, so maybe we should use something like r_promote after all.

I'm conflicted on whether we should widen a sum of, say UInt16, though.

stevengj · 2016-09-05T19:33:01Z

Okay, updated it to use the same type as sum. It just seems easier to understand if we are consistent.

timholy · 2016-09-05T21:01:40Z

base/arraymath.jl

@@ -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)


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

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.

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.)

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.

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".

…onsistent with non-empty case for size mismatches

StefanKarpinski · 2016-09-05T21:50:10Z

Slightly tangential but I really think we should probably make Bool
non-numeric like Char.

On Monday, September 5, 2016, Steven G. Johnson notifications@github.com
wrote:

As @andreasnoack https://github.com/andreasnoack pointed out, this does
not produce the expected result for Bool arrays, so maybe we should use
something like r_promote after all.

I'm conflicted on whether we should widen a sum of, say UInt16, though.

—
You are receiving this because you are subscribed to this thread.
Reply to this email directly, view it on GitHub
#18364 (comment),
or mute the thread
https://github.com/notifications/unsubscribe-auth/AAJX_Pb7xDAqkRfWioqfIiN4dBg8RDGpks5qnGgbgaJpZM4J1QuX
.

TotalVerb · 2016-09-05T23:15:37Z

I agree with @StefanKarpinski. There really ought to be a UInt1 type to model the integers mod 2 (a field that's small enough to promote upward to any other type) to replace the current Bool hackery.

I opened #18367 for that discussion.

stevengj · 2016-09-06T11:18:04Z

Seems like a random Travis timeout on OSX. I'll try restarting CI. (Is there a way to just restart Travis?)

stevengj · 2016-09-06T14:09:32Z

Darn it, now Travis is succeeding but one of the AppVeyor builds is timing out. These random CI failures are really frustrating.

tkelman · 2016-09-06T14:34:18Z

You can restart specific CI services through their UI. I also request that you please make a backup of travis failure logs to a gist before restarting, as restarted builds there overwrite previous logs.

tkelman · 2016-09-06T15:58:46Z

I'm conflicted on whether we should widen a sum of, say UInt16, though.

I'd much rather not widen UInt16 or UInt8 or Float16 arrays in cumsum, if I'm using an array of that type I am probably doing so because I want to reduce the memory usage. For reductions to a smaller output size this is less important (and the overflow prevention more worth the tradeoff) than something that returns an array of the same size as its input.

stevengj · 2016-09-06T17:15:01Z

@tkelman, okay, so it should just use typeof(zero(T)+zero(T)) for T<:Number?

tkelman · 2016-09-06T17:30:50Z

I think that sounds sane. I'm not sure which corner cases the current typeof(+zero(T)) for Number and
typeof(v[1]+v[1]) otherwise are intended to address. There might be non-Number types where addition results in a different type, but those could be tricky to deal with.

andreasnoack · 2016-09-06T18:13:14Z

typeof(zero(T)+zero(T)) is not correct. See #14237. Why not promote_op(+,T,S)?

stevengj · 2016-09-06T18:20:24Z

@tkelman, the main current corner case seems to be Bool, which gets promoted to Int for addition.

@andreasnoack, promote_op(+,T,S) seems fine, although promote_op seems like it might go away at some point.

stevengj · 2016-09-06T18:28:14Z

Hmm, no, it looks like promote_op is too conservative about falling back to Any. e.g. Base.promote_op(+, Vector{Int}, Vector{Int}) returns Any. I guess I could only use promote_op for T<:Number.

stevengj · 2016-09-06T18:39:43Z

Hmm, also need to handle cumsum of Vector{Vector{Bool}} and similar, grr.

stevengj · 2016-09-06T19:42:20Z

Okay, the new behavior should be more consistent with the old cumsum in not widening numeric types (except Bool).

stevengj · 2016-09-07T21:13:41Z

Better?

timholy · 2016-09-07T22:00:05Z

base/arraymath.jl

@@ -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)
+        li != linearindices(result) && throw(BoundsError())


Should this be an DimensionMismatch and include a helpful message about the nature of the problem?

Great. Test needs to change too.

whoops, right.

timholy · 2016-09-07T22:05:39Z

Aside from one small comment, LGTM.

pabloferz · 2016-09-08T14:17:03Z

base/arraymath.jl

+
+# 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}


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

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}.

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.

StefanKarpinski · 2016-09-08T19:04:45Z

This is a (minor) breaking change, but the previous behavior could be argued to be buggy. Should we consider backporting this fix?

stevengj · 2016-09-08T21:20:12Z

I would say the odds of this breaking someone's code is pretty low. About the only cases that should actually break would be cumsum of an array of non-Array collections of Bool, or cumsum of an array of some other weird non-Number type for which + produces a different type. I doubt there are any actual examples of such code in the wild.

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

a0f30e6

stevengj mentioned this pull request Sep 5, 2016

problem with cumsum #13244

Closed

use r_promote_type, similar to r_promote, in cumsum

4263ffb

stevengj added the breaking This change will break code label Sep 5, 2016

stevengj mentioned this pull request Sep 5, 2016

cumsum! does not check array sizes #18363

Closed

stevengj added 2 commits September 5, 2016 15:52

combine r_promote and r_promote_type

4fca7c1

NEWS for breaking change

deef68d

timholy reviewed Sep 5, 2016
View reviewed changes

kshyatt added the maths Mathematical functions label Sep 5, 2016

use length(linearindices), and also throw error on empty case to be c…

1980212

…onsistent with non-empty case for size mismatches

specify cumsum return type for sparse matrix constructor

6821c02

stevengj closed this Sep 6, 2016

stevengj reopened this Sep 6, 2016

stevengj closed this Sep 6, 2016

stevengj reopened this Sep 6, 2016

make rcum_promote_type less eager to widen than r_promote_type

f3b745f

rm obsolete NEWS link

e4a0988

timholy reviewed Sep 7, 2016
View reviewed changes

stevengj and others added 2 commits September 7, 2016 19:26

change BoundsError to DimensionMismatch

fc6643e

fix test to match DimensionMismatch exception type

9ed4912

pabloferz reviewed Sep 8, 2016
View reviewed changes

timholy merged commit c7a4897 into JuliaLang:master Sep 8, 2016

stevengj deleted the cumsum_fixes branch September 8, 2016 18:02

jw3126 mentioned this pull request Oct 21, 2016

Added accumulate, accumulate! #18931

Merged

jw3126 mentioned this pull request Jan 16, 2018

fix non numeric accumulate(op, v0, x) (#25506) #25515

Closed

simonbyrne mentioned this pull request Jan 30, 2018

create Accumulate iterator #25766

Open

2 tasks

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

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

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

stevengj commented Sep 5, 2016 •

edited

Loading

stevengj commented Sep 5, 2016

stevengj commented Sep 5, 2016

stevengj commented Sep 5, 2016

stevengj commented Sep 5, 2016

timholy Sep 5, 2016

stevengj Sep 5, 2016

stevengj Sep 5, 2016

timholy Sep 6, 2016

JeffBezanson Sep 7, 2016

StefanKarpinski commented Sep 5, 2016

TotalVerb commented Sep 5, 2016 •

edited

Loading

stevengj commented Sep 6, 2016

stevengj commented Sep 6, 2016

tkelman commented Sep 6, 2016

tkelman commented Sep 6, 2016

stevengj commented Sep 6, 2016

tkelman commented Sep 6, 2016

andreasnoack commented Sep 6, 2016

stevengj commented Sep 6, 2016 •

edited

Loading

stevengj commented Sep 6, 2016

stevengj commented Sep 6, 2016

stevengj commented Sep 6, 2016

stevengj commented Sep 7, 2016

timholy Sep 7, 2016

stevengj Sep 7, 2016

timholy Sep 8, 2016

stevengj Sep 8, 2016

timholy commented Sep 7, 2016

pabloferz Sep 8, 2016 •

edited

Loading

stevengj Sep 8, 2016 •

edited

Loading

pabloferz Sep 8, 2016 •

edited

Loading

StefanKarpinski commented Sep 8, 2016

stevengj commented Sep 8, 2016 •

edited

Loading

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

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

Conversation

stevengj commented Sep 5, 2016 • edited Loading

stevengj commented Sep 5, 2016

stevengj commented Sep 5, 2016

stevengj commented Sep 5, 2016

stevengj commented Sep 5, 2016

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

StefanKarpinski commented Sep 5, 2016

TotalVerb commented Sep 5, 2016 • edited Loading

stevengj commented Sep 6, 2016

stevengj commented Sep 6, 2016

tkelman commented Sep 6, 2016

tkelman commented Sep 6, 2016

stevengj commented Sep 6, 2016

tkelman commented Sep 6, 2016

andreasnoack commented Sep 6, 2016

stevengj commented Sep 6, 2016 • edited Loading

stevengj commented Sep 6, 2016

stevengj commented Sep 6, 2016

stevengj commented Sep 6, 2016

stevengj commented Sep 7, 2016

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

timholy commented Sep 7, 2016

pabloferz Sep 8, 2016 • edited Loading

Choose a reason for hiding this comment

stevengj Sep 8, 2016 • edited Loading

Choose a reason for hiding this comment

pabloferz Sep 8, 2016 • edited Loading

Choose a reason for hiding this comment

StefanKarpinski commented Sep 8, 2016

stevengj commented Sep 8, 2016 • edited Loading

stevengj commented Sep 5, 2016 •

edited

Loading

TotalVerb commented Sep 5, 2016 •

edited

Loading

stevengj commented Sep 6, 2016 •

edited

Loading

pabloferz Sep 8, 2016 •

edited

Loading

stevengj Sep 8, 2016 •

edited

Loading

pabloferz Sep 8, 2016 •

edited

Loading

stevengj commented Sep 8, 2016 •

edited

Loading