Comment inverser un tableau de permutation dans numpy

Tri est inutile ici. il s'agit simplement d'un algorithme de temps linéaire en un seul passage avec un besoin constant de mémoire:

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)

Le code ci-dessus imprime

 s = [2 3 1 0]

selon les besoins.

le reste de La réponse est préoccupé par l'efficacité de la vectorisation de la ci-dessus for boucle. si vous voulez juste connaître la solution, sautez à la fin de ceci réponse.

(l'original de La réplique de Aug 27, 2014; les horaires sont valables pour NumPy 1.8. Une mise à jour avec NumPy 1.11 suit plus tard.)

un algorithme de temps linéaire en un seul passage devrait être plus rapide que np.argsort; fait intéressant, la vectorisation triviale (s[p] = xrange(p.size), voir tableaux index) ci-dessus for la boucle est en fait légèrement plus lente que np.argsort tant que p.size < 700 000 (bien sur ma machine, votre kilométrage varier):

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)

à Partir de mon 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

eventuellement la complexite asymptotique entre en jeu (O(n log n)argsort vs. O(n) pour l'algorithme à passage unique) et l'algorithme à passage unique sera toujours plus rapide après un n = p.size (seuil est d'environ 700k sur ma machine).

cependant, il y a une façon moins directe de vectoriser ce qui précède for boucle avec 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

ce Qui donne n = 700 000 (de la même taille que ci-dessus):

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

C'est une belle 5,6 x vitesse pour presque rien!

Pour être juste, np.argsort bat np.put petite n (le point de bascule est d'environ n = 1210 sur ma machine):

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

c'est très probablement parce que nous attribuons et remplissons un tableau supplémentaire (à np.arange() appel) avec le np_put approche.

bien que vous n'ayez pas demandé une solution Cython, juste par curiosité, j'ai aussi chronométré la solution Cython suivante avec taped 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

Horaires:

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

Donc,np.put la solution n'est toujours pas aussi rapide que possible (ran 12.8 ms pour cette taille d'entrée; argsort a pris 72.7 ms).

mise à jour le 3 février 2017 avec NumPy 1.11

Jamie, Andris et Paul ont souligné dans les commentaires ci-dessous que la question de la performance avec l'indexation de fantaisie a été résolue. Jamie dit que C'était déjà résolu à num Py 1.9. Je l'ai testé avec Python 3.5 et NumPy 1.11 sur la machine que j'utilisais en 2014.

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

Horaires:

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

une amélioration significative en effet!

Conclusion

dans l'Ensemble, j'irais avec le

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

approche pour le code de la clarté. À mon avis, il est moins obscur que argsort, et aussi plus rapide pour les grands d'entrée tailles. Si la vitesse devient un problème, j'opterais pour la solution Cython.

je voudrais offrir un peu plus de contexte à larsmans réponse correcte. raison pourquoi argsort est correct peut être trouvé lorsque vous utilisez la représentation d'un permutation par une matrice. L'avantage mathématique d'une permutation matriceP est que la matrice "fonctionne sur des vecteurs", c'est-à-dire qu'une matrice de permutation multipliée par un vecteur Permute le vecteur.

Votre permutation ressemble à:

import numpy as np
a   = np.array([3,2,0,1])
N   = a.size
rows = np.arange(N)
P   = np.zeros((N,N),dtype=int)
P[rows,a] = 1

[[0 0 0 1]
 [0 0 1 0]
 [1 0 0 0]
 [0 1 0 0]]

Donné un matrice de permutation, nous pouvons "annuler" la multiplication en la multipliant par son inverse P^-1. La beauté des matrices de permutation est qu'elles sont orthogonales, donc P*P^(-1)=I, ou en d'autres termes P(-1)=P^T, l'inverse est la transposition. Cela signifie que nous pouvons prendre les indices de la matrice de transposition pour trouver votre vecteur de permutation inversée:

inv_a = np.where(P.T)[1]
[2 3 1 0]

ce qui, si vous y pensez, est exactement la même chose que de trouver les indices qui trient les colonnes de P!