構建列表中的點的3D立方體

Question

我有一個列表pts包含N點（Python浮點數）。我希望構造尺寸的NumPy的陣列N*N*N*3使得陣列等價於：

for i in xrange(0, N): 
    for j in xrange(0, N): 
     for k in xrange(0, N): 
      arr[i,j,k,0] = pts[i] 
      arr[i,j,k,1] = pts[j] 
      arr[i,j,k,2] = pts[k] 

我想知道我怎麼能利用與NumPy以及諸如tile功能的陣列廣播規則進行簡化。

Answer 1

我認爲以下應工作：

pts = np.array(pts) #Skip if pts is a numpy array already 
lp = len(pts) 
arr = np.zeros((lp,lp,lp,3)) 
arr[:,:,:,0] = pts[:,None,None] #None is the same as np.newaxis 
arr[:,:,:,1] = pts[None,:,None] 
arr[:,:,:,2] = pts[None,None,:] 

簡單的測試：

import numpy as np 
import timeit 

def meth1(pts): 
    pts = np.array(pts) #Skip if pts is a numpy array already 
    lp = len(pts) 
    arr = np.zeros((lp,lp,lp,3)) 
    arr[:,:,:,0] = pts[:,None,None] #None is the same as np.newaxis 
    arr[:,:,:,1] = pts[None,:,None] 
    arr[:,:,:,2] = pts[None,None,:] 
    return arr 

def meth2(pts): 
    lp = len(pts) 
    N = lp 
    arr = np.zeros((lp,lp,lp,3)) 
    for i in xrange(0, N): 
     for j in xrange(0, N): 
     for k in xrange(0, N): 
      arr[i,j,k,0] = pts[i] 
      arr[i,j,k,1] = pts[j] 
      arr[i,j,k,2] = pts[k] 

    return arr 

pts = range(10) 
a1 = meth1(pts) 
a2 = meth2(pts) 

print np.all(a1 == a2) 

NREPEAT = 10000 
print timeit.timeit('meth1(pts)','from __main__ import meth1,pts',number=NREPEAT) 
print timeit.timeit('meth2(pts)','from __main__ import meth2,pts',number=NREPEAT) 

結果：

True 
0.873255968094 #my way 
11.4249279499 #original 

所以這種新方法是一個數量級的速度更快以及。

Answer 2

import numpy as np 
N = 10 
pts = xrange(0,N) 
l = [ [ [ [ pts[i],pts[j],pts[k] ] for k in xrange(0,N) ] for j in xrange(0,N) ] for i in xrange(0,N) ] 
x = np.array(l, np.int32) 
print x.shape # (10,10,10,3) 

Answer 3

這可以在兩條線進行：

def meth3(pts): 
    arrs = np.broadcast_arrays(*np.ix_(pts, pts, pts)) 
    return np.concatenate([a[...,None] for a in arrs], axis=3) 

然而，這種方法是不一樣快mgilson的答案，因爲concatenate是煩人緩慢。儘管他的答案的一般化版本也大致如此，並且可以爲任何一組數組生成你想要的結果（即，包含在n維網格中的n維笛卡爾積）。

def meth4(arrs):  # or meth4(*arrs) for a simplified interface 
    arr = np.empty([len(a) for a in arrs] + [len(arrs)]) 
    for i, a in enumerate(np.ix_(*arrs)): 
     arr[...,i] = a 
    return arr 

該接受序列中的任何序列，只要它可以被轉換成numpy的陣列序列：

>>> meth4([[0, 1], [2, 3]]) 
array([[[ 0., 2.], 
     [ 0., 3.]], 

     [[ 1., 2.], 
     [ 1., 3.]]]) 

而這個一般性的成本不太高 - 它只是兩次小pts數組作爲慢：

>>> (meth4([pts, pts, pts]) == meth1(pts)).all() 
True 
>>> %timeit meth4([pts, pts, pts]) 
10000 loops, best of 3: 27.4 us per loop 
>>> %timeit meth1(pts) 
100000 loops, best of 3: 13.1 us per loop 

，它實際上有點快是爲較大的（雖然速度增加可能是由於我使用的empty我而不是zeros）：

>>> pts = np.linspace(0, 1, 100) 
>>> %timeit meth4([pts, pts, pts]) 
100 loops, best of 3: 13.4 ms per loop 
>>> %timeit meth1(pts) 
100 loops, best of 3: 16.7 ms per loop