revising sorting APIs (2.0)

[StefanKarpinski]

@LilithHafner, I'd love to get your take on revising the sorting APIs for Julia 2.0 while you're still very fresh in the weeds with it. It was originally designed before keyword arguments were fast and even before higher order functions were fast, which is why it leans so heavily on constructing ordering types and sort algorithm types and passing those to the core sorting logic and dispatching on them since that has basically always been fast. At some point pre-1.0 we considered changing the Ordering stuff, but we decided (I believe it was @JeffBezanson and myself, probably on a rushed triage call—the triage leading up to 1.0 was an intense time) that it was kind of useful to have a first class representation of a sorting order. So I'd love your take on whether you feel that's true or not. I've also often considered it a bit of a mistake that we expose all the gnarly internal methods of sort! externally and it seems like we should probably have _sort! or sortcore! with all the internal methods and then only have the official API methods for sort! and sort.

[LilithHafner]

Yes. Only one or two methods of sort!, sort, sortperm, etc. should be exported, and most should be internal. I like _sort.

I have mixed feelings about orderings. It's nice/necessary to have a way to make keyword arguments disappear in the implementation of e.g QuickSort, but supporting orderings is redundant with supporting AbstractVectors of arbitrary eltype + multiple dispatch, which we already have. If we could somehow attach the ordering to the array without sacrificing performance, then we could eliminate the need for

julia/base/sort.jl

Line 750 in badad9d

function _extrema(v::AbstractArray, lo::Integer, hi::Integer, o::Ordering)

And let our sorting algorithms simply call isless.

A toy (broken) implementation using reinterpret:

struct WithIsless{T, L}
    data::T
    isless::L
end
Base.isless(x::T, y::T) where {T <: WithIsless} = x.isless(x.data, y.data)
function sort!(x, lt=identity)
    _sort!(reinterpret(WithIsless{eltype(x), typeof(lt)}, x))
    x
end

[StefanKarpinski]

I'm not sure I see the advantage of that over what we do now. Can you elaborate?

[LilithHafner]

One advantage (and it's not particularly large) is that sorting algorithms could be written without any reference to a specific ordering and that all functions (e.g. Statistics.quantiles) would support arbitrary orderings by default if you pass in an AbstractVector with an attached ordering. This way one could call StatsBase.quantiles(ReveseOrdering(rand(100)), [0, .5, 1]) and get [.991023808, 0.5091720398, 0.01239874].

I think sorting internals can be substantially improved, though I'm not always sure exactly how, but for the most part, the API seems pretty good. One gap is providing an equivalent to R's order that sorts and produces a permutation at the same time (see #939 (comment) and subsequent discussion), but that need not be breaking.

[LilithHafner]

I think we could have a (mostly) linear pipeline that begins with any exported functions (sort!, partialsortperm, etc.) which all call

_sort!(x::AbstractArray, alg::Algorithm, phase::Phase, order::Ordering; kw...)

The linear pipeline is

1. Public interface and dynamic dispatch
2. constructing a new AbstractArray or Ordering to handle sortperm! and sortperm
3. copymutable to handle sort
4. sorting an AbstractArray on a dim
5. Dispatch to insertion sort for small vectors (1)
6. Optimizations for unions with Missing
7. Optimizations for floating points
8. Dispatch to insertion sort for small vectors (2)
9. Presorted check
10. Reverse sorted check
11. Dispatch to QuickSort when radix-sort is not viable (partialsort, sortperm, or no uint-mapping; in the future this may exclusively be for when there is no uint-mappping)
12. Compute extrema (store the results in kw...)
13. Dispatch to countsort for short int range
14. Dispatch to insertion sort/quick sort for small vectors with many bits
15. Radix sort

Which can be implemented like so

function _sort!(x::AbstractArray, alg::Algorithm, ::Copy, o::Order; maybecopy, kw...)
    _sort!(maybecopy(x), alg, ArrayDims(), kw...)
end

function _sort!(x::AbstractArray, alg::Algorithm, phase::ArrayDims, o::Order; dims, kw...)
    for part in slices_of(x, dims)
        _sort!(part, alg, phase, o; kw...)
    end
end
function _sort!(x::AbstractVector, alg::Algorithm, ::ArrayDims, o::Order; kw...)
    _sort!(x, alg, ::SmallCheck1(), o; kw...)
end

function _sort!(x::AbstractVector, alg::Algorithm, ::SmallCheck1, o::Order; kw...)
    length(x) < 10 && return _sort!(x, alg, Small(), o; kw...)
    _sort!(x, alg, PartitionMissings(), o; kw...)
end

function _sort!(x::AbstractVector, alg::Algorithm, ::PartitionMissings, o::Order; kw...)
    _sort!(try_to_do_missing_optimization!(x, o), alg, PartitionFloats(), o; kw...)
end

function _sort!(x::AbstractVector, alg::Algorithm, ::PartitionFloats, o::Order; kw...)
    _sort!(try_to_do_fpsort!(x, o), alg, SmallCheck2(), o; kw...)
end

...

This is extensible with

struct SmallCheckBetweenMissingAndFloat <: Algorithm end
function Base._sort!(x::AbstractVector, alg::SmallCheckBetweenMissingAndFloat, ::PartitionFloats, o::Order; kw...)
    length(x) < 15 && return _sort!(x, alg, Small(), o, kw...)
    _sort!(x, Default(), PartitionFloats(), o; kw...)
end

Linearization maximizes code reuse.

Extensibility might make it easier to turn on algorithms and optimizations with high-level dependencies (e.g. view, rand, OffsetVector) after stdlibs are loaded.

[LilithHafner]

(This is all 1.x)

[LilithHafner]

Here is a road map to get there (not totally on topic, but I want to write it down and this seems like a reasonable place)

• Revise the workspace/buffer system
• Probably rename workspace to buffer pending conversation
• Upgrade quicksort (Stabilize, optimize, and increase robustness of QuickSort #45222)
• Implement the pipeline system above, and fix performance issues with sortperm. (Refactor and document sorting dispatch #47383)
• Move as much as possible out of base and into a stdlib (Move sorting from base to stdlib #46679)

[LilithHafner]

Some links to keep track of relevant to API decisions:

• A user story about thinking sort(x; dims) would behave as sortslices(x; dims) does: https://discourse.julialang.org/t/cant-get-sortperm-to-return-the-sorting-order-of-columns-of-a-matrix/84081/6?u=lilith
• My thoughts on sortslicesperm: Adding perm version of sortslices #45787 (comment)

I like the orthogonality of sort and eachslice and I think they have the potential to combine to cover a lot of what sortslices and dims do now.

[atteson]

I just underwent an effort to write a custom sortperm! to implement a counting sort which backs off to MergeSort when there are too many unique elements. The reason for choosing this function is because I'm doing multi-key sorts so using sortperm! with initialized=true and a stable underlying sort can faciliate this. Unfortunately, since sortperm! doesn't dispatch on alg, I found this impossible and reverted to calling my function something else. My two cents is that all the useful functions that sort.jl exports should dispatch on alg so as to facilitate user extension.

Feel free to clarify if I'm missing something here.

[atteson]

In my real use case, I'm using counting sort for elements of type Date and Int, at the moment, but would be interested in having it available for all types.

I like the idea of adding that layer of indirection that you mention.

Thank you for taking the time to reply. I have 2 remaining questions:

1. I noticed some specialized sort being used for floats in base/sort.jl though didn't spend enough time to understand it fully. Is that sort in-place? I know that might sound like a strange request but in some data, e.g. prices, rounded float point numbers are common.

2. I would be interested if you think there is currently support for sorting by multiple keys, i.e. finding a permutation that sorts by vector v1 and sorts equal elements of v1 by vector v2.

[LilithHafner]

1. It is a way of sorting floating point numbers as if they were integers. It is in place if and only if the underlying algorithm is in place. (incidentally, it is also stable if and only if the underlying algorithm is stable).

2. sort!(x, by = x -> (get_key_1(x), get_key_2(x))) is the most straightforward way, but not always the most efficient. You might be better off writing your own comparison function if get_key_2 is expensive.

[atteson]

Yes. I was using the term in-place when I meant stable. Thank you for correcting me.

ZippedArrays looks promising. Thanks for mentioning it.

[LilithHafner]

Also thanks @emmt for writing ZippedArrays.

[atteson]

FYI, in the end, my counting sort was 3x faster on the real data that I'm using at the moment than using julia sort with a ZippedArray. I'm going to stick with my counting sort for the moment.

[LilithHafner]

Should we make the buffer keyword argument a documented part of the API or strictly internal? Right now it is undocumented and therefore implicitly internal and I am okay with that because whether a pre-allocated workspace is helpful depends heavily on implementation details. (c.f. #45330 (comment))

[LilithHafner]

2.0 proposal to replace all keyword arguments with first-class Orders (sounds appealing to me)

[atteson]

I have a few more suggestions:

1. I find the "merge" operation, that is merging two sorted vectors, useful outside of merge sort. In fact, I've written this routine too many times myself at too many places. If it's possible to separate out this function without impacting the performance of merge sort, I'd vote for it.

2. Personally, I like to carry around the fact that an object is sorted in the type so that it can be dispatched off of. I know this one will be a broader change. This would mean that functions like the merge operation mentioned above and searchsorted could be implemented to only work on sorted vectors. searchsorted currently gives incorrect results if the vector is not sorted. It could be changed to only work on sorted vectors or changed to "search" and dispatch to linear search for non-sorted vectors. Let me know if there is a better place to put this suggestion such as somewhere on julia discourse.

[LilithHafner]

1. That's totally possible, but I'm not interested in doing it in base. At the moment, I don't like mergesort because quicksort is usually faster and just as stable. Unless performance characteristics change before 2.0, I'd like to consider removing MergeSort entirely. If you've written it many times, perhaps consider making your own micro-package or a PR to SortingAlgorithms.jl. I'd be happy to see that method in SortingAlgorithms.jl.

1.5) Not what you asked for, but QuickSort's partition method is also quite useful on its own. I think fpsort!'s subroutines may benefit from sharing a partition! method with QuickSort.

2. I don't think that is practical with the current type system, but with abstract multiple inheritance (abstract multiple inheritance #5 or similar) it might be possible to have x::Vector{Int} and x::Sorted and dispatch on search(x::Intersection{Sorted, AbstractVector}) rather than searchsorted(x::AbstractVector). Anything is possible in 2.0, and multiple inheritance may well be possible in 1.x.

[atteson]

I don't quite follow why abstract multiple inheritance would be required to implement SortedVector. I have in mind something like the below (though sortv would be the sort method and many other methods would need to be defined for SorteVector, etc). I've been using something very similar for some time. The viewpoint is that a SortedVector is a data structure different from a vector with operations that take different amounts of time similar to DataStructures.jl.

Either way, I recognize that this is a fairly major change and this was more of a thought for 2.0. Given that this is the 2.0 sorting API, everything I've suggested here was meant for 2.0 though it would be great to have some of these features earlier.

struct SortedVector{T,V <: AbstractVector{T}} <: AbstractVector{T}
    v::V
end

search( v::AbstractVector{T}, x::T ) where {T} = findall( v .== x )

search( v::SortedVector{T}, x::T ) where {T} = searchsorted( v.v, x )

v = rand( 1:100_000, 1_000_000 );

@time l = search( v, 5 );
@time l = search( v, 5 );

sortv( v::AbstractVector{T} ) where {T} = SortedVector( sort(v) )

sv = sortv(v);

@time r = search( sv, 5 );
@time r = search( sv, 5 );

@assert( length(r) == length(l) )

[LilithHafner]

This is nice for many operations but comes at a cost. Currently, some methods already specialize on specific AbstractVectors such as SVector from StaticArrays. A SortedVector{T, SVector{5, T}} is no longer an SVector{5, T}. If you care about sortedness, this is good, but if you care about staticness, this is bad. Multiple inheritance might be able to solve this, but I haven't thought through all the details.

[LilithHafner]

The Perm ordering is inherently bad for caching. I believe that it should be removed entirely in 2.0. We could use something similar internally in special cases like sortperm of vectors of large tuples, but it is so infrequently the best approach that it should not be part of the public API.

[StefanKarpinski]

The Perm stuff is definitely something I just concocted without knowing any better, so please delete/deprecate at will.