Analyzing GW170817¶
We will demonstrate how to use jim to analyze the binary neutron star GW170817 using the IMRPhenomD waveform.
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
import numpy as np
import jax.numpy as jnp
import jax
from gwpy.timeseries import TimeSeries
from gwpy.frequencyseries import FrequencySeries
import requests
from astropy.time import Time
from scipy.signal.windows import tukey
from scipy.interpolate import interp1d
from ripple.waveforms.IMRPhenomD import gen_IMRPhenomD_polar
from jaxgw.PE.detector_preset import *
from jaxgw.PE.heterodyneLikelihood import make_heterodyne_likelihood_mutliple_detector
from jaxgw.PE.detector_projection import make_detector_response
from flowMC.nfmodel.rqSpline import RQSpline
from flowMC.sampler.MALA import MALA
from flowMC.sampler.Sampler import Sampler
from flowMC.utils.PRNG_keys import initialize_rng_keys
from flowMC.nfmodel.utils import *
2023-03-30 14:13:39.741712: W external/org_tensorflow/tensorflow/tsl/platform/default/dso_loader.cc:66] Could not load dynamic library 'libcudart.so.11.0'; dlerror: libcudart.so.11.0: cannot open shared object file: No such file or directory; LD_LIBRARY_PATH: /mnt/sw/nix/store/wxp5xscxcqq0l1nlrv8k136qs5wqaln6-vscode-1.73.1/lib:/mnt/sw/nix/store/hayjz1l94cb2ky37bhcv71aygjzq7fci-openblas-0.3.21/lib:/cm/shared/apps/slurm/current/lib64:/run/opengl-driver/lib 2023-03-30 14:13:40.065011: W external/org_tensorflow/tensorflow/tsl/platform/default/dso_loader.cc:66] Could not load dynamic library 'libcudart.so.11.0'; dlerror: libcudart.so.11.0: cannot open shared object file: No such file or directory; LD_LIBRARY_PATH: /mnt/sw/nix/store/wxp5xscxcqq0l1nlrv8k136qs5wqaln6-vscode-1.73.1/lib:/mnt/sw/nix/store/hayjz1l94cb2ky37bhcv71aygjzq7fci-openblas-0.3.21/lib:/cm/shared/apps/slurm/current/lib64:/run/opengl-driver/lib 2023-03-30 14:13:40.076833: W external/org_tensorflow/tensorflow/tsl/platform/default/dso_loader.cc:66] Could not load dynamic library 'libcudart.so.11.0'; dlerror: libcudart.so.11.0: cannot open shared object file: No such file or directory; LD_LIBRARY_PATH: /mnt/sw/nix/store/wxp5xscxcqq0l1nlrv8k136qs5wqaln6-vscode-1.73.1/lib:/mnt/sw/nix/store/hayjz1l94cb2ky37bhcv71aygjzq7fci-openblas-0.3.21/lib:/cm/shared/apps/slurm/current/lib64:/run/opengl-driver/lib
The Kernel crashed while executing code in the the current cell or a previous cell. Please review the code in the cell(s) to identify a possible cause of the failure. Click <a href='https://aka.ms/vscodeJupyterKernelCrash'>here</a> for more info. View Jupyter <a href='command:jupyter.viewOutput'>log</a> for further details.
Canceled future for execute_request message before replies were done
Strain¶
To do so, we need to know the GPS time associated with the event (in this case, $t = 1187008882.43 s$). We also need to prescribe how much data we wish to analyze around the event (in this case, $T = 128 s$, aka, the segment length or seglen). We will place the trigger $2 s$ before the end of the analysis segment, following the LVK convention.
👉 NOTE: if you don't know the tigger GPS time, you may obtain it from the event name using the
datasets.event_gps
utility from the gwosc package, e.g.,event_gps("GW170817")
.
trigger_time = 1187008882.43
seglen = 128
# determine segment bounds, placing trigger 2s before the end
post_trigger_duration = 2
start = trigger_time - seglen + post_trigger_duration
end = trigger_time + post_trigger_duration
With those parameters, we can now fetch the data from GWOSC using fetch_open_data()
. For GW170817, We make sure to specify version=2
to get the version of data without the glitch in Livingston (see GWOSC docs for this release).
ifos = ['H1', 'L1', 'V1']
data_td_dict = {i: TimeSeries.fetch_open_data(i, start, end, version=2)
for i in ifos}
For the likelihood computation, we will want frequency domain data. We can IFFT the above data after applying a window function; following common LVK practice for this event, we apply a Tukey window with a slope parameter alpha=0.00625
.
👉 NOTE: different
alpha
values may be appropriate for different events, e.g.,alpha = 0.4
is standard for shorter binary black holes.
tukey_alpha = 0.00625
data_fd_dict = {}
for ifo, d in data_td_dict.items():
w = tukey(len(d), tukey_alpha)
f = np.fft.rfftfreq(len(d), d=d.dt)
data_fd_dict[ifo] = FrequencySeries(np.fft.rfft(d*w)/d.dt, frequencies=f)
Power spectral densities (PSDs)¶
Besides the strain, to compute the likelihood we will need a PSDs characterizing the noise at each detector. Although we could estimate this oursevles directly from the data (e.g., arXiv:1907.06540), we will forgo that step and download precomputed PSDs made available by the LVK collaboration in LIGO-P1800061.
👉 NOTE: you may load any PSD you wish for this step, whether from disk or computed on the fly.
psd_url = "https://dcc.ligo.org/public/0150/P1800061/011/GW170817_PSDs.dat"
with requests.get(psd_url) as r:
psd_data = np.genfromtxt(r.iter_lines())
The psd_data
object is a 2D array where the first column is frequency and the rest are the corresponding PSD values for H1, L1 and V1, in that order. For convenience, and because these PSD data are not uniformly sampled, we will turn this into interpolants that we can evaluate over any frequency bins for each detector.
psd_dict = {}
for i, (ifo, d) in enumerate(data_fd_dict.items()):
p = interp1d(psd_data[:,0], psd_data[:,i+1], bounds_error=False,
fill_value=np.inf)
psd_dict[ifo] = FrequencySeries(p(d.frequencies), frequencies=d.frequencies)
Forming the likelihood¶
from jimgw.PE.detector_preset import *
from jimgw.PE.heterodyneLikelihood import make_heterodyne_likelihood_mutliple_detector
from jimgw.PE.detector_projection import make_detector_response
H1 = get_H1()
H1_response = make_detector_response(H1[0], H1[1])
L1 = get_L1()
L1_response = make_detector_response(L1[0], L1[1])
V1 = get_V1()
V1_response = make_detector_response(V1[0], V1[1])
def LogLikelihood(theta):
theta = theta.at[1].set(theta[1]/(1+theta[1])**2) # convert q to eta
theta = theta.at[7].set(jnp.arccos(theta[7])) # convert cos iota to iota
theta = theta.at[10].set(jnp.arcsin(theta[10])) # convert cos dec to dec
theta_waveform = theta[:8]
theta_waveform = theta_waveform.at[5].set(0)
ra = theta[9]
dec = theta[10]
hp_test, hc_test = gen_IMRPhenomD_polar(H1_frequency, theta_waveform, f_ref)
align_time = jnp.exp(-1j*2*jnp.pi*H1_frequency*(epoch+theta[5]))
h_test_H1 = H1_response(H1_frequency, hp_test, hc_test, ra, dec, gmst, theta[8]) * align_time
h_test_L1 = L1_response(L1_frequency, hp_test, hc_test, ra, dec, gmst, theta[8]) * align_time
h_test_V1 = V1_response(V1_frequency, hp_test, hc_test, ra, dec, gmst, theta[8]) * align_time
df = H1_frequency[1] - H1_frequency[0]
match_filter_SNR_H1 = 4*jnp.sum((jnp.conj(h_test_H1)*H1_data)/H1_psd*df).real
match_filter_SNR_L1 = 4*jnp.sum((jnp.conj(h_test_L1)*L1_data)/L1_psd*df).real
match_filter_SNR_V1 = 4*jnp.sum((jnp.conj(h_test_V1)*V1_data)/V1_psd*df).real
optimal_SNR_H1 = 4*jnp.sum((jnp.conj(h_test_H1)*h_test_H1)/H1_psd*df).real
optimal_SNR_L1 = 4*jnp.sum((jnp.conj(h_test_L1)*h_test_L1)/L1_psd*df).real
optimal_SNR_V1 = 4*jnp.sum((jnp.conj(h_test_V1)*h_test_V1)/V1_psd*df).
<ufunc 'degrees'>