【发布时间】：2021-11-09 20:42:14
【问题描述】：
我一直在使用 Anaconda 执行一些简单的引力计算，并且我想使用 Numba 加快计算过程。因此，我将原来的 Python 对象（列表）结果修改为仅适用于 numpy 数组的结果，然后使用 @jit 或 @njit 来加快计算速度。它似乎适用于相当短的计算，但它突然因更长的计算以及它如何与其他代码混合而崩溃。我使用的模块是；
%matplotlib notebook
import numpy as np
from matplotlib import pyplot as plt
import pandas as pd
from timeit import default_timer as timer
from matplotlib import colors
from matplotlib.widgets import Slider, Button
import re
import os
from numba import jit, njit
python 物体的引力代码是；
def verlet_int(u,p,v,t_set,N,k,d,p_num):
t = t_set[0]
T = t_set[1]
dt = t_set[2]
p_prime = []
v_prime = []
o = []
a_ik = lambda u,p,k,i,j : (-1)*u[k]*((((p[i][0] - p[k][0])**2) + ((p[i][1] - p[k][1])**2) + ((p[i][2] - p[k][2])**2))**(-1.5))*(p[i][j] - p[k][j])
while t <= T:
a_1 = []
a_2 = []
for i in np.arange(N):
b_1 = []
b_2 = []
for j in [0,1,2]:
if p_num == -1:
b_1.append(p[i][j])
b_2.append(v[i][j])
elif i == p_num:
a_1.append(p[i][j])
a_2.append(v[i][j])
v_half = v[i][j] + 0.5*sum([a_ik(u,p,k,i,j) for k in np.delete(np.arange(N),i)])*dt
p[i][j] = p[i][j] + v_half*dt
v[i][j] = v_half + 0.5*sum([a_ik(u,p,k,i,j) for k in np.delete(np.arange(N),i)])*dt
if p_num == -1:
a_1.append(b_1)
a_2.append(b_2)
elif i == p_num:
pass
if round(t/dt)%k == 0:
p_prime.append(a_1)
v_prime.append(a_2)
else:
pass
t += dt
t = np.round(t, decimals = d)
o.append(p_prime)
o.append(v_prime)
return o
在哪里，o[positions (0) or velocities (1)][pick your recorded time step][pick your particle 0...N-1][x (0), y(1), or z(2)]。也就是说，如果p_num = -1，如果不是，那么它只记录该粒子p_num = 0...N-1 的位置和速度。
我jit的修改后的代码是；
@jit(debug = True)
def verlet_int_numpy(u,p_v_1,t_i,T,dt,k,dec,p_num):
times = np.arange(t_i,T,dt)
p_v = np.copy(p_v_1)
if p_num == -1:
storage = np.zeros((round(times.size/k),2,len(u),3))
else:
storage = np.zeros((round(times.size/k),2,1,3))
for tau in range(times.size):
if tau%k == 0:
if p_num == -1:
storage[tau] = p_v
else:
storage[tau][0][0] = p_v[0][p_num]
storage[tau][0][1] = p_v[1][p_num]
else:
pass
for j in range(u.size):
a_g = np.array([0.0,0.0,0.0])
for i in np.delete(np.arange(u.size),j):
a_g = a_g + u[i]*(((p_v[0][i] - p_v[0][j])@(p_v[0][i] - p_v[0][j]))**(-1.5))*(p_v[0][i] - p_v[0][j])
v_half = p_v[1][j] + 0.5*a_g*dt
p_v[0][j] = p_v[0][j] + v_half*dt
a_g = np.array([0.0,0.0,0.0])
for i in np.delete(np.arange(u.size),j):
a_g = a_g + u[i]*(((p_v[0][i] - p_v[0][j])@(p_v[0][i] - p_v[0][j]))**(-1.5))*(p_v[0][i] - p_v[0][j])
p_v[1][j] = v_half + 0.5*a_g*dt
return storage
然后我在另一个函数中使用了它，它将每个粒子的位置和速度分别记录为熊猫数据帧到我的外部驱动器。这个存在；
def setup_csv(i_s,i_v,strings):
r_cm = (1/np.sum(i_s[0]))*np.dot(i_s[0],i_v[0])
v_cm = (1/np.sum(i_s[0]))*np.dot(i_s[0],i_v[1])
i_v[0] = i_v[0] - r_cm
i_v[1] = i_v[1] - v_cm
O = verlet_int_numpy(np.array(i_s[0]),i_v,i_s[1],i_s[2],i_s[3],i_s[6],i_s[4],i_s[7])
if i_s[7] == -1:
for j in np.arange(i_s[0].size):
D_1 = pd.DataFrame(data = O[:,0,j])
D_2 = pd.DataFrame(data = O[:,1,j])
D_1.to_csv(r"/Volumes/Storage/Physics Simulations and Research/Gravitational Research/D_test/" + "particle_{}_{}_{}_".format(j,i_s[6],i_s[3]) + strings[0], mode = 'w', header = False, index = None)
D_2.to_csv(r"/Volumes/Storage/Physics Simulations and Research/Gravitational Research/D_test/" + "particle_{}_{}_{}_".format(j,i_s[6],i_s[3]) + strings[1], mode = 'w', header = False, index = None)
else:
D_1 = pd.DataFrame(data = O[:,0,i_s[7]])
D_2 = pd.DataFrame(data = O[:,1,i_s[7]])
D_1.to_csv(r"/Volumes/Storage/Physics Simulations and Research/Gravitational Research/D_test/" + "particle_{}_{}_{}_".format(i_s[7],i_s[6],i_s[3]) + strings[0], mode = 'w', header = False, index = None)
D_2.to_csv(r"/Volumes/Storage/Physics Simulations and Research/Gravitational Research/D_test/" + "particle_{}_{}_{}_".format(i_s[7],i_s[6],i_s[3]) + strings[1], mode = 'w', header = False, index = None)
关于术语含义的更多信息在下面给出了python列表一；
#p_pick gives the particle to be recorded from orbit_int_verlet(). It takes values 0,1,...,N-1. It can also take on
#. . . -1 for recording all the possible particle positions/velocites.
#s_scale gives the magnitude of the random unit vector.
#t_i is the initial time.
#Usually assume t_o = 0.
#dt is the time step.
#T = total length of time.
#k decides the number of stored values by whether it evenly divides the number of time steps gone by.
#u are the 'N' number of masses for the bodies.
#d_places is the number of decimal places to round the time step every loop so as to not accumlate much numerical error.
#i_s = inital_scalars = [u,t_i,T,dt,d_places,s_scale,k,p_pick]
#strings = [] . . . set of string titles. The first two will be the comparison orbits; position then velocity.
#. . . the last two will be where the output files for position then velocity will go
#i_v = initial_vectors = np.array([ [[],[],...,[]], [[],[],...,[]] ]) this is 2xNx3 array.
numpy 版本有它，所以 'u' 也是一个 numpy 数组而不是一个集合。结合初始条件时；
u = np.array([0.97,0.02,0.01])
orbits = 0.05
dec = 3
T = round( (2*np.pi*(80**(1.5)))/(np.sqrt(0.98))*orbits)
dt = 0.005
p_num = -1
s_scale = 0.1
k = 950
t_i = 0
i_s = [u,t_i,T,dt,dec,s_scale,k,p_num]
i_v = np.array([[[20,0,0],[-40,0,0],[80,0,0]],[[0,0,0],[0,np.sqrt(0.85/40),0],[0,np.sqrt(0.85/80),0]]])
strings = ["comparison_p.csv","comparison_v.csv"]
然后当它运行时，它宁愿很快有一个内核已经死了，然后说它会自动重启。我认为 jitting 是有效的，这样做时没有例外，但是关于运行它的一些事情就是不合适。我不知道我做错了什么或此时发生了什么，并且非常感谢这里的一些方向。
编辑：
在之前的运行中出现了这个错误；
---------------------------------------------------------------------------
AttributeError Traceback (most recent call last)
/opt/anaconda3/lib/python3.7/site-packages/IPython/core/async_helpers.py in _pseudo_sync_runner(coro)
66 """
67 try:
---> 68 coro.send(None)
69 except StopIteration as exc:
70 return exc.value
/opt/anaconda3/lib/python3.7/site-packages/IPython/core/interactiveshell.py in run_cell_async(self, raw_cell, store_history, silent, shell_futures, transformed_cell, preprocessing_exc_tuple)
3166 # Write output to the database. Does nothing unless
3167 # history output logging is enabled.
-> 3168 self.history_manager.store_output(self.execution_count)
3169 # Each cell is a *single* input, regardless of how many lines it has
3170 self.execution_count += 1
AttributeError: 'ZMQInteractiveShell' object has no attribute 'history_manager'
【问题讨论】：