diff --git a/pyfstat/core.py b/pyfstat/core.py index a5087659a588916ea08773d86ed5485e93f532b6..0c41b5c06c1b6b62e82de8bc02e6b7150c9bed09 100755 --- a/pyfstat/core.py +++ b/pyfstat/core.py @@ -330,7 +330,8 @@ class ComputeFstat(BaseSearchClass): @helper_functions.initializer 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, injectSources=None, injectSqrtSX=None, assumeSqrtSX=None, SSBprec=None): @@ -347,8 +348,11 @@ class ComputeFstat(BaseSearchClass): this epoch binary : bool If true, search of binary parameters. - transient : bool - If true, allow for the Fstat to be computed over a transient range. + transientWindowType: str + 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 If true, compute the BSGL rather than the twoF value. detectors : str @@ -477,7 +481,7 @@ class ComputeFstat(BaseSearchClass): logging.info('Initialising FstatInput') dFreq = 0 - if self.transient: + if self.transientWindowType: self.whatToCompute = lalpulsar.FSTATQ_ATOMS_PER_DET else: self.whatToCompute = lalpulsar.FSTATQ_2F @@ -593,14 +597,36 @@ class ComputeFstat(BaseSearchClass): self.whatToCompute = (self.whatToCompute + lalpulsar.FSTATQ_2F_PER_DET) - if self.transient: + if self.transientWindowType: logging.info('Initialising transient parameters') self.windowRange = lalpulsar.transientWindowRange_t() - self.windowRange.type = lalpulsar.TRANSIENT_RECTANGULAR - self.windowRange.t0Band = 0 - self.windowRange.dt0 = 1 - self.windowRange.tauBand = 0 - self.windowRange.dtau = 1 + transientWindowTypes = {'none': lalpulsar.TRANSIENT_NONE, + 'rect': lalpulsar.TRANSIENT_RECTANGULAR, + 'exp': lalpulsar.TRANSIENT_EXPONENTIAL} + if self.transientWindowType in transientWindowTypes: + 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, asini=None, period=None, ecc=None, tp=None, @@ -624,7 +650,7 @@ class ComputeFstat(BaseSearchClass): self.whatToCompute ) - if self.transient is False: + if not self.transientWindowType: if self.BSGL is False: return self.FstatResults.twoF[0] @@ -636,7 +662,13 @@ class ComputeFstat(BaseSearchClass): return log10_BSGL/np.log10(np.exp(1)) 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( self.FstatResults.multiFatoms[0], self.windowRange, False) @@ -696,8 +728,9 @@ class ComputeFstat(BaseSearchClass): max_tau = SFTmaxStartTime - tstart taus = np.linspace(min_tau, max_tau, npoints) twoFs = [] - if self.transient is False: - self.transient = True + if not self.transientWindowType: + # still call the transient-Fstat-map function, but using the full range + self.transientWindowType = 'none' self.init_computefstatistic_single_point() for tau in taus: twoFs.append(self.get_fullycoherent_twoF( @@ -868,7 +901,9 @@ class SemiCoherentSearch(ComputeFstat): self.fs_file_name = "{}/{}_FS.dat".format(self.outdir, self.label) self.set_ephemeris_files() - self.transient = True + self.transientWindowType = 'rect' + self.t0Band = None + self.tauBand = None self.init_computefstatistic_single_point() self.init_semicoherent_parameters() @@ -876,7 +911,7 @@ class SemiCoherentSearch(ComputeFstat): logging.info(('Initialising semicoherent parameters from {} to {} in' ' {} segments').format( self.minStartTime, self.maxStartTime, self.nsegs)) - self.transient = True + self.transientWindowType = 'rect' self.whatToCompute = lalpulsar.FSTATQ_2F+lalpulsar.FSTATQ_ATOMS_PER_DET self.tboundaries = np.linspace(self.minStartTime, self.maxStartTime, self.nsegs+1) @@ -915,7 +950,7 @@ class SemiCoherentSearch(ComputeFstat): self.whatToCompute ) - #if self.transient is False: + #if not self.transientWindowType: # if self.BSGL is False: # return self.FstatResults.twoF[0] # twoF = np.float(self.FstatResults.twoF[0]) @@ -1005,8 +1040,10 @@ class SemiCoherentGlitchSearch(ComputeFstat): self.fs_file_name = "{}/{}_FS.dat".format(self.outdir, self.label) self.set_ephemeris_files() - self.transient = True - self.binary = False + self.transientWindowType = 'rect' + self.t0Band = None + self.tauBand = None + self.binary = False self.init_computefstatistic_single_point() def get_semicoherent_nglitch_twoF(self, F0, F1, F2, Alpha, Delta, *args): diff --git a/pyfstat/grid_based_searches.py b/pyfstat/grid_based_searches.py index 4a38cbaab11544041e118d9da7ecda07aed11338..daad3b110333bd3b20c585d8cef3336cb3d25154 100644 --- a/pyfstat/grid_based_searches.py +++ b/pyfstat/grid_based_searches.py @@ -31,7 +31,8 @@ class GridSearch(BaseSearchClass): Deltas, tref=None, minStartTime=None, maxStartTime=None, nsegs=1, BSGL=False, minCoverFreq=None, maxCoverFreq=None, detectors=None, SSBprec=None, injectSources=None, - input_arrays=False, assumeSqrtSX=None): + input_arrays=False, assumeSqrtSX=None, + transientWindowType=None, t0Band=None, tauBand=None): """ Parameters ---------- @@ -48,6 +49,17 @@ class GridSearch(BaseSearchClass): GPS seconds of the reference time, start time and end time input_arrays: bool 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 """ @@ -66,7 +78,9 @@ class GridSearch(BaseSearchClass): self.search = ComputeFstat( tref=self.tref, sftfilepattern=self.sftfilepattern, 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, BSGL=self.BSGL, SSBprec=self.SSBprec, injectSources=self.injectSources, @@ -563,7 +577,7 @@ class FrequencySlidingWindow(GridSearch): self.search = ComputeFstat( tref=self.tref, sftfilepattern=self.sftfilepattern, minCoverFreq=self.minCoverFreq, maxCoverFreq=self.maxCoverFreq, - detectors=self.detectors, transient=True, + detectors=self.detectors, transientWindowType=self.transientWindowType, minStartTime=self.minStartTime, maxStartTime=self.maxStartTime, BSGL=self.BSGL, SSBprec=self.SSBprec, injectSources=self.injectSources) diff --git a/pyfstat/make_sfts.py b/pyfstat/make_sfts.py index bf4e258faafca5a30ed042b1da221de3d4a9a963..fda5b60b8d4c8458922ab87d58e7281ac0e0fd4c 100644 --- a/pyfstat/make_sfts.py +++ b/pyfstat/make_sfts.py @@ -25,7 +25,8 @@ class Writer(BaseSearchClass): 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, outdir=".", sqrtSX=1, Band=4, detectors='H1', - minStartTime=None, maxStartTime=None, add_noise=True): + minStartTime=None, maxStartTime=None, add_noise=True, + transientWindowType='none'): """ Parameters ---------- @@ -92,7 +93,7 @@ class Writer(BaseSearchClass): self.run_makefakedata() 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 = ( """[TS{}] Alpha = {:1.18e} @@ -105,11 +106,11 @@ Freq = {:1.18e} f1dot = {:1.18e} f2dot = {:1.18e} refTime = {:10.6f} -transientWindowType=rect -transientStartTime={:10.3f} -transientTauDays={:1.3f}\n""") +transientWindowType = {:s} +transientStartTime = {:10.3f} +transientTauDays = {:1.3f}\n""") 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): """ @@ -119,7 +120,7 @@ transientTauDays={:1.3f}\n""") content = self.get_single_config_line( 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) if self.check_if_cff_file_needs_rewritting(content): @@ -247,7 +248,8 @@ class GlitchWriter(Writer): 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, outdir=".", sqrtSX=1, Band=4, detectors='H1', - minStartTime=None, maxStartTime=None, add_noise=True): + minStartTime=None, maxStartTime=None, add_noise=True, + transientWindowType='rect'): """ Parameters ---------- @@ -317,7 +319,7 @@ class GlitchWriter(Writer): self.tbounds[:-1])): line = self.get_single_config_line( 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 diff --git a/pyfstat/mcmc_based_searches.py b/pyfstat/mcmc_based_searches.py index fa13b688b96061e5e3d00f671e777763f127a4b2..d663874c4e80f7c077c6a6929681a89346dc0e20 100644 --- a/pyfstat/mcmc_based_searches.py +++ b/pyfstat/mcmc_based_searches.py @@ -76,6 +76,12 @@ class MCMCSearch(core.BaseSearchClass): the search assumeSqrtSX: float, optional 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 ---------- @@ -108,7 +114,8 @@ class MCMCSearch(core.BaseSearchClass): log10beta_min=-5, theta_initial=None, rhohatmax=1000, binary=False, BSGL=False, SSBprec=None, minCoverFreq=None, maxCoverFreq=None, - injectSources=None, assumeSqrtSX=None): + injectSources=None, assumeSqrtSX=None, + transientWindowType=None): if os.path.isdir(outdir) is False: os.mkdir(outdir) @@ -150,7 +157,8 @@ class MCMCSearch(core.BaseSearchClass): self.search = core.ComputeFstat( tref=self.tref, sftfilepattern=self.sftfilepattern, 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, binary=self.binary, injectSources=self.injectSources, assumeSqrtSX=self.assumeSqrtSX, SSBprec=self.SSBprec) @@ -2195,10 +2203,12 @@ class MCMCTransientSearch(MCMCSearch): def _initiate_search_object(self): logging.info('Setting up search object') + if not self.transientWindowType: + self.transientWindowType = 'rect' self.search = core.ComputeFstat( tref=self.tref, sftfilepattern=self.sftfilepattern, minCoverFreq=self.minCoverFreq, maxCoverFreq=self.maxCoverFreq, - detectors=self.detectors, transient=True, + detectors=self.detectors, transientWindowType=self.transientWindowType, minStartTime=self.minStartTime, maxStartTime=self.maxStartTime, BSGL=self.BSGL, binary=self.binary, injectSources=self.injectSources)