Skip to content
Snippets Groups Projects
Commit aa1502b4 authored by Colm Talbot's avatar Colm Talbot
Browse files

Merge branch 'master' into 'master'

change "source frame" to "detector frame"

See merge request !745
parents b88c19ca dfb35135
No related branches found
No related tags found
No related merge requests found
Pipeline #117669 passed
%% Cell type:markdown id: tags:
# Compare samplers
In this notebook, we'll compare the different samplers implemented in `bilby`. As of this version, we don't compare the outputs, only how to run them and the timings for their default setup.
## Setup
%% Cell type:code id: tags:
``` python
import numpy as np
import pylab as plt
%load_ext autoreload
%autoreload 2
import bilby
bilby.utils.setup_logger()
time_duration = 1. # set the signal duration (seconds)
sampling_frequency = 4096. # set the data sampling frequency (Hz)
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, # primary dimensionless spin magnitude
a_2=0, # secondary dimensionless spin magnitude
tilt_1=0, # polar angle between primary spin and the orbital angular momentum (radians)
tilt_2=0, # polar angle between secondary spin and the orbital angular momentum
phi_12=0, # azimuthal angle between primary and secondary spin (radians)
phi_jl=0, # azimuthal angle between total angular momentum and orbital angular momentum (radians)
luminosity_distance=100., # luminosity distance to source (Mpc)
iota=0.4, # inclination angle between line of sight and orbital angular momentum (radians)
phase=1.3, # phase (radians)
waveform_approximant='IMRPhenomPv2', # waveform approximant name
reference_frequency=50., # gravitational waveform reference frequency (Hz)
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
)
# initialise 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)
# generate a frequency-domain waveform
hf_signal = waveform_generator.frequency_domain_strain()
# initialise a single interferometer representing LIGO Hanford
H1 = bilby.gw.detector.get_empty_interferometer('H1')
# set the strain data at the interferometer
H1.set_strain_data_from_power_spectral_density(sampling_frequency=sampling_frequency, duration=time_duration)
# inject the gravitational wave signal into the interferometer model
H1.inject_signal(injection_polarizations=hf_signal, parameters=injection_parameters)
IFOs = [H1]
# compute the likelihood on each of the signal parameters
likelihood = bilby.gw.likelihood.GravitationalWaveTransient(IFOs, waveform_generator)
```
%% Cell type:markdown id: tags:
## Prior
For this test, we will simply search of the sky position, setting the other parameters to their simulated values.
%% Cell type:code id: tags:
``` python
# set the priors on each of the injection parameters to be a delta function at their given value
priors = bilby.gw.prior.BBHPriorDict()
for key in injection_parameters.keys():
priors[key] = injection_parameters[key]
# now reset the priors on the sky position co-ordinates in order to conduct a sky position search
priors['ra'] = bilby.prior.Uniform(0, 2*np.pi, 'ra')
priors['dec'] = bilby.prior.Cosine(name='dec', minimum=-np.pi/2, maximum=np.pi/2)
```
%% Cell type:markdown id: tags:
## PyMultinest
%% Cell type:code id: tags:
``` python
%%time
result = bilby.core.sampler.run_sampler(
likelihood, priors=priors, sampler='pymultinest', label='pymultinest',
npoints=200, verbose=False, resume=False)
fig = result.plot_corner(save=False)
# show the corner plot
plt.show()
print(result)
```
%% Cell type:markdown id: tags:
# dynesty
%% Cell type:code id: tags:
``` python
%%time
result = bilby.core.sampler.run_sampler(
likelihood, priors=priors, sampler='dynesty', label='dynesty',
bound='multi', sample='rwalk', npoints=200, walks=1, verbose=False,
update_interval=100)
fig = result.plot_corner(save=False)
# show the corner plot
plt.show()
print(result)
```
%% Cell type:markdown id: tags:
# Dynamic Nested Sampling (Dynesty)
See [the dynesty docs](http://dynesty.readthedocs.io/en/latest/dynamic.html#). Essentially, this methods improves the posterior estimation over that of standard nested sampling.
%% Cell type:code id: tags:
``` python
%%time
result = bilby.core.sampler.run_sampler(
likelihood, priors=priors, sampler='dynesty', label='dynesty_dynamic',
bound='multi', nlive=250, sample='unif', verbose=True,
update_interval=100, dynamic=True)
fig = result.plot_corner(save=False)
# show the corner plot
plt.show()
print(result)
```
%% Cell type:markdown id: tags:
# ptemcee
%% Cell type:code id: tags:
``` python
%%time
result = bilby.core.sampler.run_sampler(
likelihood, priors=priors, sampler='ptemcee', label='ptemcee',
nwalkers=100, nsteps=200, nburn=100, ntemps=2,
tqdm='tqdm_notebook')
fig = result.plot_corner(save=False)
# show the corner plot
plt.show()
print(result)
```
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment