Skip to content
Snippets Groups Projects
Commit e4a38195 authored by Gregory Ashton's avatar Gregory Ashton
Browse files

Adds BGSL functionality to all searches 

Adds a switch to calculate (or optimize over) BSGL rather than 2F. This
should help to quiten noise events which exist only in one detector.
parent a04d15fd
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Showing
with 62 additions and 14 deletions
pyfstat.py
+ 62
− 14
View file @ e4a38195
...@@ -165,7 +165,7 @@ class ComputeFstat(object): ...@@ -165,7 +165,7 @@ class ComputeFstat(object):
minStartTime=None, maxStartTime=None, minStartTime=None, maxStartTime=None,
minCoverFreq=None, maxCoverFreq=None, minCoverFreq=None, maxCoverFreq=None,
detector=None, earth_ephem=None, sun_ephem=None, detector=None, earth_ephem=None, sun_ephem=None,
binary=False, transient=True): binary=False, transient=True, BSGL=False):
""" """
Parameters Parameters
---------- ----------
...@@ -188,6 +188,8 @@ class ComputeFstat(object): ...@@ -188,6 +188,8 @@ class ComputeFstat(object):
If true, search of binary parameters. If true, search of binary parameters.
transient: bool transient: bool
If true, allow for the Fstat to be computed over a transient range. If true, allow for the Fstat to be computed over a transient range.
BSGL: bool
If true, compute the BSGL rather than the twoF value.
""" """
...@@ -258,7 +260,24 @@ class ComputeFstat(object): ...@@ -258,7 +260,24 @@ class ComputeFstat(object):
logging.info('Initialising FstatResults') logging.info('Initialising FstatResults')
self.FstatResults = lalpulsar.FstatResults() self.FstatResults = lalpulsar.FstatResults()
if self.BSGL:
logging.info('Initialising BSGL: this will fail if numDet < 2')
# Tuning parameters - to be reviewed
numDetectors = 2
Fstar0sc = 40.
oLGX = np.zeros(10)
oLGX[:numDetectors] = 0.5
self.BSGLSetup = lalpulsar.CreateBSGLSetup(numDetectors,
Fstar0sc,
oLGX,
False,
1)
self.twoFX = np.zeros(10)
self.whatToCompute = (lalpulsar.FSTATQ_2F +
lalpulsar.FSTATQ_2F_PER_DET)
if self.transient: if self.transient:
logging.info('Initialising transient parameters')
self.windowRange = lalpulsar.transientWindowRange_t() self.windowRange = lalpulsar.transientWindowRange_t()
self.windowRange.type = lalpulsar.TRANSIENT_RECTANGULAR self.windowRange.type = lalpulsar.TRANSIENT_RECTANGULAR
self.windowRange.t0Band = 0 self.windowRange.t0Band = 0
...@@ -290,14 +309,41 @@ class ComputeFstat(object): ...@@ -290,14 +309,41 @@ class ComputeFstat(object):
) )
if self.transient is False: if self.transient is False:
if self.BSGL is False:
return self.FstatResults.twoF[0] return self.FstatResults.twoF[0]
twoF = np.float(self.FstatResults.twoF[0])
self.twoFX[0] = self.FstatResults.twoFPerDet(0)
self.twoFX[1] = self.FstatResults.twoFPerDet(1)
BSGL = lalpulsar.ComputeBSGL(twoF, self.twoFX,
self.BSGLSetup)
return BSGL
self.windowRange.t0 = int(tstart) # TYPE UINT4 self.windowRange.t0 = int(tstart) # TYPE UINT4
self.windowRange.tau = int(tend - tstart) # TYPE UINT4 self.windowRange.tau = int(tend - tstart) # TYPE UINT4
FS = lalpulsar.ComputeTransientFstatMap( FS = lalpulsar.ComputeTransientFstatMap(
self.FstatResults.multiFatoms[0], self.windowRange, False) self.FstatResults.multiFatoms[0], self.windowRange, False)
if self.BSGL is False:
return 2*FS.F_mn.data[0][0] return 2*FS.F_mn.data[0][0]
FstatResults_single = copy.copy(self.FstatResults)
FstatResults_single.lenth = 1
FstatResults_single.data = self.FstatResults.multiFatoms[0].data[0]
FS0 = lalpulsar.ComputeTransientFstatMap(
FstatResults_single.multiFatoms[0], self.windowRange, False)
FstatResults_single.data = self.FstatResults.multiFatoms[0].data[1]
FS1 = lalpulsar.ComputeTransientFstatMap(
FstatResults_single.multiFatoms[0], self.windowRange, False)
self.twoFX[0] = 2*FS0.F_mn.data[0][0]
self.twoFX[1] = 2*FS1.F_mn.data[0][0]
BSGL = lalpulsar.ComputeBSGL(2*FS.F_mn.data[0][0], self.twoFX,
self.BSGLSetup)
return BSGL
class SemiCoherentGlitchSearch(BaseSearchClass, ComputeFstat): class SemiCoherentGlitchSearch(BaseSearchClass, ComputeFstat):
""" A semi-coherent glitch search """ A semi-coherent glitch search
...@@ -310,9 +356,10 @@ class SemiCoherentGlitchSearch(BaseSearchClass, ComputeFstat): ...@@ -310,9 +356,10 @@ class SemiCoherentGlitchSearch(BaseSearchClass, ComputeFstat):
@initializer @initializer
def __init__(self, label, outdir, tref, tstart, tend, nglitch=0, def __init__(self, label, outdir, tref, tstart, tend, nglitch=0,
sftlabel=None, sftdir=None, theta0_idx=0, minCoverFreq=None, sftlabel=None, sftdir=None, theta0_idx=0, BSGL=False,
maxCoverFreq=None, minStartTime=None, maxStartTime=None, minCoverFreq=None, maxCoverFreq=None, minStartTime=None,
detector=None, earth_ephem=None, sun_ephem=None): maxStartTime=None, detector=None, earth_ephem=None,
sun_ephem=None):
""" """
Parameters Parameters
---------- ----------
...@@ -384,7 +431,7 @@ class SemiCoherentGlitchSearch(BaseSearchClass, ComputeFstat): ...@@ -384,7 +431,7 @@ class SemiCoherentGlitchSearch(BaseSearchClass, ComputeFstat):
if np.isfinite(twoFSum): if np.isfinite(twoFSum):
return twoFSum return twoFSum
else: else:
return 0 return -np.inf
def compute_glitch_fstat_single(self, F0, F1, F2, Alpha, Delta, delta_F0, def compute_glitch_fstat_single(self, F0, F1, F2, Alpha, Delta, delta_F0,
delta_F1, tglitch): delta_F1, tglitch):
...@@ -423,8 +470,8 @@ class MCMCSearch(BaseSearchClass): ...@@ -423,8 +470,8 @@ class MCMCSearch(BaseSearchClass):
def __init__(self, label, outdir, sftlabel, sftdir, theta_prior, tref, def __init__(self, label, outdir, sftlabel, sftdir, theta_prior, tref,
tstart, tend, nsteps=[100, 100, 100], nwalkers=100, ntemps=1, tstart, tend, nsteps=[100, 100, 100], nwalkers=100, ntemps=1,
log10temperature_min=-5, theta_initial=None, scatter_val=1e-4, log10temperature_min=-5, theta_initial=None, scatter_val=1e-4,
binary=False, minCoverFreq=None, maxCoverFreq=None, binary=False, BSGL=False, minCoverFreq=None, maxCoverFreq=None,
detector=None, earth_ephem=None, sun_ephem=None): detector=None, earth_ephem=None, sun_ephem=None, theta0_idx=0):
""" """
Parameters Parameters
label, outdir: str label, outdir: str
...@@ -501,7 +548,7 @@ class MCMCSearch(BaseSearchClass): ...@@ -501,7 +548,7 @@ class MCMCSearch(BaseSearchClass):
sftdir=self.sftdir, minCoverFreq=self.minCoverFreq, sftdir=self.sftdir, minCoverFreq=self.minCoverFreq,
maxCoverFreq=self.maxCoverFreq, earth_ephem=self.earth_ephem, maxCoverFreq=self.maxCoverFreq, earth_ephem=self.earth_ephem,
sun_ephem=self.sun_ephem, detector=self.detector, sun_ephem=self.sun_ephem, detector=self.detector,
transient=False, BSGL=self.BSGL, transient=False,
minStartTime=self.minStartTime, maxStartTime=self.maxStartTime) minStartTime=self.minStartTime, maxStartTime=self.maxStartTime)
def logp(self, theta_vals, theta_prior, theta_keys, search): def logp(self, theta_vals, theta_prior, theta_keys, search):
...@@ -1095,7 +1142,7 @@ class MCMCGlitchSearch(MCMCSearch): ...@@ -1095,7 +1142,7 @@ class MCMCGlitchSearch(MCMCSearch):
tstart, tend, nglitch=1, nsteps=[100, 100, 100], nwalkers=100, tstart, tend, nglitch=1, nsteps=[100, 100, 100], nwalkers=100,
ntemps=1, log10temperature_min=-5, theta_initial=None, ntemps=1, log10temperature_min=-5, theta_initial=None,
scatter_val=1e-4, dtglitchmin=1*86400, theta0_idx=0, scatter_val=1e-4, dtglitchmin=1*86400, theta0_idx=0,
detector=None, detector=None, BSGL=False,
minCoverFreq=None, maxCoverFreq=None, earth_ephem=None, minCoverFreq=None, maxCoverFreq=None, earth_ephem=None,
sun_ephem=None): sun_ephem=None):
""" """
...@@ -1181,7 +1228,7 @@ class MCMCGlitchSearch(MCMCSearch): ...@@ -1181,7 +1228,7 @@ class MCMCGlitchSearch(MCMCSearch):
sftdir=self.sftdir, tref=self.tref, tstart=self.tstart, sftdir=self.sftdir, tref=self.tref, tstart=self.tstart,
tend=self.tend, minCoverFreq=self.minCoverFreq, tend=self.tend, minCoverFreq=self.minCoverFreq,
maxCoverFreq=self.maxCoverFreq, earth_ephem=self.earth_ephem, maxCoverFreq=self.maxCoverFreq, earth_ephem=self.earth_ephem,
sun_ephem=self.sun_ephem, detector=self.detector, sun_ephem=self.sun_ephem, detector=self.detector, BSGL=self.BSGL,
nglitch=self.nglitch, theta0_idx=self.theta0_idx, nglitch=self.nglitch, theta0_idx=self.theta0_idx,
minStartTime=self.minStartTime, maxStartTime=self.maxStartTime) minStartTime=self.minStartTime, maxStartTime=self.maxStartTime)
...@@ -1272,7 +1319,7 @@ class GridSearch(BaseSearchClass): ...@@ -1272,7 +1319,7 @@ class GridSearch(BaseSearchClass):
def __init__(self, label, outdir, sftlabel=None, sftdir=None, F0s=[0], def __init__(self, label, outdir, sftlabel=None, sftdir=None, F0s=[0],
F1s=[0], F2s=[0], Alphas=[0], Deltas=[0], tref=None, F1s=[0], F2s=[0], Alphas=[0], Deltas=[0], tref=None,
tstart=None, tend=None, minCoverFreq=None, maxCoverFreq=None, tstart=None, tend=None, minCoverFreq=None, maxCoverFreq=None,
earth_ephem=None, sun_ephem=None, detector=None): earth_ephem=None, sun_ephem=None, detector=None, BSGL=False):
""" """
Parameters Parameters
label, outdir: str label, outdir: str
...@@ -1311,7 +1358,8 @@ class GridSearch(BaseSearchClass): ...@@ -1311,7 +1358,8 @@ class GridSearch(BaseSearchClass):
sftdir=self.sftdir, minCoverFreq=self.minCoverFreq, sftdir=self.sftdir, minCoverFreq=self.minCoverFreq,
maxCoverFreq=self.maxCoverFreq, earth_ephem=self.earth_ephem, maxCoverFreq=self.maxCoverFreq, earth_ephem=self.earth_ephem,
sun_ephem=self.sun_ephem, detector=self.detector, transient=False, sun_ephem=self.sun_ephem, detector=self.detector, transient=False,
minStartTime=minStartTime, maxStartTime=maxStartTime) minStartTime=minStartTime, maxStartTime=maxStartTime,
BSGL=BSGL)
if os.path.isdir(outdir) is False: if os.path.isdir(outdir) is False:
os.mkdir(outdir) os.mkdir(outdir)
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment

AEI Privacy Policy Imprint