Numpy/Python 表现糟糕 vs. Matlab
我是个新手程序员。我正在写一个程序，用来分析一些点（细胞）在空间中的相对位置。这个程序从一个数组中获取边界和细胞类型，数组的第一列是x坐标，第二列是y坐标，第三列是细胞类型。然后，它会检查每个细胞的类型和它与边界的距离。如果符合条件，程序就会计算它与数组中其他细胞的距离，如果这个距离在指定的分析范围内，就把它添加到一个输出数组中，记录这个距离。
我的细胞标记程序是用wxpython写的，所以我希望这个程序也能用python开发，最后把它放到图形界面里。不过现在的问题是，我的机器上python运行这个核心循环大约需要20秒，而MATLAB能做到每秒大约15次循环。因为我计划要进行1000次循环（带有随机比较条件），在大约30个案例和几种探索性分析类型下，这个时间差可不是小事。
我尝试运行了一个性能分析工具，发现数组调用占用了四分之一的时间，剩下几乎全部是未指定的循环时间。
以下是主要循环的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
谢谢！
4 个回答
analysisdist_squared = analysis_dist * analysis_dist
你可以通过把以下几行
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)
来减少一些math.sqrt的调用。
因为math.sqrt是在最里面的循环中被调用的，所以你应该在模块的最上面加上from math import sqrt，然后直接用sqrt来调用这个函数。
我还会尝试把
dist=(xcomploc-xloc)**2+(ycomploc-yloc)**2
替换成
dist=(xcomploc-xloc)*(xcomploc-xloc)+(ycomploc-yloc)*(ycomploc-yloc)
这样做可能会让乘法的字节码运行得更快，而不是用指数运算。
我不确定这些改动能否让你的代码达到MATLAB的性能，但它们应该能帮助减少一些额外的开销。
我不太确定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循环，也能明显加快速度。不过，如果不知道你所有变量的形状，就很难进行向量化处理。
这里有一些加快你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
来替换。
第三： 正如Philip提到的，你可以通过提前将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中的某些元素可能是相同的。此外，我并没有实际运行所有的代码，所以里面可能有bug或者拼写错误。希望这能给你一个好的思路。哦，还有一件事。你可以把analysisdist/intervalnumber单独做成一个变量，以避免重复进行这个除法。