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

docs/likelihood.txt

@@ -167,11 +167,55 @@ In the last example, we considered only cases with known noise (e.g., a
 prespecified standard deviation. We now present a general function which can
 handle unknown noise (in which case you need to specify a prior for
 :math:`\sigma`, or known noise (in which case you pass the known noise in when
-instatiating the likelihood
-
-.. literalinclude:: ../examples/other_examples/linear_regression_unknown_noise.py
-   :lines: 52-94
-
+instatiating the likelihood::
+
+  class GaussianLikelihood(tupak.Likelihood):
+      def __init__(self, x, y, function, sigma=None):
+          """
+          A general Gaussian likelihood for known or unknown noise - the model
+          parameters are inferred from the arguments of function
+
+          Parameters
+          ----------
+          x, y: array_like
+              The data to analyse
+          function:
+              The python function to fit to the data. Note, this must take the
+              dependent variable as its first argument. The other arguments are
+              will require a prior and will be sampled over (unless a fixed
+              value is given).
+          sigma: None, float, array_like
+              If None, the standard deviation of the noise is unknown and will be
+              estimated (note: this requires a prior to be given for sigma). If
+              not None, this defined the standard-deviation of the data points.
+              This can either be a single float, or an array with length equal
+              to that for `x` and `y`.
+          """
+          self.x = x
+          self.y = y
+          self.N = len(x)
+          self.sigma = sigma
+          self.function = function
+
+          # These lines of code infer the parameters from the provided function
+          parameters = inspect.getargspec(function).args
+          parameters.pop(0)
+          self.parameters = dict.fromkeys(parameters)
+          self.function_keys = self.parameters.keys()
+          if self.sigma is None:
+              self.parameters['sigma'] = None
+
+      def log_likelihood(self):
+          sigma = self.parameters.get('sigma', self.sigma)
+          model_parameters = {k: self.parameters[k] for k in self.function_keys}
+          res = self.y - self.function(self.x, **model_parameters)
+          return -0.5 * (np.sum((res / sigma)**2)
+                         + self.N*np.log(2*np.pi*sigma**2))
+
+We provide this general-purpose class as part of tupak:
+
+.. autoclass:: tupak.core.likelihood.GaussianLikelihood
+   :members:
 
 An example using this likelihood can be found `on this page <https://git.ligo.org/Monash/tupak/blob/master/examples/other_examples/linear_regression_unknown_noise.py>`_.
 

examples/injection_examples/basic_tutorial.py

@@ -74,4 +74,4 @@ result = tupak.run_sampler(likelihood=likelihood, priors=priors, sampler='dynest
 
 # make some plots of the outputs
 result.plot_corner()
-print(result)
+

examples/injection_examples/basic_tutorial_dist_time_phase_marg.py

@@ -71,4 +71,4 @@ result = tupak.run_sampler(likelihood=likelihood, priors=priors, sampler='dynest
 
 # make some plots of the outputs
 result.plot_corner()
-print(result)
+

examples/injection_examples/basic_tutorial_time_phase_marg.py

@@ -70,4 +70,4 @@ result = tupak.run_sampler(likelihood=likelihood, priors=priors, sampler='dynest
 
 # make some plots of the outputs
 result.plot_corner()
-print(result)
+

examples/injection_examples/change_sampled_parameters.py

@@ -58,4 +58,4 @@ result = tupak.core.sampler.run_sampler(likelihood=likelihood, priors=priors, sa
                                         injection_parameters=injection_parameters, label='DifferentParameters',
                                         outdir=outdir, conversion_function=tupak.gw.conversion.generate_all_bbh_parameters)
 result.plot_corner()
-print(result)
+

examples/injection_examples/create_your_own_source_model.py

@@ -51,4 +51,4 @@ result = tupak.core.sampler.run_sampler(
     likelihood, prior, sampler='dynesty', outdir=outdir, label=label,
     resume=False, sample='unif', injection_parameters=injection_parameters)
 result.plot_corner()
-print(result)
+

examples/injection_examples/create_your_own_time_domain_source_model.py

@@ -72,4 +72,4 @@ result = tupak.core.sampler.run_sampler(likelihood, prior, sampler='dynesty', np
                                         outdir=outdir, label=label)
 
 result.plot_corner()
-print(result)
+

examples/injection_examples/how_to_specify_the_prior.py

@@ -64,4 +64,4 @@ likelihood = tupak.GravitationalWaveTransient(interferometers=IFOs, waveform_gen
 result = tupak.run_sampler(likelihood=likelihood, priors=priors, sampler='dynesty',
                            injection_parameters=injection_parameters, outdir=outdir, label='specify_prior')
 result.plot_corner()
-print(result)
+

examples/injection_examples/marginalized_likelihood.py

@@ -50,4 +50,4 @@ likelihood = tupak.GravitationalWaveTransient(
 result = tupak.run_sampler(likelihood=likelihood, priors=priors, sampler='dynesty',
                            injection_parameters=injection_parameters, outdir=outdir, label='BasicTutorial')
 result.plot_corner()
-print(result)
+

examples/injection_examples/sine_gaussian_example.py

@@ -63,7 +63,7 @@ result = tupak.core.sampler.run_sampler(likelihood=likelihood, priors=priors, sa
 
 # make some plots of the outputs
 result.plot_corner()
-print(result)
+
 
 
 

examples/open_data_examples/GW150914.py

@@ -51,4 +51,4 @@ likelihood = tupak.gw.likelihood.GravitationalWaveTransient(interferometers, wav
 result = tupak.run_sampler(likelihood, prior, sampler='dynesty',
                            outdir=outdir, label=label)
 result.plot_corner()
-print(result)
+

examples/other_examples/gaussian_example.py

@@ -47,10 +47,10 @@ class SimpleGaussianLikelihood(tupak.Likelihood):
 likelihood = SimpleGaussianLikelihood(data)
 priors = dict(mu=tupak.core.prior.Uniform(0, 5, 'mu'),
               sigma=tupak.core.prior.Uniform(0, 10, 'sigma'))
+priors['mu'] = 1
 
 # And run sampler
 result = tupak.run_sampler(
     likelihood=likelihood, priors=priors, sampler='dynesty', npoints=500,
     walks=10, outdir=outdir, label=label)
 result.plot_corner()
-print(result)

examples/other_examples/hyper_parameter_example.py

@@ -108,7 +108,7 @@ fig2.savefig('outdir/hyper_parameter_combined_posteriors.png')
 run_prior = tupak.core.prior.Uniform(minimum=-10, maximum=10, name='mu_c0')
 hyper_prior = tupak.core.prior.Gaussian(mu=0, sigma=1, name='hyper')
 
-hp_likelihood = tupak.gw.likelihood.HyperparameterLikelihood(
+hp_likelihood = tupak.core.likelihood.HyperparameterLikelihood(
         samples, hyper_prior, run_prior)
 
 hp_priors = dict(

examples/other_examples/linear_regression.py

@@ -101,4 +101,3 @@ result = tupak.run_sampler(
     walks=10, injection_parameters=injection_parameters, outdir=outdir,
     label=label)
 result.plot_corner()
-print(result)

examples/other_examples/linear_regression_unknown_noise.py

@@ -9,7 +9,6 @@ from __future__ import division
 import tupak
 import numpy as np
 import matplotlib.pyplot as plt
-import inspect
 
 # A few simple setup steps
 tupak.core.utils.setup_logger()
@@ -45,60 +44,11 @@ ax.set_ylabel('y')
 ax.legend()
 fig.savefig('{}/{}_data.png'.format(outdir, label))
 
+# Now lets instantiate the built-in GaussianLikelihood, giving it
+# the time, data and signal model. Note that, because we do not give it the
+# parameter, sigma is unknown and marginalised over during the sampling
+likelihood = tupak.core.likelihood.GaussianLikelihood(time, data, model)
 
-# Parameter estimation: we now define a Gaussian Likelihood class relevant for
-# our model.
-
-class GaussianLikelihood(tupak.Likelihood):
-    def __init__(self, x, y, function, sigma=None):
-        """
-        A general Gaussian likelihood for known or unknown noise - the model
-        parameters are inferred from the arguments of function
-
-        Parameters
-        ----------
-        x, y: array_like
-            The data to analyse
-        function:
-            The python function to fit to the data. Note, this must take the
-            dependent variable as its first argument. The other arguments are
-            will require a prior and will be sampled over (unless a fixed
-            value is given).
-        sigma: None, float, array_like
-            If None, the standard deviation of the noise is unknown and will be
-            estimated (note: this requires a prior to be given for sigma). If
-            not None, this defined the standard-deviation of the data points.
-            This can either be a single float, or an array with length equal
-            to that for `x` and `y`.
-        """
-        self.x = x
-        self.y = y
-        self.N = len(x)
-        self.sigma = sigma
-        self.function = function
-
-        # These lines of code infer the parameters from the provided function
-        parameters = inspect.getargspec(function).args
-        parameters.pop(0)
-        self.parameters = dict.fromkeys(parameters)
-        self.function_keys = self.parameters.keys()
-        if self.sigma is None:
-            self.parameters['sigma'] = None
-
-    def log_likelihood(self):
-        sigma = self.parameters.get('sigma', self.sigma)
-        model_parameters = {k: self.parameters[k] for k in self.function_keys}
-        res = self.y - self.function(self.x, **model_parameters)
-        return -0.5 * (np.sum((res / sigma)**2)
-                       + self.N*np.log(2*np.pi*sigma**2))
-
-
-# Now lets instantiate a version of our GaussianLikelihood, giving it
-# the time, data and signal model
-likelihood = GaussianLikelihood(time, data, model)
-
-# From hereon, the syntax is exactly equivalent to other tupak examples
-# We make a prior
 priors = {}
 priors['m'] = tupak.core.prior.Uniform(0, 5, 'm')
 priors['c'] = tupak.core.prior.Uniform(-2, 2, 'c')
@@ -110,4 +60,3 @@ result = tupak.run_sampler(
     walks=10, injection_parameters=injection_parameters, outdir=outdir,
     label=label)
 result.plot_corner()
-print(result)

test/detector_tests.py

@@ -305,16 +305,16 @@ class TestDetector(unittest.TestCase):
             self.ifo.xarm_tilt = 0
             self.ifo.xarm_azimuth = 0
             self.ifo.yarm_tilt = 0
-            self.ifo.yarm_azimuth = np.pi
+            self.ifo.yarm_azimuth = 90
             self.assertAlmostEqual(self.ifo.unit_vector_along_arm('x'), 1)
 
     def test_unit_vector_along_arm_y(self):
         with mock.patch('numpy.array') as m:
             m.return_value = 1
             self.ifo.xarm_tilt = 0
-            self.ifo.xarm_azimuth = 0
+            self.ifo.xarm_azimuth = 90
             self.ifo.yarm_tilt = 0
-            self.ifo.yarm_azimuth = np.pi
+            self.ifo.yarm_azimuth = 180
             self.assertAlmostEqual(self.ifo.unit_vector_along_arm('y'), -1)
 
     def test_set_data_raises_value_error(self):

tupak/core/likelihood.py

 from __future__ import division, print_function
 
+import inspect
+import logging
 import numpy as np
 
 try:
@@ -24,3 +26,96 @@ class Likelihood(object):
         return self.log_likelihood() - self.noise_log_likelihood()
 
 
+class GaussianLikelihood(Likelihood):
+    def __init__(self, x, y, function, sigma=None):
+        """
+        A general Gaussian likelihood for known or unknown noise - the model
+        parameters are inferred from the arguments of function
+
+        Parameters
+        ----------
+        x, y: array_like
+            The data to analyse
+        function:
+            The python function to fit to the data. Note, this must take the
+            dependent variable as its first argument. The other arguments are
+            will require a prior and will be sampled over (unless a fixed
+            value is given).
+        sigma: None, float, array_like
+            If None, the standard deviation of the noise is unknown and will be
+            estimated (note: this requires a prior to be given for sigma). If
+            not None, this defined the standard-deviation of the data points.
+            This can either be a single float, or an array with length equal
+            to that for `x` and `y`.
+        """
+        self.x = x
+        self.y = y
+        self.N = len(x)
+        self.sigma = sigma
+        self.function = function
+
+        # These lines of code infer the parameters from the provided function
+        parameters = inspect.getargspec(function).args
+        parameters.pop(0)
+        self.parameters = dict.fromkeys(parameters)
+        self.function_keys = self.parameters.keys()
+        if self.sigma is None:
+            self.parameters['sigma'] = None
+
+    def log_likelihood(self):
+        sigma = self.parameters.get('sigma', self.sigma)
+        model_parameters = {k: self.parameters[k] for k in self.function_keys}
+        res = self.y - self.function(self.x, **model_parameters)
+        return -0.5 * (np.sum((res / sigma)**2)
+                       + self.N*np.log(2*np.pi*sigma**2))
+
+
+class HyperparameterLikelihood(Likelihood):
+    """ A likelihood for infering hyperparameter posterior distributions
+
+    See Eq. (1) of https://arxiv.org/abs/1801.02699 for a definition.
+
+    Parameters
+    ----------
+    samples: list
+        An N-dimensional list of individual sets of samples. Each set may have
+        a different size.
+    hyper_prior: `tupak.prior.Prior`
+        A prior distribution with a `parameters` argument pointing to the
+        hyperparameters to infer from the samples. These may need to be
+        initialized to any arbitrary value, but this will not effect the
+        result.
+    run_prior: `tupak.prior.Prior`
+        The prior distribution used in the inidivudal inferences which resulted
+        in the set of samples.
+
+    """
+
+    def __init__(self, samples, hyper_prior, run_prior):
+        Likelihood.__init__(self, parameters=hyper_prior.__dict__)
+        self.samples = samples
+        self.hyper_prior = hyper_prior
+        self.run_prior = run_prior
+        if hasattr(hyper_prior, 'lnprob') and hasattr(run_prior, 'lnprob'):
+            logging.info("Using log-probabilities in likelihood")
+            self.log_likelihood = self.log_likelihood_using_lnprob
+        else:
+            logging.info("Using probabilities in likelihood")
+            self.log_likelihood = self.log_likelihood_using_prob
+
+    def log_likelihood_using_lnprob(self):
+        L = []
+        self.hyper_prior.__dict__.update(self.parameters)
+        for samp in self.samples:
+            f = self.hyper_prior.lnprob(samp) - self.run_prior.lnprob(samp)
+            L.append(logsumexp(f))
+        return np.sum(L)
+
+    def log_likelihood_using_prob(self):
+        L = []
+        self.hyper_prior.__dict__.update(self.parameters)
+        for samp in self.samples:
+            L.append(
+                np.sum(self.hyper_prior.prob(samp) /
+                       self.run_prior.prob(samp)))
+        return np.sum(np.log(L))

tupak/core/prior.py

@@ -94,10 +94,10 @@ class PriorSet(dict):
                 continue
             elif isinstance(self[key], float) or isinstance(self[key], int):
                 self[key] = DeltaFunction(self[key])
-                logging.info(
+                logging.debug(
                     "{} converted to delta function prior.".format(key))
             else:
-                logging.info(
+                logging.debug(
                     "{} cannot be converted to delta function prior."
                     .format(key))
 
@@ -139,7 +139,7 @@ class PriorSet(dict):
             if isinstance(self[key], Prior):
                 samples[key] = self[key].sample(size=size)
             else:
-                logging.info('{} not a known prior.'.format(key))
+                logging.debug('{} not a known prior.'.format(key))
         return samples
 
     def prob(self, sample):
@@ -180,7 +180,7 @@ def create_default_prior(name, default_priors_file=None):
     """
 
     if default_priors_file is None:
-        logging.info(
+        logging.debug(
             "No prior file given.")
         prior = None
     else:
@@ -188,7 +188,7 @@ def create_default_prior(name, default_priors_file=None):
         if name in default_priors.keys():
             prior = default_priors[name]
         else:
-            logging.info(
+            logging.debug(
                 "No default prior found for variable {}.".format(name))
             prior = None
     return prior
@@ -579,7 +579,7 @@ class Interped(Prior):
 
     def __initialize_attributes(self):
         if np.trapz(self.yy, self.xx) != 1:
-            logging.info('Supplied PDF for {} is not normalised, normalising.'.format(self.name))
+            logging.debug('Supplied PDF for {} is not normalised, normalising.'.format(self.name))
         self.yy /= np.trapz(self.yy, self.xx)
         self.YY = cumtrapz(self.yy, self.xx, initial=0)
         # Need last element of cumulative distribution to be exactly one.

tupak/core/result.py

@@ -91,12 +91,12 @@ class Result(dict):
         if os.path.isdir(self.outdir) is False:
             os.makedirs(self.outdir)
         if os.path.isfile(file_name):
-            logging.info(
+            logging.debug(
                 'Renaming existing file {} to {}.old'.format(file_name,
                                                              file_name))
             os.rename(file_name, file_name + '.old')
 
-        logging.info("Saving result to {}".format(file_name))
+        logging.debug("Saving result to {}".format(file_name))
         try:
             deepdish.io.save(file_name, self.get_result_dictionary())
         except Exception as e:
@@ -183,7 +183,7 @@ class Result(dict):
 
         if save:
             filename = '{}/{}_corner.png'.format(self.outdir, self.label)
-            logging.info('Saving corner plot to {}'.format(filename))
+            logging.debug('Saving corner plot to {}'.format(filename))
             fig.savefig(filename, dpi=dpi)
 
         return fig

tupak/core/sampler.py

@@ -130,7 +130,7 @@ class Sampler(object):
 
         logging.info("Search parameters:")
         for key in self.__search_parameter_keys:
-            logging.info('  {} ~ {}'.format(key, self.priors[key]))
+            logging.info('  {} = {}'.format(key, self.priors[key]))
         for key in self.__fixed_parameter_keys:
             logging.info('  {} = {}'.format(key, self.priors[key].peak))
 
@@ -156,7 +156,8 @@ class Sampler(object):
             t1 = datetime.datetime.now()
             self.likelihood.log_likelihood()
             logging.info(
-                "Single likelihood eval. took {:.3e} s".format((datetime.datetime.now() - t1).total_seconds()))
+                "Single likelihood evaluation took {:.3e} s"
+                .format((datetime.datetime.now() - t1).total_seconds()))
         except TypeError as e:
             raise TypeError(
                 "Likelihood evaluation failed with message: \n'{}'\n"
@@ -212,9 +213,9 @@ class Sampler(object):
         """
         draw = np.array(list(self.priors.sample_subset(self.__search_parameter_keys).values()))
         if np.isinf(self.log_likelihood(draw)):
-            logging.info('Prior draw {} has inf likelihood'.format(draw))
+            logging.warning('Prior draw {} has inf likelihood'.format(draw))
         if np.isinf(self.log_prior(draw)):
-            logging.info('Prior draw {} has inf prior'.format(draw))
+            logging.warning('Prior draw {} has inf prior'.format(draw))
         return draw
 
     def _run_external_sampler(self):
@@ -396,7 +397,7 @@ class Dynesty(Sampler):
 
     def generate_trace_plots(self, dynesty_results):
         filename = '{}/{}_trace.png'.format(self.outdir, self.label)
-        logging.info("Writing trace plot to {}".format(filename))
+        logging.debug("Writing trace plot to {}".format(filename))
         from dynesty import plotting as dyplot
         fig, axes = dyplot.traceplot(dynesty_results,
                                      labels=self.result.parameter_labels)
@@ -504,7 +505,7 @@ class Ptemcee(Sampler):
 
         fig.tight_layout()
         filename = '{}/{}_walkers.png'.format(self.outdir, self.label)
-        logging.info('Saving walkers plot to {}'.format('filename'))
+        logging.debug('Saving walkers plot to {}'.format('filename'))
         fig.savefig(filename)
 
 
@@ -580,7 +581,7 @@ def run_sampler(likelihood, priors=None, label='label', outdir='outdir',
                                 **kwargs)
 
         if sampler.cached_result:
-            logging.info("Using cached result")
+            logging.warning("Using cached result")
             return sampler.cached_result
 
         start_time = datetime.datetime.now()
@@ -592,7 +593,6 @@ def run_sampler(likelihood, priors=None, label='label', outdir='outdir',
 
         end_time = datetime.datetime.now()
         result.sampling_time = (end_time - start_time).total_seconds()
-        logging.info('')
         logging.info('Sampling time: {}'.format(end_time - start_time))
 
         if sampler.use_ratio:
@@ -615,6 +615,8 @@ def run_sampler(likelihood, priors=None, label='label', outdir='outdir',
         result.save_to_file()
         if plot:
             result.plot_corner()
+        logging.info("Sampling finished, results saved to {}/".format(outdir))
+        logging.info("Summary of results:\n{}".format(result))
         return result
     else:
         raise ValueError(