GBGPU Tutorial

GBGPU is a GPU-accelerated version of the FastGB waveform which has been developed by Neil Cornish, Tyson Littenberg, Travis Robson, and Stas Babak. It computes gravitational waveforms for Galactic binary systems observable by LISA using a fast/slow-type decomposition. For more details on the original construction of FastGB see arXiv:0704.1808.

The current version of the code is very closely related to the implementation of FastGB in the LISA Data Challenges’ Python code package. The waveform code is entirely Python-based. It is about 1/2 the speed of the full C version, but much simpler in Python for right now. There are also many additional functions including fast likelihood computations for individual Galactic binaries, as well as fast C-based methods to combine waveforms into global fitting templates.

The code is CPU/GPU agnostic. CUDA and NVIDIA GPUs are required to run these codes for GPUs.

See the documentation for more details. This code was designed for arXiv:2205.03461. If you use any part of this code, please cite arXiv:2205.03461, its Zenodo page, arXiv:0704.1808, and arXiv:1806.00500.

import numpy as np
import time

import matplotlib.pyplot as plt
%matplotlib inline

from gbgpu.gbgpu import GBGPU
from gbgpu.thirdbody import GBGPUThirdBody

from gbgpu.utils.constants import *
from gbgpu.utils.utility import *

Generating Galactic binary waveforms

Initialize the waveform class.

gb = GBGPU(use_gpu=False)

Setup all the binary information. GBGPU operates in a vectorized manner, so it takes arrays of parameters as inputs.

dt = 10.0
Tobs = 1.0 * YEAR


# number of points in waveform
# if None, will determine inside the code based on amp, f0 (and P2 if running third-body waveform)
N = None

# number of binaries to batch
num_bin = 10

# parameters
amp = 2e-23  # amplitude
f0 = 2e-3  # f0
fdot = 7.538331e-18  # fdot
fddot = 0.0 # fddot
phi0 = 0.1  # initial phase
iota = 0.2  # inclination
psi = 0.3  # polarization angle
lam = 0.4  # ecliptic longitude
beta_sky = 0.5  # ecliptic latitude


# for batching
amp_in = np.full(num_bin, amp)
f0_in = np.full(num_bin, f0)
fdot_in = np.full(num_bin, fdot)
fddot_in = np.full(num_bin, fddot)
phi0_in = np.full(num_bin, phi0)
iota_in = np.full(num_bin, iota)
psi_in = np.full(num_bin, psi)
lam_in = np.full(num_bin, lam)
beta_sky_in = np.full(num_bin, beta_sky)

Generate the waveforms.

params = np.array(
    [amp_in, f0_in, fdot_in, fddot_in, phi0_in, iota_in, psi_in, lam_in, beta_sky_in,]
)

gb.run_wave(*params, N=N, dt=dt, T=Tobs, oversample=2)

# signal from first binary
A = gb.A[0]
freqs = gb.freqs[0]
print("signal length:", A.shape)
plt.plot(freqs, np.abs(A))
plt.ylabel("TDI-A Channel", fontsize=16)
plt.xlabel("freqs (Hz)", fontsize=16)
dx = 7e-7
plt.xlim(f0 - dx, f0 + dx)

signal length: (128,)

(0.0019993, 0.0020007000000000002)

Adding additional GB astrophysics

It is possible in GBGPU to inherit a special class (``gbgpu.gbgpu.InheritGBGPU`) <https://mikekatz04.github.io/GBGPU/html/user/derivedwaves.html#gbgpu.gbgpu.InheritGBGPU>`__ that allows users to add other types of astrophysics to the FastGB waveform. This requires that the astrophysical effects vary slowly.

The methods that need to be written when adding astrophysics are prepare_additional_args, special_get_N, shift_frequency, and add_to_argS.

prepare_additional_args: Prepares all arguments beyond the base GB parameters.

special_get_N: Implemented if the new setup puts limitations on the sampling rate in the time-domain for the slow part of the waveform.

shift_frequency: Shifts the frequency in the slow computation.

add_to_argS: Adjusts the phasing in the transfer function of the slow waveform.

See `gbgpu.thirdbody.GBGPUThirdBody <https://mikekatz04.github.io/GBGPU/html/user/derivedwaves.html#gbgpu.thirdbody.ThirdBody>`__ for an example and more information.

Example: Third-body in orbit around the inner binary

gb_third = GBGPUThirdBody(use_gpu=False)

A2 = 400.0  # third body amplitude parameter
varpi = 0.0  # varpi phase parameter
e2 = 0.3  # eccentricity of third body
P2 = 1.2  # period of third body
T2 = 0.5 * P2  # time of periapsis passage of third body


A2_in = np.full(num_bin, A2)
P2_in = np.full(num_bin, P2)
varpi_in = np.full(num_bin, varpi)
e2_in = np.full(num_bin, e2)
T2_in = np.full(num_bin, T2)

params = np.array(
    [amp_in, f0_in, fdot_in, fddot_in, phi0_in, iota_in, psi_in, lam_in, beta_sky_in, A2_in, varpi_in, e2_in, P2_in, T2_in]
)

gb_third.run_wave(*params, N=N, dt=dt, T=Tobs, oversample=2)

# signal from first binary
A_third = gb_third.A[0]
freqs = gb_third.freqs[0]
print("Third-body signal length:", A_third.shape)
plt.plot(freqs, np.abs(A), label="No third body")
plt.plot(freqs, np.abs(A_third), label="No third body")
plt.ylabel("TDI-A Channel", fontsize=16)
plt.xlabel("freqs (Hz)", fontsize=16)
dx = 7e-7
plt.xlim(f0 - dx, f0 + dx)

Third-body signal length: (128,)

(0.0019993, 0.0020007000000000002)

Calculating the Information Matrix

params = np.array(
    [amp_in, f0_in, fdot_in, fddot_in, phi0_in, iota_in, psi_in, lam_in, beta_sky_in,]
)

inds = np.array([0, 1, 2, 4, 5, 6, 7, 8])

info_matrix = gb.information_matrix(
    params,
    easy_central_difference=False,
    eps=1e-9,
    inds=inds,
    N=1024,
    dt=dt,
    T=Tobs,
)

cov = np.linalg.pinv(info_matrix)

cov.shape

(10, 8, 8)

Covariance Matrix for first binary:

cov[0]

array([[ 1.03670758e-02, -3.03751670e-09, -4.34853085e-06,
         9.30287319e-05, -4.16250330e-04,  5.58219344e-04,
        -5.55445555e-05,  4.37703212e-04],
       [-3.03751670e-09,  2.75957296e-12,  4.60004509e-09,
        -1.30903471e-07,  2.27487105e-10, -7.85411493e-07,
         1.13683606e-07, -6.31724576e-08],
       [-4.34853085e-06,  4.60004509e-09,  8.78628478e-06,
        -2.50761137e-04,  3.73944799e-07, -1.50454902e-03,
         1.72593472e-04, -7.75682309e-05],
       [ 9.30287318e-05, -1.30903471e-07, -2.50761137e-04,
         7.17103859e-03, -9.43217244e-06,  4.30257258e-02,
        -5.08042272e-03,  1.88810467e-03],
       [-4.16250330e-04,  2.27487092e-10,  3.73944776e-07,
        -9.43217177e-06,  1.67174837e-05, -5.65945132e-05,
         6.51025824e-06, -1.93069353e-05],
       [ 5.58219343e-04, -7.85411493e-07, -1.50454902e-03,
         4.30257258e-02, -5.65945172e-05,  2.58151320e-01,
        -3.04822260e-02,  1.13284197e-02],
       [-5.55445553e-05,  1.13683606e-07,  1.72593472e-04,
        -5.08042272e-03,  6.51025872e-06, -3.04822260e-02,
         1.08018147e-02, -8.85702255e-04],
       [ 4.37703212e-04, -6.31724576e-08, -7.75682309e-05,
         1.88810467e-03, -1.93069355e-05,  1.13284197e-02,
        -8.85702255e-04,  9.57210495e-03]])

Standard deviation on the marginalized parameters:

params[inds, 0] * cov[0].diagonal() ** (1/2)

array([2.03637676e-24, 3.32239249e-09, 2.23448704e-20, 8.46819850e-03,
       8.17740391e-04, 1.52425781e-01, 4.15727117e-02, 4.89185674e-02])

Plot the Information Matrix ellipse for the intial frequency and frequency derivative:

ind1 = 1
ind2 = 2
inds_get = np.array([
    [ind1, ind1],
    [ind1, ind2],
    [ind2, ind1],
    [ind2, ind2]
]).T
sub_mat = cov[0][tuple(inds_get)].reshape(2, 2)

# calculate and draw covariance ellipse
a, b, b, c = sub_mat.flatten()

lam1 = (a + c) / 2. + np.sqrt(((a-c)/2) ** 2 + b ** 2)
lam2 = (a + c) / 2. - np.sqrt(((a-c)/2) ** 2 + b ** 2)

if b == 0. and a >= c:
    theta = 0.0
elif b == 0. and a < c:
    theta = np.pi / 2.
else:
    theta = -np.arctan2(lam1 - a, b)

t_vals = np.linspace(0., 2 * np.pi, 1000)
x = np.sqrt(lam1) * np.cos(theta) * np.cos(t_vals) - np.sqrt(lam2) * np.sin(theta) * np.sin(t_vals)
y = np.sqrt(lam1) * np.sin(theta) * np.cos(t_vals) + np.sqrt(lam2) * np.cos(theta) * np.sin(t_vals)

x_in = params[1, 0] * (1 + x)
y_in = params[2, 0] * (1 + y)
plt.plot(x_in, y_in)

[<matplotlib.lines.Line2D at 0x154fa1220>]

Utility functions

GBGPU provides many utility functions. Below are some examples of GBGPU utility functions. This may not include all utility functions. See the utility documentation for all included utility functions.

Get the instantaneous gravitational wave frequency

Given \(f_0\), \(\dot{f}_0\), and \(\ddot{f}_0\) calculate the instantaneous frequency of the gravitational waves approximating the frequency evolution to quadratic order:

num = 15
f0 = np.random.uniform(0.001, 0.002, num)
fdot = np.random.uniform(1e-17, 2e-17, num)
fddot = np.random.uniform(1e-28, 2e-28, num)

t = np.arange(0.0, 10 * YEAR, YEAR/200)

# cast t to be for all binaries
t = t[None, :]

f = get_fGW(f0, fdot, fddot, t)

plt.plot(t[0], f[0])
plt.xlabel(r"$t$ (s)")
plt.ylabel(r"$f_{gw}$ (Hz)")

Text(0, 0.5, '$f_{gw}$ (Hz)')

Get amplitude (for slowly evolving source)

m1 = np.full(num, 0.2)
m2 = np.full(num, 0.1)
f0 = np.full(num, 0.1)
d = np.full(num, 6.0)  # kpc
amp = get_amplitude(m1, m2, f0, d)

print(amp)

[6.37190427e-23 6.37190427e-23 6.37190427e-23 6.37190427e-23
 6.37190427e-23 6.37190427e-23 6.37190427e-23 6.37190427e-23
 6.37190427e-23 6.37190427e-23 6.37190427e-23 6.37190427e-23
 6.37190427e-23 6.37190427e-23 6.37190427e-23]

Get \(\dot{f}\)

m1 = np.full(num, 0.2)
m2 = np.full(num, 0.1)
f0 = np.full(num, 0.001)

fdot = get_fdot(f0, m1=m1, m2=m2)

print(fdot)

[1.73056526e-19 1.73056526e-19 1.73056526e-19 1.73056526e-19
 1.73056526e-19 1.73056526e-19 1.73056526e-19 1.73056526e-19
 1.73056526e-19 1.73056526e-19 1.73056526e-19 1.73056526e-19
 1.73056526e-19 1.73056526e-19 1.73056526e-19]

Determine necessary sampling rate in the time-domin

get_N(amp, f0, Tobs, oversample=1)

array([64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64])

Citations

You can access the necessary citations for the specific waveforms with the citation property.

print(gb.citation)

@software{michael_l_katz_2022_6500434,
  author       = {Michael L. Katz},
  title        = {mikekatz04/GBGPU: First official public release!},
  month        = apr,
  year         = 2022,
  publisher    = {Zenodo},
  version      = {v1.0.0},
  doi          = {10.5281/zenodo.6500434},
  url          = {https://doi.org/10.5281/zenodo.6500434}
}

@article{Cornish:2007if,
    author = "Cornish, Neil J. and Littenberg, Tyson B.",
    title = "{Tests of Bayesian Model Selection Techniques for Gravitational Wave Astronomy}",
    eprint = "0704.1808",
    archivePrefix = "arXiv",
    primaryClass = "gr-qc",
    doi = "10.1103/PhysRevD.76.083006",
    journal = "Phys. Rev. D",
    volume = "76",
    pages = "083006",
    year = "2007"
}

@article{Robson:2018svj,
    author = "Robson, Travis and Cornish, Neil J. and Tamanini, Nicola and Toonen, Silvia",
    title = "{Detecting hierarchical stellar systems with LISA}",
    eprint = "1806.00500",
    archivePrefix = "arXiv",
    primaryClass = "gr-qc",
    doi = "10.1103/PhysRevD.98.064012",
    journal = "Phys. Rev. D",
    volume = "98",
    number = "6",
    pages = "064012",
    year = "2018"
}