在numpy中执行多个均值计算

Question

早上好，

我正在用Numpy实现一个Cressman滤波器，用来做距离加权平均。我使用了一种叫做Ball Tree的实现（感谢Jake VanderPlas），它可以为请求数组中的每个点返回一系列位置。我的查询数组（q）的形状是[n,3]，每个点都有x、y、z坐标，我想对存储在树中的点做加权平均。包裹在树周围的代码会返回一定距离内的点，所以我得到了长度不一的数组。

我使用where函数来找到非空的条目（也就是在影响半径内至少有一些点的位置），从而创建了isgood数组……

然后，我对所有查询点进行循环，返回self.z的加权平均值（注意，这可以是dims=1或dims=2，以允许多个共同网格）。

所以，使用map或其他更快的方法的复杂之处在于self.distances和self.locations中数组长度的不均匀性……我对numpy/python还不太熟悉，但我想不出一种方法可以做到数组级别的操作（也就是说，不想回到使用循环）。

self.locations, self.distances = self.tree.query_radius( q, r, return_distance=True)
t2=time()
if debug: print "Removing voids"
isgood=np.where( np.array([len(x) for x in self.locations])!=0)[0]
interpol = np.zeros( (len(self.locations),) + np.shape(self.z[0]) )
interpol.fill(np.nan)
for dist, ix, posn, roi in zip(self.distances[isgood], self.locations[isgood], isgood, r[isgood]):
    interpol[isgood[jinterpol]] = np.average(self.z[ix], weights=(roi**2-dist**2) / (roi**2 + dist**2), axis=0)
    jinterpol += 1

那么……有没有什么提示可以加速这个循环呢？

对于一个典型的映射，比如将天气雷达数据从范围、方位、高度网格映射到笛卡尔网格，我有240x240x34个点和4个变量，查询树的时间是99秒（这个树是Jake用C和cython写的……这是个难点，因为你需要搜索数据！），计算的时间是100秒……在我看来，这个速度太慢了？我的瓶颈在哪里？np.mean效率高吗？因为它被调用了数百万次，是否可以在这里获得加速？如果我使用float32而不是默认的float64，会有帮助吗……或者甚至将其缩放为整数（这在加权时会很难避免溢出……任何提示都非常感谢！

Answer 1

你可以在这里找到关于Cressman方案和高斯权重函数的优缺点讨论：

http://www.flame.org/~cdoswell/publications/radar_oa_00.pdf

关键在于要把平滑参数和数据匹配起来（我建议使用一个接近数据点平均间距的值）。一旦你确定了平滑参数，就可以设置一个“影响半径”，这个半径是权重函数降到0.01（或者其他值）时的半径。

速度有多重要？如果你愿意，可以不调用指数函数来计算权重，而是为一些固定的半径增量制作一个离散的权重表，这样计算速度会快很多。理想情况下，你应该有一些在网格边界外的数据，这些数据可以用来映射网格点周围的值（甚至是网格的边界点）。需要注意的是，这并不是一个真正的插值方案——它不会准确返回数据点的观测值。就像Cressman方案一样，它是一种低通滤波器。