min(NaN,5) = NaN (why no @fastmath?)

General Usage
21

Sure, but it allocates the filtered array - unlike skip. Also, skip can be used with functions like findmax - unlike filter.

1 Like
22

Thank you. What do you think of these package-free alternatives that do not allocate:

maximum(u for u in x if !isnan(u))
extrema(u for u in x if !isnan(u))
findmax(ifelse(isnan(u),-Inf,u) for u in x)
1 Like
23

A few drawbacks:
Only works for “plain” numbers, eg not Unitful units.
Only type-stable for Float64, eg not for Float32.

24

Thanks, for the type stability this helps:

findmax(ifelse(isnan(u),typemin(eltype(x)),u) for u in x)
25

Yes, it does.
Also, when all values are NaNs and the result is undefined, this would return -infinity and silently continue, potentially leading to incorrect results further downstream. Meanwhile, skip + findmax would correctly throw an error.

27

Ok, this non-allocating function should be more robust:

function findmax_nan(x)
    all(isnan, x) && throw("Error: all NaN!")
    i = findfirst(!isnan, x)
    m, j = findmax(ifelse(isnan(u), typemin(eltype(x)), u) for u in x)
    return m, max(i,j)
end
1 Like
28

The plot thickens, that doesn’t apply for qNaN as “missing data” vs sNaN “propagate”.

See Table I. in the given link, I think we want to support what (some) standard does (most do not propagate, at least always, though Java does according to that doc), and preferentially what hardware does (I’m guessing these are most performant instructions), e.g. “Intel® AVX-512 VRAGESS”, “ARM® ARMv8.2 FMINNM / FMAXNM”, conformig to this (C) standard:

3.2 TS18661-1: sNaN Propagate but qNaN Are Missing Data […]

As a result, fmax will have the following new behavior:
fmax (1, qNaN) → 1 (5)
fmax (1, sNaN) → qNaN (6)

Seems in compliance with this standard (or is strengthening it):

3.1 C11: qNaN Are Missing Data […]

fmax (1, qNaN) → 1 (3)
fmax (1, sNaN) → undefined (4)

This standard was last reviewed and confirmed in 2022. Therefore this version remains current.

Abstract

ISO/IEC TS 18661-1:2014 extends programming language C to support binary floating-point arithmetic conforming to ISO/IEC/IEEE 60559:2011

29

Yeah, one can write individual functions like this, and even properly take care of all types and edgecases. Or just use nicely composable skip with regular aggregations (:
It’s like Base.skipmissing, but for all you nothings, nans, negatives, etc.

30

This feels simultaneously more ergonomic and more general/composable:

maximum(skip(isnan, A))

also allows you to write things like

findmax(skip(isnan, A))
sum(skip(isodd, B))
prod(skip(iszero, C))

It’d be nice if this was in Base.

1 Like
31

See Iterators.filter. For example,

julia> x = collect(Float64,1:10); x[3] = NaN;

julia> maximum(Iterators.filter(!isnan,x))
10.0
1 Like
32

That works ok, until you try this:

julia> findmax(Iterators.filter(!isnan,x))
ERROR: MethodError: no method matching keys(::Base.Iterators.Filter{ComposedFunction{typeof(!), typeof(isnan)}, Vector{Float64}})
33

I guess I consider that a weird case. Should it return the index into the original x or the index of the x with the skips removed? I’m not sure this generalizes very well. One can use Base.filter for the skips-removed case. For the skips-retained case, I guess one could use a monstrosity like

julia> argmax(first,Iterators.filter(!isnan∘first,zip(x,keys(x))))
(10.0, 10)

but I’ll admit the ergonomics leave something to be desired.

EDIT: maybe slightly nicer as

julia> argmax(last,Iterators.filter(!isnan∘last,pairs(x)))
10 => 10.0

but ultimately the same.

34

Well, I think “Base.filter should be fungible with Iterators.filter when possible” generalizes pretty well

35

That’s where skip brings benefits: like skipmissing, it’s understood to return just a view with the original x’s indexing retained but iterating skips over certain elements. This is unlike filter, which is understood to give a new collection altogether, even if Iterators.filter does so lazily.

36

We have the freedom to choose according to the application.
Fyi, some benchmarking:

x = rand(Float32[NaN,1,2,3.0], 1_000_000)
@btime findmax_nan($x)              # 1.344 ms (0 allocs: 0 bytes)
@btime argmax(first, Iterators.filter(!isnan∘first, zip($x, keys($x))))     # 3.391 ms (0 allocs: 0 bytes)
@btime findmax(skip(isnan, $x))     # 8.464 ms (0 allocs: 0 bytes)
37

Looks like Skipper.jl has some room for performance improvement @aplavin :eyes:

@rafael.guerra do you know how much of the performance differential is due to this?

38

Looks different to me, latest stable Julia 1.8:

julia> x = rand(Float32[NaN,1,2,3.0], 1_000_000);

julia> @btime findmax(skip(isnan, $x))
  5.855 ms (0 allocations: 0 bytes)

julia> @btime argmax(last,Iterators.filter(!isnan∘last,pairs($x)))
  6.854 ms (0 allocations: 0 bytes)

julia> @btime argmax(first, Iterators.filter(!isnan∘first, zip($x, keys($x))))
  7.034 ms (0 allocations: 0 bytes)
39

Sorry, I forgot to indicate the julia version:

versioninfo() 
julia> versioninfo()
Julia Version 1.9.0-beta3
Commit 24204a7344 (2023-01-18 07:20 UTC)
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: 8 × Intel(R) Core(TM) i7-1065G7 CPU @ 1.30GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-14.0.6 (ORCJIT, icelake-client)
  Threads: 8 on 8 virtual cores
40

This is what I get:

julia> @btime findmax_nan($x)
  1.179 ms (0 allocations: 0 bytes)
(3.0f0, 1)

julia> @btime argmax(first, Iterators.filter(!isnan∘first, zip($x, keys($x))))
  3.121 ms (0 allocations: 0 bytes)
(3.0f0, 1)

julia> @btime findmax(skip(isnan, $x))
  3.156 ms (0 allocations: 0 bytes)
(3.0f0, 1)
Julia 1.9.0-beta3 
julia> versioninfo()
Julia Version 1.9.0-beta3
Commit 24204a7344 (2023-01-18 07:20 UTC)
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: 20 × 12th Gen Intel(R) Core(TM) i9-12900HK
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-14.0.6 (ORCJIT, alderlake)
  Threads: 20 on 20 virtual cores

I suspect this.

Edit: This seems to be a case where the structure of the problem allows for findmax_nan to be a more optimal solution than using Iterators.filter or skip. For example, if x has less NaN values, then findmax(skip(isnan, x)) runs faster, and in the case where there are no NaNs, it runs nearly as fast (1.45ms here).

Personally I’d rather drop findmax(skip(isnan, x)) into my code than write my own findmax_nan function, unless it was a very hot loop that justified such an optimization, considering that the former pattern is less error-prone and more generally applicable. But the context of the OP (about fastmath) makes this off-topic.

2 Likes
41

Other language R:

min(NaN, 5.0)
min(5.0, NaN)
min(NaN, 5.0, na.rm=TRUE)
min(5.0, NaN, na.rm=TRUE)

The values returned are NaN, NaN, 5.0 and 5.0 for the above.

PS. I’m not saying that R is really analogous with Julia, but just answering the question as posed.

2 Likes