Skip to content
Snippets Groups Projects
Commit 6334cda6 authored by Colm Talbot's avatar Colm Talbot
Browse files

Make sure core examples run without error

parent e04d1190
No related branches found
No related tags found
1 merge request!1081Make sure core examples run without error
......@@ -406,7 +406,7 @@ class Multinomial(Likelihood):
Likelihood for system with N discrete possibilities.
"""
def __init__(self, data, n_dimensions, label="parameter_"):
def __init__(self, data, n_dimensions, base="parameter_"):
"""
Parameters
......@@ -415,19 +415,21 @@ class Multinomial(Likelihood):
The number of objects in each class
n_dimensions: int
The number of classes
base: str
The base of the parameter labels
"""
self.data = np.array(data)
self._total = np.sum(self.data)
super(Multinomial, self).__init__(dict())
self.n = n_dimensions
self.label = label
self.base = base
self._nll = None
def log_likelihood(self):
"""
Since n - 1 parameters are sampled, the last parameter is 1 - the rest
"""
probs = [self.parameters[self.label + str(ii)]
probs = [self.parameters[self.base + str(ii)]
for ii in range(self.n - 1)]
probs.append(1 - sum(probs))
return self._multinomial_ln_pdf(probs=probs)
......
......@@ -8,8 +8,6 @@ would give equivalent results as using the pre-defined 'Gaussian Likelihood'
"""
import inspect
import bilby
import matplotlib.pyplot as plt
import numpy as np
......@@ -52,8 +50,8 @@ fig.savefig("{}/{}_data.png".format(outdir, label))
# Parameter estimation: we now define a Gaussian Likelihood class relevant for
# our model.
class GaussianLikelihoodPyMC3(bilby.Likelihood):
def __init__(self, x, y, sigma, function):
class GaussianLikelihoodPyMC3(bilby.core.likelihood.GaussianLikelihood):
def __init__(self, x, y, sigma, func):
"""
A general Gaussian likelihood - the parameters are inferred from the
arguments of function
......@@ -64,22 +62,13 @@ class GaussianLikelihoodPyMC3(bilby.Likelihood):
The data to analyse
sigma: float
The standard deviation of the noise
function:
func:
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).
"""
self.x = x
self.y = y
self.sigma = sigma
self.N = len(x)
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)
super(GaussianLikelihoodPyMC3, self).__init__(x=x, y=y, func=func, sigma=sigma)
def log_likelihood(self, sampler=None):
"""
......@@ -88,12 +77,15 @@ class GaussianLikelihoodPyMC3(bilby.Likelihood):
sampler: :class:`bilby.core.sampler.Pymc3`
A Sampler object must be passed containing the prior distributions
and PyMC3 :class:`~pymc3.Model` to use as a context manager.
If this is not passed, the super class is called and the regular
likelihood is evaluated.
"""
from bilby.core.sampler import Pymc3
if not isinstance(sampler, Pymc3):
raise ValueError("Sampler is not a bilby Pymc3 sampler object")
print(sampler, type(sampler))
return super(GaussianLikelihoodPyMC3, self).log_likelihood()
if not hasattr(sampler, "pymc3_model"):
raise AttributeError("Sampler has not PyMC3 model attribute")
......@@ -114,12 +106,11 @@ likelihood = GaussianLikelihoodPyMC3(time, data, sigma, model)
# Define a custom prior for one of the parameter for use with PyMC3
class PriorPyMC3(bilby.core.prior.Prior):
class PyMC3UniformPrior(bilby.core.prior.Uniform):
def __init__(self, minimum, maximum, name=None, latex_label=None):
"""
Uniform prior with bounds (should be equivalent to bilby.prior.Uniform)
"""
bilby.core.prior.Prior.__init__(
self, name, latex_label, minimum=minimum, maximum=maximum
)
......@@ -128,12 +119,15 @@ class PriorPyMC3(bilby.core.prior.Prior):
"""
Change ln_prob method to take in a Sampler and return a PyMC3
distribution.
If the passed argument is not a `Pymc3` sampler, assume that it is a
float or array to be passed to the superclass.
"""
from bilby.core.sampler import Pymc3
if not isinstance(sampler, Pymc3):
raise ValueError("Sampler is not a bilby Pymc3 sampler object")
return super(PyMC3UniformPrior, self).ln_prob(sampler)
return pm.Uniform(self.name, lower=self.minimum, upper=self.maximum)
......@@ -142,7 +136,7 @@ class PriorPyMC3(bilby.core.prior.Prior):
# We make a prior
priors = dict()
priors["m"] = bilby.core.prior.Uniform(0, 5, "m")
priors["c"] = PriorPyMC3(-2, 2, "c")
priors["c"] = PyMC3UniformPrior(-2, 2, "c")
# And run sampler
result = bilby.run_sampler(
......
"""
This tutorial demonstrates how we can sample a prior in the shape of a ball.
Note that this will not end up sampling uniformly in that space, constraint
priors are more suitable for that.
This implementation will draw a value for the x-coordinate from p(x), and
given that draw a value for the
y-coordinate from p(y|x), and given that draw a value for the z-coordinate
from p(z|x,y).
Only the x-coordinate will end up being uniform for this
"""
import bilby
import numpy as np
# This tutorial demonstrates how we can sample a prior in the shape of a ball
# Note that this will not end up sampling uniformly in that space, constraint priors are more suitable for that.
# This implementation will draw a value for the x-coordinate from p(x), and given that draw a value for the
# y-coordinate from p(y|x), and given that draw a value for the z-coordinate from p(z|x,y).
# Only the x-coordinate will end up being uniform for this
class ZeroLikelihood(bilby.core.likelihood.Likelihood):
"""Flat likelihood. This always returns 0.
......
......@@ -15,7 +15,7 @@ injection_parameters = dict(
)
data = [injection_parameters[label + str(ii)] * 1000 for ii in range(n_dim)]
likelihood = Multinomial(data=data, n_dimensions=n_dim, label=label)
likelihood = Multinomial(data=data, n_dimensions=n_dim, base=label)
result = run_sampler(
likelihood=likelihood,
......
......@@ -83,16 +83,13 @@ def hyper_prior(dataset, mu, sigma):
)
def run_prior(dataset):
return 1 / 20
samples = [result.posterior for result in results]
for sample in samples:
sample["prior"] = 1 / 20
evidences = [result.log_evidence for result in results]
hp_likelihood = HyperparameterLikelihood(
posteriors=samples,
hyper_prior=hyper_prior,
sampling_prior=run_prior,
log_evidences=evidences,
max_samples=500,
)
......@@ -107,6 +104,8 @@ result = run_sampler(
likelihood=hp_likelihood,
priors=hp_priors,
sampler="dynesty",
walks=10,
nact=3,
nlive=1000,
use_ratio=False,
outdir=outdir,
......
......@@ -75,12 +75,12 @@ grid_evidence = grid.log_evidence
axes = fig.get_axes()
axes[0].plot(
grid.sample_points["c"],
np.exp(grid.marginalize_ln_posterior(not_parameter="c") - grid_evidence),
np.exp(grid.marginalize_ln_posterior(not_parameters="c") - grid_evidence),
"k--",
)
axes[3].plot(
grid.sample_points["m"],
np.exp(grid.marginalize_ln_posterior(not_parameter="m") - grid_evidence),
np.exp(grid.marginalize_ln_posterior(not_parameters="m") - grid_evidence),
"k--",
)
axes[2].contour(
......
......@@ -99,7 +99,7 @@ def fit(n):
result = bilby.run_sampler(
likelihood=likelihood,
priors=priors,
nlive=1000,
nlive=100,
outdir=outdir,
label=label,
sampler="nestle",
......
......@@ -42,15 +42,9 @@ def decay_rate(delta_t, halflife, n_init):
times = np.insert(times, 0, 0.0)
n_atoms = n_init * atto * n_avogadro
counts = np.exp(-np.log(2) * times / halflife)
rates = (
(
np.exp(-np.log(2) * (times[:-1] / halflife))
- np.exp(-np.log(2) * (times[1:] / halflife))
)
* n_atoms
/ delta_t
)
rates = (counts[:-1] - counts[1:]) * n_atoms / delta_t
return rates
......@@ -95,6 +89,7 @@ result = bilby.run_sampler(
likelihood=likelihood,
priors=priors,
sampler="dynesty",
sample="unif",
nlive=1000,
injection_parameters=injection_parameters,
outdir=outdir,
......
......@@ -5,6 +5,12 @@ In this example we look at a simple example with the sum
of two Gaussian distributions, and we try to fit with
up to three Gaussians.
We will use the `PyMultiNest` sampler which is fast but can be unreliable when
significant correlations exist in the likelihood.
To install `PyMultiNest` call
$ conda install -c conda-forge pymultinest
"""
import matplotlib.pyplot as plt
......@@ -130,7 +136,7 @@ result = bilby.run_sampler(
priors=priors,
outdir=outdir,
label=label,
sampler="dynesty",
sampler="pymultinest",
nlive=400,
)
......
import matplotlib
matplotlib.use("Agg")
matplotlib.use("Agg") # noqa
import glob
import importlib.util
import unittest
import os
import parameterized
import pytest
import shutil
import logging
# Required to run the tests
from past.builtins import execfile
import bilby.core.utils
# Imported to ensure the examples run
import numpy as np # noqa: F401
import inspect # noqa: F401
bilby.core.utils.command_line_args.clean = True
core_examples = glob.glob("examples/core_examples/*.py")
core_examples += glob.glob("examples/core_examples/*/*.py")
core_args = [(fname.split("/")[-1][:-3], fname) for fname in core_examples]
bilby.core.utils.command_line_args.bilby_test_mode = True
gw_examples = [
"examples/gw_examples/injection_examples/fast_tutorial.py",
"examples/gw_examples/data_examples/GW150914.py",
]
gw_args = [(fname.split("/")[-1][:-3], fname) for fname in gw_examples]
def _execute_file(name, fname):
print("Running {}".format(fname))
spec = importlib.util.spec_from_file_location(name, fname)
foo = importlib.util.module_from_spec(spec)
spec.loader.exec_module(foo)
class ExampleTest(unittest.TestCase):
......@@ -39,25 +53,20 @@ class ExampleTest(unittest.TestCase):
except OSError:
logging.warning("{} not removed prior to tests".format(cls.outdir))
def test_examples(self):
@parameterized.parameterized.expand(core_args)
def test_core_examples(self, name, fname):
""" Loop over examples to check they run """
examples = [
"examples/core_examples/linear_regression.py",
"examples/core_examples/linear_regression_unknown_noise.py",
]
for filename in examples:
print("Testing {}".format(filename))
execfile(filename)
bilby.core.utils.command_line_args.bilby_test_mode = False
ignore = ["15d_gaussian"]
if any([item in fname for item in ignore]):
pytest.skip()
_execute_file(name, fname)
def test_gw_examples(self):
@parameterized.parameterized.expand(gw_args)
def test_gw_examples(self, name, fname):
""" Loop over examples to check they run """
examples = [
"examples/gw_examples/injection_examples/fast_tutorial.py",
"examples/gw_examples/data_examples/GW150914.py",
]
for filename in examples:
print("Testing {}".format(filename))
execfile(filename)
bilby.core.utils.command_line_args.bilby_test_mode = True
_execute_file(name, fname)
if __name__ == "__main__":
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment