In [3]:
import logging
import pycbc
pycbc.init_logging(True)
%matplotlib notebook
from matplotlib import pyplot
import numpy
from pycbc.inference import io, models

2022-06-03 13:57:03,313 Setting verbosity to 0

Problematic 221 - without highpass¶

In [4]:
# load the model
fn = '/work/cdcapano/projects/gated_gaussian/gw190521/injections/imrdraws/margpol-inj100/inj03/samples-220_221-18.hdf'
fp = io.loadfile(fn, 'r')
cp = fp.read_config_file()
psds = fp.read_psds()
data = fp.read_data()
model = models.read_from_config(cp, data=data, psds=psds, highpass_waveforms=False)

2022-06-03 13:57:27,761 Setting up priors for each parameter
2022-06-03 13:57:27,764 No sampling_params section read from config file
2022-06-03 13:57:27,765 Loading waveform transforms
2022-06-03 13:57:27,770 Determining analysis times to use
2022-06-03 13:57:27,770 Padding H1 analysis start and end times by 4 (= psd-inverse-length/2) seconds to account for PSD wrap around effects.
2022-06-03 13:57:27,771 Padding L1 analysis start and end times by 4 (= psd-inverse-length/2) seconds to account for PSD wrap around effects.
2022-06-03 13:57:27,771 Padding V1 analysis start and end times by 4 (= psd-inverse-length/2) seconds to account for PSD wrap around effects.
In [5]:
# get the maxL point
samples = fp.read_samples(list(fp['samples'].keys()))
maxidx = samples.loglikelihood.argmax()
params = {p: samples[p][maxidx] for p in model.variable_params}
model.update(**params)
print(model.loglikelihood)
print(model.loglr)

-499356.0809044743
327.5070887125912
In [6]:
# get the waveforms
waveforms = model.get_waveforms()
gated_waveforms = model.get_gated_waveforms()
# whiten the gated waveforms
ghps = {ifo: wf[0] for ifo, wf in gated_waveforms.items()}
ghcs = {ifo: wf[1] for ifo, wf in gated_waveforms.items()}
ghps = model.whiten(ghps, 1)
ghcs = model.whiten(ghcs, 1)
white_gated_waveforms = {ifo: (ghps[ifo], ghcs[ifo]) for ifo in ghps}
In [7]:
# the maxL waveform
hp, hc = waveforms['L1']
hp = hp.to_timeseries()
fig, ax = pyplot.subplots()
ax.plot(hp.sample_times-model._current_params['tc'], hp)
ax.set_ylabel('$h_{+}$')
ax.set_xlabel('time - tc')
ax.set_xlim(-0.2, 0.4)
fig.show()
In [8]:
# frequency domain
fig, ax = pyplot.subplots()
ax.loglog()
hp, hc = waveforms['L1']
ax.plot(hp.sample_frequencies, abs(hp), label='220+221')
psd = model.psds['L1']
ax.plot(psd.sample_frequencies, psd**0.5, label='ASD')
ax.set_ylabel('$h_{+}$')
ax.set_xlabel('frequency')
ax.legend()
fig.show()
#fig.savefig('fdomain-220_221-nohighpass.png', dpi=200)
In [9]:
# the whitened gated maxL waveform
hp, hc = white_gated_waveforms['L1']
hp = hp.to_timeseries()
fig, ax = pyplot.subplots()
ax.plot(hp.sample_times-model._current_params['tc'], hp)
ax.set_ylabel('whitened gated $h_{+}$')
ax.set_xlabel('time - tc')
ax.set_xlim(-0.2, 0.4)
fig.show()

Problematic 221 - with highpass¶

In [10]:
# load the model
fn = '/work/cdcapano/projects/gated_gaussian/gw190521/injections/imrdraws/margpol-inj100/inj03/samples-220_221-18.hdf'
fp = io.loadfile(fn, 'r')
cp = fp.read_config_file()
psds = fp.read_psds()
data = fp.read_data()
model = models.read_from_config(cp, data=data, psds=psds)

2022-06-03 13:59:07,286 Setting up priors for each parameter
2022-06-03 13:59:07,288 No sampling_params section read from config file
2022-06-03 13:59:07,289 Loading waveform transforms
2022-06-03 13:59:07,295 Determining analysis times to use
2022-06-03 13:59:07,296 Padding H1 analysis start and end times by 4 (= psd-inverse-length/2) seconds to account for PSD wrap around effects.
2022-06-03 13:59:07,297 Padding L1 analysis start and end times by 4 (= psd-inverse-length/2) seconds to account for PSD wrap around effects.
2022-06-03 13:59:07,298 Padding V1 analysis start and end times by 4 (= psd-inverse-length/2) seconds to account for PSD wrap around effects.
2022-06-03 13:59:07,298 Will highpass waveforms at 10.000000 Hz
In [11]:
# get the maxL point
samples = fp.read_samples(list(fp['samples'].keys()))
maxidx = samples.loglikelihood.argmax()
params = {p: samples[p][maxidx] for p in model.variable_params}
model.update(**params)
print(model.loglikelihood)
print(model.loglr)

-499745.52640257764
-61.938398667261936
In [12]:
# get the waveforms
waveforms = model.get_waveforms()
gated_waveforms = model.get_gated_waveforms()
# whiten the gated waveforms
ghps = {ifo: wf[0] for ifo, wf in gated_waveforms.items()}
ghcs = {ifo: wf[1] for ifo, wf in gated_waveforms.items()}
ghps = model.whiten(ghps, 1)
ghcs = model.whiten(ghcs, 1)
white_gated_waveforms = {ifo: (ghps[ifo], ghcs[ifo]) for ifo in ghps}
In [13]:
# the maxL waveform
hp, hc = waveforms['L1']
hp = hp.to_timeseries()
fig, ax = pyplot.subplots()
ax.plot(hp.sample_times-model._current_params['tc'], hp)
ax.set_ylabel('$h_{+}$')
ax.set_xlabel('time - tc')
ax.set_xlim(-0.2, 0.4)
fig.show()
In [14]:
# frequency domain
fig, ax = pyplot.subplots()
ax.loglog()
hp, hc = waveforms['L1']
ax.plot(hp.sample_frequencies, abs(hp), label='220+221')
psd = model.psds['L1']
ax.plot(psd.sample_frequencies, psd**0.5, label='ASD')
ax.set_ylabel('$h_{+}$')
ax.set_xlabel('frequency')
ax.legend()
fig.show()
#fig.savefig('fdomain-220_221-withhighpass.png', dpi=200)
In [15]:
# the whitened gated maxL waveform
hp, hc = white_gated_waveforms['L1']
hp = hp.to_timeseries()
fig, ax = pyplot.subplots()
ax.plot(hp.sample_times-model._current_params['tc'], hp)
ax.set_ylabel('whitened gated $h_{+}$')
ax.set_xlabel('time - tc')
ax.set_xlim(-0.2, 0.4)
fig.show()

Ok 220+330 - without highpass¶

In [16]:
# load the model
fn = '/work/cdcapano/projects/gated_gaussian/gw190521/injections/imrdraws/margpol-inj100/inj03/samples-220_330-18.hdf'
fp = io.loadfile(fn, 'r')
cp = fp.read_config_file()
psds = fp.read_psds()
data = fp.read_data()
model = models.read_from_config(cp, data=data, psds=psds, highpass_waveforms=False)

2022-06-03 13:59:51,026 Setting up priors for each parameter
2022-06-03 13:59:51,030 No sampling_params section read from config file
2022-06-03 13:59:51,031 Loading waveform transforms
2022-06-03 13:59:51,037 Determining analysis times to use
2022-06-03 13:59:51,038 Padding H1 analysis start and end times by 4 (= psd-inverse-length/2) seconds to account for PSD wrap around effects.
2022-06-03 13:59:51,040 Padding L1 analysis start and end times by 4 (= psd-inverse-length/2) seconds to account for PSD wrap around effects.
2022-06-03 13:59:51,041 Padding V1 analysis start and end times by 4 (= psd-inverse-length/2) seconds to account for PSD wrap around effects.
In [17]:
# get the maxL point
samples = fp.read_samples(list(fp['samples'].keys()))
maxidx = samples.loglikelihood.argmax()
params = {p: samples[p][maxidx] for p in model.variable_params}
model.update(**params)
print(model.loglikelihood)
print(model.loglr)

-499665.78527988895
17.802699715364724
In [18]:
# get the waveforms
waveforms = model.get_waveforms()
gated_waveforms = model.get_gated_waveforms()
# whiten the gated waveforms
ghps = {ifo: wf[0] for ifo, wf in gated_waveforms.items()}
ghcs = {ifo: wf[1] for ifo, wf in gated_waveforms.items()}
ghps = model.whiten(ghps, 1)
ghcs = model.whiten(ghcs, 1)
white_gated_waveforms = {ifo: (ghps[ifo], ghcs[ifo]) for ifo in ghps}
In [19]:
# the maxL waveform
hp, hc = waveforms['L1']
hp = hp.to_timeseries()
fig, ax = pyplot.subplots()
ax.plot(hp.sample_times-model._current_params['tc'], hp)
ax.set_ylabel('$h_{+}$')
ax.set_xlabel('time - tc')
ax.set_xlim(-0.2, 0.4)
fig.show()
In [20]:
# frequency domain
fig, ax = pyplot.subplots()
ax.loglog()
hp, hc = waveforms['L1']
ax.plot(hp.sample_frequencies, abs(hp), label='220+330')
psd = model.psds['L1']
ax.plot(psd.sample_frequencies, psd**0.5, label='ASD')
ax.set_ylabel('$h_{+}$')
ax.set_xlabel('frequency')
ax.legend()
fig.show()
In [21]:
# the whitened gated maxL waveform
hp, hc = white_gated_waveforms['L1']
hp = hp.to_timeseries()
fig, ax = pyplot.subplots()
ax.plot(hp.sample_times-model._current_params['tc'], hp)
ax.set_ylabel('whitened gated $h_{+}$')
ax.set_xlabel('time - tc')
ax.set_xlim(-0.2, 0.4)
fig.show()

Ok 220+330 - with highpass¶

In [22]:
# load the model
fn = '/work/cdcapano/projects/gated_gaussian/gw190521/injections/imrdraws/margpol-inj100/inj03/samples-220_330-18.hdf'
fp = io.loadfile(fn, 'r')
cp = fp.read_config_file()
psds = fp.read_psds()
data = fp.read_data()
model = models.read_from_config(cp, data=data, psds=psds)

2022-06-03 14:00:14,654 Setting up priors for each parameter
2022-06-03 14:00:14,656 No sampling_params section read from config file
2022-06-03 14:00:14,656 Loading waveform transforms
2022-06-03 14:00:14,667 Determining analysis times to use
2022-06-03 14:00:14,668 Padding H1 analysis start and end times by 4 (= psd-inverse-length/2) seconds to account for PSD wrap around effects.
2022-06-03 14:00:14,669 Padding L1 analysis start and end times by 4 (= psd-inverse-length/2) seconds to account for PSD wrap around effects.
2022-06-03 14:00:14,670 Padding V1 analysis start and end times by 4 (= psd-inverse-length/2) seconds to account for PSD wrap around effects.
2022-06-03 14:00:14,671 Will highpass waveforms at 10.000000 Hz
In [23]:
# get the maxL point
samples = fp.read_samples(list(fp['samples'].keys()))
maxidx = samples.loglikelihood.argmax()
params = {p: samples[p][maxidx] for p in model.variable_params}
model.update(**params)
print(model.loglikelihood)
print(model.loglr)

-499665.7821276057
17.805856003251392
In [24]:
# get the waveforms
waveforms = model.get_waveforms()
gated_waveforms = model.get_gated_waveforms()
# whiten the gated waveforms
ghps = {ifo: wf[0] for ifo, wf in gated_waveforms.items()}
ghcs = {ifo: wf[1] for ifo, wf in gated_waveforms.items()}
ghps = model.whiten(ghps, 1)
ghcs = model.whiten(ghcs, 1)
white_gated_waveforms = {ifo: (ghps[ifo], ghcs[ifo]) for ifo in ghps}
In [25]:
# the maxL waveform
hp, hc = waveforms['L1']
hp = hp.to_timeseries()
fig, ax = pyplot.subplots()
ax.plot(hp.sample_times-model._current_params['tc'], hp)
ax.set_ylabel('$h_{+}$')
ax.set_xlabel('time - tc')
ax.set_xlim(-0.2, 0.4)
fig.show()
In [26]:
# frequency domain
fig, ax = pyplot.subplots()
ax.loglog()
hp, hc = waveforms['L1']
ax.plot(hp.sample_frequencies, abs(hp), label='220+330')
psd = model.psds['L1']
ax.plot(psd.sample_frequencies, psd**0.5, label='ASD')
ax.set_ylabel('$h_{+}$')
ax.set_xlabel('frequency')
ax.legend()
fig.show()
In [27]:
# the whitened gated maxL waveform
hp, hc = white_gated_waveforms['L1']
hp = hp.to_timeseries()
fig, ax = pyplot.subplots()
ax.plot(hp.sample_times-model._current_params['tc'], hp)
ax.set_ylabel('whitened gated $h_{+}$')
ax.set_xlabel('time - tc')
ax.set_xlim(-0.2, 0.4)
fig.show()

Timing¶

In [28]:
def timemodel(model, params):
    model.update(**params)
    model.loglikelihood
In [29]:
# Baseline: no highpass
model = models.read_from_config(cp, data=data, psds=psds, highpass_waveforms=False)

2022-06-03 14:00:28,762 Setting up priors for each parameter
2022-06-03 14:00:28,764 No sampling_params section read from config file
2022-06-03 14:00:28,765 Loading waveform transforms
2022-06-03 14:00:28,771 Determining analysis times to use
2022-06-03 14:00:28,772 Padding H1 analysis start and end times by 4 (= psd-inverse-length/2) seconds to account for PSD wrap around effects.
2022-06-03 14:00:28,773 Padding L1 analysis start and end times by 4 (= psd-inverse-length/2) seconds to account for PSD wrap around effects.
2022-06-03 14:00:28,773 Padding V1 analysis start and end times by 4 (= psd-inverse-length/2) seconds to account for PSD wrap around effects.
In [30]:
%timeit timemodel(model, params)

314 ms ± 13.6 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
In [31]:
# with highpass
model = models.read_from_config(cp, data=data, psds=psds)

2022-06-03 14:00:45,562 Setting up priors for each parameter
2022-06-03 14:00:45,565 No sampling_params section read from config file
2022-06-03 14:00:45,567 Loading waveform transforms
2022-06-03 14:00:45,572 Determining analysis times to use
2022-06-03 14:00:45,573 Padding H1 analysis start and end times by 4 (= psd-inverse-length/2) seconds to account for PSD wrap around effects.
2022-06-03 14:00:45,574 Padding L1 analysis start and end times by 4 (= psd-inverse-length/2) seconds to account for PSD wrap around effects.
2022-06-03 14:00:45,575 Padding V1 analysis start and end times by 4 (= psd-inverse-length/2) seconds to account for PSD wrap around effects.
2022-06-03 14:00:45,576 Will highpass waveforms at 10.000000 Hz
In [32]:
%timeit timemodel(model, params)

353 ms ± 21.5 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
In [ ]: