Revert "Merge branch '176-add-neutron-star-merger-to-examples' into 'master'"

This reverts merge request !171

examples/injection_examples/binary_neutron_star_example.py

-#!/bin/python
-"""
-Tutorial to demonstrate running parameter estimation on a binary neutron star
-system taking into account 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
-"""
-
-from __future__ import division, print_function
-
-import numpy as np
-
-import tupak
-
-# Specify the output directory and the name of the simulation.
-outdir = 'outdir'
-label = 'bns_example'
-tupak.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 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. 
-injection_parameters = dict(
-    mass_1=1.5, mass_2=1.3, a_1=0.0, a_2=0.0, luminosity_distance=50.,
-    iota=0.4, psi=2.659, phase=1.3, geocent_time=1126259642.413,
-    ra=1.375, dec=-1.2108, lambda1=400, lambda2=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 = 8
-sampling_frequency = 2*1570.
-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='TaylorF2',
-                          reference_frequency=50., minimum_frequency=40.0)
-
-# Create the waveform_generator using a LAL Binary Neutron Star source function
-waveform_generator = tupak.gw.WaveformGenerator(
-    duration=duration, sampling_frequency=sampling_frequency,
-    frequency_domain_source_model=tupak.gw.source.lal_binary_neutron_star,
-    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 = tupak.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)
-
-priors = tupak.gw.prior.BNSPriorSet()
-for key in ['a_1', 'a_2', 'psi', 'geocent_time', 'ra', 'dec',
-            'iota', 'luminosity_distance', 'phase']:
-    priors[key] = injection_parameters[key]
-    
-# Initialise the likelihood by passing in the interferometer data (IFOs)
-# and the waveoform generator
-likelihood = tupak.gw.GravitationalWaveTransient(
-    interferometers=interferometers, 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 `nestle` sampler
-result = tupak.run_sampler(
-    likelihood=likelihood, priors=priors, sampler='nestle', npoints=1000,
-    injection_parameters=injection_parameters, outdir=outdir, label=label)
-
-result.plot_corner()
-

test/conversion_tests.py

@@ -8,35 +8,14 @@ import numpy as np
 class TestBasicConversions(unittest.TestCase):
 
     def setUp(self):
-        self.mass_1 = 1.4
-        self.mass_2 = 1.3
-        self.mass_ratio = 13/14
-        self.total_mass = 2.7
-        self.chirp_mass = (1.4 * 1.3)**0.6 / 2.7**0.2
-        self.symmetric_mass_ratio = (1.4 * 1.3) / 2.7**2
+        self.mass_1 = 20
+        self.mass_2 = 10
+        self.mass_ratio = 0.5
+        self.total_mass = 30
+        self.chirp_mass = 200**0.6 / 30**0.2
+        self.symmetric_mass_ratio = 2/9
         self.cos_angle = -1
         self.angle = np.pi
-        self.lambda_1 = 300
-        self.lambda_2 = 300 * (14 / 13)**5
-        self.lambda_tilde = 8 / 13 * (
-            (1 + 7 * self.symmetric_mass_ratio
-             - 31 * self.symmetric_mass_ratio**2)
-            * (self.lambda_1 + self.lambda_2)
-            + (1 - 4 * self.symmetric_mass_ratio)**0.5
-            * (1 + 9 * self.symmetric_mass_ratio
-               - 11 * self.symmetric_mass_ratio**2)
-            * (self.lambda_1 - self.lambda_2)
-        )
-        self.delta_lambda = 1 / 2 * (
-                (1 - 4 * self.symmetric_mass_ratio)**0.5
-                * (1 - 13272 / 1319 * self.symmetric_mass_ratio
-                   + 8944 / 1319 * self.symmetric_mass_ratio**2)
-                * (self.lambda_1 + self.lambda_2)
-                + (1 - 15910 / 1319 * self.symmetric_mass_ratio
-                   + 32850 / 1319 * self.symmetric_mass_ratio**2
-                   + 3380 / 1319 * self.symmetric_mass_ratio**3)
-                * (self.lambda_1 - self.lambda_2)
-        )
 
     def tearDown(self):
         del self.mass_1
@@ -48,8 +27,7 @@ class TestBasicConversions(unittest.TestCase):
 
     def test_total_mass_and_mass_ratio_to_component_masses(self):
         mass_1, mass_2 = tupak.gw.conversion.total_mass_and_mass_ratio_to_component_masses(self.mass_ratio, self.total_mass)
-        self.assertTrue(all([abs(mass_1 - self.mass_1) < 1e-5,
-                             abs(mass_2 - self.mass_2) < 1e-5]))
+        self.assertTupleEqual((mass_1, mass_2), (self.mass_1, self.mass_2))
 
     def test_symmetric_mass_ratio_to_mass_ratio(self):
         mass_ratio = tupak.gw.conversion.symmetric_mass_ratio_to_mass_ratio(self.symmetric_mass_ratio)
@@ -83,20 +61,6 @@ class TestBasicConversions(unittest.TestCase):
         mass_ratio = tupak.gw.conversion.mass_1_and_chirp_mass_to_mass_ratio(self.mass_1, self.chirp_mass)
         self.assertAlmostEqual(self.mass_ratio, mass_ratio)
 
-    def test_lambda_tilde_to_lambda_1_lambda_2(self):
-        lambda_1, lambda_2 =\
-            tupak.gw.conversion.lambda_tilde_to_lambda_1_lambda_2(
-                self.lambda_tilde, self.mass_1, self.mass_2)
-        self.assertTrue(all([abs(self.lambda_1 - lambda_1) < 1e-5,
-                             abs(self.lambda_2 - lambda_2) < 1e-5]))
-
-    def test_lambda_tilde_delta_lambda_to_lambda_1_lambda_2(self):
-        lambda_1, lambda_2 =\
-            tupak.gw.conversion.lambda_tilde_delta_lambda_to_lambda_1_lambda_2(
-                self.lambda_tilde, self.delta_lambda, self.mass_1, self.mass_2)
-        self.assertTrue(all([abs(self.lambda_1 - lambda_1) < 1e-5,
-                             abs(self.lambda_2 - lambda_2) < 1e-5]))
-
 
 class TestConvertToLALBBHParams(unittest.TestCase):
 

tupak/gw/conversion.py

@@ -171,71 +171,6 @@ def convert_to_lal_binary_black_hole_parameters(parameters, search_keys, remove=
     return converted_parameters, added_keys
 
 
-def convert_to_lal_binary_neutron_star_parameters(parameters, search_keys, remove=True):
-    """
-    Convert parameters we have into parameters we need.
-
-    This is defined by the parameters of tupak.source.lal_binary_black_hole()
-
-
-    Mass: mass_1, mass_2
-    Spin: a_1, a_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.
-    The keys in added_keys should be popped after evaluating the waveform.
-
-    Parameters
-    ----------
-    parameters: dict
-        dictionary of parameter values to convert into the required parameters
-    search_keys: list
-        parameters which are needed for the waveform generation
-    remove: bool, optional
-        Whether or not to remove the extra key, necessary for sampling, default=True.
-
-    Return
-    ------
-    converted_parameters: dict
-        dict of the required parameters
-    added_keys: list
-        keys which are added to parameters during function call
-    """
-    converted_parameters = parameters.copy()
-    converted_parameters, added_keys = convert_to_lal_binary_black_hole_parameters(
-        converted_parameters, search_keys, remove=remove)
-
-    if 'lambda_1' not in search_keys and 'lambda_2' not in search_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(
-                    converted_parameters['lambda_tilde'], parameters['delta_lambda'],
-                    converted_parameters['mass_1'], converted_parameters['mass_2'])
-            added_keys.append('lambda_1')
-            added_keys.append('lambda_2')
-        elif 'lambda_tilde' in converted_parameters.keys():
-            converted_parameters['lambda_1'], converted_parameters['lambda_2'] =\
-                lambda_tilde_to_lambda_1_lambda_2(
-                    converted_parameters['lambda_tilde'],
-                    converted_parameters['mass_1'], converted_parameters['mass_2'])
-            added_keys.append('lambda_1')
-            added_keys.append('lambda_2')
-    if 'lambda_2' not in converted_parameters.keys():
-        converted_parameters['lambda_2'] =\
-            converted_parameters['lambda_1']\
-            * converted_parameters['mass_1']**5\
-            / converted_parameters['mass_2']**5
-        added_keys.append('lambda_2')
-    elif converted_parameters['lambda_2'] is None:
-        converted_parameters['lambda_2'] =\
-            converted_parameters['lambda_1']\
-            * converted_parameters['mass_1']**5\
-            / converted_parameters['mass_2']**5
-        added_keys.append('lambda_2')
-
-    return converted_parameters, added_keys
-
-
 def total_mass_and_mass_ratio_to_component_masses(mass_ratio, total_mass):
     """
     Convert total mass and mass ratio of a binary to its component masses.
@@ -425,85 +360,6 @@ def mass_1_and_chirp_mass_to_mass_ratio(mass_1, chirp_mass):
     return mass_ratio
 
 
-def lambda_tilde_delta_lambda_to_lambda_1_lambda_2(
-        lambda_tilde, delta_lambda, mass_1, mass_2):
-    """
-    Convert from dominant tidal terms to individual tidal parameters.
-
-    See, e.g., Wade et al., https://arxiv.org/pdf/1402.5156.pdf.
-
-    Parameters
-    ----------
-    lambda_tilde: float
-        Dominant tidal term.
-    delta_lambda: float
-        Secondary tidal term.
-    mass_1: float
-        Mass of more massive neutron star.
-    mass_2: float
-        Mass of less massive neutron star.
-
-    Return
-    ------
-    lambda_1: float
-        Tidal parameter of more massive neutron star.
-    lambda_2: float
-        Tidal parameter of less massive neutron star.
-    """
-    eta = component_masses_to_symmetric_mass_ratio(mass_1, mass_2)
-    coefficient_1 = (1 + 7 * eta - 31 * eta**2)
-    coefficient_2 = (1 - 4 * eta)**0.5 * (1 + 9 * eta - 11 * eta**2)
-    coefficient_3 = (1 - 4 * eta)**0.5\
-                    * (1 - 13272 / 1319 * eta + 8944 / 1319 * eta**2)
-    coefficient_4 = (1 - 15910 / 1319 * eta + 32850 / 1319 * eta**2
-                     + 3380 / 1319 * eta**3)
-    lambda_1 =\
-        (13 * lambda_tilde / 8 * (coefficient_3 - coefficient_4)
-         - 2 * delta_lambda * (coefficient_1 - coefficient_2))\
-        / ((coefficient_1 + coefficient_2) * (coefficient_3 - coefficient_4)
-           - (coefficient_1 - coefficient_2) * (coefficient_3 + coefficient_4))
-    lambda_2 =\
-        (13 * lambda_tilde / 8 * (coefficient_3 + coefficient_4)
-         - 2 * delta_lambda * (coefficient_1 + coefficient_2)) \
-        / ((coefficient_1 - coefficient_2) * (coefficient_3 + coefficient_4)
-           - (coefficient_1 + coefficient_2) * (coefficient_3 - coefficient_4))
-    return lambda_1, lambda_2
-
-
-def lambda_tilde_to_lambda_1_lambda_2(
-        lambda_tilde, mass_1, mass_2):
-    """
-    Convert from dominant tidal term to individual tidal parameters
-    assuming lambda_1 * mass_1**5 = lambda_2 * mass_2**5.
-
-    See, e.g., Wade et al., https://arxiv.org/pdf/1402.5156.pdf.
-
-    Parameters
-    ----------
-    lambda_tilde: float
-        Dominant tidal term.
-    mass_1: float
-        Mass of more massive neutron star.
-    mass_2: float
-        Mass of less massive neutron star.
-
-    Return
-    ------
-    lambda_1: float
-        Tidal parameter of more massive neutron star.
-    lambda_2: float
-        Tidal parameter of less massive neutron star.
-    """
-    eta = component_masses_to_symmetric_mass_ratio(mass_1, mass_2)
-    q = mass_2 / mass_1
-    lambda_1 = 13 / 8 * lambda_tilde / (
-            (1 + 7 * eta - 31 * eta**2) * (1 + q**-5)
-            + (1 - 4 * eta)**0.5 * (1 + 9 * eta - 11 * eta**2) * (1 - q**-5)
-        )
-    lambda_2 = lambda_1 / q**5
-    return lambda_1, lambda_2
-
-
 def generate_all_bbh_parameters(sample, likelihood=None, priors=None):
     """
     From either a single sample or a set of samples fill in all missing BBH parameters, in place.

tupak/gw/prior.py

@@ -96,45 +96,6 @@ class BBHPriorSet(PriorSet):
         return redundant
 
 
-class BNSPriorSet(PriorSet):
-
-    def __init__(self, dictionary=None, filename=None):
-        """ Initialises a Prior set for Binary Neutron Stars
-
-        Parameters
-        ----------
-        dictionary: dict, optional
-            See superclass
-        filename: str, optional
-            See superclass
-        """
-        if dictionary is None and filename is None:
-            filename = os.path.join(os.path.dirname(__file__), 'prior_files', 'binary_neutron_stars.prior')
-            logger.info('No prior given, using default BNS priors in {}.'.format(filename))
-        elif filename is not None:
-            if not os.path.isfile(filename):
-                filename = os.path.join(os.path.dirname(__file__), 'prior_files', filename)
-        PriorSet.__init__(self, dictionary=dictionary, filename=filename)
-
-    def test_redundancy(self, key):
-        bbh_redundancy = BBHPriorSet().test_redundancy(key)
-        if bbh_redundancy:
-            return True
-        redundant = False
-
-        tidal_parameters =\
-            {'lambda_1', 'lambda_2', 'lambda_tilde', 'delta_lambda'}
-
-        if key in tidal_parameters:
-            if len(tidal_parameters.intersection(self)) > 2:
-                redundant = True
-                logger.warning('{} in prior. This may lead to unexpected behaviour.'.format(
-                    tidal_parameters.intersection(self)))
-            elif len(tidal_parameters.intersection(self)) == 2:
-                redundant = True
-        return redundant
-
-
 Prior._default_latex_labels = {
     'mass_1': '$m_1$',
     'mass_2': '$m_2$',
@@ -157,11 +118,7 @@ Prior._default_latex_labels = {
     'cos_iota': '$\cos\iota$',
     'psi': '$\psi$',
     'phase': '$\phi$',
-    'geocent_time': '$t_c$',
-    'lambda_1': '$\\Lambda_1$',
-    'lambda_2': '$\\Lambda_2$',
-    'lambda_tilde': '$\\tilde{\\Lambda}$',
-    'delta_lambda': '$\\delta\\Lambda$'}
+    'geocent_time': '$t_c$'}
 
 
 class CalibrationPriorSet(PriorSet):

tupak/gw/prior_files/binary_neutron_stars.prior

-# These are the default priors we use for BNS systems.
-# Note that you may wish to use more specific mass and distance parameters.
-# These commands are all known to tupak.gw.prior.
-# Lines beginning "#" are ignored.
-mass_1 = Uniform(name='mass_1', minimum=1, maximum=2)
-mass_2 = Uniform(name='mass_2', minimum=1, maximum=2)
-# chirp_mass = Uniform(name='chirp_mass', minimum=0.87, maximum=1.74)
-# total_mass =  Uniform(name='total_mass', minimum=2, maximum=4)
-# 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)
-luminosity_distance =  tupak.gw.prior.UniformComovingVolume(name='luminosity_distance', minimum=10, maximum=500)
-dec =  Cosine(name='dec')
-ra =  Uniform(name='ra', minimum=0, maximum=2 * np.pi)
-iota =  Sine(name='iota')
-# cos_iota =  Uniform(name='cos_iota', minimum=-1, maximum=1)
-psi =  Uniform(name='psi', minimum=0, maximum=np.pi)
-phase =  Uniform(name='phase', minimum=0, maximum=2 * np.pi)
-lambda_1 = Uniform(name='lambda_1', minimum=0, maximum=3000 )
-lambda_2 = Uniform(name='lambda_2', minimum=0, maximum=3000 )
-# lambda_tilde = Uniform(name='lambda_tilde', minimum=0, maximum=5000)
-# delta_lambda = Uniform(name='delta_lambda', minimum=-5000, maximum=5000)

tupak/gw/source.py

@@ -7,7 +7,6 @@ from tupak.core import utils
 
 try:
     import lalsimulation as lalsim
-    import lal
 except ImportError:
     logger.warning("You do not have lalsuite installed currently. You will "
                    " not be able to use some of the prebuilt functions.")
@@ -252,96 +251,3 @@ def supernova_pca_model(
                          pc_coeff4 * pc4 + pc_coeff5 * pc5)
 
     return {'plus': h_plus, 'cross': h_cross}
-
-def lal_binary_neutron_star(
-        frequency_array, mass_1, mass_2, luminosity_distance, a_1, a_2,
-        iota, phase, lambda_1, lambda_2, ra, dec, geocent_time, psi, **kwargs):
-    """ A Binary Neutron Star waveform model using lalsimulation
-
-    Parameters
-    ----------
-    frequency_array: array_like
-        The frequencies at which we want to calculate the strain
-    mass_1: float
-        The mass of the heavier object in solar masses
-    mass_2: float
-        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. 
-    iota: float
-        Orbital inclination
-    phase: float
-        The phase at coalescence
-    ra: float
-        The right ascension of the binary
-    dec: float
-        The declination of the object
-    geocent_time: float
-        The time at coalescence
-    psi: float
-        Orbital polarisation
-    lambda_1: float
-        Dimensionless tidal deformability of mass_1
-    lambda_2: float
-        Dimensionless tidal deformability of mass_2
-        
-    kwargs: dict
-        Optional keyword arguments
-
-    Returns
-    -------
-    dict: A dictionary with the plus and cross polarisation strain modes
-    """
-
-    waveform_kwargs = dict(waveform_approximant='TaylorF2', reference_frequency=50.0,
-                           minimum_frequency=20.0)
-    waveform_kwargs.update(kwargs)
-    waveform_approximant = waveform_kwargs['waveform_approximant']
-    reference_frequency = waveform_kwargs['reference_frequency']
-    minimum_frequency = waveform_kwargs['minimum_frequency']
-
-    if mass_2 > mass_1:
-        return None
-
-    luminosity_distance = luminosity_distance * 1e6 * utils.parsec
-    mass_1 = mass_1 * utils.solar_mass
-    mass_2 = mass_2 * utils.solar_mass
-
-    spin_1x = 0
-    spin_1y = 0
-    spin_1z = a_1
-    spin_2x = 0
-    spin_2y = 0
-    spin_2z = a_2
-
-    longitude_ascending_nodes = 0.0
-    eccentricity = 0.0
-    mean_per_ano = 0.0
-    
-    waveform_dictionary = lal.CreateDict()
-    lalsim.SimInspiralWaveformParamsInsertTidalLambda1(waveform_dictionary, lambda_1)
-    lalsim.SimInspiralWaveformParamsInsertTidalLambda2(waveform_dictionary, lambda_2)
-
-    approximant = lalsim.GetApproximantFromString(waveform_approximant)
-
-    maximum_frequency = frequency_array[-1]
-    delta_frequency = frequency_array[1] - frequency_array[0]
-
-    hplus, hcross = lalsim.SimInspiralChooseFDWaveform(
-        mass_1, mass_2, spin_1x, spin_1y, spin_1z, spin_2x, spin_2y,
-        spin_2z, luminosity_distance, iota, phase,
-        longitude_ascending_nodes, eccentricity, mean_per_ano, delta_frequency,
-        minimum_frequency, maximum_frequency, reference_frequency,
-        waveform_dictionary, approximant)
-
-    h_plus = hplus.data.data
-    h_cross = hcross.data.data
-    
-    h_plus = h_plus[:len(frequency_array)]
-    h_cross = h_cross[:len(frequency_array)]
-
-    return {'plus': h_plus, 'cross': h_cross}