Commit 019c1fd1 authored by Meg Millhouse's avatar Meg Millhouse
Browse files

update megaplot

git-svn-id: https://svn.ligo.caltech.edu/svn/bayeswave/tags/O1_online@511 c56465c9-8126-4a4f-9d7d-ac845eff4865
parent 823c0efd
......@@ -26,6 +26,8 @@ import os
import pwd
import subprocess
import sys
print sys.argv
from scipy import integrate
######################################################################################################################
#
......@@ -176,7 +178,7 @@ def readbwb():
print '\nThis run was executed with the --signalOnly flag\n'
restrictModel = 'signal'
elif arg=='--inj':
injFlag = True
injFlag = False #True
mdc = True
elif arg=='--MDC-cache':
mdc = True
......@@ -272,9 +274,40 @@ def get_mode(n, bins):
# Plotting functions
#
######################################################################################################################
#....................................................
# New function to get SNR^2(t) for x-axis
#....................................................
def snrfunctime(median_waveform):
powerList = []
wave_integral = median_waveform
time_integrate = time
#integral at each
dt = 0.0009765
h_t_2 = []
for line, line_1, time_i in zip(wave_integral, wave_integral[1:], time_integrate):
t = time_i
w0 = line
h_t_2.append(line**2)
w1 = line_1
snr_t = (((w0**2)+(w1**2))/2)*dt
powerList.append([snr_t,t])
#plot snr^2(t)
power, time_int = zip(*powerList)
plt.plot(time_int, power)
plt.title("SNR^2(t) for axes determination")
plt.xlabel('Time (s)')
plt.ylabel('{0} SNR^2'.format(model))
plt.title('IFO {0}'.format(ifo))
plt.axis([-1.0, 1.0, -0.1, np.max(power)+15])
plt.savefig(plotsDir+'{1}_SNR_TIME_{0}.png'.format(ifo, model))
return(powerList)
# -------------------------------------------------
# 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):
axisList = []
for ifo in ifoList:
......@@ -286,19 +319,32 @@ def get_axes(jobName, postDir, ifoList, worc, model, time):
filename = str(jobName)+postDir+'{1}_recovered_waveform.dat.{0}'.format(ifo, model)
median_waveform, up_vec, down_vec = get_waveform_bigfile(filename, ifo)
# -- Get axis info
snr_t = snrfunctime(median_waveform)
#y
wave = up_vec
wave_max = wave.max()
#x
wave = median_waveform
sig_times = time[wave > 0.05*wave.max()]
axisList.append([sig_times.min(), sig_times.max()+0.1*(sig_times.max()-sig_times.min()), -wave_max*1.1, wave_max*1.1])
power, and_time = zip(*snr_t)
power = list(power)
sig_times = []
#Fixed ~90$ rep. from 5% to 95%
for row in snr_t:
if (row[0] >= 0.001*np.max(power) and row[0] <= 0.999*np.max(power)):
sig_times.append(row[1])
axisList.append([np.min(sig_times), np.max(sig_times)+0.1*(np.max(sig_times)-np.min(sig_times)), -wave_max*1.1, wave_max*1.1])
# -- Select Axis with biggest y-value
ymax = 0
for axcand in axisList:
if axcand[3] > ymax:
ymax = axcand[3]
axwinner = axcand
axwinner= axcand
return(axwinner)
# --------------
# Plot Evidence
# --------------
......@@ -405,7 +451,7 @@ def plot_waveform(jobName, postDir, ifo, plotsDir, worc, mdc, model, axwinner, t
#axisList.append([sig_times.min(), sig_times.max(), -wave_max*1.1, wave_max*1.1])
plt.xlabel('Time (s)')
plt.ylabel('Best Fit {0} Waveform (blue) and injected waveform (red)'.format(model))
plt.title('IFO {0}'.format(ifo))
plt.title('{0}'.format(ifoNames[int(ifo)]))
# -- Save the full versions of the plot
if worc == 'whitened':
......@@ -442,7 +488,7 @@ def plot_cross_waveform(jobName, postDir, ifo, plotsDir, worc, mdc, model, axwin
#axisList.append([sig_times.min(), sig_times.max(), -wave_max*1.1, wave_max*1.1])
plt.xlabel('Time (s)')
plt.ylabel('Best Fit {0} Waveform (blue) and injected waveform (red)'.format(model))
plt.title('IFO {0}'.format(ifo))
plt.title('{0}'.format(ifoNames[int(ifo)]))
# -- Save the full versions of the plot
if worc == 'whitened':
......@@ -479,7 +525,7 @@ def plot_plus_waveform(jobName, postDir, ifo, plotsDir, worc, mdc, model, axwinn
#axisList.append([sig_times.min(), sig_times.max(), -wave_max*1.1, wave_max*1.1])
plt.xlabel('Time (s)')
plt.ylabel('Best Fit {0} Waveform (blue) and injected waveform (red)'.format(model))
plt.title('IFO {0}'.format(ifo))
plt.title('{0}'.format(ifoNames[int(ifo)]))
# -- Save the full versions of the plot
if worc == 'whitened':
......@@ -521,7 +567,7 @@ def bayesogram(filename, plotsDir, worc, time, twodrange, median_waveform, axwin
plt.plot(time, median_waveform, 'w', linewidth=1, alpha=0.5)
plt.xlabel('Time (s)')
plt.ylabel('h(t)')
plt.title('IFO {0} - {1} model'.format(ifo, model))
plt.title('{0} - {1} model'.format(ifoNames[int(ifo)],model))
plt.axis(axwinner)
if worc == 'whitened':
plt.savefig(plotsDir+'{1}_t_bgram_{0}.png'.format(ifo, model))
......@@ -565,7 +611,7 @@ def waveform_spectrogram(ifo, plotsDir, worc, time, median_waveform, model):
# Adjust time axis: trigger ends up at segment length-2.0s, so we must shift back to place the trigger at t=0s
tshift = np.abs(time[0])+binwidth
cax = ax.pcolormesh(bins-tshift, freqs, np.log10(spec_power), vmin=vmin, vmax=vmax, cmap='OrRd')
plt.title('IFO {0} - {1} model'.format(ifo, model))
plt.title('{0} - {1} model'.format(ifoNames[int(ifo)],model))
plt.xlabel('Time (s)')
plt.ylabel('Frequency (Hz)')
plt.axis([time[0], time[-1], 0, fs/2])
......@@ -625,7 +671,7 @@ def make_histogram(moments, moment, ifo, plotsDir, worc, outfile, xtype, ytype):
plt.legend()
plt.ylim(1, 2*rn.max())
if not mom == 'network_overlap':
plt.title('IFO {0}'.format(ifo))
plt.title('{0}'.format(ifoNames[int(ifo)]))
plt.xlabel(moments_dict[moment][0])
plt.grid()
if worc == 'whitened':
......@@ -736,27 +782,31 @@ def plot_likelihood_1(modelList, plotsDir):
# plt.xlim([1e-4, 1])
for mod in modelList:
if mod == 'glitch':
colour = 'blue'
elif mod == 'signal':
colour = 'green'
else:
colour = 'red'
colour = 'blue'
elif mod == 'signal':
colour = 'green'
else:
colour = 'red'
try:
try:
temp, likehood, error, acl = np.recfromtxt(str(jobName)+"{0}_evidence.dat".format(mod), unpack=True)
names = ['temp','likehood','error','acl']
data = np.recfromtxt(str(jobName)+"{0}_evidence.dat".format(mod), names=names)
except:
temp, likehood, error, acl = np.recfromtxt(str(jobName)+"{0}_1_evidence.dat".format(mod), unpack=True)
names = ['temp','likehood','error','acl']
data = np.recfromtxt(str(jobName)+"{0}_1_evidence.dat".format(mod), names=names)
except:
try:
try:
temp, likehood, error = np.recfromtxt(str(jobName)+"{0}_evidence.dat".format(mod), unpack=True)
names = ['temp','likehood','error']
data = np.recfromtxt(str(jobName)+"{0}_evidence.dat".format(mod), names=names)
except:
temp, likehood, error = np.recfromtxt(str(jobName)+"{0}_1_evidence.dat".format(mod), unpack=True)
except:
names = ['temp','likehood','error']
data = np.recfromtxt(str(jobName)+"{0}_1_evidence.dat".format(mod), names=names)
except:
continue
error = np.zeros(likehood.shape)
plt.semilogx(temp, likehood, label=mod, linewidth=2, color=colour)
plt.errorbar(temp, likehood, 2*error, color=colour)
#error = np.zeros(likehood.shape)
plt.semilogx(data['temp'], data['likehood'], label=mod, linewidth=2, color=colour)
plt.errorbar(data['temp'], data['likehood'], 2*data['error'], color=colour)
plt.xlabel( '1/Temp' )
plt.ylabel('log(L)')
plt.grid()
......@@ -769,27 +819,30 @@ def plot_likelihood_2(modelList, plotsDir):
plt.title('Likelihood')
for mod in modelList:
if mod == 'glitch':
colour = 'blue'
elif mod == 'signal':
colour = 'green'
else:
colour = 'red'
colour = 'blue'
elif mod == 'signal':
colour = 'green'
else:
colour = 'red'
try:
try:
temp, likehood, error, acl = np.recfromtxt(str(jobName)+"{0}_evidence.dat".format(mod), unpack=True)
names = ['temp','likehood','error','acl']
data = np.recfromtxt(str(jobName)+"{0}_evidence.dat".format(mod), names=names)
except:
temp, likehood, error, acl = np.recfromtxt(str(jobName)+"{0}_1_evidence.dat".format(mod), unpack=True)
names = ['temp','likehood','error','acl']
data = np.recfromtxt(str(jobName)+"{0}_1_evidence.dat".format(mod), names=names)
except:
try:
try:
temp, likehood, error = np.recfromtxt(str(jobName)+"{0}_evidence.dat".format(mod), unpack=True)
names = ['temp','likehood','error']
data = np.recfromtxt(str(jobName)+"{0}_evidence.dat".format(mod), names=names)
except:
temp, likehood, error = np.recfromtxt(str(jobName)+"{0}_1_evidence.dat".format(mod), unpack=True)
names = ['temp','likehood','error']
data = np.recfromtxt(str(jobName)+"{0}_1_evidence.dat".format(mod), names=names)
except:
continue
error = np.zeros(likehood.shape)
plt.semilogx(temp, likehood*temp, label=mod, linewidth=2, color=colour)
plt.errorbar(temp, likehood*temp, 2*error, color=colour)
plt.semilogx(data['temp'], data['likehood']*data['temp'], label=mod, linewidth=2, color=colour)
plt.errorbar(data['temp'], data['likehood']*data['temp'], 2*data['error'], color=colour)
plt.grid()
plt.xlabel('1/Temp')
plt.ylabel('log(L) X 1/Temp')
......@@ -873,10 +926,10 @@ 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='m', log=True)
n,bins,patches = ax1.hist(signalChains[2], bins=range(int(min(signalChains[2])), int(max(signalChains[2]))+1, 1), histtype='bar', color='m', log=False)
if not glitchChains == []:
data = np.dstack(glitchChains[3:3+len(ifoList)])[0]
n,bins,patches = ax2.hist(data, bins=range(int(data.min()), int(data.max())+1, 1), label=ifoNames, histtype='bar', log=True)
n,bins,patches = ax2.hist(data, bins=range(int(data.min()), int(data.max())+1, 1), label=ifoNames, histtype='bar', log=False)
# -- Legend placement outside the plot
box = ax1.get_position()
ax1.set_position([box.x0, box.y0, box.width*0.8, box.height])
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment