Numpy/Python 与 Matlab 相比表现非常糟糕

Question

这里是新手程序员。我正在编写一个程序来分析点（单元格）的相对空间位置。程序从数组中获取边界和单元格类型，其中 x 坐标位于第 1 列，y 坐标位于第 2 列，单元格类型位于第 3 列。然后它检查每个单元格的单元格类型和与边界的适当距离。如果它通过，它会计算它与数组中每个其他单元格的距离，如果距离在指定的分析范围内，它会将其添加到该距离处的输出数组中。

我的单元格标记程序在 wxpython 中，所以我希望也可以在 python 中开发这个程序，并最终将它粘贴到 GUI 中。不幸的是，现在 python 在我的机器上运行核心循环需要大约 20 秒，而 MATLAB 可以每秒执行大约 15 个循环。由于我计划对约 30 个案例乘以几种探索性分析类型进行 1000 个循环（具有随机比较条件），这不是一个微不足道的差异。

我尝试运行分析器，数组调用的时间是 1/4，其余几乎都是未指定的循环时间。

这是主循环的python代码：

for basecell in range (0, cellnumber-1):
    if firstcelltype == np.array((cellrecord[basecell,2])):
        xloc=np.array((cellrecord[basecell,0]))
        yloc=np.array((cellrecord[basecell,1]))
        xedgedist=(xbound-xloc)
        yedgedist=(ybound-yloc)
        if xloc>excludedist and xedgedist>excludedist and yloc>excludedist and    yedgedist>excludedist:
            for comparecell in range (0, cellnumber-1):
                if secondcelltype==np.array((cellrecord[comparecell,2])):
                    xcomploc=np.array((cellrecord[comparecell,0]))
                    ycomploc=np.array((cellrecord[comparecell,1]))
                    dist=math.sqrt((xcomploc-xloc)**2+(ycomploc-yloc)**2)
                    dist=round(dist)
                    if dist>=1 and dist<=analysisdist:
                         arraytarget=round(dist*analysisdist/intervalnumber)
                         addone=np.array((spatialraw[arraytarget-1]))
                         addone=addone+1
                         targetcell=arraytarget-1
                         np.put(spatialraw,[targetcell,targetcell],addone)

这是主循环的matlab代码：

for basecell = 1:cellnumber;
    if firstcelltype==cellrecord(basecell,3);
         xloc=cellrecord(basecell,1);
         yloc=cellrecord(basecell,2);
         xedgedist=(xbound-xloc);
         yedgedist=(ybound-yloc);
         if (xloc>excludedist) && (yloc>excludedist) && (xedgedist>excludedist) && (yedgedist>excludedist);
             for comparecell = 1:cellnumber;
                 if secondcelltype==cellrecord(comparecell,3);
                     xcomploc=cellrecord(comparecell,1);
                     ycomploc=cellrecord(comparecell,2);
                     dist=sqrt((xcomploc-xloc)^2+(ycomploc-yloc)^2);
                     if (dist>=1) && (dist<=100.4999);
                         arraytarget=round(dist*analysisdist/intervalnumber);
                         spatialsum(1,arraytarget)=spatialsum(1,arraytarget)+1;
                    end
                end
            end            
        end
    end
end

谢谢！

Answer 1

以下是一些加快 Python 代码速度的方法。

第一：当你只存储一个值时不要创建 np 数组。您在代码中多次执行此操作。例如，

if firstcelltype == np.array((cellrecord[basecell,2])):

只能是

 if firstcelltype == cellrecord[basecell,2]:

我会用一些 timeit 语句告诉你原因：

>>> timeit.Timer('x = 111.1').timeit()
0.045882196294822819
>>> t=timeit.Timer('x = np.array(111.1)','import numpy as np').timeit()
0.55774970267830071

这些调用之间的差异是一个数量级。

第二：以下代码：

arraytarget=round(dist*analysisdist/intervalnumber)
addone=np.array((spatialraw[arraytarget-1]))
addone=addone+1
targetcell=arraytarget-1
np.put(spatialraw,[targetcell,targetcell],addone)

可以替换为

arraytarget=round(dist*analysisdist/intervalnumber)-1
spatialraw[arraytarget] += 1

第三：您可以通过事先平方analysisdist 来摆脱菲利普提到的 sqrt。但是，由于您使用 analysisdist 来获得 arraytarget，因此您可能需要创建一个单独的变量 analysisdist2，它是 analysisdist 的平方并将其用于比较。

第四：每次到达该点时，您都​​在寻找与 secondcelltype 匹配的单元格，而不是一次又一次地找到这些单元格并一遍又一遍地使用该列表。你可以定义一个数组：

comparecells = np.where(cellrecord[:,2]==secondcelltype)[0]

然后替换

for comparecell in range (0, cellnumber-1):
    if secondcelltype==np.array((cellrecord[comparecell,2])):

与

for comparecell in comparecells:

第五：使用psyco。它是一个 JIT 编译器。如果您使用的是较新的版本，Matlab 有一个内置的 JIT 编译器。这应该会加快您的代码速度。

第六：如果在前面的所有步骤之后代码仍然不够快，那么您应该尝试对代码进行矢量化。应该不会太难。基本上，你可以在 numpy 数组中拥有的东西越多越好。这是我在矢量化方面的尝试：

basecells = np.where(cellrecord[:,2]==firstcelltype)[0]
xlocs = cellrecord[basecells, 0]
ylocs = cellrecord[basecells, 1]
xedgedists = xbound - xloc
yedgedists = ybound - yloc
whichcells = np.where((xlocs>excludedist) & (xedgedists>excludedist) & (ylocs>excludedist) & (yedgedists>excludedist))[0]
selectedcells = basecells[whichcells]
comparecells = np.where(cellrecord[:,2]==secondcelltype)[0]
xcomplocs = cellrecords[comparecells,0]
ycomplocs = cellrecords[comparecells,1]
analysisdist2 = analysisdist**2
for basecell in selectedcells:
    dists = np.round((xcomplocs-xlocs[basecell])**2 + (ycomplocs-ylocs[basecell])**2)
    whichcells = np.where((dists >= 1) & (dists <= analysisdist2))[0]
    arraytargets = np.round(dists[whichcells]*analysisdist/intervalnumber) - 1
    for target in arraytargets:
        spatialraw[target] += 1

您可能可以取出内部 for 循环，但您必须小心，因为 arraytargets 的某些元素可能是相同的。另外，我实际上并没有尝试所有的代码，所以那里可能存在错误或错字。希望它能让您很好地了解如何执行此操作。哦，还有一件事。您将 analysisdist/intervalnumber 设为一个单独的变量，以避免一遍又一遍地进行除法。

Answer 2

不太确定 python 的速度，但您的 Matlab 代码可以进行高度优化。嵌套的 for 循环往往存在可怕的性能问题。您可以用矢量化函数替换内部循环...如下：

for basecell = 1:cellnumber;
    if firstcelltype==cellrecord(basecell,3);
         xloc=cellrecord(basecell,1);
         yloc=cellrecord(basecell,2);
         xedgedist=(xbound-xloc);
         yedgedist=(ybound-yloc);
         if (xloc>excludedist) && (yloc>excludedist) && (xedgedist>excludedist) && (yedgedist>excludedist);
%             for comparecell = 1:cellnumber;
%                 if secondcelltype==cellrecord(comparecell,3);
%                     xcomploc=cellrecord(comparecell,1);
%                     ycomploc=cellrecord(comparecell,2);
%                     dist=sqrt((xcomploc-xloc)^2+(ycomploc-yloc)^2);
%                     if (dist>=1) && (dist<=100.4999);
%                         arraytarget=round(dist*analysisdist/intervalnumber);
%                         spatialsum(1,arraytarget)=spatialsum(1,arraytarget)+1;
%                    end
%                end
%            end
         %replace with:
        secondcelltype_mask = secondcelltype == cellrecord(:,3);
        xcomploc_vec = cellrecord(secondcelltype_mask ,1);
                ycomploc_vec = cellrecord(secondcelltype_mask ,2);
                dist_vec = sqrt((xcomploc_vec-xloc)^2+(ycomploc_vec-yloc)^2);
                dist_mask = dist>=1 & dist<=100.4999
                arraytarget_vec = round(dist_vec(dist_mask)*analysisdist/intervalnumber);
                count = accumarray(arraytarget_vec,1, [size(spatialsum,1),1]);
                spatialsum(:,1) = spatialsum(:,1)+count;
        end
    end
end

那里可能有一些小错误，因为我没有任何数据可以用来测试代码，但它应该可以在 Matlab 代码上获得约 10 倍的速度。

根据我使用 numpy 的经验，我注意到将 for 循环换成基于矢量化/矩阵的算术也有明显的加速效果。但是，如果没有形状，所有变量的形状就很难矢量化。

Answer 3

您可以通过替换这些行来避免一些math.sqrt 调用

                dist=math.sqrt((xcomploc-xloc)**2+(ycomploc-yloc)**2)
                dist=round(dist)
                if dist>=1 and dist<=analysisdist:
                     arraytarget=round(dist*analysisdist/intervalnumber)

与

                dist=(xcomploc-xloc)**2+(ycomploc-yloc)**2
                dist=round(dist)
                if dist>=1 and dist<=analysisdist_squared:
                     arraytarget=round(math.sqrt(dist)*analysisdist/intervalnumber)

你的线路在哪里

 analysisdist_squared = analysis_dist * analysis_dist

在函数的主循环之外。

由于math.sqrt在最里面的循环中被调用，你应该在模块的顶部有from math import sqrt，然后将函数调用为sqrt。

我也会尝试替换

                dist=(xcomploc-xloc)**2+(ycomploc-yloc)**2

与

                dist=(xcomploc-xloc)*(xcomploc-xloc)+(ycomploc-yloc)*(ycomploc-yloc)

它有可能产生更快的字节码来进行乘法而不是求幂。

我怀疑这些能否让您获得 MATLAB 的性能，但它们应该有助于减少一些开销。

Answer 4

如果您有一个多核，您可以尝试使用多处理模块并使用多个进程来利用所有核。

您可以使用 x**0.5 代替 sqrt，如果我没记错的话，它会稍微快一些。