Programming

x 및 y 배열 포인트의 직교 곱을 2D 포인트의 단일 배열로

procodes 2020. 6. 29. 21:21
반응형

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.productnumpy 배열을 사용 하고 다시 돌려주는 것이 가장 효율적인 형태가 아니라고 가정합니다.


>>> 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를 사용하여 두 배열의 모든 조합으로 구성된 배열참조하십시오 .


정식 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.

주목할만한 대안

연속적인 메모리 블록을 포트란 순서로 쓰는 것이 더 빠릅니다. 이것이이 대안의 기초이며 cartesian_product_transpose일부 하드웨어에서보다 빠른 것으로 입증되었습니다 cartesian_product(아래 참조). 그러나 동일한 원칙을 사용하는 Paul Panzer의 답변은 훨씬 빠릅니다. 여전히 관심있는 독자들을 위해 여기에 포함시킵니다.

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

다른 하드웨어 및 운영 체제를 사용하는 사용자는 다른 결과를 볼 수 있습니다. 예를 들어, unutbu 보고서는 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

현재 응답 사용을 허용 tile하고 repeat두 배열 함께 방송한다. 그러나 meshgrid기능은 실질적으로 동일합니다. 다음의 출력의 tilerepeat전치에 전달되기 전에는 :

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+ dstackvs. 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)
...     

적당한 크기의 입력의 경우 속도가 크게 향상되었습니다. 그러나 numpy최신 컴퓨터에서 최신 버전의 Python (3.6.1) 및 (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위.) 비교하기 위해 cartesiandstack_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이 답변의 맨 처음에 정의 된 것보다 느립니다.이 질문에 대한 답변 중 가장 빠른 구현보다 약간 느립니다.


파이썬에서 정상적인 목록 이해를 할 수 있습니다.

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

두 배열의 경우 (고전 사례) :

여기에 이미지 설명을 입력하십시오

4 개의 배열 :

여기에 이미지 설명을 입력하십시오

(여기서 배열의 길이는 수십 항목에 불과합니다.)


줄거리를 재현하는 코드 :

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_ppcartesian_product_transpose다른 전략을 모두 사용하는 @senderle과 달리 - cartesion_product더 유리한 전치 메모리 레이아웃을 사용 하는 버전과 약간의 최적화가 있습니다.
  • cartesian_product_pp원래 메모리 레이아웃을 고수합니다. 빠른 복사는 연속 복사를 사용하는 것입니다. 연속 복사는 훨씬 더 빠른 것으로 판명되었으므로 메모리의 일부만 유효한 데이터를 포함하더라도 전체 메모리 블록을 복사하는 것이 유효한 비트 만 복사하는 것이 좋습니다.

일부 퍼 플롯. C와 Fortran 레이아웃에 대해 별도의 작업을 만들었습니다.이 작업은 서로 다른 작업 IMO이기 때문입니다.

'pp'로 끝나는 이름은 나의 접근 방식입니다.

1) 많은 작은 요소들 (각 2 요소)

여기에 이미지 설명을 입력하십시오여기에 이미지 설명을 입력하십시오

2) 많은 작은 요소 (각 4 개 요소)

여기에 이미지 설명을 입력하십시오여기에 이미지 설명을 입력하십시오

3) 길이가 같은 3 가지 요소

여기에 이미지 설명을 입력하십시오여기에 이미지 설명을 입력하십시오

4) 길이가 같은 두 가지 요인

여기에 이미지 설명을 입력하십시오여기에 이미지 설명을 입력하십시오

코드 (각 플롯마다 별도의 실행이 필요합니다. 재설정 방법을 알 수 없었으며 적절하게 편집 / 주석 입력 / 출력해야합니다) :

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)

axis=-1매개 변수 에 유의하십시오 . 결과에서 마지막 (가장 안쪽) 축입니다. 사용하는 것과 같습니다 axis=ndim.

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


@kennytm 답변 을 잠시 동안 사용 했지만 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 배우기 패키지는 정확히의 빠른 구현이 :

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

가장 빠른 방법은 생성기 표현식을 map 함수와 결합하는 것입니다.

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 matrix product로 해결하는 것입니다.

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.000155s


각 쌍에 추가와 같은 간단한 작업을 수행해야하는 특정 경우에는 추가 차원을 도입하고 브로드 캐스트가 작업을 수행하도록 할 수 있습니다.

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

실제로 쌍 자체를 얻는 비슷한 방법이 있는지 확실하지 않습니다.

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

반응형