Numpy/Python prestazioni terribilmente vs. Matlab

Question

Programmatore novizio qui. Sto scrivendo un programma che analizza le posizioni spaziali relative dei punti (cellule). Il programma ottiene i confini e il tipo di cella da una matrice con la coordinata x nella colonna 1, la coordinata y nella colonna 2 e il tipo di cella nella colonna 3. Quindi controlla ogni cella per il tipo di cella e la distanza appropriata dai limiti. Se passa, calcola la distanza tra le celle della matrice e se la distanza si trova all'interno di un intervallo di analisi specificato, la aggiunge a una matrice di output a quella distanza.

Il mio programma di marcatura celle è in wxpython, quindi speravo di sviluppare questo programma anche in python ed eventualmente inserirlo nella GUI. Sfortunatamente in questo momento python impiega circa 20 secondi per eseguire il core loop sulla mia macchina mentre MATLAB può fare ~ 15 loop/secondo. Dato che sto pianificando di fare 1000 loop (con una condizione di confronto randomizzato) su ~ 30 casi volte diversi tipi di analisi esplorative questa non è una differenza banale.

Ho provato a eseguire un profiler e le chiamate di array sono 1/4 del tempo, quasi tutto il resto è un loop non specificato.

ecco il codice Python per il ciclo 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) 

Ecco il codice MATLAB per il ciclo 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 

Grazie!

Answer 1

Ecco alcuni modi per velocizzare il codice Python.

Primo: Non creare array np quando si memorizza solo un valore. Lo fai più volte nel tuo codice. Per esempio,

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

può essere solo

if firstcelltype == cellrecord[basecell,2]: 

te lo mostrerò perché con alcune dichiarazioni timeit:

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

Questo è un ordine di grandezza nella differenza tra quelle chiamate.

Secondo: il seguente codice:

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

possono essere sostituiti con

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

Terzo: si può sbarazzarsi del sqrt come Philip menzionato da quadratura analysisdist anticipo. Tuttavia, poiché si utilizza analysisdist per ottenere arraytarget, è possibile creare una variabile separata, analysisdist2 che corrisponde al quadrato di analysisdist e utilizzarla per il confronto.

Quarto: Siete alla ricerca di cellule che corrispondono secondcelltype ogni volta che si arriva a quel punto piuttosto che trovare quelli di una volta e utilizzando la lista più e più volte. Si potrebbe definire un array:

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

e quindi sostituire

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

con

for comparecell in comparecells: 

Quinto: Usa psyco. È un compilatore JIT. Matlab ha un compilatore JIT incorporato se stai usando una versione piuttosto recente. Questo dovrebbe velocizzare il tuo codice un po '.

Sesto: Se il codice non è ancora abbastanza veloce dopo tutti i passaggi precedenti, è necessario provare a vettorializzare il codice. Non dovrebbe essere troppo difficile. Fondamentalmente, più cose si possono avere in array di numpy, meglio è. Ecco il mio tentativo per vettorializzazione:

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 

probabilmente si può prendere che interna per ciclo, ma bisogna stare attenti perché alcuni degli elementi di arraytargets potrebbe essere lo stesso. Inoltre, in realtà non ho provato tutto il codice, quindi potrebbe esserci un errore o un errore di battitura. Spero che ti dia una buona idea di come farlo. Oh, un'altra cosa. Si rende analysisdist/intervalnumber una variabile separata per evitare di ripetere questa divisione più e più volte.

Answer 2

È possibile evitare alcune delle math.sqrt chiamate sostituendo le linee

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

con

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

dove si ha la linea

analysisdist_squared = analysis_dist * analysis_dist 

al di fuori del ciclo principale del vostro funzione.

Poiché math.sqrt viene chiamato nel ciclo più interno, è necessario avere from math import sqrt nella parte superiore del modulo e chiamare la funzione come sqrt.

Vorrei anche provare a sostituire

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

con

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

C'è una possibilità che produrrà bytecode più veloce per fare la moltiplicazione invece di elevamento a potenza.

Dubito che questi ti consentiranno di ottenere prestazioni MATLAB, ma dovrebbero contribuire a ridurre alcuni costi generali.

Answer 3

Non troppo sicuro della lentezza di Python ma il codice Matlab può essere ALTAMENTE ottimizzato. I cicli for nidificati tendono ad avere problemi di prestazioni orribili. Puoi sostituire il ciclo interno con una funzione vettoriale ...come di seguito:

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 

Ci possono essere alcuni piccoli errori in là dal momento che non ho i dati per testare il codice con ma dovrebbe ottenere ~ velocità di 10 volte sul codice Matlab.

Dalla mia esperienza con Numpy ho notato che lo scambio di cicli for per aritmetica vettoriale/matrice ha anche notevoli accelerazioni. Tuttavia, senza le forme le forme di tutte le variabili sono difficili da vettorializzare.

Answer 4

Se si dispone di un multicore, è possibile provare il modulo di multiprocessing e utilizzare più processi per utilizzare tutti i core.

Invece di sqrt è possibile utilizzare x ** 0.5, che è, se ricordo corretto, leggermente più veloce.