diff --git a/pyfstat/__init__.py b/pyfstat/__init__.py index 852df7b9940a2b55f790d09d31ea477cd7084e17..2d8f4d545a705c658ba215cf645a064e12292b5d 100644 --- a/pyfstat/__init__.py +++ b/pyfstat/__init__.py @@ -1,6 +1,6 @@ from __future__ import division as _division from .core import BaseSearchClass, ComputeFstat, SemiCoherentSearch, SemiCoherentGlitchSearch -from .make_sfts import Writer +from .make_sfts import Writer, GlitchWriter from .mcmc_based_searches import MCMCSearch, MCMCGlitchSearch, MCMCSemiCoherentSearch, MCMCFollowUpSearch, MCMCTransientSearch from .grid_based_searches import GridSearch, GridUniformPriorSearch, GridGlitchSearch, FrequencySlidingWindow, DMoff_NO_SPIN diff --git a/pyfstat/make_sfts.py b/pyfstat/make_sfts.py index 96ac3faaa9da8fd41cd9db2afa8222cb6d52fd17..53dc124a8934140f207f358b0765c64927e7a3cc 100644 --- a/pyfstat/make_sfts.py +++ b/pyfstat/make_sfts.py @@ -11,11 +11,10 @@ import helper_functions class Writer(BaseSearchClass): - """ Instance object for generating SFTs containing glitch signals """ + """ Instance object for generating SFTs """ @helper_functions.initializer def __init__(self, label='Test', tstart=700000000, duration=100*86400, - dtglitch=None, delta_phi=0, delta_F0=0, delta_F1=0, - delta_F2=0, 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, outdir=".", sqrtSX=1, Band=4, detectors='H1', minStartTime=None, maxStartTime=None, add_noise=True): @@ -26,11 +25,6 @@ class Writer(BaseSearchClass): a human-readable label to be used in naming the output files tstart, tend : float start and end times (in gps seconds) of the total observation span - dtglitch: float - time (in gps seconds) of the glitch after tstart. To create data - without a glitch, set dtglitch=None - delta_phi, delta_F0, delta_F1: float - instanteneous glitch magnitudes in rad, Hz, and Hz/s respectively tref: float or None reference time (default is None, which sets the reference time to tstart) @@ -45,21 +39,12 @@ class Writer(BaseSearchClass): see `lalapps_Makefakedata_v5 --help` for help with the other paramaters """ - for d in self.delta_phi, self.delta_F0, self.delta_F1, self.delta_F2: - if np.size(d) == 1: - d = np.atleast_1d(d) self.tend = self.tstart + self.duration if self.minStartTime is None: self.minStartTime = self.tstart if self.maxStartTime is None: self.maxStartTime = self.tend - if self.dtglitch is None: - self.tbounds = [self.tstart, self.tend] - else: - self.dtglitch = np.atleast_1d(self.dtglitch) - self.tglitch = self.tstart + self.dtglitch - self.tbounds = np.concatenate(( - [self.tstart], self.tglitch, [self.tend])) + self.tbounds = [self.tstart, self.tend] logging.info('Using segment boundaries {}'.format(self.tbounds)) self.check_inputs() @@ -69,14 +54,9 @@ class Writer(BaseSearchClass): if self.tref is None: self.tref = self.tstart self.tend = self.tstart + self.duration - tbs = np.array(self.tbounds) - self.durations_days = (tbs[1:] - tbs[:-1]) / 86400 + self.duration_days = (self.tend-self.tstart) / 86400 self.config_file_name = "{}/{}.cff".format(outdir, label) - self.theta = np.array([phi, F0, F1, F2]) - self.delta_thetas = np.atleast_2d( - np.array([delta_phi, delta_F0, delta_F1, delta_F2]).T) - self.data_duration = self.maxStartTime - self.minStartTime numSFTs = int(float(self.data_duration) / self.Tsft) self.sftfilenames = [ @@ -91,12 +71,6 @@ class Writer(BaseSearchClass): def check_inputs(self): self.minStartTime = int(self.minStartTime) self.maxStartTime = int(self.maxStartTime) - shapes = np.array([np.shape(x) for x in [self.delta_phi, self.delta_F0, - self.delta_F1, self.delta_F2]] - ) - if not np.all(shapes == shapes[0]): - raise ValueError('all delta_* must be the same shape: {}'.format( - shapes)) def make_data(self): ''' A convienience wrapper to generate a cff file then sfts ''' @@ -125,21 +99,14 @@ transientTauDays={:1.3f}\n""") def make_cff(self): """ - Generates an .cff file for a 'glitching' signal + Generates a .cff file """ - thetas = self._calculate_thetas(self.theta, self.delta_thetas, - self.tbounds) - - content = '' - for i, (t, d, ts) in enumerate(zip(thetas, self.durations_days, - 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) - - content += line + 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.duration_days) if self.check_if_cff_file_needs_rewritting(content): config_file = open(self.config_file_name, "w+") @@ -272,3 +239,115 @@ transientTauDays={:1.3f}\n""") output = helper_functions.run_commandline(cl_pfs) twoF = float(output.split('\n')[-2]) return float(twoF) + + +class GlitchWriter(Writer): + """ Instance object for generating SFTs containing glitch signals """ + @helper_functions.initializer + def __init__(self, label='Test', tstart=700000000, duration=100*86400, + dtglitch=None, delta_phi=0, delta_F0=0, delta_F1=0, + 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): + """ + Parameters + ---------- + label: string + a human-readable label to be used in naming the output files + tstart, tend : float + start and end times (in gps seconds) of the total observation span + dtglitch: float + time (in gps seconds) of the glitch after tstart. To create data + without a glitch, set dtglitch=None + delta_phi, delta_F0, delta_F1: float + instanteneous glitch magnitudes in rad, Hz, and Hz/s respectively + tref: float or None + reference time (default is None, which sets the reference time to + tstart) + F0, F1, F2, Alpha, Delta, h0, cosi, psi, phi: float + frequency, sky-position, and amplitude parameters + Tsft: float + the sft duration + minStartTime, maxStartTime: float + if not None, the total span of data, this can be used to generate + transient signals + + see `lalapps_Makefakedata_v5 --help` for help with the other paramaters + """ + + for d in self.delta_phi, self.delta_F0, self.delta_F1, self.delta_F2: + if np.size(d) == 1: + d = np.atleast_1d(d) + self.tend = self.tstart + self.duration + if self.minStartTime is None: + self.minStartTime = self.tstart + if self.maxStartTime is None: + self.maxStartTime = self.tend + if self.dtglitch is None: + self.tbounds = [self.tstart, self.tend] + else: + self.dtglitch = np.atleast_1d(self.dtglitch) + self.tglitch = self.tstart + self.dtglitch + self.tbounds = np.concatenate(( + [self.tstart], self.tglitch, [self.tend])) + logging.info('Using segment boundaries {}'.format(self.tbounds)) + + self.check_inputs() + + if os.path.isdir(self.outdir) is False: + os.makedirs(self.outdir) + if self.tref is None: + self.tref = self.tstart + self.tend = self.tstart + self.duration + tbs = np.array(self.tbounds) + self.durations_days = (tbs[1:] - tbs[:-1]) / 86400 + self.config_file_name = "{}/{}.cff".format(outdir, label) + + self.theta = np.array([phi, F0, F1, F2]) + self.delta_thetas = np.atleast_2d( + np.array([delta_phi, delta_F0, delta_F1, delta_F2]).T) + + self.data_duration = self.maxStartTime - self.minStartTime + numSFTs = int(float(self.data_duration) / self.Tsft) + self.sftfilenames = [ + lalpulsar.OfficialSFTFilename( + dets[0], dets[1], numSFTs, self.Tsft, self.minStartTime, + self.data_duration, self.label) + for dets in self.detectors.split(',')] + self.sftfilepath = ';'.join([ + '{}/{}'.format(self.outdir, fn) for fn in self.sftfilenames]) + self.calculate_fmin_Band() + + def check_inputs(self): + self.minStartTime = int(self.minStartTime) + self.maxStartTime = int(self.maxStartTime) + shapes = np.array([np.shape(x) for x in [self.delta_phi, self.delta_F0, + self.delta_F1, self.delta_F2]] + ) + if not np.all(shapes == shapes[0]): + raise ValueError('all delta_* must be the same shape: {}'.format( + shapes)) + + def make_cff(self): + """ + Generates an .cff file for a 'glitching' signal + + """ + + thetas = self._calculate_thetas(self.theta, self.delta_thetas, + self.tbounds) + + content = '' + for i, (t, d, ts) in enumerate(zip(thetas, self.durations_days, + 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) + + content += line + + if self.check_if_cff_file_needs_rewritting(content): + config_file = open(self.config_file_name, "w+") + config_file.write(content) + config_file.close() diff --git a/tests.py b/tests.py index 86d2a9934012b1be9088de86a6fe80bca0b0e66d..0859151bbdb53a039580c7df0055a34215e91549 100644 --- a/tests.py +++ b/tests.py @@ -154,9 +154,9 @@ class TestSemiCoherentGlitchSearch(Test): duration = 100*86400 dtglitch = .5*100*86400 delta_F0 = 0 - Writer = pyfstat.Writer(self.label, outdir=outdir, - duration=duration, dtglitch=dtglitch, - delta_F0=delta_F0) + Writer = pyfstat.GlitchWriter( + self.label, outdir=outdir, duration=duration, dtglitch=dtglitch, + delta_F0=delta_F0) Writer.make_data() @@ -204,13 +204,12 @@ class TestMCMCSearch(Test): Alpha = 5e-3 Delta = 1.2 tref = minStartTime - delta_F0 = 0 Writer = pyfstat.Writer(F0=F0, F1=F1, F2=F2, label=self.label, h0=h0, sqrtSX=sqrtSX, outdir=outdir, tstart=minStartTime, Alpha=Alpha, Delta=Delta, tref=tref, duration=duration, - delta_F0=delta_F0, Band=4) + Band=4) Writer.make_data() predicted_FS = Writer.predict_fstat()