plot pec now makes an output file

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

postprocess/s6/condor_pp.py

@@ -3,8 +3,8 @@ import glob
 import os
 
 # -- Parameters
-mdclog = '/home/jkanner/baysewave/svn/trunk/burstinj/s6/BurstMDC-BRST_S6-Log.txt'
-
+# mdclog = '/home/jkanner/baysewave/svn/trunk/burstinj/s6/BurstMDC-BRST_S6-Log.txt'
+mdclog = '/home/jkanner/baysewave/svn/trunk/burstinj/s6/BurstMDC-ELPTC_S6-Log.txt'
 runfile = open('run_pp.sh', 'w')
 dirlist = glob.glob('*/job*')
 dirlist = [os.path.abspath(directory) for directory in dirlist]

postprocess/s6/plot_pec.py

@@ -65,6 +65,34 @@ def get_median(strain_master, ifo):
         up_vec[sample] = up_quart            
     return (median_waveform, up_vec, down_vec)
 
+
+# ------------------------------------------------------------
+# Calculate the mode and confidence interval from a histogram
+# ------------------------------------------------------------
+def get_mode(n, bins):
+    
+    # -- Take average of bins to get rid of size mismatch
+    if (bins.size) == (n.size + 1):
+        new_bins = [0.5*(bins[i] + bins[i+1]) for i in range(0,n.size)]
+        bins = new_bins
+
+    # -- Set the index of the mode to the tallest peak
+    modeindx = n.argmax()
+    samps_90 = 0.90*n.sum()
+    
+    for i in np.arange(n.size):
+        lowindx = np.max([0, modeindx-i])
+        highindx = np.min([n.size-1, modeindx+i])
+        n_inrange = n[lowindx:highindx]
+        if n_inrange.sum() >= samps_90:
+            low = bins[lowindx]
+            high = bins[highindx]
+            mode = bins[modeindx]
+            break
+    return(mode, low, high)
+    
+
+    
 # -- Create the Bayesogram
 def bayesogram(filename, time, twodrange):
     strain_master = np.loadtxt(filename)
@@ -105,7 +133,6 @@ parser.add_argument('--mdc', default='True', help='Is this this an MDC? True/Fal
 args = parser.parse_args()
 mdc = args.mdc
 
-print mdc
 # ---------------------------------------
 # Extract information from bayeswave.run
 # ---------------------------------------
@@ -316,8 +343,10 @@ for ifo in ifoList:
         # -- Read moments file
         if white:
             moments = np.recfromtxt('w_moments.dat.{0}'.format(ifo), names=True)
+            outfile = open('w_mode_{0}.txt'.format(ifo), 'w')
         else:
             moments = np.recfromtxt('moments.dat.{0}'.format(ifo), names=True)
+            outfile = open('mode_{0}.txt'.format(ifo), 'w')
 
         # -- Histogram overlaps
         plt.figure()
@@ -347,9 +376,13 @@ for ifo in ifoList:
             plt.savefig('energy_{0}.png'.format(ifo))
         plt.close()
 
+        mode, low, high = get_mode(rn, rbins)
+        outfile.write("{0} {1} {2} \n".format(mode, low, high))
+        print "Got log10(energy) of {0}, with range {1} to {2}".format(mode, low, high)
+
         # -- Histogram central time
         plt.figure()
-        rn,bins,patches = plt.hist(moments['time_rec'] - 2, bins=100, label='Recovered', alpha=0.5, log=True, linewidth=0)
+        rn,rbins,patches = plt.hist(moments['time_rec'] - 2, bins=100, label='Recovered', alpha=0.5, log=True, linewidth=0)
         n,bins,patches = plt.hist(moments['time_inj'] - 2, bins=100, label='Injected', alpha=0.5, color='yellow', log=True)
         plt.xlabel("Central Time")
         plt.title('IFO {0}'.format(ifo))
@@ -361,91 +394,102 @@ for ifo in ifoList:
         else:
            plt.savefig('tzero_{0}.png'.format(ifo))
         plt.close()
+        mode, low, high = get_mode(rn, rbins)
+        outfile.write("{0} {1} {2} \n".format(mode, low, high))
 
-    # -- Histogram of duration
-    plt.figure()
-    rn,bins,patches = plt.hist(np.log10(moments['duration_rec']), bins=100, label='Recovered', alpha=0.5)
-    n,bins,patches = plt.hist(np.log10(moments['duration_inj']), bins=100, label='Injected', alpha=0.5, color='yellow')
-    plt.xlabel("log10(Duration)")
-    plt.title('IFO {0}'.format(ifo))
-    plt.ylim(0, rn.max())
-    plt.legend()
-    if white:
-        plt.savefig('wduration_{0}.png'.format(ifo))
-    else:
-        plt.savefig('wduration_{0}.png'.format(ifo))
-    plt.close()
+        # -- Histogram of duration
+        plt.figure()
+        rn,rbins,patches = plt.hist(np.log10(moments['duration_rec']), bins=100, label='Recovered', alpha=0.5)
+        n,bins,patches = plt.hist(np.log10(moments['duration_inj']), bins=100, label='Injected', alpha=0.5, color='yellow')
+        plt.xlabel("log10(Duration)")
+        plt.title('IFO {0}'.format(ifo))
+        plt.ylim(0, rn.max())
+        plt.legend()
+        if white:
+            plt.savefig('wduration_{0}.png'.format(ifo))
+        else:
+            plt.savefig('duration_{0}.png'.format(ifo))
+        plt.close()
+        mode, low, high = get_mode(rn, rbins)
+        outfile.write("{0} {1} {2} \n".format(mode, low, high))
 
-    # -- Histogram of central frequency
-    plt.figure()
-    rn,bins,patches = plt.hist(moments['f_rec'], bins=100, label='Recovered', alpha=0.5, linewidth=0)
-    n,bins,patches = plt.hist(moments['f_inj'], bins=100, label='Injected', alpha=0.5, color='yellow')
-    plt.xlabel("Central Freq (Hz)")
-    plt.title('IFO {0}'.format(ifo))
-    plt.ylim(0, rn.max())
-    plt.legend()
-    if white:
-        plt.savefig('wfreq_{0}.png'.format(ifo))
-    else:
-        plt.savefig('freq_{0}.png'.format(ifo))
-    plt.close()
+        # -- Histogram of central frequency
+        plt.figure()
+        rn,rbins,patches = plt.hist(moments['f_rec'], bins=100, label='Recovered', alpha=0.5, linewidth=0)
+        n,bins,patches = plt.hist(moments['f_inj'], bins=100, label='Injected', alpha=0.5, color='yellow')
+        plt.xlabel("Central Freq (Hz)")
+        plt.title('IFO {0}'.format(ifo))
+        plt.ylim(0, rn.max())
+        plt.legend()
+        if white:
+            plt.savefig('wfreq_{0}.png'.format(ifo))
+        else:
+            plt.savefig('freq_{0}.png'.format(ifo))
+        plt.close()
+    
+        mode, low, high = get_mode(rn, rbins)
+        outfile.write("{0} {1} {2} \n".format(mode, low, high))
+        print "Got Central Freq. {0}, with range {1} to {2}".format(mode, low, high)
 
-    # -- Histogram of bandwidth
-    plt.figure()
-    rn,bins,patches = plt.hist(moments['ban_rec'], bins=100, label='Recovered', alpha=0.5, linewidth=0)
-    n,bins,patches = plt.hist(moments['band_inj'], bins=100, label='Injected', alpha=0.5, color='yellow')
-    plt.xlabel("Bandwidth (Hz)")
-    plt.title('IFO {0}'.format(ifo))
-    plt.ylim(0, rn.max())
-    plt.legend()
-    if white:
-        plt.savefig('wband_{0}.png'.format(ifo))
-    else:
-        plt.savefig('band_{0}.png'.format(ifo))
-    plt.close()
+        # -- Histogram of bandwidth
+        plt.figure()
+        rn,rbins,patches = plt.hist(moments['ban_rec'], bins=100, label='Recovered', alpha=0.5, linewidth=0)
+        n,bins,patches = plt.hist(moments['band_inj'], bins=100, label='Injected', alpha=0.5, color='yellow')
+        plt.xlabel("Bandwidth (Hz)")
+        plt.title('IFO {0}'.format(ifo))
+        plt.ylim(0, rn.max())
+        plt.legend()
+        if white:
+            plt.savefig('wband_{0}.png'.format(ifo))
+        else:
+            plt.savefig('band_{0}.png'.format(ifo))
+        plt.close()
+        mode, low, high = get_mode(rn, rbins)
+        outfile.write("{0} {1} {2} \n".format(mode, low, high))
  
-    # --------------------------
-    # For debugging, plot signal energy
-    # ----------------------
-    plt.figure()
-    e_freq, energy = np.loadtxt('rho.txt', unpack=True)
-    plt.plot(e_freq,energy)
-    plt.grid()
-    plt.ylabel('Signal energy')
-    plt.xlabel('Frequency')
-    plt.savefig('sig_energy.png')
-
-
-    # --------------------------
-    # For debugging, plot signal energy (time domain)
-    # ----------------------
-    plt.figure()
-    t, e = np.loadtxt('t_rho.txt', unpack=True)
-    plt.plot(t,e)
-    plt.grid()
-    plt.ylabel('Signal energy')
-    plt.xlabel('Time')
-    plt.xlim(1.98, 2.02)
-    plt.savefig('t_sig_energy.png')
-
-
-    # --------------------------
-    # For debugging, plot time domain waveform
-    # ----------------------
-    plt.figure()
-    t, h = np.loadtxt('sanity_time.txt', unpack=True)
-    plt.plot(t,h)
-    plt.grid()
-    plt.xlim(1.98, 2.02)
-    plt.ylabel('Waveform')
-    plt.xlabel('Time')
-    plt.savefig('hoft.png')
-   
+        # --------------------------
+        # For debugging, plot signal energy
+        # ----------------------
+        plt.figure()
+        e_freq, energy = np.loadtxt('rho.txt', unpack=True)
+        plt.plot(e_freq,energy)
+        plt.grid()
+        plt.ylabel('Signal energy')
+        plt.xlabel('Frequency')
+        plt.savefig('sig_energy.png')
+
+
+        # --------------------------
+        # For debugging, plot signal energy (time domain)
+        # ----------------------
+        plt.figure()
+        t, e = np.loadtxt('t_rho.txt', unpack=True)
+        plt.plot(t,e)
+        plt.grid()
+        plt.ylabel('Signal energy')
+        plt.xlabel('Time')
+        plt.xlim(1.98, 2.02)
+        plt.savefig('t_sig_energy.png')
+
+
+        # --------------------------
+        # For debugging, plot time domain waveform
+        # ----------------------
+        plt.figure()
+        t, h = np.loadtxt('sanity_time.txt', unpack=True)
+        plt.plot(t,h)
+        plt.grid()
+        plt.xlim(1.98, 2.02)
+        plt.ylabel('Waveform')
+        plt.xlabel('Time')
+        plt.savefig('hoft.png')
+        
+        # -- End the loop over white/colored moments
+        outfile.close()
 
     # -----------------------
     # Plot the median PSDs
     # -----------------------
-
     fileList = glob.glob('waveforms/'+str(jobname)+'clean_1_psd_*.dat.{0}'.format(ifo))
 
     if not len(fileList) == 0:
@@ -714,9 +758,9 @@ for ifo in ifoList:
 #    web.write("<img src='./glitch_{0}.png' alt='glitch_{0}.png' width=500>".format(ifo))
 
 # -- Write out spectrogram
-web.write("<h3>Spectrogram of median reconstructed signal model waveform</h3>")
-for ifo in ifoList:
-    web.write("<img src='./self_spec_{0}.png' alt='self_spec_{0}.png' width=500>".format(ifo) )
+# web.write("<h3>Spectrogram of median reconstructed signal model waveform</h3>")
+# for ifo in ifoList:
+#    web.write("<img src='./self_spec_{0}.png' alt='self_spec_{0}.png' width=500>".format(ifo) )
 
 # -- Write out the PSDs
 web.write("<h3>PSDs and Reconstructed PSDs</h3>")