首页 / 技术问答 / Python问答列表 / 问答详情

Python 为什么部署 Numba 并使用数组而不是类会在我的程序中为相同的参数提供不同的结果

总之，我的代码有完全相同的参数，但有三种不同的方法:

  1. 为具有三个实例/属性的粒子SP列表定义一个类.(约1分钟)

  2. 将SP定义为array.(约1分钟)

  3. 将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)))

OUTPUT是:

第二种情况:

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)))

OUTPUT是:

每个粒子的碰撞率=0.0134

最快的麻木病例；

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)))

OUTPUT人:

费率=0.00814

可以看出，即使参数相同，这三种情况下的速率也不同.

由于参数相同，我希望三种不同的方法可以给出非常接近的速率.

为什么会这样？

100

我通过简单地运行代码缩短了代码的时间.由于这个原因，我还删除了在如此短的时间尺度下看起来非常奇怪的图表.我想强调的是，这个问题似乎有两个方面，但又是相关的:许多次运行的速率相差一个数量级，而数组情况(有和没有Numba)的图形似乎比类情况更不稳定(这在长时间运行中变得非常突出).

推荐答案

TL;DR:一个问题来自天真的第SP[j], SP[j + 1] = SP[j + 1], SP[j]行，由于SP的类型不同，程序之间的行为不同.另一种是随机数生成，在Numba中表现不同(如果不考虑再现性，很难跟踪).

调试

在深入研究该问题之前，首先要做的是更改代码，以获得deterministic results，以便能够在调试器中消除错误，然后使用它跟踪差异.然后，您需要adapt the code，这样它可以是run quickly，也可以是easy to debug，因为程序通常会在找到bug之前执行几次.然后，您可以log影响counter的变量值(这是导致不同结果的变量).然后，您可以判断差异和locate步骤，其中counter是程序之间的差异.最后，您可以一步一步地找到bug的来源.这是一种非常有效地调试此类程序的通用方法.请注意，如果您可以使用reverse debugger，那么它的效率可能会更高(但事实证明，我还没有找到一个可以在我的Python机器上工作的).

获得确定性可再现结果的关键是使用以下代码:

np.random.seed(0)
seed(0)

默认情况下，Numba的随机数生成器与Numpy的随机数生成器不同步，因此需要将种子放入Numba函数中，以便实际影响Numba随机数生成器的种子(而不是CPython的Numpy的种子).此外，Numba的random()也没有与CPython的random()同步，但使用相同的种子是不够的，因为Numba似乎在内部没有使用相同的算法.使用np.random.rand()而不是random()解决了再现性问题.

引擎盖下

以下是了解问题根源的回溯步骤:

  • counter在特定步骤上不同(当(i,j)=(2,20)用于n=100和iter=3时)
  • collchk和prel在第二个程序中设置为0，而在第一个程序中它们都不同于0
  • SP[j]和SP[j+1]在第二个程序中相等，但在第一个程序中不相等
  • 100 causes the value to be equal only in the second program

这一行在前两个程序中的行为不同，因为前者处理对象的1D数组，后者处理本机值的2Darray.实际上，在前者中，对象被正确交换，但在后者中，operates on Numpy views和implicit copy of the view does not copy the array content会导致覆盖.以下是其内部工作原理:

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的视图的内容也被覆盖.因此，最终，该行被复制，从而产生一个鬼祟祟的bug.

解决此问题的一个解决方案是显式调用交换行上的np.copy.Swap two values in a numpy array.还提供了另一个有趣的解决方案(请注意，第二个经过验证的答案目前是在您的 case 中导致虚假结果的答案，并且 comments 指出这不是一个安全的解决方案).使用SP[[j, j+1]] = SP[[j+1, j]]修复了这个问题，因为SP[[j+1, j]]创建了一个临时数组(它不是SP的视图).请注意，Nuba不支持此解决方案.Numba确实支持拷贝，但创建临时数组的成本很高.对于Nuba代码，需要使用以下代码:

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中获得相同的结果.

Python相关问答推荐

如何使用函数正确索引收件箱？

根据多列和一些条件创建新列

带有pandas的分区列上的过滤器的多个条件read_parquet

Python中的函数中是否有充分的理由接受float而不接受int？

如何在Deliveryter笔记本中从同步上下文正确地安排和等待Delivercio代码中的结果？

Class_weight参数不影响RandomForestClassifier不平衡数据集中的结果

将numpy数组存储在原始二进制文件中

Python会扔掉未使用的表情吗？

使用FASTCGI在IIS上运行Django频道

Python json.转储包含一些UTF-8字符的二元组，要么失败，要么转换它们.我希望编码字符按原样保留

Python中的嵌套Ruby哈希

如何在Windows上用Python提取名称中带有逗号的文件？

如何使用LangChain和AzureOpenAI在Python中解决AttribeHelp和BadPressMessage错误？

如何记录脚本输出

使用@ guardlasses. guardlass和注释的Python继承

计算组中唯一值的数量

Python键入协议默认值

调用decorator返回原始函数的输出

搜索按钮不工作，Python tkinter

交替字符串位置的正则表达式

实用课程推荐
相关教程推荐
休闲君评测 |  TOOLFK工具网 |  古诗文网 |  良许Linux教程网 |  Python之禅 |  红色石头 |  入党申请书 | 

© 2020 Hi LearnFk 意见反馈 免责声明 关于我们 桂ICP备11002319号-3 友链&广告位+QQ: 1963612630