python – Why deploying Numba and using arrays instead of class gives different results for the same parameters in my program?

So, to summarize I have the exact same parameters for my code but with three different methods:

  1. Defining a class for my list of particles SP with three instances/attributes. (around 1 min)

  2. Defining SP as an array. (around 1 min)

  3. Defining SP as an array and using Numba (around 5 seconds)

Consider the first case:

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
        Z=[];
        Up=[];
        Down=[];
        c,s=cos
        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)))

 

The OUTPUT is: Rate of collision per particle = 0.0722

The second case:

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
        Z=[];
        Up=[];
        Down=[];
        c,s=cos
        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)))

  

The OUTPUT is :

Rate of collision per particle = 0.0134

And the quickest Numba case;

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

The OUTPUT:

rate=0.00814

As can be seen, the rates are different for all three cases even when the parameters are the same. The one with the Numba has a rate different by a factor of 10.

Since the parameters are the same I would expect the three different methods to give very close rates.

Why is this happening?

EDIT:

I have shortened the time of the code by simply running it for a much smaller time. Due to this I have also removed the graphs which look super weird for such short time scales. I want to stress that the problem seems to be two fold but related: The rates differ by an order for many runs and the graphs for the array case (with and without Numba) appear more jittered (which becomes really prominent for long time runs) than the class case.

Leave a Comment