From aadc12f69e9d89b141472f2e47ae2a24e2ae1286 Mon Sep 17 00:00:00 2001
From: Rayne Liu <rl746@cornell.edu>
Date: Sun, 16 Aug 2020 22:08:27 -0400
Subject: [PATCH] Adjustments of ranges - some new phenomena about b1 need to
be looked at
---
code/RDGW150914_ptemcee1.py | 17 ++++++++++-------
code/RDGW150914_ptemcee2.py | 25 +++++++++----------------
code/RDGW150914_ptemcee3.py | 12 ++++++------
code/RDGW150914_ptemcee4.py | 18 +++++++++---------
4 files changed, 34 insertions(+), 38 deletions(-)
diff --git a/code/RDGW150914_ptemcee1.py b/code/RDGW150914_ptemcee1.py
index 6025d8f..d31fa93 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 19f8eea..e24ae8e 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 f50cf32..4183c17 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 848d6a7..6124cdc 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):
--
GitLab