From 30b86eea9401af6baf296faee3e25c7adb4f7ce2 Mon Sep 17 00:00:00 2001 From: Yifan Wang Date: Thu, 1 Jun 2023 15:09:03 +0200 Subject: [PATCH 01/21] apply time shift to account for the gating start time is a subsampling data point. fix a few typos --- pycbc/inference/models/gated_gaussian_noise.py | 10 +++++++++- pycbc/inference/models/marginalized_gaussian_noise.py | 4 ++-- 2 files changed, 11 insertions(+), 3 deletions(-) diff --git a/pycbc/inference/models/gated_gaussian_noise.py b/pycbc/inference/models/gated_gaussian_noise.py index 8852ecc8583..bd55e6f1d39 100644 --- a/pycbc/inference/models/gated_gaussian_noise.py +++ b/pycbc/inference/models/gated_gaussian_noise.py @@ -27,7 +27,7 @@ from pycbc.types import FrequencySeries from pycbc.detector import Detector from pycbc.pnutils import hybrid_meco_frequency -from pycbc.waveform.utils import time_from_frequencyseries +from pycbc.waveform.utils import (time_from_frequencyseries, apply_fd_time_shift) from pycbc.waveform import generator from pycbc.filter import highpass from .gaussian_noise import (BaseGaussianNoise, create_waveform_generator) @@ -559,6 +559,14 @@ def _loglikelihood(self): ht = h.to_timeseries() res = data - ht rtilde = res.to_frequencyseries() + #get the subsampling offset time + gateenddelay = gatestartdelay + dgatedelay + gateidx = numpy.floor(float(gateenddelay - res.start_time) * res.sample_rate) + shifttime = gateenddelay - res.get_sample_times[gateidx] + #apply time shift such that the gating start time lands on a data sample + rtilde = apply_fd_time_shift(rtilde, -shifttime, copy=False) + res = rtilde.to_timeseries() + #gating gated_res = res.gate(gatestartdelay + dgatedelay/2, window=dgatedelay/2, copy=True, invpsd=invpsd, method='paint') diff --git a/pycbc/inference/models/marginalized_gaussian_noise.py b/pycbc/inference/models/marginalized_gaussian_noise.py index 6db6ade80e4..0a1342d0432 100644 --- a/pycbc/inference/models/marginalized_gaussian_noise.py +++ b/pycbc/inference/models/marginalized_gaussian_noise.py @@ -229,7 +229,7 @@ def __init__(self, variable_params, generator_class=generator.FDomainDetFrameTwoPolNoRespGenerator, gates=self.gates, **kwargs['static_params']) else: - # create a waveform generator for each ifo respestively + # create a waveform generator for each ifo respectively self.waveform_generator = {} for det in self.data: self.waveform_generator[det] = create_waveform_generator( @@ -386,7 +386,7 @@ def __init__(self, variable_params, data, low_frequency_cutoff, psds=None, generator_class=generator.FDomainDetFrameTwoPolGenerator, gates=self.gates, **kwargs['static_params']) else: - # create a waveform generator for each ifo respestively + # create a waveform generator for each ifo respectively self.waveform_generator = {} for det in self.data: self.waveform_generator[det] = create_waveform_generator( From eb40555183cf4952ff215430ab3d1b626b5a1655 Mon Sep 17 00:00:00 2001 From: Yifan Wang Date: Fri, 2 Jun 2023 11:57:02 +0000 Subject: [PATCH 02/21] fix int type error --- pycbc/inference/models/gated_gaussian_noise.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pycbc/inference/models/gated_gaussian_noise.py b/pycbc/inference/models/gated_gaussian_noise.py index bd55e6f1d39..b9ddf1a5772 100644 --- a/pycbc/inference/models/gated_gaussian_noise.py +++ b/pycbc/inference/models/gated_gaussian_noise.py @@ -561,8 +561,8 @@ def _loglikelihood(self): rtilde = res.to_frequencyseries() #get the subsampling offset time gateenddelay = gatestartdelay + dgatedelay - gateidx = numpy.floor(float(gateenddelay - res.start_time) * res.sample_rate) - shifttime = gateenddelay - res.get_sample_times[gateidx] + gateidx = numpy.floor(float(gateenddelay - res.start_time) * res.sample_rate).astype(int) + shifttime = gateenddelay - res.sample_times[gateidx] #apply time shift such that the gating start time lands on a data sample rtilde = apply_fd_time_shift(rtilde, -shifttime, copy=False) res = rtilde.to_timeseries() From d12d8b2a3c910e8b073e1754ec727f52aa2e4bc7 Mon Sep 17 00:00:00 2001 From: Yifan Wang Date: Thu, 15 Jun 2023 15:31:05 +0200 Subject: [PATCH 03/21] correct a misunderstanding for apply_fd_time_shift --- pycbc/inference/models/gated_gaussian_noise.py | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/pycbc/inference/models/gated_gaussian_noise.py b/pycbc/inference/models/gated_gaussian_noise.py index b9ddf1a5772..a9667475a4f 100644 --- a/pycbc/inference/models/gated_gaussian_noise.py +++ b/pycbc/inference/models/gated_gaussian_noise.py @@ -561,17 +561,19 @@ def _loglikelihood(self): rtilde = res.to_frequencyseries() #get the subsampling offset time gateenddelay = gatestartdelay + dgatedelay - gateidx = numpy.floor(float(gateenddelay - res.start_time) * res.sample_rate).astype(int) - shifttime = gateenddelay - res.sample_times[gateidx] + gateidx = numpy.floor( + float(gateenddelay - res.start_time) * res.sample_rate).astype(int) #apply time shift such that the gating start time lands on a data sample - rtilde = apply_fd_time_shift(rtilde, -shifttime, copy=False) - res = rtilde.to_timeseries() + rtilde_shift = apply_fd_time_shift(rtilde, res.sample_times[gateidx], copy=True) + res_shift = rtilde_shift.to_timeseries() #gating - gated_res = res.gate(gatestartdelay + dgatedelay/2, + gated_res_shift = res_shift.gate(gatestartdelay + dgatedelay/2, window=dgatedelay/2, copy=True, invpsd=invpsd, method='paint') self.current_proj[det] = (gated_res.proj, gated_res.projslc) - gated_rtilde = gated_res.to_frequencyseries() + gated_rtilde_shift = gated_res_shift.to_frequencyseries() + # shift back to the actual starting time + gated_rtilde = apply_fd_time_shift(gated_rtilde_shift, gateenddelay, copy=True) # overwhiten gated_rtilde *= invpsd rr = 4 * invpsd.delta_f * rtilde[slc].inner(gated_rtilde[slc]).real From f3aa2577a6e25144bb600eaab493cbebfd16691b Mon Sep 17 00:00:00 2001 From: Yifan Wang Date: Thu, 15 Jun 2023 15:52:23 +0200 Subject: [PATCH 04/21] correct a unused variable --- pycbc/inference/models/gated_gaussian_noise.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pycbc/inference/models/gated_gaussian_noise.py b/pycbc/inference/models/gated_gaussian_noise.py index a9667475a4f..0e7733f4560 100644 --- a/pycbc/inference/models/gated_gaussian_noise.py +++ b/pycbc/inference/models/gated_gaussian_noise.py @@ -570,7 +570,7 @@ def _loglikelihood(self): gated_res_shift = res_shift.gate(gatestartdelay + dgatedelay/2, window=dgatedelay/2, copy=True, invpsd=invpsd, method='paint') - self.current_proj[det] = (gated_res.proj, gated_res.projslc) + self.current_proj[det] = (gated_res_shift.proj, gated_res_shift.projslc) gated_rtilde_shift = gated_res_shift.to_frequencyseries() # shift back to the actual starting time gated_rtilde = apply_fd_time_shift(gated_rtilde_shift, gateenddelay, copy=True) From 6c45dca1c0b5b73e530a95df550059b3b0c631e0 Mon Sep 17 00:00:00 2001 From: Yifan Wang Date: Sat, 17 Jun 2023 19:11:51 +0200 Subject: [PATCH 05/21] correct the shift time --- pycbc/inference/models/gated_gaussian_noise.py | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/pycbc/inference/models/gated_gaussian_noise.py b/pycbc/inference/models/gated_gaussian_noise.py index 0e7733f4560..6146c57da12 100644 --- a/pycbc/inference/models/gated_gaussian_noise.py +++ b/pycbc/inference/models/gated_gaussian_noise.py @@ -563,8 +563,9 @@ def _loglikelihood(self): gateenddelay = gatestartdelay + dgatedelay gateidx = numpy.floor( float(gateenddelay - res.start_time) * res.sample_rate).astype(int) - #apply time shift such that the gating start time lands on a data sample - rtilde_shift = apply_fd_time_shift(rtilde, res.sample_times[gateidx], copy=True) + offset = res.sample_times[gateidx] - gateenddelay + #apply time shift such that the gating time range lands on a data sample + rtilde_shift = apply_fd_time_shift(rtilde, offset + rtilde.epoch, copy=True) res_shift = rtilde_shift.to_timeseries() #gating gated_res_shift = res_shift.gate(gatestartdelay + dgatedelay/2, @@ -572,11 +573,9 @@ def _loglikelihood(self): invpsd=invpsd, method='paint') self.current_proj[det] = (gated_res_shift.proj, gated_res_shift.projslc) gated_rtilde_shift = gated_res_shift.to_frequencyseries() - # shift back to the actual starting time - gated_rtilde = apply_fd_time_shift(gated_rtilde_shift, gateenddelay, copy=True) # overwhiten - gated_rtilde *= invpsd - rr = 4 * invpsd.delta_f * rtilde[slc].inner(gated_rtilde[slc]).real + gated_rtilde_shift *= invpsd + rr = 4 * invpsd.delta_f * rtilde_shift[slc].inner(gated_rtilde_shift[slc]).real logl += norm - 0.5*rr return float(logl) From 46a48eb2a3a75369c145a94941cc679c8692550e Mon Sep 17 00:00:00 2001 From: Yifan Wang Date: Thu, 22 Jun 2023 16:08:27 +0200 Subject: [PATCH 06/21] address cc issues --- pycbc/inference/models/gated_gaussian_noise.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/pycbc/inference/models/gated_gaussian_noise.py b/pycbc/inference/models/gated_gaussian_noise.py index 6146c57da12..4c89bf7cfda 100644 --- a/pycbc/inference/models/gated_gaussian_noise.py +++ b/pycbc/inference/models/gated_gaussian_noise.py @@ -559,18 +559,18 @@ def _loglikelihood(self): ht = h.to_timeseries() res = data - ht rtilde = res.to_frequencyseries() - #get the subsampling offset time + # get the subsampling offset time gateenddelay = gatestartdelay + dgatedelay - gateidx = numpy.floor( - float(gateenddelay - res.start_time) * res.sample_rate).astype(int) + gateidx = numpy.floor(float( + gateenddelay - res.start_time) * res.sample_rate).astype(int) offset = res.sample_times[gateidx] - gateenddelay - #apply time shift such that the gating time range lands on a data sample + # shift such that the end gating time lands on a data sample rtilde_shift = apply_fd_time_shift(rtilde, offset + rtilde.epoch, copy=True) res_shift = rtilde_shift.to_timeseries() - #gating + # gating gated_res_shift = res_shift.gate(gatestartdelay + dgatedelay/2, - window=dgatedelay/2, copy=True, - invpsd=invpsd, method='paint') + window=dgatedelay/2, copy=True, + invpsd=invpsd, method='paint') self.current_proj[det] = (gated_res_shift.proj, gated_res_shift.projslc) gated_rtilde_shift = gated_res_shift.to_frequencyseries() # overwhiten From 2cc658c39a2d7ffe34fa907d5401bbb72b0182b4 Mon Sep 17 00:00:00 2001 From: Yifan Wang Date: Sun, 2 Jul 2023 00:13:29 +0200 Subject: [PATCH 07/21] apply static method --- .../inference/models/gated_gaussian_noise.py | 28 ++++++++++--------- 1 file changed, 15 insertions(+), 13 deletions(-) diff --git a/pycbc/inference/models/gated_gaussian_noise.py b/pycbc/inference/models/gated_gaussian_noise.py index 4c89bf7cfda..1f44854441b 100644 --- a/pycbc/inference/models/gated_gaussian_noise.py +++ b/pycbc/inference/models/gated_gaussian_noise.py @@ -438,6 +438,17 @@ def _fd_data_from_strain_dict(opts, strain_dict, psd_strain_dict): if lfs is not None: opts.low_frequency_cutoff = lfs return out + + @staticmethod + def shift_to_integer_sample(fs, gate_end_time): + """Shift the gating end time to the cloest (floor) data sample + """ + gateidx = numpy.floor(float( + gateenddelay - res.start_time) * res.sample_rate).astype(int) + offset = res.sample_times[gateidx] - gateenddelay + # shift such that the end gating time lands on a data sample + fs_shift = apply_fd_time_shift(fs, offset + rtilde.epoch, copy=True) + return fs_shift def write_metadata(self, fp, group=None): """Adds writing the psds, and analyzed detectors. @@ -534,9 +545,9 @@ def _loglikelihood(self): float The value of the log likelihood. """ - # generate the template waveform + # generate the residuals try: - wfs = self.get_waveforms() + rtildes = self.get_residuals() except NoWaveformError: return self._nowaveform_logl() except FailedWaveformError as e: @@ -547,25 +558,16 @@ def _loglikelihood(self): gate_times = self.get_gate_times() logl = 0. self.current_proj.clear() - for det, h in wfs.items(): + for det, rtilde in rtiles.items(): invpsd = self._invpsds[det] norm = self.det_lognorm(det) gatestartdelay, dgatedelay = gate_times[det] # we always filter the entire segment starting from kmin, since the # gated series may have high frequency components slc = slice(self._kmin[det], self._kmax[det]) - # calculate the residual - data = self.td_data[det] - ht = h.to_timeseries() - res = data - ht - rtilde = res.to_frequencyseries() # get the subsampling offset time gateenddelay = gatestartdelay + dgatedelay - gateidx = numpy.floor(float( - gateenddelay - res.start_time) * res.sample_rate).astype(int) - offset = res.sample_times[gateidx] - gateenddelay - # shift such that the end gating time lands on a data sample - rtilde_shift = apply_fd_time_shift(rtilde, offset + rtilde.epoch, copy=True) + rtilde_shift = self.shift_to_integer_sample(rtilde,gateenddelay) res_shift = rtilde_shift.to_timeseries() # gating gated_res_shift = res_shift.gate(gatestartdelay + dgatedelay/2, From 8dae8bdf6070bbd45ea686113625c97d691b5288 Mon Sep 17 00:00:00 2001 From: Yifan Wang Date: Sun, 2 Jul 2023 21:21:13 +0200 Subject: [PATCH 08/21] complete the gated_gaussian_noise --- .../inference/models/gated_gaussian_noise.py | 28 ++++++++++++++----- 1 file changed, 21 insertions(+), 7 deletions(-) diff --git a/pycbc/inference/models/gated_gaussian_noise.py b/pycbc/inference/models/gated_gaussian_noise.py index 1f44854441b..e044e1e4a5d 100644 --- a/pycbc/inference/models/gated_gaussian_noise.py +++ b/pycbc/inference/models/gated_gaussian_noise.py @@ -441,13 +441,27 @@ def _fd_data_from_strain_dict(opts, strain_dict, psd_strain_dict): @staticmethod def shift_to_integer_sample(fs, gate_end_time): - """Shift the gating end time to the cloest (floor) data sample + """Shift the frequency series by an amount of time such that the + gate_end_time lands on a specific data sample (chosen to be the + nearest floor data sample). + + Parameters + ---------- + fs : FrequencySeries + Frequency series to be shifted. + gate_end_time : float + Gating end time. + + Returns + ------- + fs_shift : FrequencySeries + Frequency series after being shifted """ - gateidx = numpy.floor(float( - gateenddelay - res.start_time) * res.sample_rate).astype(int) - offset = res.sample_times[gateidx] - gateenddelay - # shift such that the end gating time lands on a data sample - fs_shift = apply_fd_time_shift(fs, offset + rtilde.epoch, copy=True) + floor_idx = numpy.floor( + (gate_end_delay - float(fs.start_time)) * fs.sample_rate).astype(int) + # offset is less than 0 to shift to an earlier time + offset = fs.sample_times[floor_idx] - gate_end_delay + fs_shift = apply_fd_time_shift(fs, offset + fs.epoch, copy=True) return fs_shift def write_metadata(self, fp, group=None): @@ -565,7 +579,7 @@ def _loglikelihood(self): # we always filter the entire segment starting from kmin, since the # gated series may have high frequency components slc = slice(self._kmin[det], self._kmax[det]) - # get the subsampling offset time + # get the end time of gating mask gateenddelay = gatestartdelay + dgatedelay rtilde_shift = self.shift_to_integer_sample(rtilde,gateenddelay) res_shift = rtilde_shift.to_timeseries() From 603462bdd543dd477af09f445ec939e2a6c1dc44 Mon Sep 17 00:00:00 2001 From: Yifan Wang Date: Sun, 2 Jul 2023 22:46:44 +0200 Subject: [PATCH 09/21] typo --- pycbc/inference/models/gated_gaussian_noise.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pycbc/inference/models/gated_gaussian_noise.py b/pycbc/inference/models/gated_gaussian_noise.py index e044e1e4a5d..b9ef61ed4e8 100644 --- a/pycbc/inference/models/gated_gaussian_noise.py +++ b/pycbc/inference/models/gated_gaussian_noise.py @@ -572,7 +572,7 @@ def _loglikelihood(self): gate_times = self.get_gate_times() logl = 0. self.current_proj.clear() - for det, rtilde in rtiles.items(): + for det, rtilde in rtildes.items(): invpsd = self._invpsds[det] norm = self.det_lognorm(det) gatestartdelay, dgatedelay = gate_times[det] From 793f1a6b4db6eb370936cc3e3138ba0c16c9793c Mon Sep 17 00:00:00 2001 From: Yifan Wang Date: Sun, 2 Jul 2023 22:49:13 +0200 Subject: [PATCH 10/21] typo --- pycbc/inference/models/gated_gaussian_noise.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pycbc/inference/models/gated_gaussian_noise.py b/pycbc/inference/models/gated_gaussian_noise.py index b9ef61ed4e8..3aa077510e5 100644 --- a/pycbc/inference/models/gated_gaussian_noise.py +++ b/pycbc/inference/models/gated_gaussian_noise.py @@ -458,9 +458,9 @@ def shift_to_integer_sample(fs, gate_end_time): Frequency series after being shifted """ floor_idx = numpy.floor( - (gate_end_delay - float(fs.start_time)) * fs.sample_rate).astype(int) + (gate_end_time - float(fs.start_time)) * fs.sample_rate).astype(int) # offset is less than 0 to shift to an earlier time - offset = fs.sample_times[floor_idx] - gate_end_delay + offset = fs.sample_times[floor_idx] - gate_end_time fs_shift = apply_fd_time_shift(fs, offset + fs.epoch, copy=True) return fs_shift From 0eeff1ab2ee244f89ba9c53f46601f5e69718aed Mon Sep 17 00:00:00 2001 From: Yifan Wang Date: Sun, 2 Jul 2023 23:04:58 +0200 Subject: [PATCH 11/21] re express time shift --- pycbc/inference/models/gated_gaussian_noise.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/pycbc/inference/models/gated_gaussian_noise.py b/pycbc/inference/models/gated_gaussian_noise.py index 3aa077510e5..aa478ff76f8 100644 --- a/pycbc/inference/models/gated_gaussian_noise.py +++ b/pycbc/inference/models/gated_gaussian_noise.py @@ -457,10 +457,12 @@ def shift_to_integer_sample(fs, gate_end_time): fs_shift : FrequencySeries Frequency series after being shifted """ + dt = 1 / fs.sample_rate floor_idx = numpy.floor( - (gate_end_time - float(fs.start_time)) * fs.sample_rate).astype(int) - # offset is less than 0 to shift to an earlier time - offset = fs.sample_times[floor_idx] - gate_end_time + (gate_end_time - float(fs.start_time)) / dt).astype(int) + floor_time = float(fs.start_time) + floor_idx * dt + # offset is less than 0, so shifting to an earlier time + offset = floor_time - gate_end_time fs_shift = apply_fd_time_shift(fs, offset + fs.epoch, copy=True) return fs_shift From 99520d87580a914ce1ab01545f1bd14f5550b944 Mon Sep 17 00:00:00 2001 From: Yifan Wang Date: Mon, 3 Jul 2023 00:18:52 +0200 Subject: [PATCH 12/21] space --- pycbc/inference/models/gated_gaussian_noise.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pycbc/inference/models/gated_gaussian_noise.py b/pycbc/inference/models/gated_gaussian_noise.py index aa478ff76f8..7955c106d2a 100644 --- a/pycbc/inference/models/gated_gaussian_noise.py +++ b/pycbc/inference/models/gated_gaussian_noise.py @@ -583,7 +583,7 @@ def _loglikelihood(self): slc = slice(self._kmin[det], self._kmax[det]) # get the end time of gating mask gateenddelay = gatestartdelay + dgatedelay - rtilde_shift = self.shift_to_integer_sample(rtilde,gateenddelay) + rtilde_shift = self.shift_to_integer_sample(rtilde, gateenddelay) res_shift = rtilde_shift.to_timeseries() # gating gated_res_shift = res_shift.gate(gatestartdelay + dgatedelay/2, From fc589268a40689d63b41d44743aea70872f285a0 Mon Sep 17 00:00:00 2001 From: Yifan Wang Date: Tue, 4 Jul 2023 15:37:44 +0200 Subject: [PATCH 13/21] start to work on pol marg model --- pycbc/inference/models/gated_gaussian_noise.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/pycbc/inference/models/gated_gaussian_noise.py b/pycbc/inference/models/gated_gaussian_noise.py index 7955c106d2a..5ca499690ae 100644 --- a/pycbc/inference/models/gated_gaussian_noise.py +++ b/pycbc/inference/models/gated_gaussian_noise.py @@ -270,11 +270,16 @@ def get_gated_data(self): cache = self._gated_data[det] = {} invpsd = self._invpsds[det] gatestartdelay, dgatedelay = gate_times[det] + gateenddelay = gatestartdelay + dgatedel try: dtilde = cache[gatestartdelay, dgatedelay] except KeyError: # doesn't exist yet, or the gate times changed cache.clear() + dtilde = d.to_frequencyseries() + dtilde = self.shift_to_integer_sample( + dtilde, gateenddelay) + d = dtilde.to_timeseries() d = d.gate(gatestartdelay + dgatedelay/2, window=dgatedelay/2, copy=True, invpsd=invpsd, method='paint') @@ -698,9 +703,11 @@ def get_gated_waveforms(self): for det in wfs: invpsd = self._invpsds[det] gatestartdelay, dgatedelay = gate_times[det] + gateenddelay = gatedstartdelay + dgatedelay # the waveforms are a dictionary of (hp, hc) pols = [] for h in wfs[det]: + h = self.shift_to_integer_sample(h, gateenddelay) ht = h.to_timeseries() ht = ht.gate(gatestartdelay + dgatedelay/2, window=dgatedelay/2, copy=False, From d29e92c5ff728619f23e3cc7536a19343f0da414 Mon Sep 17 00:00:00 2001 From: Yifan Wang Date: Thu, 6 Jul 2023 16:02:22 +0200 Subject: [PATCH 14/21] complete shifting for marg pol model --- pycbc/inference/models/gated_gaussian_noise.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/pycbc/inference/models/gated_gaussian_noise.py b/pycbc/inference/models/gated_gaussian_noise.py index 5ca499690ae..7432570279b 100644 --- a/pycbc/inference/models/gated_gaussian_noise.py +++ b/pycbc/inference/models/gated_gaussian_noise.py @@ -749,6 +749,7 @@ def _loglikelihood(self): # get the gated waveforms and data gated_wfs = self.get_gated_waveforms() gated_data = self.get_gated_data() + gate_times = self.get_gate_times() # cycle over loglr = 0. lognl = 0. @@ -764,12 +765,16 @@ def _loglikelihood(self): # we always filter the entire segment starting from kmin, since the # gated series may have high frequency components slc = slice(self._kmin[det], self._kmax[det]) + # get gated end time + gatestartdelay, dgatedelay = gate_times[det] + gateenddelay = gatedstartdelay + dgatedelay # get the gated values gated_hp, gated_hc = gated_wfs[det] gated_d = gated_data[det] # we'll overwhiten the ungated data and waveforms for computing # inner products d = self._overwhitened_data[det] + d = self.shift_to_integer_sample(d, gateenddelay) # overwhiten the hp and hc # we'll do this in place for computational efficiency, but as a # result we'll clear the current waveforms cache so a repeated call @@ -778,6 +783,8 @@ def _loglikelihood(self): invpsd = self._invpsds[det] hp *= invpsd hc *= invpsd + hp = self.shift_to_integer_sample(hp, gateenddelay) + hc = self.shift_to_integer_sample(hc, gateenddelay) # get the various gated inner products hpd = hp[slc].inner(gated_d[slc]).real # hcd = hc[slc].inner(gated_d[slc]).real # From 10f28c8e61b2f210b5f3fa11fc6b0ccfa8b94da9 Mon Sep 17 00:00:00 2001 From: Yifan Wang Date: Tue, 11 Jul 2023 22:40:32 +0200 Subject: [PATCH 15/21] instead of modifying models, now implementing time shifting only in gate-and-inpaint session --- .../inference/models/gated_gaussian_noise.py | 61 ++++++------------- pycbc/strain/gate.py | 6 +- pycbc/types/timeseries.py | 15 ++++- 3 files changed, 33 insertions(+), 49 deletions(-) diff --git a/pycbc/inference/models/gated_gaussian_noise.py b/pycbc/inference/models/gated_gaussian_noise.py index 7955c106d2a..4c284fdfbd2 100644 --- a/pycbc/inference/models/gated_gaussian_noise.py +++ b/pycbc/inference/models/gated_gaussian_noise.py @@ -27,7 +27,7 @@ from pycbc.types import FrequencySeries from pycbc.detector import Detector from pycbc.pnutils import hybrid_meco_frequency -from pycbc.waveform.utils import (time_from_frequencyseries, apply_fd_time_shift) +from pycbc.waveform.utils import time_from_frequencyseries from pycbc.waveform import generator from pycbc.filter import highpass from .gaussian_noise import (BaseGaussianNoise, create_waveform_generator) @@ -438,33 +438,6 @@ def _fd_data_from_strain_dict(opts, strain_dict, psd_strain_dict): if lfs is not None: opts.low_frequency_cutoff = lfs return out - - @staticmethod - def shift_to_integer_sample(fs, gate_end_time): - """Shift the frequency series by an amount of time such that the - gate_end_time lands on a specific data sample (chosen to be the - nearest floor data sample). - - Parameters - ---------- - fs : FrequencySeries - Frequency series to be shifted. - gate_end_time : float - Gating end time. - - Returns - ------- - fs_shift : FrequencySeries - Frequency series after being shifted - """ - dt = 1 / fs.sample_rate - floor_idx = numpy.floor( - (gate_end_time - float(fs.start_time)) / dt).astype(int) - floor_time = float(fs.start_time) + floor_idx * dt - # offset is less than 0, so shifting to an earlier time - offset = floor_time - gate_end_time - fs_shift = apply_fd_time_shift(fs, offset + fs.epoch, copy=True) - return fs_shift def write_metadata(self, fp, group=None): """Adds writing the psds, and analyzed detectors. @@ -561,9 +534,9 @@ def _loglikelihood(self): float The value of the log likelihood. """ - # generate the residuals + # generate the template waveform try: - rtildes = self.get_residuals() + wfs = self.get_waveforms() except NoWaveformError: return self._nowaveform_logl() except FailedWaveformError as e: @@ -574,26 +547,26 @@ def _loglikelihood(self): gate_times = self.get_gate_times() logl = 0. self.current_proj.clear() - for det, rtilde in rtildes.items(): + for det, h in wfs.items(): invpsd = self._invpsds[det] norm = self.det_lognorm(det) gatestartdelay, dgatedelay = gate_times[det] # we always filter the entire segment starting from kmin, since the # gated series may have high frequency components slc = slice(self._kmin[det], self._kmax[det]) - # get the end time of gating mask - gateenddelay = gatestartdelay + dgatedelay - rtilde_shift = self.shift_to_integer_sample(rtilde, gateenddelay) - res_shift = rtilde_shift.to_timeseries() - # gating - gated_res_shift = res_shift.gate(gatestartdelay + dgatedelay/2, - window=dgatedelay/2, copy=True, - invpsd=invpsd, method='paint') - self.current_proj[det] = (gated_res_shift.proj, gated_res_shift.projslc) - gated_rtilde_shift = gated_res_shift.to_frequencyseries() + # calculate the residual + data = self.td_data[det] + ht = h.to_timeseries() + res = data - ht + rtilde = res.to_frequencyseries() + gated_res = res.gate(gatestartdelay + dgatedelay/2, + window=dgatedelay/2, copy=True, + invpsd=invpsd, method='paint') + self.current_proj[det] = (gated_res.proj, gated_res.projslc) + gated_rtilde = gated_res.to_frequencyseries() # overwhiten - gated_rtilde_shift *= invpsd - rr = 4 * invpsd.delta_f * rtilde_shift[slc].inner(gated_rtilde_shift[slc]).real + gated_rtilde *= invpsd + rr = 4 * invpsd.delta_f * rtilde[slc].inner(gated_rtilde[slc]).real logl += norm - 0.5*rr return float(logl) @@ -799,4 +772,4 @@ def _loglikelihood(self): setattr(self._current_stats, 'maxl_logl', loglr[idx] + lognl) # compute the marginalized log likelihood marglogl = special.logsumexp(loglr) + lognl - numpy.log(len(self.pol)) - return float(marglogl) + return float(marglogl) \ No newline at end of file diff --git a/pycbc/strain/gate.py b/pycbc/strain/gate.py index e294beeae87..650c93409dd 100644 --- a/pycbc/strain/gate.py +++ b/pycbc/strain/gate.py @@ -168,15 +168,15 @@ def gate_and_paint(data, lindex, rindex, invpsd, copy=True): # Copy the data and zero inside the hole if copy: data = data.copy() - data[lindex:rindex] = 0 + data[lindex:rindex+1] = 0 # get the over-whitened gated data tdfilter = invpsd.astype('complex').to_timeseries() * invpsd.delta_t owhgated_data = (data.to_frequencyseries() * invpsd).to_timeseries() # remove the projection into the null space - proj = linalg.solve_toeplitz(tdfilter[:(rindex - lindex)], - owhgated_data[lindex:rindex]) + proj = linalg.solve_toeplitz(tdfilter[:(rindex + 1 - lindex)], + owhgated_data[lindex:rindex+1]) data[lindex:rindex] -= proj data.projslc = (lindex, rindex) data.proj = proj diff --git a/pycbc/types/timeseries.py b/pycbc/types/timeseries.py index 910922065ad..ad854395034 100644 --- a/pycbc/types/timeseries.py +++ b/pycbc/types/timeseries.py @@ -595,15 +595,26 @@ def gate(self, time, window=0.25, method='taper', copy=True, # Uses the hole-filling method of # https://arxiv.org/pdf/1908.05644.pdf from pycbc.strain.gate import gate_and_paint + from pycbc.waveform.utils import apply_fd_time_shift if invpsd is None: # These are some bare minimum settings, normally you # should probably provide a psd invpsd = 1. / self.filter_psd(self.duration/32, self.delta_f, 0) lindex = int((time - window - self.start_time) / self.delta_t) - rindex = lindex + int(2 * window / self.delta_t) + rindex = int((time + window - self.start_time) / self.delta_t) lindex = lindex if lindex >= 0 else 0 rindex = rindex if rindex <= len(self) else len(self) - return gate_and_paint(data, lindex, rindex, invpsd, copy=False) + # shift data such that gate_end_time lands on a specific data sample + fdata = data.to_frequencyseries() + rindex_time = data.sample_times[rindex] + # offset is less than 0, so shifting to an earlier time + offset = rindex_time - (time + window) + fdata = apply_fd_time_shift(fdata, offset + fdata.epoch, copy=False) + # gate and paint in time domain + data = gate_and_paint(fdata.to_timeseries(), lindex, rindex, invpsd, copy=False) + # shift back to the original time + return apply_fd_time_shift(data.to_frequencyseries(), -offset + fdata.epoch, + copy=False).to_timeseries() elif method == 'hard': tslice = data.time_slice(time - window, time + window) tslice[:] = 0 From bd96f2b45651a903771ddeca3320c54c05c672a7 Mon Sep 17 00:00:00 2001 From: Yifan Wang Date: Tue, 11 Jul 2023 22:47:33 +0200 Subject: [PATCH 16/21] remove all changes in the model --- pycbc/inference/models/gated_gaussian_noise.py | 14 -------------- 1 file changed, 14 deletions(-) diff --git a/pycbc/inference/models/gated_gaussian_noise.py b/pycbc/inference/models/gated_gaussian_noise.py index 815399df922..4c284fdfbd2 100644 --- a/pycbc/inference/models/gated_gaussian_noise.py +++ b/pycbc/inference/models/gated_gaussian_noise.py @@ -270,16 +270,11 @@ def get_gated_data(self): cache = self._gated_data[det] = {} invpsd = self._invpsds[det] gatestartdelay, dgatedelay = gate_times[det] - gateenddelay = gatestartdelay + dgatedel try: dtilde = cache[gatestartdelay, dgatedelay] except KeyError: # doesn't exist yet, or the gate times changed cache.clear() - dtilde = d.to_frequencyseries() - dtilde = self.shift_to_integer_sample( - dtilde, gateenddelay) - d = dtilde.to_timeseries() d = d.gate(gatestartdelay + dgatedelay/2, window=dgatedelay/2, copy=True, invpsd=invpsd, method='paint') @@ -676,11 +671,9 @@ def get_gated_waveforms(self): for det in wfs: invpsd = self._invpsds[det] gatestartdelay, dgatedelay = gate_times[det] - gateenddelay = gatedstartdelay + dgatedelay # the waveforms are a dictionary of (hp, hc) pols = [] for h in wfs[det]: - h = self.shift_to_integer_sample(h, gateenddelay) ht = h.to_timeseries() ht = ht.gate(gatestartdelay + dgatedelay/2, window=dgatedelay/2, copy=False, @@ -722,7 +715,6 @@ def _loglikelihood(self): # get the gated waveforms and data gated_wfs = self.get_gated_waveforms() gated_data = self.get_gated_data() - gate_times = self.get_gate_times() # cycle over loglr = 0. lognl = 0. @@ -738,16 +730,12 @@ def _loglikelihood(self): # we always filter the entire segment starting from kmin, since the # gated series may have high frequency components slc = slice(self._kmin[det], self._kmax[det]) - # get gated end time - gatestartdelay, dgatedelay = gate_times[det] - gateenddelay = gatedstartdelay + dgatedelay # get the gated values gated_hp, gated_hc = gated_wfs[det] gated_d = gated_data[det] # we'll overwhiten the ungated data and waveforms for computing # inner products d = self._overwhitened_data[det] - d = self.shift_to_integer_sample(d, gateenddelay) # overwhiten the hp and hc # we'll do this in place for computational efficiency, but as a # result we'll clear the current waveforms cache so a repeated call @@ -756,8 +744,6 @@ def _loglikelihood(self): invpsd = self._invpsds[det] hp *= invpsd hc *= invpsd - hp = self.shift_to_integer_sample(hp, gateenddelay) - hc = self.shift_to_integer_sample(hc, gateenddelay) # get the various gated inner products hpd = hp[slc].inner(gated_d[slc]).real # hcd = hc[slc].inner(gated_d[slc]).real # From 76d5462d14cb0658bf96f6de61c453a2c1a0d3c8 Mon Sep 17 00:00:00 2001 From: Yifan Wang Date: Tue, 11 Jul 2023 22:52:38 +0200 Subject: [PATCH 17/21] remove changes in model --- pycbc/inference/models/gated_gaussian_noise.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pycbc/inference/models/gated_gaussian_noise.py b/pycbc/inference/models/gated_gaussian_noise.py index 4c284fdfbd2..8852ecc8583 100644 --- a/pycbc/inference/models/gated_gaussian_noise.py +++ b/pycbc/inference/models/gated_gaussian_noise.py @@ -772,4 +772,4 @@ def _loglikelihood(self): setattr(self._current_stats, 'maxl_logl', loglr[idx] + lognl) # compute the marginalized log likelihood marglogl = special.logsumexp(loglr) + lognl - numpy.log(len(self.pol)) - return float(marglogl) \ No newline at end of file + return float(marglogl) From 0870e914c0be08ef5cde2a44960efdfb0285da2e Mon Sep 17 00:00:00 2001 From: Yifan Wang Date: Wed, 12 Jul 2023 14:01:13 +0000 Subject: [PATCH 18/21] fix some bugs when fft/ifft converting, now productive --- pycbc/strain/gate.py | 2 +- pycbc/types/timeseries.py | 15 ++++++++++----- 2 files changed, 11 insertions(+), 6 deletions(-) diff --git a/pycbc/strain/gate.py b/pycbc/strain/gate.py index 650c93409dd..f060c55b835 100644 --- a/pycbc/strain/gate.py +++ b/pycbc/strain/gate.py @@ -177,7 +177,7 @@ def gate_and_paint(data, lindex, rindex, invpsd, copy=True): # remove the projection into the null space proj = linalg.solve_toeplitz(tdfilter[:(rindex + 1 - lindex)], owhgated_data[lindex:rindex+1]) - data[lindex:rindex] -= proj + data[lindex:rindex+1] -= proj data.projslc = (lindex, rindex) data.proj = proj return data diff --git a/pycbc/types/timeseries.py b/pycbc/types/timeseries.py index ad854395034..082fbccc390 100644 --- a/pycbc/types/timeseries.py +++ b/pycbc/types/timeseries.py @@ -607,14 +607,19 @@ def gate(self, time, window=0.25, method='taper', copy=True, # shift data such that gate_end_time lands on a specific data sample fdata = data.to_frequencyseries() rindex_time = data.sample_times[rindex] - # offset is less than 0, so shifting to an earlier time offset = rindex_time - (time + window) - fdata = apply_fd_time_shift(fdata, offset + fdata.epoch, copy=False) + fdata = apply_fd_time_shift(fdata, offset + fdata.epoch, copy=copy) # gate and paint in time domain - data = gate_and_paint(fdata.to_timeseries(), lindex, rindex, invpsd, copy=False) + data = fdata.to_timeseries() + data = gate_and_paint(data, lindex, rindex, invpsd, copy=copy) # shift back to the original time - return apply_fd_time_shift(data.to_frequencyseries(), -offset + fdata.epoch, - copy=False).to_timeseries() + fdata = data.to_frequencyseries() + fdata = apply_fd_time_shift(fdata, -offset + fdata.epoch, copy=copy) + tdata = fdata.to_timeseries() + # fill in projslc information + tdata.projslc = data.projslc + tdata.proj = data.proj + return tdata elif method == 'hard': tslice = data.time_slice(time - window, time + window) tslice[:] = 0 From 1186e5bac7e9c5f1c055e578ff77b9346983bc6e Mon Sep 17 00:00:00 2001 From: Yifan Wang Date: Thu, 13 Jul 2023 21:28:47 +0200 Subject: [PATCH 19/21] remove rindex+1 in gate and inpainting, add if condition for time shifting --- pycbc/strain/gate.py | 9 +++++---- pycbc/types/timeseries.py | 31 +++++++++++++++++-------------- 2 files changed, 22 insertions(+), 18 deletions(-) diff --git a/pycbc/strain/gate.py b/pycbc/strain/gate.py index f060c55b835..73f3fa112c9 100644 --- a/pycbc/strain/gate.py +++ b/pycbc/strain/gate.py @@ -168,16 +168,17 @@ def gate_and_paint(data, lindex, rindex, invpsd, copy=True): # Copy the data and zero inside the hole if copy: data = data.copy() - data[lindex:rindex+1] = 0 + # Here's ambiguity about when gate end time exactly is, rindex-1 or rindex? + data[lindex:rindex] = 0 # get the over-whitened gated data tdfilter = invpsd.astype('complex').to_timeseries() * invpsd.delta_t owhgated_data = (data.to_frequencyseries() * invpsd).to_timeseries() # remove the projection into the null space - proj = linalg.solve_toeplitz(tdfilter[:(rindex + 1 - lindex)], - owhgated_data[lindex:rindex+1]) - data[lindex:rindex+1] -= proj + proj = linalg.solve_toeplitz(tdfilter[:(rindex - lindex)], + owhgated_data[lindex:rindex]) + data[lindex:rindex] -= proj data.projslc = (lindex, rindex) data.proj = proj return data diff --git a/pycbc/types/timeseries.py b/pycbc/types/timeseries.py index 082fbccc390..4c5be4a75b1 100644 --- a/pycbc/types/timeseries.py +++ b/pycbc/types/timeseries.py @@ -604,22 +604,25 @@ def gate(self, time, window=0.25, method='taper', copy=True, rindex = int((time + window - self.start_time) / self.delta_t) lindex = lindex if lindex >= 0 else 0 rindex = rindex if rindex <= len(self) else len(self) - # shift data such that gate_end_time lands on a specific data sample - fdata = data.to_frequencyseries() rindex_time = data.sample_times[rindex] offset = rindex_time - (time + window) - fdata = apply_fd_time_shift(fdata, offset + fdata.epoch, copy=copy) - # gate and paint in time domain - data = fdata.to_timeseries() - data = gate_and_paint(data, lindex, rindex, invpsd, copy=copy) - # shift back to the original time - fdata = data.to_frequencyseries() - fdata = apply_fd_time_shift(fdata, -offset + fdata.epoch, copy=copy) - tdata = fdata.to_timeseries() - # fill in projslc information - tdata.projslc = data.projslc - tdata.proj = data.proj - return tdata + if offset == 0: + return gate_and_paint(data, lindex, rindex, invpsd, copy=False) + else: + # time shift such that gate end time lands on a specific data sample + fdata = data.to_frequencyseries() + fdata = apply_fd_time_shift(fdata, offset + fdata.epoch, copy=False) + # gate and paint in time domain + data = fdata.to_timeseries() + data = gate_and_paint(data, lindex, rindex, invpsd, copy=False) + # shift back to the original time + fdata = data.to_frequencyseries() + fdata = apply_fd_time_shift(fdata, -offset + fdata.epoch, copy=False) + tdata = fdata.to_timeseries() + # fill in projslc information + tdata.projslc = data.projslc + tdata.proj = data.proj + return tdata elif method == 'hard': tslice = data.time_slice(time - window, time + window) tslice[:] = 0 From b7534af35a044db4034ea14196cdfcb1155772f7 Mon Sep 17 00:00:00 2001 From: Yifan Wang Date: Fri, 15 Sep 2023 22:44:12 +0200 Subject: [PATCH 20/21] get rindex_time via actually doing the math, this is faster than retrieve the data sample time array --- pycbc/types/timeseries.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pycbc/types/timeseries.py b/pycbc/types/timeseries.py index 4c5be4a75b1..c05bdf1fb62 100644 --- a/pycbc/types/timeseries.py +++ b/pycbc/types/timeseries.py @@ -604,7 +604,7 @@ def gate(self, time, window=0.25, method='taper', copy=True, rindex = int((time + window - self.start_time) / self.delta_t) lindex = lindex if lindex >= 0 else 0 rindex = rindex if rindex <= len(self) else len(self) - rindex_time = data.sample_times[rindex] + rindex_time = float(self.start_time + rindex * self.delta_t) offset = rindex_time - (time + window) if offset == 0: return gate_and_paint(data, lindex, rindex, invpsd, copy=False) From 5250673801da446886578308079be3866dd9d704 Mon Sep 17 00:00:00 2001 From: Yifan Wang Date: Fri, 15 Sep 2023 23:07:19 +0200 Subject: [PATCH 21/21] remove the proj and projslc attributes in gated_gaussian_model, it's used for debugging before --- pycbc/inference/models/gated_gaussian_noise.py | 7 ------- pycbc/strain/gate.py | 2 -- pycbc/types/timeseries.py | 3 --- 3 files changed, 12 deletions(-) diff --git a/pycbc/inference/models/gated_gaussian_noise.py b/pycbc/inference/models/gated_gaussian_noise.py index 8852ecc8583..33e7ecc27e1 100644 --- a/pycbc/inference/models/gated_gaussian_noise.py +++ b/pycbc/inference/models/gated_gaussian_noise.py @@ -46,9 +46,6 @@ def __init__(self, variable_params, data, low_frequency_cutoff, psds=None, static_params=None, highpass_waveforms=False, **kwargs): # we'll want the time-domain data, so store that self._td_data = {} - # cache the current projection for debugging - self.current_proj = {} - self.current_nproj = {} # cache the overwhitened data self._overwhitened_data = {} # cache the current gated data @@ -383,7 +380,6 @@ def _lognl(self): self._det_lognls.clear() # get the times of the gates gate_times = self.get_gate_times() - self.current_nproj.clear() for det, invpsd in self._invpsds.items(): norm = self.det_lognorm(det) gatestartdelay, dgatedelay = gate_times[det] @@ -395,7 +391,6 @@ def _lognl(self): gated_dt = data.gate(gatestartdelay + dgatedelay/2, window=dgatedelay/2, copy=True, invpsd=invpsd, method='paint') - self.current_nproj[det] = (gated_dt.proj, gated_dt.projslc) # convert to the frequency series gated_d = gated_dt.to_frequencyseries() # overwhiten @@ -546,7 +541,6 @@ def _loglikelihood(self): # get the times of the gates gate_times = self.get_gate_times() logl = 0. - self.current_proj.clear() for det, h in wfs.items(): invpsd = self._invpsds[det] norm = self.det_lognorm(det) @@ -562,7 +556,6 @@ def _loglikelihood(self): gated_res = res.gate(gatestartdelay + dgatedelay/2, window=dgatedelay/2, copy=True, invpsd=invpsd, method='paint') - self.current_proj[det] = (gated_res.proj, gated_res.projslc) gated_rtilde = gated_res.to_frequencyseries() # overwhiten gated_rtilde *= invpsd diff --git a/pycbc/strain/gate.py b/pycbc/strain/gate.py index 73f3fa112c9..6c2a70f59c3 100644 --- a/pycbc/strain/gate.py +++ b/pycbc/strain/gate.py @@ -179,6 +179,4 @@ def gate_and_paint(data, lindex, rindex, invpsd, copy=True): proj = linalg.solve_toeplitz(tdfilter[:(rindex - lindex)], owhgated_data[lindex:rindex]) data[lindex:rindex] -= proj - data.projslc = (lindex, rindex) - data.proj = proj return data diff --git a/pycbc/types/timeseries.py b/pycbc/types/timeseries.py index c05bdf1fb62..3bf94065473 100644 --- a/pycbc/types/timeseries.py +++ b/pycbc/types/timeseries.py @@ -619,9 +619,6 @@ def gate(self, time, window=0.25, method='taper', copy=True, fdata = data.to_frequencyseries() fdata = apply_fd_time_shift(fdata, -offset + fdata.epoch, copy=False) tdata = fdata.to_timeseries() - # fill in projslc information - tdata.projslc = data.projslc - tdata.proj = data.proj return tdata elif method == 'hard': tslice = data.time_slice(time - window, time + window)