Creazione di un array NumPy di ​​3D coordinate da tre matrici 1D

Question

Supponiamo di avere tre array 1D arbitrari, ad esempio:

x_p = np.array((1.0, 2.0, 3.0, 4.0, 5.0)) 
y_p = np.array((2.0, 3.0, 4.0)) 
z_p = np.array((8.0, 9.0)) 

Questi tre matrici rappresentano intervalli di campionamento in una griglia 3D, e voglio costruire un 1D serie di vettori tridimensionali per tutte le intersezioni, qualcosa come

points = np.array([[1.0, 2.0, 8.0], 
        [1.0, 2.0, 9.0], 
        [1.0, 3.0, 8.0], 
        ... 
        [5.0, 4.0, 9.0]]) 

Ordinare in realtà non importa per questo. Il modo più ovvio per generarli:

npoints = len(x_p) * len(y_p) * len(z_p) 
points = np.zeros((npoints, 3)) 
i = 0 
for x in x_p: 
    for y in y_p: 
     for z in z_p: 
      points[i, :] = (x, y, z) 
      i += 1 

Quindi la domanda è ... c'è un modo più veloce? Ho cercato ma non trovato (forse non sono riuscito a trovare le parole chiave di Google corrette).

Attualmente sto usando questo:

npoints = len(x_p) * len(y_p) * len(z_p) 
points = np.zeros((npoints, 3)) 
i = 0 
nz = len(z_p) 
for x in x_p: 
    for y in y_p: 
     points[i:i+nz, 0] = x 
     points[i:i+nz, 1] = y 
     points[i:i+nz, 2] = z_p 
     i += nz 

ma mi sento come mi manca un po 'di fantasia modo Numpy intelligente?

Answer 1

Per utilizzare griglia maglia numpy sul suddetto esempio seguente funziona:

np.vstack(np.meshgrid(x_p,y_p,z_p)).reshape(3,-1).T 

Numpy meshgrid per reti di più di due dimensioni richiedono NumPy 1.7. Per aggirare questo e tirare i dati rilevanti dal source code.

def ndmesh(*xi,**kwargs): 
    if len(xi) < 2: 
     msg = 'meshgrid() takes 2 or more arguments (%d given)' % int(len(xi) > 0) 
     raise ValueError(msg) 

    args = np.atleast_1d(*xi) 
    ndim = len(args) 
    copy_ = kwargs.get('copy', True) 

    s0 = (1,) * ndim 
    output = [x.reshape(s0[:i] + (-1,) + s0[i + 1::]) for i, x in enumerate(args)] 

    shape = [x.size for x in output] 

    # Return the full N-D matrix (not only the 1-D vector) 
    if copy_: 
     mult_fact = np.ones(shape, dtype=int) 
     return [x * mult_fact for x in output] 
    else: 
     return np.broadcast_arrays(*output) 

Controllo del risultato:

print np.vstack((ndmesh(x_p,y_p,z_p))).reshape(3,-1).T 

[[ 1. 2. 8.] 
[ 1. 2. 9.] 
[ 1. 3. 8.] 
.... 
[ 5. 3. 9.] 
[ 5. 4. 8.] 
[ 5. 4. 9.]] 

Per l'esempio precedente:

%timeit sol2() 
10000 loops, best of 3: 56.1 us per loop 

%timeit np.vstack((ndmesh(x_p,y_p,z_p))).reshape(3,-1).T 
10000 loops, best of 3: 55.1 us per loop 

Per cui ogni dimensione è 100:

%timeit sol2() 
1 loops, best of 3: 655 ms per loop 
In [10]: 

%timeit points = np.vstack((ndmesh(x_p,y_p,z_p))).reshape(3,-1).T 
10 loops, best of 3: 21.8 ms per loop 

A seconda di ciò che si vuole fare con i dati, è possibile restituire av IEW:

%timeit np.vstack((ndmesh(x_p,y_p,z_p,copy=False))).reshape(3,-1).T 
100 loops, best of 3: 8.16 ms per loop 

Answer 2

È possibile utilizzare itertools.product:

def sol1(): 
    points = np.array(list(product(x_p, y_p, z_p))) 

Un'altra soluzione utilizzando iteratori e np.take dimostrato di essere di circa 3 volte più veloce per il vostro caso:

from itertools import izip, product 

def sol2(): 
    points = np.empty((len(x_p)*len(y_p)*len(z_p),3)) 

    xi,yi,zi = izip(*product(xrange(len(x_p)), 
           xrange(len(y_p)), 
           xrange(len(z_p)) )) 

    points[:,0] = np.take(x_p,xi) 
    points[:,1] = np.take(y_p,yi) 
    points[:,2] = np.take(z_p,zi) 
    return points 

risultati Timing:

In [3]: timeit sol1() 
10000 loops, best of 3: 126 µs per loop 

In [4]: timeit sol2() 
10000 loops, best of 3: 42.9 µs per loop 

In [6]: timeit ops() 
10000 loops, best of 3: 59 µs per loop 

In [11]: timeit joekingtons() # with your permission Joe... 
10000 loops, best of 3: 56.2 µs per loop 

Così , solo la mia seconda soluzione e la soluzione di Joe Kington ti avrebbe dato qualche guadagno di prestazioni ...

Answer 3

Per esempio specifico, mgrid è molto utile .:

In [1]: import numpy as np 
In [2]: points = np.mgrid[1:6, 2:5, 8:10] 
In [3]: points.reshape(3, -1).T 
Out[3]: 
array([[1, 2, 8], 
     [1, 2, 9], 
     [1, 3, 8], 
     [1, 3, 9], 
     [1, 4, 8], 
     [1, 4, 9], 
     [2, 2, 8], 
     [2, 2, 9], 
     [2, 3, 8], 
     [2, 3, 9], 
     [2, 4, 8], 
     [2, 4, 9], 
     [3, 2, 8], 
     [3, 2, 9], 
     [3, 3, 8], 
     [3, 3, 9], 
     [3, 4, 8], 
     [3, 4, 9], 
     [4, 2, 8], 
     [4, 2, 9], 
     [4, 3, 8], 
     [4, 3, 9], 
     [4, 4, 8], 
     [4, 4, 9], 
     [5, 2, 8], 
     [5, 2, 9], 
     [5, 3, 8], 
     [5, 3, 9], 
     [5, 4, 8], 
     [5, 4, 9]]) 

Più in generale, se si dispone di array specifici che si desidera di fare questo con, devi usare invece meshgrid di mgrid. Tuttavia, avrai bisogno di numpy 1.7 o versioni successive affinché funzioni in più di due dimensioni.

Answer 4

Per chi doveva restare con NumPy < 1.7.x, ecco una semplice soluzione a due-liner:

g_p=np.zeros((x_p.size, y_p.size, z_p.size)) 
array_you_want=array(zip(*[item.flatten() for item in \ 
           [g_p+x_p[...,np.newaxis,np.newaxis],\ 
            g_p+y_p[...,np.newaxis],\ 
            g_p+z_p]])) 

Molto facile da espandere da addirittura superiore dimenision pure. Cheers!