From fe566338acbfeb344dee82605449553d2f61c194 Mon Sep 17 00:00:00 2001
From: Rayne Liu <rl746@cornell.edu>
Date: Fri, 21 Aug 2020 01:14:34 -0400
Subject: [PATCH] Included the nmax=0 scenario
---
code/RDGW150914_ptemcee1.py | 102 +++++++++++++++++-------------------
1 file changed, 49 insertions(+), 53 deletions(-)
diff --git a/code/RDGW150914_ptemcee1.py b/code/RDGW150914_ptemcee1.py
index f5489fa..6f36bc1 100755
--- a/code/RDGW150914_ptemcee1.py
+++ b/code/RDGW150914_ptemcee1.py
@@ -7,6 +7,10 @@ warnings.simplefilter("ignore") #Used to suppress the RuntimeWarnings. But if yo
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
+
+Adapted to include the nmax = 0 case.
+
+---R^a_{yne} L^i_u, 08/20/2020, don't know what to do now
'''
import numpy as np
@@ -20,14 +24,12 @@ plt.rcParams.update({'font.size': 16.5})
import ptemcee
from pycbc.pool import choose_pool
-import math
import h5py
import inspect
import pandas as pd
import json
import qnm
import random
-from scipy.optimize import minimize
#Remember to change the following global variables
#rootpath: root path to nr data
@@ -36,17 +38,18 @@ 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"
-nmax=1
+rootpath= "/Users/RayneLiu"#"/work/rayne.liu"
+nmax=0
tshift=0
vary_fund = True
#sampler parameters
-npoints=19000
+npoints=2400
nwalkers = 1200
-ntemps=137
-ndim = int(4*(nmax+1))
-burnin = 9000 #How many points do you burn before doing the corner plot. You need to watch the convergence of the chain plot a bit.
+ntemps=16
+dim = nmax+1
+ndim = 4*dim
+burnin = 1200 #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'
@@ -94,8 +97,9 @@ gw_sxs_bbh_0305rd=gw_sxs_bbh_0305[position:-1]
timesrd=gw_sxs_bbh_0305[position:-1][:,0]
# Depending on nmax, you load nmax number of freqs. and damping times from the qnm package
-omegas = [qnm.modes_cache(s=-2,l=2,m=2,n=i)(a=af)[0] for i in range (0,nmax+1)]
-
+omegas = [qnm.modes_cache(s=-2,l=2,m=2,n=i)(a=af)[0] for i in range (0,dim)]
+w = (np.real(omegas))/mf
+tau=-1/(np.imag(omegas))*mf
@@ -104,15 +108,12 @@ omegas = [qnm.modes_cache(s=-2,l=2,m=2,n=i)(a=af)[0] for i in range (0,nmax+1)]
def model_dv(theta):
#x0, y0= theta
#Your nmax might not align with the dim of theta. Better check it here.
- assert int(len(theta)/4) == nmax + 1, 'Please recheck your n and parameters'
- w = (np.real(omegas))/mf
- tau=-1/(np.imag(omegas))*mf
- dim =int(len(theta)/4)
+ assert int(len(theta)/4) == dim, 'Please recheck your n and parameters'
- avars = theta[ : (nmax+1)]
- bvars = theta[(nmax+1) : 2*(nmax+1)]
- xvars = theta[2*(nmax+1) : 3*(nmax+1)]
- yvars = theta[3*(nmax+1) : ]
+ avars = theta[ : (dim)]
+ bvars = theta[(dim) : 2*(dim)]
+ xvars = theta[2*(dim) : 3*(dim)]
+ yvars = theta[3*(dim) : ]
if vary_fund == False:
avars[0]=0
@@ -131,25 +132,25 @@ def model_dv(theta):
#It works for the (xn*Exp[iyn]) version.
def log_prior(theta):
#Warning: we are specifically working with nmax=1 so here individual prior to the parameters are manually adjusted. This does not apply to all other nmax's.
- a_0 = theta[0]
- a_1 = theta[1]
- b_0 = theta[2]
- b_1 = theta[3]
- x_0 = theta[4]
- x_1 = theta[5]
- 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, -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, -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
+ avars = theta[ : (dim)]
+ bvars = theta[(dim) : 2*(dim)]
+ xvars = theta[2*(dim) : 3*(dim)]
+ yvars = theta[3*(dim) : ]
-#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, -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
+ if nmax == 0:
+ if all([0 <= tshift <= 5, vary_fund == True, -0.45 <= avars[0] <= 0.05, -0.95 <= bvars[0] <= 0.95, 0 <= xvars[0] <= 2.4, -np.pi <= yvars[0] <= np.pi]):
+ return 0.0
+ elif nmax == 1:
+ if all([0 <= tshift <= 5, vary_fund == True, -1.2 <= avars[0] <= 1.2, -1.5 <= avars[1] <= 1.5, -4.0 <= bvars[0] <= 9.0, -4.0 <= bvars[1] <= 21.0, 0 <= xvars[0] <= 1.6, 0 <= xvars[1] <= 1.4, 0 <= yvars[0] <= 2*np.pi, 0 <= yvars[1] <= 2*np.pi]):
+ return 0.0
+ elif all([tshift == 19, vary_fund == True, -10.0 <= avars[0] <= 10.0, -10.0 <= avars[1] <= 10.0, -20.0 <= bvars[0] <= 30.0, -25.0 <= bvars[1] <= 30.0, 0 <= xvars[0] <= 0.6, 0 <= xvars[1] <= 0.9, 0 <= yvars[0] <= 2*np.pi, 0 <= yvars[1] <= 2*np.pi]):
+ return 0.0
+
+ #PAY EXTRA ATTENTION TO THESE TWO CASES, SINCE WE SHIFTED y's. THE CORNER PLOTS LABELS NEED TO BE TREATED WITH CARE.
+ elif all([0 <= tshift <= 5, vary_fund == False, -10.0 <= avars[0] <= 10.0, -1.5 <= avars[1] <= 1.5, -9.0 <= bvars[0] <= 9.0, -6.0 <= bvars[1] <= 20.0, 0 <= xvars[0] <= 2.4, 0 <= xvars[1] <= 2.5, 0 <= yvars[0]-np.pi <= 2*np.pi, 0 <= yvars[1]-np.pi <= 2*np.pi,]):
+ return 0.0
+ elif all([tshift == 19, vary_fund == False, -10.0 <= avars[0] <= 10.0, -8.0 <= avars[1] <= 8.0, -9.0 <= bvars[0] <= 9.0, -10.0 <= bvars[1] <= 42.0, 0 <= xvars[0] <= 0.6, 0 <= xvars[1] <= 0.7, 0 <= yvars[0]-np.pi <= 2*np.pi, 0 <= yvars[1]-np.pi <= 2*np.pi]):
+ return 0.0
return -np.inf
@@ -167,24 +168,21 @@ def log_likelihood(theta):
# The evidence is just a normalization factor
def log_probability(theta):
lp = log_prior(theta)
- print('lp:')
- print(lp)
if not np.isfinite(lp):
return -np.inf
return lp + log_likelihood(theta)
-
#Fit with ptemcee
#Set the number of cores of your processors
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]])
-for i in range (1,nmax+1):
- pos_aux = np.array([[random.uniform(-0.1,0.1), random.uniform(-0.1,0.1), random.uniform(0.3,0.4), random.uniform(2.01, 2.02) + (1-vary_param) * np.pi]])
+pos = np.array([[random.uniform(-0.1,0.1), random.uniform(-0.1,0.1), 4.28313743e-01, random.uniform(0.5, 0.6) + (1-vary_param) * np.pi]])
+for i in range (1,dim):
+ pos_aux = np.array([[random.uniform(-0.1,0.1), random.uniform(-0.1,0.1), random.uniform(0.3,0.4), random.uniform(0.01, 0.02) + (1-vary_param) * np.pi]])
pos = np.concatenate((pos, pos_aux), axis = 0)
pos = pos.T.flatten()
@@ -198,10 +196,10 @@ 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)]
-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_'+str(i)+'-\pi$' for i in range (nmax+1)]
+paramlabels_a = [r'$\alpha_'+str(i)+'$' for i in range (dim)]
+paramlabels_b = [r'$\beta_'+str(i)+'$' for i in range (dim)]
+paramlabels_x = [r'$x_'+str(i)+'$' for i in range (dim)]
+paramlabels_y = [r'$y_'+str(i)+'$' for i in range (dim)] if vary_fund == True else ['$y_'+str(i)+'-\pi$' for i in range (dim)]
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
@@ -212,8 +210,8 @@ 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)):
+fig, axes = plt.subplots(ndim, 1, sharex=True, figsize=(12, 9*(dim)))
+for i in range(ndim):
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])
@@ -236,10 +234,8 @@ pk = samples[np.argmax(lglk)]
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
+ samples_corn.T[-dim:] -= np.pi #This indeed changes samples_corn itself
+ pk[-dim:] -= np.pi
#print('pkFalse:')
#print(pk)
@@ -298,7 +294,7 @@ df4.to_hdf(rootpath+'/git/rdstackingproject/plotsmc/vary'+str(vary_fund)+'nmax='
#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:
+if nmax == 1 and vary_fund == True:
samplea0 = samplesT[0]
samplea1 = samplesT[1]
fighist1 = plt.figure(figsize = (8, 6))
--
GitLab