2011-10-12 90 views
4

我試圖通過計算身體各側的影響來計算埋藏物體的重力影響,然後總結一個站點的重複測量貢獻,以重複測量多個站點。代碼如下(身體是一個廣場和代碼計算順時針方向圍繞它,這就是爲什麼它的-x追溯到-x座標)我如何加速python嵌套循環?

grav = [] 
x=si.arange(-30.0,30.0,0.5) 

#-9.79742526  9.78716693 22.32153704 27.07382349 2138.27146193 
xcorn = (-9.79742526,9.78716693 ,9.78716693 ,-9.79742526,-9.79742526) 
zcorn = (22.32153704,22.32153704,27.07382349,27.07382349,22.32153704) 
gamma = (6.672*(10**-11))#'N m^2/Kg^2' 
rho = 2138.27146193#'Kg/m^3' 
grav = [] 
iter_time=[] 
def procedure(): 
    for i in si.arange(len(x)):# cycles position 
     t0=time.clock() 
     sum_lines = 0.0 

     for n in si.arange(len(xcorn)-1):#cycles corners 
      x1 = xcorn[n]-x[i] 
      x2 = xcorn[n+1]-x[i] 
      z1 = zcorn[n]-0.0 #just depth to corner since all observations are on the surface. 
      z2 = zcorn[n+1]-0.0 
      r1 = ((z1**2) + (x1**2))**0.5 
      r2 = ((z2**2) + (x2**2))**0.5 
      O1 = si.arctan2(z1,x1) 
      O2 = si.arctan2(z2,x2) 
      denom = z2-z1 
      if denom == 0.0: 
       denom = 1.0e-6 

      alpha = (x2-x1)/denom 

      beta = ((x1*z2)-(x2*z1))/denom 
      factor = (beta/(1.0+(alpha**2))) 
      term1 = si.log(r2/r1)#log base 10 
      term2 = alpha*(O2-O1) 
      sum_lines = sum_lines + (factor*(term1-term2)) 
     sum_lines = sum_lines*2*gamma*rho 
     grav.append(sum_lines) 
     t1 = time.clock() 
     dt = t1-t0 
     iter_time.append(dt) 

在加速這個循環中的任何幫助,將不勝感激謝謝。

+4

暗示http://codereview.stackexchange.com/對於這種類型的問題,更好的論壇。 – mjv

+0

@aix:'si.arange()'可能暗示他已經使用'numpy'數組。 – jfs

+0

@ J.F.Sebastian:很好。我刪除了我的評論,以免引起任何進一步的混淆。 – NPE

回答

1

重複您的xcorn和zcorn值,因此請考慮緩存某些計算的結果。

查看timeit和profile模塊以獲取有關計算時間最多的信息。

1

在Python循環中訪問numpy數組的單個元素效率非常低。例如,這條巨蟒循環:

for i in xrange(0, len(a), 2): 
    a[i] = i 

將遠慢於:

a[::2] = np.arange(0, len(a), 2) 

你可以使用一個更好的算法(更少的時間複雜度)或numpy陣列使用矢量操作,如上面的例子。但是更快的方法可能只是使用Cython編譯代碼:

#cython: boundscheck=False, wraparound=False 
#procedure_module.pyx 
import numpy as np 
cimport numpy as np 

ctypedef np.float64_t dtype_t 

def procedure(np.ndarray[dtype_t,ndim=1] x, 
       np.ndarray[dtype_t,ndim=1] xcorn): 
    cdef: 
     Py_ssize_t i, j 
     dtype_t x1, x2, z1, z2, r1, r2, O1, O2 
     np.ndarray[dtype_t,ndim=1] grav = np.empty_like(x) 
    for i in range(x.shape[0]): 
     for j in range(xcorn.shape[0]-1): 
      x1 = xcorn[j]-x[i] 
      x2 = xcorn[j+1]-x[i] 
      ... 
     grav[i] = ... 
    return grav 

沒有必要定義所有的類型,但如果你需要一個顯著加快相比,Python中,你至少應該類型的數組和循環的定義索引。

您可以使用cProfile(Cython支持它)而不是手動撥打time.clock()

要調用procedure()

#!/usr/bin/env python 
import pyximport; pyximport.install() # pip install cython 
import numpy as np 

from procedure_module import procedure 

x = np.arange(-30.0,30.0,0.5) 
xcorn = np.array((-9.79742526,9.78716693 ,9.78716693 ,-9.79742526,-9.79742526)) 
grav = procedure(x, xcorn) 
+0

感謝您的幫助,我看看cython,看看是否有幫助,我不擔心time.clock(),因爲我只有在那裏看到每次迭代的速度有多快,而我試着加快速度,它不是代碼的必要部分。 – user991926