Skip to content
Snippets Groups Projects
Commit 363535ab authored by Gregory Ashton's avatar Gregory Ashton
Browse files

Merge branch '545-tutorial-has-incorrect-documentation' into 'master' 

Resolve "Tutorial has incorrect documentation"

Closes #545

See merge request !994
parents a7c0009b b066db67
No related branches found
No related tags found
1 merge request!994Resolve "Tutorial has incorrect documentation"
Pipeline #256981 passed with warnings
Stage: initial
Stage: test
Stage: docs
Stage: deploy
Hide whitespace changes
Inline Side-by-side
examples/tutorials/visualising_the_results.ipynb
+ 22
− 7
View file @ 363535ab
%% Cell type:code id: tags:
```
``` python
! rm visualising_the_results/*
```
%% Cell type:markdown id: tags:
# Visualising the results
In this tutorial, we demonstrate the plotting tools built-in to `bilby` and how to extend them. First, we run a simple injection study and return the `result` object.
%% Cell type:code id: tags:
```
``` python
import bilby
import matplotlib.pyplot as plt
%matplotlib inline
time_duration = 4. # time duration (seconds)
sampling_frequency = 2048. # sampling frequency (Hz)
outdir = 'visualising_the_results' # directory in which to store output
label = 'example' # identifier to apply to output files
# specify injection parameters
injection_parameters = dict(
mass_1=36., # source frame (non-redshifted) primary mass (solar masses)
mass_2=29., # source frame (non-redshifted) secondary mass (solar masses)
mass_1=36., # detector frame (redshifted) primary mass (solar masses)
mass_2=29., # detector frame (redshifted) secondary mass (solar masses)
a_1=0.4, # primary dimensionless spin magnitude
a_2=0.3, # secondary dimensionless spin magnitude
tilt_1=0.5, # polar angle between primary spin and the orbital angular momentum (radians)
tilt_2=1.0, # polar angle between secondary spin and the orbital angular momentum
phi_12=1.7, # azimuthal angle between primary and secondary spin (radians)
phi_jl=0.3, # azimuthal angle between total angular momentum and orbital angular momentum (radians)
luminosity_distance=200., # luminosity distance to source (Mpc)
theta_jn=0.4, # inclination angle between line of sight and orbital angular momentum (radians)
phase=1.3, # phase (radians)
ra=1.375, # source right ascension (radians)
dec=-1.2108, # source declination (radians)
geocent_time=1126259642.413, # reference time at geocentre (time of coalescence or peak amplitude) (GPS seconds)
psi=2.659 # gravitational wave polarisation angle
)
# specify waveform arguments
waveform_arguments = dict(
waveform_approximant='IMRPhenomPv2', # waveform approximant name
reference_frequency=50., # gravitational waveform reference frequency (Hz)
)
# set up the waveform generator
waveform_generator = bilby.gw.waveform_generator.WaveformGenerator(
sampling_frequency=sampling_frequency, duration=time_duration,
frequency_domain_source_model=bilby.gw.source.lal_binary_black_hole,
parameters=injection_parameters, waveform_arguments=waveform_arguments)
# create the frequency domain signal
hf_signal = waveform_generator.frequency_domain_strain()
# initialise an interferometer based on LIGO Hanford, complete with simulated noise and injected signal
IFOs = [bilby.gw.detector.get_interferometer_with_fake_noise_and_injection(
'H1', injection_polarizations=hf_signal, injection_parameters=injection_parameters, duration=time_duration,
sampling_frequency=sampling_frequency, outdir=outdir)]
# first, set up all priors to be equal to a delta function at their designated value
priors = bilby.gw.prior.BBHPriorDict(injection_parameters.copy())
# then, reset the priors on the masses and luminosity distance to conduct a search over these parameters
priors['mass_1'] = bilby.core.prior.Uniform(20, 50, 'mass_1')
priors['mass_2'] = bilby.core.prior.Uniform(20, 50, 'mass_2')
priors['luminosity_distance'] = bilby.core.prior.Uniform(100, 300, 'luminosity_distance')
# compute the likelihoods
likelihood = bilby.gw.likelihood.GravitationalWaveTransient(interferometers=IFOs, waveform_generator=waveform_generator)
result = bilby.core.sampler.run_sampler(likelihood=likelihood, priors=priors, sampler='dynesty', npoints=100,
injection_parameters=injection_parameters, outdir=outdir, label=label,
walks=5)
# display the corner plot
plt.show()
```
%% Cell type:markdown id: tags:
In running this code, we already made the first plot! In the function `bilby.detector.get_interferometer_with_fake_noise_and_injection`, the ASD, detector data, and signal are plotted together. This figure is saved under `visualsing_the_results/H1_frequency_domain_data.png`. Note that `visualising_the_result` is our `outdir` where all the output of the run is stored. Let's take a quick look at that directory now:
%% Cell type:code id: tags:
```
``` python
!ls visualising_the_results/
```
%% Cell type:markdown id: tags:
## Corner plots
Now lets make some corner plots. You can easily generate a corner plot using `result.plot_corner()` like this:
%% Cell type:code id: tags:
```
``` python
result.plot_corner()
plt.show()
```
%% Cell type:markdown id: tags:
In a notebook, this figure will display. But by default the file is also saved to `visualising_the_result/example_corner.png`. If you change the label to something more descriptive then the `example` here will of course be replaced.
%% Cell type:markdown id: tags:
You may also want to plot a subset of the parameters, or perhaps add the `injection_paramters` as lines to check if you recovered them correctly. All this can be done through `plot_corner`. Under the hood, `plot_corner` uses
[chain consumer](https://samreay.github.io/ChainConsumer/index.html), and all the keyword arguments passed to `plot_corner` are passed through to [the `plot` function of chain consumer](https://samreay.github.io/ChainConsumer/chain_api.html#chainconsumer.plotter.Plotter.plot).
### Adding injection parameters to the plot
In the previous plot, you'll notice `bilby` added the injection parameters to the plot by default. You can switch this off by setting `truth=None` when you call `plot_corner`. Or to add different injection parameters to the plot, just pass this as a keyword argument for `truth`. In this example, we just add a line for the luminosity distance by passing a dictionary of the value we want to display.
%% Cell type:code id: tags:
```
``` python
result.plot_corner(truth=dict(luminosity_distance=201))
plt.show()
```
%% Cell type:markdown id: tags:
### Plot a subset of the corner plot
Or, to plot just a subset of parameters, just pass a list of the names you want.
%% Cell type:code id: tags:
```
``` python
result.plot_corner(parameters=['mass_1', 'mass_2'], filename='{}/subset.png'.format(outdir))
plt.show()
```
%% Cell type:markdown id: tags:
Notice here, we also passed in a keyword argument `filename=`, this overwrites the default filename and instead saves the file as `visualising_the_results/subset.png`. Useful if you want to create lots of different plots. Let's check what the outdir looks like now
%% Cell type:code id: tags:
```
``` python
!ls visualising_the_results/
```
%% Cell type:markdown id: tags:
## Alternative
If you would prefer to do the plotting yourself, you can get hold of the samples and the ordering as follows and then plot with a different module. Here is an example using the [`corner`](http://corner.readthedocs.io/en/latest/) package
%% Cell type:code id: tags:
```
``` python
import corner
samples = result.samples
labels = result.parameter_labels
fig = corner.corner(samples, labels=labels)
plt.show()
```
%% Cell type:markdown id: tags:
## Other plots
We also include some other types of plots which may be useful. Again, these are built on chain consumer so you may find it useful to check the [documentation](https://samreay.github.io/ChainConsumer/chain_api.html#plotter-class) to see how these plots can be extended. Below, we show just one example of these.
#### Distribution plots
These plots just show the 1D histograms for each parameter
%% Cell type:code id: tags:
```
``` python
result.plot_marginals()
plt.show()
```
%% Cell type:markdown id: tags:
#### Best-Fit Time Domain Waveform plot
Some plots specific to compact binary coalescence parameter estimation results can
be created by re-loading the result as a `CBCResult`:
%% Cell type:code id: tags:
```
``` python
from bilby.gw.result import CBCResult
cbc_result = CBCResult.from_json("visualising_the_results/example_result.json")
cbc_result.plot_waveform_posterior()
plt.show()
```
%% Cell type:markdown id: tags:
Again, notice that the plot is saved as a "waveform.png" in the output dir.
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment