Produit cartésien des points du tableau x et y en un seul tableau de points 2D

J'ai deux tableaux numpy qui définissent les axes x et y d'une grille. Par exemple:

x = numpy.array([1,2,3])
y = numpy.array([4,5])

Je voudrais générer le produit cartésien de ces tableaux pour générer:

array([[1,4],[2,4],[3,4],[1,5],[2,5],[3,5]])

D'une manière qui n'est pas terriblement inefficace puisque je dois le faire plusieurs fois en boucle. Je suppose que les convertir en une liste Python et utiliser itertools.productet revenir en un tableau numpy n'est pas la forme la plus efficace.

>>> numpy.transpose([numpy.tile(x, len(y)), numpy.repeat(y, len(x))])
array([[1, 4],
       [2, 4],
       [3, 4],
       [1, 5],
       [2, 5],
       [3, 5]])

Voir Utilisation de numpy pour créer un tableau de toutes les combinaisons de deux tableaux pour une solution générale pour le calcul du produit cartésien de N tableaux.

Un canonique cartesian_product(presque)

Il existe de nombreuses approches à ce problème avec des propriétés différentes. Certains sont plus rapides que d'autres et certains sont plus polyvalents. Après de nombreux tests et ajustements, j'ai constaté que la fonction suivante, qui calcule une n-dimension cartesian_product, est plus rapide que la plupart des autres pour de nombreuses entrées. Pour une paire d'approches qui sont légèrement plus complexes, mais qui sont même un peu plus rapides dans de nombreux cas, voir la réponse de Paul Panzer .

Compte tenu de cette réponse, ce n'est plus l' implémentation la plus rapide du produit cartésien à numpyma connaissance. Cependant, je pense que sa simplicité continuera à en faire une référence utile pour les améliorations futures:

def cartesian_product(*arrays):
    la = len(arrays)
    dtype = numpy.result_type(*arrays)
    arr = numpy.empty([len(a) for a in arrays] + [la], dtype=dtype)
    for i, a in enumerate(numpy.ix_(*arrays)):
        arr[...,i] = a
    return arr.reshape(-1, la)

Il convient de mentionner que cette fonction utilise ix_de manière inhabituelle; alors que l'utilisation documentée de ix_consiste à générer des indices dans un tableau, il se trouve que des tableaux de la même forme peuvent être utilisés pour l'affectation diffusée. Un grand merci à mgilson , qui m'a inspiré à essayer d'utiliser ix_cette méthode, et à unutbu , qui a fourni des commentaires extrêmement utiles sur cette réponse, y compris la suggestion à utiliser numpy.result_type.

Alternatives notables

Il est parfois plus rapide d'écrire des blocs de mémoire contigus dans l'ordre Fortran. C'est la base de cette alternative, cartesian_product_transposequi s'est avérée plus rapide sur certains matériels que cartesian_product(voir ci-dessous). Cependant, la réponse de Paul Panzer, qui utilise le même principe, est encore plus rapide. Pourtant, j'inclus ceci ici pour les lecteurs intéressés:

def cartesian_product_transpose(*arrays):
    broadcastable = numpy.ix_(*arrays)
    broadcasted = numpy.broadcast_arrays(*broadcastable)
    rows, cols = numpy.prod(broadcasted[0].shape), len(broadcasted)
    dtype = numpy.result_type(*arrays)

    out = numpy.empty(rows * cols, dtype=dtype)
    start, end = 0, rows
    for a in broadcasted:
        out[start:end] = a.reshape(-1)
        start, end = end, end + rows
    return out.reshape(cols, rows).T

Après avoir compris l'approche de Panzer, j'ai écrit une nouvelle version presque aussi rapide que la sienne, et presque aussi simple que cartesian_product:

def cartesian_product_simple_transpose(arrays):
    la = len(arrays)
    dtype = numpy.result_type(*arrays)
    arr = numpy.empty([la] + [len(a) for a in arrays], dtype=dtype)
    for i, a in enumerate(numpy.ix_(*arrays)):
        arr[i, ...] = a
    return arr.reshape(la, -1).T

Cela semble avoir une surcharge en temps constant qui le rend plus lent que celui de Panzer pour les petites entrées. Mais pour des entrées plus importantes, dans tous les tests que j'ai exécutés, il fonctionne aussi bien que son implémentation la plus rapide ( cartesian_product_transpose_pp).

Dans les sections suivantes, j'inclus quelques tests d'autres alternatives. Celles-ci sont maintenant quelque peu dépassées, mais plutôt que de dupliquer les efforts, j'ai décidé de les laisser ici par intérêt historique. Pour des tests à jour, voir la réponse de Panzer, ainsi que celle de Nico Schlömer .

Tests par rapport aux alternatives

Voici une batterie de tests qui montrent l'amélioration des performances que certaines de ces fonctions fournissent par rapport à un certain nombre d'alternatives. Tous les tests présentés ici ont été effectués sur une machine quadricœur, exécutant Mac OS 10.12.5, Python 3.6.1 et numpy1.12.1. Les variations sur le matériel et les logiciels sont connues pour produire des résultats différents, donc YMMV. Exécutez ces tests pour vous-même pour être sûr!

Définitions:

import numpy
import itertools
from functools import reduce

### Two-dimensional products ###

def repeat_product(x, y):
    return numpy.transpose([numpy.tile(x, len(y)), 
                            numpy.repeat(y, len(x))])

def dstack_product(x, y):
    return numpy.dstack(numpy.meshgrid(x, y)).reshape(-1, 2)

### Generalized N-dimensional products ###

def cartesian_product(*arrays):
    la = len(arrays)
    dtype = numpy.result_type(*arrays)
    arr = numpy.empty([len(a) for a in arrays] + [la], dtype=dtype)
    for i, a in enumerate(numpy.ix_(*arrays)):
        arr[...,i] = a
    return arr.reshape(-1, la)

def cartesian_product_transpose(*arrays):
    broadcastable = numpy.ix_(*arrays)
    broadcasted = numpy.broadcast_arrays(*broadcastable)
    rows, cols = numpy.prod(broadcasted[0].shape), len(broadcasted)
    dtype = numpy.result_type(*arrays)

    out = numpy.empty(rows * cols, dtype=dtype)
    start, end = 0, rows
    for a in broadcasted:
        out[start:end] = a.reshape(-1)
        start, end = end, end + rows
    return out.reshape(cols, rows).T

# from https://stackoverflow.com/a/1235363/577088

def cartesian_product_recursive(*arrays, out=None):
    arrays = [numpy.asarray(x) for x in arrays]
    dtype = arrays[0].dtype

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

    m = n // arrays[0].size
    out[:,0] = numpy.repeat(arrays[0], m)
    if arrays[1:]:
        cartesian_product_recursive(arrays[1:], out=out[0:m,1:])
        for j in range(1, arrays[0].size):
            out[j*m:(j+1)*m,1:] = out[0:m,1:]
    return out

def cartesian_product_itertools(*arrays):
    return numpy.array(list(itertools.product(*arrays)))

### Test code ###

name_func = [('repeat_product',                                                 
              repeat_product),                                                  
             ('dstack_product',                                                 
              dstack_product),                                                  
             ('cartesian_product',                                              
              cartesian_product),                                               
             ('cartesian_product_transpose',                                    
              cartesian_product_transpose),                                     
             ('cartesian_product_recursive',                           
              cartesian_product_recursive),                            
             ('cartesian_product_itertools',                                    
              cartesian_product_itertools)]

def test(in_arrays, test_funcs):
    global func
    global arrays
    arrays = in_arrays
    for name, func in test_funcs:
        print('{}:'.format(name))
        %timeit func(*arrays)

def test_all(*in_arrays):
    test(in_arrays, name_func)

# `cartesian_product_recursive` throws an 
# unexpected error when used on more than
# two input arrays, so for now I've removed
# it from these tests.

def test_cartesian(*in_arrays):
    test(in_arrays, name_func[2:4] + name_func[-1:])

x10 = [numpy.arange(10)]
x50 = [numpy.arange(50)]
x100 = [numpy.arange(100)]
x500 = [numpy.arange(500)]
x1000 = [numpy.arange(1000)]

Résultats de test:

In [2]: test_all(*(x100 * 2))
repeat_product:
67.5 µs ± 633 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
dstack_product:
67.7 µs ± 1.09 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
cartesian_product:
33.4 µs ± 558 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
cartesian_product_transpose:
67.7 µs ± 932 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
cartesian_product_recursive:
215 µs ± 6.01 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
cartesian_product_itertools:
3.65 ms ± 38.7 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

In [3]: test_all(*(x500 * 2))
repeat_product:
1.31 ms ± 9.28 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
dstack_product:
1.27 ms ± 7.5 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
cartesian_product:
375 µs ± 4.5 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
cartesian_product_transpose:
488 µs ± 8.88 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
cartesian_product_recursive:
2.21 ms ± 38.4 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
cartesian_product_itertools:
105 ms ± 1.17 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

In [4]: test_all(*(x1000 * 2))
repeat_product:
10.2 ms ± 132 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
dstack_product:
12 ms ± 120 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
cartesian_product:
4.75 ms ± 57.1 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
cartesian_product_transpose:
7.76 ms ± 52.7 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
cartesian_product_recursive:
13 ms ± 209 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
cartesian_product_itertools:
422 ms ± 7.77 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

Dans tous les cas, cartesian_producttel que défini au début de cette réponse est le plus rapide.

Pour les fonctions qui acceptent un nombre arbitraire de tableaux d'entrée, cela vaut également la peine de vérifier les performances len(arrays) > 2. (Jusqu'à ce que je puisse déterminer pourquoi cartesian_product_recursivegénère une erreur dans ce cas, je l'ai supprimée de ces tests.)

In [5]: test_cartesian(*(x100 * 3))
cartesian_product:
8.8 ms ± 138 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
cartesian_product_transpose:
7.87 ms ± 91.5 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
cartesian_product_itertools:
518 ms ± 5.5 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

In [6]: test_cartesian(*(x50 * 4))
cartesian_product:
169 ms ± 5.1 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
cartesian_product_transpose:
184 ms ± 4.32 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
cartesian_product_itertools:
3.69 s ± 73.5 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

In [7]: test_cartesian(*(x10 * 6))
cartesian_product:
26.5 ms ± 449 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
cartesian_product_transpose:
16 ms ± 133 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
cartesian_product_itertools:
728 ms ± 16 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

In [8]: test_cartesian(*(x10 * 7))
cartesian_product:
650 ms ± 8.14 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
cartesian_product_transpose:
518 ms ± 7.09 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
cartesian_product_itertools:
8.13 s ± 122 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

Comme le montrent ces tests, cartesian_productreste compétitif jusqu'à ce que le nombre de tableaux d'entrée dépasse (environ) quatre. Après cela, cartesian_product_transposea un léger avantage.

Il convient de rappeler que les utilisateurs disposant d'autres matériels et systèmes d'exploitation peuvent voir des résultats différents. Par exemple, unutbu rapporte avoir vu les résultats suivants pour ces tests utilisant Ubuntu 14.04, Python 3.4.3 et numpy1.14.0.dev0 + b7050a9:

>>> %timeit cartesian_product_transpose(x500, y500) 
1000 loops, best of 3: 682 µs per loop
>>> %timeit cartesian_product(x500, y500)
1000 loops, best of 3: 1.55 ms per loop

Ci-dessous, je vais dans quelques détails sur les tests précédents que j'ai exécutés dans ce sens. Les performances relatives de ces approches ont changé au fil du temps, pour différents matériels et différentes versions de Python et numpy. Bien que ce ne soit pas immédiatement utile pour les personnes utilisant des versions à jour de numpy, cela illustre comment les choses ont changé depuis la première version de cette réponse.

Une alternative simple: meshgrid+dstack

La réponse actuellement acceptée utilise tileet repeatpour diffuser deux tableaux ensemble. Mais la meshgridfonction fait pratiquement la même chose. Voici la sortie de tileet repeatavant d'être passée à transposer:

In [1]: import numpy
In [2]: x = numpy.array([1,2,3])
   ...: y = numpy.array([4,5])
   ...: 

In [3]: [numpy.tile(x, len(y)), numpy.repeat(y, len(x))]
Out[3]: [array([1, 2, 3, 1, 2, 3]), array([4, 4, 4, 5, 5, 5])]

Et voici la sortie de meshgrid:

In [4]: numpy.meshgrid(x, y)
Out[4]: 
[array([[1, 2, 3],
        [1, 2, 3]]), array([[4, 4, 4],
        [5, 5, 5]])]

Comme vous pouvez le voir, c'est presque identique. Il suffit de remodeler le résultat pour obtenir exactement le même résultat.

In [5]: xt, xr = numpy.meshgrid(x, y)
   ...: [xt.ravel(), xr.ravel()]
Out[5]: [array([1, 2, 3, 1, 2, 3]), array([4, 4, 4, 5, 5, 5])]

Plutôt que de remodeler à ce stade, cependant, nous pourrions passer la sortie de meshgridà dstacket remodeler ensuite, ce qui économise du travail:

In [6]: numpy.dstack(numpy.meshgrid(x, y)).reshape(-1, 2)
Out[6]: 
array([[1, 4],
       [2, 4],
       [3, 4],
       [1, 5],
       [2, 5],
       [3, 5]])

Contrairement à ce que prétend ce commentaire , je n'ai vu aucune preuve que différentes entrées produiraient des sorties de forme différente, et comme le montre ce qui précède, elles font des choses très similaires, il serait donc assez étrange qu'elles le fassent. Veuillez me faire savoir si vous trouvez un contre-exemple.

Test meshgrid+ dstackvs repeat+transpose

La performance relative de ces deux approches a évolué au fil du temps. Dans une version antérieure de Python (2.7), le résultat en utilisant meshgrid+ dstackétait nettement plus rapide pour les petites entrées. (Notez que ces tests proviennent d'une ancienne version de cette réponse.) Définitions:

>>> def repeat_product(x, y):
...     return numpy.transpose([numpy.tile(x, len(y)), 
                                numpy.repeat(y, len(x))])
...
>>> def dstack_product(x, y):
...     return numpy.dstack(numpy.meshgrid(x, y)).reshape(-1, 2)
...     

Pour une entrée de taille moyenne, j'ai vu une accélération significative. Mais j'ai réessayé ces tests avec des versions plus récentes de Python (3.6.1) et numpy(1.12.1), sur une machine plus récente. Les deux approches sont désormais presque identiques.

Ancien test

>>> x, y = numpy.arange(500), numpy.arange(500)
>>> %timeit repeat_product(x, y)
10 loops, best of 3: 62 ms per loop
>>> %timeit dstack_product(x, y)
100 loops, best of 3: 12.2 ms per loop

Nouveau test

In [7]: x, y = numpy.arange(500), numpy.arange(500)
In [8]: %timeit repeat_product(x, y)
1.32 ms ± 24.7 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
In [9]: %timeit dstack_product(x, y)
1.26 ms ± 8.47 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

Comme toujours, YMMV, mais cela suggère que dans les versions récentes de Python et de numpy, ceux-ci sont interchangeables.

Fonctions produit généralisées

En général, on peut s'attendre à ce que l'utilisation des fonctions intégrées soit plus rapide pour les petites entrées, tandis que pour les grandes entrées, une fonction spécialement conçue peut être plus rapide. De plus, pour un produit généralisé à n dimensions, tileet repeatn'aidera pas, car ils n'ont pas d'analogues clairs de dimension supérieure. Il vaut donc la peine d'étudier le comportement des fonctions spécialement conçues à cet effet.

La plupart des tests pertinents apparaissent au début de cette réponse, mais voici quelques-uns des tests effectués sur les versions antérieures de Python et numpyà titre de comparaison.

La cartesianfonction définie dans une autre réponse fonctionnait plutôt bien pour des entrées plus importantes. (C'est la même que la fonction appelée cartesian_product_recursiveci - dessus.) Afin de comparer cartesianà dstack_prodct, nous utilisons seulement deux dimensions.

Là encore, l'ancien test montrait une différence significative, tandis que le nouveau test n'en montrait presque aucune.

Ancien test

>>> x, y = numpy.arange(1000), numpy.arange(1000)
>>> %timeit cartesian([x, y])
10 loops, best of 3: 25.4 ms per loop
>>> %timeit dstack_product(x, y)
10 loops, best of 3: 66.6 ms per loop

Nouveau test

In [10]: x, y = numpy.arange(1000), numpy.arange(1000)
In [11]: %timeit cartesian([x, y])
12.1 ms ± 199 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
In [12]: %timeit dstack_product(x, y)
12.7 ms ± 334 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

Comme avant, dstack_productbat toujours cartesianà des échelles plus petites.

Nouveau test ( ancien test redondant non illustré )

In [13]: x, y = numpy.arange(100), numpy.arange(100)
In [14]: %timeit cartesian([x, y])
215 µs ± 4.75 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
In [15]: %timeit dstack_product(x, y)
65.7 µs ± 1.15 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

Ces distinctions sont, je pense, intéressantes et méritent d'être enregistrées; mais ils sont finalement académiques. Comme l'ont montré les tests au début de cette réponse, toutes ces versions sont presque toujours plus lentes que cartesian_product, défini au tout début de cette réponse - qui est elle-même un peu plus lente que les implémentations les plus rapides parmi les réponses à cette question.

Vous pouvez simplement faire une compréhension de liste normale en python

x = numpy.array([1,2,3])
y = numpy.array([4,5])
[[x0, y0] for x0 in x for y0 in y]

qui devrait vous donner

[[1, 4], [1, 5], [2, 4], [2, 5], [3, 4], [3, 5]]

Cela m'intéressait également et j'ai fait une petite comparaison des performances, peut-être un peu plus claire que dans la réponse de @ senderle.

Pour deux tableaux (le cas classique):

Pour quatre baies:

(Notez que la longueur des tableaux n'est que de quelques dizaines d'entrées ici.)

Code pour reproduire les graphiques:

from functools import reduce
import itertools
import numpy
import perfplot


def dstack_product(arrays):
    return numpy.dstack(numpy.meshgrid(*arrays, indexing="ij")).reshape(-1, len(arrays))


# Generalized N-dimensional products
def cartesian_product(arrays):
    la = len(arrays)
    dtype = numpy.find_common_type([a.dtype for a in arrays], [])
    arr = numpy.empty([len(a) for a in arrays] + [la], dtype=dtype)
    for i, a in enumerate(numpy.ix_(*arrays)):
        arr[..., i] = a
    return arr.reshape(-1, la)


def cartesian_product_transpose(arrays):
    broadcastable = numpy.ix_(*arrays)
    broadcasted = numpy.broadcast_arrays(*broadcastable)
    rows, cols = reduce(numpy.multiply, broadcasted[0].shape), len(broadcasted)
    dtype = numpy.find_common_type([a.dtype for a in arrays], [])

    out = numpy.empty(rows * cols, dtype=dtype)
    start, end = 0, rows
    for a in broadcasted:
        out[start:end] = a.reshape(-1)
        start, end = end, end + rows
    return out.reshape(cols, rows).T


# from https://stackoverflow.com/a/1235363/577088
def cartesian_product_recursive(arrays, out=None):
    arrays = [numpy.asarray(x) for x in arrays]
    dtype = arrays[0].dtype

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

    m = n // arrays[0].size
    out[:, 0] = numpy.repeat(arrays[0], m)
    if arrays[1:]:
        cartesian_product_recursive(arrays[1:], out=out[0:m, 1:])
        for j in range(1, arrays[0].size):
            out[j * m : (j + 1) * m, 1:] = out[0:m, 1:]
    return out


def cartesian_product_itertools(arrays):
    return numpy.array(list(itertools.product(*arrays)))


perfplot.show(
    setup=lambda n: 2 * (numpy.arange(n, dtype=float),),
    n_range=[2 ** k for k in range(13)],
    # setup=lambda n: 4 * (numpy.arange(n, dtype=float),),
    # n_range=[2 ** k for k in range(6)],
    kernels=[
        dstack_product,
        cartesian_product,
        cartesian_product_transpose,
        cartesian_product_recursive,
        cartesian_product_itertools,
    ],
    logx=True,
    logy=True,
    xlabel="len(a), len(b)",
    equality_check=None,
)

En me basant sur le travail de terrain exemplaire de @ senderle, j'ai proposé deux versions - une pour les mises en page C et une pour les mises en page Fortran - qui sont souvent un peu plus rapides.

• cartesian_product_transpose_ppest - contrairement à @ senderle cartesian_product_transposequi utilise une stratégie complètement différente - une version de cartesion_productqui utilise la disposition de mémoire de transposition la plus favorable + quelques optimisations très mineures.
• cartesian_product_ppcolle avec la disposition de la mémoire d'origine. Ce qui le rend rapide, c'est son utilisation de la copie contiguë. Les copies contiguës s'avèrent être tellement plus rapides que copier un bloc complet de mémoire même si seulement une partie de celui-ci contient des données valides est préférable à ne copier que les bits valides.

Quelques perfplots. J'en ai fait des séparés pour les mises en page C et Fortran, car ce sont des tâches différentes de l'OMI.

Les noms se terminant par «pp» sont mes approches.

1) de nombreux petits facteurs (2 éléments chacun)

2) de nombreux petits facteurs (4 éléments chacun)

3) trois facteurs d'égale longueur

4) deux facteurs d'égale longueur

Code (besoin de faire des exécutions séparées pour chaque tracé b / c je ne pouvais pas comprendre comment réinitialiser; également besoin d'éditer / commenter correctement):

import numpy
import numpy as np
from functools import reduce
import itertools
import timeit
import perfplot

def dstack_product(arrays):
    return numpy.dstack(
        numpy.meshgrid(*arrays, indexing='ij')
        ).reshape(-1, len(arrays))

def cartesian_product_transpose_pp(arrays):
    la = len(arrays)
    dtype = numpy.result_type(*arrays)
    arr = numpy.empty((la, *map(len, arrays)), dtype=dtype)
    idx = slice(None), *itertools.repeat(None, la)
    for i, a in enumerate(arrays):
        arr[i, ...] = a[idx[:la-i]]
    return arr.reshape(la, -1).T

def cartesian_product(arrays):
    la = len(arrays)
    dtype = numpy.result_type(*arrays)
    arr = numpy.empty([len(a) for a in arrays] + [la], dtype=dtype)
    for i, a in enumerate(numpy.ix_(*arrays)):
        arr[...,i] = a
    return arr.reshape(-1, la)

def cartesian_product_transpose(arrays):
    broadcastable = numpy.ix_(*arrays)
    broadcasted = numpy.broadcast_arrays(*broadcastable)
    rows, cols = numpy.prod(broadcasted[0].shape), len(broadcasted)
    dtype = numpy.result_type(*arrays)

    out = numpy.empty(rows * cols, dtype=dtype)
    start, end = 0, rows
    for a in broadcasted:
        out[start:end] = a.reshape(-1)
        start, end = end, end + rows
    return out.reshape(cols, rows).T

from itertools import accumulate, repeat, chain

def cartesian_product_pp(arrays, out=None):
    la = len(arrays)
    L = *map(len, arrays), la
    dtype = numpy.result_type(*arrays)
    arr = numpy.empty(L, dtype=dtype)
    arrs = *accumulate(chain((arr,), repeat(0, la-1)), np.ndarray.__getitem__),
    idx = slice(None), *itertools.repeat(None, la-1)
    for i in range(la-1, 0, -1):
        arrs[i][..., i] = arrays[i][idx[:la-i]]
        arrs[i-1][1:] = arrs[i]
    arr[..., 0] = arrays[0][idx]
    return arr.reshape(-1, la)

def cartesian_product_itertools(arrays):
    return numpy.array(list(itertools.product(*arrays)))


# from https://stackoverflow.com/a/1235363/577088
def cartesian_product_recursive(arrays, out=None):
    arrays = [numpy.asarray(x) for x in arrays]
    dtype = arrays[0].dtype

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

    m = n // arrays[0].size
    out[:, 0] = numpy.repeat(arrays[0], m)
    if arrays[1:]:
        cartesian_product_recursive(arrays[1:], out=out[0:m, 1:])
        for j in range(1, arrays[0].size):
            out[j*m:(j+1)*m, 1:] = out[0:m, 1:]
    return out

### Test code ###
if False:
  perfplot.save('cp_4el_high.png',
    setup=lambda n: n*(numpy.arange(4, dtype=float),),
                n_range=list(range(6, 11)),
    kernels=[
        dstack_product,
        cartesian_product_recursive,
        cartesian_product,
#        cartesian_product_transpose,
        cartesian_product_pp,
#        cartesian_product_transpose_pp,
        ],
    logx=False,
    logy=True,
    xlabel='#factors',
    equality_check=None
    )
else:
  perfplot.save('cp_2f_T.png',
    setup=lambda n: 2*(numpy.arange(n, dtype=float),),
    n_range=[2**k for k in range(5, 11)],
    kernels=[
#        dstack_product,
#        cartesian_product_recursive,
#        cartesian_product,
        cartesian_product_transpose,
#        cartesian_product_pp,
        cartesian_product_transpose_pp,
        ],
    logx=True,
    logy=True,
    xlabel='length of each factor',
    equality_check=None
    )

Depuis octobre 2017, numpy dispose désormais d'une np.stackfonction générique qui prend un paramètre d'axe. En l'utilisant, on peut avoir un "produit cartésien généralisé" en utilisant la technique "dstack and meshgrid":

import numpy as np
def cartesian_product(*arrays):
    ndim = len(arrays)
    return np.stack(np.meshgrid(*arrays), axis=-1).reshape(-1, ndim)

Remarque sur le axis=-1paramètre. Il s'agit du dernier axe (le plus à l'intérieur) du résultat. Cela équivaut à utiliseraxis=ndim .

Un autre commentaire, puisque les produits cartésiens explosent très rapidement, à moins que nous ayons besoin de réaliser le tableau en mémoire pour une raison quelconque, si le produit est très grand, nous pouvons vouloir utiliser itertoolset utiliser les valeurs à la volée.

J'ai utilisé la réponse @kennytm pendant un certain temps, mais en essayant de faire la même chose dans TensorFlow, mais j'ai trouvé que TensorFlow n'a pas d'équivalent denumpy.repeat() . Après un peu d'expérimentation, je pense avoir trouvé une solution plus générale pour les vecteurs de points arbitraires.

Pour numpy:

import numpy as np

def cartesian_product(*args: np.ndarray) -> np.ndarray:
    """
    Produce the cartesian product of arbitrary length vectors.

    Parameters
    ----------
    np.ndarray args
        vector of points of interest in each dimension

    Returns
    -------
    np.ndarray
        the cartesian product of size [m x n] wherein:
            m = prod([len(a) for a in args])
            n = len(args)
    """
    for i, a in enumerate(args):
        assert a.ndim == 1, "arg {:d} is not rank 1".format(i)
    return np.concatenate([np.reshape(xi, [-1, 1]) for xi in np.meshgrid(*args)], axis=1)

et pour TensorFlow:

import tensorflow as tf

def cartesian_product(*args: tf.Tensor) -> tf.Tensor:
    """
    Produce the cartesian product of arbitrary length vectors.

    Parameters
    ----------
    tf.Tensor args
        vector of points of interest in each dimension

    Returns
    -------
    tf.Tensor
        the cartesian product of size [m x n] wherein:
            m = prod([len(a) for a in args])
            n = len(args)
    """
    for i, a in enumerate(args):
        tf.assert_rank(a, 1, message="arg {:d} is not rank 1".format(i))
    return tf.concat([tf.reshape(xi, [-1, 1]) for xi in tf.meshgrid(*args)], axis=1)

Le package Scikit-learn a une implémentation rapide de exactement ceci:

from sklearn.utils.extmath import cartesian
product = cartesian((x,y))

Notez que la convention de cette implémentation est différente de ce que vous voulez, si vous vous souciez de l'ordre de la sortie. Pour votre commande exacte, vous pouvez faire

product = cartesian((y,x))[:, ::-1]

Plus généralement, si vous avez deux tableaux 2d numpy a et b, et que vous souhaitez concaténer chaque ligne de a à chaque ligne de b (un produit cartésien de lignes, un peu comme une jointure dans une base de données), vous pouvez utiliser cette méthode :

import numpy
def join_2d(a, b):
    assert a.dtype == b.dtype
    a_part = numpy.tile(a, (len(b), 1))
    b_part = numpy.repeat(b, len(a), axis=0)
    return numpy.hstack((a_part, b_part))

Le plus rapide que vous pouvez obtenir est soit de combiner une expression de générateur avec la fonction de carte:

import numpy
import datetime
a = np.arange(1000)
b = np.arange(200)

start = datetime.datetime.now()

foo = (item for sublist in [list(map(lambda x: (x,i),a)) for i in b] for item in sublist)

print (list(foo))

print ('execution time: {} s'.format((datetime.datetime.now() - start).total_seconds()))

Sorties (en fait, toute la liste résultante est imprimée):

[(0, 0), (1, 0), ...,(998, 199), (999, 199)]
execution time: 1.253567 s

ou en utilisant une expression de double générateur:

a = np.arange(1000)
b = np.arange(200)

start = datetime.datetime.now()

foo = ((x,y) for x in a for y in b)

print (list(foo))

print ('execution time: {} s'.format((datetime.datetime.now() - start).total_seconds()))

Sorties (liste complète imprimée):

[(0, 0), (1, 0), ...,(998, 199), (999, 199)]
execution time: 1.187415 s

Tenez compte du fait que la majeure partie du temps de calcul est consacrée à la commande d'impression. Les calculs du générateur sont par ailleurs assez efficaces. Sans impression, les temps de calcul sont:

execution time: 0.079208 s

pour l'expression du générateur + la fonction de carte et:

execution time: 0.007093 s

pour l'expression du double générateur.

Si vous voulez réellement calculer le produit réel de chacune des paires de coordonnées, le plus rapide est de le résoudre comme un produit matriciel numpy:

a = np.arange(1000)
b = np.arange(200)

start = datetime.datetime.now()

foo = np.dot(np.asmatrix([[i,0] for i in a]), np.asmatrix([[i,0] for i in b]).T)

print (foo)

print ('execution time: {} s'.format((datetime.datetime.now() - start).total_seconds()))

Les sorties:

 [[     0      0      0 ...,      0      0      0]
 [     0      1      2 ...,    197    198    199]
 [     0      2      4 ...,    394    396    398]
 ..., 
 [     0    997   1994 ..., 196409 197406 198403]
 [     0    998   1996 ..., 196606 197604 198602]
 [     0    999   1998 ..., 196803 197802 198801]]
execution time: 0.003869 s

et sans impression (dans ce cas, cela ne permet pas d'économiser beaucoup puisque seule une petite partie de la matrice est réellement imprimée):

execution time: 0.003083 s

Cela peut également être facilement fait en utilisant la méthode itertools.product

from itertools import product
import numpy as np

x = np.array([1, 2, 3])
y = np.array([4, 5])
cart_prod = np.array(list(product(*[x, y])),dtype='int32')

Résultat: tableau ([
[1, 4],
[1, 5],
[2, 4],
[2, 5],
[3, 4],
[3, 5]], dtype = int32)

Temps d'exécution: 0,000155 s

Dans le cas spécifique où vous devez effectuer des opérations simples telles que l'ajout sur chaque paire, vous pouvez introduire une dimension supplémentaire et laisser la diffusion faire le travail:

>>> a, b = np.array([1,2,3]), np.array([10,20,30])
>>> a[None,:] + b[:,None]
array([[11, 12, 13],
       [21, 22, 23],
       [31, 32, 33]])

Je ne sais pas s'il existe un moyen similaire d'obtenir les paires elles-mêmes.