为什么部署 Numba 并使用数组而不是 class 会为我的程序中的相同参数提供不同的结果？
Why deploying Numba and using arrays instead of class gives different results for the same parameters in my program?
总而言之，我的代码具有完全相同的参数，但使用了三种不同的方法：
正在为我的粒子 SP 列表定义一个 class，其中包含三个 instances/attributes。 （大约 1 分钟）
正在将 SP 定义为数组。 （大约 1 分钟）
将 SP 定义为数组并使用 Numba（大约 5 秒）
考虑第一种情况：
n=1000
mu=np.random.uniform(0,1,n)
r=[sqrt(-2*log(1-i)) for i in mu]
eta=np.random.uniform(0,1,n)
theta=2*pi*eta;
cuz=[cos(i) for i in theta]
suz=[sin(i) for i in theta]
Zinitial=[a*b for a,b in zip(r,cuz)];
Pinitial=[a*b for a,b in zip(r,suz)];
class Particle:
def __init__(self, pos, mom, spin):
self.pos = pos
self.mom = mom
self.spin = spin
SP = sorted([Particle(pos = i, mom = j, spin = choice([1, 0])) for i,j in zip(Zinitial,Pinitial)],key=lambda x:x.pos)
Upi=[];
Downi=[];
count_plot=[];
for j in range(len(SP)):
if SP[j].spin == 1:
Upi.append(SP[j].pos)
else:
Downi.append(SP[j].pos)
Zavgi=sum(Zinitial)/len(Zinitial)
Zreli=sum(Upi)/len(Upi)-sum(Downi)/len(Downi)
"Observables"
Zavg=[Zavgi];
Zrelm=[Zreli];
T_plot=[0];
"Time"
iter=10**(4);
dt=1/(2*n);
alf=sqrt(n);
"Dynamics"
counter=0;
sum1,sum2=0,0;
for i in range(1,iter+1):
t=i*dt;
T_plot.append(t)
Z=[];
Up=[];
Down=[];
c,s=cos(t),sin(t);
c1,s1=cos(t-dt),sin(t-dt);
for j in range(n-1):
collchk=((c*(SP[j].pos)+s*(SP[j].mom))-(c*(SP[j+1].pos)+s*(SP[j+1].mom)))*(c1*(SP[j].pos)+s1*(SP[j].mom)-(c1*(SP[j+1].pos)+s1*(SP[j+1].mom)));
prel=((c*(SP[j].mom)-s*(SP[j].pos))-(c*(SP[j+1].mom)-s*(SP[j+1].pos)))/2;
rcoeff=1/(1+(prel*alf)**2);
rand_value=random();
if collchk<0:
SP[j], SP[j+1]=SP[j+1],SP[j];
if rcoeff>rand_value:
counter=counter+1
SP[j].spin,SP[j+1].spin=SP[j+1].spin,SP[j].spin;
if SP[j].spin == 1:
Up.append(c*(SP[j].pos)+s*(SP[j].mom))
else:
Down.append(c*(SP[j].pos)+s*(SP[j].mom))
Z.append(c*(SP[j].pos)+s*(SP[j].mom))
Zrel=sum(Up[0:])/len(Up) - sum(Down[0:])/len(Down);
Zrelm.append(Zrel)
Zm=sum(Z)/len(Z)
Zavg.append(Zm)
print("Rate of collision per particle = ",counter/(n*(dt*iter)))
输出是：
每个粒子的碰撞率 = 0.0722
第二种情况：
n=1000
mu=np.random.uniform(0,1,n)
r=np.array([sqrt(-2*log(1-i)) for i in mu])
eta=np.random.uniform(0,1,n)
theta=2*pi*eta;
cuz=np.array([cos(i) for i in theta])
suz=np.array([sin(i) for i in theta])
Zinitial=np.multiply(r,cuz);
Pinitial=np.multiply(r,suz);
SP = np.array(sorted(np.array([ np.array([i,j,choice([1,0])]) for i, j in zip(Zinitial, Pinitial)]),
key=lambda x: x[0]))
Upi=[];
Downi=[];
count_plot=[];
for j in range(len(SP)):
if SP[j][2] == 1:
Upi.append(SP[j][0])
else:
Downi.append(SP[j][0])
Zavgi=sum(Zinitial)/len(Zinitial)
Zreli=sum(Upi)/len(Upi)-sum(Downi)/len(Downi)
"Observables"
Zavg=[Zavgi];
Zrelm=[Zreli];
T_plot=[0];
"Time"
iter=10**(4);
dt=1/(2*n);
alf=sqrt(n);
"Dynamics"
counter=0;
sum1,sum2=0,0;
for i in range(1,iter+1):
t=i*dt;
T_plot.append(t)
Z=[];
Up=[];
Down=[];
c,s=cos(t),sin(t);
c1,s1=cos(t-dt),sin(t-dt);
for j in range(n-1):
collchk=((c*(SP[j][0])+s*(SP[j][1]))-(c*(SP[j+1][0])+s*(SP[j+1][1])))*(c1*(SP[j][0])+s1*(SP[j][1])-(c1*(SP[j+1][0])+s1*(SP[j+1][1])));
prel=((c*(SP[j][1])-s*(SP[j][0]))-(c*(SP[j+1][1])-s*(SP[j+1][0])))/2;
rcoeff=1/(1+(prel*alf)**2);
rand_value=random();
if collchk<0:
SP[j], SP[j+1]=SP[j+1],SP[j];
if rcoeff>rand_value:
counter=counter+1
SP[j][2],SP[j+1][2]=SP[j+1][2],SP[j][2];
if SP[j][2] == 1:
Up.append(c*(SP[j][0])+s*(SP[j][1]))
else:
Down.append(c*(SP[j][0])+s*(SP[j][1]))
Z.append(c*(SP[j][0])+s*(SP[j][1]))
Zrel=sum(Up[0:])/len(Up) - sum(Down[0:])/len(Down);
Zrelm.append(Zrel)
Zm=sum(Z)/len(Z)
Zavg.append(Zm)
print("Rate of collision per particle = ",counter/(n*(dt*iter)))
输出是：
每个粒子的碰撞率 = 0.0134
以及最快的 Numba 案例；
import numpy as np
import matplotlib.pyplot as plt
import random
from math import *
from random import *
from numba import jit
"Dynamics"
@jit(nopython=True)
def f(SP, Zavgi, Zreli, alf, dt, n):
"Time"
counter = 0;
Zavg = np.array([Zavgi]);
Zrelm = np.array([Zreli]);
T_plot = np.array([0]);
for i in range(1, iter + 1):
t = i * dt;
np.append(T_plot,t)
Z = [];
Up = [];
Down = [];
c, s = cos(t), sin(t);
c1, s1 = cos(t - dt), sin(t - dt);
for j in range(n - 1):
collchk = ((c * (SP[j][0]) + s * (SP[j][1])) - (c * (SP[j + 1][0]) + s * (SP[j + 1][1]))) * (
c1 * (SP[j][0]) + s1 * (SP[j][1]) - (c1 * (SP[j + 1][0]) + s1 * (SP[j + 1][1])));
prel = ((c * (SP[j][1]) - s * (SP[j][0])) - (c * (SP[j + 1][1]) - s * (SP[j + 1][0]))) / 2;
rcoeff = 1 / (1 + (prel * alf) ** 2);
rand_value = random();
if collchk < 0:
SP[j], SP[j + 1] = SP[j + 1], SP[j];
if rcoeff > rand_value:
counter = counter + 1
SP[j][2], SP[j + 1][2] = SP[j + 1][2], SP[j][2];
if SP[j][2] == 1:
Up.append(c * (SP[j][0]) + s * (SP[j][1]))
else:
Down.append(c * (SP[j][0]) + s * (SP[j][1]))
Z.append(c * (SP[j][0]) + s * (SP[j][1]))
Zrel = np.sum(np.array(Up)) / len(Up) - np.sum(np.array(Down)) / len(Down);
Zrelm = np.append(Zrelm, Zrel)
Zm = np.sum(np.array(Z)) / len(Z)
Zavg = np.append(Zavg, Zm)
return Zavg, Zrelm, counter, T_plot
if __name__ == '__main__':
n = 1000
mu = np.random.uniform(0, 1, n)
r = [sqrt(-2 * log(1 - i)) for i in mu]
eta = np.random.uniform(0, 1, n)
theta = 2 * pi * eta;
cuz = [cos(i) for i in theta]
suz = [sin(i) for i in theta]
Zinitial = [a * b for a, b in zip(r, cuz)];
Pinitial = [a * b for a, b in zip(r, suz)];
iter = 10 ** (4);
dt = 1 / (2 * n);
alf = sqrt(n);
SP = np.array(sorted(np.array([ np.array([i,j,choice([1,0])]) for i, j in zip(Zinitial, Pinitial)]),
key=lambda x: x[0]))
Upi = [];
Downi = [];
count_plot = [];
for j in range(len(SP)):
if SP[j][2] == 1:
Upi.append(SP[j][0])
else:
Downi.append(SP[j][0])
Zavgi = np.sum(Zinitial) / len(Zinitial)
Zreli = np.sum(Upi) / len(Upi) - np.sum(Downi) / len(Downi)
Zavg, Zrelm, counter, T_plot = f(SP, Zavgi, Zreli, alf, dt, n)
print("rate= ", counter/(n*(iter*dt)))
输出:
比率= 0.00814
可以看出，即使参数相同，三种情况下的速率也不同。
带有 Numba 的速率相差 10 倍。
由于参数相同，我希望三种不同的方法给出非常接近的比率。
为什么会这样？
编辑：
我通过简单地 运行 缩短了代码的时间，缩短了很多时间。因此，我还删除了在如此短的时间范围内看起来非常奇怪的图表。我想强调的是，这个问题似乎有两方面但相关：多次运行的速率因顺序而异，并且阵列情况（有和没有 Numba）的图表显得更加紧张（这在长时间运行时变得非常突出）比 class 案例。
TL;DR: 一个问题来自天真的行 SP[j], SP[j + 1] = SP[j + 1], SP[j] 由于 SP 的类型不同，程序之间的行为也不同。另一个来自在 Numba 中表现不同的随机数生成（如果不考虑再现性就很难跟踪）。
调试
在深入研究问题之前，首先要做的是更改代码以获得确定性结果以便能够重现调试器中的错误，然后跟踪与此的差异。然后，您需要调整代码，以便它可以运行快速并且易于调试 因为程序通常要执行几次才能找到错误。然后，您可以 log 影响 counter 的变量值（即导致不同结果的变量）。然后，您可以检查差异并 找到 程序之间 counter 不同的步骤。最后，你可以一步步回溯到bug的源头。这是一种非常有效地调试此类程序的通用方法。请注意，如果您可以使用 反向调试器 ，这可能会更有效（但事实证明我还没有找到在我的机器上适用于 Python 的调试器）。
获得确定性可重现结果的关键是使用以下代码：
np.random.seed(0)
seed(0)
默认情况下，Numba 的随机数生成器与 Numpy 的随机数生成器不同步，因此需要将种子放入 Numba 函数中，以便实际影响 Numba 的随机数生成器的种子（而不是来自 CPython 的 Numpy 之一）。此外，Numba 的 random() 也不与 CPython 同步，但使用相同的种子是不够的，因为 Numba 似乎并不在内部使用相同的算法。使用 np.random.rand() 而不是 random() 修复了重现性问题。
引擎盖下
以下是了解问题根源的回溯步骤：
counter 在特定步骤不同（当 (i,j)=(2,20) 用于 n=100 和 iter=3）collchk和prel在第二个程序中设置为0，而在第一个SP[j] 和 SP[j+1] 在第二个程序中相等，但在第一个程序中不相等 SP[j], SP[j + 1] = SP[j + 1], SP[j]导致值只在第二个程序中相等
此行在前两个程序中的行为不同，因为前者处理一维对象数组，而后者处理原生值的二维数组。事实上，在前者中，对象被正确交换，但后来 在 Numpy 视图上运行 并且视图的 隐式副本不复制数组内容 导致覆盖。以下是它在内部的运作方式：
SP[j], SP[j + 1] = SP[j + 1], SP[j] 被隐含地翻译成等价于：
tmp1 = SP[j+1,:]
tmp2 = SP[j,:]
SP[j,:] = tmp1
SP[j+1,:] = tmp2
问题在于 tmp1 和 tmp2 是两个 Numpy 视图的引用：SP[j,:] 和 SP[j+1,:]。因此，SP[j] = tmp1 导致 tmp2 的视图内容也被覆盖。因此，最后，该行被复制，导致一个偷偷摸摸的错误。
解决此问题的一种方法是在交换的行上显式调用 np.copy。 Swap two values in a numpy array. 还提供了另一种有趣的解决方案（请注意，目前已验证的第二个答案是导致您的案例出现虚假结果的答案，并且评论指出这不是一个安全的解决方案）。使用 SP[[j, j+1]] = SP[[j+1, j]] 解决了这个问题，因为 SP[[j+1, j]] 创建了一个临时数组（这不是 SP 的视图）。请注意，Numba 不支持此解决方案。 Numba 确实支持复制，但创建临时数组的成本很高。对于 Numba 代码，您需要使用以下代码：
for k in range(SP.shape[1]):
SP[j, k], SP[j + 1, k] = SP[j + 1, k], SP[j, k]
最后，请注意，使用 np.random.rand 进行上述修复足以在 Numba 中获得相同的结果。
总而言之，我的代码具有完全相同的参数，但使用了三种不同的方法：
正在为我的粒子 SP 列表定义一个 class，其中包含三个 instances/attributes。 （大约 1 分钟）
正在将 SP 定义为数组。 （大约 1 分钟）
将 SP 定义为数组并使用 Numba（大约 5 秒）
考虑第一种情况：
n=1000
mu=np.random.uniform(0,1,n)
r=[sqrt(-2*log(1-i)) for i in mu]
eta=np.random.uniform(0,1,n)
theta=2*pi*eta;
cuz=[cos(i) for i in theta]
suz=[sin(i) for i in theta]
Zinitial=[a*b for a,b in zip(r,cuz)];
Pinitial=[a*b for a,b in zip(r,suz)];
class Particle:
def __init__(self, pos, mom, spin):
self.pos = pos
self.mom = mom
self.spin = spin
SP = sorted([Particle(pos = i, mom = j, spin = choice([1, 0])) for i,j in zip(Zinitial,Pinitial)],key=lambda x:x.pos)
Upi=[];
Downi=[];
count_plot=[];
for j in range(len(SP)):
if SP[j].spin == 1:
Upi.append(SP[j].pos)
else:
Downi.append(SP[j].pos)
Zavgi=sum(Zinitial)/len(Zinitial)
Zreli=sum(Upi)/len(Upi)-sum(Downi)/len(Downi)
"Observables"
Zavg=[Zavgi];
Zrelm=[Zreli];
T_plot=[0];
"Time"
iter=10**(4);
dt=1/(2*n);
alf=sqrt(n);
"Dynamics"
counter=0;
sum1,sum2=0,0;
for i in range(1,iter+1):
t=i*dt;
T_plot.append(t)
Z=[];
Up=[];
Down=[];
c,s=cos(t),sin(t);
c1,s1=cos(t-dt),sin(t-dt);
for j in range(n-1):
collchk=((c*(SP[j].pos)+s*(SP[j].mom))-(c*(SP[j+1].pos)+s*(SP[j+1].mom)))*(c1*(SP[j].pos)+s1*(SP[j].mom)-(c1*(SP[j+1].pos)+s1*(SP[j+1].mom)));
prel=((c*(SP[j].mom)-s*(SP[j].pos))-(c*(SP[j+1].mom)-s*(SP[j+1].pos)))/2;
rcoeff=1/(1+(prel*alf)**2);
rand_value=random();
if collchk<0:
SP[j], SP[j+1]=SP[j+1],SP[j];
if rcoeff>rand_value:
counter=counter+1
SP[j].spin,SP[j+1].spin=SP[j+1].spin,SP[j].spin;
if SP[j].spin == 1:
Up.append(c*(SP[j].pos)+s*(SP[j].mom))
else:
Down.append(c*(SP[j].pos)+s*(SP[j].mom))
Z.append(c*(SP[j].pos)+s*(SP[j].mom))
Zrel=sum(Up[0:])/len(Up) - sum(Down[0:])/len(Down);
Zrelm.append(Zrel)
Zm=sum(Z)/len(Z)
Zavg.append(Zm)
print("Rate of collision per particle = ",counter/(n*(dt*iter)))
输出是： 每个粒子的碰撞率 = 0.0722
第二种情况：
n=1000
mu=np.random.uniform(0,1,n)
r=np.array([sqrt(-2*log(1-i)) for i in mu])
eta=np.random.uniform(0,1,n)
theta=2*pi*eta;
cuz=np.array([cos(i) for i in theta])
suz=np.array([sin(i) for i in theta])
Zinitial=np.multiply(r,cuz);
Pinitial=np.multiply(r,suz);
SP = np.array(sorted(np.array([ np.array([i,j,choice([1,0])]) for i, j in zip(Zinitial, Pinitial)]),
key=lambda x: x[0]))
Upi=[];
Downi=[];
count_plot=[];
for j in range(len(SP)):
if SP[j][2] == 1:
Upi.append(SP[j][0])
else:
Downi.append(SP[j][0])
Zavgi=sum(Zinitial)/len(Zinitial)
Zreli=sum(Upi)/len(Upi)-sum(Downi)/len(Downi)
"Observables"
Zavg=[Zavgi];
Zrelm=[Zreli];
T_plot=[0];
"Time"
iter=10**(4);
dt=1/(2*n);
alf=sqrt(n);
"Dynamics"
counter=0;
sum1,sum2=0,0;
for i in range(1,iter+1):
t=i*dt;
T_plot.append(t)
Z=[];
Up=[];
Down=[];
c,s=cos(t),sin(t);
c1,s1=cos(t-dt),sin(t-dt);
for j in range(n-1):
collchk=((c*(SP[j][0])+s*(SP[j][1]))-(c*(SP[j+1][0])+s*(SP[j+1][1])))*(c1*(SP[j][0])+s1*(SP[j][1])-(c1*(SP[j+1][0])+s1*(SP[j+1][1])));
prel=((c*(SP[j][1])-s*(SP[j][0]))-(c*(SP[j+1][1])-s*(SP[j+1][0])))/2;
rcoeff=1/(1+(prel*alf)**2);
rand_value=random();
if collchk<0:
SP[j], SP[j+1]=SP[j+1],SP[j];
if rcoeff>rand_value:
counter=counter+1
SP[j][2],SP[j+1][2]=SP[j+1][2],SP[j][2];
if SP[j][2] == 1:
Up.append(c*(SP[j][0])+s*(SP[j][1]))
else:
Down.append(c*(SP[j][0])+s*(SP[j][1]))
Z.append(c*(SP[j][0])+s*(SP[j][1]))
Zrel=sum(Up[0:])/len(Up) - sum(Down[0:])/len(Down);
Zrelm.append(Zrel)
Zm=sum(Z)/len(Z)
Zavg.append(Zm)
print("Rate of collision per particle = ",counter/(n*(dt*iter)))
输出是：
每个粒子的碰撞率 = 0.0134
以及最快的 Numba 案例；
import numpy as np
import matplotlib.pyplot as plt
import random
from math import *
from random import *
from numba import jit
"Dynamics"
@jit(nopython=True)
def f(SP, Zavgi, Zreli, alf, dt, n):
"Time"
counter = 0;
Zavg = np.array([Zavgi]);
Zrelm = np.array([Zreli]);
T_plot = np.array([0]);
for i in range(1, iter + 1):
t = i * dt;
np.append(T_plot,t)
Z = [];
Up = [];
Down = [];
c, s = cos(t), sin(t);
c1, s1 = cos(t - dt), sin(t - dt);
for j in range(n - 1):
collchk = ((c * (SP[j][0]) + s * (SP[j][1])) - (c * (SP[j + 1][0]) + s * (SP[j + 1][1]))) * (
c1 * (SP[j][0]) + s1 * (SP[j][1]) - (c1 * (SP[j + 1][0]) + s1 * (SP[j + 1][1])));
prel = ((c * (SP[j][1]) - s * (SP[j][0])) - (c * (SP[j + 1][1]) - s * (SP[j + 1][0]))) / 2;
rcoeff = 1 / (1 + (prel * alf) ** 2);
rand_value = random();
if collchk < 0:
SP[j], SP[j + 1] = SP[j + 1], SP[j];
if rcoeff > rand_value:
counter = counter + 1
SP[j][2], SP[j + 1][2] = SP[j + 1][2], SP[j][2];
if SP[j][2] == 1:
Up.append(c * (SP[j][0]) + s * (SP[j][1]))
else:
Down.append(c * (SP[j][0]) + s * (SP[j][1]))
Z.append(c * (SP[j][0]) + s * (SP[j][1]))
Zrel = np.sum(np.array(Up)) / len(Up) - np.sum(np.array(Down)) / len(Down);
Zrelm = np.append(Zrelm, Zrel)
Zm = np.sum(np.array(Z)) / len(Z)
Zavg = np.append(Zavg, Zm)
return Zavg, Zrelm, counter, T_plot
if __name__ == '__main__':
n = 1000
mu = np.random.uniform(0, 1, n)
r = [sqrt(-2 * log(1 - i)) for i in mu]
eta = np.random.uniform(0, 1, n)
theta = 2 * pi * eta;
cuz = [cos(i) for i in theta]
suz = [sin(i) for i in theta]
Zinitial = [a * b for a, b in zip(r, cuz)];
Pinitial = [a * b for a, b in zip(r, suz)];
iter = 10 ** (4);
dt = 1 / (2 * n);
alf = sqrt(n);
SP = np.array(sorted(np.array([ np.array([i,j,choice([1,0])]) for i, j in zip(Zinitial, Pinitial)]),
key=lambda x: x[0]))
Upi = [];
Downi = [];
count_plot = [];
for j in range(len(SP)):
if SP[j][2] == 1:
Upi.append(SP[j][0])
else:
Downi.append(SP[j][0])
Zavgi = np.sum(Zinitial) / len(Zinitial)
Zreli = np.sum(Upi) / len(Upi) - np.sum(Downi) / len(Downi)
Zavg, Zrelm, counter, T_plot = f(SP, Zavgi, Zreli, alf, dt, n)
print("rate= ", counter/(n*(iter*dt)))
输出:
比率= 0.00814
可以看出，即使参数相同，三种情况下的速率也不同。 带有 Numba 的速率相差 10 倍。
由于参数相同，我希望三种不同的方法给出非常接近的比率。
为什么会这样？
编辑：
我通过简单地 运行 缩短了代码的时间，缩短了很多时间。因此，我还删除了在如此短的时间范围内看起来非常奇怪的图表。我想强调的是，这个问题似乎有两方面但相关：多次运行的速率因顺序而异，并且阵列情况（有和没有 Numba）的图表显得更加紧张（这在长时间运行时变得非常突出）比 class 案例。
TL;DR: 一个问题来自天真的行 SP[j], SP[j + 1] = SP[j + 1], SP[j] 由于 SP 的类型不同，程序之间的行为也不同。另一个来自在 Numba 中表现不同的随机数生成（如果不考虑再现性就很难跟踪）。
调试
在深入研究问题之前，首先要做的是更改代码以获得确定性结果以便能够重现调试器中的错误，然后跟踪与此的差异。然后，您需要调整代码，以便它可以运行快速并且易于调试 因为程序通常要执行几次才能找到错误。然后，您可以 log 影响 counter 的变量值（即导致不同结果的变量）。然后，您可以检查差异并 找到 程序之间 counter 不同的步骤。最后，你可以一步步回溯到bug的源头。这是一种非常有效地调试此类程序的通用方法。请注意，如果您可以使用 反向调试器 ，这可能会更有效（但事实证明我还没有找到在我的机器上适用于 Python 的调试器）。
获得确定性可重现结果的关键是使用以下代码：
np.random.seed(0)
seed(0)
默认情况下，Numba 的随机数生成器与 Numpy 的随机数生成器不同步，因此需要将种子放入 Numba 函数中，以便实际影响 Numba 的随机数生成器的种子（而不是来自 CPython 的 Numpy 之一）。此外，Numba 的 random() 也不与 CPython 同步，但使用相同的种子是不够的，因为 Numba 似乎并不在内部使用相同的算法。使用 np.random.rand() 而不是 random() 修复了重现性问题。
引擎盖下
以下是了解问题根源的回溯步骤：
counter在特定步骤不同（当(i,j)=(2,20)用于n=100和iter=3）collchk和prel在第二个程序中设置为0，而在第一个 中它们都不同为0
SP[j]和SP[j+1]在第二个程序中相等，但在第一个程序中不相等SP[j], SP[j + 1] = SP[j + 1], SP[j]导致值只在第二个程序中相等
此行在前两个程序中的行为不同，因为前者处理一维对象数组，而后者处理原生值的二维数组。事实上，在前者中，对象被正确交换，但后来 在 Numpy 视图上运行 并且视图的 隐式副本不复制数组内容 导致覆盖。以下是它在内部的运作方式：
SP[j], SP[j + 1] = SP[j + 1], SP[j] 被隐含地翻译成等价于：
tmp1 = SP[j+1,:]
tmp2 = SP[j,:]
SP[j,:] = tmp1
SP[j+1,:] = tmp2
问题在于 tmp1 和 tmp2 是两个 Numpy 视图的引用：SP[j,:] 和 SP[j+1,:]。因此，SP[j] = tmp1 导致 tmp2 的视图内容也被覆盖。因此，最后，该行被复制，导致一个偷偷摸摸的错误。
解决此问题的一种方法是在交换的行上显式调用 np.copy。 Swap two values in a numpy array. 还提供了另一种有趣的解决方案（请注意，目前已验证的第二个答案是导致您的案例出现虚假结果的答案，并且评论指出这不是一个安全的解决方案）。使用 SP[[j, j+1]] = SP[[j+1, j]] 解决了这个问题，因为 SP[[j+1, j]] 创建了一个临时数组（这不是 SP 的视图）。请注意，Numba 不支持此解决方案。 Numba 确实支持复制，但创建临时数组的成本很高。对于 Numba 代码，您需要使用以下代码：
for k in range(SP.shape[1]):
SP[j, k], SP[j + 1, k] = SP[j + 1, k], SP[j, k]
最后，请注意，使用 np.random.rand 进行上述修复足以在 Numba 中获得相同的结果。