bbhx Tutorial
bbhx is a software package that produces black hole binary waveforms. It focuses on LISA and provides the proper LISA response function for MBHBs (arXiv:1806.10734, arXiv:2003.00357). bbhx also provides fast likelihood functions. The package is GPU-accelerated for fast analysis.
If you use this software please cite arXiv:2005.01827, arXiv:2111.01064, and the associated Zenodo page Please also cite any consituent parts used like the response function or waveforms. See the citation attribute for each class or docstring for functions for more information.
[1]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from bbhx.waveformbuild import BBHWaveformFD
from bbhx.waveforms.phenomhm import PhenomHMAmpPhase
from bbhx.response.fastfdresponse import LISATDIResponse
from bbhx.likelihood import Likelihood, HeterodynedLikelihood
from bbhx.utils.constants import *
from bbhx.utils.transform import *
from lisatools.sensitivity import get_sensitivity
np.random.seed(111222)
No CuPy
No CuPy or GPU PhenomHM module.
No CuPy or GPU interpolation available.
No CuPy or GPU response available.
No CuPy
GPU accelerated MBHB waveforms for LISA
First, we will detail how to produce full waveforms for MBHBs to be detected by LISA. We will look at how to generate PhenomHM waveforms (arXiv:1708.00404, arXiv:1508.07253, arXiv:1508.07250) put through the LISA response function ((arXiv:1806.10734, arXiv:2003.00357). More information on generating the waveforms and response separately can be found below.
[2]:
wave_gen = BBHWaveformFD(amp_phase_kwargs=dict(run_phenomd=False))
[3]:
# set parameters
f_ref = 0.0 # let phenom codes set f_ref -> fmax = max(f^2A(f))
phi_ref = 0.0 # phase at f_ref
m1 = 1e6
m2 = 5e5
a1 = 0.2
a2 = 0.4
dist = 18e3 * PC_SI * 1e6 # 3e3 in Mpc
inc = np.pi/3.
beta = np.pi/4. # ecliptic latitude
lam = np.pi/5. # ecliptic longitude
psi = np.pi/6. # polarization angle
t_ref = 0.5 * YRSID_SI # t_ref (in the SSB reference frame)
# frequencies to interpolate to
freq_new = np.logspace(-4, 0, 10000)
modes = [(2,2), (2,1), (3,3), (3,2), (4,4), (4,3)]
[4]:
wave = wave_gen(m1, m2, a1, a2,
dist, phi_ref, f_ref, inc, lam,
beta, psi, t_ref, freqs=freq_new,
modes=modes, direct=False, fill=True, squeeze=True, length=1024)[0]
for i, let in enumerate(["A", "E", "T"]):
plt.loglog(freq_new, np.abs(wave[i]), label=let)
plt.legend()
plt.xlabel("Frequency (Hz)")
plt.ylabel(r"$\tilde{h}(f)$ (Hz$^{-1/2}$)")
[4]:
Text(0, 0.5, '$\\tilde{h}(f)$ (Hz$^{-1/2}$)')
Adjust inspiral timing
[5]:
# start at the beginning of observation
t_obs_start = 0.0
t_obs_end = (t_ref - 24 * 3600) / YRSID_SI # 1 day before merger
wave = wave_gen(m1, m2, a1, a2,
dist, phi_ref, f_ref, inc, lam,
beta, psi, t_ref, freqs=freq_new,
modes=modes, direct=False, fill=True, squeeze=True, length=1024, t_obs_start=t_obs_start, t_obs_end=t_obs_end)[0]
for i, let in enumerate(["A", "E", "T"]):
plt.loglog(freq_new, np.abs(wave[i]), label=let)
plt.legend()
plt.xlabel("Frequency (Hz)")
plt.ylabel(r"$\tilde{h}(f)$ (Hz$^{-1/2}$)")
[5]:
Text(0, 0.5, '$\\tilde{h}(f)$ (Hz$^{-1/2}$)')
Controlling interpolation
The user has the option to control the interpolation by providing the number of frequencies to interpolate over with the length keyword argument.
[6]:
# number of points is too small
length = 128
wave1 = wave_gen(m1, m2, a1, a2,
dist, phi_ref, f_ref, inc, lam,
beta, psi, t_ref, freqs=freq_new,
modes=modes, direct=False, fill=True, squeeze=True, length=length,
)[0]
# high number of points is good
length = 16384
wave2 = wave_gen(m1, m2, a1, a2,
dist, phi_ref, f_ref, inc, lam,
beta, psi, t_ref, freqs=freq_new,
modes=modes, direct=False, fill=True, squeeze=True, length=length,
)[0]
fig, ax = plt.subplots(1, 2)
plt.subplots_adjust(wspace=0.4)
for i, let in enumerate(["A", "E", "T"]):
ax[0].loglog(freq_new, np.abs(wave1[i]), label=let)
ax[1].loglog(freq_new, np.abs(wave2[i]), label=let)
ax[0].legend()
ax[0].set_ylabel(r"$\tilde{h}(f)$ (Hz$^{-1/2}$)")
ax[0].set_xlabel("Frequency (Hz)")
ax[1].set_xlabel("Frequency (Hz)")
[6]:
Text(0.5, 0, 'Frequency (Hz)')
Check number of points that is good for interpolation:
Generally, most reasonable numbers of points are fine for interpolation. We usually put 1024 to be conservative. It can slightly affect the waveform timing.
[7]:
# compute directly first (see below)
wave = wave_gen(m1, m2, a1, a2,
dist, phi_ref, f_ref, inc, lam,
beta, psi, t_ref, freqs=freq_new,
modes=modes, direct=True, compress=True, squeeze=False,
) [0]
lengths_in = 2**np.arange(5, 15)[::-1]
mismatch = np.zeros_like(lengths_in, dtype=float)
for j, length in enumerate(lengths_in):
# compute directly first (see below)
wave_temp = wave_gen(m1, m2, a1, a2,
dist, phi_ref, f_ref, inc, lam,
beta, psi, t_ref, freqs=freq_new,
modes=modes,
direct=False, fill=True, squeeze=False, length=length
)[0]
num = np.sum([np.dot(wave_temp[i].conj(), wave[i]) for i in range(3)])
den1 = np.sum([np.dot(wave[i].conj(), wave[i]) for i in range(3)])
den2 = np.sum([np.dot(wave_temp[i].conj(), wave_temp[i]) for i in range(3)])
overlap = num / np.sqrt(den1 * den2)
mismatch[j] = 1 - overlap.real
plt.loglog(lengths_in, mismatch)
plt.xlabel("Number of initial frequencies")
plt.ylabel("Mismatch")
[7]:
Text(0, 0.5, 'Mismatch')
Generate waveforms without interpolation
When interpolating, frequencies between \(10^{-4}/M\) and \(0.6/M\) are evaluated. Therefore, when not interpolating, the evaluation of the waveform will stretch outside these limits if frequencies outside those limits are given by the user.
With compress=True, all modes will be directly combined.
[8]:
wave = wave_gen(m1, m2, a1, a2,
dist, phi_ref, f_ref, inc, lam,
beta, psi, t_ref, freqs=freq_new,
modes=modes, direct=True, compress=True, squeeze=True,
)
for i, let in enumerate(["A", "E", "T"]):
plt.loglog(freq_new, np.abs(wave[i]), label=let)
plt.legend()
plt.xlabel("Frequency (Hz)")
plt.ylabel(r"$\tilde{h}(f)$ (Hz$^{-1/2}$)")
[8]:
Text(0, 0.5, '$\\tilde{h}(f)$ (Hz$^{-1/2}$)')
With compress=False, it will keep the modes separate.
[9]:
wave = wave_gen(m1, m2, a1, a2,
dist, phi_ref, f_ref, inc, lam,
beta, psi, t_ref, freqs=freq_new,
modes=modes, direct=True, compress=False, squeeze=True,
)
for i, mode in enumerate(wave_gen.amp_phase_gen.modes):
plt.loglog(freq_new, np.abs(wave[0][i]), label=mode)
plt.legend()
plt.xlabel("Frequency (Hz)")
plt.ylabel(r"$\tilde{h}(f)$ (Hz$^{-1/2}$)")
plt.xlim(1e-4, 1e-1)
plt.ylim(1e-27, 1e-14)
[9]:
(1e-27, 1e-14)
Adjusting mode content
You can adjust the mode content with the modes keyword argument. You will receive an error if the requested mode indices \((l,m)\) are not available for the given waveform model chosen.
[10]:
wave = wave_gen(m1, m2, a1, a2,
dist, phi_ref, f_ref, inc, lam,
beta, psi, t_ref, freqs=freq_new,
modes=[(3,3),(4,4)], direct=True, compress=False, squeeze=True,
)
for i, mode in enumerate(wave_gen.amp_phase_gen.modes):
plt.loglog(freq_new, np.abs(wave[0][i]), label=mode)
plt.legend()
plt.xlabel("Frequency (Hz)")
plt.ylabel(r"$\tilde{h}(f)$ (Hz$^{-1/2}$)")
plt.xlim(1e-4, 1e-1)
plt.ylim(1e-27, 1e-14)
[10]:
(1e-27, 1e-14)
All of these functions can take arrays for parameters
[11]:
# set parameters
f_ref = np.array([0.0, 0.0]) # let phenom codes set fRef -> fmax = max(f^2A(f))
phi_ref = np.array([0.0, 1.0])
m1 = np.array([1e6, 4e5])
m2 = np.array([5e5, 1e5])
a1 = np.array([0.2, 0.8])
a2 = np.array([0.4, 0.7])
dist = np.array([10e3, 1e3]) * PC_SI * 1e6 # 3e3 in Mpc
inc = np.array([np.pi/3, np.pi/4])
beta = np.array([np.pi/4., np.pi/5])
lam = np.array([np.pi/5., np.pi/6])
psi = np.array([np.pi/6., np.pi/7])
t_ref = np.array([0.5, 1.2]) * YRSID_SI # in the SSB reference frame
# frequencies to interpolate to
freq_new = np.logspace(-4, 0, 10000)
modes = [(2,2), (2,1), (3,3), (3,2), (4,4), (4,3)]
waves = wave_gen(m1, m2, a1, a2,
dist, phi_ref, f_ref, inc, lam,
beta, psi, t_ref, freqs=freq_new,
modes=modes, direct=False, fill=True, squeeze=False, length=1024)
for j, ls in enumerate(["solid", "dashed"]):
for i, let in enumerate(["A", "E", "T"]):
plt.loglog(freq_new, np.abs(waves[j][i]), label=let + f" bin {j}", ls=ls)
plt.legend()
plt.xlabel("Frequency (Hz)")
plt.ylabel(r"$\tilde{h}(f)$ (Hz$^{-1/2}$)")
[11]:
Text(0, 0.5, '$\\tilde{h}(f)$ (Hz$^{-1/2}$)')
The fill keyword argument
The fill keyword argument allows the user to fill the data stream as in the first example above. However, it can give two other returns. (It is only available when direct=False.
If fill=False, it will return special information used in the fast likelihood function. This includes templates as shortened arrays (only the values where \(\tilde{h}(f)\) is non-zero), the information on which frequency each array starts at and its length.
[12]:
templates, start_inds, lengths = wave_gen(m1, m2, a1, a2,
dist, phi_ref, f_ref, inc, lam,
beta, psi, t_ref, freqs=freq_new,
modes=modes, direct=False, fill=False, squeeze=False, length=1024)
print(templates, start_inds, lengths)
# plot the first binary
for i, let in enumerate(["A", "E", "T"]):
plt.loglog(freq_new[start_inds[0]:start_inds[0] + lengths[0]], np.abs(templates[0][i]), label=let)
plt.legend()
plt.xlabel("Frequency (Hz)")
plt.ylabel(r"$\tilde{h}(f)$ (Hz$^{-1/2}$)")
[array([[-3.14154099e-19+2.93930145e-19j, -4.11910745e-19-2.31073392e-19j,
7.18639789e-20-4.30372868e-19j, ...,
-1.21330575e-43-4.29769485e-44j, 6.36570618e-46-1.07049893e-43j,
-1.51012483e-45-8.87380319e-44j],
[-3.77142996e-19+6.09479474e-19j, -7.23835376e-19-1.45368093e-19j,
-9.26662932e-20-7.14586769e-19j, ...,
2.57629108e-43-8.19014428e-44j, 1.33780103e-43+1.81800785e-43j,
1.11242618e-43+1.49366094e-43j],
[-2.82108635e-28+4.25544336e-28j, -5.89161226e-28+6.21183363e-30j,
-9.93313342e-29-5.74008671e-28j, ...,
-3.89107774e-44-4.92386953e-44j, 2.82041850e-44-4.37426151e-44j,
2.24209609e-44-3.62845654e-44j]]), array([[-3.07245461e-20+6.26998385e-21j, 2.18054463e-20-2.25976697e-20j,
-4.29271950e-21+3.11543715e-20j, ...,
-2.78259785e-20-8.22146447e-20j, 7.77498602e-20+3.89050092e-20j,
-7.88818912e-20+3.68982312e-20j],
[ 1.30420115e-20+4.61283100e-22j, -1.09311648e-20+7.26757282e-21j,
4.32441893e-21-1.24795359e-20j, ...,
9.23185723e-20+2.51389976e-20j, -8.34510053e-20+4.74818772e-20j,
1.43365804e-20-9.52734595e-20j],
[ 3.26882976e-30+8.96343560e-29j, -5.51509029e-29-7.21135974e-29j,
8.84898291e-29+2.47851396e-29j, ...,
-2.76775148e-26+4.58288676e-26j, -1.27450933e-26-5.22741814e-26j,
4.91910732e-26+2.24536166e-26j]])] [0 0] [8026 1681]
[12]:
Text(0, 0.5, '$\\tilde{h}(f)$ (Hz$^{-1/2}$)')
The binaries can also be combined into one data stream with fill=True and combine=True.
[13]:
waves = wave_gen(m1, m2, a1, a2,
dist, phi_ref, f_ref, inc, lam,
beta, psi, t_ref, freqs=freq_new,
modes=modes, direct=False, fill=True, combine=True, squeeze=False, length=1024)
for i, let in enumerate(["A", "E", "T"]):
plt.loglog(freq_new, np.abs(waves[i]), label=let)
plt.legend()
plt.xlabel("Frequency (Hz)")
plt.ylabel(r"$\tilde{h}(f)$ (Hz$^{-1/2}$)")
[13]:
Text(0, 0.5, '$\\tilde{h}(f)$ (Hz$^{-1/2}$)')
Generating PhenomHM waveforms
Now we will create PhenomHM waveforms in the source frame, but scaled for the distance.
[14]:
phenomhm = PhenomHMAmpPhase(use_gpu=False, run_phenomd=False)
[15]:
f_ref = 0.0 # let phenom codes set f_ref -> fmax = max(f^2A(f))
phi_ref = 0.0 # phase at f_ref
m1 = 1e6
m2 = 5e5
a1 = 0.2
a2 = 0.4
dist = 18e3 * PC_SI * 1e6 # 3e3 in Mpc
t_ref = 0.95 * YRSID_SI
modes = [(2,2), (3, 3), (4, 4), (2, 1), (3, 2), (4, 3)]
phenomhm(
m1,
m2,
a1,
a2,
dist,
phi_ref,
f_ref,
t_ref,
1024,
modes=modes
)
# get important quantities
freqs = phenomhm.freqs_shaped # shape (num_bin_all, num_modes, length)
amps = phenomhm.amp # shape (num_bin_all, num_modes, length)
phase = phenomhm.phase # shape (num_bin_all, num_modes, length)
tf = phenomhm.tf # shape (num_bin_all, num_modes, length)
Plot Amplitudes
[16]:
for i, mode in enumerate(modes):
plt.loglog(freqs[0, i], np.sqrt(freqs[0, i]) * amps[0, i], label=mode)
plt.ylabel("Characteristic Strain")
plt.xlabel("Frequency (Hz)")
plt.legend()
[16]:
<matplotlib.legend.Legend at 0x13660aa20>
Fast FD Response
A large amount of the Fast FD Response codes were provided by Sylvain Marsat, so thank you to him! Detailed information on the response can be found in arXiv:1806.10734 and arXiv:2003.00357. The response works by applying an effective transfer function, \(\mathcal{T}(f, t(f))\), to the waveform to get the TDI channels:
where \(t(f)\) is the time-frequency correspondence. This determines where the LISA constellation is in its orbit. This response function assumes the arms of the constellation are constant in time.
The response function operates directly on information coming from the base waveform generation in PhenomHM when you run BBHWaveformFD. To run the response separately you need to provide a phase and tf vector. The phase is adjusted inside the response code. If you do not want that to happen, you can add the keyword argument adjust_phase=False. If you are just running the response function to use it itself and not using to adjust phase, you can just put in zeros for the phase, it will not
affect the computation. All parameters for the response are in the SSB frame. This function can take 1D arrays for the parameters as well. In this case, you need freqs.shape=(num_bin_all, length), phases.shape=tf.shape=(num_bin_all, num_modes, length) with num_bin_all the total number of binaries, num_modes the number of harmonics, and length number of frequencies to evaluate per binary.
[17]:
# use phase/tf information from last waveform run
freqs = phenomhm.freqs.copy()
phase = phenomhm.phase.copy()
tf = phenomhm.tf.copy()
modes = phenomhm.modes
phi_ref = 0.0
inc = np.pi/4
beta = np.pi/5
lam = np.pi/6
psi = np.pi/7
length = phenomhm.freqs_shaped.shape[-1]
freqs_shaped = phenomhm.freqs_shaped
response = LISATDIResponse()
response(freqs,
inc,
lam,
beta,
psi,
phi_ref,
length,
phase=phase,
tf=tf,
modes=modes)
# plot parts of the response
for i in range(1, 4):
# (2, 2) mode
index = response.modes.index((2,2))
# response.transferL1
transfer = getattr(response, f"transferL{i}")
plt.loglog(freqs_shaped[0, index], np.abs(transfer[0, index]), label=f"L{i} Real", ls="solid", color=f"C{i-1}")
#plt.loglog(freqs, transfer.imag[0, index], label=f"L{i} Imag", ls="dashed", color=f"C{i-1}")
plt.legend()
plt.xlabel("Frequency (Hz)")
plt.ylabel(r"$\mathcal{T}(f, t(f))$")
[17]:
Text(0, 0.5, '$\\mathcal{T}(f, t(f))$')
[18]:
# use phase/tf information from last waveform run
num_bin_all = 10
freqs = np.tile(phenomhm.freqs_shaped.copy(), (num_bin_all, 1, 1))
phase = np.tile(phenomhm.phase.copy(), (num_bin_all, 1, 1))
tf = np.tile(phenomhm.tf.copy(), (num_bin_all, 1, 1))
phi_ref = np.full(num_bin_all, 0.0)
inc = np.full(num_bin_all, np.pi/4)
beta = np.full(num_bin_all, np.pi/5)
lam = np.full(num_bin_all, np.pi/6)
psi = np.full(num_bin_all, np.pi/7)
length = phenomhm.freqs_shaped.shape[-1]
response(
freqs,
inc,
lam,
beta,
psi,
phi_ref,
length,
phase=phase,
tf=tf,
modes=None, # defaults to phenomhm modes
adjust_phase=True, # if you want to keep the phase array you input the same, set this to false.
) # You can get adjusted phases with response.phase
print(freqs.shape, phase.shape, tf.shape, response.transferL1.shape)
# get phase adjustment by response function
phase_diff = phase - phenomhm.phase
plt.semilogx(freqs[0, 0], phase_diff[0, 0])
plt.xlabel("Frequency (Hz)")
plt.ylabel("Phase Difference from Response Function")
(10, 6, 1024) (10, 6, 1024) (10, 6, 1024) (10, 6, 1024)
[18]:
Text(0, 0.5, 'Phase Difference from Response Function')
If building a waveform, you can also put a larger buffer object directly into the response code. This larger buffer object would be carrying amplitudes, phases, and tf from the source-frame waveform generator in a flattened 1D array.
[19]:
phenomhm.tf.min()
[19]:
2074.3644391335547
[20]:
num_bin_all = 10
length = 1024
num_modes = len(response.modes)
buffer_shape = (9, num_bin_all, num_modes, length)
out_buffer = np.zeros(buffer_shape).flatten()
f_ref = np.full(num_bin_all, 0.0) # let phenom codes set f_ref -> fmax = max(f^2A(f))
phi_ref = np.full(num_bin_all, 0.0) # phase at f_ref
t_ref = np.full(num_bin_all, 3e6) # time at f_ref
m1 = np.full(num_bin_all, 1e6)
m2 = np.full(num_bin_all, 5e5)
a1 = np.full(num_bin_all, 0.2)
a2 = np.full(num_bin_all, 0.4)
dist = np.full(num_bin_all, 3e3) * PC_SI * 1e6 # 3e3 in Mpc
inc = np.full(num_bin_all, np.pi/4)
beta = np.full(num_bin_all, np.pi/5)
lam = np.full(num_bin_all, np.pi/6)
psi = np.full(num_bin_all, np.pi/7)
phenomhm(
m1,
m2,
a1,
a2,
dist,
phi_ref,
f_ref,
t_ref,
length,
out_buffer=out_buffer
)
print(out_buffer.reshape(buffer_shape)[:, 0, 0, 0])
response(
freqs,
inc,
lam,
beta,
psi,
phi_ref,
length,
out_buffer=out_buffer,
modes=None, # defaults to phenomhm modes
)
print(out_buffer.reshape(buffer_shape)[:, 0, 0, 0])
# these are all numbers at the first frequency for (2,2) mode for
# amp, phase, tf, transferL1re, transferL1im, transferL2re, transferL2im, transferL3re, transferL3im
[1.70824294e-13 4.14100197e+03 1.65814682e+01 0.00000000e+00
0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00
0.00000000e+00]
[ 1.70824294e-13 4.14087439e+03 1.65814682e+01 -3.09383403e-06
2.34941055e-06 4.04192751e-06 1.44557752e-06 2.79724168e-16
-4.52904052e-15]
Computing the Likelihood
The gravitational-wave likelihood, \(\mathcal{L}\), is given by
where \(\langle a|b\rangle\) is the noise-weighted inner product between time-domain strain signals \(a(t)\) and \(b(t)\). This inner product is given by
where \(S_n(f)\) is the one-sided Power Spectral Density in the noise.
Computing the log-Likelihood for MBHBs directly against the Fourier transform of data
One option in bbhx for computing the likelihood is a direct Likeihood computation: bbhx.likelihood.Likelihood. You will need lisatools if you want to compute the sensitivity curves. Otherwise, put in your own sensitivity functions or np.ones to get it to work.
[22]:
######## generate data
# set parameters
f_ref = 0.0 # let phenom codes set f_ref -> fmax = max(f^2A(f))
phi_ref = 0.0 # phase at f_ref
m1 = 1e6
m2 = 5e5
a1 = 0.2
a2 = 0.4
dist = 18e3 * PC_SI * 1e6 # 3e3 in Mpc
inc = np.pi/3.
beta = np.pi/4. # ecliptic latitude
lam = np.pi/5. # ecliptic longitude
psi = np.pi/6. # polarization angle
t_ref = 0.9 * YRSID_SI # t_ref (in the SSB reference frame)
T_obs = 1.2 # years
dt = 10.0
n = int(T_obs * YRSID_SI / dt)
data_freqs = np.fft.rfftfreq(n, dt)[1:] # remove DC
# frequencies to interpolate to
modes = [(2,2), (2,1), (3,3), (3,2), (4,4), (4,3)]
waveform_kwargs = dict(modes=modes, direct=False, fill=True, squeeze=True, length=1024)
data_channels = wave_gen(m1, m2, a1, a2,
dist, phi_ref, f_ref, inc, lam,
beta, psi, t_ref, freqs=data_freqs,
**waveform_kwargs)[0]
######## get noise information (need lisatools)
PSD_A = get_sensitivity(data_freqs, sens_fn="A1TDISens")
PSD_E = get_sensitivity(data_freqs, sens_fn="E1TDISens")
PSD_T = get_sensitivity(data_freqs, sens_fn="T1TDISens")
df = data_freqs[1] - data_freqs[0]
psd = np.array([PSD_A, PSD_E, PSD_T])
# initialize Likelihood
like = Likelihood(
wave_gen,
data_freqs,
data_channels,
psd,
use_gpu=False,
)
# get params
num_bins = 10
params_in = np.tile(np.array([m1, m2, a1, a2, dist, phi_ref, f_ref, inc, lam, beta, psi, t_ref]), (num_bins, 1))
# change masses for test
params_in[:, 0] *= (1 + 1e-4 * np.random.randn(num_bins))
# get_ll and not __call__ to work with lisatools
ll = like.get_ll(params_in.T, **waveform_kwargs)
print(ll, like.d_h)
[ -1.07304677 -9.67441448 -0.82049942 -6.35227414 -10.27292833
-0.18575894 -0.35926129 -2.60369419 -1.14904283 -5.69420711] [764982.31202969 -591.06049369j 764946.22578225-1778.38208222j
765008.20318771 +516.83052761j 764957.38520581-1439.90214805j
764944.36630292-1832.8289125j 764991.18762188 -245.87126849j
765004.61719784 +341.9503657j 764973.14459642 -921.03288148j
765010.06945941 +611.6637902j 764959.82211007-1363.06563164j]
You can also phase marginalize. This occurs in the \(\langle d|h \rangle\) term in the Likelihood.
[23]:
# get_ll and not __call__ to work with lisatools
ll_marg = like.get_ll(params_in.T, phase_marginalize=True, **waveform_kwargs)
print(np.array([ll, ll_marg]).T)
[[ -1.07304677 -0.84470656]
[ -9.67441448 -7.60718515]
[ -0.82049942 -0.64591713]
[ -6.35227414 -4.99708991]
[-10.27292833 -8.07717584]
[ -0.18575894 -0.14624693]
[ -0.35926129 -0.28283688]
[ -2.60369419 -2.04922946]
[ -1.14904283 -0.90451498]
[ -5.69420711 -4.47979925]]
You can also return the signal-to-noise ratio: \(\langle d|h \rangle / \sqrt{\langle h|h \rangle}\).
[24]:
# get_ll and not __call__ to work with lisatools
out = like.get_ll(params_in.T, return_extracted_snr=True, **waveform_kwargs)
# first column is log-Likelihood, second is snr
print(out)
[[-1.07304677e+00 8.74639879e+02]
[-9.67441448e+00 8.74630045e+02]
[-8.20499422e-01 8.74640167e+02]
[-6.35227414e+00 8.74633843e+02]
[-1.02729283e+01 8.74629361e+02]
[-1.85758940e-01 8.74640893e+02]
[-3.59261290e-01 8.74640695e+02]
[-2.60369419e+00 8.74638129e+02]
[-1.14904283e+00 8.74639792e+02]
[-5.69420711e+00 8.74634595e+02]]
Heterodyned Likelihoods
Heterodyned Likelihoods are also available. See arXiv:2109.02728 for details on Heterodyning in general. We specifically implement the method used in arXiv:1806.08792. The practical reason for using a Heterodyned Likelihood method is the computation can be performed along a sparse array of frequencies, rather than a dense array comming from a Fourier transform.
It is reasonably accurate with higher modes and produces the same posteriors as the full Likelihood. It uses a direct sum of modes prior to residual computation. It gets more accurate as modes are removed.
[25]:
reference_index = ll.argmax()
reference_params = params_in[reference_index]
# how many frequencies to use
length_f_het = 128
# initialize Likelihood
like_het = HeterodynedLikelihood(
wave_gen,
data_freqs,
data_channels,
reference_params,
length_f_het,
use_gpu=False,
)
# get_ll and not __call__ to work with lisatools
ll_het = like_het.get_ll(params_in.T, **waveform_kwargs)
print(np.array([ll, ll_het]).T)
[[ -1.07304677 -1.00832284]
[ -9.67441448 -8.90280807]
[ -0.82049942 -0.74911253]
[ -6.35227414 -5.85909352]
[-10.27292833 -9.45085245]
[ -0.18575894 -0.18590962]
[ -0.35926129 -0.32960139]
[ -2.60369419 -2.41820795]
[ -1.14904283 -1.04894937]
[ -5.69420711 -5.25570331]]
Utility functions
There are some utility functions in bbhx than can be very useful.
Transformation functions
When sampling on MBHBs, it can be useful to sample in the LISA constellation reference frame where the posterior modes in the sky have a generally predictable structure. In the LISA constellation referenence frame, the sky localization for a source generally has 8 sky modes: 4 longitudinal and 2 latitudinal. Longitudinal modes involve \(\lambda + (0, 1, 2, 3)\times\pi/2\) and \(\psi + (0, 1, 2, 3)\times\pi/2\), where \(\lambda\) is the ecliptic longitude and \(\psi\) is the polarization angle. The latitudinal modes are at \((\pm\beta, \pm\cos{\iota}, \pm\cos{\psi})\) where \(\beta\) and \(\iota\) are the ecliptic latitude and inclination, respectively.
Transformations to and from LISA and SSB reference frames
First, we will examine transformations of coordinates to and from the LISA and SSB reference frames. All of these functions can take arrays as input as well.
[26]:
# start in the LISA frame
tL = YRSID_SI # time
lambdaL = 0.5
betaL = 0.8
psiL = 0.3
tSSB, lambdaSSB, betaSSB, psiSSB = LISA_to_SSB(tL, lambdaL, betaL, psiL)
print(tSSB, lambdaSSB, betaSSB, psiSSB)
31557992.305078838 2.3277015675208466 1.0933520437673625 2.3134690745168953
From SSB to LISA:
[27]:
tL, lambdaL, betaL, psiL = SSB_to_LISA(tSSB, lambdaSSB, betaSSB, psiSSB)
print(tL, lambdaL, betaL, psiL)
31558149.7635457 0.4999999999999887 0.8000000000000038 0.29999999999999005
Posterior sky-mode transformation
You can transform your sky mode parameters to all longitudinal modes. If inplace==False, the first dimension will be 4 times as long as the input array.
[28]:
inc = 0.4
lam = 3.42
beta = 0.6
psi = 1.0
coords = np.array([inc, lam, beta, psi])
# if we had the whole mbh parameter array like used in BBHWaveformFD,
# it would be dict(inc=7, lam=8, beta=9, psi=10) included f_ref as a parameter
ind_map = dict(inc=0, lam=1, beta=2, psi=3)
print(f"start params: {coords}")
out_coords = mbh_sky_mode_transform(coords, ind_map=ind_map, inplace=False, kind="long", cos_i=False)
print(f"end params: {out_coords}")
start params: [0.4 3.42 0.6 1. ]
end params: [[0.4 3.42 0.6 1. ]
[0.4 4.99079633 0.6 2.57079633]
[0.4 0.27840735 0.6 1. ]
[0.4 1.84920367 0.6 2.57079633]]
Or latitudinal transforms. If inplace==False, the first dimension will be 2 times as long as the input array.
[29]:
print(f"start params: {coords}")
out_coords = mbh_sky_mode_transform(coords, ind_map=ind_map, inplace=False, kind="lat", cos_i=False)
print(f"end params: {out_coords}")
start params: [0.4 3.42 0.6 1. ]
end params: [[ 0.4 3.42 0.6 1. ]
[ 2.74159265 3.42 -0.6 2.14159265]]
Or both transforms. If inplace==False, the first dimension will be 8 times as long as the input array.
[30]:
print(f"start params: {coords}")
out_coords = mbh_sky_mode_transform(coords, ind_map=ind_map, inplace=False, kind="both", cos_i=False)
print(f"end params: {out_coords}")
start params: [0.4 3.42 0.6 1. ]
end params: [[ 0.4 3.42 0.6 1. ]
[ 2.74159265 3.42 -0.6 2.14159265]
[ 0.4 4.99079633 0.6 2.57079633]
[ 2.74159265 4.99079633 -0.6 0.57079633]
[ 0.4 0.27840735 0.6 1. ]
[ 2.74159265 0.27840735 -0.6 2.14159265]
[ 0.4 1.84920367 0.6 2.57079633]
[ 2.74159265 1.84920367 -0.6 0.57079633]]
You can also do these transformations in-place. When you do this, the input binaries will be assigned to new sky modes based on the transformation. Every 2, 4, or 8 binaries will be given the same sky mode transformation if latitudinal, longitudinal, or both transformations are applied, respectively.
[31]:
in_coords = np.tile(np.array([inc, lam, beta, psi]), (16, 1))
print(f"start params: {in_coords}")
out_coords = mbh_sky_mode_transform(in_coords, ind_map=ind_map, inplace=True, kind="both", cos_i=False)
print(f"end params: {out_coords}")
start params: [[0.4 3.42 0.6 1. ]
[0.4 3.42 0.6 1. ]
[0.4 3.42 0.6 1. ]
[0.4 3.42 0.6 1. ]
[0.4 3.42 0.6 1. ]
[0.4 3.42 0.6 1. ]
[0.4 3.42 0.6 1. ]
[0.4 3.42 0.6 1. ]
[0.4 3.42 0.6 1. ]
[0.4 3.42 0.6 1. ]
[0.4 3.42 0.6 1. ]
[0.4 3.42 0.6 1. ]
[0.4 3.42 0.6 1. ]
[0.4 3.42 0.6 1. ]
[0.4 3.42 0.6 1. ]
[0.4 3.42 0.6 1. ]]
end params: [[ 0.4 3.42 0.6 1. ]
[ 2.74159265 3.42 -0.6 2.14159265]
[ 0.4 4.99079633 0.6 2.57079633]
[ 2.74159265 4.99079633 -0.6 0.57079633]
[ 0.4 0.27840735 0.6 1. ]
[ 2.74159265 0.27840735 -0.6 2.14159265]
[ 0.4 1.84920367 0.6 2.57079633]
[ 2.74159265 1.84920367 -0.6 0.57079633]
[ 0.4 3.42 0.6 1. ]
[ 2.74159265 3.42 -0.6 2.14159265]
[ 0.4 4.99079633 0.6 2.57079633]
[ 2.74159265 4.99079633 -0.6 0.57079633]
[ 0.4 0.27840735 0.6 1. ]
[ 2.74159265 0.27840735 -0.6 2.14159265]
[ 0.4 1.84920367 0.6 2.57079633]
[ 2.74159265 1.84920367 -0.6 0.57079633]]
Citations
You can access citations for each class with the citation property.
[32]:
print(wave_gen.citation)
@article{Katz:2020hku,
author = "Katz, Michael L. and Marsat, Sylvain and Chua, Alvin J. K. and Babak, Stanislav and Larson, Shane L.",
title = "{GPU-accelerated massive black hole binary parameter estimation with LISA}",
eprint = "2005.01827",
archivePrefix = "arXiv",
primaryClass = "gr-qc",
doi = "10.1103/PhysRevD.102.023033",
journal = "Phys. Rev. D",
volume = "102",
number = "2",
pages = "023033",
year = "2020"
}
@article{Katz:2021uax,
author = "Katz, Michael L.",
title = "{A fully-automated end-to-end pipeline for massive black hole binary signal extraction from LISA data}",
eprint = "2111.01064",
archivePrefix = "arXiv",
primaryClass = "gr-qc",
month = "11",
year = "2021"
}
@article{michael katz_2021,
title={mikekatz04/BBHx: First official public release},
DOI={10.5281/zenodo.5730688},
abstractNote={No description provided.},
publisher={Zenodo},
author={Michael Katz},
year={2021},
month={Nov}
}
@article{Marsat:2018oam,
author = "Marsat, Sylvain and Baker, John G.",
title = "{Fourier-domain modulations and delays of gravitational-wave signals}",
eprint = "1806.10734",
archivePrefix = "arXiv",
primaryClass = "gr-qc",
month = "6",
year = "2018"
}
@article{Marsat:2020rtl,
author = "Marsat, Sylvain and Baker, John G. and Dal Canton, Tito",
title = "{Exploring the Bayesian parameter estimation of binary black holes with LISA}",
eprint = "2003.00357",
archivePrefix = "arXiv",
primaryClass = "gr-qc",
doi = "10.1103/PhysRevD.103.083011",
journal = "Phys. Rev. D",
volume = "103",
number = "8",
pages = "083011",
year = "2021"
}
@article{London2018,
archivePrefix = {arXiv},
arxivId = {gr-qc/1708.00404},
author = {London, Lionel and Khan, Sebastian and Fauchon-Jones, Edward and Garc'ia, Cecilio and Hannam, Mark and Husa, Sascha and Jim{'{e}}nez-Forteza, Xisco and Kalaghatgi, Chinmay and Ohme, Frank and Pannarale, Francesco},
doi = {10.1103/PhysRevLett.120.161102},
eprint = {1708.00404},
journal = {\prl},
keywords = {General Relativity and Quantum Cosmology},
month = {apr},
number = {16},
pages = {161102},
primaryClass = {gr-qc},
title = {{First Higher-Multipole Model of Gravitational Waves from Spinning and Coalescing Black-Hole Binaries}},
volume = {120},
year = {2018}
}
@article{Husa2016,
archivePrefix = {arXiv},
arxivId = {gr-qc/1508.07250},
author = {Husa, S and Khan, S and Hannam, M and P{"{u}}rrer, M and Ohme, F and Forteza, X.$\sim$J. and Boh{'{e}}, A},
doi = {10.1103/PhysRevD.93.044006},
eprint = {1508.07250},
journal = {\prd},
month = {feb},
number = {4},
pages = {44006},
primaryClass = {gr-qc},
title = {{Frequency-domain gravitational waves from nonprecessing black-hole binaries. I. New numerical waveforms and anatomy of the signal}},
volume = {93},
year = {2016}
}
@article{Khan2016,
archivePrefix = {arXiv},
arxivId = {gr-qc/1508.07253},
author = {Khan, S and Husa, S and Hannam, M and Ohme, F and P{"{u}}rrer, M and Forteza, X.$\sim$J. and Boh{'{e}}, A},
doi = {10.1103/PhysRevD.93.044007},
eprint = {1508.07253},
journal = {\prd},
month = {feb},
number = {4},
pages = {44007},
primaryClass = {gr-qc},
title = {{Frequency-domain gravitational waves from nonprecessing black-hole binaries. II. A phenomenological model for the advanced detector era}},
volume = {93},
year = {2016}
}