Merge branch 'master' of git.ligo.org:Monash/tupak

.gitlab-ci.yml

@@ -22,11 +22,11 @@ exitcode-jessie:
     - pip install coverage
     - pip install coverage-badge
     - coverage erase
-    - coverage run -a test/prior_tests.py
-    - coverage run -a test/tests.py
-    - coverage run -a test/waveform_generator_tests.py
-    - coverage run -a test/noise_realisation_tests.py
-    - coverage html --include=tupak/*
+    - coverage run --include=tupak/* -a test/prior_tests.py
+    - coverage run --include=tupak/* -a test/tests.py
+    - coverage run --include=tupak/* -a test/waveform_generator_tests.py
+    - coverage run --include=tupak/* -a test/noise_realisation_tests.py
+    - coverage html
     - coverage-badge -o coverage_badge.svg -f
   artifacts:
     paths:

README.md

 [![pipeline status](https://git.ligo.org/Monash/tupak/badges/master/pipeline.svg)](https://git.ligo.org/Monash/tupak/commits/master)
-[![coverage report](https://monash.docs.ligo.org/tupak/coverage_report.svg)](
+[![coverage report](https://monash.docs.ligo.org/tupak/coverage_badge.svg)](
 https://monash.docs.ligo.org/tupak/)
 
 # Tupak

examples/other_examples/alternative_likelihoods.py

@@ -20,7 +20,7 @@ outdir = 'outdir'
 # use of the `tupak` waveform_generator to make the signal (more on this later)
 # But, one could make this work without the waveform generator.
 
-class GaussianGravitationalWaveTransient():
+class GaussianLikelihood():
     def __init__(self, x, y, waveform_generator):
         self.x = x
         self.y = y
@@ -76,7 +76,7 @@ waveform_generator = tupak.waveform_generator.WaveformGenerator(time_duration=ti
 
 # Now lets instantiate a version of out GravitationalWaveTransient, giving it the time, data
 # and waveform_generator
-likelihood = GaussianGravitationalWaveTransient(time, data, waveform_generator)
+likelihood = GaussianLikelihood(time, data, waveform_generator)
 
 # From hereon, the syntax is exactly equivalent to other tupak examples
 # We make a prior

test/prior_tests.py

@@ -198,6 +198,7 @@ class TestFillPrior(unittest.TestCase):
     def setUp(self):
         self.likelihood = Mock()
         self.likelihood.parameters = dict(a=0, b=0, c=0, d=0, asdf=0, ra=1)
+        self.likelihood.non_standard_sampling_parameter_keys = dict(t=8)
         self.priors = dict(a=1, b=1.1, c='string', d=tupak.prior.Uniform(0, 1))
         self.priors = tupak.prior.fill_priors(self.priors, self.likelihood)
 

test/tests.py

@@ -57,8 +57,10 @@ class Test(unittest.TestCase):
         priors['luminosity_distance'] = tupak.prior.Uniform(
             name='luminosity_distance', minimum=dL - 10, maximum=dL + 10)
 
-        result = tupak.sampler.run_sampler(likelihood, priors, sampler='nestle')
-        self.assertAlmostEqual(np.mean(result.samples), dL, delta=np.std(result.samples))
+        result = tupak.sampler.run_sampler(
+            likelihood, priors, sampler='nestle', verbose=False, npoints=100)
+        self.assertAlmostEqual(np.mean(result.samples), dL,
+                               delta=np.std(result.samples))
 
 
 if __name__ == '__main__':

tupak/likelihood.py

@@ -11,7 +11,20 @@ import tupak
 import logging
 
 
-class GravitationalWaveTransient(object):
+class Likelihood(object):
+    """ Empty likelihood class to be subclassed by other likelihoods """
+
+    def log_likelihood():
+        return np.nan
+
+    def noise_log_likelihood():
+        return np.nan
+
+    def log_likelihood_ratio(self):
+        return self.log_likelihood() - self.noise_log_likelihood()
+
+
+class GravitationalWaveTransient(Likelihood):
     """ A gravitational-wave transient likelihood object
 
     This is the usual likelihood object to use for transient gravitational
@@ -114,7 +127,7 @@ class GravitationalWaveTransient(object):
         return log_l.real
 
     def log_likelihood(self):
-        return self.log_likelihood() + self.noise_log_likelihood()
+        return self.log_likelihood_ratio() + self.noise_log_likelihood()
 
     def setup_distance_marginalization(self):
         if 'luminosity_distance' not in self.prior.keys():
@@ -136,6 +149,62 @@ class GravitationalWaveTransient(object):
                                                  bounds_error=False, fill_value=-np.inf)
 
 
+class BasicGravitationalWaveTransient(Likelihood):
+    """ A basic gravitaitonal wave transient likelihood
+
+    The simplest frequency-domain gravitational wave transient likelihood. Does
+    not include distance/phase marginalization.
+
+    Parameters
+    ----------
+    interferometers: list
+        A list of `tupak.detector.Interferometer` instances - contains the
+        detector data and power spectral densities
+    waveform_generator: `tupak.waveform_generator.WaveformGenerator`
+        An object which computes the frequency-domain strain of the signal,
+        given some set of parameters
+
+    Returns
+    -------
+    Likelihood: `tupak.likelihood.Likelihood`
+        A likehood object, able to compute the likelihood of the data given
+        some model parameters
+
+    """
+    def __init__(self, interferometers, waveform_generator):
+        self.interferometers = interferometers
+        self.waveform_generator = waveform_generator
+
+    def noise_log_likelihood(self):
+        log_l = 0
+        for interferometer in self.interferometers:
+            log_l -= 2. / self.waveform_generator.time_duration * np.sum(
+                abs(interferometer.data) ** 2 /
+                interferometer.power_spectral_density_array)
+        return log_l.real
+
+    def log_likelihood(self):
+        log_l = 0
+        waveform_polarizations = self.waveform_generator.frequency_domain_strain()
+        if waveform_polarizations is None:
+            return np.nan_to_num(-np.inf)
+        for interferometer in self.interferometers:
+            log_l += self.log_likelihood_interferometer(
+                waveform_polarizations, interferometer)
+        return log_l.real
+
+    def log_likelihood_interferometer(self, waveform_polarizations,
+                                      interferometer):
+        signal_ifo = interferometer.get_detector_response(
+            waveform_polarizations, self.waveform_generator.parameters)
+
+        log_l = - 2. / self.waveform_generator.time_duration * np.vdot(
+            interferometer.data - signal_ifo,
+            (interferometer.data - signal_ifo)
+            / interferometer.power_spectral_density_array)
+        return log_l.real
+
+
 def get_binary_black_hole_likelihood(interferometers):
     """ A rapper to quickly set up a likelihood for BBH parameter estimation
 
@@ -157,3 +226,4 @@ def get_binary_black_hole_likelihood(interferometers):
         parameters={'waveform_approximant': 'IMRPhenomPv2', 'reference_frequency': 50})
     likelihood = tupak.likelihood.GravitationalWaveTransient(interferometers, waveform_generator)
     return likelihood
+

tupak/prior.py

@@ -433,7 +433,7 @@ def parse_keys_to_parameters(keys):
     return parameters
 
 
-def fill_priors(prior, likelihood, parameters=None):
+def fill_priors(prior, likelihood):
     """
     Fill dictionary of priors based on required parameters of likelihood
 
@@ -446,9 +446,9 @@ def fill_priors(prior, likelihood, parameters=None):
         dictionary of prior objects and floats
     likelihood: tupak.likelihood.GravitationalWaveTransient instance
         Used to infer the set of parameters to fill the prior with
-    parameters: list
-        list of parameters to be sampled in, this can override the default
-        priors for the waveform generator
+
+    Note: if `likelihood` has `non_standard_sampling_parameter_keys`, then this
+    will set-up default priors for those as well.
 
     Returns
     -------
@@ -470,8 +470,8 @@ def fill_priors(prior, likelihood, parameters=None):
 
     missing_keys = set(likelihood.parameters) - set(prior.keys())
 
-    if parameters is not None:
-        for parameter in parameters:
+    if getattr(likelihood, 'non_standard_sampling_parameter_keys', None) is not None:
+        for parameter in likelihood.non_standard_sampling_parameter_keys:
             prior[parameter] = create_default_prior(parameter)
 
     for missing_key in missing_keys:

tupak/result.py

@@ -67,6 +67,9 @@ class Result(dict):
         else:
             return ''
 
+    def get_result_dictionary(self):
+        return dict(self)
+
     def save_to_file(self, outdir, label):
         """ Writes the Result to a deepdish h5 file """
         file_name = result_file_name(outdir, label)
@@ -80,7 +83,7 @@ class Result(dict):
 
         logging.info("Saving result to {}".format(file_name))
         try:
-            deepdish.io.save(file_name, self)
+            deepdish.io.save(file_name, self.get_result_dictionary())
         except Exception as e:
             logging.error(
                 "\n\n Saving the data has failed with the following message:\n {} \n\n"

tupak/sampler.py

@@ -245,8 +245,9 @@ class Nestle(Sampler):
     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
+        if 'verbose' in self.kwargs:
+            if self.kwargs['verbose']:
+                self.kwargs['callback'] = nestle.print_progress
             self.kwargs.pop('verbose')
         self.verify_kwargs_against_external_sampler_function()
 
@@ -444,7 +445,7 @@ def run_sampler(likelihood, priors=None, label='label', outdir='outdir',
 
     if priors is None:
         priors = dict()
-    priors = fill_priors(priors, likelihood, parameters=likelihood.non_standard_sampling_parameter_keys)
+    priors = fill_priors(priors, likelihood)
     tupak.prior.write_priors_to_file(priors, outdir)
 
     if implemented_samplers.__contains__(sampler.title()):
@@ -468,7 +469,7 @@ def run_sampler(likelihood, priors=None, label='label', outdir='outdir',
             result.injection_parameters = injection_parameters
             if conversion_function is not None:
                 conversion_function(result.injection_parameters)
-        result.fixed_parameter_keys = [key for key in priors if isinstance(key, prior.DeltaFunction)]
+        result.fixed_parameter_keys = sampler.fixed_parameter_keys
         # result.prior = prior  # Removed as this breaks the saving of the data
         result.samples_to_data_frame(likelihood=likelihood, priors=priors, conversion_function=conversion_function)
         result.kwargs = sampler.kwargs

tupak/waveform_generator.py

@@ -89,7 +89,7 @@ class WaveformGenerator(object):
         if self.parameter_conversion is not None:
             for key in added_keys:
                 self.parameters.pop(key)
-        return model_time_series[key]
+        return model_time_series
 
     @property
     def frequency_array(self):