Merge branch 'develop-DK' into 'master'

add transient F(t0,tau) file output

See merge request GregAshton/PyFstat!8

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()
 

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)
 

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
         """
@@ -124,7 +128,7 @@ class GridSearch(BaseSearchClass):
         if args.clean:
             return False
         if os.path.isfile(self.out_file) is False:
-            logging.info('No old data found, continuing with grid search')
+            logging.info('No old data found in "{:s}", continuing with grid search'.format(self.out_file))
             return False
         if self.sftfilepattern is not None:
             oldest_sft = min([os.path.getmtime(f) for f in
@@ -135,13 +139,13 @@ class GridSearch(BaseSearchClass):
                 return False
 
         data = np.atleast_2d(np.genfromtxt(self.out_file, delimiter=' '))
-        if np.all(data[:, 0:-1] == self.input_data):
+        if np.all(data[:,0:len(self.coord_arrays)] == self.input_data[:,0:len(self.coord_arrays)]):
             logging.info(
-                'Old data found with matching input, no search performed')
+                'Old data found in "{:s}" with matching input, no search performed'.format(self.out_file))
             return data
         else:
             logging.info(
-                'Old data found, input differs, continuing with grid search')
+                'Old data found in "{:s}", input differs, continuing with grid search'.format(self.out_file))
             return False
         return False
 
@@ -157,8 +161,21 @@ 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)
+            thisCand = list(vals) + [detstat]
+            if self.transientWindowType:
+                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)
         if return_data: