英文:
Increasing the efficiency of a python code by using arrays
问题
以下是您提供的Python代码的中文翻译：
import matplotlib.pyplot as plt
import numpy as np
class bifurcation_diagram(object):
def __init__(self):
self.omega = []
self.theta = []
self.dt = (2 * np.pi / (2.0 / 3.0)) / 600
self.time = []
self.theta_in_bifurcation_diagram = []
self.F_D = np.arange(1.35, 1.5, 0.001)
self.theta_versus_F_D = []
def calculate(self):
l = 9.8
g = 9.8
q = 0.5
Omega_D = 2.0 / 3.0
for f_d in self.F_D:
self.omega.append([0])
self.theta.append([0.2])
self.time.append([0])
for i in range(600000):
k1_theta = self.dt * self.omega[-1][-1]
k1_omega = self.dt * ((-g / l) * np.sin(self.theta[-1][-1]) - q * self.omega[-1][-1] + f_d * np.sin(Omega_D * self.time[-1][-1]))
k2_theta = self.dt * (self.omega[-1][-1] + 0.5 * k1_omega)
k2_omega = self.dt * ((-g / l) * np.sin(self.theta[-1][-1] + 0.5 * k1_theta) - q * (self.omega[-1][-1] + 0.5 * k1_omega) + f_d * np.sin(Omega_D * (self.time[-1][-1] + 0.5 * self.dt)))
k3_theta = self.dt * (self.omega[-1][-1] + 0.5 * k2_omega)
k3_omega = self.dt * ((-g / l) * np.sin(self.theta[-1][-1] + 0.5 * k2_theta) - q * (self.omega[-1][-1] + 0.5 * k2_omega) + f_d * np.sin(Omega_D * (self.time[-1][-1] + 0.5 * self.dt)))
k4_theta = self.dt * (self.omega[-1][-1] + k3_omega)
k4_omega = self.dt * ((-g / l) * np.sin(self.theta[-1][-1] + k3_theta) - q * (self.omega[-1][-1] + k3_omega) + f_d * np.sin(Omega_D * (self.time[-1][-1] + self.dt)))
temp_theta = self.theta[-1][-1] + (1.0 / 6.0) * (k1_theta + 2 * k2_theta + 2 * k3_theta + k4_theta)
temp_omega = self.omega[-1][-1] + (1.0 / 6.0) * (k1_omega + 2 * k2_omega + 2 * k3_omega + k4_omega)
while temp_theta > np.pi:
temp_theta -= 2 * np.pi
while temp_theta < -np.pi:
temp_theta += 2 * np.pi
self.omega[-1].append(temp_omega)
self.theta[-1].append(temp_theta)
self.time[-1].append(self.dt * i)
for i in range(500, 1000):
n = i * 600
self.theta_in_bifurcation_diagram.append(self.theta[-1][n])
self.theta_versus_F_D.append(f_d)
def show_results(self):
plt.plot(self.theta_versus_F_D, self.theta_in_bifurcation_diagram, '.')
plt.title('分岔图' + '\n' + r'$\theta$ versus $F_D$')
plt.xlabel(r'$F_D$')
plt.ylabel(r'$\theta$ (弧度)')
plt.xlim(1.35, 1.5)
plt.ylim(1, 3)
plt.show()
bifurcation = bifurcation_diagram()
bifurcation.calculate()
bifurcation.show_results()
如果您希望缩短代码并提高效率，可以考虑使用NumPy数组操作以减少循环的使用，从而提高计算速度。要添加颜色，您可以在plt.plot函数中使用不同的颜色参数来绘制彩色图形。此外，您可以考虑使用函数或方法来更好地组织代码。
英文:
I have the following python code:
import matplotlib.pyplot as plt
import numpy as np
class bifurcation_diagram(object):
def __init__(self):
self.omega = []
self.theta = []
self.dt = (2 * np.pi / (2.0 / 3.0)) / 600
self.time = []
self.theta_in_bifurcation_diagram = []
self.F_D = np.arange(1.35,1.5,0.001)
self.theta_versus_F_D = []
def calculate(self):
l = 9.8
g = 9.8
q = 0.5
Omega_D = 2.0 / 3.0
for f_d in self.F_D:
self.omega.append([0])
self.theta.append([0.2])
self.time.append([0])
for i in range(600000):
k1_theta = self.dt * self.omega[-1][-1]
k1_omega = self.dt * ((-g / l) * np.sin(self.theta[-1][-1]) - q * self.omega[-1][-1] + f_d * np.sin(Omega_D * self.time[-1][-1]))
k2_theta = self.dt * (self.omega[-1][-1] + 0.5 * k1_omega)
k2_omega = self.dt * ((-g / l) * np.sin(self.theta[-1][-1] + 0.5 * k1_theta) - q * (self.omega[-1][-1] + 0.5 * k1_omega) + f_d * np.sin(Omega_D * (self.time[-1][-1] + 0.5 * self.dt)))
k3_theta = self.dt * (self.omega[-1][-1] + 0.5 * k2_omega)
k3_omega = self.dt * ((-g / l) * np.sin(self.theta[-1][-1] + 0.5 * k2_theta) - q * (self.omega[-1][-1] + 0.5 * k2_omega) + f_d * np.sin(Omega_D * (self.time[-1][-1] + 0.5 * self.dt)))
k4_theta = self.dt * (self.omega[-1][-1] + k3_omega)
k4_omega = self.dt * ((-g / l) * np.sin(self.theta[-1][-1] + k3_theta) - q * (self.omega[-1][-1] + k3_omega) + f_d * np.sin(Omega_D * (self.time[-1][-1] + self.dt)))
temp_theta = self.theta[-1][-1] + (1.0 / 6.0) * (k1_theta + 2 * k2_theta + 2 * k3_theta + k4_theta)
temp_omega = self.omega[-1][-1] + (1.0 / 6.0) * (k1_omega + 2 * k2_omega + 2 * k3_omega + k4_omega)
while temp_theta > np.pi:
temp_theta -= 2 * np.pi
while temp_theta < -np.pi:
temp_theta += 2 * np.pi
self.omega[-1].append(temp_omega)
self.theta[-1].append(temp_theta)
self.time[-1].append(self.dt * i)
for i in range(500,1000):
n = i * 600
self.theta_in_bifurcation_diagram.append(self.theta[-1][n])
self.theta_versus_F_D.append(f_d)
def show_results(self):
plt.plot(self.theta_versus_F_D,self.theta_in_bifurcation_diagram,'.')
plt.title('Bifurcation diagram' + '\n' + r'$\theta$ versus $F_D$')
plt.xlabel(r'$F_D$')
plt.ylabel(r'$\theta$ (radians)')
plt.xlim(1.35,1.5)
plt.ylim(1,3)
plt.show()
bifurcation = bifurcation_diagram()
bifurcation.calculate()
bifurcation.show_results()
I wish to make it more compact and make a coloured plot, along with increasing its efficiency. Any help in this regard would be truly beneficial.
I was expecting the code to run fast but it takes more than fifteen minutes to run.
答案1
得分: 2
以下是您提供的代码的翻译部分：
Python is not made for such a code because it is generally **interpreted** using CPython (it is meant to be done for glue codes and scripts) and CPython performs (nearly) no optimizations of the code so repeated expressions are recomputed. The usual way to speed up such a code is to use **Numpy** functions operating on large arrays (not lists). THis is called **vectorization**. That being said, this code is pretty expensive and using Numpy may not be enough and it is also a non-negligible work to use Numpy here. An alternative solution is to use **Numba** so to compile this code (thanks to a JIT compiler). Note that Numba is meant to speed-up codes using mostly Numpy and basic numerical computations (not general-purpose codes like plotting).
首先要做的是摆脱列表，并将它们替换为Numpy数组。数组可以预先分配到正确的大小，以避免昂贵的`append`调用。然后，计算函数应该从类中移出，以便轻松使用Numba。然后，两个`while`循环可以替换为`np.remainder`。最后，可以使用多线程并行计算数组的每一行。
这是最终的代码（几乎未经测试）：
```python
import matplotlib.pyplot as plt
import numpy as np
import numba as nb
@nb.njit('(float64[::1], float64)', parallel=True)
def compute_numba(F_D, dt):
# ...
# (代码中的其余部分)
class bifurcation_diagram(object):
# ...
# (代码中的其余部分)
bifurcation = bifurcation_diagram()
bifurcation.calculate()
bifurcation.show_results()
在我的i5-9600KF处理器上，这段代码执行时间为1.07秒，而初始代码执行时间约为19分20秒！因此，这个新代码快约1150倍！
此外，我还发现最终代码更容易阅读。
我建议您检查结果，尽管乍看之下看起来不错。以下是我迄今为止获得的结果图表：
请注意，这个翻译是您提供的代码的中文版本，不包括代码中的注释或解释性文本。
<details>
<summary>英文:</summary>
Python is not made for such a code because it is generally **interpreted** using CPython (it is meant to be done for glue codes and scripts) and CPython performs (nearly) no optimizations of the code so repeated expressions are recomputed. The usual way to speed up such a code is to use **Numpy** functions operating on large arrays (not lists). THis is called **vectorization**. That being said, this code is pretty expensive and using Numpy may not be enough and it is also a non-negligible work to use Numpy here. An alternative solution is to use **Numba** so to compile this code (thanks to a JIT compiler). Note that Numba is meant to speed-up codes using mostly Numpy and basic numerical computations (not general-purpose codes like plotting).
The first thing to do is to get ride of lists and replace them by Numpy arrays. Arrays can be preallocated at the right size so to avoid the expensive calls to `append`. Then, the computing function should be moved away from the class so for Numba to be easily used. Then, the two `while` loops can be replaced by `np.remainder`. Finally, you can compute each row of the arrays in **parallel** using multiple threads.
Here is the resulting code (barely tested):
```python
import matplotlib.pyplot as plt
import numpy as np
import numba as nb
@nb.njit('(float64[::1], float64)', parallel=True)
def compute_numba(F_D, dt):
l = 9.8
g = 9.8
q = 0.5
Omega_D = 2.0 / 3.0
size = 600000
loop2_start, loop2_end = 500, 1000
loop2_stride = 600
omega = np.empty((F_D.size, size+1), dtype=np.float64)
theta = np.empty((F_D.size, size+1), dtype=np.float64)
time = np.empty((F_D.size, size+1), dtype=np.float64)
theta_in_bifurcation_diagram = np.empty((loop2_end-loop2_start) * F_D.size, dtype=np.float64)
theta_versus_F_D = np.empty((loop2_end-loop2_start) * F_D.size, dtype=np.float64)
for i in nb.prange(F_D.size):
f_d = F_D[i]
omega[i, 0] = 0.0
theta[i, 0] = 0.2
time[i, 0] = 0.0
for j in range(size):
cur_omega = omega[i,j]
cur_theta = theta[i,j]
cur_time = time[i,j]
tmp = np.sin(Omega_D * (cur_time + 0.5 * dt))
k1_theta = dt * cur_omega
k1_omega = dt * ((-g / l) * np.sin(cur_theta) - q * cur_omega + f_d * np.sin(Omega_D * cur_time))
k2_theta = dt * (cur_omega + 0.5 * k1_omega)
k2_omega = dt * ((-g / l) * np.sin(cur_theta + 0.5 * k1_theta) - q * (cur_omega + 0.5 * k1_omega) + f_d * tmp)
k3_theta = dt * (cur_omega + 0.5 * k2_omega)
k3_omega = dt * ((-g / l) * np.sin(cur_theta + 0.5 * k2_theta) - q * (cur_omega + 0.5 * k2_omega) + f_d * tmp)
k4_theta = dt * (cur_omega + k3_omega)
k4_omega = dt * ((-g / l) * np.sin(cur_theta + k3_theta) - q * (cur_omega + k3_omega) + f_d * np.sin(Omega_D * (cur_time + dt)))
temp_theta = cur_theta + (1.0 / 6.0) * (k1_theta + 2 * k2_theta + 2 * k3_theta + k4_theta)
temp_omega = cur_omega + (1.0 / 6.0) * (k1_omega + 2 * k2_omega + 2 * k3_omega + k4_omega)
temp_theta = np.remainder(temp_theta + np.pi, 2 * np.pi) - np.pi
omega[i,j+1] = temp_omega
theta[i,j+1] = temp_theta
time[i,j+1] = dt * j
for k, j in enumerate(range(loop2_start, loop2_end)):
n = j * loop2_stride
offset = (loop2_end - loop2_start) * i + k
theta_in_bifurcation_diagram[offset] = theta[i,n]
theta_versus_F_D[offset] = f_d
return (omega, theta, time, theta_in_bifurcation_diagram, theta_versus_F_D)
class bifurcation_diagram(object):
def __init__(self):
self.omega = np.array([])
self.theta = np.array([])
self.dt = (2 * np.pi / (2.0 / 3.0)) / 600
self.time = np.array([])
self.theta_in_bifurcation_diagram = np.array([])
self.F_D = np.arange(1.35,1.5,0.001)
self.theta_versus_F_D = np.array([])
def calculate(self):
self.omega, self.theta, self.time, self.theta_in_bifurcation_diagram, self.theta_versus_F_D = compute_numba(self.F_D, self.dt)
def show_results(self):
plt.plot(self.theta_versus_F_D,self.theta_in_bifurcation_diagram,'.')
plt.title('Bifurcation diagram' + '\n' + r'$\theta$ versus $F_D$')
plt.xlabel(r'$F_D$')
plt.ylabel(r'$\theta$ (radians)')
plt.xlim(1.35,1.5)
plt.ylim(1,3)
plt.show()
bifurcation = bifurcation_diagram()
bifurcation.calculate()
bifurcation.show_results()
On my i5-9600KF processor with 6 cores, this code takes 1.07 secondes rather than about 19 min 20 for the initial one! As a result, this new code is about 1150 times faster!
Besides, I also find the resulting code easier to read.
I advise you to check the results, although they look fine at first glance.
Indeed, here is the resulting graph I got so far:
通过集体智慧和协作来改善编程学习和解决问题的方式。致力于成为全球开发者共同参与的知识库，让每个人都能够通过互相帮助和分享经验来进步。
评论