From 9354cce33a5c092bd57c602e3f6e0b3fa4bea3b3 Mon Sep 17 00:00:00 2001 From: Alex Correia Date: Tue, 17 Dec 2024 15:21:32 -0500 Subject: [PATCH] add support for injection module --- pycbc/detector/__init__.py | 2 + pycbc/detector/ground.py | 18 +- pycbc/detector/space.py | 376 ++++++++++++++------------------ pycbc/distributions/__init__.py | 10 +- pycbc/inject/inject.py | 79 +++++-- 5 files changed, 240 insertions(+), 245 deletions(-) diff --git a/pycbc/detector/__init__.py b/pycbc/detector/__init__.py index 635932fb9ff..0562bc1c881 100644 --- a/pycbc/detector/__init__.py +++ b/pycbc/detector/__init__.py @@ -1,2 +1,4 @@ from .ground import * +from .ground import _ground_detectors from .space import * +from .space import _space_detectors diff --git a/pycbc/detector/ground.py b/pycbc/detector/ground.py index 31637c04aec..de02edf66d9 100644 --- a/pycbc/detector/ground.py +++ b/pycbc/detector/ground.py @@ -227,7 +227,7 @@ def __init__(self, detector_name, reference_time=1126259462.0): using a slower but higher precision method. """ self.name = str(detector_name) - + lal_detectors = [pfx for pfx, name in get_available_lal_detectors()] if detector_name in _ground_detectors: self.info = _ground_detectors[detector_name] @@ -247,6 +247,12 @@ def __init__(self, detector_name, reference_time=1126259462.0): self.sday = None self.gmst_reference = None + @property + def sky_coords(self): + """ List the sky coordinates used in response function. + """ + return 'ra', 'dec' + def set_gmst_reference(self): if self.reference_time is not None: self.sday = float(sday.si.scale) @@ -477,7 +483,8 @@ def time_delay_from_detector(self, other_detector, right_ascension, def project_wave(self, hp, hc, ra, dec, polarization, method='lal', - reference_time=None): + reference_time=None, + **kwargs): """Return the strain of a waveform as measured by the detector. Apply the time shift for the given detector relative to the assumed geocentric frame and apply the antenna patterns to the plus and cross @@ -663,10 +670,3 @@ def ppdets(ifos, separator=', '): if ifos: return separator.join(sorted(ifos)) return 'no detectors' - -__all__ = ['Detector', 'get_available_detectors', - 'get_available_lal_detectors', - 'gmst_accurate', 'add_detector_on_earth', - 'single_arm_frequency_response', 'ppdets', - 'overhead_antenna_pattern', 'load_detector_config', - '_ground_detectors',] diff --git a/pycbc/detector/space.py b/pycbc/detector/space.py index ac1233cb8e7..9fdce92e20a 100644 --- a/pycbc/detector/space.py +++ b/pycbc/detector/space.py @@ -67,12 +67,13 @@ def apply_polarization(hp, hc, polarization): return hp_ssb, hc_ssb -def check_signal_times(hp, hc, orbit_start_time, orbit_end_time, - offset=TIME_OFFSET_20_DEGREES, pad_data=False, t0=1e4): +def preprocess(hp, hc, orbit_start_time, orbit_end_time, + polarization=0, offset=TIME_OFFSET_20_DEGREES, + pad_data=False, t0=1e4): """ - Ensure that input signal lies within the provided orbital window. This - assumes that the start times of hp and hc are relative to the detector - mission start time. + Apply polarization and ensure that input signal lies within the + provided orbital window. This assumes that the start times of + hp and hc are relative to the detector mission start time. Parameters ---------- @@ -125,28 +126,31 @@ def check_signal_times(hp, hc, orbit_start_time, orbit_end_time, hc.prepend_zeros(pad_idx) hc.append_zeros(pad_idx) + # rotate GW from radiation frame to SSB using polarization angle + hp, hc = apply_polarization(hp, hc, polarization) + # make sure signal lies within orbit length if hp.duration + hp.start_time > orbit_end_time: - logging.warning('Time of signal end is greater than end of orbital ' + - f'data. Cutting signal at {orbit_end_time}.') + logging.warning("Time of signal end is greater than end of orbital data. " + + "Cutting signal at orbit end time.") # cut off data succeeding orbit end time - end_idx = numpy.argwhere(hp.sample_times.numpy() <= orbit_end_time)[-1][0] - hp = hp[:end_idx] - hc = hc[:end_idx] + orbit_end_idx = numpy.argwhere(hp.sample_times.numpy() <= orbit_end_time)[-1][0] + hp = hp[:orbit_end_idx] + hc = hc[:orbit_end_idx] if hp.start_time < orbit_start_time: - logging.warning('Time of signal start is less than start of orbital ' + - f'data. Cutting signal at {orbit_start_time}.') + logging.warning("Time of signal start is less than start of orbital data. " + + "Cutting signal at orbit start time.") # cut off data preceding orbit start time - start_idx = numpy.argwhere(hp.sample_times.numpy() >= orbit_start_time)[0][0] - hp = hp[start_idx:] - hc = hc[start_idx:] + orbit_start_idx = numpy.argwhere(hp.sample_times.numpy() >= orbit_start_time)[0][0] + hp = hp[orbit_start_idx:] + hc = hc[orbit_start_idx:] return hp, hc -def cut_channels(tdi_dict, remove_garbage=False, t0=1e4): +def postprocess(tdi_dict, remove_garbage=False, t0=1e4): """ - Cut TDI channels if needed. + Apply start times to TDI channels and cut if needed. Parameters ---------- @@ -177,100 +181,90 @@ def cut_channels(tdi_dict, remove_garbage=False, t0=1e4): slc = slice(pad_idx, -pad_idx) tdi_dict[chan] = tdi_dict[chan][slc] else: - raise ValueError('remove_garbage arg must be a bool or "zero"') + raise ValueError('remove_garbage arg must be a bool ' + + 'or "zero"') return tdi_dict -_space_detectors = {'LISA': {'armlength': 2.5e9, - 'aliases': ['LISA_A', 'LISA_E', 'LISA_T', - 'LISA_X', 'LISA_Y', 'LISA_Z'], - }, - } - class AbsSpaceDet(ABC): """ Abstract base class to set structure for space detector classes. Parameters ---------- - detector_name : str - The name of the detector. Accepts any output from - `get_available_space_detectors`. - reference_time : float (optional) - The reference time in seconds of the signal in the SSB frame. This is - defined such that the detector mission start time corresponds to 0. - Default None. + The reference time in seconds of the signal in the SSB frame. This + will correspond to the zero point of the input signal. By default, + a reference time is automatically assigned such that the input signal + starts at zero. Default None. + + apply_offset : bool (optional) + Flag whether to shift the times of the input waveforms by + a given value. Some backends require this such that the + detector is oriented correctly at t = 0. Default False. + + offset : float (optional) + The time in seconds by which to offset the input waveform if + apply_offset is True. Default 7365189.431698299. """ - def __init__(self, detector_name, reference_time=None, **kwargs): - self.det, self.chan = parse_det_name(detector_name) - if detector_name not in get_available_space_detectors(): - raise NotImplementedError('Unrecognized detector. ', - 'Currently accepts: ', - f'{get_available_space_detectors()}') - self.reference_time = reference_time + def __init__(self, apply_offset=False, + offset=TIME_OFFSET_20_DEGREES, **kwargs): + # specify whether to apply offsets to GPS times + if apply_offset: + self.offset = offset + else: + self.offset = 0. @property @abstractmethod def sky_coords(self): """ List the sky coordinate names for the detector class. + Raises a NotImplementedError if not specified in child class. + """ + return + + @abstractmethod + def orbits_init(self): + """ + Placeholder for initializing constellation orbital data. + Raises a NotImplementedError if not specified in child class. + """ + return + + @abstractmethod + def get_links(self): + """ + Placeholder for calculating GW projections onto detector links. + Raises a NotImplementedError if not specified in child class. """ return @abstractmethod - def project_wave(self, hp, hc, lamb, beta): + def project_wave(self): """ Placeholder for evaluating the TDI channels from the GW projections. + Raises a NotImplementedError if not specified in child class. """ return class _LDC_detector(AbsSpaceDet): """ - LISA detector modeled using LDC software. Constellation orbits are - generated using LISA Orbits (https://pypi.org/project/lisaorbits/). - Link projections are generated using LISA GW Response - (10.5281/zenodo.6423435). TDI channels are generated using pyTDI - (10.5281/zenodo.6351736). + LISA detector modeled using LDC software. Constellation orbits are generated + using LISA Orbits (https://pypi.org/project/lisaorbits/). Link projections + are generated using LISA GW Response (10.5281/zenodo.6423435). TDI channels + are generated using pyTDI (10.5281/zenodo.6351736). Parameters ---------- - detector_name : str - The name of the detector. Accepts any output from - `get_available_space_detectors`. - - reference_time : float (optional) - The reference time in seconds of the signal in the SSB frame. This is - defined such that the detector mission start time corresponds to 0. - Default None. - - apply_offset : bool (optional) - Flag whether to shift the times of the input waveforms by - a given value. Some backends require this such that the - detector is oriented correctly at t = 0. Default False. - - offset : float (optional) - The time in seconds by which to offset the input waveform if - apply_offset is True. Default 7365189.431698299. - orbits : str (optional) The constellation orbital data used for generating projections and TDI. See self.orbits_init for accepted inputs. Default 'EqualArmlength'. """ - def __init__(self, detector_name, reference_time=None, apply_offset=False, - offset=TIME_OFFSET_20_DEGREES, - orbits='EqualArmlength', **kwargs): - super().__init__(detector_name, reference_time, **kwargs) - assert self.det == 'LISA', 'LDC backend only works with LISA detector' - - # specify whether to apply offsets to GPS times - if apply_offset: - self.offset = offset - else: - self.offset = 0. - + def __init__(self, orbits='EqualArmlength', *args, **kwargs): + super().__init__(*args, **kwargs) # orbits properties self.orbits = orbits self.orbits_start_time = None @@ -290,8 +284,9 @@ def __init__(self, detector_name, reference_time=None, apply_offset=False, self.proj_init = None self.tdi_init = None self.tdi_chan = 'AET' - if self.chan is not None and self.chan in 'XYZ': - self.tdi_chan = 'XYZ' + if 'tdi_chan' in kwargs.keys(): + if kwargs['tdi_chan'] is not None and kwargs['tdi_chan'] in 'XYZ': + self.tdi_chan = 'XYZ' @property def sky_coords(self): @@ -299,9 +294,7 @@ def sky_coords(self): def orbits_init(self, orbits, size=316, dt=100000.0, t_init=0.0): """ - Initialize the orbital information for the constellation. Defualt args - generate roughly 1 year worth of data starting at LISA mission - start time. + Initialize the orbital information for the constellation. Parameters ---------- @@ -313,26 +306,26 @@ def orbits_init(self, orbits, size=316, dt=100000.0, t_init=0.0): length : int (optional) The number of samples to generate if creating a new orbit file. - Default 316. + Default 316; generates ~1 year worth of data. dt : float (optional) The time step in seconds to use if generating a new orbit file. - Default 100000. + Default 100000; generates ~1 year worth of data. t_init : float (optional) The start time in seconds to use if generating a new orbit file. - Default 0. + Default 0; generates data starting at the LISA mission start. """ defaults = ['EqualArmlength', 'Keplerian'] - assert type(orbits) == str, ('Must input either a file path as ', - 'str, "EqualArmlength", or "Keplerian"') + assert type(orbits) == str, ("Must input either a file path as a string, " + + "'EqualArmlength', or 'Keplerian'") # generate a new file if orbits in defaults: try: import lisaorbits except ImportError: - raise ImportError('lisaorbits not found') + raise ImportError('lisaorbits required if generating an orbits file') if orbits == 'EqualArmlength': o = lisaorbits.EqualArmlengthOrbits() if orbits == 'Keplerian': @@ -349,8 +342,7 @@ def orbits_init(self, orbits, size=316, dt=100000.0, t_init=0.0): ofile = orbits with h5py.File(ofile, 'r') as f: self.orbits_start_time = f.attrs['t0'] - self.orbits_end_time = self.orbit_start_time + \ - f.attrs['dt']*f.attrs['size'] + self.orbits_end_time = self.orbit_start_time + f.attrs['dt']*f.attrs['size'] # add light travel buffer times lisa_arm = _space_detectors['LISA']['armlength'] @@ -361,23 +353,23 @@ def orbits_init(self, orbits, size=316, dt=100000.0, t_init=0.0): def strain_container(self, response, orbits=None): """ - Read in the necessary link and orbit information for generating TDI - channels. Replicates the functionality of `pyTDI.Data.from_gws()`. + Read in the necessary link and orbit information for generating TDI channels. + Replicates the functionality of pyTDI.Data.from_gws(). Parameters ---------- response : array - The laser link projections of the GW. Uses get_links output format. + The laser link projections of the GW. Uses format of self.get_links() output. orbits : str, optional - The path to the file containing orbital information for the LISA - constellation. Default to orbits class attribute. + The path to the file containing orbital information for the LISA constellation. + Default to class attribute self.orbits. Returns ------- dict, array - The arguments and measurements associated with the link and orbital - data. + The arguments and measurements associated with the link and orbital data. + Can be passed into pyTDI.michelson.X1.build(**args)(measurements). """ try: from pytdi import Data @@ -433,22 +425,19 @@ def get_links(self, hp, hc, lamb, beta, polarization): try: from lisagwresponse import ReadStrain except ImportError: - raise ImportError('LISA GW Response not found') + raise ImportError('LISA GW Response required to generate projections') if self.dt is None: self.dt = hp.delta_t # configure orbits and signal self.orbits_init(orbits=self.orbits) - hp, hc = check_signal_times(hp, hc, self.orbits_start_time, - self.orbits_end_time, offset=self.offset, - pad_data=self.pad_data, t0=self.t0) + hp, hc = preprocess(hp, hc, self.orbits_start_time, self.orbits_end_time, + polarization=polarization, offset=self.offset, + pad_data=self.pad_data, t0=self.t0) self.start_time = hp.start_time - self.offset self.sample_times = hp.sample_times.numpy() - # apply polarization - hp, hc = apply_polarization(hp, hc, polarization) - if self.proj_init is None: # initialize the class self.proj_init = ReadStrain(self.sample_times, hp, hc, @@ -467,19 +456,18 @@ def get_links(self, hp, hc, lamb, beta, polarization): return wf_proj def project_wave(self, hp, hc, lamb, beta, polarization=0, - tdi=1.5, tdi_chan=None, pad_data=False, + tdi=1, tdi_chan=None, pad_data=False, remove_garbage=False, t0=1e4, **kwargs): """ Evaluate the TDI observables. - The TDI generation requires some startup time at the start and end of - the waveform, creating erroneous ringing or "garbage" at the edges of - the signal. By default, this method will cut off a time length t0 from - the start and end to remove this garbage, which may delete sensitive - data at the edges of the input strains (e.g., the late inspiral and - ringdown of a binary merger). Thus, the default output will be shorter - than the input by (2*t0) seconds. See pad_data and remove_garbage to - modify this behavior. + The TDI generation requires some startup time at the start and end of the + waveform, creating erroneous ringing or "garbage" at the edges of the signal. + By default, this method will cut off a time length t0 from the start and end + to remove this garbage, which may delete sensitive data at the edges of the + input data (e.g., the late inspiral and ringdown of a binary merger). Thus, + the default output will be shorter than the input by (2*t0) seconds. See + pad_data and remove_garbage to modify this behavior. Parameters ---------- @@ -498,9 +486,9 @@ def project_wave(self, hp, hc, lamb, beta, polarization=0, polarization : float The polarization angle of the GW in radians. - tdi : float (optional) - TDI channel configuration. Accepts 1.5 for 1st generation TDI or - 2 for 2nd generation TDI. Default 1.5. + tdi : int (optional) + TDI channel configuration. Accepts 1 for 1st generation TDI or + 2 for 2nd generation TDI. Default 2. tdi_chan : str (optional) The TDI observables to calculate. Accepts 'XYZ', 'AET', or 'AE'. @@ -530,10 +518,10 @@ def project_wave(self, hp, hc, lamb, beta, polarization=0, try: from pytdi import michelson except ImportError: - raise ImportError('pyTDI not found') + raise ImportError('pyTDI required for TDI combinations') # set TDI generation - if tdi == 1.5: + if tdi == 1: X, Y, Z = michelson.X1, michelson.Y1, michelson.Z1 elif tdi == 2: X, Y, Z = michelson.X2, michelson.Y2, michelson.Z2 @@ -549,8 +537,7 @@ def project_wave(self, hp, hc, lamb, beta, polarization=0, self.pad_data = pad_data self.remove_garbage = remove_garbage self.t0 = t0 - response = self.get_links(hp, hc, lamb, beta, - polarization=polarization) + response = self.get_links(hp, hc, lamb, beta, polarization=polarization) # load in data using response measurements self.tdi_init = self.strain_container(response, self.orbits) @@ -562,59 +549,33 @@ def project_wave(self, hp, hc, lamb, beta, polarization=0, # convert to AET if specified if tdi_chan == 'XYZ': - tdi_dict = {'LISA_X': TimeSeries(chanx, delta_t=self.dt, - epoch=self.start_time), - 'LISA_Y': TimeSeries(chany, delta_t=self.dt, - epoch=self.start_time), - 'LISA_Z': TimeSeries(chanz, delta_t=self.dt, - epoch=self.start_time)} + tdi_dict = {'LISA_X': TimeSeries(chanx, delta_t=self.dt, epoch=self.start_time), + 'LISA_Y': TimeSeries(chany, delta_t=self.dt, epoch=self.start_time), + 'LISA_Z': TimeSeries(chanz, delta_t=self.dt, epoch=self.start_time)} elif tdi_chan == 'AET': chana = (chanz - chanx)/numpy.sqrt(2) chane = (chanx - 2*chany + chanz)/numpy.sqrt(6) chant = (chanx + chany + chanz)/numpy.sqrt(3) - tdi_dict = {'LISA_A': TimeSeries(chana, delta_t=self.dt, - epoch=self.start_time), - 'LISA_E': TimeSeries(chane, delta_t=self.dt, - epoch=self.start_time), - 'LISA_T': TimeSeries(chant, delta_t=self.dt, - epoch=self.start_time)} + tdi_dict = {'LISA_A': TimeSeries(chana, delta_t=self.dt, epoch=self.start_time), + 'LISA_E': TimeSeries(chane, delta_t=self.dt, epoch=self.start_time), + 'LISA_T': TimeSeries(chant, delta_t=self.dt, epoch=self.start_time)} else: raise ValueError('Unrecognized TDI channel input. ' + 'Please input either "XYZ" or "AET".') # processing - tdi_dict = cut_channels(tdi_dict, remove_garbage=self.remove_garbage, - t0=self.t0) + tdi_dict = postprocess(tdi_dict, remove_garbage=self.remove_garbage, t0=self.t0) return tdi_dict class _FLR_detector(AbsSpaceDet): """ - LISA detector modeled using FastLISAResponse. Constellation orbits are - generated using LISA Analysis Tools (10.5281/zenodo.10930979). Link - projections and TDI channels are generated using FastLISAResponse - (https://arxiv.org/abs/2204.06633). + LISA detector modeled using FastLISAResponse. Constellation orbits are generated + using LISA Analysis Tools (10.5281/zenodo.10930979). Link projections and TDI + channels are generated using FastLISAResponse (https://arxiv.org/abs/2204.06633). Parameters ---------- - detector_name : str - The name of the detector. Accepts any output from - `get_available_space_detectors`. - - reference_time : float (optional) - The reference time in seconds of the signal in the SSB frame. This is - defined such that the detector mission start time corresponds to 0. - Default None. - - apply_offset : bool (optional) - Flag whether to shift the times of the input waveforms by - a given value. Some backends require this such that the - detector is oriented correctly at t = 0. Default False. - - offset : float (optional) - The time in seconds by which to offset the input waveform if - apply_offset is True. Default 7365189.431698299. - use_gpu : bool (optional) Specify whether to run class on GPU support via CuPy. Default False. @@ -623,21 +584,11 @@ class _FLR_detector(AbsSpaceDet): and TDI. See self.orbits_init for accepted inputs. Default 'EqualArmlength'. """ - def __init__(self, detector_name, reference_time=None, apply_offset=False, - offset=TIME_OFFSET_20_DEGREES, - orbits='EqualArmlength', use_gpu=False, **kwargs): - logging.warning('WARNING: FastLISAResponse TDI implementation is a ', - 'work in progress. Currently unable to reproduce LDC ', - 'or BBHx waveforms.') + def __init__(self, orbits='EqualArmlength', use_gpu=False, *args, **kwargs): + logging.warning('WARNING: FastLISAResponse TDI implementation is a work in progress. ', + 'Currently unable to reproduce LDC or BBHx waveforms.') + super().__init__(*args, **kwargs) self.use_gpu = use_gpu - super().__init__(detector_name, reference_time, **kwargs) - assert self.det == 'LISA', 'FLR backend only works with LISA detector' - - # specify whether to apply offsets to GPS times - if apply_offset: - self.offset = offset - else: - self.offset = 0. # orbits properties self.orbits = orbits @@ -660,7 +611,7 @@ def __init__(self, detector_name, reference_time=None, apply_offset=False, if 'tdi_chan' in kwargs.keys(): if kwargs['tdi_chan'] is not None and kwargs['tdi_chan'] in 'XYZ': self.tdi_chan = 'XYZ' - + @property def sky_coords(self): return 'eclipticlongitude', 'eclipticlatitude' @@ -728,7 +679,8 @@ def get_links(self, hp, hc, lamb, beta, polarization=0, use_gpu=None): use_gpu : bool (optional) Flag whether to use GPU support. Default to class input. - CuPy is required if use_gpu is True. + CuPy is required if use_gpu is True; an ImportError will be raised + if CuPy could not be imported. Returns ------- @@ -739,21 +691,18 @@ def get_links(self, hp, hc, lamb, beta, polarization=0, use_gpu=None): try: from fastlisaresponse import pyResponseTDI except ImportError: - raise ImportError('FastLISAResponse not found') + raise ImportError('FastLISAResponse required for LISA projection/TDI') if self.dt is None: self.dt = hp.delta_t # configure the orbit and signal self.orbits_init(orbits=self.orbits) - hp, hc = check_signal_times(hp, hc, self.orbits_start_time, - self.orbits_end_time, offset=self.offset, - pad_data=self.pad_data, t0=self.t0) + hp, hc = preprocess(hp, hc, self.orbits_start_time, self.orbits_end_time, + polarization=polarization, offset=self.offset, + pad_data=self.pad_data, t0=self.t0) self.start_time = hp.start_time - self.offset self.sample_times = hp.sample_times.numpy() - - # apply polarization - hp, hc = apply_polarization(hp, hc, polarization) # interpolate orbital data to signal sample times self.orbits.configure(t_arr=self.sample_times) @@ -773,8 +722,7 @@ def get_links(self, hp, hc, lamb, beta, polarization=0, use_gpu=None): if self.tdi_init is None: # initialize the class - self.tdi_init = pyResponseTDI(1/self.dt, len(wf), - orbits=self.orbits, + self.tdi_init = pyResponseTDI(1/self.dt, len(wf), orbits=self.orbits, use_gpu=use_gpu) else: # update params in the initialized class @@ -790,19 +738,18 @@ def get_links(self, hp, hc, lamb, beta, polarization=0, use_gpu=None): return wf_proj def project_wave(self, hp, hc, lamb, beta, polarization=0, - tdi=1.5, tdi_chan=None, use_gpu=None, pad_data=False, + tdi=1, tdi_chan=None, use_gpu=None, pad_data=False, remove_garbage=False, t0=1e4, **kwargs): """ Evaluate the TDI observables. - The TDI generation requires some startup time at the start and end of - the waveform, creating erroneous ringing or "garbage" at the edges of - the signal. By default, this method will cut off a time length t0 from - the start and end to remove this garbage, which may delete sensitive - data at the edges of the input strains (e.g., the late inspiral and - ringdown of a binary merger). Thus, the default output will be shorter - than the input by (2*t0) seconds. See pad_data and remove_garbage to - modify this behavior. + The TDI generation requires some startup time at the start and end of the + waveform, creating erroneous ringing or "garbage" at the edges of the signal. + By default, this method will cut off a time length t0 from the start and end + to remove this garbage, which may delete sensitive data at the edges of the + input data (e.g., the late inspiral and ringdown of a binary merger). Thus, + the default output will be shorter than the input by (2*t0) seconds. See + pad_data and remove_garbage to modify this behavior. Parameters ---------- @@ -821,9 +768,9 @@ def project_wave(self, hp, hc, lamb, beta, polarization=0, polarization : float (optional) The polarization angle of the GW in radians. - tdi : float(optional) - TDI channel configuration. Accepts 1.5 for 1st generation TDI or - 2 for 2nd generation TDI. Default 1.5. + tdi : int (optional) + TDI channel configuration. Accepts 1 for 1st generation TDI or + 2 for 2nd generation TDI. Default 1. tdi_chan : str (optional) The TDI observables to calculate. Accepts 'XYZ', 'AET', or 'AE'. @@ -841,7 +788,7 @@ def project_wave(self, hp, hc, lamb, beta, polarization=0, time length t0 worth of data at the start and end of the waveform will be cut from TDI channels. If 'zero', time length t0 worth of edge data will be zeroed. If False, TDI channels will not be - modified. Default False. + modified. Default True. t0 : float (optional) Time length in seconds to pad/cut from the start and end of @@ -865,7 +812,7 @@ def project_wave(self, hp, hc, lamb, beta, polarization=0, use_gpu=use_gpu) # set TDI configuration (let FLR handle if not 1 or 2) - if tdi == 1.5: + if tdi == 1: tdi_opt = '1st generation' elif tdi == 2: tdi_opt = '2nd generation' @@ -896,10 +843,15 @@ def project_wave(self, hp, hc, lamb, beta, polarization=0, tdi_dict[f'LISA_{chan}'] = TimeSeries(tdi_obs[i], delta_t=self.dt, epoch=self.start_time) - tdi_dict = cut_channels(tdi_dict, remove_garbage=self.remove_garbage, - t0=self.t0) + tdi_dict = postprocess(tdi_dict, remove_garbage=self.remove_garbage, + t0=self.t0) return tdi_dict +_space_detectors = {'LISA': {'armlength': 2.5e9, + 'aliases': ['LISA_A', 'LISA_E', 'LISA_T', + 'LISA_X', 'LISA_Y', 'LISA_Z'], + }, + } _backends = {'LISA': {'LDC': _LDC_detector, 'FLR': _FLR_detector, @@ -913,39 +865,35 @@ class SpaceDetector(AbsSpaceDet): Parameters ---------- detector_name : str - The name of the detector. Accepts any output from - `get_available_space_detectors`. - - reference_time : float (optional) - The reference time in seconds of the signal in the SSB frame. This is - defined such that the detector mission start time corresponds to 0. - Default None. + The detector name. Accepts 'LISA' or one of its aliases in _space_detectors. backend : str (optional) The backend architecture to use for generating TDI. Accepts 'LDC' or 'FLR'. Default 'LDC'. """ - def __init__(self, detector_name, reference_time=None, backend='LDC', - **kwargs): - super().__init__(detector_name, reference_time, **kwargs) - if backend in _backends[self.det].keys(): - c = _backends[self.det][backend] - self.backend = c(detector_name, reference_time, **kwargs) + def __init__(self, detector_name, backend='LDC', *args, **kwargs): + det, chan = parse_det_name(detector_name) + if detector_name in get_available_space_detectors(): + if backend in _backends[det].keys(): + c = _backends[det][backend] + kwargs['tdi_chan'] = chan + self.backend = c(*args, **kwargs) + else: + raise ValueError(f'Detector {det} does not support backend {backend}', + f'This detector accepts: {_backends[det].keys()}') else: - raise ValueError(f'Detector {self.det} does not support backend ', - f'{backend}.This detector accepts: ' - f'{_backends[self.det].keys()}') + raise NotImplementedError('Unrecognized detector. Currently accepts: ', + f'{get_available_space_detectors()}') @property def sky_coords(self): return self.backend.sky_coords + def orbits_init(self, *args, **kwargs): + return self.backend.orbits_init(*args, **kwargs) + def get_links(self, hp, hc, lamb, beta, *args, **kwargs): return self.backend.get_links(hp, hc, lamb, beta, *args, **kwargs) def project_wave(self, hp, hc, lamb, beta, *args, **kwargs): return self.backend.project_wave(hp, hc, lamb, beta, *args, **kwargs) - - -__all__ = ['get_available_space_detectors', 'SpaceDetector', - '_space_detectors',] \ No newline at end of file diff --git a/pycbc/distributions/__init__.py b/pycbc/distributions/__init__.py index 1ab45c41f24..c1430da9a08 100644 --- a/pycbc/distributions/__init__.py +++ b/pycbc/distributions/__init__.py @@ -170,9 +170,13 @@ def read_params_from_config(cp, prior_section='prior', # float (as we would expect for string arguments) static_args[key] = float(val) except ValueError: - # try converting to a list of strings; this function will just - # return val if it does not begin (end) with [ (]) - static_args[key] = _convert_liststring_to_list(val) + if val is None or val == '': + # convert None to True in cases where specifying `key = ` + static_args[key] = True + else: + # try converting to a list of strings; this function will just + # return val if it does not begin (end) with [ (]) + static_args[key] = _convert_liststring_to_list(val) return variable_args, static_args diff --git a/pycbc/inject/inject.py b/pycbc/inject/inject.py index 1e15c692df8..5801f21cf0e 100644 --- a/pycbc/inject/inject.py +++ b/pycbc/inject/inject.py @@ -41,8 +41,12 @@ from pycbc.waveform import utils as wfutils from pycbc.waveform import ringdown_td_approximants from pycbc.types import float64, float32, TimeSeries, load_timeseries -from pycbc.detector import Detector -from pycbc.conversions import tau0_from_mass1_mass2 +from pycbc.detector import (Detector, SpaceDetector, _space_detectors, + get_available_space_detectors, parse_det_name) +from pycbc.conversions import (tau0_from_mass1_mass2, + get_final_from_initial, + tau_from_final_mass_spin) +from pycbc.coordinates.space import TIME_OFFSET_20_DEGREES from pycbc.filter import resample_to_delta_t import pycbc.io from pycbc.io.ligolw import LIGOLWContentHandler @@ -85,19 +89,33 @@ def projector(detector_name, inj, hp, hc, distance_scale=1): """ Use the injection row to project the polarizations into the detector frame """ - detector = Detector(detector_name) + if detector_name in get_available_space_detectors(): + # FIXME : This should probably be moved to project_wave; + # wait on code reviews to see if this is OK here + apply_offset = False + offset = TIME_OFFSET_20_DEGREES + if 'apply_offset' in inj.dtype.names: + apply_offset = inj.apply_offset + if 'offset' in inj.dtype.names: + offset = inj.offset + detector = SpaceDetector(detector_name, + apply_offset=apply_offset, + offset=offset) + else: + detector = Detector(detector_name) hp /= distance_scale hc /= distance_scale try: tc = inj.tc - ra = inj.ra - dec = inj.dec + sky1_name, sky2_name = detector.sky_coords + sky1 = inj[sky1_name] + sky2 = inj[sky2_name] except: tc = inj.time_geocent - ra = inj.longitude - dec = inj.latitude + sky1 = inj.longitude + sky2 = inj.latitude hp.start_time += tc hc.start_time += tc @@ -117,10 +135,15 @@ def projector(detector_name, inj, hp, hc, distance_scale=1): logger.info('Injecting at %s, method is %s', tc, projection_method) # compute the detector response and add it to the strain + extra_args = dict(method = projection_method, + reference_time = tc) + pos_args = [sky1_name, sky2_name, 'polarization'] + for field in inj.dtype.names: + if field not in pos_args: + extra_args[field] = inj[field] signal = detector.project_wave(hp_tapered, hc_tapered, - ra, dec, inj.polarization, - method=projection_method, - reference_time=tc,) + sky1, sky2, inj.polarization, + **extra_args) return signal def legacy_approximant_name(apx): @@ -541,6 +564,14 @@ def apply(self, strain, detector_name, f_lower=None, distance_scale=1, if self.table[0]['approximant'] in fd_det: t0 = float(strain.start_time) t1 = float(strain.end_time) + elif detector_name in get_available_space_detectors(): + au_travel_time = lal.AU_SI / lal.C_SI + try: + arm_travel_time = _space_detectors[detector_name]['armlength'] / lal.C_SI + except: + arm_travel_time = 2.5e9 / lal.C_SI # default to LISA + t0 = float(strain.start_time) - au_travel_time - arm_travel_time + t1 = float(strain.end_time) + au_travel_time + arm_travel_time else: earth_travel_time = lal.REARTH_SI / lal.C_SI t0 = float(strain.start_time) - earth_travel_time @@ -561,12 +592,15 @@ def apply(self, strain, detector_name, f_lower=None, distance_scale=1, for ii, inj in enumerate(injections): f_l = inj.f_lower if f_lower is None else f_lower # roughly estimate if the injection may overlap with the segment - # Add 2s to end_time to account for ringdown and light-travel delay - end_time = inj.tc + 2 + # Add to end_time to account for ringdown and light-travel delay + # The ringdown buffer corresponds to ~2s for a final mass of 70M + mf, chif = get_final_from_initial(inj.mass1, inj.mass2) + rd_buffer = 500*tau_from_final_mass_spin(mf, chif) + end_time = inj.tc + rd_buffer inj_length = tau0_from_mass1_mass2(inj.mass1, inj.mass2, f_l) # Start time is taken as twice approx waveform length with a 1s # safety buffer - start_time = inj.tc - 2 * (inj_length + 1) + start_time = inj.tc - 2 * (inj_length + rd_buffer/2) if end_time < t0 or start_time > t1: continue signal = self.make_strain_from_inj_object(inj, delta_t, @@ -638,8 +672,13 @@ def make_strain_from_inj_object(self, inj, delta_t, detector_name, # compute the waveform time series hp, hc = get_td_waveform(inj, delta_t=delta_t, f_lower=f_l, **self.extra_args) - strain = projector(detector_name, - inj, hp, hc, distance_scale=distance_scale) + if detector_name in get_available_space_detectors(): + # call a specific TDI channel + strain = projector(detector_name, inj, hp, hc, + distance_scale=distance_scale)[detector_name] + else: + strain = projector(detector_name, inj, hp, hc, + distance_scale=distance_scale) return strain def end_times(self): @@ -757,13 +796,15 @@ def make_strain_from_inj_object(self, inj, delta_t, detector_name, def end_times(self): """Return the approximate end times of all injections. - - Currently, this just assumes all ringdowns are 2 seconds long. """ # the start times are the tcs tcs = self.table.tc - # FIXME: this could be figured out using the ringdown module - return tcs + 2 + mfs = self.table.final_mass + chifs = self.table.final_spin + # end time is calculated from damping time such that a 70M final mass + # returns a ringdown roughly 2s long + rd_buffs = 500*tau_from_final_mass_spin(mfs, chifs) + return tcs + rd_buffs @staticmethod def supported_approximants():