Numpy/Python performant terriblement vs Matlab

Question

Programmeur novice ici. J'écris un programme qui analyse les emplacements spatiaux relatifs des points (cellules). Le programme obtient les limites et le type de cellule d'un tableau avec la coordonnée x dans la colonne 1, la coordonnée y dans la colonne 2 et le type de cellule dans la colonne 3. Il vérifie ensuite chaque cellule pour le type de cellule et la distance appropriée des limites. S'il passe, il calcule ensuite sa distance à partir de chaque cellule du tableau et si la distance est dans une plage d'analyse spécifiée, il l'ajoute à un tableau de sortie à cette distance. Mon programme de marquage de cellules est en wxpython. J'espérais donc développer ce programme en python et le coller dans l'interface graphique. Malheureusement en ce moment python prend ~ 20 secondes pour exécuter la boucle de base sur ma machine tandis que MATLAB peut faire ~ 15 boucles/seconde. Puisque je prévois de faire 1000 boucles (avec une condition de comparaison aléatoire) sur ~ 30 cas fois plusieurs types d'analyse exploratoire, ce n'est pas une différence triviale.

J'ai essayé d'exécuter un profileur et les appels de groupe sont 1/4 du temps, presque tout le reste est temps de boucle non spécifié.

Voici le code python pour la boucle principale:

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) 

Voici le code Matlab pour la boucle principale:

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 

Merci!

Answer 1

Voici quelques façons d'accélérer votre code python.

Première: Ne faites pas de tableaux np lorsque vous ne stockez qu'une seule valeur. Vous faites cela plusieurs fois dans votre code. Par exemple,

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

peut juste être

if firstcelltype == cellrecord[basecell,2]: 

Je vais vous montrer pourquoi avec quelques déclarations TimeIt:

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

C'est un ordre de grandeur de la différence entre ces appels.

Deuxième: Le code suivant:

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

peut être remplacé par

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

Troisième: Vous pouvez vous débarrasser de l'sqrt comme Philippe mentionné par élévation au carré analysisdist au préalable. Cependant, puisque vous utilisez analysisdist pour obtenir arraytarget, vous pouvez créer une variable distincte, analysisdist2 qui est le carré de analysisdist et l'utiliser pour votre comparaison.

Quatrième: Vous êtes à la recherche de cellules qui correspondent à secondcelltype chaque fois que vous arrivez à ce point plutôt que de trouver les une seule fois et en utilisant la liste encore et encore. Vous pouvez définir un tableau:

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

puis remplacez

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

avec

for comparecell in comparecells: 

Cinquième: Utilisez psyco. C'est un compilateur JIT. Matlab a un compilateur JIT intégré si vous utilisez une version assez récente. Cela devrait accélérer votre code.

Sixième: Si le code n'est toujours pas assez rapide après toutes les étapes précédentes, vous devriez essayer de vectoriser votre code. Cela ne devrait pas être trop difficile. Fondamentalement, le plus de choses que vous pouvez avoir dans des tableaux numpy le mieux. Voici mon essai à vectorisation:

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 

Vous pouvez probablement prendre que pour la boucle intérieure, mais vous devez faire attention parce que certains des éléments de arraytargets pourrait être le même. De plus, je n'ai pas vraiment essayé tout le code, donc il pourrait y avoir un bug ou une faute de frappe. J'espère que cela vous donne une bonne idée de la façon de le faire. Oh, encore une chose. Vous faites analysisdist/intervalnumber une variable distincte pour éviter de faire cette division encore et encore.

Answer 2

Vous pouvez éviter certaines des math.sqrt appels en remplaçant les lignes

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

avec

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

où vous avez la ligne

analysisdist_squared = analysis_dist * analysis_dist 

en dehors de la boucle principale de votre fonction.

Puisque math.sqrt est appelée dans la boucle la plus interne, vous devez avoir from math import sqrt en haut du module et appeler simplement la fonction en tant que sqrt.

Je voudrais aussi essayer de remplacer

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

avec

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

Il y a une chance qu'il produira un code plus rapide d'octets à faire la multiplication plutôt que exponentiation.

Je doute que ceux-ci vous mèneront aux performances de MATLAB, mais ils devraient aider à réduire certains frais généraux.

Answer 3

Pas trop sûr de la lenteur de python mais le code Matlab peut être FORTEMENT optimisé. Les boucles for-nested ont tendance à avoir d'horribles problèmes de performance. Vous pouvez remplacer la boucle interne par une fonction vectorisée ...comme ci-dessous:

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 

Il peut y avoir quelques petites erreurs là-dedans depuis que je n'ai pas de données pour tester le code avec mais il devrait obtenir ~ vitesse 10X sur le code Matlab. D'après mon expérience avec numpy, j'ai remarqué que l'échange de boucles pour l'arithmétique vectorisée/matricielle a également des accélérations visibles. Cependant, sans les formes, les formes de toutes vos variables sont difficiles à vectoriser.

Answer 4

Si vous avez un multicœur, vous pouvez essayer le module multitraitement et utiliser plusieurs processus pour utiliser tous les cœurs. Au lieu de sqrt, vous pouvez utiliser x ** 0.5, ce qui est, si je me souviens bien, légèrement plus rapide.