From eff154bc39d763b1174b7e64d020b73536f447b7 Mon Sep 17 00:00:00 2001
From: Rayne Liu <rl746@cornell.edu>
Date: Wed, 7 Oct 2020 23:00:39 +0000
Subject: [PATCH] Added waveform band plotting
---
code/NR_Interpolate-0001_t_10M_wandt.py | 32 +++++++++++++++++++++----
1 file changed, 27 insertions(+), 5 deletions(-)
diff --git a/code/NR_Interpolate-0001_t_10M_wandt.py b/code/NR_Interpolate-0001_t_10M_wandt.py
index 63002fc..045ee44 100755
--- a/code/NR_Interpolate-0001_t_10M_wandt.py
+++ b/code/NR_Interpolate-0001_t_10M_wandt.py
@@ -42,12 +42,12 @@ tshift=10
vary_fund = True
#sampler parameters
-npoints = 1200
-nwalkers = 256
+npoints = 100#1200
+nwalkers = 30#256
ntemps=12
dim = nmax+1
ndim = 4*dim
-burnin = 600 #How many points do you burn before doing the corner plot. You need to watch the convergence of the chain plot a bit.
+burnin = 50#600 #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! I usually use half the points.
numbins = 32 #corner plot parameter - how many bins you want
datacolor = '#105670' #'#4fa3a7'
@@ -180,8 +180,8 @@ def log_prior(theta):
return 0.0
elif tshift == 10:
- if all([0.54 <= omega0 <= 0.6, 0.4 <= omega1 <= 0.58, 8.8 <= tau0 <= 14.2, 3.7 <= tau1 <= 10., \
- 0. <= xvar0 <= 1.0, 0. <= xvar1 <= 1.2, 0 <= yvar0 <= 2*np.pi, -np.pi <= yvar1 <= np.pi]):
+ if all([0.54 <= omega0 <= 0.6, 0.35 <= omega1 <= 0.64, 9.1 <= tau0 <= 15.7, 0. <= tau1 <= 9., \
+ 0. <= xvar0 <= 1.0, 0. <= xvar1 <= 1.2, -np.pi <= yvar0 <= np.pi, -np.pi <= yvar1 <= np.pi]):
return 0.0
return -np.inf
@@ -288,3 +288,25 @@ for yi in range(naxes):
ax.axhline(median[yi], color=mediancolor)
ax.plot(median[xi], median[yi], color = mediancolor, marker = 's')
figcorn.savefig(rootpath + '/plotsmc/0001_10M_interpolated_cornerplot_wandt_'+'nmax'+str(nmax)+'_tshift'+str(tshift)+'_'+str(nwalkers)+'walkers_'+str(npoints)+'pts.png', format='png', bbox_inches='tight', dpi=300)
+
+
+
+#Now plot the NR data against the mcmc fit data, together with the 1-sigma varying error data
+onesig_bounds = np.array([np.percentile(samples[:, i], [16, 84]) for i in range(len(samples[0]))]).T
+modelfitpk = model_dv(pk)
+figband = plt.figure(figsize = (12, 9))
+#Plot the 1-sigma_percentile
+for j in range(len(samples)):
+ sample = samples[j]
+ if np.all(onesig_bounds[0] <= sample) and np.all(sample <= onesig_bounds[1]):
+ plt.plot(timespan_new, model_dv(sample).real, "#79CAF2", alpha=0.3)
+
+plt.plot(timespan_new, gwdatanew_re, "k", alpha=0.7, lw=2, label=r'NR_re')
+plt.plot(timespan_new, modelfitpk.real, "r", alpha=0.7, lw=2, label=r'FitMCmax_re')
+plt.title(r'Comparison of the MC fit data and the $1-\sigma$ error band')
+plt.legend()
+plt.xlabel("t")
+plt.ylabel("h")
+
+figband.savefig(rootpath + '/plotsmc/0001_10M_interpolated_waveform_wandt_'+'nmax'+str(nmax)+'_tshift'+str(tshift)+'_'+str(nwalkers)+'walkers_'+str(npoints)+'pts.png', format = 'png', dpi = 384, bbox_inches = 'tight')
+
--
GitLab