4

这是一个小代码，它演示了我得到的错误。

import numpy as np

r=4.0
L=20.0
ratio = 4*np.pi / 3.0 * (r/L)**3

for i in range(5, 16):
    n = 10**i
    m = int(ratio * n)
    print i,n,m

    ip = np.random.random_integers(100, size=(n,3))    
    jp = np.random.random_integers(100, size=(m,3))

    a = np.expand_dims(ip, -1) == jp.T
    b = np.where( a.all(axis=1).any(axis=1) )[0]

我得到以下输出：

5 100000 3351
6 1000000 33510 
Traceback (most recent call last):
  File "example.py", line 16, in <module>
    b = np.where( a.all(axis=1).any(axis=1) )[0]
AttributeError: 'bool' object has no attribute 'all'

有人知道这里发生了什么吗？

或者，一种在 ip 中索引 jp 元素位置的相当快速的方法也可以工作。我可以从这里选择第二种解决方案

4

1 回答 1

1

您正在广播ip反对jp创建非常大的数组。当i==6你有一个 100GB 的阵列时。

一种解决方案是遍历数组：

for i in range(2,6):
    t=time.time()
    n = 10**i+1
    m = int(ratio * n)
    print i,n,m

    ip = np.random.random_integers(10, size=(n,3))
    jp = np.random.random_integers(10, size=(m,3))

    chunksize=10000
    if chunksize>ip.shape[0]:
        chunksize=ip.shape[0]
    span=ip.shape[0]/chunksize
    remainder=(ip.shape[0]-span*chunksize)

    out=[]
    start=0
    for n in xrange(span):
        end=start+chunksize
        a = np.expand_dims(ip[start:end], -1) == jp.T
        b = np.where( a.all(axis=1).any(axis=1) )[0]
        out.append(b+start)
        start+=chunksize

    if remainder!=0:
        a = np.expand_dims(ip[-remainder:], -1) == jp.T
        b = np.where( a.all(axis=1).any(axis=1) )[0]
        out.append(b+end)

    end=np.sort(np.concatenate(out))
    print time.time()-t,end.shape

时间约为 10 秒，i==6因此i==7大约需要 20 分钟。

于 2013-08-14T23:59:06.827 回答