Merge branch 'master' into 'add-ptemcee-sampler'

# Conflicts:
#   peyote/utils.py

.gitlab-ci.yml

@@ -24,6 +24,7 @@ exitcode-jessie:
     - coverage run test/prior_tests.py
     - coverage run test/tests.py
     - coverage run test/waveform_generator_tests.py
+    - coverage run test/noise_realisation_tests.py
     - coverage html --include=peyote/*
     - coverage-badge -o coverage.svg
   artifacts:

peyote/detector.py

@@ -495,12 +495,13 @@ def get_inteferometer(
 
     utils.check_directory_exists_and_if_not_mkdir(outdir)
 
-    strain = TimeSeries.fetch_open_data(
-            name, center_time-T/2, center_time+T/2, cache=cache, **kwargs)
+    strain = get_open_strain_data(
+            name, center_time-T/2, center_time+T/2, outdir=outdir, cache=cache,
+            **kwargs)
 
-    strain_psd = TimeSeries.fetch_open_data(
+    strain_psd = get_open_strain_data(
             name, center_time+psd_offset, center_time+psd_offset+psd_duration,
-            cache=cache, **kwargs)
+            outdir=outdir, cache=cache, **kwargs)
 
     sampling_frequency = int(strain.sample_rate.value)
 
@@ -552,3 +553,17 @@ def get_inteferometer(
         fig.savefig('{}/{}_frequency_domain_data.png'.format(outdir, name))
 
     return interferometer, sampling_frequency, time_duration
+
+
+def get_open_strain_data(name, t1, t2, outdir, cache=False, **kwargs):
+    filename = '{}/{}_{}_{}.txt'.format(outdir, name, t1, t2)
+    if os.path.isfile(filename) and cache:
+        logging.info('Using cached data from {}'.format(filename))
+        strain = TimeSeries.read(filename)
+    else:
+        logging.info('Fetching open data ...')
+        strain = TimeSeries.fetch_open_data(name, t1, t2, **kwargs)
+        logging.info('Saving data to {}'.format(filename))
+        strain.write(filename)
+    return strain
+

peyote/prior.py

 #!/bin/python
+from __future__ import division
 
 import numpy as np
 from scipy.interpolate import interp1d
@@ -57,15 +58,15 @@ class Prior(object):
         elif self.name == 'mchirp':
             return '$\mathcal{M}$'
         elif self.name == 'q':
-            return 'q'
+            return '$q$'
         elif self.name == 'a1':
-            return 'a_1'
+            return '$a_1$'
         elif self.name == 'a2':
-            return 'a_2'
+            return '$a_2$'
         elif self.name == 'tilt1':
-            return 'tilt_1'
+            return '$\\text{tilt}_1$'
         elif self.name == 'tilt2':
-            return 'tilt_2'
+            return '$\\text{tilt}_2$'
         elif self.name == 'phi1':
             return '$\phi_1$'
         elif self.name == 'phi2':
@@ -262,17 +263,17 @@ def create_default_prior(name):
         prior = PowerLaw(name=name, alpha=0, bounds=(5, 100))
     elif name == 'q':
         prior = PowerLaw(name=name, alpha=0, bounds=(0, 1))
-    elif name == 'a1':
-        prior = PowerLaw(name=name, alpha=0, bounds=(0, 1))
-    elif name == 'a2':
-        prior = PowerLaw(name=name, alpha=0, bounds=(0, 1))
-    elif name == 'tilt1':
+    elif name == 'a_1':
+        prior = PowerLaw(name=name, alpha=0, bounds=(0, 0.8))
+    elif name == 'a_2':
+        prior = PowerLaw(name=name, alpha=0, bounds=(0, 0.8))
+    elif name == 'tilt_1':
         prior = Sine(name=name)
-    elif name == 'tilt2':
+    elif name == 'tilt_2':
         prior = Sine(name=name)
-    elif name == 'phi1':
+    elif name == 'phi_1':
         prior = PowerLaw(name=name, alpha=0, bounds=(0, 2 * np.pi))
-    elif name == 'phi2':
+    elif name == 'phi_2':
         prior = PowerLaw(name=name, alpha=0, bounds=(0, 2 * np.pi))
     elif name == 'luminosity_distance':
         prior = PowerLaw(name=name, alpha=2, bounds=(1e2, 5e3))

peyote/sampler.py

@@ -100,6 +100,15 @@ class Sampler(object):
     def kwargs(self, kwargs):
         self.__kwargs = kwargs
 
+    def verify_kwargs_against_external_sampler_function(self):
+        args = inspect.getargspec(self.external_sampler_function).args
+        for user_input in self.kwargs.keys():
+            if user_input not in args:
+                logging.warning(
+                    "Supplied argument '{}' not an argument of '{}', removing."
+                    .format(user_input, self.external_sampler_function))
+                self.kwargs.pop(user_input)
+
     def initialise_parameters(self):
 
         for key in self.priors:
@@ -195,14 +204,22 @@ class Nestle(Sampler):
 
     @kwargs.setter
     def kwargs(self, kwargs):
-        self.__kwargs = kwargs
+        self.__kwargs = dict(verbose=True, method='multi')
+        self.__kwargs.update(kwargs)
+
+        if 'npoints' not in self.__kwargs:
+            for equiv in ['nlive', 'nlives', 'n_live_points']:
+                if equiv in self.__kwargs:
+                    self.__kwargs['npoints'] = self.__kwargs.pop(equiv)
 
     def run_sampler(self):
         nestle = self.external_sampler
+        self.external_sampler_function = nestle.sample
         if self.kwargs.get('verbose', True):
             self.kwargs['callback'] = nestle.print_progress
+        self.verify_kwargs_against_external_sampler_function()
 
-        out = nestle.sample(
+        out = self.external_sampler_function(
             loglikelihood=self.log_likelihood,
             prior_transform=self.prior_transform,
             ndim=self.ndim, **self.kwargs)
@@ -250,9 +267,18 @@ class Pymultinest(Sampler):
         if self.__kwargs['outputfiles_basename'].endswith('/') is False:
             self.__kwargs['outputfiles_basename'] = '{}/'.format(
                 self.__kwargs['outputfiles_basename'])
+        if 'n_live_points' not in self.__kwargs:
+            for equiv in ['nlive', 'nlives', 'npoints', 'npoint']:
+                if equiv in self.__kwargs:
+                    self.__kwargs['n_live_points'] = self.__kwargs.pop(equiv)
 
     def run_sampler(self):
         pymultinest = self.external_sampler
+        self.external_sampler_function = pymultinest.run
+        self.verify_kwargs_against_external_sampler_function()
+        # Note: pymultinest.solve adds some extra steps, but underneath
+        # we are calling pymultinest.run - hence why it is used in checking
+        # the arguments.
         out = pymultinest.solve(
             LogLikelihood=self.log_likelihood, Prior=self.prior_transform,
             n_dims=self.ndim, **self.kwargs)
@@ -339,7 +365,11 @@ def run_sampler(likelihood, priors, label='label', outdir='outdir',
                                 **sampler_kwargs)
         result = sampler.run_sampler()
         result.noise_logz = likelihood.noise_log_likelihood()
-        result.log_bayes_factor = result.logz - result.noise_logz
+        if use_ratio:
+            result.log_bayes_factor = result.logz
+            result.logz = result.log_bayes_factor + result.noise_logz
+        else:
+            result.log_bayes_factor = result.logz - result.noise_logz
         result.injection_parameters = injection_parameters
         result.save_to_file(outdir=outdir, label=label)
         return result, sampler

peyote/source.py

 from __future__ import division, print_function
 
-try:
-    import lal
-except ImportError:
-    raise ImportWarning("You do not have lalsuite installed currently. You will not be able to use some of the "
-                        "prebuilt functions.")
 try:
     import lalsimulation as lalsim
 except ImportError:
     raise ImportWarning("You do not have lalsuite installed currently. You will not be able to use some of the "
                         "prebuilt functions.")
 
+from . import utils
+
 
 def lal_binary_black_hole(
-        frequency_array, mass_1, mass_2, luminosity_distance, spin11, spin12, spin13, spin21, spin22, spin23,
+        frequency_array, mass_1, mass_2, luminosity_distance, a_1, tilt_1, phi_1, a_2, tilt_2, phi_2,
         iota, phase, waveform_approximant, reference_frequency, ra, dec,
         geocent_time, psi):
     """ A Binary Black Hole waveform model using lalsimulation """
-    luminosity_distance = luminosity_distance * 1e6 * lal.PC_SI
     if mass_2 > mass_1:
         return {'plus': 0, 'cross': 0}
-    mass_1 = mass_1 * lal.MSUN_SI
-    mass_2 = mass_2 * lal.MSUN_SI
+
+    luminosity_distance = luminosity_distance * 1e6 * utils.parsec
+    mass_1 = mass_1 * utils.solar_mass
+    mass_2 = mass_2 * utils.solar_mass
+
+    spin_1x, spin_1y, spin_1z = utils.spherical_to_cartesian(a_1, tilt_1, phi_1)
+    spin_2x, spin_2y, spin_2z = utils.spherical_to_cartesian(a_2, tilt_2, phi_2)
 
     longitude_ascending_nodes = 0.0
     eccentricity = 0.0
     mean_per_ano = 0.0
 
-    waveform_dictionary = lal.CreateDict()
+    waveform_dictionary = None
 
     approximant = lalsim.GetApproximantFromString(waveform_approximant)
 
@@ -36,8 +37,8 @@ def lal_binary_black_hole(
     delta_frequency = frequency_array[1] - frequency_array[0]
 
     hplus, hcross = lalsim.SimInspiralChooseFDWaveform(
-        mass_1, mass_2, spin11, spin12, spin13, spin21, spin22,
-        spin23, luminosity_distance, iota, phase,
+        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,
         frequency_minimum, frequency_maximum, reference_frequency,
         waveform_dictionary, approximant)

peyote/utils.py

@@ -2,11 +2,12 @@ from __future__ import division
 import logging
 import os
 import numpy as np
-from astropy.time import Time
+from math import fmod
 
 # Constants
 speed_of_light = 299792458.0  # speed of light in m/s
-
+parsec = 3.085677581 * 1e16
+solar_mass = 1.98855 * 1e30
 
 def get_sampling_frequency(time_series):
     """
@@ -39,23 +40,32 @@ def ra_dec_to_theta_phi(ra, dec, gmst):
     return theta, phi
 
 
-def gps_time_to_gmst(time):
+def gps_time_to_gmst(gps_time):
     """
     Convert gps time to Greenwich mean sidereal time in radians
+
+    This method assumes a constant rotation rate of earth since 00:00:00, 1 Jan. 2000
+    A correction has been applied to give the exact correct value for 00:00:00, 1 Jan. 2018
+    Error accumulates at a rate of ~0.0001 radians/decade.
+
     Input:
     time - gps time
     Output:
     gmst - Greenwich mean sidereal time in radians
     """
-    gps_time = Time(time, format='gps', scale='utc')
-    gmst = gps_time.sidereal_time('mean', 'greenwich').value * np.pi / 12
+    omega_earth = 2 * np.pi * (1 / 365.2425 + 1) / 86400.
+    gps_2000 = 630720013.
+    gmst_2000 = (6 + 39. / 60 + 51.251406103947375 / 3600) * np.pi / 12
+    correction_2018 = -0.00017782487379358614
+    sidereal_time = omega_earth * (gps_time - gps_2000) + gmst_2000 + correction_2018
+    gmst = fmod(sidereal_time, 2 * np.pi)
     return gmst
 
 
 def create_fequency_series(sampling_frequency, duration):
     """
-    Create a frequency series with the correct length and spacing. 
-    
+    Create a frequency series with the correct length and spacing.
+
     :param sampling_frequency: sampling frequency
     :param duration: duration of data
     :return: frequencies, frequency series
@@ -81,8 +91,8 @@ def create_fequency_series(sampling_frequency, duration):
 def create_white_noise(sampling_frequency, duration):
     """
     Create white_noise which is then coloured by a given PSD
-    
-    
+
+
     :param sampling_frequency: sampling frequency
     :param duration: duration of data
     """
@@ -291,6 +301,19 @@ def get_progress_bar(module='tqdm'):
         return tqdm
 
 
+def spherical_to_cartesian(radius, theta, phi):
+    """
+    Convert from spherical coordinates to cartesian.
+
+    :param radius: radial coordinate
+    :param theta: axial coordinate
+    :param phi: azimuthal coordinate
+    :return cartesian: cartesian vector
+    """
+    cartesian = [radius * np.sin(theta) * np.cos(phi), radius * np.sin(theta) * np.sin(phi), radius * np.cos(theta)]
+    return cartesian
+
+
 def check_directory_exists_and_if_not_mkdir(directory):
     if not os.path.exists(directory):
         os.makedirs(directory)
@@ -298,3 +321,26 @@ def check_directory_exists_and_if_not_mkdir(directory):
     else:
         logging.debug('Directory {} exists'.format(directory))
 
+def inner_product(aa, bb, frequency, PSD):
+    '''
+    Calculate the inner product defined in the matched filter statistic
+
+    arguments:
+    aai, bb: single-sided Fourier transform, created, e.g., by the nfft function above
+    frequency: an array of frequencies associated with aa, bb, also returned by nfft
+    PSD: PSD object
+
+    Returns:
+    The matched filter inner product for aa and bb
+    '''
+    PSD_interp = PSD.power_spectral_density_interpolated(frequency)
+
+    # caluclate the inner product
+    integrand = np.conj(aa) * bb / PSD_interp
+
+    df = frequency[1] - frequency[0]
+    integral = np.sum(integrand) * df
+
+    product = 4. * np.real(integral)
+
+    return product

peyote/waveform_generator.py

@@ -23,11 +23,11 @@ class WaveformGenerator(object):
 
     """
 
-    def __init__(self, source_model, sampling_frequency=4096, time_duration=1,
+    def __init__(self, frequency_domain_source_model, sampling_frequency=4096, time_duration=1,
                  parameters=None):
         self.time_duration = time_duration
         self.sampling_frequency = sampling_frequency
-        self.source_model = source_model
+        self.source_model = frequency_domain_source_model
         self.parameters = parameters
         self.__frequency_array_updated = False
         self.__time_array_updated = False

test/noise_realisation_tests.py

+import numpy as np
+import unittest
+from context import peyote
+
+
+class TestNoiseRealisation(unittest.TestCase):
+    def test_averaged_noise(self):
+        time_duration = 1.
+        sampling_frequency = 4096.
+        factor = np.sqrt(2./time_duration)
+        navg = 1000
+        psd_avg = 0
+        H1 = peyote.detector.H1
+        for x in range(0, navg):
+            H1_hf_noise, frequencies = H1.power_spectral_density.get_noise_realisation(sampling_frequency,
+                                                                                       time_duration)
+            H1.set_data(sampling_frequency, time_duration, frequency_domain_strain=H1_hf_noise)
+            hf_tmp = H1.data
+            psd_avg += abs(hf_tmp)**2
+
+        psd_avg = psd_avg/navg
+        asd_avg = np.sqrt(abs(psd_avg))
+
+        a = H1.amplitude_spectral_density_array/factor
+        b = asd_avg
+        self.assertTrue(np.isclose(a[2]/b[2], 1.00, atol=1e-1))
+
+    def test_noise_normalisation(self):
+        time_duration = 1.
+        sampling_frequency = 4096.
+        number_of_samples = int(np.round(time_duration*sampling_frequency))
+        time_array = (1./sampling_frequency) * np.linspace(0, number_of_samples, number_of_samples)
+
+        # generate some toy-model signal for matched filtering SNR testing
+        navg = 10000
+        snr = np.zeros(navg)
+        mu = np.exp(-(time_array-time_duration/2.)**2 / (2.*0.1**2)) * np.sin(2*np.pi*100*time_array)
+        muf, frequency_array = peyote.utils.nfft(mu, sampling_frequency)
+        for x in range(0, navg):
+            H1 = peyote.detector.H1
+            H1_hf_noise, frequencies = H1.power_spectral_density.get_noise_realisation(sampling_frequency, time_duration)
+            H1.set_data(sampling_frequency, time_duration,frequency_domain_strain=H1_hf_noise)
+            hf_tmp = H1.data
+            Sh = H1.power_spectral_density
+            snr[x] = peyote.utils.inner_product(hf_tmp, muf, frequency_array, Sh) / np.sqrt(peyote.utils.inner_product(muf, muf, frequency_array, Sh))
+
+        self.assertTrue(np.isclose(np.std(snr), 1.00, atol=1e-2))

test/waveform_generator_tests.py

@@ -17,7 +17,8 @@ def gaussian_frequency_domain_strain_2(frequency_array, a, m, s, ra, dec, geocen
 class TestWaveformGeneratorInstantiationWithoutOptionalParameters(unittest.TestCase):
 
     def setUp(self):
-        self.waveform_generator = peyote.waveform_generator.WaveformGenerator(source_model=gaussian_frequency_domain_strain)
+        self.waveform_generator = peyote.waveform_generator.WaveformGenerator(
+            frequency_domain_source_model=gaussian_frequency_domain_strain)
         self.simulation_parameters = dict(amplitude=1e-21, mu=100, sigma=1,
                                      ra=1.375,
                                      dec=-1.2108,
@@ -51,7 +52,8 @@ class TestWaveformGeneratorInstantiationWithoutOptionalParameters(unittest.TestC
 class TestParameterSetter(unittest.TestCase):
 
     def setUp(self):
-        self.waveform_generator = peyote.waveform_generator.WaveformGenerator(source_model=gaussian_frequency_domain_strain)
+        self.waveform_generator = peyote.waveform_generator.WaveformGenerator(
+            frequency_domain_source_model=gaussian_frequency_domain_strain)
         self.simulation_parameters = dict(amplitude=1e-21, mu=100, sigma=1,
                                      ra=1.375,
                                      dec=-1.2108,
@@ -89,7 +91,8 @@ class TestParameterSetter(unittest.TestCase):
 class TestSourceModelSetter(unittest.TestCase):
 
     def setUp(self):
-        self.waveform_generator = peyote.waveform_generator.WaveformGenerator(source_model=gaussian_frequency_domain_strain)
+        self.waveform_generator = peyote.waveform_generator.WaveformGenerator(
+            frequency_domain_source_model=gaussian_frequency_domain_strain)
         self.waveform_generator.source_model = gaussian_frequency_domain_strain_2
         self.simulation_parameters = dict(a=1e-21, m=100, s=1,
                                      ra=1.375,

tutorials/BasicTutorial.py

@@ -11,12 +11,12 @@ sampling_frequency = 4096.
 injection_parameters = dict(
     mass_1=36.,
     mass_2=29.,
-    spin11=0,
-    spin12=0,
-    spin13=0,
-    spin21=0,
-    spin22=0,
-    spin23=0,
+    a_1=0,
+    a_2=0,
+    tilt_1=0,
+    tilt_2=0,
+    phi_1=0,
+    phi_2=0,
     luminosity_distance=100.,
     iota=0.4,
     phase=1.3,
@@ -32,7 +32,7 @@ injection_parameters = dict(
 waveform_generator = peyote.waveform_generator.WaveformGenerator(
     sampling_frequency=sampling_frequency,
     time_duration=time_duration,
-    source_model=peyote.source.lal_binary_black_hole,
+    frequency_domain_source_model=peyote.source.lal_binary_black_hole,
     parameters=injection_parameters)
 hf_signal = waveform_generator.frequency_domain_strain()
 

tutorials/GW150914.py

 from __future__ import division
+import matplotlib
+matplotlib.use('AGG')
 import peyote
+import numpy as np
 
 peyote.utils.setup_logger()
 
@@ -12,7 +15,7 @@ L1, sampling_frequency, time_duration = peyote.detector.get_inteferometer('L1',
 IFOs = [H1, L1]
 
 maximum_posterior_estimates = dict(
-    spin11=0, spin12=0, spin13=0, spin21=0, spin22=0, spin23=0,
+    a_1=0, a_2=0, tilt_1=0, tilt_2=0, phi_1=0, phi_2=0,
     luminosity_distance=410., iota=2.97305, phase=1.145,
     waveform_approximant='IMRPhenomPv2', reference_frequency=50., ra=1.375,
     dec=-1.2108, geocent_time=time_of_event, psi=2.659, mass_1=36, mass_2=29)
@@ -21,8 +24,8 @@ maximum_posterior_estimates = dict(
 prior = peyote.prior.parse_floats_to_fixed_priors(maximum_posterior_estimates)
 prior['ra'] = peyote.prior.create_default_prior(name='ra')
 prior['dec'] = peyote.prior.create_default_prior(name='dec')
-prior['psi'] = peyote.prior.create_default_prior(name='psi')
-prior['phase'] = peyote.prior.create_default_prior(name='phase')
+prior['psi'] = peyote.prior.Uniform(0, np.pi/2, 'psi')
+prior['phase'] = peyote.prior.Uniform(0, np.pi/2, 'phase')
 prior['iota'] = peyote.prior.create_default_prior(name='iota')
 prior['mass_1'] = peyote.prior.Uniform(10, 80, 'mass_1')
 prior['mass_2'] = peyote.prior.Uniform(10, 80, 'mass_2')
@@ -42,7 +45,6 @@ likelihood = peyote.likelihood.Likelihood(IFOs, waveformgenerator)
 # Run the sampler
 result, sampler = peyote.sampler.run_sampler(
     likelihood, prior, label='GW150914', sampler='pymultinest',
-    n_live_points=512, verbose=True, resume=False, outdir=outdir,
-    use_ratio=True)
+    npoints=1024, resume=False, outdir=outdir, use_ratio=True)
 truth = [maximum_posterior_estimates[x] for x in result.search_parameter_keys]
 sampler.plot_corner(truth=truth)

tutorials/custom_source_tutorial.py

@@ -73,7 +73,7 @@ injection_parameters = dict(amplitude=1e-21,
 sampling_parameters = peyote.prior.parse_floats_to_fixed_priors(injection_parameters)
 
 wg = peyote.waveform_generator.WaveformGenerator(
-     source_model=gaussian_frequency_domain_strain,
+     frequency_domain_source_model=gaussian_frequency_domain_strain,
      parameters=injection_parameters
      )
 

tutorials/make_standard_data.py

@@ -15,12 +15,12 @@ sampling_frequency = 4096.
 simulation_parameters = dict(
     mass_1=36.,
     mass_2=29.,
-    spin11=0,
-    spin12=0,
-    spin13=0,
-    spin21=0,
-    spin22=0,
-    spin23=0,
+    a_1=0,
+    a_2=0,
+    tilt_1=0,
+    tilt_2=0,
+    phi_1=0,
+    phi_2=0,
     luminosity_distance=100.,
     iota=0.4,
     phase=1.3,
@@ -32,7 +32,7 @@ simulation_parameters = dict(
     psi=2.659
 )
 #sampling_parameters = peyote.parameter.Parameter.parse_floats_to_parameters(simulation_parameters)
-waveform_generator = WaveformGenerator(source_model=peyote.source.lal_binary_black_hole,
+waveform_generator = WaveformGenerator(frequency_domain_source_model=peyote.source.lal_binary_black_hole,
                                        sampling_frequency=sampling_frequency,
                                        time_duration=time_duration,
                                        parameters=simulation_parameters)