У меня есть два массива numpy, которые определяют оси x и y сетки. Например:
x = numpy.array([1,2,3])
y = numpy.array([4,5])
Я хотел бы сгенерировать декартово произведение этих массивов для генерации:
array([[1,4],[2,4],[3,4],[1,5],[2,5],[3,5]])
В некотором смысле это не ужасно неэффективно, так как мне нужно делать это много раз в цикле. Я предполагаю, что преобразование их в список Python и использование itertools.product и обратно в массив numpy не является самой эффективной формой.
>>> 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]])
См. Использование numpy для построения массива всех комбинаций из двух массивов для общего решения для вычисления декартова произведения N массивов.
cartesian_product (почти)Существует много подходов к этой проблеме с разными свойствами. Некоторые быстрее, чем другие, а некоторые более универсальны. После большого тестирования и настройки, я обнаружил, что следующая функция, которая вычисляет n-мерный cartesian_product, быстрее, чем большинство других для многих входов. Для пары подходов, которые немного сложнее, но во многих случаях даже немного быстрее, см. Ответ Paul Panzer.
Учитывая этот ответ, это уже не самая быстрая реализация декартова продукта в numpy, о котором я знаю. Однако, я думаю, его простота будет продолжать делать его полезным ориентиром для дальнейшего улучшения:
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)
Стоит отметить, что эта функция использует ix_ необычным образом; тогда как документированное использование ix_ заключается в генерировать индексы в массив, так бывает, что массивы с одинаковой формой могут использоваться для трансляции назначение. Большое спасибо mgilson, который вдохновил меня попробовать использовать ix_ таким образом, и unutbu, который предоставил некоторые чрезвычайно полезные отзывы об этом ответе, включая предложение использовать numpy.result_type.
Иногда быстрее записывать непрерывные блоки памяти в порядке Fortran. Это основа этой альтернативы, cartesian_product_transpose, которая оказалась быстрее на каком-то оборудовании, чем cartesian_product (см. Ниже). Однако ответ Павла Панзера, который использует тот же принцип, еще быстрее. Тем не менее, я включаю это здесь для заинтересованных читателей:
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
После того, как я понял подход Panzer, я написал новую версию, которая почти так же быстро, как и его, и почти такая же простая, как 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
Похоже, что накладные расходы постоянны, что заставляет его работать медленнее, чем Panzer для небольших входов. Но для больших входных данных во всех тестах, которые я выполнял, он работает так же хорошо, как и его самая быстрая реализация (cartesian_product_transpose_pp).
В следующих разделах я включаю некоторые тесты других альтернатив. Сейчас они несколько устарели, но вместо того, чтобы дублировать усилия, я решил оставить их здесь из исторических интересов. Для получения последних тестов см. Ответ Panzer, а также Nico Schlömer.
Вот батарея тестов, которые показывают повышение производительности, которое некоторые из этих функций предоставляют по сравнению с несколькими альтернативами. Все тесты, показанные здесь, выполнялись на четырехъядерной машине, работающей под управлением Mac OS 10.12.5, Python 3.6.1 и numpy 1.12.1. Известно, что вариации на аппаратное и программное обеспечение дают разные результаты, поэтому YMMV. Запустите эти тесты для себя, чтобы быть уверенным!
Определения:
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 /questions/64378/using-numpy-to-build-an-array-of-all-combinations-of-two-arrays/441983#441983
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)]
Результаты тестирования:
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)
Во всех случаях cartesian_product, как определено в начале этого ответа, является самым быстрым.
Для тех функций, которые принимают произвольное количество входных массивов, стоит проверить производительность при len(arrays) > 2. (Пока я не могу определить, почему cartesian_product_recursive выдает ошибку в этом случае, я удалил ее из этих тестов.)
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)
Как показывают эти тесты, cartesian_product остается конкурентоспособным до тех пор, пока количество входных массивов не будет превышать (примерно) четыре. После этого cartesian_product_transpose имеет небольшое ребро.
Стоит повторить, что пользователи с другим оборудованием и операционными системами могут видеть разные результаты. Например, unutbu сообщает о следующих результатах для этих тестов с использованием Ubuntu 14.04, Python 3.4.3 и numpy 1.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
Ниже я расскажу о нескольких более ранних тестах, которые я провел в этих строках. Относительная производительность этих подходов со временем изменилась, для разных аппаратных и различных версий Python и numpy. Хотя это не сразу полезно для людей, использующих современные версии numpy, это иллюстрирует, как все изменилось со времени первой версии этого ответа.
meshgrid + dstackВ принятом в настоящее время ответе используются tile и repeat для трансляции двух массивов. Но функция meshgrid делает практически одно и то же. Здесь вывод tile и repeat перед передачей в транспонирование:
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])]
И здесь вывод meshgrid:
In [4]: numpy.meshgrid(x, y)
Out[4]:
[array([[1, 2, 3],
[1, 2, 3]]), array([[4, 4, 4],
[5, 5, 5]])]
Как вы можете видеть, он почти идентичен. Нам нужно только изменить результат, чтобы получить точно такой же результат.
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])]
Вместо того, чтобы переделывать в этот момент, мы могли бы передать вывод meshgrid в dstack и впоследствии изменить его, что экономит некоторую работу:
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]])
В отличие от претензии в этом комментарии, я не видел никаких доказательств того, что разные входы будут вызывать результаты по-разному, и, как показано выше, они очень похожи вещи, так что было бы довольно странно, если бы они это сделали. Пожалуйста, дайте мне знать, если вы найдете контрпример.
meshgrid + dstack против repeat + transposeОтносительная производительность этих двух подходов со временем изменилась. В более ранней версии Python (2.7) результат с использованием meshgrid + dstack был заметно быстрее для небольших входов. (Обратите внимание, что эти тесты взяты из старой версии этого ответа.) Определения:
>>> 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)
...
Для ввода среднего размера я видел значительное ускорение. Но я повторил эти тесты с более поздними версиями Python (3.6.1) и numpy (1.12.1) на более новой машине. Эти два подхода почти идентичны.
Старый тест
>>> 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
Новый тест
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)
Как всегда, YMMV, но это говорит о том, что в последних версиях Python и numpy они взаимозаменяемы.
В общем, мы можем ожидать, что использование встроенных функций будет быстрее для небольших входов, в то время как для больших входов целевая функция может быть быстрее. Кроме того, для обобщенного n-мерного произведения tile и repeat не помогут, поскольку у них нет четких высокоразмерных аналогов. Поэтому стоит исследовать поведение специально созданных функций.
Большинство соответствующих тестов появляются в начале этого ответа, но вот несколько тестов, выполненных для более ранних версий Python и numpy для сравнения.
Функция cartesian, определенная в другом ответе, использовалась для выполнения больших задач. (Это то же самое, что и функция, названная выше cartesian_product_recursive.) Чтобы сравнить cartesian с dstack_prodct, мы используем только два измерения.
Здесь снова старый тест показал значительную разницу, в то время как новый тест почти не отображается.
Старый тест
>>> 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
Новый тест
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)
Как и раньше, dstack_product по-прежнему бьет cartesian в меньших масштабах.
Новый тест (избыточный старый тест не показан)
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)
Эти различия, я думаю, интересны и достойны записи; но в конце они академичны. Как показали тесты в начале этого ответа, все эти версии почти всегда медленнее, чем cartesian_product, определенные в самом начале этого ответа, что само по себе немного медленнее, чем самые быстрые реализации среди ответов на этот вопрос.
Вы можете просто выполнить обычное понимание списка в python
x = numpy.array([1,2,3])
y = numpy.array([4,5])
[[x0, y0] for x0 in x for y0 in y]
который должен дать вам
[[1, 4], [1, 5], [2, 4], [2, 5], [3, 4], [3, 5]]
Я также заинтересовался этим и провел небольшое сравнение производительности, возможно, несколько понятнее, чем в ответе @senderle.
Для двух массивов (классический случай):
Для четырех массивов:
(Обратите внимание, что длина массивов здесь составляет всего несколько десятков записей.)
Код для воспроизведения участков:
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: 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
)
Основываясь на образцовой работе на основе @senderle, я придумал две версии: одну для C и одну для компоновки Fortran, - которые часто бывают немного быстрее.
cartesian_product_transpose_pp - в отличие от @senderle cartesian_product_transpose который использует совсем другую стратегию - версию cartesion_product которая использует более благоприятную компоновку памяти транспонирования + некоторые очень незначительные оптимизации.cartesian_product_pp использует оригинальную макет памяти. То, что делает это быстро, - это использование непрерывного копирования. Смежные копии оказываются намного быстрее, чем копирование полного блока памяти, хотя только часть из них содержит достоверные данные, предпочтительнее только копировать действительные биты.Некоторые перфорации. Я сделал отдельные для макетов C и Fortran, потому что это разные задачи IMO.
Имена, оканчивающиеся на "pp", - мои подходы.
1) множество крошечных факторов (по 2 элемента)
2) множество мелких факторов (по 4 элемента)
3) три фактора равной длины
4) два фактора равной длины
Код (нужно делать отдельные прогоны для каждого графика b/c. Я не мог понять, как сбросить, а также отредактировать/прокомментировать в/из соответственно):
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
)
По состоянию на октябрь 2017 года numpy теперь имеет общую функцию np.stack которая принимает параметр оси. Используя его, мы можем иметь "обобщенное декартово произведение" с использованием техники "dstack и meshgrid":
import numpy as np
def cartesian_product(*arrays):
ndim = len(arrays)
return np.stack(np.meshgrid(*arrays), axis=-1).reshape(-1, ndim)
Обратите внимание на параметр axis=-1. Это последняя (внутренняя) ось в результате. Это эквивалентно использованию axis=ndim.
Еще один комментарий, так как декартовы продукты взрываются очень быстро, если нам не нужно каким-то образом реализовать массив в памяти, если продукт очень большой, мы можем захотеть использовать itertools и использовать значения "на лету".
Я использовал @kennytm answer какое-то время, но при попытке сделать то же самое в TensorFlow, но я обнаружил, что TensorFlow не имеет эквивалента numpy.repeat(). После небольшого эксперимента, я думаю, я нашел более общее решение для произвольных векторов точек.
Для 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)
и для 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)
Пакет Scikit-learn имеет быструю реализацию именно этого:
from sklearn.utils.extmath import cartesian
product = cartesian((x,y))
Обратите внимание, что соглашение этой реализации отличается от того, что вы хотите, если вам нужен порядок вывода. Для вашего точного заказа вы можете сделать
product = cartesian((y,x))[:, ::-1]
В более общем плане, если у вас есть два массива 2d numpy a и b, и вы хотите объединить каждую строку a в каждую строку b (декартово произведение строк, вроде объединения в базу данных), вы можете используйте этот метод:
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))
Самое быстрое, что вы можете получить, это либо комбинировать выражение генератора с функцией карты:
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()))
Выходы (фактически распечатывается весь итоговый список):
[(0, 0), (1, 0), ...,(998, 199), (999, 199)]
execution time: 1.253567 s
или с использованием выражения с двойным генератором:
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()))
Выходы (напечатан весь список):
[(0, 0), (1, 0), ...,(998, 199), (999, 199)]
execution time: 1.187415 s
Учтите, что большая часть времени вычислений входит в команду печати. Расчеты генератора в остальном прилично эффективны. Без печати время вычисления:
execution time: 0.079208 s
для выражения генератора + функция отображения и:
execution time: 0.007093 s
для выражения двойного генератора.
Если вы действительно хотите, чтобы вычислить фактический продукт каждой из пар координат, наиболее быстрым является его решение в виде матричного продукта:
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()))
Выходы:
[[ 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
и без печати (в этом случае он мало экономит, поскольку на самом деле распечатывается только малая часть матрицы):
execution time: 0.003083 s
Это также легко сделать с помощью метода 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')
Результат: массив ([
[1, 4],
[1, 5],
[2, 4],
[2, 5],
[3, 4],
[3, 5]], dtype = int32)
Время выполнения: 0,000155 с
В конкретном случае, когда вам нужно выполнить простые операции, такие как сложение для каждой пары, вы можете ввести дополнительное измерение и позволить вещанию делать свою работу:
>>> 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]])
Я не уверен, есть ли подобный способ на самом деле получить сами пары.