From 2834bb6e91ff85a295c54143bb922259819afac1 Mon Sep 17 00:00:00 2001
From: David Keitel <david.keitel@ligo.org>
Date: Wed, 16 Aug 2017 12:48:10 +0100
Subject: [PATCH] gridsearch: optionally write full transient-Fstat maps to
file(s)
-new output: full (t0,tau) Fstat map,
if bands in those params are searched
-one file per Doppler grid point
---
.../short_transient_search_gridded.py | 6 ++++--
pyfstat/core.py | 13 +++++++------
pyfstat/grid_based_searches.py | 17 ++++++++++++++---
3 files changed, 25 insertions(+), 11 deletions(-)
diff --git a/examples/transient_examples/short_transient_search_gridded.py b/examples/transient_examples/short_transient_search_gridded.py
index fd9b83b..64e6423 100644
--- a/examples/transient_examples/short_transient_search_gridded.py
+++ b/examples/transient_examples/short_transient_search_gridded.py
@@ -35,7 +35,8 @@ search1 = pyfstat.GridSearch(
sftfilepattern=os.path.join(datadir,'*simulated_transient_signal*sft'),
F0s=F0s, F1s=F1s, F2s=F2s, Alphas=Alphas, Deltas=Deltas, tref=tref,
minStartTime=minStartTime, maxStartTime=maxStartTime,
- BSGL=False)
+ BSGL=False,
+ outputTransientFstatMap=True)
search1.run()
search1.print_max_twoF()
@@ -49,7 +50,8 @@ search2 = pyfstat.GridSearch(
F0s=F0s, F1s=F1s, F2s=F2s, Alphas=Alphas, Deltas=Deltas, tref=tref,
minStartTime=minStartTime, maxStartTime=maxStartTime,
transientWindowType='rect', t0Band=Tspan-2*Tsft, tauBand=Tspan,
- BSGL=False)
+ BSGL=False,
+ outputTransientFstatMap=True)
search2.run()
search2.print_max_twoF()
diff --git a/pyfstat/core.py b/pyfstat/core.py
index 5d06781..e336acc 100755
--- a/pyfstat/core.py
+++ b/pyfstat/core.py
@@ -642,7 +642,6 @@ class ComputeFstat(BaseSearchClass):
asini=None, period=None, ecc=None, tp=None,
argp=None):
""" Returns twoF or ln(BSGL) fully-coherently at a single point """
-
self.PulsarDopplerParams.fkdot = np.array([F0, F1, F2, 0, 0, 0, 0])
self.PulsarDopplerParams.Alpha = Alpha
self.PulsarDopplerParams.Delta = Delta
@@ -681,10 +680,11 @@ class ComputeFstat(BaseSearchClass):
# F-stat computation
self.windowRange.tau = int(2*self.Tsft)
- FS = lalpulsar.ComputeTransientFstatMap(
+ self.FstatMap = lalpulsar.ComputeTransientFstatMap(
self.FstatResults.multiFatoms[0], self.windowRange, False)
+ F_mn = self.FstatMap.F_mn.data
- twoF = 2*np.max(FS.F_mn.data)
+ twoF = 2*np.max(F_mn)
if self.BSGL is False:
if np.isnan(twoF):
return 0
@@ -705,7 +705,7 @@ class ComputeFstat(BaseSearchClass):
# to return BSGL
# FIXME: should we instead compute BSGL over the whole F_mn
# and return the maximum of that?
- idx_maxTwoF = np.argmax(FS.F_mn.data)
+ idx_maxTwoF = np.argmax(F_mn)
self.twoFX[0] = 2*FS0.F_mn.data[idx_maxTwoF]
self.twoFX[1] = 2*FS1.F_mn.data[idx_maxTwoF]
@@ -751,10 +751,11 @@ class ComputeFstat(BaseSearchClass):
self.transientWindowType = 'none'
self.init_computefstatistic_single_point()
for tau in taus:
- twoFs.append(self.get_fullycoherent_twoF(
+ detstat = self.get_fullycoherent_twoF(
tstart=tstart, tend=tstart+tau, F0=F0, F1=F1, F2=F2,
Alpha=Alpha, Delta=Delta, asini=asini, period=period, ecc=ecc,
- tp=tp, argp=argp))
+ tp=tp, argp=argp)
+ twoFs.append(detstat)
return taus, np.array(twoFs)
diff --git a/pyfstat/grid_based_searches.py b/pyfstat/grid_based_searches.py
index f0e67a3..4617936 100644
--- a/pyfstat/grid_based_searches.py
+++ b/pyfstat/grid_based_searches.py
@@ -32,7 +32,8 @@ class GridSearch(BaseSearchClass):
nsegs=1, BSGL=False, minCoverFreq=None, maxCoverFreq=None,
detectors=None, SSBprec=None, injectSources=None,
input_arrays=False, assumeSqrtSX=None,
- transientWindowType=None, t0Band=None, tauBand=None):
+ transientWindowType=None, t0Band=None, tauBand=None,
+ outputTransientFstatMap=False):
"""
Parameters
----------
@@ -59,6 +60,9 @@ class GridSearch(BaseSearchClass):
and tau in (2*Tsft,2*Tsft+tauBand).
if =0, only compute CW Fstat with t0=minStartTime,
tau=maxStartTime-minStartTime.
+ outputTransientFstatMap: bool
+ if true, write output files for (t0,tau) Fstat maps
+ (one file for each doppler grid point!)
For all other parameters, see `pyfstat.ComputeFStat` for details
"""
@@ -157,8 +161,15 @@ class GridSearch(BaseSearchClass):
data = []
for vals in tqdm(self.input_data):
- FS = self.search.get_det_stat(*vals)
- data.append(list(vals) + [FS])
+ detstat = self.search.get_det_stat(*vals)
+ windowRange = getattr(self.search, 'windowRange', None)
+ FstatMap = getattr(self.search, 'FstatMap', None)
+ data.append(list(vals) + [detstat])
+ if self.outputTransientFstatMap and self.transientWindowType:
+ tCWfile = os.path.splitext(self.out_file)[0]+'_tCW_%.16f_%.16f_%.16f_%.16g_%.16g.dat' % (vals[2],vals[5],vals[6],vals[3],vals[4]) # freq alpha delta f1dot f2dot
+ fo = lal.FileOpen(tCWfile, 'w')
+ lalpulsar.write_transientFstatMap_to_fp ( fo, FstatMap, windowRange, None )
+ del fo # instead of lal.FileClose() which is not SWIG-exported
data = np.array(data, dtype=np.float)
if return_data:
--
GitLab