Why am I not getting the same result every time I run through my loop?

Question

I'm doing some signal processing, and I have a loop in Python in which I'm trying to optimize a result, ln(B), over a parameter alpha. When I ran it, I got an optimal alpha value of -0.8, but when I used it in my regular code it gave me a different result than the loop had given me. For clarity's sake, the regular code's result was ~4 and the loop's result for alpha = -0.8 was ~5.

I thought that was weird, so after double checking that the process was identical in each case I ran the loop code over one instance of alpha =-0.8 and got the same result as the regular code, 4, as opposed to the result the loop had given me at alpha = -0.8 when iterating over a range. I then iterated the code over the list [-0.8,-0.8,-0.8, -0.8, -0.8, -0.8], and found that for the first -0.8, I got 4, and then for each subsequent one I got 5.

Is there some reason a loop would give me a different answer or process for the first item than the rest of them? I've been staring at it for hours and I can't find my mistake. I thought maybe I wasn't resetting a variable, but the difference only happens between the first iteration and the rest, not for each iteration.

This is my loop:

alpharange = [-0.8,-0.8,-0.8, -0.8, -0.8, -0.8]

for alpha in alpharange:
  run += 1
  print("Run {0}, alpha = {1}".format(run, alpha))

  #find var_n
  denominator = 0 
  r = np.arange(1,16)
  for n in r:
    if n*fecho <= fmax: 
        denominator += (n*fecho)**(2-2*alpha)/(var_f[n*nf])**2
  var_n = numerator/denominator 
  print("var_n = ", var_n)

  #find lnB
  YB_ln = np.zeros(spec_H1.shape)
  B_t = np.arange(0,len(YB_ln[0,:]),1) #length of a row
  B_f = np.arange(0,len(YB_ln[:,0]),1) #length of a column 

  for nt in tzoomindx:
    for nf in fzoomindx:
      for n in r:
          if n*nf<realfreqs.size-1:
              f = realfreqs[n*nf]
              YB_ln[nf,nt] += (np.abs(2*np.real(spec_H1[n*nf,nt]*np.conj(spec_L1[n*nf,nt])) + np.abs(spec_H1[n*nf,nt])**2 + np.abs(spec_L1[n*nf,nt])**2 ))/(1/(var_n*(1/f**alpha)) + 1/var_f[n*nf]) - np.log(1 + (var_n*(1/f**alpha))/var_f[n*nf])    
  peak = np.max(YB_ln)
  peakcoord = np.where((YB_ln==peak))

  if peak > lnB_max:
    alpha_max = alpha
    lnB_max = peak
    time_max = t[peakcoord[1]]-tinterval
    freq_max = realfreqs[peakcoord[0]]
    var_n_max = var_n
    print("!!!!!NEW MAX:  alpha = {0} with peak {1} at time {2} and frequency {3}".format(alpha_max,lnB_max,time_max,freq_max))

  alphalist.append(alpha)
  lnBlist.append(peak)
  timelist.append(t[peakcoord[1]]-tinterval)
  freqlist.append(realfreqs[peakcoord[0]])
  print("Peak lnB = {0} at time {1} and frequency {2}. Time: {3} min {4} s".format(peak,t[peakcoord[1]]-tinterval,realfreqs[peakcoord[0]],(time.time()-start_time)//60,(time.time()-start_time)%60))

print("lnB is optimized at alpha = {0} with peak {1} at time {2} and frequency {3}, var_n = {4}".format(alpha_max,lnB_max,time_max,freq_max,var_n_max))
print("Run took {0} min, {1} s".format((time.time()-start_time)//60,(time.time()-start_time)%60)) 

and the output:

Run 1, alpha = -0.8
var_n =  (1.1668471858083481e-14+0j)

/usr/local/lib/python3.6/dist-packages/ipykernel_launcher.py:72: ComplexWarning: Casting complex values to real discards the imaginary part

!!!!!NEW MAX:  alpha = -0.8 with peak 4.115523906554817 at time [1.0625] and frequency [72.]
Peak lnB = 4.115523906554817 at time [1.0625] and frequency [72.]. Time: 0.0 min 0.3648514747619629 s
Run 2, alpha = -0.8
var_n =  (3.664163403845884e-14+0j)
!!!!!NEW MAX:  alpha = -0.8 with peak 5.330720524005124 at time [1.0625] and frequency [72.]
Peak lnB = 5.330720524005124 at time [1.0625] and frequency [72.]. Time: 0.0 min 0.702958345413208 s
Run 3, alpha = -0.8
var_n =  (3.664163403845884e-14+0j)
Peak lnB = 5.330720524005124 at time [1.0625] and frequency [72.]. Time: 0.0 min 1.0434083938598633 s
Run 4, alpha = -0.8
var_n =  (3.664163403845884e-14+0j)
Peak lnB = 5.330720524005124 at time [1.0625] and frequency [72.]. Time: 0.0 min 1.375929832458496 s
Run 5, alpha = -0.8
var_n =  (3.664163403845884e-14+0j)
Peak lnB = 5.330720524005124 at time [1.0625] and frequency [72.]. Time: 0.0 min 1.7248213291168213 s
Run 6, alpha = -0.8
var_n =  (3.664163403845884e-14+0j)
Peak lnB = 5.330720524005124 at time [1.0625] and frequency [72.]. Time: 0.0 min 2.0683481693267822 s
lnB is optimized at alpha = -0.8 with peak 5.330720524005124 at time [1.0625] and frequency [72.], var_n = (3.664163403845884e-14+0j)
Run took 0.0 min, 2.069751739501953 s

(I ignore that error because the complex values only have a real part).

It's clear to me that it's the variable var_n that's changing and that's what's throwing off my result. I'm just not sure why it's different only for the first iteration.

My full code with all the variable definitions is as follows, I just didn't want to clutter the one above:

start_time = time.time()

#find the energy over the leading seconds

dt = t[1] - t[0]

inspiral = np.where((t-tinterval >= -2) & (t-tinterval <= 0))
Einspiral_sum = np.sum(E[inspiral])*dt

print("Summed E_inspiral = ", Einspiral_sum)

#find the 'constant'

noise = np.where((t-tinterval <= -5) & (t-tinterval >= -50))
Enoise = np.sum(E[noise])/(E[noise].size)

print("E_noise = ", Enoise)

#finding the variance 

var_f = var_f_same

fecho = 72 
nf = np.where((realfreqs==fecho))
nf = int(nf[0])
fmax = 300

tzoomindx = np.where((t-tinterval>=0.5)&(t-tinterval<=1.5))
tzoomindx = np.array(tzoomindx[0])
fzoomindx = np.where((realfreqs>=63)&(realfreqs<=92))
fzoomindx=np.array(fzoomindx[0])

upperbound = 0
lowerbound = -2
numerator = (Einspiral_sum - (Enoise*upperbound-Enoise*lowerbound))

run = 0
lnB_max = -1e5
lnBlist = []
alphalist = []
timelist = []
freqlist = []

#alpharange = np.arange(-1.5,1,0.1)
#alpharange = np.array([-0.5,-0.4,-0.3,-0.2,-0.1,0,0.1,0.2,0.3,0.4])
alpharange = [-0.8,-0.8,-0.8, -0.8, -0.8, -0.8]

for alpha in alpharange:
  run += 1
  print("Run {0}, alpha = {1}".format(run, alpha))

  #find var_n
  denominator = 0 
  r = np.arange(1,16)
  for n in r:
    if n*fecho <= fmax: 
        denominator += (n*fecho)**(2-2*alpha)/(var_f[n*nf])**2
  var_n = numerator/denominator #for Einspiral = sum

  #find lnB
  YB_ln = np.zeros(spec_H1.shape)
  B_t = np.arange(0,len(YB_ln[0,:]),1) #length of a row
  B_f = np.arange(0,len(YB_ln[:,0]),1) #length of a column 

  for nt in tzoomindx:
    for nf in fzoomindx:
      for n in r:
          if n*nf<realfreqs.size-1:
              f = realfreqs[n*nf]
              YB_ln[nf,nt] += (np.abs(2*np.real(spec_H1[n*nf,nt]*np.conj(spec_L1[n*nf,nt])) + np.abs(spec_H1[n*nf,nt])**2 + np.abs(spec_L1[n*nf,nt])**2 ))/(1/(var_n*(1/f**alpha)) + 1/var_f[n*nf]) - np.log(1 + (var_n*(1/f**alpha))/var_f[n*nf])    
  peak = np.max(YB_ln)
  peakcoord = np.where((YB_ln==peak))

  if peak > lnB_max:
    alpha_max = alpha
    lnB_max = peak
    time_max = t[peakcoord[1]]-tinterval
    freq_max = realfreqs[peakcoord[0]]
    var_n_max = var_n
    print("!!!!!NEW MAX:  alpha = {0} with peak {1} at time {2} and frequency {3}".format(alpha_max,lnB_max,time_max,freq_max))

  alphalist.append(alpha)
  lnBlist.append(peak)
  timelist.append(t[peakcoord[1]]-tinterval)
  freqlist.append(realfreqs[peakcoord[0]])
  print("Peak lnB = {0} at time {1} and frequency {2}. Time: {3} min {4} s".format(peak,t[peakcoord[1]]-tinterval,realfreqs[peakcoord[0]],(time.time()-start_time)//60,(time.time()-start_time)%60))

print("lnB is optimized at alpha = {0} with peak {1} at time {2} and frequency {3}, var_n = {4}".format(alpha_max,lnB_max,time_max,freq_max,var_n_max))
print("Run took {0} min, {1} s".format((time.time()-start_time)//60,(time.time()-start_time)%60))                 

E is the energy of my signal, spec_H1/spec_L1 are the spectrograms of my signal, and all in all I am trying to find the Bayes factor of a signal in a certain time and frequency range.

Answer 1

Possibly, as you iterate over nf in fzoomindx, you change nf, which is used in calculation of denominator in denominator += (n*fecho)**(2-2*alpha)/(var_f[n*nf])**2 ...