diff --git a/CHANGELOG.md b/CHANGELOG.md index 83db3fefce47a6e937343425ef63c3f67bf48d3a..00db6a52e5d281f4f58fdca7adfc4025c90fa838 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,9 +4,10 @@ ### Added - `emcee` now writes all progress to disk and can resume from a previous run. -- ### Changed +- Cosmology generalised, users can now specify the cosmology used, default is astropy Planck15 +- UniformComovingVolume prior *requires* the name to be one of "luminosity_distance", "comoving_distance", "redshift" - Time/frequency array generation/conversion improved. We now impose `duration` is an integer multiple of `sampling_frequency`. Converting back and forth between time/frequency arrays now works for all valid arrays. - Updates the bilby.core.utils constants to match those of Astropy v3.0.4 diff --git a/bilby/core/prior.py b/bilby/core/prior.py index 8bc6e4e5daf9324868b700d96ea03443dd9f5218..228148d6c11dae4671d6b424c09e7959a0fe469b 100644 --- a/bilby/core/prior.py +++ b/bilby/core/prior.py @@ -512,19 +512,19 @@ class Prior(object): @property def minimum(self): - return self.__minimum + return self._minimum @minimum.setter def minimum(self, minimum): - self.__minimum = minimum + self._minimum = minimum @property def maximum(self): - return self.__maximum + return self._maximum @maximum.setter def maximum(self, maximum): - self.__maximum = maximum + self._maximum = maximum @property def __default_latex_label(self): @@ -1587,7 +1587,7 @@ class Interped(Prior): Prior.__init__(self, name=name, latex_label=latex_label, unit=unit, minimum=np.nanmax(np.array((min(xx), minimum))), maximum=np.nanmin(np.array((max(xx), maximum)))) - self.__initialize_attributes() + self.__update_instance() def __eq__(self, other): if self.__class__ != other.__class__: @@ -1632,12 +1632,12 @@ class Interped(Prior): float: Minimum of the prior distribution """ - return self.__minimum + return self._minimum @minimum.setter def minimum(self, minimum): - self.__minimum = minimum - if '_Interped__maximum' in self.__dict__ and self.__maximum < np.inf: + self._minimum = minimum + if '_maximum' in self.__dict__ and self._maximum < np.inf: self.__update_instance() @property @@ -1651,12 +1651,12 @@ class Interped(Prior): float: Maximum of the prior distribution """ - return self.__maximum + return self._maximum @maximum.setter def maximum(self, maximum): - self.__maximum = maximum - if '_Interped__minimum' in self.__dict__ and self.__minimum < np.inf: + self._maximum = maximum + if '_minimum' in self.__dict__ and self._minimum < np.inf: self.__update_instance() def __update_instance(self): diff --git a/bilby/gw/__init__.py b/bilby/gw/__init__.py index 75e3a3805a9c3829f450994d4742c17c7caaaa9d..9facdc6ee7c56c2d4fb7d5d2f61ef087cd623453 100644 --- a/bilby/gw/__init__.py +++ b/bilby/gw/__init__.py @@ -1,5 +1,8 @@ -from . import (calibration, conversion, detector, likelihood, prior, source, - utils, waveform_generator, result) +from . import (calibration, conversion, cosmology, detector, likelihood, prior, + result, source, utils, waveform_generator) + from .waveform_generator import WaveformGenerator from .likelihood import GravitationalWaveTransient + + diff --git a/bilby/gw/conversion.py b/bilby/gw/conversion.py index 65d9b12d9294689e8b1047aa1fbc0dd572d74ec2..ea5853d55eeac0440eae8cae539898fe9030ac10 100644 --- a/bilby/gw/conversion.py +++ b/bilby/gw/conversion.py @@ -6,42 +6,48 @@ from pandas import DataFrame from ..core.utils import logger, solar_mass from ..core.prior import DeltaFunction, Interped from .utils import lalsim_SimInspiralTransformPrecessingNewInitialConditions - +from .cosmology import get_cosmology try: - from astropy.cosmology import z_at_value, Planck15 - import astropy.units as u + from astropy import units + from astropy.cosmology import z_at_value except ImportError: - logger.warning("You do not have astropy installed currently. You will" - " not be able to use some of the prebuilt functions.") + logger.debug("You do not have astropy installed currently. You will" + " not be able to use some of the prebuilt functions.") -def redshift_to_luminosity_distance(redshift): - return Planck15.luminosity_distance(redshift).value +def redshift_to_luminosity_distance(redshift, cosmology=None): + cosmology = get_cosmology(cosmology) + return cosmology.luminosity_distance(redshift).value -def redshift_to_comoving_distance(redshift): - return Planck15.comoving_distance(redshift).value +def redshift_to_comoving_distance(redshift, cosmology=None): + cosmology = get_cosmology(cosmology) + return cosmology.comoving_distance(redshift).value @np.vectorize -def luminosity_distance_to_redshift(distance): - return z_at_value(Planck15.luminosity_distance, distance * u.Mpc) +def luminosity_distance_to_redshift(distance, cosmology=None): + cosmology = get_cosmology(cosmology) + return z_at_value(cosmology.luminosity_distance, distance * units.Mpc) @np.vectorize -def comoving_distance_to_redshift(distance): - return z_at_value(Planck15.comoving_distance, distance * u.Mpc) +def comoving_distance_to_redshift(distance, cosmology=None): + cosmology = get_cosmology(cosmology) + return z_at_value(cosmology.comoving_distance, distance * units.Mpc) -def comoving_distance_to_luminosity_distance(distance): - redshift = comoving_distance_to_redshift(distance) - return redshift_to_luminosity_distance(redshift) +def comoving_distance_to_luminosity_distance(distance, cosmology=None): + cosmology = get_cosmology(cosmology) + redshift = comoving_distance_to_redshift(distance, cosmology) + return redshift_to_luminosity_distance(redshift, cosmology) -def luminosity_distance_to_comoving_distance(distance): - redshift = luminosity_distance_to_redshift(distance) - return redshift_to_comoving_distance(redshift) +def luminosity_distance_to_comoving_distance(distance, cosmology=None): + cosmology = get_cosmology(cosmology) + redshift = luminosity_distance_to_redshift(distance, cosmology) + return redshift_to_comoving_distance(redshift, cosmology) @np.vectorize @@ -123,6 +129,23 @@ def convert_to_lal_binary_black_hole_parameters(parameters): converted_parameters = parameters.copy() original_keys = list(converted_parameters.keys()) + if 'redshift' in converted_parameters.keys(): + converted_parameters['luminosity_distance'] = \ + redshift_to_luminosity_distance(parameters['redshift']) + elif 'comoving_distance' in converted_parameters.keys(): + converted_parameters['luminosity_distance'] = \ + comoving_distance_to_luminosity_distance( + parameters['comoving_distance']) + + for key in original_keys: + if key[-7:] == '_source': + if 'redshift' not in converted_parameters.keys(): + converted_parameters['redshift'] =\ + luminosity_distance_to_redshift( + parameters['luminosity_distance']) + converted_parameters[key[:-7]] = converted_parameters[key] * ( + 1 + converted_parameters['redshift']) + if 'chirp_mass' in converted_parameters.keys(): if 'total_mass' in converted_parameters.keys(): converted_parameters['symmetric_mass_ratio'] =\ @@ -198,14 +221,6 @@ def convert_to_lal_binary_black_hole_parameters(parameters): converted_parameters[angle] =\ np.arccos(converted_parameters[cos_angle]) - if 'redshift' in converted_parameters.keys(): - converted_parameters['luminosity_distance'] =\ - redshift_to_luminosity_distance(parameters['redshift']) - elif 'comoving_distance' in converted_parameters.keys(): - converted_parameters['luminosity_distance'] = \ - comoving_distance_to_luminosity_distance( - parameters['comoving_distance']) - added_keys = [key for key in converted_parameters.keys() if key not in original_keys] @@ -243,6 +258,18 @@ def convert_to_lal_binary_neutron_star_parameters(parameters): converted_parameters, added_keys =\ convert_to_lal_binary_black_hole_parameters(converted_parameters) + for idx in ['1', '2']: + mag = 'a_{}'.format(idx) + if mag in original_keys: + tilt = 'tilt_{}'.format(idx) + if tilt in original_keys: + converted_parameters['chi_{}'.format(idx)] = ( + converted_parameters[mag] * + np.cos(converted_parameters[tilt])) + else: + converted_parameters['chi_{}'.format(idx)] = ( + converted_parameters[mag]) + # 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']]): @@ -273,18 +300,6 @@ def convert_to_lal_binary_neutron_star_parameters(parameters): * converted_parameters['mass_1']**5\ / converted_parameters['mass_2']**5 - for idx in ['1', '2']: - mag = 'a_{}'.format(idx) - if mag in original_keys: - tilt = 'tilt_{}'.format(idx) - if tilt in original_keys: - converted_parameters['chi_{}'.format(idx)] = ( - converted_parameters[mag] * - np.cos(converted_parameters[tilt])) - else: - converted_parameters['chi_{}'.format(idx)] = ( - converted_parameters[mag]) - added_keys = [key for key in converted_parameters.keys() if key not in original_keys] diff --git a/bilby/gw/cosmology.py b/bilby/gw/cosmology.py new file mode 100644 index 0000000000000000000000000000000000000000..df82acf7177fbfca7b5137c633997a078435339d --- /dev/null +++ b/bilby/gw/cosmology.py @@ -0,0 +1,62 @@ +from ..core.utils import logger + +try: + from astropy import cosmology as cosmo + DEFAULT_COSMOLOGY = cosmo.Planck15 + COSMOLOGY = [DEFAULT_COSMOLOGY, DEFAULT_COSMOLOGY.name] +except ImportError: + logger.debug("You do not have astropy installed currently. You will" + " not be able to use some of the prebuilt functions.") + DEFAULT_COSMOLOGY = None + COSMOLOGY = [None, str(None)] + + +def get_cosmology(cosmology=None): + if cosmology is None: + cosmology = COSMOLOGY[0] + elif isinstance(cosmology, str): + cosmology = cosmo.__dict__[cosmology] + return cosmology + + +def set_cosmology(cosmology=None): + """ + Get an instance of a astropy.cosmology.FLRW subclass. + + To avoid repeatedly instantiating the same class, test if it is the same + as the last used cosmology. + + Parameters + ---------- + cosmology: astropy.cosmology.FLRW, str, dict + Description of cosmology, one of: + None - Use DEFAULT_COSMOLOGY + Instance of astropy.cosmology.FLRW subclass + String with name of known Astropy cosmology, e.g., "Planck13" + Dictionary with arguments required to instantiate the cosmology + class. + + Returns + ------- + cosmo: astropy.cosmology.FLRW + Cosmology instance + """ + if cosmology is None: + cosmology = DEFAULT_COSMOLOGY + elif isinstance(cosmology, cosmo.FLRW): + cosmology = cosmology + elif isinstance(cosmology, str): + cosmology = cosmo.__dict__[cosmology] + elif isinstance(cosmology, dict): + if 'Ode0' in cosmology.keys(): + if 'w0' in cosmology.keys(): + cosmology = cosmo.wCDM(**cosmology) + else: + cosmology = cosmo.LambdaCDM(**cosmology) + else: + cosmology = cosmo.FlatLambdaCDM(**cosmology) + COSMOLOGY[0] = cosmology + if cosmology.name is not None: + COSMOLOGY[1] = cosmology.name + else: + COSMOLOGY[1] = repr(cosmology) diff --git a/bilby/gw/prior.py b/bilby/gw/prior.py index 03a4fb4db679ce677115ba9378e889b453e9fd5b..375124147a1dbbbb36e972dc0f0cda15adb969c1 100644 --- a/bilby/gw/prior.py +++ b/bilby/gw/prior.py @@ -3,31 +3,150 @@ import os import numpy as np from scipy.interpolate import UnivariateSpline -from ..core.prior import (PriorDict, Uniform, FromFile, Prior, DeltaFunction, - Gaussian, Interped) -from ..core.utils import logger - - -class UniformComovingVolume(FromFile): - - def __init__(self, minimum=None, maximum=None, - name='luminosity distance', latex_label='$d_L$', unit='Mpc'): - """ - - Parameters - ---------- - minimum: float, optional - See superclass - maximum: float, optional - See superclass - name: str, optional - See superclass - latex_label: str, optional - See superclass - """ - file_name = os.path.join(os.path.dirname(__file__), 'prior_files', 'comoving.txt') - FromFile.__init__(self, file_name=file_name, minimum=minimum, maximum=maximum, name=name, - latex_label=latex_label, unit=unit) +from ..core.prior import (PriorDict, Uniform, Prior, DeltaFunction, Gaussian, + Interped) +from ..core.utils import infer_args_from_method, logger +from .cosmology import get_cosmology + +try: + from astropy import cosmology as cosmo, units +except ImportError: + logger.debug("You do not have astropy installed currently. You will" + " not be able to use some of the prebuilt functions.") + + +class Cosmological(Interped): + + _default_args_dict = dict( + redshift=dict(name='redshift', latex_label='$z$', unit=None), + luminosity_distance=dict( + name='luminosity_distance', latex_label='$d_L$', unit=units.Mpc), + comoving_distance=dict( + name='comoving_distance', latex_label='$d_C$', unit=units.Mpc)) + + def __init__(self, minimum, maximum, cosmology=None, name=None, + latex_label=None, unit=None): + self.cosmology = get_cosmology(cosmology) + if name not in self._default_args_dict: + raise ValueError( + "Name {} not recognised. Must be one of luminosity_distance, " + "comoving_distance, redshift".format(name)) + self.name = name + label_args = self._default_args_dict[self.name] + if latex_label is not None: + label_args['latex_label'] = latex_label + if unit is not None: + if isinstance(unit, str): + unit = units.__dict__[unit] + label_args['unit'] = unit + self.unit = label_args['unit'] + self._minimum = dict() + self._maximum = dict() + self.minimum = minimum + self.maximum = maximum + if name == 'redshift': + xx, yy = self._get_redshift_arrays() + elif name == 'comoving_distance': + xx, yy = self._get_comoving_distance_arrays() + elif name == 'luminosity_distance': + xx, yy = self._get_luminosity_distance_arrays() + else: + raise ValueError('Name {} not recognized.'.format(name)) + Interped.__init__(self, xx=xx, yy=yy, minimum=minimum, maximum=maximum, + **label_args) + + @property + def minimum(self): + return self._minimum[self.name] + + @minimum.setter + def minimum(self, minimum): + cosmology = get_cosmology(self.cosmology) + self._minimum[self.name] = minimum + if self.name == 'redshift': + self._minimum['luminosity_distance'] =\ + cosmology.luminosity_distance(minimum).value + self._minimum['comoving_distance'] =\ + cosmology.comoving_distance(minimum).value + elif self.name == 'luminosity_distance': + if minimum == 0: + self._minimum['redshift'] = 0 + else: + self._minimum['redshift'] = cosmo.z_at_value( + cosmology.luminosity_distance, minimum * self.unit) + self._minimum['comoving_distance'] = self._minimum['redshift'] + elif self.name == 'comoving_distance': + if minimum == 0: + self._minimum['redshift'] = 0 + else: + self._minimum['redshift'] = cosmo.z_at_value( + cosmology.comoving_distance, minimum * self.unit) + self._minimum['luminosity_distance'] = self._minimum['redshift'] + if getattr(self._maximum, self.name, np.inf) < np.inf: + self.__update_instance() + + @property + def maximum(self): + return self._maximum[self.name] + + @maximum.setter + def maximum(self, maximum): + cosmology = get_cosmology(self.cosmology) + self._maximum[self.name] = maximum + if self.name == 'redshift': + self._maximum['luminosity_distance'] = \ + cosmology.luminosity_distance(maximum).value + self._maximum['comoving_distance'] = \ + cosmology.comoving_distance(maximum).value + elif self.name == 'luminosity_distance': + self._maximum['redshift'] = cosmo.z_at_value( + cosmology.luminosity_distance, maximum * self.unit) + self._maximum['comoving_distance'] = self._maximum['redshift'] + elif self.name == 'comoving_distance': + self._maximum['redshift'] = cosmo.z_at_value( + cosmology.comoving_distance, maximum * self.unit) + self._maximum['luminosity_distance'] = self._maximum['redshift'] + if getattr(self._minimum, self.name, np.inf) < np.inf: + self.__update_instance() + + def get_corresponding_prior(self, name=None, unit=None): + subclass_args = infer_args_from_method(self.__init__) + args_dict = {key: getattr(self, key) for key in subclass_args} + self._convert_to(new=name, args_dict=args_dict) + if unit is not None: + args_dict['unit'] = unit + return self.__class__(**args_dict) + + def _convert_to(self, new, args_dict): + args_dict.update(self._default_args_dict[new]) + args_dict['minimum'] = self._minimum[args_dict['name']] + args_dict['maximum'] = self._maximum[args_dict['name']] + + def _get_comoving_distance_arrays(self): + zs, p_dz = self._get_redshift_arrays() + dc_of_z = self.cosmology.comoving_distance(zs).value + ddc_dz = np.gradient(dc_of_z, zs) + p_dc = p_dz / ddc_dz + return dc_of_z, p_dc + + def _get_luminosity_distance_arrays(self): + zs, p_dz = self._get_redshift_arrays() + dl_of_z = self.cosmology.luminosity_distance(zs).value + ddl_dz = np.gradient(dl_of_z, zs) + p_dl = p_dz / ddl_dz + return dl_of_z, p_dl + + def _get_redshift_arrays(self): + raise NotImplementedError + + +class UniformComovingVolume(Cosmological): + + def _get_redshift_arrays(self): + zs = np.linspace(self._minimum['redshift'] * 0.99, + self._maximum['redshift'] * 1.01, 1000) + p_dz = self.cosmology.differential_comoving_volume(zs).value + return zs, p_dz class AlignedSpin(Interped): diff --git a/examples/injection_examples/fast_tutorial.py b/examples/injection_examples/fast_tutorial.py index cb9831e2a2e4f4b5b011b9d92dcad7cf15ead732..c0e5cee868c71522abab5fc655361c2b7e05020a 100644 --- a/examples/injection_examples/fast_tutorial.py +++ b/examples/injection_examples/fast_tutorial.py @@ -5,7 +5,7 @@ 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. +between luminosity distances of 100Mpc and 5Gpc, the cosmology is Planck15. """ from __future__ import division, print_function diff --git a/test/conversion_test.py b/test/conversion_test.py index 3d2d698582db814f0b58bf34e32b6f7d75a27132..7ee6bda1d60d0572061473d00fb7bb42ffa4af67 100644 --- a/test/conversion_test.py +++ b/test/conversion_test.py @@ -1,9 +1,12 @@ from __future__ import division, absolute_import import unittest import mock -from bilby.gw import conversion + import numpy as np +import bilby +from bilby.gw import conversion + class TestBasicConversions(unittest.TestCase): @@ -105,7 +108,7 @@ class TestBasicConversions(unittest.TestCase): def test_lambda_1_lambda_2_to_delta_lambda(self): delta_lambda = \ - conversion.lambda_1_lambda_2_to_lambda_tilde( + conversion.lambda_1_lambda_2_to_delta_lambda( self.lambda_1, self.lambda_2, self.mass_1, self.mass_2) self.assertTrue((self.delta_lambda - delta_lambda) < 1e-5) @@ -115,31 +118,165 @@ class TestConvertToLALParams(unittest.TestCase): def setUp(self): self.search_keys = [] self.parameters = dict() - self.remove = True + self.component_mass_pars = dict(mass_1=1.4, mass_2=1.4) + self.mass_parameters = self.component_mass_pars.copy() + self.mass_parameters['mass_ratio'] =\ + conversion.component_masses_to_mass_ratio(**self.component_mass_pars) + self.mass_parameters['symmetric_mass_ratio'] = \ + conversion.component_masses_to_symmetric_mass_ratio(**self.component_mass_pars) + self.mass_parameters['chirp_mass'] = \ + conversion.component_masses_to_chirp_mass(**self.component_mass_pars) + self.mass_parameters['total_mass'] = \ + conversion.component_masses_to_total_mass(**self.component_mass_pars) + self.component_tidal_parameters = dict(lambda_1=300, lambda_2=300) + self.all_component_pars = self.component_tidal_parameters.copy() + self.all_component_pars.update(self.component_mass_pars) + self.tidal_parameters = self.component_tidal_parameters.copy() + self.tidal_parameters['lambda_tilde'] = \ + conversion.lambda_1_lambda_2_to_lambda_tilde(**self.all_component_pars) + self.tidal_parameters['delta_lambda'] = \ + conversion.lambda_1_lambda_2_to_delta_lambda(**self.all_component_pars) def tearDown(self): del self.search_keys del self.parameters - del self.remove + + def bbh_convert(self): + self.parameters, self.added_keys =\ + conversion.convert_to_lal_binary_black_hole_parameters( + self.parameters) + + def bns_convert(self): + self.parameters, self.added_keys = \ + conversion.convert_to_lal_binary_neutron_star_parameters( + self.parameters) + + def test_redshift_to_luminosity_distance(self): + self.parameters['redshift'] = 1 + dl = conversion.redshift_to_luminosity_distance( + self.parameters['redshift']) + self.bbh_convert() + self.assertEqual(self.parameters['luminosity_distance'], dl) + + def test_comoving_to_luminosity_distance(self): + self.parameters['comoving_distance'] = 1 + dl = conversion.comoving_distance_to_luminosity_distance( + self.parameters['comoving_distance']) + self.bbh_convert() + self.assertEqual(self.parameters['luminosity_distance'], dl) + + def test_source_to_lab_frame(self): + self.parameters['test_source'] = 1 + self.parameters['redshift'] = 1 + lab = self.parameters['test_source'] * (1 + self.parameters['redshift']) + self.bbh_convert() + self.assertEqual(self.parameters['test'], lab) + + def _conversion_to_component_mass(self, keys): + for key in keys: + self.parameters[key] = self.mass_parameters[key] + self.bbh_convert() + self.assertAlmostEqual( + max([abs(self.parameters[key] - self.component_mass_pars[key]) + for key in ['mass_1', 'mass_2']]), 0) + + def test_chirp_mass_total_mass(self): + self._conversion_to_component_mass(['chirp_mass', 'total_mass']) + + def test_chirp_mass_sym_mass_ratio(self): + self._conversion_to_component_mass(['chirp_mass', 'symmetric_mass_ratio']) + + def test_chirp_mass_mass_ratio(self): + self._conversion_to_component_mass(['chirp_mass', 'mass_ratio']) + + def test_total_mass_sym_mass_ratio(self): + self._conversion_to_component_mass(['total_mass', 'symmetric_mass_ratio']) + + def test_total_mass_mass_ratio(self): + self._conversion_to_component_mass(['total_mass', 'mass_ratio']) + + def test_total_mass_mass_1(self): + self._conversion_to_component_mass(['total_mass', 'mass_1']) + + def test_total_mass_mass_2(self): + self._conversion_to_component_mass(['total_mass', 'mass_2']) + + def test_sym_mass_ratio_mass_1(self): + self._conversion_to_component_mass(['symmetric_mass_ratio', 'mass_1']) + + def test_sym_mass_ratio_mass_2(self): + self._conversion_to_component_mass(['symmetric_mass_ratio', 'mass_2']) + + def test_mass_ratio_mass_1(self): + self._conversion_to_component_mass(['mass_ratio', 'mass_1']) + + def test_mass_ratio_mass_2(self): + self._conversion_to_component_mass(['mass_ratio', 'mass_2']) + + def test_bbh_aligned_spin_to_spherical(self): + self.parameters['chi_1'] = -0.5 + a_1 = abs(self.parameters['chi_1']) + tilt_1 = np.arccos(np.sign(self.parameters['chi_1'])) + phi_jl = 0.0 + phi_12 = 0.0 + self.bbh_convert() + self.assertDictEqual( + {key: self.parameters[key] + for key in ['a_1', 'tilt_1', 'phi_12', 'phi_jl']}, + dict(a_1=a_1, tilt_1=tilt_1, phi_jl=phi_jl, phi_12=phi_12)) 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 - self.parameters, _ =\ - conversion.convert_to_lal_binary_black_hole_parameters( - 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)) + self.parameters['cos_tilt_1'] = 1 + t1 = np.arccos(self.parameters['cos_tilt_1']) + self.bbh_convert() + self.assertEqual(self.parameters['tilt_1'], t1) + + def _conversion_to_component_tidal(self, keys): + for key in keys: + self.parameters[key] = self.tidal_parameters[key] + for key in ['mass_1', 'mass_2']: + self.parameters[key] = self.mass_parameters[key] + self.bns_convert() + component_dict = { + key: self.parameters[key] for key in ['lambda_1', 'lambda_2']} + self.assertDictEqual(component_dict, self.component_tidal_parameters) + + def test_lambda_tilde_delta_lambda(self): + self._conversion_to_component_tidal(['lambda_tilde', 'delta_lambda']) + + def test_lambda_tilde(self): + self._conversion_to_component_tidal(['lambda_tilde']) + + def test_lambda_1(self): + self._conversion_to_component_tidal(['lambda_1']) + + def test_bns_spherical_spin_to_aligned(self): + self.parameters['a_1'] = -1 + self.parameters['tilt_1'] = np.arccos(0.5) + chi_1 = self.parameters['a_1'] * np.cos(self.parameters['tilt_1']) + self.bns_convert() + self.assertEqual(self.parameters['chi_1'], chi_1) + + +class TestDistanceTransformations(unittest.TestCase): + + def setUp(self): + self.distances = np.linspace(1, 1000, 100) + + def test_luminosity_redshift_with_cosmology(self): + z = conversion.luminosity_distance_to_redshift(self.distances, cosmology='WMAP9') + dl = conversion.redshift_to_luminosity_distance(z, cosmology='WMAP9') + self.assertAlmostEqual(max(abs(dl - self.distances)), 0, 4) + + def test_comoving_redshift_with_cosmology(self): + z = conversion.comoving_distance_to_redshift(self.distances, cosmology='WMAP9') + dc = conversion.redshift_to_comoving_distance(z, cosmology='WMAP9') + self.assertAlmostEqual(max(abs(dc - self.distances)), 0, 4) + + def test_comoving_luminosity_with_cosmology(self): + dc = conversion.comoving_distance_to_luminosity_distance(self.distances, cosmology='WMAP9') + dl = conversion.luminosity_distance_to_comoving_distance(dc, cosmology='WMAP9') + self.assertAlmostEqual(max(abs(dl - self.distances)), 0, 4) if __name__ == '__main__': diff --git a/test/cosmology_test.py b/test/cosmology_test.py new file mode 100644 index 0000000000000000000000000000000000000000..ca8ef65693497ea6f1c71affabea391e9084202d --- /dev/null +++ b/test/cosmology_test.py @@ -0,0 +1,56 @@ +from __future__ import division, absolute_import +import unittest + +from astropy.cosmology import WMAP9, Planck15 +from bilby.gw import cosmology + + +class TestSetCosmology(unittest.TestCase): + + def setUp(self): + pass + + def test_setting_cosmology_with_string(self): + cosmology.set_cosmology('WMAP9') + self.assertEqual(cosmology.COSMOLOGY[1], 'WMAP9') + cosmology.set_cosmology('Planck15') + + def test_setting_cosmology_with_astropy_object(self): + cosmology.set_cosmology(WMAP9) + self.assertEqual(cosmology.COSMOLOGY[1], 'WMAP9') + cosmology.set_cosmology(Planck15) + + def test_setting_cosmology_with_default(self): + cosmology.set_cosmology() + self.assertEqual(cosmology.COSMOLOGY[1], cosmology.DEFAULT_COSMOLOGY.name) + + def test_setting_cosmology_with_flat_lambda_cdm_dict(self): + cosmo_dict = dict(H0=67.7, Om0=0.3) + cosmology.set_cosmology(cosmo_dict) + self.assertEqual(cosmology.COSMOLOGY[1][:13], 'FlatLambdaCDM') + + def test_setting_cosmology_with_lambda_cdm_dict(self): + cosmo_dict = dict(H0=67.7, Om0=0.3, Ode0=0.7) + cosmology.set_cosmology(cosmo_dict) + self.assertEqual(cosmology.COSMOLOGY[1][:9], 'LambdaCDM') + + def test_setting_cosmology_with_w_cdm_dict(self): + cosmo_dict = dict(H0=67.7, Om0=0.3, Ode0=0.7, w0=-1.0) + cosmology.set_cosmology(cosmo_dict) + self.assertEqual(cosmology.COSMOLOGY[1][:4], 'wCDM') + + +class TestGetCosmology(unittest.TestCase): + + def setUp(self): + pass + + def test_getting_cosmology_with_string(self): + self.assertEqual(cosmology.get_cosmology('WMAP9').name, 'WMAP9') + + def test_getting_cosmology_with_default(self): + self.assertEqual(cosmology.get_cosmology(), cosmology.COSMOLOGY[0]) + + +if __name__ == '__main__': + unittest.main() diff --git a/test/gw_prior_test.py b/test/gw_prior_test.py index dcab53d9e0bf1fd829e32c312f55563d9968bbf0..f7c7454e9472f496b048045e1719b23f08321cfa 100644 --- a/test/gw_prior_test.py +++ b/test/gw_prior_test.py @@ -1,9 +1,13 @@ from __future__ import division, absolute_import import unittest -import bilby import os import sys +import numpy as np +from astropy import cosmology + +import bilby + class TestBBHPriorDict(unittest.TestCase): @@ -69,5 +73,73 @@ class TestCalibrationPrior(unittest.TestCase): self.assertEqual(len(test), n_nodes * 3) +class TestUniformComovingVolumePrior(unittest.TestCase): + + def setUp(self): + pass + + def test_minimum(self): + prior = bilby.gw.prior.UniformComovingVolume( + minimum=10, maximum=10000, name='luminosity_distance') + self.assertEqual(prior.minimum, 10) + + def test_maximum(self): + prior = bilby.gw.prior.UniformComovingVolume( + minimum=10, maximum=10000, name='luminosity_distance') + self.assertEqual(prior.maximum, 10000) + + def test_zero_minimum_works(self): + prior = bilby.gw.prior.UniformComovingVolume( + minimum=0, maximum=10000, name='luminosity_distance') + self.assertEqual(prior.minimum, 0) + + def test_specify_cosmology(self): + prior = bilby.gw.prior.UniformComovingVolume( + minimum=10, maximum=10000, name='luminosity_distance', + cosmology='Planck13') + self.assertEqual(repr(prior.cosmology), repr(cosmology.Planck13)) + + def test_comoving_prior_creation(self): + prior = bilby.gw.prior.UniformComovingVolume( + minimum=10, maximum=1000, name='comoving_distance') + self.assertEqual(prior.latex_label, '$d_C$') + + def test_redshift_prior_creation(self): + prior = bilby.gw.prior.UniformComovingVolume( + minimum=0.1, maximum=1, name='redshift') + self.assertEqual(prior.latex_label, '$z$') + + def test_redshift_to_luminosity_distance(self): + prior = bilby.gw.prior.UniformComovingVolume( + minimum=0.1, maximum=1, name='redshift') + new_prior = prior.get_corresponding_prior('luminosity_distance') + self.assertEqual(new_prior.name, 'luminosity_distance') + + def test_luminosity_distance_to_redshift(self): + prior = bilby.gw.prior.UniformComovingVolume( + minimum=10, maximum=10000, name='luminosity_distance') + new_prior = prior.get_corresponding_prior('redshift') + self.assertEqual(new_prior.name, 'redshift') + + def test_luminosity_distance_to_comoving_distance(self): + prior = bilby.gw.prior.UniformComovingVolume( + minimum=10, maximum=10000, name='luminosity_distance') + new_prior = prior.get_corresponding_prior('comoving_distance') + self.assertEqual(new_prior.name, 'comoving_distance') + + +class TestAlignedSpin(unittest.TestCase): + + def setUp(self): + pass + + def test_default_prior_matches_analytic(self): + prior = bilby.gw.prior.AlignedSpin() + chis = np.linspace(-1, 1, 20) + analytic = - np.log(np.abs(chis)) / 2 + max_difference = max(abs(analytic - prior.prob(chis))) + self.assertAlmostEqual(max_difference, 0, 2) + + if __name__ == '__main__': unittest.main() diff --git a/test/prior_test.py b/test/prior_test.py index 5d78d919bcff173cb077387006c4fd644c76d309..d8d4041de3a719e4a4238945f60e522b3a641388 100644 --- a/test/prior_test.py +++ b/test/prior_test.py @@ -133,7 +133,7 @@ class TestPriorClasses(unittest.TestCase): bilby.core.prior.PowerLaw(name='test', unit='unit', alpha=2, minimum=1, maximum=1e2), bilby.core.prior.Uniform(name='test', unit='unit', minimum=0, maximum=1), bilby.core.prior.LogUniform(name='test', unit='unit', minimum=5e0, maximum=1e2), - bilby.gw.prior.UniformComovingVolume(name='test', minimum=2e2, maximum=5e3), + bilby.gw.prior.UniformComovingVolume(name='redshift', minimum=0.1, maximum=1.0), bilby.core.prior.Sine(name='test', unit='unit'), bilby.core.prior.Cosine(name='test', unit='unit'), bilby.core.prior.Interped(name='test', unit='unit', xx=np.linspace(0, 10, 1000), @@ -153,7 +153,7 @@ class TestPriorClasses(unittest.TestCase): bilby.core.prior.Lorentzian(name='test', unit='unit', alpha=0, beta=1), bilby.core.prior.Gamma(name='test', unit='unit', k=1, theta=1), bilby.core.prior.ChiSquared(name='test', unit='unit', nu=2), - bilby.gw.prior.AlignedSpin(name='test', unit='unit') + bilby.gw.prior.AlignedSpin(name='test', unit='unit'), ] def test_minimum_rescaling(self): @@ -253,8 +253,8 @@ class TestPriorClasses(unittest.TestCase): def test_unit_setting(self): for prior in self.priors: - if isinstance(prior, bilby.gw.prior.UniformComovingVolume): - self.assertEqual('Mpc', prior.unit) + if isinstance(prior, bilby.gw.prior.Cosmological): + self.assertEqual(None, prior.unit) else: self.assertEqual('unit', prior.unit)