From eac5426b354cfaf772112139675a08ce67993cbb Mon Sep 17 00:00:00 2001
From: Moritz Huebner <moritz.huebner@ligo.org>
Date: Mon, 1 Apr 2019 18:38:19 -0500
Subject: [PATCH] Merge remote-tracking branch 'origin/master' into
245-proposal-library
# Conflicts:
# bilby/core/sampler/cpnest.py
# examples/injection_examples/fast_tutorial.py
---
bilby/core/prior.py | 183 +++++---
bilby/core/sampler/__init__.py | 3 +-
bilby/core/sampler/cpnest.py | 84 +++-
bilby/core/sampler/proposal.py | 339 ++++++++++++++
bilby/gw/prior.py | 21 +-
bilby/gw/prior_files/GW150914.prior | 32 +-
.../aligned_spin_binary_black_holes.prior | 32 +-
bilby/gw/prior_files/binary_black_holes.prior | 44 +-
.../gw/prior_files/binary_neutron_stars.prior | 40 +-
bilby/gw/sampler/__init__.py | 2 +
bilby/gw/sampler/proposal.py | 49 ++
.../custom_proposal_example.py | 68 +++
test/gw_prior_test.py | 8 +-
test/prior_files/binary_black_holes.prior | 44 +-
test/prior_files/binary_neutron_stars.prior | 40 +-
test/prior_test.py | 137 +++---
test/proposal_test.py | 433 ++++++++++++++++++
test/sampler_test.py | 4 +-
18 files changed, 1318 insertions(+), 245 deletions(-)
create mode 100644 bilby/core/sampler/proposal.py
create mode 100644 bilby/gw/sampler/__init__.py
create mode 100644 bilby/gw/sampler/proposal.py
create mode 100644 examples/injection_examples/custom_proposal_example.py
create mode 100644 test/proposal_test.py
diff --git a/bilby/core/prior.py b/bilby/core/prior.py
index 85ed0fdff..2bac06180 100644
--- a/bilby/core/prior.py
+++ b/bilby/core/prior.py
@@ -407,7 +407,7 @@ class Prior(object):
_default_latex_labels = dict()
def __init__(self, name=None, latex_label=None, unit=None, minimum=-np.inf,
- maximum=np.inf):
+ maximum=np.inf, periodic_boundary=False):
""" Implements a Prior object
Parameters
@@ -422,13 +422,15 @@ class Prior(object):
Minimum of the domain, default=-np.inf
maximum: float, optional
Maximum of the domain, default=np.inf
-
+ periodic_boundary: bool, optional
+ Whether or not the boundary condition is periodic. Not available in all samplers.
"""
self.name = name
self.latex_label = latex_label
self.unit = unit
self.minimum = minimum
self.maximum = maximum
+ self.periodic_boundary = periodic_boundary
def __call__(self):
"""Overrides the __call__ special method. Calls the sample method.
@@ -632,6 +634,16 @@ class Prior(object):
def maximum(self, maximum):
self._maximum = maximum
+ @property
+ def periodic_boundary(self):
+ return self._periodic_boundary
+
+ @periodic_boundary.setter
+ def periodic_boundary(self, periodic_boundary):
+ if type(periodic_boundary) is not bool:
+ raise ValueError('{} is not a valid setting for prior boundaries'.format(periodic_boundary))
+ self._periodic_boundary = periodic_boundary
+
@property
def __default_latex_label(self):
if self.name in self._default_latex_labels.keys():
@@ -709,7 +721,7 @@ class DeltaFunction(Prior):
class PowerLaw(Prior):
def __init__(self, alpha, minimum, maximum, name=None, latex_label=None,
- unit=None):
+ unit=None, periodic_boundary=False):
"""Power law with bounds and alpha, spectral index
Parameters
@@ -726,9 +738,12 @@ class PowerLaw(Prior):
See superclass
unit: str
See superclass
+ periodic_boundary: bool
+ See superclass
"""
Prior.__init__(self, name=name, latex_label=latex_label,
- minimum=minimum, maximum=maximum, unit=unit)
+ minimum=minimum, maximum=maximum, unit=unit,
+ periodic_boundary=periodic_boundary)
self.alpha = alpha
def rescale(self, val):
@@ -789,13 +804,14 @@ class PowerLaw(Prior):
normalising = (1 + self.alpha) / (self.maximum ** (1 + self.alpha) -
self.minimum ** (1 + self.alpha))
- return (self.alpha * np.nan_to_num(np.log(val)) + np.log(normalising)) + np.log(1. * self.is_in_prior_range(val))
+ return (self.alpha * np.nan_to_num(np.log(val)) + np.log(normalising)) + np.log(
+ 1. * self.is_in_prior_range(val))
class Uniform(Prior):
def __init__(self, minimum, maximum, name=None, latex_label=None,
- unit=None):
+ unit=None, periodic_boundary=False):
"""Uniform prior with bounds
Parameters
@@ -810,9 +826,12 @@ class Uniform(Prior):
See superclass
unit: str
See superclass
+ periodic_boundary: bool
+ See superclass
"""
Prior.__init__(self, name=name, latex_label=latex_label,
- minimum=minimum, maximum=maximum, unit=unit)
+ minimum=minimum, maximum=maximum, unit=unit,
+ periodic_boundary=periodic_boundary)
def rescale(self, val):
Prior.test_valid_for_rescaling(val)
@@ -850,7 +869,7 @@ class Uniform(Prior):
class LogUniform(PowerLaw):
def __init__(self, minimum, maximum, name=None, latex_label=None,
- unit=None):
+ unit=None, periodic_boundary=False):
"""Log-Uniform prior with bounds
Parameters
@@ -865,9 +884,11 @@ class LogUniform(PowerLaw):
See superclass
unit: str
See superclass
+ periodic_boundary: bool
+ See superclass
"""
PowerLaw.__init__(self, name=name, latex_label=latex_label, unit=unit,
- minimum=minimum, maximum=maximum, alpha=-1)
+ minimum=minimum, maximum=maximum, alpha=-1, periodic_boundary=periodic_boundary)
if self.minimum <= 0:
logger.warning('You specified a uniform-in-log prior with minimum={}'.format(self.minimum))
@@ -875,7 +896,7 @@ class LogUniform(PowerLaw):
class SymmetricLogUniform(Prior):
def __init__(self, minimum, maximum, name=None, latex_label=None,
- unit=None):
+ unit=None, periodic_boundary=False):
"""Symmetric Log-Uniform distribtions with bounds
This is identical to a Log-Uniform distribition, but mirrored about
@@ -895,9 +916,12 @@ class SymmetricLogUniform(Prior):
See superclass
unit: str
See superclass
+ periodic_boundary: bool
+ See superclass
"""
Prior.__init__(self, name=name, latex_label=latex_label,
- minimum=minimum, maximum=maximum, unit=unit)
+ minimum=minimum, maximum=maximum, unit=unit,
+ periodic_boundary=periodic_boundary)
def rescale(self, val):
"""
@@ -955,7 +979,7 @@ class SymmetricLogUniform(Prior):
class Cosine(Prior):
def __init__(self, name=None, latex_label=None, unit=None,
- minimum=-np.pi / 2, maximum=np.pi / 2):
+ minimum=-np.pi / 2, maximum=np.pi / 2, periodic_boundary=False):
"""Cosine prior with bounds
Parameters
@@ -970,9 +994,11 @@ class Cosine(Prior):
See superclass
unit: str
See superclass
+ periodic_boundary: bool
+ See superclass
"""
Prior.__init__(self, name=name, latex_label=latex_label, unit=unit,
- minimum=minimum, maximum=maximum)
+ minimum=minimum, maximum=maximum, periodic_boundary=periodic_boundary)
def rescale(self, val):
"""
@@ -1001,7 +1027,7 @@ class Cosine(Prior):
class Sine(Prior):
def __init__(self, name=None, latex_label=None, unit=None, minimum=0,
- maximum=np.pi):
+ maximum=np.pi, periodic_boundary=False):
"""Sine prior with bounds
Parameters
@@ -1016,9 +1042,11 @@ class Sine(Prior):
See superclass
unit: str
See superclass
+ periodic_boundary: bool
+ See superclass
"""
Prior.__init__(self, name=name, latex_label=latex_label, unit=unit,
- minimum=minimum, maximum=maximum)
+ minimum=minimum, maximum=maximum, periodic_boundary=periodic_boundary)
def rescale(self, val):
"""
@@ -1046,7 +1074,7 @@ class Sine(Prior):
class Gaussian(Prior):
- def __init__(self, mu, sigma, name=None, latex_label=None, unit=None):
+ def __init__(self, mu, sigma, name=None, latex_label=None, unit=None, periodic_boundary=False):
"""Gaussian prior with mean mu and width sigma
Parameters
@@ -1061,8 +1089,10 @@ class Gaussian(Prior):
See superclass
unit: str
See superclass
+ periodic_boundary: bool
+ See superclass
"""
- Prior.__init__(self, name=name, latex_label=latex_label, unit=unit)
+ Prior.__init__(self, name=name, latex_label=latex_label, unit=unit, periodic_boundary=periodic_boundary)
self.mu = mu
self.sigma = sigma
@@ -1094,7 +1124,7 @@ class Gaussian(Prior):
class Normal(Gaussian):
- def __init__(self, mu, sigma, name=None, latex_label=None, unit=None):
+ def __init__(self, mu, sigma, name=None, latex_label=None, unit=None, periodic_boundary=False):
"""A synonym for the Gaussian distribution.
Parameters
@@ -1109,15 +1139,17 @@ class Normal(Gaussian):
See superclass
unit: str
See superclass
+ periodic_boundary: bool
+ See superclass
"""
- Gaussian.__init__(self, mu=mu, sigma=sigma, name=name,
- latex_label=latex_label, unit=unit)
+ Gaussian.__init__(self, mu=mu, sigma=sigma, name=name, latex_label=latex_label,
+ unit=unit, periodic_boundary=periodic_boundary)
class TruncatedGaussian(Prior):
def __init__(self, mu, sigma, minimum, maximum, name=None,
- latex_label=None, unit=None):
+ latex_label=None, unit=None, periodic_boundary=False):
"""Truncated Gaussian prior with mean mu and width sigma
https://en.wikipedia.org/wiki/Truncated_normal_distribution
@@ -1138,9 +1170,11 @@ class TruncatedGaussian(Prior):
See superclass
unit: str
See superclass
+ periodic_boundary: bool
+ See superclass
"""
Prior.__init__(self, name=name, latex_label=latex_label, unit=unit,
- minimum=minimum, maximum=maximum)
+ minimum=minimum, maximum=maximum, periodic_boundary=periodic_boundary)
self.mu = mu
self.sigma = sigma
@@ -1176,14 +1210,14 @@ class TruncatedGaussian(Prior):
-------
float: Prior probability of val
"""
- return np.exp(-(self.mu - val) ** 2 / (2 * self.sigma ** 2)) / (
- 2 * np.pi) ** 0.5 / self.sigma / self.normalisation * self.is_in_prior_range(val)
+ return np.exp(-(self.mu - val) ** 2 / (2 * self.sigma ** 2)) / \
+ (2 * np.pi) ** 0.5 / self.sigma / self.normalisation * self.is_in_prior_range(val)
class TruncatedNormal(TruncatedGaussian):
def __init__(self, mu, sigma, minimum, maximum, name=None,
- latex_label=None, unit=None):
+ latex_label=None, unit=None, periodic_boundary=False):
"""A synonym for the TruncatedGaussian distribution.
Parameters
@@ -1202,14 +1236,16 @@ class TruncatedNormal(TruncatedGaussian):
See superclass
unit: str
See superclass
+ periodic_boundary: bool
+ See superclass
"""
TruncatedGaussian.__init__(self, mu=mu, sigma=sigma, minimum=minimum,
- maximum=maximum, name=name,
- latex_label=latex_label, unit=unit)
+ maximum=maximum, name=name, latex_label=latex_label,
+ unit=unit, periodic_boundary=periodic_boundary)
class HalfGaussian(TruncatedGaussian):
- def __init__(self, sigma, name=None, latex_label=None, unit=None):
+ def __init__(self, sigma, name=None, latex_label=None, unit=None, periodic_boundary=False):
"""A Gaussian with its mode at zero, and truncated to only be positive.
Parameters
@@ -1222,14 +1258,16 @@ class HalfGaussian(TruncatedGaussian):
See superclass
unit: str
See superclass
+ periodic_boundary: bool
+ See superclass
"""
TruncatedGaussian.__init__(self, 0., sigma, minimum=0., maximum=np.inf,
name=name, latex_label=latex_label,
- unit=unit)
+ unit=unit, periodic_boundary=periodic_boundary)
class HalfNormal(HalfGaussian):
- def __init__(self, sigma, name=None, latex_label=None, unit=None):
+ def __init__(self, sigma, name=None, latex_label=None, unit=None, periodic_boundary=False):
"""A synonym for the HalfGaussian distribution.
Parameters
@@ -1242,14 +1280,16 @@ class HalfNormal(HalfGaussian):
See superclass
unit: str
See superclass
-
+ periodic_boundary: bool
+ See superclass
"""
HalfGaussian.__init__(self, sigma=sigma, name=name,
- latex_label=latex_label, unit=unit)
+ latex_label=latex_label, unit=unit,
+ periodic_boundary=periodic_boundary)
class LogNormal(Prior):
- def __init__(self, mu, sigma, name=None, latex_label=None, unit=None):
+ def __init__(self, mu, sigma, name=None, latex_label=None, unit=None, periodic_boundary=False):
"""Log-normal prior with mean mu and width sigma
https://en.wikipedia.org/wiki/Log-normal_distribution
@@ -1266,10 +1306,11 @@ class LogNormal(Prior):
See superclass
unit: str
See superclass
-
+ periodic_boundary: bool
+ See superclass
"""
Prior.__init__(self, name=name, minimum=0., latex_label=latex_label,
- unit=unit)
+ unit=unit, periodic_boundary=periodic_boundary)
if sigma <= 0.:
raise ValueError("For the LogGaussian prior the standard deviation must be positive")
@@ -1305,7 +1346,7 @@ class LogNormal(Prior):
class LogGaussian(LogNormal):
- def __init__(self, mu, sigma, name=None, latex_label=None, unit=None):
+ def __init__(self, mu, sigma, name=None, latex_label=None, unit=None, periodic_boundary=False):
"""Synonym of LogNormal prior
https://en.wikipedia.org/wiki/Log-normal_distribution
@@ -1322,14 +1363,15 @@ class LogGaussian(LogNormal):
See superclass
unit: str
See superclass
-
+ periodic_boundary: bool
+ See superclass
"""
LogNormal.__init__(self, mu=mu, sigma=sigma, name=name,
- latex_label=latex_label, unit=unit)
+ latex_label=latex_label, unit=unit, periodic_boundary=periodic_boundary)
class Exponential(Prior):
- def __init__(self, mu, name=None, latex_label=None, unit=None):
+ def __init__(self, mu, name=None, latex_label=None, unit=None, periodic_boundary=False):
"""Exponential prior with mean mu
Parameters
@@ -1342,10 +1384,11 @@ class Exponential(Prior):
See superclass
unit: str
See superclass
-
+ periodic_boundary: bool
+ See superclass
"""
Prior.__init__(self, name=name, minimum=0., latex_label=latex_label,
- unit=unit)
+ unit=unit, periodic_boundary=periodic_boundary)
self.mu = mu
def rescale(self, val):
@@ -1377,7 +1420,7 @@ class Exponential(Prior):
class StudentT(Prior):
def __init__(self, df, mu=0., scale=1., name=None, latex_label=None,
- unit=None):
+ unit=None, periodic_boundary=False):
"""Student's t-distribution prior with number of degrees of freedom df,
mean mu and scale
@@ -1397,8 +1440,10 @@ class StudentT(Prior):
See superclass
unit: str
See superclass
+ periodic_boundary: bool
+ See superclass
"""
- Prior.__init__(self, name=name, latex_label=latex_label, unit=unit)
+ Prior.__init__(self, name=name, latex_label=latex_label, unit=unit, periodic_boundary=periodic_boundary)
if df <= 0. or scale <= 0.:
raise ValueError("For the StudentT prior the number of degrees of freedom and scale must be positive")
@@ -1437,7 +1482,7 @@ class StudentT(Prior):
class Beta(Prior):
def __init__(self, alpha, beta, minimum=0, maximum=1, name=None,
- latex_label=None, unit=None):
+ latex_label=None, unit=None, periodic_boundary=False):
"""Beta distribution
https://en.wikipedia.org/wiki/Beta_distribution
@@ -1461,7 +1506,8 @@ class Beta(Prior):
See superclass
unit: str
See superclass
-
+ periodic_boundary: bool
+ See superclass
"""
if alpha <= 0. or beta <= 0.:
raise ValueError("alpha and beta must both be positive values")
@@ -1471,7 +1517,7 @@ class Beta(Prior):
self._minimum = minimum
self._maximum = maximum
Prior.__init__(self, minimum=minimum, maximum=maximum, name=name,
- latex_label=latex_label, unit=unit)
+ latex_label=latex_label, unit=unit, periodic_boundary=periodic_boundary)
self._set_dist()
def rescale(self, val):
@@ -1565,7 +1611,7 @@ class Beta(Prior):
class Logistic(Prior):
- def __init__(self, mu, scale, name=None, latex_label=None, unit=None):
+ def __init__(self, mu, scale, name=None, latex_label=None, unit=None, periodic_boundary=False):
"""Logistic distribution
https://en.wikipedia.org/wiki/Logistic_distribution
@@ -1582,8 +1628,10 @@ class Logistic(Prior):
See superclass
unit: str
See superclass
+ periodic_boundary: bool
+ See superclass
"""
- Prior.__init__(self, name=name, latex_label=latex_label, unit=unit)
+ Prior.__init__(self, name=name, latex_label=latex_label, unit=unit, periodic_boundary=periodic_boundary)
if scale <= 0.:
raise ValueError("For the Logistic prior the scale must be positive")
@@ -1620,7 +1668,7 @@ class Logistic(Prior):
class Cauchy(Prior):
- def __init__(self, alpha, beta, name=None, latex_label=None, unit=None):
+ def __init__(self, alpha, beta, name=None, latex_label=None, unit=None, periodic_boundary=False):
"""Cauchy distribution
https://en.wikipedia.org/wiki/Cauchy_distribution
@@ -1637,8 +1685,10 @@ class Cauchy(Prior):
See superclass
unit: str
See superclass
+ periodic_boundary: bool
+ See superclass
"""
- Prior.__init__(self, name=name, latex_label=latex_label, unit=unit)
+ Prior.__init__(self, name=name, latex_label=latex_label, unit=unit, periodic_boundary=periodic_boundary)
if beta <= 0.:
raise ValueError("For the Cauchy prior the scale must be positive")
@@ -1675,7 +1725,7 @@ class Cauchy(Prior):
class Lorentzian(Cauchy):
- def __init__(self, alpha, beta, name=None, latex_label=None, unit=None):
+ def __init__(self, alpha, beta, name=None, latex_label=None, unit=None, periodic_boundary=False):
"""Synonym for the Cauchy distribution
https://en.wikipedia.org/wiki/Cauchy_distribution
@@ -1692,13 +1742,15 @@ class Lorentzian(Cauchy):
See superclass
unit: str
See superclass
+ periodic_boundary: bool
+ See superclass
"""
Cauchy.__init__(self, alpha=alpha, beta=beta, name=name,
- latex_label=latex_label, unit=unit)
+ latex_label=latex_label, unit=unit, periodic_boundary=periodic_boundary)
class Gamma(Prior):
- def __init__(self, k, theta=1., name=None, latex_label=None, unit=None):
+ def __init__(self, k, theta=1., name=None, latex_label=None, unit=None, periodic_boundary=False):
"""Gamma distribution
https://en.wikipedia.org/wiki/Gamma_distribution
@@ -1715,9 +1767,11 @@ class Gamma(Prior):
See superclass
unit: str
See superclass
+ periodic_boundary: bool
+ See superclass
"""
Prior.__init__(self, name=name, minimum=0., latex_label=latex_label,
- unit=unit)
+ unit=unit, periodic_boundary=periodic_boundary)
if k <= 0 or theta <= 0:
raise ValueError("For the Gamma prior the shape and scale must be positive")
@@ -1755,7 +1809,7 @@ class Gamma(Prior):
class ChiSquared(Gamma):
- def __init__(self, nu, name=None, latex_label=None, unit=None):
+ def __init__(self, nu, name=None, latex_label=None, unit=None, periodic_boundary=False):
"""Chi-squared distribution
https://en.wikipedia.org/wiki/Chi-squared_distribution
@@ -1770,13 +1824,15 @@ class ChiSquared(Gamma):
See superclass
unit: str
See superclass
+ periodic_boundary: bool
+ See superclass
"""
if nu <= 0 or not isinstance(nu, int):
raise ValueError("For the ChiSquared prior the number of degrees of freedom must be a positive integer")
Gamma.__init__(self, name=name, k=nu / 2., theta=2.,
- latex_label=latex_label, unit=unit)
+ latex_label=latex_label, unit=unit, periodic_boundary=periodic_boundary)
@property
def nu(self):
@@ -1790,7 +1846,7 @@ class ChiSquared(Gamma):
class Interped(Prior):
def __init__(self, xx, yy, minimum=np.nan, maximum=np.nan, name=None,
- latex_label=None, unit=None):
+ latex_label=None, unit=None, periodic_boundary=False):
"""Creates an interpolated prior function from arrays of xx and yy=p(xx)
Parameters
@@ -1809,6 +1865,8 @@ class Interped(Prior):
See superclass
unit: str
See superclass
+ periodic_boundary: bool
+ See superclass
Attributes
-------
@@ -1827,7 +1885,8 @@ class Interped(Prior):
self.__all_interpolated = interp1d(x=xx, y=yy, bounds_error=False, fill_value=0)
Prior.__init__(self, name=name, latex_label=latex_label, unit=unit,
minimum=np.nanmax(np.array((min(xx), minimum))),
- maximum=np.nanmin(np.array((max(xx), maximum))))
+ maximum=np.nanmin(np.array((max(xx), maximum))),
+ periodic_boundary=periodic_boundary)
self._update_instance()
def __eq__(self, other):
@@ -1920,7 +1979,7 @@ class Interped(Prior):
class FromFile(Interped):
def __init__(self, file_name, minimum=None, maximum=None, name=None,
- latex_label=None, unit=None):
+ latex_label=None, unit=None, periodic_boundary=False):
"""Creates an interpolated prior function from arrays of xx and yy=p(xx) extracted from a file
Parameters
@@ -1937,6 +1996,8 @@ class FromFile(Interped):
See superclass
unit: str
See superclass
+ periodic_boundary: bool
+ See superclass
Attributes
-------
@@ -1948,8 +2009,8 @@ class FromFile(Interped):
self.id = file_name
xx, yy = np.genfromtxt(self.id).T
Interped.__init__(self, xx=xx, yy=yy, minimum=minimum,
- maximum=maximum, name=name,
- latex_label=latex_label, unit=unit)
+ maximum=maximum, name=name, latex_label=latex_label,
+ unit=unit, periodic_boundary=periodic_boundary)
except IOError:
logger.warning("Can't load {}.".format(self.id))
logger.warning("Format should be:")
diff --git a/bilby/core/sampler/__init__.py b/bilby/core/sampler/__init__.py
index 028c5f5aa..7feed4e80 100644
--- a/bilby/core/sampler/__init__.py
+++ b/bilby/core/sampler/__init__.py
@@ -17,6 +17,7 @@ from .ptmcmc import PTMCMCSampler
from .pymc3 import Pymc3
from .pymultinest import Pymultinest
from .fake_sampler import FakeSampler
+from . import proposal
IMPLEMENTED_SAMPLERS = {
'cpnest': Cpnest, 'dynesty': Dynesty, 'emcee': Emcee, 'nestle': Nestle,
@@ -89,7 +90,7 @@ def run_sampler(likelihood, priors=None, label='label', outdir='outdir',
If true, save the priors and results to disk.
If hdf5, save as an hdf5 file instead of json.
gzip: bool
- If true, and save is true, gzip the saved results file.
+ If true, and save is true, gzip the saved results file.
result_class: bilby.core.result.Result, or child of
The result class to use. By default, `bilby.core.result.Result` is used,
but objects which inherit from this class can be given providing
diff --git a/bilby/core/sampler/cpnest.py b/bilby/core/sampler/cpnest.py
index eec05efc3..6acd33689 100644
--- a/bilby/core/sampler/cpnest.py
+++ b/bilby/core/sampler/cpnest.py
@@ -1,9 +1,12 @@
from __future__ import absolute_import
+import copy
+
import numpy as np
from pandas import DataFrame
from .base_sampler import NestedSampler
+from .proposal import Sample, JumpProposalCycle
from ..utils import logger, check_directory_exists_and_if_not_mkdir
@@ -37,7 +40,7 @@ class Cpnest(NestedSampler):
"""
default_kwargs = dict(verbose=1, nthreads=1, nlive=500, maxmcmc=1000,
seed=None, poolsize=100, nhamiltonian=0, resume=True,
- output=None)
+ output=None, proposals=None)
def _translate_kwargs(self, kwargs):
if 'nlive' not in kwargs:
@@ -78,8 +81,19 @@ class Cpnest(NestedSampler):
for name in self.names])
return point
+ self._resolve_proposal_functions()
model = Model(self.search_parameter_keys, self.priors)
- out = CPNest(model, **self.kwargs)
+ try:
+ out = CPNest(model, **self.kwargs)
+ except TypeError as e:
+ if 'proposals' in self.kwargs.keys():
+ logger.warning('YOU ARE TRYING TO USE PROPOSALS IN A VERSION OF CPNEST THAT DOES'
+ 'NOT ACCEPT CUSTOM PROPOSALS. SAMPLING WILL COMMENCE WITH THE DEFAULT'
+ 'PROPOSALS.')
+ del self.kwargs['proposals']
+ out = CPNest(model, **self.kwargs)
+ else:
+ raise TypeError(e)
out.run()
if self.plot:
@@ -103,3 +117,69 @@ class Cpnest(NestedSampler):
self.kwargs['output'] = '{}/'.format(self.kwargs['output'])
check_directory_exists_and_if_not_mkdir(self.kwargs['output'])
NestedSampler._verify_kwargs_against_default_kwargs(self)
+
+ def _resolve_proposal_functions(self):
+ from cpnest.proposal import ProposalCycle
+ if 'proposals' in self.kwargs:
+ if self.kwargs['proposals'] is None:
+ return
+ if type(self.kwargs['proposals']) == JumpProposalCycle:
+ self.kwargs['proposals'] = dict(mhs=self.kwargs['proposals'], hmc=self.kwargs['proposals'])
+ for key, proposal in self.kwargs['proposals'].items():
+ if isinstance(proposal, JumpProposalCycle):
+ self.kwargs['proposals'][key] = cpnest_proposal_cycle_factory(proposal)
+ elif isinstance(proposal, ProposalCycle):
+ pass
+ else:
+ raise TypeError("Unknown proposal type")
+
+
+def cpnest_proposal_factory(jump_proposal):
+ import cpnest.proposal
+
+ class CPNestEnsembleProposal(cpnest.proposal.EnsembleProposal):
+
+ def __init__(self, jp):
+ self.jump_proposal = jp
+ self.ensemble = None
+
+ def __call__(self, sample, **kwargs):
+ return self.get_sample(sample, **kwargs)
+
+ def get_sample(self, cpnest_sample, **kwargs):
+ sample = Sample.from_cpnest_live_point(cpnest_sample)
+ self.ensemble = kwargs.get('coordinates', self.ensemble)
+ sample = self.jump_proposal(sample=sample, sampler_name='cpnest', **kwargs)
+ self.log_J = self.jump_proposal.log_j
+ return self._update_cpnest_sample(cpnest_sample, sample)
+
+ @staticmethod
+ def _update_cpnest_sample(cpnest_sample, sample):
+ cpnest_sample.names = list(sample.keys())
+ for i, value in enumerate(sample.values()):
+ cpnest_sample.values[i] = value
+ return cpnest_sample
+
+ return CPNestEnsembleProposal(jump_proposal)
+
+
+def cpnest_proposal_cycle_factory(jump_proposals):
+ import cpnest.proposal
+
+ class CPNestProposalCycle(cpnest.proposal.ProposalCycle):
+ def __init__(self):
+ self.jump_proposals = copy.deepcopy(jump_proposals)
+ for i, prop in enumerate(self.jump_proposals.proposal_functions):
+ self.jump_proposals.proposal_functions[i] = cpnest_proposal_factory(prop)
+ self.jump_proposals.update_cycle()
+ super(CPNestProposalCycle, self).__init__(proposals=self.jump_proposals.proposal_functions,
+ weights=self.jump_proposals.weights,
+ cyclelength=self.jump_proposals.cycle_length)
+
+ def get_sample(self, old, **kwargs):
+ return self.jump_proposals(sample=old, coordinates=self.ensemble, **kwargs)
+
+ def set_ensemble(self, ensemble):
+ self.ensemble = ensemble
+
+ return CPNestProposalCycle
diff --git a/bilby/core/sampler/proposal.py b/bilby/core/sampler/proposal.py
new file mode 100644
index 000000000..a928292b3
--- /dev/null
+++ b/bilby/core/sampler/proposal.py
@@ -0,0 +1,339 @@
+from collections import OrderedDict
+from inspect import isclass
+
+import numpy as np
+import random
+
+from bilby.core.prior import Uniform
+
+
+class Sample(OrderedDict):
+
+ def __init__(self, dictionary=None):
+ if dictionary is None:
+ dictionary = dict()
+ super(Sample, self).__init__(dictionary)
+
+ def __add__(self, other):
+ return Sample({key: self[key] + other[key] for key in self.keys()})
+
+ def __sub__(self, other):
+ return Sample({key: self[key] - other[key] for key in self.keys()})
+
+ def __mul__(self, other):
+ return Sample({key: self[key] * other for key in self.keys()})
+
+ @classmethod
+ def from_cpnest_live_point(cls, cpnest_live_point):
+ res = cls(dict())
+ for i, key in enumerate(cpnest_live_point.names):
+ res[key] = cpnest_live_point.values[i]
+ return res
+
+ @classmethod
+ def from_external_type(cls, external_sample, sampler_name):
+ if sampler_name == 'cpnest':
+ return cls.from_cpnest_live_point(external_sample)
+ return external_sample
+
+
+class JumpProposal(object):
+
+ def __init__(self, priors=None):
+ """ A generic class for jump proposals
+
+ Parameters
+ ----------
+ priors: bilby.core.prior.PriorDict
+ Dictionary of priors used in this sampling run
+
+ Attributes
+ ----------
+ log_j: float
+ Log Jacobian of the proposal. Characterises whether or not detailed balance
+ is preserved. If not, log_j needs to be adjusted accordingly.
+ """
+ self.priors = priors
+ self.log_j = 0.0
+
+ def __call__(self, sample, **kwargs):
+ """ A generic wrapper for the jump proposal function
+
+ Parameters
+ ----------
+ args: Arguments that are going to be passed into the proposal function
+ kwargs: Keyword arguments that are going to be passed into the proposal function
+
+ Returns
+ -------
+ dict: A dictionary with the new samples. Boundary conditions are applied.
+
+ """
+ return self._apply_boundaries(sample)
+
+ def _move_reflecting_keys(self, sample):
+ keys = [key for key in sample.keys() if not self.priors[key].periodic_boundary]
+ for key in keys:
+ if sample[key] > self.priors[key].maximum or sample[key] < self.priors[key].minimum:
+ r = self.priors[key].maximum - self.priors[key].minimum
+ delta = (sample[key] - self.priors[key].minimum) % (2 * r)
+ if delta > r:
+ sample[key] = 2 * self.priors[key].maximum - self.priors[key].minimum - delta
+ elif delta < r:
+ sample[key] = self.priors[key].minimum + delta
+ return sample
+
+ def _move_periodic_keys(self, sample):
+ keys = [key for key in sample.keys() if self.priors[key].periodic_boundary]
+ for key in keys:
+ if sample[key] > self.priors[key].maximum or sample[key] < self.priors[key].minimum:
+ sample[key] = (self.priors[key].minimum +
+ ((sample[key] - self.priors[key].minimum) %
+ (self.priors[key].maximum - self.priors[key].minimum)))
+ return sample
+
+ def _apply_boundaries(self, sample):
+ sample = self._move_periodic_keys(sample)
+ sample = self._move_reflecting_keys(sample)
+ return sample
+
+
+class JumpProposalCycle(object):
+
+ def __init__(self, proposal_functions, weights, cycle_length=100):
+ """ A generic wrapper class for proposal cycles
+
+ Parameters
+ ----------
+ proposal_functions: list
+ A list of callable proposal functions/objects
+ weights: list
+ A list of integer weights for the respective proposal functions
+ cycle_length: int, optional
+ Length of the proposal cycle
+ """
+ self.proposal_functions = proposal_functions
+ self.weights = weights
+ self.cycle_length = cycle_length
+ self._index = 0
+ self._cycle = np.array([])
+ self.update_cycle()
+
+ def __call__(self, **kwargs):
+ proposal = self._cycle[self.index]
+ self._index = (self.index + 1) % self.cycle_length
+ return proposal(**kwargs)
+
+ def __len__(self):
+ return len(self.proposal_functions)
+
+ def update_cycle(self):
+ self._cycle = np.random.choice(self.proposal_functions, size=self.cycle_length,
+ p=self.weights, replace=True)
+
+ @property
+ def proposal_functions(self):
+ return self._proposal_functions
+
+ @proposal_functions.setter
+ def proposal_functions(self, proposal_functions):
+ for i, proposal in enumerate(proposal_functions):
+ if isclass(proposal):
+ proposal_functions[i] = proposal()
+ self._proposal_functions = proposal_functions
+
+ @property
+ def index(self):
+ return self._index
+
+ @property
+ def weights(self):
+ """
+
+ Returns
+ -------
+ Normalised proposal weights
+
+ """
+ return np.array(self._weights) / np.sum(np.array(self._weights))
+
+ @weights.setter
+ def weights(self, weights):
+ assert len(weights) == len(self.proposal_functions)
+ self._weights = weights
+
+ @property
+ def unnormalised_weights(self):
+ return self._weights
+
+
+class NormJump(JumpProposal):
+ def __init__(self, step_size, priors=None):
+ """
+ A normal distributed step centered around the old sample
+
+ Parameters
+ ----------
+ step_size: float
+ The scalable step size
+ priors:
+ See superclass
+ """
+ super(NormJump, self).__init__(priors)
+ self.step_size = step_size
+
+ def __call__(self, sample, **kwargs):
+ for key in sample.keys():
+ sample[key] = np.random.normal(sample[key], self.step_size)
+ return super(NormJump, self).__call__(sample)
+
+
+class EnsembleWalk(JumpProposal):
+
+ def __init__(self, random_number_generator=random.random, n_points=3, priors=None,
+ **random_number_generator_args):
+ """
+ An ensemble walk
+ Parameters
+ ----------
+ random_number_generator: func, optional
+ A random number generator. Default is random.random
+ n_points: int, optional
+ Number of points in the ensemble to average over. Default is 3.
+ priors:
+ See superclass
+ random_number_generator_args:
+ Additional keyword arguments for the random number generator
+ """
+ super(EnsembleWalk, self).__init__(priors)
+ self.random_number_generator = random_number_generator
+ self.n_points = n_points
+ self.random_number_generator_args = random_number_generator_args
+
+ def __call__(self, sample, **kwargs):
+ subset = random.sample(kwargs['coordinates'], self.n_points)
+ for i in range(len(subset)):
+ subset[i] = Sample.from_external_type(subset[i], kwargs.get('sampler_name', None))
+ center_of_mass = self.get_center_of_mass(subset)
+ for x in subset:
+ sample += (x - center_of_mass) * self.random_number_generator(**self.random_number_generator_args)
+ return super(EnsembleWalk, self).__call__(sample)
+
+ @staticmethod
+ def get_center_of_mass(subset):
+ return {key: np.mean([c[key] for c in subset]) for key in subset[0].keys()}
+
+
+class EnsembleStretch(JumpProposal):
+
+ def __init__(self, scale=2.0, priors=None):
+ """
+ Stretch move. Calculates the log Jacobian which can be used in cpnest to bias future moves.
+
+ Parameters
+ ----------
+ scale: float, optional
+ Stretching scale. Default is 2.0.
+ """
+ super(EnsembleStretch, self).__init__(priors)
+ self.scale = scale
+
+ def __call__(self, sample, **kwargs):
+ second_sample = random.choice(kwargs['coordinates'])
+ second_sample = Sample.from_external_type(second_sample, kwargs.get('sampler_name', None))
+ step = random.uniform(-1, 1) * np.log(self.scale)
+ sample = second_sample + (sample - second_sample) * np.exp(step)
+ self.log_j = len(sample) * step
+ return super(EnsembleStretch, self).__call__(sample)
+
+
+class DifferentialEvolution(JumpProposal):
+
+ def __init__(self, sigma=1e-4, mu=1.0, priors=None):
+ """
+ Differential evolution step. Takes two elements from the existing coordinates and differentially evolves the
+ old sample based on them using some Gaussian randomisation in the step.
+
+ Parameters
+ ----------
+ sigma: float, optional
+ Random spread in the evolution step. Default is 1e-4
+ mu: float, optional
+ Scale of the randomization. Default is 1.0
+ """
+ super(DifferentialEvolution, self).__init__(priors)
+ self.sigma = sigma
+ self.mu = mu
+
+ def __call__(self, sample, **kwargs):
+ a, b = random.sample(kwargs['coordinates'], 2)
+ sample = sample + (b - a) * random.gauss(self.mu, self.sigma)
+ return super(DifferentialEvolution, self).__call__(sample)
+
+
+class EnsembleEigenVector(JumpProposal):
+
+ def __init__(self, priors=None):
+ """
+ Ensemble step based on the ensemble eigenvectors.
+
+ Parameters
+ ----------
+ priors:
+ See superclass
+ """
+ super(EnsembleEigenVector, self).__init__(priors)
+ self.eigen_values = None
+ self.eigen_vectors = None
+ self.covariance = None
+
+ def update_eigenvectors(self, coordinates):
+ if coordinates is None:
+ return
+ elif len(coordinates[0]) == 1:
+ self._set_1_d_eigenvectors(coordinates)
+ else:
+ self._set_n_d_eigenvectors(coordinates)
+
+ def _set_1_d_eigenvectors(self, coordinates):
+ n_samples = len(coordinates)
+ key = list(coordinates[0].keys())[0]
+ variance = np.var([coordinates[j][key] for j in range(n_samples)])
+ self.eigen_values = np.atleast_1d(variance)
+ self.covariance = self.eigen_values
+ self.eigen_vectors = np.eye(1)
+
+ def _set_n_d_eigenvectors(self, coordinates):
+ n_samples = len(coordinates)
+ dim = len(coordinates[0])
+ cov_array = np.zeros((dim, n_samples))
+ for i, key in enumerate(coordinates[0].keys()):
+ for j in range(n_samples):
+ cov_array[i, j] = coordinates[j][key]
+ self.covariance = np.cov(cov_array)
+ self.eigen_values, self.eigen_vectors = np.linalg.eigh(self.covariance)
+
+ def __call__(self, sample, **kwargs):
+ self.update_eigenvectors(kwargs['coordinates'])
+ i = random.randrange(len(sample))
+ jump_size = np.sqrt(np.fabs(self.eigen_values[i])) * random.gauss(0, 1)
+ for j, key in enumerate(sample.keys()):
+ sample[key] += jump_size * self.eigen_vectors[j, i]
+ return super(EnsembleEigenVector, self).__call__(sample)
+
+
+class DrawFlatPrior(JumpProposal):
+ """
+ Draws a proposal from the flattened prior distribution.
+ """
+
+ def __call__(self, sample, **kwargs):
+ sample = _draw_from_flat_priors(sample, self.priors)
+ return super(DrawFlatPrior, self).__call__(sample)
+
+
+def _draw_from_flat_priors(sample, priors):
+ for key in sample.keys():
+ flat_prior = Uniform(priors[key].minimum, priors[key].maximum, priors[key].name)
+ sample[key] = flat_prior.sample()
+ return sample
diff --git a/bilby/gw/prior.py b/bilby/gw/prior.py
index 58790e41c..c751303d2 100644
--- a/bilby/gw/prior.py
+++ b/bilby/gw/prior.py
@@ -31,7 +31,7 @@ class Cosmological(Interped):
name='comoving_distance', latex_label='$d_C$', unit=units.Mpc))
def __init__(self, minimum, maximum, cosmology=None, name=None,
- latex_label=None, unit=None):
+ latex_label=None, unit=None, periodic_boundary=False):
self.cosmology = get_cosmology(cosmology)
if name not in self._default_args_dict:
raise ValueError(
@@ -59,7 +59,7 @@ class Cosmological(Interped):
else:
raise ValueError('Name {} not recognized.'.format(name))
Interped.__init__(self, xx=xx, yy=yy, minimum=minimum, maximum=maximum,
- **label_args)
+ periodic_boundary=periodic_boundary, **label_args)
@property
def minimum(self):
@@ -162,7 +162,7 @@ class UniformComovingVolume(Cosmological):
class AlignedSpin(Interped):
def __init__(self, a_prior=Uniform(0, 1), z_prior=Uniform(-1, 1),
- name=None, latex_label=None, unit=None):
+ name=None, latex_label=None, unit=None, periodic_boundary=False):
"""
Prior distribution for the aligned (z) component of the spin.
@@ -193,7 +193,8 @@ class AlignedSpin(Interped):
yy = [np.trapz(np.nan_to_num(a_prior.prob(aas) / aas *
z_prior.prob(x / aas)), aas) for x in xx]
Interped.__init__(self, xx=xx, yy=yy, name=name,
- latex_label=latex_label, unit=unit)
+ latex_label=latex_label, unit=unit,
+ periodic_boundary=periodic_boundary)
class BBHPriorDict(PriorDict):
@@ -509,13 +510,15 @@ class CalibrationPriorDict(PriorDict):
latex_label = "$A^{}_{}$".format(label, ii)
prior[name] = Gaussian(mu=amplitude_mean_nodes[ii],
sigma=amplitude_sigma_nodes[ii],
- name=name, latex_label=latex_label)
+ name=name, latex_label=latex_label,
+ periodic_boundary=False)
for ii in range(n_nodes):
name = "recalib_{}_phase_{}".format(label, ii)
latex_label = "$\\phi^{}_{}$".format(label, ii)
prior[name] = Gaussian(mu=phase_mean_nodes[ii],
sigma=phase_sigma_nodes[ii],
- name=name, latex_label=latex_label)
+ name=name, latex_label=latex_label,
+ periodic_boundary=False)
for ii in range(n_nodes):
name = "recalib_{}_frequency_{}".format(label, ii)
latex_label = "$f^{}_{}$".format(label, ii)
@@ -568,13 +571,15 @@ class CalibrationPriorDict(PriorDict):
latex_label = "$A^{}_{}$".format(label, ii)
prior[name] = Gaussian(mu=amplitude_mean_nodes[ii],
sigma=amplitude_sigma_nodes[ii],
- name=name, latex_label=latex_label)
+ name=name, latex_label=latex_label,
+ periodic_boundary=False)
for ii in range(n_nodes):
name = "recalib_{}_phase_{}".format(label, ii)
latex_label = "$\\phi^{}_{}$".format(label, ii)
prior[name] = Gaussian(mu=phase_mean_nodes[ii],
sigma=phase_sigma_nodes[ii],
- name=name, latex_label=latex_label)
+ name=name, latex_label=latex_label,
+ periodic_boundary=False)
for ii in range(n_nodes):
name = "recalib_{}_frequency_{}".format(label, ii)
latex_label = "$f^{}_{}$".format(label, ii)
diff --git a/bilby/gw/prior_files/GW150914.prior b/bilby/gw/prior_files/GW150914.prior
index 410fc9276..9641b77d9 100644
--- a/bilby/gw/prior_files/GW150914.prior
+++ b/bilby/gw/prior_files/GW150914.prior
@@ -1,20 +1,20 @@
# These are the default priors for analysing GW150914.
-mass_1 = Uniform(name='mass_1', minimum=30, maximum=50, unit='$M_{\\odot}$')
-mass_2 = Uniform(name='mass_2', minimum=20, maximum=40, unit='$M_{\\odot}$')
-mass_ratio = Constraint(name='mass_ratio', minimum=0.125, maximum=1)
-a_1 = Uniform(name='a_1', minimum=0, maximum=0.8)
-a_2 = Uniform(name='a_2', minimum=0, maximum=0.8)
-tilt_1 = Sine(name='tilt_1')
-tilt_2 = Sine(name='tilt_2')
-phi_12 = Uniform(name='phi_12', minimum=0, maximum=2 * np.pi)
-phi_jl = Uniform(name='phi_jl', minimum=0, maximum=2 * np.pi)
-luminosity_distance = bilby.gw.prior.UniformComovingVolume(name='luminosity_distance', minimum=1e2, maximum=1e3, unit='Mpc')
-dec = Cosine(name='dec')
-ra = Uniform(name='ra', minimum=0, maximum=2 * np.pi)
-theta_jn = Sine(name='theta_jn')
-psi = Uniform(name='psi', minimum=0, maximum=np.pi)
-phase = Uniform(name='phase', minimum=0, maximum=2 * np.pi)
-geocent_time = Uniform(1126259462.322, 1126259462.522, name='geocent_time', unit='$s$')
+mass_1 = Uniform(name='mass_1', minimum=30, maximum=50, unit='$M_{\\odot}$', periodic_boundary=False)
+mass_2 = Uniform(name='mass_2', minimum=20, maximum=40, unit='$M_{\\odot}$', periodic_boundary=False)
+mass_ratio = Constraint(name='mass_ratio', minimum=0.125, maximum=1)
+a_1 = Uniform(name='a_1', minimum=0, maximum=0.8, periodic_boundary=False)
+a_2 = Uniform(name='a_2', minimum=0, maximum=0.8, periodic_boundary=False)
+tilt_1 = Sine(name='tilt_1', periodic_boundary=False)
+tilt_2 = Sine(name='tilt_2', periodic_boundary=False)
+phi_12 = Uniform(name='phi_12', minimum=0, maximum=2 * np.pi, periodic_boundary=True)
+phi_jl = Uniform(name='phi_jl', minimum=0, maximum=2 * np.pi, periodic_boundary=True)
+luminosity_distance = bilby.gw.prior.UniformComovingVolume(name='luminosity_distance', minimum=1e2, maximum=1e3, unit='Mpc', periodic_boundary=False)
+dec = Cosine(name='dec', periodic_boundary=False)
+ra = Uniform(name='ra', minimum=0, maximum=2 * np.pi, periodic_boundary=True)
+theta_jn = Sine(name='theta_jn', periodic_boundary=False)
+psi = Uniform(name='psi', minimum=0, maximum=np.pi, periodic_boundary=True)
+phase = Uniform(name='phase', minimum=0, maximum=2 * np.pi, periodic_boundary=True)
+geocent_time = Uniform(1126259462.322, 1126259462.522, name='geocent_time', unit='$s$', periodic_boundary=False)
# These are the calibration parameters as described in
# https://journals.aps.org/prx/abstract/10.1103/PhysRevX.6.041015
# recalib_H1_frequency_0 = 20
diff --git a/bilby/gw/prior_files/aligned_spin_binary_black_holes.prior b/bilby/gw/prior_files/aligned_spin_binary_black_holes.prior
index ca49a506b..8be6dcd1e 100644
--- a/bilby/gw/prior_files/aligned_spin_binary_black_holes.prior
+++ b/bilby/gw/prior_files/aligned_spin_binary_black_holes.prior
@@ -2,19 +2,19 @@
# Note that you may wish to use more specific mass and distance parameters.
# These commands are all known to bilby.gw.prior.
# Lines beginning "#" are ignored.
-mass_1 = Uniform(name='mass_1', minimum=5, maximum=100, unit='$M_{\\odot}$')
-mass_2 = Uniform(name='mass_2', minimum=5, maximum=100, unit='$M_{\\odot}$')
-mass_ratio = Constraint(name='mass_ratio', minimum=0.125, maximum=1)
-# chirp_mass = Uniform(name='chirp_mass', minimum=25, maximum=100, unit='$M_{\\odot}$')
-# total_mass = Uniform(name='total_mass', minimum=10, maximum=200, unit='$M_{\\odot}$')
-# mass_ratio = Uniform(name='mass_ratio', minimum=0.125, maximum=1)
-# symmetric_mass_ratio = Uniform(name='symmetric_mass_ratio', minimum=8 / 81, maximum=0.25)
-chi_1 = bilby.gw.prior.AlignedSpin(name='chi_1', a_prior=Uniform(minimum=0, maximum=0.8))
-chi_2 = bilby.gw.prior.AlignedSpin(name='chi_2', a_prior=Uniform(minimum=0, maximum=0.8))
-luminosity_distance = bilby.gw.prior.UniformComovingVolume(name='luminosity_distance', minimum=1e2, maximum=5e3, unit='Mpc')
-dec = Cosine(name='dec')
-ra = Uniform(name='ra', minimum=0, maximum=2 * np.pi)
-theta_jn = Sine(name='theta_jn')
-# cos_theta_jn = Uniform(name='cos_theta_jn', minimum=-1, maximum=1)
-psi = Uniform(name='psi', minimum=0, maximum=np.pi)
-phase = Uniform(name='phase', minimum=0, maximum=2 * np.pi)
+mass_1 = Uniform(name='mass_1', minimum=5, maximum=100, unit='$M_{\\odot}$', periodic_boundary=False)
+mass_2 = Uniform(name='mass_2', minimum=5, maximum=100, unit='$M_{\\odot}$', periodic_boundary=False)
+mass_ratio = Constraint(name='mass_ratio', minimum=0.125, maximum=1)
+# chirp_mass = Uniform(name='chirp_mass', minimum=25, maximum=100, unit='$M_{\\odot}$', periodic_boundary=False)
+# total_mass = Uniform(name='total_mass', minimum=10, maximum=200, unit='$M_{\\odot}$', periodic_boundary=False)
+# mass_ratio = Uniform(name='mass_ratio', minimum=0.125, maximum=1, periodic_boundary=False)
+# symmetric_mass_ratio = Uniform(name='symmetric_mass_ratio', minimum=8 / 81, maximum=0.25, periodic_boundary=False)
+chi_1 = bilby.gw.prior.AlignedSpin(name='chi_1', a_prior=Uniform(minimum=0, maximum=0.8), periodic_boundary=False)
+chi_2 = bilby.gw.prior.AlignedSpin(name='chi_2', a_prior=Uniform(minimum=0, maximum=0.8), periodic_boundary=False)
+luminosity_distance = bilby.gw.prior.UniformComovingVolume(name='luminosity_distance', minimum=1e2, maximum=5e3, unit='Mpc', periodic_boundary=False)
+dec = Cosine(name='dec', periodic_boundary=False)
+ra = Uniform(name='ra', minimum=0, maximum=2 * np.pi, periodic_boundary=True)
+theta_jn = Sine(name='theta_jn', periodic_boundary=False)
+# cos_theta_jn = Uniform(name='cos_theta_jn', minimum=-1, maximum=1, periodic_boundary=False)
+psi = Uniform(name='psi', minimum=0, maximum=np.pi, periodic_boundary=True)
+phase = Uniform(name='phase', minimum=0, maximum=2 * np.pi, periodic_boundary=True)
diff --git a/bilby/gw/prior_files/binary_black_holes.prior b/bilby/gw/prior_files/binary_black_holes.prior
index e79bd5baa..6df97a1dd 100644
--- a/bilby/gw/prior_files/binary_black_holes.prior
+++ b/bilby/gw/prior_files/binary_black_holes.prior
@@ -2,25 +2,25 @@
# Note that you may wish to use more specific mass and distance parameters.
# These commands are all known to bilby.gw.prior.
# Lines beginning "#" are ignored.
-mass_1 = Uniform(name='mass_1', minimum=5, maximum=100, unit='$M_{\\odot}$')
-mass_2 = Uniform(name='mass_2', minimum=5, maximum=100, unit='$M_{\\odot}$')
-mass_ratio = Constraint(name='mass_ratio', minimum=0.125, maximum=1)
-# chirp_mass = Uniform(name='chirp_mass', minimum=25, maximum=100, unit='$M_{\\odot}$')
-# total_mass = Uniform(name='total_mass', minimum=10, maximum=200, unit='$M_{\\odot}$')
-# mass_ratio = Uniform(name='mass_ratio', minimum=0.125, maximum=1)
-# symmetric_mass_ratio = Uniform(name='symmetric_mass_ratio', minimum=8 / 81, maximum=0.25)
-a_1 = Uniform(name='a_1', minimum=0, maximum=0.8)
-a_2 = Uniform(name='a_2', minimum=0, maximum=0.8)
-tilt_1 = Sine(name='tilt_1')
-tilt_2 = Sine(name='tilt_2')
-# cos_tilt_1 = Uniform(name='cos_tilt_1', minimum=-1, maximum=1)
-# cos_tilt_2 = Uniform(name='cos_tilt_2', minimum=-1, maximum=1)
-phi_12 = Uniform(name='phi_12', minimum=0, maximum=2 * np.pi)
-phi_jl = Uniform(name='phi_jl', minimum=0, maximum=2 * np.pi)
-luminosity_distance = bilby.gw.prior.UniformComovingVolume(name='luminosity_distance', minimum=1e2, maximum=5e3, unit='Mpc')
-dec = Cosine(name='dec')
-ra = Uniform(name='ra', minimum=0, maximum=2 * np.pi)
-theta_jn = Sine(name='theta_jn')
-# cos_theta_jn = Uniform(name='cos_theta_jn', minimum=-1, maximum=1)
-psi = Uniform(name='psi', minimum=0, maximum=np.pi)
-phase = Uniform(name='phase', minimum=0, maximum=2 * np.pi)
+mass_1 = Uniform(name='mass_1', minimum=5, maximum=100, unit='$M_{\\odot}$', periodic_boundary=False)
+mass_2 = Uniform(name='mass_2', minimum=5, maximum=100, unit='$M_{\\odot}$', periodic_boundary=False)
+mass_ratio = Constraint(name='mass_ratio', minimum=0.125, maximum=1)
+# chirp_mass = Uniform(name='chirp_mass', minimum=25, maximum=100, unit='$M_{\\odot}$', periodic_boundary=False)
+# total_mass = Uniform(name='total_mass', minimum=10, maximum=200, unit='$M_{\\odot}$', periodic_boundary=False)
+# mass_ratio = Uniform(name='mass_ratio', minimum=0.125, maximum=1, periodic_boundary=False)
+# symmetric_mass_ratio = Uniform(name='symmetric_mass_ratio', minimum=8 / 81, maximum=0.25, periodic_boundary=False)
+a_1 = Uniform(name='a_1', minimum=0, maximum=0.8, periodic_boundary=False)
+a_2 = Uniform(name='a_2', minimum=0, maximum=0.8, periodic_boundary=False)
+tilt_1 = Sine(name='tilt_1', periodic_boundary=False)
+tilt_2 = Sine(name='tilt_2', periodic_boundary=False)
+# cos_tilt_1 = Uniform(name='cos_tilt_1', minimum=-1, maximum=1, periodic_boundary=False)
+# cos_tilt_2 = Uniform(name='cos_tilt_2', minimum=-1, maximum=1, periodic_boundary=False)
+phi_12 = Uniform(name='phi_12', minimum=0, maximum=2 * np.pi, periodic_boundary=True)
+phi_jl = Uniform(name='phi_jl', minimum=0, maximum=2 * np.pi, periodic_boundary=True)
+luminosity_distance = bilby.gw.prior.UniformComovingVolume(name='luminosity_distance', minimum=1e2, maximum=5e3, unit='Mpc', periodic_boundary=False)
+dec = Cosine(name='dec', periodic_boundary=False)
+ra = Uniform(name='ra', minimum=0, maximum=2 * np.pi, periodic_boundary=True)
+theta_jn = Sine(name='theta_jn', periodic_boundary=False)
+# cos_theta_jn = Uniform(name='cos_theta_jn', minimum=-1, maximum=1, periodic_boundary=False)
+psi = Uniform(name='psi', minimum=0, maximum=np.pi, periodic_boundary=True)
+phase = Uniform(name='phase', minimum=0, maximum=2 * np.pi, periodic_boundary=True)
diff --git a/bilby/gw/prior_files/binary_neutron_stars.prior b/bilby/gw/prior_files/binary_neutron_stars.prior
index eff1a2f2a..09acbdd14 100644
--- a/bilby/gw/prior_files/binary_neutron_stars.prior
+++ b/bilby/gw/prior_files/binary_neutron_stars.prior
@@ -2,23 +2,23 @@
# Note that you may wish to use more specific mass and distance parameters.
# These commands are all known to bilby.gw.prior.
# Lines beginning "#" are ignored.
-mass_1 = Uniform(name='mass_1', minimum=1, maximum=2, unit='$M_{\\odot}$')
-mass_2 = Uniform(name='mass_2', minimum=1, maximum=2, unit='$M_{\\odot}$')
-mass_ratio = Constraint(name='mass_ratio', minimum=0.125, maximum=1)
-# chirp_mass = Uniform(name='chirp_mass', minimum=0.87, maximum=1.74, unit='$M_{\\odot}$')
-# total_mass = Uniform(name='total_mass', minimum=2, maximum=4, unit='$M_{\\odot}$')
-# mass_ratio = Uniform(name='mass_ratio', minimum=0.5, maximum=1)
-# symmetric_mass_ratio = Uniform(name='symmetric_mass_ratio', minimum=0.22, maximum=0.25)
-chi_1 = bilby.gw.prior.AlignedSpin(a_prior=Uniform(0, 0.05), z_prior=Uniform(-1, 1), name='chi_1', latex_label='$\\chi_1$')
-chi_2 = bilby.gw.prior.AlignedSpin(a_prior=Uniform(0, 0.05), z_prior=Uniform(-1, 1), name='chi_2', latex_label='$\\chi_2$')
-luminosity_distance = bilby.gw.prior.UniformComovingVolume(name='luminosity_distance', minimum=10, maximum=500, unit='Mpc')
-dec = Cosine(name='dec')
-ra = Uniform(name='ra', minimum=0, maximum=2 * np.pi)
-theta_jn = Sine(name='theta_jn')
-# cos_theta_jn = Uniform(name='cos_theta_jn', minimum=-1, maximum=1)
-psi = Uniform(name='psi', minimum=0, maximum=np.pi)
-phase = Uniform(name='phase', minimum=0, maximum=2 * np.pi)
-lambda_1 = Uniform(name='lambda_1', minimum=0, maximum=3000 )
-lambda_2 = Uniform(name='lambda_2', minimum=0, maximum=3000 )
-# lambda_tilde = Uniform(name='lambda_tilde', minimum=0, maximum=5000)
-# delta_lambda = Uniform(name='delta_lambda', minimum=-5000, maximum=5000)
+mass_1 = Uniform(name='mass_1', minimum=1, maximum=2, unit='$M_{\\odot}$', periodic_boundary=False)
+mass_2 = Uniform(name='mass_2', minimum=1, maximum=2, unit='$M_{\\odot}$', periodic_boundary=False)
+mass_ratio = Constraint(name='mass_ratio', minimum=0.125, maximum=1)
+# chirp_mass = Uniform(name='chirp_mass', minimum=0.87, maximum=1.74, unit='$M_{\\odot}$', periodic_boundary=False)
+# total_mass = Uniform(name='total_mass', minimum=2, maximum=4, unit='$M_{\\odot}$', periodic_boundary=False)
+# mass_ratio = Uniform(name='mass_ratio', minimum=0.5, maximum=1, periodic_boundary=False)
+# symmetric_mass_ratio = Uniform(name='symmetric_mass_ratio', minimum=0.22, maximum=0.25, periodic_boundary=False)
+chi_1 = bilby.gw.prior.AlignedSpin(a_prior=Uniform(0, 0.05), z_prior=Uniform(-1, 1), name='chi_1', latex_label='$\\chi_1$', periodic_boundary=False)
+chi_2 = bilby.gw.prior.AlignedSpin(a_prior=Uniform(0, 0.05), z_prior=Uniform(-1, 1), name='chi_2', latex_label='$\\chi_2$', periodic_boundary=False)
+luminosity_distance = bilby.gw.prior.UniformComovingVolume(name='luminosity_distance', minimum=10, maximum=500, unit='Mpc', periodic_boundary=False)
+dec = Cosine(name='dec', periodic_boundary=False)
+ra = Uniform(name='ra', minimum=0, maximum=2 * np.pi, periodic_boundary=True)
+theta_jn = Sine(name='theta_jn', periodic_boundary=False)
+# cos_theta_jn = Uniform(name='cos_theta_jn', minimum=-1, maximum=1, periodic_boundary=False)
+psi = Uniform(name='psi', minimum=0, maximum=np.pi, periodic_boundary=True)
+phase = Uniform(name='phase', minimum=0, maximum=2 * np.pi, periodic_boundary=True)
+lambda_1 = Uniform(name='lambda_1', minimum=0, maximum=3000, periodic_boundary=False)
+lambda_2 = Uniform(name='lambda_2', minimum=0, maximum=3000, periodic_boundary=False)
+# lambda_tilde = Uniform(name='lambda_tilde', minimum=0, maximum=5000, periodic_boundary=False)
+# delta_lambda = Uniform(name='delta_lambda', minimum=-5000, maximum=5000, periodic_boundary=False)
diff --git a/bilby/gw/sampler/__init__.py b/bilby/gw/sampler/__init__.py
new file mode 100644
index 000000000..25ad2ca38
--- /dev/null
+++ b/bilby/gw/sampler/__init__.py
@@ -0,0 +1,2 @@
+from __future__ import absolute_import
+from . import proposal
diff --git a/bilby/gw/sampler/proposal.py b/bilby/gw/sampler/proposal.py
new file mode 100644
index 000000000..328fd98a8
--- /dev/null
+++ b/bilby/gw/sampler/proposal.py
@@ -0,0 +1,49 @@
+import random
+
+import numpy as np
+
+from bilby.core.sampler.proposal import JumpProposal
+
+
+class SkyLocationWanderJump(JumpProposal):
+ """
+ Jump proposal for wandering over the sky location. Does a Gaussian step in
+ RA and DEC depending on the temperature.
+ """
+
+ def __call__(self, sample, **kwargs):
+ temperature = 1 / kwargs.get('inverse_temperature', 1.0)
+ sigma = np.sqrt(temperature) / 2 / np.pi
+ sample['ra'] += random.gauss(0, sigma)
+ sample['dec'] += random.gauss(0, sigma)
+ return super(SkyLocationWanderJump, self).__call__(sample)
+
+
+class CorrelatedPolarisationPhaseJump(JumpProposal):
+ """
+ Correlated polarisation/phase jump proposal. Jumps between degenerate phi/psi regions.
+ """
+
+ def __call__(self, sample, **kwargs):
+ alpha = sample['psi'] + sample['phase']
+ beta = sample['psi'] - sample['phase']
+
+ draw = random.random()
+ if draw < 0.5:
+ alpha = 3.0 * np.pi * random.random()
+ else:
+ beta = 3.0 * np.pi * random.random() - 2 * np.pi
+ sample['psi'] = (alpha + beta) * 0.5
+ sample['phase'] = (alpha - beta) * 0.5
+ return super(CorrelatedPolarisationPhaseJump, self).__call__(sample)
+
+
+class PolarisationPhaseJump(JumpProposal):
+ """
+ Correlated polarisation/phase jump proposal. Jumps between degenerate phi/psi regions.
+ """
+
+ def __call__(self, sample, **kwargs):
+ sample['phase'] += np.pi
+ sample['psi'] += np.pi / 2
+ return super(PolarisationPhaseJump, self).__call__(sample)
diff --git a/examples/injection_examples/custom_proposal_example.py b/examples/injection_examples/custom_proposal_example.py
new file mode 100644
index 000000000..8480dd361
--- /dev/null
+++ b/examples/injection_examples/custom_proposal_example.py
@@ -0,0 +1,68 @@
+#!/usr/bin/env python
+"""
+Tutorial for running cpnest with custom jump proposals.
+"""
+from __future__ import division, print_function
+
+import numpy as np
+import bilby.gw.sampler.proposal
+from bilby.core.sampler import proposal
+
+
+# The set up here is the same as in fast_tutorial.py. Look there for descriptive explanations.
+
+duration = 4.
+sampling_frequency = 2048.
+
+outdir = 'outdir'
+label = 'custom_jump_proposals'
+bilby.core.utils.setup_logger(outdir=outdir, label=label)
+
+np.random.seed(88170235)
+
+injection_parameters = dict(
+ mass_1=36., mass_2=29., a_1=0.4, a_2=0.3, tilt_1=0.5, tilt_2=1.0,
+ phi_12=1.7, phi_jl=0.3, luminosity_distance=2000., theta_jn=0.4, psi=2.659,
+ phase=1.3, geocent_time=1126259642.413, ra=1.375, dec=-1.2108)
+waveform_arguments = dict(waveform_approximant='IMRPhenomPv2',
+ reference_frequency=50., minimum_frequency=20.)
+waveform_generator = bilby.gw.WaveformGenerator(
+ duration=duration, sampling_frequency=sampling_frequency,
+ frequency_domain_source_model=bilby.gw.source.lal_binary_black_hole,
+ waveform_arguments=waveform_arguments)
+ifos = bilby.gw.detector.InterferometerList(['H1', 'L1'])
+ifos.set_strain_data_from_power_spectral_densities(
+ sampling_frequency=sampling_frequency, duration=duration,
+ start_time=injection_parameters['geocent_time'] - 3)
+ifos.inject_signal(waveform_generator=waveform_generator,
+ parameters=injection_parameters)
+priors = bilby.gw.prior.BBHPriorDict()
+priors['geocent_time'] = bilby.core.prior.Uniform(
+ minimum=injection_parameters['geocent_time'] - 1,
+ maximum=injection_parameters['geocent_time'] + 1,
+ name='geocent_time', latex_label='$t_c$', unit='$s$')
+for key in ['a_1', 'a_2', 'tilt_1', 'tilt_2', 'phi_12', 'phi_jl', 'geocent_time']:
+ priors[key] = injection_parameters[key]
+likelihood = bilby.gw.GravitationalWaveTransient(
+ interferometers=ifos, waveform_generator=waveform_generator)
+
+# Definition of the custom jump proposals. Define a JumpProposalCycle. The first argument is a list
+# of all allowed jump proposals. The second argument is a list of weights for the respective jump
+# proposal.
+
+jump_proposals = proposal.JumpProposalCycle(
+ [proposal.EnsembleWalk(priors=priors), proposal.EnsembleStretch(priors=priors),
+ proposal.DifferentialEvolution(priors=priors), proposal.EnsembleEigenVector(priors=priors),
+ bilby.gw.sampler.proposal.SkyLocationWanderJump(priors=priors), bilby.gw.sampler.proposal.CorrelatedPolarisationPhaseJump(priors=priors),
+ bilby.gw.sampler.proposal.PolarisationPhaseJump(priors=priors), proposal.DrawFlatPrior(priors=priors)],
+ weights=[2, 2, 5, 1, 1, 1, 1, 1])
+
+# Run cpnest with the proposals kwarg specified.
+# Make sure to have a version of cpnest installed that supports custom proposals.
+result = bilby.run_sampler(
+ likelihood=likelihood, priors=priors, sampler='cpnest', npoints=1000,
+ injection_parameters=injection_parameters, outdir=outdir, label=label,
+ proposals=jump_proposals)
+
+# Make a corner plot.
+result.plot_corner()
diff --git a/test/gw_prior_test.py b/test/gw_prior_test.py
index 893d4f3de..485341dbf 100644
--- a/test/gw_prior_test.py
+++ b/test/gw_prior_test.py
@@ -36,8 +36,10 @@ class TestBBHPriorDict(unittest.TestCase):
for key in default.keys()])
names = all([self.bbh_prior_dict[key].name == default[key].name
for key in default.keys()])
+ boundaries = all([self.bbh_prior_dict[key].periodic_boundary is default[key].periodic_boundary
+ for key in default.keys()])
- self.assertTrue(all([minima, maxima, names]))
+ self.assertTrue(all([minima, maxima, names, boundaries]))
def test_create_from_dict(self):
new_dict = bilby.gw.prior.BBHPriorDict(dictionary=self.prior_dict)
@@ -135,8 +137,10 @@ class TestBNSPriorDict(unittest.TestCase):
for key in default.keys()])
names = all([self.bns_prior_dict[key].name == default[key].name
for key in default.keys()])
+ boundaries = all([self.bns_prior_dict[key].periodic_boundary == default[key].periodic_boundary
+ for key in default.keys()])
- self.assertTrue(all([minima, maxima, names]))
+ self.assertTrue(all([minima, maxima, names, boundaries]))
def test_create_from_dict(self):
new_dict = bilby.gw.prior.BNSPriorDict(dictionary=self.prior_dict)
diff --git a/test/prior_files/binary_black_holes.prior b/test/prior_files/binary_black_holes.prior
index e79bd5baa..b9abc7a45 100644
--- a/test/prior_files/binary_black_holes.prior
+++ b/test/prior_files/binary_black_holes.prior
@@ -2,25 +2,25 @@
# Note that you may wish to use more specific mass and distance parameters.
# These commands are all known to bilby.gw.prior.
# Lines beginning "#" are ignored.
-mass_1 = Uniform(name='mass_1', minimum=5, maximum=100, unit='$M_{\\odot}$')
-mass_2 = Uniform(name='mass_2', minimum=5, maximum=100, unit='$M_{\\odot}$')
-mass_ratio = Constraint(name='mass_ratio', minimum=0.125, maximum=1)
-# chirp_mass = Uniform(name='chirp_mass', minimum=25, maximum=100, unit='$M_{\\odot}$')
-# total_mass = Uniform(name='total_mass', minimum=10, maximum=200, unit='$M_{\\odot}$')
-# mass_ratio = Uniform(name='mass_ratio', minimum=0.125, maximum=1)
-# symmetric_mass_ratio = Uniform(name='symmetric_mass_ratio', minimum=8 / 81, maximum=0.25)
-a_1 = Uniform(name='a_1', minimum=0, maximum=0.8)
-a_2 = Uniform(name='a_2', minimum=0, maximum=0.8)
-tilt_1 = Sine(name='tilt_1')
-tilt_2 = Sine(name='tilt_2')
-# cos_tilt_1 = Uniform(name='cos_tilt_1', minimum=-1, maximum=1)
-# cos_tilt_2 = Uniform(name='cos_tilt_2', minimum=-1, maximum=1)
-phi_12 = Uniform(name='phi_12', minimum=0, maximum=2 * np.pi)
-phi_jl = Uniform(name='phi_jl', minimum=0, maximum=2 * np.pi)
-luminosity_distance = bilby.gw.prior.UniformComovingVolume(name='luminosity_distance', minimum=1e2, maximum=5e3, unit='Mpc')
-dec = Cosine(name='dec')
-ra = Uniform(name='ra', minimum=0, maximum=2 * np.pi)
-theta_jn = Sine(name='theta_jn')
-# cos_theta_jn = Uniform(name='cos_theta_jn', minimum=-1, maximum=1)
-psi = Uniform(name='psi', minimum=0, maximum=np.pi)
-phase = Uniform(name='phase', minimum=0, maximum=2 * np.pi)
+mass_1 = Uniform(name='mass_1', minimum=5, maximum=100, unit='$M_{\\odot}$', periodic_boundary=False)
+mass_2 = Uniform(name='mass_2', minimum=5, maximum=100, unit='$M_{\\odot}$', periodic_boundary=False)
+mass_ratio = Constraint(name='mass_ratio', minimum=0.125, maximum=1)
+# chirp_mass = Uniform(name='chirp_mass', minimum=25, maximum=100, unit='$M_{\\odot}$', periodic_boundary=False)
+# total_mass = Uniform(name='total_mass', minimum=10, maximum=200, unit='$M_{\\odot}$', periodic_boundary=False)
+# mass_ratio = Uniform(name='mass_ratio', minimum=0.125, maximum=1, periodic_boundary=False)
+# symmetric_mass_ratio = Uniform(name='symmetric_mass_ratio', minimum=8 / 81, maximum=0.25, periodic_boundary=False)
+a_1 = Uniform(name='a_1', minimum=0, maximum=0.8, periodic_boundary=False)
+a_2 = Uniform(name='a_2', minimum=0, maximum=0.8, periodic_boundary=False)
+tilt_1 = Sine(name='tilt_1', periodic_boundary=False)
+tilt_2 = Sine(name='tilt_2', periodic_boundary=False)
+# cos_tilt_1 = Uniform(name='cos_tilt_1', minimum=-1, maximum=1, periodic_boundary=False)
+# cos_tilt_2 = Uniform(name='cos_tilt_2', minimum=-1, maximum=1, periodic_boundary=False)
+phi_12 = Uniform(name='phi_12', minimum=0, maximum=2 * np.pi, periodic_boundary=True)
+phi_jl = Uniform(name='phi_jl', minimum=0, maximum=2 * np.pi, periodic_boundary=True)
+luminosity_distance = bilby.gw.prior.UniformComovingVolume(name='luminosity_distance', minimum=1e2, maximum=5e3, unit='Mpc', periodic_boundary=False)
+dec = Cosine(name='dec', periodic_boundary=False)
+ra = Uniform(name='ra', minimum=0, maximum=2 * np.pi, periodic_boundary=True)
+theta_jn = Sine(name='theta_jn', periodic_boundary=False)
+# cos_theta_jn = Uniform(name='cos_theta_jn', minimum=-1, maximum=1, periodic_boundary=False)
+psi = Uniform(name='psi', minimum=0, maximum=np.pi, periodic_boundary=True)
+phase = Uniform(name='phase', minimum=0, maximum=2 * np.pi, periodic_boundary=True)
diff --git a/test/prior_files/binary_neutron_stars.prior b/test/prior_files/binary_neutron_stars.prior
index eff1a2f2a..09acbdd14 100644
--- a/test/prior_files/binary_neutron_stars.prior
+++ b/test/prior_files/binary_neutron_stars.prior
@@ -2,23 +2,23 @@
# Note that you may wish to use more specific mass and distance parameters.
# These commands are all known to bilby.gw.prior.
# Lines beginning "#" are ignored.
-mass_1 = Uniform(name='mass_1', minimum=1, maximum=2, unit='$M_{\\odot}$')
-mass_2 = Uniform(name='mass_2', minimum=1, maximum=2, unit='$M_{\\odot}$')
-mass_ratio = Constraint(name='mass_ratio', minimum=0.125, maximum=1)
-# chirp_mass = Uniform(name='chirp_mass', minimum=0.87, maximum=1.74, unit='$M_{\\odot}$')
-# total_mass = Uniform(name='total_mass', minimum=2, maximum=4, unit='$M_{\\odot}$')
-# mass_ratio = Uniform(name='mass_ratio', minimum=0.5, maximum=1)
-# symmetric_mass_ratio = Uniform(name='symmetric_mass_ratio', minimum=0.22, maximum=0.25)
-chi_1 = bilby.gw.prior.AlignedSpin(a_prior=Uniform(0, 0.05), z_prior=Uniform(-1, 1), name='chi_1', latex_label='$\\chi_1$')
-chi_2 = bilby.gw.prior.AlignedSpin(a_prior=Uniform(0, 0.05), z_prior=Uniform(-1, 1), name='chi_2', latex_label='$\\chi_2$')
-luminosity_distance = bilby.gw.prior.UniformComovingVolume(name='luminosity_distance', minimum=10, maximum=500, unit='Mpc')
-dec = Cosine(name='dec')
-ra = Uniform(name='ra', minimum=0, maximum=2 * np.pi)
-theta_jn = Sine(name='theta_jn')
-# cos_theta_jn = Uniform(name='cos_theta_jn', minimum=-1, maximum=1)
-psi = Uniform(name='psi', minimum=0, maximum=np.pi)
-phase = Uniform(name='phase', minimum=0, maximum=2 * np.pi)
-lambda_1 = Uniform(name='lambda_1', minimum=0, maximum=3000 )
-lambda_2 = Uniform(name='lambda_2', minimum=0, maximum=3000 )
-# lambda_tilde = Uniform(name='lambda_tilde', minimum=0, maximum=5000)
-# delta_lambda = Uniform(name='delta_lambda', minimum=-5000, maximum=5000)
+mass_1 = Uniform(name='mass_1', minimum=1, maximum=2, unit='$M_{\\odot}$', periodic_boundary=False)
+mass_2 = Uniform(name='mass_2', minimum=1, maximum=2, unit='$M_{\\odot}$', periodic_boundary=False)
+mass_ratio = Constraint(name='mass_ratio', minimum=0.125, maximum=1)
+# chirp_mass = Uniform(name='chirp_mass', minimum=0.87, maximum=1.74, unit='$M_{\\odot}$', periodic_boundary=False)
+# total_mass = Uniform(name='total_mass', minimum=2, maximum=4, unit='$M_{\\odot}$', periodic_boundary=False)
+# mass_ratio = Uniform(name='mass_ratio', minimum=0.5, maximum=1, periodic_boundary=False)
+# symmetric_mass_ratio = Uniform(name='symmetric_mass_ratio', minimum=0.22, maximum=0.25, periodic_boundary=False)
+chi_1 = bilby.gw.prior.AlignedSpin(a_prior=Uniform(0, 0.05), z_prior=Uniform(-1, 1), name='chi_1', latex_label='$\\chi_1$', periodic_boundary=False)
+chi_2 = bilby.gw.prior.AlignedSpin(a_prior=Uniform(0, 0.05), z_prior=Uniform(-1, 1), name='chi_2', latex_label='$\\chi_2$', periodic_boundary=False)
+luminosity_distance = bilby.gw.prior.UniformComovingVolume(name='luminosity_distance', minimum=10, maximum=500, unit='Mpc', periodic_boundary=False)
+dec = Cosine(name='dec', periodic_boundary=False)
+ra = Uniform(name='ra', minimum=0, maximum=2 * np.pi, periodic_boundary=True)
+theta_jn = Sine(name='theta_jn', periodic_boundary=False)
+# cos_theta_jn = Uniform(name='cos_theta_jn', minimum=-1, maximum=1, periodic_boundary=False)
+psi = Uniform(name='psi', minimum=0, maximum=np.pi, periodic_boundary=True)
+phase = Uniform(name='phase', minimum=0, maximum=2 * np.pi, periodic_boundary=True)
+lambda_1 = Uniform(name='lambda_1', minimum=0, maximum=3000, periodic_boundary=False)
+lambda_2 = Uniform(name='lambda_2', minimum=0, maximum=3000, periodic_boundary=False)
+# lambda_tilde = Uniform(name='lambda_tilde', minimum=0, maximum=5000, periodic_boundary=False)
+# delta_lambda = Uniform(name='delta_lambda', minimum=-5000, maximum=5000, periodic_boundary=False)
diff --git a/test/prior_test.py b/test/prior_test.py
index 6a2afc750..de5e49fd8 100644
--- a/test/prior_test.py
+++ b/test/prior_test.py
@@ -35,8 +35,10 @@ class TestPriorInstantiationWithoutOptionalPriors(unittest.TestCase):
self.assertIsNone(self.prior.rescale(1))
def test_base_repr(self):
- self.prior = bilby.core.prior.Prior(name='test_name', latex_label='test_label', minimum=0, maximum=1)
- expected_string = "Prior(name='test_name', latex_label='test_label', unit=None, minimum=0, maximum=1)"
+ self.prior = bilby.core.prior.Prior(name='test_name', latex_label='test_label', minimum=0, maximum=1,
+ periodic_boundary=False)
+ expected_string = "Prior(name='test_name', latex_label='test_label', unit=None, minimum=0, maximum=1, " \
+ "periodic_boundary=False)"
self.assertEqual(expected_string, self.prior.__repr__())
def test_base_prob(self):
@@ -60,6 +62,9 @@ class TestPriorInstantiationWithoutOptionalPriors(unittest.TestCase):
self.assertFalse(self.prior.is_in_prior_range(val_below))
self.assertFalse(self.prior.is_in_prior_range(val_above))
+ def test_periodic_boundary_is_false(self):
+ self.assertFalse(self.prior.periodic_boundary)
+
class TestPriorName(unittest.TestCase):
@@ -105,7 +110,7 @@ class TestPriorIsFixed(unittest.TestCase):
pass
def tearDown(self):
- pass
+ del self.prior
def test_is_fixed_parent_class(self):
self.prior = bilby.core.prior.Prior()
@@ -120,6 +125,23 @@ class TestPriorIsFixed(unittest.TestCase):
self.assertFalse(self.prior.is_fixed)
+class TestPriorBoundary(unittest.TestCase):
+
+ def setUp(self):
+ self.prior = bilby.core.prior.Prior(periodic_boundary=False)
+
+ def tearDown(self):
+ del self.prior
+
+ def test_set_boundary_valid(self):
+ self.prior.periodic_boundary = True
+ self.assertTrue(self.prior.periodic_boundary)
+
+ def test_set_boundary_invalid(self):
+ with self.assertRaises(ValueError):
+ self.prior.periodic_boundary = 'else'
+
+
class TestPriorClasses(unittest.TestCase):
def setUp(self):
@@ -153,7 +175,6 @@ class TestPriorClasses(unittest.TestCase):
bilby.core.prior.Lorentzian(name='test', unit='unit', alpha=0, beta=1),
bilby.core.prior.Gamma(name='test', unit='unit', k=1, theta=1),
bilby.core.prior.ChiSquared(name='test', unit='unit', nu=2),
- bilby.core.prior.FermiDirac(name='test', unit='unit', sigma=1., r=10.),
bilby.gw.prior.AlignedSpin(name='test', unit='unit'),
]
@@ -357,8 +378,9 @@ class TestPriorClasses(unittest.TestCase):
class TestPriorDict(unittest.TestCase):
def setUp(self):
- self.first_prior = bilby.core.prior.Uniform(name='a', minimum=0, maximum=1, unit='kg')
- self.second_prior = bilby.core.prior.PowerLaw(name='b', alpha=3, minimum=1, maximum=2, unit='m/s')
+ self.first_prior = bilby.core.prior.Uniform(name='a', minimum=0, maximum=1, unit='kg', periodic_boundary=False)
+ self.second_prior = bilby.core.prior.PowerLaw(name='b', alpha=3, minimum=1, maximum=2, unit='m/s',
+ periodic_boundary=False)
self.third_prior = bilby.core.prior.DeltaFunction(name='c', peak=42, unit='m')
self.priors = dict(mass=self.first_prior,
speed=self.second_prior,
@@ -398,43 +420,46 @@ class TestPriorDict(unittest.TestCase):
def test_read_from_file(self):
expected = dict(
mass_1=bilby.core.prior.Uniform(
- name='mass_1', minimum=5, maximum=100, unit='$M_{\\odot}$'),
+ name='mass_1', minimum=5, maximum=100, unit='$M_{\\odot}$', periodic_boundary=False),
mass_2=bilby.core.prior.Uniform(
- name='mass_2', minimum=5, maximum=100, unit='$M_{\\odot}$'),
- mass_ratio=bilby.core.prior.Constraint(
- name='mass_ratio', minimum=0.125, maximum=1.0),
- a_1=bilby.core.prior.Uniform(name='a_1', minimum=0, maximum=0.8),
- a_2=bilby.core.prior.Uniform(name='a_2', minimum=0, maximum=0.8),
- tilt_1=bilby.core.prior.Sine(name='tilt_1'),
- tilt_2=bilby.core.prior.Sine(name='tilt_2'),
+ name='mass_2', minimum=5, maximum=100, unit='$M_{\\odot}$', periodic_boundary=False),
+ mass_ratio=bilby.core.prior.Constraint(name='mass_ratio', minimum=0.125, maximum=1, latex_label='$q$',
+ unit=None),
+ a_1=bilby.core.prior.Uniform(name='a_1', minimum=0, maximum=0.8, periodic_boundary=False),
+ a_2=bilby.core.prior.Uniform(name='a_2', minimum=0, maximum=0.8, periodic_boundary=False),
+ tilt_1=bilby.core.prior.Sine(name='tilt_1', periodic_boundary=False),
+ tilt_2=bilby.core.prior.Sine(name='tilt_2', periodic_boundary=False),
phi_12=bilby.core.prior.Uniform(
- name='phi_12', minimum=0, maximum=2 * np.pi),
+ name='phi_12', minimum=0, maximum=2 * np.pi, periodic_boundary=True),
phi_jl=bilby.core.prior.Uniform(
- name='phi_jl', minimum=0, maximum=2 * np.pi),
+ name='phi_jl', minimum=0, maximum=2 * np.pi, periodic_boundary=True),
luminosity_distance=bilby.gw.prior.UniformComovingVolume(
name='luminosity_distance', minimum=1e2,
- maximum=5e3, unit='Mpc'),
- dec=bilby.core.prior.Cosine(name='dec'),
+ maximum=5e3, unit='Mpc', periodic_boundary=False),
+ dec=bilby.core.prior.Cosine(name='dec', periodic_boundary=False),
ra=bilby.core.prior.Uniform(
- name='ra', minimum=0, maximum=2 * np.pi),
- theta_jn=bilby.core.prior.Sine(name='theta_jn'),
- psi=bilby.core.prior.Uniform(name='psi', minimum=0, maximum=np.pi),
+ name='ra', minimum=0, maximum=2 * np.pi, periodic_boundary=True),
+ theta_jn=bilby.core.prior.Sine(name='theta_jn', periodic_boundary=False),
+ psi=bilby.core.prior.Uniform(name='psi', minimum=0, maximum=np.pi, periodic_boundary=True),
phase=bilby.core.prior.Uniform(
- name='phase', minimum=0, maximum=2 * np.pi)
+ name='phase', minimum=0, maximum=2 * np.pi, periodic_boundary=True)
)
self.assertDictEqual(expected, self.prior_set_from_file)
def test_to_file(self):
expected = ["length = DeltaFunction(peak=42, name='c', latex_label='c', unit='m')\n",
- "speed = PowerLaw(alpha=3, minimum=1, maximum=2, name='b', latex_label='b', unit='m/s')\n",
- "mass = Uniform(minimum=0, maximum=1, name='a', latex_label='a', unit='kg')\n"]
+ "speed = PowerLaw(alpha=3, minimum=1, maximum=2, name='b', latex_label='b', "
+ "unit='m/s', periodic_boundary=False)\n",
+ "mass = Uniform(minimum=0, maximum=1, name='a', latex_label='a', "
+ "unit='kg', periodic_boundary=False)\n"]
self.prior_set_from_dict.to_file(outdir='prior_files', label='to_file_test')
with open('prior_files/to_file_test.prior') as f:
for i, line in enumerate(f.readlines()):
self.assertTrue(line in expected)
def test_from_dict_with_string(self):
- string_prior = "bilby.core.prior.PowerLaw(name='b', alpha=3, minimum=1, maximum=2, unit='m/s')"
+ string_prior = "bilby.core.prior.PowerLaw(name='b', alpha=3, minimum=1, maximum=2, unit='m/s', " \
+ "periodic_boundary=False)"
self.priors['speed'] = string_prior
from_dict = bilby.core.prior.PriorDict(dictionary=self.priors)
self.assertDictEqual(self.prior_set_from_dict, from_dict)
@@ -444,8 +469,10 @@ class TestPriorDict(unittest.TestCase):
self.prior_set_from_dict['e'] = 7.3
self.prior_set_from_dict['f'] = 'unconvertable'
self.prior_set_from_dict.convert_floats_to_delta_functions()
- expected = dict(mass=bilby.core.prior.Uniform(name='a', minimum=0, maximum=1, unit='kg'),
- speed=bilby.core.prior.PowerLaw(name='b', alpha=3, minimum=1, maximum=2, unit='m/s'),
+ expected = dict(mass=bilby.core.prior.Uniform(name='a', minimum=0, maximum=1, unit='kg',
+ periodic_boundary=False),
+ speed=bilby.core.prior.PowerLaw(name='b', alpha=3, minimum=1, maximum=2, unit='m/s',
+ periodic_boundary=False),
length=bilby.core.prior.DeltaFunction(name='c', peak=42, unit='m'),
d=bilby.core.prior.DeltaFunction(peak=5),
e=bilby.core.prior.DeltaFunction(peak=7.3),
@@ -454,33 +481,37 @@ class TestPriorDict(unittest.TestCase):
def test_prior_set_from_dict_but_using_a_string(self):
prior_set = bilby.core.prior.PriorDict(dictionary=self.default_prior_file)
- expected = dict(
- mass_1=bilby.core.prior.Uniform(
- name='mass_1', minimum=5, maximum=100, unit='$M_{\\odot}$'),
- mass_2=bilby.core.prior.Uniform(
- name='mass_2', minimum=5, maximum=100, unit='$M_{\\odot}$'),
- mass_ratio=bilby.core.prior.Constraint(
- name='mass_ratio', minimum=0.125, maximum=1.0),
- a_1=bilby.core.prior.Uniform(name='a_1', minimum=0, maximum=0.8),
- a_2=bilby.core.prior.Uniform(name='a_2', minimum=0, maximum=0.8),
- tilt_1=bilby.core.prior.Sine(name='tilt_1'),
- tilt_2=bilby.core.prior.Sine(name='tilt_2'),
- phi_12=bilby.core.prior.Uniform(
- name='phi_12', minimum=0, maximum=2 * np.pi),
- phi_jl=bilby.core.prior.Uniform(
- name='phi_jl', minimum=0, maximum=2 * np.pi),
- luminosity_distance=bilby.gw.prior.UniformComovingVolume(
- name='luminosity_distance', minimum=1e2,
- maximum=5e3, unit='Mpc'),
- dec=bilby.core.prior.Cosine(name='dec'),
- ra=bilby.core.prior.Uniform(
- name='ra', minimum=0, maximum=2 * np.pi),
- theta_jn=bilby.core.prior.Sine(name='theta_jn'),
- psi=bilby.core.prior.Uniform(name='psi', minimum=0, maximum=np.pi),
- phase=bilby.core.prior.Uniform(
- name='phase', minimum=0, maximum=2 * np.pi)
+ expected = bilby.core.prior.PriorDict(
+ dict(
+ mass_1=bilby.core.prior.Uniform(
+ name='mass_1', minimum=5, maximum=100, unit='$M_{\\odot}$', periodic_boundary=False),
+ mass_2=bilby.core.prior.Uniform(
+ name='mass_2', minimum=5, maximum=100, unit='$M_{\\odot}$', periodic_boundary=False),
+ mass_ratio=bilby.core.prior.Constraint(name='mass_ratio', minimum=0.125, maximum=1, latex_label='$q$',
+ unit=None),
+ a_1=bilby.core.prior.Uniform(name='a_1', minimum=0, maximum=0.8, periodic_boundary=False),
+ a_2=bilby.core.prior.Uniform(name='a_2', minimum=0, maximum=0.8, periodic_boundary=False),
+ tilt_1=bilby.core.prior.Sine(name='tilt_1', periodic_boundary=False),
+ tilt_2=bilby.core.prior.Sine(name='tilt_2', periodic_boundary=False),
+ phi_12=bilby.core.prior.Uniform(
+ name='phi_12', minimum=0, maximum=2 * np.pi, periodic_boundary=True),
+ phi_jl=bilby.core.prior.Uniform(
+ name='phi_jl', minimum=0, maximum=2 * np.pi, periodic_boundary=True),
+ luminosity_distance=bilby.gw.prior.UniformComovingVolume(
+ name='luminosity_distance', minimum=1e2,
+ maximum=5e3, unit='Mpc', periodic_boundary=False),
+ dec=bilby.core.prior.Cosine(name='dec', periodic_boundary=False),
+ ra=bilby.core.prior.Uniform(
+ name='ra', minimum=0, maximum=2 * np.pi, periodic_boundary=True),
+ theta_jn=bilby.core.prior.Sine(name='theta_jn', periodic_boundary=False),
+ psi=bilby.core.prior.Uniform(name='psi', minimum=0, maximum=np.pi, periodic_boundary=True),
+ phase=bilby.core.prior.Uniform(
+ name='phase', minimum=0, maximum=2 * np.pi, periodic_boundary=True)
+ )
)
- self.assertDictEqual(expected, prior_set)
+ all_keys = set(prior_set.keys()).union(set(expected.keys()))
+ for key in all_keys:
+ self.assertEqual(expected[key], prior_set[key])
def test_dict_argument_is_not_string_or_dict(self):
with self.assertRaises(ValueError):
diff --git a/test/proposal_test.py b/test/proposal_test.py
new file mode 100644
index 000000000..dc86444c4
--- /dev/null
+++ b/test/proposal_test.py
@@ -0,0 +1,433 @@
+import unittest
+import mock
+import random
+
+import numpy as np
+
+import bilby.gw.sampler.proposal
+from bilby.core import prior
+from bilby.core.sampler import proposal
+
+
+class TestSample(unittest.TestCase):
+
+ def setUp(self):
+ self.sample = proposal.Sample(dict(a=1, c=2))
+
+ def tearDown(self):
+ del self.sample
+
+ def test_add_sample(self):
+ other = proposal.Sample(dict(a=2, c=5))
+ expected = proposal.Sample(dict(a=3, c=7))
+ self.assertDictEqual(expected, self.sample + other)
+
+ def test_subtract_sample(self):
+ other = proposal.Sample(dict(a=2, c=5))
+ expected = proposal.Sample(dict(a=-1, c=-3))
+ self.assertDictEqual(expected, self.sample - other)
+
+ def test_multiply_sample(self):
+ other = 2
+ expected = proposal.Sample(dict(a=2, c=4))
+ self.assertDictEqual(expected, self.sample * other)
+
+
+class TestJumpProposal(unittest.TestCase):
+
+ def setUp(self):
+ self.priors = prior.PriorDict(dict(reflecting=prior.Uniform(minimum=-0.5, maximum=1, periodic_boundary=False),
+ periodic=prior.Uniform(minimum=-0.5, maximum=1, periodic_boundary=True),
+ default=prior.Uniform(minimum=-0.5, maximum=1)))
+ self.sample_above = dict(reflecting=1.1, periodic=1.1, default=1.1)
+ self.sample_below = dict(reflecting=-0.6, periodic=-0.6, default=-0.6)
+ self.sample_way_above_case1 = dict(reflecting=272, periodic=272, default=272)
+ self.sample_way_above_case2 = dict(reflecting=270.1, periodic=270.1, default=270.1)
+ self.sample_way_below_case1 = dict(reflecting=-274, periodic=-274.1, default=-274)
+ self.sample_way_below_case2 = dict(reflecting=-273.1, periodic=-273.1, default=-273.1)
+ self.jump_proposal = proposal.JumpProposal(priors=self.priors)
+
+ def tearDown(self):
+ del self.priors
+ del self.sample_above
+ del self.sample_below
+ del self.sample_way_above_case1
+ del self.sample_way_above_case2
+ del self.sample_way_below_case1
+ del self.sample_way_below_case2
+ del self.jump_proposal
+
+ def test_set_get_log_j(self):
+ self.jump_proposal.log_j = 2.3
+ self.assertEqual(2.3, self.jump_proposal.log_j)
+
+ def test_boundary_above_reflecting(self):
+ new_sample = self.jump_proposal(self.sample_above)
+ self.assertAlmostEqual(0.9, new_sample['reflecting'])
+
+ def test_boundary_above_periodic(self):
+ new_sample = self.jump_proposal(self.sample_above)
+ self.assertAlmostEqual(-0.4, new_sample['periodic'])
+
+ def test_boundary_above_default(self):
+ new_sample = self.jump_proposal(self.sample_above)
+ self.assertAlmostEqual(0.9, new_sample['default'])
+
+ def test_boundary_below_reflecting(self):
+ new_sample = self.jump_proposal(self.sample_below)
+ self.assertAlmostEqual(-0.4, new_sample['reflecting'])
+
+ def test_boundary_below_periodic(self):
+ new_sample = self.jump_proposal(self.sample_below)
+ self.assertAlmostEqual(0.9, new_sample['periodic'])
+
+ def test_boundary_below_default(self):
+ new_sample = self.jump_proposal(self.sample_below)
+ self.assertAlmostEqual(-0.4, new_sample['default'])
+
+ def test_boundary_way_below_reflecting_case1(self):
+ new_sample = self.jump_proposal(self.sample_way_below_case1)
+ self.assertAlmostEqual(0.0, new_sample['reflecting'])
+
+ def test_boundary_way_below_reflecting_case2(self):
+ new_sample = self.jump_proposal(self.sample_way_below_case2)
+ self.assertAlmostEqual(-0.1, new_sample['reflecting'])
+
+ def test_boundary_way_below_periodic(self):
+ new_sample = self.jump_proposal(self.sample_way_below_case2)
+ self.assertAlmostEqual(-0.1, new_sample['periodic'])
+
+ def test_boundary_way_above_reflecting_case1(self):
+ new_sample = self.jump_proposal(self.sample_way_above_case1)
+ self.assertAlmostEqual(0.0, new_sample['reflecting'])
+
+ def test_boundary_way_above_reflecting_case2(self):
+ new_sample = self.jump_proposal(self.sample_way_above_case2)
+ self.assertAlmostEqual(0.1, new_sample['reflecting'])
+
+ def test_boundary_way_above_periodic(self):
+ new_sample = self.jump_proposal(self.sample_way_above_case2)
+ self.assertAlmostEqual(0.1, new_sample['periodic'])
+
+ def test_priors(self):
+ self.assertEqual(self.priors, self.jump_proposal.priors)
+
+
+class TestNormJump(unittest.TestCase):
+
+ def setUp(self):
+ self.priors = prior.PriorDict(dict(reflecting=prior.Uniform(minimum=-0.5, maximum=1, periodic_boundary=True),
+ periodic=prior.Uniform(minimum=-0.5, maximum=1, periodic_boundary=False),
+ default=prior.Uniform(minimum=-0.5, maximum=1)))
+ self.jump_proposal = proposal.NormJump(step_size=3.0, priors=self.priors)
+
+ def tearDown(self):
+ del self.priors
+ del self.jump_proposal
+
+ def test_step_size_init(self):
+ self.assertEqual(3.0, self.jump_proposal.step_size)
+
+ def test_set_step_size(self):
+ self.jump_proposal.step_size = 1.0
+ self.assertEqual(1.0, self.jump_proposal.step_size)
+
+ def test_jump_proposal_call(self):
+ with mock.patch("numpy.random.normal") as m:
+ m.return_value = 0.5
+ sample = proposal.Sample(dict(reflecting=0.0, periodic=0.0, default=0.0))
+ new_sample = self.jump_proposal(sample)
+ expected = proposal.Sample(dict(reflecting=0.5, periodic=0.5, default=0.5))
+ self.assertDictEqual(expected, new_sample)
+
+
+class TestEnsembleWalk(unittest.TestCase):
+
+ def setUp(self):
+ self.priors = prior.PriorDict(dict(reflecting=prior.Uniform(minimum=-0.5, maximum=1, periodic_boundary=False),
+ periodic=prior.Uniform(minimum=-0.5, maximum=1, periodic_boundary=True),
+ default=prior.Uniform(minimum=-0.5, maximum=1)))
+ self.jump_proposal = proposal.EnsembleWalk(random_number_generator=random.random,
+ n_points=4, priors=self.priors)
+
+ def tearDown(self):
+ del self.priors
+ del self.jump_proposal
+
+ def test_n_points_init(self):
+ self.assertEqual(4, self.jump_proposal.n_points)
+
+ def test_set_n_points(self):
+ self.jump_proposal.n_points = 3
+ self.assertEqual(3, self.jump_proposal.n_points)
+
+ def test_random_number_generator_init(self):
+ self.assertEqual(random.random, self.jump_proposal.random_number_generator)
+
+ def test_get_center_of_mass(self):
+ samples = [proposal.Sample(dict(reflecting=0.1*i, periodic=0.1*i, default=0.1*i)) for i in range(3)]
+ expected = proposal.Sample(dict(reflecting=0.1, periodic=0.1, default=0.1))
+ actual = self.jump_proposal.get_center_of_mass(samples)
+ for key in samples[0].keys():
+ self.assertAlmostEqual(expected[key], actual[key])
+
+ def test_jump_proposal_call(self):
+ with mock.patch('random.sample') as m:
+ self.jump_proposal.random_number_generator = lambda: 2
+ m.return_value = [proposal.Sample(dict(periodic=0.3, reflecting=0.3, default=0.3)),
+ proposal.Sample(dict(periodic=0.1, reflecting=0.1, default=0.1))]
+ sample = proposal.Sample(dict(periodic=0.1, reflecting=0.1, default=0.1))
+ new_sample = self.jump_proposal(sample, coordinates=None)
+ expected = proposal.Sample(dict(periodic=0.1, reflecting=0.1, default=0.1))
+ for key, value in new_sample.items():
+ self.assertAlmostEqual(expected[key], value)
+
+
+class TestEnsembleEnsembleStretch(unittest.TestCase):
+
+ def setUp(self):
+ self.priors = prior.PriorDict(dict(reflecting=prior.Uniform(minimum=-0.5, maximum=1, periodic_boundary=False),
+ periodic=prior.Uniform(minimum=-0.5, maximum=1, periodic_boundary=True),
+ default=prior.Uniform(minimum=-0.5, maximum=1)))
+ self.jump_proposal = proposal.EnsembleStretch(scale=3.0, priors=self.priors)
+
+ def tearDown(self):
+ del self.priors
+ del self.jump_proposal
+
+ def test_scale_init(self):
+ self.assertEqual(3.0, self.jump_proposal.scale)
+
+ def test_set_get_scale(self):
+ self.jump_proposal.scale = 5.0
+ self.assertEqual(5.0, self.jump_proposal.scale)
+
+ def test_jump_proposal_call(self):
+ with mock.patch('random.choice') as m:
+ with mock.patch('random.uniform') as n:
+ second_sample = proposal.Sample(dict(periodic=0.3, reflecting=0.3, default=0.3))
+ random_number = 0.5
+ m.return_value = second_sample
+ n.return_value = random_number
+ sample = proposal.Sample(dict(periodic=0.1, reflecting=0.1, default=0.1))
+ new_sample = self.jump_proposal(sample, coordinates=None)
+ coords = 0.3 - 0.2 * np.exp(random_number * np.log(self.jump_proposal.scale))
+ expected = proposal.Sample(dict(periodic=coords, reflecting=coords, default=coords))
+ for key, value in new_sample.items():
+ self.assertAlmostEqual(expected[key], value)
+
+ def test_log_j_after_call(self):
+ with mock.patch('random.uniform') as m1:
+ with mock.patch('numpy.log') as m2:
+ with mock.patch('numpy.exp') as m3:
+ m1.return_value = 1
+ m2.return_value = 1
+ m3.return_value = 1
+ coordinates = [proposal.Sample(dict(periodic=0.3, reflecting=0.3, default=0.3)),
+ proposal.Sample(dict(periodic=0.3, reflecting=0.3, default=0.3))]
+ sample = proposal.Sample(dict(periodic=0.2, reflecting=0.2, default=0.2))
+ self.jump_proposal(sample=sample,
+ coordinates=coordinates)
+ self.assertEqual(3, self.jump_proposal.log_j)
+
+
+class TestDifferentialEvolution(unittest.TestCase):
+
+ def setUp(self):
+ self.priors = prior.PriorDict(dict(reflecting=prior.Uniform(minimum=-0.5, maximum=1, periodic_boundary=False),
+ periodic=prior.Uniform(minimum=-0.5, maximum=1, periodic_boundary=True),
+ default=prior.Uniform(minimum=-0.5, maximum=1)))
+ self.jump_proposal = proposal.DifferentialEvolution(sigma=1e-3, mu=0.5, priors=self.priors)
+
+ def tearDown(self):
+ del self.priors
+ del self.jump_proposal
+
+ def test_mu_init(self):
+ self.assertEqual(0.5, self.jump_proposal.mu)
+
+ def test_set_get_mu(self):
+ self.jump_proposal.mu = 1
+ self.assertEqual(1, self.jump_proposal.mu)
+
+ def test_set_get_sigma(self):
+ self.jump_proposal.sigma = 2
+ self.assertEqual(2, self.jump_proposal.sigma)
+
+ def test_jump_proposal_call(self):
+ with mock.patch('random.sample') as m:
+ with mock.patch('random.gauss') as n:
+ m.return_value = proposal.Sample(dict(periodic=0.2, reflecting=0.2, default=0.2)),\
+ proposal.Sample(dict(periodic=0.3, reflecting=0.3, default=0.3))
+ n.return_value = 1
+ sample = proposal.Sample(dict(periodic=0.1, reflecting=0.1, default=0.1))
+ expected = proposal.Sample(dict(periodic=0.2, reflecting=0.2, default=0.2))
+ new_sample = self.jump_proposal(sample, coordinates=None)
+ for key, value in new_sample.items():
+ self.assertAlmostEqual(expected[key], value)
+
+
+class TestEnsembleEigenVector(unittest.TestCase):
+
+ def setUp(self):
+ self.priors = prior.PriorDict(dict(reflecting=prior.Uniform(minimum=-0.5, maximum=1, periodic_boundary=False),
+ periodic=prior.Uniform(minimum=-0.5, maximum=1, periodic_boundary=True),
+ default=prior.Uniform(minimum=-0.5, maximum=1)))
+ self.jump_proposal = proposal.EnsembleEigenVector(priors=self.priors)
+
+ def tearDown(self):
+ del self.priors
+ del self.jump_proposal
+
+ def test_init_eigen_values(self):
+ self.assertIsNone(self.jump_proposal.eigen_values)
+
+ def test_init_eigen_vectors(self):
+ self.assertIsNone(self.jump_proposal.eigen_vectors)
+
+ def test_init_covariance(self):
+ self.assertIsNone(self.jump_proposal.covariance)
+
+ def test_jump_proposal_update_eigenvectors_none(self):
+ self.assertIsNone(self.jump_proposal.update_eigenvectors(coordinates=None))
+
+ def test_jump_proposal_update_eigenvectors_1_d(self):
+ coordinates = [proposal.Sample(dict(periodic=0.3)), proposal.Sample(dict(periodic=0.1))]
+ with mock.patch('numpy.var') as m:
+ m.return_value = 1
+ self.jump_proposal.update_eigenvectors(coordinates)
+ self.assertTrue(np.equal(np.array([1]), self.jump_proposal.eigen_values))
+ self.assertTrue(np.equal(np.array([1]), self.jump_proposal.covariance))
+ self.assertTrue(np.equal(np.array([[1.]]), self.jump_proposal.eigen_vectors))
+
+ def test_jump_proposal_update_eigenvectors_n_d(self):
+ coordinates = [proposal.Sample(dict(periodic=0.3, reflecting=0.3, default=0.3)),
+ proposal.Sample(dict(periodic=0.1, reflecting=0.1, default=0.1))]
+ with mock.patch('numpy.cov') as m:
+ with mock.patch('numpy.linalg.eigh') as n:
+ m.side_effect = lambda x: x
+ n.return_value = 1, 2
+ self.jump_proposal.update_eigenvectors(coordinates)
+ self.assertTrue(np.array_equal(np.array([[0.3, 0.1], [0.3, 0.1], [0.3, 0.1]]), self.jump_proposal.covariance))
+ self.assertEqual(1, self.jump_proposal.eigen_values)
+ self.assertEqual(2, self.jump_proposal.eigen_vectors)
+
+ def test_jump_proposal_call(self):
+ self.jump_proposal.update_eigenvectors = lambda x: None
+ self.jump_proposal.eigen_values = np.array([1, np.nan, np.nan])
+ self.jump_proposal.eigen_vectors = np.array([[0.1, np.nan, np.nan],
+ [0.4, np.nan, np.nan],
+ [0.7, np.nan, np.nan]])
+ with mock.patch('random.randrange') as m:
+ with mock.patch('random.gauss') as n:
+ m.return_value = 0
+ n.return_value = 1
+ expected = proposal.Sample()
+ expected['periodic'] = 0.2
+ expected['reflecting'] = 0.5
+ expected['default'] = 0.8
+ sample = proposal.Sample()
+ sample['periodic'] = 0.1
+ sample['reflecting'] = 0.1
+ sample['default'] = 0.1
+ new_sample = self.jump_proposal(sample, coordinates=None)
+ for key, value in new_sample.items():
+ self.assertAlmostEqual(expected[key], value)
+
+
+class TestSkyLocationWanderJump(unittest.TestCase):
+
+ def setUp(self):
+ self.priors = prior.PriorDict(dict(ra=prior.Uniform(minimum=0.0, maximum=2*np.pi, periodic_boundary=True),
+ dec=prior.Uniform(minimum=0.0, maximum=np.pi, periodic_boundary=False)))
+ self.jump_proposal = bilby.gw.sampler.proposal.SkyLocationWanderJump(priors=self.priors)
+
+ def tearDown(self):
+ del self.priors
+ del self.jump_proposal
+
+ def test_jump_proposal_call_without_inverse_temperature(self):
+ with mock.patch('random.gauss') as m:
+ m.return_value = 1
+ sample = proposal.Sample(dict(ra=0.2, dec=-0.5))
+ expected = proposal.Sample(dict(ra=1.2, dec=0.5))
+ new_sample = self.jump_proposal(sample)
+ for key, value in new_sample.items():
+ self.assertAlmostEqual(expected[key], value)
+ m.assert_called_with(0, 1.0 / 2 / np.pi)
+
+ def test_jump_proposal_call_with_inverse_temperature(self):
+ with mock.patch('random.gauss') as m:
+ m.return_value = 1
+ sample = proposal.Sample(dict(ra=0.2, dec=-0.5))
+ expected = proposal.Sample(dict(ra=1.2, dec=0.5))
+ new_sample = self.jump_proposal(sample, inverse_temperature=2.0)
+ for key, value in new_sample.items():
+ self.assertAlmostEqual(expected[key], value)
+ m.assert_called_with(0, np.sqrt(1 / 2.0) / 2 / np.pi)
+
+
+class TestCorrelatedPolarisationPhaseJump(unittest.TestCase):
+
+ def setUp(self):
+ self.priors = prior.PriorDict(dict(phase=prior.Uniform(minimum=0.0, maximum=2*np.pi),
+ psi=prior.Uniform(minimum=0.0, maximum=np.pi)))
+ self.jump_proposal = bilby.gw.sampler.proposal.CorrelatedPolarisationPhaseJump(priors=self.priors)
+
+ def tearDown(self):
+ del self.priors
+ del self.jump_proposal
+
+ def test_jump_proposal_call_case_1(self):
+ with mock.patch('random.random') as m:
+ m.return_value = 0.3
+ sample = proposal.Sample(dict(phase=0.2, psi=0.5))
+ alpha = 3.0 * np.pi * 0.3
+ beta = 0.3
+ expected = proposal.Sample(dict(phase=0.5*(alpha-beta), psi=0.5*(alpha+beta)))
+ self.assertEqual(expected, self.jump_proposal(sample, coordinates=None))
+
+ def test_jump_proposal_call_case_2(self):
+ with mock.patch('random.random') as m:
+ m.return_value = 0.7
+ sample = proposal.Sample(dict(phase=0.2, psi=0.5))
+ alpha = 0.7
+ beta = 3.0 * np.pi * 0.7 - 2 * np.pi
+ expected = proposal.Sample(dict(phase=0.5*(alpha-beta), psi=0.5*(alpha+beta)))
+ self.assertEqual(expected, self.jump_proposal(sample))
+
+
+class TestPolarisationPhaseJump(unittest.TestCase):
+
+ def setUp(self):
+ self.priors = prior.PriorDict(dict(phase=prior.Uniform(minimum=0.0, maximum=2*np.pi),
+ psi=prior.Uniform(minimum=0.0, maximum=np.pi)))
+ self.jump_proposal = bilby.gw.sampler.proposal.PolarisationPhaseJump(priors=self.priors)
+
+ def tearDown(self):
+ del self.priors
+ del self.jump_proposal
+
+ def test_jump_proposal_call(self):
+ sample = proposal.Sample(dict(phase=0.2, psi=0.5))
+ expected = proposal.Sample(dict(phase=0.2+np.pi, psi=0.5+np.pi/2))
+ self.assertEqual(expected, self.jump_proposal(sample))
+
+
+class TestDrawFlatPrior(unittest.TestCase):
+
+ def setUp(self):
+ self.priors = prior.PriorDict(dict(phase=prior.Uniform(minimum=0.0, maximum=2*np.pi),
+ psi=prior.Cosine(minimum=0.0, maximum=np.pi)))
+ self.jump_proposal = proposal.DrawFlatPrior(priors=self.priors)
+
+ def tearDown(self):
+ del self.priors
+ del self.jump_proposal
+
+ def test_jump_proposal_call(self):
+ with mock.patch('bilby.core.prior.Uniform.sample') as m:
+ m.return_value = 0.3
+ sample = proposal.Sample(dict(phase=0.2, psi=0.5))
+ expected = proposal.Sample(dict(phase=0.3, psi=0.3))
+ self.assertEqual(expected, self.jump_proposal(sample))
diff --git a/test/sampler_test.py b/test/sampler_test.py
index 2500fab13..314605470 100644
--- a/test/sampler_test.py
+++ b/test/sampler_test.py
@@ -113,13 +113,13 @@ class TestCPNest(unittest.TestCase):
def test_default_kwargs(self):
expected = dict(verbose=1, nthreads=1, nlive=500, maxmcmc=1000,
seed=None, poolsize=100, nhamiltonian=0, resume=True,
- output='outdir/cpnest_label/')
+ output='outdir/cpnest_label/', proposals=None)
self.assertDictEqual(expected, self.sampler.kwargs)
def test_translate_kwargs(self):
expected = dict(verbose=1, nthreads=1, nlive=250, maxmcmc=1000,
seed=None, poolsize=100, nhamiltonian=0, resume=True,
- output='outdir/cpnest_label/')
+ output='outdir/cpnest_label/', proposals=None)
for equiv in bilby.core.sampler.base_sampler.NestedSampler.npoints_equiv_kwargs:
new_kwargs = self.sampler.kwargs.copy()
del new_kwargs['nlive']
--
GitLab