In-Place-Shuffling von mehrdimensionalen Arrays

Question

Ich versuche, ein NaN-sicheres Shuffling-Verfahren in Cython zu implementieren, das entlang mehrerer Achsen einer mehrdimensionalen Matrix beliebiger Dimension mische.

Im einfachen Fall einer 1D-Matrix kann man einfach alle Indizes mit nicht-NaN Shuffle über Werte, die Fisher-Yates-Algorithmus:

def shuffle1D(np.ndarray[double, ndim=1] x): 
    cdef np.ndarray[long, ndim=1] idx = np.where(~np.isnan(x))[0] 
    cdef unsigned int i,j,n,m 

    randint = np.random.randint 
    for i in xrange(len(idx)-1, 0, -1): 
     j = randint(i+1) 
     n,m = idx[i], idx[j] 
     x[n], x[m] = x[m], x[n] 

Ich möchte diesen Algorithmus erweitern multidimensional zu handhaben große Arrays ohne Umformung (wodurch eine Kopie für kompliziertere Fälle ausgelöst wird, die hier nicht berücksichtigt werden). Zu diesem Zweck müsste ich die feste Eingabedimension loswerden, die mit numpy arrays oder Memoryviews in Cython nicht möglich ist. Gibt es eine Problemumgehung?

Vielen Dank im Voraus!

Answer 1

Dank der Kommentare von @Veedrac diese Antwort mehr von Cython Fähigkeiten nutzt.

• Ein Zeiger-Array speichert die Speicheradresse der Werte entlang axis
• Ihr Algorithmus mit einer Modifikation verwendet wird that checks for nan values, sie zu verhindern, die sortiert werden
• Es wird keine Kopie für C geordnete Anordnungen erstellen. Bei Fortran geordneten Arrays gibt der Befehl ravel() eine Kopie zurück. Dies könnte durch die Schaffung eine andere Anordnung von Doppelzeigern verbessert werden, um die Werte von x, wahrscheinlich mit einiger Cache Strafe ...
ist

mindestens eine Größenordnung Dieser Code zu tragen schneller als die anderen basierend auf Scheiben.

from libc.stdlib cimport malloc, free 

cimport numpy as np 
import numpy as np 
from numpy.random import randint 

cdef extern from "numpy/npy_math.h": 
    bint npy_isnan(double x) 

def shuffleND(x, int axis=-1): 
    cdef np.ndarray[double, ndim=1] v # view of x 
    cdef np.ndarray[int, ndim=1] strides 
    cdef int i, j 
    cdef int num_axis, pos, stride 
    cdef double tmp 
    cdef double **v_axis 

    if axis==-1: 
     axis = x.ndim-1 

    shape = list(x.shape) 
    num_axis = shape.pop(axis) 

    v_axis = <double **>malloc(num_axis*sizeof(double *)) 
    for i in range(num_axis): 
     v_axis[i] = <double *>malloc(1*sizeof(double)) 

    try: 
     tmp_strides = [s//x.itemsize for s in x.strides] 
     stride = tmp_strides.pop(axis) 
     strides = np.array(tmp_strides, dtype=np.int32) 
     v = x.ravel() 
     for indices in np.ndindex(*shape): 
      pos = (strides*indices).sum() 
      for i in range(num_axis): 
       v_axis[i] = &v[pos + i*stride] 
      for i in range(num_axis-1, 0, -1): 
       j = randint(i+1) 
       if npy_isnan(v_axis[i][0]) or npy_isnan(v_axis[j][0]): 
        continue 
       tmp = v_axis[i][0] 
       v_axis[i][0] = v_axis[j][0] 
       v_axis[j][0] = tmp 
    finally: 
     free(v_axis) 

    return x 

Answer 2

Der folgende Algorithmus basiert auf Schichten, in denen keine Kopie erstellt wird, und es sollte für alle np.ndarray funktionieren. Die wichtigsten Schritte sind:

• np.ndindex() verwendet wird throught die verschiedenen mehrdimensionalen Indizes zu laufen, die eine Ausnahme von der Achse gehören, Sie
• die Shuffle bereits entwickelt von Ihnen für die 1-D Fall mischen wollen angewendet .

Code:

def shuffleND(np.ndarray x, axis=-1): 
    cdef np.ndarray[long long, ndim=1] idx 
    cdef unsigned int i, j, n, m 
    if axis==-1: 
     axis = x.ndim-1 
    all_shape = list(np.shape(x)) 
    shape = all_shape[:] 
    shape.pop(axis) 
    for slices in np.ndindex(*shape): 
     slices = list(slices) 
     axis_slice = slices[:] 
     axis_slice.insert(axis, slice(None)) 
     idx = np.where(~np.isnan(x[tuple(axis_slice)]))[0] 
     for i in range(idx.shape[0]-1, 0, -1): 
      j = randint(i+1) 
      n, m = idx[i], idx[j] 
      slice1 = slices[:] 
      slice1.insert(axis, n) 
      slice2 = slices[:] 
      slice2.insert(axis, m) 
      slice1 = tuple(slice1) 
      slice2 = tuple(slice2) 
      x[slice1], x[slice2] = x[slice2], x[slice1] 
    return x