diff --git a/CHANGELOG.md b/CHANGELOG.md
index 58c91df81b115d2d65941fce1d6ca732001779ca..169ab1ac37537eb82430a965d5b658cb33442a5f 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -5,6 +5,8 @@
Changes currently on master, but not under a tag.
### Changes
+- Make BBH/BNS parameter conversion more logical
+- Source frame masses/spins included in posterior
- Make filling in posterior with fixed parameters work
## [0.3] 2018-01-02
diff --git a/bilby/core/result.py b/bilby/core/result.py
index 3ca06858412997c83f0f78fde58e79b0598c2839..86faa4aa732370ad89a82f8f3d14e4351c7eef52 100644
--- a/bilby/core/result.py
+++ b/bilby/core/result.py
@@ -234,8 +234,9 @@ class Result(dict):
elif k in self.parameter_labels:
latex_labels.append(k)
else:
- raise ValueError('key {} not a parameter label or latex label'
- .format(k))
+ logger.info(
+ 'key {} not a parameter label or latex label'.format(k))
+ latex_labels.append(' '.join(k.split('_')))
return latex_labels
@property
@@ -588,39 +589,6 @@ class Result(dict):
self.prior_values[key]\
= priors[key].prob(self.posterior[key].values)
- def construct_cbc_derived_parameters(self):
- """ Construct widely used derived parameters of CBCs """
- self.posterior['mass_chirp'] = (
- (self.posterior.mass_1 * self.posterior.mass_2) ** 0.6 / (
- self.posterior.mass_1 + self.posterior.mass_2) ** 0.2)
- self.search_parameter_keys.append('mass_chirp')
- self.parameter_labels.append('$\mathcal{M}$')
-
- self.posterior['q'] = self.posterior.mass_2 / self.posterior.mass_1
- self.search_parameter_keys.append('q')
- self.parameter_labels.append('$q$')
-
- self.posterior['eta'] = (
- (self.posterior.mass_1 * self.posterior.mass_2) / (
- self.posterior.mass_1 + self.posterior.mass_2) ** 2)
- self.search_parameter_keys.append('eta')
- self.parameter_labels.append('$\eta$')
-
- self.posterior['chi_eff'] = (
- (self.posterior.a_1 * np.cos(self.posterior.tilt_1) +
- self.posterior.q * self.posterior.a_2 *
- np.cos(self.posterior.tilt_2)) / (1 + self.posterior.q))
- self.search_parameter_keys.append('chi_eff')
- self.parameter_labels.append('$\chi_{\mathrm eff}$')
-
- self.posterior['chi_p'] = (
- np.maximum(self.posterior.a_1 * np.sin(self.posterior.tilt_1),
- (4 * self.posterior.q + 3) / (3 * self.posterior.q + 4) *
- self.posterior.q * self.posterior.a_2 *
- np.sin(self.posterior.tilt_2)))
- self.search_parameter_keys.append('chi_p')
- self.parameter_labels.append('$\chi_{\mathrm p}$')
-
def check_attribute_match_to_other_object(self, name, other_object):
""" Check attribute name exists in other_object and is the same
diff --git a/bilby/gw/conversion.py b/bilby/gw/conversion.py
index 793fecb0b8862fd074aa469444e54929c3099634..ac40efe0cf323ae2c290d1a881da2251d6d3a8fa 100644
--- a/bilby/gw/conversion.py
+++ b/bilby/gw/conversion.py
@@ -1,6 +1,5 @@
from __future__ import division
import numpy as np
-import pandas as pd
from ..core.utils import logger, solar_mass
from ..core.prior import DeltaFunction
@@ -47,6 +46,53 @@ def luminosity_distance_to_comoving_distance(distance):
return redshift_to_comoving_distance(redshift)
+@np.vectorize
+def transform_precessing_spins(theta_jn, phi_jl, tilt_1, tilt_2, phi_12, a_1,
+ a_2, mass_1, mass_2, reference_frequency, phase):
+ """
+ Vectorized version of
+ lalsimulation.SimInspiralTransformPrecessingNewInitialConditions
+
+ All parameters are defined at the reference frequency
+
+ Parameters
+ ----------
+ theta_jn: float
+ Inclination angle
+ phi_jl: float
+ Spin phase angle
+ tilt_1: float
+ Primary object tilt
+ tilt_2: float
+ Secondary object tilt
+ phi_12: float
+ Relative spin azimuthal angle
+ a_1: float
+ Primary dimensionless spin magnitude
+ a_2: float
+ Secondary dimensionless spin magnitude
+ mass_1: float
+ Primary mass _in SI units_
+ mass_2: float
+ Secondary mass _in SI units_
+ reference_frequency: float
+ phase: float
+ Orbital phase
+
+ Returns
+ -------
+ iota: float
+ Transformed inclination
+ spin_1x, spin_1y, spin_1z, spin_2x, spin_2y, spin_2z: float
+ Cartesian spin components
+ """
+ iota, spin_1x, spin_1y, spin_1z, spin_2x, spin_2y, spin_2z = \
+ lalsim.SimInspiralTransformPrecessingNewInitialConditions(
+ theta_jn, phi_jl, tilt_1, tilt_2, phi_12, a_1, a_2, mass_1,
+ mass_2, reference_frequency, phase)
+ return iota, spin_1x, spin_1y, spin_1z, spin_2x, spin_2y, spin_2z
+
+
def convert_to_lal_binary_black_hole_parameters(parameters):
"""
Convert parameters we have into parameters we need.
@@ -164,7 +210,7 @@ def convert_to_lal_binary_neutron_star_parameters(parameters):
Mass: mass_1, mass_2
- Spin: a_1, a_2, tilt_1, tilt_2, phi_12, phi_jl
+ Spin: chi_1, chi_2, tilt_1, tilt_2, phi_12, phi_jl
Extrinsic: luminosity_distance, theta_jn, phase, ra, dec, geocent_time, psi
This involves popping a lot of things from parameters.
@@ -187,6 +233,14 @@ def convert_to_lal_binary_neutron_star_parameters(parameters):
converted_parameters, added_keys =\
convert_to_lal_binary_black_hole_parameters(converted_parameters)
+ # catch if tidal parameters aren't present
+ if not any([key in converted_parameters for key in
+ ['lambda_1', 'lambda_2', 'lambda_tilde', 'delta_lambda']]):
+ converted_parameters['lambda_1'] = 0
+ converted_parameters['lambda_2'] = 0
+ added_keys = added_keys + ['lambda_1', 'lambda_2']
+ return converted_parameters, added_keys
+
if 'delta_lambda' in converted_parameters.keys():
converted_parameters['lambda_1'], converted_parameters['lambda_2'] =\
lambda_tilde_delta_lambda_to_lambda_1_lambda_2(
@@ -406,6 +460,71 @@ def mass_1_and_chirp_mass_to_mass_ratio(mass_1, chirp_mass):
return mass_ratio
+def lambda_1_lambda_2_to_lambda_tilde(lambda_1, lambda_2, mass_1, mass_2):
+ """
+ Convert from individual tidal parameters to domainant tidal term.
+
+ See, e.g., Wade et al., https://arxiv.org/pdf/1402.5156.pdf.
+
+ Parameters
+ ----------
+ lambda_1: float
+ Tidal parameter of more massive neutron star.
+ lambda_2: float
+ Tidal parameter of less massive neutron star.
+ mass_1: float
+ Mass of more massive neutron star.
+ mass_2: float
+ Mass of less massive neutron star.
+
+ Return
+ ------
+ lambda_tilde: float
+ Dominant tidal term.
+ """
+ eta = component_masses_to_symmetric_mass_ratio(mass_1, mass_2)
+ lambda_plus = lambda_1 + lambda_2
+ lambda_minus = lambda_1 - lambda_2
+ lambda_tilde = 8 / 13 * (
+ (1 + 7 * eta - 31 * eta**2) * lambda_plus +
+ (1 - 4 * eta)**0.5 * (1 + 9 * eta - 11 * eta**2) * lambda_minus)
+
+ return lambda_tilde
+
+
+def lambda_1_lambda_2_to_delta_lambda(lambda_1, lambda_2, mass_1, mass_2):
+ """
+ Convert from individual tidal parameters to second domainant tidal term.
+
+ See, e.g., Wade et al., https://arxiv.org/pdf/1402.5156.pdf.
+
+ Parameters
+ ----------
+ lambda_1: float
+ Tidal parameter of more massive neutron star.
+ lambda_2: float
+ Tidal parameter of less massive neutron star.
+ mass_1: float
+ Mass of more massive neutron star.
+ mass_2: float
+ Mass of less massive neutron star.
+
+ Return
+ ------
+ delta_lambda: float
+ Second dominant tidal term.
+ """
+ eta = component_masses_to_symmetric_mass_ratio(mass_1, mass_2)
+ lambda_plus = lambda_1 + lambda_2
+ lambda_minus = lambda_1 - lambda_2
+ delta_lambda = 1 / 2 * (
+ (1 - 4 * eta) ** 0.5 * (1 - 13272 / 1319 * eta + 8944 / 1319 * eta**2) *
+ lambda_plus + (1 - 15910 / 1319 * eta + 32850 / 1319 * eta**2 +
+ 3380 / 1319 * eta**3) * lambda_minus)
+
+ return delta_lambda
+
+
def lambda_tilde_delta_lambda_to_lambda_1_lambda_2(
lambda_tilde, delta_lambda, mass_1, mass_2):
"""
@@ -484,6 +603,29 @@ def lambda_tilde_to_lambda_1_lambda_2(
return lambda_1, lambda_2
+def _generate_all_cbc_parameters(sample, defaults, base_conversion,
+ likelihood=None, priors=None):
+ """Generate all cbc parameters, helper function for BBH/BNS"""
+ output_sample = sample.copy()
+ waveform_defaults = defaults
+ for key in waveform_defaults:
+ try:
+ output_sample[key] = \
+ likelihood.waveform_generator.waveform_arguments[key]
+ except (KeyError, AttributeError):
+ default = waveform_defaults[key]
+ output_sample[key] = default
+ logger.warning('Assuming {} = {}'.format(key, default))
+
+ output_sample = fill_from_fixed_priors(output_sample, priors)
+ output_sample, _ = base_conversion(output_sample)
+ output_sample = generate_mass_parameters(output_sample)
+ output_sample = generate_spin_parameters(output_sample)
+ output_sample = generate_source_frame_parameters(output_sample)
+ compute_snrs(output_sample, likelihood)
+ return output_sample
+
+
def generate_all_bbh_parameters(sample, likelihood=None, priors=None):
"""
From either a single sample or a set of samples fill in all missing
@@ -500,21 +642,43 @@ def generate_all_bbh_parameters(sample, likelihood=None, priors=None):
priors: dict, optional
Dictionary of prior objects, used to fill in non-sampled parameters.
"""
- output_sample = sample.copy()
- if likelihood is not None:
- output_sample['reference_frequency'] =\
- likelihood.waveform_generator.waveform_arguments[
- 'reference_frequency']
- output_sample['waveform_approximant'] =\
- likelihood.waveform_generator.waveform_arguments[
- 'waveform_approximant']
+ waveform_defaults = {
+ 'reference_frequency': 50.0, 'waveform_approximant': 'IMRPhenomPv2',
+ 'minimum_frequency': 20.0}
+ output_sample = _generate_all_cbc_parameters(
+ sample, defaults=waveform_defaults,
+ base_conversion=convert_to_lal_binary_black_hole_parameters,
+ likelihood=likelihood, priors=priors)
+ return output_sample
- output_sample = fill_from_fixed_priors(output_sample, priors)
- output_sample, _ =\
- convert_to_lal_binary_black_hole_parameters(output_sample)
- output_sample = generate_non_standard_parameters(output_sample)
- output_sample = generate_component_spins(output_sample)
- compute_snrs(output_sample, likelihood)
+
+def generate_all_bns_parameters(sample, likelihood=None, priors=None):
+ """
+ From either a single sample or a set of samples fill in all missing
+ BNS parameters, in place.
+
+ Since we assume BNS waveforms are aligned, component spins won't be
+ calculated.
+
+ Parameters
+ ----------
+ sample: dict or pandas.DataFrame
+ Samples to fill in with extra parameters, this may be either an
+ injection or posterior samples.
+ likelihood: bilby.gw.likelihood.GravitationalWaveTransient, optional
+ GravitationalWaveTransient used for sampling, used for waveform and
+ likelihood.interferometers.
+ priors: dict, optional
+ Dictionary of prior objects, used to fill in non-sampled parameters.
+ """
+ waveform_defaults = {
+ 'reference_frequency': 50.0, 'waveform_approximant': 'TaylorF2',
+ 'minimum_frequency': 20.0}
+ output_sample = _generate_all_cbc_parameters(
+ sample, defaults=waveform_defaults,
+ base_conversion=convert_to_lal_binary_neutron_star_parameters,
+ likelihood=likelihood, priors=priors)
+ output_sample = generate_tidal_parameters(output_sample)
return output_sample
@@ -540,13 +704,12 @@ def fill_from_fixed_priors(sample, priors):
return output_sample
-def generate_non_standard_parameters(sample):
+def generate_mass_parameters(sample):
"""
- Add the known non-standard parameters to the data frame/dictionary.
+ Add the known mass parameters to the data frame/dictionary.
We add:
- chirp mass, total mass, symmetric mass ratio, mass ratio, cos tilt 1,
- cos tilt 2, cos iota
+ chirp mass, total mass, symmetric mass ratio, mass ratio
Parameters
----------
@@ -569,12 +732,47 @@ def generate_non_standard_parameters(sample):
output_sample['mass_ratio'] =\
component_masses_to_mass_ratio(sample['mass_1'], sample['mass_2'])
- output_sample['cos_tilt_1'] = np.cos(output_sample['tilt_1'])
- output_sample['cos_tilt_2'] = np.cos(output_sample['tilt_2'])
- output_sample['cos_iota'] = np.cos(output_sample['iota'])
+ return output_sample
+
+
+def generate_spin_parameters(sample):
+ """
+ Add all spin parameters to the data frame/dictionary.
+
+ We add:
+ cartestian spin components, chi_eff, chi_p cos tilt 1, cos tilt 2
+
+ Parameters
+ ----------
+ sample : dict, pandas.DataFrame
+ The input dictionary with some spin parameters
+
+ Returns
+ -------
+ dict: The updated dictionary
+
+ """
+ output_sample = sample.copy()
+
+ output_sample = generate_component_spins(output_sample)
+
+ output_sample['chi_eff'] = (output_sample['spin_1z'] +
+ output_sample['spin_2z'] *
+ output_sample['mass_ratio']) /\
+ (1 + output_sample['mass_ratio'])
+
+ output_sample['chi_p'] = np.maximum(
+ (output_sample['spin_1x'] ** 2 + output_sample['spin_1y']**2)**0.5,
+ (4 * output_sample['mass_ratio'] + 3) /
+ (3 * output_sample['mass_ratio'] + 4) * output_sample['mass_ratio'] *
+ (output_sample['spin_2x'] ** 2 + output_sample['spin_2y']**2)**0.5)
+
+ try:
+ output_sample['cos_tilt_1'] = np.cos(output_sample['tilt_1'])
+ output_sample['cos_tilt_2'] = np.cos(output_sample['tilt_2'])
+ except KeyError:
+ pass
- output_sample['redshift'] =\
- luminosity_distance_to_redshift(sample['luminosity_distance'])
return output_sample
@@ -599,13 +797,12 @@ def generate_component_spins(sample):
spin_conversion_parameters =\
['iota', 'phi_jl', 'tilt_1', 'tilt_2', 'phi_12', 'a_1', 'a_2', 'mass_1',
'mass_2', 'reference_frequency', 'phase']
- if all(key in output_sample.keys() for key in spin_conversion_parameters)\
- and isinstance(output_sample, dict):
+ if all(key in output_sample.keys() for key in spin_conversion_parameters):
output_sample['iota'], output_sample['spin_1x'],\
output_sample['spin_1y'], output_sample['spin_1z'], \
output_sample['spin_2x'], output_sample['spin_2y'],\
output_sample['spin_2z'] =\
- lalsim.SimInspiralTransformPrecessingNewInitialConditions(
+ transform_precessing_spins(
output_sample['iota'], output_sample['phi_jl'],
output_sample['tilt_1'], output_sample['tilt_2'],
output_sample['phi_12'], output_sample['a_1'],
@@ -618,40 +815,73 @@ def generate_component_spins(sample):
np.arctan(output_sample['spin_1y'] / output_sample['spin_1x'])
output_sample['phi_2'] =\
np.arctan(output_sample['spin_2y'] / output_sample['spin_2x'])
+ elif 'chi_1' in output_sample and 'chi_2' in output_sample:
+ output_sample['spin_1x'] = 0
+ output_sample['spin_1y'] = 0
+ output_sample['spin_1z'] = output_sample['chi_1']
+ output_sample['spin_2x'] = 0
+ output_sample['spin_2y'] = 0
+ output_sample['spin_2z'] = output_sample['chi_2']
+ else:
+ logger.warning("Component spin extraction failed.")
+ logger.warning(output_sample.keys())
- elif all(key in output_sample.keys() for key in spin_conversion_parameters)\
- and isinstance(output_sample, pd.DataFrame):
- logger.debug('Extracting component spins.')
- new_spin_parameters =\
- ['spin_1x', 'spin_1y', 'spin_1z', 'spin_2x', 'spin_2y', 'spin_2z']
- new_spins =\
- {name: np.zeros(len(output_sample)) for name in new_spin_parameters}
-
- for ii in range(len(output_sample)):
- new_spins['iota'], new_spins['spin_1x'][ii],\
- new_spins['spin_1y'][ii], new_spins['spin_1z'][ii], \
- new_spins['spin_2x'][ii], new_spins['spin_2y'][ii],\
- new_spins['spin_2z'][ii] =\
- lalsim.SimInspiralTransformPrecessingNewInitialConditions(
- output_sample['iota'][ii], output_sample['phi_jl'][ii],
- output_sample['tilt_1'][ii], output_sample['tilt_2'][ii],
- output_sample['phi_12'][ii], output_sample['a_1'][ii],
- output_sample['a_2'][ii],
- output_sample['mass_1'][ii] * solar_mass,
- output_sample['mass_2'][ii] * solar_mass,
- output_sample['reference_frequency'][ii],
- output_sample['phase'][ii])
-
- for name in new_spin_parameters:
- output_sample[name] = new_spins[name]
+ return output_sample
- output_sample['phi_1'] =\
- np.arctan(output_sample['spin_1y'] / output_sample['spin_1x'])
- output_sample['phi_2'] =\
- np.arctan(output_sample['spin_2y'] / output_sample['spin_2x'])
- else:
- logger.warning("Component spin extraction failed.")
+def generate_tidal_parameters(sample):
+ """
+ Generate all tidal parameters
+
+ lambda_tilde, delta_lambda
+
+ Parameters
+ ----------
+ sample: dict, pandas.DataFrame
+ Should contain lambda_1, lambda_2
+
+ Returns
+ -------
+ output_sample: dict, pandas.DataFrame
+ Updated sample
+ """
+ output_sample = sample.copy()
+
+ output_sample['lambda_tilde'] =\
+ lambda_1_lambda_2_to_lambda_tilde(
+ output_sample['lambda_1'], output_sample['lambda_2'],
+ output_sample['mass_1'], output_sample['mass_2'])
+ output_sample['delta_lambda'] = \
+ lambda_1_lambda_2_to_delta_lambda(
+ output_sample['lambda_1'], output_sample['lambda_2'],
+ output_sample['mass_1'], output_sample['mass_2'])
+
+ return output_sample
+
+
+def generate_source_frame_parameters(sample):
+ """
+ Generate source frame masses along with redshifts and comoving distance.
+
+ Parameters
+ ----------
+ sample: dict, pandas.DataFrame
+
+ Returns
+ -------
+ output_sample: dict, pandas.DataFrame
+ """
+ output_sample = sample.copy()
+
+ output_sample['redshift'] =\
+ luminosity_distance_to_redshift(output_sample['luminosity_distance'])
+ output_sample['comoving_distance'] =\
+ redshift_to_comoving_distance(output_sample['redshift'])
+
+ for key in ['mass_1', 'mass_2', 'chirp_mass', 'total_mass']:
+ if key in output_sample:
+ output_sample['{}_source'.format(key)] =\
+ output_sample[key] / (1 + output_sample['redshift'])
return output_sample
@@ -674,7 +904,6 @@ def compute_snrs(sample, likelihood):
if isinstance(temp_sample, dict):
temp = dict()
for key in likelihood.waveform_generator.parameters.keys():
- # likelihood.waveform_generator.parameters[key] = temp_sample[key]
temp[key] = temp_sample[key]
signal_polarizations =\
likelihood.waveform_generator.frequency_domain_strain(temp)
@@ -696,8 +925,6 @@ def compute_snrs(sample, likelihood):
temp = dict()
for key in set(temp_sample.keys()).intersection(
likelihood.waveform_generator.parameters.keys()):
- # likelihood.waveform_generator.parameters[key] =\
- # temp_sample[key][ii]
temp[key] = temp_sample[key][ii]
signal_polarizations =\
likelihood.waveform_generator.frequency_domain_strain(temp)
diff --git a/bilby/gw/prior_files/binary_neutron_stars.prior b/bilby/gw/prior_files/binary_neutron_stars.prior
index c30fa00bca63d5d64675c361f5e64b42fe45409d..1bc4d485c1cfdf89ddee3ab4b7b33d6cabd80654 100644
--- a/bilby/gw/prior_files/binary_neutron_stars.prior
+++ b/bilby/gw/prior_files/binary_neutron_stars.prior
@@ -8,8 +8,8 @@ mass_2 = Uniform(name='mass_2', minimum=1, maximum=2, unit='$M_{\\odot}$')
# total_mass = Uniform(name='total_mass', minimum=2, maximum=4, unit='$M_{\\odot}$')
# mass_ratio = Uniform(name='mass_ratio', minimum=0.5, maximum=1)
# symmetric_mass_ratio = Uniform(name='symmetric_mass_ratio', minimum=0.22, maximum=0.25)
-a_1 = Uniform(name='a_1', minimum= -0.05, maximum= 0.05)
-a_2 = Uniform(name='a_2', minimum= -0.05, maximum= 0.05)
+chi_1 = bilby.gw.prior.AlignedSpin(a_prior=Uniform(0, 0.05), z_prior=Uniform(-1, 1), name='chi_1', latex_label='$\\chi_1$')
+chi_2 = bilby.gw.prior.AlignedSpin(a_prior=Uniform(0, 0.05), z_prior=Uniform(-1, 1), name='chi_2', latex_label='$\\chi_2$')
luminosity_distance = bilby.gw.prior.UniformComovingVolume(name='luminosity_distance', minimum=10, maximum=500, unit='Mpc')
dec = Cosine(name='dec')
ra = Uniform(name='ra', minimum=0, maximum=2 * np.pi)
diff --git a/bilby/gw/source.py b/bilby/gw/source.py
index 6ea32a002b47d9db62b0e1faf2ed77ca0af6662b..75bf22febf755b10732b21d057d7a738826f44bc 100644
--- a/bilby/gw/source.py
+++ b/bilby/gw/source.py
@@ -255,7 +255,7 @@ def supernova_pca_model(
def lal_binary_neutron_star(
- frequency_array, mass_1, mass_2, luminosity_distance, a_1, a_2,
+ frequency_array, mass_1, mass_2, luminosity_distance, chi_1, chi_2,
iota, phase, lambda_1, lambda_2, ra, dec, geocent_time, psi, **kwargs):
""" A Binary Neutron Star waveform model using lalsimulation
@@ -269,10 +269,10 @@ def lal_binary_neutron_star(
The mass of the lighter object in solar masses
luminosity_distance: float
The luminosity distance in megaparsec
- a_1: float
- Dimensionless spin magnitude
- a_2: float
- Dimensionless secondary spin magnitude
+ chi_1: float
+ Dimensionless aligned spin
+ chi_2: float
+ Dimensionless aligned spin
iota: float
Orbital inclination
phase: float
@@ -314,10 +314,10 @@ def lal_binary_neutron_star(
spin_1x = 0
spin_1y = 0
- spin_1z = a_1
+ spin_1z = chi_1
spin_2x = 0
spin_2y = 0
- spin_2z = a_2
+ spin_2z = chi_2
longitude_ascending_nodes = 0.0
eccentricity = 0.0
diff --git a/examples/injection_examples/basic_tutorial.py b/examples/injection_examples/basic_tutorial.py
index 1e90bcddcda848fb378792bb6ff7f788ba10a41f..2f3766f8adfd656b63e22b706fe80b5fc54d3bf9 100644
--- a/examples/injection_examples/basic_tutorial.py
+++ b/examples/injection_examples/basic_tutorial.py
@@ -1,9 +1,12 @@
#!/usr/bin/env python
"""
-Tutorial to demonstrate running parameter estimation on a reduced parameter space for an injected signal.
+Tutorial to demonstrate running parameter estimation on a reduced parameter
+space for an injected signal.
-This example estimates the masses using a uniform prior in both component masses and distance using a uniform in
-comoving volume prior on luminosity distance between luminosity distances of 100Mpc and 5Gpc, the cosmology is WMAP7.
+This example estimates the masses using a uniform prior in both component masses
+and distance using a uniform in
+comoving volume prior on luminosity distance between luminosity distances of
+100Mpc and 5Gpc, the cosmology is WMAP7.
"""
from __future__ import division, print_function
@@ -11,7 +14,8 @@ import numpy as np
import bilby
-# Set the duration and sampling frequency of the data segment that we're going to inject the signal into
+# Set the duration and sampling frequency of the data segment that we're going
+# to inject the signal into
duration = 4.
sampling_frequency = 2048.
@@ -24,12 +28,14 @@ bilby.core.utils.setup_logger(outdir=outdir, label=label)
# Set up a random seed for result reproducibility. This is optional!
np.random.seed(88170235)
-# We are going to inject a binary black hole waveform. We first establish a dictionary of parameters that
-# includes all of the different waveform parameters, including masses of the two black holes (mass_1, mass_2),
+# We are going to inject a binary black hole waveform. We first establish a
+# dictionary of parameters that includes all of the different waveform
+# parameters, including masses of the two black holes (mass_1, mass_2),
# spins of both black holes (a, tilt, phi), etc.
-injection_parameters = dict(mass_1=36., mass_2=29., a_1=0.4, a_2=0.3, tilt_1=0.5, tilt_2=1.0, phi_12=1.7, phi_jl=0.3,
- luminosity_distance=2000., iota=0.4, psi=2.659, phase=1.3, geocent_time=1126259642.413,
- ra=1.375, dec=-1.2108)
+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., iota=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',
@@ -51,29 +57,37 @@ ifos.set_strain_data_from_power_spectral_densities(
ifos.inject_signal(waveform_generator=waveform_generator,
parameters=injection_parameters)
-# Set up prior, which is a dictionary
-# By default we will sample all terms in the signal models. However, this will take a long time for the calculation,
-# so for this example we will set almost all of the priors to be equall to their injected values. This implies the
-# prior is a delta function at the true, injected value. In reality, the sampler implementation is smart enough to
-# not sample any parameter that has a delta-function prior.
-# The above list does *not* include mass_1, mass_2, iota and luminosity_distance, which means those are the parameters
-# that will be included in the sampler. If we do nothing, then the default priors get used.
+# Set up PriorSet, which inherits from dict
+# By default we will sample all terms in the signal models. However, this will
+# take a long time for the calculation, so for this example we will set almost
+# all of the priors to be equal to their injected values. This implies the
+# prior is a delta function at the true, injected value. In reality, the
+# sampler implementation is smart enough to not sample any parameter that has a
+# delta-function prior. The above list does *not* include mass_1, mass_2, iota
+# and luminosity_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.BBHPriorSet()
-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', 'psi', 'ra', 'dec', 'geocent_time', 'phase']:
+# 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', 'psi', 'ra',
+ 'dec', 'geocent_time', 'phase']:
priors[key] = injection_parameters[key]
-# Initialise the likelihood by passing in the interferometer data (IFOs) and the waveoform generator
-likelihood = bilby.gw.GravitationalWaveTransient(interferometers=ifos, waveform_generator=waveform_generator,
- time_marginalization=False, phase_marginalization=False,
- distance_marginalization=False, prior=priors)
+# Initialise the likelihood by passing in the interferometer data (IFOs)
+# and the waveoform generator
+likelihood = bilby.gw.GravitationalWaveTransient(
+ interferometers=ifos, waveform_generator=waveform_generator,
+ time_marginalization=False, phase_marginalization=False,
+ distance_marginalization=False, prior=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=1000,
- injection_parameters=injection_parameters, outdir=outdir, label=label)
+result = bilby.run_sampler(
+ likelihood=likelihood, priors=priors, sampler='nestle', npoints=1000,
+ injection_parameters=injection_parameters, outdir=outdir, label=label,
+ conversion_function=bilby.gw.conversion.generate_all_bbh_parameters)
+
# make some plots of the outputs
result.plot_corner()
diff --git a/examples/injection_examples/binary_neutron_star_example.py b/examples/injection_examples/binary_neutron_star_example.py
index 4f826a060037129c891286c28e2a022d626a6edd..76cf0c1ec9d27f79099f2798df6847a963070ee1 100644
--- a/examples/injection_examples/binary_neutron_star_example.py
+++ b/examples/injection_examples/binary_neutron_star_example.py
@@ -25,9 +25,9 @@ np.random.seed(88170235)
# We are going to inject a binary neutron star waveform. We first establish a
# dictionary of parameters that includes all of the different waveform
# parameters, including masses of the two black holes (mass_1, mass_2),
-# spins of both black holes (a_1,a_2) , etc.
+# aligned spins of both black holes (chi_1, chi_2), etc.
injection_parameters = dict(
- mass_1=1.5, mass_2=1.3, a_1=0.0, a_2=0.0, luminosity_distance=50.,
+ mass_1=1.5, mass_2=1.3, chi_1=0.02, chi_2=0.02, luminosity_distance=50.,
iota=0.4, psi=2.659, phase=1.3, geocent_time=1126259642.413,
ra=1.375, dec=-1.2108, lambda_1=400, lambda_2=450)
@@ -46,6 +46,7 @@ waveform_arguments = dict(waveform_approximant='TaylorF2',
waveform_generator = bilby.gw.WaveformGenerator(
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)
# Set up interferometers. In this case we'll use three interferometers
@@ -60,11 +61,25 @@ interferometers.set_strain_data_from_power_spectral_densities(
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
+# delta_lambda rather than mass_1, mass_2, lambda_1, and lambda_2.
priors = bilby.gw.prior.BNSPriorSet()
-for key in ['a_1', 'a_2', 'psi', 'geocent_time', 'ra', 'dec',
+for key in ['psi', 'geocent_time', 'ra', 'dec', 'chi_1', 'chi_2',
'iota', '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['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 waveoform generator
likelihood = bilby.gw.GravitationalWaveTransient(
@@ -75,7 +90,8 @@ likelihood = bilby.gw.GravitationalWaveTransient(
# Run sampler. In this case we're going to use the `nestle` sampler
result = bilby.run_sampler(
likelihood=likelihood, priors=priors, sampler='nestle', npoints=1000,
- injection_parameters=injection_parameters, outdir=outdir, label=label)
+ injection_parameters=injection_parameters, outdir=outdir, label=label,
+ conversion_function=bilby.gw.conversion.generate_all_bns_parameters)
result.plot_corner()
diff --git a/test/conversion_test.py b/test/conversion_test.py
index b578237877fac820ba2d7952f68048a8b2b11f9b..3d2d698582db814f0b58bf34e32b6f7d75a27132 100644
--- a/test/conversion_test.py
+++ b/test/conversion_test.py
@@ -97,8 +97,20 @@ class TestBasicConversions(unittest.TestCase):
self.assertTrue(all([abs(self.lambda_1 - lambda_1) < 1e-5,
abs(self.lambda_2 - lambda_2) < 1e-5]))
+ def test_lambda_1_lambda_2_to_lambda_tilde(self):
+ lambda_tilde = \
+ conversion.lambda_1_lambda_2_to_lambda_tilde(
+ self.lambda_1, self.lambda_2, self.mass_1, self.mass_2)
+ self.assertTrue((self.lambda_tilde - lambda_tilde) < 1e-5)
-class TestConvertToLALBBHParams(unittest.TestCase):
+ def test_lambda_1_lambda_2_to_delta_lambda(self):
+ delta_lambda = \
+ conversion.lambda_1_lambda_2_to_lambda_tilde(
+ self.lambda_1, self.lambda_2, self.mass_1, self.mass_2)
+ self.assertTrue((self.delta_lambda - delta_lambda) < 1e-5)
+
+
+class TestConvertToLALParams(unittest.TestCase):
def setUp(self):
self.search_keys = []
@@ -110,7 +122,7 @@ class TestConvertToLALBBHParams(unittest.TestCase):
del self.parameters
del self.remove
- def test_cos_angle_to_angle_conversion(self):
+ def test_bbh_cos_angle_to_angle_conversion(self):
with mock.patch('numpy.arccos') as m:
m.return_value = 42
self.parameters['cos_tilt_1'] = 1
@@ -119,6 +131,16 @@ class TestConvertToLALBBHParams(unittest.TestCase):
self.parameters)
self.assertDictEqual(self.parameters, dict(tilt_1=42, cos_tilt_1=1))
+ def test_bns_cos_angle_to_angle_conversion(self):
+ with mock.patch('numpy.arccos') as m:
+ m.return_value = 42
+ self.parameters['cos_tilt_1'] = 1
+ self.parameters, _ = \
+ conversion.convert_to_lal_binary_neutron_star_parameters(
+ self.parameters)
+ self.assertDictEqual(self.parameters, dict(
+ tilt_1=42, cos_tilt_1=1, lambda_1=0, lambda_2=0))
+
if __name__ == '__main__':
unittest.main()