From ad1d6a2466793f875e8d8f9daa22d77972e6a46f Mon Sep 17 00:00:00 2001
From: Gregory Ashton <gregory.ashton@ligo.org>
Date: Tue, 8 Nov 2016 16:20:16 +0100
Subject: [PATCH] Several minor improvements
1) Adds function to produce a .loudest file
2) Improves log output of data, including a histogram of SFT timestamps
3) Generalised automated estimation of the F0s when many SFTs are loaded
4) Adds tref to the par file
---
pyfstat.py | 46 +++++++++++++++++++++++++++++++++++++---------
1 file changed, 37 insertions(+), 9 deletions(-)
diff --git a/pyfstat.py b/pyfstat.py
index ca8ef3f..511fef4 100755
--- a/pyfstat.py
+++ b/pyfstat.py
@@ -188,6 +188,20 @@ class BaseSearchClass(object):
post_theta_at_ith_glitch, self.tref - tbounds[i+1]))
return thetas
+ def generate_loudest(self):
+ params = read_par(self.label, self.outdir)
+ for key in ['Alpha', 'Delta', 'F0', 'F1']:
+ if key not in params:
+ params[key] = self.theta_prior[key]
+ cmd = ('lalapps_ComputeFstatistic_v2 -a {} -d {} -f {} -s {} -D "{}"'
+ ' --refTime={} --outputLoudest="{}/{}.loudest" '
+ '--minStartTime={} --maxStartTime={}').format(
+ params['Alpha'], params['Delta'], params['F0'],
+ params['F1'], self.sftfilepath, params['tref'],
+ self.outdir, self.label, self.minStartTime,
+ self.maxStartTime)
+ subprocess.call([cmd], shell=True)
+
class ComputeFstat(object):
""" Base class providing interface to `lalpulsar.ComputeFstat` """
@@ -255,10 +269,23 @@ class ComputeFstat(object):
SFTCatalog = lalpulsar.SFTdataFind(self.sftfilepath, constraints)
names = list(set([d.header.name for d in SFTCatalog.data]))
SFT_timestamps = [d.header.epoch for d in SFTCatalog.data]
- logging.info(
- 'Loaded {} data files from detectors {} spanning {} to {}'.format(
- len(SFT_timestamps), names, int(SFT_timestamps[0]),
- int(SFT_timestamps[-1])))
+ try:
+ from bashplotlib.histogram import plot_hist
+ print('Data timestamps histogram:')
+ plot_hist(SFT_timestamps, height=5, bincount=50)
+ except IOError:
+ pass
+ if len(names) == 0:
+ raise ValueError('No data loaded.')
+ logging.info('Loaded {} data files from detectors {}'.format(
+ len(SFT_timestamps), names))
+ logging.info('Data spans from {} ({}) to {} ({})'.format(
+ int(SFT_timestamps[0]),
+ subprocess.check_output('lalapps_tconvert {}'.format(
+ int(SFT_timestamps[0])), shell=True).rstrip('\n'),
+ int(SFT_timestamps[-1]),
+ subprocess.check_output('lalapps_tconvert {}'.format(
+ int(SFT_timestamps[-1])), shell=True)).rstrip('\n'))
logging.info('Initialising ephems')
ephems = lalpulsar.InitBarycenter(self.earth_ephem, self.sun_ephem)
@@ -303,11 +330,11 @@ class ComputeFstat(object):
FstatOAs.injectSources = lalpulsar.FstatOptionalArgsDefaults.injectSources
if self.minCoverFreq is None or self.maxCoverFreq is None:
- fA = SFTCatalog.data[0].header.f0
- numBins = SFTCatalog.data[0].numBins
- fB = fA + (numBins-1)*SFTCatalog.data[0].header.deltaF
- self.minCoverFreq = fA + 0.5
- self.maxCoverFreq = fB - 0.5
+ fAs = [d.header.f0 for d in SFTCatalog.data]
+ fBs = [d.header.f0 + (d.numBins-1)*d.header.deltaF
+ for d in SFTCatalog.data]
+ self.minCoverFreq = np.min(fAs) + 0.5
+ self.maxCoverFreq = np.max(fBs) - 0.5
logging.info('Min/max cover freqs not provided, using '
'{} and {}, est. from SFTs'.format(
self.minCoverFreq, self.maxCoverFreq))
@@ -1488,6 +1515,7 @@ class MCMCSearch(BaseSearchClass):
filename = '{}/{}.par'.format(self.outdir, self.label)
with open(filename, 'w+') as f:
f.write('MaxtwoF = {}\n'.format(max_twoF))
+ f.write('tref = {}\n'.format(self.tref))
if hasattr(self, 'theta0_index'):
f.write('theta0_index = {}\n'.format(self.theta0_idx))
if method == 'med':
--
GitLab