diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 316093ccc974d296bacbb4833b91cb3a30830f21..a1066d505b40b22fd13392a3278068881d86fd80 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -1,11 +1,11 @@
test_app:
script:
- - . /home/greg/lalsuite-install/etc/lalapps-user-env.sh
- - /home/greg/anaconda2/bin/python tests.py Writer
- - /home/greg/anaconda2/bin/python tests.py par
- - /home/greg/anaconda2/bin/python tests.py BaseSearchClass
- - /home/greg/anaconda2/bin/python tests.py ComputeFstat
- - /home/greg/anaconda2/bin/python tests.py SemiCoherentSearch
- - /home/greg/anaconda2/bin/python tests.py SemiCoherentGlitchSearch
- - /home/greg/anaconda2/bin/python tests.py MCMCSearch
- - /home/greg/anaconda2/bin/python tests.py GridSearch
+ - . /home/user1/lalsuite-install/etc/lalapps-user-env.sh
+ - /home/user1/anaconda2/bin/python tests.py Writer
+ - /home/user1/anaconda2/bin/python tests.py par
+ - /home/user1/anaconda2/bin/python tests.py BaseSearchClass
+ - /home/user1/anaconda2/bin/python tests.py ComputeFstat
+ - /home/user1/anaconda2/bin/python tests.py SemiCoherentSearch
+ - /home/user1/anaconda2/bin/python tests.py SemiCoherentGlitchSearch
+ - /home/user1/anaconda2/bin/python tests.py MCMCSearch
+ - /home/user1/anaconda2/bin/python tests.py GridSearch
diff --git a/examples/transient_examples/short_transient_search_gridded.py b/examples/transient_examples/short_transient_search_gridded.py
index 64e642342f40b088a5c68c1d44c7b69e68be20e9..dea9c26745e968bc6af825b61aaeb78b72fef22a 100644
--- a/examples/transient_examples/short_transient_search_gridded.py
+++ b/examples/transient_examples/short_transient_search_gridded.py
@@ -35,8 +35,7 @@ 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,
- outputTransientFstatMap=True)
+ BSGL=False)
search1.run()
search1.print_max_twoF()
@@ -44,7 +43,7 @@ search1.plot_1D(xkey='F0',
xlabel='freq [Hz]', ylabel='$2\mathcal{F}$')
print('with t0,tau bands:')
-search2 = pyfstat.GridSearch(
+search2 = pyfstat.TransientGridSearch(
label='tCW', outdir=datadir,
sftfilepattern=os.path.join(datadir,'*simulated_transient_signal*sft'),
F0s=F0s, F1s=F1s, F2s=F2s, Alphas=Alphas, Deltas=Deltas, tref=tref,
diff --git a/pyfstat/__init__.py b/pyfstat/__init__.py
index 5c447084af048abdbb789435a6a48e65511cc713..3516c4a8fdffbedb1697a40dbed4ba2f37dd835c 100644
--- a/pyfstat/__init__.py
+++ b/pyfstat/__init__.py
@@ -3,4 +3,4 @@ from __future__ import division as _division
from .core import BaseSearchClass, ComputeFstat, SemiCoherentSearch, SemiCoherentGlitchSearch
from .make_sfts import Writer, GlitchWriter, FrequencyModulatedArtifactWriter, FrequencyAmplitudeModulatedArtifactWriter
from .mcmc_based_searches import MCMCSearch, MCMCGlitchSearch, MCMCSemiCoherentSearch, MCMCFollowUpSearch, MCMCTransientSearch
-from .grid_based_searches import GridSearch, GridUniformPriorSearch, GridGlitchSearch, FrequencySlidingWindow, DMoff_NO_SPIN, SliceGridSearch
+from .grid_based_searches import GridSearch, GridUniformPriorSearch, GridGlitchSearch, FrequencySlidingWindow, DMoff_NO_SPIN, SliceGridSearch, TransientGridSearch
diff --git a/pyfstat/core.py b/pyfstat/core.py
index 1e7452257b53260eff2ac69c34024f73cc3b579e..d0731a87a4032db9b397180ac765b22a0fa0f7b9 100755
--- a/pyfstat/core.py
+++ b/pyfstat/core.py
@@ -331,6 +331,7 @@ class ComputeFstat(BaseSearchClass):
def __init__(self, tref, sftfilepattern=None, minStartTime=None,
maxStartTime=None, binary=False, BSGL=False,
transientWindowType=None, t0Band=None, tauBand=None,
+ dt0=None, dtau=None,
detectors=None, minCoverFreq=None, maxCoverFreq=None,
injectSources=None, injectSqrtSX=None, assumeSqrtSX=None,
SSBprec=None):
@@ -359,6 +360,9 @@ class ComputeFstat(BaseSearchClass):
and tau in (2*Tsft,2*Tsft+tauBand).
if =0, only compute CW Fstat with t0=minStartTime,
tau=maxStartTime-minStartTime.
+ dt0, dtau: int
+ grid resolutions in transient start-time and duration,
+ both default to Tsft
detectors : str
Two character reference to the data to use, specify None for no
contraint. If multiple-separate by comma.
@@ -615,27 +619,38 @@ class ComputeFstat(BaseSearchClass):
.format(self.transientWindowType,
', '.join(transientWindowTypes)))
+ # default spacing
self.Tsft = int(1.0/SFTCatalog.data[0].header.deltaF)
- if self.t0Band is None:
+ self.windowRange.dt0 = self.Tsft
+ self.windowRange.dtau = self.Tsft
+
+ # special treatment of window_type = none ==> replace by rectangular window spanning all the data
+ if self.windowRange.type == lalpulsar.TRANSIENT_NONE:
+ self.windowRange.t0 = int(self.minStartTime)
self.windowRange.t0Band = 0
- self.windowRange.dt0 = 1
- else:
- if not isinstance(self.t0Band, int):
- logging.warn('Casting non-integer t0Band={} to int...'
- .format(self.t0Band))
- self.t0Band = int(self.t0Band)
- self.windowRange.t0Band = self.t0Band
- self.windowRange.dt0 = self.Tsft
- if self.tauBand is None:
+ self.windowRange.tau = int(self.maxStartTime-self.minStartTime)
self.windowRange.tauBand = 0
- self.windowRange.dtau = 1
- else:
- if not isinstance(self.tauBand, int):
- logging.warn('Casting non-integer tauBand={} to int...'
- .format(self.tauBand))
- self.tauBand = int(self.tauBand)
- self.windowRange.tauBand = self.tauBand
- self.windowRange.dtau = self.Tsft
+ else: # user-set bands and spacings
+ if self.t0Band is None:
+ self.windowRange.t0Band = 0
+ else:
+ if not isinstance(self.t0Band, int):
+ logging.warn('Casting non-integer t0Band={} to int...'
+ .format(self.t0Band))
+ self.t0Band = int(self.t0Band)
+ self.windowRange.t0Band = self.t0Band
+ if self.dt0:
+ self.windowRange.dt0 = self.dt0
+ if self.tauBand is None:
+ self.windowRange.tauBand = 0
+ else:
+ if not isinstance(self.tauBand, int):
+ logging.warn('Casting non-integer tauBand={} to int...'
+ .format(self.tauBand))
+ self.tauBand = int(self.tauBand)
+ self.windowRange.tauBand = self.tauBand
+ if self.dtau:
+ self.windowRange.dtau = self.dtau
def get_fullycoherent_twoF(self, tstart, tend, F0, F1, F2, Alpha, Delta,
asini=None, period=None, ecc=None, tp=None,
@@ -891,6 +906,18 @@ class ComputeFstat(BaseSearchClass):
else:
return ax
+ def write_atoms_to_file(self, fnamebase=''):
+ multiFatoms = getattr(self.FstatResults, 'multiFatoms', None)
+ if multiFatoms and multiFatoms[0]:
+ dopplerName = lalpulsar.PulsarDopplerParams2String ( self.PulsarDopplerParams )
+ #fnameAtoms = os.path.join(self.outdir,'Fstatatoms_%s.dat' % dopplerName)
+ fnameAtoms = fnamebase + '_Fstatatoms_%s.dat' % dopplerName
+ fo = lal.FileOpen(fnameAtoms, 'w')
+ lalpulsar.write_MultiFstatAtoms_to_fp ( fo, multiFatoms[0] )
+ del fo # instead of lal.FileClose() which is not SWIG-exported
+ else:
+ raise RuntimeError('Cannot print atoms vector to file: no FstatResults.multiFatoms, or it is None!')
+
class SemiCoherentSearch(ComputeFstat):
""" A semi-coherent search """
diff --git a/pyfstat/grid_based_searches.py b/pyfstat/grid_based_searches.py
index de51d0fa03936a72ca2d3fdbe1a852671b49ef0f..2126906e194d9b6784366362cd42d4104510c5cc 100644
--- a/pyfstat/grid_based_searches.py
+++ b/pyfstat/grid_based_searches.py
@@ -36,9 +36,7 @@ class GridSearch(BaseSearchClass):
Deltas, tref=None, minStartTime=None, maxStartTime=None,
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,
- outputTransientFstatMap=False):
+ input_arrays=False, assumeSqrtSX=None):
"""
Parameters
----------
@@ -55,19 +53,6 @@ class GridSearch(BaseSearchClass):
GPS seconds of the reference time, start time and end time
input_arrays: bool
if true, use the F0s, F1s, etc as is
- transientWindowType: str
- If 'rect' or 'exp', compute atoms so that a transient (t0,tau) map
- can later be computed. ('none' instead of None explicitly calls
- the transient-window function, but with the full range, for
- debugging). Currently only supported for nsegs=1.
- t0Band, tauBand: int
- if >0, search t0 in (minStartTime,minStartTime+t0Band)
- 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
@@ -91,8 +76,6 @@ class GridSearch(BaseSearchClass):
tref=self.tref, sftfilepattern=self.sftfilepattern,
minCoverFreq=self.minCoverFreq, maxCoverFreq=self.maxCoverFreq,
detectors=self.detectors,
- transientWindowType=self.transientWindowType,
- t0Band=self.t0Band, tauBand=self.tauBand,
minStartTime=self.minStartTime, maxStartTime=self.maxStartTime,
BSGL=self.BSGL, SSBprec=self.SSBprec,
injectSources=self.injectSources,
@@ -185,19 +168,7 @@ class GridSearch(BaseSearchClass):
for vals in tqdm(iterable,
total=getattr(self, 'total_iterations', None)):
detstat = self.search.get_det_stat(*vals)
- windowRange = getattr(self.search, 'windowRange', None)
- FstatMap = getattr(self.search, 'FstatMap', None)
thisCand = list(vals) + [detstat]
- if getattr(self, 'transientWindowType', None):
- if self.outputTransientFstatMap:
- 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
- Fmn = FstatMap.F_mn.data
- maxidx = np.unravel_index(Fmn.argmax(), Fmn.shape)
- thisCand += [windowRange.t0+maxidx[0]*windowRange.dt0,
- windowRange.tau+maxidx[1]*windowRange.dtau]
data.append(thisCand)
data = np.array(data, dtype=np.float)
@@ -387,6 +358,117 @@ class GridSearch(BaseSearchClass):
self.outdir, self.label, dets, type(self).__name__)
+class TransientGridSearch(GridSearch):
+ """ Gridded transient-continous search using ComputeFstat """
+
+ @helper_functions.initializer
+ def __init__(self, label, outdir, sftfilepattern, F0s, F1s, F2s, Alphas,
+ Deltas, tref=None, minStartTime=None, maxStartTime=None,
+ BSGL=False, minCoverFreq=None, maxCoverFreq=None,
+ detectors=None, SSBprec=None, injectSources=None,
+ input_arrays=False, assumeSqrtSX=None,
+ transientWindowType=None, t0Band=None, tauBand=None,
+ dt0=None, dtau=None,
+ outputTransientFstatMap=False,
+ outputAtoms=False):
+ """
+ Parameters
+ ----------
+ label, outdir: str
+ A label and directory to read/write data from/to
+ sftfilepattern: str
+ Pattern to match SFTs using wildcards (*?) and ranges [0-9];
+ mutiple patterns can be given separated by colons.
+ F0s, F1s, F2s, delta_F0s, delta_F1s, tglitchs, Alphas, Deltas: tuple
+ Length 3 tuple describing the grid for each parameter, e.g
+ [F0min, F0max, dF0], for a fixed value simply give [F0]. Unless
+ input_arrays == True, then these are the values to search at.
+ tref, minStartTime, maxStartTime: int
+ GPS seconds of the reference time, start time and end time
+ input_arrays: bool
+ if true, use the F0s, F1s, etc as is
+ transientWindowType: str
+ If 'rect' or 'exp', compute atoms so that a transient (t0,tau) map
+ can later be computed. ('none' instead of None explicitly calls
+ the transient-window function, but with the full range, for
+ debugging). Currently only supported for nsegs=1.
+ t0Band, tauBand: int
+ if >0, search t0 in (minStartTime,minStartTime+t0Band)
+ and tau in (2*Tsft,2*Tsft+tauBand).
+ if =0, only compute CW Fstat with t0=minStartTime,
+ tau=maxStartTime-minStartTime.
+ dt0, dtau: int
+ grid resolutions in transient start-time and duration,
+ both default to Tsft
+ 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
+ """
+
+ self.nsegs = 1
+ if os.path.isdir(outdir) is False:
+ os.mkdir(outdir)
+ self.set_out_file()
+ self.keys = ['_', '_', 'F0', 'F1', 'F2', 'Alpha', 'Delta']
+ self.search_keys = [x+'s' for x in self.keys[2:]]
+ for k in self.search_keys:
+ setattr(self, k, np.atleast_1d(getattr(self, k)))
+
+ def inititate_search_object(self):
+ logging.info('Setting up search object')
+ self.search = ComputeFstat(
+ tref=self.tref, sftfilepattern=self.sftfilepattern,
+ minCoverFreq=self.minCoverFreq, maxCoverFreq=self.maxCoverFreq,
+ detectors=self.detectors,
+ transientWindowType=self.transientWindowType,
+ t0Band=self.t0Band, tauBand=self.tauBand,
+ dt0=self.dt0, dtau=self.dtau,
+ minStartTime=self.minStartTime, maxStartTime=self.maxStartTime,
+ BSGL=self.BSGL, SSBprec=self.SSBprec,
+ injectSources=self.injectSources,
+ assumeSqrtSX=self.assumeSqrtSX)
+ self.search.get_det_stat = self.search.get_fullycoherent_twoF
+
+ def run(self, return_data=False):
+ self.get_input_data_array()
+ old_data = self.check_old_data_is_okay_to_use()
+ if old_data is not False:
+ self.data = old_data
+ return
+
+ if hasattr(self, 'search') is False:
+ self.inititate_search_object()
+
+ data = []
+ for vals in tqdm(self.input_data):
+ detstat = self.search.get_det_stat(*vals)
+ windowRange = getattr(self.search, 'windowRange', None)
+ FstatMap = getattr(self.search, 'FstatMap', None)
+ thisCand = list(vals) + [detstat]
+ if getattr(self, 'transientWindowType', None):
+ if self.outputTransientFstatMap:
+ 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
+ Fmn = FstatMap.F_mn.data
+ maxidx = np.unravel_index(Fmn.argmax(), Fmn.shape)
+ thisCand += [windowRange.t0+maxidx[0]*windowRange.dt0,
+ windowRange.tau+maxidx[1]*windowRange.dtau]
+ data.append(thisCand)
+ if self.outputAtoms:
+ self.search.write_atoms_to_file(os.path.splitext(self.out_file)[0])
+
+ data = np.array(data, dtype=np.float)
+ if return_data:
+ return data
+ else:
+ self.save_array_to_disk(data)
+ self.data = data
+
+
class SliceGridSearch(GridSearch):
""" Slice gridded search using ComputeFstat """
@helper_functions.initializer
diff --git a/pyfstat/helper_functions.py b/pyfstat/helper_functions.py
index 5c011fd91fda3a2785ee596d7c5b85a14e00baf0..556503e347080d1ab785de5bd8048be537ed5505 100644
--- a/pyfstat/helper_functions.py
+++ b/pyfstat/helper_functions.py
@@ -11,12 +11,19 @@ import inspect
import peakutils
from functools import wraps
from scipy.stats.distributions import ncx2
+import numpy as np
import lal
import lalpulsar
-import matplotlib.pyplot as plt
-import numpy as np
-
+# workaround for matplotlib on X-less remote logins
+if 'DISPLAY' in os.environ:
+ import matplotlib.pyplot as plt
+else:
+ logging.info('No $DISPLAY environment variable found, so importing \
+ matplotlib.pyplot with non-interactive "Agg" backend.')
+ import matplotlib
+ matplotlib.use('Agg')
+ import matplotlib.pyplot as plt
def set_up_optional_tqdm():
try:
diff --git a/tests.py b/tests.py
index 4b853e4fd0d61247fe68a75c481ffc935a542ba8..bf92c9ffb3aa9b417d2aefce64111128f83f329b 100644
--- a/tests.py
+++ b/tests.py
@@ -4,6 +4,7 @@ import os
import shutil
import pyfstat
import lalpulsar
+import logging
class Test(unittest.TestCase):
@@ -12,7 +13,11 @@ class Test(unittest.TestCase):
@classmethod
def setUpClass(self):
if os.path.isdir(self.outdir):
- shutil.rmtree(self.outdir)
+ try:
+ shutil.rmtree(self.outdir)
+ except OSError:
+ logging.warning(
+ "{} not removed prior to tests".format(self.outdir))
h0 = 1
sqrtSX = 1
F0 = 30
@@ -38,7 +43,11 @@ class Test(unittest.TestCase):
@classmethod
def tearDownClass(self):
if os.path.isdir(self.outdir):
- shutil.rmtree(self.outdir)
+ try:
+ shutil.rmtree(self.outdir)
+ except OSError:
+ logging.warning(
+ "{} not removed prior to tests".format(self.outdir))
class Writer(Test):