Change accumulation promotion behaviour

[simonbyrne]

This is an interim fix until we can get the performance of #25766 to be acceptable.

It changes the promotion behaviour of small integers in cumsum/cumprod, and includes the fix from #25515.

[mbauman]

I verified that the performance of the dims keyword methods are equivalent — but cumsum(::Vector) is broken:

julia> A = rand(100);

julia> cumsum(A);
ERROR: MethodError: getfield(Base, Symbol("#kw##cumsum!"))()(::NamedTuple{(:dims,),Tuple{Int64}}, ::typeof(cumsum!), ::Array{Float64,1}, ::Array{Float64,1}) is ambiguous. Candidates:
  (::getfield(Base, Symbol("#kw##cumsum!")))(::Any, ::typeof(cumsum!), out, v::AbstractArray{T,1} where T) in Base
  (::getfield(Base, Symbol("#kw##cumsum!")))(::Any, ::typeof(cumsum!), B::AbstractArray{T,N} where N, A) where T in Base
Possible fix, define
  (::getfield(Base, Symbol("#kw##cumsum!")))(::Any, ::typeof(cumsum!), ::AbstractArray{T,N} where N, ::AbstractArray{T,1} where T)
Stacktrace:
 [1] #cumsum#602(::Int64, ::Function, ::Array{Float64,1}) at ./accumulate.jl:93
 [2] #cumsum at ./<missing>:0 [inlined]
 [3] cumsum(::Array{Float64,1}) at ./accumulate.jl:118

[StefanKarpinski]

Triage accepts and looks forward to a future Accumulate type that is fast :)

[JeffBezanson]

👍 I think it makes sense for sum and cumsum to match (one can maybe debate what they should do, but they should certainly match).

[StefanKarpinski]

This already had a CI straight flush, but we wanted to squash the last three commits but not the first one, so I rebased on top of master and repushed it. Merge as soon as CI passes.

[tkelman]

I don't think this is a good idea now, as it wasn't a good idea before. cumsum returns an array, so it's much more sensitive to its element type than sum. It's very telling that of the few uses of this in base, at least one of them (that tests identified) wanted the accumulate non-promoting behavior. In packages, this would come across as a no-deprecation breaking change where they may suddenly consume more memory and get a different element type than they had been expecting, in return for somewhat lower chance of overflow.

Why does cumsum still exist, seems it could easily be deprecated in favor of the more general and no-longer accumulate(+, ...)

[mbauman]

I'm not sure that's so telling — that line you identify (colptrB = cumsum(colptrB)) would be better as simply cumsum!. Any time we construct a new array as the result of some computation we're faced with this difficulty. We've converged on using the scalar behaviors to choose the result. I agree this is problematic in cases like sparse array index computations, but I find this to be the most sensible overarching rule.

I do find the disconnect between cumsum and accumulate(+, …) that this introduces a bit disconcerting. Would this prevent us from introducing the general Accumulate iterator in 1.x? Or if we did, would Accumulate and accumulate differ slightly in their array types? That part seems more worrisome.

[StefanKarpinski]

We've already accepted the exact same difference between sum and reduce(+, ...) and that has worked well; the difference between cumsum and accumulate(+, ...) exactly mirrors that. Making ad hoc distinctions between functions that produce scalars and arrays just leads to unpredictable, inconsistent behaviors. There is simply no choice here that will work well for everybody in every situation: sometimes you want to widen when reducing/accumulating, sometimes you don't. We now make both behaviors easy to express. That is the best that can be done in such situations.

I find the argument that this causes a dangerous silent breakage uncompelling since this change goes in the safe direction: from non-widening, which may overflow and give incorrect values, to widening which is certain to give correct answers if it did before. In some cases this change will be a bug fix. In some cases, code may have been relying on overflow for correctness, in which case it will now get an InexactError when too-large values are extracted from the accumulated array and used in an unwidedened location. The only other problem is the potential memory blow up from a widened type, which is also relatively safe and discoverable with standard performance analysis tools.

[simonbyrne]

I do find the disconnect between cumsum and accumulate(+, …) that this introduces a bit disconcerting. Would this prevent us from introducing the general Accumulate iterator in 1.x? Or if we did, would Accumulate and accumulate differ slightly in their array types? That part seems more worrisome.

It won't make a difference to the Accumulate iterator: accumulate(op,X) will just be collect(Accumulate(op, X)); we will still have cumsum(X) = accumulate(add_sum, X) (in the same way sum(X) = reduce(add_sum, X)).

[tkelman]

We've converged on using the scalar behaviors to choose the result.

Right, but which scalar behavior is appropriate for cumsum? Is it accumulate(+, ...) or accumulate(add_sum, ...) ? I don't think silently switching from the former, which it has always been, to the latter, is a hands-down improvement.

The only other problem is the potential memory blow up from a widened type, which is also relatively safe and discoverable with standard performance analysis tools.

It's not merely a memory concern, it's also going to break API's that aren't prepared to be getting a different element type than their inputs had. Functions that work with arrays will be pickier about their element types a bit more often than functions that work with scalars, thanks to covariance and the frequent desire or requirement for element types to be homogenous in arrays. Widening of element types in a handful of privileged array functions seems more ad-hoc to me than predictably going through the standard accumulate code paths using the intuitively expected addition function (vs the implementation detail add_sum function that most people won't commonly use or necessarily want).

There is simply no choice here that will work well for everybody in every situation: sometimes you want to widen when reducing/accumulating, sometimes you don't. We now make both behaviors easy to express. That is the best that can be done in such situations.

a) which behavior is all of the code in the wild that has been using cumsum expecting, and
b) which is actually more intuitive and expected?

I don't think the answer to either of those questions is clear-cut and unanimous enough to say the promoting behavior suddenly deserves the cumsum name when previously only scalar functions have behaved in that way. sum throws away the intermediate steps in its computation, the point of cumsum is to keep them and return them - promoting them along the way means copying cumsum would return a different element type and different result than in-place cumsum! which is odd to me - many of the copying versions of functions can be thought of, and implemented, as in_place!(copy(x)) which this change would move away from here. I'm not sure cumsum is even a notable, widely-used enough function to need to preserve the export and make that judgement call within base rather than allowing users to specify exactly what behavior they want via the function they choose to pass to accumulate.

[StefanKarpinski]

which behavior is all of the code in the wild that has been using cumsum expecting,

This argument can always be made for not changing anything. An argument that can always be made in all circumstances doesn't actually help decide anything.

cumsum would return a different element type and different result than in-place cumsum! which is odd to me

This is the most compelling argument against this change so far. The choice is between cumsum being consistent with sum or cumsum!. Another consideration is that if cumsum and accumulate(+, ...) do the same thing then how does one write the often-useful widening version?

[mbauman]

cumsum would return a different element type and different result than in-place cumsum! which is odd to me

I don't find this nearly as problematic because it'll error if the values are different.

# using accumulate(add_sum, …) as a standin for cumsum since I haven't locally built this PR
julia> A = UInt8[254, 2]
       accumulate(Base.add_sum, A)
2-element Array{UInt64,1}:
 0x00000000000000fe
 0x0000000000000100

julia> accumulate!(Base.add_sum, A, A)
ERROR: InexactError: trunc(UInt8, 256)

[simonbyrne]

actually that's not true: cumsum! doesn't use add_sum (it should have the same behaviour it has now), which is why we could switch out the sparse change for cumsum!

[StefanKarpinski]

If we do this, shouldn't cumsum! use add_sum and error like @mbauman showed?

[mbauman]

We actually do have a sum! function (for reducing sums along a dimension), and it does behave as I expected:

julia> sum!(UInt8[0], UInt8[254, 2])
ERROR: InexactError: trunc(UInt8, 256)

[simonbyrne]

I figured inheriting the type of the output array would probably be the most useful behaviour.

[simonbyrne]

But maybe that would make more sense to be consistent with sum!

[JeffBezanson]

I think cumsum! should compute the same values as cumsum, and just try to store them into whatever array you supply. That way results can't silently change when you optimize code to use an in-place operation.

[tkelman]

An errors-by-default cumsum! seems way less useful than a might-overflow one.

[JeffBezanson]

Ok, but the right way to get that would be to make cumsum not widen, not make cumsum! different in an ad-hoc way.

[JeffBezanson]

Looks like this is identical to the test 2 lines above?

[JeffBezanson]

Is this unused now?

[tkelman]

the right way to get that would be to make cumsum not widen

Precisely. What real problem does making cumsum suddenly widen, with no deprecation or opt-in to that change, solve? Would be better IMO to deprecate cumsum in favor of accumulate(+, ...) now, which gives exactly the same behavior it has always had of not widening. Bringing it back in a different widening form could be done later, but again I don't see what problem that's solving - sum! is pretty obscure and not used often, and one of very few functions I can think of that errors on overflow.

[JeffBezanson]

I guess I don't have a strong opinion about whether cumsum should widen. Much more important to me is that f and f! compute the same values whenever possible. Deprecating cumsum to accumulate(+) is also kind of appealing since I like deleting things :)

[mbauman]

I just don't see a consistent rationale in the status quo. Why does sum widen? And why should cumsum behave differently?

[JeffBezanson]

Well, ok, maybe we should get rid of all implicit widening.

We could add a wideadd to go with widemul as well.

[JeffBezanson]

Which reminds me, sum widens Int8 to Int, but widen(Int8) is Int32. They should probably be the same (to the extent sum widens).

[JeffBezanson]

See also #9665.

[simonbyrne]

One other option is to reduce with checked_add in all cases?

[tkelman]

Much more important to me is that f and f! compute the same values whenever possible.

this now gives different element types for cumsum and cumsum!, so is not behaving as cumsum(x) = cumsum!(copy(x)) either.

[StefanKarpinski]

This identity always holds, however: cumsum(v) == cumsum!(copy(v)).

[tkelman]

Except that the right side now errors in cases it never did before, and the element types differ which matters quite a bit for array functions.

[kmsquire]

How about adding a method which takes the eltype or container type (including eltype) as a parameter (and not widening in the general case)? This would need to be implemented for many more functions, of course, but I think that making array widening less magical and more under user control would be a good thing.

[martinholters]

Before this PR:

julia> cumprod!(zeros(2,2), [1 2; 3 4], dims=1)
2×2 Array{Float64,2}:
 1.0  2.0
 3.0  8.0

After:

julia> cumprod!(zeros(2,2), [1 2; 3 4], dims=1)
2×2 Array{Float64,2}:
 1.0  2.0
 4.0  6.0

That's just plain wrong, is it?

EDIT: bisection specifically points at 82c8f45