Contour regions on PP plots

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

postprocess/s6/plot_evidence.py

@@ -15,6 +15,7 @@ import acor
 from pylal import bayespputils
 import sky_area.sky_area_clustering as sky
 import skyview
+from scipy import stats
 
 # --------------
 # Parameters
@@ -30,7 +31,7 @@ plt.close('all')
 # Read S6 MDC log file
 # -----------------
 print "Reading MDC log (this may take a few minutes) ..."
-mdc = ft.Mdc(filename=mdcfile)
+# mdc = ft.Mdc(filename=mdcfile)
 print "I read the mdc log!"
 
 
@@ -149,7 +150,23 @@ for job in dirList:
     # --------------------------------
     # Process skymap information
     # --------------------------------
-    skyview.make_skyview(mdc=mdc, NSIDE=NSIDE, results=skydict, sim=True)
+    # -- Check if skystats exists
+    if os.path.isfile('skystats.txt'):
+        # area50, area90, searcharea, injcontour = np.loadtxt('skystats.txt', unpack=True)
+        print "Found skystats file!"
+        # skydata = np.recfromtxt('skystats.txt', names='True')
+        # for key in ['area50', 'area90', 'searcharea', 'injcontour']:
+        #    skydict[key].append(skydata[key])
+
+        skydata = open('skystats.txt', 'r')
+        keylist = skydata.readline().split()[1:]
+        valuelist = skydata.readline().split()
+        for key, value in zip(keylist, valuelist):
+            if key not in skydict: skydict[key] = []
+            skydict[key].append(float(value))
+    else:
+        print "No skystats found :(  Running skyview.  You have time for a coffee break ..."
+        skyview.make_skyview(mdc=mdc, NSIDE=NSIDE, results=skydict, sim=True)
 
     # -- Return to top level directory, continue through list of jobs
     plt.close('all')
@@ -162,15 +179,52 @@ plt.rcParams['figure.figsize'] = 8, 6
 # Make PP plot
 # ------------
 injconArray = np.array(skydict['injcontour'])
+print injconArray
 xtics = np.arange(0,1.01,0.01)
 numFound = np.array( [ (injconArray <= xvalue).sum() for xvalue in xtics] )
 fracFound = 1.0*numFound / injconArray.size
+
 plt.close('all')
 plt.figure()
+
+# -- Make 1 sigma confidence region
+print "attempting to calculate PP plots"
+bi = stats.binom
+N = len(injconArray)
+
+# -- 1 sigma
+con = 0.68
+edge = (1 - con)/2.0
+uppery = bi.ppf(1-edge, N, xtics)
+lowery = bi.ppf(edge, N, xtics)
+uppery[0]=0
+lowery[0]=0
+plt.fill_between(xtics, lowery/N, uppery/N, facecolor='gray', alpha=0.2)
+
+# -- 2 sigma
+con = 0.95
+edge = (1 - con)/2.0
+uppery = bi.ppf(1-edge, N, xtics)
+lowery = bi.ppf(edge, N, xtics)
+uppery[0]=0
+lowery[0]=0
+plt.fill_between(xtics, lowery/N, uppery/N, facecolor='gray', alpha=0.2)
+
+
+# -- 3 sigma
+con = 0.997
+edge = (1 - con)/2.0
+uppery = bi.ppf(1-edge, N, xtics)
+lowery = bi.ppf(edge, N, xtics)
+uppery[0]=0
+lowery[0]=0
+plt.fill_between(xtics, lowery/N, uppery/N, facecolor='gray', alpha=0.2)
+
+
 plt.plot(xtics, fracFound)
 plt.xlabel('Confidence Region')
 plt.ylabel('Fraction of injections found')
-plt.plot( [0, 1], [0, 1], 'k--' )
+plt.plot( [0, 1], [0, 1], 'k--', linewidth=2 )
 plt.savefig('pp.png')
 plt.close()