网格/where/索引解决方案已经非常快速。我使它快了约65%。虽然这不算太多,但我还是会逐步解释一下:
对于我来说,最简单的方法是将网格中所有3D向量作为一个大的2D 3 x M
数组中的列来处理。 meshgrid
是创建所有组合的正确工具(请注意,需要numpy版本>= 1.7才能使用3D meshgrid),而vstack
+ reshape
则将数据转换为所需的形式。例如:
>>> np.vstack(np.meshgrid(*[np.arange(0, 2)]*3)).reshape(3,-1)
array([[0, 0, 1, 1, 0, 0, 1, 1],
[0, 0, 0, 0, 1, 1, 1, 1],
[0, 1, 0, 1, 0, 1, 0, 1]])
每一列都是一个三维向量。这八个向量中的每一个代表一个
1x1x1
的立方体顶点(在所有维度上步长和长度均为1的三维网格)。
我们将这个数组称为
vectors
(其中包含表示网格中所有点的所有三维向量)。然后,准备一个布尔掩码以选择符合您的mod2条件的那些向量:
mod2bool = np.sum(vectors, axis=0) % 2 == 0
np.sum(vectors, axis=0)
会创建一个包含每个列向量元素总和的1 x M
数组。因此,mod2bool
是一个1 x M
数组,其中每个列向量都有一个布尔值。现在使用这个布尔掩码:
vectorsubset = vectors[:,mod2bool]
这个代码选择所有行(:),并使用布尔索引来过滤列,这两个操作在numpy中非常快。使用原生的numpy方法计算剩余向量的长度:
lengths = np.sqrt(np.sum(vectorsubset**2, axis=0))
这很快 - 但是,scipy.stats.ss
和bottleneck.ss
甚至可以比此更快地执行平方和计算。
按照您的说明转换长度:
with np.errstate(divide='ignore'):
p = (r/lengths)**n
这涉及到有限数除以零,导致输出数组中出现
Inf
。这是完全可以接受的。我们使用numpy的
errstate
上下文管理器来确保这些零除法不会引发异常或运行时警告。
现在将有限元素求和(忽略无穷大),并返回总和:
return np.sum(p[np.isfinite(p)])
我已经实现了以下两种方法。一种完全像刚才解释的那样,另一种涉及到瓶颈的ss
和nansum
函数。我还添加了您的方法进行比较,以及跳过np.where((x*x+y*y+z*z)!=0)
索引的修改版本,而是创建Inf
,最后用isfinite
方式求和。
import sys
import numpy as np
import bottleneck as bn
N = 100
n = 12
r = np.sqrt(2)
x,y,z = np.meshgrid(*[np.arange(-N, N+1)]*3)
gridvectors = np.vstack((x,y,z)).reshape(3, -1)
def measure_time(func):
import time
def modified_func(*args, **kwargs):
t0 = time.time()
result = func(*args, **kwargs)
duration = time.time() - t0
print("%s duration: %.3f s" % (func.__name__, duration))
return result
return modified_func
@measure_time
def method_columnvecs(vectors):
mod2bool = np.sum(vectors, axis=0) % 2 == 0
vectorsubset = vectors[:,mod2bool]
lengths = np.sqrt(np.sum(vectorsubset**2, axis=0))
with np.errstate(divide='ignore'):
p = (r/lengths)**n
return np.sum(p[np.isfinite(p)])
@measure_time
def method_columnvecs_opt(vectors):
mod2bool = bn.nansum(vectors, axis=0) % 2 == 0
lengths = np.sqrt(bn.ss(vectors[:,mod2bool]))
with np.errstate(divide='ignore'):
p = (r/lengths)**n
return bn.nansum(p[np.isfinite(p)])
@measure_time
def method_original(x,y,z):
ind = np.where((x+y+z)%2==0)
x = x[ind]
y = y[ind]
z = z[ind]
ind = np.where((x*x+y*y+z*z)!=0)
x = x[ind]
y = y[ind]
z = z[ind]
p=np.sqrt(x*x+y*y+z*z)/r
return np.sum((1/p)**n)
@measure_time
def method_original_finitesum(x,y,z):
ind = np.where((x+y+z)%2==0)
x = x[ind]
y = y[ind]
z = z[ind]
lengths = np.sqrt(x*x+y*y+z*z)
with np.errstate(divide='ignore'):
p = (r/lengths)**n
return np.sum(p[np.isfinite(p)])
print method_columnvecs(gridvectors)
print method_columnvecs_opt(gridvectors)
print method_original(x,y,z)
print method_original_finitesum(x,y,z)
这是输出结果:
$ python test.py
method_columnvecs duration: 1.295 s
12.1318801965
method_columnvecs_opt duration: 1.162 s
12.1318801965
method_original duration: 1.936 s
12.1318801965
method_original_finitesum duration: 1.714 s
12.1318801965
所有方法产生的结果都是相同的。当使用isfinite
风格求和时,您的方法会稍微快一些。我的方法更快,但我认为这是一种学术性质的练习,而不是重要的改进 :-)
我还有一个问题:您说对于N=3,计算应该产生12。即使您的方法也没有做到这一点。以上所有方法对于N=3都产生12.1317530867的结果。这是预期的吗?
pdist
,هڈ¯èƒ½è؟کو¶‰هڈٹbottleneckçڑ„ss
(ه¹³و–¹ه’Œï¼Œن½ 需è¦په†چو¬،وڈگé«که®ƒçڑ„و½œهٹ›ï¼‰م€‚ - Dr. Jan-Philip GehrckeC
应该等于什么?我认为你可以轻松地使用x,y,z = np.meshgrid(* [np.arange(-N,N + 1)] * 3)
来完成这个问题。 - wflynny