Working With Bilby Components¶
Written by Rhiannon Udall, with editing from Jacob Golomb
This tutorial is oriented around the skills required to work with some of the most important classes in bilby, before any analysis has been performed. These skills are often most important in post-processing, and we'll be building on these foundations in future examples. However, they also have utility for people who are learning to extend and modify bilby for their own purposes, since there are a lot of gotchas and unintuitive tricks to be aware of.
It assumes that you have the basic knowledge of PE and bilby as introduced in the Open Data Workshop (Section 8). If you've never used bilby before, start there!
It also includes a few challenge problems, implementing features which are present in bilby and ptfp. Obviously, if you want the feature yourself go ahead and use them! For those who want to learn how to modify bilby, though, these are useful excercises to complete.
0. Standard Setup¶
First, some imports that will almost universally be needed
import copy
import numpy as np
import matplotlib.pyplot as plt
import bilby
from typing import Tuple
%matplotlib inline
Although ptfp is not a focus of this tutorial, we'll also import it, and use it to set some pleasant defaults (credit to Derek Davis, Sophie Hourihane, and Simona Miller for these)
from ptfp.configuration import set_better_matplotlib_defaults
set_better_matplotlib_defaults()
/home/rhiannon.udall/.conda/envs/bilby-MRs-testing/lib/python3.11/site-packages/gwpy/time/__init__.py:36: UserWarning: Wswiglal-redir-stdio:
SWIGLAL standard output/error redirection is enabled in IPython.
This may lead to performance penalties. To disable locally, use:
with lal.no_swig_redirect_standard_output_error():
...
To disable globally, use:
lal.swig_redirect_standard_output_error(False)
Note however that this will likely lead to error messages from
LAL functions being either misdirected or lost when called from
Jupyter notebooks.
To suppress this warning, use:
import warnings
warnings.filterwarnings("ignore", "Wswiglal-redir-stdio")
import lal
from lal import LIGOTimeGPS
1. Priors and PriorDicts¶
We'll begin with the simplest and (subective statement here) single best component of bilby: the prior class. bilby makes priors very easy to work with. You were briefly introduced to these in ODW-3.1, but let's take a closer look. The generic parent class of all priors in bilby is bilby.core.prior.Prior, which you can view in a web browser here. This is not usable on its own - it's just an abstract base class - but it allows us to describe how priors actually get used in bilby.
The Prior¶
One of the most common things you'll do with a prior is draw samples. This is very simple in bilby, since for some prior it's just prior.sample(). The methods for doing this efficiently are clever, and they relate to how new priors can be created in practice, but we'll save that discussion for Appendix A. For now it suffices to say that for any prior, you can sample from it arbitrarily. Let's see this in practice with some example priors:
a_prior = bilby.core.prior.Uniform(minimum=-5, maximum=5, name="a", latex_label="$a$")
b_prior = bilby.core.prior.PowerLaw(alpha=2, minimum=1, maximum=100, name='b', latex_label='$b$', boundary='periodic')
c_prior = bilby.core.prior.Gaussian(mu=0, sigma=5, name="c", latex_label="$c$")
The first of these is a Uniform Prior, which is exactly what it sounds like:
a_samples = a_prior.sample(50000)
a_domain = np.linspace(-5, 5, 100)
fig, ax = plt.subplots()
ax.hist(a_samples, bins=50, histtype='step', density=True, label="Probability from 50000 Samples")
ax.plot(a_domain, a_prior.prob(a_domain), label="Analytic Probability")
ax.set_xlabel(a_prior.latex_label)
ax.set_ylabel("$p(a)$")
ax.legend()
fig
As you can see, we called sample to produce many samples distributed according to prob, and then we can verify by calling prob directly. Next is a power law prior, in which
$$ p(x) \propto x^{\alpha} $$
That ends up looking like this for $\alpha=2$:
b_samples = b_prior.sample(50000)
b_domain = np.linspace(1, 100, 100)
fig, ax = plt.subplots()
ax.hist(b_samples, bins=50, histtype='step', density=True, label="Probability from 50000 Samples")
ax.plot(b_domain, b_prior.prob(b_domain), label="Analytic Probability")
ax.set_xlabel(b_prior.latex_label)
ax.set_ylabel("$p(b)$")
ax.legend()
fig
Finally we have the classic Gaussian prior:
c_samples = c_prior.sample(50000)
c_domain = np.linspace(-25, 25, 100)
fig, ax = plt.subplots()
ax.hist(c_samples, bins=50, histtype='step', density=True, label="Probability from 50000 Samples")
ax.plot(c_domain, c_prior.prob(c_domain), label="Analytic Probability")
ax.set_xlabel(c_prior.latex_label)
ax.set_ylabel("$p(c)$")
ax.legend(loc='upper left')
fig
As you run this, you may notice that this is all going very fast. That's because actual sampling is from a unit line (very fast), and then the transformations are just (vectorized) analytic math in this case.
Now, let's put these all together into something more complicated: the prior dict.
PriorDicts¶
PriorDicts are exactly what they sound like: dictionaries of priors. They're easily instantiated as such:
prior_dict = bilby.core.prior.PriorDict(
{'a':a_prior,
'b':b_prior,
'c':c_prior}
)
Sampling from a PriorDict only takes one call, rather than the three we needed before:
prior_dict_samples = prior_dict.sample(5000)
prior_dict_samples
{'a': array([-0.11063755, 4.84514728, 1.70398471, ..., 4.90045904,
-3.29989104, -1.34539395]),
'b': array([82.3856183 , 27.06464433, 66.46781613, ..., 87.45323611,
96.8585779 , 94.53205941]),
'c': array([ 2.75739375, -0.62899639, -5.90322004, ..., -9.03435585,
-2.57534025, 0.66528773])}
You can see this now returns a dictionary of the samples. The useful thing is that the ith index of each array corresponds to the others with that same index. Since our priors are independent, this doesn't make a difference, but it's very important when working with joint priors that are not independent: a slice of the ith element in each array forms a correctly distributed sample.
We can do some interesting things like looking at a 2-dimensional slice of the data:
fig, ax = plt.subplots()
# For *reasons* the prior_dict.prob() is not well behaved, so we still use the individual priors .prob() methods
# We can get this by then computing the product of the two
prior_dict_samples['prob'] = prior_dict['b'].prob(prior_dict_samples['b']) * prior_dict['c'].prob(prior_dict_samples['c'])
# Scatter the two
# Passing 'c=prior_dict_samples['prob']' lets us set a color gradient based on the probability
scatter = ax.scatter(prior_dict_samples['b'], prior_dict_samples['c'], c=prior_dict_samples['prob'], label="Samples", alpha=0.5)
fig.colorbar(scatter, label="$p(b, c)$")
ax.set_xlabel("$b$")
ax.set_ylabel("$c$")
ax.legend()
fig
Challenge:¶
One common situation involving priors is setting priors on some quantities (say, binary masses) and then looking at plots of a derived quantity (say, chirp mass and mass ratio). These quantities are defined as:
$$ q = \frac{m_2}{m_1} $$
$$ \mathcal{M}_c = \frac{(m_1 m_2)^{3/5}}{(m_1 + m_2)^{1/5}} = \frac{m_1 q^{3/5}}{(1+q)^{1/5}} $$
Let's do exactly that: Make a prior which is uniform in mass_1 and mass_2, then draw these samples and convert to chirp_mass and mass_ratio. Make a scatter plot like above in the chirp_mass and mass_ratio slice.
Hint: It's possible to sample chirp_mass and mass_ratio under this prior (that is, uniform in components) by use of bilby.gw.prior.UniformInComponentsChirpMass and bilby.gw.prior.UniformInComponentsMassRatio. Try drawing some samples and plotting them, to make sure they line up with the answer you got. These use an analytic transformation and compute the appropriate Jacobian, which we'll discuss more in Appendix A.
# Solution to the challenge goes here!
2. The WaveformGenerator Class¶
Now, let's go to something gravitational wave specific: how do you produce a waveform with bilby? The answer is the WaveformGenerator, which you interacted with briefly in the ODW tutorial. The point of the WaveformGenerator is to handle the steps between a pure waveform model and having a waveform that you can actually use for analysis. Luckily they aren't too complicated, but it can be useful to see what's needed to actually make one from scratch:
class WaveformGenerator(object):
def __init__(self, duration=None, sampling_frequency=None, start_time=0, frequency_domain_source_model=None,
time_domain_source_model=None, parameters=None, parameter_conversion=None, waveform_arguments=None):
We can group the arguments into 4 different groups:
- Parameters that describe the data:
duration,sampling_frequency, andstart_time - The model (time or frequency domain) which will be used to do the computation. This will usually be
lal_binary_black_holefor any BBH model, andlal_binary_neutron_starfor any BNS model. - The
parametersof the binary itself: generally you shouldn't actually pass these, because there's a place to put them more generically down the line (as we shall see later). - The
parameter_conversionand thewaveform_arguments. The former tells generator how to convert some parameters (say chirp mass and mass ratio) into others (say the component masses) if necessary, and when it is passed as None (the default) it will use standard CBC conversions. The latter tells the generator what special (non-parameter) arguments to put into the model. The most prominent of these are the minimum/maximum/reference frequencies and the waveform approximant, but others may apply.
In order to explore this further, let's make a reasonable example instance.
# Remember to always set configurations like these at the start and reference them as variables!
# You never know when you'll need them again later, and if you then change them you'll have to go through
# the whole notebook changing things, which is no fun at all and a recipe for weird bugs
duration = 4
sampling_frequency=4096
start_time=0
post_trigger_duration=2
trigger_time = start_time + duration - post_trigger_duration
minimum_frequency = 20
reference_frequency = 20
waveform_approximant = 'IMRPhenomXPHM'
waveform_generator = bilby.gw.waveform_generator.WaveformGenerator(
duration=duration,
sampling_frequency=sampling_frequency,
start_time=start_time,
frequency_domain_source_model=bilby.gw.source.lal_binary_black_hole,
waveform_arguments={
'minimum_frequency':minimum_frequency,
'maximum_frequency':sampling_frequency / 2, # Nyquist frequency is sampling frequency / 2
'reference_frequency':reference_frequency,
'waveform_approximant':waveform_approximant
}
)
15:03 bilby INFO : Waveform generator initiated with frequency_domain_source_model: bilby.gw.source.lal_binary_black_hole time_domain_source_model: None parameter_conversion: bilby.gw.conversion.convert_to_lal_binary_black_hole_parameters
The main thing to do with a waveform generator is, well, generate waveforms. So, let's do that. There are two options: time and frequency domain. Both of these will work no matter what the domain of the source model we put in was, because the waveform generator will take care of the fft/ifft under the hood to convert back and forth. That being said, bilby is primarily oriented around frequency domain waveforms - in fact, bilby_pipe doesn't even allow waveforms which are explicitly time domain. So, we're going to plot a TD waveform, but remember that all the actual analysis for this waveform is done in the frequency domain:
waveform_parameters = {
'mass_1':30, # Note these are *detector frame* masses
'mass_2':30,
'a_1':0,
'a_2':0,
'tilt_1':0,
'tilt_2':0,
'phi_12':0,
'phi_jl':0,
'luminosity_distance':400,
'theta_jn':np.pi/3,
'phase':np.pi
}
waveform = waveform_generator.time_domain_strain(waveform_parameters)
fig, ax = plt.subplots()
ax.plot(waveform_generator.time_array, waveform['plus'], label='$+$ polarization')
ax.plot(waveform_generator.time_array, waveform['cross'], label='$\\times$ polarization')
ax.set_xlabel("Time [s]")
ax.set_ylabel("h")
ax.legend()
fig
Now, this actually looks kind of weird! The reason is that bilby defines the end of the waveform at the peak strain (a proxy for the merger). All the information that come over then wraps around if you aren't careful. It's important to note that the analysis itself will take care of this (as we'll see shortly). But if you want to check the plus and cross polarizations directly, you'll need to deal with this.
This is a reason to generally not use this output directly, but if you must you can implement the bilby standard formatting by rolling the result so that the post-trigger duration is the standard bilby value of 2 seconds:
# Rolling it backwards by 2 seconds, so 2 * sampling_frequency elements of the array
waveform['plus'] = np.roll(waveform['plus'], -post_trigger_duration * waveform_generator.sampling_frequency)
waveform['cross'] = np.roll(waveform['cross'], -post_trigger_duration * waveform_generator.sampling_frequency)
fig, ax = plt.subplots()
ax.plot(waveform_generator.time_array, waveform['plus'], label='$+$ polarization')
ax.plot(waveform_generator.time_array, waveform['cross'], label='$\\times$ polarization')
ax.set_xlabel("Time [s]")
ax.set_ylabel("h")
ax.legend()
fig
This is still not great but it is better. Among other issues, you can see oscillations well after the merger, where the ringdown should already have died away. Unfortunately, that's not physics, it's just an issue with data processing: if you count the number of peaks carefully, you'll see there are 16 in 2 seconds, corresponding to 8 Hz: the minimum frequency we put into the waveform generator. This came out due to the process of using ifft on the frequency domain waveform. This is not an issue in analysis, because bilby computes all the likelihoods in the frequency domain, meaning no effects like this are accrued. But, it's something to remember when making plots. One solution to this could be to whiten the data, which we'll talk about in a later tutorial, as the PSD would be very high at these frequencies (in analysis it will actually be artificially high below the minimum analysis frequency), such that these completely disappear.
Challenge¶
Create a negative waveform generator, which takes the arguments of a waveform generator and produces the negative result. This is useful for plotting (as we'll explore shortly). It's also a useful example of exploiting inheritance. I've doing the somewhat tricky part in the initialization already: this class either takes normal arguments or it takes an existing waveform generator and uses its arguments. What you need to do now is modify the methods of the waveform generator to make the plus and cross strain negative.
This is implemented in ptfp as ptfp.waveform_generators.NegativeWaveformGenerator.
hint: the super() function will call the class that was inherited from, which in this case is the normal waveform generator.
class NegativeWaveformGenerator(bilby.gw.waveform_generator.WaveformGenerator):
def __init__(
self,
parent=None,
duration=None,
sampling_frequency=None,
start_time=0,
frequency_domain_source_model=None,
time_domain_source_model=None,
parameters=None,
parameter_conversion=None,
waveform_arguments=None,
):
if parent is not None:
self.__dict__ = parent.__dict__.copy()
else:
super(NegativeWaveformGenerator, self).__init__(
duration,
sampling_frequency,
start_time,
frequency_domain_source_model,
time_domain_source_model,
parameters,
parameter_conversion,
waveform_arguments,
)
# Modify this function to get negative strain!
def frequency_domain_strain(self, parameters=None):
pass
3. The Interferometer Class¶
Now we'll discuss the use of the Interferometer in bilby. This class is what bilby uses to interact with two things:
- The geometry of the detectors (so, how gravitational waves project onto it)
- The data itself (if it is written)
We'll explore both of these functionalities a bit.
Using the detector geometry¶
One thing that the detector does in bilby is handle the projection of the waveform, which makes sense since that is a detector dependent process. First, let's see how to get the antenna response functions $F_{+/\times}$. You may recall these show up in the projection as:
$$ h = F_+ h_+ + F_\times h_\times $$
First we'll need an interferometer. These are actually annoyingly difficult to instantiate, requiring a multistage process. For this though we only need to start with an empty interferometer (that is, not adding any data yet). We can do this with the helper function bilby.gw.detector.get_empty_interferometer(). We just pass the name of the detector ("H1", "L1", etc.) and it will fetch the correct geometry automatically.
lho_detector = bilby.gw.detector.get_empty_interferometer("H1")
lho_detector.sampling_frequency = sampling_frequency
lho_detector.duration = duration
lho_detector.start_time = start_time
# copy.deepcopy creates a full copy of the object you pass it. This is valuable if you then want to modify
# a version of the object, while keeping the original in place
# we'll be using this PSD later, so I'm stashing it away now to prevent cell rerun shenanigans
default_psd = copy.deepcopy(lho_detector.power_spectral_density)
Now we can compute the antenna response to a specific part of the sky at a specific time. We'll do $\alpha=\delta=\psi=t_c=0$ for simplicity
# The standard bilby segment has 2 seconds post merger, so the merger time is the segment start time + the duration - 2
antenna_response_parameters = dict(ra=0, dec=0, psi=0, time=trigger_time)
# ** is a double splat operator: it takes the dictionary and expands it as a bunch of key=value arguments to the function
# The * is a single splat operator and works similarly on lists and tuples, expanding them out as a bunch of position arguments
plus_response = lho_detector.antenna_response(**antenna_response_parameters, mode="plus")
cross_response = lho_detector.antenna_response(**antenna_response_parameters, mode="cross")
print(f"The plus response at this time/sky location is :{plus_response}")
print(f"The cross response at this time/sky location is :{cross_response}")
# antenna response calls it time which get_detector_response calls it geocent_time
# while this makes a certain amount of sense, note that for the antenna response fractions of a second don't matter
# (only timescales such that the Earth has moved significantly do), so this is just something to keep track of
antenna_response_parameters['geocent_time'] = antenna_response_parameters.pop("time")
The plus response at this time/sky location is :-0.38131980492882117 The cross response at this time/sky location is :0.5945726842640582
It's pretty rare that you'll use these directly (only really if you are doing sensitivity studies or something over the entire sky, in which case bilby might not be the best option anyways). Much more common is to use the function which calls this and then projects the waveform: get_detector_response. Let's use our waveform again, but this time we need it in the frequency domain, and we project it onto the detector:
waveform_fd = waveform_generator.frequency_domain_strain(waveform_parameters)
lho_projected_waveform_fd = lho_detector.get_detector_response(
waveform_polarizations=waveform_fd,
parameters=antenna_response_parameters
)
fig, ax = plt.subplots()
# We plot the absolute value of the projected waveform, because it's a complex values frequency vector
ax.plot(lho_detector.frequency_array, np.abs(lho_projected_waveform_fd))
# When plotting FD waveforms, it's preferable to use a log-log plot
ax.set_xscale("log")
ax.set_yscale("log")
ax.set_xlabel("Frequency [Hz]")
ax.set_ylabel("Strain [$1/\sqrt{\\mathrm{Hz}}$]")
ax.set_xlim(20, 2048)
fig
ptfp implements these steps together as ptfp.waveforms.get_projected_frequency_domain_waveform.
We can also get the TD waveform now with ifft. Note that now we don't have to roll anything, because get_detector_response manages that part for us.
lho_projected_waveform_td = bilby.core.utils.infft(lho_projected_waveform_fd, sampling_frequency)
fig, ax = plt.subplots()
ax.plot(lho_detector.time_array, lho_projected_waveform_td, label="TD Waveform Projection")
ax.legend()
ax.set_xlabel("Time [s]")
ax.set_ylabel("Strain")
fig
Setting and Plotting the Data¶
There are multiple different kinds of "detector" data, depending on what you're trying to do. Obviously, there's real detector data, which can be read in a variety of ways. There's also fake Gaussian noise, in which one uses a model of the Gaussian component of the noise to produce idealized detector noise, and zero noise, in which one uses the same noise model but generates no noise at all (producing an unbiased but also unrealistic realization of the data).
The noise model in question is the PSD. You can look these up to understand where they come from, but for our purposes it's most important for how it shows up in the likelihood:
$$ \ln \mathcal{L} = -\frac{1}{2} \sum_k \left( \frac{|\tilde{d}_k - \tilde{h}_k|^2}{\sigma_k^2} + \ln(2 \pi \sigma_k^2)\right) $$
For this demonstration we'll work with the middle of those two: generating Gaussian noise. Now, our detector actually has a power spectral density set by default, but it's not a particularly realistic one. To get a more realistic one, we'll need to load it in from somewhere. Here I am using the GW150914 PSD (which will also be used in later examples).
lho_loaded_psd = bilby.gw.detector.PowerSpectralDensity.from_power_spectral_density_file("/home/rhiannon.udall/PlottingToolsForPEExampleRun/GW150914_LHO_PSD.dat")
lho_detector.power_spectral_density = lho_loaded_psd
fig, ax = plt.subplots()
ax.plot(default_psd.frequency_array, default_psd.psd_array, label="Original PSD")
ax.plot(lho_detector.frequency_array, lho_detector.power_spectral_density_array, label="Loaded PSD Written to Detector")
ax.legend()
ax.set_xscale("log")
ax.set_yscale("log")
ax.set_xlabel("Frequency [Hz]")
ax.set_ylabel("$S_n$ [$1/\sqrt{\mathrm{Hz}}$]")
ax.set_xlim(20, 1024)
ax.set_ylim(1e-50, 1e-38)
fig
You can see that the real PSD is less sensitive than the idealized PSD, because it's from O1. If we wanted to use an O4 PSD we'd find that the higher frequency sensitivity is now actually better than the idealized prediction, but the low frequency noise is still much worse. Now, let's generate some Gaussian noise from this PSD. This can be done with:
# By setting this seed we guarantee the data will show the same noise realization each time
# This is good for reproducibility, but obviously in real analyses you want to allow more randomness
np.random.seed(88170235)
lho_detector.set_strain_data_from_power_spectral_density(
sampling_frequency=sampling_frequency,
duration=duration,
start_time=start_time
)
fig, ax = plt.subplots()
# asd_from_freq_series() normalizes to the proper units of an ASD, which is 4 * \Delta f = 4 / duration
ax.plot(
lho_detector.frequency_array,
bilby.gw.utils.asd_from_freq_series(
lho_detector.strain_data.frequency_domain_strain, 1 / lho_detector.duration
),
label="Noise Realization")
ax.plot(lho_detector.frequency_array, lho_detector.amplitude_spectral_density_array, label="ASD")
ax.set_xscale("log")
ax.set_yscale("log")
ax.set_xlim(20, 1024)
ax.set_ylim(1e-24, 1e-18)
ax.set_xlabel("Frequency [Hz]")
ax.set_ylabel("Strain [$1/\sqrt{\mathrm{Hz}}$]")
ax.legend()
fig
fig, ax = plt.subplots()
ax.plot(lho_detector.time_array, lho_detector.strain_data.time_domain_strain, label="TD Strain")
ax.legend()
ax.set_xlabel("Time [Hz]")
ax.set_ylabel("Strain")
fig
Alright, so this looks good in the frequency domain, but still looks quite ugly in the time domain. We can make this better by looking at the whitened data instead:
fig, ax = plt.subplots()
ax.plot(
lho_detector.frequency_array,
bilby.gw.utils.asd_from_freq_series(
lho_detector.whitened_frequency_domain_strain, 1 / lho_detector.duration
),
label="Whitened Noise Realization")
ax.set_xscale("log")
ax.set_yscale("log")
ax.set_xlim(20, 1024)
ax.set_xlabel("Frequency [Hz]")
ax.set_ylabel("Whitened Strain [$1/\sqrt{\mathrm{Hz}}$]")
ax.legend()
fig
And generally more interesting, the whitened time domain strain:
whitened_td_strain = lho_detector.whitened_time_domain_strain
fig, ax = plt.subplots()
ax.plot(lho_detector.time_array, whitened_td_strain, label="Whitened TD Strain")
ax.legend()
ax.set_xlabel("Time [s]")
ax.set_ylabel("Whitened Strain")
fig
This looks a lot better! No weird low frequency stuff, and a normalized scale instead of $10^{-21}$.
Injecting a Signal¶
Injecting a signal with bilby is actually pretty easy with the inject_signal method:
all_parameters = copy.deepcopy(waveform_parameters)
all_parameters.update(antenna_response_parameters)
lho_detector.inject_signal(parameters=all_parameters, waveform_generator=waveform_generator)
lho_whitened_fd_waveform = lho_detector.whiten_frequency_series(lho_projected_waveform_fd)
lho_whitened_td_waveform = lho_detector.get_whitened_time_series_from_whitened_frequency_series(lho_whitened_fd_waveform)
fig, ax = plt.subplots()
ax.plot(lho_detector.time_array, whitened_td_strain, label="Whitened TD Strain")
ax.plot(lho_detector.time_array, lho_whitened_td_waveform, label="Whitened Injected Waveform")
ax.legend()
ax.set_xlabel("Time [s]")
ax.set_ylabel("Strain")
fig
15:03 bilby INFO : Injected signal in H1: 15:03 bilby INFO : optimal SNR = 12.56 15:03 bilby INFO : matched filter SNR = 13.03-1.60j 15:03 bilby INFO : mass_1 = 30 15:03 bilby INFO : mass_2 = 30 15:03 bilby INFO : a_1 = 0 15:03 bilby INFO : a_2 = 0 15:03 bilby INFO : tilt_1 = 0 15:03 bilby INFO : tilt_2 = 0 15:03 bilby INFO : phi_12 = 0 15:03 bilby INFO : phi_jl = 0 15:03 bilby INFO : luminosity_distance = 400 15:03 bilby INFO : theta_jn = 1.0471975511965976 15:03 bilby INFO : phase = 3.141592653589793 15:03 bilby INFO : ra = 0 15:03 bilby INFO : dec = 0 15:03 bilby INFO : psi = 0 15:03 bilby INFO : geocent_time = 2
And now we can see our waveform looks much better as well. This is why it is always best too look at whitened time domain data, because it's much more informative than the raw data would be.
Challenge¶
Now, a few things:
- First, make a spectrogram from the strain data we have (hint: you can get a
gwpytimeseries from the bilby strain data withinterferometer.strain_data.to_gwpy_timeseries()) - Next, use the
NegativeWaveformGeneratorto subtract the signal from the strain data, producing the residual. Since we know what we injected, the residual should be perfect, but as we shall see in tutorial 3 residuals won't necessarily be perfect in real analyses. (hint: if you do this correctly then thebilby INFOthat gets printed out should show an optimal SNR like the one above, but a matched filter SNR of around 0) - Once you have the residual data, make the spectrogram and timeseries plots to confirm that the signal has been removed.
# Make a copy to mess around with the residual, so you don't have to go back and remake the original interferometer all the time
residual_lho_detector = copy.deepcopy(lho_detector)
# Challenge solution goes here
4. The Likelihood Class¶
Finally, we come to our last subject (for now): the Likelihood. The Likelihood is one of the central components of bilby, but interestingly enough it's mostly just a container for all the stuff we've already done. This is reflected in its instantiation:
# This is the high mass prior from bilby_pipe
# https://git.ligo.org/lscsoft/bilby_pipe/-/blob/master/bilby_pipe/data_files/high_mass.prior?ref_type=heads
from bilby.core.prior import Uniform, Sine, Cosine, Constraint
high_mass_prior = bilby.core.prior.PriorDict(
dict(
chirp_mass = Uniform(name='chirp_mass', minimum=25, maximum=175, unit='$M_{\odot}$'),
mass_ratio = Uniform(name='mass_ratio', minimum=0.125, maximum=1),
mass_1 = Constraint(name='mass_1', minimum=1.001398, maximum=200),
mass_2 = Constraint(name='mass_2', minimum=1.001398, maximum=200),
a_1 = Uniform(name='a_1', minimum=0, maximum=0.99),
a_2 = Uniform(name='a_2', minimum=0, maximum=0.99),
tilt_1 = Sine(name='tilt_1'),
tilt_2 = Sine(name='tilt_2'),
phi_12 = Uniform(name='phi_12', minimum=0, maximum=2 * np.pi, boundary='periodic'),
phi_jl = Uniform(name='phi_jl', minimum=0, maximum=2 * np.pi, boundary='periodic'),
luminosity_distance = bilby.gw.prior.UniformSourceFrame(name='luminosity_distance', minimum=1e2, maximum=7e3, unit='Mpc'),
dec = Cosine(name='dec'),
ra = Uniform(name='ra', minimum=0, maximum=2 * np.pi, boundary='periodic'),
theta_jn = Sine(name='theta_jn'),
psi = Uniform(name='psi', minimum=0, maximum=np.pi, boundary='periodic'),
phase = Uniform(name='phase', minimum=0, maximum=2 * np.pi, boundary='periodic'),
)
)
likelihood = bilby.gw.likelihood.GravitationalWaveTransient(
interferometers=[lho_detector],
waveform_generator=waveform_generator,
priors=high_mass_prior
)
Now, we can use the likelihood to compute lots of things all at once, for example:
likelihood.parameters = all_parameters
waveform_polarizations = likelihood.waveform_generator.frequency_domain_strain(all_parameters)
likelihood.calculate_snrs(waveform_polarizations=waveform_polarizations, interferometer=likelihood.interferometers[0])
GravitationalWaveTransient._CalculatedSNRs(d_inner_h=(163.66882428148253-20.03536514491011j), optimal_snr_squared=157.77014047964337, complex_matched_filter_snr=(13.030274252422862-1.5950887637380655j), d_inner_h_array=None, optimal_snr_squared_array=None)
The real usefulness though will be explored in tutorial 3: if you can construct the likelihood used in an analysis (which is harder than it should be, but we'll get to that), then you have what you need to to mess with all the outputs of that analysis.
Appendix A: How Priors are Written¶
Under construction for now
Earlier we discussed some of the practical ways that priors are used. Now, let's take a closer look at the code: we see that what it actually does is draw a sample from a uniform distribution between 0 and 1 (which numpy can do), and then rescale this by computing the inverse CDF. This method is called "Inverse Transform Sampling." This also explains how multidimensional priors work: they just sample over multidimensional cubes of this type.
What this means is that, for any given prior the methods that need to be written are:
rescale, which is just the inverse CDF.probwhich is the probability for a given input value.