diff --git a/bilby/core/sampler/emcee.py b/bilby/core/sampler/emcee.py index b52c2bee052e10e2e144bef3f5d0aabcfff5e1ea..3117a6c7a2cdd11266bce26fda6e9b06feb9a36d 100644 --- a/bilby/core/sampler/emcee.py +++ b/bilby/core/sampler/emcee.py @@ -51,6 +51,11 @@ class Emcee(MCMCSampler): use_ratio=False, plot=False, skip_import_verification=False, pos0=None, nburn=None, burn_in_fraction=0.25, resume=True, burn_in_act=3, **kwargs): + import emcee + if LooseVersion(emcee.__version__) > LooseVersion('2.2.1'): + self.prerelease = True + else: + self.prerelease = False MCMCSampler.__init__( self, likelihood=likelihood, priors=priors, outdir=outdir, label=label, use_ratio=use_ratio, plot=plot, @@ -82,7 +87,6 @@ class Emcee(MCMCSampler): @property def sampler_function_kwargs(self): import emcee - keys = ['lnprob0', 'rstate0', 'blobs0', 'iterations', 'thin', 'storechain', 'mh_proposal'] # updated function keywords for emcee > v2.2.1 @@ -93,7 +97,7 @@ class Emcee(MCMCSampler): function_kwargs = {key: self.kwargs[key] for key in keys if key in self.kwargs} function_kwargs['p0'] = self.pos0 - if LooseVersion(emcee.__version__) > LooseVersion('2.2.1'): + if self.prerelease: if function_kwargs['mh_proposal'] is not None: logger.warning("The 'mh_proposal' option is no longer used " "in emcee v{}, and will be ignored.".format(emcee.__version__)) @@ -109,8 +113,6 @@ class Emcee(MCMCSampler): @property def sampler_init_kwargs(self): - import emcee - init_kwargs = {key: value for key, value in self.kwargs.items() if key not in self.sampler_function_kwargs} @@ -122,7 +124,7 @@ class Emcee(MCMCSampler): updatekeys = {'dim': 'ndim', 'lnpostfn': 'log_prob_fn'} - if LooseVersion(emcee.__version__) > LooseVersion('2.2.1'): + if self.prerelease: for key in updatekeys: if key in init_kwargs: init_kwargs[updatekeys[key]] = init_kwargs.pop(key) @@ -193,8 +195,10 @@ class Emcee(MCMCSampler): for sample in tqdm(sampler.sample(**self.sampler_function_kwargs), total=self.nsteps): - points = np.hstack([sample[0], np.array(sample[3])]) - # import IPython; IPython.embed() + if self.prerelease: + points = np.hstack([sample.coords, sample.blobs]) + else: + points = np.hstack([sample[0], np.array(sample[3])]) with open(out_file, "a") as ff: for ii, point in enumerate(points): ff.write(template.format(ii, *point)) diff --git a/bilby/gw/detector.py b/bilby/gw/detector.py index 6da387ab93d592bd360b06dff4028b7ebc4d5c3d..9b01722dc445c8c0c1cea58b12ead3a864045cde 100644 --- a/bilby/gw/detector.py +++ b/bilby/gw/detector.py @@ -147,38 +147,8 @@ class InterferometerList(list): if utils.command_line_args.test: return - fig = plt.figure() - for ii, interferometer in enumerate(self): - ax = fig.add_subplot(len(self) // 2, 2, ii + 1) - ax.loglog(interferometer.frequency_array, - gwutils.asd_from_freq_series(freq_data=interferometer.frequency_domain_strain, - df=(interferometer.frequency_array[1] - - interferometer.frequency_array[0])), - color='C0', label=interferometer.name) - ax.loglog(interferometer.frequency_array, - interferometer.amplitude_spectral_density_array, - color='C1', lw=0.5, label=interferometer.name + ' ASD') - ax.grid('on') - ax.set_ylabel(r'strain [strain/$\sqrt{\rm Hz}$]') - ax.set_xlabel(r'frequency [Hz]') - ax.set_xlim(20, 2000) - ax.legend(loc='best') - if signal is not None: - ax.loglog(self.frequency_array, - gwutils.asd_from_freq_series(freq_data=signal, - df=(self.frequency_array[1] - - self.frequency_array[0]) - ), - color='C2', - label='Signal') - fig.tight_layout() - if label is None: - fig.savefig( - '{}/frequency_domain_data.png'.format(outdir)) - else: - fig.savefig( - '{}/{}_frequency_domain_data.png'.format( - outdir, label)) + for interferometer in self: + interferometer.plot_data(signal=signal, outdir=outdir, label=label) @property def number_of_interferometers(self): diff --git a/bilby/gw/source.py b/bilby/gw/source.py index 6d4071776db7f9f9f361b7eb04b8ad354ab728cf..22cc8f2353401e6258b7a79324c281331d776a00 100644 --- a/bilby/gw/source.py +++ b/bilby/gw/source.py @@ -4,12 +4,14 @@ import numpy as np from ..core import utils from ..core.utils import logger, spherical_to_cartesian +from .conversion import transform_precessing_spins from .utils import (lalsim_GetApproximantFromString, lalsim_SimInspiralFD, lalsim_SimInspiralChooseFDWaveform, lalsim_SimInspiralWaveformParamsInsertTidalLambda1, - lalsim_SimInspiralWaveformParamsInsertTidalLambda2) -from .conversion import transform_precessing_spins + lalsim_SimInspiralWaveformParamsInsertTidalLambda2, + lalsim_SimIMRPhenomPCalculateModelParametersFromSourceFrame, + lalsim_SimIMRPhenomPFrequencySequence) try: import lal @@ -402,20 +404,20 @@ def roq(frequency_array, mass_1, mass_2, luminosity_distance, a_1, tilt_1, reference_frequency, phase) chi_1_l, chi_2_l, chi_p, theta_jn, alpha, phase_aligned, zeta =\ - lalsim.SimIMRPhenomPCalculateModelParametersFromSourceFrame( + lalsim_SimIMRPhenomPCalculateModelParametersFromSourceFrame( mass_1, mass_2, reference_frequency, phase, iota, spin_1x, spin_1y, spin_1z, spin_2x, spin_2y, spin_2z, version) waveform_polarizations = dict() - h_linear_plus, h_linear_cross = lalsim.SimIMRPhenomPFrequencySequence( + h_linear_plus, h_linear_cross = lalsim_SimIMRPhenomPFrequencySequence( frequency_nodes_linear, chi_1_l, chi_2_l, chi_p, theta_jn, mass_1, mass_2, luminosity_distance, - alpha, phase_aligned, reference_frequency, version, None) - h_quadratic_plus, h_quadratic_cross = lalsim.SimIMRPhenomPFrequencySequence( + alpha, phase_aligned, reference_frequency, version) + h_quadratic_plus, h_quadratic_cross = lalsim_SimIMRPhenomPFrequencySequence( frequency_nodes_quadratic, chi_1_l, chi_2_l, chi_p, theta_jn, mass_1, mass_2, luminosity_distance, - alpha, phase_aligned, reference_frequency, version, None) + alpha, phase_aligned, reference_frequency, version) waveform_polarizations['linear'] = dict( plus=(np.cos(2 * zeta) * h_linear_plus.data.data + diff --git a/bilby/gw/utils.py b/bilby/gw/utils.py index e4c2007342b9f7387d2d28d821bbecfe74558639..94c3cafb2c0cdb63a73e72f6cef8db7c9757974a 100644 --- a/bilby/gw/utils.py +++ b/bilby/gw/utils.py @@ -799,6 +799,33 @@ def lalsim_SimInspiralChooseFDWaveform( waveform_dictionary, approximant) +def lalsim_SimIMRPhenomPCalculateModelParametersFromSourceFrame( + mass_1, mass_2, reference_frequency, phase, iota, spin_1x, + spin_1y, spin_1z, spin_2x, spin_2y, spin_2z, version): + [mass_1, mass_2, reference_frequency, phase, iota, spin_1x, + spin_1y, spin_1z, spin_2x, spin_2y, spin_2z] = convert_args_list_to_float( + mass_1, mass_2, reference_frequency, phase, iota, spin_1x, + spin_1y, spin_1z, spin_2x, spin_2y, spin_2z) + return lalsim.SimIMRPhenomPCalculateModelParametersFromSourceFrame( + mass_1, mass_2, reference_frequency, phase, iota, spin_1x, + spin_1y, spin_1z, spin_2x, spin_2y, spin_2z, version) + + +def lalsim_SimIMRPhenomPFrequencySequence( + frequency_nodes, chi_1_l, chi_2_l, chi_p, theta_jn, + mass_1, mass_2, luminosity_distance, + alpha, phase_aligned, reference_frequency, version): + [chi_1_l, chi_2_l, chi_p, theta_jn, mass_1, mass_2, luminosity_distance, + alpha, phase_aligned, reference_frequency] = convert_args_list_to_float( + chi_1_l, chi_2_l, chi_p, theta_jn, mass_1, mass_2, luminosity_distance, + alpha, phase_aligned, reference_frequency) + + return lalsim.SimIMRPhenomPFrequencySequence( + frequency_nodes, chi_1_l, chi_2_l, chi_p, theta_jn, mass_1, mass_2, + luminosity_distance, alpha, phase_aligned, reference_frequency, version, + None) + + def lalsim_SimInspiralWaveformParamsInsertTidalLambda1( waveform_dictionary, lambda_1): try: diff --git a/test/detector_test.py b/test/detector_test.py index 95b2f5f3247d2c604b5200c8ac0f4c82b5ebf3ef..edaa139fc74c12f694254986771decf16512ceee 100644 --- a/test/detector_test.py +++ b/test/detector_test.py @@ -39,7 +39,8 @@ class TestInterferometer(unittest.TestCase): xarm_tilt=self.xarm_tilt, yarm_tilt=self.yarm_tilt) self.ifo.strain_data.set_from_frequency_domain_strain( np.linspace(0, 4096, 4097), sampling_frequency=4096, duration=2) - bilby.core.utils.check_directory_exists_and_if_not_mkdir('outdir') + self.outdir = 'outdir' + bilby.core.utils.check_directory_exists_and_if_not_mkdir(self.outdir) def tearDown(self): del self.name @@ -55,7 +56,7 @@ class TestInterferometer(unittest.TestCase): del self.xarm_tilt del self.yarm_tilt del self.ifo - rmtree('outdir') + rmtree(self.outdir) def test_name_setting(self): self.assertEqual(self.ifo.name, self.name) @@ -798,7 +799,8 @@ class TestInterferometerList(unittest.TestCase): self.ifo2.strain_data.set_from_frequency_domain_strain( self.frequency_arrays, sampling_frequency=4096, duration=2) self.ifo_list = bilby.gw.detector.InterferometerList([self.ifo1, self.ifo2]) - bilby.core.utils.check_directory_exists_and_if_not_mkdir('outdir') + self.outdir = 'outdir' + bilby.core.utils.check_directory_exists_and_if_not_mkdir(self.outdir) def tearDown(self): del self.frequency_arrays @@ -829,7 +831,7 @@ class TestInterferometerList(unittest.TestCase): del self.ifo1 del self.ifo2 del self.ifo_list - rmtree('outdir') + rmtree(self.outdir) def test_init_with_string(self): with self.assertRaises(TypeError): @@ -986,6 +988,15 @@ class TestInterferometerList(unittest.TestCase): with self.assertRaises(TypeError): bilby.gw.detector.InterferometerList.from_hdf5(filename) + def test_plot_data(self): + ifos = bilby.gw.detector.InterferometerList(['H1', 'L1']) + ifos.set_strain_data_from_power_spectral_densities(2048, 4) + ifos.plot_data(outdir=self.outdir) + + ifos = bilby.gw.detector.InterferometerList(['H1', 'L1', 'V1']) + ifos.set_strain_data_from_power_spectral_densities(2048, 4) + ifos.plot_data(outdir=self.outdir) + class TestPowerSpectralDensityWithoutFiles(unittest.TestCase): diff --git a/test/likelihood_test.py b/test/likelihood_test.py index 595c1f9affea1e0ef62ac1373765ba767d9bd2a0..cf878ad85d4407bcceabd6a0d641cace791a34b6 100644 --- a/test/likelihood_test.py +++ b/test/likelihood_test.py @@ -31,6 +31,18 @@ class TestLikelihoodBase(unittest.TestCase): def test_base_log_likelihood_ratio(self): self.assertTrue(np.isnan(self.likelihood.log_likelihood_ratio())) + def test_meta_data_unset(self): + self.assertEqual(self.likelihood.meta_data, None) + + def test_meta_data_set_fail(self): + with self.assertRaises(ValueError): + self.likelihood.meta_data = 10 + + def test_meta_data(self): + meta_data = dict(x=1, y=2) + self.likelihood.meta_data = meta_data + self.assertEqual(self.likelihood.meta_data, meta_data) + class TestAnalytical1DLikelihood(unittest.TestCase): diff --git a/test/prior_test.py b/test/prior_test.py index d8d4041de3a719e4a4238945f60e522b3a641388..4b85a66c625b43f1194fa79213278e807b65f399 100644 --- a/test/prior_test.py +++ b/test/prior_test.py @@ -4,7 +4,7 @@ import unittest from mock import Mock import numpy as np import os -import shutil +import copy from collections import OrderedDict @@ -130,6 +130,7 @@ class TestPriorClasses(unittest.TestCase): bilby.core.prior.Gaussian(name='test', unit='unit', mu=0, sigma=1), bilby.core.prior.Normal(name='test', unit='unit', mu=0, sigma=1), bilby.core.prior.PowerLaw(name='test', unit='unit', alpha=0, minimum=0, maximum=1), + bilby.core.prior.PowerLaw(name='test', unit='unit', alpha=-1, minimum=0.5, maximum=1), 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), @@ -205,6 +206,28 @@ class TestPriorClasses(unittest.TestCase): outside_domain = np.linspace(prior.minimum - 1e4, prior.minimum - 1, 1000) self.assertTrue(all(prior.prob(outside_domain) == 0)) + def test_prob_and_ln_prob(self): + for prior in self.priors: + sample = prior.sample() + self.assertAlmostEqual(np.log(prior.prob(sample)), prior.ln_prob(sample), 12) + + def test_log_normal_fail(self): + with self.assertRaises(ValueError): + bilby.core.prior.LogNormal(name='test', unit='unit', mu=0, sigma=-1) + + def test_studentt_fail(self): + with self.assertRaises(ValueError): + bilby.core.prior.StudentT(name='test', unit='unit', df=3, mu=0, scale=-1) + with self.assertRaises(ValueError): + bilby.core.prior.StudentT(name='test', unit='unit', df=0, mu=0, scale=1) + + def test_beta_fail(self): + with self.assertRaises(ValueError): + bilby.core.prior.Beta(name='test', unit='unit', alpha=-2.0, beta=2.0), + + with self.assertRaises(ValueError): + bilby.core.prior.Beta(name='test', unit='unit', alpha=2.0, beta=-2.0), + def test_probability_in_domain(self): """Test that the prior probability is non-negative in domain of validity and zero outside.""" for prior in self.priors: @@ -319,6 +342,15 @@ class TestPriorDict(unittest.TestCase): del self.default_prior_file del self.prior_set_from_file + def test_copy(self): + priors = bilby.core.prior.PriorDict(self.priors) + self.assertEqual(priors, priors.copy()) + + def test_prior_set(self): + priors_dict = bilby.core.prior.PriorDict(self.priors) + priors_set = bilby.core.prior.PriorSet(self.priors) + self.assertEqual(priors_dict, priors_set) + def test_prior_set_is_ordered_dict(self): self.assertIsInstance(self.prior_set_from_dict, OrderedDict)