diff --git a/code/RDGW150914_ptemcee1.py b/code/RDGW150914_ptemcee1.py
index 6025d8f806e56fc1cf5b90ce32965a4ce100ea24..d31fa93dd0e9ca92996c7c1d6cbfae0ca4840498 100755
--- a/code/RDGW150914_ptemcee1.py
+++ b/code/RDGW150914_ptemcee1.py
@@ -33,17 +33,17 @@ from scipy.optimize import minimize
#tshift: time shift after the strain peak
#vary_fund: whether you vary the fundamental frequency. Works in the model_dv function.
-rootpath= "/Users/RayneLiu"#"/work/rayne.liu"
+rootpath= "/work/rayne.liu"#"/Users/RayneLiu"
nmax=1
tshift=3.6
vary_fund = True
#sampler parameters
-npoints=3000
+npoints=60002
nwalkers = 42
ntemps=12
ndim = int(4*(nmax+1))
-burnin = 1500 #How many points do you burn before doing the corner plot. You need to watch the convergence of the chain plot a bit.
+burnin = 5000 #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'
@@ -137,15 +137,15 @@ def log_prior(theta):
y_0 = theta[6]
y_1 = theta[7]
- if all([nmax == 1, 0 <= tshift <= 5, vary_fund == True, -0.6 <= a_0 <= 0.6, -2.0 <= a_1 <= 2.0, -1.0 <= b_0 <= 3.0, -1.0 <= b_1 <= 2.0, 0 <= x_0 <= 1.6, 0 <= x_1 <= 1.6, 0 <= y_0 <= 2*np.pi, 0 <= y_1 <= 2*np.pi]):
+ if all([nmax == 1, 0 <= tshift <= 5, vary_fund == True, -0.8 <= a_0 <= 0.8, -1.0 <= a_1 <= 1.0, -1.0 <= b_0 <= 9.0, -1.0 <= b_1 <= 21.0, 0 <= x_0 <= 1.6, 0 <= x_1 <= 1.6, 0 <= y_0 <= 2*np.pi, 0 <= y_1 <= 2*np.pi]):
return 0.0
- elif all([nmax == 1, tshift == 19, vary_fund == True, 0 <= x_0 <= 1.5, 0 <= y_0 <= 2*np.pi, -1.0 <= a_0 <= 1.0, -1.0 <= b_0 <= 2.0, 0 <= x_1 <= 1.8, 0 <= y_1 <= 2*np.pi, -1.0 <= a_1 <= 1.0, -1.0 <= b_1 <= 2.0]):
+ elif all([nmax == 1, tshift == 19, vary_fund == True, -10.0 <= a_0 <= 10.0, -10.0 <= a_1 <= 10.0, -1.0 <= b_0 <= 10.0, -1.0 <= b_1 <= 10.0, 0 <= x_0 <= 1.0, 0 <= x_1 <= 1.2, 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, 0 <= x_0 <= 1.6, 0 <= y_0-np.pi <= 2*np.pi, -1.0 <= a_0 <= 1.0, -1.0 <= b_0 <= 1.0, 0 <= x_1 <= 2.0, 0 <= y_1 <= 2*np.pi, -1.0 <= a_1 <= 1.2, -1.0 <= b_1 <= 4.2]):
+ elif all([nmax == 1, 0 <= tshift <= 5, vary_fund == False, -1.0 <= a_0 <= 1.0, -1.5 <= a_1 <= 1.5, -1.0 <= b_0 <= 1.0, -1.0 <= b_1 <= 12.0, 0 <= x_0 <= 2.0, 0 <= x_1 <= 2.4, 0 <= y_0-np.pi <= 2*np.pi, 0 <= y_1 <= 2*np.pi,]):
return 0.0
- elif all([nmax == 1, tshift == 19, vary_fund == False, 0 <= x_0 <= 0.6, 0 <= y_0-np.pi <= 2*np.pi, -1.0 <= a_0 <= 1.0, -1.0 <= b_0 <= 1.0, 0 <= x_1 <= 1.2, 0 <= y_1 <= 2*np.pi, -1.0 <= a_1 <= 1.0, -1.0 <= b_1 <= 2.0]):
+ elif all([nmax == 1, tshift == 19, vary_fund == False, -1.0 <= a_0 <= 1.0, -10.0 <= a_1 <= 10.0, -1.0 <= b_0 <= 1.0, -1.0 <= b_1 <= 10.0, 0 <= x_0 <= 0.6, 0 <= x_1 <= 1.2, 0 <= y_0-np.pi <= 2*np.pi, 0 <= y_1 <= 2*np.pi]):
return 0.0
return -np.inf
@@ -294,3 +294,6 @@ 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')
+
+#import os
+#os.system('afplay /System/Library/Sounds/Submarine.aiff')
\ No newline at end of file
diff --git a/code/RDGW150914_ptemcee2.py b/code/RDGW150914_ptemcee2.py
index 19f8eea566cd166e7118de15e4896eb1b8b0625f..e24ae8ed0109f1b8b84d6a548f2337cf3adc705b 100755
--- a/code/RDGW150914_ptemcee2.py
+++ b/code/RDGW150914_ptemcee2.py
@@ -33,17 +33,17 @@ from scipy.optimize import minimize
#tshift: time shift after the strain peak
#vary_fund: whether you vary the fundamental frequency. Works in the model_dv function.
-rootpath="/work/rayne.liu" #"/Users/RayneLiu"
+rootpath="/work/rayne.liu"#"/Users/RayneLiu"
nmax=1
tshift=19
vary_fund = True
#sampler parameters
-npoints=1000000
+npoints=60001
nwalkers = 42
ntemps=12
ndim = int(4*(nmax+1))
-burnin = 200000 #How many points do you burn before doing the corner plot. You need to watch the convergence of the chain plot a bit.
+burnin = 5000 #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'
@@ -137,15 +137,15 @@ def log_prior(theta):
y_0 = theta[6]
y_1 = theta[7]
- if all([nmax == 1, 0 <= tshift <= 5, vary_fund == True, 0 <= x_0 <= 2.0, 0 <= y_0 <= 2*np.pi, -0.4 <= a_0 <= 0.4, -1.0 <= b_0 <= 1.6, 0 <= x_1 <= 1.8, 0 <= y_1 <= 2*np.pi, -1.0 <= a_1 <= 1.6, -1.0 <= b_1 <= 2.0]):
+ if all([nmax == 1, 0 <= tshift <= 5, vary_fund == True, -0.8 <= a_0 <= 0.8, -1.0 <= a_1 <= 1.0, -1.0 <= b_0 <= 9.0, -1.0 <= b_1 <= 21.0, 0 <= x_0 <= 1.6, 0 <= x_1 <= 1.6, 0 <= y_0 <= 2*np.pi, 0 <= y_1 <= 2*np.pi]):
return 0.0
- elif all([nmax == 1, tshift == 19, vary_fund == True, 0 <= x_0 <= 1.5, 0 <= y_0 <= 2*np.pi, -1.0 <= a_0 <= 3.0, -1.0 <= b_0 <= 2.0, 0 <= x_1 <= 1.8, 0 <= y_1 <= 2*np.pi, -1.0 <= a_1 <= 3.0, -1.0 <= b_1 <= 2.0]):
+ elif all([nmax == 1, tshift == 19, vary_fund == True, -10.0 <= a_0 <= 10.0, -10.0 <= a_1 <= 10.0, -1.0 <= b_0 <= 10.0, -1.0 <= b_1 <= 10.0, 0 <= x_0 <= 1.0, 0 <= x_1 <= 1.2, 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, 0 <= x_0 <= 2.0, 0 <= y_0-np.pi <= 2*np.pi, -1.0 <= a_0 <= 1.0, -1.0 <= b_0 <= 1.0, 0 <= x_1 <= 1.6, 0 <= y_1 <= 2*np.pi, -1.0 <= a_1 <= 1.2, -1.0 <= b_1 <= 2.8]):
+ elif all([nmax == 1, 0 <= tshift <= 5, vary_fund == False, -1.0 <= a_0 <= 1.0, -1.5 <= a_1 <= 1.5, -1.0 <= b_0 <= 1.0, -1.0 <= b_1 <= 12.0, 0 <= x_0 <= 2.0, 0 <= x_1 <= 2.4, 0 <= y_0-np.pi <= 2*np.pi, 0 <= y_1 <= 2*np.pi,]):
return 0.0
- elif all([nmax == 1, tshift == 19, vary_fund == False, 0 <= x_0 <= 0.6, 0 <= y_0-np.pi <= 2*np.pi, -1.0 <= a_0 <= 1.0, -1.0 <= b_0 <= 1.0, 0 <= x_1 <= 1.2, 0 <= y_1 <= 2*np.pi, -1.0 <= a_1 <= 3.0, -1.0 <= b_1 <= 2.0]):
+ elif all([nmax == 1, tshift == 19, vary_fund == False, -1.0 <= a_0 <= 1.0, -10.0 <= a_1 <= 10.0, -1.0 <= b_0 <= 1.0, -1.0 <= b_1 <= 10.0, 0 <= x_0 <= 0.6, 0 <= x_1 <= 1.2, 0 <= y_0-np.pi <= 2*np.pi, 0 <= y_1 <= 2*np.pi]):
return 0.0
return -np.inf
@@ -295,13 +295,6 @@ 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')
-"""
-
-
-
-
-
-
-
-"""
+import os
+os.system('afplay /System/Library/Sounds/Submarine.aiff')
\ No newline at end of file
diff --git a/code/RDGW150914_ptemcee3.py b/code/RDGW150914_ptemcee3.py
index f50cf3250ed3e0d02c43e2ef70cb324b8317e16c..4183c172c9398ff6dd1fdd3db8e8afb0359d2057 100755
--- a/code/RDGW150914_ptemcee3.py
+++ b/code/RDGW150914_ptemcee3.py
@@ -39,11 +39,11 @@ tshift=0
vary_fund = False
#sampler parameters
-npoints=1002
+npoints=60000
nwalkers = 42
ntemps=12
ndim = int(4*(nmax+1))
-burnin = 200 #How many points do you burn before doing the corner plot. You need to watch the convergence of the chain plot a bit.
+burnin = 5000 #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'
@@ -137,15 +137,15 @@ def log_prior(theta):
y_0 = theta[6]
y_1 = theta[7]
- if all([nmax == 1, 0 <= tshift <= 5, vary_fund == True, 0 <= x_0 <= 2.0, 0 <= y_0 <= 2*np.pi, -0.4 <= a_0 <= 0.4, -1.0 <= b_0 <= 1.6, 0 <= x_1 <= 1.8, 0 <= y_1 <= 2*np.pi, -1.0 <= a_1 <= 1.6, -1.0 <= b_1 <= 2.0]):
+ if all([nmax == 1, 0 <= tshift <= 5, vary_fund == True, -0.8 <= a_0 <= 0.8, -1.0 <= a_1 <= 1.0, -1.0 <= b_0 <= 9.0, -1.0 <= b_1 <= 21.0, 0 <= x_0 <= 1.6, 0 <= x_1 <= 1.6, 0 <= y_0 <= 2*np.pi, 0 <= y_1 <= 2*np.pi]):
return 0.0
- elif all([nmax == 1, tshift == 19, vary_fund == True, 0 <= x_0 <= 1.5, 0 <= y_0 <= 2*np.pi, -1.0 <= a_0 <= 1.0, -1.0 <= b_0 <= 2.0, 0 <= x_1 <= 1.8, 0 <= y_1 <= 2*np.pi, -1.0 <= a_1 <= 1.0, -1.0 <= b_1 <= 2.0]):
+ elif all([nmax == 1, tshift == 19, vary_fund == True, -10.0 <= a_0 <= 10.0, -10.0 <= a_1 <= 10.0, -1.0 <= b_0 <= 10.0, -1.0 <= b_1 <= 10.0, 0 <= x_0 <= 1.0, 0 <= x_1 <= 1.2, 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, 0 <= x_0 <= 1.6, 0 <= y_0-np.pi <= 2*np.pi, -1.0 <= a_0 <= 1.0, -1.0 <= b_0 <= 1.0, 0 <= x_1 <= 2.0, 0 <= y_1 <= 2*np.pi, -1.0 <= a_1 <= 1.2, -1.0 <= b_1 <= 4.2]):
+ elif all([nmax == 1, 0 <= tshift <= 5, vary_fund == False, -1.0 <= a_0 <= 1.0, -1.5 <= a_1 <= 1.5, -1.0 <= b_0 <= 1.0, -1.0 <= b_1 <= 12.0, 0 <= x_0 <= 2.0, 0 <= x_1 <= 2.4, 0 <= y_0-np.pi <= 2*np.pi, 0 <= y_1 <= 2*np.pi,]):
return 0.0
- elif all([nmax == 1, tshift == 19, vary_fund == False, 0 <= x_0 <= 0.6, 0 <= y_0-np.pi <= 2*np.pi, -1.0 <= a_0 <= 1.0, -1.0 <= b_0 <= 1.0, 0 <= x_1 <= 1.2, 0 <= y_1 <= 2*np.pi, -1.0 <= a_1 <= 1.0, -1.0 <= b_1 <= 2.0]):
+ elif all([nmax == 1, tshift == 19, vary_fund == False, -1.0 <= a_0 <= 1.0, -10.0 <= a_1 <= 10.0, -1.0 <= b_0 <= 1.0, -1.0 <= b_1 <= 10.0, 0 <= x_0 <= 0.6, 0 <= x_1 <= 1.2, 0 <= y_0-np.pi <= 2*np.pi, 0 <= y_1 <= 2*np.pi]):
return 0.0
return -np.inf
diff --git a/code/RDGW150914_ptemcee4.py b/code/RDGW150914_ptemcee4.py
index 848d6a72b4a9e7e41fc76d67d9585b153e880bb6..6124cdc8f3cc1344e3e05743f13e40ccc067d5db 100755
--- a/code/RDGW150914_ptemcee4.py
+++ b/code/RDGW150914_ptemcee4.py
@@ -33,17 +33,17 @@ from scipy.optimize import minimize
#tshift: time shift after the strain peak
#vary_fund: whether you vary the fundamental frequency. Works in the model_dv function.
-rootpath="/work/rayne.liu"# "/Users/RayneLiu"
+rootpath= "/work/rayne.liu"#"/Users/RayneLiu"
nmax=1
tshift=19
vary_fund = False
#sampler parameters
-npoints=1000000
+npoints=1002
nwalkers = 42
ntemps=12
ndim = int(4*(nmax+1))
-burnin = 420000 #How many points do you burn before doing the corner plot. You need to watch the convergence of the chain plot a bit.
+burnin = 20 #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'
@@ -137,15 +137,15 @@ def log_prior(theta):
y_0 = theta[6]
y_1 = theta[7]
- if all([nmax == 1, 0 <= tshift <= 5, vary_fund == True, 0 <= x_0 <= 2.0, 0 <= y_0 <= 2*np.pi, -0.4 <= a_0 <= 0.4, -1.0 <= b_0 <= 1.6, 0 <= x_1 <= 1.8, 0 <= y_1 <= 2*np.pi, -1.0 <= a_1 <= 1.6, -1.0 <= b_1 <= 2.0]):
+ if all([nmax == 1, 0 <= tshift <= 5, vary_fund == True, -0.8 <= a_0 <= 0.8, -1.0 <= a_1 <= 1.0, -1.0 <= b_0 <= 9.0, -1.0 <= b_1 <= 21.0, 0 <= x_0 <= 1.6, 0 <= x_1 <= 1.6, 0 <= y_0 <= 2*np.pi, 0 <= y_1 <= 2*np.pi]):
return 0.0
- elif all([nmax == 1, tshift == 19, vary_fund == True, 0 <= x_0 <= 1.5, 0 <= y_0 <= 2*np.pi, -1.0 <= a_0 <= 1.0, -1.0 <= b_0 <= 2.0, 0 <= x_1 <= 1.8, 0 <= y_1 <= 2*np.pi, -1.0 <= a_1 <= 1.0, -1.0 <= b_1 <= 2.0]):
+ elif all([nmax == 1, tshift == 19, vary_fund == True, -10.0 <= a_0 <= 10.0, -10.0 <= a_1 <= 10.0, -1.0 <= b_0 <= 10.0, -1.0 <= b_1 <= 10.0, 0 <= x_0 <= 1.0, 0 <= x_1 <= 1.2, 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, 0 <= x_0 <= 1.6, 0 <= y_0-np.pi <= 2*np.pi, -1.0 <= a_0 <= 1.0, -1.0 <= b_0 <= 1.0, 0 <= x_1 <= 2.0, 0 <= y_1 <= 2*np.pi, -1.0 <= a_1 <= 1.2, -1.0 <= b_1 <= 4.2]):
+ elif all([nmax == 1, 0 <= tshift <= 5, vary_fund == False, -1.0 <= a_0 <= 1.0, -1.5 <= a_1 <= 1.5, -1.0 <= b_0 <= 1.0, -1.0 <= b_1 <= 12.0, 0 <= x_0 <= 2.0, 0 <= x_1 <= 2.4, 0 <= y_0-np.pi <= 2*np.pi, 0 <= y_1 <= 2*np.pi,]):
return 0.0
- elif all([nmax == 1, tshift == 19, vary_fund == False, 0 <= x_0 <= 0.6, 0 <= y_0-np.pi <= 2*np.pi, -1.0 <= a_0 <= 1.0, -1.0 <= b_0 <= 1.0, 0 <= x_1 <= 1.2, 0 <= y_1 <= 2*np.pi, -1.0 <= a_1 <= 1.0, -1.0 <= b_1 <= 2.0]):
+ elif all([nmax == 1, tshift == 19, vary_fund == False, -1.0 <= a_0 <= 1.0, -10.0 <= a_1 <= 10.0, -1.0 <= b_0 <= 1.0, -1.0 <= b_1 <= 10.0, 0 <= x_0 <= 0.6, 0 <= x_1 <= 1.2, 0 <= y_0-np.pi <= 2*np.pi, 0 <= y_1 <= 2*np.pi]):
return 0.0
return -np.inf
@@ -172,8 +172,8 @@ def log_probability(theta):
#Fit with ptemcee
#Set the number of cores of your processors
-pool = choose_pool(12)
-pool.size = 12
+pool = choose_pool(6)
+pool.size = 6
vary_param = float(vary_fund)
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]])
for i in range (1,nmax+1):