diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index 14215c6f183db9d846f38f82067ecd0f7878f1f4..c2c49cea108d3c1ae687cabc8bbc479cd9b6c596 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -9,7 +9,7 @@ repos:
hooks:
- id: black
language_version: python3
- files: '(^bilby/bilby_mcmc/|^examples/core_examples/|examples/gw_examples/data_examples)'
+ files: '(^bilby/bilby_mcmc/|^examples/)'
- repo: https://github.com/codespell-project/codespell
rev: v2.1.0
hooks:
@@ -20,7 +20,7 @@ repos:
hooks:
- id: isort # sort imports alphabetically and separates import into sections
args: [-w=88, -m=3, -tc, -sp=setup.cfg ]
- files: '(^bilby/bilby_mcmc/|^examples/core_examples/examples/gw_examples/data_examples)'
+ files: '(^bilby/bilby_mcmc/|^examples/)'
- repo: https://github.com/datarootsio/databooks
rev: 0.1.14
hooks:
diff --git a/bilby/gw/source.py b/bilby/gw/source.py
index 5ea9415e98fb0bb83bc8e3fb214c2d4654d8c549..3c335aaeade8f40df71bacdd1ef2d00d368ada9a 100644
--- a/bilby/gw/source.py
+++ b/bilby/gw/source.py
@@ -885,8 +885,7 @@ def sinegaussian(frequency_array, hrss, Q, frequency, **kwargs):
return{'plus': h_plus, 'cross': h_cross}
-def supernova(
- frequency_array, realPCs, imagPCs, file_path, luminosity_distance, **kwargs):
+def supernova(frequency_array, luminosity_distance, **kwargs):
"""
A source model that reads a simulation from a text file.
@@ -897,8 +896,6 @@ def supernova(
----------
frequency_array: array-like
Unused
- realPCs: UNUSED
- imagPCs: UNUSED
file_path: str
Path to the file containing the NR simulation. The format of this file
should be readable by :code:`numpy.loadtxt` and have four columns
@@ -907,7 +904,8 @@ def supernova(
luminosity_distance: float
The distance to the source in kpc, this scales the amplitude of the
signal. The simulation is assumed to be at 10kpc.
- kwargs: UNUSED
+ kwargs:
+ extra keyword arguments, this should include the :code:`file_path`
Returns
-------
@@ -915,20 +913,20 @@ def supernova(
A dictionary containing the plus and cross components of the signal.
"""
- realhplus, imaghplus, realhcross, imaghcross = np.loadtxt(
- file_path, usecols=(0, 1, 2, 3), unpack=True)
+ file_path = kwargs["file_path"]
+ data = np.genfromtxt(file_path)
# waveform in file at 10kpc
scaling = 1e-3 * (10.0 / luminosity_distance)
- h_plus = scaling * (realhplus + 1.0j * imaghplus)
- h_cross = scaling * (realhcross + 1.0j * imaghcross)
+ h_plus = scaling * (data[:, 0] + 1j * data[:, 1])
+ h_cross = scaling * (data[:, 2] + 1j * data[:, 3])
return {'plus': h_plus, 'cross': h_cross}
def supernova_pca_model(
- frequency_array, pc_coeff1, pc_coeff2, pc_coeff3, pc_coeff4, pc_coeff5,
- luminosity_distance, **kwargs):
+ frequency_array, pc_coeff1, pc_coeff2, pc_coeff3, pc_coeff4, pc_coeff5, luminosity_distance, **kwargs
+):
r"""
Signal model based on a five-component principal component decomposition
of a model.
@@ -966,22 +964,19 @@ def supernova_pca_model(
The plus and cross polarizations of the signal
"""
- realPCs = kwargs['realPCs']
- imagPCs = kwargs['imagPCs']
+ principal_components = kwargs["realPCs"] + 1j * kwargs["imagPCs"]
+ coefficients = [pc_coeff1, pc_coeff2, pc_coeff3, pc_coeff4, pc_coeff5]
- pc1 = realPCs[:, 0] + 1.0j * imagPCs[:, 0]
- pc2 = realPCs[:, 1] + 1.0j * imagPCs[:, 1]
- pc3 = realPCs[:, 2] + 1.0j * imagPCs[:, 2]
- pc4 = realPCs[:, 3] + 1.0j * imagPCs[:, 3]
- pc5 = realPCs[:, 4] + 1.0j * imagPCs[:, 5]
+ strain = np.sum(
+ [coeff * principal_components[:, ii] for ii, coeff in enumerate(coefficients)],
+ axis=0
+ )
# file at 10kpc
scaling = 1e-23 * (10.0 / luminosity_distance)
- h_plus = scaling * (pc_coeff1 * pc1 + pc_coeff2 * pc2 + pc_coeff3 * pc3 +
- pc_coeff4 * pc4 + pc_coeff5 * pc5)
- h_cross = scaling * (pc_coeff1 * pc1 + pc_coeff2 * pc2 + pc_coeff3 * pc3 +
- pc_coeff4 * pc4 + pc_coeff5 * pc5)
+ h_plus = scaling * strain
+ h_cross = scaling * strain
return {'plus': h_plus, 'cross': h_cross}
diff --git a/examples/core_examples/gaussian_process_celerite_example.py b/examples/core_examples/gaussian_process_celerite_example.py
index 417493f57745ba211b333d70094d731bed0291cf..574566a9af00c4208667f62b11572ff3892afc76 100644
--- a/examples/core_examples/gaussian_process_celerite_example.py
+++ b/examples/core_examples/gaussian_process_celerite_example.py
@@ -1,13 +1,11 @@
-import matplotlib.pyplot as plt
-import numpy as np
from pathlib import Path
-import celerite.terms
-
import bilby
+import celerite.terms
+import matplotlib.pyplot as plt
+import numpy as np
from bilby.core.prior import Uniform
-
# In this example we show how we can use the `celerite` package within `bilby`.
# We begin by synthesizing some data and then use a simple Gaussian Process
# model to fit and interpolate the data.
diff --git a/examples/core_examples/gaussian_process_george_example.py b/examples/core_examples/gaussian_process_george_example.py
index 2c2946583cbe7ab32aa4ddd106e1b2f44559e76e..969a56a72f5182962a8f7087041b57e791c48248 100644
--- a/examples/core_examples/gaussian_process_george_example.py
+++ b/examples/core_examples/gaussian_process_george_example.py
@@ -1,13 +1,11 @@
-import matplotlib.pyplot as plt
-import numpy as np
from pathlib import Path
-import george
-
import bilby
+import george
+import matplotlib.pyplot as plt
+import numpy as np
from bilby.core.prior import Uniform
-
# In this example we show how we can use the `george` package within
# `bilby`. We begin by synthesizing some data and then use a simple Gaussian
# Process model to fit and interpolate the data. `bilby` implements a
diff --git a/examples/core_examples/hyper_parameter_example.py b/examples/core_examples/hyper_parameter_example.py
index 638d47edabaa280881bfc04703f5504402cbb5f8..288e9c21b07aa254da6725f57d95b7e3eb05fa99 100644
--- a/examples/core_examples/hyper_parameter_example.py
+++ b/examples/core_examples/hyper_parameter_example.py
@@ -4,7 +4,6 @@ An example of how to use bilby to perform parameter estimation for hyper params
"""
import matplotlib.pyplot as plt
import numpy as np
-
from bilby.core.likelihood import GaussianLikelihood
from bilby.core.prior import Uniform
from bilby.core.result import make_pp_plot
diff --git a/examples/core_examples/slabspike_example.py b/examples/core_examples/slabspike_example.py
index 21df0790ff55d88a183686ec76c0f5c96b08d3d2..2c42c174bda64e8535017cfce190850e3e9faad5 100644
--- a/examples/core_examples/slabspike_example.py
+++ b/examples/core_examples/slabspike_example.py
@@ -13,11 +13,10 @@ To install `PyMultiNest` call
$ conda install -c conda-forge pymultinest
"""
+import bilby
import matplotlib.pyplot as plt
import numpy as np
-import bilby
-
outdir = "outdir"
label = "slabspike"
bilby.utils.check_directory_exists_and_if_not_mkdir(outdir)
diff --git a/examples/gw_examples/injection_examples/australian_detector.py b/examples/gw_examples/injection_examples/australian_detector.py
index 6b903a9fc0695a279ef777e0d0fa166aa7fae140..85c5cf03e7cb640207e27e731b95c1b3acae7820 100644
--- a/examples/gw_examples/injection_examples/australian_detector.py
+++ b/examples/gw_examples/injection_examples/australian_detector.py
@@ -3,52 +3,54 @@
Tutorial to demonstrate a new interferometer
We place a new instrument in Gingin, with an A+ sensitivity in a network of A+
-interferometers at Hanford and Livingston
-"""
+interferometers at Hanford and Livingston.
-import numpy as np
+This requires :code:`gwinc` to be installed. This is available via conda-forge.
+"""
import bilby
import gwinc
+import numpy as np
# Set the duration and sampling frequency of the data segment that we're going
# to inject the signal into
-duration = 4.
-sampling_frequency = 2048.
+duration = 4
+sampling_frequency = 1024
# Specify the output directory and the name of the simulation.
-outdir = 'outdir'
-label = 'australian_detector'
+outdir = "outdir"
+label = "australian_detector"
bilby.core.utils.setup_logger(outdir=outdir, label=label)
# Set up a random seed for result reproducibility. This is optional!
np.random.seed(88170232)
# create a new detector using a PyGwinc sensitivity curve
-frequencies = np.logspace(0, 3, 1000)
-budget, gwinc_ifo, _, _ = gwinc.load_ifo('Aplus')
-gwinc_ifo = gwinc.precompIFO(frequencies, gwinc_ifo)
-gwinc_traces = budget(frequencies, ifo=gwinc_ifo).calc_trace()
-gwinc_noises = {n: d[0] for n, d in gwinc_traces.items()}
-
-Aplus_psd = gwinc_noises['Total']
+curve = gwinc.load_budget("Aplus").run()
# Set up the detector as a four-kilometer detector in Gingin
# The location of this detector is not defined in Bilby, so we need to add it
AusIFO = bilby.gw.detector.Interferometer(
power_spectral_density=bilby.gw.detector.PowerSpectralDensity(
- frequency_array=frequencies, psd_array=Aplus_psd),
- name='AusIFO', length=4,
- minimum_frequency=min(frequencies), maximum_frequency=max(frequencies),
- latitude=-31.34, longitude=115.91,
- elevation=0., xarm_azimuth=2., yarm_azimuth=125.)
+ frequency_array=curve.freq, psd_array=curve.psd
+ ),
+ name="AusIFO",
+ length=4,
+ minimum_frequency=20,
+ maximum_frequency=sampling_frequency / 2,
+ latitude=-31.34,
+ longitude=115.91,
+ elevation=0.0,
+ xarm_azimuth=2.0,
+ yarm_azimuth=125.0,
+)
# Set up two other detectors at Hanford and Livingston
-interferometers = bilby.gw.detector.InterferometerList(['H1', 'L1'])
-for interferometer in interferometers:
- interferometer.power_spectral_density =\
- bilby.gw.detector.PowerSpectralDensity(
- frequency_array=frequencies, psd_array=Aplus_psd)
+interferometers = bilby.gw.detector.InterferometerList(["H1", "L1"])
+for ifo in interferometers:
+ ifo.power_spectral_density = bilby.gw.detector.PowerSpectralDensity(
+ frequency_array=curve.freq, psd_array=curve.psd
+ )
# append the Australian detector to the list of other detectors
interferometers.append(AusIFO)
@@ -58,53 +60,88 @@ interferometers.append(AusIFO)
# as we're using a three-detector network of A+, we inject a GW150914-like
# signal at 4 Gpc
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=4000., theta_jn=0.4, psi=2.659,
- phase=1.3, geocent_time=1126259642.413, ra=1.375, dec=0.2108)
+ mass_1=36.0,
+ mass_2=29.0,
+ 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=4000.0,
+ theta_jn=0.4,
+ psi=2.659,
+ phase=1.3,
+ geocent_time=1126259642.413,
+ ra=1.375,
+ dec=0.2108,
+)
# Fixed arguments passed into the source model
-waveform_arguments = dict(waveform_approximant='IMRPhenomPv2',
- reference_frequency=50.)
+waveform_arguments = dict(waveform_approximant="IMRPhenomXP", reference_frequency=50.0)
# Create the waveform_generator using a LAL BinaryBlackHole source function
waveform_generator = bilby.gw.WaveformGenerator(
- duration=duration, sampling_frequency=sampling_frequency,
+ duration=duration,
+ sampling_frequency=sampling_frequency,
frequency_domain_source_model=bilby.gw.source.lal_binary_black_hole,
- waveform_arguments=waveform_arguments)
+ waveform_arguments=waveform_arguments,
+)
-start_time = injection_parameters['geocent_time'] + 2 - duration
+start_time = injection_parameters["geocent_time"] + 2 - duration
# inject the signal into the interferometers
-for interferometer in interferometers:
- interferometer.set_strain_data_from_power_spectral_density(
- sampling_frequency=sampling_frequency, duration=duration)
- interferometer.inject_signal(
- parameters=injection_parameters, waveform_generator=waveform_generator)
+for ifo in interferometers:
+ ifo.set_strain_data_from_power_spectral_density(
+ sampling_frequency=sampling_frequency, duration=duration
+ )
+ ifo.inject_signal(
+ parameters=injection_parameters, waveform_generator=waveform_generator
+ )
# plot the data for sanity
- signal = interferometer.get_detector_response(
- waveform_generator.frequency_domain_strain(), injection_parameters)
- interferometer.plot_data(signal=signal, outdir=outdir, label=label)
+ signal = ifo.get_detector_response(
+ waveform_generator.frequency_domain_strain(), injection_parameters
+ )
+ ifo.plot_data(signal=signal, outdir=outdir, label=label)
# set up priors
priors = bilby.gw.prior.BBHPriorDict()
-for key in ['a_1', 'a_2', 'tilt_1', 'tilt_2', 'phi_12', 'phi_jl', 'psi',
- 'geocent_time', 'phase']:
+for key in [
+ "a_1",
+ "a_2",
+ "tilt_1",
+ "tilt_2",
+ "phi_12",
+ "phi_jl",
+ "psi",
+ "geocent_time",
+ "phase",
+ "ra",
+ "dec",
+ "luminosity_distance",
+ "theta_jn",
+]:
priors[key] = injection_parameters[key]
# Initialise the likelihood by passing in the interferometer data (IFOs)
# and the waveoform generator
likelihood = bilby.gw.GravitationalWaveTransient(
- interferometers=interferometers, waveform_generator=waveform_generator,
- time_marginalization=False, phase_marginalization=False,
- distance_marginalization=False, priors=priors)
+ interferometers=interferometers,
+ waveform_generator=waveform_generator,
+)
result = bilby.run_sampler(
- likelihood=likelihood, priors=priors, npoints=1000,
- injection_parameters=injection_parameters, outdir=outdir, label=label)
+ likelihood=likelihood,
+ priors=priors,
+ npoints=1000,
+ injection_parameters=injection_parameters,
+ outdir=outdir,
+ label=label,
+)
# make some plots of the outputs
result.plot_corner()
diff --git a/examples/gw_examples/injection_examples/binary_neutron_star_example.py b/examples/gw_examples/injection_examples/binary_neutron_star_example.py
index e6aa18a9a47c9594f05a446d47d6294a31fc6032..90308c3c5b3d2bc973335eb1658647af685a6380 100644
--- a/examples/gw_examples/injection_examples/binary_neutron_star_example.py
+++ b/examples/gw_examples/injection_examples/binary_neutron_star_example.py
@@ -9,13 +9,12 @@ tidal deformabilities
"""
-import numpy as np
-
import bilby
+import numpy as np
# Specify the output directory and the name of the simulation.
-outdir = 'outdir'
-label = 'bns_example'
+outdir = "outdir"
+label = "bns_example"
bilby.core.utils.setup_logger(outdir=outdir, label=label)
# Set up a random seed for result reproducibility. This is optional!
@@ -26,39 +25,56 @@ np.random.seed(88170235)
# parameters, including masses of the two black holes (mass_1, mass_2),
# aligned spins of both black holes (chi_1, chi_2), etc.
injection_parameters = dict(
- mass_1=1.5, mass_2=1.3, chi_1=0.02, chi_2=0.02, luminosity_distance=50.,
- theta_jn=0.4, psi=2.659, phase=1.3, geocent_time=1126259642.413,
- ra=1.375, dec=-1.2108, lambda_1=400, lambda_2=450)
+ mass_1=1.5,
+ mass_2=1.3,
+ chi_1=0.02,
+ chi_2=0.02,
+ luminosity_distance=50.0,
+ theta_jn=0.4,
+ psi=2.659,
+ phase=1.3,
+ geocent_time=1126259642.413,
+ ra=1.375,
+ dec=-1.2108,
+ lambda_1=400,
+ lambda_2=450,
+)
# Set the duration and sampling frequency of the data segment that we're going
# to inject the signal into. For the
# TaylorF2 waveform, we cut the signal close to the isco frequency
duration = 32
-sampling_frequency = 2 * 1024
-start_time = injection_parameters['geocent_time'] + 2 - duration
+sampling_frequency = 2048
+start_time = injection_parameters["geocent_time"] + 2 - duration
# Fixed arguments passed into the source model. The analysis starts at 40 Hz.
-waveform_arguments = dict(waveform_approximant='IMRPhenomPv2_NRTidal',
- reference_frequency=50., minimum_frequency=40.0)
+waveform_arguments = dict(
+ waveform_approximant="IMRPhenomPv2_NRTidal",
+ reference_frequency=50.0,
+ minimum_frequency=40.0,
+)
# Create the waveform_generator using a LAL Binary Neutron Star source function
waveform_generator = bilby.gw.WaveformGenerator(
- duration=duration, sampling_frequency=sampling_frequency,
+ duration=duration,
+ sampling_frequency=sampling_frequency,
frequency_domain_source_model=bilby.gw.source.lal_binary_neutron_star,
parameter_conversion=bilby.gw.conversion.convert_to_lal_binary_neutron_star_parameters,
- waveform_arguments=waveform_arguments)
+ waveform_arguments=waveform_arguments,
+)
# Set up interferometers. In this case we'll use three interferometers
# (LIGO-Hanford (H1), LIGO-Livingston (L1), and Virgo (V1)).
# These default to their design sensitivity and start at 40 Hz.
-interferometers = bilby.gw.detector.InterferometerList(['H1', 'L1', 'V1'])
+interferometers = bilby.gw.detector.InterferometerList(["H1", "L1", "V1"])
for interferometer in interferometers:
interferometer.minimum_frequency = 40
interferometers.set_strain_data_from_power_spectral_densities(
- sampling_frequency=sampling_frequency, duration=duration,
- start_time=start_time)
-interferometers.inject_signal(parameters=injection_parameters,
- waveform_generator=waveform_generator)
+ sampling_frequency=sampling_frequency, duration=duration, start_time=start_time
+)
+interferometers.inject_signal(
+ parameters=injection_parameters, waveform_generator=waveform_generator
+)
# Load the default prior for binary neutron stars.
# We're going to sample in chirp_mass, symmetric_mass_ratio, lambda_tilde, and
@@ -66,31 +82,45 @@ interferometers.inject_signal(parameters=injection_parameters,
# BNS have aligned spins by default, if you want to allow precessing spins
# pass aligned_spin=False to the BNSPriorDict
priors = bilby.gw.prior.BNSPriorDict()
-for key in ['psi', 'geocent_time', 'ra', 'dec', 'chi_1', 'chi_2',
- 'theta_jn', 'luminosity_distance', 'phase']:
+for key in [
+ "psi",
+ "geocent_time",
+ "ra",
+ "dec",
+ "chi_1",
+ "chi_2",
+ "theta_jn",
+ "luminosity_distance",
+ "phase",
+]:
priors[key] = injection_parameters[key]
-priors.pop('mass_ratio')
-priors.pop('lambda_1')
-priors.pop('lambda_2')
-priors['chirp_mass'] = bilby.core.prior.Gaussian(
- 1.215, 0.1, name='chirp_mass', unit='$M_{\\odot}$')
-priors['symmetric_mass_ratio'] = bilby.core.prior.Uniform(
- 0.1, 0.25, name='symmetric_mass_ratio')
-priors['lambda_tilde'] = bilby.core.prior.Uniform(0, 5000, name='lambda_tilde')
-priors['delta_lambda'] = bilby.core.prior.Uniform(
- -5000, 5000, name='delta_lambda')
+del priors["mass_ratio"], priors["lambda_1"], priors["lambda_2"]
+priors["chirp_mass"] = bilby.core.prior.Gaussian(
+ 1.215, 0.1, name="chirp_mass", unit="$M_{\\odot}$"
+)
+priors["symmetric_mass_ratio"] = bilby.core.prior.Uniform(
+ 0.1, 0.25, name="symmetric_mass_ratio"
+)
+priors["lambda_tilde"] = bilby.core.prior.Uniform(0, 5000, name="lambda_tilde")
+priors["delta_lambda"] = bilby.core.prior.Uniform(-5000, 5000, name="delta_lambda")
# Initialise the likelihood by passing in the interferometer data (IFOs)
# and the waveform generator
likelihood = bilby.gw.GravitationalWaveTransient(
- interferometers=interferometers, waveform_generator=waveform_generator,
- time_marginalization=False, phase_marginalization=False,
- distance_marginalization=False, priors=priors)
+ interferometers=interferometers,
+ waveform_generator=waveform_generator,
+)
# Run sampler. In this case we're going to use the `nestle` sampler
result = bilby.run_sampler(
- likelihood=likelihood, priors=priors, sampler='nestle', npoints=100,
- injection_parameters=injection_parameters, outdir=outdir, label=label,
- conversion_function=bilby.gw.conversion.generate_all_bns_parameters)
+ likelihood=likelihood,
+ priors=priors,
+ sampler="nestle",
+ npoints=100,
+ injection_parameters=injection_parameters,
+ outdir=outdir,
+ label=label,
+ conversion_function=bilby.gw.conversion.generate_all_bns_parameters,
+)
result.plot_corner()
diff --git a/examples/gw_examples/injection_examples/bns_eos_example.py b/examples/gw_examples/injection_examples/bns_eos_example.py
index 5e32138b94bf6b6f402f7618fc07a0d12bcc3c61..881cc5ae8f09f02e2d2c2b705023a5fe5ca4207b 100644
--- a/examples/gw_examples/injection_examples/bns_eos_example.py
+++ b/examples/gw_examples/injection_examples/bns_eos_example.py
@@ -1,22 +1,20 @@
#!/usr/bin/env python
"""
Tutorial to demonstrate running parameter estimation on a binary neutron star
-system taking into account tidal deformabilities.
+system taking into account tidal deformabilities with a physically motivated
+model for the tidal deformabilities.
-This example estimates the masses using a uniform prior in both component masses
-and also estimates the tidal deformabilities using a uniform prior in both
-tidal deformabilities
+WARNING: The code is extremely slow.
"""
-import numpy as np
-
import bilby
-from bilby.gw.eos import TabularEOS, EOSFamily
+import numpy as np
+from bilby.gw.eos import EOSFamily, TabularEOS
# Specify the output directory and the name of the simulation.
-outdir = 'outdir'
-label = 'bns_eos_example'
+outdir = "outdir"
+label = "bns_eos_example"
bilby.core.utils.setup_logger(outdir=outdir, label=label)
# Set up a random seed for result reproducibility. This is optional!
@@ -32,7 +30,7 @@ np.random.seed(88170235)
# assuming a specific equation of state, and calculating
# corresponding tidal deformability parameters from the EoS and
# masses.
-mpa1_eos = TabularEOS('MPA1')
+mpa1_eos = TabularEOS("MPA1")
mpa1_fam = EOSFamily(mpa1_eos)
mass_1 = 1.5
@@ -42,41 +40,58 @@ lambda_2 = mpa1_fam.lambda_from_mass(mass_2)
injection_parameters = dict(
- mass_1=mass_1, mass_2=mass_2, chi_1=0.02, chi_2=0.02, luminosity_distance=50.,
- theta_jn=0.4, psi=2.659, phase=1.3, geocent_time=1126259642.413,
- ra=1.375, dec=-1.2108, lambda_1=lambda_1, lambda_2=lambda_2)
+ mass_1=mass_1,
+ mass_2=mass_2,
+ chi_1=0.02,
+ chi_2=0.02,
+ luminosity_distance=50.0,
+ theta_jn=0.4,
+ psi=2.659,
+ phase=1.3,
+ geocent_time=1126259642.413,
+ ra=1.375,
+ dec=-1.2108,
+ lambda_1=lambda_1,
+ lambda_2=lambda_2,
+)
# Set the duration and sampling frequency of the data segment that we're going
# to inject the signal into. For the
# TaylorF2 waveform, we cut the signal close to the isco frequency
duration = 32
-sampling_frequency = 2 * 1024
-start_time = injection_parameters['geocent_time'] + 2 - duration
+sampling_frequency = 2048
+start_time = injection_parameters["geocent_time"] + 2 - duration
# Fixed arguments passed into the source model. The analysis starts at 40 Hz.
# Note that the EoS sampling is agnostic to waveform model as long as the approximant
# can include tides.
-waveform_arguments = dict(waveform_approximant='IMRPhenomPv2_NRTidal',
- reference_frequency=50., minimum_frequency=40.0)
+waveform_arguments = dict(
+ waveform_approximant="IMRPhenomPv2_NRTidal",
+ reference_frequency=50.0,
+ minimum_frequency=40.0,
+)
# Create the waveform_generator using a LAL Binary Neutron Star source function
waveform_generator = bilby.gw.WaveformGenerator(
- duration=duration, sampling_frequency=sampling_frequency,
+ duration=duration,
+ sampling_frequency=sampling_frequency,
frequency_domain_source_model=bilby.gw.source.lal_binary_neutron_star,
parameter_conversion=bilby.gw.conversion.convert_to_lal_binary_neutron_star_parameters,
- waveform_arguments=waveform_arguments)
+ waveform_arguments=waveform_arguments,
+)
# Set up interferometers. In this case we'll use three interferometers
# (LIGO-Hanford (H1), LIGO-Livingston (L1), and Virgo (V1)).
# These default to their design sensitivity and start at 40 Hz.
-interferometers = bilby.gw.detector.InterferometerList(['H1', 'L1', 'V1'])
+interferometers = bilby.gw.detector.InterferometerList(["H1", "L1", "V1"])
for interferometer in interferometers:
interferometer.minimum_frequency = 40
interferometers.set_strain_data_from_power_spectral_densities(
- sampling_frequency=sampling_frequency, duration=duration,
- start_time=start_time)
-interferometers.inject_signal(parameters=injection_parameters,
- waveform_generator=waveform_generator)
+ sampling_frequency=sampling_frequency, duration=duration, start_time=start_time
+)
+interferometers.inject_signal(
+ parameters=injection_parameters, waveform_generator=waveform_generator
+)
# We're going to sample in chirp_mass, symmetric_mass_ratio, and
# specific EoS model parameters. We're using a 4-parameter
@@ -84,38 +99,62 @@ interferometers.inject_signal(parameters=injection_parameters,
# BNS have aligned spins by default, if you want to allow precessing spins
# pass aligned_spin=False to the BNSPriorDict
priors = bilby.gw.prior.BNSPriorDict()
-for key in ['psi', 'geocent_time', 'ra', 'dec', 'chi_1', 'chi_2',
- 'theta_jn', 'luminosity_distance', 'phase']:
+for key in [
+ "psi",
+ "geocent_time",
+ "ra",
+ "dec",
+ "chi_1",
+ "chi_2",
+ "theta_jn",
+ "luminosity_distance",
+ "phase",
+]:
priors[key] = injection_parameters[key]
-priors.pop('mass_1')
-priors.pop('mass_2')
-priors.pop('lambda_1')
-priors.pop('lambda_2')
-priors.pop('mass_ratio')
-priors['chirp_mass'] = bilby.core.prior.Gaussian(1.215, 0.1, name='chirp_mass', unit='$M_{\\odot}$')
-priors['symmetric_mass_ratio'] = bilby.core.prior.Uniform(0.1, 0.25, name='symmetric_mass_ratio')
-priors['eos_spectral_gamma_0'] = bilby.core.prior.Uniform(0.2, 2.0, name='gamma0', latex_label='$\\gamma_0')
-priors['eos_spectral_gamma_1'] = bilby.core.prior.Uniform(-1.6, 1.7, name='gamma1', latex_label='$\\gamma_1')
-priors['eos_spectral_gamma_2'] = bilby.core.prior.Uniform(-0.6, 0.6, name='gamma2', latex_label='$\\gamma_2')
-priors['eos_spectral_gamma_3'] = bilby.core.prior.Uniform(-0.02, 0.02, name='gamma3', latex_label='$\\gamma_3')
+for key in ["mass_1", "mass_2", "lambda_1", "lambda_2", "mass_ratio"]:
+ del priors[key]
+priors["chirp_mass"] = bilby.core.prior.Gaussian(
+ 1.215, 0.1, name="chirp_mass", unit="$M_{\\odot}$"
+)
+priors["symmetric_mass_ratio"] = bilby.core.prior.Uniform(
+ 0.1, 0.25, name="symmetric_mass_ratio"
+)
+priors["eos_spectral_gamma_0"] = bilby.core.prior.Uniform(
+ 0.2, 2.0, name="gamma0", latex_label="$\\gamma_0"
+)
+priors["eos_spectral_gamma_1"] = bilby.core.prior.Uniform(
+ -1.6, 1.7, name="gamma1", latex_label="$\\gamma_1"
+)
+priors["eos_spectral_gamma_2"] = bilby.core.prior.Uniform(
+ -0.6, 0.6, name="gamma2", latex_label="$\\gamma_2"
+)
+priors["eos_spectral_gamma_3"] = bilby.core.prior.Uniform(
+ -0.02, 0.02, name="gamma3", latex_label="$\\gamma_3"
+)
# The eos_check prior imposes several hard physical constraints on samples like
# enforcing causality and monotinicity of the EoSs. In almost ever conceivable
# sampling scenario, this should be enabled.
-priors['eos_check'] = bilby.gw.prior.EOSCheck()
+priors["eos_check"] = bilby.gw.prior.EOSCheck()
# Initialise the likelihood by passing in the interferometer data (IFOs)
# and the waveform generator
likelihood = bilby.gw.GravitationalWaveTransient(
- interferometers=interferometers, waveform_generator=waveform_generator,
- time_marginalization=False, phase_marginalization=False,
- distance_marginalization=False, priors=priors)
+ interferometers=interferometers,
+ waveform_generator=waveform_generator,
+)
# Run sampler. In this case we're going to use the `dynesty` sampler
result = bilby.run_sampler(
- likelihood=likelihood, priors=priors, sampler='dynesty', npoints=1000,
- injection_parameters=injection_parameters, outdir=outdir, label=label,
+ likelihood=likelihood,
+ priors=priors,
+ sampler="dynesty",
+ npoints=1000,
+ injection_parameters=injection_parameters,
+ outdir=outdir,
+ label=label,
conversion_function=bilby.gw.conversion.generate_all_bns_parameters,
- resume=True)
+ resume=True,
+)
result.plot_corner()
diff --git a/examples/gw_examples/injection_examples/calibration_example.py b/examples/gw_examples/injection_examples/calibration_example.py
index 8b9556f48ca43a2419629cc1d4f8fdb1daa25b96..cd5c4da87b9ef79acc53af50167a67c1aa343c19 100644
--- a/examples/gw_examples/injection_examples/calibration_example.py
+++ b/examples/gw_examples/injection_examples/calibration_example.py
@@ -2,20 +2,23 @@
"""
Tutorial to demonstrate running parameter estimation with calibration
uncertainties included.
+
+We set up the full problem as is required and then just sample over a small
+number of calibration parameters.
"""
-import numpy as np
import bilby
+import numpy as np
# Set the duration and sampling frequency of the data segment
# that we're going to create and inject the signal into.
-duration = 4.
-sampling_frequency = 2048.
+duration = 4
+sampling_frequency = 1024
# Specify the output directory and the name of the simulation.
-outdir = 'outdir'
-label = 'calibration'
+outdir = "outdir"
+label = "calibration"
bilby.core.utils.setup_logger(outdir=outdir, label=label)
# Set up a random seed for result reproducibility. This is optional!
@@ -26,37 +29,58 @@ np.random.seed(88170235)
# 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., theta_jn=0.4, psi=2.659,
- phase=1.3, geocent_time=1126259642.413, ra=1.375, dec=-1.2108)
+ mass_1=36.0,
+ mass_2=29.0,
+ 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.0,
+ 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='IMRPhenomPv2',
- reference_frequency=50.)
+waveform_arguments = dict(waveform_approximant="IMRPhenomXP", reference_frequency=50.0)
# Create the waveform_generator using a LAL BinaryBlackHole source function
waveform_generator = bilby.gw.WaveformGenerator(
- duration=duration, sampling_frequency=sampling_frequency,
+ duration=duration,
+ sampling_frequency=sampling_frequency,
frequency_domain_source_model=bilby.gw.source.lal_binary_black_hole,
- parameters=injection_parameters, waveform_arguments=waveform_arguments)
+ parameters=injection_parameters,
+ waveform_arguments=waveform_arguments,
+)
# Set up interferometers. In this case we'll use three interferometers
# (LIGO-Hanford (H1), LIGO-Livingston (L1), and Virgo (V1)).
# These default to their design sensitivity
-ifos = bilby.gw.detector.InterferometerList(['H1', 'L1', 'V1'])
+ifos = bilby.gw.detector.InterferometerList(["H1", "L1", "V1"])
for ifo in ifos:
- injection_parameters.update({
- 'recalib_{}_amplitude_{}'.format(ifo.name, ii): 0.1 for ii in range(5)})
- injection_parameters.update({
- 'recalib_{}_phase_{}'.format(ifo.name, ii): 0.01 for ii in range(5)})
+ injection_parameters.update(
+ {f"recalib_{ifo.name}_amplitude_{ii}": 0.1 for ii in range(5)}
+ )
+ injection_parameters.update(
+ {f"recalib_{ifo.name}_phase_{ii}": 0.01 for ii in range(5)}
+ )
ifo.calibration_model = bilby.gw.calibration.CubicSpline(
- prefix='recalib_{}_'.format(ifo.name),
+ prefix=f"recalib_{ifo.name}_",
minimum_frequency=ifo.minimum_frequency,
- maximum_frequency=ifo.maximum_frequency, n_points=5)
+ maximum_frequency=ifo.maximum_frequency,
+ n_points=5,
+ )
ifos.set_strain_data_from_power_spectral_densities(
- sampling_frequency=sampling_frequency, duration=duration)
-ifos.inject_signal(parameters=injection_parameters,
- waveform_generator=waveform_generator)
+ sampling_frequency=sampling_frequency, duration=duration
+)
+ifos.inject_signal(
+ parameters=injection_parameters, waveform_generator=waveform_generator
+)
# Set up prior, which is a dictionary
# Here we fix the injected cbc parameters and most of the calibration parameters
@@ -64,21 +88,30 @@ ifos.inject_signal(parameters=injection_parameters,
# We allow a subset of the calibration parameters to vary.
priors = injection_parameters.copy()
for key in injection_parameters:
- if 'recalib' in key:
+ if "recalib" in key:
priors[key] = injection_parameters[key]
-for name in ['recalib_H1_amplitude_0', 'recalib_H1_amplitude_1']:
- priors[name] = bilby.prior.Gaussian(
- mu=0, sigma=0.2, name=name, latex_label='H1 $A_{}$'.format(name[-1]))
+for name in ["recalib_H1_amplitude_0", "recalib_H1_amplitude_1"]:
+ priors[name] = bilby.core.prior.Gaussian(
+ mu=0, sigma=0.2, name=name, latex_label=f"H1 $A_{name[-1]}$"
+ )
# Initialise the likelihood by passing in the interferometer data (IFOs) and
# the waveform generator
likelihood = bilby.gw.GravitationalWaveTransient(
- interferometers=ifos, waveform_generator=waveform_generator)
+ interferometers=ifos, waveform_generator=waveform_generator
+)
# Run sampler. In this case we're going to use the `dynesty` sampler
result = bilby.run_sampler(
- likelihood=likelihood, priors=priors, sampler='dynesty', npoints=1000,
- injection_parameters=injection_parameters, outdir=outdir, label=label)
+ likelihood=likelihood,
+ priors=priors,
+ sampler="dynesty",
+ npoints=1000,
+ sample="unif",
+ injection_parameters=injection_parameters,
+ outdir=outdir,
+ label=label,
+)
# make some plots of the outputs
result.plot_corner()
diff --git a/examples/gw_examples/injection_examples/calibration_marginalization_example.py b/examples/gw_examples/injection_examples/calibration_marginalization_example.py
index 73279aadbf4435d59d51ce08b33d59569e475028..2f8dc326a04299df5dba5168097f67dcf04723fb 100644
--- a/examples/gw_examples/injection_examples/calibration_marginalization_example.py
+++ b/examples/gw_examples/injection_examples/calibration_marginalization_example.py
@@ -4,19 +4,20 @@ Tutorial to demonstrate running parameter estimation with calibration
uncertainties marginalized over using a finite set of realizations.
"""
-import numpy as np
+from copy import deepcopy
+
import bilby
-from copy import copy
import matplotlib.pyplot as plt
+import numpy as np
# Set the duration and sampling frequency of the data segment
# that we're going to create and inject the signal into.
-duration = 4.
-sampling_frequency = 2048.
+duration = 4
+sampling_frequency = 1024
# Specify the output directory and the name of the simulation.
-outdir = 'outdir'
-label = 'calibration_marginalization'
+outdir = "outdir"
+label = "calibration_marginalization"
bilby.core.utils.setup_logger(outdir=outdir, label=label)
# Set up a random seed for result reproducibility. This is optional!
@@ -27,110 +28,169 @@ np.random.seed(170817)
# 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., theta_jn=0.4, psi=2.659,
- phase=1.3, geocent_time=1126259642.413, ra=1.375, dec=-1.2108)
-start_time = injection_parameters['geocent_time'] - duration + 2
+ mass_1=36.0,
+ mass_2=29.0,
+ 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.0,
+ theta_jn=0.4,
+ psi=2.659,
+ phase=1.3,
+ geocent_time=1126259642.413,
+ ra=1.375,
+ dec=-1.2108,
+)
+start_time = injection_parameters["geocent_time"] - duration + 2
# Fixed arguments passed into the source model
-waveform_arguments = dict(waveform_approximant='IMRPhenomPv2',
- reference_frequency=50.)
+waveform_arguments = dict(waveform_approximant="IMRPhenomXP", reference_frequency=50.0)
# Create the waveform_generator using a LAL BinaryBlackHole source function
waveform_generator = bilby.gw.WaveformGenerator(
- duration=duration, sampling_frequency=sampling_frequency,
+ duration=duration,
+ sampling_frequency=sampling_frequency,
frequency_domain_source_model=bilby.gw.source.lal_binary_black_hole,
- parameters=injection_parameters, waveform_arguments=waveform_arguments)
-waveform_generator_rew = copy(waveform_generator)
+ parameters=injection_parameters,
+ waveform_arguments=waveform_arguments,
+)
+waveform_generator_rew = deepcopy(waveform_generator)
# Set up interferometers. In this case we'll use three interferometers
# (LIGO-Hanford (H1), LIGO-Livingston (L1), and Virgo (V1)).
# These default to their design sensitivity
-ifos = bilby.gw.detector.InterferometerList(['H1', 'L1', 'V1'])
+ifos = bilby.gw.detector.InterferometerList(["H1", "L1", "V1"])
for ifo in ifos:
- injection_parameters.update({
- 'recalib_{}_amplitude_{}'.format(ifo.name, ii): 0.0 for ii in range(10)})
- injection_parameters.update({
- 'recalib_{}_phase_{}'.format(ifo.name, ii): 0.0 for ii in range(10)})
+ injection_parameters.update(
+ {f"recalib_{ifo.name}_amplitude_{ii}": 0.0 for ii in range(10)}
+ )
+ injection_parameters.update(
+ {f"recalib_{ifo.name}_phase_{ii}": 0.0 for ii in range(10)}
+ )
ifo.calibration_model = bilby.gw.calibration.CubicSpline(
- prefix='recalib_{}_'.format(ifo.name),
+ prefix=f"recalib_{ifo.name}_",
minimum_frequency=ifo.minimum_frequency,
- maximum_frequency=ifo.maximum_frequency, n_points=10)
+ maximum_frequency=ifo.maximum_frequency,
+ n_points=10,
+ )
ifos.set_strain_data_from_power_spectral_densities(
- sampling_frequency=sampling_frequency, duration=duration, start_time=start_time)
-ifos.inject_signal(parameters=injection_parameters,
- waveform_generator=waveform_generator)
-ifos_rew = copy(ifos)
+ sampling_frequency=sampling_frequency, duration=duration, start_time=start_time
+)
+ifos.inject_signal(
+ parameters=injection_parameters, waveform_generator=waveform_generator
+)
+ifos_rew = deepcopy(ifos)
# Set up prior, which is a dictionary
# Here we fix the injected cbc parameters (except the distance)
# to the injected values.
priors = injection_parameters.copy()
-priors['luminosity_distance'] = bilby.prior.Uniform(
- injection_parameters['luminosity_distance'] - 1000, injection_parameters['luminosity_distance'] + 1000,
- name='luminosity_distance', latex_label='$d_L$')
+priors["luminosity_distance"] = bilby.prior.Uniform(
+ injection_parameters["luminosity_distance"] - 1000,
+ injection_parameters["luminosity_distance"] + 1000,
+ name="luminosity_distance",
+ latex_label="$d_L$",
+)
for key in injection_parameters:
- if 'recalib' in key:
+ if "recalib" in key:
priors[key] = injection_parameters[key]
# Convert to prior dictionary to replace the floats with delta function priors
priors = bilby.core.prior.PriorDict(priors)
-priors_rew = copy(priors)
+priors_rew = deepcopy(priors)
# Initialise the likelihood by passing in the interferometer data (IFOs) and
# the waveform generator. Here we assume there is no calibration uncertainty
likelihood = bilby.gw.GravitationalWaveTransient(
- interferometers=ifos, waveform_generator=waveform_generator, priors=priors)
+ interferometers=ifos, waveform_generator=waveform_generator, priors=priors
+)
# Run sampler. In this case we're going to use the `dynesty` sampler
result = bilby.run_sampler(
- likelihood=likelihood, priors=priors, sampler='dynesty', npoints=500,
- injection_parameters=injection_parameters, outdir=outdir, label=label)
+ likelihood=likelihood,
+ priors=priors,
+ sampler="dynesty",
+ npoints=500,
+ walks=20,
+ nact=3,
+ injection_parameters=injection_parameters,
+ outdir=outdir,
+ label=label,
+)
# Setting the log likelihood to actually be the log likelihood and not the log likelihood ratio...
# This is used the for reweighting
-result.posterior['log_likelihood'] = result.posterior['log_likelihood'] + result.log_noise_evidence
+result.posterior["log_likelihood"] = (
+ result.posterior["log_likelihood"] + result.log_noise_evidence
+)
# Setting the priors we want on the calibration response curve parameters - as an example.
-for name in ['recalib_H1_amplitude_1', 'recalib_H1_amplitude_4']:
+for name in ["recalib_H1_amplitude_1", "recalib_H1_amplitude_4"]:
priors_rew[name] = bilby.prior.Gaussian(
- mu=0, sigma=0.03, name=name, latex_label='H1 $A_{}$'.format(name[-1]))
+ mu=0, sigma=0.03, name=name, latex_label=f"H1 $A_{name[-1]}$"
+ )
# Setting up the calibration marginalized likelihood.
# We save the calibration response curve files into the output directory under {ifo.name}_calibration_file.h5
cal_likelihood = bilby.gw.GravitationalWaveTransient(
- interferometers=ifos_rew, waveform_generator=waveform_generator_rew,
- calibration_marginalization=True, priors=priors_rew,
+ interferometers=ifos_rew,
+ waveform_generator=waveform_generator_rew,
+ calibration_marginalization=True,
+ priors=priors_rew,
number_of_response_curves=100,
- calibration_lookup_table={ifos[i].name: f'{outdir}/{ifos[i].name}_calibration_file.h5' for i in range(len(ifos))})
+ calibration_lookup_table={
+ ifos[i].name: f"{outdir}/{ifos[i].name}_calibration_file.h5"
+ for i in range(len(ifos))
+ },
+)
# Plot the magnitude of the curves to be used in the marginalization
-plt.semilogx(ifos[0].frequency_array[ifos[0].frequency_mask][0:-1],
- cal_likelihood.calibration_draws[ifos[0].name][:, 0:-1].T)
+plt.semilogx(
+ ifos[0].frequency_array[ifos[0].frequency_mask][0:-1],
+ cal_likelihood.calibration_draws[ifos[0].name][:, 0:-1].T,
+)
plt.xlim(20, 1024)
-plt.ylabel('Magnitude')
-plt.xlabel('Frequency [Hz]')
+plt.ylabel("Magnitude")
+plt.xlabel("Frequency [Hz]")
plt.savefig(f"{outdir}/calibration_draws.pdf")
plt.clf()
# Reweight the posterior samples from a distribution with no calibration uncertainty to one with uncertainty.
# This method utilizes rejection sampling which can be inefficient at drawing samples at higher SNRs.
-result_rew = bilby.core.result.reweight(result,
- new_likelihood=cal_likelihood,
- conversion_function=bilby.gw.conversion.generate_all_bbh_parameters)
+result_rew = bilby.core.result.reweight(
+ result,
+ new_likelihood=cal_likelihood,
+ conversion_function=bilby.gw.conversion.generate_all_bbh_parameters,
+)
# Plot distance posterior with and without the calibration
-for res, label in zip([result, result_rew], ["No calibration uncertainty", "Calibration uncertainty"]):
- plt.hist(res.posterior['luminosity_distance'], label=label, bins=50, histtype='step', density=True)
+for res, label in zip(
+ [result, result_rew], ["No calibration uncertainty", "Calibration uncertainty"]
+):
+ plt.hist(
+ res.posterior["luminosity_distance"],
+ label=label,
+ bins=50,
+ histtype="step",
+ density=True,
+ )
plt.legend()
-plt.xlabel('Luminosity distance [Mpc]')
+plt.xlabel("Luminosity distance [Mpc]")
plt.savefig(f"{outdir}/luminosity_distance_posterior.pdf")
plt.clf()
plt.hist(
- result_rew.posterior['recalib_index'],
- bins=np.linspace(0, cal_likelihood.number_of_response_curves - 1, cal_likelihood.number_of_response_curves),
- density=True)
+ result_rew.posterior["recalib_index"],
+ bins=np.linspace(
+ 0,
+ cal_likelihood.number_of_response_curves - 1,
+ cal_likelihood.number_of_response_curves,
+ ),
+ density=True,
+)
plt.xlim(0, cal_likelihood.number_of_response_curves - 1)
-plt.xlabel('Calibration index')
+plt.xlabel("Calibration index")
plt.savefig(f"{outdir}/calibration_index_histogram.pdf")
diff --git a/examples/gw_examples/injection_examples/change_sampled_parameters.py b/examples/gw_examples/injection_examples/change_sampled_parameters.py
index 4e0266d08b6290395d2a5e8dc47bdf0dc6ca9432..10eedb7a6d0c19c3e9fae60a214bd6093bf170ef 100644
--- a/examples/gw_examples/injection_examples/change_sampled_parameters.py
+++ b/examples/gw_examples/injection_examples/change_sampled_parameters.py
@@ -10,71 +10,107 @@ The cosmology is according to the Planck 2015 data release.
import bilby
import numpy as np
-
bilby.core.utils.setup_logger(log_level="info")
-duration = 4.
-sampling_frequency = 2048.
-outdir = 'outdir'
+duration = 4
+sampling_frequency = 2048
+outdir = "outdir"
+label = "different_parameters"
np.random.seed(151226)
injection_parameters = dict(
- total_mass=66., mass_ratio=0.9, 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, theta_jn=0.4, psi=2.659,
- phase=1.3, geocent_time=1126259642.413, ra=1.375, dec=-1.2108)
+ total_mass=66.0,
+ mass_ratio=0.9,
+ 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,
+ theta_jn=0.4,
+ psi=2.659,
+ phase=1.3,
+ geocent_time=1126259642.413,
+ ra=1.375,
+ dec=-1.2108,
+)
-waveform_arguments = dict(waveform_approximant='IMRPhenomPv2',
- reference_frequency=50.)
+waveform_arguments = dict(waveform_approximant="IMRPhenomXP", reference_frequency=50.0)
# Create the waveform_generator using a LAL BinaryBlackHole source function
# We specify a function which transforms a dictionary of parameters into the
# appropriate parameters for the source model.
waveform_generator = bilby.gw.waveform_generator.WaveformGenerator(
- sampling_frequency=sampling_frequency, duration=duration,
+ sampling_frequency=sampling_frequency,
+ duration=duration,
frequency_domain_source_model=bilby.gw.source.lal_binary_black_hole,
parameter_conversion=bilby.gw.conversion.convert_to_lal_binary_black_hole_parameters,
- waveform_arguments=waveform_arguments)
+ waveform_arguments=waveform_arguments,
+)
# Set up interferometers.
-ifos = bilby.gw.detector.InterferometerList(['H1', 'L1', 'V1', 'K1'])
+ifos = bilby.gw.detector.InterferometerList(["H1", "L1", "V1", "K1"])
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)
+ sampling_frequency=sampling_frequency,
+ duration=duration,
+ start_time=injection_parameters["geocent_time"] - 2,
+)
+ifos.inject_signal(
+ waveform_generator=waveform_generator, parameters=injection_parameters
+)
# Set up prior
# Note it is possible to sample in different parameters to those that were
# injected.
priors = bilby.gw.prior.BBHPriorDict()
-priors.pop('mass_1')
-priors.pop('mass_2')
-priors.pop('luminosity_distance')
-priors['chirp_mass'] = bilby.prior.Uniform(
- name='chirp_mass', latex_label='$m_c$', minimum=13, maximum=45,
- unit='$M_{\\odot}$')
-priors['symmetric_mass_ratio'] = bilby.prior.Uniform(
- name='symmetric_mass_ratio', latex_label='q', minimum=0.1, maximum=0.25)
-priors['redshift'] = bilby.prior.Uniform(
- name='redshift', latex_label='$z$', minimum=0, maximum=0.5)
+
+del priors["mass_ratio"]
+priors["symmetric_mass_ratio"] = bilby.prior.Uniform(
+ name="symmetric_mass_ratio", latex_label="q", minimum=0.1, maximum=0.25
+)
+
+del priors["luminosity_distance"]
+priors["redshift"] = bilby.prior.Uniform(
+ name="redshift", latex_label="$z$", minimum=0, maximum=0.5
+)
# These parameters will not be sampled
-for key in ['a_1', 'a_2', 'tilt_1', 'tilt_2', 'phi_12', 'phi_jl', 'psi',
- 'ra', 'dec', 'geocent_time', 'phase']:
+for key in [
+ "a_1",
+ "a_2",
+ "tilt_1",
+ "tilt_2",
+ "phi_12",
+ "phi_jl",
+ "psi",
+ "ra",
+ "dec",
+ "geocent_time",
+ "phase",
+]:
priors[key] = injection_parameters[key]
-priors.pop('theta_jn')
-priors['cos_theta_jn'] = np.cos(injection_parameters['theta_jn'])
+del priors["theta_jn"]
+priors["cos_theta_jn"] = np.cos(injection_parameters["theta_jn"])
print(priors)
# Initialise GravitationalWaveTransient
likelihood = bilby.gw.likelihood.GravitationalWaveTransient(
- interferometers=ifos, waveform_generator=waveform_generator)
+ interferometers=ifos, waveform_generator=waveform_generator
+)
# Run sampler
# Note we've added a post-processing conversion function, this will generate
# many useful additional parameters, e.g., source-frame masses.
result = bilby.core.sampler.run_sampler(
- likelihood=likelihood, priors=priors, sampler='dynesty', outdir=outdir,
- injection_parameters=injection_parameters, label='DifferentParameters',
- conversion_function=bilby.gw.conversion.generate_all_bbh_parameters)
+ likelihood=likelihood,
+ priors=priors,
+ sampler="dynesty",
+ walks=25,
+ nact=5,
+ outdir=outdir,
+ injection_parameters=injection_parameters,
+ label=label,
+ conversion_function=bilby.gw.conversion.generate_all_bbh_parameters,
+)
result.plot_corner()
diff --git a/examples/gw_examples/injection_examples/create_your_own_source_model.py b/examples/gw_examples/injection_examples/create_your_own_source_model.py
index f755363f091a12f11df425ccaea3d7d7831090c2..77f576a0bd80cc38aaaef1050c04ac4a63fb7eb4 100644
--- a/examples/gw_examples/injection_examples/create_your_own_source_model.py
+++ b/examples/gw_examples/injection_examples/create_your_own_source_model.py
@@ -6,48 +6,97 @@ import bilby
import numpy as np
# First set up logging and some output directories and labels
-outdir = 'outdir'
-label = 'create_your_own_source_model'
+outdir = "outdir"
+label = "create_your_own_source_model"
sampling_frequency = 4096
duration = 1
# Here we define out source model - this is the sine-Gaussian model in the
# frequency domain.
-def sine_gaussian(f, A, f0, tau, phi0, geocent_time, ra, dec, psi):
- arg = -(np.pi * tau * (f - f0))**2 + 1j * phi0
- plus = np.sqrt(np.pi) * A * tau * np.exp(arg) / 2.
+def gaussian(frequency_array, amplitude, f0, tau, phi0):
+ r"""
+ Our custom source model, this is just a Gaussian in frequency with
+ variable global phase.
+
+ .. math::
+
+ \tilde{h}_{\plus}(f) = \frac{A \tau}{2\sqrt{\pi}}}
+ e^{- \pi \tau (f - f_{0})^2 + i \phi_{0}} \\
+ \tilde{h}_{\times}(f) = \tilde{h}_{\plus}(f) e^{i \pi / 2}
+
+
+ Parameters
+ ----------
+ frequency_array: array-like
+ The frequencies to evaluate the model at. This is required for all
+ Bilby source models.
+ amplitude: float
+ An overall amplitude prefactor.
+ f0: float
+ The central frequency.
+ tau: float
+ The damping rate.
+ phi0: float
+ The reference phase.
+
+ Returns
+ -------
+ dict:
+ A dictionary containing "plus" and "cross" entries.
+
+ """
+ arg = -((np.pi * tau * (frequency_array - f0)) ** 2) + 1j * phi0
+ plus = np.sqrt(np.pi) * amplitude * tau * np.exp(arg) / 2.0
cross = plus * np.exp(1j * np.pi / 2)
- return {'plus': plus, 'cross': cross}
+ return {"plus": plus, "cross": cross}
# We now define some parameters that we will inject
-injection_parameters = dict(A=1e-23, f0=100, tau=1, phi0=0, geocent_time=0,
- ra=0, dec=0, psi=0)
+injection_parameters = dict(
+ amplitude=1e-23, f0=100, tau=1, phi0=0, geocent_time=0, ra=0, dec=0, psi=0
+)
# Now we pass our source function to the WaveformGenerator
waveform_generator = bilby.gw.waveform_generator.WaveformGenerator(
- duration=duration, sampling_frequency=sampling_frequency,
- frequency_domain_source_model=sine_gaussian)
+ duration=duration,
+ sampling_frequency=sampling_frequency,
+ frequency_domain_source_model=gaussian,
+)
# Set up interferometers.
-ifos = bilby.gw.detector.InterferometerList(['H1', 'L1'])
+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)
+ sampling_frequency=sampling_frequency,
+ duration=duration,
+ start_time=injection_parameters["geocent_time"] - 0.5,
+)
+ifos.inject_signal(
+ waveform_generator=waveform_generator,
+ parameters=injection_parameters,
+ raise_error=False,
+)
# Here we define the priors for the search. We use the injection parameters
# except for the amplitude, f0, and geocent_time
prior = injection_parameters.copy()
-prior['A'] = bilby.core.prior.LogUniform(minimum=1e-25, maximum=1e-21, name='A')
-prior['f0'] = bilby.core.prior.Uniform(90, 110, 'f')
+prior["amplitude"] = bilby.core.prior.LogUniform(
+ minimum=1e-25, maximum=1e-21, latex_label="$\\mathcal{A}$"
+)
+prior["f0"] = bilby.core.prior.Uniform(90, 110, latex_label="$f_{0}$")
likelihood = bilby.gw.likelihood.GravitationalWaveTransient(
- interferometers=ifos, waveform_generator=waveform_generator)
+ interferometers=ifos, waveform_generator=waveform_generator
+)
result = bilby.core.sampler.run_sampler(
- likelihood, prior, sampler='dynesty', outdir=outdir, label=label,
- resume=False, sample='unif', injection_parameters=injection_parameters)
+ likelihood,
+ prior,
+ sampler="dynesty",
+ outdir=outdir,
+ label=label,
+ resume=False,
+ sample="unif",
+ injection_parameters=injection_parameters,
+)
result.plot_corner()
diff --git a/examples/gw_examples/injection_examples/create_your_own_time_domain_source_model.py b/examples/gw_examples/injection_examples/create_your_own_time_domain_source_model.py
index 11290fde30e52f9d1851654d3ffa4349f697fbe8..1427f8ceb2813a262db031a9e494a0ce580c012f 100644
--- a/examples/gw_examples/injection_examples/create_your_own_time_domain_source_model.py
+++ b/examples/gw_examples/injection_examples/create_your_own_time_domain_source_model.py
@@ -6,62 +6,114 @@ noise in two interferometers (LIGO Livingston and Hanford at design
sensitivity), and then recovered.
"""
-import numpy as np
import bilby
+import numpy as np
# define the time-domain model
-def time_domain_damped_sinusoid(
- time, amplitude, damping_time, frequency, phase, t0):
- """
+def time_domain_damped_sinusoid(time, amplitude, damping_time, frequency, phase, t0):
+ r"""
This example only creates a linearly polarised signal with only plus
polarisation.
+
+ .. math::
+
+ h_{\plus}(t) =
+ \Theta(t - t_{0}) A
+ e^{-(t - t_{0}) / \tau}
+ \sin \left( 2 \pi f t + \phi \right)
+
+ Parameters
+ ----------
+ time: array-like
+ The times at which to evaluate the model. This is required for all
+ time-domain models.
+ amplitude: float
+ The peak amplitude.
+ damping_time: float
+ The damping time of the exponential.
+ frequency: float
+ The frequency of the oscillations.
+ phase: float
+ The initial phase of the signal.
+ t0: float
+ The offset of the start of the signal from the start time.
+
+ Returns
+ -------
+ dict:
+ A dictionary containing "plus" and "cross" entries.
+
"""
plus = np.zeros(len(time))
tidx = time >= t0
- plus[tidx] = amplitude * np.exp(-(time[tidx] - t0) / damping_time) *\
- np.sin(2 * np.pi * frequency * (time[tidx] - t0) + phase)
+ plus[tidx] = (
+ amplitude
+ * np.exp(-(time[tidx] - t0) / damping_time)
+ * np.sin(2 * np.pi * frequency * (time[tidx] - t0) + phase)
+ )
cross = np.zeros(len(time))
- return {'plus': plus, 'cross': cross}
+ return {"plus": plus, "cross": cross}
# define parameters to inject.
-injection_parameters = dict(amplitude=5e-22, damping_time=0.1, frequency=50,
- phase=0, ra=0, dec=0, psi=0, t0=0., geocent_time=0.)
+injection_parameters = dict(
+ amplitude=5e-22,
+ damping_time=0.1,
+ frequency=50,
+ phase=0,
+ ra=0,
+ dec=0,
+ psi=0,
+ t0=0.0,
+ geocent_time=0.0,
+)
-duration = 1.0
+duration = 1
sampling_frequency = 1024
-outdir = 'outdir'
-label = 'time_domain_source_model'
+outdir = "outdir"
+label = "time_domain_source_model"
# call the waveform_generator to create our waveform model.
waveform = bilby.gw.waveform_generator.WaveformGenerator(
- duration=duration, sampling_frequency=sampling_frequency,
+ duration=duration,
+ sampling_frequency=sampling_frequency,
time_domain_source_model=time_domain_damped_sinusoid,
- start_time=injection_parameters['geocent_time'] - 0.5)
+ start_time=injection_parameters["geocent_time"] - 0.5,
+)
# inject the signal into three interferometers
-ifos = bilby.gw.detector.InterferometerList(['H1', 'L1'])
+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'] - 0.5)
-ifos.inject_signal(waveform_generator=waveform,
- parameters=injection_parameters)
+ sampling_frequency=sampling_frequency,
+ duration=duration,
+ start_time=injection_parameters["geocent_time"] - 0.5,
+)
+ifos.inject_signal(
+ waveform_generator=waveform, parameters=injection_parameters, raise_error=False
+)
# create the priors
prior = injection_parameters.copy()
-prior['amplitude'] = bilby.core.prior.LogUniform(1e-23, 1e-21, r'$h_0$')
-prior['damping_time'] = bilby.core.prior.Uniform(
- 0.01, 1, r'damping time', unit='$s$')
-prior['frequency'] = bilby.core.prior.Uniform(0, 200, r'frequency', unit='Hz')
-prior['phase'] = bilby.core.prior.Uniform(-np.pi / 2, np.pi / 2, r'$\phi$')
+prior["amplitude"] = bilby.core.prior.LogUniform(1e-23, 1e-21, r"$h_0$")
+prior["damping_time"] = bilby.core.prior.Uniform(0.01, 1, r"damping time", unit="$s$")
+prior["frequency"] = bilby.core.prior.Uniform(0, 200, r"frequency", unit="Hz")
+prior["phase"] = bilby.core.prior.Uniform(-np.pi / 2, np.pi / 2, r"$\phi$")
# define likelihood
likelihood = bilby.gw.likelihood.GravitationalWaveTransient(ifos, waveform)
# launch sampler
result = bilby.core.sampler.run_sampler(
- likelihood, prior, sampler='dynesty', npoints=1000,
- injection_parameters=injection_parameters, outdir=outdir, label=label)
+ likelihood,
+ prior,
+ sampler="dynesty",
+ npoints=500,
+ walks=5,
+ nact=3,
+ injection_parameters=injection_parameters,
+ outdir=outdir,
+ label=label,
+)
result.plot_corner()
diff --git a/examples/gw_examples/injection_examples/custom_proposal_example.py b/examples/gw_examples/injection_examples/custom_proposal_example.py
index 2a19f4ee6944ff7aceff31591dad3c837ce0a73c..47019dc858da3d9cd640395a5174955bb40ac012 100644
--- a/examples/gw_examples/injection_examples/custom_proposal_example.py
+++ b/examples/gw_examples/injection_examples/custom_proposal_example.py
@@ -1,71 +1,105 @@
#!/usr/bin/env python
"""
Tutorial for running cpnest with custom jump proposals.
+
+This example takes longer than most to run.
+
+Due to how cpnest creates parallel processes, the multiprocessing start method
+needs to be set on some operating systems.
"""
+import multiprocessing
+
+multiprocessing.set_start_method("fork") # noqa
-import numpy as np
import bilby.gw.sampler.proposal
+import numpy as np
from bilby.core.sampler import proposal
-
# The set up here is the same as in fast_tutorial.py. Look there for descriptive explanations.
-duration = 4.
-sampling_frequency = 2048.
+duration = 4
+sampling_frequency = 1024
-outdir = 'outdir'
-label = 'custom_jump_proposals'
+outdir = "outdir"
+label = "custom_jump_proposals"
bilby.core.utils.setup_logger(outdir=outdir, label=label)
np.random.seed(88170235)
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., theta_jn=0.4, psi=2.659,
- phase=1.3, geocent_time=1126259642.413, ra=1.375, dec=-1.2108)
-waveform_arguments = dict(waveform_approximant='IMRPhenomPv2',
- reference_frequency=50., minimum_frequency=20.)
+ mass_1=36.0,
+ mass_2=29.0,
+ 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.0,
+ theta_jn=0.4,
+ psi=2.659,
+ phase=1.3,
+ geocent_time=1126259642.413,
+ ra=1.375,
+ dec=-1.2108,
+)
+waveform_arguments = dict(
+ waveform_approximant="IMRPhenomPv2",
+ reference_frequency=50.0,
+ minimum_frequency=20.0,
+)
waveform_generator = bilby.gw.WaveformGenerator(
- duration=duration, sampling_frequency=sampling_frequency,
+ duration=duration,
+ sampling_frequency=sampling_frequency,
frequency_domain_source_model=bilby.gw.source.lal_binary_black_hole,
- waveform_arguments=waveform_arguments)
-ifos = bilby.gw.detector.InterferometerList(['H1', 'L1'])
+ waveform_arguments=waveform_arguments,
+)
+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)
+ sampling_frequency=sampling_frequency,
+ duration=duration,
+ start_time=injection_parameters["geocent_time"] - 2,
+)
+ifos.inject_signal(
+ waveform_generator=waveform_generator, parameters=injection_parameters
+)
priors = bilby.gw.prior.BBHPriorDict()
-priors['geocent_time'] = bilby.core.prior.Uniform(
- minimum=injection_parameters['geocent_time'] - 1,
- maximum=injection_parameters['geocent_time'] + 1,
- name='geocent_time', latex_label='$t_c$', unit='$s$')
-for key in ['a_1', 'a_2', 'tilt_1', 'tilt_2', 'phi_12', 'phi_jl', 'geocent_time']:
+for key in ["a_1", "a_2", "tilt_1", "tilt_2", "phi_12", "phi_jl", "geocent_time"]:
priors[key] = injection_parameters[key]
likelihood = bilby.gw.GravitationalWaveTransient(
- interferometers=ifos, waveform_generator=waveform_generator)
+ interferometers=ifos, waveform_generator=waveform_generator
+)
# Definition of the custom jump proposals. Define a JumpProposalCycle. The first argument is a list
# of all allowed jump proposals. The second argument is a list of weights for the respective jump
# proposal.
jump_proposals = proposal.JumpProposalCycle(
- [proposal.EnsembleWalk(priors=priors),
- proposal.EnsembleStretch(priors=priors),
- proposal.DifferentialEvolution(priors=priors),
- proposal.EnsembleEigenVector(priors=priors),
- bilby.gw.sampler.proposal.SkyLocationWanderJump(priors=priors),
- bilby.gw.sampler.proposal.CorrelatedPolarisationPhaseJump(priors=priors),
- bilby.gw.sampler.proposal.PolarisationPhaseJump(priors=priors),
- proposal.DrawFlatPrior(priors=priors)],
- weights=[2, 2, 5, 1, 1, 1, 1, 1])
+ [
+ proposal.EnsembleWalk(priors=priors),
+ proposal.EnsembleStretch(priors=priors),
+ proposal.DifferentialEvolution(priors=priors),
+ proposal.EnsembleEigenVector(priors=priors),
+ bilby.gw.sampler.proposal.SkyLocationWanderJump(priors=priors),
+ bilby.gw.sampler.proposal.CorrelatedPolarisationPhaseJump(priors=priors),
+ bilby.gw.sampler.proposal.PolarisationPhaseJump(priors=priors),
+ proposal.DrawFlatPrior(priors=priors),
+ ],
+ weights=[2, 2, 5, 1, 1, 1, 1, 1],
+)
# Run cpnest with the proposals kwarg specified.
# Make sure to have a version of cpnest installed that supports custom proposals.
result = bilby.run_sampler(
- likelihood=likelihood, priors=priors, sampler='cpnest', npoints=1000,
- injection_parameters=injection_parameters, outdir=outdir, label=label,
- proposals=jump_proposals)
+ likelihood=likelihood,
+ priors=priors,
+ sampler="cpnest",
+ npoints=1000,
+ injection_parameters=injection_parameters,
+ outdir=outdir,
+ label=label,
+ proposals=jump_proposals,
+)
# Make a corner plot.
result.plot_corner()
diff --git a/examples/gw_examples/injection_examples/eccentric_inspiral.py b/examples/gw_examples/injection_examples/eccentric_inspiral.py
index d081cb575bc15855f1e05d06e5ec0426c62ac82d..55a5853ea4e7f135ef5bb37b7aa1e3c79e775cdd 100644
--- a/examples/gw_examples/injection_examples/eccentric_inspiral.py
+++ b/examples/gw_examples/injection_examples/eccentric_inspiral.py
@@ -6,84 +6,121 @@ similar to GW150914.
This uses the same binary parameters that were used to make Figures 1, 2 & 5 in
Lower et al. (2018) -> arXiv:1806.05350.
-
-For a more comprehensive look at what goes on in each step, refer to the
-"basic_tutorial.py" example.
"""
-import numpy as np
import bilby
+import numpy as np
-duration = 64.
-sampling_frequency = 256.
+duration = 64
+sampling_frequency = 256
-outdir = 'outdir'
-label = 'eccentric_GW140914'
+outdir = "outdir"
+label = "eccentric_GW150914"
bilby.core.utils.setup_logger(outdir=outdir, label=label)
# Set up a random seed for result reproducibility.
np.random.seed(150914)
injection_parameters = dict(
- mass_1=35., mass_2=30., eccentricity=0.1, luminosity_distance=440.,
- theta_jn=0.4, psi=0.1, phase=1.2, geocent_time=1180002601.0, ra=45, dec=5.73)
-
-waveform_arguments = dict(waveform_approximant='EccentricFD',
- reference_frequency=10., minimum_frequency=10.)
+ mass_1=35.0,
+ mass_2=30.0,
+ eccentricity=0.1,
+ luminosity_distance=440.0,
+ theta_jn=0.4,
+ psi=0.1,
+ phase=1.2,
+ geocent_time=1180002601.0,
+ ra=45,
+ dec=5.73,
+)
+
+waveform_arguments = dict(
+ waveform_approximant="EccentricFD", reference_frequency=10.0, minimum_frequency=10.0
+)
# Create the waveform_generator using the LAL eccentric black hole no spins
# source function
waveform_generator = bilby.gw.WaveformGenerator(
- duration=duration, sampling_frequency=sampling_frequency,
+ duration=duration,
+ sampling_frequency=sampling_frequency,
frequency_domain_source_model=bilby.gw.source.lal_eccentric_binary_black_hole_no_spins,
- parameters=injection_parameters, waveform_arguments=waveform_arguments)
+ parameters=injection_parameters,
+ waveform_arguments=waveform_arguments,
+)
# Setting up three interferometers (LIGO-Hanford (H1), LIGO-Livingston (L1), and
# Virgo (V1)) at their design sensitivities. The maximum frequency is set just
# prior to the point at which the waveform model terminates. This is to avoid
# any biases introduced from using a sharply terminating waveform model.
-minimum_frequency = 10.
-maximum_frequency = 128.
+minimum_frequency = 10.0
+maximum_frequency = 128.0
-ifos = bilby.gw.detector.InterferometerList(['H1', 'L1'])
+ifos = bilby.gw.detector.InterferometerList(["H1", "L1"])
for ifo in ifos:
ifo.minimum_frequency = minimum_frequency
ifo.maximum_frequency = maximum_frequency
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)
+ sampling_frequency=sampling_frequency,
+ duration=duration,
+ start_time=injection_parameters["geocent_time"] + 2 - duration,
+)
+ifos.inject_signal(
+ waveform_generator=waveform_generator, parameters=injection_parameters
+)
# Now we set up the priors on each of the binary parameters.
priors = bilby.core.prior.PriorDict()
priors["mass_1"] = bilby.core.prior.Uniform(
- name='mass_1', minimum=5, maximum=60, unit='$M_{\\odot}$')
+ name="mass_1", minimum=5, maximum=60, unit="$M_{\\odot}$", latex_label="$m_1$"
+)
priors["mass_2"] = bilby.core.prior.Uniform(
- name='mass_2', minimum=5, maximum=60, unit='$M_{\\odot}$')
+ name="mass_2", minimum=5, maximum=60, unit="$M_{\\odot}$", latex_label="$m_2$"
+)
priors["eccentricity"] = bilby.core.prior.LogUniform(
- name='eccentricity', latex_label='$e$', minimum=1e-4, maximum=0.4)
-priors["luminosity_distance"] = bilby.gw.prior.UniformComovingVolume(
- name='luminosity_distance', minimum=1e2, maximum=2e3)
-priors["dec"] = bilby.core.prior.Cosine(name='dec')
+ name="eccentricity", latex_label="$e$", minimum=1e-4, maximum=0.4
+)
+priors["luminosity_distance"] = bilby.gw.prior.UniformSourceFrame(
+ name="luminosity_distance", minimum=1e2, maximum=2e3
+)
+priors["dec"] = bilby.core.prior.Cosine(name="dec")
priors["ra"] = bilby.core.prior.Uniform(
- name='ra', minimum=0, maximum=2 * np.pi)
-priors["theta_jn"] = bilby.core.prior.Sine(name='theta_jn')
-priors["psi"] = bilby.core.prior.Uniform(name='psi', minimum=0, maximum=np.pi)
+ name="ra", minimum=0, maximum=2 * np.pi, boundary="periodic"
+)
+priors["theta_jn"] = bilby.core.prior.Sine(name="theta_jn")
+priors["psi"] = bilby.core.prior.Uniform(
+ name="psi", minimum=0, maximum=np.pi, boundary="periodic"
+)
priors["phase"] = bilby.core.prior.Uniform(
- name='phase', minimum=0, maximum=2 * np.pi)
+ name="phase", minimum=0, maximum=2 * np.pi, boundary="periodic"
+)
priors["geocent_time"] = bilby.core.prior.Uniform(
- 1180002600.9, 1180002601.1, name='geocent_time', unit='s')
+ injection_parameters["geocent_time"] - 0.1,
+ injection_parameters["geocent_time"] + 0.1,
+ name="geocent_time",
+ unit="s",
+)
# Initialising the likelihood function.
likelihood = bilby.gw.likelihood.GravitationalWaveTransient(
- interferometers=ifos, waveform_generator=waveform_generator)
+ interferometers=ifos,
+ waveform_generator=waveform_generator,
+ priors=priors,
+ time_marginalization=True,
+ distance_marginalization=True,
+ phase_marginalization=True,
+)
# Now we run sampler (PyMultiNest in our case).
result = bilby.run_sampler(
- likelihood=likelihood, priors=priors, sampler='pymultinest', npoints=1000,
- injection_parameters=injection_parameters, outdir=outdir, label=label)
+ likelihood=likelihood,
+ priors=priors,
+ sampler="pymultinest",
+ npoints=500,
+ injection_parameters=injection_parameters,
+ outdir=outdir,
+ label=label,
+)
# And finally we make some plots of the output posteriors.
result.plot_corner()
diff --git a/examples/gw_examples/injection_examples/fake_sampler_example.py b/examples/gw_examples/injection_examples/fake_sampler_example.py
index e1e7d600d88be737405d12dcdd21be3d5fbcba90..e4632ccb633d6306a8cae177cbe26b7c17a6252f 100755
--- a/examples/gw_examples/injection_examples/fake_sampler_example.py
+++ b/examples/gw_examples/injection_examples/fake_sampler_example.py
@@ -1,119 +1,143 @@
#!/usr/bin/env python
"""
-Read ROQ posterior and calculate full likelihood at same parameter space points.
+Demonstrate using the FakeSampler to reweight a result to include higher-order
+emission modes. This is a simplified version of the method presented in
+arXiv:1905.05477, however, the method can be applied to a much wider range of
+initial and more complex likelihoods.
"""
-import numpy as np
-import deepdish as dd
import bilby
import matplotlib.pyplot as plt
+import numpy as np
-def make_comparison_histograms(file_full, file_roq):
- # Returns a dictionary
- data_full = dd.io.load(file_full)
- data_roq = dd.io.load(file_roq)
-
- # These are pandas dataframes
- pos_full = data_full['posterior']
- pos_roq = data_roq['posterior']
+def make_comparison_histograms(result_1, result_2):
+ pos_full = result_1.posterior
+ pos_simple = result_2.posterior
plt.figure()
- plt.hist(pos_full['log_likelihood_evaluations'], 50, label='full', histtype='step')
- plt.hist(pos_roq['log_likelihood_evaluations'], 50, label='roq', histtype='step')
- plt.xlabel(r'delta_logl')
+ plt.hist(
+ pos_full["log_likelihood"],
+ 50,
+ label=result_1.label,
+ histtype="step",
+ density=True,
+ )
+ plt.hist(
+ pos_simple["log_likelihood"],
+ 50,
+ label=result_2.label,
+ histtype="step",
+ density=True,
+ )
+ plt.xlabel(r"delta_logl")
plt.legend(loc=2)
- plt.savefig('delta_logl.pdf')
- plt.close()
-
- plt.figure()
- delta_dlogl = pos_full['log_likelihood_evaluations'] - pos_roq['log_likelihood_evaluations']
- plt.hist(delta_dlogl, 50)
- plt.xlabel(r'delta_logl_full - delta_logl_roq')
- plt.savefig('delta_delta_logl.pdf')
- plt.close()
-
- plt.figure()
- delta_dlogl = np.abs(pos_full['log_likelihood_evaluations'] - pos_roq['log_likelihood_evaluations'])
- bins = np.logspace(np.log10(delta_dlogl.min()), np.log10(delta_dlogl.max()), 25)
- plt.hist(delta_dlogl, bins=bins)
- plt.xscale('log')
- plt.xlabel(r'|delta_logl_full - delta_logl_roq|')
- plt.savefig('log_abs_delta_delta_logl.pdf')
+ plt.savefig(f"{result_1.outdir}/delta_logl.pdf")
plt.close()
def main():
- outdir = 'outdir_full'
- label = 'full'
+ outdir = "outdir"
np.random.seed(170808)
duration = 4
- sampling_frequency = 2048
- noise = 'zero'
-
- sampler = 'fake_sampler'
- # This example assumes that the following posterior file exists.
- # It comes from a run using the full likelihood using the same
- # injection and sampling parameters, but the ROQ likelihood.
- # See roq_example.py for such an example.
- sample_file = 'outdir_dynesty_zero_noise_SNR22/roq_result.h5'
+ sampling_frequency = 1024
injection_parameters = dict(
- chirp_mass=36., mass_ratio=0.9, a_1=0.4, a_2=0.3, tilt_1=0.0, tilt_2=0.0,
- phi_12=1.7, phi_jl=0.3, luminosity_distance=2000., iota=0.4, psi=0.659,
- phase=1.3, geocent_time=1126259642.413, ra=1.375, dec=-1.2108)
-
- waveform_arguments = dict(waveform_approximant='IMRPhenomPv2',
- reference_frequency=20.0, minimum_frequency=20.0)
+ chirp_mass=36.0,
+ mass_ratio=0.2,
+ chi_1=0.4,
+ chi_2=0.3,
+ luminosity_distance=2000.0,
+ theta_jn=0.4,
+ psi=0.659,
+ phase=1.3,
+ geocent_time=1126259642.413,
+ ra=1.375,
+ dec=-1.2108,
+ )
+
+ waveform_arguments = dict(
+ waveform_approximant="IMRPhenomXAS",
+ reference_frequency=20.0,
+ minimum_frequency=20.0,
+ )
waveform_generator = bilby.gw.WaveformGenerator(
- duration=duration, sampling_frequency=sampling_frequency,
+ duration=duration,
+ sampling_frequency=sampling_frequency,
frequency_domain_source_model=bilby.gw.source.lal_binary_black_hole,
waveform_arguments=waveform_arguments,
- parameter_conversion=bilby.gw.conversion.convert_to_lal_binary_black_hole_parameters)
-
- ifos = bilby.gw.detector.InterferometerList(['H1', 'L1', 'V1'])
-
- if noise == 'Gaussian':
- ifos.set_strain_data_from_power_spectral_densities(
- sampling_frequency=sampling_frequency, duration=duration,
- start_time=injection_parameters['geocent_time'] - 3)
- elif noise == 'zero':
- ifos.set_strain_data_from_zero_noise(
- sampling_frequency=sampling_frequency, duration=duration,
- start_time=injection_parameters['geocent_time'] - 3)
-
- ifos.inject_signal(waveform_generator=waveform_generator,
- parameters=injection_parameters)
-
- priors = bilby.gw.prior.BBHPriorDict()
- for key in ['a_1', 'a_2', 'tilt_1', 'tilt_2', 'iota', 'psi', 'ra',
- 'dec', 'phi_12', 'phi_jl', 'luminosity_distance']:
+ parameter_conversion=bilby.gw.conversion.convert_to_lal_binary_black_hole_parameters,
+ )
+
+ ifos = bilby.gw.detector.InterferometerList(["H1", "L1", "V1"])
+
+ ifos.set_strain_data_from_zero_noise(
+ sampling_frequency=sampling_frequency,
+ duration=duration,
+ start_time=injection_parameters["geocent_time"] - 2,
+ )
+
+ ifos.inject_signal(
+ waveform_generator=waveform_generator,
+ parameters=injection_parameters,
+ )
+
+ priors = bilby.gw.prior.BBHPriorDict(aligned_spin=True)
+ for key in [
+ "luminosity_distance",
+ "theta_jn",
+ "phase",
+ "psi",
+ "ra",
+ "dec",
+ "geocent_time",
+ ]:
priors[key] = injection_parameters[key]
- priors.pop('mass_1')
- priors.pop('mass_2')
- priors['chirp_mass'] = bilby.core.prior.Uniform(
- 15, 40, latex_label='$\\mathcal{M}$')
- priors['mass_ratio'] = bilby.core.prior.Uniform(0.5, 1, latex_label='$q$')
- priors['geocent_time'] = bilby.core.prior.Uniform(
- injection_parameters['geocent_time'] - 0.1,
- injection_parameters['geocent_time'] + 0.1, latex_label='$t_c$', unit='s')
-
- likelihood = bilby.gw.GravitationalWaveTransient(
- interferometers=ifos, waveform_generator=waveform_generator)
-
- result = bilby.run_sampler(
- likelihood=likelihood, priors=priors, sampler=sampler, sample_file=sample_file,
- injection_parameters=injection_parameters, outdir=outdir, label=label)
-
- # Make a corner plot.
- result.plot_corner()
-
- # Compare full and ROQ likelihoods
- make_comparison_histograms(outdir + '/%s_result.h5' % label, sample_file)
-
-
-if __name__ == '__main__':
+ priors["chirp_mass"].maximum = 45
+
+ likelihood = bilby.gw.likelihood.GravitationalWaveTransient(
+ interferometers=ifos, waveform_generator=waveform_generator
+ )
+
+ # perform the initial sampling with our simple model
+ original_result = bilby.run_sampler(
+ likelihood=likelihood,
+ priors=priors,
+ sampler="dynesty",
+ walks=10,
+ nact=3,
+ bound="single",
+ injection_parameters=injection_parameters,
+ outdir=outdir,
+ label="primary_mode_only",
+ save="hdf5",
+ )
+
+ # update the waveform generator to use our higher-order mode waveform
+ likelihood.waveform_generator.waveform_arguments[
+ "waveform_approximant"
+ ] = "IMRPhenomXHM"
+
+ # call the FakeSampler to compute the new likelihoods
+ new_result = bilby.run_sampler(
+ likelihood=likelihood,
+ priors=priors,
+ sampler="fake_sampler",
+ sample_file=f"{outdir}/{original_result.label}_result.hdf5",
+ injection_parameters=injection_parameters,
+ outdir=outdir,
+ verbose=False,
+ label="higher_order_mode",
+ save="hdf5",
+ )
+
+ # make some comparison plots
+ bilby.core.result.plot_multiple([original_result, new_result])
+ make_comparison_histograms(new_result, original_result)
+
+
+if __name__ == "__main__":
main()
diff --git a/examples/gw_examples/injection_examples/fast_tutorial.py b/examples/gw_examples/injection_examples/fast_tutorial.py
index 3a2e18fc135b6225e369e28643d5749d0b9db50d..f462008a70994982d47c9cbd69179161c96aae3c 100644
--- a/examples/gw_examples/injection_examples/fast_tutorial.py
+++ b/examples/gw_examples/injection_examples/fast_tutorial.py
@@ -8,18 +8,18 @@ and distance using a uniform in comoving volume prior on luminosity distance
between luminosity distances of 100Mpc and 5Gpc, the cosmology is Planck15.
"""
-import numpy as np
import bilby
+import numpy as np
# Set the duration and sampling frequency of the data segment that we're
# going to inject the signal into
-duration = 4.
-sampling_frequency = 2048.
+duration = 4.0
+sampling_frequency = 2048.0
minimum_frequency = 20
# Specify the output directory and the name of the simulation.
-outdir = 'outdir'
-label = 'fast_tutorial'
+outdir = "outdir"
+label = "fast_tutorial"
bilby.core.utils.setup_logger(outdir=outdir, label=label)
# Set up a random seed for result reproducibility. This is optional!
@@ -30,31 +30,51 @@ np.random.seed(88170235)
# 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., theta_jn=0.4, psi=2.659,
- phase=1.3, geocent_time=1126259642.413, ra=1.375, dec=-1.2108)
+ mass_1=36.0,
+ mass_2=29.0,
+ 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.0,
+ 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='IMRPhenomPv2',
- reference_frequency=50.,
- minimum_frequency=minimum_frequency)
+waveform_arguments = dict(
+ waveform_approximant="IMRPhenomPv2",
+ reference_frequency=50.0,
+ minimum_frequency=minimum_frequency,
+)
# Create the waveform_generator using a LAL BinaryBlackHole source function
waveform_generator = bilby.gw.WaveformGenerator(
- duration=duration, sampling_frequency=sampling_frequency,
+ duration=duration,
+ sampling_frequency=sampling_frequency,
frequency_domain_source_model=bilby.gw.source.lal_binary_black_hole,
parameter_conversion=bilby.gw.conversion.convert_to_lal_binary_black_hole_parameters,
- waveform_arguments=waveform_arguments)
+ 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 = 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)
+ sampling_frequency=sampling_frequency,
+ duration=duration,
+ start_time=injection_parameters["geocent_time"] - 2,
+)
+ifos.inject_signal(
+ waveform_generator=waveform_generator, parameters=injection_parameters
+)
# Set up a PriorDict, which inherits from dict.
# By default we will sample all terms in the signal models. However, this will
@@ -67,12 +87,19 @@ ifos.inject_signal(waveform_generator=waveform_generator,
# distance, which means those are the parameters that will be included in the
# sampler. If we do nothing, then the default priors get used.
priors = bilby.gw.prior.BBHPriorDict()
-priors['geocent_time'] = bilby.core.prior.Uniform(
- minimum=injection_parameters['geocent_time'] - 0.1,
- maximum=injection_parameters['geocent_time'] + 0.1,
- name='geocent_time', latex_label='$t_c$', unit='$s$')
-for key in ['a_1', 'a_2', 'tilt_1', 'tilt_2', 'phi_12', 'phi_jl', 'psi', 'ra',
- 'dec', 'geocent_time', 'phase']:
+for key in [
+ "a_1",
+ "a_2",
+ "tilt_1",
+ "tilt_2",
+ "phi_12",
+ "phi_jl",
+ "psi",
+ "ra",
+ "dec",
+ "geocent_time",
+ "phase",
+]:
priors[key] = injection_parameters[key]
# Perform a check that the prior does not extend to a parameter space longer than the data
@@ -81,12 +108,19 @@ priors.validate_prior(duration, minimum_frequency)
# Initialise the likelihood by passing in the interferometer data (ifos) and
# the waveform generator
likelihood = bilby.gw.GravitationalWaveTransient(
- interferometers=ifos, waveform_generator=waveform_generator)
+ interferometers=ifos, waveform_generator=waveform_generator
+)
# Run sampler. In this case we're going to use the `dynesty` sampler
result = bilby.run_sampler(
- likelihood=likelihood, priors=priors, sampler='dynesty', npoints=1000,
- injection_parameters=injection_parameters, outdir=outdir, label=label)
+ likelihood=likelihood,
+ priors=priors,
+ sampler="dynesty",
+ npoints=1000,
+ injection_parameters=injection_parameters,
+ outdir=outdir,
+ label=label,
+)
# Make a corner plot.
result.plot_corner()
diff --git a/examples/gw_examples/injection_examples/how_to_specify_the_prior.py b/examples/gw_examples/injection_examples/how_to_specify_the_prior.py
index 851a3aab37cfbc1ce915af5e80e104dddf4b3a80..816288c5a727dc482c1524f7357e638004271228 100644
--- a/examples/gw_examples/injection_examples/how_to_specify_the_prior.py
+++ b/examples/gw_examples/injection_examples/how_to_specify_the_prior.py
@@ -4,58 +4,94 @@ Tutorial to demonstrate how to specify the prior distributions used for
parameter estimation.
"""
-import numpy as np
import bilby
+import numpy as np
-
-duration = 4.
-sampling_frequency = 2048.
-outdir = 'outdir'
+duration = 4
+sampling_frequency = 1024
+outdir = "outdir"
np.random.seed(151012)
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=4000., theta_jn=0.4, psi=2.659,
- phase=1.3, geocent_time=1126259642.413, ra=1.375, dec=-1.2108)
+ mass_1=36.0,
+ mass_2=29.0,
+ 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=4000.0,
+ theta_jn=0.4,
+ psi=2.659,
+ phase=1.3,
+ geocent_time=1126259642.413,
+ ra=1.375,
+ dec=-1.2108,
+)
-waveform_arguments = dict(waveform_approximant='IMRPhenomPv2',
- reference_frequency=50., minimum_frequency=20.)
+waveform_arguments = dict(
+ waveform_approximant="IMRPhenomXPHM",
+ reference_frequency=50.0,
+ minimum_frequency=20.0,
+)
# Create the waveform_generator using a LAL BinaryBlackHole source function
waveform_generator = bilby.gw.WaveformGenerator(
- duration=duration, sampling_frequency=sampling_frequency,
+ duration=duration,
+ sampling_frequency=sampling_frequency,
frequency_domain_source_model=bilby.gw.source.lal_binary_black_hole,
- waveform_arguments=waveform_arguments)
+ waveform_arguments=waveform_arguments,
+)
# Set up interferometers.
-ifos = bilby.gw.detector.InterferometerList(['H1', 'L1'])
+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)
+ sampling_frequency=sampling_frequency,
+ duration=duration,
+ start_time=injection_parameters["geocent_time"] - 2,
+)
+ifos.inject_signal(
+ waveform_generator=waveform_generator, parameters=injection_parameters
+)
# Set up prior
# This loads in a predefined set of priors for BBHs.
priors = bilby.gw.prior.BBHPriorDict()
# These parameters will not be sampled
-for key in ['tilt_1', 'tilt_2', 'phi_12', 'phi_jl', 'phase', 'theta_jn', 'ra',
- 'dec', 'geocent_time', 'psi']:
+for key in [
+ "tilt_1",
+ "tilt_2",
+ "phi_12",
+ "phi_jl",
+ "phase",
+ "theta_jn",
+ "ra",
+ "dec",
+ "geocent_time",
+ "psi",
+]:
priors[key] = injection_parameters[key]
# We can make uniform distributions.
-priors['mass_2'] = bilby.core.prior.Uniform(
- name='mass_2', minimum=20, maximum=40, unit='$M_{\\odot}$')
+del priors["chirp_mass"], priors["mass_ratio"]
+# We can make uniform distributions.
+priors["mass_1"] = bilby.core.prior.Uniform(
+ name="mass_1", minimum=20, maximum=40, unit="$M_{\\odot}$"
+)
+priors["mass_2"] = bilby.core.prior.Uniform(
+ name="mass_2", minimum=20, maximum=40, unit="$M_{\\odot}$"
+)
# We can make a power-law distribution, p(x) ~ x^{alpha}
# Note: alpha=0 is a uniform distribution, alpha=-1 is uniform-in-log
-priors['a_1'] = bilby.core.prior.PowerLaw(
- name='a_1', alpha=-1, minimum=1e-2, maximum=1)
+priors["a_1"] = bilby.core.prior.PowerLaw(name="a_1", alpha=-1, minimum=1e-2, maximum=1)
# We can define a prior from an array as follows.
# Note: this doesn't have to be properly normalised.
a_2 = np.linspace(0, 1, 1001)
-p_a_2 = a_2 ** 4
-priors['a_2'] = bilby.core.prior.Interped(
- name='a_2', xx=a_2, yy=p_a_2, minimum=0, maximum=0.5)
+p_a_2 = a_2**4
+priors["a_2"] = bilby.core.prior.Interped(
+ name="a_2", xx=a_2, yy=p_a_2, minimum=0, maximum=0.5
+)
# Additionally, we have Gaussian, TruncatedGaussian, Sine and Cosine.
# It's also possible to load an interpolate a prior from a file.
# Finally, if you don't specify any necessary parameters it will be filled in
@@ -64,10 +100,16 @@ priors['a_2'] = bilby.core.prior.Interped(
# Initialise GravitationalWaveTransient
likelihood = bilby.gw.GravitationalWaveTransient(
- interferometers=ifos, waveform_generator=waveform_generator)
+ interferometers=ifos, waveform_generator=waveform_generator
+)
# Run sampler
result = bilby.run_sampler(
- likelihood=likelihood, priors=priors, sampler='dynesty', outdir=outdir,
- injection_parameters=injection_parameters, label='specify_prior')
+ likelihood=likelihood,
+ priors=priors,
+ sampler="dynesty",
+ outdir=outdir,
+ injection_parameters=injection_parameters,
+ label="specify_prior",
+)
result.plot_corner()
diff --git a/examples/gw_examples/injection_examples/marginalized_likelihood.py b/examples/gw_examples/injection_examples/marginalized_likelihood.py
index c28efca13d4d9ddda064d4e6eef14c72b8086407..55cd2d4c7774613b5f520cbadccbad915804aa03 100644
--- a/examples/gw_examples/injection_examples/marginalized_likelihood.py
+++ b/examples/gw_examples/injection_examples/marginalized_likelihood.py
@@ -9,55 +9,90 @@ parameter can be recovered in post-processing.
import bilby
import numpy as np
-
-duration = 4.
-sampling_frequency = 2048.
-outdir = 'outdir'
-label = 'marginalized_likelihood'
+duration = 4
+sampling_frequency = 1024
+outdir = "outdir"
+label = "marginalized_likelihood"
np.random.seed(170608)
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=4000., theta_jn=0.4, psi=2.659,
- phase=1.3, geocent_time=1126259642.413, ra=1.375, dec=-1.2108)
+ mass_1=36.0,
+ mass_2=29.0,
+ 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=4000.0,
+ theta_jn=0.4,
+ psi=2.659,
+ phase=1.3,
+ geocent_time=1126259642.413,
+ ra=1.375,
+ dec=-1.2108,
+)
-waveform_arguments = dict(waveform_approximant='IMRPhenomPv2',
- reference_frequency=50.)
+waveform_arguments = dict(waveform_approximant="IMRPhenomXP", reference_frequency=50)
# Create the waveform_generator using a LAL BinaryBlackHole source function
waveform_generator = bilby.gw.WaveformGenerator(
- duration=duration, sampling_frequency=sampling_frequency,
+ duration=duration,
+ sampling_frequency=sampling_frequency,
frequency_domain_source_model=bilby.gw.source.lal_binary_black_hole,
- waveform_arguments=waveform_arguments)
+ waveform_arguments=waveform_arguments,
+)
# Set up interferometers.
-ifos = bilby.gw.detector.InterferometerList(['H1', 'L1'])
+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)
+ sampling_frequency=sampling_frequency,
+ duration=duration,
+ start_time=injection_parameters["geocent_time"] - 2,
+)
+ifos.inject_signal(
+ waveform_generator=waveform_generator, parameters=injection_parameters
+)
# Set up prior
priors = bilby.gw.prior.BBHPriorDict()
-priors['geocent_time'] = bilby.core.prior.Uniform(
- minimum=injection_parameters['geocent_time'] - 1,
- maximum=injection_parameters['geocent_time'] + 1,
- name='geocent_time', latex_label='$t_c$', unit='$s$')
+priors["geocent_time"] = bilby.core.prior.Uniform(
+ minimum=injection_parameters["geocent_time"] - 0.1,
+ maximum=injection_parameters["geocent_time"] + 0.1,
+ name="geocent_time",
+ latex_label="$t_c$",
+ unit="$s$",
+)
# These parameters will not be sampled
-for key in ['a_1', 'a_2', 'tilt_1', 'tilt_2', 'phi_12', 'phi_jl', 'theta_jn', 'ra',
- 'dec', 'mass_1', 'mass_2']:
+for key in [
+ "a_1",
+ "a_2",
+ "tilt_1",
+ "tilt_2",
+ "phi_12",
+ "phi_jl",
+ "theta_jn",
+ "ra",
+ "dec",
+ "mass_1",
+ "mass_2",
+]:
priors[key] = injection_parameters[key]
+del priors["chirp_mass"], priors["mass_ratio"]
# Initialise GravitationalWaveTransient
# Note that we now need to pass the: priors and flags for each thing that's
# being marginalised. A lookup table is used for distance marginalisation which
# takes a few minutes to build.
likelihood = bilby.gw.GravitationalWaveTransient(
- interferometers=ifos, waveform_generator=waveform_generator, priors=priors,
- distance_marginalization=True, phase_marginalization=True,
- time_marginalization=True)
+ interferometers=ifos,
+ waveform_generator=waveform_generator,
+ priors=priors,
+ distance_marginalization=True,
+ phase_marginalization=True,
+ time_marginalization=True,
+)
# Run sampler
# Note that we've added an additional argument `conversion_function`, this is
@@ -66,7 +101,12 @@ likelihood = bilby.gw.GravitationalWaveTransient(
# reconstructs posterior distributions for the parameters which were
# marginalised over in the likelihood.
result = bilby.run_sampler(
- likelihood=likelihood, priors=priors, sampler='dynesty',
- injection_parameters=injection_parameters, outdir=outdir, label=label,
- conversion_function=bilby.gw.conversion.generate_all_bbh_parameters)
+ likelihood=likelihood,
+ priors=priors,
+ sampler="dynesty",
+ injection_parameters=injection_parameters,
+ outdir=outdir,
+ label=label,
+ conversion_function=bilby.gw.conversion.generate_all_bbh_parameters,
+)
result.plot_corner()
diff --git a/examples/gw_examples/injection_examples/multiband_example.py b/examples/gw_examples/injection_examples/multiband_example.py
index ce6c9f436319a96695d0c9337b3e69701d454bec..e43b6ca2190a399442723639b9fd79ef06c418d6 100644
--- a/examples/gw_examples/injection_examples/multiband_example.py
+++ b/examples/gw_examples/injection_examples/multiband_example.py
@@ -4,72 +4,104 @@ Example of how to use the multi-banding method (see Morisaki, (2021) arXiv:
2104.07813) for a binary neutron star simulated signal in Gaussian noise.
"""
-import numpy as np
-
import bilby
+import numpy as np
-outdir = 'outdir'
-label = 'multibanding'
+outdir = "outdir"
+label = "multibanding"
np.random.seed(170808)
-minimum_frequency = 20.
-reference_frequency = 100.
-duration = 256.
-sampling_frequency = 4096.
-approximant = 'IMRPhenomD'
+minimum_frequency = 20
+reference_frequency = 100
+duration = 256
+sampling_frequency = 4096
+approximant = "IMRPhenomD"
injection_parameters = dict(
- chirp_mass=1.2, mass_ratio=0.8, chi_1=0., chi_2=0.,
- ra=3.44616, dec=-0.408084, luminosity_distance=200.,
- theta_jn=0.4, psi=0.659, phase=1.3, geocent_time=1187008882)
+ chirp_mass=1.2,
+ mass_ratio=0.8,
+ chi_1=0.0,
+ chi_2=0.0,
+ ra=3.44616,
+ dec=-0.408084,
+ luminosity_distance=200.0,
+ theta_jn=0.4,
+ psi=0.659,
+ phase=1.3,
+ geocent_time=1187008882,
+)
# inject signal
waveform_generator = bilby.gw.WaveformGenerator(
- duration=duration, sampling_frequency=sampling_frequency,
+ duration=duration,
+ sampling_frequency=sampling_frequency,
frequency_domain_source_model=bilby.gw.source.lal_binary_black_hole,
- waveform_arguments=dict(waveform_approximant=approximant,
- reference_frequency=reference_frequency),
- parameter_conversion=bilby.gw.conversion.convert_to_lal_binary_black_hole_parameters)
-ifos = bilby.gw.detector.InterferometerList(['H1', 'L1', 'V1'])
+ waveform_arguments=dict(
+ waveform_approximant=approximant, reference_frequency=reference_frequency
+ ),
+ parameter_conversion=bilby.gw.conversion.convert_to_lal_binary_black_hole_parameters,
+)
+ifos = bilby.gw.detector.InterferometerList(["H1", "L1", "V1"])
ifos.set_strain_data_from_power_spectral_densities(
- sampling_frequency=sampling_frequency, duration=duration,
- start_time=injection_parameters['geocent_time'] - duration + 2.)
-ifos.inject_signal(waveform_generator=waveform_generator,
- parameters=injection_parameters)
+ sampling_frequency=sampling_frequency,
+ duration=duration,
+ start_time=injection_parameters["geocent_time"] - duration + 2,
+)
+ifos.inject_signal(
+ waveform_generator=waveform_generator, parameters=injection_parameters
+)
for ifo in ifos:
ifo.minimum_frequency = minimum_frequency
# make waveform generator for likelihood evaluations
search_waveform_generator = bilby.gw.waveform_generator.WaveformGenerator(
- duration=duration, sampling_frequency=sampling_frequency,
+ duration=duration,
+ sampling_frequency=sampling_frequency,
frequency_domain_source_model=bilby.gw.source.binary_black_hole_frequency_sequence,
- waveform_arguments=dict(waveform_approximant=approximant,
- reference_frequency=reference_frequency),
- parameter_conversion=bilby.gw.conversion.convert_to_lal_binary_black_hole_parameters)
+ waveform_arguments=dict(
+ waveform_approximant=approximant, reference_frequency=reference_frequency
+ ),
+ parameter_conversion=bilby.gw.conversion.convert_to_lal_binary_black_hole_parameters,
+)
# make prior
priors = bilby.gw.prior.BNSPriorDict()
-priors['chi_1'] = 0.
-priors['chi_2'] = 0.
-priors.pop('lambda_1')
-priors.pop('lambda_2')
-priors['chirp_mass'] = bilby.core.prior.Uniform(name='chirp_mass', minimum=1.15, maximum=1.25)
-priors['geocent_time'] = bilby.core.prior.Uniform(
- injection_parameters['geocent_time'] - 0.1,
- injection_parameters['geocent_time'] + 0.1, latex_label='$t_c$', unit='s')
+priors["chi_1"] = 0
+priors["chi_2"] = 0
+del priors["lambda_1"], priors["lambda_2"]
+priors["chirp_mass"] = bilby.core.prior.Uniform(
+ name="chirp_mass", minimum=1.15, maximum=1.25
+)
+priors["geocent_time"] = bilby.core.prior.Uniform(
+ injection_parameters["geocent_time"] - 0.1,
+ injection_parameters["geocent_time"] + 0.1,
+ latex_label="$t_c$",
+ unit="s",
+)
# make multi-banded likelihood
likelihood = bilby.gw.likelihood.MBGravitationalWaveTransient(
- interferometers=ifos, waveform_generator=search_waveform_generator, priors=priors,
- reference_chirp_mass=priors['chirp_mass'].minimum,
- distance_marginalization=True, phase_marginalization=True
+ interferometers=ifos,
+ waveform_generator=search_waveform_generator,
+ priors=priors,
+ reference_chirp_mass=priors["chirp_mass"].minimum,
+ distance_marginalization=True,
+ phase_marginalization=True,
)
# sampling
result = bilby.run_sampler(
- likelihood=likelihood, priors=priors, sampler='dynesty',
- nlive=500, walks=100, maxmcmc=5000, nact=5,
- injection_parameters=injection_parameters, outdir=outdir, label=label)
+ likelihood=likelihood,
+ priors=priors,
+ sampler="dynesty",
+ nlive=500,
+ walks=100,
+ maxmcmc=5000,
+ nact=5,
+ injection_parameters=injection_parameters,
+ outdir=outdir,
+ label=label,
+)
# Make a corner plot.
result.plot_corner()
diff --git a/examples/gw_examples/injection_examples/non_tensor.py b/examples/gw_examples/injection_examples/non_tensor.py
index 6c44216954f92aab978cac1be612d3ebbc17b042..1020a10002a869afffb9d82384844649b67f8c61 100644
--- a/examples/gw_examples/injection_examples/non_tensor.py
+++ b/examples/gw_examples/injection_examples/non_tensor.py
@@ -28,69 +28,99 @@ def vector_tensor_sine_gaussian(frequency_array, hrss, Q, frequency, epsilon):
Relative size of the vector modes compared to the tensor modes.
"""
waveform_polarizations = bilby.gw.source.sinegaussian(
- frequency_array, hrss, Q, frequency)
+ frequency_array, hrss, Q, frequency
+ )
- waveform_polarizations['x'] = epsilon * waveform_polarizations['plus']
- waveform_polarizations['y'] = epsilon * waveform_polarizations['cross']
+ waveform_polarizations["x"] = epsilon * waveform_polarizations["plus"]
+ waveform_polarizations["y"] = epsilon * waveform_polarizations["cross"]
return waveform_polarizations
-duration = 4.
-sampling_frequency = 2048.
+duration = 1
+sampling_frequency = 512
-outdir = 'outdir'
-label = 'vector_tensor'
+outdir = "outdir"
+label = "vector_tensor"
bilby.core.utils.setup_logger(outdir=outdir, label=label)
np.random.seed(170801)
injection_parameters = dict(
- hrss=1e-22, Q=5.0, frequency=200.0, ra=1.375, dec=-1.2108,
- geocent_time=1126259642.413, psi=2.659, epsilon=0.2)
-
-waveform_generator =\
- bilby.gw.waveform_generator.WaveformGenerator(
- duration=duration, sampling_frequency=sampling_frequency,
- frequency_domain_source_model=vector_tensor_sine_gaussian)
-
-ifos = bilby.gw.detector.InterferometerList(['H1', 'L1', 'V1'])
+ hrss=1e-22,
+ Q=5.0,
+ frequency=200.0,
+ ra=1.375,
+ dec=-1.2108,
+ geocent_time=1126259642.413,
+ psi=2.659,
+ epsilon=0.2,
+)
+
+waveform_generator = bilby.gw.waveform_generator.WaveformGenerator(
+ duration=duration,
+ sampling_frequency=sampling_frequency,
+ frequency_domain_source_model=vector_tensor_sine_gaussian,
+)
+
+ifos = bilby.gw.detector.InterferometerList(["H1", "L1", "V1"])
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)
-
-priors = dict()
-for key in ['psi', 'geocent_time', 'hrss', 'Q', 'frequency']:
+ sampling_frequency=sampling_frequency,
+ duration=duration,
+ start_time=injection_parameters["geocent_time"] - 0.5,
+)
+ifos.inject_signal(
+ waveform_generator=waveform_generator,
+ parameters=injection_parameters,
+ raise_error=False,
+)
+
+priors = bilby.core.prior.PriorDict()
+for key in ["psi", "geocent_time", "hrss", "Q", "frequency"]:
priors[key] = injection_parameters[key]
-priors['ra'] = bilby.core.prior.Uniform(0, 2 * np.pi, latex_label='$\\alpha$')
-priors['dec'] = bilby.core.prior.Cosine(latex_label='$\\delta$')
-priors['epsilon'] = bilby.core.prior.Uniform(0, 1, latex_label='$\\epsilon$')
+priors["ra"] = bilby.core.prior.Uniform(0, 2 * np.pi, latex_label="$\\alpha$")
+priors["dec"] = bilby.core.prior.Cosine(latex_label="$\\delta$")
+priors["epsilon"] = bilby.core.prior.Uniform(0, 1, latex_label="$\\epsilon$")
vector_tensor_likelihood = bilby.gw.likelihood.GravitationalWaveTransient(
- interferometers=ifos, waveform_generator=waveform_generator)
+ interferometers=ifos, waveform_generator=waveform_generator
+)
# Run sampler. In this case we're going to use the `dynesty` sampler
vector_tensor_result = bilby.core.sampler.run_sampler(
- likelihood=vector_tensor_likelihood, priors=priors, sampler='nestle',
- npoints=1000, injection_parameters=injection_parameters,
- outdir=outdir, label='vector_tensor')
+ likelihood=vector_tensor_likelihood,
+ priors=priors,
+ sampler="nestle",
+ nlive=1000,
+ injection_parameters=injection_parameters,
+ outdir=outdir,
+ label="vector_tensor",
+)
vector_tensor_result.plot_corner()
tensor_likelihood = bilby.gw.likelihood.GravitationalWaveTransient(
- interferometers=ifos, waveform_generator=waveform_generator)
+ interferometers=ifos, waveform_generator=waveform_generator
+)
-priors['epsilon'] = 0
+priors["epsilon"] = 0
-# Run sampler. In this case we're going to use the `dynesty` sampler
+# Run sampler. In this case we're going to use the `nestle` sampler
tensor_result = bilby.core.sampler.run_sampler(
- likelihood=tensor_likelihood, priors=priors, sampler='nestle', npoints=1000,
- injection_parameters=injection_parameters, outdir=outdir, label='tensor')
+ likelihood=tensor_likelihood,
+ priors=priors,
+ sampler="nestle",
+ nlive=1000,
+ injection_parameters=injection_parameters,
+ outdir=outdir,
+ label="tensor",
+)
# make some plots of the outputs
tensor_result.plot_corner()
bilby.result.plot_multiple(
- [tensor_result, vector_tensor_result], labels=['Tensor', 'Vector + Tensor'],
- parameters=dict(ra=1.375, dec=-1.2108), evidences=True)
+ [tensor_result, vector_tensor_result],
+ labels=["Tensor", "Vector + Tensor"],
+ parameters=dict(ra=1.375, dec=-1.2108),
+ evidences=True,
+)
diff --git a/examples/gw_examples/injection_examples/plot_skymap.py b/examples/gw_examples/injection_examples/plot_skymap.py
index 4c7ecb9a3e3abbaa4bab6099330b1332d5846909..7d0dbb98eca15dab0f19ea5e83506614e527182c 100644
--- a/examples/gw_examples/injection_examples/plot_skymap.py
+++ b/examples/gw_examples/injection_examples/plot_skymap.py
@@ -5,47 +5,85 @@ skymap
"""
import bilby
-duration = 4.
-sampling_frequency = 2048.
-outdir = 'outdir'
-label = 'plot_skymap'
+duration = 4
+sampling_frequency = 1024
+outdir = "outdir"
+label = "plot_skymap"
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=4000., theta_jn=0.4, psi=2.659,
- phase=1.3, geocent_time=1126259642.413, ra=1.375, dec=-0.2108)
+ mass_1=36.0,
+ mass_2=29.0,
+ 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=4000.0,
+ theta_jn=0.4,
+ psi=2.659,
+ phase=1.3,
+ geocent_time=1126259642.413,
+ ra=1.375,
+ dec=-0.2108,
+)
-waveform_arguments = dict(waveform_approximant='IMRPhenomPv2',
- reference_frequency=50.)
+waveform_arguments = dict(waveform_approximant="IMRPhenomXP", reference_frequency=50.0)
waveform_generator = bilby.gw.WaveformGenerator(
- duration=duration, sampling_frequency=sampling_frequency,
+ duration=duration,
+ sampling_frequency=sampling_frequency,
frequency_domain_source_model=bilby.gw.source.lal_binary_black_hole,
- parameters=injection_parameters, waveform_arguments=waveform_arguments)
+ parameters=injection_parameters,
+ waveform_arguments=waveform_arguments,
+)
-ifos = bilby.gw.detector.InterferometerList(['H1', 'L1'])
+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)
+ sampling_frequency=sampling_frequency,
+ duration=duration,
+ start_time=injection_parameters["geocent_time"] - 2,
+)
+ifos.inject_signal(
+ waveform_generator=waveform_generator, parameters=injection_parameters
+)
priors = bilby.gw.prior.BBHPriorDict()
-for key in ['a_1', 'a_2', 'tilt_1', 'tilt_2', 'phi_12', 'phi_jl', 'psi',
- 'mass_1', 'mass_2', 'phase', 'geocent_time', 'luminosity_distance',
- 'theta_jn']:
+for key in [
+ "a_1",
+ "a_2",
+ "tilt_1",
+ "tilt_2",
+ "phi_12",
+ "phi_jl",
+ "psi",
+ "mass_1",
+ "mass_2",
+ "phase",
+ "geocent_time",
+ "theta_jn",
+]:
priors[key] = injection_parameters[key]
+del priors["chirp_mass"], priors["mass_ratio"]
likelihood = bilby.gw.GravitationalWaveTransient(
- interferometers=ifos, waveform_generator=waveform_generator,
- time_marginalization=True, phase_marginalization=True,
- distance_marginalization=False, priors=priors)
+ interferometers=ifos, waveform_generator=waveform_generator
+)
result = bilby.run_sampler(
- likelihood=likelihood, priors=priors, sampler='dynesty', npoints=1000,
- injection_parameters=injection_parameters, outdir=outdir, label=label, result_class=bilby.gw.result.CBCResult)
+ likelihood=likelihood,
+ priors=priors,
+ sampler="nestle",
+ npoints=250,
+ injection_parameters=injection_parameters,
+ outdir=outdir,
+ label=label,
+ result_class=bilby.gw.result.CBCResult,
+)
# make some plots of the outputs
result.plot_corner()
# will require installation of ligo.skymap (pip install ligo.skymap)
-result.plot_skymap(maxpts=5000)
+# the skymap generation code is fairly slow when using many points so limit
+# ourselves to 500 points in the fit
+result.plot_skymap(maxpts=500)
diff --git a/examples/gw_examples/injection_examples/plot_time_domain_data.py b/examples/gw_examples/injection_examples/plot_time_domain_data.py
index 4e572e8a53a330a0be47a9bfc5f2b50b6879cd0a..707c0f18c101a3d7ac59e3c143a699a7ed1c88db 100644
--- a/examples/gw_examples/injection_examples/plot_time_domain_data.py
+++ b/examples/gw_examples/injection_examples/plot_time_domain_data.py
@@ -1,37 +1,64 @@
#!/usr/bin/env python
"""
+This example demonstrates how to simulate some data, add an injected signal
+and plot the data.
"""
import numpy as np
-import bilby
+from bilby.gw.detector import get_empty_interferometer
+from bilby.gw.source import lal_binary_black_hole
+from bilby.gw.waveform_generator import WaveformGenerator
np.random.seed(1)
duration = 4
-sampling_frequency = 2048.
+sampling_frequency = 2048.0
-outdir = 'outdir'
-label = 'example'
+outdir = "outdir"
+label = "example"
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=1000., theta_jn=0.4, psi=2.659,
- phase=1.3, geocent_time=1126259642.413, ra=1.375, dec=-1.2108)
+ mass_1=36.0,
+ mass_2=29.0,
+ 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=1000.0,
+ theta_jn=0.4,
+ psi=2.659,
+ phase=1.3,
+ geocent_time=1126259642.413,
+ ra=1.375,
+ dec=-1.2108,
+)
-waveform_arguments = dict(waveform_approximant='IMRPhenomPv2',
- reference_frequency=50.)
+waveform_arguments = dict(
+ waveform_approximant="IMRPhenomTPHM", reference_frequency=50.0
+)
-waveform_generator = bilby.gw.WaveformGenerator(
- duration=duration, sampling_frequency=sampling_frequency,
- frequency_domain_source_model=bilby.gw.source.lal_binary_black_hole,
- parameters=injection_parameters, waveform_arguments=waveform_arguments)
+waveform_generator = WaveformGenerator(
+ duration=duration,
+ sampling_frequency=sampling_frequency,
+ frequency_domain_source_model=lal_binary_black_hole,
+ parameters=injection_parameters,
+ waveform_arguments=waveform_arguments,
+)
hf_signal = waveform_generator.frequency_domain_strain(injection_parameters)
-H1 = bilby.gw.detector.get_interferometer_with_fake_noise_and_injection(
- 'H1', injection_polarizations=hf_signal,
- injection_parameters=injection_parameters, duration=duration,
- sampling_frequency=sampling_frequency, outdir=outdir)
+ifo = get_empty_interferometer("H1")
+ifo.set_strain_data_from_power_spectral_density(
+ duration=duration, sampling_frequency=sampling_frequency
+)
+ifo.inject_signal(injection_polarizations=hf_signal, parameters=injection_parameters)
-t0 = injection_parameters['geocent_time']
-H1.plot_time_domain_data(outdir=outdir, label=label, notches=[50],
- bandpass_frequencies=(50, 200), start_end=(-0.5, 0.5),
- t0=t0)
+t0 = injection_parameters["geocent_time"]
+ifo.plot_time_domain_data(
+ outdir=outdir,
+ label=label,
+ notches=[50],
+ bandpass_frequencies=(50, 200),
+ start_end=(-0.5, 0.5),
+ t0=t0,
+)
diff --git a/examples/gw_examples/injection_examples/reproduce_mpa1_eos.py b/examples/gw_examples/injection_examples/reproduce_mpa1_eos.py
index 699b56c102d6b4fb462ce1bd2ba16d1f73787899..6d2ee4fc2f1e0bf83391b41f5087ed016b80bfdc 100644
--- a/examples/gw_examples/injection_examples/reproduce_mpa1_eos.py
+++ b/examples/gw_examples/injection_examples/reproduce_mpa1_eos.py
@@ -9,19 +9,25 @@ from bilby.gw import eos
# First, we specify the spectral parameter values for the MPA1 EoS.
MPA1_gammas = [1.0215, 0.1653, -0.0235, -0.0004]
-MPA1_p0 = 1.51e33 # Pressure in CGS
-MPA1_e0_c2 = 2.04e14 # \epsilon_0 / c^2 in CGS
-MPA1_xmax = 6.63 # Dimensionless ending pressure
+MPA1_p0 = 1.51e33 # Pressure in CGS
+MPA1_e0_c2 = 2.04e14 # \epsilon_0 / c^2 in CGS
+MPA1_xmax = 6.63 # Dimensionless ending pressure
# Create the spectral decomposition EoS class
-MPA1_spectral = eos.SpectralDecompositionEOS(MPA1_gammas, p0=MPA1_p0, e0=MPA1_e0_c2, xmax=MPA1_xmax, npts=100)
+MPA1_spectral = eos.SpectralDecompositionEOS(
+ MPA1_gammas, p0=MPA1_p0, e0=MPA1_e0_c2, xmax=MPA1_xmax, npts=100
+)
# And create another from tabulated data
-MPA1_tabulated = eos.TabularEOS('MPA1')
+MPA1_tabulated = eos.TabularEOS("MPA1")
# Now let's plot them
# To do so, we specify a representation and plot ranges.
-MPA1_spectral_plot = MPA1_spectral.plot('pressure-energy_density', xlim=[1e22, 1e36], ylim=[1e9, 1e36])
-MPA1_tabular_plot = MPA1_tabulated.plot('pressure-energy_density', xlim=[1e22, 1e36], ylim=[1e9, 1e36])
-MPA1_spectral_plot.savefig('spectral_mpa1.pdf')
-MPA1_tabular_plot.savefig('tabular_mpa1.pdf')
+MPA1_spectral_plot = MPA1_spectral.plot(
+ "pressure-energy_density", xlim=[1e22, 1e36], ylim=[1e9, 1e36]
+)
+MPA1_tabular_plot = MPA1_tabulated.plot(
+ "pressure-energy_density", xlim=[1e22, 1e36], ylim=[1e9, 1e36]
+)
+MPA1_spectral_plot.savefig("spectral_mpa1.pdf")
+MPA1_tabular_plot.savefig("tabular_mpa1.pdf")
diff --git a/examples/gw_examples/injection_examples/roq_example.py b/examples/gw_examples/injection_examples/roq_example.py
index 82f87d3e418620c5b86f136b2a81b0da7bfec3b7..254e29aad0e15e27075186cdf40ac5ae579b09e1 100644
--- a/examples/gw_examples/injection_examples/roq_example.py
+++ b/examples/gw_examples/injection_examples/roq_example.py
@@ -7,14 +7,16 @@ Gaussian noise.
This requires files specifying the appropriate basis weights.
These aren't shipped with Bilby, but are available on LDG clusters and
from the public repository https://git.ligo.org/lscsoft/ROQ_data.
-"""
-import numpy as np
+We also reweight the result using the regular waveform model to check how
+correct the ROQ result is.
+"""
import bilby
+import numpy as np
-outdir = 'outdir'
-label = 'roq'
+outdir = "outdir"
+label = "roq"
# The ROQ bases can be given an overall scaling.
# Note how this is applied to durations, frequencies and masses.
@@ -33,9 +35,9 @@ freq_nodes_quadratic = np.load("fnodes_quadratic.npy") * scale_factor
params = np.genfromtxt("params.dat", names=True)
# Get scaled ROQ quantities
-minimum_chirp_mass = params['chirpmassmin'] / scale_factor
-maximum_chirp_mass = params['chirpmassmax'] / scale_factor
-minimum_component_mass = params['compmin'] / scale_factor
+minimum_chirp_mass = params["chirpmassmin"] / scale_factor
+maximum_chirp_mass = params["chirpmassmax"] / scale_factor
+minimum_component_mass = params["compmin"] / scale_factor
np.random.seed(170808)
@@ -43,73 +45,144 @@ duration = 4 / scale_factor
sampling_frequency = 2048 * scale_factor
injection_parameters = dict(
- mass_1=36.0, mass_2=29.0, a_1=0.4, a_2=0.3, tilt_1=0.0, tilt_2=0.0,
- phi_12=1.7, phi_jl=0.3, luminosity_distance=1000., theta_jn=0.4, psi=0.659,
- phase=1.3, geocent_time=1126259642.413, ra=1.375, dec=-1.2108)
-
-waveform_arguments = dict(waveform_approximant='IMRPhenomPv2',
- reference_frequency=20. * scale_factor)
+ mass_1=36.0,
+ mass_2=29.0,
+ a_1=0.4,
+ a_2=0.3,
+ tilt_1=0.0,
+ tilt_2=0.0,
+ phi_12=1.7,
+ phi_jl=0.3,
+ luminosity_distance=1000.0,
+ theta_jn=0.4,
+ psi=0.659,
+ phase=1.3,
+ geocent_time=1126259642.413,
+ ra=1.375,
+ dec=-1.2108,
+)
+
+waveform_arguments = dict(
+ waveform_approximant="IMRPhenomPv2", reference_frequency=20.0 * scale_factor
+)
waveform_generator = bilby.gw.WaveformGenerator(
- duration=duration, sampling_frequency=sampling_frequency,
+ duration=duration,
+ sampling_frequency=sampling_frequency,
frequency_domain_source_model=bilby.gw.source.lal_binary_black_hole,
waveform_arguments=waveform_arguments,
- parameter_conversion=bilby.gw.conversion.convert_to_lal_binary_black_hole_parameters)
+ parameter_conversion=bilby.gw.conversion.convert_to_lal_binary_black_hole_parameters,
+)
-ifos = bilby.gw.detector.InterferometerList(['H1', 'L1', 'V1'])
+ifos = bilby.gw.detector.InterferometerList(["H1", "L1", "V1"])
ifos.set_strain_data_from_power_spectral_densities(
- sampling_frequency=sampling_frequency, duration=duration,
- start_time=injection_parameters['geocent_time'] - 3 / scale_factor)
-ifos.inject_signal(waveform_generator=waveform_generator,
- parameters=injection_parameters)
+ sampling_frequency=sampling_frequency,
+ duration=duration,
+ start_time=injection_parameters["geocent_time"] - 2 / scale_factor,
+)
+ifos.inject_signal(
+ waveform_generator=waveform_generator, parameters=injection_parameters
+)
for ifo in ifos:
ifo.minimum_frequency = 20 * scale_factor
# make ROQ waveform generator
search_waveform_generator = bilby.gw.waveform_generator.WaveformGenerator(
- duration=duration, sampling_frequency=sampling_frequency,
+ duration=duration,
+ sampling_frequency=sampling_frequency,
frequency_domain_source_model=bilby.gw.source.binary_black_hole_roq,
waveform_arguments=dict(
frequency_nodes_linear=freq_nodes_linear,
frequency_nodes_quadratic=freq_nodes_quadratic,
- reference_frequency=20. * scale_factor, waveform_approximant='IMRPhenomPv2'),
- parameter_conversion=bilby.gw.conversion.convert_to_lal_binary_black_hole_parameters)
+ reference_frequency=20.0 * scale_factor,
+ waveform_approximant="IMRPhenomPv2",
+ ),
+ parameter_conversion=bilby.gw.conversion.convert_to_lal_binary_black_hole_parameters,
+)
# Here we add constraints on chirp mass and mass ratio to the prior, these are
# determined by the domain of validity of the ROQ basis.
priors = bilby.gw.prior.BBHPriorDict()
-for key in ['a_1', 'a_2', 'tilt_1', 'tilt_2', 'theta_jn', 'phase', 'psi', 'ra',
- 'dec', 'phi_12', 'phi_jl', 'luminosity_distance']:
+for key in [
+ "a_1",
+ "a_2",
+ "tilt_1",
+ "tilt_2",
+ "theta_jn",
+ "phase",
+ "psi",
+ "ra",
+ "dec",
+ "phi_12",
+ "phi_jl",
+ "luminosity_distance",
+]:
priors[key] = injection_parameters[key]
-for key in ['mass_1', 'mass_2']:
+for key in ["mass_1", "mass_2"]:
priors[key].minimum = max(priors[key].minimum, minimum_component_mass)
-priors['chirp_mass'] = bilby.core.prior.Uniform(
- name='chirp_mass', minimum=float(minimum_chirp_mass),
- maximum=float(maximum_chirp_mass))
-priors['mass_ratio'] = bilby.core.prior.Uniform(0.125, 1, name='mass_ratio')
-priors['geocent_time'] = bilby.core.prior.Uniform(
- injection_parameters['geocent_time'] - 0.1,
- injection_parameters['geocent_time'] + 0.1, latex_label='$t_c$', unit='s')
+priors["chirp_mass"] = bilby.core.prior.Uniform(
+ name="chirp_mass",
+ minimum=float(minimum_chirp_mass),
+ maximum=float(maximum_chirp_mass),
+)
+# The roq parameters typically store the mass ratio bounds as m1/m2 not m2/m1 as in the
+# Bilby convention.
+priors["mass_ratio"] = bilby.core.prior.Uniform(
+ 1 / params["qmax"], 1, name="mass_ratio"
+)
+priors["geocent_time"] = bilby.core.prior.Uniform(
+ injection_parameters["geocent_time"] - 0.1,
+ injection_parameters["geocent_time"] + 0.1,
+ latex_label="$t_c$",
+ unit="s",
+)
likelihood = bilby.gw.likelihood.ROQGravitationalWaveTransient(
- interferometers=ifos, waveform_generator=search_waveform_generator,
- linear_matrix=basis_matrix_linear, quadratic_matrix=basis_matrix_quadratic,
- priors=priors, roq_params=params, roq_scale_factor=scale_factor)
+ interferometers=ifos,
+ waveform_generator=search_waveform_generator,
+ linear_matrix=basis_matrix_linear,
+ quadratic_matrix=basis_matrix_quadratic,
+ priors=priors,
+ roq_params=params,
+ roq_scale_factor=scale_factor,
+)
# write the weights to file so they can be loaded multiple times
-likelihood.save_weights('weights.npz')
+likelihood.save_weights("weights.npz")
# remove the basis matrices as these are big for longer bases
del basis_matrix_linear, basis_matrix_quadratic
# load the weights from the file
likelihood = bilby.gw.likelihood.ROQGravitationalWaveTransient(
- interferometers=ifos, waveform_generator=search_waveform_generator,
- weights='weights.npz', priors=priors)
+ interferometers=ifos,
+ waveform_generator=search_waveform_generator,
+ weights="weights.npz",
+ priors=priors,
+)
result = bilby.run_sampler(
- likelihood=likelihood, priors=priors, sampler='dynesty', npoints=500,
- injection_parameters=injection_parameters, outdir=outdir, label=label)
-
-# Make a corner plot.
-result.plot_corner()
+ likelihood=likelihood,
+ priors=priors,
+ sampler="dynesty",
+ npoints=500,
+ injection_parameters=injection_parameters,
+ outdir=outdir,
+ label=label,
+)
+
+# Resample the result using the full waveform model with the FakeSampler.
+# This will give us an idea of how good a job the ROQ does.
+full_likelihood = bilby.gw.likelihood.GravitationalWaveTransient(
+ interferometers=ifos, waveform_generator=waveform_generator
+)
+resampled_result = bilby.run_sampler(
+ likelihood=full_likelihood,
+ priors=priors,
+ sampler="fake_sampler",
+ label="roq_resampled",
+ outdir=outdir,
+)
+
+# Make a comparison corner plot with the two likelihoods.
+bilby.core.result.plot_multiple([result, resampled_result], labels=["ROQ", "Regular"])
diff --git a/examples/gw_examples/injection_examples/sine_gaussian_example.py b/examples/gw_examples/injection_examples/sine_gaussian_example.py
index d40cb68e89e37a28409c8a292dad4176117a9c20..d253826ea6129ec56a9619b64d097a147b179504 100644
--- a/examples/gw_examples/injection_examples/sine_gaussian_example.py
+++ b/examples/gw_examples/injection_examples/sine_gaussian_example.py
@@ -8,12 +8,12 @@ import numpy as np
# Set the duration and sampling frequency of the data segment that we're going
# to inject the signal into
-duration = 4.
-sampling_frequency = 2048.
+duration = 1
+sampling_frequency = 512
# Specify the output directory and the name of the simulation.
-outdir = 'outdir'
-label = 'sine_gaussian'
+outdir = "outdir"
+label = "sine_gaussian"
bilby.core.utils.setup_logger(outdir=outdir, label=label)
# Set up a random seed for result reproducibility. This is optional!
@@ -23,42 +23,64 @@ np.random.seed(170801)
# dictionary of parameters that includes all of the different waveform
# parameters
injection_parameters = dict(
- hrss=1e-22, Q=5.0, frequency=200.0, ra=1.375, dec=-1.2108,
- geocent_time=1126259642.413, psi=2.659)
+ hrss=1e-22,
+ Q=5.0,
+ frequency=200.0,
+ ra=1.375,
+ dec=-1.2108,
+ geocent_time=1126259642.413,
+ psi=2.659,
+)
# Create the waveform_generator using a sine Gaussian source function
waveform_generator = bilby.gw.waveform_generator.WaveformGenerator(
- duration=duration, sampling_frequency=sampling_frequency,
- frequency_domain_source_model=bilby.gw.source.sinegaussian)
+ duration=duration,
+ sampling_frequency=sampling_frequency,
+ frequency_domain_source_model=bilby.gw.source.sinegaussian,
+)
# Set up interferometers. In this case we'll use three interferometers
# (LIGO-Hanford (H1), LIGO-Livingston (L1), and Virgo (V1)). These default to
# their design sensitivity
-ifos = bilby.gw.detector.InterferometerList(['H1', 'L1', 'V1'])
+ifos = bilby.gw.detector.InterferometerList(["H1", "L1", "V1"])
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)
+ sampling_frequency=sampling_frequency,
+ duration=duration,
+ start_time=injection_parameters["geocent_time"] - 0.5,
+)
+ifos.inject_signal(
+ waveform_generator=waveform_generator,
+ parameters=injection_parameters,
+ raise_error=False,
+)
-# Set up prior, which is a dictionary
-priors = dict()
-for key in ['psi', 'ra', 'dec', 'geocent_time']:
+# Set up the prior. We will fix the "extrinsic" parameters to their true values.
+priors = bilby.core.prior.PriorDict()
+for key in ["psi", "ra", "dec", "geocent_time"]:
priors[key] = injection_parameters[key]
-priors['Q'] = bilby.core.prior.Uniform(2, 50, 'Q')
-priors['frequency'] = bilby.core.prior.Uniform(30, 1000, 'frequency', unit='Hz')
-priors['hrss'] = bilby.core.prior.Uniform(1e-23, 1e-21, 'hrss')
+priors["Q"] = bilby.core.prior.Uniform(2, 50, "Q")
+priors["frequency"] = bilby.core.prior.Uniform(160, 240, "frequency", unit="Hz")
+priors["hrss"] = bilby.core.prior.Uniform(1e-23, 1e-21, "hrss")
# Initialise the likelihood by passing in the interferometer data (IFOs) and
# the waveoform generator
likelihood = bilby.gw.likelihood.GravitationalWaveTransient(
- interferometers=ifos, waveform_generator=waveform_generator)
+ interferometers=ifos, waveform_generator=waveform_generator
+)
# Run sampler. In this case we're going to use the `dynesty` sampler
result = bilby.core.sampler.run_sampler(
- likelihood=likelihood, priors=priors, sampler='dynesty', npoints=1000,
- injection_parameters=injection_parameters, outdir=outdir, label=label)
+ likelihood=likelihood,
+ priors=priors,
+ sampler="dynesty",
+ nlive=1000,
+ walks=10,
+ nact=5,
+ injection_parameters=injection_parameters,
+ outdir=outdir,
+ label=label,
+)
# make some plots of the outputs
result.plot_corner()
diff --git a/examples/gw_examples/injection_examples/standard_15d_cbc_tutorial.py b/examples/gw_examples/injection_examples/standard_15d_cbc_tutorial.py
index 0e084c051a740537621b2c1e55dfeb88c62708fe..db8222e1e7b86425905225b5483d4b298546b2ed 100644
--- a/examples/gw_examples/injection_examples/standard_15d_cbc_tutorial.py
+++ b/examples/gw_examples/injection_examples/standard_15d_cbc_tutorial.py
@@ -3,18 +3,20 @@
Tutorial to demonstrate running parameter estimation on a full 15 parameter
space for an injected cbc signal. This is the standard injection analysis script
one can modify for the study of injected CBC events.
+
+This will take many hours to run.
"""
-import numpy as np
import bilby
+import numpy as np
# Set the duration and sampling frequency of the data segment that we're
# going to inject the signal into
-duration = 4.
-sampling_frequency = 2048.
+duration = 4.0
+sampling_frequency = 1024.0
# Specify the output directory and the name of the simulation.
-outdir = 'outdir'
-label = 'full_15_parameters'
+outdir = "outdir"
+label = "full_15_parameters"
bilby.core.utils.setup_logger(outdir=outdir, label=label)
# Set up a random seed for result reproducibility. This is optional!
@@ -25,75 +27,120 @@ np.random.seed(88170235)
# 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., theta_jn=0.4, psi=2.659,
- phase=1.3, geocent_time=1126259642.413, ra=1.375, dec=-1.2108)
+ mass_1=36.0,
+ mass_2=29.0,
+ 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.0,
+ 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='IMRPhenomPv2',
- reference_frequency=50., minimum_frequency=20.)
+waveform_arguments = dict(
+ waveform_approximant="IMRPhenomXPHM",
+ reference_frequency=50.0,
+ minimum_frequency=20.0,
+)
# Create the waveform_generator using a LAL BinaryBlackHole source function
# the generator will convert all the parameters
waveform_generator = bilby.gw.WaveformGenerator(
- duration=duration, sampling_frequency=sampling_frequency,
+ duration=duration,
+ sampling_frequency=sampling_frequency,
frequency_domain_source_model=bilby.gw.source.lal_binary_black_hole,
parameter_conversion=bilby.gw.conversion.convert_to_lal_binary_black_hole_parameters,
- waveform_arguments=waveform_arguments)
+ 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 = 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)
+ sampling_frequency=sampling_frequency,
+ duration=duration,
+ start_time=injection_parameters["geocent_time"] - 2,
+)
+ifos.inject_signal(
+ waveform_generator=waveform_generator, parameters=injection_parameters
+)
# For this analysis, we implement the standard BBH priors defined, except for
# the definition of the time prior, which is defined as uniform about the
# injected value.
-# Furthermore, we decide to sample in chirp mass and mass ratio, due to the
-# preferred shape for the associated posterior distributions.
+# We change the mass boundaries to be more targeted for the source we
+# injected.
+# We define priors in the time at the Hanford interferometer and two
+# parameters (zenith, azimuth) defining the sky position wrt the two
+# interferometers.
priors = bilby.gw.prior.BBHPriorDict()
-priors.pop('mass_1')
-priors.pop('mass_2')
-priors['chirp_mass'] = bilby.prior.Uniform(
- name='chirp_mass', latex_label='$M$', minimum=10.0, maximum=100.0,
- unit='$M_{\\odot}$')
-
-priors['mass_ratio'] = bilby.prior.Uniform(
- name='mass_ratio', latex_label='$q$', minimum=0.5, maximum=1.0)
-
-priors['geocent_time'] = bilby.core.prior.Uniform(
- minimum=injection_parameters['geocent_time'] - 0.1,
- maximum=injection_parameters['geocent_time'] + 0.1,
- name='geocent_time', latex_label='$t_c$', unit='$s$')
+time_delay = ifos[0].time_delay_from_geocenter(
+ injection_parameters["ra"],
+ injection_parameters["dec"],
+ injection_parameters["geocent_time"],
+)
+priors["H1_time"] = bilby.core.prior.Uniform(
+ minimum=injection_parameters["geocent_time"] + time_delay - 0.1,
+ maximum=injection_parameters["geocent_time"] + time_delay + 0.1,
+ name="H1_time",
+ latex_label="$t_H$",
+ unit="$s$",
+)
+del priors["ra"], priors["dec"]
+priors["zenith"] = bilby.core.prior.Sine(latex_label="$\\kappa$")
+priors["azimuth"] = bilby.core.prior.Uniform(
+ minimum=0, maximum=2 * np.pi, latex_label="$\\epsilon$", boundary="periodic"
+)
# Initialise the likelihood by passing in the interferometer data (ifos) and
# the waveoform generator, as well the priors.
-# The explicit time, distance, and phase marginalizations are turned on to
-# improve convergence, and the parameters are recovered by the conversion
-# function.
+# The explicit distance marginalization is turned on to improve
+# convergence, and the posterior is recovered by the conversion function.
likelihood = bilby.gw.GravitationalWaveTransient(
- interferometers=ifos, waveform_generator=waveform_generator, priors=priors,
- distance_marginalization=True, phase_marginalization=True, time_marginalization=True)
-
-# Run sampler. In this case we're going to use the `cpnest` sampler
-# Note that the maxmcmc parameter is increased so that between each iteration of
-# the nested sampler approach, the walkers will move further using an mcmc
-# approach, searching the full parameter space.
-# The conversion function will determine the distance, phase and coalescence
-# time posteriors in post processing.
+ interferometers=ifos,
+ waveform_generator=waveform_generator,
+ priors=priors,
+ distance_marginalization=True,
+ phase_marginalization=False,
+ time_marginalization=False,
+ reference_frame="H1L1",
+ time_reference="H1",
+)
+
+# Run sampler. In this case we're going to use the `dynesty` sampler
+# Note that the `walks`, `nact`, and `maxmcmc` parameter are specified
+# to ensure sufficient convergence of the analysis.
+# We set `npool=4` to parallelize the analysis over four cores.
+# The conversion function will determine the distance posterior in post processing
result = bilby.run_sampler(
- likelihood=likelihood, priors=priors, sampler='cpnest', npoints=2000,
- injection_parameters=injection_parameters, outdir=outdir,
- label=label, maxmcmc=2000,
- conversion_function=bilby.gw.conversion.generate_all_bbh_parameters)
+ likelihood=likelihood,
+ priors=priors,
+ sampler="dynesty",
+ nlive=1000,
+ walks=20,
+ nact=50,
+ maxmcmc=2000,
+ npool=4,
+ injection_parameters=injection_parameters,
+ outdir=outdir,
+ label=label,
+ conversion_function=bilby.gw.conversion.generate_all_bbh_parameters,
+ result_class=bilby.gw.result.CBCResult,
+)
+
+# Plot the inferred waveform superposed on the actual data.
+result.plot_waveform_posterior(n_samples=1000)
# Make a corner plot.
result.plot_corner()
diff --git a/examples/gw_examples/injection_examples/using_gwin.py b/examples/gw_examples/injection_examples/using_gwin.py
deleted file mode 100644
index de68a82886ed8f262d0f5bef5ddf98629be13b8d..0000000000000000000000000000000000000000
--- a/examples/gw_examples/injection_examples/using_gwin.py
+++ /dev/null
@@ -1,92 +0,0 @@
-#!/usr/bin/env python
-"""
-An example of how to use bilby 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 bilby is that while fixed parameters
-are specified via the prior in bilby, 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.
-
-"""
-import numpy as np
-import bilby
-
-import gwin
-from pycbc import psd as pypsd
-from pycbc.waveform.generator import (FDomainDetFrameGenerator,
- FDomainCBCGenerator)
-
-label = 'using_gwin'
-
-# Search priors
-priors = dict()
-priors['distance'] = bilby.core.prior.Uniform(500, 2000, 'distance')
-priors['polarization'] = bilby.core.prior.Uniform(0, np.pi, 'theta_jn')
-
-# 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 bilby.likelihood object
-class GWINLikelihood(bilby.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 = bilby.run_sampler(
- likelihood=likelihood, priors=priors, sampler='dynesty', npoints=500,
- label=label)
-result.plot_corner()
diff --git a/examples/gw_examples/supernova_example/supernova_example.py b/examples/gw_examples/supernova_example/supernova_example.py
index b5de9f40775a91f27dce75a2a3fe549d6ca4a2c8..1d79fc650bcb6d63d5734d8715bccab1dc67ef0a 100644
--- a/examples/gw_examples/supernova_example/supernova_example.py
+++ b/examples/gw_examples/supernova_example/supernova_example.py
@@ -6,18 +6,22 @@ supernova injected signal. Signal model is made by applying PCA to a set of
supernova waveforms. The first few PCs are then linearly combined with a scale
factor. (See https://arxiv.org/pdf/1202.3256.pdf)
+For this example we use `PyMultiNest`, this can be installed using
+
+conda install -c conda-forge pymultinest
"""
-import numpy as np
import bilby
+import numpy as np
# Set the duration and sampling frequency of the data segment that we're going
-# to inject the signal into
-duration = 3.
-sampling_frequency = 4096.
+# to inject the signal into.
+# These are fixed by the resolution in the injection file that we are using.
+duration = 3
+sampling_frequency = 4096
# Specify the output directory and the name of the simulation.
-outdir = 'outdir'
-label = 'supernova'
+outdir = "outdir"
+label = "supernova"
bilby.core.utils.setup_logger(outdir=outdir, label=label)
# Set up a random seed for result reproducibility. This is optional!
@@ -27,70 +31,92 @@ np.random.seed(170801)
# of parameters that includes all of the different waveform parameters. It will
# read in a signal to inject from a txt file.
-injection_parameters = dict(file_path='MuellerL15_example_inj.txt',
- luminosity_distance=7.0, ra=4.6499,
- dec=-0.5063, geocent_time=1126259642.413,
- psi=2.659)
+injection_parameters = dict(
+ luminosity_distance=7.0,
+ ra=4.6499,
+ dec=-0.5063,
+ geocent_time=1126259642.413,
+ psi=2.659,
+)
# Create the waveform_generator using a supernova source function
waveform_generator = bilby.gw.waveform_generator.WaveformGenerator(
- duration=duration, sampling_frequency=sampling_frequency,
+ duration=duration,
+ sampling_frequency=sampling_frequency,
frequency_domain_source_model=bilby.gw.source.supernova,
- parameters=injection_parameters)
+ parameters=injection_parameters,
+ parameter_conversion=lambda parameters: (parameters, list()),
+ waveform_arguments=dict(file_path="MuellerL15_example_inj.txt"),
+)
# Set up interferometers. In this case we'll use three interferometers
# (LIGO-Hanford (H1), LIGO-Livingston (L1), and Virgo (V1)). These default to
# their design sensitivity
-ifos = bilby.gw.detector.InterferometerList(['H1', 'L1'])
+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)
+ sampling_frequency=sampling_frequency,
+ duration=duration,
+ start_time=injection_parameters["geocent_time"] - 1.5,
+)
+ifos.inject_signal(
+ waveform_generator=waveform_generator,
+ parameters=injection_parameters,
+ raise_error=False,
+)
# read in from a file the PCs used to create the signal model.
-realPCs = np.loadtxt('SupernovaRealPCs.txt')
-imagPCs = np.loadtxt('SupernovaImagPCs.txt')
+realPCs = np.genfromtxt("SupernovaRealPCs.txt")
+imagPCs = np.genfromtxt("SupernovaImagPCs.txt")
# Now we make another waveform_generator because the signal model is
# not the same as the injection in this case.
-simulation_parameters = dict(
- realPCs=realPCs, imagPCs=imagPCs)
+simulation_parameters = dict(realPCs=realPCs, imagPCs=imagPCs)
search_waveform_generator = bilby.gw.waveform_generator.WaveformGenerator(
- duration=duration, sampling_frequency=sampling_frequency,
+ duration=duration,
+ sampling_frequency=sampling_frequency,
frequency_domain_source_model=bilby.gw.source.supernova_pca_model,
- waveform_arguments=simulation_parameters)
+ waveform_arguments=simulation_parameters,
+)
# Set up prior
-priors = dict()
-for key in ['psi', 'geocent_time']:
+priors = bilby.core.prior.PriorDict()
+for key in ["psi", "geocent_time"]:
priors[key] = injection_parameters[key]
-priors['luminosity_distance'] = bilby.core.prior.Uniform(
- 2, 20, 'luminosity_distance', unit='$kpc$')
-priors['pc_coeff1'] = bilby.core.prior.Uniform(-1, 1, 'pc_coeff1')
-priors['pc_coeff2'] = bilby.core.prior.Uniform(-1, 1, 'pc_coeff2')
-priors['pc_coeff3'] = bilby.core.prior.Uniform(-1, 1, 'pc_coeff3')
-priors['pc_coeff4'] = bilby.core.prior.Uniform(-1, 1, 'pc_coeff4')
-priors['pc_coeff5'] = bilby.core.prior.Uniform(-1, 1, 'pc_coeff5')
-priors['ra'] = bilby.core.prior.Uniform(minimum=0, maximum=2 * np.pi,
- name='ra')
-priors['dec'] = bilby.core.prior.Sine(name='dec')
-priors['geocent_time'] = bilby.core.prior.Uniform(
- injection_parameters['geocent_time'] - 1,
- injection_parameters['geocent_time'] + 1,
- 'geocent_time', unit='$s$')
+priors["luminosity_distance"] = bilby.core.prior.Uniform(
+ 2, 20, "luminosity_distance", unit="$kpc$"
+)
+priors["pc_coeff1"] = bilby.core.prior.Uniform(-1, 1, "pc_coeff1")
+priors["pc_coeff2"] = bilby.core.prior.Uniform(-1, 1, "pc_coeff2")
+priors["pc_coeff3"] = bilby.core.prior.Uniform(-1, 1, "pc_coeff3")
+priors["pc_coeff4"] = bilby.core.prior.Uniform(-1, 1, "pc_coeff4")
+priors["pc_coeff5"] = bilby.core.prior.Uniform(-1, 1, "pc_coeff5")
+priors["ra"] = bilby.core.prior.Uniform(
+ minimum=0, maximum=2 * np.pi, name="ra", boundary="periodic"
+)
+priors["dec"] = bilby.core.prior.Sine(name="dec")
+priors["geocent_time"] = bilby.core.prior.Uniform(
+ injection_parameters["geocent_time"] - 1,
+ injection_parameters["geocent_time"] + 1,
+ "geocent_time",
+ unit="$s$",
+)
# Initialise the likelihood by passing in the interferometer data (IFOs) and
# the waveoform generator
likelihood = bilby.gw.likelihood.GravitationalWaveTransient(
- interferometers=ifos, waveform_generator=search_waveform_generator)
+ interferometers=ifos, waveform_generator=search_waveform_generator
+)
# Run sampler.
result = bilby.run_sampler(
- likelihood=likelihood, priors=priors, sampler='dynesty', npoints=500,
- outdir=outdir, label=label)
+ likelihood=likelihood,
+ priors=priors,
+ sampler="pymultinest",
+ npoints=500,
+ outdir=outdir,
+ label=label,
+)
# make some plots of the outputs
result.plot_corner()
-print(result)