Commit 73fee787 authored by David Keitel's avatar David Keitel
Browse files

configurable transient windows:

 -replaces boolean 'transient' flag
 -window type [none,rect,exp]
  for both injections and searches
 -t0, tau bands for gridsearch,
  single-point use (e.g. MCMC) should be unchanged
parent 5216d198
...@@ -330,7 +330,8 @@ class ComputeFstat(BaseSearchClass): ...@@ -330,7 +330,8 @@ class ComputeFstat(BaseSearchClass):
@helper_functions.initializer @helper_functions.initializer
def __init__(self, tref, sftfilepattern=None, minStartTime=None, def __init__(self, tref, sftfilepattern=None, minStartTime=None,
maxStartTime=None, binary=False, transient=True, BSGL=False, maxStartTime=None, binary=False, BSGL=False,
transientWindowType=None, t0Band=None, tauBand=None,
detectors=None, minCoverFreq=None, maxCoverFreq=None, detectors=None, minCoverFreq=None, maxCoverFreq=None,
injectSources=None, injectSqrtSX=None, assumeSqrtSX=None, injectSources=None, injectSqrtSX=None, assumeSqrtSX=None,
SSBprec=None): SSBprec=None):
...@@ -347,8 +348,11 @@ class ComputeFstat(BaseSearchClass): ...@@ -347,8 +348,11 @@ class ComputeFstat(BaseSearchClass):
this epoch this epoch
binary : bool binary : bool
If true, search of binary parameters. If true, search of binary parameters.
transient : bool transientWindowType: str
If true, allow for the Fstat to be computed over a transient range. If 'rect' or 'exp',
allow for the Fstat to be computed over a transient range.
('none' instead of None explicitly calls the transient-window function,
but with the full range, for debugging)
BSGL : bool BSGL : bool
If true, compute the BSGL rather than the twoF value. If true, compute the BSGL rather than the twoF value.
detectors : str detectors : str
...@@ -477,7 +481,7 @@ class ComputeFstat(BaseSearchClass): ...@@ -477,7 +481,7 @@ class ComputeFstat(BaseSearchClass):
logging.info('Initialising FstatInput') logging.info('Initialising FstatInput')
dFreq = 0 dFreq = 0
if self.transient: if self.transientWindowType:
self.whatToCompute = lalpulsar.FSTATQ_ATOMS_PER_DET self.whatToCompute = lalpulsar.FSTATQ_ATOMS_PER_DET
else: else:
self.whatToCompute = lalpulsar.FSTATQ_2F self.whatToCompute = lalpulsar.FSTATQ_2F
...@@ -593,14 +597,36 @@ class ComputeFstat(BaseSearchClass): ...@@ -593,14 +597,36 @@ class ComputeFstat(BaseSearchClass):
self.whatToCompute = (self.whatToCompute + self.whatToCompute = (self.whatToCompute +
lalpulsar.FSTATQ_2F_PER_DET) lalpulsar.FSTATQ_2F_PER_DET)
if self.transient: if self.transientWindowType:
logging.info('Initialising transient parameters') logging.info('Initialising transient parameters')
self.windowRange = lalpulsar.transientWindowRange_t() self.windowRange = lalpulsar.transientWindowRange_t()
self.windowRange.type = lalpulsar.TRANSIENT_RECTANGULAR transientWindowTypes = {'none': lalpulsar.TRANSIENT_NONE,
self.windowRange.t0Band = 0 'rect': lalpulsar.TRANSIENT_RECTANGULAR,
self.windowRange.dt0 = 1 'exp': lalpulsar.TRANSIENT_EXPONENTIAL}
self.windowRange.tauBand = 0 if self.transientWindowType in transientWindowTypes:
self.windowRange.dtau = 1 self.windowRange.type = transientWindowTypes[self.transientWindowType]
else:
raise ValueError('Unknown window-type ({}) passed as input. Allowed are: [{}].'.format(self.transientWindowType, ', '.join(transientWindowTypes)))
self.Tsft = int(1.0/SFTCatalog.data[0].header.deltaF)
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.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
def get_fullycoherent_twoF(self, tstart, tend, F0, F1, F2, Alpha, Delta, def get_fullycoherent_twoF(self, tstart, tend, F0, F1, F2, Alpha, Delta,
asini=None, period=None, ecc=None, tp=None, asini=None, period=None, ecc=None, tp=None,
...@@ -624,7 +650,7 @@ class ComputeFstat(BaseSearchClass): ...@@ -624,7 +650,7 @@ class ComputeFstat(BaseSearchClass):
self.whatToCompute self.whatToCompute
) )
if self.transient is False: if not self.transientWindowType:
if self.BSGL is False: if self.BSGL is False:
return self.FstatResults.twoF[0] return self.FstatResults.twoF[0]
...@@ -636,7 +662,13 @@ class ComputeFstat(BaseSearchClass): ...@@ -636,7 +662,13 @@ class ComputeFstat(BaseSearchClass):
return log10_BSGL/np.log10(np.exp(1)) return log10_BSGL/np.log10(np.exp(1))
self.windowRange.t0 = int(tstart) # TYPE UINT4 self.windowRange.t0 = int(tstart) # TYPE UINT4
self.windowRange.tau = int(tend - tstart) # TYPE UINT4 if self.windowRange.tauBand == 0:
# true single-template search also in transient params:
# actual (t0,tau) window was set with tstart, tend before
self.windowRange.tau = int(tend - tstart) # TYPE UINT4
else:
# grid search: start at minimum tau required for nondegenerate F-stat computation
self.windowRange.tau = int(2*self.Tsft)
FS = lalpulsar.ComputeTransientFstatMap( FS = lalpulsar.ComputeTransientFstatMap(
self.FstatResults.multiFatoms[0], self.windowRange, False) self.FstatResults.multiFatoms[0], self.windowRange, False)
...@@ -696,8 +728,9 @@ class ComputeFstat(BaseSearchClass): ...@@ -696,8 +728,9 @@ class ComputeFstat(BaseSearchClass):
max_tau = SFTmaxStartTime - tstart max_tau = SFTmaxStartTime - tstart
taus = np.linspace(min_tau, max_tau, npoints) taus = np.linspace(min_tau, max_tau, npoints)
twoFs = [] twoFs = []
if self.transient is False: if not self.transientWindowType:
self.transient = True # still call the transient-Fstat-map function, but using the full range
self.transientWindowType = 'none'
self.init_computefstatistic_single_point() self.init_computefstatistic_single_point()
for tau in taus: for tau in taus:
twoFs.append(self.get_fullycoherent_twoF( twoFs.append(self.get_fullycoherent_twoF(
...@@ -868,7 +901,9 @@ class SemiCoherentSearch(ComputeFstat): ...@@ -868,7 +901,9 @@ class SemiCoherentSearch(ComputeFstat):
self.fs_file_name = "{}/{}_FS.dat".format(self.outdir, self.label) self.fs_file_name = "{}/{}_FS.dat".format(self.outdir, self.label)
self.set_ephemeris_files() self.set_ephemeris_files()
self.transient = True self.transientWindowType = 'rect'
self.t0Band = None
self.tauBand = None
self.init_computefstatistic_single_point() self.init_computefstatistic_single_point()
self.init_semicoherent_parameters() self.init_semicoherent_parameters()
...@@ -876,7 +911,7 @@ class SemiCoherentSearch(ComputeFstat): ...@@ -876,7 +911,7 @@ class SemiCoherentSearch(ComputeFstat):
logging.info(('Initialising semicoherent parameters from {} to {} in' logging.info(('Initialising semicoherent parameters from {} to {} in'
' {} segments').format( ' {} segments').format(
self.minStartTime, self.maxStartTime, self.nsegs)) self.minStartTime, self.maxStartTime, self.nsegs))
self.transient = True self.transientWindowType = 'rect'
self.whatToCompute = lalpulsar.FSTATQ_2F+lalpulsar.FSTATQ_ATOMS_PER_DET self.whatToCompute = lalpulsar.FSTATQ_2F+lalpulsar.FSTATQ_ATOMS_PER_DET
self.tboundaries = np.linspace(self.minStartTime, self.maxStartTime, self.tboundaries = np.linspace(self.minStartTime, self.maxStartTime,
self.nsegs+1) self.nsegs+1)
...@@ -915,7 +950,7 @@ class SemiCoherentSearch(ComputeFstat): ...@@ -915,7 +950,7 @@ class SemiCoherentSearch(ComputeFstat):
self.whatToCompute self.whatToCompute
) )
#if self.transient is False: #if not self.transientWindowType:
# if self.BSGL is False: # if self.BSGL is False:
# return self.FstatResults.twoF[0] # return self.FstatResults.twoF[0]
# twoF = np.float(self.FstatResults.twoF[0]) # twoF = np.float(self.FstatResults.twoF[0])
...@@ -1005,8 +1040,10 @@ class SemiCoherentGlitchSearch(ComputeFstat): ...@@ -1005,8 +1040,10 @@ class SemiCoherentGlitchSearch(ComputeFstat):
self.fs_file_name = "{}/{}_FS.dat".format(self.outdir, self.label) self.fs_file_name = "{}/{}_FS.dat".format(self.outdir, self.label)
self.set_ephemeris_files() self.set_ephemeris_files()
self.transient = True self.transientWindowType = 'rect'
self.binary = False self.t0Band = None
self.tauBand = None
self.binary = False
self.init_computefstatistic_single_point() self.init_computefstatistic_single_point()
def get_semicoherent_nglitch_twoF(self, F0, F1, F2, Alpha, Delta, *args): def get_semicoherent_nglitch_twoF(self, F0, F1, F2, Alpha, Delta, *args):
......
...@@ -31,7 +31,8 @@ class GridSearch(BaseSearchClass): ...@@ -31,7 +31,8 @@ class GridSearch(BaseSearchClass):
Deltas, tref=None, minStartTime=None, maxStartTime=None, Deltas, tref=None, minStartTime=None, maxStartTime=None,
nsegs=1, BSGL=False, minCoverFreq=None, maxCoverFreq=None, nsegs=1, BSGL=False, minCoverFreq=None, maxCoverFreq=None,
detectors=None, SSBprec=None, injectSources=None, detectors=None, SSBprec=None, injectSources=None,
input_arrays=False, assumeSqrtSX=None): input_arrays=False, assumeSqrtSX=None,
transientWindowType=None, t0Band=None, tauBand=None):
""" """
Parameters Parameters
---------- ----------
...@@ -48,6 +49,17 @@ class GridSearch(BaseSearchClass): ...@@ -48,6 +49,17 @@ class GridSearch(BaseSearchClass):
GPS seconds of the reference time, start time and end time GPS seconds of the reference time, start time and end time
input_arrays: bool input_arrays: bool
if true, use the F0s, F1s, etc as is 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.
For all other parameters, see `pyfstat.ComputeFStat` for details For all other parameters, see `pyfstat.ComputeFStat` for details
""" """
...@@ -66,7 +78,9 @@ class GridSearch(BaseSearchClass): ...@@ -66,7 +78,9 @@ class GridSearch(BaseSearchClass):
self.search = ComputeFstat( self.search = ComputeFstat(
tref=self.tref, sftfilepattern=self.sftfilepattern, tref=self.tref, sftfilepattern=self.sftfilepattern,
minCoverFreq=self.minCoverFreq, maxCoverFreq=self.maxCoverFreq, minCoverFreq=self.minCoverFreq, maxCoverFreq=self.maxCoverFreq,
detectors=self.detectors, transient=False, detectors=self.detectors,
transientWindowType=self.transientWindowType,
t0Band=self.t0Band, tauBand=self.tauBand,
minStartTime=self.minStartTime, maxStartTime=self.maxStartTime, minStartTime=self.minStartTime, maxStartTime=self.maxStartTime,
BSGL=self.BSGL, SSBprec=self.SSBprec, BSGL=self.BSGL, SSBprec=self.SSBprec,
injectSources=self.injectSources, injectSources=self.injectSources,
...@@ -563,7 +577,7 @@ class FrequencySlidingWindow(GridSearch): ...@@ -563,7 +577,7 @@ class FrequencySlidingWindow(GridSearch):
self.search = ComputeFstat( self.search = ComputeFstat(
tref=self.tref, sftfilepattern=self.sftfilepattern, tref=self.tref, sftfilepattern=self.sftfilepattern,
minCoverFreq=self.minCoverFreq, maxCoverFreq=self.maxCoverFreq, minCoverFreq=self.minCoverFreq, maxCoverFreq=self.maxCoverFreq,
detectors=self.detectors, transient=True, detectors=self.detectors, transientWindowType=self.transientWindowType,
minStartTime=self.minStartTime, maxStartTime=self.maxStartTime, minStartTime=self.minStartTime, maxStartTime=self.maxStartTime,
BSGL=self.BSGL, SSBprec=self.SSBprec, BSGL=self.BSGL, SSBprec=self.SSBprec,
injectSources=self.injectSources) injectSources=self.injectSources)
......
...@@ -25,7 +25,8 @@ class Writer(BaseSearchClass): ...@@ -25,7 +25,8 @@ class Writer(BaseSearchClass):
tref=None, F0=30, F1=1e-10, F2=0, Alpha=5e-3, tref=None, F0=30, F1=1e-10, F2=0, Alpha=5e-3,
Delta=6e-2, h0=0.1, cosi=0.0, psi=0.0, phi=0, Tsft=1800, Delta=6e-2, h0=0.1, cosi=0.0, psi=0.0, phi=0, Tsft=1800,
outdir=".", sqrtSX=1, Band=4, detectors='H1', outdir=".", sqrtSX=1, Band=4, detectors='H1',
minStartTime=None, maxStartTime=None, add_noise=True): minStartTime=None, maxStartTime=None, add_noise=True,
transientWindowType='none'):
""" """
Parameters Parameters
---------- ----------
...@@ -92,7 +93,7 @@ class Writer(BaseSearchClass): ...@@ -92,7 +93,7 @@ class Writer(BaseSearchClass):
self.run_makefakedata() self.run_makefakedata()
def get_single_config_line(self, i, Alpha, Delta, h0, cosi, psi, phi, F0, def get_single_config_line(self, i, Alpha, Delta, h0, cosi, psi, phi, F0,
F1, F2, tref, tstart, duration_days): F1, F2, tref, window, tstart, duration_days):
template = ( template = (
"""[TS{}] """[TS{}]
Alpha = {:1.18e} Alpha = {:1.18e}
...@@ -105,11 +106,11 @@ Freq = {:1.18e} ...@@ -105,11 +106,11 @@ Freq = {:1.18e}
f1dot = {:1.18e} f1dot = {:1.18e}
f2dot = {:1.18e} f2dot = {:1.18e}
refTime = {:10.6f} refTime = {:10.6f}
transientWindowType=rect transientWindowType = {:s}
transientStartTime={:10.3f} transientStartTime = {:10.3f}
transientTauDays={:1.3f}\n""") transientTauDays = {:1.3f}\n""")
return template.format(i, Alpha, Delta, h0, cosi, psi, phi, F0, F1, return template.format(i, Alpha, Delta, h0, cosi, psi, phi, F0, F1,
F2, tref, tstart, duration_days) F2, tref, window, tstart, duration_days)
def make_cff(self): def make_cff(self):
""" """
...@@ -119,7 +120,7 @@ transientTauDays={:1.3f}\n""") ...@@ -119,7 +120,7 @@ transientTauDays={:1.3f}\n""")
content = self.get_single_config_line( content = self.get_single_config_line(
0, self.Alpha, self.Delta, self.h0, self.cosi, self.psi, 0, self.Alpha, self.Delta, self.h0, self.cosi, self.psi,
self.phi, self.F0, self.F1, self.F2, self.tref, self.tstart, self.phi, self.F0, self.F1, self.F2, self.tref, self.transientWindowType, self.tstart,
self.duration_days) self.duration_days)
if self.check_if_cff_file_needs_rewritting(content): if self.check_if_cff_file_needs_rewritting(content):
...@@ -247,7 +248,8 @@ class GlitchWriter(Writer): ...@@ -247,7 +248,8 @@ class GlitchWriter(Writer):
delta_F2=0, tref=None, F0=30, F1=1e-10, F2=0, Alpha=5e-3, delta_F2=0, tref=None, F0=30, F1=1e-10, F2=0, Alpha=5e-3,
Delta=6e-2, h0=0.1, cosi=0.0, psi=0.0, phi=0, Tsft=1800, Delta=6e-2, h0=0.1, cosi=0.0, psi=0.0, phi=0, Tsft=1800,
outdir=".", sqrtSX=1, Band=4, detectors='H1', outdir=".", sqrtSX=1, Band=4, detectors='H1',
minStartTime=None, maxStartTime=None, add_noise=True): minStartTime=None, maxStartTime=None, add_noise=True,
transientWindowType='rect'):
""" """
Parameters Parameters
---------- ----------
...@@ -317,7 +319,7 @@ class GlitchWriter(Writer): ...@@ -317,7 +319,7 @@ class GlitchWriter(Writer):
self.tbounds[:-1])): self.tbounds[:-1])):
line = self.get_single_config_line( line = self.get_single_config_line(
i, self.Alpha, self.Delta, self.h0, self.cosi, self.psi, i, self.Alpha, self.Delta, self.h0, self.cosi, self.psi,
t[0], t[1], t[2], t[3], self.tref, ts, d) t[0], t[1], t[2], t[3], self.tref, self.transientWindowType, ts, d)
content += line content += line
......
...@@ -76,6 +76,12 @@ class MCMCSearch(core.BaseSearchClass): ...@@ -76,6 +76,12 @@ class MCMCSearch(core.BaseSearchClass):
the search the search
assumeSqrtSX: float, optional assumeSqrtSX: float, optional
Don't estimate noise-floors, but assume (stationary) per-IFO sqrt{SX} Don't estimate noise-floors, but assume (stationary) per-IFO sqrt{SX}
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.
Attributes Attributes
---------- ----------
...@@ -108,7 +114,8 @@ class MCMCSearch(core.BaseSearchClass): ...@@ -108,7 +114,8 @@ class MCMCSearch(core.BaseSearchClass):
log10beta_min=-5, theta_initial=None, log10beta_min=-5, theta_initial=None,
rhohatmax=1000, binary=False, BSGL=False, rhohatmax=1000, binary=False, BSGL=False,
SSBprec=None, minCoverFreq=None, maxCoverFreq=None, SSBprec=None, minCoverFreq=None, maxCoverFreq=None,
injectSources=None, assumeSqrtSX=None): injectSources=None, assumeSqrtSX=None,
transientWindowType=None):
if os.path.isdir(outdir) is False: if os.path.isdir(outdir) is False:
os.mkdir(outdir) os.mkdir(outdir)
...@@ -150,7 +157,8 @@ class MCMCSearch(core.BaseSearchClass): ...@@ -150,7 +157,8 @@ class MCMCSearch(core.BaseSearchClass):
self.search = core.ComputeFstat( self.search = core.ComputeFstat(
tref=self.tref, sftfilepattern=self.sftfilepattern, tref=self.tref, sftfilepattern=self.sftfilepattern,
minCoverFreq=self.minCoverFreq, maxCoverFreq=self.maxCoverFreq, minCoverFreq=self.minCoverFreq, maxCoverFreq=self.maxCoverFreq,
detectors=self.detectors, BSGL=self.BSGL, transient=False, detectors=self.detectors, BSGL=self.BSGL,
transientWindowType=self.transientWindowType,
minStartTime=self.minStartTime, maxStartTime=self.maxStartTime, minStartTime=self.minStartTime, maxStartTime=self.maxStartTime,
binary=self.binary, injectSources=self.injectSources, binary=self.binary, injectSources=self.injectSources,
assumeSqrtSX=self.assumeSqrtSX, SSBprec=self.SSBprec) assumeSqrtSX=self.assumeSqrtSX, SSBprec=self.SSBprec)
...@@ -2195,10 +2203,12 @@ class MCMCTransientSearch(MCMCSearch): ...@@ -2195,10 +2203,12 @@ class MCMCTransientSearch(MCMCSearch):
def _initiate_search_object(self): def _initiate_search_object(self):
logging.info('Setting up search object') logging.info('Setting up search object')
if not self.transientWindowType:
self.transientWindowType = 'rect'
self.search = core.ComputeFstat( self.search = core.ComputeFstat(
tref=self.tref, sftfilepattern=self.sftfilepattern, tref=self.tref, sftfilepattern=self.sftfilepattern,
minCoverFreq=self.minCoverFreq, maxCoverFreq=self.maxCoverFreq, minCoverFreq=self.minCoverFreq, maxCoverFreq=self.maxCoverFreq,
detectors=self.detectors, transient=True, detectors=self.detectors, transientWindowType=self.transientWindowType,
minStartTime=self.minStartTime, maxStartTime=self.maxStartTime, minStartTime=self.minStartTime, maxStartTime=self.maxStartTime,
BSGL=self.BSGL, binary=self.binary, BSGL=self.BSGL, binary=self.binary,
injectSources=self.injectSources) injectSources=self.injectSources)
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment