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

Merge branch 'add-gwin-example' into 'master'

Adds a simple example of using gwin for an injection

See merge request Monash/tupak!166
parents 9223a8b7 a33ffc79
No related branches found
No related tags found
1 merge request!166Adds a simple example of using gwin for an injection
Pipeline #30264 passed with warnings
#!/bin/python
"""
An example of how to use tupak with `gwin` (https://github.com/gwastro/gwin) to
perform CBC parameter estimation.
To run this example, it is sufficient to use the pip-installable pycbc package,
but the source installation of gwin. You can install this by cloning the
repository (https://github.com/gwastro/gwin) and running
$ python setup.py install
A practical difference between gwin and tupak is that while fixed parameters
are specified via the prior in tupak, in gwin, these are fixed at instantiation
of the model. So, in the following, we only create priors for the parameters
to be searched over.
"""
from __future__ import division, print_function
import numpy as np
import tupak
import gwin
from pycbc import psd as pypsd
from pycbc.waveform.generator import (FDomainDetFrameGenerator,
FDomainCBCGenerator)
label = 'using_gwin'
# Search priors
priors = dict()
priors['distance'] = tupak.core.prior.Uniform(500, 2000, 'distance')
priors['polarization'] = tupak.core.prior.Uniform(0, np.pi, 'iota')
# Data variables
seglen = 4
sample_rate = 2048
N = seglen * sample_rate / 2 + 1
fmin = 30.
# Injected signal variables
injection_parameters = dict(mass1=38.6, mass2=29.3, spin1z=0, spin2z=0,
tc=0, ra=3.1, dec=1.37, polarization=2.76,
distance=1500)
# These lines figure out which parameters are to be searched over
variable_parameters = priors.keys()
fixed_parameters = injection_parameters.copy()
for key in priors:
fixed_parameters.pop(key)
# These lines generate the `model` object - see https://gwin.readthedocs.io/en/latest/api/gwin.models.gaussian_noise.html
generator = FDomainDetFrameGenerator(
FDomainCBCGenerator, 0.,
variable_args=variable_parameters, detectors=['H1', 'L1'],
delta_f=1. / seglen, f_lower=fmin,
approximant='IMRPhenomPv2', **fixed_parameters)
signal = generator.generate(**injection_parameters)
psd = pypsd.aLIGOZeroDetHighPower(int(N), 1. / seglen, 20.)
psds = {'H1': psd, 'L1': psd}
model = gwin.models.GaussianNoise(
variable_parameters, signal, generator, fmin, psds=psds)
model.update(**injection_parameters)
# This create a dummy class to convert the model into a tupak.likelihood object
class GWINLikelihood(tupak.core.likelihood.Likelihood):
def __init__(self, model):
""" A likelihood to wrap around GWIN model objects
Parameters
----------
model: gwin.model.GaussianNoise
A gwin model
"""
self.model = model
self.parameters = {x: None for x in self.model.variable_params}
def log_likelihood(self):
self.model.update(**self.parameters)
return self.model.loglikelihood
likelihood = GWINLikelihood(model)
# Now run the inference
result = tupak.run_sampler(
likelihood=likelihood, priors=priors, sampler='dynesty', npoints=500,
label=label)
result.plot_corner()
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