Commit 80b7bff8 authored by Tyson Littenberg's avatar Tyson Littenberg
Browse files

review ready postprocessing


git-svn-id: https://svn.ligo.caltech.edu/svn/bayeswave/trunk@494 c56465c9-8126-4a4f-9d7d-ac845eff4865
parent c4a898d6
......@@ -176,7 +176,7 @@ def readbwb():
print '\nThis run was executed with the --signalOnly flag\n'
restrictModel = 'signal'
elif arg=='--inj':
injFlag = True
injFlag = False #True
mdc = True
elif arg=='--MDC-cache':
mdc = True
......@@ -733,7 +733,6 @@ def mode_values(moment, ifo, mdc, momentsPlus, momentsCross, ifoList, tablesDir,
def plot_likelihood_1(modelList, plotsDir):
plt.figure()
plt.title('Likelihood')
# plt.xlim([1e-4, 1])
for mod in modelList:
if mod == 'glitch':
colour = 'blue'
......@@ -755,8 +754,9 @@ def plot_likelihood_1(modelList, plotsDir):
except:
continue
error = np.zeros(likehood.shape)
plt.semilogx(temp, likehood, label=mod, linewidth=2, color=colour)
plt.errorbar(temp, likehood, 2*error, color=colour)
mask=(temp>=1e-3)
plt.semilogx(temp[mask], likehood[mask], label=mod, linewidth=2, color=colour)
plt.errorbar(temp[mask], likehood[mask], 2*error[mask], color=colour)
plt.xlabel( '1/Temp' )
plt.ylabel('log(L)')
plt.grid()
......@@ -873,10 +873,10 @@ def plot_model_dims(modelList, ifoList, ifoNames, plotsDir):
fig, (ax1, ax2) = plt.subplots(2, sharex=True)
ax1.set_title('Model Dimension Histogram')
if not signalChains == []:
n,bins,patches = ax1.hist(signalChains[2], bins=range(int(min(signalChains[2])), int(max(signalChains[2]))+1, 1), histtype='bar', color='m', log=True)
n,bins,patches = ax1.hist(signalChains[2], bins=range(int(min(signalChains[2])), int(max(signalChains[2]))+1, 1), histtype='bar', color='m', log=False)
if not glitchChains == []:
data = np.dstack(glitchChains[3:3+len(ifoList)])[0]
n,bins,patches = ax2.hist(data, bins=range(int(data.min()), int(data.max())+1, 1), label=ifoNames, histtype='bar', log=True)
n,bins,patches = ax2.hist(data, bins=range(int(data.min()), int(data.max())+1, 1), label=ifoNames, histtype='bar', log=False)
# -- Legend placement outside the plot
box = ax1.get_position()
ax1.set_position([box.x0, box.y0, box.width*0.8, box.height])
......
......@@ -17,12 +17,16 @@ import sky_area.sky_area_clustering as sky
import getopt
from readbwb import BwbParams
import lalinference.fits as lf
import lalinference.plot as lp
import lalinference.cmap
# -------------------------------------------
# Module to make skymaps, skyview webpage, etc.
# Works with single output directory from
# BayesWaveBurst
# -------------------------------------------
def make_skyview(directory='.', mdc=None, NSIDE=256, ra=None, dec=None, results=None):
def make_skyview(directory='.', mdc=None, NSIDE=64, ra=None, dec=None, results=None, npost=5000):
# -- Close any open figures
plt.close('all')
......@@ -91,7 +95,7 @@ def make_skyview(directory='.', mdc=None, NSIDE=256, ra=None, dec=None, results=
print "Total samples are {0}".format(ralist.size)
# -- Remove burn in samples
burnin = ralist.size/2
burnin = ralist.size/4
ralist = ralist[burnin:]
sin_dec = sin_dec[burnin:]
print "After removing burn-in samples are {0}".format(ralist.size)
......@@ -99,30 +103,12 @@ def make_skyview(directory='.', mdc=None, NSIDE=256, ra=None, dec=None, results=
declist = np.arcsin(sin_dec)
thetalist = np.pi/2.0 - declist
# -- Determin auto-correlation length
print "The size of the RA list is {0}".format(ralist.size)
tau, mean, sigma = acor.acor(ralist)
acl = int(tau)
print "The ACL is {0}".format(acl)
# -- Thin the chains
ralist = ralist[::acl]
declist = declist[::acl]
thetalist = thetalist[::acl]
radec = np.column_stack((ralist, declist))
# -- Make the KDE map
redo = True
count = 0
while redo:
try:
#kde = sky.ClusteredKDEPosterior(radec)
kde = sky.ClusteredSkyKDEPosterior(radec)
redo = False
except:
count += 1
print "KDE failure number {0}".format(count)
print "Will attempt to try again"
if radec.shape[0] > npost:
radec = np.random.permutation(radec)[:npost,:]
kde = sky.ClusteredSkyKDEPosterior(radec)
# -- Get the sidereal time
print "Finding sidereal time"
......@@ -150,7 +136,7 @@ def make_skyview(directory='.', mdc=None, NSIDE=256, ra=None, dec=None, results=
injra = params.ra
injtheta = np.pi/2 - injdec
mdc = True
# -- Make plots directory, if needed
plotsDir = './plots'
if not os.path.exists(plotsDir):
......@@ -178,17 +164,25 @@ def make_skyview(directory='.', mdc=None, NSIDE=256, ra=None, dec=None, results=
rarad = rarad + 2*np.pi
# Shift the declination to match hp.projscatter conventions
decrad = np.pi*0.5 - decrad
# -- Plot the skymap and injection location
skymap = kde.as_healpix(NSIDE, nest=True)
hp.mollview(skymap, title='Skymap from KDE clusters', nest=True)
#fig = plt.figure(frameon=False)
fig = plt.figure(figsize=(8,6), frameon=False)
ax = plt.subplot(111, projection='astro mollweide')
ax.cla()
ax.grid()
lp.healpix_heatmap(skymap, nest=True, vmin=0.0, vmax=np.max(skymap), cmap=plt.get_cmap('cylon'))
if (ra is not None and dec is not None):
hp.projscatter(decrad, rarad, s=200, c='yellow', marker=(5,1), alpha=0.5)
plt.plot(rarad, dec, 'kx', ms=30, mew=1)
if mdc is not None:
hp.projscatter(injtheta, injra, s=200, c='yellow', marker=(5,1), alpha=0.5)
plt.plot(injra, injdec, 'kx', ms=30, mew=1)
plt.savefig(plotsDir+'/skymap.png')
plt.close()
hp.write_map("skymap_{0}.fits".format(gps), skymap, nest=True)
lf.write_sky_map('skymap_{0}.fits'.format(gps), skymap, nest=True, gps_time=trigtime)
# -- Calculate the 50 and 90% credible intervals
sq_deg = (180.0/np.pi)**2
......@@ -210,25 +204,12 @@ def make_skyview(directory='.', mdc=None, NSIDE=256, ra=None, dec=None, results=
searcharea = 0
pixarea = hp.nside2pixarea(NSIDE, degrees=True)
print("here")
if results is not None:
for name in ['injcontour', 'area50', 'area90', 'searcharea']:
if name not in results: results[name] = []
results[name].append(eval(name))
# -- Compute autocorrelation functions
acf = bayespputils.autocorrelation(ralist)
plt.figure()
plt.plot(acf,'o', markersize=2)
plt.title('RA autocorrelation', fontsize=14)
plt.savefig(plotsDir+'/ra_acf.png')
plt.close()
# -- Try scatter plot of ra
plt.figure()
plt.plot(ralist, 'o', markersize=2)
plt.title('RA chain: {0} samples'.format(ralist.size), fontsize=14)
plt.savefig(plotsDir+'/ra_chain.png')
# -- Make an output file with key statistics
outfile = open('skystats.txt', 'w')
......@@ -236,38 +217,6 @@ def make_skyview(directory='.', mdc=None, NSIDE=256, ra=None, dec=None, results=
outfile.write("{0} {1} {2} {3}".format(area50, area90, searcharea, injcontour))
outfile.close()
# -- Since we've done all the sky stuff, let's try making the skyview web page right here
skyview = open('skyview.html','w')
skyview.write('<html><head><title>Skyview</title></head><body>\n')
skyview.write('<h1>Sky View for trigger at {0}</h1>\n'.format(trigtime))
skyview.write('<h2>Links</h2>\n')
skyview.write('<a href="./index.html">Event Overview</a> for this trigger<br/>\n')
skyview.write('<a href="../../top.html">Top level page</a> for this injection set<br/>\n')
skyview.write('<h2>Statistics</h2>\n')
skyview.write('The area per pixel is {0:.2f} sq. degrees.<br/>\n'.format(pixarea))
skyview.write('The area of the 50% region is {0:.2f} sq. degrees.<br/>\n'.format(area50))
skyview.write('The area of the 90% region is {0:.2f} sq. degrees.<br/>\n'.format(area90))
skyview.write('The searched area is {0:.2f} sq. degrees.<br/>\n'.format(searcharea))
skyview.write('The injection was found at the {0:.2f} percent area contour.<br/>\n'.format(injcontour * 100.0))
skyview.write('<h2>Plots</h2>\n')
skyview.write('<img src="'+plotsDir+'/skymap.png" width=600>\n')
skyview.write('<img src="'+plotsDir+'/ra_chain.png" width=600>\n')
skyview.write('<img src="'+plotsDir+'/ra_acf.png" width=600>\n')
skyview.write('</body></html>')
skyview.close()
# -- Let's also save the main results on a txt file for megaplot.py to read in
datoutput = open('megasky_results.dat', 'w')
datoutput.write('# Megasky main results\n')
datoutput.write('# Rows are: pixarea [sq. degrees.], area50 [sq. degrees.], area90 [sq. degrees.], searcharea [sq. degrees.], injcontour%\n')
datoutput.write('{0}\n'.format(pixarea))
datoutput.write('{0}\n'.format(area50))
datoutput.write('{0}\n'.format(area90))
datoutput.write('{0}\n'.format(searcharea))
datoutput.write('{0}'.format(format(injcontour * 100.0)))
datoutput.close
os.chdir(topdir)
# -- Write main script for command line running
if __name__ == "__main__":
......@@ -277,7 +226,7 @@ if __name__ == "__main__":
# -- Set default argument values
directory='.'
mdc=None
NSIDE=256
NSIDE=64
ra=None
dec=None
geo=False
......@@ -311,6 +260,6 @@ if __name__ == "__main__":
dec = float(arg)
if opt == '--geo':
geo = 'True'
make_skyview(directory, mdc, NSIDE, ra, dec, results)
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment