Add full only Bayes factors

BayesWaveUtils/scripts/megaplot.py

@@ -101,7 +101,7 @@ ncolor = 'darkgrey'
 gcolor = 'darkgoldenrod'
 scolor = 'darkorchid'
 
-ifoColors = ['darkgoldenrod','darkkhaki','darkseagreen','olive','cadetblue']
+ifoColors = ['darkgoldenrod','darkkhaki','darkseagreen','olive','cadetblue','green','slategray','darkcyan']
 
 injcolor = 'teal'
 
@@ -132,11 +132,12 @@ def readbwb():
     bayeswaverunfile = bayeswaverunfile[0]
     bayeswave = open(bayeswaverunfile, 'r')
     # -- Parse command line arguments
-    lines = bayeswave.readlines()
-    if(len(lines[8].strip())!=0):
-        cmdline = lines[8].rstrip().split(' ')
-    else:
-        cmdline = lines[7].rstrip().split(' ')
+    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
@@ -173,12 +174,7 @@ def readbwb():
     # -- Determine number of IFOs
     ifoNum = len(ifoNames)
     # -- Convert the IFO name list into a list of numbers useful for opening datafiles
-    if ifoNum == 1:
-        ifoList = ['0']
-    elif ifoNum == 2:
-        ifoList = ['0', '1']
-    elif ifoNum == 3: 
-        ifoList = ['0', '1', '2']
+    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
@@ -992,6 +988,40 @@ def plot_model_dims(modelList, ifoList, ifoNames, plotsDir, fullOnly_flag=0):
         ax1.grid()
         plt.savefig(plotsDir+'model_dimensions_histo.png')
 
+    # -- 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
+        
+        ev_sig = 0
+        ev_glitch = 0
+        ev_full = 0
+        for s, g in zip(signal_on, glitch_on):
+            if s and not g:
+                ev_sig += 1
+            if g and not s:
+                ev_glitch += 1
+            if s and g:
+                ev_full += 1
+    else:
+        ev_sig = 0
+        ev_glitch = 0
+        ev_full = 0
+
+    # Actually full vs. signal only
+    try:
+        sig_noise = ev_full/ev_sig
+    except:
+        sig_noise = np.inf
+
+    # Actually full vs. glitch only
+    try:
+        sig_gl = ev_full/ev_gl
+    except:
+        sig_gl = np.inf
+
+    return sig_noise, sig_gl
+
 
 
 ######################################################################################################################
@@ -1127,7 +1157,10 @@ def make_skymap_page(htmlDir, plotsDir):
 
 def splash_page(plotsDir, ifoNames, snrList, bayeswaverunfile, sig_gl, sig_noise, fullOnly_flag=0):
     html_string = ''
-    html_string += '    <a href="./'+plotsDir+'odds.png"><img src="./'+plotsDir+'odds.png" style="float: right;" width=600 alt="odds.png"></a>\n'
+    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))
     if mdc:
         html_string += '    <h3>Matched Filter SNRs of Injections</h3>\n'
@@ -1143,8 +1176,11 @@ def splash_page(plotsDir, ifoNames, snrList, bayeswaverunfile, sig_gl, sig_noise
         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)
     else:
-        html_string +='    <h3>This is full model run.</h3>\n'
-        html_string +='    This was run that simultaneously fits for signals and glitches, therefore there are no Bayes factors for this run.\n'
+        html_string +='    <h3>Bayes factor for signal+glitch vs signal only:</h3>\n'
+        html_string +='    {0} <br/>\n'.format(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 +='  </p>\n'
 
     return html_string
@@ -1680,7 +1716,10 @@ except:
 
 
 if 'signal' in modelList or 'glitch' in modelList:
-    plot_model_dims(modelList, ifoList, ifoNames, plotsDir, fullOnly_flag)
+    ev = plot_model_dims(modelList, ifoList, ifoNames, plotsDir, fullOnly_flag)
+    if fullOnly_flag:
+        sig_noise = ev[0]
+        sig_gl = ev[1]
 
 
 # -- Posterior predictive checks: