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

Fix theta jn definition

parent dded5f6c
No related branches found
No related tags found
1 merge request!777Resolve "Double check `theta_jn` and `iota` consistency in the docs"
Pipeline #123469 passed
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
# Compare samplers # 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. 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 ## Setup
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
import numpy as np import numpy as np
import pylab as plt import pylab as plt
%load_ext autoreload %load_ext autoreload
%autoreload 2 %autoreload 2
import bilby import bilby
bilby.utils.setup_logger() bilby.utils.setup_logger()
time_duration = 1. # set the signal duration (seconds) time_duration = 1. # set the signal duration (seconds)
sampling_frequency = 4096. # set the data sampling frequency (Hz) sampling_frequency = 4096. # set the data sampling frequency (Hz)
injection_parameters = dict( injection_parameters = dict(
mass_1=36., # detector frame (redshifted) primary mass (solar masses) mass_1=36., # detector frame (redshifted) primary mass (solar masses)
mass_2=29., # detector frame (redshifted) secondary mass (solar masses) mass_2=29., # detector frame (redshifted) secondary mass (solar masses)
a_1=0, # primary dimensionless spin magnitude a_1=0, # primary dimensionless spin magnitude
a_2=0, # secondary 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_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 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_12=0, # azimuthal angle between primary and secondary spin (radians)
phi_jl=0, # azimuthal angle between total angular momentum and orbital angular momentum (radians) phi_jl=0, # azimuthal angle between total angular momentum and orbital angular momentum (radians)
luminosity_distance=100., # luminosity distance to source (Mpc) luminosity_distance=100., # luminosity distance to source (Mpc)
theta_jn=0.4, # inclination angle between line of sight and orbital angular momentum (radians) theta_jn=0.4, # angle between the total angular momentum (both spin and orbital) and the line of sight
phase=1.3, # phase (radians) phase=1.3, # phase (radians)
waveform_approximant='IMRPhenomPv2', # waveform approximant name waveform_approximant='IMRPhenomPv2', # waveform approximant name
reference_frequency=50., # gravitational waveform reference frequency (Hz) reference_frequency=50., # gravitational waveform reference frequency (Hz)
ra=1.375, # source right ascension (radians) ra=1.375, # source right ascension (radians)
dec=-1.2108, # source declination (radians) dec=-1.2108, # source declination (radians)
geocent_time=1126259642.413, # reference time at geocentre (time of coalescence or peak amplitude) (GPS seconds) geocent_time=1126259642.413, # reference time at geocentre (time of coalescence or peak amplitude) (GPS seconds)
psi=2.659 # gravitational wave polarisation angle psi=2.659 # gravitational wave polarisation angle
) )
# initialise the waveform generator # initialise the waveform generator
waveform_generator = bilby.gw.waveform_generator.WaveformGenerator( waveform_generator = bilby.gw.waveform_generator.WaveformGenerator(
sampling_frequency=sampling_frequency, sampling_frequency=sampling_frequency,
duration=time_duration, duration=time_duration,
frequency_domain_source_model=bilby.gw.source.lal_binary_black_hole, frequency_domain_source_model=bilby.gw.source.lal_binary_black_hole,
parameters=injection_parameters) parameters=injection_parameters)
# generate a frequency-domain waveform # generate a frequency-domain waveform
hf_signal = waveform_generator.frequency_domain_strain() hf_signal = waveform_generator.frequency_domain_strain()
# initialise a single interferometer representing LIGO Hanford # initialise a single interferometer representing LIGO Hanford
H1 = bilby.gw.detector.get_empty_interferometer('H1') H1 = bilby.gw.detector.get_empty_interferometer('H1')
# set the strain data at the interferometer # set the strain data at the interferometer
H1.set_strain_data_from_power_spectral_density(sampling_frequency=sampling_frequency, duration=time_duration) H1.set_strain_data_from_power_spectral_density(sampling_frequency=sampling_frequency, duration=time_duration)
# inject the gravitational wave signal into the interferometer model # inject the gravitational wave signal into the interferometer model
H1.inject_signal(injection_polarizations=hf_signal, parameters=injection_parameters) H1.inject_signal(injection_polarizations=hf_signal, parameters=injection_parameters)
IFOs = [H1] IFOs = [H1]
# compute the likelihood on each of the signal parameters # compute the likelihood on each of the signal parameters
likelihood = bilby.gw.likelihood.GravitationalWaveTransient(IFOs, waveform_generator) likelihood = bilby.gw.likelihood.GravitationalWaveTransient(IFOs, waveform_generator)
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
## Prior ## Prior
For this test, we will simply search of the sky position, setting the other parameters to their simulated values. For this test, we will simply search of the sky position, setting the other parameters to their simulated values.
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
# set the priors on each of the injection parameters to be a delta function at their given value # set the priors on each of the injection parameters to be a delta function at their given value
priors = bilby.gw.prior.BBHPriorDict() priors = bilby.gw.prior.BBHPriorDict()
for key in injection_parameters.keys(): for key in injection_parameters.keys():
priors[key] = injection_parameters[key] priors[key] = injection_parameters[key]
# now reset the priors on the sky position co-ordinates in order to conduct a sky position search # 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['ra'] = bilby.prior.Uniform(0, 2*np.pi, 'ra')
priors['dec'] = bilby.prior.Cosine(name='dec', minimum=-np.pi/2, maximum=np.pi/2) priors['dec'] = bilby.prior.Cosine(name='dec', minimum=-np.pi/2, maximum=np.pi/2)
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
## PyMultinest ## PyMultinest
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
%%time %%time
result = bilby.core.sampler.run_sampler( result = bilby.core.sampler.run_sampler(
likelihood, priors=priors, sampler='pymultinest', label='pymultinest', likelihood, priors=priors, sampler='pymultinest', label='pymultinest',
npoints=200, verbose=False, resume=False) npoints=200, verbose=False, resume=False)
fig = result.plot_corner(save=False) fig = result.plot_corner(save=False)
# show the corner plot # show the corner plot
plt.show() plt.show()
print(result) print(result)
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
# dynesty # dynesty
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
%%time %%time
result = bilby.core.sampler.run_sampler( result = bilby.core.sampler.run_sampler(
likelihood, priors=priors, sampler='dynesty', label='dynesty', likelihood, priors=priors, sampler='dynesty', label='dynesty',
bound='multi', sample='rwalk', npoints=200, walks=1, verbose=False, bound='multi', sample='rwalk', npoints=200, walks=1, verbose=False,
update_interval=100) update_interval=100)
fig = result.plot_corner(save=False) fig = result.plot_corner(save=False)
# show the corner plot # show the corner plot
plt.show() plt.show()
print(result) print(result)
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
# Dynamic Nested Sampling (Dynesty) # 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. 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: %% Cell type:code id: tags:
``` python ``` python
%%time %%time
result = bilby.core.sampler.run_sampler( result = bilby.core.sampler.run_sampler(
likelihood, priors=priors, sampler='dynesty', label='dynesty_dynamic', likelihood, priors=priors, sampler='dynesty', label='dynesty_dynamic',
bound='multi', nlive=250, sample='unif', verbose=True, bound='multi', nlive=250, sample='unif', verbose=True,
update_interval=100, dynamic=True) update_interval=100, dynamic=True)
fig = result.plot_corner(save=False) fig = result.plot_corner(save=False)
# show the corner plot # show the corner plot
plt.show() plt.show()
print(result) print(result)
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
# ptemcee # ptemcee
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
%%time %%time
result = bilby.core.sampler.run_sampler( result = bilby.core.sampler.run_sampler(
likelihood, priors=priors, sampler='ptemcee', label='ptemcee', likelihood, priors=priors, sampler='ptemcee', label='ptemcee',
nwalkers=100, nsteps=200, nburn=100, ntemps=2, nwalkers=100, nsteps=200, nburn=100, ntemps=2,
tqdm='tqdm_notebook') tqdm='tqdm_notebook')
fig = result.plot_corner(save=False) fig = result.plot_corner(save=False)
# show the corner plot # show the corner plot
plt.show() plt.show()
print(result) 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