Skip to content
Snippets Groups Projects
Commit aadc12f6 authored by Rayne Liu's avatar Rayne Liu
Browse files

Adjustments of ranges - some new phenomena about b1 need to be looked at

parent 131d0fcb
Branches
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
code/RDGW150914_ptemcee1.py
+ 10
− 7
View file @ aadc12f6
......@@ -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
code/RDGW150914_ptemcee2.py
+ 9
− 16
View file @ aadc12f6
......@@ -39,11 +39,11 @@ 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
code/RDGW150914_ptemcee3.py
+ 6
− 6
View file @ aadc12f6
......@@ -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
......
code/RDGW150914_ptemcee4.py
+ 9
− 9
View file @ aadc12f6
......@@ -39,11 +39,11 @@ 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):
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment

AEI Privacy Policy Imprint