Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Improve quantile performance #86

Open
wants to merge 4 commits into
base: master
Choose a base branch
from
Open
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
23 changes: 17 additions & 6 deletions src/Statistics.jl
Original file line number Diff line number Diff line change
Expand Up @@ -937,22 +937,33 @@ function quantile!(q::AbstractArray, v::AbstractVector, p::AbstractArray;
end
isempty(q) && return q

minp, maxp = extrema(p)
_quantilesort!(v, sorted, minp, maxp)
if length(p) == 2
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Special-casing 2 is kind of weird. For example, for [0.25, 0.5, 0.75] this branch would probably also be faster, right? Actually, isn't this approach faster than the other in most cases?

A possible optimization would be to call partialsort! on the full array for the first quantile, then call it on a view from the first quantile to the end of the array for the second, and so on. Not sure that would make a big difference, but at least it shouldn't be really slower than sorting everything between the two extreme quantiles, right?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, it is faster:

julia> using Statistics

julia> using BenchmarkTools

julia> x = rand(10^4);

julia> @btime quantile($x, [0.25, 0.5, 0.75]);
  410.400 μs (4 allocations: 78.42 KiB)

julia> @btime [quantile($x, 0.25), quantile($x, 0.5), quantile($x, 0.75)];
  314.000 μs (7 allocations: 234.72 KiB)

The problem is that quantiles do not need to be sorted, so this complicates the code (but of course is doable as I guess sorting requested quantiles should not be problematic for performance).

We do not need views I think, as sort! supports passing start and stop ranges.

However, as it seems to be a bigger rework I will leave it for after DataFrames.jl 1.3 releaes when we work on general Statistics update.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

OK. Yes, sorting quantiles should have a negligible cost (and most of the time they will already be sorted).

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do you feel like finishing this?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I have opened https://github.com/JuliaLang/Statistics.jl/pull/91 to perform easy comparison of both.

for (i, j) in zip(eachindex(p), eachindex(q))
@inbounds q[j] = quantile!(v, p[i], sorted=sorted, alpha=alpha, beta=beta)
end
else
minp, maxp = extrema(p)
_quantilesort!(v, sorted, minp, maxp)

for (i, j) in zip(eachindex(p), eachindex(q))
@inbounds q[j] = _quantile(v,p[i], alpha=alpha, beta=beta)
for (i, j) in zip(eachindex(p), eachindex(q))
@inbounds q[j] = _quantile(v, p[i], alpha=alpha, beta=beta)
end
end

return q
end

function quantile!(v::AbstractVector, p::Union{AbstractArray, Tuple{Vararg{Real}}};
sorted::Bool=false, alpha::Real=1., beta::Real=alpha)
if length(p) == 2
return map(x -> quantile!(v, x, sorted=sorted, alpha=alpha, beta=beta), p)
end

if !isempty(p)
minp, maxp = extrema(p)
_quantilesort!(v, sorted, minp, maxp)
end
return map(x->_quantile(v, x, alpha=alpha, beta=beta), p)
return map(x -> _quantile(v, x, alpha=alpha, beta=beta), p)
end

quantile!(v::AbstractVector, p::Real; sorted::Bool=false, alpha::Real=1., beta::Real=alpha) =
Expand All @@ -965,10 +976,10 @@ function _quantilesort!(v::AbstractArray, sorted::Bool, minp::Real, maxp::Real)

if !sorted
lv = length(v)
# only need to perform partial sort
lo = floor(Int,minp*(lv))
hi = ceil(Int,1+maxp*(lv))

# only need to perform partial sort
sort!(v, 1, lv, Base.Sort.PartialQuickSort(lo:hi), Base.Sort.Forward)
end
if (sorted && (ismissing(v[end]) || (v[end] isa Number && isnan(v[end])))) ||
Expand Down