加快對numpy數組的分析

Question

我有一個python代碼，它導入4列txt文件，其中第 前三列是x，y，z座標，第四列是該座標下的密度。

下面是讀取代碼，轉換爲ndarray，傅立葉變換該字段，計算原點距離（k =（0,0,0））和變換後的座標，並取平均值並繪製它們。 感謝大熊貓（用於數據分析的python庫）和python FFT，加載256^3行並對它們進行傅里葉變換非常快，並在幾秒鐘內完成。

但是，將加載的txt轉換爲numpy ndarray，計算平均密度（每個座標的平均值）以及計算與原點的距離（k =（0,0,0））需要很長時間。

我認爲問題是在最後np.around部分，但我不知道如何優化它。

我有一個32核心機器的資源。

有人可以教我如何加速，使其成爲一個多進程代碼或類似的東西，以便這些可以很快完成？謝謝。

（如果你是宇宙學家和以往任何時候都需要這個代碼，你可以，如果你可以的。謝謝使用它，但請與我聯繫，）

from __future__ import division 
import numpy as np 

ngridx = 128 
ngridy = 128  
ngridz = 128 

maxK = max(ngridx,ngridy,ngridz) 

#making input file 
f = np.zeros((ngridx*ngridy*ngridz,4)) 

i = 0 
for i in np.arange(len(f)): 
    f[i][0] = int(i/(ngridy*ngridz)) 
    f[i][1] = int((i/ngridz))%ngridy 
    f[i][2] = int(i%ngridz) 
    f[i][3] = np.random.rand(1) 
    if i%1000000 ==0: 
     print i 
#This takes forever 
#end making input file 

#Thanks to Mike, 
a = f[:,3].reshape(ngridx,ngridy,ngridz) 

avg =np.sum(f[:,3])/len(f) 
a /= avg 
p = np.fft.fftn(a) 
#This part is much much faster than before (Original Post). 

#Keeping track of corresponding wavenumbers (k_x, k_y,k_z) for each element in p 
#This is just a convension on fourier transformation so you can ignore this part 
kValx = np.fft.fftfreq(ngridx , (1.0/ngridx)) 
kValy = np.fft.fftfreq(ngridy , (1.0/ngridy)) 
kValz = np.fft.fftfreq(ngridz , (1.0/ngridz)) 
kx = np.zeros((ngridx,ngridy,ngridz)) 
ky = np.zeros((ngridx,ngridy,ngridz)) 
kz = np.zeros((ngridx,ngridy,ngridz)) 
rangecolx = np.arange(ngridx) 
rangecoly = np.arange(ngridy) 
rangecolz = np.arange(ngridz) 
for row in np.arange(ngridx): 
    for column in np.arange(ngridy): 
     for height in np.arange(ngridz): 
      kx[row][column][height] = (kValx[row]) 
      ky[row][column][height] = (kValy[column]) 
      kz[row][column][height] = (kValz[height]) 
    if row%10 == 0: 
     print row 
print 'wavenumber generate complete!' 

#Calculating the average powerspectrum in terms of fixed K (Distance from origin to a point in fourier space) 
#by taking the spherical shell of thickness 1 and averaging out the values inside it. 
#I am sure that this process can be optimised somehow, but I gave up. 

qlen = maxK/2 #Nyquist frequency 
q = np.zeros(((qlen),4),dtype=complex) 
#q is a four column array with length maxK/2. 
#q[:,0] is integer wavenumber (K, which is the distance from the origin = sqrt(kx^2+ky^2+kz^2)) 
#q[:,1] is the sum of square of the fourier transformed value 
#q[:,2] is the sum of the fourier transformed value, 
#and q[:,3] is the total number of samples with K=q[:,0] 

for i in np.arange(len(q)): 
    q[i][0] = i 
i = 0 
for i in np.arange(len(p)): 
    for r in np.arange(len(p[0])): 
     for s in np.arange(len(p[0,0])): 
      K = np.around(np.sqrt(kx[i,r,s]**2+ky[i,r,s]**2+kz[i,r,s]**2)) 
      if K < qlen: 
       q[K][1]=q[K][1]+np.abs(p[i,r,s])**2 
       q[K][2]=q[K][2]+p[i,r,s] 
       q[K][3]=q[K][3]+1 
    if i%10 ==0: 
     print 'i = ',i,' !' 
print q 

Answer 1

numpy的通常可以做的事情數百時間比普通的Python快，只需很少的額外努力。你只需要知道編寫代碼的正確方法。僅舉的第一件事情我想：

索引

普通的Python是經常的事情，電腦應該是很大的很慢。一個例子是索引，所以像

a[f[i,0]][f[i,1]][f[i,2]]=f[i,3] 

讓我非常懷疑。當你說「將加載的txt轉換爲numpy ndarray」需要很長時間時，這就是你指的那個嗎？這並不會讓我感到意外，因爲每次python看到a[f[i,0]]時，它必須首先索引f，確保i是一個整數，並且您還沒有跑出f的邊緣;那麼它必須確保f[i,0]是一個整數，並且您還沒有跑出a的邊緣。然後在它知道你要設置哪個元素之前，必須重複這兩次。

一個改進是使用a[f[i,0],f[i,1],f[i,2]]，因爲在這種索引下numpy速度更快。

但我假設你的數據實際上是某種順序。例如，f[i,2]是否從0循環到256，然後f[i,1]增加1，並且f [i，2]從0開始？如果是這樣，你真正需要做的是這樣說

a = f[:,3].reshape(ngridx,ngridy,ngridz) 

這是一個可笑的快速操作，以一毫秒的一小部分。形狀可能是錯誤的，所以你可能不得不改變參數的順序，對轉置做些事情，但基本的想法肯定是存在的。您可以在the documentation中閱讀。

複製數據是很糟糕

你並不需要複製的一切，當你需要複製（數組或部分）的數組，你應該讓numpy的爲你做它。例如，您可以使用a[1:]而不是Firstdel函數。或者，如果你真的需要使數據的副本（你只是爲了繪製不）使用：

def Firstdel(a): 
    return numpy.copy(a[1:]) 

但在一般情況下，你可以只使用numpy的陣列的「切片」，而不是複製它們。閱讀關於它here。

循環

循環也臭名昭着的時間浪費。首先，while在Python的簡單循環中並不常見。因此，而不是

while i < len(f): 
    # do stuff 
    i = i+1 

你應該使用

for i in range(len(f)): 
    # do stuff 

擺脫儘可能多的循環，你可以。要設置kx，ky，並kz，該代碼比嵌套循環快約10倍，但秤如N代替N-^3（其中N = ngridx ngridy ngridz）：

for row in range(ngridx): 
    kx[row,:,:] = kValx[row] 
for column in range(ngridy): 
    ky[:,column,:] = kValy[column] 
for height in range(ngridz): 
    kz[:,:,height] = kValz[height] 

切片可以對於設置值也很有用，因爲循環進入了numpy。取而代之的是代碼

i = 0 
while i < len(q): 
    q[i][0] = i 
    i = i + 1 

只使用

q[:,0] = range(len(q)) 

在這裏，我只是設置的q「切片」等於另一個數組。

該循環之後的嵌套循環也可以加速，但它們會更復雜一些。

但是你也希望儘可能避免循環。這給我們帶來...

使用內置numpy的功能

numpy的存在的原因是把這些蟒蛇緩慢環路和成等快速C代碼（我們並不需要知道其存在）。所以有很多功能可以完成你想要做的事情，而這些功能已經內置在numpy中。

而不是

while i < len(f): 
    masstot = masstot + f[i,3] 
    i = i+1 

你應該使用類似

masstot = np.sum(f[:,3]) 

所以，很簡單閱讀，但也很可能是方式更快，因爲numpy的必須在該數據的直接訪問計算機的內存，並可以使用快速的C函數來查找總和，而不是使用慢python函數。 （再一次，你不需要知道任何關於C函數的東西，他們只會完成這項工作。）

代替K表明大嵌套循環計算值通過循環每一次，只是使K用適當的值的數組：

K = np.around(np.sqrt(kx**2+ky**2+kz**2)) 

K將是大小相同kx等。然後，您可以使用advanced indexing來設置q的值。這是我會怎麼做，最後一節：

# Again, we get rid of nested loops, to get a large improvement in speed and scaling 
K = np.around(np.sqrt(kx**2+ky**2+kz**2)).astype(int) 
for i in range(qlen): 
    indices = (K==i) # This will be an array of True/False values, 
        # which will be used for "advanced indexing" of p 
    q[i,0] = i 
    q[i,1] = sum(np.abs(p[indices])**2) 
    q[i,2] = sum(p[indices]) 
    q[i,3] = sum(indices) 
print q 

把這些都在一起，我得到的35提高的因素在代碼目前在你的問題。

Answer 2

還不如加快輸入文件創建以及：

size = ngridx*ngridy*ngridz 
f = np.zeros((size,4)) 
a = np.arange(size) 
f[:, 0] = np.floor_divide(a, ngridy*ngridz) 
f[:, 1] = np.fmod(np.floor_divide(a, ngridz), ngridy) 
f[:, 2] = np.fmod(a, ngridz) 
f[:, 3] = np.random.rand(size) 

爲了kx，ky，並kz，你可以擺脫循環使用broadcasting：

kx += kValx[:, np.newaxis, np.newaxis] 
ky += kValy[np.newaxis, :, np.newaxis] 
kz += kValz[np.newaxis, np.newaxis, :]