From db216728b75971709dc4786a7a8de975f01a5b2b Mon Sep 17 00:00:00 2001 From: acorreia61201 Date: Thu, 4 Apr 2024 18:01:19 +0000 Subject: [PATCH 01/24] Added LISA detector class with GPU support --- pycbc/detector.py | 260 ++++++++++++++++++++++++++++------------------ 1 file changed, 160 insertions(+), 100 deletions(-) diff --git a/pycbc/detector.py b/pycbc/detector.py index 261b0eb5feb..70f9ae027b4 100644 --- a/pycbc/detector.py +++ b/pycbc/detector.py @@ -649,132 +649,192 @@ def overhead_antenna_pattern(right_ascension, declination, polarization): """ LISA class """ +try: + import cupy +except ImportError: + cupy = None -class LISA(object): - """For LISA detector +class LISA_detector(object): """ - def __init__(self): - None + LISA-like GW detector. Applies detector response from FastLISAResponse. + """ + def __init__(self, orbit_kwargs, use_gpu=False): + """ + Parameters + ---------- + orbit_kwargs : dict + List of keywords for generating LISA orbital data. Passed into pyResponseTDI + instance. Keywords and defaults as set by lisaresponse are: + + 'orbit_module' = None + 'order' = 0 + 'max_t_orbits' = 3.15576e7 + 'orbit_file' = None + + Either an orbit file or orbit module is required for projection. + If both a file and module are provided, the file will take priority. + If neither are provided, a ValueError is raised. + + use_gpu : bool (optional) + Specify whether to run class on GPU support. Default False. + """ + # intialize orbit kwargs to FLR defaults + self.orbit_kwargs = {'orbit_module': None, 'order': 0, 'max_t_orbits': 3.15576e7, + 'orbit_file': None} + + # get kwargs from class input + for key, val in orbit_kwargs.items(): + self.orbit_kwargs[key] = val + + # handle if both/neither module and file are specified + if self.orbit_kwargs['orbit_module'] is not None: + if self.orbit_kwargs['orbit_file'] is not None: + # FLR defaults to the module; override this and default to the file + self.orbit_kwargs['orbit_module'] = None + else: + if self.orbit_kwargs['orbit_file'] is None: + raise ValueError('Either orbit file or orbit module required to initialize LISA detector') + + # cache the FLR instance along with dt and n + self.dt = None + self.n = None + self.response_init = None + self.t0 = None + + # initialize whether to use gpu; FLR has handles if this cannot be done + self.use_gpu = use_gpu + + def project_wf(self, hp, hc, lamb, beta, t0=10000., use_gpu=None): + """ + Project a radiation frame (SSB) waveform to the LISA constellation. + This *does not* perform TDI; this only produces a detector frame + representation of the SSB frame signal. - def get_pos(self, ref_time): - """Return the position of LISA detector for a given reference time Parameters ---------- - ref_time : numpy.ScalarType + hp : pycbc.types.TimeSeries + The plus polarization of the GW in the radiation frame. + + hc : pycbc.types.TimeSeries + The cross polarization of the GW in the radiation frame. + + lambda : float + The ecleptic latitude in the SSB frame. + + beta : float + The ecleptic longitude in the SSB frame. + + t0 : float (optional) + Number of seconds to omit from start and end of waveform. Defaults + to FastLISAResponse preset of 10000 s. + + use_gpu : bool (optional) + Flag whether to use GPU support. Default to class input. + CuPy is required if use_gpu == True; an ImportError will be raised + if CuPy could not be imported. Returns ------- - location : numpy.ndarray of shape (3,3) - Returns the position of all 3 sattelites with each row - correspoding to a single axis. + ndarray + The waveform projected to the LISA laser links. """ - ref_time = Time(val=ref_time, format='gps', scale='utc').jyear - n = np.array(range(1, 4)) - kappa, _lambda_ = 0, 0 - alpha = 2. * np.pi * ref_time/1 + kappa - beta_n = (n - 1) * 2.0 * pi / 3.0 + _lambda_ - a, L = 1., 0.03342293561 - e = L/(2. * a * np.sqrt(3)) + from fastlisaresponse import pyResponseTDI - x = a*cos(alpha) + a*e*(sin(alpha)*cos(alpha)*sin(beta_n) - (1 + sin(alpha)**2)*cos(beta_n)) - y = a*sin(alpha) + a*e*(sin(alpha)*cos(alpha)*cos(beta_n) - (1 + cos(alpha)**2)*sin(beta_n)) - z = -np.sqrt(3)*a*e*cos(alpha - beta_n) - self.location = np.array([x, y, z]) + s_to_yr = 3.168808781402895e-08 # seconds to years conversion factor + dt = hp.delta_t # timestep (s) + n = len(hp) # number of points - return self.location + # format wf to FLR standard hp + i*hc as ndarray + hp = hp.numpy() + hc = hc.numpy() + wf = hp + 1j*hc - def get_gcrs_pos(self, location): - """ Transforms ICRS frame to GCRS frame + # set use_gpu to class input if not specified + if use_gpu is None: + use_gpu = self.use_gpu - Parameters - ---------- - loc : numpy.ndarray shape (3,1) units: AU - Cartesian Coordinates of the location - in ICRS frame + # convert to cupy if needed + if use_gpu: + if cupy is None: + raise ImportError('CuPy not imported but is required for GPU usage') + else: + wf = cupy.asarray(wf) - Returns - ---------- - loc : numpy.ndarray shape (3,1) units: meters - GCRS coordinates in cartesian system - """ - loc = location - loc = coordinates.SkyCoord(x=loc[0], y=loc[1], z=loc[2], unit=units.AU, - frame='icrs', representation_type='cartesian').transform_to('gcrs') - loc.representation_type = 'cartesian' - conv = np.float32(((loc.x.unit/units.m).decompose()).to_string()) - loc = np.array([np.float32(loc.x), np.float32(loc.y), - np.float32(loc.z)])*conv - return loc + # initialize the FLR class + response_init = pyResponseTDI(dt, n, self.orbit_kwargs, use_gpu=use_gpu) - def time_delay_from_location(self, other_location, right_ascension, - declination, t_gps): - """Return the time delay from the LISA detector to detector for - a signal with the given sky location. In other words return - `t1 - t2` where `t1` is the arrival time in this detector and - `t2` is the arrival time in the other location. Units(AU) + # cache n, dt, response init + self.dt = dt + self.n = n + self.response_init = response_init + self.t0 = t0 - Parameters - ---------- - other_location : numpy.ndarray of coordinates in ICRS frame - A detector instance. - right_ascension : float - The right ascension (in rad) of the signal. - declination : float - The declination (in rad) of the signal. - t_gps : float - The GPS time (in s) of the signal. + # project the signal + response_init.get_projections(wf, lamb, beta, t0) + wf_proj = response_init.y_gw + + return wf_proj - Returns - ------- - numpy.ndarray - The arrival time difference between the detectors. + def get_tdi(self, hp, hc, lamb, beta, tdi='XYZ', tdi_orbit_kwargs=None, use_gpu=None): """ - dx = self.location - other_location - cosd = cos(declination) - e0 = cosd * cos(right_ascension) - e1 = cosd * -sin(right_ascension) - e2 = sin(declination) - ehat = np.array([e0, e1, e2]) - return dx.dot(ehat) / constants.c.value - - def time_delay_from_detector(self, det, right_ascension, - declination, t_gps): - """Return the time delay from the LISA detector for a signal with - the given sky location in ICRS frame; i.e. return `t1 - t2` where - `t1` is the arrival time in this detector and `t2` is the arrival - time in the other detector. + Calculate the TDI variables given a waveform and LISA orbit data. Parameters ---------- - other_detector : detector.Detector - A detector instance. - right_ascension : float - The right ascension (in rad) of the signal. - declination : float - The declination (in rad) of the signal. - t_gps : float - The GPS time (in s) of the signal. + hp : pycbc.types.TimeSeries + The plus polarization of the GW in the radiation frame. + + hc : pycbc.types.TimeSeries + The cross polarization of the GW in the radiation frame. + + lambda : float + The ecleptic latitude in the SSB frame. + + beta : float + The ecleptic longitude in the SSB frame. + + tdi : str (optional) + The TDI observables to calculate. Accepts 'XYZ', 'AET', or 'AE'. + Default 'XYZ'. + + tdi_orbit_kwargs : dict (optional) + Orbit keywords specifically for TDI generation. If None, use + orbit_kwargs provided on class initialization. Default None. + + use_gpu : bool (optional) + Flag whether to use GPU support. Default to class input. Returns ------- - numpy.ndarray - The arrival time difference between the detectors. + dict + The TDI observables keyed by their corresponding TDI channel. """ - loc = Detector(det, t_gps).get_icrs_pos() - return self.time_delay_from_location(loc, right_ascension, - declination, t_gps) - def time_delay_from_earth_center(self, right_ascension, declination, t_gps): - """Return the time delay from the earth center in ICRS frame - """ - t_gps = Time(val=t_gps, format='gps', scale='utc') - earth = coordinates.get_body('earth', t_gps, - location=None).transform_to('icrs') - earth.representation_type = 'cartesian' - return self.time_delay_from_location( - np.array([np.float32(earth.x), np.float32(earth.y), - np.float32(earth.z)]), right_ascension, - declination, t_gps) + # set use_gpu + if use_gpu is None: + use_gpu = self.use_gpu + + # generate the Doppler time series + links = self.project_wf(hp, hc, lamb, beta, t0=self.t0, use_gpu=use_gpu) + + # if TDI kwargs are provided, update the response_init + if type(tdi_orbit_kwargs) == dict: + self.response_init.tdi_orbit_kwargs = tdi_orbit_kwargs + + # set TDI channels + if tdi in ['XYZ', 'AET', 'AE']: + self.response_init.tdi_chan = tdi + else: + raise ValueError('TDI channels must be one of: XYZ, AET, AE') + + tdi_obs = self.response_init.get_tdi_delays() + + # convert to dict + tdi_dict = {} + for i in range(len(tdi)): + tdi_dict[tdi[i]] = tdi_obs[i] + + return tdi_dict def ppdets(ifos, separator=', '): From 5d5e16b1448eb78aa4bd87416faa0df3e8b99c31 Mon Sep 17 00:00:00 2001 From: acorreia61201 Date: Thu, 4 Apr 2024 20:04:59 +0000 Subject: [PATCH 02/24] Fixed bug when calling get_tdi w/out setting t0 --- pycbc/detector.py | 28 ++++++++++++++++++---------- 1 file changed, 18 insertions(+), 10 deletions(-) diff --git a/pycbc/detector.py b/pycbc/detector.py index 70f9ae027b4..68e75e1e6a6 100644 --- a/pycbc/detector.py +++ b/pycbc/detector.py @@ -699,12 +699,12 @@ def __init__(self, orbit_kwargs, use_gpu=False): self.dt = None self.n = None self.response_init = None - self.t0 = None + self.t0 = 10000 # initialize whether to use gpu; FLR has handles if this cannot be done self.use_gpu = use_gpu - def project_wf(self, hp, hc, lamb, beta, t0=10000., use_gpu=None): + def project_wf(self, hp, hc, lamb, beta, t0=None, use_gpu=None): """ Project a radiation frame (SSB) waveform to the LISA constellation. This *does not* perform TDI; this only produces a detector frame @@ -725,8 +725,8 @@ def project_wf(self, hp, hc, lamb, beta, t0=10000., use_gpu=None): The ecleptic longitude in the SSB frame. t0 : float (optional) - Number of seconds to omit from start and end of waveform. Defaults - to FastLISAResponse preset of 10000 s. + Number of seconds to omit from start and end of waveform. Default + to FLR default of 10000 s. use_gpu : bool (optional) Flag whether to use GPU support. Default to class input. @@ -740,7 +740,6 @@ def project_wf(self, hp, hc, lamb, beta, t0=10000., use_gpu=None): """ from fastlisaresponse import pyResponseTDI - s_to_yr = 3.168808781402895e-08 # seconds to years conversion factor dt = hp.delta_t # timestep (s) n = len(hp) # number of points @@ -759,23 +758,24 @@ def project_wf(self, hp, hc, lamb, beta, t0=10000., use_gpu=None): raise ImportError('CuPy not imported but is required for GPU usage') else: wf = cupy.asarray(wf) - + # initialize the FLR class response_init = pyResponseTDI(dt, n, self.orbit_kwargs, use_gpu=use_gpu) - # cache n, dt, response init + # cache response init self.dt = dt self.n = n self.response_init = response_init - self.t0 = t0 + if t0 is not None: + self.t0 = t0 # project the signal - response_init.get_projections(wf, lamb, beta, t0) + response_init.get_projections(wf, lamb, beta, self.t0) wf_proj = response_init.y_gw return wf_proj - def get_tdi(self, hp, hc, lamb, beta, tdi='XYZ', tdi_orbit_kwargs=None, use_gpu=None): + def get_tdi(self, hp, hc, lamb, beta, t0=None, tdi='XYZ', tdi_orbit_kwargs=None, use_gpu=None): """ Calculate the TDI variables given a waveform and LISA orbit data. @@ -792,6 +792,10 @@ def get_tdi(self, hp, hc, lamb, beta, tdi='XYZ', tdi_orbit_kwargs=None, use_gpu= beta : float The ecleptic longitude in the SSB frame. + + t0 : float (optional) + Number of seconds to omit from start and end of waveform. Default + to FLR default of 10000 s. tdi : str (optional) The TDI observables to calculate. Accepts 'XYZ', 'AET', or 'AE'. @@ -813,6 +817,10 @@ def get_tdi(self, hp, hc, lamb, beta, tdi='XYZ', tdi_orbit_kwargs=None, use_gpu= # set use_gpu if use_gpu is None: use_gpu = self.use_gpu + + # update t0 if provided + if t0 is not None: + self.t0 = t0 # generate the Doppler time series links = self.project_wf(hp, hc, lamb, beta, t0=self.t0, use_gpu=use_gpu) From 9dd0b839a63ea485d18c04f3ee01682b0a0cfb96 Mon Sep 17 00:00:00 2001 From: acorreia61201 Date: Mon, 29 Apr 2024 19:30:46 +0000 Subject: [PATCH 03/24] Fixed sampling frequency bug in pyResponseTDI class initialization; fixed bug where orbit file was truncated by default init params --- pycbc/detector.py | 12 +++--------- 1 file changed, 3 insertions(+), 9 deletions(-) diff --git a/pycbc/detector.py b/pycbc/detector.py index 68e75e1e6a6..7f40bd8002d 100644 --- a/pycbc/detector.py +++ b/pycbc/detector.py @@ -664,12 +664,7 @@ def __init__(self, orbit_kwargs, use_gpu=False): ---------- orbit_kwargs : dict List of keywords for generating LISA orbital data. Passed into pyResponseTDI - instance. Keywords and defaults as set by lisaresponse are: - - 'orbit_module' = None - 'order' = 0 - 'max_t_orbits' = 3.15576e7 - 'orbit_file' = None + instance. Either an orbit file or orbit module is required for projection. If both a file and module are provided, the file will take priority. @@ -679,8 +674,7 @@ def __init__(self, orbit_kwargs, use_gpu=False): Specify whether to run class on GPU support. Default False. """ # intialize orbit kwargs to FLR defaults - self.orbit_kwargs = {'orbit_module': None, 'order': 0, 'max_t_orbits': 3.15576e7, - 'orbit_file': None} + self.orbit_kwargs = {'orbit_module': None, 'orbit_file': None} # get kwargs from class input for key, val in orbit_kwargs.items(): @@ -760,7 +754,7 @@ def project_wf(self, hp, hc, lamb, beta, t0=None, use_gpu=None): wf = cupy.asarray(wf) # initialize the FLR class - response_init = pyResponseTDI(dt, n, self.orbit_kwargs, use_gpu=use_gpu) + response_init = pyResponseTDI(1/dt, n, self.orbit_kwargs, use_gpu=use_gpu) # cache response init self.dt = dt From 40e83d356bf3414a534bfe5a598ea32bdbd1a97d Mon Sep 17 00:00:00 2001 From: acorreia61201 Date: Fri, 31 May 2024 16:28:40 +0000 Subject: [PATCH 04/24] TDI configuration now actually updates when providing tdi input in get_tdi() --- pycbc/detector.py | 153 +++++++++++++++++++++++++++------------------- 1 file changed, 89 insertions(+), 64 deletions(-) diff --git a/pycbc/detector.py b/pycbc/detector.py index 7f40bd8002d..ecd38e7f970 100644 --- a/pycbc/detector.py +++ b/pycbc/detector.py @@ -654,51 +654,41 @@ def overhead_antenna_pattern(right_ascension, declination, polarization): except ImportError: cupy = None +try: + from lisatools.detector import ESAOrbits +except ImportError: + ESAOrbits = None + class LISA_detector(object): """ LISA-like GW detector. Applies detector response from FastLISAResponse. """ - def __init__(self, orbit_kwargs, use_gpu=False): + def __init__(self, orbits=ESAOrbits(), use_gpu=False): """ Parameters ---------- - orbit_kwargs : dict - List of keywords for generating LISA orbital data. Passed into pyResponseTDI - instance. - - Either an orbit file or orbit module is required for projection. - If both a file and module are provided, the file will take priority. - If neither are provided, a ValueError is raised. + orbits: lisatools.detector.Orbits, optional + Orbital information to pass into pyResponseTDI. Accepts LISA Analysis Tools + format. Default ESAOrbits. use_gpu : bool (optional) Specify whether to run class on GPU support. Default False. """ - # intialize orbit kwargs to FLR defaults - self.orbit_kwargs = {'orbit_module': None, 'orbit_file': None} - - # get kwargs from class input - for key, val in orbit_kwargs.items(): - self.orbit_kwargs[key] = val - - # handle if both/neither module and file are specified - if self.orbit_kwargs['orbit_module'] is not None: - if self.orbit_kwargs['orbit_file'] is not None: - # FLR defaults to the module; override this and default to the file - self.orbit_kwargs['orbit_module'] = None - else: - if self.orbit_kwargs['orbit_file'] is None: - raise ValueError('Either orbit file or orbit module required to initialize LISA detector') + # intialize orbit information + if orbits is None: + raise ImportError('LISAanalysistools required for inputting orbital data') + self.orbits = orbits # cache the FLR instance along with dt and n self.dt = None self.n = None self.response_init = None - self.t0 = 10000 - + self.t0 = 10000. + # initialize whether to use gpu; FLR has handles if this cannot be done self.use_gpu = use_gpu - def project_wf(self, hp, hc, lamb, beta, t0=None, use_gpu=None): + def project_wf(self, hp, hc, lamb, beta, t0=10000., use_gpu=None): """ Project a radiation frame (SSB) waveform to the LISA constellation. This *does not* perform TDI; this only produces a detector frame @@ -719,8 +709,8 @@ def project_wf(self, hp, hc, lamb, beta, t0=None, use_gpu=None): The ecleptic longitude in the SSB frame. t0 : float (optional) - Number of seconds to omit from start and end of waveform. Default - to FLR default of 10000 s. + Number of seconds to omit from start and end of waveform. Defaults + to FastLISAResponse preset of 10000 s. use_gpu : bool (optional) Flag whether to use GPU support. Default to class input. @@ -734,9 +724,16 @@ def project_wf(self, hp, hc, lamb, beta, t0=None, use_gpu=None): """ from fastlisaresponse import pyResponseTDI + # get dt and n from waveform data dt = hp.delta_t # timestep (s) n = len(hp) # number of points + # configure the orbits file according to wf times + ### FIXME: This currently doesn't work. There seems to be a bug in + ### LISAanalysistools that breaks the code when specifying a time array. + # self.sample_times = hp.sample_times.numpy() + # self.orbits.configure(t_arr=self.sample_times) + # format wf to FLR standard hp + i*hc as ndarray hp = hp.numpy() hc = hc.numpy() @@ -752,24 +749,33 @@ def project_wf(self, hp, hc, lamb, beta, t0=None, use_gpu=None): raise ImportError('CuPy not imported but is required for GPU usage') else: wf = cupy.asarray(wf) - - # initialize the FLR class - response_init = pyResponseTDI(1/dt, n, self.orbit_kwargs, use_gpu=use_gpu) - # cache response init + if self.response_init is None: + # initialize the class + ### TDI is set to '2nd generation' here to match default value in get_tdi() + # see FIXME in get_tdi() + response_init = pyResponseTDI(1/dt, n, orbits=self.orbits, use_gpu=use_gpu, tdi='2nd generation') + self.response_init = response_init + else: + # update params in the initialized class + self.response_init.sampling_frequency = 1/dt + self.response_init.num_pts = n + self.response_init.orbits = self.orbits + self.response_init.use_gpu = use_gpu + + # cache n, dt, response class self.dt = dt self.n = n - self.response_init = response_init - if t0 is not None: - self.t0 = t0 + self.t0 = t0 # project the signal - response_init.get_projections(wf, lamb, beta, self.t0) - wf_proj = response_init.y_gw - + self.response_init.get_projections(wf, lamb, beta, t0) + wf_proj = self.response_init.y_gw + return wf_proj - def get_tdi(self, hp, hc, lamb, beta, t0=None, tdi='XYZ', tdi_orbit_kwargs=None, use_gpu=None): + def get_tdi(self, hp, hc, lamb, beta, tdi='2nd generation', tdi_chan='XYZ', tdi_orbits=None, + use_gpu=None, remove_garbage=True): """ Calculate the TDI variables given a waveform and LISA orbit data. @@ -781,27 +787,31 @@ def get_tdi(self, hp, hc, lamb, beta, t0=None, tdi='XYZ', tdi_orbit_kwargs=None, hc : pycbc.types.TimeSeries The cross polarization of the GW in the radiation frame. - lambda : float + lamb : float The ecleptic latitude in the SSB frame. beta : float The ecleptic longitude in the SSB frame. - - t0 : float (optional) - Number of seconds to omit from start and end of waveform. Default - to FLR default of 10000 s. - tdi : str (optional) + tdi : string (optional) + TDI channel configuration. Accepts "1st generation" or "2nd generation" + as inputs. Default "2nd generation". + + tdi_chan : str (optional) The TDI observables to calculate. Accepts 'XYZ', 'AET', or 'AE'. Default 'XYZ'. - tdi_orbit_kwargs : dict (optional) - Orbit keywords specifically for TDI generation. If None, use - orbit_kwargs provided on class initialization. Default None. + tdi_orbits : lisatools.detector.Orbits (optional) + Orbit keywords specifically for TDI generation. Default to class input. use_gpu : bool (optional) Flag whether to use GPU support. Default to class input. + remove_garbage : bool (optional) + Flag whether to excise data from start and end of TDI projections. + This will remove erroneous ringing at start and end by excising + time length t0. Default True. + Returns ------- dict @@ -811,31 +821,46 @@ def get_tdi(self, hp, hc, lamb, beta, t0=None, tdi='XYZ', tdi_orbit_kwargs=None, # set use_gpu if use_gpu is None: use_gpu = self.use_gpu - - # update t0 if provided - if t0 is not None: - self.t0 = t0 - - # generate the Doppler time series - links = self.project_wf(hp, hc, lamb, beta, t0=self.t0, use_gpu=use_gpu) - # if TDI kwargs are provided, update the response_init - if type(tdi_orbit_kwargs) == dict: - self.response_init.tdi_orbit_kwargs = tdi_orbit_kwargs + tdi_args = {} + + # set TDI configuration + if tdi in ['1st generation', '2nd generation']: + original_tdi = self.response_init.tdi + self.response_init.tdi = tdi + else: + raise ValueError('TDI configuration must be either "1st generation" or "2nd generation"') # set TDI channels - if tdi in ['XYZ', 'AET', 'AE']: - self.response_init.tdi_chan = tdi + if tdi_chan in ['XYZ', 'AET', 'AE']: + self.response_init.tdi_chan = tdi_chan else: raise ValueError('TDI channels must be one of: XYZ, AET, AE') + # if TDI orbit class is provided, update the response_init + # tdi_orbits are set to class input automatically by FLR otherwise + if tdi_orbits is not None: + ### FIXME: bug in LISAanalysistools when specifying time array + # tdi_orbits.configure(t_arr=self.sample_times) + self.response_init.tdi_orbits = tdi_orbits + + ### the TDI config is not automatically updated by setting tdi attributes + ### for now, update the TDI config manually if different from default + ### FIXME: find a better, faster way to do this + if self.response_init.tdi != original_tdi: + self.response_init._init_tdi_delays() + + # generate the Doppler time series + links = self.project_wf(hp, hc, lamb, beta, t0=self.t0, use_gpu=use_gpu) + + # generate the TDI channels tdi_obs = self.response_init.get_tdi_delays() - + # convert to dict tdi_dict = {} - for i in range(len(tdi)): - tdi_dict[tdi[i]] = tdi_obs[i] - + for i in range(len(tdi_chan)): + tdi_dict[tdi_chan[i]] = tdi_obs[i] + return tdi_dict From 9ca5a865e13839caecf8726d296d4c0cfd40cd56 Mon Sep 17 00:00:00 2001 From: acorreia61201 Date: Mon, 3 Jun 2024 18:57:16 +0000 Subject: [PATCH 05/24] Added reference_time to methods as kwarg; revised methods to be more consistent with LIGO detector class --- pycbc/detector.py | 25 ++++++++++++++++++------- 1 file changed, 18 insertions(+), 7 deletions(-) diff --git a/pycbc/detector.py b/pycbc/detector.py index ecd38e7f970..bc1fbb70ec9 100644 --- a/pycbc/detector.py +++ b/pycbc/detector.py @@ -688,7 +688,7 @@ def __init__(self, orbits=ESAOrbits(), use_gpu=False): # initialize whether to use gpu; FLR has handles if this cannot be done self.use_gpu = use_gpu - def project_wf(self, hp, hc, lamb, beta, t0=10000., use_gpu=None): + def get_links(self, hp, hc, lamb, beta, t0=10000., use_gpu=None, reference_time=None): """ Project a radiation frame (SSB) waveform to the LISA constellation. This *does not* perform TDI; this only produces a detector frame @@ -716,6 +716,10 @@ def project_wf(self, hp, hc, lamb, beta, t0=10000., use_gpu=None): Flag whether to use GPU support. Default to class input. CuPy is required if use_gpu == True; an ImportError will be raised if CuPy could not be imported. + + reference_time : float (optional) + The GPS start time of the GW signal in the SSB frame. Default to + value in input signals hp and hc. Returns ------- @@ -729,9 +733,12 @@ def project_wf(self, hp, hc, lamb, beta, t0=10000., use_gpu=None): n = len(hp) # number of points # configure the orbits file according to wf times + if reference_time is not None: + hp.start_time = reference_time + self.sample_times = hp.sample_times.numpy() + ### FIXME: This currently doesn't work. There seems to be a bug in ### LISAanalysistools that breaks the code when specifying a time array. - # self.sample_times = hp.sample_times.numpy() # self.orbits.configure(t_arr=self.sample_times) # format wf to FLR standard hp + i*hc as ndarray @@ -753,7 +760,7 @@ def project_wf(self, hp, hc, lamb, beta, t0=10000., use_gpu=None): if self.response_init is None: # initialize the class ### TDI is set to '2nd generation' here to match default value in get_tdi() - # see FIXME in get_tdi() + ### see FIXME in get_tdi() response_init = pyResponseTDI(1/dt, n, orbits=self.orbits, use_gpu=use_gpu, tdi='2nd generation') self.response_init = response_init else: @@ -774,10 +781,10 @@ def project_wf(self, hp, hc, lamb, beta, t0=10000., use_gpu=None): return wf_proj - def get_tdi(self, hp, hc, lamb, beta, tdi='2nd generation', tdi_chan='XYZ', tdi_orbits=None, + def project_wave(self, hp, hc, lamb, beta, reference_time=None, tdi='2nd generation', tdi_chan='XYZ', tdi_orbits=None, use_gpu=None, remove_garbage=True): """ - Calculate the TDI variables given a waveform and LISA orbit data. + Evaluate the TDI observables Parameters ---------- @@ -792,7 +799,11 @@ def get_tdi(self, hp, hc, lamb, beta, tdi='2nd generation', tdi_chan='XYZ', tdi_ beta : float The ecleptic longitude in the SSB frame. - + + reference_time : float (optional) + The GPS start time of the GW signal in the SSB frame. Default to + value in input signals hp and hc. + tdi : string (optional) TDI channel configuration. Accepts "1st generation" or "2nd generation" as inputs. Default "2nd generation". @@ -851,7 +862,7 @@ def get_tdi(self, hp, hc, lamb, beta, tdi='2nd generation', tdi_chan='XYZ', tdi_ self.response_init._init_tdi_delays() # generate the Doppler time series - links = self.project_wf(hp, hc, lamb, beta, t0=self.t0, use_gpu=use_gpu) + links = self.get_links(hp, hc, lamb, beta, t0=self.t0, use_gpu=use_gpu, reference_time=reference_time) # generate the TDI channels tdi_obs = self.response_init.get_tdi_delays() From bfd83b50ef051312e2a44505df6793cbc34247c5 Mon Sep 17 00:00:00 2001 From: acorreia61201 Date: Tue, 4 Jun 2024 19:58:19 +0000 Subject: [PATCH 06/24] Revised init signature to be more similar to LIGO detector; rebase to most recent PyCBC patch (6/4/24) --- pycbc/detector.py | 61 ++++++++++++++++++++++++++++++----------------- 1 file changed, 39 insertions(+), 22 deletions(-) diff --git a/pycbc/detector.py b/pycbc/detector.py index bc1fbb70ec9..d513dbd1a06 100644 --- a/pycbc/detector.py +++ b/pycbc/detector.py @@ -663,21 +663,39 @@ class LISA_detector(object): """ LISA-like GW detector. Applies detector response from FastLISAResponse. """ - def __init__(self, orbits=ESAOrbits(), use_gpu=False): + def __init__(self, detector='LISA', reference_time=None, orbits=ESAOrbits(), use_gpu=False): """ Parameters ---------- - orbits: lisatools.detector.Orbits, optional - Orbital information to pass into pyResponseTDI. Accepts LISA Analysis Tools - format. Default ESAOrbits. + detector: str (optional) + String specifying space-borne detector. Currently only accepts 'LISA', + which is the default setting. + + reference_time: float (optional) + The GPS start time of signal in the SSB frame. Default to start time of + orbits input. + + orbits: lisatools.detector.Orbits (optional) + Orbital information to pass into pyResponseTDI. Default + lisatools.detector.ESAOrbits. use_gpu : bool (optional) - Specify whether to run class on GPU support. Default False. + Specify whether to run class on GPU support via CuPy. Default False. """ + # initialize detector; for now we only accept LISA + assert (detector=='LISA'), 'Currently only supports LISA detector' + # intialize orbit information if orbits is None: raise ImportError('LISAanalysistools required for inputting orbital data') + orbits.configure(linear_interp_setup=True) self.orbits = orbits + + # specify and cache the start time and orbital time series + if reference_time is None: + reference_time = self.orbits.t_base[0] + self.ref_time = reference_time + self.sample_times = None # cache the FLR instance along with dt and n self.dt = None @@ -719,7 +737,7 @@ def get_links(self, hp, hc, lamb, beta, t0=10000., use_gpu=None, reference_time= reference_time : float (optional) The GPS start time of the GW signal in the SSB frame. Default to - value in input signals hp and hc. + input on class initialization. Returns ------- @@ -732,11 +750,18 @@ def get_links(self, hp, hc, lamb, beta, t0=10000., use_gpu=None, reference_time= dt = hp.delta_t # timestep (s) n = len(hp) # number of points - # configure the orbits file according to wf times + # cache n, dt + self.dt = dt + self.n = n + self.t0 = t0 + + # set waveform start time and cache the resultant time series if reference_time is not None: - hp.start_time = reference_time + self.ref_time = reference_time + hp.start_time = self.ref_time self.sample_times = hp.sample_times.numpy() - + + # rescale the orbital time series to match waveform ### FIXME: This currently doesn't work. There seems to be a bug in ### LISAanalysistools that breaks the code when specifying a time array. # self.orbits.configure(t_arr=self.sample_times) @@ -759,8 +784,8 @@ def get_links(self, hp, hc, lamb, beta, t0=10000., use_gpu=None, reference_time= if self.response_init is None: # initialize the class - ### TDI is set to '2nd generation' here to match default value in get_tdi() - ### see FIXME in get_tdi() + ### TDI is set to '2nd generation' here to match default value in project_wave() + ### see FIXME in project_wave() response_init = pyResponseTDI(1/dt, n, orbits=self.orbits, use_gpu=use_gpu, tdi='2nd generation') self.response_init = response_init else: @@ -770,11 +795,6 @@ def get_links(self, hp, hc, lamb, beta, t0=10000., use_gpu=None, reference_time= self.response_init.orbits = self.orbits self.response_init.use_gpu = use_gpu - # cache n, dt, response class - self.dt = dt - self.n = n - self.t0 = t0 - # project the signal self.response_init.get_projections(wf, lamb, beta, t0) wf_proj = self.response_init.y_gw @@ -828,12 +848,12 @@ def project_wave(self, hp, hc, lamb, beta, reference_time=None, tdi='2nd generat dict The TDI observables keyed by their corresponding TDI channel. """ - # set use_gpu if use_gpu is None: use_gpu = self.use_gpu - - tdi_args = {} + + # generate the Doppler time series + links = self.get_links(hp, hc, lamb, beta, t0=self.t0, use_gpu=use_gpu, reference_time=reference_time) # set TDI configuration if tdi in ['1st generation', '2nd generation']: @@ -861,9 +881,6 @@ def project_wave(self, hp, hc, lamb, beta, reference_time=None, tdi='2nd generat if self.response_init.tdi != original_tdi: self.response_init._init_tdi_delays() - # generate the Doppler time series - links = self.get_links(hp, hc, lamb, beta, t0=self.t0, use_gpu=use_gpu, reference_time=reference_time) - # generate the TDI channels tdi_obs = self.response_init.get_tdi_delays() From c9c19ea148236ad61cc38451b19a48842acb7ab6 Mon Sep 17 00:00:00 2001 From: acorreia61201 Date: Thu, 6 Jun 2024 18:49:26 +0000 Subject: [PATCH 07/24] Added polarization to project_wave signature; allow for shifting of orbit file start time --- pycbc/detector.py | 121 +++++++++++++++++++++++++++++++--------------- 1 file changed, 82 insertions(+), 39 deletions(-) diff --git a/pycbc/detector.py b/pycbc/detector.py index d513dbd1a06..0c3ad703b22 100644 --- a/pycbc/detector.py +++ b/pycbc/detector.py @@ -705,12 +705,41 @@ def __init__(self, detector='LISA', reference_time=None, orbits=ESAOrbits(), use # initialize whether to use gpu; FLR has handles if this cannot be done self.use_gpu = use_gpu + + def source_to_ssb(self, hp, hc, polarization): + """ + Project a source-frame signal to SSB coordinates. + + Parameters + ---------- + hp : pycbc.types.TimeSeries + The plus polarization of the GW in the radiation frame. + + hc : pycbc.types.TimeSeries + The cross polarization of the GW in the radiation frame. + + polarization : float + The polarization angle of the GW in radians. + + Returns + ------- + (pycbc.types.TimeSeries, pycbc.types.TimeSeries) + The plus and cross polarizations of the GW in the SSB frame. + """ + cphi = cos(2*polarization) + sphi = sin(2*polarization) + + hp_ssb = hp*cphi - hc*sphi + hc_ssb = hp*sphi + hc*cphi + + return hp_ssb, hc_ssb - def get_links(self, hp, hc, lamb, beta, t0=10000., use_gpu=None, reference_time=None): + def get_links(self, hp, hc, lamb, beta, polarization=0, reference_time=None, + t0=10000., use_gpu=None, tdi=2): """ - Project a radiation frame (SSB) waveform to the LISA constellation. - This *does not* perform TDI; this only produces a detector frame - representation of the SSB frame signal. + Project a radiation frame waveform to the LISA constellation. + This does not perform TDI; this only produces a detector frame + representation of the signal. Parameters ---------- @@ -720,11 +749,20 @@ def get_links(self, hp, hc, lamb, beta, t0=10000., use_gpu=None, reference_time= hc : pycbc.types.TimeSeries The cross polarization of the GW in the radiation frame. - lambda : float - The ecleptic latitude in the SSB frame. + lamb : float + The ecleptic latitude of the source in the SSB frame. beta : float - The ecleptic longitude in the SSB frame. + The ecleptic longitude of the source in the SSB frame. + + polarization : float (optional) + The polarization angle of the GW in radians. If polarization is + 0, the hp and hc inputs are treated as SSB frame projections. + Default 0. + + reference_time : float (optional) + The GPS start time of the GW signal in the SSB frame. Default to + input on class initialization. t0 : float (optional) Number of seconds to omit from start and end of waveform. Defaults @@ -735,9 +773,9 @@ def get_links(self, hp, hc, lamb, beta, t0=10000., use_gpu=None, reference_time= CuPy is required if use_gpu == True; an ImportError will be raised if CuPy could not be imported. - reference_time : float (optional) - The GPS start time of the GW signal in the SSB frame. Default to - input on class initialization. + tdi : int (optional) + TDI channel configuration. Accepts 1 for 1st generation TDI or + 2 for 2nd generation TDI. Default 2. Returns ------- @@ -759,12 +797,15 @@ def get_links(self, hp, hc, lamb, beta, t0=10000., use_gpu=None, reference_time= if reference_time is not None: self.ref_time = reference_time hp.start_time = self.ref_time + hc.start_time = self.ref_time self.sample_times = hp.sample_times.numpy() # rescale the orbital time series to match waveform - ### FIXME: This currently doesn't work. There seems to be a bug in - ### LISAanalysistools that breaks the code when specifying a time array. - # self.orbits.configure(t_arr=self.sample_times) + self.orbits.configure(t_arr=self.sample_times) + + # convert to SSB frame (if pol = 0, the signal is already in SSB frame) + if polarization != 0.: + hp, hc = self.source_to_ssb(hp, hc, polarization) # format wf to FLR standard hp + i*hc as ndarray hp = hp.numpy() @@ -781,12 +822,18 @@ def get_links(self, hp, hc, lamb, beta, t0=10000., use_gpu=None, reference_time= raise ImportError('CuPy not imported but is required for GPU usage') else: wf = cupy.asarray(wf) + + # set TDI configuration + if tdi == 1: + tdi_opt = '1st generation' + if tdi == 2: + tdi_opt = '2nd generation' if self.response_init is None: # initialize the class - ### TDI is set to '2nd generation' here to match default value in project_wave() - ### see FIXME in project_wave() - response_init = pyResponseTDI(1/dt, n, orbits=self.orbits, use_gpu=use_gpu, tdi='2nd generation') + print(self.orbits.pycppdetector_args) + response_init = pyResponseTDI(1/dt, n, orbits=self.orbits, + use_gpu=use_gpu, tdi=tdi_opt) self.response_init = response_init else: # update params in the initialized class @@ -794,6 +841,10 @@ def get_links(self, hp, hc, lamb, beta, t0=10000., use_gpu=None, reference_time= self.response_init.num_pts = n self.response_init.orbits = self.orbits self.response_init.use_gpu = use_gpu + if tdi_opt != self.response_init.tdi: + # re-initialize TDI in existing response_init class + self.response_init.tdi = tdi_opt + self.response_init._init_TDI_delays() # project the signal self.response_init.get_projections(wf, lamb, beta, t0) @@ -801,10 +852,10 @@ def get_links(self, hp, hc, lamb, beta, t0=10000., use_gpu=None, reference_time= return wf_proj - def project_wave(self, hp, hc, lamb, beta, reference_time=None, tdi='2nd generation', tdi_chan='XYZ', tdi_orbits=None, - use_gpu=None, remove_garbage=True): + def project_wave(self, hp, hc, lamb, beta, polarization=0, reference_time=None, tdi=2, + tdi_chan='AET', tdi_orbits=None, use_gpu=None, remove_garbage=True): """ - Evaluate the TDI observables + Evaluate the TDI observables. Parameters ---------- @@ -820,17 +871,22 @@ def project_wave(self, hp, hc, lamb, beta, reference_time=None, tdi='2nd generat beta : float The ecleptic longitude in the SSB frame. + polarization : float (optional) + The polarization angle of the GW in radians. If polarization is + 0, the hp and hc inputs are treated as SSB frame projections. + Default 0. + reference_time : float (optional) The GPS start time of the GW signal in the SSB frame. Default to value in input signals hp and hc. - tdi : string (optional) - TDI channel configuration. Accepts "1st generation" or "2nd generation" - as inputs. Default "2nd generation". + 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'. - Default 'XYZ'. + Default 'AET'. tdi_orbits : lisatools.detector.Orbits (optional) Orbit keywords specifically for TDI generation. Default to class input. @@ -853,14 +909,8 @@ def project_wave(self, hp, hc, lamb, beta, reference_time=None, tdi='2nd generat use_gpu = self.use_gpu # generate the Doppler time series - links = self.get_links(hp, hc, lamb, beta, t0=self.t0, use_gpu=use_gpu, reference_time=reference_time) - - # set TDI configuration - if tdi in ['1st generation', '2nd generation']: - original_tdi = self.response_init.tdi - self.response_init.tdi = tdi - else: - raise ValueError('TDI configuration must be either "1st generation" or "2nd generation"') + self.get_links(hp, hc, lamb, beta, polarization=polarization, reference_time=reference_time, + t0=self.t0, use_gpu=use_gpu, tdi=tdi) # set TDI channels if tdi_chan in ['XYZ', 'AET', 'AE']: @@ -871,16 +921,9 @@ def project_wave(self, hp, hc, lamb, beta, reference_time=None, tdi='2nd generat # if TDI orbit class is provided, update the response_init # tdi_orbits are set to class input automatically by FLR otherwise if tdi_orbits is not None: - ### FIXME: bug in LISAanalysistools when specifying time array - # tdi_orbits.configure(t_arr=self.sample_times) + tdi_orbits.configure(t_arr=self.sample_times) self.response_init.tdi_orbits = tdi_orbits - ### the TDI config is not automatically updated by setting tdi attributes - ### for now, update the TDI config manually if different from default - ### FIXME: find a better, faster way to do this - if self.response_init.tdi != original_tdi: - self.response_init._init_tdi_delays() - # generate the TDI channels tdi_obs = self.response_init.get_tdi_delays() From a2fb6a440616a39d78efac392f7d78a483d352aa Mon Sep 17 00:00:00 2001 From: acorreia61201 Date: Sat, 15 Jun 2024 15:46:39 +0000 Subject: [PATCH 08/24] Add functionality to remove or zero edge data --- pycbc/detector.py | 76 +++++++++++++++++++++++++++++++---------------- 1 file changed, 51 insertions(+), 25 deletions(-) diff --git a/pycbc/detector.py b/pycbc/detector.py index 0c3ad703b22..4b8cec9386f 100644 --- a/pycbc/detector.py +++ b/pycbc/detector.py @@ -700,6 +700,7 @@ def __init__(self, detector='LISA', reference_time=None, orbits=ESAOrbits(), use # cache the FLR instance along with dt and n self.dt = None self.n = None + self.pad_idx = None self.response_init = None self.t0 = 10000. @@ -708,7 +709,7 @@ def __init__(self, detector='LISA', reference_time=None, orbits=ESAOrbits(), use def source_to_ssb(self, hp, hc, polarization): """ - Project a source-frame signal to SSB coordinates. + Project a source-frame signal to SSB. Parameters ---------- @@ -735,7 +736,7 @@ def source_to_ssb(self, hp, hc, polarization): return hp_ssb, hc_ssb def get_links(self, hp, hc, lamb, beta, polarization=0, reference_time=None, - t0=10000., use_gpu=None, tdi=2): + use_gpu=None, tdi=2): """ Project a radiation frame waveform to the LISA constellation. This does not perform TDI; this only produces a detector frame @@ -764,13 +765,9 @@ def get_links(self, hp, hc, lamb, beta, polarization=0, reference_time=None, The GPS start time of the GW signal in the SSB frame. Default to input on class initialization. - t0 : float (optional) - Number of seconds to omit from start and end of waveform. Defaults - to FastLISAResponse preset of 10000 s. - use_gpu : bool (optional) Flag whether to use GPU support. Default to class input. - CuPy is required if use_gpu == True; an ImportError will be raised + CuPy is required if use_gpu is True; an ImportError will be raised if CuPy could not be imported. tdi : int (optional) @@ -783,17 +780,17 @@ def get_links(self, hp, hc, lamb, beta, polarization=0, reference_time=None, The waveform projected to the LISA laser links. """ from fastlisaresponse import pyResponseTDI - + # get dt and n from waveform data dt = hp.delta_t # timestep (s) n = len(hp) # number of points - # cache n, dt + # cache n, dt, pad length self.dt = dt self.n = n - self.t0 = t0 + self.pad_idx = int(self.t0/self.dt) - # set waveform start time and cache the resultant time series + # set waveform start time and cache time series if reference_time is not None: self.ref_time = reference_time hp.start_time = self.ref_time @@ -801,6 +798,8 @@ def get_links(self, hp, hc, lamb, beta, polarization=0, reference_time=None, self.sample_times = hp.sample_times.numpy() # rescale the orbital time series to match waveform + ### requires bugfix in LISAanalysistools to function + ### awaiting response from dev team for formal release self.orbits.configure(t_arr=self.sample_times) # convert to SSB frame (if pol = 0, the signal is already in SSB frame) @@ -823,15 +822,16 @@ def get_links(self, hp, hc, lamb, beta, polarization=0, reference_time=None, else: wf = cupy.asarray(wf) - # set TDI configuration + # set TDI configuration (let FLR handle if not 1 or 2) if tdi == 1: tdi_opt = '1st generation' - if tdi == 2: + elif tdi == 2: tdi_opt = '2nd generation' + else: + tdi_opt = tdi if self.response_init is None: # initialize the class - print(self.orbits.pycppdetector_args) response_init = pyResponseTDI(1/dt, n, orbits=self.orbits, use_gpu=use_gpu, tdi=tdi_opt) self.response_init = response_init @@ -841,13 +841,14 @@ def get_links(self, hp, hc, lamb, beta, polarization=0, reference_time=None, self.response_init.num_pts = n self.response_init.orbits = self.orbits self.response_init.use_gpu = use_gpu + if tdi_opt != self.response_init.tdi: - # re-initialize TDI in existing response_init class + # update TDI in existing response_init class self.response_init.tdi = tdi_opt self.response_init._init_TDI_delays() # project the signal - self.response_init.get_projections(wf, lamb, beta, t0) + self.response_init.get_projections(wf, lamb, beta, t0=self.t0) wf_proj = self.response_init.y_gw return wf_proj @@ -889,28 +890,32 @@ def project_wave(self, hp, hc, lamb, beta, polarization=0, reference_time=None, Default 'AET'. tdi_orbits : lisatools.detector.Orbits (optional) - Orbit keywords specifically for TDI generation. Default to class input. + Orbit keywords specifically for TDI generation. Default to class + input. use_gpu : bool (optional) Flag whether to use GPU support. Default to class input. - - remove_garbage : bool (optional) - Flag whether to excise data from start and end of TDI projections. - This will remove erroneous ringing at start and end by excising - time length t0. Default True. + + remove_garbage : bool, str (optional) + Flag whether to remove gaps in TDI from start and end. If True, + edge data will be cut from TDI channels. If 'zero', edge data + will be zeroed. If False, TDI channels will not be modified. + Default True. Returns ------- dict The TDI observables keyed by their corresponding TDI channel. """ + start_time = time.time() + # set use_gpu if use_gpu is None: use_gpu = self.use_gpu # generate the Doppler time series self.get_links(hp, hc, lamb, beta, polarization=polarization, reference_time=reference_time, - t0=self.t0, use_gpu=use_gpu, tdi=tdi) + use_gpu=use_gpu, tdi=tdi) # set TDI channels if tdi_chan in ['XYZ', 'AET', 'AE']: @@ -921,16 +926,37 @@ def project_wave(self, hp, hc, lamb, beta, polarization=0, reference_time=None, # if TDI orbit class is provided, update the response_init # tdi_orbits are set to class input automatically by FLR otherwise if tdi_orbits is not None: + ### requires bugfix in LISAanalysistools to function + ### awaiting response from dev team for formal release tdi_orbits.configure(t_arr=self.sample_times) self.response_init.tdi_orbits = tdi_orbits # generate the TDI channels tdi_obs = self.response_init.get_tdi_delays() - - # convert to dict + + # processing tdi_dict = {} + slc = slice(self.pad_idx, -self.pad_idx) + for i in range(len(tdi_chan)): + # add to dict tdi_dict[tdi_chan[i]] = tdi_obs[i] + + # treat start and end gaps (if remove_garbage is not False) + if remove_garbage: + if remove_garbage == 'zero': + # zero the edge data + tdi_dict[tdi_chan[i]][self.pad_idx:] = 0 + tdi_dict[tdi_chan[i]][:-self.pad_idx] = 0 + elif type(remove_garbage) == bool: + # cut the edge data + tdi_dict[tdi_chan[i]] = tdi_dict[tdi_chan[i]][slc] + if i == 0: + # cut the sample times (only do once) + self.sample_times = self.sample_times[slc] + else: + # raise error + raise ValueError('remove_garbage arg must be a bool or "zero"') return tdi_dict From 5fe2702f9f66c23046b32bb5a28d4051fbebba95 Mon Sep 17 00:00:00 2001 From: acorreia61201 Date: Thu, 11 Jul 2024 20:01:29 +0000 Subject: [PATCH 09/24] edited polarization rotation function; final touches for first draft --- pycbc/detector.py | 110 +++++++++++++++++++++++++++++++--------------- 1 file changed, 74 insertions(+), 36 deletions(-) diff --git a/pycbc/detector.py b/pycbc/detector.py index 4b8cec9386f..6b6ef7e06f4 100644 --- a/pycbc/detector.py +++ b/pycbc/detector.py @@ -648,7 +648,6 @@ def overhead_antenna_pattern(right_ascension, declination, polarization): """ LISA class """ - try: import cupy except ImportError: @@ -694,7 +693,7 @@ def __init__(self, detector='LISA', reference_time=None, orbits=ESAOrbits(), use # specify and cache the start time and orbital time series if reference_time is None: reference_time = self.orbits.t_base[0] - self.ref_time = reference_time + self.com_ref_time = reference_time self.sample_times = None # cache the FLR instance along with dt and n @@ -707,25 +706,28 @@ def __init__(self, detector='LISA', reference_time=None, orbits=ESAOrbits(), use # initialize whether to use gpu; FLR has handles if this cannot be done self.use_gpu = use_gpu - def source_to_ssb(self, hp, hc, polarization): + def apply_polarization(self, hp, hc, polarization): """ - Project a source-frame signal to SSB. + Apply polarization rotation matrix. + + .. math:: + \bmatrix{} Parameters ---------- - hp : pycbc.types.TimeSeries - The plus polarization of the GW in the radiation frame. + hp : array + The plus polarization of the GW. - hc : pycbc.types.TimeSeries - The cross polarization of the GW in the radiation frame. + hc : array + The cross polarization of the GW. polarization : float The polarization angle of the GW in radians. Returns ------- - (pycbc.types.TimeSeries, pycbc.types.TimeSeries) - The plus and cross polarizations of the GW in the SSB frame. + (array, array) + The plus and cross polarizations of the GW rotated by the polarization angle. """ cphi = cos(2*polarization) sphi = sin(2*polarization) @@ -734,6 +736,44 @@ def source_to_ssb(self, hp, hc, polarization): hc_ssb = hp*sphi + hc*cphi return hp_ssb, hc_ssb + + def time_delay_from_ssb(self, orbits, reference_time): + """ + Calculate the time delay from the SSB frame to the LISA COM frame. + + Parameters + ---------- + orbits: lisatools.detector.Orbits + The orbital information of the satellites. + + reference_time: float + The time in seconds in the SSB frame. + + Returns + ------- + float + The additive time delay factor between the SSB and LISA COM frame + at the given time. + """ + C_SI = constants.c.value + + # configure orbits if not already + try: + orbits.t + except ValueError: + orbits.configure(linear_interp_setup=True) + + # get positions of spacecraft at reference time + sc_x = [] + for i in [1, 2, 3]: + sc_x.append(orbits.get_pos(reference_time, i)) + + # average sc positions to get COM position + com_vec = sum(sc_x)/len(sc_x) + com_dist = sum(com_vec*com_vec)**0.5 + + # time delay from SSB to LISA COM is distance/c + return com_dist/C_SI def get_links(self, hp, hc, lamb, beta, polarization=0, reference_time=None, use_gpu=None, tdi=2): @@ -757,9 +797,7 @@ def get_links(self, hp, hc, lamb, beta, polarization=0, reference_time=None, The ecleptic longitude of the source in the SSB frame. polarization : float (optional) - The polarization angle of the GW in radians. If polarization is - 0, the hp and hc inputs are treated as SSB frame projections. - Default 0. + The polarization angle of the GW in radians. Default 0. reference_time : float (optional) The GPS start time of the GW signal in the SSB frame. Default to @@ -781,35 +819,35 @@ def get_links(self, hp, hc, lamb, beta, polarization=0, reference_time=None, """ from fastlisaresponse import pyResponseTDI - # get dt and n from waveform data - dt = hp.delta_t # timestep (s) - n = len(hp) # number of points - - # cache n, dt, pad length + # get dt from waveform data + dt = hp.delta_t self.dt = dt - self.n = n self.pad_idx = int(self.t0/self.dt) # set waveform start time and cache time series if reference_time is not None: - self.ref_time = reference_time - hp.start_time = self.ref_time - hc.start_time = self.ref_time + frame_delay = self.time_delay_from_ssb(self.orbits, reference_time) + self.com_ref_time = reference_time + frame_delay + hp.start_time = self.com_ref_time + hc.start_time = self.com_ref_time self.sample_times = hp.sample_times.numpy() - + # rescale the orbital time series to match waveform ### requires bugfix in LISAanalysistools to function ### awaiting response from dev team for formal release self.orbits.configure(t_arr=self.sample_times) - # convert to SSB frame (if pol = 0, the signal is already in SSB frame) - if polarization != 0.: - hp, hc = self.source_to_ssb(hp, hc, polarization) + # rotate GW from radiation frame to SSB using polarization angle + hp, hc = self.apply_polarization(hp, hc, polarization) # format wf to FLR standard hp + i*hc as ndarray hp = hp.numpy() hc = hc.numpy() wf = hp + 1j*hc + + # get n from waveform data + n = len(wf) + self.n = n # set use_gpu to class input if not specified if use_gpu is None: @@ -853,8 +891,8 @@ def get_links(self, hp, hc, lamb, beta, polarization=0, reference_time=None, return wf_proj - def project_wave(self, hp, hc, lamb, beta, polarization=0, reference_time=None, tdi=2, - tdi_chan='AET', tdi_orbits=None, use_gpu=None, remove_garbage=True): + def project_wave(self, hp, hc, lamb, beta, polarization, reference_time=None, tdi=2, + tdi_chan='AET', tdi_orbits=None, use_gpu=None, remove_garbage=False): """ Evaluate the TDI observables. @@ -872,10 +910,8 @@ def project_wave(self, hp, hc, lamb, beta, polarization=0, reference_time=None, beta : float The ecleptic longitude in the SSB frame. - polarization : float (optional) - The polarization angle of the GW in radians. If polarization is - 0, the hp and hc inputs are treated as SSB frame projections. - Default 0. + polarization : float + The polarization angle of the GW in radians. reference_time : float (optional) The GPS start time of the GW signal in the SSB frame. Default to @@ -907,8 +943,6 @@ def project_wave(self, hp, hc, lamb, beta, polarization=0, reference_time=None, dict The TDI observables keyed by their corresponding TDI channel. """ - start_time = time.time() - # set use_gpu if use_gpu is None: use_gpu = self.use_gpu @@ -916,6 +950,9 @@ def project_wave(self, hp, hc, lamb, beta, polarization=0, reference_time=None, # generate the Doppler time series self.get_links(hp, hc, lamb, beta, polarization=polarization, reference_time=reference_time, use_gpu=use_gpu, tdi=tdi) + + print(len(self.response_init.y_gw[0])) + print(len(self.sample_times)) # set TDI channels if tdi_chan in ['XYZ', 'AET', 'AE']: @@ -946,14 +983,15 @@ def project_wave(self, hp, hc, lamb, beta, polarization=0, reference_time=None, if remove_garbage: if remove_garbage == 'zero': # zero the edge data - tdi_dict[tdi_chan[i]][self.pad_idx:] = 0 - tdi_dict[tdi_chan[i]][:-self.pad_idx] = 0 + tdi_dict[tdi_chan[i]][:self.pad_idx] = 0 + tdi_dict[tdi_chan[i]][-self.pad_idx:] = 0 elif type(remove_garbage) == bool: # cut the edge data tdi_dict[tdi_chan[i]] = tdi_dict[tdi_chan[i]][slc] if i == 0: # cut the sample times (only do once) self.sample_times = self.sample_times[slc] + print('check') else: # raise error raise ValueError('remove_garbage arg must be a bool or "zero"') From fa7ee73c1ae70e1c3a3ea9bc212008d147ae8c6b Mon Sep 17 00:00:00 2001 From: acorreia61201 Date: Thu, 11 Jul 2024 20:31:29 +0000 Subject: [PATCH 10/24] changed t0 to kwarg in init --- pycbc/detector.py | 61 ++++++++++------------------------------------- 1 file changed, 13 insertions(+), 48 deletions(-) diff --git a/pycbc/detector.py b/pycbc/detector.py index 6b6ef7e06f4..f9b42da7465 100644 --- a/pycbc/detector.py +++ b/pycbc/detector.py @@ -662,7 +662,7 @@ class LISA_detector(object): """ LISA-like GW detector. Applies detector response from FastLISAResponse. """ - def __init__(self, detector='LISA', reference_time=None, orbits=ESAOrbits(), use_gpu=False): + def __init__(self, detector='LISA', reference_time=None, orbits=ESAOrbits(), use_gpu=False, t0=10000.): """ Parameters ---------- @@ -680,6 +680,11 @@ def __init__(self, detector='LISA', reference_time=None, orbits=ESAOrbits(), use use_gpu : bool (optional) Specify whether to run class on GPU support via CuPy. Default False. + + t0 : float (optional) + Time length (in seconds) to cut from start and end of TDI observables. By + default, TDI channels will have erroneous noise at edges that will be cut + according to this argument. Default 10000 s. """ # initialize detector; for now we only accept LISA assert (detector=='LISA'), 'Currently only supports LISA detector' @@ -693,7 +698,7 @@ def __init__(self, detector='LISA', reference_time=None, orbits=ESAOrbits(), use # specify and cache the start time and orbital time series if reference_time is None: reference_time = self.orbits.t_base[0] - self.com_ref_time = reference_time + self.ref_time = reference_time self.sample_times = None # cache the FLR instance along with dt and n @@ -701,7 +706,9 @@ def __init__(self, detector='LISA', reference_time=None, orbits=ESAOrbits(), use self.n = None self.pad_idx = None self.response_init = None - self.t0 = 10000. + + if t0 is None: + self.t0 = 10000. # initialize whether to use gpu; FLR has handles if this cannot be done self.use_gpu = use_gpu @@ -710,9 +717,6 @@ def apply_polarization(self, hp, hc, polarization): """ Apply polarization rotation matrix. - .. math:: - \bmatrix{} - Parameters ---------- hp : array @@ -736,44 +740,6 @@ def apply_polarization(self, hp, hc, polarization): hc_ssb = hp*sphi + hc*cphi return hp_ssb, hc_ssb - - def time_delay_from_ssb(self, orbits, reference_time): - """ - Calculate the time delay from the SSB frame to the LISA COM frame. - - Parameters - ---------- - orbits: lisatools.detector.Orbits - The orbital information of the satellites. - - reference_time: float - The time in seconds in the SSB frame. - - Returns - ------- - float - The additive time delay factor between the SSB and LISA COM frame - at the given time. - """ - C_SI = constants.c.value - - # configure orbits if not already - try: - orbits.t - except ValueError: - orbits.configure(linear_interp_setup=True) - - # get positions of spacecraft at reference time - sc_x = [] - for i in [1, 2, 3]: - sc_x.append(orbits.get_pos(reference_time, i)) - - # average sc positions to get COM position - com_vec = sum(sc_x)/len(sc_x) - com_dist = sum(com_vec*com_vec)**0.5 - - # time delay from SSB to LISA COM is distance/c - return com_dist/C_SI def get_links(self, hp, hc, lamb, beta, polarization=0, reference_time=None, use_gpu=None, tdi=2): @@ -826,10 +792,9 @@ def get_links(self, hp, hc, lamb, beta, polarization=0, reference_time=None, # set waveform start time and cache time series if reference_time is not None: - frame_delay = self.time_delay_from_ssb(self.orbits, reference_time) - self.com_ref_time = reference_time + frame_delay - hp.start_time = self.com_ref_time - hc.start_time = self.com_ref_time + self.ref_time = reference_time + hp.start_time = self.ref_time + hc.start_time = self.ref_time self.sample_times = hp.sample_times.numpy() # rescale the orbital time series to match waveform From 8cc2d5984a3fd99e879dabe2fa3c79e676e6fff3 Mon Sep 17 00:00:00 2001 From: acorreia61201 Date: Fri, 12 Jul 2024 16:15:50 +0000 Subject: [PATCH 11/24] fix TypeError in unittest --- pycbc/detector.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/pycbc/detector.py b/pycbc/detector.py index f9b42da7465..155f692e12e 100644 --- a/pycbc/detector.py +++ b/pycbc/detector.py @@ -697,8 +697,10 @@ def __init__(self, detector='LISA', reference_time=None, orbits=ESAOrbits(), use # specify and cache the start time and orbital time series if reference_time is None: - reference_time = self.orbits.t_base[0] - self.ref_time = reference_time + ref_time = self.orbits.t_base[0] + else: + ref_time = reference_time + self.ref_time = ref_time self.sample_times = None # cache the FLR instance along with dt and n From 90ed9a8f64638e588df85b8358d8bb50f3dc260c Mon Sep 17 00:00:00 2001 From: acorreia61201 Date: Fri, 12 Jul 2024 16:46:08 +0000 Subject: [PATCH 12/24] fix bug when calling ESAOrbits without lisatools --- pycbc/detector.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/pycbc/detector.py b/pycbc/detector.py index 155f692e12e..d1da98788d1 100644 --- a/pycbc/detector.py +++ b/pycbc/detector.py @@ -662,7 +662,7 @@ class LISA_detector(object): """ LISA-like GW detector. Applies detector response from FastLISAResponse. """ - def __init__(self, detector='LISA', reference_time=None, orbits=ESAOrbits(), use_gpu=False, t0=10000.): + def __init__(self, detector='LISA', reference_time=None, orbits=None, use_gpu=False, t0=10000.): """ Parameters ---------- @@ -691,7 +691,10 @@ def __init__(self, detector='LISA', reference_time=None, orbits=ESAOrbits(), use # intialize orbit information if orbits is None: - raise ImportError('LISAanalysistools required for inputting orbital data') + # set ESAOrbits as default; raise error if ESAOrbits cannot be imported + if ESAOrbits is None: + raise ImportError('LISAanalysistools required for inputting orbital data') + orbits = ESAOrbits() orbits.configure(linear_interp_setup=True) self.orbits = orbits From b36b378b8c3ad998087f84cc70d730c89e228143 Mon Sep 17 00:00:00 2001 From: acorreia61201 Date: Fri, 30 Aug 2024 12:57:55 -0400 Subject: [PATCH 13/24] fix case where t0 is not None; remove debug statements/white space; edit companion (temp pypmc fix, add FLR) --- companion.txt | 4 +++ pycbc/detector.py | 66 +++++++++++++++++++++++------------------------ 2 files changed, 36 insertions(+), 34 deletions(-) diff --git a/companion.txt b/companion.txt index 769c574d315..1eb5d5988b4 100644 --- a/companion.txt +++ b/companion.txt @@ -18,6 +18,7 @@ https://github.com/willvousden/ptemcee/archive/master.tar.gz --extra-index-url https://download.pytorch.org/whl/cpu torch nessai>=0.11.0 +pypmc==1.2.2 snowline # useful to look at PyCBC Live with htop @@ -29,3 +30,6 @@ sympy>=1.9 # Needed for KDE trigger statistics git+https://github.com/mennthor/awkde.git@master scikit-learn + +# LISA response +fastlisaresponse diff --git a/pycbc/detector.py b/pycbc/detector.py index d1da98788d1..805a0271780 100644 --- a/pycbc/detector.py +++ b/pycbc/detector.py @@ -669,18 +669,18 @@ def __init__(self, detector='LISA', reference_time=None, orbits=None, use_gpu=Fa detector: str (optional) String specifying space-borne detector. Currently only accepts 'LISA', which is the default setting. - + reference_time: float (optional) The GPS start time of signal in the SSB frame. Default to start time of orbits input. - + orbits: lisatools.detector.Orbits (optional) Orbital information to pass into pyResponseTDI. Default lisatools.detector.ESAOrbits. - + use_gpu : bool (optional) Specify whether to run class on GPU support via CuPy. Default False. - + t0 : float (optional) Time length (in seconds) to cut from start and end of TDI observables. By default, TDI channels will have erroneous noise at edges that will be cut @@ -688,7 +688,7 @@ def __init__(self, detector='LISA', reference_time=None, orbits=None, use_gpu=Fa """ # initialize detector; for now we only accept LISA assert (detector=='LISA'), 'Currently only supports LISA detector' - + # intialize orbit information if orbits is None: # set ESAOrbits as default; raise error if ESAOrbits cannot be imported @@ -697,7 +697,7 @@ def __init__(self, detector='LISA', reference_time=None, orbits=None, use_gpu=Fa orbits = ESAOrbits() orbits.configure(linear_interp_setup=True) self.orbits = orbits - + # specify and cache the start time and orbital time series if reference_time is None: ref_time = self.orbits.t_base[0] @@ -711,17 +711,19 @@ def __init__(self, detector='LISA', reference_time=None, orbits=None, use_gpu=Fa self.n = None self.pad_idx = None self.response_init = None - + if t0 is None: self.t0 = 10000. + else: + self.t0 = t0 # initialize whether to use gpu; FLR has handles if this cannot be done self.use_gpu = use_gpu - + def apply_polarization(self, hp, hc, polarization): """ Apply polarization rotation matrix. - + Parameters ---------- hp : array @@ -729,10 +731,10 @@ def apply_polarization(self, hp, hc, polarization): hc : array The cross polarization of the GW. - + polarization : float The polarization angle of the GW in radians. - + Returns ------- (array, array) @@ -740,10 +742,10 @@ def apply_polarization(self, hp, hc, polarization): """ cphi = cos(2*polarization) sphi = sin(2*polarization) - + hp_ssb = hp*cphi - hc*sphi hc_ssb = hp*sphi + hc*cphi - + return hp_ssb, hc_ssb def get_links(self, hp, hc, lamb, beta, polarization=0, reference_time=None, @@ -766,10 +768,10 @@ def get_links(self, hp, hc, lamb, beta, polarization=0, reference_time=None, beta : float The ecleptic longitude of the source in the SSB frame. - + polarization : float (optional) The polarization angle of the GW in radians. Default 0. - + reference_time : float (optional) The GPS start time of the GW signal in the SSB frame. Default to input on class initialization. @@ -778,7 +780,7 @@ def get_links(self, hp, hc, lamb, beta, polarization=0, reference_time=None, Flag whether to use GPU support. Default to class input. CuPy is required if use_gpu is True; an ImportError will be raised if CuPy could not be imported. - + tdi : int (optional) TDI channel configuration. Accepts 1 for 1st generation TDI or 2 for 2nd generation TDI. Default 2. @@ -789,7 +791,7 @@ def get_links(self, hp, hc, lamb, beta, polarization=0, reference_time=None, The waveform projected to the LISA laser links. """ from fastlisaresponse import pyResponseTDI - + # get dt from waveform data dt = hp.delta_t self.dt = dt @@ -801,12 +803,12 @@ def get_links(self, hp, hc, lamb, beta, polarization=0, reference_time=None, hp.start_time = self.ref_time hc.start_time = self.ref_time self.sample_times = hp.sample_times.numpy() - + # rescale the orbital time series to match waveform ### requires bugfix in LISAanalysistools to function ### awaiting response from dev team for formal release self.orbits.configure(t_arr=self.sample_times) - + # rotate GW from radiation frame to SSB using polarization angle hp, hc = self.apply_polarization(hp, hc, polarization) @@ -814,7 +816,7 @@ def get_links(self, hp, hc, lamb, beta, polarization=0, reference_time=None, hp = hp.numpy() hc = hc.numpy() wf = hp + 1j*hc - + # get n from waveform data n = len(wf) self.n = n @@ -829,7 +831,7 @@ def get_links(self, hp, hc, lamb, beta, polarization=0, reference_time=None, raise ImportError('CuPy not imported but is required for GPU usage') else: wf = cupy.asarray(wf) - + # set TDI configuration (let FLR handle if not 1 or 2) if tdi == 1: tdi_opt = '1st generation' @@ -849,7 +851,7 @@ def get_links(self, hp, hc, lamb, beta, polarization=0, reference_time=None, self.response_init.num_pts = n self.response_init.orbits = self.orbits self.response_init.use_gpu = use_gpu - + if tdi_opt != self.response_init.tdi: # update TDI in existing response_init class self.response_init.tdi = tdi_opt @@ -879,14 +881,14 @@ def project_wave(self, hp, hc, lamb, beta, polarization, reference_time=None, td beta : float The ecleptic longitude in the SSB frame. - + polarization : float The polarization angle of the GW in radians. - + reference_time : float (optional) The GPS start time of the GW signal in the SSB frame. Default to value in input signals hp and hc. - + tdi : int (optional) TDI channel configuration. Accepts 1 for 1st generation TDI or 2 for 2nd generation TDI. Default 2. @@ -901,7 +903,7 @@ def project_wave(self, hp, hc, lamb, beta, polarization, reference_time=None, td use_gpu : bool (optional) Flag whether to use GPU support. Default to class input. - + remove_garbage : bool, str (optional) Flag whether to remove gaps in TDI from start and end. If True, edge data will be cut from TDI channels. If 'zero', edge data @@ -916,13 +918,10 @@ def project_wave(self, hp, hc, lamb, beta, polarization, reference_time=None, td # set use_gpu if use_gpu is None: use_gpu = self.use_gpu - + # generate the Doppler time series self.get_links(hp, hc, lamb, beta, polarization=polarization, reference_time=reference_time, use_gpu=use_gpu, tdi=tdi) - - print(len(self.response_init.y_gw[0])) - print(len(self.sample_times)) # set TDI channels if tdi_chan in ['XYZ', 'AET', 'AE']: @@ -940,15 +939,15 @@ def project_wave(self, hp, hc, lamb, beta, polarization, reference_time=None, td # generate the TDI channels tdi_obs = self.response_init.get_tdi_delays() - + # processing tdi_dict = {} slc = slice(self.pad_idx, -self.pad_idx) - + for i in range(len(tdi_chan)): # add to dict tdi_dict[tdi_chan[i]] = tdi_obs[i] - + # treat start and end gaps (if remove_garbage is not False) if remove_garbage: if remove_garbage == 'zero': @@ -961,7 +960,6 @@ def project_wave(self, hp, hc, lamb, beta, polarization, reference_time=None, td if i == 0: # cut the sample times (only do once) self.sample_times = self.sample_times[slc] - print('check') else: # raise error raise ValueError('remove_garbage arg must be a bool or "zero"') From 3b99cd35b25f267b99daa2d7b832ee0860968154 Mon Sep 17 00:00:00 2001 From: acorreia61201 Date: Fri, 30 Aug 2024 16:31:09 -0400 Subject: [PATCH 14/24] Indentation fix in init --- pycbc/detector.py | 38 +++++++++++++++++--------------------- 1 file changed, 17 insertions(+), 21 deletions(-) diff --git a/pycbc/detector.py b/pycbc/detector.py index 805a0271780..e4302846fc8 100644 --- a/pycbc/detector.py +++ b/pycbc/detector.py @@ -662,28 +662,28 @@ class LISA_detector(object): """ LISA-like GW detector. Applies detector response from FastLISAResponse. """ - def __init__(self, detector='LISA', reference_time=None, orbits=None, use_gpu=False, t0=10000.): + def __init__(self, detector='LISA', reference_time=None, orbits=None, use_gpu=False, t0=None): """ Parameters ---------- detector: str (optional) - String specifying space-borne detector. Currently only accepts 'LISA', + String specifying space-borne detector. Currently only accepts 'LISA', which is the default setting. reference_time: float (optional) - The GPS start time of signal in the SSB frame. Default to start time of + The GPS start time of signal in the SSB frame. Default to start time of orbits input. orbits: lisatools.detector.Orbits (optional) - Orbital information to pass into pyResponseTDI. Default + Orbital information to pass into pyResponseTDI. Default lisatools.detector.ESAOrbits. use_gpu : bool (optional) Specify whether to run class on GPU support via CuPy. Default False. t0 : float (optional) - Time length (in seconds) to cut from start and end of TDI observables. By - default, TDI channels will have erroneous noise at edges that will be cut + Time length (in seconds) to cut from start and end of TDI observables. By + default, TDI channels will have erroneous noise at edges that will be cut according to this argument. Default 10000 s. """ # initialize detector; for now we only accept LISA @@ -711,11 +711,11 @@ def __init__(self, detector='LISA', reference_time=None, orbits=None, use_gpu=Fa self.n = None self.pad_idx = None self.response_init = None - + if t0 is None: self.t0 = 10000. - else: - self.t0 = t0 + else: + self.t0 = t0 # initialize whether to use gpu; FLR has handles if this cannot be done self.use_gpu = use_gpu @@ -748,7 +748,7 @@ def apply_polarization(self, hp, hc, polarization): return hp_ssb, hc_ssb - def get_links(self, hp, hc, lamb, beta, polarization=0, reference_time=None, + def get_links(self, hp, hc, lamb, beta, polarization=0, reference_time=None, use_gpu=None, tdi=2): """ Project a radiation frame waveform to the LISA constellation. @@ -773,12 +773,12 @@ def get_links(self, hp, hc, lamb, beta, polarization=0, reference_time=None, The polarization angle of the GW in radians. Default 0. reference_time : float (optional) - The GPS start time of the GW signal in the SSB frame. Default to + The GPS start time of the GW signal in the SSB frame. Default to input on class initialization. use_gpu : bool (optional) Flag whether to use GPU support. Default to class input. - CuPy is required if use_gpu is True; an ImportError will be raised + CuPy is required if use_gpu is True; an ImportError will be raised if CuPy could not be imported. tdi : int (optional) @@ -805,8 +805,6 @@ def get_links(self, hp, hc, lamb, beta, polarization=0, reference_time=None, self.sample_times = hp.sample_times.numpy() # rescale the orbital time series to match waveform - ### requires bugfix in LISAanalysistools to function - ### awaiting response from dev team for formal release self.orbits.configure(t_arr=self.sample_times) # rotate GW from radiation frame to SSB using polarization angle @@ -842,7 +840,7 @@ def get_links(self, hp, hc, lamb, beta, polarization=0, reference_time=None, if self.response_init is None: # initialize the class - response_init = pyResponseTDI(1/dt, n, orbits=self.orbits, + response_init = pyResponseTDI(1/dt, n, orbits=self.orbits, use_gpu=use_gpu, tdi=tdi_opt) self.response_init = response_init else: @@ -863,7 +861,7 @@ def get_links(self, hp, hc, lamb, beta, polarization=0, reference_time=None, return wf_proj - def project_wave(self, hp, hc, lamb, beta, polarization, reference_time=None, tdi=2, + def project_wave(self, hp, hc, lamb, beta, polarization, reference_time=None, tdi=2, tdi_chan='AET', tdi_orbits=None, use_gpu=None, remove_garbage=False): """ Evaluate the TDI observables. @@ -886,7 +884,7 @@ def project_wave(self, hp, hc, lamb, beta, polarization, reference_time=None, td The polarization angle of the GW in radians. reference_time : float (optional) - The GPS start time of the GW signal in the SSB frame. Default to + The GPS start time of the GW signal in the SSB frame. Default to value in input signals hp and hc. tdi : int (optional) @@ -898,7 +896,7 @@ def project_wave(self, hp, hc, lamb, beta, polarization, reference_time=None, td Default 'AET'. tdi_orbits : lisatools.detector.Orbits (optional) - Orbit keywords specifically for TDI generation. Default to class + Orbit keywords specifically for TDI generation. Default to class input. use_gpu : bool (optional) @@ -920,7 +918,7 @@ def project_wave(self, hp, hc, lamb, beta, polarization, reference_time=None, td use_gpu = self.use_gpu # generate the Doppler time series - self.get_links(hp, hc, lamb, beta, polarization=polarization, reference_time=reference_time, + self.get_links(hp, hc, lamb, beta, polarization=polarization, reference_time=reference_time, use_gpu=use_gpu, tdi=tdi) # set TDI channels @@ -932,8 +930,6 @@ def project_wave(self, hp, hc, lamb, beta, polarization, reference_time=None, td # if TDI orbit class is provided, update the response_init # tdi_orbits are set to class input automatically by FLR otherwise if tdi_orbits is not None: - ### requires bugfix in LISAanalysistools to function - ### awaiting response from dev team for formal release tdi_orbits.configure(t_arr=self.sample_times) self.response_init.tdi_orbits = tdi_orbits From c1123ffb2dbbc49701f03119f3a8d4f18119a87a Mon Sep 17 00:00:00 2001 From: acorreia61201 Date: Tue, 3 Sep 2024 16:14:43 -0400 Subject: [PATCH 15/24] add padding (for signals that taper to zero); more cleanup --- pycbc/detector.py | 30 +++++++++++++++++++++--------- 1 file changed, 21 insertions(+), 9 deletions(-) diff --git a/pycbc/detector.py b/pycbc/detector.py index e4302846fc8..bcf6d8be10d 100644 --- a/pycbc/detector.py +++ b/pycbc/detector.py @@ -711,7 +711,7 @@ def __init__(self, detector='LISA', reference_time=None, orbits=None, use_gpu=Fa self.n = None self.pad_idx = None self.response_init = None - + if t0 is None: self.t0 = 10000. else: @@ -795,9 +795,20 @@ def get_links(self, hp, hc, lamb, beta, polarization=0, reference_time=None, # get dt from waveform data dt = hp.delta_t self.dt = dt + + # index corresponding to time length t0 self.pad_idx = int(self.t0/self.dt) - # set waveform start time and cache time series + # pad the data with zeros + ### this assumes that the signal naturally tapers to zero + ### this will not work with e.g. GBs or sinusoids + ### i.g. we should be tapering to prevent discontinuities + hp.prepend_zeros(self.pad_idx) + hp.append_zeros(self.pad_idx) + hc.prepend_zeros(self.pad_idx) + hc.append_zeros(self.pad_idx) + + # set waveform start time and cache time samples if reference_time is not None: self.ref_time = reference_time hp.start_time = self.ref_time @@ -815,7 +826,7 @@ def get_links(self, hp, hc, lamb, beta, polarization=0, reference_time=None, hc = hc.numpy() wf = hp + 1j*hc - # get n from waveform data + # save length of wf n = len(wf) self.n = n @@ -862,7 +873,7 @@ def get_links(self, hp, hc, lamb, beta, polarization=0, reference_time=None, return wf_proj def project_wave(self, hp, hc, lamb, beta, polarization, reference_time=None, tdi=2, - tdi_chan='AET', tdi_orbits=None, use_gpu=None, remove_garbage=False): + tdi_chan='AET', tdi_orbits=None, use_gpu=None, remove_garbage=True): """ Evaluate the TDI observables. @@ -904,9 +915,9 @@ def project_wave(self, hp, hc, lamb, beta, polarization, reference_time=None, td remove_garbage : bool, str (optional) Flag whether to remove gaps in TDI from start and end. If True, - edge data will be cut from TDI channels. If 'zero', edge data - will be zeroed. If False, TDI channels will not be modified. - Default True. + time length t0 worth of edge data 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 True. Returns ------- @@ -918,8 +929,9 @@ def project_wave(self, hp, hc, lamb, beta, polarization, reference_time=None, td use_gpu = self.use_gpu # generate the Doppler time series - self.get_links(hp, hc, lamb, beta, polarization=polarization, reference_time=reference_time, - use_gpu=use_gpu, tdi=tdi) + self.get_links(hp, hc, lamb, beta, polarization=polarization, + reference_time=reference_time, use_gpu=use_gpu, + tdi=tdi) # set TDI channels if tdi_chan in ['XYZ', 'AET', 'AE']: From 0d3f8e62e0d02d370bc466ebf5d9d0c88d1bdf53 Mon Sep 17 00:00:00 2001 From: Alex Correia Date: Thu, 5 Sep 2024 13:51:42 -0400 Subject: [PATCH 16/24] add controls for zero padding data to project_wave; add controls for specifying orbit start times for numerical data --- pycbc/detector.py | 131 ++++++++++++++++++++++++++++++---------------- 1 file changed, 87 insertions(+), 44 deletions(-) diff --git a/pycbc/detector.py b/pycbc/detector.py index bcf6d8be10d..14ea1db03e8 100644 --- a/pycbc/detector.py +++ b/pycbc/detector.py @@ -662,7 +662,8 @@ class LISA_detector(object): """ LISA-like GW detector. Applies detector response from FastLISAResponse. """ - def __init__(self, detector='LISA', reference_time=None, orbits=None, use_gpu=False, t0=None): + def __init__(self, detector='LISA', reference_time=None, orbits=None, + orbits_reference_time=None, use_gpu=False, t0=None): """ Parameters ---------- @@ -678,6 +679,10 @@ def __init__(self, detector='LISA', reference_time=None, orbits=None, use_gpu=Fa Orbital information to pass into pyResponseTDI. Default lisatools.detector.ESAOrbits. + orbits_reference_time: float (optional) + The GPS start time of the orbital data in the SSB frame. Default to the + start time given in the orbital data. + use_gpu : bool (optional) Specify whether to run class on GPU support via CuPy. Default False. @@ -695,29 +700,37 @@ def __init__(self, detector='LISA', reference_time=None, orbits=None, use_gpu=Fa if ESAOrbits is None: raise ImportError('LISAanalysistools required for inputting orbital data') orbits = ESAOrbits() - orbits.configure(linear_interp_setup=True) self.orbits = orbits - # specify and cache the start time and orbital time series + # specify a start time for the orbit data + if orbits_reference_time is None: + self.orb_ref_time = self.orbits.t_base[0] + else: + self.orb_ref_time = orbits_reference_time + + # specify and cache the start time if reference_time is None: - ref_time = self.orbits.t_base[0] + self.ref_time = self.orbits.t_base[0] else: - ref_time = reference_time - self.ref_time = ref_time - self.sample_times = None + assert (reference_time >= self.orb_ref_time) and ( + reference_time <= self.orbits.t_base[-1]-self.orbits.t_base[0]+self.orb_ref_time), ( + "Signal reference time is not in time domain of orbital data") + self.ref_time = reference_time - # cache the FLR instance along with dt and n + # allocate caches self.dt = None self.n = None self.pad_idx = None + self.sample_times = None self.response_init = None + # initialize padding/cutting time length if t0 is None: self.t0 = 10000. else: self.t0 = t0 - # initialize whether to use gpu; FLR has handles if this cannot be done + # initialize whether to use gpu; FLR handles if this cannot be done self.use_gpu = use_gpu def apply_polarization(self, hp, hc, polarization): @@ -764,10 +777,10 @@ def get_links(self, hp, hc, lamb, beta, polarization=0, reference_time=None, The cross polarization of the GW in the radiation frame. lamb : float - The ecleptic latitude of the source in the SSB frame. + The ecliptic latitude of the source in the SSB frame. beta : float - The ecleptic longitude of the source in the SSB frame. + The ecliptic longitude of the source in the SSB frame. polarization : float (optional) The polarization angle of the GW in radians. Default 0. @@ -792,21 +805,9 @@ def get_links(self, hp, hc, lamb, beta, polarization=0, reference_time=None, """ from fastlisaresponse import pyResponseTDI - # get dt from waveform data - dt = hp.delta_t - self.dt = dt - - # index corresponding to time length t0 - self.pad_idx = int(self.t0/self.dt) - - # pad the data with zeros - ### this assumes that the signal naturally tapers to zero - ### this will not work with e.g. GBs or sinusoids - ### i.g. we should be tapering to prevent discontinuities - hp.prepend_zeros(self.pad_idx) - hp.append_zeros(self.pad_idx) - hc.prepend_zeros(self.pad_idx) - hc.append_zeros(self.pad_idx) + # get dt from waveform data (skipped if running project_wave) + if self.dt is None: + self.dt = hp.delta_t # set waveform start time and cache time samples if reference_time is not None: @@ -815,20 +816,28 @@ def get_links(self, hp, hc, lamb, beta, polarization=0, reference_time=None, hc.start_time = self.ref_time self.sample_times = hp.sample_times.numpy() + # For simulated orbital data that starts at t=0, we shift the + # signal times back by the given orbital ref time. This + # converts the GPS time inputs to be relative to the orbital start + # time. For data with associated GPS times, this won't be necessary. + if self.orbits.t_base[0] == 0.: + t_arr = self.sample_times - self.orb_ref_time + else: + t_arr = self.sample_times + # rescale the orbital time series to match waveform - self.orbits.configure(t_arr=self.sample_times) + self.orbits.configure(t_arr=t_arr) # rotate GW from radiation frame to SSB using polarization angle hp, hc = self.apply_polarization(hp, hc, polarization) - # format wf to FLR standard hp + i*hc as ndarray + # format wf to hp + i*hc hp = hp.numpy() hc = hc.numpy() wf = hp + 1j*hc # save length of wf - n = len(wf) - self.n = n + self.n = len(wf) # set use_gpu to class input if not specified if use_gpu is None: @@ -851,13 +860,12 @@ def get_links(self, hp, hc, lamb, beta, polarization=0, reference_time=None, if self.response_init is None: # initialize the class - response_init = pyResponseTDI(1/dt, n, orbits=self.orbits, - use_gpu=use_gpu, tdi=tdi_opt) - self.response_init = response_init + self.response_init = pyResponseTDI(1/self.dt, self.n, orbits=self.orbits, + use_gpu=use_gpu, tdi=tdi_opt) else: # update params in the initialized class - self.response_init.sampling_frequency = 1/dt - self.response_init.num_pts = n + self.response_init.sampling_frequency = 1/self.dt + self.response_init.num_pts = self.n self.response_init.orbits = self.orbits self.response_init.use_gpu = use_gpu @@ -872,11 +880,20 @@ def get_links(self, hp, hc, lamb, beta, polarization=0, reference_time=None, return wf_proj - def project_wave(self, hp, hc, lamb, beta, polarization, reference_time=None, tdi=2, - tdi_chan='AET', tdi_orbits=None, use_gpu=None, remove_garbage=True): + def project_wave(self, hp, hc, lamb, beta, polarization, reference_time=None, + tdi=2, tdi_chan='AET', tdi_orbits=None, use_gpu=None, + pad_data=False, remove_garbage=True): """ 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 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 ---------- hp : pycbc.types.TimeSeries @@ -886,10 +903,10 @@ def project_wave(self, hp, hc, lamb, beta, polarization, reference_time=None, td The cross polarization of the GW in the radiation frame. lamb : float - The ecleptic latitude in the SSB frame. + The ecliptic latitude in the SSB frame. beta : float - The ecleptic longitude in the SSB frame. + The ecliptic longitude in the SSB frame. polarization : float The polarization angle of the GW in radians. @@ -913,17 +930,39 @@ def project_wave(self, hp, hc, lamb, beta, polarization, reference_time=None, td use_gpu : bool (optional) Flag whether to use GPU support. Default to class input. + pad_data : bool (optional) + Flag whether to pad the data with time length t0 of zeros at the + start and end. If True, remove_garbage will interact with the + pads rather than the input data. Default False. + remove_garbage : bool, str (optional) Flag whether to remove gaps in TDI from start and end. If True, - time length t0 worth of edge data 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 True. + 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 True. Returns ------- dict The TDI observables keyed by their corresponding TDI channel. """ + # get dt from waveform data + self.dt = hp.delta_t + + # index corresponding to time length t0 + self.pad_idx = int(self.t0/self.dt) + + # pad the data with zeros + ### this assumes that the signal naturally tapers to zero + ### this will not work with e.g. GBs or sinusoids + ### we should be tapering in general to prevent discontinuities + if pad_data: + hp.prepend_zeros(self.pad_idx) + hp.append_zeros(self.pad_idx) + hc.prepend_zeros(self.pad_idx) + hc.append_zeros(self.pad_idx) + # set use_gpu if use_gpu is None: use_gpu = self.use_gpu @@ -942,7 +981,12 @@ def project_wave(self, hp, hc, lamb, beta, polarization, reference_time=None, td # if TDI orbit class is provided, update the response_init # tdi_orbits are set to class input automatically by FLR otherwise if tdi_orbits is not None: - tdi_orbits.configure(t_arr=self.sample_times) + # handle for analytic orbits + if tdi_orbits.t_base[0] == 0: + t_arr = self.sample_times - self.orbit_ref_time + else: + t_arr = self.sample_times + tdi_orbits.configure(t_arr=t_arr) self.response_init.tdi_orbits = tdi_orbits # generate the TDI channels @@ -969,7 +1013,6 @@ def project_wave(self, hp, hc, lamb, beta, polarization, reference_time=None, td # cut the sample times (only do once) self.sample_times = self.sample_times[slc] else: - # raise error raise ValueError('remove_garbage arg must be a bool or "zero"') return tdi_dict From 7ff76660621854427dc9f11d0d7fbbe3739de7ee Mon Sep 17 00:00:00 2001 From: Alex Correia Date: Mon, 9 Sep 2024 18:31:33 -0400 Subject: [PATCH 17/24] remove orbit reference time; move TDI arguments to project_wave --- pycbc/detector.py | 78 ++++++++++++++++------------------------------- 1 file changed, 26 insertions(+), 52 deletions(-) diff --git a/pycbc/detector.py b/pycbc/detector.py index 14ea1db03e8..e490eff1897 100644 --- a/pycbc/detector.py +++ b/pycbc/detector.py @@ -663,7 +663,7 @@ class LISA_detector(object): LISA-like GW detector. Applies detector response from FastLISAResponse. """ def __init__(self, detector='LISA', reference_time=None, orbits=None, - orbits_reference_time=None, use_gpu=False, t0=None): + use_gpu=False, t0=None): """ Parameters ---------- @@ -679,10 +679,6 @@ def __init__(self, detector='LISA', reference_time=None, orbits=None, Orbital information to pass into pyResponseTDI. Default lisatools.detector.ESAOrbits. - orbits_reference_time: float (optional) - The GPS start time of the orbital data in the SSB frame. Default to the - start time given in the orbital data. - use_gpu : bool (optional) Specify whether to run class on GPU support via CuPy. Default False. @@ -701,20 +697,15 @@ def __init__(self, detector='LISA', reference_time=None, orbits=None, raise ImportError('LISAanalysistools required for inputting orbital data') orbits = ESAOrbits() self.orbits = orbits - - # specify a start time for the orbit data - if orbits_reference_time is None: - self.orb_ref_time = self.orbits.t_base[0] - else: - self.orb_ref_time = orbits_reference_time + orbit_start = self.orbits.t_base[0] + orbit_end = self.orbits.t_base[-1] # specify and cache the start time if reference_time is None: - self.ref_time = self.orbits.t_base[0] + self.ref_time = orbit_start else: - assert (reference_time >= self.orb_ref_time) and ( - reference_time <= self.orbits.t_base[-1]-self.orbits.t_base[0]+self.orb_ref_time), ( - "Signal reference time is not in time domain of orbital data") + assert (reference_time >= orbit_start) and (reference_time <= orbit_end), ( + "Reference time is not in time domain of orbital data") self.ref_time = reference_time # allocate caches @@ -762,7 +753,7 @@ def apply_polarization(self, hp, hc, polarization): return hp_ssb, hc_ssb def get_links(self, hp, hc, lamb, beta, polarization=0, reference_time=None, - use_gpu=None, tdi=2): + use_gpu=None): """ Project a radiation frame waveform to the LISA constellation. This does not perform TDI; this only produces a detector frame @@ -794,10 +785,6 @@ def get_links(self, hp, hc, lamb, beta, polarization=0, reference_time=None, CuPy is required if use_gpu is True; an ImportError will be raised if CuPy could not be imported. - tdi : int (optional) - TDI channel configuration. Accepts 1 for 1st generation TDI or - 2 for 2nd generation TDI. Default 2. - Returns ------- ndarray @@ -816,17 +803,10 @@ def get_links(self, hp, hc, lamb, beta, polarization=0, reference_time=None, hc.start_time = self.ref_time self.sample_times = hp.sample_times.numpy() - # For simulated orbital data that starts at t=0, we shift the - # signal times back by the given orbital ref time. This - # converts the GPS time inputs to be relative to the orbital start - # time. For data with associated GPS times, this won't be necessary. - if self.orbits.t_base[0] == 0.: - t_arr = self.sample_times - self.orb_ref_time - else: - t_arr = self.sample_times - # rescale the orbital time series to match waveform - self.orbits.configure(t_arr=t_arr) + assert hp.duration + self.ref_time <= self.orbits.t_base[-1], ( + "Time of signal end is greater than end of orbital data") + self.orbits.configure(t_arr=self.sample_times) # rotate GW from radiation frame to SSB using polarization angle hp, hc = self.apply_polarization(hp, hc, polarization) @@ -850,18 +830,10 @@ def get_links(self, hp, hc, lamb, beta, polarization=0, reference_time=None, else: wf = cupy.asarray(wf) - # set TDI configuration (let FLR handle if not 1 or 2) - if tdi == 1: - tdi_opt = '1st generation' - elif tdi == 2: - tdi_opt = '2nd generation' - else: - tdi_opt = tdi - if self.response_init is None: # initialize the class self.response_init = pyResponseTDI(1/self.dt, self.n, orbits=self.orbits, - use_gpu=use_gpu, tdi=tdi_opt) + use_gpu=use_gpu) else: # update params in the initialized class self.response_init.sampling_frequency = 1/self.dt @@ -869,11 +841,6 @@ def get_links(self, hp, hc, lamb, beta, polarization=0, reference_time=None, self.response_init.orbits = self.orbits self.response_init.use_gpu = use_gpu - if tdi_opt != self.response_init.tdi: - # update TDI in existing response_init class - self.response_init.tdi = tdi_opt - self.response_init._init_TDI_delays() - # project the signal self.response_init.get_projections(wf, lamb, beta, t0=self.t0) wf_proj = self.response_init.y_gw @@ -969,8 +936,20 @@ def project_wave(self, hp, hc, lamb, beta, polarization, reference_time=None, # generate the Doppler time series self.get_links(hp, hc, lamb, beta, polarization=polarization, - reference_time=reference_time, use_gpu=use_gpu, - tdi=tdi) + reference_time=reference_time, use_gpu=use_gpu) + + # set TDI configuration (let FLR handle if not 1 or 2) + if tdi == 1: + tdi_opt = '1st generation' + elif tdi == 2: + tdi_opt = '2nd generation' + else: + tdi_opt = tdi + + if tdi_opt != self.response_init.tdi: + # update TDI in existing response_init class + self.response_init.tdi = tdi_opt + self.response_init._init_TDI_delays() # set TDI channels if tdi_chan in ['XYZ', 'AET', 'AE']: @@ -981,12 +960,7 @@ def project_wave(self, hp, hc, lamb, beta, polarization, reference_time=None, # if TDI orbit class is provided, update the response_init # tdi_orbits are set to class input automatically by FLR otherwise if tdi_orbits is not None: - # handle for analytic orbits - if tdi_orbits.t_base[0] == 0: - t_arr = self.sample_times - self.orbit_ref_time - else: - t_arr = self.sample_times - tdi_orbits.configure(t_arr=t_arr) + tdi_orbits.configure(t_arr=self.sample_times) self.response_init.tdi_orbits = tdi_orbits # generate the TDI channels From 533d253aa7de7747a704f718277989ad7ffa785f Mon Sep 17 00:00:00 2001 From: Alex Correia Date: Tue, 10 Sep 2024 13:34:33 -0400 Subject: [PATCH 18/24] add time conversions (GPS time offset, SSB to LISA); moved more padding code from get_links to project_wave --- pycbc/detector.py | 141 ++++++++++++++++++++++++++++++++++------------ 1 file changed, 104 insertions(+), 37 deletions(-) diff --git a/pycbc/detector.py b/pycbc/detector.py index e490eff1897..bbdaf14df3e 100644 --- a/pycbc/detector.py +++ b/pycbc/detector.py @@ -648,6 +648,9 @@ def overhead_antenna_pattern(right_ascension, declination, polarization): """ LISA class """ + +from pycbc.coordinates.space import ssb_to_lisa, TIME_OFFSET_20_DEGREES + try: import cupy except ImportError: @@ -657,13 +660,18 @@ def overhead_antenna_pattern(right_ascension, declination, polarization): from lisatools.detector import ESAOrbits except ImportError: ESAOrbits = None + +try: + from fastlisaresponse import pyResponseTDI +except ImportError: + pyResponseTDI = None class LISA_detector(object): """ LISA-like GW detector. Applies detector response from FastLISAResponse. """ def __init__(self, detector='LISA', reference_time=None, orbits=None, - use_gpu=False, t0=None): + use_gpu=False, apply_offset=False, offset=TIME_OFFSET_20_DEGREES): """ Parameters ---------- @@ -682,10 +690,16 @@ def __init__(self, detector='LISA', reference_time=None, orbits=None, use_gpu : bool (optional) Specify whether to run class on GPU support via CuPy. Default False. - t0 : float (optional) - Time length (in seconds) to cut from start and end of TDI observables. By - default, TDI channels will have erroneous noise at edges that will be cut - according to this argument. Default 10000 s. + apply_offset : bool (optional) + Specify whether to add a time offset to input times. If True, + times are treated as GPS times, to which an offset is added to ensure + LISA is ~20 degrees behind the Earth at t=0 (the LDC convention for the + LISA start date). If False, no offset is applied. Default False. + + offset : float (optional) + Time offset in seconds to apply to GPS times if apply_offset = True. + Default pycbc.coordinates.space.TIME_OFFSET_20_DEGREES (places LISA + ~20 deg behind Earth). """ # initialize detector; for now we only accept LISA assert (detector=='LISA'), 'Currently only supports LISA detector' @@ -700,10 +714,17 @@ def __init__(self, detector='LISA', reference_time=None, orbits=None, orbit_start = self.orbits.t_base[0] orbit_end = self.orbits.t_base[-1] + # specify whether to apply offsets to GPS times + if apply_offset: + self.offset = offset + else: + self.offset = 0. + # specify and cache the start time if reference_time is None: - self.ref_time = orbit_start + self.ref_time = orbit_start + self.offset else: + reference_time += self.offset assert (reference_time >= orbit_start) and (reference_time <= orbit_end), ( "Reference time is not in time domain of orbital data") self.ref_time = reference_time @@ -711,19 +732,21 @@ def __init__(self, detector='LISA', reference_time=None, orbits=None, # allocate caches self.dt = None self.n = None - self.pad_idx = None + self.pad_data = False # don't pad by default if only projecting self.sample_times = None self.response_init = None # initialize padding/cutting time length - if t0 is None: - self.t0 = 10000. - else: - self.t0 = t0 + self.t0 = 10000. # initialize whether to use gpu; FLR handles if this cannot be done self.use_gpu = use_gpu + # extrinsic params for LISA time conversion + self.beta = None + self.lamb = None + self.pol = None + def apply_polarization(self, hp, hc, polarization): """ Apply polarization rotation matrix. @@ -737,7 +760,7 @@ def apply_polarization(self, hp, hc, polarization): The cross polarization of the GW. polarization : float - The polarization angle of the GW in radians. + The SSB polarization angle of the GW in radians. Returns ------- @@ -756,8 +779,6 @@ def get_links(self, hp, hc, lamb, beta, polarization=0, reference_time=None, use_gpu=None): """ Project a radiation frame waveform to the LISA constellation. - This does not perform TDI; this only produces a detector frame - representation of the signal. Parameters ---------- @@ -790,22 +811,36 @@ def get_links(self, hp, hc, lamb, beta, polarization=0, reference_time=None, ndarray The waveform projected to the LISA laser links. """ - from fastlisaresponse import pyResponseTDI + if pyResponseTDI is None: + raise ImportError('FastLISAResponse required for LISA projection/TDI') - # get dt from waveform data (skipped if running project_wave) + # get dt from waveform data if self.dt is None: self.dt = hp.delta_t - # set waveform start time and cache time samples + # set waveform start time if specified if reference_time is not None: self.ref_time = reference_time - hp.start_time = self.ref_time - hc.start_time = self.ref_time + + # specify and cache sample times + start = self.ref_time + base_dur = hp.duration + if self.pad_data: + start -= self.t0 # start should correspond to input signal, not pad + hp.start_time = start + hc.start_time = start self.sample_times = hp.sample_times.numpy() - # rescale the orbital time series to match waveform - assert hp.duration + self.ref_time <= self.orbits.t_base[-1], ( - "Time of signal end is greater than end of orbital data") + # make sure signal still lies within orbit length + assert hp.duration + start <= self.orbits.t_base[-1], ( + "Time of signal end is greater than end of orbital data.") + if self.pad_data: + # specify that the padding is causing the issue + assert start >= self.orbits.t_base[0], ( + "Starting pad extends before start of orbital data. " + + "Consider decreasing t0 or increasing reference time.") + + # configure the orbit to match signal self.orbits.configure(t_arr=self.sample_times) # rotate GW from radiation frame to SSB using polarization angle @@ -826,7 +861,8 @@ def get_links(self, hp, hc, lamb, beta, polarization=0, reference_time=None, # convert to cupy if needed if use_gpu: if cupy is None: - raise ImportError('CuPy not imported but is required for GPU usage') + raise ImportError('CuPy not imported but is required for GPU usage. ' + + 'Ensure use_gpu = False if not using GPU.') else: wf = cupy.asarray(wf) @@ -849,7 +885,7 @@ def get_links(self, hp, hc, lamb, beta, polarization=0, reference_time=None, def project_wave(self, hp, hc, lamb, beta, polarization, reference_time=None, tdi=2, tdi_chan='AET', tdi_orbits=None, use_gpu=None, - pad_data=False, remove_garbage=True): + pad_data=False, remove_garbage=True, t0=None): """ Evaluate the TDI observables. @@ -909,26 +945,34 @@ def project_wave(self, hp, hc, lamb, beta, polarization, reference_time=None, edge data will be zeroed. If False, TDI channels will not be modified. Default True. + t0 : float (optional) + Time length in seconds to pad/cut from the start and end of + the data if pad_data/remove_garbage is True. Default 10000. + Returns ------- dict The TDI observables keyed by their corresponding TDI channel. """ + self.pad_data = pad_data + # get dt from waveform data self.dt = hp.delta_t - # index corresponding to time length t0 - self.pad_idx = int(self.t0/self.dt) + # get index corresponding to time length t0 + if t0 is not None: + self.t0 = t0 + if pad_data or remove_garbage: + pad_idx = int(self.t0/self.dt) # pad the data with zeros ### this assumes that the signal naturally tapers to zero ### this will not work with e.g. GBs or sinusoids - ### we should be tapering in general to prevent discontinuities if pad_data: - hp.prepend_zeros(self.pad_idx) - hp.append_zeros(self.pad_idx) - hc.prepend_zeros(self.pad_idx) - hc.append_zeros(self.pad_idx) + hp.prepend_zeros(pad_idx) + hp.append_zeros(pad_idx) + hc.prepend_zeros(pad_idx) + hc.append_zeros(pad_idx) # set use_gpu if use_gpu is None: @@ -968,18 +1012,17 @@ def project_wave(self, hp, hc, lamb, beta, polarization, reference_time=None, # processing tdi_dict = {} - slc = slice(self.pad_idx, -self.pad_idx) + slc = slice(pad_idx, -pad_idx) for i in range(len(tdi_chan)): - # add to dict tdi_dict[tdi_chan[i]] = tdi_obs[i] - # treat start and end gaps (if remove_garbage is not False) + # treat start and end gaps if remove_garbage: if remove_garbage == 'zero': # zero the edge data - tdi_dict[tdi_chan[i]][:self.pad_idx] = 0 - tdi_dict[tdi_chan[i]][-self.pad_idx:] = 0 + tdi_dict[tdi_chan[i]][:pad_idx] = 0 + tdi_dict[tdi_chan[i]][-pad_idx:] = 0 elif type(remove_garbage) == bool: # cut the edge data tdi_dict[tdi_chan[i]] = tdi_dict[tdi_chan[i]][slc] @@ -987,10 +1030,34 @@ def project_wave(self, hp, hc, lamb, beta, polarization, reference_time=None, # cut the sample times (only do once) self.sample_times = self.sample_times[slc] else: - raise ValueError('remove_garbage arg must be a bool or "zero"') + raise ValueError('remove_garbage arg must be a bool ' + + 'or "zero"') + + # unapply offset for GPS SSB times + self.sample_times -= self.offset + self.ref_time -= self.offset + + # save params for LISA time conversion + self.beta = beta + self.lamb = lamb + self.pol = polarization return tdi_dict + @property + def sample_times_LISA(self): + """ + Sample times of signal converted to LISA frame. + """ + params = [self.sample_times, self.beta, self.lamb, self.pol] + assert all(i is not None for i in params), ("Need to run project_wave for conversion") + + # calculate first time converted to LISA; shift sample times to match LISA start time + lisa_start, _, _, _ = ssb_to_lisa(self.sample_times[0], self.beta, self.lamb, + self.pol, t0=0) + diff = lisa_start - self.sample_times[0] + return self.sample_times + diff + def ppdets(ifos, separator=', '): """Pretty-print a list (or set) of detectors: return a string listing From 0e6c7d37565d13cf62fe45228ee3ecb00b3a37bf Mon Sep 17 00:00:00 2001 From: Alex Correia Date: Mon, 16 Sep 2024 10:06:44 -0400 Subject: [PATCH 19/24] add debug statements, adjust SSB to LISA time conversions, fix sky coord labels --- pycbc/detector.py | 79 +++++++++++++++++++++++++++++++---------------- 1 file changed, 52 insertions(+), 27 deletions(-) diff --git a/pycbc/detector.py b/pycbc/detector.py index bbdaf14df3e..5df010c295a 100644 --- a/pycbc/detector.py +++ b/pycbc/detector.py @@ -660,7 +660,7 @@ def overhead_antenna_pattern(right_ascension, declination, polarization): from lisatools.detector import ESAOrbits except ImportError: ESAOrbits = None - + try: from fastlisaresponse import pyResponseTDI except ImportError: @@ -789,10 +789,10 @@ def get_links(self, hp, hc, lamb, beta, polarization=0, reference_time=None, The cross polarization of the GW in the radiation frame. lamb : float - The ecliptic latitude of the source in the SSB frame. + The ecliptic longitude of the source in the SSB frame. beta : float - The ecliptic longitude of the source in the SSB frame. + The ecliptic latitude of the source in the SSB frame. polarization : float (optional) The polarization angle of the GW in radians. Default 0. @@ -820,7 +820,7 @@ def get_links(self, hp, hc, lamb, beta, polarization=0, reference_time=None, # set waveform start time if specified if reference_time is not None: - self.ref_time = reference_time + self.ref_time = reference_time + self.offset # specify and cache sample times start = self.ref_time @@ -830,6 +830,7 @@ def get_links(self, hp, hc, lamb, beta, polarization=0, reference_time=None, hp.start_time = start hc.start_time = start self.sample_times = hp.sample_times.numpy() + print('cached wf inputs and times') # make sure signal still lies within orbit length assert hp.duration + start <= self.orbits.t_base[-1], ( @@ -842,14 +843,17 @@ def get_links(self, hp, hc, lamb, beta, polarization=0, reference_time=None, # configure the orbit to match signal self.orbits.configure(t_arr=self.sample_times) + print('configured orbits') # rotate GW from radiation frame to SSB using polarization angle hp, hc = self.apply_polarization(hp, hc, polarization) + print('applied polarization') # format wf to hp + i*hc hp = hp.numpy() hc = hc.numpy() wf = hp + 1j*hc + print('converted to numpy') # save length of wf self.n = len(wf) @@ -868,14 +872,17 @@ def get_links(self, hp, hc, lamb, beta, polarization=0, reference_time=None, if self.response_init is None: # initialize the class + print('fresh init') self.response_init = pyResponseTDI(1/self.dt, self.n, orbits=self.orbits, use_gpu=use_gpu) else: # update params in the initialized class + print('update init') self.response_init.sampling_frequency = 1/self.dt self.response_init.num_pts = self.n self.response_init.orbits = self.orbits self.response_init.use_gpu = use_gpu + print('response initialized') # project the signal self.response_init.get_projections(wf, lamb, beta, t0=self.t0) @@ -906,10 +913,10 @@ def project_wave(self, hp, hc, lamb, beta, polarization, reference_time=None, The cross polarization of the GW in the radiation frame. lamb : float - The ecliptic latitude in the SSB frame. + The ecliptic longitude in the SSB frame. beta : float - The ecliptic longitude in the SSB frame. + The ecliptic latitude in the SSB frame. polarization : float The polarization angle of the GW in radians. @@ -959,10 +966,16 @@ def project_wave(self, hp, hc, lamb, beta, polarization, reference_time=None, # get dt from waveform data self.dt = hp.delta_t + # save params for LISA time conversion + self.beta = beta + self.lamb = lamb + self.pol = polarization + # get index corresponding to time length t0 if t0 is not None: self.t0 = t0 if pad_data or remove_garbage: + global pad_idx pad_idx = int(self.t0/self.dt) # pad the data with zeros @@ -981,6 +994,7 @@ def project_wave(self, hp, hc, lamb, beta, polarization, reference_time=None, # generate the Doppler time series self.get_links(hp, hc, lamb, beta, polarization=polarization, reference_time=reference_time, use_gpu=use_gpu) + print('get links') # set TDI configuration (let FLR handle if not 1 or 2) if tdi == 1: @@ -1009,13 +1023,19 @@ def project_wave(self, hp, hc, lamb, beta, polarization, reference_time=None, # generate the TDI channels tdi_obs = self.response_init.get_tdi_delays() + print('tdi complete') # processing tdi_dict = {} - slc = slice(pad_idx, -pad_idx) + self.sample_times -= self.offset + self.ref_time -= self.offset + print('start preprocessing') for i in range(len(tdi_chan)): - tdi_dict[tdi_chan[i]] = tdi_obs[i] + # save as TimeSeries with LISA frame times + tdi_dict[tdi_chan[i]] = TimeSeries(tdi_obs[i], delta_t=self.dt, + epoch=self.ref_time_LISA) + print(f'saved {i} to TimeSeries') # treat start and end gaps if remove_garbage: @@ -1025,37 +1045,42 @@ def project_wave(self, hp, hc, lamb, beta, polarization, reference_time=None, tdi_dict[tdi_chan[i]][-pad_idx:] = 0 elif type(remove_garbage) == bool: # cut the edge data + slc = slice(pad_idx, -pad_idx) tdi_dict[tdi_chan[i]] = tdi_dict[tdi_chan[i]][slc] if i == 0: - # cut the sample times (only do once) - self.sample_times = self.sample_times[slc] + # update sample times once + self.sample_times = tdi_dict[tdi_chan[i]].sample_times else: - raise ValueError('remove_garbage arg must be a bool ' + + raise ValueError('remove_garbage arg must be a bool ' + 'or "zero"') - - # unapply offset for GPS SSB times - self.sample_times -= self.offset - self.ref_time -= self.offset - - # save params for LISA time conversion - self.beta = beta - self.lamb = lamb - self.pol = polarization + print(f'finished postprocessing {i}') return tdi_dict @property - def sample_times_LISA(self): + def ref_time_LISA(self): """ - Sample times of signal converted to LISA frame. + Signal start time converted to LISA frame. """ - params = [self.sample_times, self.beta, self.lamb, self.pol] + params = [self.ref_time, self.lamb, self.beta, self.pol] assert all(i is not None for i in params), ("Need to run project_wave for conversion") - # calculate first time converted to LISA; shift sample times to match LISA start time - lisa_start, _, _, _ = ssb_to_lisa(self.sample_times[0], self.beta, self.lamb, - self.pol, t0=0) - diff = lisa_start - self.sample_times[0] + # convert ref time to LISA + ssb_start = self.ref_time.gpsSeconds + 1e-9*self.ref_time.gpsNanoSeconds + lisa_start, _, _, _ = ssb_to_lisa(t_ssb = ssb_start, + longitude_ssb = self.lamb, + latitude_ssb = self.beta, + polarization_ssb = self.pol, + t0=self.offset) + + return lisa_start + + @property + def sample_times_LISA(self): + """ + Signal sample times converted to LISA frame. + """ + diff = self.ref_time_LISA - self.ref_time return self.sample_times + diff From 8750bf05a18216b05e0ee8dbe78b6a1ef4356ca0 Mon Sep 17 00:00:00 2001 From: alcorrei Date: Mon, 16 Sep 2024 17:47:05 -0400 Subject: [PATCH 20/24] ref time accepts anything castable to float, tweak how TDI chans are processed and saved --- pycbc/detector.py | 46 ++++++++++++++++++++++++++++------------------ 1 file changed, 28 insertions(+), 18 deletions(-) diff --git a/pycbc/detector.py b/pycbc/detector.py index 5df010c295a..65aed4a22d7 100644 --- a/pycbc/detector.py +++ b/pycbc/detector.py @@ -679,9 +679,10 @@ def __init__(self, detector='LISA', reference_time=None, orbits=None, String specifying space-borne detector. Currently only accepts 'LISA', which is the default setting. - reference_time: float (optional) - The GPS start time of signal in the SSB frame. Default to start time of - orbits input. + reference_time: float-like (optional) + The start time of signal in the SSB frame. Accepts any type that is + castable to float (e.g. LIGOTimeGPS). Default to start time of orbits + input. orbits: lisatools.detector.Orbits (optional) Orbital information to pass into pyResponseTDI. Default @@ -724,14 +725,17 @@ def __init__(self, detector='LISA', reference_time=None, orbits=None, if reference_time is None: self.ref_time = orbit_start + self.offset else: - reference_time += self.offset - assert (reference_time >= orbit_start) and (reference_time <= orbit_end), ( + try: + rt = float(reference_time) + except ValueError: + raise ValueError('reference time input must be castable to float') + rt += self.offset + assert (rt >= orbit_start) and (rt <= orbit_end), ( "Reference time is not in time domain of orbital data") - self.ref_time = reference_time + self.ref_time = rt # allocate caches self.dt = None - self.n = None self.pad_data = False # don't pad by default if only projecting self.sample_times = None self.response_init = None @@ -820,7 +824,11 @@ def get_links(self, hp, hc, lamb, beta, polarization=0, reference_time=None, # set waveform start time if specified if reference_time is not None: - self.ref_time = reference_time + self.offset + try: + rt = float(reference_time) + except ValueError: + raise ValueError('reference time input must be castable to float') + self.ref_time = rt + self.offset # specify and cache sample times start = self.ref_time @@ -856,7 +864,7 @@ def get_links(self, hp, hc, lamb, beta, polarization=0, reference_time=None, print('converted to numpy') # save length of wf - self.n = len(wf) + n = len(wf) # set use_gpu to class input if not specified if use_gpu is None: @@ -873,13 +881,13 @@ def get_links(self, hp, hc, lamb, beta, polarization=0, reference_time=None, if self.response_init is None: # initialize the class print('fresh init') - self.response_init = pyResponseTDI(1/self.dt, self.n, orbits=self.orbits, + self.response_init = pyResponseTDI(1/self.dt, n, orbits=self.orbits, use_gpu=use_gpu) else: # update params in the initialized class print('update init') self.response_init.sampling_frequency = 1/self.dt - self.response_init.num_pts = self.n + self.response_init.num_pts = n self.response_init.orbits = self.orbits self.response_init.use_gpu = use_gpu print('response initialized') @@ -1032,10 +1040,8 @@ def project_wave(self, hp, hc, lamb, beta, polarization, reference_time=None, print('start preprocessing') for i in range(len(tdi_chan)): - # save as TimeSeries with LISA frame times - tdi_dict[tdi_chan[i]] = TimeSeries(tdi_obs[i], delta_t=self.dt, - epoch=self.ref_time_LISA) - print(f'saved {i} to TimeSeries') + # save as numpy arrays + tdi_dict[tdi_chan[i]] = tdi_obs[i] # treat start and end gaps if remove_garbage: @@ -1049,11 +1055,16 @@ def project_wave(self, hp, hc, lamb, beta, polarization, reference_time=None, tdi_dict[tdi_chan[i]] = tdi_dict[tdi_chan[i]][slc] if i == 0: # update sample times once - self.sample_times = tdi_dict[tdi_chan[i]].sample_times + self.sample_times = self.sample_times[slc] else: raise ValueError('remove_garbage arg must be a bool ' + 'or "zero"') print(f'finished postprocessing {i}') + + # save as TimeSeries with LISA frame times + tdi_dict[tdi_chan[i]] = TimeSeries(tdi_dict[tdi_chan[i]], delta_t=self.dt, + epoch=self.ref_time_LISA) + print(f'saved {i} to TimeSeries') return tdi_dict @@ -1066,8 +1077,7 @@ def ref_time_LISA(self): assert all(i is not None for i in params), ("Need to run project_wave for conversion") # convert ref time to LISA - ssb_start = self.ref_time.gpsSeconds + 1e-9*self.ref_time.gpsNanoSeconds - lisa_start, _, _, _ = ssb_to_lisa(t_ssb = ssb_start, + lisa_start, _, _, _ = ssb_to_lisa(t_ssb = self.ref_time, longitude_ssb = self.lamb, latitude_ssb = self.beta, polarization_ssb = self.pol, From 880ed4fa2fdc82ef1911b782013377946cfe386f Mon Sep 17 00:00:00 2001 From: alcorrei Date: Tue, 17 Sep 2024 16:30:41 -0400 Subject: [PATCH 21/24] start reworking reference time to replicate bbhx (i.e. reference time inputs are no longer necessarily start time) --- pycbc/detector.py | 116 +++++++++++++++++++++++++--------------------- 1 file changed, 63 insertions(+), 53 deletions(-) diff --git a/pycbc/detector.py b/pycbc/detector.py index 65aed4a22d7..78ddf59a641 100644 --- a/pycbc/detector.py +++ b/pycbc/detector.py @@ -661,11 +661,6 @@ def overhead_antenna_pattern(right_ascension, declination, polarization): except ImportError: ESAOrbits = None -try: - from fastlisaresponse import pyResponseTDI -except ImportError: - pyResponseTDI = None - class LISA_detector(object): """ LISA-like GW detector. Applies detector response from FastLISAResponse. @@ -680,9 +675,9 @@ def __init__(self, detector='LISA', reference_time=None, orbits=None, which is the default setting. reference_time: float-like (optional) - The start time of signal in the SSB frame. Accepts any type that is - castable to float (e.g. LIGOTimeGPS). Default to start time of orbits - input. + The reference time of signal in the SSB frame. Accepts any type that is + castable to float (e.g. LIGOTimeGPS). By default, the reference time is + set to the midpoint of the time series. orbits: lisatools.detector.Orbits (optional) Orbital information to pass into pyResponseTDI. Default @@ -721,24 +716,13 @@ def __init__(self, detector='LISA', reference_time=None, orbits=None, else: self.offset = 0. - # specify and cache the start time - if reference_time is None: - self.ref_time = orbit_start + self.offset - else: - try: - rt = float(reference_time) - except ValueError: - raise ValueError('reference time input must be castable to float') - rt += self.offset - assert (rt >= orbit_start) and (rt <= orbit_end), ( - "Reference time is not in time domain of orbital data") - self.ref_time = rt - # allocate caches self.dt = None self.pad_data = False # don't pad by default if only projecting self.sample_times = None self.response_init = None + self.ref_time = reference_time + self.start_time = None # initialize padding/cutting time length self.t0 = 10000. @@ -802,8 +786,8 @@ def get_links(self, hp, hc, lamb, beta, polarization=0, reference_time=None, The polarization angle of the GW in radians. Default 0. reference_time : float (optional) - The GPS start time of the GW signal in the SSB frame. Default to - input on class initialization. + The reference time of the GW signal in the SSB frame. Default + behavior places start time of the signal at GPS time t=0. use_gpu : bool (optional) Flag whether to use GPU support. Default to class input. @@ -815,41 +799,53 @@ def get_links(self, hp, hc, lamb, beta, polarization=0, reference_time=None, ndarray The waveform projected to the LISA laser links. """ - if pyResponseTDI is None: + try: + from fastlisaresponse import pyResponseTDI + except ImportError: raise ImportError('FastLISAResponse required for LISA projection/TDI') # get dt from waveform data if self.dt is None: self.dt = hp.delta_t - # set waveform start time if specified - if reference_time is not None: - try: - rt = float(reference_time) - except ValueError: - raise ValueError('reference time input must be castable to float') - self.ref_time = rt + self.offset - - # specify and cache sample times - start = self.ref_time - base_dur = hp.duration + if reference_time is None: + reference_time = self.ref_time + + # get start time of waveform + if reference_time is None: + start_time = self.offset # t = 0 scaled by offset + else: + base_dur = hp.duration + if self.pad_data: + base_dur -= 2*self.t0 # subtract off pads for base wf length + reference_time += self.offset + self.ref_time = float(reference_time) + ### this assumes the wf is generated such that tref = 0 + ### is this generally true? + ### if not we need to specify the time value of the input wf tref + start_time = float(reference_time + hp.start_time) # start time of unpadded wf + if self.pad_data: - start -= self.t0 # start should correspond to input signal, not pad - hp.start_time = start - hc.start_time = start - self.sample_times = hp.sample_times.numpy() - print('cached wf inputs and times') + start_time -= self.t0 # start time of padded wf + # set start times; save reference time as midpoint if not specified + hp.start_time = start_time + hc.start_time = start_time + self.start_time = start_time + if reference_time is None: + self.ref_time = float((hp.end_time - hp.start_time)/2) + # make sure signal still lies within orbit length - assert hp.duration + start <= self.orbits.t_base[-1], ( + assert hp.duration + start_time <= self.orbits.t_base[-1], ( "Time of signal end is greater than end of orbital data.") if self.pad_data: # specify that the padding is causing the issue - assert start >= self.orbits.t_base[0], ( + assert start_time >= self.orbits.t_base[0], ( "Starting pad extends before start of orbital data. " + "Consider decreasing t0 or increasing reference time.") # configure the orbit to match signal + self.sample_times = hp.sample_times.numpy() self.orbits.configure(t_arr=self.sample_times) print('configured orbits') @@ -930,8 +926,8 @@ def project_wave(self, hp, hc, lamb, beta, polarization, reference_time=None, The polarization angle of the GW in radians. reference_time : float (optional) - The GPS start time of the GW signal in the SSB frame. Default to - value in input signals hp and hc. + The reference time of the GW signal in the SSB frame. Default + behavior places start time of the signal at GPS time t=0. tdi : int (optional) TDI channel configuration. Accepts 1 for 1st generation TDI or @@ -970,6 +966,8 @@ def project_wave(self, hp, hc, lamb, beta, polarization, reference_time=None, The TDI observables keyed by their corresponding TDI channel. """ self.pad_data = pad_data + if reference_time is None: + reference_time = self.ref_time # get dt from waveform data self.dt = hp.delta_t @@ -1035,8 +1033,8 @@ def project_wave(self, hp, hc, lamb, beta, polarization, reference_time=None, # processing tdi_dict = {} - self.sample_times -= self.offset self.ref_time -= self.offset + self.start_time -= self.offset print('start preprocessing') for i in range(len(tdi_chan)): @@ -1063,7 +1061,11 @@ def project_wave(self, hp, hc, lamb, beta, polarization, reference_time=None, # save as TimeSeries with LISA frame times tdi_dict[tdi_chan[i]] = TimeSeries(tdi_dict[tdi_chan[i]], delta_t=self.dt, - epoch=self.ref_time_LISA) + epoch=self.start_time_LISA) + if i == 0: + # reset sample times to LISA + self.sample_times = tdi_dict[tdi_chan[i]].sample_times.numpy() + print(f'saved {i} to TimeSeries') return tdi_dict @@ -1071,28 +1073,36 @@ def project_wave(self, hp, hc, lamb, beta, polarization, reference_time=None, @property def ref_time_LISA(self): """ - Signal start time converted to LISA frame. + Reference time converted to LISA frame. """ params = [self.ref_time, self.lamb, self.beta, self.pol] assert all(i is not None for i in params), ("Need to run project_wave for conversion") # convert ref time to LISA - lisa_start, _, _, _ = ssb_to_lisa(t_ssb = self.ref_time, + lisa_ref, _, _, _ = ssb_to_lisa(t_ssb = self.ref_time, longitude_ssb = self.lamb, latitude_ssb = self.beta, polarization_ssb = self.pol, t0=self.offset) - return lisa_start + return lisa_ref @property - def sample_times_LISA(self): + def start_time_LISA(self): """ - Signal sample times converted to LISA frame. + Start time converted to LISA frame. """ - diff = self.ref_time_LISA - self.ref_time - return self.sample_times + diff + params = [self.start_time, self.lamb, self.beta, self.pol] + assert all(i is not None for i in params), ("Need to run project_wave for conversion") + # convert ref time to LISA + lisa_start, _, _, _ = ssb_to_lisa(t_ssb = self.start_time, + longitude_ssb = self.lamb, + latitude_ssb = self.beta, + polarization_ssb = self.pol, + t0=self.offset) + + return lisa_start def ppdets(ifos, separator=', '): """Pretty-print a list (or set) of detectors: return a string listing From 859ad1991f82fc861debd38ea550c77a83cd17d3 Mon Sep 17 00:00:00 2001 From: Alex Correia Date: Wed, 18 Sep 2024 17:11:23 -0400 Subject: [PATCH 22/24] fix reference time to correctly calculate signal start time; fix sample times conversion to be based on reference time --- pycbc/detector.py | 76 +++++++++++++++++------------------------------ 1 file changed, 28 insertions(+), 48 deletions(-) diff --git a/pycbc/detector.py b/pycbc/detector.py index 78ddf59a641..222e702fc0a 100644 --- a/pycbc/detector.py +++ b/pycbc/detector.py @@ -804,37 +804,36 @@ def get_links(self, hp, hc, lamb, beta, polarization=0, reference_time=None, except ImportError: raise ImportError('FastLISAResponse required for LISA projection/TDI') + # save params for LISA time conversion + self.lamb = lamb + self.beta = beta + self.pol = polarization + # get dt from waveform data if self.dt is None: self.dt = hp.delta_t - if reference_time is None: - reference_time = self.ref_time - # get start time of waveform if reference_time is None: start_time = self.offset # t = 0 scaled by offset else: - base_dur = hp.duration - if self.pad_data: - base_dur -= 2*self.t0 # subtract off pads for base wf length reference_time += self.offset self.ref_time = float(reference_time) - ### this assumes the wf is generated such that tref = 0 - ### is this generally true? - ### if not we need to specify the time value of the input wf tref - start_time = float(reference_time + hp.start_time) # start time of unpadded wf - - if self.pad_data: - start_time -= self.t0 # start time of padded wf + # assume tref = 0 in input waveforms (default from get_td_waveform) + start_time = float(self.ref_time + hp.start_time) - # set start times; save reference time as midpoint if not specified + # set start times hp.start_time = start_time hc.start_time = start_time - self.start_time = start_time + if self.pad_data: + self.start_time = start_time + self.t0 # signal starts t0 after pad start + else: + self.start_time = start_time + + # save reference time as midpoint if not specified if reference_time is None: - self.ref_time = float((hp.end_time - hp.start_time)/2) - + self.ref_time = float((hp.end_time + hp.start_time)/2) + # make sure signal still lies within orbit length assert hp.duration + start_time <= self.orbits.t_base[-1], ( "Time of signal end is greater than end of orbital data.") @@ -972,22 +971,17 @@ def project_wave(self, hp, hc, lamb, beta, polarization, reference_time=None, # get dt from waveform data self.dt = hp.delta_t - # save params for LISA time conversion - self.beta = beta - self.lamb = lamb - self.pol = polarization - # get index corresponding to time length t0 if t0 is not None: self.t0 = t0 - if pad_data or remove_garbage: + if self.pad_data or remove_garbage: global pad_idx pad_idx = int(self.t0/self.dt) # pad the data with zeros ### this assumes that the signal naturally tapers to zero ### this will not work with e.g. GBs or sinusoids - if pad_data: + if self.pad_data: hp.prepend_zeros(pad_idx) hp.append_zeros(pad_idx) hc.prepend_zeros(pad_idx) @@ -997,6 +991,10 @@ def project_wave(self, hp, hc, lamb, beta, polarization, reference_time=None, if use_gpu is None: use_gpu = self.use_gpu + # get reference time from class if not supplied + if reference_time is None: + reference_time = self.ref_time + # generate the Doppler time series self.get_links(hp, hc, lamb, beta, polarization=polarization, reference_time=reference_time, use_gpu=use_gpu) @@ -1051,21 +1049,18 @@ def project_wave(self, hp, hc, lamb, beta, polarization, reference_time=None, # cut the edge data slc = slice(pad_idx, -pad_idx) tdi_dict[tdi_chan[i]] = tdi_dict[tdi_chan[i]][slc] - if i == 0: - # update sample times once - self.sample_times = self.sample_times[slc] else: raise ValueError('remove_garbage arg must be a bool ' + 'or "zero"') print(f'finished postprocessing {i}') - + # save as TimeSeries with LISA frame times + loc = self.ref_time - self.start_time tdi_dict[tdi_chan[i]] = TimeSeries(tdi_dict[tdi_chan[i]], delta_t=self.dt, - epoch=self.start_time_LISA) + epoch=self.ref_time_LISA - loc) if i == 0: - # reset sample times to LISA + # convert sample times to LISA self.sample_times = tdi_dict[tdi_chan[i]].sample_times.numpy() - print(f'saved {i} to TimeSeries') return tdi_dict @@ -1076,7 +1071,8 @@ def ref_time_LISA(self): Reference time converted to LISA frame. """ params = [self.ref_time, self.lamb, self.beta, self.pol] - assert all(i is not None for i in params), ("Need to run project_wave for conversion") + assert all(i is not None for i in params), ( + "Need to run get_links (or project_wave) for conversion") # convert ref time to LISA lisa_ref, _, _, _ = ssb_to_lisa(t_ssb = self.ref_time, @@ -1087,22 +1083,6 @@ def ref_time_LISA(self): return lisa_ref - @property - def start_time_LISA(self): - """ - Start time converted to LISA frame. - """ - params = [self.start_time, self.lamb, self.beta, self.pol] - assert all(i is not None for i in params), ("Need to run project_wave for conversion") - - # convert ref time to LISA - lisa_start, _, _, _ = ssb_to_lisa(t_ssb = self.start_time, - longitude_ssb = self.lamb, - latitude_ssb = self.beta, - polarization_ssb = self.pol, - t0=self.offset) - - return lisa_start def ppdets(ifos, separator=', '): """Pretty-print a list (or set) of detectors: return a string listing From cade855d96973c461c00d54701f6a0b15bf27239 Mon Sep 17 00:00:00 2001 From: Alex Correia Date: Thu, 19 Sep 2024 19:14:51 -0400 Subject: [PATCH 23/24] fixes to start time calculation --- pycbc/detector.py | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/pycbc/detector.py b/pycbc/detector.py index 222e702fc0a..f438ae4d433 100644 --- a/pycbc/detector.py +++ b/pycbc/detector.py @@ -649,7 +649,7 @@ def overhead_antenna_pattern(right_ascension, declination, polarization): """ LISA class """ -from pycbc.coordinates.space import ssb_to_lisa, TIME_OFFSET_20_DEGREES +from pycbc.coordinates.space import ssb_to_lisa, TIME_OFFSET_20_DEGREES, t_ssb_from_t_lisa try: import cupy @@ -820,11 +820,14 @@ def get_links(self, hp, hc, lamb, beta, polarization=0, reference_time=None, reference_time += self.offset self.ref_time = float(reference_time) # assume tref = 0 in input waveforms (default from get_td_waveform) - start_time = float(self.ref_time + hp.start_time) + start_time_LISA = float(self.ref_time_LISA + hp.start_time) + # convert the LISA start time to SSB + start_time = t_ssb_from_t_lisa(start_time_LISA, self.lamb, self.beta, t0=0) # set start times hp.start_time = start_time hc.start_time = start_time + if self.pad_data: self.start_time = start_time + self.t0 # signal starts t0 after pad start else: @@ -1076,10 +1079,10 @@ def ref_time_LISA(self): # convert ref time to LISA lisa_ref, _, _, _ = ssb_to_lisa(t_ssb = self.ref_time, - longitude_ssb = self.lamb, - latitude_ssb = self.beta, - polarization_ssb = self.pol, - t0=self.offset) + longitude_ssb = self.lamb, + latitude_ssb = self.beta, + polarization_ssb = self.pol, + t0=self.offset) return lisa_ref From 68b61c78571f7bfdbea5aea9a95dfc417341b789 Mon Sep 17 00:00:00 2001 From: Alex Correia Date: Thu, 26 Sep 2024 10:13:05 -0400 Subject: [PATCH 24/24] set TDI chans to SSB times, clean up unnecessary attributes --- pycbc/detector.py | 143 ++++++++++++++++------------------------------ 1 file changed, 49 insertions(+), 94 deletions(-) diff --git a/pycbc/detector.py b/pycbc/detector.py index f438ae4d433..81ba7a755a8 100644 --- a/pycbc/detector.py +++ b/pycbc/detector.py @@ -649,21 +649,17 @@ def overhead_antenna_pattern(right_ascension, declination, polarization): """ LISA class """ -from pycbc.coordinates.space import ssb_to_lisa, TIME_OFFSET_20_DEGREES, t_ssb_from_t_lisa +from pycbc.coordinates.space import TIME_OFFSET_20_DEGREES try: import cupy except ImportError: cupy = None -try: - from lisatools.detector import ESAOrbits -except ImportError: - ESAOrbits = None - class LISA_detector(object): """ - LISA-like GW detector. Applies detector response from FastLISAResponse. + LISA-like GW detector. Applies detector response from FastLISAResponse + (arXiv:2204.06633). """ def __init__(self, detector='LISA', reference_time=None, orbits=None, use_gpu=False, apply_offset=False, offset=TIME_OFFSET_20_DEGREES): @@ -681,19 +677,18 @@ def __init__(self, detector='LISA', reference_time=None, orbits=None, orbits: lisatools.detector.Orbits (optional) Orbital information to pass into pyResponseTDI. Default - lisatools.detector.ESAOrbits. + lisatools.detector.EqualArmlengthOrbits. use_gpu : bool (optional) Specify whether to run class on GPU support via CuPy. Default False. apply_offset : bool (optional) Specify whether to add a time offset to input times. If True, - times are treated as GPS times, to which an offset is added to ensure - LISA is ~20 degrees behind the Earth at t=0 (the LDC convention for the - LISA start date). If False, no offset is applied. Default False. + an offset is added to ensure LISA is oriented correctly relative to the + Earth at t=0. If False, no offset is applied. Default False. offset : float (optional) - Time offset in seconds to apply to GPS times if apply_offset = True. + Time offset in seconds to apply to input times if apply_offset = True. Default pycbc.coordinates.space.TIME_OFFSET_20_DEGREES (places LISA ~20 deg behind Earth). """ @@ -702,10 +697,12 @@ def __init__(self, detector='LISA', reference_time=None, orbits=None, # intialize orbit information if orbits is None: - # set ESAOrbits as default; raise error if ESAOrbits cannot be imported - if ESAOrbits is None: + try: + # use equal armlengths by default + from lisatools.detector import EqualArmlengthOrbits + orbits = EqualArmlengthOrbits() + except ImportError: raise ImportError('LISAanalysistools required for inputting orbital data') - orbits = ESAOrbits() self.orbits = orbits orbit_start = self.orbits.t_base[0] orbit_end = self.orbits.t_base[-1] @@ -721,6 +718,8 @@ def __init__(self, detector='LISA', reference_time=None, orbits=None, self.pad_data = False # don't pad by default if only projecting self.sample_times = None self.response_init = None + if reference_time is not None: + reference_time = float(reference_time) self.ref_time = reference_time self.start_time = None @@ -730,11 +729,6 @@ def __init__(self, detector='LISA', reference_time=None, orbits=None, # initialize whether to use gpu; FLR handles if this cannot be done self.use_gpu = use_gpu - # extrinsic params for LISA time conversion - self.beta = None - self.lamb = None - self.pol = None - def apply_polarization(self, hp, hc, polarization): """ Apply polarization rotation matrix. @@ -752,7 +746,7 @@ def apply_polarization(self, hp, hc, polarization): Returns ------- - (array, array) + (pycbc.types.TimeSeries, pycbc.types.TimeSeries) The plus and cross polarizations of the GW rotated by the polarization angle. """ cphi = cos(2*polarization) @@ -762,7 +756,7 @@ def apply_polarization(self, hp, hc, polarization): hc_ssb = hp*sphi + hc*cphi return hp_ssb, hc_ssb - + def get_links(self, hp, hc, lamb, beta, polarization=0, reference_time=None, use_gpu=None): """ @@ -797,54 +791,23 @@ def get_links(self, hp, hc, lamb, beta, polarization=0, reference_time=None, Returns ------- ndarray - The waveform projected to the LISA laser links. + The waveform projected to the LISA laser links. Shape is (6, N) + for input waveforms with N total samples. """ try: from fastlisaresponse import pyResponseTDI except ImportError: raise ImportError('FastLISAResponse required for LISA projection/TDI') - # save params for LISA time conversion - self.lamb = lamb - self.beta = beta - self.pol = polarization - # get dt from waveform data if self.dt is None: self.dt = hp.delta_t - # get start time of waveform - if reference_time is None: - start_time = self.offset # t = 0 scaled by offset - else: - reference_time += self.offset - self.ref_time = float(reference_time) - # assume tref = 0 in input waveforms (default from get_td_waveform) - start_time_LISA = float(self.ref_time_LISA + hp.start_time) - # convert the LISA start time to SSB - start_time = t_ssb_from_t_lisa(start_time_LISA, self.lamb, self.beta, t0=0) - - # set start times - hp.start_time = start_time - hc.start_time = start_time - - if self.pad_data: - self.start_time = start_time + self.t0 # signal starts t0 after pad start - else: - self.start_time = start_time - - # save reference time as midpoint if not specified - if reference_time is None: - self.ref_time = float((hp.end_time + hp.start_time)/2) - # make sure signal still lies within orbit length - assert hp.duration + start_time <= self.orbits.t_base[-1], ( + assert hp.duration + hp.start_time <= self.orbits.t_base[-1], ( "Time of signal end is greater than end of orbital data.") - if self.pad_data: - # specify that the padding is causing the issue - assert start_time >= self.orbits.t_base[0], ( - "Starting pad extends before start of orbital data. " + - "Consider decreasing t0 or increasing reference time.") + assert hp.start_time >= self.orbits.t_base[0], ( + "Time of signal start is less than start of orbital data.") # configure the orbit to match signal self.sample_times = hp.sample_times.numpy() @@ -964,12 +927,11 @@ def project_wave(self, hp, hc, lamb, beta, polarization, reference_time=None, Returns ------- - dict - The TDI observables keyed by their corresponding TDI channel. + dict ({str: pycbc.types.TimeSeries}) + The TDI observables as TimeSeries objects keyed by their + corresponding TDI channel name. """ self.pad_data = pad_data - if reference_time is None: - reference_time = self.ref_time # get dt from waveform data self.dt = hp.delta_t @@ -981,7 +943,23 @@ def project_wave(self, hp, hc, lamb, beta, polarization, reference_time=None, global pad_idx pad_idx = int(self.t0/self.dt) - # pad the data with zeros + # get reference time from class + if reference_time is None: + if self.ref_time is None: + # take ref time as midpoint of signal, start time as t=0 + self.ref_time = float(hp.duration/2) + reference_time = self.ref_time + + # get times corresponding to unpadded data + self.ref_time = reference_time + # assume tref = 0 in input waveforms (default from get_td_waveform) + self.start_time = float(self.ref_time + hp.start_time) + + # apply times to wfs + hp.start_time = self.start_time + self.offset + hc.start_time = self.start_time + self.offset + + # pad the data with zeros in the SSB frame ### this assumes that the signal naturally tapers to zero ### this will not work with e.g. GBs or sinusoids if self.pad_data: @@ -989,15 +967,13 @@ def project_wave(self, hp, hc, lamb, beta, polarization, reference_time=None, hp.append_zeros(pad_idx) hc.prepend_zeros(pad_idx) hc.append_zeros(pad_idx) - + # set use_gpu if use_gpu is None: use_gpu = self.use_gpu - # get reference time from class if not supplied - if reference_time is None: - reference_time = self.ref_time - + self.hp = hp + # generate the Doppler time series self.get_links(hp, hc, lamb, beta, polarization=polarization, reference_time=reference_time, use_gpu=use_gpu) @@ -1034,8 +1010,6 @@ def project_wave(self, hp, hc, lamb, beta, polarization, reference_time=None, # processing tdi_dict = {} - self.ref_time -= self.offset - self.start_time -= self.offset print('start preprocessing') for i in range(len(tdi_chan)): @@ -1057,35 +1031,16 @@ def project_wave(self, hp, hc, lamb, beta, polarization, reference_time=None, 'or "zero"') print(f'finished postprocessing {i}') - # save as TimeSeries with LISA frame times - loc = self.ref_time - self.start_time + # save as TimeSeries with SSB times tdi_dict[tdi_chan[i]] = TimeSeries(tdi_dict[tdi_chan[i]], delta_t=self.dt, - epoch=self.ref_time_LISA - loc) - if i == 0: - # convert sample times to LISA - self.sample_times = tdi_dict[tdi_chan[i]].sample_times.numpy() + epoch=self.start_time) + if pad_data and (not remove_garbage or remove_garbage == 'zero'): + # scale the start since the pads haven't been removed + tdi_dict[tdi_chan[i]].start_time -= self.t0 print(f'saved {i} to TimeSeries') return tdi_dict - @property - def ref_time_LISA(self): - """ - Reference time converted to LISA frame. - """ - params = [self.ref_time, self.lamb, self.beta, self.pol] - assert all(i is not None for i in params), ( - "Need to run get_links (or project_wave) for conversion") - - # convert ref time to LISA - lisa_ref, _, _, _ = ssb_to_lisa(t_ssb = self.ref_time, - longitude_ssb = self.lamb, - latitude_ssb = self.beta, - polarization_ssb = self.pol, - t0=self.offset) - - return lisa_ref - def ppdets(ifos, separator=', '): """Pretty-print a list (or set) of detectors: return a string listing