加快对 numpy 中数组的分析
Speeding up analysis on arrays in numpy
我有一个 python 代码，它导入包含数字的 4 列 txt 文件
前三列是 x、y、z 坐标，第四列是该坐标处的密度。
下面是读取的代码，转换为 ndarray，对该字段进行傅里叶变换，计算距原点的距离 (k=(0,0,0)) 和变换后的坐标，取平均值，然后绘制它们。
多亏了 pandas（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
Numpy 通常可以比普通 python 快数百倍的速度，只需很少的额外努力（有时甚至会自动使用多个内核，而无需您的工作）。您只需要知道编写代码的正确方法。仅举出我想到的第一件事：
- 考虑一下您的方法 — 简化它并确保它在编码之前是正确的
- 学习使用 numpy 的索引
- 尽可能不要复制数据
- 不惜一切代价尽量避免循环
- 使用 numpy 的内置函数
- 处理整个数组（或大块数组），而不是一次处理单个元素
- 分析您的代码以找出慢点
- 使用numba，尤其是对于您无法避免的循环或 numpy 无法帮助解决的循环
你的方法
我认为我们都会犯的一个错误是试图在代码中只写，即使我们仍在试图弄清楚代码应该做什么——基本上，说一种我们不太了解的语言，它不是真正的描述性语言，而且单词数量非常有限。我发现最好首先用简单的英语将我的方法概括为一系列步骤。这让我清楚地看到，用我自己的话来说，正是我想要代码做什么。因此，很容易获得概览并思考正在发生的事情，但如果我意识到某些东西需要返工，也很容易重写。
事实上，我将这些步骤作为注释写在代码文件的不同行上。一旦我对这种方法感到满意，我就开始在这些评论之间编写代码本身。这给我留下了结构良好且注释正确的代码。
在 OP 中原始代码的特定示例中，整个 f 数组可能不需要存在。如果您只是用简单的英语编写程序，我想您会看到这一点，并且能够在编写代码之前轻松地 in plain English 更改您的方法。就目前而言，如此多的代码依赖于 f，需要完全重写才能更改它。
索引
Plain python 通常在计算机应该擅长的事情上非常慢。一个例子是索引，所以一行
a[f[i,0]][f[i,1]][f[i,2]]=f[i,3]
让我很怀疑。当您说 "transforming the loaded txt to 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.
中阅读相关信息
更具体地说，如果我正确理解代码，前三部分代码可以完全替换为
a = np.random.rand(ngridx,ngridy,ngridz)
也就是说，据我所知，f根本不需要存在。
复制数据错误
你不需要复制所有东西，当你确实需要复制一个数组（或数组的一部分）时，你应该让 numpy 为你做。例如，不用 Firstdel 函数，只需使用 a[1:]:
def Firstdel(a):
return a[1:]
这根本不会复制数据；它 returns 只是忽略了另一个数字存储在它之前的 RAM 中的事实。或者，如果你真的需要复制数据（你不只是为了绘图）使用这个：
def Firstdel(a):
return numpy.copy(a[1:])
但通常情况下，您可以只使用 "slices" 个 numpy 数组，而不是复制它们。了解它 here。
但是，还要注意，如果您只获取一个切片（并且没有复制发生），这意味着如果您写入切片，原始数组也会被写入。
循环
循环也是出了名的浪费时间——特别是当这些循环在 python 代码中，而不是某些 C 代码中时。事实上，我认为 numpy 的全部意义在于用他们用 C 语言为你编写的循环替换你将在 python 中编写的循环。
首先，对于简单循环，while在python中并不常见。所以而不是
while i < len(f):
# do stuff
i = i+1
你可能应该使用
for i in range(len(f)):
# do stuff
但这只是风格和可读性的问题。为了提高速度，您需要尽可能多地摆脱循环——尤其是嵌套循环。要设置 kx、ky 和 kz，此代码比 OP 中参数的嵌套循环快大约 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] = np.arange(len(q))
在这里，我只是将 q 的 "slice" 设置为另一个数组。
该循环之后的嵌套循环也可以加速，但它们会稍微复杂一些。
但您还希望尽可能避免循环。这将我们带到...
使用内置的 numpy 函数
同样，numpy 存在的主要原因是将这些慢速 python 循环转换为快速 C 代码（我们不需要知道它的存在）。所以有很多函数可以做你想做的事情已经内置到 numpy 中。
而不是
while i < len(f):
masstot = masstot + f[i,3]
i = i+1
你应该使用像
这样的东西
masstot = np.sum(f[:,3])
阅读起来更简单，但也可能 方式 更快，因为 numpy 可以直接访问计算机内存中的数据，并且可以使用快速 C求和的函数，而不是使用缓慢的 python 函数。 （同样，您无需了解 C 函数；它们会完成工作。）
处理整个数组（或切片）
一个常见的错误是使用 numpy 函数，但仍然自己执行循环。这段代码中有一个很好的例子：
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))
与其每次通过循环计算 K 的大嵌套循环不同，只需将 K 设为具有适当值的数组：
K = np.around(np.sqrt(kx**2+ky**2+kz**2))
K 将与 kx、ky 和 kz 的大小和形状相同。然后，您可以使用 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):
p_i = p[K==i] # This will be just the elements of p where K==i
q[i,0] = i
q[i,1] = np.sum(np.abs(p_i)**2)
q[i,2] = np.sum(p_i)
q[i,3] = len(p_i)
print q
将这些技巧放在一起，我发现您的问题中的代码比当前代码提高了大约 70 倍。如果它仍然太慢，您可以使用分析器找到慢的部分。如果您仍然无法改进它们，可能值得尝试 numba，它会根据需要编译代码，运行 几乎与等效的 C 代码一样快。
还可以加快 输入文件 的创建速度：
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, :]
我有一个 python 代码，它导入包含数字的 4 列 txt 文件 前三列是 x、y、z 坐标，第四列是该坐标处的密度。
下面是读取的代码，转换为 ndarray，对该字段进行傅里叶变换，计算距原点的距离 (k=(0,0,0)) 和变换后的坐标，取平均值，然后绘制它们。 多亏了 pandas（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
Numpy 通常可以比普通 python 快数百倍的速度，只需很少的额外努力（有时甚至会自动使用多个内核，而无需您的工作）。您只需要知道编写代码的正确方法。仅举出我想到的第一件事：
- 考虑一下您的方法 — 简化它并确保它在编码之前是正确的
- 学习使用 numpy 的索引
- 尽可能不要复制数据
- 不惜一切代价尽量避免循环
- 使用 numpy 的内置函数
- 处理整个数组（或大块数组），而不是一次处理单个元素
- 分析您的代码以找出慢点
- 使用numba，尤其是对于您无法避免的循环或 numpy 无法帮助解决的循环
你的方法
我认为我们都会犯的一个错误是试图在代码中只写，即使我们仍在试图弄清楚代码应该做什么——基本上，说一种我们不太了解的语言，它不是真正的描述性语言，而且单词数量非常有限。我发现最好首先用简单的英语将我的方法概括为一系列步骤。这让我清楚地看到，用我自己的话来说，正是我想要代码做什么。因此，很容易获得概览并思考正在发生的事情，但如果我意识到某些东西需要返工，也很容易重写。
事实上，我将这些步骤作为注释写在代码文件的不同行上。一旦我对这种方法感到满意，我就开始在这些评论之间编写代码本身。这给我留下了结构良好且注释正确的代码。
在 OP 中原始代码的特定示例中，整个 f 数组可能不需要存在。如果您只是用简单的英语编写程序，我想您会看到这一点，并且能够在编写代码之前轻松地 in plain English 更改您的方法。就目前而言，如此多的代码依赖于 f，需要完全重写才能更改它。
索引
Plain python 通常在计算机应该擅长的事情上非常慢。一个例子是索引，所以一行
a[f[i,0]][f[i,1]][f[i,2]]=f[i,3]
让我很怀疑。当您说 "transforming the loaded txt to 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.
中阅读相关信息更具体地说，如果我正确理解代码，前三部分代码可以完全替换为
a = np.random.rand(ngridx,ngridy,ngridz)
也就是说，据我所知，f根本不需要存在。
复制数据错误
你不需要复制所有东西，当你确实需要复制一个数组（或数组的一部分）时，你应该让 numpy 为你做。例如，不用 Firstdel 函数，只需使用 a[1:]:
def Firstdel(a):
return a[1:]
这根本不会复制数据；它 returns 只是忽略了另一个数字存储在它之前的 RAM 中的事实。或者，如果你真的需要复制数据（你不只是为了绘图）使用这个：
def Firstdel(a):
return numpy.copy(a[1:])
但通常情况下，您可以只使用 "slices" 个 numpy 数组，而不是复制它们。了解它 here。
但是，还要注意，如果您只获取一个切片（并且没有复制发生），这意味着如果您写入切片，原始数组也会被写入。
循环
循环也是出了名的浪费时间——特别是当这些循环在 python 代码中，而不是某些 C 代码中时。事实上，我认为 numpy 的全部意义在于用他们用 C 语言为你编写的循环替换你将在 python 中编写的循环。
首先，对于简单循环，while在python中并不常见。所以而不是
while i < len(f):
# do stuff
i = i+1
你可能应该使用
for i in range(len(f)):
# do stuff
但这只是风格和可读性的问题。为了提高速度，您需要尽可能多地摆脱循环——尤其是嵌套循环。要设置 kx、ky 和 kz，此代码比 OP 中参数的嵌套循环快大约 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] = np.arange(len(q))
在这里，我只是将 q 的 "slice" 设置为另一个数组。
该循环之后的嵌套循环也可以加速，但它们会稍微复杂一些。
但您还希望尽可能避免循环。这将我们带到...
使用内置的 numpy 函数
同样，numpy 存在的主要原因是将这些慢速 python 循环转换为快速 C 代码（我们不需要知道它的存在）。所以有很多函数可以做你想做的事情已经内置到 numpy 中。
而不是
while i < len(f):
masstot = masstot + f[i,3]
i = i+1
你应该使用像
这样的东西masstot = np.sum(f[:,3])
阅读起来更简单，但也可能 方式 更快，因为 numpy 可以直接访问计算机内存中的数据，并且可以使用快速 C求和的函数，而不是使用缓慢的 python 函数。 （同样，您无需了解 C 函数；它们会完成工作。）
处理整个数组（或切片）
一个常见的错误是使用 numpy 函数，但仍然自己执行循环。这段代码中有一个很好的例子：
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))
与其每次通过循环计算 K 的大嵌套循环不同，只需将 K 设为具有适当值的数组：
K = np.around(np.sqrt(kx**2+ky**2+kz**2))
K 将与 kx、ky 和 kz 的大小和形状相同。然后，您可以使用 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):
p_i = p[K==i] # This will be just the elements of p where K==i
q[i,0] = i
q[i,1] = np.sum(np.abs(p_i)**2)
q[i,2] = np.sum(p_i)
q[i,3] = len(p_i)
print q
将这些技巧放在一起，我发现您的问题中的代码比当前代码提高了大约 70 倍。如果它仍然太慢，您可以使用分析器找到慢的部分。如果您仍然无法改进它们，可能值得尝试 numba，它会根据需要编译代码，运行 几乎与等效的 C 代码一样快。
还可以加快 输入文件 的创建速度：
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, :]