Numpy/Python con un rendimiento terrible en comparación con Matlab

Question

Programador novato aquí. Estoy escribiendo un programa que analiza las ubicaciones espaciales relativas de puntos (células). El programa obtiene los límites y el tipo de celda de una matriz con la coordenada x en la columna 1, la coordenada y en la columna 2 y el tipo de celda en la columna 3. A continuación, verifica el tipo de celda y la distancia apropiada desde los límites. Si pasa, entonces calcula su distancia de cada otra celda en la matriz y si la distancia está dentro de un rango de análisis especificado, lo agrega a una matriz de salida a esa distancia.

Mi programa de marcado de células está en wxpython, así que esperaba desarrollar este programa también en python y eventualmente incluirlo en la GUI. Lamentablemente, en este momento, Python tarda ~ 20 segundos en ejecutar el bucle central en mi máquina, mientras que MATLAB puede hacer ~ 15 bucles/segundo. Como estoy planeando hacer 1000 bucles (con una condición de comparación aleatorizada) en ~ 30 casos por varios tipos de análisis exploratorios, esta no es una diferencia trivial.

He intentado ejecutar un generador de perfiles y las llamadas de conjunto son 1/4 de las veces, casi todo el resto es tiempo de ciclo no especificado.

Este es el código Python para el bucle principal:

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) 

Aquí está el código de MATLAB para el bucle principal:

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 

Gracias!

Answer 1

Aquí hay algunas maneras de acelerar su código python.

Primero: No haga matrices np cuando solo está almacenando un valor. Usted hace esto muchas veces en su código. Por ejemplo,

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

puede ser sólo

if firstcelltype == cellrecord[basecell,2]: 

te muestro lo que con algunas declaraciones TimeIt:

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

Eso es un orden de magnitud de diferencia entre esas llamadas.

Segundo: El siguiente código:

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

pueden ser sustituidos por

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

Tercero: usted puede deshacerse de la raíz cuadrada como Philip mencionado elevando al cuadrado analysisdist de antemano. Sin embargo, dado que usa analysisdist para obtener arraytarget, es posible que desee crear una variable separada, analysisdist2, que sea el cuadrado de la barra de análisis y usarla para su comparación.

Cuarto: Usted está buscando células que responden a secondcelltype cada vez que se llega a ese punto en lugar de encontrar los que una vez y utilizando la lista una y otra vez. Se podría definir una matriz:

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

y luego vuelva a colocar

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

con

for comparecell in comparecells: 

Quinto: Uso psyco. Es un compilador JIT. Matlab tiene un compilador JIT incorporado si está utilizando una versión algo reciente. Esto debería acelerar tu código un poco.

Sexto: Si el código aún no es lo suficientemente rápido después de todos los pasos anteriores, intente vectorizar su código. No debería ser demasiado difícil. Básicamente, cuantas más cosas puedas tener en matrices numpy, mejor. Aquí está mi intento en la vectorización:

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 

probablemente le puede sacar ese interior de bucle, pero hay que tener cuidado porque algunos de los elementos de arraytargets podría ser el mismo. Además, en realidad no probé todo el código, por lo que podría haber un error o un error tipográfico allí. Con suerte, te da una buena idea de cómo hacer esto. Oh, una cosa más. Usted hace analysisdist/intervalnumber una variable separada para evitar hacer esa división una y otra vez.

Answer 2

Usted puede evitar algunos de los math.sqrt llamadas mediante la sustitución de las líneas

   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) 

donde se tiene la línea

analysisdist_squared = analysis_dist * analysis_dist 

fuera del bucle principal de su función.

Dado que math.sqrt se llama en el ciclo más interno, debe tener from math import sqrt en la parte superior del módulo y simplemente llamar a la función como sqrt.

También me intente reemplazar

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

con

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

Hay una posibilidad de que se produzca el código de bytes más rápida de hacer la multiplicación en lugar de exponenciación.

Dudo que estos te lleven al rendimiento de MATLAB, pero deberían ayudar a reducir algunos gastos generales.

Answer 3

No estoy seguro de la lentitud de Python pero su código Matlab puede ser ALTAMENTE optimizado. Los bucles foráneos anidados tienden a tener problemas de rendimiento horribles. Puedes reemplazar el lazo interno con una función vectorizada ...de la siguiente manera:

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 

Puede haber algunos pequeños errores en allí ya que no tengo datos para probar el código con pero debe obtener la velocidad de ~ 10 veces más arriba en el código Matlab.

Según mi experiencia con Numpy, he notado que el intercambio de bucles for vectorizados/aritméticos basados ​​en matrices también tiene una notable aceleración. Sin embargo, sin las formas de todas las variables, es difícil vectorizar las cosas.

Answer 4

Si tiene una multinúcleo, quizás podría probar el módulo de multiprocesamiento y usar procesos múltiples para hacer uso de todos los núcleos.

En lugar de sqrt puede usar x ** 0.5, que es, si mal no recuerdo, un poco más rápido.