Modo Python efficiente per equazioni ricorsive

Question

Sto cercando di ottimizzare un loop che ho in un pezzo del mio codice. Ho pensato che scriverlo in un modo più intorpidito avrebbe reso più veloce, ma ora è più lento! le equazioni prende in input una vec numpy.array di lunghezza n:

from numpy import * 

def f(vec): 
    n=len(vec) 
    aux=0 
    for i in range(n): 
     aux = aux + (1- aux)*vec[i] 
    return aux 

def f2(vec): 
    n=len(vec) 
    G=tril(array([-vec]*n),-1)+1    #numpy way! 
    aux=dot(G.prod(1),vec) 
    return aux 


if __name__ == '__main__': 
    import timeit 
    print(timeit.timeit("f(ones(225)+4)", setup="from __main__ import f\nfrom numpy import ones",number=1000)) 
    print(timeit.timeit("f2(ones(225)+4)", setup="from __main__ import f2\nfrom numpy import ones,tril,dot",number=1000)) 

0,429496049881 [s]

5,66514706612 [s]

ho finalmente deciso di inserte l'intera funzione nel mio ciclo, ottenere un incremento delle prestazioni 3 volte. Sto davvero cercando un aumento delle prestazioni di 100 volte, ma non so cos'altro fare. Questa è la mia ultima funzione:

def CALC_PROB_LOC2(int nSectors, int nZones,double[:] beta, double[:] thetaLoc,np.ndarray[double, ndim=2] h, np.ndarray[double, ndim=2] p, np.ndarray[np.float64_t, ndim=3] U_nij, np.ndarray[double, ndim=2] A_ni): 
    cdef np.ndarray[double, ndim=3] Pr_nij =np.zeros((nSectors,nZones,nZones),dtype="d") 
    cdef np.ndarray[double, ndim=2] U_ni =np.zeros((nSectors,nZones),dtype="d") 
    #cdef np.ndarray[np.float64_t, ndim=1] A_ni_pos 
    cdef Py_ssize_t n,i,opt 
    cdef int aux_bool,options 
    cdef np.ndarray[np.float64_t, ndim=1] prob,attractor,optionCosts 
    cdef np.ndarray[np.float64_t, ndim=1] eq23,utilities 
    cdef double disu 
    cdef double eq22 
    cdef double aux17 
    for n in range(nSectors): 
     aux_bool=1 
     if n in [0,2,9,10,11,12,13,14,18,19,20]: 
      for i in xrange(nZones): 
       U_ni[n,i]=p[n,i]+h[n,i] 
       Pr_nij[n,i,i]=1 
      aux_bool=0 
     if aux_bool==1: 
      if beta[n]<=0: 
       for i in xrange(nZones): 
        U_ni[n,i]=U_nij[n,i,i] 
      else: 
       A_ni_pos=A_ni[n,:]>0 
       options=len(A_ni[n,:][A_ni_pos]) 
       attractor=A_ni[n,:][A_ni_pos] 
       if options>0: 
        for i in xrange(nZones): 
         optionCosts=U_nij[n,i,A_ni_pos] 
         disu=0 
         eq22=0 
         aux17=0 
         prob=np.ones(options)/options #default value 
         if beta[n]==0: 
          Pr_nij[n,i,A_ni_pos],U_ni[n,i]= prob,0 
         if options==1: 
          Pr_nij[n,i,A_ni_pos],U_ni[n,i]= prob,optionCosts 
         else: 
          if thetaLoc[n]<=0: 
           cmin=1 
          else: 
           cmin=(optionCosts**thetaLoc[n]).min() 
           if cmin==0: 
            cmin=100 
          utilities=optionCosts/cmin 
          eq23=-beta[n]*utilities 
          eq23=np.exp(eq23) 
          aux17=np.dot(attractor,eq23) 
          if aux17==0: 
           Pr_nij[n,i,A_ni_pos],U_ni[n,i]= 0*prob,0 
          else: 
           for opt in range(options): 
            eq22=eq22+(1-eq22)*eq23[opt] 
           prob=attractor*eq23/aux17 
           disu=cmin*(-np.log(eq22)/beta[n]) 
           Pr_nij[n,i,A_ni_pos],U_ni[n,i]= prob,disu 


    return Pr_nij,U_ni 

Answer 1

è quello che succede quando un algoritmi lineari viene sostituito da un quadratica uno: Non importa quanto velocemente è eseguito, l'algoritmo migliore vince sempre (per un problema abbastanza grande).

È abbastanza chiaro che f viene eseguito in tempo lineare e f2 viene eseguito in tempo quadratorio perché è la complessità temporale di un prodotto con punti a matrice.

Un grafico log-log mostra chiaramente la differenza di tempo di esecuzione (lineare si riferisce a f, quadractic a f2):

più a destra parte della linea verde (cioè, quando non si comporta come una linea retta) può essere spiegato perché le funzioni di numpy di solito hanno un overhead elevato che è trascurabile per gli array che non sono minuscoli ma domina il tempo di esecuzione quando sono piccoli.

---

Il modo "standard" per accelerare il codice in Python che è già utilizzando un algoritmo veloce è quello di raggiungere per il codice compilato e scrivere una proroga. Cython consente di farlo annotando il codice sorgente di Python con alcune annotazioni di tipo e comprende gli array numpy.

Raccontando Cython che vec è una serie di doppie, aux è un doppio e i un intero, è in grado di generare un'estensione C che è 400x più veloce per me.

def f(double[:] vec): 
    n = len(vec) 
    cdef double aux = 0 
    cdef int i 
    for i in range(n): 
     aux = aux + (1- aux)*vec[i] 
    return aux 

Se vi capita di essere utilizzando IPython, si può semplicemente eseguire %load_ext cythonmagic e quindi copiare quella funzione a una cella prefissato dalla linea %%cython di provarlo. Altri metodi per compilarlo e compilarlo sono spiegati in Cython documentation. A proposito, IPython ti permette anche di scrivere il codice scrivendo %timeit prima dell'affermazione, è davvero utile.

Un'opzione completamente diversa è l'uso di PyPy, un'implementazione di Python 2.7 che viene fornito con un JIT e ha un supporto numerico di base. Può eseguire questo piccolo snippet sostituendo import numpypy per import numpy, ma è possibile che non sia in grado di eseguire l'intero programma. È un po 'più lento di Cython ma non richiede un compilatore né annota il codice.