Commit 8475e4ac authored by Jonah Kanner's avatar Jonah Kanner 🤓
Browse files

megasky now removes burnin

git-svn-id: https://svn.ligo.caltech.edu/svn/bayeswave/trunk@243 c56465c9-8126-4a4f-9d7d-ac845eff4865
parent 9ef7b471
......@@ -1123,8 +1123,12 @@ for mod in modelList:
sig_gl, sig_noise, err_sig_gl, err_sig_noise = plot_evidence(jobName, plotsDir)
# -- Plot Temp Chains
plot_likelihood_1(restrictModel, plotsDir)
plot_likelihood_2(restrictModel, plotsDir)
# WARNING: This feature appears to be broken because of a new file format!
try:
plot_likelihood_1(restrictModel, plotsDir)
plot_likelihood_2(restrictModel, plotsDir)
except:
print "Failed to plot temp vs. likelihood. This may be due to a change in file format for files like signal_evidence.dat"
# -- Plot the evolution of the dimensions of each model
plot_model_dims(modelList, ifoList, ifoNames, plotsDir)
......
import numpy as np
# ---------------
# Script to read
# BWB parameters
# --------------
# -- Read the run file
class BwbParams:
def __init__(self):
# -- Read run file
bayeswaverunfile = glob.glob('*bayeswave.run')
if not len(bayeswaverunfile):
sys.exit("\n *bayeswave.run not found! \n")
bayeswave = open(bayeswaverunfile[0], 'r')
for line in runfile:
spl = line.split()
if len(spl) > 2:
if line.split()[0].strip() == "Command":
spl = line.split()
try:
indx = spl.index('--MDC-prefactor')
self.scale = spl[indx+1]
except:
self.scale = None
try:
indx = spl.index('--trigtime')
self.gps = spl[indx+1]
except:
self.gps = None
print "WARNING! Failed to read GPS time"
if len(spl) > 3:
if spl[2] == 'detectors':
self.numifo = spl[3]
runfile.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]))
infile.close()
self.snr = np.array(snr)
import numpy as np
import glob
# ---------------
# Script to read
......@@ -12,8 +13,21 @@ import numpy as np
class BwbParams:
def __init__(self):
# -- Determine name of job
bayeswaverunfile = glob.glob('*bayeswave.run')
if not len(bayeswaverunfile):
sys.exit("\n bayeswave.run not found! \n")
bayeswaverunfile = bayeswaverunfile[0]
jobName = bayeswaverunfile.split('_')
del jobName[-1]
if not len(jobName) == 0:
self.jobname = '_'.join(jobName)+'_'
else:
self.jobname = ''
# -- Read run file
runfile = open('bayeswave.run', 'r')
runfile = open(bayeswaverunfile, 'r')
for line in runfile:
spl = line.split()
......
......@@ -15,6 +15,7 @@ import acor
from pylal import bayespputils
import sky_area.sky_area_clustering as sky
import getopt
from readbwb import BwbParams
# -------------------------------------------
# Module to make skymaps, skyview webpage, etc.
......@@ -73,19 +74,33 @@ def make_skyview(directory='.', mdc=None, NSIDE=256, results=None):
# --------------------------------
# Create skymap summary statistics
# --------------------------------
# -- Get run name
params = BwbParams()
jobname = params.jobname
# -- Input skymap data
# num, L, ralist, sin_dec, psi, e, dA, dphi, dt = np.loadtxt(filename, unpack=True)
print "Extracting RA/DEC samples"
filename = './chains/signal_extchain.dat.0'
filename = './chains/' + jobname + 'signal_extchain.dat.0'
data = np.loadtxt(filename, unpack=True)
ralist = data[2]
sin_dec = data[3]
print "Total samples are {0}".format(ralist.size)
# -- Remove burn in samples
burnin = ralist.size/2
ralist = ralist[burnin:]
sin_dec = sin_dec[burnin:]
print "After removing burn-in samples are {0}".format(ralist.size)
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]
......
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