Mescolamento in loco di array multidimensionali

Question

Sto cercando di implementare una procedura di shuffling NaN-safe in Cython che può mescolare lungo diversi assi di una matrice multidimensionale di dimensione arbitraria.

Nel caso semplice di una matrice 1D, si può semplicemente mescolare sopra tutti gli indici con non-NaN valori utilizzando l'algoritmo di Fisher-Yates:

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] 

Vorrei estendere questo algoritmo per gestire grandi multidimensionale array senza risagoma (che innesca una copia per casi più complicati non considerati qui). A tal fine, avrei bisogno di eliminare la dimensione di input fissa, che non sembra possibile con gli array numpy né con le visualizzazioni di memoria in Cython. C'è una soluzione?

Molte grazie in anticipo!

Answer 1

Grazie ai commenti di @Veedrac questa risposta utilizza più di capacità Cython.

• Una matrice puntatore memorizza l'indirizzo di memoria dei valori lungo axis
• l'algoritmo viene utilizzato con una modifica that checks for nan values, impedendo loro di essere ordinati
• non creerà una copia per C array ordinati. Nel caso degli array ordinati Fortran, il comando ravel() restituirà una copia. Questo potrebbe essere migliorata con la creazione di un altro array di doppi puntatori per portare i valori di x, probabilmente con un po 'di pena di cache ...

Questo codice è almeno un ordine di grandezza più veloce rispetto agli altri sulla base di fette.

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

Il seguente algoritmo è basato su sezioni, in cui non viene eseguita alcuna copia e dovrebbe funzionare per qualsiasi np.ndarray. Le fasi principali sono:

• np.ndindex() utilizzato per eseguire throught i diversi indici multidimensionali, escludendo quello appartenente all'asse si vuole mischiare
• shuffle già sviluppato da voi per il caso 1-D viene applicata .

Codice:

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