Merge branch 'master' into develop-GA

.gitlab-ci.yml

 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

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,

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

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)
+            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.tau = int(self.maxStartTime-self.minStartTime)
+                self.windowRange.tauBand = 0
+            else: # user-set bands and spacings
                 if self.t0Band is None:
                     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.dt0:
+                        self.windowRange.dt0 = self.dt0
                 if self.tauBand is None:
                     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
+                    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 """

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

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
 
+# 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
-import numpy as np
-
 
 def set_up_optional_tqdm():
     try:

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):
+            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):
+            try:
                 shutil.rmtree(self.outdir)
+            except OSError:
+                logging.warning(
+                    "{} not removed prior to tests".format(self.outdir))
 
 
 class Writer(Test):