Exposes assumeSqrtSX functionality

pyfstat/core.py

@@ -151,8 +151,8 @@ class ComputeFstat(object):
     def __init__(self, tref, sftfilepath=None, minStartTime=None,
                  maxStartTime=None, binary=False, transient=True, BSGL=False,
                  detector=None, minCoverFreq=None, maxCoverFreq=None,
-                 earth_ephem=None, sun_ephem=None, injectSources=None
-                 ):
+                 earth_ephem=None, sun_ephem=None, injectSources=None,
+                 assumeSqrtSX=None):
         """
         Parameters
         ----------
@@ -180,6 +180,10 @@ class ComputeFstat(object):
             Paths of the two files containing positions of Earth and Sun,
             respectively at evenly spaced times, as passed to CreateFstatInput.
             If None defaults defined in BaseSearchClass will be used.
+        assumeSqrtSX: float
+            Don't estimate noise-floors but assume (stationary) per-IFO
+            sqrt{SX} (if single value: use for all IFOs). If signal only,
+            set sqrtSX=1
 
         """
 
@@ -250,7 +254,14 @@ class ComputeFstat(object):
         FstatOAs.runningMedianWindow = lalpulsar.FstatOptionalArgsDefaults.runningMedianWindow
         FstatOAs.FstatMethod = lalpulsar.FstatOptionalArgsDefaults.FstatMethod
         FstatOAs.InjectSqrtSX = lalpulsar.FstatOptionalArgsDefaults.injectSqrtSX
+        if self.assumeSqrtSX is None:
             FstatOAs.assumeSqrtSX = lalpulsar.FstatOptionalArgsDefaults.assumeSqrtSX
+        else:
+            mnf = lalpulsar.MultiNoiseFloor()
+            assumeSqrtSX = np.atleast_1d(self.assumeSqrtSX)
+            mnf.sqrtSn[:len(assumeSqrtSX)] = assumeSqrtSX
+            mnf.length = len(assumeSqrtSX)
+            FstatOAs.assumeSqrtSX = mnf
         FstatOAs.prevInput = lalpulsar.FstatOptionalArgsDefaults.prevInput
         FstatOAs.collectTiming = lalpulsar.FstatOptionalArgsDefaults.collectTiming
 
@@ -460,7 +471,7 @@ class SemiCoherentSearch(BaseSearchClass, ComputeFstat):
                  binary=False, BSGL=False, minStartTime=None,
                  maxStartTime=None, minCoverFreq=None, maxCoverFreq=None,
                  detector=None, earth_ephem=None, sun_ephem=None,
-                 injectSources=None):
+                 injectSources=None, assumeSqrtSX=None):
         """
         Parameters
         ----------
@@ -570,7 +581,7 @@ class SemiCoherentGlitchSearch(BaseSearchClass, ComputeFstat):
     @helper_functions.initializer
     def __init__(self, label, outdir, tref, minStartTime, maxStartTime,
                  nglitch=0, sftfilepath=None, theta0_idx=0, BSGL=False,
-                 minCoverFreq=None, maxCoverFreq=None,
+                 minCoverFreq=None, maxCoverFreq=None, assumeSqrtSX=None,
                  detector=None, earth_ephem=None, sun_ephem=None):
         """
         Parameters

pyfstat/mcmc_based_searches.py

@@ -13,7 +13,8 @@ import emcee
 import corner
 import dill as pickle
 
-from core import BaseSearchClass, ComputeFstat
+from core import BaseSearchClass, ComputeFstat, SemiCoherentSearch
+from optimal_setup_functions import get_V_estimate
 from core import tqdm, args, earth_ephem, sun_ephem
 from optimal_setup_functions import get_optimal_setup
 import helper_functions
@@ -28,7 +29,7 @@ class MCMCSearch(BaseSearchClass):
                  theta_initial=None, scatter_val=1e-10,
                  binary=False, BSGL=False, minCoverFreq=None,
                  maxCoverFreq=None, detector=None, earth_ephem=None,
-                 sun_ephem=None, injectSources=None):
+                 sun_ephem=None, injectSources=None, assumeSqrtSX=None):
         """
         Parameters
         label, outdir: str
@@ -113,7 +114,8 @@ class MCMCSearch(BaseSearchClass):
             earth_ephem=self.earth_ephem, sun_ephem=self.sun_ephem,
             detector=self.detector, BSGL=self.BSGL, transient=False,
             minStartTime=self.minStartTime, maxStartTime=self.maxStartTime,
-            binary=self.binary, injectSources=self.injectSources)
+            binary=self.binary, injectSources=self.injectSources,
+            assumeSqrtSX=self.assumeSqrtSX)
 
     def logp(self, theta_vals, theta_prior, theta_keys, search):
         H = [self.generic_lnprior(**theta_prior[key])(p) for p, key in
@@ -1247,7 +1249,7 @@ class MCMCSemiCoherentSearch(MCMCSearch):
                  scatter_val=1e-10, detector=None, BSGL=False,
                  minStartTime=None, maxStartTime=None, minCoverFreq=None,
                  maxCoverFreq=None, earth_ephem=None, sun_ephem=None,
-                 injectSources=None):
+                 injectSources=None, assumeSqrtSX=None):
         """
 
         """
@@ -1284,7 +1286,7 @@ class MCMCSemiCoherentSearch(MCMCSearch):
             maxStartTime=self.maxStartTime, minCoverFreq=self.minCoverFreq,
             maxCoverFreq=self.maxCoverFreq, detector=self.detector,
             earth_ephem=self.earth_ephem, sun_ephem=self.sun_ephem,
-            injectSources=self.injectSources)
+            injectSources=self.injectSources, assumeSqrtSX=self.assumeSqrtSX)
 
     def logp(self, theta_vals, theta_prior, theta_keys, search):
         H = [self.generic_lnprior(**theta_prior[key])(p) for p, key in

tests.py

@@ -105,6 +105,23 @@ class TestComputeFstat(Test):
         print predicted_FS, FS
         self.assertTrue(np.abs(predicted_FS-FS)/FS < 0.2)
 
+    def run_computefstatistic_single_point_no_noise(self):
+        Writer = pyfstat.Writer(self.label, outdir=outdir, add_noise=False)
+        Writer.make_data()
+        predicted_FS = Writer.predict_fstat()
+
+        search = pyfstat.ComputeFstat(
+            tref=Writer.tref, assumeSqrtSX=1,
+            sftfilepath='{}/*{}*sft'.format(Writer.outdir, Writer.label))
+        FS = search.run_computefstatistic_single_point(Writer.tstart,
+                                                       Writer.tend,
+                                                       Writer.F0,
+                                                       Writer.F1,
+                                                       Writer.F2,
+                                                       Writer.Alpha,
+                                                       Writer.Delta)
+        print predicted_FS, FS
+        self.assertTrue(np.abs(predicted_FS-FS)/FS < 0.2)
 
 class TestSemiCoherentGlitchSearch(Test):
     label = "Test"