Skip to content
Snippets Groups Projects
Commit ff94e42f authored by Meg Millhouse's avatar Meg Millhouse
Browse files

Merge branch 'skip_noise_model' into 'master'

Skip gaussian noise model, and adjustments to output pages

See merge request lscsoft/bayeswave!208
parents 2d5ea0e9 5e2fde98
Branches master
No related tags found
No related merge requests found
......@@ -124,6 +124,7 @@ def readbwb():
restrictModel = ''
injFlag = False
lalsimFlag = False
noNoiseFlag = False
mdc = False
snrList = []
# -- Open *bayeswave.run
......@@ -157,12 +158,17 @@ def readbwb():
elif arg=='--glitchOnly':
print('\nThis run was executed with the --glitchOnly flag\n')
restrictModel = 'glitch'
noNoiseFlag = True
elif arg=='--noiseOnly':
print('\nThis run was executed with the --noiseOnly flag\n')
restrictModel = 'noise'
elif arg=='--signalOnly':
print('\nThis run was executed with the --signalOnly flag\n')
restrictModel = 'signal'
noNoiseFlag = True
elif arg=='--skipNoise':
print('\nThis run was executed with the --skipNoise flag\n')
noNoiseFlag = True
elif arg=='--inj':
injFlag = True
mdc = True
......@@ -193,7 +199,7 @@ def readbwb():
bayeswave.close()
# -- Report to user
print("{0}".format(info))
return(jobName, restrictModel, mdc, injFlag, bayeswaverunfile, ifoList, ifoNames, gps, snrList, info)
return(jobName, restrictModel, mdc, injFlag, noNoiseFlag, bayeswaverunfile, ifoList, ifoNames, gps, snrList, info)
# ---------------------------------------------
......@@ -263,7 +269,10 @@ def get_axes(jobName, postDir, ifoList, model, time, fullOnly_flag=0):
for row in snr_t:
if (row[0] >= 0.0001*np.max(power) and row[0] <= 0.9999*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])
try:
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])
except:
axisList.append([time[0],time[-1],-4.,4.])
# -- Select Axis with biggest y-value
......@@ -344,6 +353,7 @@ def plot_evidence(jobName, plotsDir):
err_sig_noise=0
err_sig_gl=0
err_sig_si=0
sig_noise = sig_si - sig_noise
sig_gl = sig_si - sig_gl
......@@ -352,9 +362,29 @@ def plot_evidence(jobName, plotsDir):
err_sig_noise = math.sqrt(err_sig_noise)
err_sig_gl += err_sig_si
err_sig_gl = math.sqrt(err_sig_gl)
# -- Report to the user
print(' log( E_signal / E_noise ) =', sig_noise)
print(' log( E_signal / E_glitch ) =', sig_gl)
# -- Account for case where noise model is skipped
if noNoiseFlag:
sig_noise = 0.0
err_sig_noise = 1.0
print(' log( E_signal / E_noise ) = N/A (no noise model)')
else:
# -- Report to the user
print(' log( E_signal / E_noise ) = {0}'.format(sig_noise))
if fullOnly_flag:
print(' log( E_signal / E_glitch ) = N/A (--fullOnly mode)')
elif restrictModel == 'signal':
sig_gl = 0.0
err_sig_gl = 1.0
print(' log( E_signal / E_glitch ) = N/A (no glitch model)')
elif restrictModel == 'glitch':
sig_gl = 0.0
err_sig_gl = 1.0
print(' log( E_signal / E_glitch ) = N/A (no signal model)')
else:
print(' log( E_signal / E_glitch ) = {0}'.format(sig_gl))
# -- Plot the data point
plt.figure()
plt.errorbar(sig_gl, sig_noise, 2*err_sig_gl, 2*err_sig_noise, color='black')
......@@ -773,24 +803,25 @@ def plot_likelihood_1(modelList, plotsDir, fullOnly_flag=0):
# Now do noise
mod = "noise"
if not noNoiseFlag:
mod = "noise"
colour = ncolor
colour = ncolor
try:
names = ['temp','likehood','error','acl']
data = np.recfromtxt(str(jobName)+"{0}_evidence.dat".format(mod), names=names)
except:
try:
names = ['temp','likehood','error']
names = ['temp','likehood','error','acl']
data = np.recfromtxt(str(jobName)+"{0}_evidence.dat".format(mod), names=names)
except:
print('No noise evidence')
#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)
minlist.append(data['likehood'][np.where(data['temp']>1e-4)[0][-1]])
maxlist.append(data['likehood'][0])
try:
names = ['temp','likehood','error']
data = np.recfromtxt(str(jobName)+"{0}_evidence.dat".format(mod), names=names)
except:
print('No noise evidence')
#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)
minlist.append(data['likehood'][np.where(data['temp']>1e-4)[0][-1]])
maxlist.append(data['likehood'][0])
plt.ylim(min(minlist)-1000,max(maxlist)+500)
......@@ -992,7 +1023,7 @@ def plot_model_dims(modelList, ifoList, ifoNames, plotsDir, fullOnly_flag=0):
# -- Calculate bayes factors for signal+glitch vs signal only, glitch only
if fullOnly_flag:
signal_on = signalChains > 0
glitch_on = np.sum(glitchChains[0:len(ifoList)+1],axis=1)>0
glitch_on = np.sum(glitchChains[0:len(ifoList)+1],axis=0)>0
ev_sig = 0
ev_glitch = 0
......@@ -1011,20 +1042,24 @@ def plot_model_dims(modelList, ifoList, ifoNames, plotsDir, fullOnly_flag=0):
# Actually full vs. signal only
try:
sig_noise = ev_full/ev_sig
sig_noise = float(ev_full)/float(ev_sig)
except:
sig_noise = np.inf
# Actually full vs. glitch only
try:
sig_gl = ev_full/ev_gl
sig_gl = float(ev_full)/float(ev_glitch)
except:
sig_gl = np.inf
if fullOnly_flag:
print("\nNumber of iterations in signal only, glitch only, and signal+glitch modes:\nNsig = {0}, Ngl = {1}, Nboth = {2}\n".format(ev_sig,ev_glitch,ev_full))
return sig_noise, sig_gl
######################################################################################################################
#
# BWB webpage production
......@@ -1161,28 +1196,53 @@ def splash_page(plotsDir, ifoNames, snrList, bayeswaverunfile, sig_gl, sig_noise
html_string = ''
if not fullOnly_flag:
html_string += ' <a href="./'+plotsDir+'odds.png"><img src="./'+plotsDir+'odds.png" style="float: right;" width=600 alt="odds.png"></a>\n'
else:
html_string += ' <h3>This run was done with the --fullOnly flag, meaning glitch and signal were run simultaneously.</h3>\n'
html_string += ' <h3>Detector Names: {0}</h3>\n'.format(', '.join(ifoNames))
html_string += ' <h3>Models used:</h3>\n'
html_string += ' <b>Simultaneous signal+glitch fitting:</b> '
if fullOnly_flag:
html_string += 'On</br>\n'
else:
html_string += 'Off</br>\n'
html_string += ' <b>Signal model:</b> '
if not restrictModel == 'glitch':
html_string += 'On</br>\n'
else:
html_string += 'Off</br>\n'
html_string += ' <b>Glitch model:</b> '
if not restrictModel == 'signal':
html_string += 'On</br>\n'
else:
html_string += 'Off</br>\n'
html_string += ' <b>Noise model:</b> '
if noNoiseFlag or fullOnly_flag:
html_string += 'Off</br>\n'
else:
html_string += 'On</br>\n'
if mdc:
html_string += ' <h3>Matched Filter SNRs of Injections</h3>\n'
for ifo, snr in zip(ifoNames, snrList):
html_string += ' {0} injected with SNR {1:.1f}<br/>\n'.format(ifo, snr)
html_string += ' <h3>Log Info</h3>\n'
html_string += ' <a href="./'+bayeswaverunfile+'">See full log file</a>\n'
if not restrictModel:
html_string +=' <h3>This is a {0} only run.</h3>\n'.format(restrictModel)
html_string +=' There are no Bayes factors for this run since it\'s only a single model.\n'
html_string +=' <h3>Evidence for Signal</h3>\n'
if not fullOnly_flag:
html_string +=' <h3>Evidence for Signal</h3>\n'
html_string +=' log(Evidence_signal / Evidence_glitch) = {0:.1f} &plusmn; {1:.1f} <br/>\n'.format(sig_gl,err_sig_gl)
html_string +=' log(Evidence_signal / Evidence_noise) = {0:.1f} &plusmn; {1:.1f} <br/>\n'.format(sig_noise,err_sig_noise)
if restrictModel == 'signal':
html_string += ' log(Evidence_signal / Evidence_glitch) = N/A (no glitch model)<br/>\n'
elif restrictModel == 'glitch':
html_string += ' log(Evidence_signal / Evidence_glitch) = N/A (no signal model)<br/>\n'
else:
html_string +=' log(Evidence_signal / Evidence_glitch) = {0:.1f} &plusmn; {1:.1f} <br/>\n'.format(sig_gl,err_sig_gl)
if not noNoiseFlag:
html_string +=' log(Evidence_signal / Evidence_noise) = {0:.1f} &plusmn; {1:.1f} <br/>\n'.format(sig_noise,err_sig_noise)
else:
html_string += ' log(Evidence_signal / Evidence_noise) = N/A (no noise model)\n'
else:
html_string +=' <h3>Bayes factor for signal+glitch vs signal only:</h3>\n'
html_string +=' {0} <br/>\n'.format(sig_noise)
html_string +=' <b>log(Bayes factor for signal+glitch vs signal only):</b>\n'
html_string +=' {0} <br/>\n'.format(np.log(sig_noise))
html_string +=' <h3>Bayes factor for signal+glitch vs glitch only:</h3>\n'
html_string +=' {0} <br/>\n'.format(sig_gl)
html_string +=' <b>log(Bayes factor for signal+glitch vs glitch only):</b>\n'
html_string +=' {0} <br/>\n'.format(np.log(sig_gl))
html_string +=' </p>\n'
return html_string
......@@ -1503,7 +1563,7 @@ parser = argparse.ArgumentParser(description='Produce html page for a bayeswave
# -- Get basic info on the run
jobName, restrictModel, mdc, injFlag, bayeswaverunfile, ifoList, ifoNames, gps, snrList, info = readbwb()
jobName, restrictModel, mdc, injFlag, noNoiseFlag, bayeswaverunfile, ifoList, ifoNames, gps, snrList, info = readbwb()
if not restrictModel == '':
model = restrictModel #TODO: think about the --noiseOnly case
......
......@@ -866,6 +866,11 @@ void print_run_flags(FILE *fptr, struct Data *data, struct Prior *prior)
else if(!data->glitchFlag) fprintf(fptr, "REQUIRED");
else fprintf(fptr, "ENABLED");
fprintf(fptr, "\n");
fprintf(fptr, " noise model is ........... ");
if(!data->noiseModelFlag) fprintf(fptr, "DISABLED");
else fprintf(fptr, "ENABLED");
fprintf(fptr, "\n");
fprintf(fptr, " wavelet prior is ......... ");
if(data->waveletPriorFlag) fprintf(fptr, "ENABLED");
......@@ -1457,6 +1462,15 @@ void parse_command_line(struct Data *data, struct Chain *chain, struct Prior *pr
data->fullModelFlag = 0;
}
ppt = LALInferenceGetProcParamVal(commandLine, "--skipNoise");
if(ppt)
{
data->noiseModelFlag = 0;
data->glitchModelFlag = 1;
data->signalModelFlag = 1;
data->fullModelFlag = 0;
}
ppt = LALInferenceGetProcParamVal(commandLine, "--glitchOnly");
if(ppt)
{
......@@ -1483,7 +1497,7 @@ void parse_command_line(struct Data *data, struct Chain *chain, struct Prior *pr
data->signalModelFlag = 0;
data->fullModelFlag = 1;
}
ppt = LALInferenceGetProcParamVal(commandLine, "--stochastic");
if(ppt)
{
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment