diff --git a/BayesWaveUtils/scripts/megaplot.py b/BayesWaveUtils/scripts/megaplot.py
index f3efd4b7f680f2e03f1c81c1a0d542ecac81c9fd..5e0af486724564bd968a3ca76bea0a930d8444b1 100755
--- a/BayesWaveUtils/scripts/megaplot.py
+++ b/BayesWaveUtils/scripts/megaplot.py
@@ -119,8 +119,10 @@ htmlDir = 'html/'
#### --fullOnly ###
fullOnly = 0
-if os.path.exists('post/full'):
+if os.path.exists(postDir+'full'):
fullOnly_flag = 1
+ print("Found --fullOnly flag\n")
+
#Adopt common color scheme for different models
ncolor = 'darkgrey'
@@ -331,16 +333,18 @@ def snrfunctime(median_waveform):
# Read in data and median waveform to get plot axis
# Now using SNR^2 for time axis determination
# -------------------------------------------------
-def get_axes(jobName, postDir, ifoList, worc, model, time):
+def get_axes(jobName, postDir, ifoList, worc, model, time, fullOnly_flag=0):
axisList = []
for ifo in ifoList:
# -- Read Signal model
+ if fullOnly_flag==1:
+ filename = str(jobName)+postDir+'full/{1}_median_time_domain_waveform_{0}.dat'.format(ifoNames[int(ifo)], model)
+ else:
+ filename = str(jobName)+postDir+'{1}/{1}_median_time_domain_waveform_{0}.dat'.format(ifoNames[int(ifo)], model)
try:
- filename = str(jobName)+postDir+'{1}/{1}_recovered_{2}_waveform_{0}.dat'.format(ifoNames[int(ifo)], model, worc)
- median_waveform, up_vec, down_vec = get_waveform_bigfile(filename)
+ timesamp, median_waveform, dummy1, dummy2, up_vec, down_vec = get_waveform(filename)
except:
- filename = str(jobName)+postDir+'{1}/{1}_median_time_domain_waveform_{0}.dat'.format(ifoNames[int(ifo)], model)
- timesamp, median_waveform, up_vec, down_vec = get_waveform(filename)
+ print("Couldn't find waveform file %s" % filename)
# filename = str(jobName)+postDir+'{1}_median_time_domain_waveform_{0}.dat'.format(ifoNames[int(ifo)], model)
# timesamp, median_waveform, up_vec, down_vec, up_vec2, down_vec2 = get_waveform(filename)
@@ -838,8 +842,13 @@ def bayesogram(filename, plotsDir, worc, time, twodrange, median_waveform, axwin
# --------------------------------------
# Plot Q scans of waveforms
# --------------------------------------
-def Q_scan(subDir,model, Q, t, f, ifo, axwinner, f_axwinner, climwinner=[1,1]):
- data = np.genfromtxt(postDir+subDir+'/'+'{0}_spectrogram_{1}_{2}.dat'.format(model, Q, ifoNames[int(ifo)]))
+def Q_scan(subDir,model, Q, t, f, ifo, axwinner, f_axwinner, climwinner=[1,1], fullOnly_flag=0):
+
+ if fullOnly_flag==1:
+ filename = postDir+'/full/{0}_spectrogram_{1}_{2}.dat'.format(model, Q, ifoNames[int(ifo)])
+ else:
+ filename = postDir+subDir+'/'+'{0}_spectrogram_{1}_{2}.dat'.format(model, Q, ifoNames[int(ifo)])
+ data = np.genfromtxt(filename)
#plt.clim(np.min(data)/np.sum(data),np.max(data)/np.sum(data))
@@ -1192,31 +1201,38 @@ def plot_psd(jobName, postDir, ifoList, ifoNames, plotsDir):
-def plot_model_dims(modelList, ifoList, ifoNames, plotsDir):
+def plot_model_dims(modelList, ifoList, ifoNames, plotsDir, fullOnly_flag=0):
lineStyles = ['-', '--', ':']
lineColors = [H1color, L1color, V1color]
glitchChains = []
signalChains = []
# -- Read in data
- for mod in modelList:
+ if fullOnly_flag==1:
+ intChainFile = 'chains/'+str(jobName)+"full_model.dat.0"
+ N = len(open(intChainFile).readlines())
+ chains = np.loadtxt(intChainFile,skiprows=N/2)
+ chains = np.transpose(chains)
+
+ else:
+ for mod in modelList:
- intChainFile = 'chains/'+str(jobName)+"{0}_model.dat.0".format(mod)
- N = len(open(intChainFile).readlines())
- chains = np.loadtxt(intChainFile,skiprows=N/2)
- chains = np.transpose(chains)
- if mod == 'glitch':
- glitchChains = chains
- elif mod == 'signal':
- signalChains = chains
+ intChainFile = 'chains/'+str(jobName)+"{0}_model.dat.0".format(mod)
+ N = len(open(intChainFile).readlines())
+ chains = np.loadtxt(intChainFile,skiprows=N/2)
+ chains = np.transpose(chains)
+ if mod == 'glitch':
+ glitchChains = chains[3:3+len(ifoList)]
+ elif mod == 'signal':
+ signalChains = chains[2]
- signalOnly_flag = 0
- glitchOnly_flag = 0
+ signalOnly_flag = 0
+ glitchOnly_flag = 0
- if signalChains == []:
- signalOnly_flag == 1
- if glitchChains == []:
- glitchOnly_flag == 1
+ if signalChains == []:
+ signalOnly_flag == 1
+ if glitchChains == []:
+ glitchOnly_flag == 1
# -- Variation of the model dimension over MCMC iterations. 2 subplots for signal and glitch models
@@ -1224,11 +1240,11 @@ def plot_model_dims(modelList, ifoList, ifoNames, plotsDir):
fig, (ax1, ax2) = plt.subplots(2, sharex=True, sharey=True)
ax1.set_title('Model Dimension')
if not signalChains == []:
- ax1.plot(signalChains[2], linewidth=2, color=scolor, label='Signal')
+ ax1.plot(signalChains, linewidth=2, color=scolor, label='Signal')
if not glitchChains == []:
for ifo in ifoList:
ifoint=int(ifo)
- ax2.plot(glitchChains[3+ifoint], lineStyles[ifoint], color=lineColors[ifoint], linewidth=2, label=ifoNames[ifoint])
+ ax2.plot(glitchChains[ifoint], lineStyles[ifoint], color=lineColors[ifoint], linewidth=2, label=ifoNames[ifoint])
# -- Legend placement outside the plot
box = ax1.get_position()
@@ -1252,11 +1268,11 @@ def plot_model_dims(modelList, ifoList, ifoNames, plotsDir):
ax1 = plt.subplot(111)
ax1.set_title('Model Dimension')
if not signalChains == []:
- ax1.plot(signalChains[2], linewidth=2, color=scolor, label='Signal')
+ ax1.plot(signalChains, linewidth=2, color=scolor, label='Signal')
if not glitchChains == []:
for ifo in ifoList:
ifoint=int(ifo)
- ax1.plot(glitchChains[3+ifoint], lineStyles[ifoint], color=lineColors[ifoint], linewidth=2, label=ifoNames[ifoint])
+ ax1.plot(glitchChains[ifoint], lineStyles[ifoint], color=lineColors[ifoint], linewidth=2, label=ifoNames[ifoint])
# -- Legend placement outside the plot
box = ax1.get_position()
@@ -1281,9 +1297,9 @@ def plot_model_dims(modelList, ifoList, ifoNames, plotsDir):
fig, (ax1, ax2) = plt.subplots(2, sharex=True)
ax1.set_title('Model Dimension Histogram')
if not signalChains == []:
- n,bins,patches = ax1.hist(signalChains[2], bins=range(int(min(signalChains[2])), int(max(signalChains[2]))+1, 1), histtype='bar', color=scolor, log=False)
+ n,bins,patches = ax1.hist(signalChains, bins=range(int(min(signalChains[2])), int(max(signalChains[2]))+1, 1), histtype='bar', color=scolor, log=False)
if not glitchChains == []:
- data = np.dstack(glitchChains[3:3+len(ifoList)])[0]
+ data = np.dstack(glitchChains)[0]
#ax2.set_prop_cycle(cycler('color',['darkgoldenrod','darkkhaki','dkarsage']))
n,bins,patches = ax2.hist(data, bins=range(int(data.min()), int(data.max())+1, 1), label=ifoNames, histtype='bar', log=False, color=[lineColors[int(i)] for i in ifoList])
# -- Legend placement outside the plot
@@ -1308,9 +1324,9 @@ def plot_model_dims(modelList, ifoList, ifoNames, plotsDir):
ax1 = plt.subplot(111)
ax1.set_title('Model Dimension Histogram')
if not signalChains == []:
- n,bins,patches = ax1.hist(signalChains[2], bins=range(int(min(signalChains[2])), int(max(signalChains[2]))+1, 1), histtype='bar', color=scolor, log=False)
+ n,bins,patches = ax1.hist(signalChains, bins=range(int(min(signalChains[2])), int(max(signalChains[2]))+1, 1), histtype='bar', color=scolor, log=False)
if not glitchChains == []:
- data = np.dstack(glitchChains[3:3+len(ifoList)])[0]
+ data = np.dstack(glitchChains)[0]
#ax2.set_prop_cycle(cycler('color',['darkgoldenrod','darkkhaki','dkarsage']))
n,bins,patches = ax1.hist(data, bins=range(int(data.min()), int(data.max())+1, 1), label=ifoNames, histtype='bar', log=False, color=[lineColors[int(i)] for i in ifoList])
# -- Legend placement outside the plot
@@ -1393,7 +1409,7 @@ def make_mode_table(page, subpage, imdc, momentsPlus, momentsCross, ifoList, tab
subpage.write(' </table>\n')
-def whitened_residual_plots(model,ifoList,ifoNames):
+def whitened_residual_plots(model,ifoList,ifoNames,fullOnly_flag=0):
if model == 'glitch':
colour = gcolor
@@ -1418,8 +1434,11 @@ def whitened_residual_plots(model,ifoList,ifoNames):
psd_info = {}
for ifo in ifoList:
-
- filename = 'post/{1}/{1}_median_PSD_{0}.dat'.format(ifoNames[int(ifo)],model)
+
+ if fullOnly_flag==1:
+ filename = 'post/full/{1}_median_PSD_{0}.dat'.format(ifoNames[int(ifo)],model)
+ else:
+ filename = 'post/{1}/{1}_median_PSD_{0}.dat'.format(ifoNames[int(ifo)],model)
psd_info[ifo] = get_waveform(filename)
imin = 0
@@ -2107,6 +2126,9 @@ def make_webpage(htmlDir, model, mdc, gps, ifoList, ifoNames, modelList, snrList
#
######################################################################################################################
# -- Parse command line arguments
+
+print("Running version compatable with --fullOnly flag\n\n")
+
parser = argparse.ArgumentParser(description='Produce html page for a bayeswave run.')
@@ -2151,7 +2173,7 @@ for mod in modelList:
for worc in ['whitened']:
# -- Determine axes for waveform plots
- axwinner = get_axes(jobName, postDir, ifoList, worc, model, time)
+ axwinner = get_axes(jobName, postDir, ifoList, worc, model, time, fullOnly_flag)
# -- Read in the reconstructed waveform
try:
@@ -2171,11 +2193,16 @@ for mod in modelList:
# -- Plot the median waveform and injected waveform (if present)
plot_waveform(jobName, postDir, ifo, plotsDir, worc, mdc, mod, axwinner, time, low_50, high_50, low_90, high_90, median_waveform)
-
- filename = str(jobName)+postDir+'{1}/{1}_median_frequency_domain_waveform_spectrum_{0}.dat'.format(ifoNames[int(ifo)],mod)
+ if fullOnly_flag==1:
+ filename = str(jobName)+postDir+'full/{1}_median_frequency_domain_waveform_spectrum_{0}.dat'.format(ifoNames[int(ifo)],mod)
+ else:
+ filename = str(jobName)+postDir+'{1}/{1}_median_frequency_domain_waveform_spectrum_{0}.dat'.format(ifoNames[int(ifo)],mod)
powerspec_info = get_waveform(filename)
- filename = str(jobName)+postDir+'{1}/{1}_median_PSD_{0}.dat'.format(ifoNames[int(ifo)],mod)
+ if fullOnly_flag==1:
+ filename = str(jobName)+postDir+'full/{1}_median_PSD_{0}.dat'.format(ifoNames[int(ifo)],mod)
+ else:
+ filename = str(jobName)+postDir+'{1}/{1}_median_PSD_{0}.dat'.format(ifoNames[int(ifo)],mod)
psd_info = get_waveform(filename)
f_axwinner = get_axes_fdom(jobName, postDir, ifoList, worc, model, time)
@@ -2185,7 +2212,10 @@ for mod in modelList:
# ---- tf track
- filename = str(jobName)+postDir+'{1}/{1}_median_tf_{0}.dat'.format(ifoNames[int(ifo)],mod)
+ if fullOnly_flag==1:
+ filename = str(jobName)+postDir+'full/{1}_median_tf_{0}.dat'.format(ifoNames[int(ifo)],mod)
+ else:
+ filename = str(jobName)+postDir+'{1}/{1}_median_tf_{0}.dat'.format(ifoNames[int(ifo)],mod)
freqsamp, median_powerspec, high_50, low_50, high_90, low_90 = get_waveform(filename)
#plot_tf_tracks(jobName, postDir, ifo, plotsDir, worc, mdc, mod, axwinner, f_axwinner, freqsamp, low_50, high_50, low_90, high_90, median_powerspec);
@@ -2203,18 +2233,18 @@ for mod in modelList:
climwinner8 = [1,1]
climwinner16 = [1,1]
- Q_scan(mod,mod, 4, time, freq, ifo, axwinner, f_axwinner)
- Q_scan(mod,mod, 8, time, freq, ifo, axwinner, f_axwinner)
- Q_scan(mod,mod, 16, time, freq, ifo, axwinner, f_axwinner)
+ Q_scan(mod,mod, 4, time, freq, ifo, axwinner, f_axwinner, fullOnly_flag=fullOnly_flag)
+ Q_scan(mod,mod, 8, time, freq, ifo, axwinner, f_axwinner, fullOnly_flag=fullOnly_flag)
+ Q_scan(mod,mod, 16, time, freq, ifo, axwinner, f_axwinner, fullOnly_flag=fullOnly_flag)
climwinner4 = Q_scan('data','data', 4, time, freq, ifo, axwinner, f_axwinner)
climwinner8 = Q_scan('data','data', 8, time, freq, ifo, axwinner, f_axwinner)
climwinner16 = Q_scan('data','data', 16, time, freq, ifo, axwinner, f_axwinner)
- Q_scan(mod,'{0}_residual'.format(mod), 4, time, freq, ifo, axwinner, f_axwinner, climwinner4)
- Q_scan(mod,'{0}_residual'.format(mod), 8, time, freq, ifo, axwinner, f_axwinner, climwinner8)
- Q_scan(mod,'{0}_residual'.format(mod), 16, time, freq, ifo, axwinner, f_axwinner, climwinner16)
+ Q_scan(mod,'{0}_residual'.format(mod), 4, time, freq, ifo, axwinner, f_axwinner, climwinner4, fullOnly_flag)
+ Q_scan(mod,'{0}_residual'.format(mod), 8, time, freq, ifo, axwinner, f_axwinner, climwinner8, fullOnly_flag)
+ Q_scan(mod,'{0}_residual'.format(mod), 16, time, freq, ifo, axwinner, f_axwinner, climwinner16, fullOnly_flag)
if mdc:
@@ -2228,7 +2258,11 @@ for mod in modelList:
# ---------------------------------
for worc in ['whitened', 'colored']:
# -- Read moments file and open mode output file
- moments = np.recfromtxt(str(jobName)+postDir+'{1}/{1}_{2}_moments_{0}.dat'.format(ifoNames[int(ifo)], mod, worc), names=True)
+ if fullOnly_flag==1:
+ momentsfile = str(jobName)+postDir+'full/{1}_{2}_moments_{0}.dat'.format(ifoNames[int(ifo)], mod, worc)
+ else:
+ momentsfile = str(jobName)+postDir+'{1}/{1}_{2}_moments_{0}.dat'.format(ifoNames[int(ifo)], mod, worc)
+ moments = np.recfromtxt(momentsfile, names=True)
momentsPlus = [] # Moments relative to the + polarization only
momentsCross = [] # Moments relative to the x polarization only
if mdc:
@@ -2343,18 +2377,18 @@ except:
if 'signal' in modelList or 'glitch' in modelList:
- plot_model_dims(modelList, ifoList, ifoNames, plotsDir)
+ plot_model_dims(modelList, ifoList, ifoNames, plotsDir, fullOnly_flag)
# -- Posterior predictive checks:
for mod in modelList:
try:
- whitened_residual_plots(mod,ifoList,ifoNames)
+ whitened_residual_plots(mod,ifoList,ifoNames,fullOnly_flag)
except Exception as e:
print("Error producing %s residual histograms"%mod)
print(traceback.format_exc())
try:
- whitened_residual_plots('noise',ifoList,ifoNames)
+ whitened_residual_plots('noise',ifoList,ifoNames,fullOnly_flag)
except Exception as e:
print("Error producing noise residual histograms")
print(traceback.format_exc())