I don’t understand how this is possible. I tried it many times, targmin always wins, and by quite a significant margin (1.5.2)
using BenchmarkTools
function targmin(a)
m = Inf
ind = 0
for i in 1:length(a)
if a[i] < m
m = a[i]
ind = i
end
end
ind
end
a = rand(10 ^ 8);
println(argmin(a) == targmin(a))
@btime argmin(a)
@btime targmin(a)
true
382.492 ms (1 allocation: 16 bytes)
182.026 ms (1 allocation: 16 bytes)
I can replicate this on master. It would be great to get to the bottom of this — ideally with a fix, but at least opening an issue.
Note that your code assumes that eltype(a) is comparable to Float64, since m is Inf in the first iteration. That might not necessarily be the case. You’re also assuming 1-based indexing, eachindex would be better.
That said, I think the crux of the matter is that argmin does findmin, which itself calls findminmax! (seemingly very complex function)
Additionally, argmin allows for reduction across arbitrary dimensions, your targmin always reduces over the whole array.
The core of argmin is:
It does some extra stuff to handle NaN etc and it is generic on the type of the array. But maybe it could be optimized as well (the code in question is quite old), for example, is both the m != m and ai != ai needed in every loop iteration?
Oh wow, that code is ages old!
Those checks are needed to ensure the same behaviour as min and max for NaN, I think.
Yeah, I commented on the NaN handling but it does for exampe m != m over and over in the loop, even if m hasn’t changed.
I wonder why the “change in iteration protocol” commit didn’t change to the new for loop as well and instead did it manually?
Probably because you want to peel off the first value to use as the initial minimum.
How about:
findmax(a) = _findmax(a, :)
function _findmax(a, ::Colon)
p = pairs(a)
y = @inbounds iterate(p)
if y === nothing
throw(ArgumentError("collection must be non-empty"))
end
(mi, m), s = y
isnan(m) && return m, mi
i = mi
while true
y = @inbounds iterate(p, s)
y === nothing && break
(i, ai), s = y
isnan(ai) && return ai, i
# Neither `m` nor `ai` can be NaN here so can use `<` instead of `isless`.
# Edit: not true, consider 0.0 < -0-0 == false
if m < ai
m, mi = ai, i
end
end
return m, mi
end
does anyone see any problems with that?
Edit, the isnan should be changed back to ai != ai to handle missing…
Edit Edit: actually, findmax is already broken with respect to missing
julia> findmax([missing, 1])
ERROR: TypeError: non-boolean (Missing) used in boolean context
Edit Edit Edit: My code fails to handle signed zero I think.
You would not need to bail out early at all, isless already treats NaN the way you want
while true
y = iterate(p, s)
y === nothing && break
(i, ai), s = y
if isless(m, ai)
m = ai
mi = i
end
end
julia> findmax([1, NaN, 2])
(NaN, 2)
Not sure if you should have a branch just for finding NaNs faster at all.
But early exit is nice in case you have a NaN at the beginning, or?
There are many early exits you could imagine, e.g. typemax(T) where T <: Union{Int32, Int64, Int128} etc., the NaN early exit seems to be of less use. And here luckily
julia> isless(NaN, NaN)
false
so you actually find the first NaN as promised.
Yeah, that’s true.
PS: This behaviour is also not very convincing to me:
julia> struct LargerThanNaN
end
julia> Base.isless(_, ::LargerThanNaN) = true
julia> findmax([1, NaN, LargerThanNaN()])
(NaN, 2)
IEEE specifies that NaN is supposed to propagate to indicate an error in an earlier computation. We’re not returning NaN because it’s the maximum (or minimum), but because it flags an illegal operation in a preceding computation.
Personally, since m < NaN is false, I would never want NaN as the result of maximum (or minimum for that matter), but that’s so 2008 of me. As of 2019, we should really be specifying whether we want maximum or maximumNumber to clarify our intent as to whether we want errors to propagate or be discarded.
How fun it is that Julia predates the latest floating point standard.