maximum([fill(NaN, 255); missing]) !== maximum([fill(NaN, 256); missing]) (Should we fix reduce(max, ...) etc.?)

[tkf]

Currently, we have

julia> maximum([fill(NaN, 256); missing])
NaN

julia> maximum([fill(NaN, 255); missing])
missing

julia> maximum(x for x in [fill(NaN, 256); missing])
missing

julia> maximum(x for x in [fill(NaN, 255); missing])
missing

julia> max(missing, NaN)
missing

As you can see, maximum is incompatible with the definition of max depending on the size of the array. Furthermore, the result changes when processed with the identity iterator transformation.

I think maximum(xs) should be defined as reduce(max, xs) semantically. Currently, it is defined as such at the "syntax" level:

julia/base/reducedim.jl

Lines 716 to 725 in 13b07fc

	for (fname, _fname, op) in [(:sum, :_sum, :add_sum), (:prod, :_prod, :mul_prod),
	(:maximum, :_maximum, :max), (:minimum, :_minimum, :min)]
	@eval begin
	# User-facing methods with keyword arguments
	@inline ($fname)(a::AbstractArray; dims=:, kw...) = ($_fname)(a, dims; kw...)
	@inline ($fname)(f, a::AbstractArray; dims=:, kw...) = ($_fname)(f, a, dims; kw...)
	# Underlying implementations using dispatch
	($_fname)(a, ::Colon; kw...) = ($_fname)(identity, a, :; kw...)
	($_fname)(f, a, ::Colon; kw...) = mapreduce(f, $op, a; kw...)

However, mapreduce implementation for max/min is not compatible with the mathematical definition of mapreduce.

There were related discussions on findmax/min and argmax/min (#27613, #35316) and I've been suggesting to define the total orderings compatible with max/min and use it to define all the similar functions, at least at the API definition level. I think this is the right direction as otherwise parallelizing these reducers is hard; i.e., we need to express them as folds with an associative operator which requires transitivity in the corresponding comparison. Furthermore, rigorously defining these reducers would let us confidently leverage the commutativity of max/min and define efficient maximum(::Diagonal) etc. #30236.

So, I suggest to:

1. Use the mathematically correct definition in mapreduce(min, ...) and mapreduce(max, ...).
2. Define the total ordering compatible with min (maybe in Add 2-arg versions of findmax/min, argmax/min #35316) and use it in argmin and findmin. (Note: isless is already max-compatible)

cc @nalimilan @StefanKarpinski @mbauman @JeffBezanson

[nalimilan]

Good catch. I don't think this qualifies as a minor change though, it's just a bug that should be fixed. maximum([missing; fill(NaN, 256)]), which is very close to one of your examples, even errors in 1.4. The fix I implemented at #35989 should be improved to avoid short-circuiting for NaN when the array can contain missing.

I can't comment on your more general proposal. It sounds appealing but I'm not able to judge. :-)

[tkf]

@nalimilan What do you think about

julia> findmax([NaN, missing])
(NaN, 1)

being incompatible with maximum? That's the minor change part that I wanted to fix; i.e., return (missing, 2) (which is not clear at all from reading the OP, now that I look at it again).

[nalimilan]

Well I'd say it qualifies as a bug fix, just like changing maximum([fill(NaN, 256); missing]) to return missing. Anyway saying it's a minor change is OK too, what matters is not to wait until 2.0. :-)

[tkf]

OK. findmax explicitly says "If any data element is NaN" (even though it directly contradicts with "The result is in line with max.") so I was a bit worried:

julia/base/array.jl

Lines 2172 to 2177 in 13b07fc

	findmax(itr) -> (x, index)
	Return the maximum element of the collection `itr` and its index. If there are multiple
	maximal elements, then the first one will be returned.
	If any data element is `NaN`, this element is returned.
	The result is in line with `max`.

[nalimilan]

Woops. I wasn't aware of that. So yeah that would be a minor change.

[tkf]

I should have mentioned it :)

[ViralBShah]

This specific issue is fixed in 1.9:

julia> maximum([fill(NaN, 256); missing])
missing

julia> maximum([fill(NaN, 255); missing])
missing

julia> maximum(x for x in [fill(NaN, 256); missing])
missing

julia> maximum(x for x in [fill(NaN, 255); missing])
missing

julia> max(missing, NaN)
missing