4

我已经编写了以下代码，它完全符合我的要求，但是速度太慢了。我确信有一种方法可以让它更快，但我似乎无法找到它应该如何完成。代码的第一部分只是显示什么是什么形状。

两个测量图像（VV1和HH1）
预先计算的值，VV模拟和HH模拟，它们都依赖于 3 个参数（预先计算的(101, 31, 11)值）
索引 2 只是将VV和HH图像放在同一个 ndarray 中，而不是制作两个 3darray

VV1 = numpy.ndarray((54, 43)).flatten()
HH1 = numpy.ndarray((54, 43)).flatten()
precomp = numpy.ndarray((101, 31, 11, 2))

我们让变化的三个参数中的两个

comp = numpy.zeros((len(parameter1), len(parameter2)))

for i,(vv,hh) in enumerate(zip(VV1,HH1)):
    comp0 = numpy.zeros((len(parameter1),len(parameter2)))
    for j in range(len(parameter1)):
        for jj in range(len(parameter2)):
            comp0[j,jj] = numpy.min((vv-precomp[j,jj,:,0])**2+(hh-precomp[j,jj,:,1])**2)
    comp+=comp0

我知道我应该做的显而易见的事情是尽可能多地摆脱 for 循环，但我不知道如何numpy.min在处理更多维度时使行为正确。

我注意到的第二件事（如果它可以矢量化则不那么重要，但仍然很有趣）我注意到它主要占用 CPU 时间，而不是 RAM，但我已经搜索了很长时间，但我找不到写类似“parfor”的方法" 而不是 matlab 中的 "for"，（@parallel如果我只是将 for 循环放在单独的方法中，是否可以制作装饰器？）

编辑：回复 Janne Karila：是的，这肯定会改善很多，

for (vv,hh) in zip(VV1,HH1):
    comp+= numpy.min((vv-precomp[...,0])**2+(hh-precomp[...,1])**2, axis=2)

肯定要快很多，但是有没有可能删除外部 for 循环呢？有没有办法用一个@parallel或什么东西来使一个for循环并行？

4

3 回答 3

5

这可以替换内部循环，j并且jj

comp0 = numpy.min((vv-precomp[...,0])**2+(hh-precomp[...,1])**2, axis=2)

这可能是整个循环的替代品，尽管所有这些索引都让我有点牵强。（虽然这会创建一个大的中间数组）

comp = numpy.sum(
        numpy.min((VV1.reshape(-1,1,1,1) - precomp[numpy.newaxis,...,0])**2
                 +(HH1.reshape(-1,1,1,1) - precomp[numpy.newaxis,...,1])**2, 
                 axis=2),
        axis=0)
于 2013-05-17T11:43:23.323 回答
0

在计算机科学中，有大 O 表示法的概念，用于获得做某事所需工作量的近似值。要快速编写程序，请尽可能少做。

这就是为什么 Janne 的答案如此之快，你做的计算更少。更进一步，我们可以应用memoization的概念，因为您受 CPU 限制而不是 RAM 限制。如果需要比以下示例更复杂，您可以使用内存库。

class AutoVivification(dict):
        """Implementation of perl's autovivification feature."""
        def __getitem__(self, item):
            try:
                return dict.__getitem__(self, item)
            except KeyError:
                value = self[item] = type(self)()
                return value

memo = AutoVivification()

def memoize(n, arr, end):
    if not memo[n][arr][end]:
        memo[n][arr][end] = (n-arr[...,end])**2        
    return memo[n][arr][end]


for (vv,hh) in zip(VV1,HH1):
    first = memoize(vv, precomp, 0)
    second = memoize(hh, precomp, 1)
    comp+= numpy.min(first+second, axis=2)

任何已经计算过的东西都会保存到字典中的内存中，我们可以稍后查找它而不是重新计算它。您甚至可以将正在完成的数学分解为更小的步骤，如有必要，每个步骤都被记住。

AutoVivification 字典只是为了更容易将结果保存在 memoize 中，因为我很懒。同样，你可以记住你做的任何数学，所以如果numpy.min很慢，也记住它。

于 2013-05-17T12:30:15.830 回答
0

并行化循环的一种方法是以使用map. 在这种情况下，您可以使用multiprocessing.Pool来使用并行地图。

我会改变这个：

for (vv,hh) in zip(VV1,HH1):
    comp+= numpy.min((vv-precomp[...,0])**2+(hh-precomp[...,1])**2, axis=2)

对于这样的事情：

def buildcomp(vvhh):
    vv, hh = vvhh
    return numpy.min((vv-precomp[...,0])**2+(hh-precomp[...,1])**2, axis=2)

if __name__=='__main__':
    from multiprocessing import Pool
    nthreads = 2
    p = Pool(nthreads)
    complist = p.map(buildcomp, np.column_stack((VV1,HH1)))
    comp = np.dstack(complist).sum(-1)

请注意，dstack假设每个comp.ndim都是2，因为它将添加第三个轴，并沿它求和。这会减慢它的速度，因为您必须构建列表，将其堆叠，然后对其求和，但这些都是并行或 numpy 操作。

我还将其更改zip为 numpy operation np.column_stack，因为对于长数组zip来说要慢得多，假设它们已经是一维数组（它们在您的示例中）。

我不能轻易测试这个，所以如果有问题，请随时让我知道。

于 2013-05-17T15:01:31.813 回答