Skip to content
Snippets Groups Projects
Commit 77aa17aa authored by Moritz Huebner's avatar Moritz Huebner
Browse files

Merge branch 'add_wavform_plot_to_tutorial' into 'master' 

add waveform plot to tutorial

See merge request !817
parents 47dc99b3 678f1dbc
No related branches found
No related tags found
No related merge requests found
Pipeline #162400 passed
Stage: test
Stage: deploy
Hide whitespace changes
Inline Side-by-side
examples/tutorials/visualising_the_results.ipynb
+ 48
− 145
View file @ 77aa17aa
%% 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)
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()
```
%% Output
10:31 bilby INFO : No power spectral density provided, using aLIGO,zero detuning, high power.
10:31 bilby INFO : Injected signal in H1:
10:31 bilby INFO : optimal SNR = 120.28
10:31 bilby INFO : matched filter SNR = 120.60-0.46j
10:31 bilby INFO : mass_1 = 36.0
10:31 bilby INFO : mass_2 = 29.0
10:31 bilby INFO : a_1 = 0.4
10:31 bilby INFO : a_2 = 0.3
10:31 bilby INFO : tilt_1 = 0.5
10:31 bilby INFO : tilt_2 = 1.0
10:31 bilby INFO : phi_12 = 1.7
10:31 bilby INFO : phi_jl = 0.3
10:31 bilby INFO : luminosity_distance = 200.0
10:31 bilby INFO : theta_jn = 0.4
10:31 bilby INFO : phase = 1.3
10:31 bilby INFO : ra = 1.375
10:31 bilby INFO : dec = -1.2108
10:31 bilby INFO : geocent_time = 1126259642.413
10:31 bilby INFO : psi = 2.659
10:31 bilby WARNING : The waveform_generator start_time is not equal to that of the provided interferometers. Overwriting the waveform_generator.
10:31 bilby INFO : Running for label 'example', output will be saved to 'visualising_the_results'
10:31 bilby INFO : Using LAL version Branch: None;Tag: lalsuite-v6.60;Id: 413788d4afeff8b759d8d75abe589ae4a846120c;;Builder: Unknown User <>;Repository status: CLEAN: All modifications committed
10:31 bilby INFO : Search parameters:
10:31 bilby INFO : mass_1 = Uniform(minimum=20, maximum=50, name='mass_1', latex_label='$m_1$', unit=None, boundary=None)
10:31 bilby INFO : mass_2 = Uniform(minimum=20, maximum=50, name='mass_2', latex_label='$m_2$', unit=None, boundary=None)
10:31 bilby INFO : luminosity_distance = Uniform(minimum=100, maximum=300, name='luminosity_distance', latex_label='$d_L$', unit=None, boundary=None)
10:31 bilby INFO : a_1 = 0.4
10:31 bilby INFO : a_2 = 0.3
10:31 bilby INFO : tilt_1 = 0.5
10:31 bilby INFO : tilt_2 = 1.0
10:31 bilby INFO : phi_12 = 1.7
10:31 bilby INFO : phi_jl = 0.3
10:31 bilby INFO : theta_jn = 0.4
10:31 bilby INFO : phase = 1.3
10:31 bilby INFO : ra = 1.375
10:31 bilby INFO : dec = -1.2108
10:31 bilby INFO : geocent_time = 1126259642.413
10:31 bilby INFO : psi = 2.659
10:31 bilby INFO : Single likelihood evaluation took 2.305e-03 s
10:31 bilby INFO : Using sampler Dynesty with kwargs {'bound': 'multi', 'sample': 'rwalk', 'verbose': True, 'periodic': None, 'reflective': None, 'check_point_delta_t': 600, 'nlive': 100, 'first_update': None, 'walks': 5, 'npdim': None, 'rstate': None, 'queue_size': None, 'pool': None, 'use_pool': None, 'live_points': None, 'logl_args': None, 'logl_kwargs': None, 'ptform_args': None, 'ptform_kwargs': None, 'enlarge': None, 'bootstrap': None, 'vol_dec': 0.5, 'vol_check': 2.0, 'facc': 0.5, 'slices': 5, 'update_interval': 60, 'print_func': <bound method Dynesty._print_func of <bilby.core.sampler.dynesty.Dynesty object at 0x1236747f0>>, 'dlogz': 0.1, 'maxiter': None, 'maxcall': None, 'logl_max': inf, 'add_live': True, 'print_progress': True, 'save_bounds': False, 'n_effective': None}
10:31 bilby INFO : Checkpoint every n_check_point = 300000
10:31 bilby INFO : Using dynesty version 0.9.7
1864| logz ratio=7255.416 +/- 0.541 | dlogz: 0.118 > 0.100000
10:32 bilby WARNING : Run terminated with signal 2
10:32 bilby INFO : Writing checkpoint file visualising_the_results/example_resume.pickle
Exception while calling loglikelihood function:
params: [ 35.99600559 28.99205241 200.17592464]
args: []
kwargs: {}
exception:
Traceback (most recent call last):
File "/Users/virginiademilio/Documents/Bilby/bilby/bilby/core/sampler/dynesty.py", line 279, in _run_nested_wrapper
self.sampler.run_nested(**kwargs)
TypeError: run_nested() got an unexpected keyword argument 'n_effective'
During handling of the above exception, another exception occurred:
Traceback (most recent call last):
File "/Users/virginiademilio/anaconda3/envs/bilby_source/lib/python3.7/site-packages/dynesty/dynesty.py", line 805, in __call__
return self.func(x, *self.args, **self.kwargs)
File "/Users/virginiademilio/Documents/Bilby/bilby/bilby/core/sampler/base_sampler.py", line 583, in log_likelihood
return Sampler.log_likelihood(self, theta)
File "/Users/virginiademilio/Documents/Bilby/bilby/bilby/core/sampler/base_sampler.py", line 381, in log_likelihood
return self.likelihood.log_likelihood_ratio()
File "/Users/virginiademilio/Documents/Bilby/bilby/bilby/gw/likelihood.py", line 244, in log_likelihood_ratio
self.waveform_generator.frequency_domain_strain(self.parameters)
File "/Users/virginiademilio/Documents/Bilby/bilby/bilby/gw/waveform_generator.py", line 119, in frequency_domain_strain
transformed_model_data_points=self.time_array)
File "/Users/virginiademilio/Documents/Bilby/bilby/bilby/gw/waveform_generator.py", line 158, in _calculate_strain
model_strain = self._strain_from_model(model_data_points, model)
File "/Users/virginiademilio/Documents/Bilby/bilby/bilby/gw/waveform_generator.py", line 170, in _strain_from_model
return model(model_data_points, **self.parameters)
File "/Users/virginiademilio/Documents/Bilby/bilby/bilby/gw/source.py", line 71, in lal_binary_black_hole
phi_jl=phi_jl, **waveform_kwargs)
File "/Users/virginiademilio/Documents/Bilby/bilby/bilby/gw/source.py", line 272, in _base_lal_cbc_fd_waveform
waveform_dictionary, approximant)
File "/Users/virginiademilio/Documents/Bilby/bilby/bilby/gw/utils.py", line 798, in lalsim_SimInspiralChooseFDWaveform
waveform_dictionary, approximant)
File "/Users/virginiademilio/Documents/Bilby/bilby/bilby/core/sampler/dynesty.py", line 385, in write_current_state_and_exit
sys.exit(130)
SystemExit: 130
An exception has occurred, use %tb to see the full traceback.
SystemExit: 130
/Users/virginiademilio/anaconda3/envs/bilby_source/lib/python3.7/site-packages/IPython/core/interactiveshell.py:3334: UserWarning: To exit: use 'exit', 'quit', or Ctrl-D.
warn("To exit: use 'exit', 'quit', or Ctrl-D.", stacklevel=1)
%% 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 sepcific 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