missing sortperm(A, dim)?

[stevengj]

We have sort(A, dim) but not sortperm(A, dim). Is there any reason for this omission?

(See also #9258, which is slightly different because sortrows(A) is not the same thing as sort(A, 1).)

[stevengj]

This was just brought up on julia-users

[stevengj]

Note that it should probably return scalar indices that can be passed to A[i], regardless of the dimensions of A. For example, here is an implementation for 2d arrays:

function Base.sortperm(A::AbstractMatrix, dim::Integer)
    P = mapslices(sortperm, A, dim)
    for j = 1:size(P,2)
        offset = (j-1) * size(P,1)
        for i = 1:size(P,1)
            P[i,j] += offset
        end
    end
    return P
end

e.g. try A = rand(1:100, 5,4); i = sortperm(A, 1); A[i].

[mes9]

That example works with dim=1, but when it sorts the values across rows, it will not add the right offset amount. It's not clear to me how to keep track of the index in an efficient way... An inelegant solution is to transpose the matrix and treat it as sorting over dim=1

[stevengj]

Whoops, right. To handle both dim == 1 and dim == 2 you need something like:

function Base.sortperm(A::AbstractMatrix, dim::Integer)
    P = mapslices(sortperm, A, dim)
    if dim == 1
        for j = 1:size(P,2)
            offset = (j-1) * size(P,1)
            for i = 1:size(P,1)
                P[i,j] += offset
            end
        end
    else # if dim == 2
        for j = 1:size(P,2)
            for i = 1:size(P,1)
                P[i,j] = (P[i,j] - 1) * size(P,1) + i
            end
        end
    end
    return P
end

[timholy]

It seems better to use CartesianIndex and sub2ind for this kind of chicanery---then you can write one algorithm that works for any number of dimensions. The "exponential smoothing" example in http://julialang.org/blog/2016/02/iteration#filtering-along-a-specified-dimension-exploiting-multiple-indexes could almost be lifted directly:

for Ipost in Rpost
    for i = 1:size(P, dim)
        for Ipre in Rpre
            P[Ipre, i, Ipost] = sub2ind(P, Ipre.I..., P[Ipre, i, Ipost], Ipost.I...)
        end
    end
end

all inside an @inline (to avoid the splat penalty).

(Parenthetically, it's annoying that the blog doesn't seem to support links for specific sections.)

EDIT jwn: updated link to point directly to the specific section. it's a shame our blog layout doesn't expose these links more obviously.

[stevengj]

@timholy, I was thinking that what I wrote is probably a lot faster for the common case of 2d matrices, but it would be interesting to see a benchmark. Maybe the index-correction speed is irrelevant because it spends all of its time in the mapslices.

[timholy]

Aside from #9080 (which last I timed it was about a 20% penalty), it shouldn't be. A lovely thing about Julia is that the generalization penalty is often nonexistent. (At least until you hit MAX_TUPLETYPE_LEN or whatever...)

[mirestrepo]

Was there a consensus on a good solution? Is there a plan to incorporate this feature?

[aminya]

I edited @stevengj's solution for AbstractMatrices to be compatible with the current syntax:

function Base.sortperm(A::AbstractMatrix; dims::Integer, rev::Bool = false)
    P = mapslices(x -> sortperm(x; rev = rev), A, dims = dims)
    if dims == 1
        for j = 1:size(P, 2)
            offset = (j - 1) * size(P, 1)
            for i = 1:size(P, 1)
                P[i, j] += offset
            end
        end
    else # if dims == 2
        for j = 1:size(P, 2)
            for i = 1:size(P, 1)
                P[i, j] = (P[i, j] - 1) * size(P, 1) + i
            end
        end
    end
    return P
end

[mbauman]

Do we really want to return an array of linear indices here? That doesn't really sound right to me anymore.

Then do we really want to compute an array of CartesianIndexes? I suppose we could, but is that really any more useful than the raw output from mapslices(sortperm, ...)? I guess I don't fully understand the use-cases here.

[aminya]

Do we really want to return an array of linear indices here? That doesn't really sound right to me anymore.

Then do we really want to compute an array of CartesianIndexes? I suppose we could, but is that really any more useful than the raw output from mapslices(sortperm, ...)? I guess I don't fully understand the use-cases here.

I think the check is this:

A[sortperm(A, dims = dim)] == sort(A, dims = dim)