Skip to content
Snippets Groups Projects
Commit 4cdf9f76 authored by Matthew Carney's avatar Matthew Carney
Browse files

Changed iota to theta_jn in tutorials.

parent 3d257816
No related branches found
No related tags found
1 merge request!648updated iota to theta_jn in tutorials
%% Cell type:code id: tags:
``` python
from __future__ import division, print_function
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns
import bilby
%matplotlib inline
```
%% Output
14:35 bilby INFO : Running bilby version: 0.3.1: (CLEAN) 06cce8b 2018-10-18 05:02:54 -0500
%% Cell type:markdown id: tags:
## This notebook will show how to use the PTMCMCSampler, in particular this will highlight how to add custom jump proposals
%% Cell type:markdown id: tags:
## create 150914 like injection
%% Cell type:code id: tags:
``` python
# Set the duration and sampling frequency of the data segment that we're
# going to inject the signal into
duration = 4.
sampling_frequency = 2048.
# Specify the output directory and the name of the simulation.
outdir = 'outdir'
label = 'basic_tutorial4'
bilby.core.utils.setup_logger(outdir=outdir, label=label)
# Set up a random seed for result reproducibility. This is optional!
np.random.seed(88170235)
# We are going to inject a binary black hole waveform. We first establish a
# dictionary of parameters that includes all of the different waveform
# parameters, including masses of the two black holes (mass_1, mass_2),
# spins of both black holes (a, tilt, phi), etc.
injection_parameters = dict(
mass_1=36., mass_2=29., a_1=0.4, a_2=0.3, tilt_1=0.5, tilt_2=1.0,
phi_12=1.7, phi_jl=0.3, luminosity_distance=2000., iota=0.4, psi=2.659,
phi_12=1.7, phi_jl=0.3, luminosity_distance=2000., theta_jn=0.4, psi=2.659,
phase=1.3, geocent_time=1126259642.413, ra=1.375, dec=-1.2108)
# Fixed arguments passed into the source model
waveform_arguments = dict(waveform_approximant='IMRPhenomP',
reference_frequency=50., minimum_frequency=20.)
```
%% Cell type:markdown id: tags:
## Inject into data
%% Cell type:code id: tags:
``` python
# Create the waveform_generator using a LAL BinaryBlackHole source function
waveform_generator = bilby.gw.WaveformGenerator(
duration=duration, sampling_frequency=sampling_frequency,
frequency_domain_source_model=bilby.gw.source.lal_binary_black_hole,
waveform_arguments=waveform_arguments)
# Set up interferometers. In this case we'll use two interferometers
# (LIGO-Hanford (H1), LIGO-Livingston (L1). These default to their design
# sensitivity
ifos = bilby.gw.detector.InterferometerList(['H1', 'L1'])
ifos.set_strain_data_from_power_spectral_densities(
sampling_frequency=sampling_frequency, duration=duration,
start_time=injection_parameters['geocent_time'] - 3)
ifos.inject_signal(waveform_generator=waveform_generator,
parameters=injection_parameters)
```
%% Output
/home/c1572221/src/bilby/bilby/gw/detector.py:1986: RuntimeWarning: invalid value encountered in multiply
frequency_domain_strain = self.__power_spectral_density_interpolated(frequencies) ** 0.5 * white_noise
14:36 bilby INFO : Injected signal in H1:
14:36 bilby INFO : optimal SNR = 12.09
14:36 bilby INFO : matched filter SNR = 10.88-0.37j
14:36 bilby INFO : luminosity_distance = 2000.0
14:36 bilby INFO : psi = 2.659
14:36 bilby INFO : a_2 = 0.3
14:36 bilby INFO : a_1 = 0.4
14:36 bilby INFO : geocent_time = 1126259642.41
14:36 bilby INFO : tilt_2 = 1.0
14:36 bilby INFO : phi_jl = 0.3
14:36 bilby INFO : ra = 1.375
14:36 bilby INFO : phase = 1.3
14:36 bilby INFO : mass_2 = 29.0
14:36 bilby INFO : mass_1 = 36.0
14:36 bilby INFO : phi_12 = 1.7
14:36 bilby INFO : dec = -1.2108
14:36 bilby INFO : iota = 0.4
14:36 bilby INFO : tilt_1 = 0.5
14:36 bilby INFO : Injected signal in L1:
14:36 bilby INFO : optimal SNR = 9.79
14:36 bilby INFO : matched filter SNR = 10.00+0.07j
14:36 bilby INFO : luminosity_distance = 2000.0
14:36 bilby INFO : psi = 2.659
14:36 bilby INFO : a_2 = 0.3
14:36 bilby INFO : a_1 = 0.4
14:36 bilby INFO : geocent_time = 1126259642.41
14:36 bilby INFO : tilt_2 = 1.0
14:36 bilby INFO : phi_jl = 0.3
14:36 bilby INFO : ra = 1.375
14:36 bilby INFO : phase = 1.3
14:36 bilby INFO : mass_2 = 29.0
14:36 bilby INFO : mass_1 = 36.0
14:36 bilby INFO : phi_12 = 1.7
14:36 bilby INFO : dec = -1.2108
14:36 bilby INFO : iota = 0.4
14:36 bilby INFO : tilt_1 = 0.5
[{'cross': array([0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j]),
'plus': array([-0.-0.j, -0.-0.j, -0.-0.j, ..., -0.-0.j, -0.-0.j,
-0.-0.j])},
{'cross': array([0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j]),
'plus': array([-0.-0.j, -0.-0.j, -0.-0.j, ..., -0.-0.j, -0.-0.j,
-0.-0.j])}]
%% Cell type:markdown id: tags:
### For simplicity, we will fix all parameters here to the injected value and only vary over mass1 and mass2,
%% Cell type:code id: tags:
``` python
priors = injection_parameters.copy()
priors['mass_1'] = bilby.prior.Uniform(name='mass_1', minimum=10, maximum=80, unit=r'$M_{\\odot}$')
priors['mass_2'] = bilby.prior.Uniform(name='mass_1', minimum=10, maximum=80, unit=r'$M_{\\odot}$')
```
%% Cell type:markdown id: tags:
## Here we create arbitrary jump proposals. This will highlight the necessary features of a jump proposal in ptmcmc. That is it takes the current position, x, then outputs a new position , q, and the jump probability i.e. p(x -> q). These will then be passed to the standard metropolis hastings condition.
## The two proposals below are probably not very good ones, ideally we would use proposals based upon our kmowledge of the problem/parameter space. In general for these proposals lqxy will certainly not be 0
%% Cell type:code id: tags:
``` python
class UniformJump(object):
def __init__(self, pmin, pmax):
"""Draw random parameters from pmin, pmax"""
self.pmin = pmin
self.pmax = pmax
def unjump(self, x, it, beta):
"""
Function prototype must read in parameter vector x,
sampler iteration number it, and inverse temperature beta
"""
# log of forward-backward jump probability
lqxy = 0
# uniformly drawm parameters
q = np.random.uniform(self.pmin, self.pmax, len(x))
return q, lqxy
```
%% Cell type:code id: tags:
``` python
class NormJump(object):
def __init__(self, step_size):
"""Draw random parameters from pmin, pmax"""
self.step_size = step_size
def normjump(self, x, it, beta):
"""
Function prototype must read in parameter vector x,
sampler iteration number it, and inverse temperature beta
"""
# log of forward-backward jump probability. this is only zero for simple examples.
lqxy = 0
# uniformly drawm parameters
q = np.random.multivariate_normal(x , self.step_size * np.eye(len(x)) , 1)
return q[0], lqxy
```
%% Cell type:markdown id: tags:
### Below we create a dictionary containing our jump proposals and the relative weight of that proposal in the proposal cycle, these are then passed to bilby.run_sampler under the keyword argument custom_proposals =
%% Cell type:code id: tags:
``` python
normjump = NormJump(1)
normweight = 5
ujump = UniformJump(20, 40)
uweight = 1
custom = {'uniform': [ujump.unjump , uweight],
'normal': [normjump.normjump , normweight]}
```
%% Cell type:code id: tags:
``` python
# Initialise the likelihood by passing in the interferometer data (ifos) and
# the waveoform generator
likelihood = bilby.gw.GravitationalWaveTransient(
interferometers=ifos,waveform_generator=waveform_generator)
result = bilby.run_sampler(
likelihood=likelihood, priors=priors, sampler= 'PTMCMCsampler',custom_proposals = custom , Niter = 10**4 )
```
%% Output
14:36 bilby INFO : Running for label 'label', output will be saved to 'outdir'
14:36 bilby INFO : Search parameters:
14:36 bilby INFO : mass_2 = Uniform(minimum=10, maximum=80, name='mass_1', latex_label='$m_1$', unit='$M_{\\\\odot}$')
14:36 bilby INFO : mass_1 = Uniform(minimum=10, maximum=80, name='mass_1', latex_label='$m_1$', unit='$M_{\\\\odot}$')
14:36 bilby INFO : phi_jl = 0.3
14:36 bilby INFO : dec = -1.2108
14:36 bilby INFO : psi = 2.659
14:36 bilby INFO : a_2 = 0.3
14:36 bilby INFO : a_1 = 0.4
14:36 bilby INFO : geocent_time = 1126259642.41
14:36 bilby INFO : luminosity_distance = 2000.0
14:36 bilby INFO : ra = 1.375
14:36 bilby INFO : phase = 1.3
14:36 bilby INFO : phi_12 = 1.7
14:36 bilby INFO : tilt_2 = 1.0
14:36 bilby INFO : iota = 0.4
14:36 bilby INFO : tilt_1 = 0.5
14:36 bilby INFO : Single likelihood evaluation took 1.670e-03 s
14:36 bilby INFO : Using sampler PTMCMCSampler with kwargs {'Niter': 10000, 'Tskip': 100, 'verbose': True, 'covUpdate': 1000, 'ladder': None, 'burn': 5000, 'NUTSweight': 0, 'AMweight': 1, 'logl_grad': None, 'HMCstepsize': 0.1, 'groups': None, 'Tmax': None, 'Tmin': 1, 'HMCsteps': 300, 'p0': None, 'neff': 10000, 'logp_grad': None, 'HMCweight': 0, 'loglkwargs': {}, 'custom_proposals': {'normal': [<bound method NormJump.normjump of <__main__.NormJump object at 0x7f2ee3cf7e10>>, 5], 'uniform': [<bound method UniformJump.unjump of <__main__.UniformJump object at 0x7f2ee3cf7d50>>, 1]}, 'logpargs': {}, 'MALAweight': 0, 'loglargs': {}, 'DEweight': 1, 'thin': 1, 'isave': 1000, 'outDir': None, 'SCAMweight': 1, 'logpkwargs': {}}
14:36 bilby INFO : Adding normal to proposals with weight 5
14:36 bilby INFO : Adding uniform to proposals with weight 1
Finished 50.00 percent in 26.952440 s Acceptance rate = 0.1046667Adding DE jump with weight 1
Finished 90.00 percent in 47.835580 s Acceptance rate = 0.140667
14:37 bilby INFO : Sampling time: 0:00:53.551375
14:37 bilby ERROR :
Saving the data has failed with the following message:
Can't pickle <type 'instancemethod'>: attribute lookup __builtin__.instancemethod failed
14:37 bilby INFO : Results saved to outdir/
14:37 bilby INFO : Summary of results:
nsamples: 4000
log_noise_evidence: -8086.279
log_evidence: nan +/- nan
log_bayes_factor: nan +/- nan
Run Complete
%% Cell type:code id: tags:
``` python
result.plot_corner()
```
%% Output
<Figure size 396x396 with 4 Axes>
%% Cell type:markdown id: tags:
### PTMCMC produces the acceptance rate for each of the proposals (including the ones built in). This is taken as an average at a specified checkpoint. This is one (acceptnace rate is certainly not the only/even the best metric here. Think exploration v exploitation problem ) indicators of whether our jump proposal is a good one
%% Cell type:code id: tags:
``` python
sampler_meta = result.meta_data['sampler_meta']
jumps = sampler_meta['proposals']
```
%% Cell type:code id: tags:
``` python
plt.figure()
plt.xlabel('epoch')
plt.ylabel('acceptance rate')
for i,proposal in enumerate(jumps):
plt.plot(jumps[proposal] , label = proposal)
plt.legend(loc='best', frameon=True)
```
%% Output
<matplotlib.legend.Legend at 0x7f2ed97a0f50>
%% Cell type:markdown id: tags:
## We can generate the 1d chains for each of the parameters too and the likelihood of those points on the chain
%% Cell type:code id: tags:
``` python
m2 = result.posterior.mass_2.values
m1 = result.posterior.mass_1.values
fig, ax = plt.subplots(nrows = 2 , ncols =1 , sharex = True , figsize = (8,8))
ax[0].plot(m1 , 'o', label = 'm1' )
ax[0].plot(m2 , 'o', label = 'm2' )
ax[0].set_ylabel(r'$M_{\odot}$')
ax[0].legend(loc = 'best' , frameon = True , fontsize = 12)
ax[1].plot(result.log_likelihood_evaluations)
ax[1].set_ylabel(r'$\mathcal{L}$')
ax[1].set_xlabel('iterations')
ax[1].set_xscale('log')
```
%% Output
%% Cell type:code id: tags:
``` python
```
......
%% 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)
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)
theta_jn=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