Używanie numpy do zbudowania tablicy wszystkich kombinacji dwóch tablic

Próbuję przeanalizować przestrzeń parametrów funkcji 6-parametrowej, aby zbadać jej zachowanie numeryczne, zanim spróbuję zrobić z nią coś złożonego, więc szukam skutecznego sposobu, aby to zrobić.

Moja funkcja przyjmuje wartości zmiennoprzecinkowe, które jako dane wejściowe przyjmują 6-dim tablicę numpy. Na początku próbowałem zrobić tak:

Najpierw stworzyłem funkcję, która pobiera 2 tablice i generuje tablicę ze wszystkimi kombinacjami wartości z dwóch tablic

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

Następnie reduce()stosowałem to do m kopii tej samej tablicy:

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

A potem oceniam moją funkcję w ten sposób:

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

To działa, ale jest zbyt wolne. Wiem, że przestrzeń parametrów jest ogromna, ale nie powinno to być takie wolne. W tym przykładzie pobrałem tylko 10 ⁶ (milion) punktów, a samo utworzenie tablicy zajęło więcej niż 15 sekund values.

Czy znasz lepszy sposób na zrobienie tego z numpy?

Mogę zmodyfikować sposób, w jaki funkcja Fprzyjmuje argumenty, jeśli jest to konieczne.

W nowszej wersji numpy(> 1.8.x) numpy.meshgrid()zapewnia znacznie szybszą implementację:

Rozwiązanie @ 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()Kiedyś był tylko 2D, teraz jest zdolny do ND. W tym przypadku 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]])

Zwróć uwagę, że kolejność wyniku końcowego jest nieco inna.

Oto prosta implementacja. To około 5 razy szybciej niż przy użyciu 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

itertools.combinations to generalnie najszybszy sposób na uzyskanie kombinacji z kontenera Pythona (jeśli faktycznie potrzebujesz kombinacji, tj. układów BEZ powtórzeń i niezależnych od kolejności; wydaje się, że nie to robi twój kod, ale nie mogę powiedz, czy dzieje się tak dlatego, że twój kod zawiera błędy, czy też dlatego, że używasz złej terminologii).

Jeśli chcesz czegoś innego niż kombinacje, być może inne iteratory w itertools productlub permutationsmogą ci służyć lepiej. Na przykład wygląda na to, że Twój kod jest mniej więcej taki sam, jak:

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

Wszystkie te iteratory dają krotki, a nie listy lub tablice numpy, więc jeśli twój F jest wybredny, jeśli chodzi o uzyskanie konkretnej tablicy numpy, będziesz musiał zaakceptować dodatkowe obciążenie związane z konstruowaniem lub czyszczeniem i ponownym wypełnianiem jednej na każdym kroku.

Możesz zrobić coś takiego

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

co daje

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

Następująca implementacja numpy powinna wynosić ok. 2x szybkość udzielonej odpowiedzi:

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

Wygląda na to, że potrzebujesz siatki do oceny funkcji, w takim przypadku możesz użyć numpy.ogrid(otwórz) lub numpy.mgrid(uzupełnij):

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

możesz użyć np.array(itertools.product(a, b))

Oto jeszcze jeden sposób, używając czystego NumPy, bez rekursji, bez rozumienia list i bez jawnych pętli for. Jest około 20% wolniejsze niż oryginalna odpowiedź i opiera się na 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

Na przykład,

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

daje

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

Aby uzyskać czystą numpy implementację iloczynu kartezjańskiego tablic 1D (lub płaskich list Pythona), po prostu użyj meshgrid(), obróć osie transpose()i zmień kształt na pożądany wynik:

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

Zauważ, że ma to konwencję najszybszej zmiany ostatniej osi („styl C” lub „wiersz główny”).

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

Jeśli chcesz zmienić pierwszą oś najszybciej („styl FORTRAN” lub „kolumna główna”), po prostu zmień orderparametr w reshape()ten sposób:reshape((-1, N), order='F')

Pandy mergeoferują naiwne, szybkie rozwiązanie problemu:

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