1

我有一个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来简化这一点。

4

3 回答 3

3

我认为以下应该有效：

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

所以这种新方法也快了一个数量级。

于 2012-10-18T12:46:37.513 回答
1 
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)
于 2012-10-18T12:48:59.853 回答
1

这可以通过两行来完成：

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
于 2012-10-18T15:21:30.510 回答