From ccd68a81338590d1db5864bead1b6cc3999828f1 Mon Sep 17 00:00:00 2001 From: Sophie Hourihane Date: Wed, 4 Aug 2021 14:57:43 -0700 Subject: [PATCH 01/39] Adding cleanwindow parameter to only subtract glitches at the edge --- src/BayesWave.c | 29 ++++++++++++++++++++++++----- src/BayesWave.h | 1 + src/BayesWaveIO.c | 4 ++++ 3 files changed, 29 insertions(+), 5 deletions(-) diff --git a/src/BayesWave.c b/src/BayesWave.c index 0b4e5db..be82654 100644 --- a/src/BayesWave.c +++ b/src/BayesWave.c @@ -746,15 +746,34 @@ int main(int argc, char *argv[]) g = m->glitch[ifo]; for(i=1; i<=g->size; i++) { - if(g->intParams[i][0]range[0][0] || g->intParams[i][0]>prior->range[0][1]) + // if cleanwindow > -1 , + // then we only subtract glitches from (segment_start, segment_start + cleanwindow) + // and (segment_end - cleanwindow, segment_end) + // Otherwise subtract glitches outside of window parameter + if (data->cleanwindow == -1) { - //regress from the data & glitch model se we can compute window SNR - m->wavelet(data->s[ifo], g->intParams[i], N, -1, data->Tobs); - m->wavelet(g->templates, g->intParams[i], N, -1, data->Tobs); + // subtracts glitches outside of the window parameter + if(g->intParams[i][0]range[0][0] || g->intParams[i][0]>prior->range[0][1]) + { + //regress from the data & glitch model se we can compute window SNR + m->wavelet(data->s[ifo], g->intParams[i], N, -1, data->Tobs); + m->wavelet(g->templates, g->intParams[i], N, -1, data->Tobs); + } } - //regress all wavelets from residual array for plotting + else + { + // subtracts glitches outside of the cleanwindow parameter + if(g->intParams[i][0]< data->cleanwindow || g->intParams[i][0]>data->Tobs - data->cleanwindow) + { + //regress from the data & glitch model se we can compute window SNR + m->wavelet(data->s[ifo], g->intParams[i], N, -1, data->Tobs); + m->wavelet(g->templates, g->intParams[i], N, -1, data->Tobs); + } + } + //regress all wavelets from residual array for plotting m->wavelet(r[ifo], g->intParams[i], N, -1, data->Tobs); } + ifoSNR[ifo] += fourier_nwip(data->imin,data->imax,g->templates,g->templates,m->invSnf[ifo]); rSNR += ifoSNR[ifo]; fprintf(stdout, " Residual SNR in %s = %g\n",data->ifos[ifo],sqrt(ifoSNR[ifo])); diff --git a/src/BayesWave.h b/src/BayesWave.h index f9d7511..0ea6540 100644 --- a/src/BayesWave.h +++ b/src/BayesWave.h @@ -355,6 +355,7 @@ struct Data double fmax; // maximum frequency for inner products double Tobs; // segment length double Twin; // window length + double cleanwindow; // time at either edge of the segment that the cleaning phase subtracts double Qmin; // minimum wavelet quality factor double Qmax; // maximum wavelet quality factor double gmst; // Greenwich time diff --git a/src/BayesWaveIO.c b/src/BayesWaveIO.c index cafae48..53c828a 100644 --- a/src/BayesWaveIO.c +++ b/src/BayesWaveIO.c @@ -1540,6 +1540,10 @@ void parse_command_line(struct Data *data, struct Chain *chain, struct Prior *pr if(ppt) data->Twin = atof(ppt->value); else data->Twin = 1.0; + ppt = LALInferenceGetProcParamVal(commandLine, "--cleanwindow"); + if(ppt) data->cleanwindow = atof(ppt->value); + else data->cleanwindow = -1; + ppt = LALInferenceGetProcParamVal(commandLine, "--Nchain"); if(ppt) chain->NC = atoi(ppt->value); else chain->NC = 20; -- GitLab From fbcff50b0560c8d56464c2f53eb028e53e36edaa Mon Sep 17 00:00:00 2001 From: Sophie Hourihane Date: Fri, 13 Aug 2021 16:06:47 -0700 Subject: [PATCH 02/39] Transferring the trigdir to the remote node so that if the run fails after checkpointing, the dag can be resubmitted without restarting from the very beginning --- BayesWaveUtils/bayeswave_pipe/bayeswave_pipe_utils.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/BayesWaveUtils/bayeswave_pipe/bayeswave_pipe_utils.py b/BayesWaveUtils/bayeswave_pipe/bayeswave_pipe_utils.py index 87b8e45..8f21572 100644 --- a/BayesWaveUtils/bayeswave_pipe/bayeswave_pipe_utils.py +++ b/BayesWaveUtils/bayeswave_pipe/bayeswave_pipe_utils.py @@ -1145,7 +1145,8 @@ class bayeswaveJob(pipeline.CondorDAGJob, pipeline.AnalysisJob): self.add_condor_cmd("+PreArgs", '"$(macrooutputDir) bayeswave"') # Configure file transfers - transferstring='datafind,setupdirs.py' + transferstring='$(macrooutputDir),datafind,setupdirs.py' + print("Transfer input files", transferstring) if cp.getboolean('condor','copy-frames'): transferstring += ',$(macroframes)' -- GitLab From 3aab31da5415dd7bc3b9649a7680d6cd62dd4354 Mon Sep 17 00:00:00 2001 From: Sophie Hourihane Date: Tue, 17 Aug 2021 23:21:55 +0000 Subject: [PATCH 03/39] Revert "Transferring the trigdir to the remote node so that if the run fails after..." This reverts commit fbcff50b0560c8d56464c2f53eb028e53e36edaa --- BayesWaveUtils/bayeswave_pipe/bayeswave_pipe_utils.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/BayesWaveUtils/bayeswave_pipe/bayeswave_pipe_utils.py b/BayesWaveUtils/bayeswave_pipe/bayeswave_pipe_utils.py index 8f21572..87b8e45 100644 --- a/BayesWaveUtils/bayeswave_pipe/bayeswave_pipe_utils.py +++ b/BayesWaveUtils/bayeswave_pipe/bayeswave_pipe_utils.py @@ -1145,8 +1145,7 @@ class bayeswaveJob(pipeline.CondorDAGJob, pipeline.AnalysisJob): self.add_condor_cmd("+PreArgs", '"$(macrooutputDir) bayeswave"') # Configure file transfers - transferstring='$(macrooutputDir),datafind,setupdirs.py' - print("Transfer input files", transferstring) + transferstring='datafind,setupdirs.py' if cp.getboolean('condor','copy-frames'): transferstring += ',$(macroframes)' -- GitLab From 9cd2f60bea39f7ad05930cc3fc8faef18fb71006 Mon Sep 17 00:00:00 2001 From: "cailin.plunkett" Date: Mon, 26 Jul 2021 10:50:56 -0700 Subject: [PATCH 04/39] megaplot fix signalOnly bug --- BayesWaveUtils/scripts/megaplot.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/BayesWaveUtils/scripts/megaplot.py b/BayesWaveUtils/scripts/megaplot.py index b04e5ca..9ec8519 100755 --- a/BayesWaveUtils/scripts/megaplot.py +++ b/BayesWaveUtils/scripts/megaplot.py @@ -2614,7 +2614,7 @@ for mod in RUN_FLAGS.modelList: plot_waveform(jobName, postDir, ifo, plotsDir, worc, mdc, mod, axwinner, time, low_50, high_50, low_90, high_90, median_waveform) # make corner plots for signal or glitch data - if mod in ('signal', 'glitch'): + if mod in ('signal', 'glitch') and ('glitch' in RUN_FLAGS.modelList): if RUN_FLAGS.CBCType: # makes a glitch parameter plot filename = 'chains/cbc_params_{0}.dat.0'.format(ifoNames[int(ifo)], mod) -- GitLab From 4c6ef66948d179777c387ac9e3c56a98bd5ea46b Mon Sep 17 00:00:00 2001 From: "cailin.plunkett" Date: Tue, 17 Aug 2021 17:43:40 -0700 Subject: [PATCH 05/39] Print smooth PSDs from arbitary IFO list --- src/BayesWave.c | 19 +++++++------------ 1 file changed, 7 insertions(+), 12 deletions(-) diff --git a/src/BayesWave.c b/src/BayesWave.c index be82654..3bfac53 100644 --- a/src/BayesWave.c +++ b/src/BayesWave.c @@ -978,19 +978,14 @@ int main(int argc, char *argv[]) // Dump smooth psd for comparison debug if (bayescbc->debug) { - sprintf(filename,"smooth_psd_H1.dat"); - chain->runFile = fopen(filename,"w"); - for(i=1; iN/2; i++) { - fprintf(chain->runFile, "%f, %e, %e\n", i/data->Tobs, model[chain->index[0]]->Snf[0][i], model[chain->index[0]]->SnS[0][i]); - } - fclose(chain->runFile); - - sprintf(filename,"smooth_psd_L1.dat"); - chain->runFile = fopen(filename,"w"); - for(i=0; iN/2; i++) { - fprintf(chain->runFile, "%f, %e, %e\n", i/data->Tobs, model[chain->index[0]]->Snf[1][i], model[chain->index[0]]->SnS[1][i]); + for( ifo = 0; ifo < data->NI; ifo++ ){ + sprintf(filename,"smooth_psd_%s.dat", data->ifos[ifo]); + chain->runFile = fopen(filename,"w"); + for(i=0; iN/2; i++) { + fprintf(chain->runFile, "%f, %e, %e\n", i/data->Tobs, model[chain->index[0]]->Snf[ifo][i], model[chain->index[0]]->SnS[ifo][i]); + } + fclose(chain->runFile); } - fclose(chain->runFile); } } -- GitLab From f47f01c6ed671f5a6d17df2cd9a6274bad324ff0 Mon Sep 17 00:00:00 2001 From: "cailin.plunkett" Date: Tue, 17 Aug 2021 17:45:03 -0700 Subject: [PATCH 06/39] Making BW_Flags an executable and adding cleaning phase output page --- BayesWaveUtils/scripts/BW_Flags.py | 382 +++++++++++++++++++++++++++++ BayesWaveUtils/scripts/megaplot.py | 306 +++++++++-------------- BayesWaveUtils/setup.py | 2 +- 3 files changed, 496 insertions(+), 194 deletions(-) create mode 100644 BayesWaveUtils/scripts/BW_Flags.py diff --git a/BayesWaveUtils/scripts/BW_Flags.py b/BayesWaveUtils/scripts/BW_Flags.py new file mode 100644 index 0000000..39d4739 --- /dev/null +++ b/BayesWaveUtils/scripts/BW_Flags.py @@ -0,0 +1,382 @@ +#!/usr/bin/env python +from bayeswave_pipe import bayeswave_pipe_utils as pipe_utils +import numpy as np +import os +import sys + + +""" + Class containing run information for a BayesWave run + The parameters are kept inside of trigdir / bayeswave.run + + Must use readbwb to actually initialize values + + Example use: + # Initialize runFlags + runFlags = Flags('/home/sophie.hourihane/public_html/glitch_events/noNoise/seeds_noClean/0noise_2048_5chains_heterodyne_distance600') + # Set runFlags + +""" +class Flags: + def __init__(self, trigdir): + # Parameters describing where things live + self.trigdir = trigdir + self.XMLfname = None + self.jobName='' + + # Run settings + self.trigtime = 0 # gps time of event + self.flow = 32 + self.srate = 2048 + self.segment_start = 0 + self.seglen = 0 + self.Nchain = 20 # 20 is default number of chains + + # Windows, cleaning, and bayesLine + self.bayesLine = False + self.noClean = 0 + self.CBCwindow = -1 #window around cbc event + self.window = 0 + + # Parameters describing what run type was used + self.fullOnly_flag = 0 + self.GlitchCBC_flag = 0 + self.SignalCBC_flag = 0 + self.glitchOnly_flag = 0 + self.CBCType = 0 # ie, literally if CBC is in the modellist + self.CBCName = 0 # I think this is unused now + self.multi_type = 0 # true for full, cbcsignal, cbcglitch + self.noNoiseFlag = False + self.cleanOnly_flag = 0 + + self.modelList = () # Gets filled in readbwb, theres probably a clearer way to to this + + # Injection settings + self.injFlag = False + self.event = 0 + self.mdc = False + self.inj_dat = None #If there is an injection, this holds an injection table + + # Likelihood calculation + self.heterodyneL = False + + # Printed output + self.verbose = 0 + self.checkpoint = False + + # Lists of IFOs + self.ifoList = [] + self.ifoNames = [] + self.snrList = [] + + def run_dirname(self): + if self.fullOnly_flag: + return('full') + if self.SignalCBC_flag: + return('cbcsignal') + if self.GlitchCBC_flag: + return('cbcglitch') + if len(self.modelList) == 1: + return(self.modelList[0]) + else: + print("uh oh! what's the dirname?") + return('') + + """ + Given an xml parameter name, will return that value + Accepted names are: + param_names = ['M1_s', 'M2_s', 'Mchirp_s', 'Mtot_s', + 's1', 's2', 'chi_eff', + 'phase', 't', + 'distance', 'redshift', + 'RA', 'sindec', 'psi', 'cosiota', + 'M1_d', 'M2_d', 'Mchirp_d', 'Mtot_d'] + """ + def get_xml_parameter(self, p): + if p[-2:] == '_d': + p = p[:-2] + # if we get a value in the source frame, divide by (1 + z) and return the detector frame + if p[-2:] == '_s': + p = p[:-2] + z = redshift_from_distance(self.inj_dat.getColumnByName('distance')[self.event]) + return(self.get_xml_parameter(p) / (1.0 + z)) + if p == 'Mc' or p == 'Mchirp': + return(self.inj_dat.getColumnByName('mchirp')[self.event]) + elif p == 'Z' or p == 'redshift': + return(redshift_from_distance(self.inj_dat.getColumnByName('distance')[self.event])) + elif p == 'spin1' or p == 'spin2': + print("TODO: get_xml_paramter spin1 and spin2, returning 0") + return 0 + elif (p == 'chi_eff'): + print("FIX get_xml_parameter chi_eff and spin parameters, returning 0") + return 0 + elif p == 'Q': + return(self.inj_dat.getColumnByName('mass2')[self.event] / self.inj_dat.getColumnByName('mass1')[self.event]) + elif p == 'mtot' or p == 'Mt': + return(self.inj_dat.getColumnByName('mass2')[self.event] + self.inj_dat.getColumnByName('mass1')[self.event]) + elif p == 'coa_t' or p == 't': + return(self.inj_dat.getColumnByName('geocent_end_time')[self.event] + self.inj_dat.getColumnByName('geocent_end_time_ns')[self.event] * 1e-9 - self.segment_start) + elif p == 'psi': + return(self.inj_dat.getColumnByName('polarization')[self.event]) + elif p == 'ra' or p == 'RA': + return(self.inj_dat.getColumnByName('longitude')[self.event]) + elif p == 'sin_dec': + return(np.sin(self.inj_dat.getColumnByName('latitude')[self.event])) + elif p == 'dec': + return(self.inj_dat.getColumnByName('latitude')[self.event]) + elif p == 'cos_iota' or p == 'cosiota': + return(np.cos(self.inj_dat.getColumnByName('inclination')[self.event])) + elif p == 'iota': + return(self.inj_dat.getColumnByName('inclination')[self.event]) + elif p == 'phase' or p == 'phi' : + return(self.inj_dat.getColumnByName('coa_phase')[self.event]) + else: + return(self.inj_dat.getColumnByName(p)[self.event]) + +# ---------------------------------- +# Function for initializing runFlags +# ---------------------------------- +def readbwb(runFlags): + # -- Initialize variables with info about the job + info = "" + jobName = '' + restrictModel = '' + lalsimFlag = False + + # -- Open *bayeswave.run + bayeswaverunfile = 'bayeswave.run' + runFlags.bayeswaverunfile = bayeswaverunfile + bayeswave = open(bayeswaverunfile, 'r') + # -- Parse command line arguments + raw_lines = bayeswave.readlines() + lines = [l.lstrip() for l in raw_lines] + cmdline = '' + for l in lines: + if 'Command' in l: + cmdline = l.rstrip().split(' ') + + # Go through the different commands in bayeswave.run + for index, arg in enumerate(cmdline): + # Parameters describing where things live + if arg=='--runName': + jobName = cmdline[index+1] + print("The job name is: {0}".format(jobName)) + runFlags.jobName = jobName+'_' + + # Run settings + elif arg=='--trigtime': + gps = float(cmdline[index+1]) + runFlags.gps = gps + info = info + "The trigger time is GPS: {0}\n".format(gps) + elif arg=='--srate': + runFlags.srate = float(cmdline[index + 1]) + elif arg=='--segment-start': + runFlags.segment_start = float(cmdline[index + 1]) + elif arg=='--seglen': + runFlags.seglen = float(cmdline[index + 1]) + elif arg == '--Nchain': + runFlags.Nchain = int(cmdline[index + 1]) + + # IFOs + elif arg=='--ifo': + runFlags.ifoNames.append(cmdline[index+1]) + + # Windows, cleaning, and bayesLine + elif arg == '--bayesLine': + runFlags.bayesLine = True + elif arg == '--window': + runFlags.window = float(cmdline[index + 1]) + elif arg == '--CBCwindow': + runFlags.CBCwindow = float(cmdline[index + 1]) + elif arg == '--noClean': + print('\nThis run had no cleaning phase, --noClean\n') + runFlags.noClean = True + # Parameters describing what run type was + elif arg=='--glitchOnly': + print('\nThis run was executed with the --glitchOnly flag\n') + runFlags.modelList = ['glitch'] + runFlags.glitchOnly_flag = True + runFlags.noNoiseFlag = True + elif arg=='--signalOnly': + print('\nThis run was executed with the --signalOnly flag\n') + runFlags.modelList = ['signal'] + elif arg=='--cbcOnly': + print('\nThis run was executed with the --cbcOnly flag\n') + runFlags.modelList = ['cbc'] + runFlags.CBCType = True + elif arg=='--noiseOnly': + print('\nThis run was executed with the --noiseOnly flag\n') + runFlags.modelList = ['noise'] + elif arg=='--skipNoise': + print('\nThis run was executed with the --skipNoise flag\n') + runFlags.noNoiseFlag = True + # Combined types + elif arg=='--fullOnly': + print('\nThis run was executed with the --fullOnly flag\n') + runFlags.fullOnly_flag = True + runFlags.modelList = ('signal', 'glitch') + elif arg=='--GlitchCBC': + print('\nThis run was executed with the --GlitchCBC flag\n') + runFlags.GlitchCBC_flag = True + runFlags.CBCType = True + runFlags.modelList = ['glitch', 'cbc'] + runFlags.multi_type = True + elif arg=='--SignalCBC': + print('\nThis run was executed with the --SignalCBC flag\n') + runFlags.SignalCBC_flag = True + runFlags.CBCType = True + runFlags.multi_type = True + runFlags.modelList = ('signal', 'cbc') + elif arg=='--cleanOnly': + print('\nThis run was executed with the --cleanOnly flag\n') + runFlags.cleanOnly_flag = True + runFlags.modelList = ['clean'] + # Injection settings + elif arg=='--inj': + runFlags.injFlag = True + runFlags.mdc = True + runFlags.XMLfname = cmdline[index+1] + xmldoc = pipe_utils.ligolw_utils.load_filename(runFlags.XMLfname, contenthandler = pipe_utils.LIGOLWContentHandler, verbose = True) + runFlags.inj_dat = pipe_utils.lsctables.SimInspiralTable.get_table(xmldoc) + print('xml filename is {0}'.format(runFlags.XMLfname)) + elif arg =='--event': + # if injection then this indexes the event in the .xml file + runFlags.event = int(cmdline[index + 1]) + elif arg=='--MDC-cache': + mdc = True + elif arg=='--lalinspiralinjection': + lalsimFlag = True + # Likelihood calculation + elif arg=='--heterodyneL': + runFlags.heterodyneL = True + + # Printed output settings + elif arg == '--checkpoint': + runFlags.checkpoint = True + elif arg=='--verbose': + print('\nThis run was executed with the --verbose flag\n') + runFlags.verbose = True + runFlags.restrictModel = runFlags.modelList[0] + + + if lalsimFlag: + runFlags.injFlag = False + + # -- Determine number of IFOs + ifoNum = len(runFlags.ifoNames) + # -- Convert the IFO name list into a list of numbers useful for opening datafiles + runFlags.ifoList = [str(ii) for ii in range(0,len(runFlags.ifoNames))] + + # -- Getting flow + for index, arg in enumerate(cmdline): + for ifo in runFlags.ifoNames: + if arg == '--{0}-flow'.format(ifo): + #TODO THIS IS WRONG IF DIFFERENT IFOS HAVE DIFFERENT FLOW + runFlags.flow = float(cmdline[index+1]) + + info = info + "Detectors(s) found: {0}\n".format(len(runFlags.ifoList)) + info = info + "The detectors are: {0}\n".format(', '.join(runFlags.ifoNames)) + info = info + "Models are: {0}\n".format(', '.join(runFlags.modelList)) + info = info + "\tflow = {0}\n".format(runFlags.flow) + info = info + "\tseglen = {0}\n".format(runFlags.seglen) + info = info + "\tNchains = {0}\n".format(runFlags.Nchain) + info = info + "\tCleaning phase = {0}\n".format(not runFlags.noClean) + + if 'cbc' in runFlags.modelList: + info = info + "\theterodyneL = {0}\n".format(runFlags.heterodyneL) + + # -- If MDC, read SNR info + if runFlags.mdc: + for ifo in runFlags.ifoList: + snrFile = 'post/injection_whitened_moments_%s.dat'%(runFlags.ifoNames[int(ifo)]) + if not os.path.exists(snrFile): + sys.exit("\n {0} not found! \n".format(snrFile)) + snrdata = open(snrFile, 'r') + snrdata.readline() + snr = float(snrdata.readline().split(' ')[0]) + runFlags.snrList.append(snr) + snrdata.close() + info = info + 'Injected SNR in detector {0} = {1}\n'.format(runFlags.ifoNames[int(ifo)],runFlags.snrList[-1]) + bayeswave.close() + # -- Report to user + print("{0}".format(info)) + +""" + Functions for calculating redshift (used in Flags.get_xml_paramter) +""" +def redshift_from_distance(DL): # distance (MPC) + #double zlow, zhigh; + #double zmid; + #double Dlow, Dhigh, Dmid; + #double u, v, w; + #int i; + + zlow = 1.0; + zhigh = 1000.0; + + if(DL < 1.0): + zlow = 1.0e-10 + zhigh = 3.0e-4 + elif (DL < 100.0): + zlow = 2.0e-4 + zhigh = 5.0e-2 + elif (DL < 1000.0): + zlow = 1.0e-2 + zhigh = 3.0e-1 + elif (DL < 10000.0): + zlow = 1.0e-1 + zhigh = 2.0 + + zmid = (zhigh + zlow) / 2.0; + + Dlow = DL_fit(zlow); + Dhigh = DL_fit(zhigh); + Dmid = DL_fit(zmid); + zmid = (zhigh + zlow) / 2.0; + + Dlow = DL_fit(zlow); + Dhigh = DL_fit(zhigh); + Dmid = DL_fit(zmid); + + for i in range(30): + + u = Dlow-DL; + v = Dmid-DL; + + if(u*v < 0.0): + zhigh = zmid; + Dhigh = Dmid; + else: + zlow = zmid; + Dlow = Dmid; + + zmid = (zhigh+zlow)/2.0; + Dmid = DL_fit(zmid); + + return(zmid); + +# fits distance to redshift assuming flat cosmological constant +def DL_fit(z): + # Planck values https://arxiv.org/abs/1807.06209 + Om = 0.315; + H0 = 67.4; + #double RH, x1, x2, Om1, Om2; + + CLIGHT = 2.99792458e8 # m / s + MPC = 3.085677581e+22 # meters + + # See http://arxiv.org/pdf/1111.6396v1.pdf + + RH = CLIGHT / H0 * (1.0e-3 * MPC); + + x1 = (1.0-Om)/(Om); + + x2 = (1.0-Om)/(Om*pow((1.0+z),3.0)); + + Om1 = (1.0+1.32*x1+0.4415*x1*x1+0.02656*x1*x1*x1)/(1.0+1.392*x1+0.5121*x1*x1+0.03944*x1*x1*x1); + + Om2 = (1.0+1.32*x2+0.4415*x2*x2+0.02656*x2*x2*x2)/(1.0+1.392*x2+0.5121*x2*x2+0.03944*x2*x2*x2); + + D = 2.0*RH*(1.0+z) / np.sqrt(Om)*(Om1-Om2 / np.sqrt(1.0+z)); + + return(D / MPC); diff --git a/BayesWaveUtils/scripts/megaplot.py b/BayesWaveUtils/scripts/megaplot.py index 9ec8519..a00a568 100755 --- a/BayesWaveUtils/scripts/megaplot.py +++ b/BayesWaveUtils/scripts/megaplot.py @@ -32,7 +32,7 @@ import re import traceback from bayeswave_pipe import bayeswave_pipe_utils as pipe_utils import lalsimulation as lalsim - +import BW_Flags @@ -143,48 +143,6 @@ cbc_dict['lambda2'] = [20, 'lambda2', r'$\lambda_2$', 'Lambda 2'] cbc_dict['lambda3'] = [21, 'lambda3', r'$\lambda_3$', 'Lambda 3'] cbc_dict['lambda4'] = [22, 'lambda4', r'$\lambda_4$', 'Lambda 4'] -class Flags: - def __init__(self): - self.fullOnly_flag = 0 - self.GlitchCBC_flag = 0 - self.SignalCBC_flag = 0 - self.glitchOnly_flag = 0 - self.CBCType = 0 - self.multi_type = 0 # true for full, cbcsignal, cbcglitch - self.modelList = () # Gets filled in readbwb, theres probably a clearer way to to this - self.bayesLine = False - self.XMLfname = None - self.injFlag = False - self.event = 0 - self.gps = 0 # gps time of event - self.CBCwindow = -1 #window around cbc event - self.segment_start = 0 - self.CBCName = 0 - self.multi_type = 0 # true for full, cbcsignal, cbcglitch - self.modelList = () # Gets filled in readbwb, theres probably a clearer way to to this - - self.verbose = 0 - self.noClean = 0 - self.Nchain = 20 # default value for number of chains - - self.checkpoint = False - self.noClean = False - - def run_dirname(self): - if self.fullOnly_flag: - return('full') - if self.SignalCBC_flag: - return('cbcsignal') - if self.GlitchCBC_flag: - return('cbcglitch') - if len(self.modelList) == 1: - return(self.modelList[0]) - #if self.glitchOnly_flag: - #return('glitch') - else: - print("uh oh! what's the dirname?") - return('') - modelList = ('signal', 'glitch', 'cbc') postDir = 'post/' @@ -192,7 +150,7 @@ plotsDir = 'plots/' htmlDir = 'html/' #### --fullOnly ### -RUN_FLAGS = Flags() +RUN_FLAGS = BW_Flags.Flags(workdir) # I think RUN_FLAGS gets set twice, I dont think this codeblock is necesssary if os.path.exists(postDir+'full'): @@ -227,6 +185,8 @@ def set_color(model): return(gcolor) #glitchcbccolor) elif model == 'cbcsignal': return(signalcbccolor) + elif model == 'clean': + return('mediumseagreen') return(ncolor) @@ -240,144 +200,6 @@ injcolor = 'teal' # Read in run data and info # ###################################################################################################################### - -# ---------------------------------- -# Extract information about the run -# ---------------------------------- -# initialize RUN_FLAGS also -def readbwb(runFlags): - # -- Initialize variables with info about the job - info = "" - ifoNames = [] - ifoList = [] - jobName = '' - restrictModel = '' - injFlag = False - lalsimFlag = False - runFlags.noNoiseFlag = False - mdc = False - snrList = [] - # -- Open *bayeswave.run - bayeswaverunfile = glob.glob('*bayeswave.run') - if not len(bayeswaverunfile): - sys.exit("\n *bayeswave.run not found! \n") - bayeswaverunfile = bayeswaverunfile[0] - bayeswave = open(bayeswaverunfile, 'r') - # -- Parse command line arguments - raw_lines = bayeswave.readlines() - lines = [l.lstrip() for l in raw_lines] - cmdline = '' - for l in lines: - if 'Command' in l: - cmdline = l.rstrip().split(' ') - - for index, arg in enumerate(cmdline): - # Determine the IFOs involved - if arg=='--ifo': - ifoNames.append(cmdline[index+1]) - # Determine the trigger time - elif arg=='--trigtime': - gps = float(cmdline[index+1]) - runFlags.gps = gps - info = info + "The trigger time is GPS: {0}\n".format(gps) - # Determine the job name - elif arg=='--runName': - jobName = cmdline[index+1] - print("The job name is: {0}".format(jobName)) - jobName = jobName+'_' - elif arg=='--fullOnly': - print('\nThis run was executed with the --fullOnly flag\n') - runFlags.fullOnly_flag = True - runFlags.modelList = ('signal', 'glitch') - # Determine if the job was glitchOnly, noiseOnly, or signalOnly - elif arg=='--glitchOnly': - print('\nThis run was executed with the --glitchOnly flag\n') - restrictModel = 'glitch' - runFlags.glitchOnly_flag = True - runFlags.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' - elif arg=='--cbcOnly': - print('\nThis run was executed with the --cbcOnly flag\n') - restrictModel = 'cbc' - runFlags.modelList = ('cbc') - runFlags.CBCType = True - elif arg == '--noClean': - print('\nThis run had no cleaning phase, --noClean\n') - runFlags.noClean = True - elif arg=='--GlitchCBC': - print('\nThis run was executed with the --GlitchCBC flag\n') - runFlags.GlitchCBC_flag = True - runFlags.CBCType = True - runFlags.modelList = ('glitch', 'cbc') - runFlags.multi_type = True - elif arg == '--CBCwindow': - runFlags.CBCwindow = float(cmdline[index + 1]) - elif arg=='--SignalCBC': - print('\nThis run was executed with the --SignalCBC flag\n') - runFlags.SignalCBC_flag = True - runFlags.CBCType = True - runFlags.multi_type = True - runFlags.modelList = ('signal', 'cbc') - elif arg=='--skipNoise': - print('\nThis run was executed with the --skipNoise flag\n') - runFlags.noNoiseFlag = True - elif arg=='--verbose': - print('\nThis run was executed with the --verbose flag\n') - runFlags.verbose = True - elif arg == '--Nchain': - runFlags.Nchain = int(cmdline[index + 1]) - elif arg == '--checkpoint': - runFlags.checkpoint = True - elif arg=='--segment-start': - runFlags.segment_start = float(cmdline[index + 1]) - elif arg=='--inj': - injFlag = True - runFlags.injFlag = True - mdc = True - runFlags.XMLfname = cmdline[index+1] - print('xml filename is {0}'.format(runFlags.XMLfname)) - elif arg =='--event': - # if injection then this indexes the event in the .xml file - runFlags.event = int(cmdline[index + 1]) - elif arg=='--MDC-cache': - mdc = True - elif arg=='--lalinspiralinjection': - lalsimFlag = True - elif arg=='--bayesLine': - runFlags.bayesLine = True - - if lalsimFlag: - injFlag = False - self.injFlag = False - # -- Determine number of IFOs - ifoNum = len(ifoNames) - # -- Convert the IFO name list into a list of numbers useful for opening datafiles - ifoList = [str(ii) for ii in range(0,len(ifoNames))] - info = info + "Detectors(s) found: {0}\n".format(len(ifoList)) - info = info + "The detectors are: {0}\n".format(', '.join(ifoNames)) - # -- If MDC, read SNR info - if mdc: - for ifo in ifoList: - snrFile = jobName+'post/injection_whitened_moments_%s.dat'%(ifoNames[int(ifo)]) - if not os.path.exists(snrFile): - sys.exit("\n {0} not found! \n".format(snrFile)) - snrdata = open(snrFile, 'r') - snrdata.readline() - snr = float(snrdata.readline().split(' ')[0]) - snrList.append(snr) - snrdata.close() - info = info + 'Injected SNR in detector {0} = {1}\n'.format(ifoNames[int(ifo)],snrList[-1]) - bayeswave.close() - # -- Report to user - print("{0}".format(info)) - return(jobName, restrictModel, mdc, injFlag, bayeswaverunfile, ifoList, ifoNames, gps, snrList, info, runFlags) - - # --------------------------------------------- # Get Median waveform with CIs # --------------------------------------------- @@ -1663,11 +1485,81 @@ def plot_model_dims(modelList, ifoList, ifoNames, plotsDir, runFlags = RUN_FLAGS return sig_noise, sig_gl +#-------------------- +# plot cleaning phase +#-------------------- +def plot_cleaning_phase_time_domain(runFlags, plotsDir='plots/'): + + if not runFlags.noClean: + print('making cleaning phase plot(s)') + + #find window ends + window_start = runFlags.gps - runFlags.segment_start + runFlags.window / 2 + window_end = runFlags.gps - runFlags.segment_start - runFlags.window / 2 + + color = set_color(runFlags.run_dirname()) + cleancolor = set_color('clean') + + for ifo in ifoNames: + name = 'cleaning_phase_time_domain_plot_{0}'.format(ifo) + f, ax = plt.subplots(dpi=200) + + #get glitches found in cleaning phase & those found in run + clean = get_waveform('post/clean/glitch_median_time_domain_waveform_{0}.dat'.format(ifo)) + run = get_waveform('post/{1}/{1}_median_time_domain_waveform_{0}.dat'.format(ifo, runFlags.run_dirname())) + #plot cleaning phase waveform + ax.plot(clean[0], clean[1],color = cleancolor, label = 'clean', alpha=0.7,zorder=10,linewidth=1) + ax.fill_between(clean[0], clean[4], clean[5],color = cleancolor, alpha=0.3,edgecolor=None,zorder=8) + #plot run waveform + ax.plot(run[0], run[1],color = color, label = 'run',linewidth=1) + ax.fill_between(run[0], run[4], run[5],color = color, alpha=0.3,edgecolor=None) + + #plot window + ax.axvspan(window_start, window_end, color='gray',alpha=0.2,zorder=1,label='window') + + #other plot elements + ax.legend(fontsize=6); ax.set_xlim(0, runFlags.seglen) + ax.set_xlabel('Time (s)'); ax.set_ylabel('Whitened Strain') + ax.set_title(runFlags.run_dirname()+', '+ifo) + plt.savefig(plotsDir+name); plt.close() +def plot_cleaning_phase_frequency_domain(runFlags, plotsDir='plots/'): + + if not runFlags.noClean: + color = set_color(runFlags.run_dirname()) + for ifo in ifoNames: + name = 'cleaning_phase_frequency_domain_plot_{0}'.format(ifo) + f, ax = plt.subplots(dpi=200) + #get glitches found in cleaning phase & those found in run + clean = get_waveform('post/clean/glitch_median_frequency_domain_waveform_spectrum_{0}.dat'.format(ifo)) + run = get_waveform('post/{1}/{1}_median_frequency_domain_waveform_spectrum_{0}.dat'.format(ifo, runFlags.run_dirname())) + + #get psd in cleaning phase and full run + clean_psd = get_waveform('post/clean/glitch_median_PSD_{0}.dat'.format(ifo)) + run_psd = get_waveform('post/{1}/{1}_median_PSD_{0}.dat'.format(ifo, runFlags.run_dirname())) + + #plot cleaning phase waveform + ax.plot(clean[0], clean[1],color = set_color('clean'), label = 'clean', alpha=0.7,zorder=10,linewidth=1) + ax.fill_between(clean[0], clean[4], clean[5],color = set_color('clean'), alpha=0.3,edgecolor=None,zorder=8) + + #plot run waveform + ax.plot(run[0], run[1],color = color, label = 'run',linewidth=1) + ax.fill_between(run[0], run[4], run[5],color = color, alpha=0.3,edgecolor=None) + + #plot PSDs + ax.plot(clean_psd[0],clean_psd[1],color='darkred',linewidth=1,label='clean psd') + ax.plot(run_psd[0],run_psd[1],color='k',linewidth=1,label='run psd') + + #other plot elements + ax.legend(fontsize=6,ncol=2) + ax.set_xlim(runFlags.flow, runFlags.srate/2); ax.set_ylim(1e-50,1e-42) + ax.set_xscale('log'); ax.set_yscale('log') + ax.set_xlabel('Frequency (Hz)'); ax.set_ylabel('Power') + ax.set_title(ifo) + plt.savefig(plotsDir+name); plt.close() ###################################################################################################################### -# # BWB webpage production # ###################################################################################################################### @@ -1961,9 +1853,6 @@ def splash_page(plotsDir, ifoNames, snrList, bayeswaverunfile, sig_gl, sig_noise html_string += '\n' - - - html_string +='

Evidence for Signal

\n' if not runFlags.fullOnly_flag: if restrictModel == 'signal': @@ -2075,6 +1964,9 @@ def make_index(htmlDir, plotsDir, model, gps, ifoList, ifoNames, snrList, bayesw index.write('