Wie zum invertieren einer permutation array in numpy

Sortierung ist ein übermaß hier. Dies ist nur ein single-pass-linearzeit-Algorithmus mit konstanter Speicherbedarf:

from __future__ import print_function
import numpy as np

p = np.array([3, 2, 0, 1])
s = np.empty(p.size, dtype=np.int32)
for i in np.arange(p.size):
    s[p[i]] = i

print('s =', s)

Der obige code druckt

 s = [2 3 1 0]

als erforderlich.

Den rest der Antwort befasst sich mit der effizienten Vektorisierung der oben for Schleife. , Wenn Sie nur wollen, die Lösung zu wissen, springen Sie an das Ende dieser Antwort.

(Die original-Antwort von Aug 27, 2014; die timings sind gültig für NumPy 1.8. Ein update mit NumPy 1.11 folgt später.)

Single-pass-linearzeit-Algorithmus wird erwartet, schneller zu sein als np.argsort; interessanterweise die trivial-Vektorisierung (s[p] = xrange(p.size), siehe index-arrays) der oben genannten for Schleife ist tatsächlich etwas langsamer, als np.argsort solange p.size < 700 000 (gut, auf meiner Maschine, Ihre Laufleistung wird variieren):

import numpy as np

def np_argsort(p):
    return np.argsort(p)

def np_fancy(p):
    s = np.zeros(p.size, p.dtype) # np.zeros is better than np.empty here, at least on Linux
    s[p] = xrange(p.size) 
    return s

def create_input(n):
    np.random.seed(31)
    indices = np.arange(n, dtype = np.int32)
    return np.random.permutation(indices)

Aus meiner IPython notebook:

p = create_input(700000)
%timeit np_argsort(p)
10 loops, best of 3: 72.7 ms per loop
%timeit np_fancy(p)
10 loops, best of 3: 70.2 ms per loop

Schließlich die asymptotische Komplexität kicks (O(n log n) für argsort vs. O(n) für die single-pass-Algorithmus) und die single-pass-Algorithmus, werden konsequent schneller nach einem ausreichend großen n = p.size (Schwelle etwa 700k auf meinem Rechner).

Allerdings gibt es eine weniger einfache Weise zu Vektorisieren die oben for Schleife mit np.put:

def np_put(p):
    n = p.size
    s = np.zeros(n, dtype = np.int32)
    i = np.arange(n, dtype = np.int32)
    np.put(s, p, i) # s[p[i]] = i 
    return s

Gibt für n = 700 000 (die gleiche Größe wie oben):

p = create_input(700000)
%timeit np_put(p)
100 loops, best of 3: 12.8 ms per loop

Dies ist ein schönes 5,6 x speed-up für den nächsten zu nichts!

Um fair zu sein, np.argsort schlägt immer noch die np.put Ansatz für kleinere n (der Wendepunkt ist um n = 1210 auf meiner Maschine):

p = create_input(1210)
%timeit np_argsort(p)
10000 loops, best of 3: 25.1 µs per loop
%timeit np_fancy(p)
10000 loops, best of 3: 118 µs per loop
%timeit np_put(p)
10000 loops, best of 3: 25 µs per loop

Dies ist wahrscheinlich, weil wir reservieren und füllen Sie in ein extra array (an der np.arange() call) mit der np_put Ansatz.

Obwohl Sie nicht Fragen, für eine Cython-Lösung, nur aus Neugier, ich auch zeitlich die folgenden Cython-Lösung mit eingegeben memoryviews:

import numpy as np
cimport numpy as np

def in_cython(np.ndarray[np.int32_t] p):    
    cdef int i
    cdef int[:] pmv
    cdef int[:] smv 
    pmv = p
    s = np.empty(p.size, dtype=np.int32)
    smv = s
    for i in xrange(p.size):
        smv[pmv[i]] = i
    return s

Timings:

p = create_input(700000)
%timeit in_cython(p)
100 loops, best of 3: 2.59 ms per loop

So, die np.put Lösung ist noch nicht so schnell wie möglich (lief 12.8 ms für diese Eingangs Größe; argsort nahm 72.7 ms).

Update am Feb 3, 2017 mit NumPy 1.11

Jamie, Andris und Paul darauf hingewiesen, in den Kommentaren unten, dass das performance-Problem mit fancy indexing gelöst wurde. Jamie sagt, es sei bereits aufgelöst NumPy 1.9. Getestet habe ich es mit Python 3.5 und NumPy 1.11 auf die Maschine, die ich mit zurück in 2014.

def invert_permutation(p):
    s = np.empty(p.size, p.dtype)
    s[p] = np.arange(p.size)
    return s

Timings:

p = create_input(880)
%timeit np_argsort(p)
100000 loops, best of 3: 11.6 µs per loop
%timeit invert_permutation(p)
100000 loops, best of 3: 11.5 µs per loop

Eine signifikante Verbesserung, in der Tat!

Abschluss

Alles in allem, würde ich mit dem

def invert_permutation(p):
    '''The argument p is assumed to be some permutation of 0, 1, ..., len(p)-1. 
    Returns an array s, where s[i] gives the index of i in p.
    '''
    s = np.empty(p.size, p.dtype)
    s[p] = np.arange(p.size)
    return s

Ansatz für die code-Klarheit. Meiner Meinung nach ist es weniger dunkel als argsort, und auch schneller für große input-Größen. Wenn Geschwindigkeit ein Problem wird, würde ich mit dem Cython-Lösung.