From 0ed304aa1d04816bcd75447958ae3c7b6ef61f7b Mon Sep 17 00:00:00 2001
From: Rayne Liu <rl746@cornell.edu>
Date: Tue, 18 Aug 2020 22:42:09 -0400
Subject: [PATCH] Found and fixed a bug on phases, adjusted plot saving, added
 autocorrelation time

---
 code/RDGW150914_ptemcee1.py          | 47 ++++++++++++++++++----------
 code/RDGW150914_ptemcee2.py          | 35 ++++++++++++++-------
 code/RDGW150914_ptemcee3.py          | 41 ++++++++++++++----------
 code/RDGW150914_ptemcee4.py          | 37 +++++++++++++---------
 code/condor_submit_RdownPtemcee1.sub |  2 +-
 code/condor_submit_RdownPtemcee3.sub |  2 +-
 6 files changed, 103 insertions(+), 61 deletions(-)

diff --git a/code/RDGW150914_ptemcee1.py b/code/RDGW150914_ptemcee1.py
index 21c204f..ff3e6ad 100755
--- a/code/RDGW150914_ptemcee1.py
+++ b/code/RDGW150914_ptemcee1.py
@@ -2,8 +2,9 @@
 # coding: utf-8
 import warnings
 warnings.simplefilter("ignore") #Used to suppress the RuntimeWarnings. But if you are changing the code for TESTING, MAKE SURE TO COMMENT IT OUT.
+
 '''
-This script calculates the RDGW150914 constraints with ptemcee - specifically, the n=1 case both varying or not varying the fundamental frequencies. It produces the chain plot, corner plot, parameter constraints, and data plotted with the 1-sigma band. Since we are working specifically for the n=1 case, we also add in the corner plot the combined chain of alpha_0 and alpha_1, in order to demonstrate that the two are indistinguishable with each other.
+This script calculates the RDGW150914 constraints with ptemcee - specifically, the n=1 case varying the fundamental frequencies, tshift = 0. It produces the chain plot, corner plot, parameter constraints, and data plotted with the 1-sigma band. Since we are working specifically for the n=1 case, we also compare alpha_0 and alpha_1 in a histogram, in order to demonstrate that the two are indistinguishable with each other.
 
 --- R^a_{yne} L^i_u, 08/09/2020
 '''
@@ -41,11 +42,11 @@ tshift=0
 vary_fund = True
 
 #sampler parameters
-npoints=20000 
-nwalkers = 1000
+npoints=200
+nwalkers = 42
 ntemps=16
 ndim = int(4*(nmax+1))
-burnin = 4200  #How many points do you burn before doing the corner plot. You need to watch the convergence of the chain plot a bit.
+burnin = 42  #How many points do you burn before doing the corner plot. You need to watch the convergence of the chain plot a bit.
             #This is trivial but often forgotten: this cannot be more than npoints! Usually 1/5~1/4 npoints is what I observe.
 numbins = 42 #corner plot parameter - how many bins you want
 datacolor = '#105670' #'#4fa3a7'
@@ -139,15 +140,15 @@ def log_prior(theta):
     y_0 = theta[6]
     y_1 = theta[7]
 
-    if all([nmax == 1, 0 <= tshift <= 5, vary_fund == True, -1.2 <= a_0 <= 1.2, -1.5 <= a_1 <= 1.5, -9.0 <= b_0 <= 9.0, -10.0 <= b_1 <= 21.0, 0 <= x_0 <= 1.6, 0 <= x_1 <= 1.4, 0 <= y_0 <= 2*np.pi, 0 <= y_1 <= 2*np.pi]):        
+    if all([nmax == 1, 0 <= tshift <= 5, vary_fund == True, -1.2 <= a_0 <= 1.2, -1.5 <= a_1 <= 1.5, -4.0 <= b_0 <= 9.0, -4.0 <= b_1 <= 21.0, 0 <= x_0 <= 1.6, 0 <= x_1 <= 1.4, 0 <= y_0 <= 2*np.pi, 0 <= y_1 <= 2*np.pi]):        
         return 0.0
-    elif all([nmax == 1, tshift == 19, vary_fund == True, -10.0 <= a_0 <= 10.0, -10.0 <= a_1 <= 10.0, -9.0 <= b_0 <= 20.0, -9.0 <= b_1 <= 25.0, 0 <= x_0 <= 0.8, 0 <= x_1 <= 0.9, 0 <= y_0 <= 2*np.pi, 0 <= y_1 <= 2*np.pi]):
+    elif all([nmax == 1, tshift == 19, vary_fund == True, -10.0 <= a_0 <= 10.0, -10.0 <= a_1 <= 10.0, -20.0 <= b_0 <= 30.0, -25.0 <= b_1 <= 30.0, 0 <= x_0 <= 0.6, 0 <= x_1 <= 0.9, 0 <= y_0 <= 2*np.pi, 0 <= y_1 <= 2*np.pi]):
         return 0.0
     
-#PAY EXTRA ATTENTION TO THESE TWO CASES, SINCE WE SHIFTED y_0. THE CORNER PLOTS LABELS NEED TO BE TREATED WITH CARE.
-    elif all([nmax == 1, 0 <= tshift <= 5, vary_fund == False, -10.0 <= a_0 <= 10.0, -1.5 <= a_1 <= 1.5, -9.0 <= b_0 <= 9.0, -10.0 <= b_1 <= 16.0, 0 <= x_0 <= 2.2, 0 <= x_1 <= 2.5, 0 <= y_0-np.pi <= 2*np.pi, 0 <= y_1 <= 2*np.pi,]):
+#PAY EXTRA ATTENTION TO THESE TWO CASES, SINCE WE SHIFTED y's. THE CORNER PLOTS LABELS NEED TO BE TREATED WITH CARE.
+    elif all([nmax == 1, 0 <= tshift <= 5, vary_fund == False, -10.0 <= a_0 <= 10.0, -1.5 <= a_1 <= 1.5, -9.0 <= b_0 <= 9.0, -6.0 <= b_1 <= 20.0, 0 <= x_0 <= 2.4, 0 <= x_1 <= 2.5, 0 <= y_0-np.pi <= 2*np.pi, 0 <= y_1-np.pi <= 2*np.pi,]):
         return 0.0
-    elif all([nmax == 1, tshift == 19, vary_fund == False, -10.0 <= a_0 <= 10.0, -10.0 <= a_1 <= 10.0, -9.0 <= b_0 <= 9.0, -10.0 <= b_1 <= 32.0, 0 <= x_0 <= 0.6, 0 <= x_1 <= 1.0, 0 <= y_0-np.pi <= 2*np.pi, 0 <= y_1 <= 2*np.pi]):
+    elif all([nmax == 1, tshift == 19, vary_fund == False, -10.0 <= a_0 <= 10.0, -8.0 <= a_1 <= 8.0, -9.0 <= b_0 <= 9.0, -10.0 <= b_1 <= 42.0, 0 <= x_0 <= 0.6, 0 <= x_1 <= 0.7, 0 <= y_0-np.pi <= 2*np.pi, 0 <= y_1-np.pi <= 2*np.pi]):
         return 0.0
 
     return -np.inf
@@ -161,6 +162,7 @@ def log_likelihood(theta):
         return -np.inf
     return result
 
+
 # Logposterior distribution for the residuals case.
 # The evidence is just a normalization factor
 def log_probability(theta):
@@ -176,8 +178,8 @@ def log_probability(theta):
 
 #Fit with ptemcee
 #Set the number of cores of your processors
-pool = choose_pool(28)
-pool.size = 28
+pool = choose_pool(24)
+pool.size = 24
 vary_param = float(vary_fund)
 np.random.seed(42)
 pos = np.array([[random.uniform(-0.1,0.1), random.uniform(-0.1,0.1), 4.28313743e-01, random.uniform(2.5, 2.6) + (1-vary_param) * np.pi]])
@@ -194,6 +196,7 @@ sampler.run_mcmc(pos,npoints)
 
 
 
+
 #Define labels and start plotting
 paramlabels_a = [r'$\alpha_'+str(i)+'$' for i in range (nmax+1)]
 paramlabels_b = [r'$\beta_'+str(i)+'$' for i in range (nmax+1)]
@@ -211,12 +214,12 @@ if vary_fund == False:
 #Chain plot
 fig, axes = plt.subplots(4*(nmax+1), 1, sharex=True, figsize=(12, 9*(nmax+1)))
 for i in range(4*(nmax+1)):
-    axes[i].plot(sampler.chain[0,:, :, i].T, color="k", alpha=0.4)
+    axes[i].plot(sampler.chain[0,:, :, i].T, color="k", alpha=0.4, rasterized=True)
     axes[i].yaxis.set_major_locator(MaxNLocator(5))
     axes[i].set_ylabel(paramlabels[i])
 axes[-1].set_xlabel('Iterations')
 #plt.show()
-fig.savefig(rootpath+'/git/rdstackingproject/plotsmc/vary'+str(vary_fund)+'nmax='+str(nmax)+'_tshift='+str(tshift)+'_'+str(npoints)+'pt_chain.pdf', format = 'pdf')
+fig.savefig(rootpath+'/git/rdstackingproject/plotsmc/vary'+str(vary_fund)+'nmax='+str(nmax)+'_tshift='+str(tshift)+'_'+str(npoints)+'pt_chain.png', format = 'png', bbox_inches = 'tight')
 
 
 
@@ -224,12 +227,20 @@ fig.savefig(rootpath+'/git/rdstackingproject/plotsmc/vary'+str(vary_fund)+'nmax=
 samples = sampler.chain[0,:, burnin:, :].reshape((-1, ndim))
 #samples for corner plot
 samples_corn = samples if vary_fund == True else np.delete(samples, np.s_[0,2], 1)
+
 #print('Values with peak likelihood:')
 lglk = np.array([log_likelihood(samples[i]) for i in range(len(samples))])
 pk = samples[np.argmax(lglk)]
 #print('pk:')
 #print(pk)
 pk_corn = pk if vary_fund == True else np.delete(pk, [0,2])
+#y_0 range needs some messaging to make the plot. But in order to make the whole picture consistent, better change the range of y_1 too.
+if vary_fund == False:
+    samples_corn.T[-1] = samples_corn.T[-1] - np.pi 
+    samples_corn.T[-2] = samples_corn.T[-2] - np.pi #This indeed changes samples_corn itself
+    pk[-1] -= np.pi
+    pk[-2] -= np.pi
+
 #print('pkFalse:')
 #print(pk)
     
@@ -241,6 +252,7 @@ median = np.median(samples_corn, axis=0)
 
 figcorn = corner.corner(samples_corn, bins = numbins, hist_bin_factor = 5, color = datacolor, truths=pk_corn, truth_color = pkcolor, plot_contours = True, labels = paramlabels_corner, quantiles=(0.05, 0.16, 0.5, 0.84, 0.95), levels=[1-np.exp(-0.5), 1-np.exp(-1.64 ** 2/2)], show_titles=True)
 
+
 #Extract the axes in order to add more important line plots
 naxes = len(pk_corn)
 axes = np.array(figcorn.axes).reshape((naxes, naxes))
@@ -260,7 +272,7 @@ for yi in range(naxes):
 figcorn.savefig(rootpath+'/git/rdstackingproject/plotsmc/vary'+str(vary_fund)+'nmax='+str(nmax)+'_tshift='+str(tshift)+'_'+str(npoints)+'pt_corner.pdf', format = 'pdf')
 
 
-#Now export the covariance and correlation matrices, as well as the logevidence
+#Now export the covariance and correlation matrices, as well as the logevidence and autocorrelation time
 #The covariant matrix
 samplesT = samples_corn.T
 cov_matr=np.cov(samplesT)
@@ -279,6 +291,10 @@ logevidence = sampler.log_evidence_estimate(fburnin=0.2)
 df3 = pd.DataFrame(logevidence,index = [r'logevidence', r'error'], columns=[r'vary='+str(vary_fund)+', nmax='+str(nmax)+', tshift='+str(tshift)+', '+str(npoints)+'pts'])
 df3.to_hdf(rootpath+'/git/rdstackingproject/plotsmc/vary'+str(vary_fund)+'nmax='+str(nmax)+'_tshift='+str(tshift)+'_'+str(npoints)+'pt_logevidence.h5', key='df')
 
+#The autocorrelation time, this time I don't intend to delete alpha0 and beta0, just keep the nan's there since we don't really use this data, only look at it. It's good to keep it for sanity check's sake.
+atcrr = sampler.get_autocorr_time()[0]
+df4 = pd.DataFrame(atcrr,index = paramlabels, columns=[r'Autocorrelation time, vary='+str(vary_fund)+', nmax='+str(nmax)+', tshift='+str(tshift)+', '+str(npoints)+'pts'])
+df4.to_hdf(rootpath+'/git/rdstackingproject/plotsmc/vary'+str(vary_fund)+'nmax='+str(nmax)+'_tshift='+str(tshift)+'_'+str(npoints)+'pt_autocrrtime.h5', key='df')    
 
 
 #Now plot alpha_0 on top of alpha_1 - only on top, not stacked, and only valid for vary_fund == True
@@ -318,7 +334,6 @@ plt.xlim(timesrd[0], timesrd[0]+80)
 plt.xlabel("t")
 plt.ylabel("h")
 
-figband.savefig(rootpath+'/git/rdstackingproject/plotsmc/vary'+str(vary_fund)+'nmax='+str(nmax)+'_tshift='+str(tshift)+'_'+str(npoints)+'pt_band.pdf', format = 'pdf')
-
+figband.savefig(rootpath+'/git/rdstackingproject/plotsmc/vary'+str(vary_fund)+'nmax='+str(nmax)+'_tshift='+str(tshift)+'_'+str(npoints)+'pt_band.png', format = 'png', bbox_inches = 'tight')
 #import os
 #os.system('afplay /System/Library/Sounds/Submarine.aiff')
diff --git a/code/RDGW150914_ptemcee2.py b/code/RDGW150914_ptemcee2.py
index 384fbe9..546c153 100755
--- a/code/RDGW150914_ptemcee2.py
+++ b/code/RDGW150914_ptemcee2.py
@@ -2,8 +2,9 @@
 # coding: utf-8
 import warnings
 warnings.simplefilter("ignore") #Used to suppress the RuntimeWarnings. But if you are changing the code for TESTING, MAKE SURE TO COMMENT IT OUT.
+
 '''
-This script calculates the RDGW150914 constraints with ptemcee - specifically, the n=1 case both varying or not varying the fundamental frequencies. It produces the chain plot, corner plot, parameter constraints, and data plotted with the 1-sigma band. Since we are working specifically for the n=1 case, we also add in the corner plot the combined chain of alpha_0 and alpha_1, in order to demonstrate that the two are indistinguishable with each other.
+This script calculates the RDGW150914 constraints with ptemcee - specifically, the n=1 case varying the fundamental frequencies, tshift = 19. It produces the chain plot, corner plot, parameter constraints, and data plotted with the 1-sigma band. Since we are working specifically for the n=1 case, we also compare alpha_0 and alpha_1 in a histogram, in order to demonstrate that the two are indistinguishable with each other.
 
 --- R^a_{yne} L^i_u, 08/09/2020
 '''
@@ -139,15 +140,15 @@ def log_prior(theta):
     y_0 = theta[6]
     y_1 = theta[7]
 
-    if all([nmax == 1, 0 <= tshift <= 5, vary_fund == True, -1.2 <= a_0 <= 1.2, -1.5 <= a_1 <= 1.5, -9.0 <= b_0 <= 9.0, -10.0 <= b_1 <= 21.0, 0 <= x_0 <= 1.6, 0 <= x_1 <= 1.4, 0 <= y_0 <= 2*np.pi, 0 <= y_1 <= 2*np.pi]):        
+    if all([nmax == 1, 0 <= tshift <= 5, vary_fund == True, -1.2 <= a_0 <= 1.2, -1.5 <= a_1 <= 1.5, -4.0 <= b_0 <= 9.0, -4.0 <= b_1 <= 21.0, 0 <= x_0 <= 1.6, 0 <= x_1 <= 1.4, 0 <= y_0 <= 2*np.pi, 0 <= y_1 <= 2*np.pi]):        
         return 0.0
-    elif all([nmax == 1, tshift == 19, vary_fund == True, -10.0 <= a_0 <= 10.0, -10.0 <= a_1 <= 10.0, -9.0 <= b_0 <= 20.0, -9.0 <= b_1 <= 25.0, 0 <= x_0 <= 0.8, 0 <= x_1 <= 0.9, 0 <= y_0 <= 2*np.pi, 0 <= y_1 <= 2*np.pi]):
+    elif all([nmax == 1, tshift == 19, vary_fund == True, -10.0 <= a_0 <= 10.0, -10.0 <= a_1 <= 10.0, -20.0 <= b_0 <= 30.0, -25.0 <= b_1 <= 30.0, 0 <= x_0 <= 0.6, 0 <= x_1 <= 0.9, 0 <= y_0 <= 2*np.pi, 0 <= y_1 <= 2*np.pi]):
         return 0.0
     
-#PAY EXTRA ATTENTION TO THESE TWO CASES, SINCE WE SHIFTED y_0. THE CORNER PLOTS LABELS NEED TO BE TREATED WITH CARE.
-    elif all([nmax == 1, 0 <= tshift <= 5, vary_fund == False, -10.0 <= a_0 <= 10.0, -1.5 <= a_1 <= 1.5, -9.0 <= b_0 <= 9.0, -10.0 <= b_1 <= 16.0, 0 <= x_0 <= 2.2, 0 <= x_1 <= 2.5, 0 <= y_0-np.pi <= 2*np.pi, 0 <= y_1 <= 2*np.pi,]):
+#PAY EXTRA ATTENTION TO THESE TWO CASES, SINCE WE SHIFTED y's. THE CORNER PLOTS LABELS NEED TO BE TREATED WITH CARE.
+    elif all([nmax == 1, 0 <= tshift <= 5, vary_fund == False, -10.0 <= a_0 <= 10.0, -1.5 <= a_1 <= 1.5, -9.0 <= b_0 <= 9.0, -6.0 <= b_1 <= 20.0, 0 <= x_0 <= 2.4, 0 <= x_1 <= 2.5, 0 <= y_0-np.pi <= 2*np.pi, 0 <= y_1-np.pi <= 2*np.pi,]):
         return 0.0
-    elif all([nmax == 1, tshift == 19, vary_fund == False, -10.0 <= a_0 <= 10.0, -10.0 <= a_1 <= 10.0, -9.0 <= b_0 <= 9.0, -10.0 <= b_1 <= 32.0, 0 <= x_0 <= 0.6, 0 <= x_1 <= 1.0, 0 <= y_0-np.pi <= 2*np.pi, 0 <= y_1 <= 2*np.pi]):
+    elif all([nmax == 1, tshift == 19, vary_fund == False, -10.0 <= a_0 <= 10.0, -8.0 <= a_1 <= 8.0, -9.0 <= b_0 <= 9.0, -10.0 <= b_1 <= 42.0, 0 <= x_0 <= 0.6, 0 <= x_1 <= 0.7, 0 <= y_0-np.pi <= 2*np.pi, 0 <= y_1-np.pi <= 2*np.pi]):
         return 0.0
 
     return -np.inf
@@ -218,7 +219,7 @@ for i in range(4*(nmax+1)):
     axes[i].set_ylabel(paramlabels[i])
 axes[-1].set_xlabel('Iterations')
 #plt.show()
-fig.savefig(rootpath+'/git/rdstackingproject/plotsmc/vary'+str(vary_fund)+'nmax='+str(nmax)+'_tshift='+str(tshift)+'_'+str(npoints)+'pt_chain.pdf', format = 'pdf')
+fig.savefig(rootpath+'/git/rdstackingproject/plotsmc/vary'+str(vary_fund)+'nmax='+str(nmax)+'_tshift='+str(tshift)+'_'+str(npoints)+'pt_chain.png', format = 'png', bbox_inches = 'tight')
 
 
 
@@ -226,12 +227,20 @@ fig.savefig(rootpath+'/git/rdstackingproject/plotsmc/vary'+str(vary_fund)+'nmax=
 samples = sampler.chain[0,:, burnin:, :].reshape((-1, ndim))
 #samples for corner plot
 samples_corn = samples if vary_fund == True else np.delete(samples, np.s_[0,2], 1)
+
 #print('Values with peak likelihood:')
 lglk = np.array([log_likelihood(samples[i]) for i in range(len(samples))])
 pk = samples[np.argmax(lglk)]
 #print('pk:')
 #print(pk)
 pk_corn = pk if vary_fund == True else np.delete(pk, [0,2])
+#y_0 range needs some messaging to make the plot. But in order to make the whole picture consistent, better change the range of y_1 too.
+if vary_fund == False:
+    samples_corn.T[-1] = samples_corn.T[-1] - np.pi 
+    samples_corn.T[-2] = samples_corn.T[-2] - np.pi #This indeed changes samples_corn itself
+    pk[-1] -= np.pi
+    pk[-2] -= np.pi
+
 #print('pkFalse:')
 #print(pk)
     
@@ -243,6 +252,7 @@ median = np.median(samples_corn, axis=0)
 
 figcorn = corner.corner(samples_corn, bins = numbins, hist_bin_factor = 5, color = datacolor, truths=pk_corn, truth_color = pkcolor, plot_contours = True, labels = paramlabels_corner, quantiles=(0.05, 0.16, 0.5, 0.84, 0.95), levels=[1-np.exp(-0.5), 1-np.exp(-1.64 ** 2/2)], show_titles=True)
 
+
 #Extract the axes in order to add more important line plots
 naxes = len(pk_corn)
 axes = np.array(figcorn.axes).reshape((naxes, naxes))
@@ -281,12 +291,16 @@ logevidence = sampler.log_evidence_estimate(fburnin=0.2)
 df3 = pd.DataFrame(logevidence,index = [r'logevidence', r'error'], columns=[r'vary='+str(vary_fund)+', nmax='+str(nmax)+', tshift='+str(tshift)+', '+str(npoints)+'pts'])
 df3.to_hdf(rootpath+'/git/rdstackingproject/plotsmc/vary'+str(vary_fund)+'nmax='+str(nmax)+'_tshift='+str(tshift)+'_'+str(npoints)+'pt_logevidence.h5', key='df')
 
+#The autocorrelation time, this time I don't intend to delete alpha0 and beta0, just keep the nan's there since we don't really use this data, only look at it. It's good to keep it for sanity check's sake.
+atcrr = sampler.get_autocorr_time()[0]
+df4 = pd.DataFrame(atcrr,index = paramlabels, columns=[r'Autocorrelation time, vary='+str(vary_fund)+', nmax='+str(nmax)+', tshift='+str(tshift)+', '+str(npoints)+'pts'])
+df4.to_hdf(rootpath+'/git/rdstackingproject/plotsmc/vary'+str(vary_fund)+'nmax='+str(nmax)+'_tshift='+str(tshift)+'_'+str(npoints)+'pt_autocrrtime.h5', key='df')    
 
 
 #Now plot alpha_0 on top of alpha_1 - only on top, not stacked, and only valid for vary_fund == True
 if vary_fund == True:
-    samplea0 = samples.T[0]
-    samplea1 = samples.T[1]
+    samplea0 = samplesT[0]
+    samplea1 = samplesT[1]
     fighist1 = plt.figure(figsize = (8, 6))
     n0, bins0, patches0 = plt.hist(samplea0, bins = numbins * 6, alpha = 0.5, label = r'$\alpha_0$')
     n1, bins1, patches1 = plt.hist(samplea1, bins = numbins * 10, alpha = 0.5, label = r'$\alpha_1$')
@@ -320,7 +334,6 @@ plt.xlim(timesrd[0], timesrd[0]+80)
 plt.xlabel("t")
 plt.ylabel("h")
 
-figband.savefig(rootpath+'/git/rdstackingproject/plotsmc/vary'+str(vary_fund)+'nmax='+str(nmax)+'_tshift='+str(tshift)+'_'+str(npoints)+'pt_band.pdf', format = 'pdf')
-
+figband.savefig(rootpath+'/git/rdstackingproject/plotsmc/vary'+str(vary_fund)+'nmax='+str(nmax)+'_tshift='+str(tshift)+'_'+str(npoints)+'pt_band.png', format = 'png', bbox_inches = 'tight')
 #import os
 #os.system('afplay /System/Library/Sounds/Submarine.aiff')
diff --git a/code/RDGW150914_ptemcee3.py b/code/RDGW150914_ptemcee3.py
index 38b61c8..1e8c2d9 100755
--- a/code/RDGW150914_ptemcee3.py
+++ b/code/RDGW150914_ptemcee3.py
@@ -4,7 +4,7 @@ import warnings
 warnings.simplefilter("ignore") #Used to suppress the RuntimeWarnings. But if you are changing the code for TESTING, MAKE SURE TO COMMENT IT OUT.
 
 '''
-This script calculates the RDGW150914 constraints with ptemcee - specifically, the n=1 case both varying or not varying the fundamental frequencies. It produces the chain plot, corner plot, parameter constraints, and data plotted with the 1-sigma band. Since we are working specifically for the n=1 case, we also add in the corner plot the combined chain of alpha_0 and alpha_1, in order to demonstrate that the two are indistinguishable with each other.
+This script calculates the RDGW150914 constraints with ptemcee - specifically, the n=1 case not varying the fundamental frequencies, tshift = 0. It produces the chain plot, corner plot, parameter constraints, and data plotted with the 1-sigma band. Since we are working specifically for the n=1 case, we also compare alpha_0 and alpha_1 in a histogram, in order to demonstrate that the two are indistinguishable with each other.
 
 --- R^a_{yne} L^i_u, 08/09/2020
 '''
@@ -42,11 +42,11 @@ tshift=0
 vary_fund = False
 
 #sampler parameters
-npoints=20000 
-nwalkers = 1000
+npoints=200 
+nwalkers = 42
 ntemps=16
 ndim = int(4*(nmax+1))
-burnin = 4200  #How many points do you burn before doing the corner plot. You need to watch the convergence of the chain plot a bit.
+burnin = 42  #How many points do you burn before doing the corner plot. You need to watch the convergence of the chain plot a bit.
             #This is trivial but often forgotten: this cannot be more than npoints! Usually 1/5~1/4 npoints is what I observe.
 numbins = 42 #corner plot parameter - how many bins you want
 datacolor = '#105670' #'#4fa3a7'
@@ -140,15 +140,15 @@ def log_prior(theta):
     y_0 = theta[6]
     y_1 = theta[7]
 
-    if all([nmax == 1, 0 <= tshift <= 5, vary_fund == True, -1.2 <= a_0 <= 1.2, -1.5 <= a_1 <= 1.5, -9.0 <= b_0 <= 9.0, -10.0 <= b_1 <= 21.0, 0 <= x_0 <= 1.6, 0 <= x_1 <= 1.4, 0 <= y_0 <= 2*np.pi, 0 <= y_1 <= 2*np.pi]):        
+    if all([nmax == 1, 0 <= tshift <= 5, vary_fund == True, -1.2 <= a_0 <= 1.2, -1.5 <= a_1 <= 1.5, -4.0 <= b_0 <= 9.0, -4.0 <= b_1 <= 21.0, 0 <= x_0 <= 1.6, 0 <= x_1 <= 1.4, 0 <= y_0 <= 2*np.pi, 0 <= y_1 <= 2*np.pi]):        
         return 0.0
-    elif all([nmax == 1, tshift == 19, vary_fund == True, -10.0 <= a_0 <= 10.0, -10.0 <= a_1 <= 10.0, -9.0 <= b_0 <= 20.0, -9.0 <= b_1 <= 25.0, 0 <= x_0 <= 0.8, 0 <= x_1 <= 0.9, 0 <= y_0 <= 2*np.pi, 0 <= y_1 <= 2*np.pi]):
+    elif all([nmax == 1, tshift == 19, vary_fund == True, -10.0 <= a_0 <= 10.0, -10.0 <= a_1 <= 10.0, -20.0 <= b_0 <= 30.0, -25.0 <= b_1 <= 30.0, 0 <= x_0 <= 0.6, 0 <= x_1 <= 0.9, 0 <= y_0 <= 2*np.pi, 0 <= y_1 <= 2*np.pi]):
         return 0.0
     
-#PAY EXTRA ATTENTION TO THESE TWO CASES, SINCE WE SHIFTED y_0. THE CORNER PLOTS LABELS NEED TO BE TREATED WITH CARE.
-    elif all([nmax == 1, 0 <= tshift <= 5, vary_fund == False, -10.0 <= a_0 <= 10.0, -1.5 <= a_1 <= 1.5, -9.0 <= b_0 <= 9.0, -10.0 <= b_1 <= 16.0, 0 <= x_0 <= 2.2, 0 <= x_1 <= 2.5, 0 <= y_0-np.pi <= 2*np.pi, 0 <= y_1 <= 2*np.pi,]):
+#PAY EXTRA ATTENTION TO THESE TWO CASES, SINCE WE SHIFTED y's. THE CORNER PLOTS LABELS NEED TO BE TREATED WITH CARE.
+    elif all([nmax == 1, 0 <= tshift <= 5, vary_fund == False, -10.0 <= a_0 <= 10.0, -1.5 <= a_1 <= 1.5, -9.0 <= b_0 <= 9.0, -6.0 <= b_1 <= 20.0, 0 <= x_0 <= 2.4, 0 <= x_1 <= 2.5, 0 <= y_0-np.pi <= 2*np.pi, 0 <= y_1-np.pi <= 2*np.pi,]):
         return 0.0
-    elif all([nmax == 1, tshift == 19, vary_fund == False, -10.0 <= a_0 <= 10.0, -10.0 <= a_1 <= 10.0, -9.0 <= b_0 <= 9.0, -10.0 <= b_1 <= 32.0, 0 <= x_0 <= 0.6, 0 <= x_1 <= 1.0, 0 <= y_0-np.pi <= 2*np.pi, 0 <= y_1 <= 2*np.pi]):
+    elif all([nmax == 1, tshift == 19, vary_fund == False, -10.0 <= a_0 <= 10.0, -8.0 <= a_1 <= 8.0, -9.0 <= b_0 <= 9.0, -10.0 <= b_1 <= 42.0, 0 <= x_0 <= 0.6, 0 <= x_1 <= 0.7, 0 <= y_0-np.pi <= 2*np.pi, 0 <= y_1-np.pi <= 2*np.pi]):
         return 0.0
 
     return -np.inf
@@ -178,8 +178,8 @@ def log_probability(theta):
 
 #Fit with ptemcee
 #Set the number of cores of your processors
-pool = choose_pool(28)
-pool.size = 28
+pool = choose_pool(24)
+pool.size = 24
 vary_param = float(vary_fund)
 np.random.seed(42)
 pos = np.array([[random.uniform(-0.1,0.1), random.uniform(-0.1,0.1), 4.28313743e-01, random.uniform(2.5, 2.6) + (1-vary_param) * np.pi]])
@@ -201,7 +201,7 @@ sampler.run_mcmc(pos,npoints)
 paramlabels_a = [r'$\alpha_'+str(i)+'$' for i in range (nmax+1)]
 paramlabels_b = [r'$\beta_'+str(i)+'$' for i in range (nmax+1)]
 paramlabels_x = [r'$x_'+str(i)+'$' for i in range (nmax+1)]
-paramlabels_y = [r'$y_'+str(i)+'$' for i in range (nmax+1)] if vary_fund == True else ['$y_0-\pi$'] + ['$y_'+str(i)+'$' for i in range (1, nmax+1)]
+paramlabels_y = [r'$y_'+str(i)+'$' for i in range (nmax+1)] if vary_fund == True else ['$y_'+str(i)+'-\pi$' for i in range (nmax+1)]
 paramlabels = paramlabels_a + paramlabels_b + paramlabels_x + paramlabels_y
 #Need to delete alpha_0 and alpha_1 for the corner plot
 paramlabels_corner = paramlabels_a + paramlabels_b + paramlabels_x + paramlabels_y
@@ -219,7 +219,7 @@ for i in range(4*(nmax+1)):
     axes[i].set_ylabel(paramlabels[i])
 axes[-1].set_xlabel('Iterations')
 #plt.show()
-fig.savefig(rootpath+'/git/rdstackingproject/plotsmc/vary'+str(vary_fund)+'nmax='+str(nmax)+'_tshift='+str(tshift)+'_'+str(npoints)+'pt_chain.pdf', format = 'pdf')
+fig.savefig(rootpath+'/git/rdstackingproject/plotsmc/vary'+str(vary_fund)+'nmax='+str(nmax)+'_tshift='+str(tshift)+'_'+str(npoints)+'pt_chain.png', format = 'png', bbox_inches = 'tight')
 
 
 
@@ -234,9 +234,11 @@ pk = samples[np.argmax(lglk)]
 #print('pk:')
 #print(pk)
 pk_corn = pk if vary_fund == True else np.delete(pk, [0,2])
-#y_0 range needs some messaging to make the plot
+#y_0 range needs some messaging to make the plot. But in order to make the whole picture consistent, better change the range of y_1 too.
 if vary_fund == False:
+    samples_corn.T[-1] = samples_corn.T[-1] - np.pi 
     samples_corn.T[-2] = samples_corn.T[-2] - np.pi #This indeed changes samples_corn itself
+    pk[-1] -= np.pi
     pk[-2] -= np.pi
 
 #print('pkFalse:')
@@ -289,11 +291,16 @@ logevidence = sampler.log_evidence_estimate(fburnin=0.2)
 df3 = pd.DataFrame(logevidence,index = [r'logevidence', r'error'], columns=[r'vary='+str(vary_fund)+', nmax='+str(nmax)+', tshift='+str(tshift)+', '+str(npoints)+'pts'])
 df3.to_hdf(rootpath+'/git/rdstackingproject/plotsmc/vary'+str(vary_fund)+'nmax='+str(nmax)+'_tshift='+str(tshift)+'_'+str(npoints)+'pt_logevidence.h5', key='df')
 
+#The autocorrelation time
+atcrr = sampler.get_autocorr_time()[0]
+df4 = pd.DataFrame(atcrr,index = paramlabels, columns=[r'Autocorrelation time, vary='+str(vary_fund)+', nmax='+str(nmax)+', tshift='+str(tshift)+', '+str(npoints)+'pts'])
+df4.to_hdf(rootpath+'/git/rdstackingproject/plotsmc/vary'+str(vary_fund)+'nmax='+str(nmax)+'_tshift='+str(tshift)+'_'+str(npoints)+'pt_autocrrtime.h5', key='df')    
+
 
 #Now plot alpha_0 on top of alpha_1 - only on top, not stacked, and only valid for vary_fund == True
 if vary_fund == True:
-    samplea0 = samples.T[0]
-    samplea1 = samples.T[1]
+    samplea0 = samplesT[0]
+    samplea1 = samplesT[1]
     fighist1 = plt.figure(figsize = (8, 6))
     n0, bins0, patches0 = plt.hist(samplea0, bins = numbins * 6, alpha = 0.5, label = r'$\alpha_0$')
     n1, bins1, patches1 = plt.hist(samplea1, bins = numbins * 10, alpha = 0.5, label = r'$\alpha_1$')
@@ -327,4 +334,4 @@ plt.xlim(timesrd[0], timesrd[0]+80)
 plt.xlabel("t")
 plt.ylabel("h")
 
-figband.savefig(rootpath+'/git/rdstackingproject/plotsmc/vary'+str(vary_fund)+'nmax='+str(nmax)+'_tshift='+str(tshift)+'_'+str(npoints)+'pt_band.pdf', format = 'pdf')
+figband.savefig(rootpath+'/git/rdstackingproject/plotsmc/vary'+str(vary_fund)+'nmax='+str(nmax)+'_tshift='+str(tshift)+'_'+str(npoints)+'pt_band.png', format = 'png', bbox_inches = 'tight')
\ No newline at end of file
diff --git a/code/RDGW150914_ptemcee4.py b/code/RDGW150914_ptemcee4.py
index 8eac5d0..e39042b 100755
--- a/code/RDGW150914_ptemcee4.py
+++ b/code/RDGW150914_ptemcee4.py
@@ -4,7 +4,7 @@ import warnings
 warnings.simplefilter("ignore") #Used to suppress the RuntimeWarnings. But if you are changing the code for TESTING, MAKE SURE TO COMMENT IT OUT.
 
 '''
-This script calculates the RDGW150914 constraints with ptemcee - specifically, the n=1 case both varying or not varying the fundamental frequencies. It produces the chain plot, corner plot, parameter constraints, and data plotted with the 1-sigma band. Since we are working specifically for the n=1 case, we also add in the corner plot the combined chain of alpha_0 and alpha_1, in order to demonstrate that the two are indistinguishable with each other.
+This script calculates the RDGW150914 constraints with ptemcee - specifically, the n=1 case not varying the fundamental frequencies, tshift = 19. It produces the chain plot, corner plot, parameter constraints, and data plotted with the 1-sigma band. Since we are working specifically for the n=1 case, we also compare alpha_0 and alpha_1 in a histogram, in order to demonstrate that the two are indistinguishable with each other.
 
 --- R^a_{yne} L^i_u, 08/09/2020
 '''
@@ -140,15 +140,15 @@ def log_prior(theta):
     y_0 = theta[6]
     y_1 = theta[7]
 
-    if all([nmax == 1, 0 <= tshift <= 5, vary_fund == True, -1.2 <= a_0 <= 1.2, -1.5 <= a_1 <= 1.5, -9.0 <= b_0 <= 9.0, -10.0 <= b_1 <= 21.0, 0 <= x_0 <= 1.6, 0 <= x_1 <= 1.4, 0 <= y_0 <= 2*np.pi, 0 <= y_1 <= 2*np.pi]):        
+    if all([nmax == 1, 0 <= tshift <= 5, vary_fund == True, -1.2 <= a_0 <= 1.2, -1.5 <= a_1 <= 1.5, -4.0 <= b_0 <= 9.0, -4.0 <= b_1 <= 21.0, 0 <= x_0 <= 1.6, 0 <= x_1 <= 1.4, 0 <= y_0 <= 2*np.pi, 0 <= y_1 <= 2*np.pi]):        
         return 0.0
-    elif all([nmax == 1, tshift == 19, vary_fund == True, -10.0 <= a_0 <= 10.0, -10.0 <= a_1 <= 10.0, -9.0 <= b_0 <= 20.0, -9.0 <= b_1 <= 25.0, 0 <= x_0 <= 0.8, 0 <= x_1 <= 0.9, 0 <= y_0 <= 2*np.pi, 0 <= y_1 <= 2*np.pi]):
+    elif all([nmax == 1, tshift == 19, vary_fund == True, -10.0 <= a_0 <= 10.0, -10.0 <= a_1 <= 10.0, -20.0 <= b_0 <= 30.0, -25.0 <= b_1 <= 30.0, 0 <= x_0 <= 0.6, 0 <= x_1 <= 0.9, 0 <= y_0 <= 2*np.pi, 0 <= y_1 <= 2*np.pi]):
         return 0.0
     
-#PAY EXTRA ATTENTION TO THESE TWO CASES, SINCE WE SHIFTED y_0. THE CORNER PLOTS LABELS NEED TO BE TREATED WITH CARE.
-    elif all([nmax == 1, 0 <= tshift <= 5, vary_fund == False, -10.0 <= a_0 <= 10.0, -1.5 <= a_1 <= 1.5, -9.0 <= b_0 <= 9.0, -10.0 <= b_1 <= 16.0, 0 <= x_0 <= 2.2, 0 <= x_1 <= 2.5, 0 <= y_0-np.pi <= 2*np.pi, 0 <= y_1 <= 2*np.pi,]):
+#PAY EXTRA ATTENTION TO THESE TWO CASES, SINCE WE SHIFTED y's. THE CORNER PLOTS LABELS NEED TO BE TREATED WITH CARE.
+    elif all([nmax == 1, 0 <= tshift <= 5, vary_fund == False, -10.0 <= a_0 <= 10.0, -1.5 <= a_1 <= 1.5, -9.0 <= b_0 <= 9.0, -6.0 <= b_1 <= 20.0, 0 <= x_0 <= 2.4, 0 <= x_1 <= 2.5, 0 <= y_0-np.pi <= 2*np.pi, 0 <= y_1-np.pi <= 2*np.pi,]):
         return 0.0
-    elif all([nmax == 1, tshift == 19, vary_fund == False, -10.0 <= a_0 <= 10.0, -10.0 <= a_1 <= 10.0, -9.0 <= b_0 <= 9.0, -10.0 <= b_1 <= 32.0, 0 <= x_0 <= 0.6, 0 <= x_1 <= 1.0, 0 <= y_0-np.pi <= 2*np.pi, 0 <= y_1 <= 2*np.pi]):
+    elif all([nmax == 1, tshift == 19, vary_fund == False, -10.0 <= a_0 <= 10.0, -8.0 <= a_1 <= 8.0, -9.0 <= b_0 <= 9.0, -10.0 <= b_1 <= 42.0, 0 <= x_0 <= 0.6, 0 <= x_1 <= 0.7, 0 <= y_0-np.pi <= 2*np.pi, 0 <= y_1-np.pi <= 2*np.pi]):
         return 0.0
 
     return -np.inf
@@ -178,8 +178,8 @@ def log_probability(theta):
 
 #Fit with ptemcee
 #Set the number of cores of your processors
-pool = choose_pool(28)
-pool.size = 28
+pool = choose_pool(24)
+pool.size = 24
 vary_param = float(vary_fund)
 np.random.seed(42)
 pos = np.array([[random.uniform(-0.1,0.1), random.uniform(-0.1,0.1), 4.28313743e-01, random.uniform(2.5, 2.6) + (1-vary_param) * np.pi]])
@@ -201,7 +201,7 @@ sampler.run_mcmc(pos,npoints)
 paramlabels_a = [r'$\alpha_'+str(i)+'$' for i in range (nmax+1)]
 paramlabels_b = [r'$\beta_'+str(i)+'$' for i in range (nmax+1)]
 paramlabels_x = [r'$x_'+str(i)+'$' for i in range (nmax+1)]
-paramlabels_y = [r'$y_'+str(i)+'$' for i in range (nmax+1)] if vary_fund == True else ['$y_0-\pi$'] + ['$y_'+str(i)+'$' for i in range (1, nmax+1)]
+paramlabels_y = [r'$y_'+str(i)+'$' for i in range (nmax+1)] if vary_fund == True else ['$y_'+str(i)+'-\pi$' for i in range (nmax+1)]
 paramlabels = paramlabels_a + paramlabels_b + paramlabels_x + paramlabels_y
 #Need to delete alpha_0 and alpha_1 for the corner plot
 paramlabels_corner = paramlabels_a + paramlabels_b + paramlabels_x + paramlabels_y
@@ -219,7 +219,7 @@ for i in range(4*(nmax+1)):
     axes[i].set_ylabel(paramlabels[i])
 axes[-1].set_xlabel('Iterations')
 #plt.show()
-fig.savefig(rootpath+'/git/rdstackingproject/plotsmc/vary'+str(vary_fund)+'nmax='+str(nmax)+'_tshift='+str(tshift)+'_'+str(npoints)+'pt_chain.pdf', format = 'pdf')
+fig.savefig(rootpath+'/git/rdstackingproject/plotsmc/vary'+str(vary_fund)+'nmax='+str(nmax)+'_tshift='+str(tshift)+'_'+str(npoints)+'pt_chain.png', format = 'png', bbox_inches = 'tight')
 
 
 
@@ -234,9 +234,11 @@ pk = samples[np.argmax(lglk)]
 #print('pk:')
 #print(pk)
 pk_corn = pk if vary_fund == True else np.delete(pk, [0,2])
-#y_0 range needs some messaging to make the plot
+#y_0 range needs some messaging to make the plot. But in order to make the whole picture consistent, better change the range of y_1 too.
 if vary_fund == False:
-    samples_corn.T[-2] = samples_corn.T[-2] - np.pi
+    samples_corn.T[-1] = samples_corn.T[-1] - np.pi 
+    samples_corn.T[-2] = samples_corn.T[-2] - np.pi #This indeed changes samples_corn itself
+    pk[-1] -= np.pi
     pk[-2] -= np.pi
 
 #print('pkFalse:')
@@ -250,6 +252,7 @@ median = np.median(samples_corn, axis=0)
 
 figcorn = corner.corner(samples_corn, bins = numbins, hist_bin_factor = 5, color = datacolor, truths=pk_corn, truth_color = pkcolor, plot_contours = True, labels = paramlabels_corner, quantiles=(0.05, 0.16, 0.5, 0.84, 0.95), levels=[1-np.exp(-0.5), 1-np.exp(-1.64 ** 2/2)], show_titles=True)
 
+
 #Extract the axes in order to add more important line plots
 naxes = len(pk_corn)
 axes = np.array(figcorn.axes).reshape((naxes, naxes))
@@ -288,12 +291,16 @@ logevidence = sampler.log_evidence_estimate(fburnin=0.2)
 df3 = pd.DataFrame(logevidence,index = [r'logevidence', r'error'], columns=[r'vary='+str(vary_fund)+', nmax='+str(nmax)+', tshift='+str(tshift)+', '+str(npoints)+'pts'])
 df3.to_hdf(rootpath+'/git/rdstackingproject/plotsmc/vary'+str(vary_fund)+'nmax='+str(nmax)+'_tshift='+str(tshift)+'_'+str(npoints)+'pt_logevidence.h5', key='df')
 
+#The autocorrelation time, this time I don't intend to delete alpha0 and beta0, just keep the nan's there since we don't really use this data, only look at it. It's good to keep it for sanity check's sake.
+atcrr = sampler.get_autocorr_time()[0]
+df4 = pd.DataFrame(atcrr,index = paramlabels, columns=[r'Autocorrelation time, vary='+str(vary_fund)+', nmax='+str(nmax)+', tshift='+str(tshift)+', '+str(npoints)+'pts'])
+df4.to_hdf(rootpath+'/git/rdstackingproject/plotsmc/vary'+str(vary_fund)+'nmax='+str(nmax)+'_tshift='+str(tshift)+'_'+str(npoints)+'pt_autocrrtime.h5', key='df')    
 
 
 #Now plot alpha_0 on top of alpha_1 - only on top, not stacked, and only valid for vary_fund == True
 if vary_fund == True:
-    samplea0 = samples.T[0]
-    samplea1 = samples.T[1]
+    samplea0 = samplesT[0]
+    samplea1 = samplesT[1]
     fighist1 = plt.figure(figsize = (8, 6))
     n0, bins0, patches0 = plt.hist(samplea0, bins = numbins * 6, alpha = 0.5, label = r'$\alpha_0$')
     n1, bins1, patches1 = plt.hist(samplea1, bins = numbins * 10, alpha = 0.5, label = r'$\alpha_1$')
@@ -327,4 +334,4 @@ plt.xlim(timesrd[0], timesrd[0]+80)
 plt.xlabel("t")
 plt.ylabel("h")
 
-figband.savefig(rootpath+'/git/rdstackingproject/plotsmc/vary'+str(vary_fund)+'nmax='+str(nmax)+'_tshift='+str(tshift)+'_'+str(npoints)+'pt_band.pdf', format = 'pdf')
+figband.savefig(rootpath+'/git/rdstackingproject/plotsmc/vary'+str(vary_fund)+'nmax='+str(nmax)+'_tshift='+str(tshift)+'_'+str(npoints)+'pt_band.png', format = 'png', bbox_inches = 'tight')
\ No newline at end of file
diff --git a/code/condor_submit_RdownPtemcee1.sub b/code/condor_submit_RdownPtemcee1.sub
index 52aa793..d7afc43 100755
--- a/code/condor_submit_RdownPtemcee1.sub
+++ b/code/condor_submit_RdownPtemcee1.sub
@@ -13,7 +13,7 @@ notify_user = rl746@cornell.edu
 notification = Complete
 arguments  = "-processid $(Process)" 
 request_memory = 232GB
-request_cpus = 28
+request_cpus = 24
 on_exit_remove = (ExitBySignal == False) || ((ExitBySignal == True) && (ExitSignal != 11))
 accounting_group = aei.dev.test_dynesty
 queue 1
diff --git a/code/condor_submit_RdownPtemcee3.sub b/code/condor_submit_RdownPtemcee3.sub
index 781ec9b..7e1fbd5 100755
--- a/code/condor_submit_RdownPtemcee3.sub
+++ b/code/condor_submit_RdownPtemcee3.sub
@@ -13,7 +13,7 @@ notify_user = rl746@cornell.edu
 notification = Complete
 arguments  = "-processid $(Process)" 
 request_memory = 232GB
-request_cpus = 28
+request_cpus = 24
 on_exit_remove = (ExitBySignal == False) || ((ExitBySignal == True) && (ExitSignal != 11))
 accounting_group = aei.dev.test_dynesty
 queue 1
-- 
GitLab