title: Why deploying Numba and using arrays instead of class gives different results for the same parameters in my program?
body: <p>So, to summarize I have the exact same parameters for my code but with three different methods:</p>
  <ol>
  <li><p>Defining a class for my list of particles SP with three instances/attributes. (around 1 mins)</p>
  </li>
  <li><p>Defining SP as an array. (around 1 mins)</p>
  </li>
  <li><p>Defining SP as an array and using Numba (around 5 seconds)</p>
  </li>
  </ol>
  <p>Consider the first case:</p>
  <pre><code>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)
  &quot;Observables&quot;
  Zavg=[Zavgi];
  Zrelm=[Zreli];
  T_plot=[0];
  
  &quot;Time&quot;
  
  iter=10**(4);
  dt=1/(2*n);
  alf=sqrt(n);
  
  &quot;Dynamics&quot;
  
  
  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&lt;0:
                 
                
                 SP[j], SP[j+1]=SP[j+1],SP[j];
                 
                
                 if rcoeff&gt;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(&quot;Rate of collision per particle = &quot;,counter/(n*(dt*iter)))
  
   
  </code></pre>
  <p>The <em>OUTPUT</em> is:
  Rate of collision per particle =  0.0722</p>
  <p>The second case:</p>
  <pre><code>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)
  &quot;Observables&quot;
  Zavg=[Zavgi];
  Zrelm=[Zreli];
  T_plot=[0];
  
  &quot;Time&quot;
  
  iter=10**(4);
  dt=1/(2*n);
  alf=sqrt(n);
  
  &quot;Dynamics&quot;
  
  
  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&lt;0:
                 
                
                 SP[j], SP[j+1]=SP[j+1],SP[j];
                 
                
                 if rcoeff&gt;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(&quot;Rate of collision per particle = &quot;,counter/(n*(dt*iter)))
  
    
  </code></pre>
  <p>The <em>OUTPUT</em> is :</p>
  <p>Rate of collision per particle =  0.0134</p>
  <p>And the quickest Numba case;</p>
  <pre><code>import numpy as np
  import matplotlib.pyplot as plt
  import random
  from math import *
  from random import *
  from numba import jit
  
  
  &quot;Dynamics&quot;
  
  @jit(nopython=True)
  def f(SP, Zavgi, Zreli, alf, dt, n):
      &quot;Time&quot;
      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 &lt; 0:
  
                  SP[j], SP[j + 1] = SP[j + 1], SP[j];
  
                  if rcoeff &gt; 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(&quot;rate= &quot;, counter/(n*(iter*dt)))
      
  </code></pre>
  <p>The <em>OUTPUT</em>:</p>
  <p>rate=  0.00814</p>
  <p>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.</p>
  <p>Since the parameters are the same I would expect the three different methods to give very close rates.</p>
  <p>Why is this happening?</p>
  <p><em><strong>EDIT:</strong></em></p>
  <p>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.</p>