programing

x 및 y 배열 점을 2D 점의 단일 배열로 데카르트 곱

telecom 2023. 6. 14. 21:40
반응형

x 및 y 배열 점을 2D 점의 단일 배열로 데카르트 곱

저는 격자의 x축과 y축을 정의하는 두 개의 numpy 배열을 가지고 있습니다.예:

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

이 어레이의 데카르트 제품을 생성하여 다음을 생성하고 싶습니다.

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

여러 번 반복해야 하기 때문에 비효율적이지는 않습니다.파이썬 목록으로 변환하여 사용하는 것으로 가정합니다.itertools.product그리고 numpy 배열로 돌아가는 것은 가장 효율적인 형태가 아닙니다.

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_은 이적인방식로으의문; 에된화사용서반례면▁in▁of.ix_ 인덱스를 배열로 생성하는 입니다. 마침 같은 모양의 배열을 브로드캐스트 할당에 사용할 수 있습니다.사용해 볼 수 있도록 영감을 준 mgilson에게 감사드립니다.ix_이런 식으로, 그리고 사용할 제안을 포함하여 이 답변에 대해 매우 도움이 되는 피드백을 제공한 unutbu에게.

주목할 만한 대안

때때로 연속적인 메모리 블록을 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

방식을 하게 된 접 방 식 이 해 을 한 후 근 의 나 그 의 의 썼 간 것 는 만 거 습 니 새 단 한 버 다 을 전 로 운 의 르 고 큼 빠 거 ▁simple ▁as ▁panzer ▁as 썼 ▁after ▁almost 습 ▁a ▁that s ▁version 니 다 ▁and ▁i ' ▁is ▁new , ' 을 ▁coming , 버 전 , 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).

다음 섹션에는 다른 대안에 대한 몇 가지 테스트가 포함되어 있습니다.이것들은 이제 다소 구식이지만, 중복적인 노력보다는 역사적인 관심에서 그것들을 여기에 남겨두기로고 결정했습니다.최신 테스트는 Nico Schlömer의 답변뿐만 아니라 Panzer의 답변을 참조하십시오.

대안에 대한 테스트

다음은 이러한 기능 중 일부가 여러 가지 대안과 비교하여 제공하는 성능 향상을 보여주는 일련의 테스트입니다.에 표시된 10 3및 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 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)]

테스트 결과:

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또한. (이 경우에 오류가 발생하는 이유를 확인할 수 있을 때까지 이 테스트에서 제거했습니다.)

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입력 어레이의 수가 4개 이상으로 증가할 때까지 경쟁력을 유지합니다.다음에 그이로후.cartesian_product_transpose약간의 가장자리가 있습니다.

다른 하드웨어 및 운영 체제를 사용하는 사용자는 다른 결과를 볼 수도 있음을 거듭 강조할 필요가 있습니다.를 들어.04, 3. 및 Ubuntu 14.04, Python 3.4.3을 합니다.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

아래에서는 이전에 수행한 테스트에 대한 몇 가지 세부 사항을 참조하십시오.이러한 접근 방식의 상대적인 성능은 시간이 지남에 따라 다양한 하드웨어 및 버전의 Python 및numpy의 최신 버전을 사용하는 사용자에게 즉시 유용하지는 않습니다.numpy이 답변의 첫 번째 버전 이후 상황이 어떻게 변했는지 설명합니다.

대안은 다음과 같습니다.meshgrid+dstack

은 현수락답변다사음용다니합을은된재다▁uses니▁the를 사용합니다.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])]

하지만, 이 시점에서 재구성하는 것보다, 우리는 출력을 전달할 수 있습니다.meshgriddstack나중에 모양을 변경하여 작업을 절약할 수 있습니다.

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(에서는 이전전버 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)
...     

중간 크기의 입력에서는 속도가 크게 향상되었습니다.하지만 저는 더 최신 버전의 파이썬(3.6.1)으로 이러한 테스트를 다시 시도했습니다.numpy 최신에서. (1..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에서 이것들이 상호 교환 가능하다는 것을 시사합니다.

일반화된 제품 기능

일반적으로, 우리는 작은 입력의 경우 내장된 기능을 사용하는 것이 더 빠를 것으로 예상할 수 있지만, 큰 입력의 경우 목적에 맞게 만들어진 기능이 더 빠를 수 있습니다. 곱의 , 또한일차된화제경우의품원반차,tile그리고.repeat도움이 되지 않을 겁니다. 왜냐하면 그들은 명확한 고차원 아날로그를 가지고 있지 않기 때문입니다.따라서 목적에 맞게 구축된 기능의 동작도 조사할 가치가 있습니다.

는 이, 이전 및 대분의관테는에이시나부작다, 버서입니 에서 수행된 몇 가지 .numpy비교를 위하여

cartesian더 큰 입력에 대해 상당히 잘 수행하기 위해 사용되는 다른 답변에 정의된 함수.(이 함수는 다음과 같습니다.)cartesian_product_recursive아래.)비교하기 는 위하기해교.cartesiandstack_prodct우리는 단지 2차원을 사용합니다.

여기서도 이전 테스트는 유의한 차이를 보인 반면 새로운 테스트는 거의 없음을 보여줍니다.

이전 테스트

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

이러한 차이는 흥미롭고 기록할 가치가 있다고 생각합니다. 하지만 결국에는 학술적인 것입니다., 이 버전들은 의대답시있테보이듯들이었여주스트다느, 이모버은거다다립니음보항상의들전이든작에는▁are▁than▁as다느▁almost니▁slower립▁all▁always▁of▁versions▁showed▁at다▁these,▁answer▁the보다 느립니다.cartesian_product이 답변의 맨 앞에 정의되어 있습니다. 이 질문에 대한 답변 중 가장 빠른 구현보다 그 자체가 약간 느립니다.

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

N 배열의 데카르트 곱을 계산하기 위한 일반 솔루션에 대한 두 배열의 모든 조합 배열을 구성하려면 numpy 사용을 참조하십시오.

그냥 파이썬으로 일반적인 목록 이해를 할 수 있습니다.

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의 답변보다 약간 더 명확한 성능 비교를 했습니다.

두 개의 어레이(일반적인 경우):

enter image description here

4개 어레이의 경우:

enter image description here

(여기서는 어레이 길이가 수십 개 항목에 불과합니다.)


그림을 재현하는 코드:

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

@senderle의 모범적인 기초 작업을 바탕으로 C용과 Fortran 레이아웃용의 두 가지 버전을 생각해 냈습니다. 이 두 가지 버전은 종종 조금 더 빠릅니다.

  • cartesian_product_transpose_pp@freelle의 경우와 다릅니다.cartesian_product_transpose완전히 다른 전략을 사용하는 - 의 버전.cartesion_product보다 유리한 전치 메모리 레이아웃 + 일부 매우 사소한 최적화를 사용합니다.
  • cartesian_product_pp원래 메모리 레이아웃이 있는 스틱.그것을 빠르게 만드는 것은 연속적인 복사를 사용하기 때문입니다.연속 복사는 메모리의 일부에만 유효한 데이터가 포함되어 있더라도 전체 블록을 복사하는 것이 유효한 비트만 복사하는 것보다 훨씬 빠르다는 것이 밝혀졌습니다.

퍼플 플롯.C와 Fortran 레이아웃은 IMO 작업이 다르기 때문에 따로 만들었습니다.

'pp'로 끝나는 이름들이 저의 접근법입니다.

많은 작은 요인(각각 2개의 요소)

enter image description hereenter image description here

많은 작은 요인(각각 4개의 요소)

enter image description hereenter image description here

길이가 같은 세 가지 요인

enter image description hereenter image description here

길이가 같은 두 인자

enter image description hereenter image description here

코드(각 플롯 b/c에 대해 별도의 런을 수행해야 함 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년 10월 기준으로 numpy는 현재 제네릭을 보유하고 있습니다.np.stack축 매개 변수를 사용하는 함수입니다.이를 사용하면 "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))

a = np.array([1,2])
b = np.array([10,20])
cartesian_product(a,b)

# output:
# array([[ 1, 10],
#        [ 2, 10],
#        [ 1, 20],
#        [ 2, 20]])  

에 대한 axis=-1은 결과의 안쪽) .이 축은 결과의 마지막 축입니다.은 사하는것같다니습을 사용하는 것과 .axis=ndim.

또 다른 의견은, 데카르트 제품이 매우 빠르게 폭발하기 때문에, 어떤 이유로 메모리의 어레이를 실현할 필요가 없는 한, 제품이 매우 크면, 우리는 사용하고 싶을 수도 있습니다.itertools값을 즉시 사용할 수 있습니다.

Scikit-learn 패키지에는 다음과 같은 기능이 신속하게 구현되어 있습니다.

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

출력 순서에 관심이 있는 경우 이 구현의 관례는 원하는 것과 다릅니다.정확한 주문은 다음과 같습니다.

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

한동안 @kennytm answer를 사용했지만 TensorFlow에서 동일한 작업을 수행하려고 할 때 TensorFlow에 해당하는 작업이 없습니다.numpy.repeat()약간의 실험 후, 저는 임의의 점 벡터에 대한 더 일반적인 해결책을 찾았다고 생각합니다.

뚱딴지의 경우:

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)

더 일반적으로, 두 개의 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

이중 발생기 식에 대해 설명합니다.

실제로 원하는 것이 각 좌표 쌍의 실제 곱을 계산하는 것이라면 가장 빠른 방법은 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()))

출력:

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

이것은 또한 iter 도구를 사용하여 쉽게 할 수 있습니다.제품 방법

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

저는 그들 스스로 실제로 그들을 얻을 수 있는 비슷한 방법이 있는지 확신할 수 없습니다.

저는 파티에 좀 늦었지만, 그 문제의 까다로운 변종을 만났습니다.여러 어레이의 데카르트 제품을 원하지만 데카르트 제품은 컴퓨터의 메모리보다 훨씬 큽니다(그러나 해당 제품으로 수행되는 계산은 빠르거나 최소한 병렬 처리가 가능합니다).

분명한 해결책은 이 데카르트 제품을 청크로 나누고 이러한 청크를 차례로 처리하는 것입니다("스트리밍" 방식).로 쉽게 할 수 있습니다.itertools.product끔찍할 정도로 느립니다.또한, 여기서 제안된 솔루션 중 어느 것도(만큼 빠른) 이러한 가능성을 제공하지 않습니다.은 Numba를 " 제가를는보솔규하안루, "적례다다빠"보다 약간 더 .cartesian_product가 할 있는 모든 에 꽤 제가 할 수 있는 모든 곳에서 최적화하려고 노력했기 때문에 꽤 긴 시간입니다.

import numba as nb
import numpy as np
from typing import List


@nb.njit(nb.types.Tuple((nb.int32[:, :],
                         nb.int32[:]))(nb.int32[:],
                                       nb.int32[:],
                                       nb.int64, nb.int64))
def cproduct(sizes: np.ndarray, current_tuple: np.ndarray, start_idx: int, end_idx: int):
    """Generates ids tuples from start_id to end_id"""
    assert len(sizes) >= 2
    assert start_idx < end_idx

    tuples = np.zeros((end_idx - start_idx, len(sizes)), dtype=np.int32)
    tuple_idx = 0
    # stores the current combination
    current_tuple = current_tuple.copy()
    while tuple_idx < end_idx - start_idx:
        tuples[tuple_idx] = current_tuple
        current_tuple[0] += 1
        # using a condition here instead of including this in the inner loop
        # to gain a bit of speed: this is going to be tested each iteration,
        # and starting a loop to have it end right away is a bit silly
        if current_tuple[0] == sizes[0]:
            # the reset to 0 and subsequent increment amount to carrying
            # the number to the higher "power"
            current_tuple[0] = 0
            current_tuple[1] += 1
            for i in range(1, len(sizes) - 1):
                if current_tuple[i] == sizes[i]:
                    # same as before, but in a loop, since this is going
                    # to get called less often
                    current_tuple[i + 1] += 1
                    current_tuple[i] = 0
                else:
                    break
        tuple_idx += 1
    return tuples, current_tuple


def chunked_cartesian_product_ids(sizes: List[int], chunk_size: int):
    """Just generates chunks of the cartesian product of the ids of each
    input arrays (thus, we just need their sizes here, not the actual arrays)"""
    prod = np.prod(sizes)

    # putting the largest number at the front to more efficiently make use
    # of the cproduct numba function
    sizes = np.array(sizes, dtype=np.int32)
    sorted_idx = np.argsort(sizes)[::-1]
    sizes = sizes[sorted_idx]
    if chunk_size > prod:
        chunk_bounds = (np.array([0, prod])).astype(np.int64)
    else:
        num_chunks = np.maximum(np.ceil(prod / chunk_size), 2).astype(np.int32)
        chunk_bounds = (np.arange(num_chunks + 1) * chunk_size).astype(np.int64)
        chunk_bounds[-1] = prod
    current_tuple = np.zeros(len(sizes), dtype=np.int32)
    for start_idx, end_idx in zip(chunk_bounds[:-1], chunk_bounds[1:]):
        tuples, current_tuple = cproduct(sizes, current_tuple, start_idx, end_idx)
        # re-arrange columns to match the original order of the sizes list
        # before yielding
        yield tuples[:, np.argsort(sorted_idx)]


def chunked_cartesian_product(*arrays, chunk_size=2 ** 25):
    """Returns chunks of the full cartesian product, with arrays of shape
    (chunk_size, n_arrays). The last chunk will obviously have the size of the
    remainder"""
    array_lengths = [len(array) for array in arrays]
    for array_ids_chunk in chunked_cartesian_product_ids(array_lengths, chunk_size):
        slices_lists = [arrays[i][array_ids_chunk[:, i]] for i in range(len(arrays))]
        yield np.vstack(slices_lists).swapaxes(0,1)


def cartesian_product(*arrays):
    """Actual cartesian product, not chunked, still fast"""
    total_prod = np.prod([len(array) for array in arrays])
    return next(chunked_cartesian_product(*arrays, total_prod))


a = np.arange(0, 3)
b = np.arange(8, 10)
c = np.arange(13, 16)
for cartesian_tuples in chunked_cartesian_product(*[a, b, c], chunk_size=5):
    print(cartesian_tuples)

이렇게 하면 데카르트 제품을 53배로 출력할 수 있습니다.

[[ 0  8 13]
 [ 0  8 14]
 [ 0  8 15]
 [ 1  8 13]
 [ 1  8 14]]
[[ 1  8 15]
 [ 2  8 13]
 [ 2  8 14]
 [ 2  8 15]
 [ 0  9 13]]
[[ 0  9 14]
 [ 0  9 15]
 [ 1  9 13]
 [ 1  9 14]
 [ 1  9 15]]
[[ 2  9 13]
 [ 2  9 14]
 [ 2  9 15]]

만약 당신이 여기서 무엇이 행해지고 있는지 이해할 의향이 있다면, 그 이면에 있는 직관.njitted함수는 각 "숫자"를 입력 배열의 크기(정규 이진법, 십진법 또는 십진법 기본에서 동일한 숫자로 구성된 요소)로 구성된 이상한 숫자 베이스로 열거하는 것입니다.

분명히, 이 솔루션은 대형 제품에 흥미롭습니다.작은 것들의 경우, 오버헤드가 약간 비쌀 수 있습니다.

참고: numba는 아직 개발이 많이 진행 중이기 때문에 python 3.6으로 이를 실행하기 위해 numba 0.50을 사용하고 있습니다.

또 하나:

>>>x1, y1 = np.meshgrid(x, y)
>>>np.c_[x1.ravel(), y1.ravel()]
array([[1, 4],
       [2, 4],
       [3, 4],
       [1, 5],
       [2, 5],
       [3, 5]])

Ashkan의 대답에 영감을 받아 다음을 시도할 수도 있습니다.

>>> x, y = np.meshgrid(x, y)
>>> np.concatenate([x.flatten().reshape(-1,1), y.flatten().reshape(-1,1)], axis=1)

이것은 당신에게 필요한 데카르트 제품을 줄 것입니다!

입니다. (을 사용한 곱)numpy.tile그리고.numpy.repeat함수)를 선택합니다.

from functors import reduce
from operator import mul

def cartesian_product(arrays):
    return np.vstack(
        np.tile(
            np.repeat(arrays[j], reduce(mul, map(len, arrays[j+1:]), 1)),
            reduce(mul, map(len, arrays[:j]), 1),
        )
        for j in range(len(arrays))
    ).T

PyTorch를 사용할 의향이 있다면 매우 효율적이라고 생각합니다.

>>> import torch

>>> torch.cartesian_prod(torch.as_tensor(x), torch.as_tensor(y))
tensor([[1, 4],
        [1, 5],
        [2, 4],
        [2, 5],
        [3, 4],
        [3, 5]])

그리고 당신은 쉽게 얻을 수 있습니다.numpy매개 변수:

>>> torch.cartesian_prod(torch.as_tensor(x), torch.as_tensor(y)).numpy()
array([[1, 4],
       [1, 5],
       [2, 4],
       [2, 5],
       [3, 4],
       [3, 5]])

언급URL : https://stackoverflow.com/questions/11144513/cartesian-product-of-x-and-y-array-points-into-single-array-of-2d-points

반응형