Allow exploiting mutability in `sum` et al

[mgkurtz]

Currently the implementations in reduce.jl which lie behind sum, prod et al do not really care for mutable input.

This means two things:

First, sum([a]) === a, and since sum(array) is never identical to an entry for longer arrays, this speciality may escape people’s notice, so they may modify the result of sum(array) leading to strange bugs in one element arrays. (On the other hand, probably, minimum(array) shall remain identical to an entry of array.)

Second, calling reduce(add!, array) modifies array[1], or if the execution goes into the mapreduce_impl branch, this will even modify several entries. This is bad. But using a modifying addition is rather useful. If, for example array is a Vector{BigInt} using Base.GMP.MPZ.add! is ten times faster than using +. Similar speedups can be achieved for various computations in the computer algebra system OSCAR, where we have more than one hundred types with an in place add! method, that could all benefit of allowing modifying operators in reduce.

So, I would be rather glad, if we could add

reduce_first(op, x) = deepcopy(x)
reduce_first(::Union{typeof(min), typeof(max)}, x) = x

always call reduce_first in the appropriate places in reduce.jl, add a note about mutable operators in reducing functions and officially document add_sum and mul_prod.

Then, I could overwrite

Base.add_sum(a::T, b::T) where T<:SuperTypeOfMutableAdditiveTypes = add!(a, a, b)
Base.mul_prod(a::T, b::T) where T<:SuperTypeOfMutableMultiplicativeTypes = mul!(a, a, b)

to get much faster sums and products for all appropriate types.

Of course, we could also choose reduce_first to only deepcopy for + and * instead of for all but min and max.

By the way, I know, there is a MutableArithmetic package, but as it can hardly hijack the sum implementation, this needs the MutableArithmetic.@rewrite macro when using, and quite more complicated setup than overwriting add_sum as above. (Or, I do not understand the docs correctly.) Especially as I want sum(f(x) for x in array) to work, where type dispatching gets complicated.

As allowing mutable data and mutating operators in reduce would allow for a big speedup in some use cases, increase consistency in some others and seems to be little effort, I hope, we can manage to explicitly support them. If we can settle to a solution in this issue, I will gladly provide a pull request.

[andrewjradcliffe]

Mutability in reductions can be surprisingly complicated, in particular due to the variability of the beholder's conception of where mutability should occur/not occur. The base implementations attempt to be as general as possible, so as to avoid imposing limitations (e.g. performance) and constraints (if I want mutability then forcing a deepcopy makes that impossible).

Largely, it is up to the programmer to provide an init for cases in which a neutral mutable element is desired as the basis for the computation. One notes that the functions ending in ! do so for a good reason -- these as "beware! I mutate and demand appropriate action!". With freedom and generality in a language, the programmer must think a little more about the what is "appropriate"; this is a price worth paying.

where we have more than one hundred types with an in place add! method, that could all benefit of allowing modifying operators in reduce.

It is worth noting that sum and prod already respect the special nature of BigFloat and BigInt. More generally, this sounds like a design problem that necessitates detailed thought. Type piracy (w.r.t. base) is never a good option, as it means that your package would introduce unexpected behavior. However, modifying base to support special cases is a slippery slope. There is no escaping the fact that programming problems demand the same careful consideration afforded to mathematical problems (as they are mathematical problems in disguise).

Then, I could overwrite [...] Base.add_sum, Base.mul_prod

Doing so would deviate from their current purpose. The further addition of mutation would lead to unexpected results (not to mention problems throughout the package ecosystem which rely on current behavior). If a mutating + or * is truly desired, one always has the option to define +ꜝ(a, b) = Base.GMP.MPZ.add!(a, b) (adding definitions as one needs).

On a practical note, the examples based on reduce can and should be dealt with by providing a neutral initial element via the init keyword. This is how one can utilize reduce with mutating functions to achieve both the "intended" behavior and high performance. e.g.

julia> v = big.(rand(Int128, 2^10));

julia> @allocated sum(v)
64

julia> @allocated reduce(+, v)
64496

julia> @allocated reduce(Base.GMP.MPZ.add!, v, init=zero(eltype(v)))
96

julia> a = [1,2,3];

julia> A = [a];

julia> reduce(+, A, init=zero(first(A))) == a
true

julia> reduce(+, A, init=zero(first(A))) === a
false

[StefanKarpinski]

I fully agree with everything @andrewjradcliffe wrote. Automatically doing deepcopy is pretty much never a good idea. More generally, numbers should always act like they are immutable values—we only allow mutating them as a dirty concession to the practicalities of BigInts. So "mutating arithmetic" operations like add! shouldn't really exist in the first place, let alone providing motivation for doing automatic deepcopy calls in fundamental functional operations like reduce. It's possible that there's some improvement to be made here, but it's definitely not that. The suggestion of always using the init form of reduction when the reducer does mutation seems like the correct one.

[mgkurtz]

Perhaps, I should have pointed out, that in the end, I want sum and prod to work quickly. So yes, perhaps the improvement shall be done at another place. As sum and prod both directly delegate to mapreduce and add_sum/mul_prod changing them seemed sensible to me. Anyway, I learn, that I should name my issues by what I really want. I ask for forgiveness.

More generally, numbers should always act like they are immutable values

I totally agree with you, that they should. I am an algebraist after all. I love functional programming with all its constraints, immutability and clean code. Only, if we write our entire code in a purely functional way, without ever mutating anything (as most of the things in a computer algebra library are numbers in some sense), the code will end up much slower and we could just as well stick to existing alternatives using Python et al. So not mutating our “numbers“ is sadly no practical option for us.

we only allow mutating them as a dirty concession to the practicalities of BigInts.

Actually summing BigInts is only rudimentary implemented, only the first call is efficient, the others take ten times as long for summing.

as = rand(BigInt, 1000, 1000)
f(x::Int) :: BigInt

sum(as)
sum(as; init=big(0))
sum(as; dims=1)
sum(f, 1:1_000_000)
sum(f(x*x) for x in 1:1_000_000)

Here, BigInt is just a place holder for all kind of mutable “numbers“ which are a lot in OSCAR and in base, this can be arrays (think of a generator yielding vectors or matrices) for example.

modifying base to support special cases is a slippery slope.

I agree here, too. Then we need to answer whether sum is intended for numbers of bitstype only. So, is it a mere side effect of Julia’s ducktyping that we can sum up BigInts and arrays? Meaning the halfhearted implementation of summing BigInts is kind of an accident? I have the feeling @StefanKarpinski argues in that direction. From a technical point of view this actually seems natural to me, from a high level programmer point of view less so.

The suggestion of always using the init form of reduction when the reducer does mutation seems like the correct one.

True enough. (Although the documentation of reduce does not guarantee that, so I should use foldl instead.) I should have stated more clearly, that I actually want sum to work as I want.