Skip to content

Commit

Permalink
Draft fix for followup indexing error
Browse files Browse the repository at this point in the history
  • Loading branch information
titodalcanton committed Mar 26, 2024
1 parent b95b260 commit 037f5cd
Show file tree
Hide file tree
Showing 2 changed files with 78 additions and 22 deletions.
41 changes: 36 additions & 5 deletions pycbc/filter/matchedfilter.py
Original file line number Diff line number Diff line change
Expand Up @@ -1831,16 +1831,42 @@ def followup_event_significance(ifo, data_reader, bank,
onsource_end += coinc_threshold

# Calculate how much time needed to calculate significance
trim_pad = (data_reader.trim_padding * data_reader.strain.delta_t)
bdur = int(lookback + 2.0 * trim_pad + length_in_time)
if bdur > data_reader.strain.duration * .75:
bdur = data_reader.strain.duration * .75
trim_pad = data_reader.trim_padding * data_reader.strain.delta_t
bdur = int(lookback + 2 * trim_pad + length_in_time)
bdur_samples = bank.round_up(bdur * bank.sample_rate)
max_safe_bdur_samples = int(
data_reader.strain.duration * .9 * bank.sample_rate
)
if bdur_samples > max_safe_bdur_samples:
bdur_samples = max_safe_bdur_samples
new_lookback = bdur_samples / bank.sample_rate - (2 * trim_pad + length_in_time)
if new_lookback > 0:
logging.warning(
'Strain buffer too short for a lookback time of %f s, '
'reducing lookback to %f s',
lookback,
new_lookback
)
else:
logging.error(
'Strain buffer too short to compute the followup SNR time '
'series for template %d, will not use %s for followup',
template_id,
ifo
)
return None
bdur = bdur_samples / bank.sample_rate

# Require all strain be valid within lookback time
if data_reader.state is not None:
state_start_time = data_reader.strain.end_time \
- data_reader.reduced_pad * data_reader.strain.delta_t - bdur
if not data_reader.state.is_extent_valid(state_start_time, bdur):
logging.info(
'%s strain buffer contains invalid data during lookback, '
'will not use for followup',
ifo
)
return None

# We won't require that all DQ checks be valid for now, except at
Expand All @@ -1849,10 +1875,15 @@ def followup_event_significance(ifo, data_reader, bank,
dq_start_time = onsource_start - duration / 2.0
dq_duration = onsource_end - onsource_start + duration
if not data_reader.dq.is_extent_valid(dq_start_time, dq_duration):
logging.info(
'%s DQ buffer indicates invalid data during onsource window, '
'will not use for followup',
ifo
)
return None

# Calculate SNR time series for this duration
htilde = bank.get_template(template_id, min_buffer=bdur)
htilde = bank.get_template(template_id, td_samples=bdur_samples)
stilde = data_reader.overwhitened_data(htilde.delta_f)

sigma2 = htilde.sigmasq(stilde.psd)
Expand Down
59 changes: 42 additions & 17 deletions pycbc/waveform/bank.py
Original file line number Diff line number Diff line change
Expand Up @@ -582,34 +582,59 @@ def __getitem__(self, index):

return self.get_template(index)

def get_template(self, index, min_buffer=None):
def get_template(self, index, td_samples=None):
"""Calculate and return the frequency-domain waveform for the template
with the given index. The frequency resolution can optionally be
controlled by giving the duration of the waveform in the time domain in
terms of number of samples.
Parameters
----------
index: int
Index of the template in the bank.
td_samples: int, optional
Duration of the waveform in the time domain, in units of samples,
used to define the length and resolution of the frequency series.
If not given, this is calculated from the template parameters.
Returns
-------
htilde: FrequencySeries
Template waveform in the frequency domain.
"""
approximant = self.approximant(index)
f_end = self.end_frequency(index)
flow = self.table[index].f_lower

# Determine the length of time of the filter, rounded up to
# nearest power of two
if min_buffer is None:
min_buffer = self.minimum_buffer
min_buffer += 0.5
if td_samples is None:
time_duration = self.minimum_buffer
time_duration += 0.5

from pycbc.waveform.waveform import props
p = props(self.table[index])
p.pop('approximant')
buff_size = pycbc.waveform.get_waveform_filter_length_in_time(approximant, **p)
if not buff_size:
raise RuntimeError('Template waveform %s not recognized!' % approximant)
from pycbc.waveform.waveform import props
p = props(self.table[index])
p.pop('approximant')
waveform_duration = pycbc.waveform.get_waveform_filter_length_in_time(
approximant, **p
)
if waveform_duration is None:
raise RuntimeError(f'Template waveform {approximant} not recognized!')
time_duration += waveform_duration
td_samples = self.round_up(time_duration * self.sample_rate)

tlen = self.round_up((buff_size + min_buffer) * self.sample_rate)
flen = int(tlen / 2 + 1)
flen = int(td_samples / 2 + 1)

delta_f = self.sample_rate / float(tlen)
delta_f = self.sample_rate / float(td_samples)

if f_end is None or f_end >= (flen * delta_f):
f_end = (flen - 1) * delta_f

logging.info("Generating %s, %ss, %i, starting from %s Hz",
approximant, 1.0 / delta_f, index, flow)
logging.info(
"Generating %s, duration %s s, index %i, starting from %s Hz",
approximant,
1.0 / delta_f,
index,
flow
)

# Get the waveform filter
distance = 1.0 / DYN_RANGE_FAC
Expand Down

0 comments on commit 037f5cd

Please sign in to comment.