Skip to content

Commit

Permalink
add controls for zero padding data to project_wave; add controls for …
Browse files Browse the repository at this point in the history 
…specifying orbit start times for numerical data
  • Loading branch information
@acorreia61201
acorreia61201 committed Sep 5, 2024
1 parent c1123ff commit 0d3f8e6
Showing 1 changed file with 87 additions and 44 deletions.
131 changes: 87 additions & 44 deletions pycbc/detector.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
----------
Expand All @@ -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.
Expand All @@ -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):
Expand Down Expand Up @@ -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.
Expand All @@ -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:
Expand All @@ -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:
Expand All @@ -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

Expand All @@ -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
Expand All @@ -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.
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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
Expand Down

0 comments on commit 0d3f8e6

Please sign in to comment.