Improvements to the template estimation (V, Vsky, Vpe)

1) Moves the actual estimation to external function since it is not
contingent on using the MCMC methods etc
2) Change notation to V, Vsky and Vpe
3) Fix some table outputs (Tcoh was forcing integer)
4) Rename of names -> detector_names
5) Adds tests
6) Writes a simpler table if the sky is not searched over

pyfstat.py

@@ -125,6 +125,71 @@ def read_par(label, outdir):
     return d
 
 
+def get_V_estimate(
+        nsegs, tref, minStartTime, maxStartTime, DeltaOmega, DeltaFs,
+        fiducial_freq, detector_names, earth_ephem, sun_ephem):
+    """ Returns V, Vsky, Vpe estimated from the super-sky metric
+
+    Parameters
+    ----------
+    nsegs: int
+        Number of semi-coherent segments
+    tref: int
+        Reference time in GPS seconds
+    minStartTime, maxStartTime: int
+        Minimum and maximum SFT timestamps
+    DeltaOmega: float
+        Solid angle of the sky-patch
+    DeltaFs: array
+        Array of [DeltaF0, DeltaF1, ...], length determines the number of
+        spin-down terms.
+    fiducial_freq: float
+        Fidicual frequency
+    detector_names: array
+        Array of detectors to average over
+    earth_ephem, sun_ephem: st
+        Paths to the ephemeris files
+
+    """
+    spindowns = len(DeltaFs) - 1
+    tboundaries = np.linspace(minStartTime, maxStartTime, nsegs+1)
+
+    ref_time = lal.LIGOTimeGPS(tref)
+    segments = lal.SegListCreate()
+    for j in range(len(tboundaries)-1):
+        seg = lal.SegCreate(lal.LIGOTimeGPS(tboundaries[j]),
+                            lal.LIGOTimeGPS(tboundaries[j+1]),
+                            j)
+        lal.SegListAppend(segments, seg)
+    detNames = lal.CreateStringVector(*detector_names)
+    detectors = lalpulsar.MultiLALDetector()
+    lalpulsar.ParseMultiLALDetector(detectors, detNames)
+    detector_weights = None
+    detector_motion = (lalpulsar.DETMOTION_SPIN
+                       + lalpulsar.DETMOTION_ORBIT)
+    ephemeris = lalpulsar.InitBarycenter(earth_ephem, sun_ephem)
+    try:
+        SSkyMetric = lalpulsar.ComputeSuperskyMetrics(
+            spindowns, ref_time, segments, fiducial_freq, detectors,
+            detector_weights, detector_motion, ephemeris)
+    except RuntimeError as e:
+        logging.debug('Encountered run-time error {}'.format(e))
+        return None, None, None
+
+    sqrtdetG_SKY = np.sqrt(np.linalg.det(
+        SSkyMetric.semi_rssky_metric.data[:2, :2]))
+    sqrtdetG_PE = np.sqrt(np.linalg.det(
+        SSkyMetric.semi_rssky_metric.data[2:, 2:]))
+
+    Vsky = .5*sqrtdetG_SKY*DeltaOmega
+    Vpe = sqrtdetG_PE * np.prod(DeltaFs)
+    if Vsky == 0:
+        Vsky = 1
+    if Vpe == 0:
+        Vpe = 1
+    return (Vsky * Vpe, Vsky, Vpe)
+
+
 class BaseSearchClass(object):
     """ The base search class, provides general functions """
 
@@ -291,8 +356,8 @@ class ComputeFstat(object):
         logging.info('Loading data matching pattern {}'.format(
                      self.sftfilepath))
         SFTCatalog = lalpulsar.SFTdataFind(self.sftfilepath, constraints)
-        names = list(set([d.header.name for d in SFTCatalog.data]))
-        self.names = names
+        detector_names = list(set([d.header.name for d in SFTCatalog.data]))
+        self.detector_names = detector_names
         SFT_timestamps = [d.header.epoch for d in SFTCatalog.data]
         try:
             from bashplotlib.histogram import plot_hist
@@ -300,10 +365,10 @@ class ComputeFstat(object):
             plot_hist(SFT_timestamps, height=5, bincount=50)
         except IOError:
             pass
-        if len(names) == 0:
+        if len(detector_names) == 0:
             raise ValueError('No data loaded.')
         logging.info('Loaded {} data files from detectors {}'.format(
-            len(SFT_timestamps), names))
+            len(SFT_timestamps), detector_names))
         logging.info('Data spans from {} ({}) to {} ({})'.format(
             int(SFT_timestamps[0]),
             subprocess.check_output('lalapps_tconvert {}'.format(
@@ -2009,74 +2074,39 @@ class MCMCFollowUpSearch(MCMCSemiCoherentSearch):
         if prior[key]['type'] == 'unif':
             return .5*(prior[key]['upper'] + prior[key]['lower'])
 
-    def get_number_of_templates_estimate(self, nsegs):
-        """ Returns V, Vsky, Vf estimated from the super-sky metric """
-        if args.no_template_counting:
-            return 'N/A'
-        tboundaries = np.linspace(self.minStartTime, self.maxStartTime,
-                                  nsegs+1)
-
+    def init_V_estimate_parameters(self):
         if 'Alpha' in self.theta_keys:
             DeltaAlpha = self.get_width_from_prior(self.theta_prior, 'Alpha')
             DeltaDelta = self.get_width_from_prior(self.theta_prior, 'Delta')
             DeltaMid = self.get_mid_from_prior(self.theta_prior, 'Delta')
             DeltaOmega = np.sin(DeltaMid) * DeltaDelta * DeltaAlpha
+        else:
+            DeltaOmega = 0
         if 'F0' in self.theta_keys:
             DeltaF0 = self.get_width_from_prior(self.theta_prior, 'F0')
         else:
-            DeltaF0 = 1
+            raise ValueError('You are searching over F0?')
+        DeltaFs = [DeltaF0]
         if 'F1' in self.theta_keys:
             DeltaF1 = self.get_width_from_prior(self.theta_prior, 'F1')
-            spindowns = 1
-        else:
-            DeltaF1 = 1
-            spindowns = 0
+            DeltaFs.append(DeltaF1)
+        if 'F2' in self.theta_keys:
+            DeltaF2 = self.get_width_from_prior(self.theta_prior, 'F2')
+            DeltaFs.append(DeltaF2)
 
-        ref_time = lal.LIGOTimeGPS(self.tref)
-        segments = lal.SegListCreate()
-        for j in range(len(tboundaries)-1):
-            seg = lal.SegCreate(lal.LIGOTimeGPS(tboundaries[j]),
-                                lal.LIGOTimeGPS(tboundaries[j+1]),
-                                j)
-            lal.SegListAppend(segments, seg)
         if type(self.theta_prior['F0']) == dict:
             fiducial_freq = self.get_mid_from_prior(self.theta_prior, 'F0')
         else:
             fiducial_freq = self.theta_prior['F0']
-        detector_names = self.search.names
-        detNames = lal.CreateStringVector(*detector_names)
-        detectors = lalpulsar.MultiLALDetector()
-        lalpulsar.ParseMultiLALDetector(detectors, detNames)
-        detector_weights = None
-        detector_motion = (lalpulsar.DETMOTION_SPIN
-                           + lalpulsar.DETMOTION_ORBIT)
-        ephemeris = lalpulsar.InitBarycenter(self.earth_ephem,
-                                             self.sun_ephem)
-        try:
-            SSkyMetric = lalpulsar.ComputeSuperskyMetrics(
-                spindowns, ref_time, segments, fiducial_freq, detectors,
-                detector_weights, detector_motion, ephemeris)
-        except RuntimeError as e:
-            print e
-            return 'N/A'
 
-        sqrtdetG_SKY = np.sqrt(np.linalg.det(
-            SSkyMetric.semi_rssky_metric.data[:2, :2]))
-        sqrtdetG_F = np.sqrt(np.linalg.det(
-            SSkyMetric.semi_rssky_metric.data[2:, 2:]))
-
-        if 'Alpha' in self.theta_keys:
-            return (.5*sqrtdetG_SKY*sqrtdetG_F*DeltaOmega*DeltaF1*DeltaF0,
-                    .5*sqrtdetG_SKY*DeltaOmega,
-                    sqrtdetG_F*DeltaF1*DeltaF0)
-        else:
-            return (sqrtdetG_F*DeltaF1*DeltaF0,
-                    'N/A',
-                    sqrtdetG_F*DeltaF1*DeltaF0)
+        return fiducial_freq, DeltaOmega, DeltaFs
 
     def init_run_setup(self, run_setup, log_table=True, gen_tex_table=True):
         logging.info('Calculating the number of templates for this setup..')
-        number_of_templates = []
+        Vs = []
+        Vskys = []
+        Vpes = []
+        fiducial_freq, DeltaOmega, DeltaFs = self.init_V_estimate_parameters()
         for i, rs in enumerate(run_setup):
             rs = list(rs)
             if len(rs) == 2:
@@ -2084,49 +2114,72 @@ class MCMCFollowUpSearch(MCMCSemiCoherentSearch):
             if np.shape(rs[0]) == ():
                 rs[0] = (rs[0], 0)
             run_setup[i] = rs
-            number_of_templates.append(
-                    self.get_number_of_templates_estimate(rs[1]))
+
+            V, Vsky, Vpe = get_V_estimate(
+                rs[1], self.tref, self.minStartTime, self.maxStartTime,
+                DeltaOmega, DeltaFs, fiducial_freq, self.search.detector_names,
+                self.earth_ephem, self.sun_ephem)
+            Vs.append(V)
+            Vskys.append(Vsky)
+            Vpes.append(Vpe)
 
         if log_table:
             logging.info('Using run-setup as follow:')
-            logging.info('Stage | nburn | nprod | nsegs | Tcoh | resetp0 |'
+            logging.info('Stage | nburn | nprod | nsegs | Tcoh d | resetp0 |'
                          '# templates = # sky x # Freq')
             for i, rs in enumerate(run_setup):
                 Tcoh = (self.maxStartTime - self.minStartTime) / rs[1] / 86400
-                if number_of_templates[i] == 'N/A':
-                    vtext = number_of_templates[i]
-                elif 'N/A' in number_of_templates[i]:
-                    vtext = '{:1.0e} = {} x {:1.0e}'.format(
-                            *number_of_templates[i])
+                if Vs[i] is None:
+                    vtext = 'N/A'
                 else:
                     vtext = '{:1.0e} = {:1.0e} x {:1.0e}'.format(
-                            *number_of_templates[i])
+                            Vs[i], Vskys[i], Vpes[i])
                 logging.info('{} | {} | {} | {} | {} | {} | {}'.format(
                     str(i).ljust(5), str(rs[0][0]).ljust(5),
                     str(rs[0][1]).ljust(5), str(rs[1]).ljust(5),
-                    '{:1.2f}'.format(Tcoh).ljust(4), str(rs[2]).ljust(7), vtext))
+                    '{:6.1f}'.format(Tcoh), str(rs[2]).ljust(7),
+                    vtext))
 
         if gen_tex_table:
             filename = '{}/{}_run_setup.tex'.format(self.outdir, self.label)
+            if DeltaOmega > 0:
                 with open(filename, 'w+') as f:
                     f.write(r'\begin{tabular}{c|cccccc}' + '\n')
                     f.write(r'Stage & $\Nseg$ & $\Tcoh^{\rm days}$ &'
-                        r'$\Nsteps$ & $\V$ & $\Vsky$ & $\Vf$ \\ \hline'
+                            r'$\Nsteps$ & $\V$ & $\Vsky$ & $\Vpe$ \\ \hline'
                             '\n')
                     for i, rs in enumerate(run_setup):
-                    Tcoh = (self.maxStartTime - self.minStartTime) / rs[1] / 86400
+                        Tcoh = float(self.maxStartTime - self.minStartTime)/rs[1]/86400
                         line = r'{} & {} & {} & {} & {} & {} & {} \\' + '\n'
-                    if number_of_templates[i] == 'N/A':
-                        V = Vsky = Vf = 'N/A'
+                        if Vs[i] is None:
+                            Vs[i] = Vskys[i] = Vpes[i] = 'N/A'
+                        if rs[0][-1] == 0:
+                            nsteps = rs[0][0]
                         else:
-                        V, Vsky, Vf = number_of_templates[i]
+                            nsteps = '{},{}'.format(*rs[0])
+                        line = line.format(i, rs[1], Tcoh, nsteps,
+                                           texify_float(Vs[i]),
+                                           texify_float(Vskys[i]),
+                                           texify_float(Vpes[i]))
+                        f.write(line)
+                    f.write(r'\end{tabular}' + '\n')
+            else:
+                with open(filename, 'w+') as f:
+                    f.write(r'\begin{tabular}{c|cccc}' + '\n')
+                    f.write(r'Stage & $\Nseg$ & $\Tcoh^{\rm days}$ &'
+                            r'$\Nsteps$ & $\Vpe$ \\ \hline'
+                            '\n')
+                    for i, rs in enumerate(run_setup):
+                        Tcoh = float(self.maxStartTime - self.minStartTime)/rs[1]/86400
+                        line = r'{} & {} & {} & {} & {} \\' + '\n'
+                        if Vs[i] is None:
+                            Vs[i] = Vskys[i] = Vpes[i] = 'N/A'
                         if rs[0][-1] == 0:
                             nsteps = rs[0][0]
                         else:
                             nsteps = '{},{}'.format(*rs[0])
                         line = line.format(i, rs[1], Tcoh, nsteps,
-                                       texify_float(V), texify_float(Vsky),
-                                       texify_float(Vf))
+                                           texify_float(Vpes[i]))
                         f.write(line)
                     f.write(r'\end{tabular}' + '\n')
 

tests.py

@@ -194,6 +194,42 @@ class TestMCMCSearch(Test):
             FS > predicted_FS or np.abs((FS-predicted_FS))/predicted_FS < 0.3)
 
 
+class TestAuxillaryFunctions(Test):
+    nsegs = 10
+    minStartTime = 1e9
+    maxStartTime = minStartTime + 100 * 86400
+    tref = .5*(minStartTime + maxStartTime)
+    DeltaOmega = 1e-2
+    DeltaFs = [1e-4, 1e-14]
+    fiducial_freq = 100
+    detector_names = ['H1', 'L1']
+    earth_ephem = pyfstat.earth_ephem
+    sun_ephem = pyfstat.sun_ephem
+
+    def test_get_V_estimate_sky_F0_F1(self):
+
+        out = pyfstat.get_V_estimate(
+            self.nsegs, self.tref, self.minStartTime, self.maxStartTime,
+            self.DeltaOmega, self.DeltaFs, self.fiducial_freq,
+            self.detector_names, self.earth_ephem, self.sun_ephem)
+        V, Vsky, Vpe = out
+        self.assertTrue(V == Vsky * Vpe)
+        self.__class__.Vpe_COMPUTED_WITH_SKY = Vpe
+
+    def test_get_V_estimate_F0_F1(self):
+        out = pyfstat.get_V_estimate(
+            self.nsegs, self.tref, self.minStartTime, self.maxStartTime,
+            self.DeltaOmega, self.DeltaFs, self.fiducial_freq,
+            self.detector_names, self.earth_ephem, self.sun_ephem)
+        V, Vsky, Vpe = out
+        self.assertTrue(V == Vsky * Vpe)
+        self.__class__.Vpe_COMPUTED_WITHOUT_SKY = Vpe
+
+    def test_the_equivalence_of_Vpe(self):
+        """Tests if the Vpe computed with and without the sky are equal """
+        self.assertEqual(self.__class__.Vpe_COMPUTED_WITHOUT_SKY,
+                         self.__class__.Vpe_COMPUTED_WITH_SKY)
+
 if __name__ == '__main__':
     outdir = 'TestData'
     if os.path.isdir(outdir):