Creación de una matriz numpy de coordenadas 3D a partir de tres matrices 1D

Supongamos que tengo tres matrices 1D arbitrarias, por ejemplo:

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)) 

Estas tres matrices representan intervalos de muestreo en una cuadrícula 3D, y quiero construir una matriz 1D de vectores tridimensionales para todas las intersecciones, algo así como

 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]]) 

La orden en realidad no importa para esto. La forma obvia de generarlos:

 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 

Entonces la pregunta es … ¿hay una manera más rápida? He buscado, pero no he encontrado (posiblemente no he podido encontrar las palabras clave de Google adecuadas).

Actualmente estoy usando esto:

 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 

¿Pero siento que me estoy perdiendo alguna manera inteligente y elegante de Numpy?

Para usar la rejilla de malla numpy en el ejemplo anterior, funcionará lo siguiente:

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

Numpy meshgrid para cuadrículas de más de dos dimensiones requieren numpy 1.7. Para evitar esto y extraer los datos relevantes del código fuente .

 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 ND 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) 

Comprobando el resultado:

 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.]] 

Para el ejemplo anterior:

 %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 

Para cuando cada dimensión es 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 

Dependiendo de lo que desee hacer con los datos, puede devolver una vista:

 %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 

Para su ejemplo específico, mgrid es bastante útil:

 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]]) 

Más generalmente, si tiene matrices específicas con las que desea hacer esto, usaría meshgrid lugar de mgrid . Sin embargo, necesitará numpy 1.7 o posterior para que funcione en más de dos dimensiones.

Puedes usar itertools.product :

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

Otra solución que usa iteradores y np.take demostró ser aproximadamente 3 veces más rápida para su 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 

Resultados de tiempo:

 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 

Por lo tanto, solo mi segunda solución y la solución de Joe Kington le darían cierta ganancia de rendimiento …

Para aquellos que tuvieron que quedarse con numpy <1.7.x, aquí hay una solución simple de dos líneas:

 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]])) 

Muy fácil de expandir a una dimensión aún mayor también. ¡Aclamaciones!