From 2b4910013cfd2abc8b5c3b4fea7b80e23b77708d Mon Sep 17 00:00:00 2001 From: Jan Weidner Date: Sat, 8 Dec 2018 22:44:11 +0100 Subject: [PATCH 1/6] faster maximum, minimum (#28849) --- base/reduce.jl | 17 ++++++++++++----- test/reduce.jl | 2 ++ 2 files changed, 14 insertions(+), 5 deletions(-) diff --git a/base/reduce.jl b/base/reduce.jl index 0fbcc5a8a2a45..15f0fe8b26eed 100644 --- a/base/reduce.jl +++ b/base/reduce.jl @@ -453,16 +453,23 @@ prod(a) = mapreduce(identity, mul_prod, a) ## maximum & minimum +# these propagate NaN correctly only in the second argument +_max_fast(x, y) = ifelse(x > y, x, y) +_min_fast(x, y) = ifelse(x < y, x, y) + +_fast(::typeof(min)) = _min_fast +_fast(::typeof(max)) = _max_fast + function mapreduce_impl(f, op::Union{typeof(max), typeof(min)}, A::AbstractArray, first::Int, last::Int) - # locate the first non NaN number + op_fast = _fast(op) @inbounds a1 = A[first] v = mapreduce_first(f, op, a1) - i = first + 1 - while (v == v) && (i <= last) + for i in (first+1):last + # short circuit NaN correctly + (v == v) || return v @inbounds ai = A[i] - v = op(v, f(ai)) - i += 1 + v = op_fast(v, f(ai)) end v end diff --git a/test/reduce.jl b/test/reduce.jl index f2bdd30d21684..66b1e646e795e 100644 --- a/test/reduce.jl +++ b/test/reduce.jl @@ -190,7 +190,9 @@ prod2(itr) = invoke(prod, Tuple{Any}, itr) @test isequal(extrema([NaN]), (NaN, NaN)) @test isnan(maximum([NaN, 2.])) +@test isnan(maximum([2., NaN])) @test isnan(minimum([NaN, 2.])) +@test isnan(minimum([2., NaN])) @test isequal(extrema([NaN, 2.]), (NaN,NaN)) @test isnan(maximum([NaN, 2., 3.])) From ae2f40ef7c6983302735aecdd0ad076bcecbef6a Mon Sep 17 00:00:00 2001 From: Jan Weidner Date: Sun, 9 Dec 2018 00:04:54 +0100 Subject: [PATCH 2/6] fix --- base/reduce.jl | 25 +++++++++++++------------ test/reduce.jl | 10 ++++++++++ 2 files changed, 23 insertions(+), 12 deletions(-) diff --git a/base/reduce.jl b/base/reduce.jl index 15f0fe8b26eed..ada69667ba68d 100644 --- a/base/reduce.jl +++ b/base/reduce.jl @@ -454,22 +454,23 @@ prod(a) = mapreduce(identity, mul_prod, a) ## maximum & minimum # these propagate NaN correctly only in the second argument -_max_fast(x, y) = ifelse(x > y, x, y) -_min_fast(x, y) = ifelse(x < y, x, y) - -_fast(::typeof(min)) = _min_fast -_fast(::typeof(max)) = _max_fast +_fast(::typeof(max), x, y) = ifelse(x > y, x, y) +_fast(::typeof(min), x, y) = ifelse(x < y, x, y) function mapreduce_impl(f, op::Union{typeof(max), typeof(min)}, A::AbstractArray, first::Int, last::Int) - op_fast = _fast(op) - @inbounds a1 = A[first] + a1 = @inbounds A[first] v = mapreduce_first(f, op, a1) - for i in (first+1):last - # short circuit NaN correctly - (v == v) || return v - @inbounds ai = A[i] - v = op_fast(v, f(ai)) + chunk_len = 256 + stop = first + while last - stop > 0 + v == v || return v + start = stop + 1 + stop = min(stop + chunk_len, last) + @simd for i in start:stop + @inbounds ai = A[i] + v = _fast(op, v, f(ai)) + end end v end diff --git a/test/reduce.jl b/test/reduce.jl index 66b1e646e795e..1e0c977f2b1f6 100644 --- a/test/reduce.jl +++ b/test/reduce.jl @@ -185,6 +185,16 @@ prod2(itr) = invoke(prod, Tuple{Any}, itr) @test minimum([4, 3, 5, 2]) == 2 @test extrema([4, 3, 5, 2]) == (2, 5) +@testset "maximum checks all elements" begin + N = 1025 + for i in 1:N + arr = fill(0., N) + truth = rand() + arr[i] = truth + @test maximum(arr) == truth + end +end + @test isnan(maximum([NaN])) @test isnan(minimum([NaN])) @test isequal(extrema([NaN]), (NaN, NaN)) From c1f5569a8316c37a4921cfda52f25f9b235dad47 Mon Sep 17 00:00:00 2001 From: Jan Weidner Date: Mon, 10 Dec 2018 21:25:42 +0100 Subject: [PATCH 3/6] fix --- base/reduce.jl | 16 ++++++++-------- test/reduce.jl | 10 +++++++++- 2 files changed, 17 insertions(+), 9 deletions(-) diff --git a/base/reduce.jl b/base/reduce.jl index ada69667ba68d..12eab79ab7c7b 100644 --- a/base/reduce.jl +++ b/base/reduce.jl @@ -452,21 +452,21 @@ julia> prod(1:20) prod(a) = mapreduce(identity, mul_prod, a) ## maximum & minimum - -# these propagate NaN correctly only in the second argument -_fast(::typeof(max), x, y) = ifelse(x > y, x, y) -_fast(::typeof(min), x, y) = ifelse(x < y, x, y) +_fast(::typeof(max), x, y) = ifelse(x == x, + ifelse(x > y, x, y), + x) +_fast(::typeof(min), x, y) = ifelse(x == x, + ifelse(x < y, x, y), + x) function mapreduce_impl(f, op::Union{typeof(max), typeof(min)}, A::AbstractArray, first::Int, last::Int) a1 = @inbounds A[first] v = mapreduce_first(f, op, a1) chunk_len = 256 - stop = first - while last - stop > 0 + for start in (first + 1):chunk_len:last v == v || return v - start = stop + 1 - stop = min(stop + chunk_len, last) + stop = min(start + chunk_len-1, last) @simd for i in start:stop @inbounds ai = A[i] v = _fast(op, v, f(ai)) diff --git a/test/reduce.jl b/test/reduce.jl index 1e0c977f2b1f6..f283e42db9eb2 100644 --- a/test/reduce.jl +++ b/test/reduce.jl @@ -185,13 +185,21 @@ prod2(itr) = invoke(prod, Tuple{Any}, itr) @test minimum([4, 3, 5, 2]) == 2 @test extrema([4, 3, 5, 2]) == (2, 5) -@testset "maximum checks all elements" begin +@testset "minimum/maximum checks all elements" begin N = 1025 for i in 1:N arr = fill(0., N) truth = rand() arr[i] = truth @test maximum(arr) == truth + + truth = -rand() + arr[i] = truth + @test minimum(arr) == truth + + arr[i] = NaN + @test isnan(maximum(arr)) + @test isnan(minimum(arr)) end end From 1ca74c5e41dffe3c7d034a81c902cc95eaa0b677 Mon Sep 17 00:00:00 2001 From: Jan Weidner Date: Mon, 10 Dec 2018 23:13:34 +0100 Subject: [PATCH 4/6] add max_maximum and min_minimum --- base/reduce.jl | 73 +++++++++++++++++++++++++++++++++++------------ base/reducedim.jl | 6 ++-- 2 files changed, 57 insertions(+), 22 deletions(-) diff --git a/base/reduce.jl b/base/reduce.jl index 12eab79ab7c7b..68030ff351e78 100644 --- a/base/reduce.jl +++ b/base/reduce.jl @@ -452,31 +452,66 @@ julia> prod(1:20) prod(a) = mapreduce(identity, mul_prod, a) ## maximum & minimum -_fast(::typeof(max), x, y) = ifelse(x == x, - ifelse(x > y, x, y), - x) -_fast(::typeof(min), x, y) = ifelse(x == x, - ifelse(x < y, x, y), - x) - -function mapreduce_impl(f, op::Union{typeof(max), typeof(min)}, + +""" + max_maximum(x,y) + +Compute the bigger of two elements. This function dose the same thing as `max(x,y)`, except that +* It is slightly faster +* The order of `-0.0` and `0.0` is undefined. + +See also [`max`](@ref), [`min_minimum`](@ref) +""" +function max_maximum(x,y) + ifelse(isnan(x), + x, + ifelse(x > y, x, y)) +end + +""" + min_minimum(x,y) + +Variant of [`max_maximum`](@ref) for computing minima. +""" +function min_minimum(x,y) + ifelse(isnan(x), + x, + ifelse(x < y, x, y)) +end + +function mapreduce_impl(f, op::Union{typeof(max_maximum), typeof(min_minimum)}, A::AbstractArray, first::Int, last::Int) a1 = @inbounds A[first] - v = mapreduce_first(f, op, a1) + v1 = mapreduce_first(f, op, a1) + v2 = v3 = v4 = v1 chunk_len = 256 - for start in (first + 1):chunk_len:last - v == v || return v - stop = min(start + chunk_len-1, last) - @simd for i in start:stop - @inbounds ai = A[i] - v = _fast(op, v, f(ai)) + start = first + stop = start + chunk_len - 4 + while stop <= last + isnan(v1) && return v1 + isnan(v2) && return v2 + isnan(v3) && return v3 + isnan(v4) && return v4 + @inbounds for i in start:4:stop + v1 = op(v1, f(A[i+1])) + v2 = op(v2, f(A[i+2])) + v3 = op(v3, f(A[i+3])) + v4 = op(v4, f(A[i+4])) end + start = stop + stop = start + chunk_len - 4 + end + v = op(op(v1,v2),op(v3,v4)) + start += 1 + for i in start:last + @inbounds ai = A[i] + v = op(v, f(A[i])) end v end -maximum(f, a) = mapreduce(f, max, a) -minimum(f, a) = mapreduce(f, min, a) +maximum(f, a) = mapreduce(f, max_maximum, a) +minimum(f, a) = mapreduce(f, min_minimum, a) """ maximum(itr) @@ -492,7 +527,7 @@ julia> maximum([1,2,3]) 3 ``` """ -maximum(a) = mapreduce(identity, max, a) +maximum(a) = mapreduce(identity, max_maximum, a) """ minimum(itr) @@ -508,7 +543,7 @@ julia> minimum([1,2,3]) 1 ``` """ -minimum(a) = mapreduce(identity, min, a) +minimum(a) = mapreduce(identity, min_minimum, a) ## all & any diff --git a/base/reducedim.jl b/base/reducedim.jl index 3a3f9f14469c8..04dde5e298172 100644 --- a/base/reducedim.jl +++ b/base/reducedim.jl @@ -125,7 +125,7 @@ function _reducedim_init(f, op, fv, fop, A, region) end # initialization when computing minima and maxima requires a little care -for (f1, f2, initval) in ((:min, :max, :Inf), (:max, :min, :(-Inf))) +for (f1, f2, initval) in ((:min_minimum, :max_maximum, :Inf), (:max_maximum, :min_minimum, :(-Inf))) @eval function reducedim_init(f, op::typeof($f1), A::AbstractArray, region) # First compute the reduce indices. This will throw an ArgumentError # if any region is invalid @@ -642,7 +642,7 @@ julia> any!([1 1], A) any!(r, A) for (fname, _fname, op) in [(:sum, :_sum, :add_sum), (:prod, :_prod, :mul_prod), - (:maximum, :_maximum, :max), (:minimum, :_minimum, :min)] + (:maximum, :_maximum, :max_maximum), (:minimum, :_minimum, :min_minimum)] @eval begin # User-facing methods with keyword arguments @inline ($fname)(a::AbstractArray; dims=:) = ($_fname)(a, dims) @@ -662,7 +662,7 @@ all(f::Function, a::AbstractArray; dims=:) = _all(f, a, dims) _all(a, ::Colon) = _all(identity, a, :) for (fname, op) in [(:sum, :add_sum), (:prod, :mul_prod), - (:maximum, :max), (:minimum, :min), + (:maximum, :max_maximum), (:minimum, :min_minimum), (:all, :&), (:any, :|)] fname! = Symbol(fname, '!') _fname = Symbol('_', fname) From d84337fa113624a97d86c4090304bd809ac8daef Mon Sep 17 00:00:00 2001 From: Jan Weidner Date: Tue, 11 Dec 2018 22:51:23 +0100 Subject: [PATCH 5/6] fix --- base/reduce.jl | 54 +++++++++++++++++++++++------------------------ base/reducedim.jl | 6 +++--- test/reduce.jl | 51 ++++++++++++++++++++++++++++++++------------ 3 files changed, 67 insertions(+), 44 deletions(-) diff --git a/base/reduce.jl b/base/reduce.jl index 68030ff351e78..2bd506f7eda3f 100644 --- a/base/reduce.jl +++ b/base/reduce.jl @@ -452,34 +452,24 @@ julia> prod(1:20) prod(a) = mapreduce(identity, mul_prod, a) ## maximum & minimum - -""" - max_maximum(x,y) - -Compute the bigger of two elements. This function dose the same thing as `max(x,y)`, except that -* It is slightly faster -* The order of `-0.0` and `0.0` is undefined. - -See also [`max`](@ref), [`min_minimum`](@ref) -""" -function max_maximum(x,y) +function _fast(::typeof(max),x,y) ifelse(isnan(x), x, ifelse(x > y, x, y)) end -""" - min_minimum(x,y) - -Variant of [`max_maximum`](@ref) for computing minima. -""" -function min_minimum(x,y) +function _fast(::typeof(min),x,y) ifelse(isnan(x), x, ifelse(x < y, x, y)) end -function mapreduce_impl(f, op::Union{typeof(max_maximum), typeof(min_minimum)}, +isbadzero(::typeof(max), x) = (x == zero(x)) & signbit(x) +isbadzero(::typeof(min), x) = (x == zero(x)) & !signbit(x) +isgoodzero(::typeof(max), x) = isbadzero(min, x) +isgoodzero(::typeof(min), x) = isbadzero(max, x) + +function mapreduce_impl(f, op::Union{typeof(max), typeof(min)}, A::AbstractArray, first::Int, last::Int) a1 = @inbounds A[first] v1 = mapreduce_first(f, op, a1) @@ -493,10 +483,10 @@ function mapreduce_impl(f, op::Union{typeof(max_maximum), typeof(min_minimum)}, isnan(v3) && return v3 isnan(v4) && return v4 @inbounds for i in start:4:stop - v1 = op(v1, f(A[i+1])) - v2 = op(v2, f(A[i+2])) - v3 = op(v3, f(A[i+3])) - v4 = op(v4, f(A[i+4])) + v1 = _fast(op, v1, f(A[i+1])) + v2 = _fast(op, v2, f(A[i+2])) + v3 = _fast(op, v3, f(A[i+3])) + v4 = _fast(op, v4, f(A[i+4])) end start = stop stop = start + chunk_len - 4 @@ -507,11 +497,21 @@ function mapreduce_impl(f, op::Union{typeof(max_maximum), typeof(min_minimum)}, @inbounds ai = A[i] v = op(v, f(A[i])) end - v + + # enforce correct order of 0.0 and -0.0 + # e.g. maximum([0.0, -0.0]) === 0.0 + # should hold + if isbadzero(op, v) + for i in first:last + x = @inbounds A[i] + isgoodzero(op,x) && return x + end + end + return v end -maximum(f, a) = mapreduce(f, max_maximum, a) -minimum(f, a) = mapreduce(f, min_minimum, a) +maximum(f, a) = mapreduce(f, max, a) +minimum(f, a) = mapreduce(f, min, a) """ maximum(itr) @@ -527,7 +527,7 @@ julia> maximum([1,2,3]) 3 ``` """ -maximum(a) = mapreduce(identity, max_maximum, a) +maximum(a) = mapreduce(identity, max, a) """ minimum(itr) @@ -543,7 +543,7 @@ julia> minimum([1,2,3]) 1 ``` """ -minimum(a) = mapreduce(identity, min_minimum, a) +minimum(a) = mapreduce(identity, min, a) ## all & any diff --git a/base/reducedim.jl b/base/reducedim.jl index 04dde5e298172..3a3f9f14469c8 100644 --- a/base/reducedim.jl +++ b/base/reducedim.jl @@ -125,7 +125,7 @@ function _reducedim_init(f, op, fv, fop, A, region) end # initialization when computing minima and maxima requires a little care -for (f1, f2, initval) in ((:min_minimum, :max_maximum, :Inf), (:max_maximum, :min_minimum, :(-Inf))) +for (f1, f2, initval) in ((:min, :max, :Inf), (:max, :min, :(-Inf))) @eval function reducedim_init(f, op::typeof($f1), A::AbstractArray, region) # First compute the reduce indices. This will throw an ArgumentError # if any region is invalid @@ -642,7 +642,7 @@ julia> any!([1 1], A) any!(r, A) for (fname, _fname, op) in [(:sum, :_sum, :add_sum), (:prod, :_prod, :mul_prod), - (:maximum, :_maximum, :max_maximum), (:minimum, :_minimum, :min_minimum)] + (:maximum, :_maximum, :max), (:minimum, :_minimum, :min)] @eval begin # User-facing methods with keyword arguments @inline ($fname)(a::AbstractArray; dims=:) = ($_fname)(a, dims) @@ -662,7 +662,7 @@ all(f::Function, a::AbstractArray; dims=:) = _all(f, a, dims) _all(a, ::Colon) = _all(identity, a, :) for (fname, op) in [(:sum, :add_sum), (:prod, :mul_prod), - (:maximum, :max_maximum), (:minimum, :min_minimum), + (:maximum, :max), (:minimum, :min), (:all, :&), (:any, :|)] fname! = Symbol(fname, '!') _fname = Symbol('_', fname) diff --git a/test/reduce.jl b/test/reduce.jl index f283e42db9eb2..317fb814fb280 100644 --- a/test/reduce.jl +++ b/test/reduce.jl @@ -185,21 +185,44 @@ prod2(itr) = invoke(prod, Tuple{Any}, itr) @test minimum([4, 3, 5, 2]) == 2 @test extrema([4, 3, 5, 2]) == (2, 5) +@test maximum([-0.,0.]) === 0.0 +@test maximum([0.,-0.]) === 0.0 +@test maximum([0.,-0.,0.]) === 0.0 +@test minimum([-0.,0.]) === -0.0 +@test minimum([0.,-0.]) === -0.0 +@test minimum([0.,-0.,0.]) === -0.0 + @testset "minimum/maximum checks all elements" begin - N = 1025 - for i in 1:N - arr = fill(0., N) - truth = rand() - arr[i] = truth - @test maximum(arr) == truth - - truth = -rand() - arr[i] = truth - @test minimum(arr) == truth - - arr[i] = NaN - @test isnan(maximum(arr)) - @test isnan(minimum(arr)) + for N in [2:20;150;300] + for i in 1:N + arr = fill(0., N) + truth = rand() + arr[i] = truth + @test maximum(arr) == truth + + truth = -rand() + arr[i] = truth + @test minimum(arr) == truth + + arr[i] = NaN + @test isnan(maximum(arr)) + @test isnan(minimum(arr)) + + arr = zeros(N) + @test minimum(arr) === 0.0 + @test maximum(arr) === 0.0 + + arr[i] = -0.0 + @test minimum(arr) === -0.0 + @test maximum(arr) === 0.0 + + arr = -zeros(N) + @test minimum(arr) === -0.0 + @test maximum(arr) === -0.0 + arr[i] = 0.0 + @test minimum(arr) === -0.0 + @test maximum(arr) === 0.0 + end end end From afbd564a3028b4e51180e49baca582676358f478 Mon Sep 17 00:00:00 2001 From: Jan Weidner Date: Wed, 12 Dec 2018 13:23:55 +0100 Subject: [PATCH 6/6] fix --- base/reduce.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/base/reduce.jl b/base/reduce.jl index 2bd506f7eda3f..8717181e80e39 100644 --- a/base/reduce.jl +++ b/base/reduce.jl @@ -464,8 +464,9 @@ function _fast(::typeof(min),x,y) ifelse(x < y, x, y)) end -isbadzero(::typeof(max), x) = (x == zero(x)) & signbit(x) -isbadzero(::typeof(min), x) = (x == zero(x)) & !signbit(x) +isbadzero(::typeof(max), x::AbstractFloat) = (x == zero(x)) & signbit(x) +isbadzero(::typeof(min), x::AbstractFloat) = (x == zero(x)) & !signbit(x) +isbadzero(op, x) = false isgoodzero(::typeof(max), x) = isbadzero(min, x) isgoodzero(::typeof(min), x) = isbadzero(max, x)