Usando numpy para construir una matriz de todas las combinaciones de dos matrices

143

Estoy tratando de ejecutar el espacio de parámetros de una función de 6 parámetros para estudiar su comportamiento numérico antes de intentar hacer algo complejo con él, así que estoy buscando una manera eficiente de hacerlo.

Mi función toma valores flotantes dada una matriz numpy de 6 dim como entrada. Lo que intenté hacer inicialmente fue esto:

Primero creé una función que toma 2 matrices y genera una matriz con todas las combinaciones de valores de las dos matrices

from numpy import *
def comb(a,b):
    c = []
    for i in a:
        for j in b:
            c.append(r_[i,j])
    return c

Luego solía reduce()aplicar eso a m copias de la misma matriz:

def combs(a,m):
    return reduce(comb,[a]*m)

Y luego evalúo mi función así:

values = combs(np.arange(0,1,0.1),6)
for val in values:
    print F(val)

Esto funciona pero es muuuuy lento. Sé que el espacio de parámetros es enorme, pero esto no debería ser tan lento. Solo he muestreado 10 6 (un millón) puntos en este ejemplo y me llevó más de 15 segundos crear la matriz values.

¿Conoces alguna forma más eficiente de hacer esto con numpy?

Puedo modificar la forma en que la función Ftoma sus argumentos si es necesario.

Rafael S. Calsaverini
fuente
Para el producto cartesiano más rápido que he encontrado, vea esta respuesta . (Dado que la pregunta se formula de manera muy diferente de éste, considero que las preguntas no son duplicados, pero la mejor solución a las dos preguntas es la misma.)
senderle

Respuestas:

127

En la versión más reciente de numpy(> 1.8.x), numpy.meshgrid()proporciona una implementación mucho más rápida:

La solución de @ pv

In [113]:

%timeit cartesian(([1, 2, 3], [4, 5], [6, 7]))
10000 loops, best of 3: 135 µs per loop
In [114]:

cartesian(([1, 2, 3], [4, 5], [6, 7]))

Out[114]:
array([[1, 4, 6],
       [1, 4, 7],
       [1, 5, 6],
       [1, 5, 7],
       [2, 4, 6],
       [2, 4, 7],
       [2, 5, 6],
       [2, 5, 7],
       [3, 4, 6],
       [3, 4, 7],
       [3, 5, 6],
       [3, 5, 7]])

numpy.meshgrid()solía ser solo 2D, ahora es capaz de ND. En este caso, 3D:

In [115]:

%timeit np.array(np.meshgrid([1, 2, 3], [4, 5], [6, 7])).T.reshape(-1,3)
10000 loops, best of 3: 74.1 µs per loop
In [116]:

np.array(np.meshgrid([1, 2, 3], [4, 5], [6, 7])).T.reshape(-1,3)

Out[116]:
array([[1, 4, 6],
       [1, 5, 6],
       [2, 4, 6],
       [2, 5, 6],
       [3, 4, 6],
       [3, 5, 6],
       [1, 4, 7],
       [1, 5, 7],
       [2, 4, 7],
       [2, 5, 7],
       [3, 4, 7],
       [3, 5, 7]])

Tenga en cuenta que el orden de la resultante final es ligeramente diferente.

CT Zhu
fuente
14
np.stack(np.meshgrid([1, 2, 3], [4, 5], [6, 7]), -1).reshape(-1, 3)dará el orden correcto
Eric
@CT Zhu ¿Hay una manera fácil de transformar esto para que la matriz que contiene las diferentes matrices como columnas se use como entrada?
Dole
2
Cabe señalar que meshgrid solo funciona para conjuntos de rango más pequeños, tengo uno grande y recibo un error: ValueError: la dimensión máxima admitida para un ndarray es 32, encontrado 69
mikkom
157

Aquí hay una implementación puramente numpy. Es aproximadamente 5 veces más rápido que usar itertools.


import numpy as np

def cartesian(arrays, out=None):
    """
    Generate a cartesian product of input arrays.

    Parameters
    ----------
    arrays : list of array-like
        1-D arrays to form the cartesian product of.
    out : ndarray
        Array to place the cartesian product in.

    Returns
    -------
    out : ndarray
        2-D array of shape (M, len(arrays)) containing cartesian products
        formed of input arrays.

    Examples
    --------
    >>> cartesian(([1, 2, 3], [4, 5], [6, 7]))
    array([[1, 4, 6],
           [1, 4, 7],
           [1, 5, 6],
           [1, 5, 7],
           [2, 4, 6],
           [2, 4, 7],
           [2, 5, 6],
           [2, 5, 7],
           [3, 4, 6],
           [3, 4, 7],
           [3, 5, 6],
           [3, 5, 7]])

    """

    arrays = [np.asarray(x) for x in arrays]
    dtype = arrays[0].dtype

    n = np.prod([x.size for x in arrays])
    if out is None:
        out = np.zeros([n, len(arrays)], dtype=dtype)

    m = n / arrays[0].size
    out[:,0] = np.repeat(arrays[0], m)
    if arrays[1:]:
        cartesian(arrays[1:], out=out[0:m, 1:])
        for j in xrange(1, arrays[0].size):
            out[j*m:(j+1)*m, 1:] = out[0:m, 1:]
    return out
pv.
fuente
46
¿Alguna vez consideró enviar esto para incluirlo en numpy? Esta no es la primera vez que busco esta funcionalidad y encuentro tu publicación.
endolito
1
Hay un error en esta implementación. Para matrices de cadenas, por ejemplo: matrices [0] .dtype = "| S3" y matrices [1] .dtype = "| S5". Por lo tanto, es necesario encontrar la cadena más larga en la entrada y usar su tipo en out = np.zeros ([n, len (arrays)], dtype = dtype)
norecces
38
FYI: parece haber llegado al paquete scikit-learn enfrom sklearn.utils.extmath import cartesian
Gus
2
Me acabo de dar cuenta: esto es ligeramente diferente de itertools.combinations, ya que esta función respeta el orden de los valores, mientras que las combinaciones no lo hacen, por lo que esta función devuelve más valores que las combinaciones. Sigue siendo muy impresionante, pero desafortunadamente no es lo que estaba buscando :(
David Marx
66
TypeError: slice indices must be integers or None or have an __index__ methodlanzado porcartesian(arrays[1:], out=out[0:m,1:])
Boern
36

itertools.combinations es en general la forma más rápida de obtener combinaciones de un contenedor de Python (si de hecho quiere combinaciones, es decir, arreglos SIN repeticiones e independientes del orden; eso no es lo que parece estar haciendo su código, pero no puedo diga si eso se debe a que su código tiene errores o porque está usando la terminología incorrecta).

Si desea algo diferente a las combinaciones, quizás otros iteradores en itertools, producto permutations, podrían servirle mejor. Por ejemplo, parece que su código es aproximadamente el mismo que:

for val in itertools.product(np.arange(0, 1, 0.1), repeat=6):
    print F(val)

Todos estos iteradores producen tuplas, no listas o matrices numpy, por lo que si su F es exigente para obtener específicamente una matriz numpy, tendrá que aceptar la sobrecarga adicional de construir o borrar y volver a llenar una en cada paso.

Alex Martelli
fuente
8

Puedes hacer algo como esto

import numpy as np

def cartesian_coord(*arrays):
    grid = np.meshgrid(*arrays)        
    coord_list = [entry.ravel() for entry in grid]
    points = np.vstack(coord_list).T
    return points

a = np.arange(4)  # fake data
print(cartesian_coord(*6*[a])

lo que da

array([[0, 0, 0, 0, 0, 0],
   [0, 0, 0, 0, 0, 1],
   [0, 0, 0, 0, 0, 2],
   ..., 
   [3, 3, 3, 3, 3, 1],
   [3, 3, 3, 3, 3, 2],
   [3, 3, 3, 3, 3, 3]])
felippe
fuente
2
¿Hay alguna manera de hacer que NumPy acepte más de 32 matrices para meshgrid? Este método funciona para mí siempre que no pase más de 32 matrices.
Joelmob
8

La siguiente implementación numpy debería ser de aprox. 2 veces la velocidad de la respuesta dada:

def cartesian2(arrays):
    arrays = [np.asarray(a) for a in arrays]
    shape = (len(x) for x in arrays)

    ix = np.indices(shape, dtype=int)
    ix = ix.reshape(len(arrays), -1).T

    for n, arr in enumerate(arrays):
        ix[:, n] = arrays[n][ix[:, n]]

    return ix
Stefan van der Walt
fuente
1
Se ve bien. Según mis pruebas rudimentarias, esto parece más rápido que la respuesta original para todos los pares, triples y 4-tuplas de {1,2, ..., 100}. Después de eso, la respuesta original gana. Además, para futuros lectores que busquen generar todas las k-tuplas de {1, ..., n}, funcionará np.indices((n,...,n)).reshape(k,-1).T.
jme
Esto solo funciona para enteros, mientras que la respuesta aceptada también funciona para flotantes.
FJC
7

Parece que desea una cuadrícula para evaluar su función, en cuyo caso puede usar numpy.ogrid(abrir) o numpy.mgrid(desarrollar):

import numpy
my_grid = numpy.mgrid[[slice(0,1,0.1)]*6]
steabert
fuente
6

puedes usar np.array(itertools.product(a, b))

William Song
fuente
np.array (list (itertools.product (l, l2)))
ZirconCode
4

Aquí hay otra manera, usando NumPy puro, sin recursión, sin comprensión de la lista y sin explícito para los bucles. Es aproximadamente un 20% más lento que la respuesta original, y se basa en np.meshgrid.

def cartesian(*arrays):
    mesh = np.meshgrid(*arrays)  # standard numpy meshgrid
    dim = len(mesh)  # number of dimensions
    elements = mesh[0].size  # number of elements, any index will do
    flat = np.concatenate(mesh).ravel()  # flatten the whole meshgrid
    reshape = np.reshape(flat, (dim, elements)).T  # reshape and transpose
    return reshape

Por ejemplo,

x = np.arange(3)
a = cartesian(x, x, x, x, x)
print(a)

da

[[0 0 0 0 0]
 [0 0 0 0 1]
 [0 0 0 0 2]
 ..., 
 [2 2 2 2 0]
 [2 2 2 2 1]
 [2 2 2 2 2]]
étale-cohomology
fuente
3

Para una implementación puramente numpy del producto cartesiano de matrices 1D (o listas planas de python), simplemente use meshgrid(), enrolle los ejes con transpose()y vuelva a dar forma a la salida deseada:

 def cartprod(*arrays):
     N = len(arrays)
     return transpose(meshgrid(*arrays, indexing='ij'), 
                      roll(arange(N + 1), -1)).reshape(-1, N)

Tenga en cuenta que esto tiene la convención de que el último eje cambia más rápido ("estilo C" o "fila mayor").

In [88]: cartprod([1,2,3], [4,8], [100, 200, 300, 400], [-5, -4])
Out[88]: 
array([[  1,   4, 100,  -5],
       [  1,   4, 100,  -4],
       [  1,   4, 200,  -5],
       [  1,   4, 200,  -4],
       [  1,   4, 300,  -5],
       [  1,   4, 300,  -4],
       [  1,   4, 400,  -5],
       [  1,   4, 400,  -4],
       [  1,   8, 100,  -5],
       [  1,   8, 100,  -4],
       [  1,   8, 200,  -5],
       [  1,   8, 200,  -4],
       [  1,   8, 300,  -5],
       [  1,   8, 300,  -4],
       [  1,   8, 400,  -5],
       [  1,   8, 400,  -4],
       [  2,   4, 100,  -5],
       [  2,   4, 100,  -4],
       [  2,   4, 200,  -5],
       [  2,   4, 200,  -4],
       [  2,   4, 300,  -5],
       [  2,   4, 300,  -4],
       [  2,   4, 400,  -5],
       [  2,   4, 400,  -4],
       [  2,   8, 100,  -5],
       [  2,   8, 100,  -4],
       [  2,   8, 200,  -5],
       [  2,   8, 200,  -4],
       [  2,   8, 300,  -5],
       [  2,   8, 300,  -4],
       [  2,   8, 400,  -5],
       [  2,   8, 400,  -4],
       [  3,   4, 100,  -5],
       [  3,   4, 100,  -4],
       [  3,   4, 200,  -5],
       [  3,   4, 200,  -4],
       [  3,   4, 300,  -5],
       [  3,   4, 300,  -4],
       [  3,   4, 400,  -5],
       [  3,   4, 400,  -4],
       [  3,   8, 100,  -5],
       [  3,   8, 100,  -4],
       [  3,   8, 200,  -5],
       [  3,   8, 200,  -4],
       [  3,   8, 300,  -5],
       [  3,   8, 300,  -4],
       [  3,   8, 400,  -5],
       [  3,   8, 400,  -4]])

Si desea cambiar el primer eje más rápido ("estilo FORTRAN" o "columna mayor"), simplemente cambie el orderparámetro de reshape()esta manera:reshape((-1, N), order='F')

RBF06
fuente
1

Pandas mergeofrece una solución ingenua y rápida al problema:

# given the lists
x, y, z = [1, 2, 3], [4, 5], [6, 7]

# get dfs with same, constant index 
x = pd.DataFrame({'x': x}, index=np.repeat(0, len(x))
y = pd.DataFrame({'y': y}, index=np.repeat(0, len(y))
z = pd.DataFrame({'z': z}, index=np.repeat(0, len(z))

# get all permutations stored in a new df
df = pd.merge(x, pd.merge(y, z, left_index=True, righ_index=True),
              left_index=True, right_index=True)
simone
fuente