Merge branch...

Merge branch '202-bns-ability-to-sample-in-and-convert-to-tilde-lambda-and-delta-tilde-lambda' into 'master'

Resolve "BNS: ability to sample in and convert to \tilde{Lambda} and \delta\tilde{\Lambda}"

Closes #172 and #202

See merge request Monash/bilby!222

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

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
 

bilby/gw/conversion.py

 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)

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)

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

examples/injection_examples/basic_tutorial.py

 #!/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()

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()
 

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()