Julia: partialsortperm! silently gives wrong result

[AzamatB]

I believe I have found a silent bug with partialsortperm! (almost went nuts debugging the huge code base where this was happening). Here is the MRE (compare the outputs of the partialsortperm! with partialsortperm and the outputs of v[ix] in both cases):

julia> v = [3, 1, 5, 2, 4,
            3, 1, 5, 2, 4]
10-element Array{Int64,1}:
 3
 1
 5
 2
 4
 3
 1
 5
 2
 4

julia> k = 1:5
1:5

julia> ix = collect(k)
5-element Array{Int64,1}:
 1
 2
 3
 4
 5

julia> partialsortperm!(ix, v, k)
5-element view(::Array{Int64,1}, 1:5) with eltype Int64:
 2
 4
 1
 5
 3

julia> v[ix]
5-element Array{Int64,1}:
 1
 2
 3
 4
 5

julia> ix = partialsortperm(v, k)
5-element view(::Array{Int64,1}, 1:5) with eltype Int64:
 2
 7
 4
 9
 1

julia> v[ix]
5-element Array{Int64,1}:
 1
 1
 2
 2
 3

julia> versioninfo()
Julia Version 1.2.0
Commit c6da87ff4b (2019-08-20 00:03 UTC)
Platform Info:
  OS: macOS (x86_64-apple-darwin18.6.0)
  CPU: Intel(R) Core(TM) i7-4770HQ CPU @ 2.20GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-6.0.1 (ORCJIT, haswell)
Environment:
  JULIA_EDITOR = atom  -a
  JULIA_NUM_THREADS = 4

[KristofferC]

ix should be collect(1:10), as in

julia> v[partialsortperm!(collect(1:10), v, k)]
5-element Array{Int64,1}:
 1
 1
 2
 2
 3

, no?

[AzamatB]

Well, from the documentation of partialsortperm!:

Like partialsortperm, but accepts a preallocated index vector ix

And given that partialsortperm returns ix of the same length as k, shouldn't the preallocated version behave the same way, i.e. also take ix of the same size as k. Otherwise, it's kinda contradicts it's documentation and also the general use of bang in other mutating-nonmutating pairs of identical functions

[KristofferC]

It mutates the passed in vector but gives back a view of length k. Agree that the docs could be improved though.

[AzamatB]

I don't think that this is just a documentation issue.
To me the only advantage of partialsortperm! over partialsortperm is that you can use it like

k = 1:5
ix = collect(k)

partialsortperm!(ix, v1, k)
# do something with v1[ix]
...
# do something with u1[ix]
...
partialsortperm!(ix, v2, k)
# do something with v2[ix]
...
# do something with u2[ix]
...

If ix is not of the same size as k it would be this would not be possible, and a workaround would have to allocate.

EDIT: see more poignant example in my next post below.

[KristofferC]

I thought the advantage is that you don't have to allocate the ix vector with the ! version (you are passing in the workspace).

[AzamatB]

@KristofferC yes, exactly. My real use case is something like

k = 1:5
ix = collect(k)
# vofvs is a vector of vectors, whose elements are not necessarily
# of the same length, but each one is longer than k
for v in vofvs
   partialsortperm!(ix, v, k)
   # do something with ix
   ...
end

which in its current form produces the wrong result, and to fix it I will have to allocate at every iteration of the loop defeating the purpose of the ! version.

[KristofferC]

You just have to use the return value instead (a view of ix) instead of ix itself.

[AzamatB]

No, that will not solve it, because with the current interface ix must match v in length.
Besides the fact that it diverges from the meaning of ! in every other mutating-nonmutating pairs of identical functions, this means that the loop code above must be rewritten as

for v in vofvs
   ix = collect(1:length(v))
   ix_view = partialsortperm!(ix, v, k)
   # do something with ix_view
   ...
end

i.e. it will require reallocatiing ix at every loop iteration in order to produce correct result. So this is not a documentation issue, but the one of interface.

[haampie]

You can still have your preallocation if you write it like this:

ix = Vector{Int}(undef, mapreduce(length, max, vs))

and then pass a view of the correct length to partialsortperm!

for v in vs
  ix_view = partialsortperm!(view(ix, eachindex(v)), v, k)
end

[AzamatB]

@haampie Sure, but I want the language to help me instead of having to bend over backwards to get the desired behavior. This is ignoring the fact that this workaround still has its own problems (like when k is short and the longest v is very long).

Granted I don't know the implementation details of partialsortperm and its ! version, but from the perspective of the user of those functions, the documented behavior of partialsortperm! is clearly better and more useful than it's actual current behavior, which is why I don't all the defensiveness about it.

[AzamatB]

@KristofferC could you remove the doc label on this issue?
I would say that this is a bug as it contradicts the documented behavior and we should fix the behavior, not the documentation for reasons outlined above.

[haampie]

The docs should definitely be improved to avoid the confusion, so no need to remove the label.

Note that sortperm!([1], [3,2,1]) throws an an ArgumentError "index vector must have the same indices as the source vector", one can only assume it's an oversight this lacks in partialsortperm!.

For your understanding, the function is named sortperm!(p, v) because it sorts a _permutation_ p of the indices of a vector v, which is equivalent* to calling sort!(p, by = i -> v[i]). Analogously partialsortperm!(p, v, k) partially sorts a _permutation_ p and is equivalent* to partialsort!(p, k, by = i -> v[i]). For convenience, you do not have to provide a permutation of the indices, you just have to provide a workspace in which the function can store a dummy permutation. As a consequence, it should have the same size as the vector.

* not always strictly equivalent in that sortperm additionally requires stable sort

[KristofferC]

i.e. it will require reallocatiing ix at every loop iteration in order to produce correct result

It doesn't, just use resize!?

[AzamatB]

just use resize!

Again, I'm arguing that this is an implementation detail, and should be hidden away from the user. That is my feedback as a user, who had a really bad experience with that function. You can, of course, take it as you wish.

@haampie in its current form the pair partialsortperm – partialsortperm! does not have the symmetry between them that every other allocating – nonallocating pair of identical functions in Base enjoys (e.g. append – append!, fill – fill!, etc.). This, together with the fact that it is documented to have that symmetry and does not throw any error when you assume it, is an indication that the current interface of partialsortperm! is broken.

[KristofferC]

It is fully possible to use the current interface to repeatedly call the function without reallocating the work buffer, which, to me, shows that there needs to be more docs of how you use it correctly. In addition there might be good to add length check on the input vector.

[haampie]

the pair partialsortperm – partialsortperm! does not have the symmetry between them that every other allocating – nonallocating pair of identical functions in Base enjoys

It's analogous to copy(x) - copy!(y, x). The docs of partialsortperm! should simply be updated to read "Like partialsortperm, but accepts a preallocated permutation vector of the indices of v", and an assertion on the length should be added. A PR is welcome.

[AzamatB]

@KristofferC I honestly don't understand this. It is fully possible to use C for every Julia use case since both languages are Turing complete, but I hardly can see how that can be an argument in favor of C. partialsortperm – partialsortperm! is an exception from the symmetry between every other mutating – nonmutating function pair in Base and there is no any reason for this exception.

@haampie How is it analogous? You don't have to do any resize!ing and can just turn y = copy(x) into copy!(y, x), so in that regard, it is analogous to every other mutating – nonmutating function pair in Base and _is at odds_ with partialsortperm!.

[haampie]

So maybe your confusion with the method has to do with the return value being a view of the mutated permutation, rather than the full mutated permutation?

Notice that partialsort! has the same behavior:

julia> partialsort!(rand(10), 5:8)
4-element view(::Array{Float64,1}, 5:8) with eltype Float64:
 0.43184474913373627
 0.60587737927092   
 0.6914000729286978 
 0.7451764940939822

[KristofferC]

and there is no any reason for this exception.

The function requires length(v) workspace but returns a result of length k. (As you have already seen), if you only pass it a workspace of length k it doesn't work.

julia> partialsortperm!(ix, v, k)'
1×5 LinearAlgebra....
 2  7  4  9  1

julia> ix'
1×10 LinearAlgebra....
 2  7  4  9  1  6  5  10  3  8

The stuff after k is also modified.

There is no getting around the fact that if you don't want this function to not allocate, you need to pass it a work buffer of length v because that is what it needs. And if you want the work buffer to be directly used as the permutation result, the function would need to resize! it but that is kinda weird and not really done.

[dkarrasch]

I spent a lot of time to understand partialsortperm! some time ago, but I simply couldn't understand what it's doing, what it is returning, and how to retrieve the correct result of partialsortperm, yet in a preallocated vector...

https://discourse.julialang.org/t/what-does-partialsortperm-do/25656

[haampie]

I'm thinking that in principle it should not be necessary to make partialsortperm! throw when the first argument is not a permutation of the indices. Maybe the names of the functions sortperm! and partialsortperm! are just off.

Both sortperm!(ix, x) and partialsortperm!(ix, x, k) operate on a vector ix of valid indices of x and do not exploit the fact that ix is a permutation anywhere. So if you drop the length check in sortperm!, the invariants that hold are

sort(x[ix]) == x[sortperm!(ix, x, initialized=true)] 
sort(x[ix])[k] == x[partialsortperm!(ix, x, k, initialized=true)]

With the current implementation of base, the second invariant always holds since there is no size check on ix. The first holds only when you pass in a vector of indices of the right size, but does not fail when the input is not a permutation:

julia> ix = [1,2,5,2,1,7,1,1,1,2];
julia> x = rand(10)
julia> isperm(ix)
false
julia> sort(x[ix]) == x[sortperm!(ix, x, initialized=true)]
true

So maybe sortindices! and partialsortindices! are better names?