Make SNR plots for S6 injections

git-svn-id: https://svn.ligo.caltech.edu/svn/bayeswave/trunk@62 c56465c9-8126-4a4f-9d7d-ac845eff4865

postprocess/s6/plot_evidence.py

+
+
+import matplotlib
+matplotlib.use('Agg')
+import numpy as np
+import matplotlib.pyplot as plt
+import glob
+import os
+import sys
+# from scipy.optimize import curve_fit
+
+# --------------
+# Plot Evidence
+# --------------
+plt.close('all')
+maxodds=50
+
+# -- Get directory list
+glitch = []
+signal = []
+topdir = os.getcwd()
+dirList = glob.glob('*/job*')
+# dirList = [os.path.abspath(directory) for directory in dirList]
+# print dirList
+
+# -- Loop over jobs
+snrList = []
+goodJobList = []
+for job in dirList:
+    os.chdir(job)
+    print "Entered", os.getcwd()
+    jobname = glob.glob('*bayeswave.run')[0]
+    bayeswave = open(jobname, 'r')
+
+    # ---------------------------------------
+    # Extract information from bayeswave.run
+    # ---------------------------------------
+    ifoNames = []
+    for line in bayeswave:
+        spl = line.split()
+        if len(spl) < 1: continue
+        del spl[-1]
+        # Determine names of IFOs
+        if spl[0] == 'Command' and spl[1] == 'line:':
+            for index, splElement in enumerate(spl):
+                if splElement == '--ifo':
+                    ifoNames.append(spl[index+1]) 
+            continue         
+        # Determine number of IFOs
+        if ' '.join(spl) == 'number of detectors':
+            spl = line.split()
+            ifoNum = spl[3]
+            continue         
+        # Determine gps trigger time
+        if ' '.join(spl) == 'trigger time (s)':
+            spl = line.split()
+            gps = spl[3]
+            break
+    ifoNum = int(ifoNum)
+    if ifoNum == 1:
+        ifoList = ['0']
+    elif ifoNum == 2:
+        ifoList = ['0', '1']
+    elif ifoNum == 3:
+        ifoList = ['0', '1', '2']
+    
+    # ----------------------------------
+    # Determine name of job 
+    # ----------------------------------
+    jobname = jobname.split('_')
+    del jobname[-1]
+    if not len(jobname) == 0:
+        jobname = '_'.join(jobname)+'_'
+    else:
+        jobname = ''
+
+    try:
+        infile = open(str(jobname)+'evidence.dat', 'r')
+    except:
+        os.chdir(topdir)
+        continue
+    
+    for line in infile:
+        spl = line.split()
+        if spl[0] == 'noise':
+            sig_noise = -float(spl[1])
+            signal.append(sig_noise)
+        if spl[0] == 'glitch':
+            sig_gl = -float(spl[1])
+            glitch.append(sig_gl)
+    infile.close()
+
+    # -- Parse log file
+    infile = open('condorOut.txt', 'r')
+    redo = False
+    snr = []
+    print "Parsing log file..."
+    for line in infile:
+        if line.find('Recompute') > 0:
+            print "Found recompute line"
+            print line
+            redo = True
+        if redo:
+            if line.find('Injected SNR') > -1:
+                print "Found line with SNR"
+                print line
+                snr.append(float(line.split()[2]))
+    if len(snr) == 0:
+        sys.exit("No injections found in this run!")
+    # snrList.append(np.min(snr))
+    # snrList.append(np.log(snr[0]/snr[1]))
+    snrList.append(np.sqrt(snr[0]**2 + snr[1]**2 ))
+    goodJobList.append(job)
+
+    infile.close()
+    os.chdir(topdir)
+
+# ----------------------
+# Make histogram of SNRs
+# ----------------------
+plt.figure()
+plt.hist(snrList)
+plt.xlabel('SNR')
+plt.ylabel('Number')
+plt.savefig('snr_hist.png')
+plt.close()
+
+# ---------------------------
+# Make plot of all injections 
+# ---------------------------
+plt.figure()
+# for sig_noise, sig_gl, snr in zip(signal, glitch, snrList):
+#      plt.scatter(limit_odds(sig_gl), limit_odds(sig_noise), c=snr, s=40, cmap=matplotlib.cm.gray)
+
+plt.scatter(glitch, signal, c=snrList, s=40)
+cbar = plt.colorbar()
+cbar.set_label('Network SNR', rotation=270)
+plt.ylabel('log( E_signal / E_noise )')
+plt.xlabel('log( E_signal / E_glitch )')
+plt.axis([np.min(glitch), np.max(glitch), np.min(signal), np.max(signal)])
+    
+# plt.fill_between([0,maxodds], [0, 0], [maxodds, maxodds], facecolor='green', interpolate=True, alpha=0.3)
+# plt.fill_between([0,maxodds], [0, 0], [-maxodds, -maxodds], facecolor='red', interpolate=True, alpha=0.3)
+# plt.fill_between([-maxodds,0], [0, 0], [maxodds, maxodds], facecolor='red', interpolate=True, alpha=0.3)
+
+plt.grid()
+    
+# plt.title('network SNR')
+plt.savefig('net_snr_all.png')
+
+print np.max(signal)
+print np.max(glitch)
+print signal
+print "Found {0} injections".format(len(glitch))
+
+plt.axis([-50, 50, -50, 50])
+plt.savefig('net_snr_zoom.png')
+plt.close()
+
+
+# --------------------
+# Make SNR table 
+# -------------------
+snrTable = open('snrTable.html', 'w')
+snrTable.write('<html><body>')
+snrTable.write('<table>')
+snrTable.write('<tr><b><td>SNR <td>Job    <td>Evidence(S/N)  <td>Evidence(S/G)</b></tr>')
+
+# -- Make tuple, sort by SNR
+raw_eventList = zip(snrList, signal, glitch, goodJobList)
+eventList =  sorted(raw_eventList, key=lambda event: event[0], reverse=True)
+
+for (snr, sn, sg, job) in eventList:
+    snrTable.write("""<tr>
+<td>{0:.2f}  
+<td><a href='{1}'>{1}</a>  
+<td>{2:.2f}  
+<td>{3:.2f}
+</tr>
+""".format(snr, job, sn, sg) )
+snrTable.write('</table></html></body>')
+snrTable.close()
+
+
+# -----------------------------
+# -- Try making efficiency plot
+# ----------------------------
+glitch = np.array(glitch)
+signal = np.array(signal)
+snrList = np.array(snrList)
+
+binsize = 2
+snrBins = np.arange(0,20,binsize)
+
+effList = np.array([])
+for bin in snrBins:
+    indx = np.logical_and(bin<snrList, snrList<bin+binsize)
+    this_gl = glitch[indx]
+    this_n  = signal[indx]
+    detected = np.logical_and(this_gl>0, this_n>0)
+    eff = detected.sum()*1.0 / this_gl.size
+    effList = np.append(effList, eff)
+
+plt.figure()
+plt.plot(snrBins + binsize/2.0, effList, 'o')
+plt.xlabel('Network SNR')
+plt.ylabel('Detection Efficiency')
+plt.savefig('eff.png')
+plt.close()
+
+# -- Make plot without high SNR signals
+snrmax = 50
+glitch = np.array(glitch)
+signal = np.array(signal)
+snrList = np.array(snrList)
+qglitch = glitch[ snrList < snrmax ]
+qsignal = signal[ snrList < snrmax ]
+qsnr    = snrList[ snrList < snrmax ] 
+
+if len(qglitch) > 0 and len(qsignal) > 0 and len(qsnr) > 0:
+    plt.figure()
+    plt.scatter(qglitch, qsignal, c=qsnr, s=40)
+    cbar = plt.colorbar()
+    cbar.set_label('Network SNR', rotation=270)
+    plt.ylabel('log( E_signal / E_noise )')
+    plt.xlabel('log( E_signal / E_glitch )')
+    plt.grid()
+    plt.axis([np.min(qglitch)-10, np.max(qglitch)+10, np.min(qsignal)-10, np.max(qsignal)+10])
+    plt.savefig('quiet.png')
+    plt.close()
+
+# -- Super low SNR
+snrmax = 18
+glitch = np.array(glitch)
+signal = np.array(signal)
+snrList = np.array(snrList)
+qglitch = glitch[ snrList < snrmax ]
+qsignal = signal[ snrList < snrmax ]
+qsnr    = snrList[ snrList < snrmax ] 
+
+if len(qglitch) > 0 and len(qsignal) > 0 and len(qsnr) > 0:
+    plt.figure()
+    plt.scatter(qglitch, qsignal, c=qsnr, s=40)
+    cbar = plt.colorbar()
+    cbar.set_label('Network SNR', rotation=270)
+    plt.ylabel('log( E_signal / E_noise )')
+    plt.xlabel('log( E_signal / E_glitch )')
+    plt.grid()
+    plt.axis([np.min(qglitch)-10, np.max(qglitch)+10, np.min(qsignal)-10, np.max(qsignal)+10])
+    plt.savefig('super_quiet.png')
+    plt.close()
+
+
+# -----------------------
+# Make top level web page
+# ----------------------
+outfile = open('top.html', 'w')
+outfile.write('<html><body>')
+outfile.write('<h1>BayesWave Burst Injection Run</h1>')
+outfile.write("<h4><a href='./snrTable.html'>Link to table of all injections</a></h4><br/>")
+outfile.write('<h3>All Injections</h3>')
+outfile.write("<img src='net_snr_all.png' alt='net_snr_all.png' width=600>")
+outfile.write("<h3>All injections, fixed axis scale</h3>")
+outfile.write("<img src='net_snr_zoom.png' alt='zoom' width=600>")
+outfile.write("<h3>Injections with moderate SNR</h3>")
+outfile.write("<img src='quiet.png' alt='quiet' width=600>")
+outfile.write("<h3>Injections with low SNR</h3>")
+outfile.write("<img src='super_quiet.png' alt='super_quiet' width=600>")
+outfile.write("</body></html>")
+outfile.close()

postprocess/skymap/skygram.py

-
 import numpy as np
 import healpy as hp
 import matplotlib
@@ -40,7 +39,12 @@ auto = autocorr[autocorr.size/2:]
 plt.figure()
 plt.plot(auto)
 
+rand_ra = 2*np.pi*np.random.rand(ralist.size)
+randcorr = np.correlate(rand_ra, rand_ra, mode='full')
+plt.plot(randcorr[autocorr.size/2:])
+
 plt.ylabel('Autocorrelation')
 plt.xlabel('Number of sample offsets')
+plt.xlim([0, 30])
 plt.savefig('RA_autocorr.png')