diff --git a/CHANGELOG.md b/CHANGELOG.md
index cdd2320735c7b47ca6c7f080a80a5b4a18f2d8eb..f34d8285fb71c74ad06d1b1e5c2e46915aefbcf7 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -3,14 +3,31 @@
## Unreleased
### Added
--
+-
### Changed
+- Make calibration work, maybe with interp1d
-
### Removed
-
+## [0.4.4] 2019-04-03
+
+### Added
+- Infrastucture for custom jump proposals (cpnest-only)
+- Evidence uncertainty estimate to cpnest
+
+### Changed
+- Bug fix to close figures after creation
+- Improved the frequency-mask to entirely remove values outside the mask rather
+ than simply set them to zero
+- Fix problem with Prior prob and ln_prob if passing multiple samples
+- Improved cpnest prior sampling
+
+### Removed
+-
+
## [0.4.3] 2019-03-21
### Added
diff --git a/bilby/core/prior.py b/bilby/core/prior.py
index aa3693e1642726dc7ca1c3610ae50bb1a497ce05..3108124e57e3737b43c4786cc5a482dd83ba7d74 100644
--- a/bilby/core/prior.py
+++ b/bilby/core/prior.py
@@ -268,12 +268,20 @@ class PriorDict(OrderedDict):
"""
prob = np.product([self[key].prob(sample[key])
for key in sample], **kwargs)
- if prob == 0:
- return 0
- elif self.evaluate_constraints(sample):
+
+ if np.all(prob == 0.):
return prob
else:
- return 0
+ if isinstance(prob, float):
+ if self.evaluate_constraints(sample):
+ return prob
+ else:
+ return 0.
+ else:
+ constrained_prob = np.zeros_like(prob)
+ keep = np.array(self.evaluate_constraints(sample), dtype=bool)
+ constrained_prob[keep] = prob[keep]
+ return constrained_prob
def ln_prob(self, sample, axis=None):
"""
@@ -293,12 +301,20 @@ class PriorDict(OrderedDict):
"""
ln_prob = np.sum([self[key].ln_prob(sample[key])
for key in sample], axis=axis)
- if np.isinf(ln_prob):
- return ln_prob
- elif self.evaluate_constraints(sample):
+
+ if np.all(np.isinf(ln_prob)):
return ln_prob
else:
- return -np.inf
+ if isinstance(ln_prob, float):
+ if self.evaluate_constraints(sample):
+ return ln_prob
+ else:
+ return -np.inf
+ else:
+ constrained_ln_prob = -np.inf * np.ones_like(ln_prob)
+ keep = np.array(self.evaluate_constraints(sample), dtype=bool)
+ constrained_ln_prob[keep] = ln_prob[keep]
+ return constrained_ln_prob
def rescale(self, keys, theta):
"""Rescale samples from unit cube to prior
@@ -392,7 +408,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
@@ -407,13 +423,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.
@@ -617,6 +635,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():
@@ -694,7 +722,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
@@ -711,9 +739,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):
@@ -774,13 +805,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
@@ -795,9 +827,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)
@@ -835,7 +870,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
@@ -850,9 +885,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))
@@ -860,7 +897,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
@@ -880,9 +917,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):
"""
@@ -940,7 +980,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
@@ -955,9 +995,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):
"""
@@ -986,7 +1028,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
@@ -1001,9 +1043,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):
"""
@@ -1031,7 +1075,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
@@ -1046,8 +1090,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
@@ -1079,7 +1125,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
@@ -1094,15 +1140,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
@@ -1123,9 +1171,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
@@ -1161,14 +1211,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
@@ -1187,14 +1237,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
@@ -1207,14 +1259,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
@@ -1227,14 +1281,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
@@ -1251,10 +1307,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")
@@ -1290,7 +1347,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
@@ -1307,14 +1364,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
@@ -1327,10 +1385,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):
@@ -1362,7 +1421,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
@@ -1382,8 +1441,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")
@@ -1422,7 +1483,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
@@ -1446,7 +1507,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")
@@ -1456,7 +1518,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):
@@ -1550,7 +1612,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
@@ -1567,8 +1629,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")
@@ -1605,7 +1669,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
@@ -1622,8 +1686,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")
@@ -1660,7 +1726,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
@@ -1677,13 +1743,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
@@ -1700,9 +1768,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")
@@ -1740,7 +1810,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
@@ -1755,13 +1825,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):
@@ -1775,7 +1847,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
@@ -1794,6 +1866,8 @@ class Interped(Prior):
See superclass
unit: str
See superclass
+ periodic_boundary: bool
+ See superclass
Attributes
-------
@@ -1812,7 +1886,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):
@@ -1905,7 +1980,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
@@ -1922,6 +1997,8 @@ class FromFile(Interped):
See superclass
unit: str
See superclass
+ periodic_boundary: bool
+ See superclass
Attributes
-------
@@ -1933,8 +2010,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/result.py b/bilby/core/result.py
index 82fb07792c4ee844a2e6ec5cb67ddbc6c271e4be..a4630475fad32a85d70dd3d0becc946d948c8e87 100644
--- a/bilby/core/result.py
+++ b/bilby/core/result.py
@@ -412,12 +412,17 @@ class Result(object):
default=False
outdir: str, optional
Path to the outdir. Default is the one stored in the result object.
- extension: str, optional {json, hdf5}
- Determines the method to use to store the data
+ extension: str, optional {json, hdf5, True}
+ Determines the method to use to store the data (if True defaults
+ to json)
gzip: bool, optional
If true, and outputing to a json file, this will gzip the resulting
file and add '.gz' to the file extension.
"""
+
+ if extension is True:
+ extension = "json"
+
outdir = self._safe_outdir_creation(outdir, self.save_to_file)
file_name = result_file_name(outdir, self.label, extension, gzip)
@@ -978,6 +983,17 @@ class Result(object):
fig.savefig(filename, dpi=dpi)
plt.close(fig)
+ @staticmethod
+ def _add_prior_fixed_values_to_posterior(posterior, priors):
+ if priors is None:
+ return posterior
+ for key in priors:
+ if isinstance(priors[key], DeltaFunction):
+ posterior[key] = priors[key].peak
+ elif isinstance(priors[key], float):
+ posterior[key] = priors[key]
+ return posterior
+
def samples_to_posterior(self, likelihood=None, priors=None,
conversion_function=None):
"""
@@ -1000,11 +1016,8 @@ class Result(object):
except ValueError:
data_frame = pd.DataFrame(
self.samples, columns=self.search_parameter_keys)
- for key in priors:
- if isinstance(priors[key], DeltaFunction):
- data_frame[key] = priors[key].peak
- elif isinstance(priors[key], float):
- data_frame[key] = priors[key]
+ data_frame = self._add_prior_fixed_values_to_posterior(
+ data_frame, priors)
data_frame['log_likelihood'] = getattr(
self, 'log_likelihood_evaluations', np.nan)
if self.log_prior_evaluations is None:
diff --git a/bilby/core/sampler/__init__.py b/bilby/core/sampler/__init__.py
index c8a5d2307e808ee6dab2d1ff4529edc7b7a04576..7feed4e8049f570e7c03c541f57e3c3cc4adfd4b 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
@@ -167,6 +168,10 @@ def run_sampler(likelihood, priors=None, label='label', outdir='outdir',
else:
result = sampler.run_sampler()
+ # Initial save of the sampler in case of failure in post-processing
+ if save:
+ result.save_to_file(extension=save, gzip=gzip)
+
end_time = datetime.datetime.now()
result.sampling_time = (end_time - start_time).total_seconds()
logger.info('Sampling time: {}'.format(end_time - start_time))
@@ -188,12 +193,11 @@ def run_sampler(likelihood, priors=None, label='label', outdir='outdir',
result.samples_to_posterior(likelihood=likelihood, priors=result.priors,
conversion_function=conversion_function)
- if save == 'hdf5':
- result.save_to_file(extension='hdf5')
- logger.info("Results saved to {}/".format(outdir))
- elif save:
- result.save_to_file(gzip=gzip)
- logger.info("Results saved to {}/".format(outdir))
+
+ if save:
+ # The overwrite here ensures we overwrite the initially stored data
+ result.save_to_file(overwrite=True, extension=save, gzip=gzip)
+
if plot:
result.plot_corner()
logger.info("Summary of results:\n{}".format(result))
diff --git a/bilby/core/sampler/cpnest.py b/bilby/core/sampler/cpnest.py
index 0f1cd8780827756b25e4bbce9bbea064f9c2301d..6acd33689f2c6dbcf6204387f52d2737bd6ef9b7 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:
@@ -49,14 +52,17 @@ class Cpnest(NestedSampler):
def run_sampler(self):
from cpnest import model as cpmodel, CPNest
+ from cpnest.parameter import LivePoint
class Model(cpmodel.Model):
""" A wrapper class to pass our log_likelihood into cpnest """
- def __init__(self, names, bounds):
+ def __init__(self, names, priors):
self.names = names
- self.bounds = bounds
- self._check_bounds()
+ self.priors = priors
+ self.bounds = [
+ [self.priors[key].minimum, self.priors[key].maximum]
+ for key in self.names]
@staticmethod
def log_likelihood(x, **kwargs):
@@ -68,16 +74,26 @@ class Cpnest(NestedSampler):
theta = [x[n] for n in self.search_parameter_keys]
return self.log_prior(theta)
- def _check_bounds(self):
- for bound in self.bounds:
- if not all(np.isfinite(bound)):
- raise ValueError(
- 'CPNest requires priors to have finite bounds.')
-
- bounds = [[self.priors[key].minimum, self.priors[key].maximum]
- for key in self.search_parameter_keys]
- model = Model(self.search_parameter_keys, bounds)
- out = CPNest(model, **self.kwargs)
+ def new_point(self):
+ """Draw a point from the prior"""
+ point = LivePoint(
+ self.names, [self.priors[name].sample()
+ for name in self.names])
+ return point
+
+ self._resolve_proposal_functions()
+ model = Model(self.search_parameter_keys, self.priors)
+ 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:
@@ -87,7 +103,7 @@ class Cpnest(NestedSampler):
self.result.posterior.rename(columns=dict(
logL='log_likelihood', logPrior='log_prior'), inplace=True)
self.result.log_evidence = out.NS.state.logZ
- self.result.log_evidence_err = np.nan
+ self.result.log_evidence_err = np.sqrt(out.NS.state.info / out.NS.state.nlive)
return self.result
def _verify_kwargs_against_default_kwargs(self):
@@ -101,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/dynesty.py b/bilby/core/sampler/dynesty.py
index a51ce406d75e2061af47270484a0f68a2d872199..e2787895f0cf220896b5e6c7aa232f046b33058e 100644
--- a/bilby/core/sampler/dynesty.py
+++ b/bilby/core/sampler/dynesty.py
@@ -189,7 +189,10 @@ class Dynesty(NestedSampler):
if self.kwargs["verbose"]:
print("")
- # self.result.sampler_output = out
+ dynesty_result = "{}/{}_dynesty.pickle".format(self.outdir, self.label)
+ with open(dynesty_result, 'wb') as file:
+ pickle.dump(out, file)
+
weights = np.exp(out['logwt'] - out['logz'][-1])
nested_samples = DataFrame(
out.samples, columns=self.search_parameter_keys)
@@ -307,8 +310,8 @@ class Dynesty(NestedSampler):
return False
def write_current_state_and_exit(self, signum=None, frame=None):
+ logger.warning("Run terminated with signal {}".format(signum))
self.write_current_state()
- logger.warning("Run terminated")
sys.exit()
def write_current_state(self):
@@ -328,6 +331,7 @@ class Dynesty(NestedSampler):
NestedSampler to write to disk.
"""
check_directory_exists_and_if_not_mkdir(self.outdir)
+ logger.info("Writing checkpoint file {}".format(self.resume_file))
current_state = dict(
unit_cube_samples=self.sampler.saved_u,
diff --git a/bilby/core/sampler/emcee.py b/bilby/core/sampler/emcee.py
index 3117a6c7a2cdd11266bce26fda6e9b06feb9a36d..6b4a0e732aa3274dbe5fee8ce4f211138964f1d5 100644
--- a/bilby/core/sampler/emcee.py
+++ b/bilby/core/sampler/emcee.py
@@ -1,10 +1,15 @@
from __future__ import absolute_import, print_function
+from collections import namedtuple
import os
+import signal
+from shutil import copyfile
+import sys
import numpy as np
from pandas import DataFrame
from distutils.version import LooseVersion
+import dill as pickle
from ..utils import (
logger, get_progress_bar, check_directory_exists_and_if_not_mkdir)
@@ -41,11 +46,11 @@ class Emcee(MCMCSampler):
"""
- default_kwargs = dict(nwalkers=500, a=2, args=[], kwargs={},
- postargs=None, pool=None, live_dangerously=False,
- runtime_sortingfn=None, lnprob0=None, rstate0=None,
- blobs0=None, iterations=100, thin=1, storechain=True,
- mh_proposal=None)
+ default_kwargs = dict(
+ nwalkers=500, a=2, args=[], kwargs={}, postargs=None, pool=None,
+ live_dangerously=False, runtime_sortingfn=None, lnprob0=None,
+ rstate0=None, blobs0=None, iterations=100, thin=1, storechain=True,
+ mh_proposal=None)
def __init__(self, likelihood, priors, outdir='outdir', label='label',
use_ratio=False, plot=False, skip_import_verification=False,
@@ -65,7 +70,9 @@ class Emcee(MCMCSampler):
self.nburn = nburn
self.burn_in_fraction = burn_in_fraction
self.burn_in_act = burn_in_act
- self._old_chain = None
+
+ signal.signal(signal.SIGTERM, self.checkpoint_and_exit)
+ signal.signal(signal.SIGINT, self.checkpoint_and_exit)
def _translate_kwargs(self, kwargs):
if 'nwalkers' not in kwargs:
@@ -136,6 +143,14 @@ class Emcee(MCMCSampler):
return init_kwargs
+ def lnpostfn(self, theta):
+ log_prior = self.log_prior(theta)
+ if np.isinf(log_prior):
+ return -np.inf, [np.nan, np.nan]
+ else:
+ log_likelihood = self.log_likelihood(theta)
+ return log_likelihood + log_prior, [log_likelihood, log_prior]
+
@property
def nburn(self):
if type(self.__nburn) in [float, int]:
@@ -166,71 +181,141 @@ class Emcee(MCMCSampler):
def nsteps(self, nsteps):
self.kwargs['iterations'] = nsteps
- def __getstate__(self):
- # In order to be picklable with dill, we need to discard the pool
- # object before trying.
- d = self.__dict__
- d["_Sampler__kwargs"]["pool"] = None
- return d
+ @property
+ def stored_chain(self):
+ """ Read the stored zero-temperature chain data in from disk """
+ return np.genfromtxt(self.checkpoint_info.chain_file, names=True)
- def run_sampler(self):
- import emcee
- tqdm = get_progress_bar()
- sampler = emcee.EnsembleSampler(**self.sampler_init_kwargs)
- out_dir = os.path.join(self.outdir, 'emcee_{}'.format(self.label))
- out_file = os.path.join(out_dir, 'chain.dat')
+ @property
+ def stored_samples(self):
+ """ Returns the samples stored on disk """
+ return self.stored_chain[self.search_parameter_keys]
- if self.resume:
- self.load_old_chain(out_file)
- else:
- self._set_pos0()
+ @property
+ def stored_loglike(self):
+ """ Returns the log-likelihood stored on disk """
+ return self.stored_chain['log_l']
+ @property
+ def stored_logprior(self):
+ """ Returns the log-prior stored on disk """
+ return self.stored_chain['log_p']
+
+ def _init_chain_file(self):
+ with open(self.checkpoint_info.chain_file, "w+") as ff:
+ ff.write('walker\t{}\tlog_l\tlog_p\n'.format(
+ '\t'.join(self.search_parameter_keys)))
+
+ @property
+ def checkpoint_info(self):
+ """ Defines various things related to checkpointing and storing data
+
+ Returns
+ -------
+ checkpoint_info: named_tuple
+ An object with attributes `sampler_file`, `chain_file`, and
+ `chain_template`. The first two give paths to where the sampler and
+ chain data is stored, the last a formatted-str-template with which
+ to write the chain data to disk
+
+ """
+ out_dir = os.path.join(
+ self.outdir, '{}_{}'.format(self.__class__.__name__.lower(),
+ self.label))
check_directory_exists_and_if_not_mkdir(out_dir)
- if not os.path.isfile(out_file):
- with open(out_file, "w") as ff:
- ff.write('walker\t{}\tlog_l'.format(
- '\t'.join(self.search_parameter_keys)))
- template =\
+
+ chain_file = os.path.join(out_dir, 'chain.dat')
+ sampler_file = os.path.join(out_dir, 'sampler.pickle')
+ chain_template =\
'{:d}' + '\t{:.9e}' * (len(self.search_parameter_keys) + 2) + '\n'
- for sample in tqdm(sampler.sample(**self.sampler_function_kwargs),
- total=self.nsteps):
- if self.prerelease:
- points = np.hstack([sample.coords, sample.blobs])
- else:
- points = np.hstack([sample[0], np.array(sample[3])])
- with open(out_file, "a") as ff:
- for ii, point in enumerate(points):
- ff.write(template.format(ii, *point))
+ CheckpointInfo = namedtuple(
+ 'CheckpointInfo', ['sampler_file', 'chain_file', 'chain_template'])
- self.result.sampler_output = np.nan
- blobs_flat = np.array(sampler.blobs).reshape((-1, 2))
- log_likelihoods, log_priors = blobs_flat.T
- if self._old_chain is not None:
- chain = np.vstack([self._old_chain[:, :-2],
- sampler.chain.reshape((-1, self.ndim))])
- log_ls = np.hstack([self._old_chain[:, -2], log_likelihoods])
- log_ps = np.hstack([self._old_chain[:, -1], log_priors])
- self.nsteps = chain.shape[0] // self.nwalkers
+ checkpoint_info = CheckpointInfo(
+ sampler_file=sampler_file, chain_file=chain_file,
+ chain_template=chain_template)
+
+ return checkpoint_info
+
+ @property
+ def sampler_chain(self):
+ nsteps = self._previous_iterations
+ return self.sampler.chain[:, :nsteps, :]
+
+ def checkpoint(self):
+ """ Writes a pickle file of the sampler to disk using dill """
+ logger.info("Checkpointing sampler to file {}"
+ .format(self.checkpoint_info.sampler_file))
+ with open(self.checkpoint_info.sampler_file, 'wb') as f:
+ # Overwrites the stored sampler chain with one that is truncated
+ # to only the completed steps
+ self.sampler._chain = self.sampler_chain
+ pickle.dump(self._sampler, f)
+
+ def checkpoint_and_exit(self, signum, frame):
+ logger.info("Recieved signal {}".format(signum))
+ self.checkpoint()
+ sys.exit()
+
+ def _initialise_sampler(self):
+ import emcee
+ self._sampler = emcee.EnsembleSampler(**self.sampler_init_kwargs)
+ self._init_chain_file()
+
+ @property
+ def sampler(self):
+ """ Returns the ptemcee sampler object
+
+ If, alrady initialized, returns the stored _sampler value. Otherwise,
+ first checks if there is a pickle file from which to load. If there is
+ not, then initialize the sampler and set the initial random draw
+
+ """
+ if hasattr(self, '_sampler'):
+ pass
+ elif self.resume and os.path.isfile(self.checkpoint_info.sampler_file):
+ logger.info("Resuming run from checkpoint file {}"
+ .format(self.checkpoint_info.sampler_file))
+ with open(self.checkpoint_info.sampler_file, 'rb') as f:
+ self._sampler = pickle.load(f)
+ self._set_pos0_for_resume()
else:
- chain = sampler.chain.reshape((-1, self.ndim))
- log_ls = log_likelihoods
- log_ps = log_priors
- self.calculate_autocorrelation(chain)
- self.print_nburn_logging_info()
- self.result.nburn = self.nburn
- n_samples = self.nwalkers * self.nburn
- if self.result.nburn > self.nsteps:
- raise SamplerError(
- "The run has finished, but the chain is not burned in: "
- "`nburn < nsteps`. Try increasing the number of steps.")
- self.result.samples = chain[n_samples:, :]
- self.result.log_likelihood_evaluations = log_ls[n_samples:]
- self.result.log_prior_evaluations = log_ps[n_samples:]
- self.result.walkers = sampler.chain
- self.result.log_evidence = np.nan
- self.result.log_evidence_err = np.nan
- return self.result
+ self._initialise_sampler()
+ self._set_pos0()
+ return self._sampler
+
+ def write_chains_to_file(self, sample):
+ chain_file = self.checkpoint_info.chain_file
+ temp_chain_file = chain_file + '.temp'
+ if os.path.isfile(chain_file):
+ copyfile(chain_file, temp_chain_file)
+
+ if self.prerelease:
+ points = np.hstack([sample.coords, sample.blobs])
+ else:
+ points = np.hstack([sample[0], np.array(sample[3])])
+ with open(temp_chain_file, "a") as ff:
+ for ii, point in enumerate(points):
+ ff.write(self.checkpoint_info.chain_template.format(ii, *point))
+ os.rename(temp_chain_file, chain_file)
+
+ @property
+ def _previous_iterations(self):
+ """ Returns the number of iterations that the sampler has saved
+
+ This is used when loading in a sampler from a pickle file to figure out
+ how much of the run has already been completed
+ """
+ return len(self.sampler.blobs)
+
+ def _draw_pos0_from_prior(self):
+ return np.array(
+ [self.get_random_draw_from_prior() for _ in range(self.nwalkers)])
+
+ @property
+ def _pos0_shape(self):
+ return (self.nwalkers, self.ndim)
def _set_pos0(self):
if self.pos0 is not None:
@@ -238,9 +323,9 @@ class Emcee(MCMCSampler):
if isinstance(self.pos0, DataFrame):
self.pos0 = self.pos0[self.search_parameter_keys].values
elif type(self.pos0) in (list, np.ndarray):
- self.pos0 = np.squeeze(self.kwargs['pos0'])
+ self.pos0 = np.squeeze(self.pos0)
- if self.pos0.shape != (self.nwalkers, self.ndim):
+ if self.pos0.shape != self._pos0_shape:
raise ValueError(
'Input pos0 should be of shape ndim, nwalkers')
logger.debug("Checking input pos0")
@@ -248,27 +333,44 @@ class Emcee(MCMCSampler):
self.check_draw(draw)
else:
logger.debug("Generating initial walker positions from prior")
- self.pos0 = [self.get_random_draw_from_prior()
- for _ in range(self.nwalkers)]
-
- def load_old_chain(self, file_name=None):
- if file_name is None:
- out_dir = os.path.join(self.outdir, 'emcee_{}'.format(self.label))
- file_name = os.path.join(out_dir, 'chain.dat')
- if os.path.isfile(file_name):
- old_chain = np.genfromtxt(file_name, skip_header=1)
- self.pos0 = [np.squeeze(old_chain[-(self.nwalkers - ii), 1:-2])
- for ii in range(self.nwalkers)]
- self._old_chain = old_chain[:-self.nwalkers + 1, 1:]
- logger.info('Resuming from {}'.format(os.path.abspath(file_name)))
- else:
- logger.warning('Failed to resume. {} not found.'.format(file_name))
- self._set_pos0()
+ self.pos0 = self._draw_pos0_from_prior()
- def lnpostfn(self, theta):
- log_prior = self.log_prior(theta)
- if np.isinf(log_prior):
- return -np.inf, [np.nan, np.nan]
- else:
- log_likelihood = self.log_likelihood(theta)
- return log_likelihood + log_prior, [log_likelihood, log_prior]
+ def _set_pos0_for_resume(self):
+ self.pos0 = self.sampler.chain[:, -1, :]
+
+ def run_sampler(self):
+ tqdm = get_progress_bar()
+ sampler_function_kwargs = self.sampler_function_kwargs
+ iterations = sampler_function_kwargs.pop('iterations')
+ iterations -= self._previous_iterations
+
+ sampler_function_kwargs['p0'] = self.pos0
+
+ # main iteration loop
+ for sample in tqdm(
+ self.sampler.sample(iterations=iterations, **sampler_function_kwargs),
+ total=iterations):
+ self.write_chains_to_file(sample)
+ self.checkpoint()
+
+ self.result.sampler_output = np.nan
+ blobs_flat = np.array(self.sampler.blobs).reshape((-1, 2))
+ log_likelihoods, log_priors = blobs_flat.T
+ chain = self.sampler.chain.reshape((-1, self.ndim))
+ log_ls = log_likelihoods
+ log_ps = log_priors
+ self.calculate_autocorrelation(chain)
+ self.print_nburn_logging_info()
+ self.result.nburn = self.nburn
+ n_samples = self.nwalkers * self.nburn
+ if self.result.nburn > self.nsteps:
+ raise SamplerError(
+ "The run has finished, but the chain is not burned in: "
+ "`nburn < nsteps`. Try increasing the number of steps.")
+ self.result.samples = chain[n_samples:, :]
+ self.result.log_likelihood_evaluations = log_ls[n_samples:]
+ self.result.log_prior_evaluations = log_ps[n_samples:]
+ self.result.walkers = self.sampler.chain
+ self.result.log_evidence = np.nan
+ self.result.log_evidence_err = np.nan
+ return self.result
diff --git a/bilby/core/sampler/proposal.py b/bilby/core/sampler/proposal.py
new file mode 100644
index 0000000000000000000000000000000000000000..a928292b32b3293d593b0080e49f316c849c3a20
--- /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/core/sampler/ptemcee.py b/bilby/core/sampler/ptemcee.py
index 8ea4165cd4a2f600700f18c81d46bfe762a07d91..31c5609a6a72ac338cf36981108bd37d4f73154b 100644
--- a/bilby/core/sampler/ptemcee.py
+++ b/bilby/core/sampler/ptemcee.py
@@ -1,8 +1,11 @@
from __future__ import absolute_import, division, print_function
+import os
+from shutil import copyfile
+
import numpy as np
-from ..utils import get_progress_bar
+from ..utils import logger, get_progress_bar
from . import Emcee
from .base_sampler import SamplerError
@@ -27,23 +30,22 @@ class Ptemcee(Emcee):
The number of temperatures used by ptemcee
"""
- default_kwargs = dict(ntemps=2, nwalkers=500,
- Tmax=None, betas=None,
- threads=1, pool=None, a=2.0,
- loglargs=[], logpargs=[],
- loglkwargs={}, logpkwargs={},
- adaptation_lag=10000, adaptation_time=100,
- random=None, iterations=100, thin=1,
- storechain=True, adapt=True,
- swap_ratios=False,
- )
-
- def __init__(self, likelihood, priors, outdir='outdir', label='label', use_ratio=False, plot=False,
- skip_import_verification=False, nburn=None, burn_in_fraction=0.25,
- burn_in_act=3, **kwargs):
- Emcee.__init__(self, likelihood=likelihood, priors=priors, outdir=outdir, label=label,
- use_ratio=use_ratio, plot=plot, skip_import_verification=skip_import_verification,
- nburn=nburn, burn_in_fraction=burn_in_fraction, burn_in_act=burn_in_act, **kwargs)
+ default_kwargs = dict(
+ ntemps=2, nwalkers=500, Tmax=None, betas=None, threads=1, pool=None,
+ a=2.0, loglargs=[], logpargs=[], loglkwargs={}, logpkwargs={},
+ adaptation_lag=10000, adaptation_time=100, random=None, iterations=100,
+ thin=1, storechain=True, adapt=True, swap_ratios=False)
+
+ def __init__(self, likelihood, priors, outdir='outdir', label='label',
+ use_ratio=False, plot=False, skip_import_verification=False,
+ nburn=None, burn_in_fraction=0.25, burn_in_act=3, resume=True,
+ **kwargs):
+ Emcee.__init__(
+ self, likelihood=likelihood, priors=priors, outdir=outdir,
+ label=label, use_ratio=use_ratio, plot=plot,
+ skip_import_verification=skip_import_verification,
+ nburn=nburn, burn_in_fraction=burn_in_fraction,
+ burn_in_act=burn_in_act, resume=resume, **kwargs)
@property
def sampler_function_kwargs(self):
@@ -56,42 +58,97 @@ class Ptemcee(Emcee):
for key, value in self.kwargs.items()
if key not in self.sampler_function_kwargs}
- def run_sampler(self):
+ @property
+ def ntemps(self):
+ return self.kwargs['ntemps']
+
+ @property
+ def sampler_chain(self):
+ nsteps = self._previous_iterations
+ return self.sampler.chain[:, :, :nsteps, :]
+
+ def _initialise_sampler(self):
import ptemcee
+ self._sampler = ptemcee.Sampler(
+ dim=self.ndim, logl=self.log_likelihood, logp=self.log_prior,
+ **self.sampler_init_kwargs)
+ self._init_chain_file()
+
+ def print_tswap_acceptance_fraction(self):
+ logger.info("Sampler per-chain tswap acceptance fraction = {}".format(
+ self.sampler.tswap_acceptance_fraction))
+
+ def write_chains_to_file(self, pos, loglike, logpost):
+ chain_file = self.checkpoint_info.chain_file
+ temp_chain_file = chain_file + '.temp'
+ if os.path.isfile(chain_file):
+ copyfile(chain_file, temp_chain_file)
+
+ with open(temp_chain_file, "a") as ff:
+ loglike = np.squeeze(loglike[0, :])
+ logprior = np.squeeze(logpost[0, :]) - loglike
+ for ii, (point, logl, logp) in enumerate(zip(pos[0, :, :], loglike, logprior)):
+ line = np.concatenate((point, [logl, logp]))
+ ff.write(self.checkpoint_info.chain_template.format(ii, *line))
+ os.rename(temp_chain_file, chain_file)
+
+ @property
+ def _previous_iterations(self):
+ """ Returns the number of iterations that the sampler has saved
+
+ This is used when loading in a sampler from a pickle file to figure out
+ how much of the run has already been completed
+ """
+ return self.sampler.time
+
+ def _draw_pos0_from_prior(self):
+ # for ptemcee, the pos0 has the shape ntemps, nwalkers, ndim
+ return [[self.get_random_draw_from_prior()
+ for _ in range(self.nwalkers)]
+ for _ in range(self.kwargs['ntemps'])]
+
+ @property
+ def _pos0_shape(self):
+ return (self.ntemps, self.nwalkers, self.ndim)
+
+ def _set_pos0_for_resume(self):
+ self.pos0 = None
+
+ def run_sampler(self):
tqdm = get_progress_bar()
- sampler = ptemcee.Sampler(dim=self.ndim, logl=self.log_likelihood,
- logp=self.log_prior, **self.sampler_init_kwargs)
- self.pos0 = [[self.get_random_draw_from_prior()
- for _ in range(self.nwalkers)]
- for _ in range(self.kwargs['ntemps'])]
-
- log_likelihood_evaluations = []
- log_prior_evaluations = []
+ sampler_function_kwargs = self.sampler_function_kwargs
+ iterations = sampler_function_kwargs.pop('iterations')
+ iterations -= self._previous_iterations
+
+ # main iteration loop
for pos, logpost, loglike in tqdm(
- sampler.sample(self.pos0, **self.sampler_function_kwargs),
- total=self.nsteps):
- log_likelihood_evaluations.append(loglike)
- log_prior_evaluations.append(logpost - loglike)
- pass
+ self.sampler.sample(self.pos0, iterations=iterations,
+ **sampler_function_kwargs),
+ total=iterations):
+ self.write_chains_to_file(pos, loglike, logpost)
+ self.checkpoint()
- self.calculate_autocorrelation(sampler.chain.reshape((-1, self.ndim)))
+ self.calculate_autocorrelation(self.sampler.chain.reshape((-1, self.ndim)))
self.result.sampler_output = np.nan
self.print_nburn_logging_info()
+ self.print_tswap_acceptance_fraction()
+
self.result.nburn = self.nburn
if self.result.nburn > self.nsteps:
raise SamplerError(
"The run has finished, but the chain is not burned in: "
"`nburn < nsteps`. Try increasing the number of steps.")
- self.result.samples = sampler.chain[0, :, self.nburn:, :].reshape(
+
+ self.result.samples = self.sampler.chain[0, :, self.nburn:, :].reshape(
(-1, self.ndim))
- self.result.log_likelihood_evaluations = np.array(
- log_likelihood_evaluations)[self.nburn:, 0, :].reshape((-1))
- self.result.log_prior_evaluations = np.array(
- log_prior_evaluations)[self.nburn:, 0, :].reshape((-1))
- self.result.betas = sampler.betas
+ self.result.walkers = self.sampler.chain[0, :, :, :]
+
+ n_samples = self.nwalkers * self.nburn
+ self.result.log_likelihood_evaluations = self.stored_loglike[n_samples:]
+ self.result.log_prior_evaluations = self.stored_logprior[n_samples:]
+ self.result.betas = self.sampler.betas
self.result.log_evidence, self.result.log_evidence_err =\
- sampler.log_evidence_estimate(
- sampler.loglikelihood, self.nburn / self.nsteps)
- self.result.walkers = sampler.chain[0, :, :, :]
+ self.sampler.log_evidence_estimate(
+ self.sampler.loglikelihood, self.nburn / self.nsteps)
return self.result
diff --git a/bilby/core/series.py b/bilby/core/series.py
index 980c6c11afe529e4fc1cfbef3b4f26e69a512a6c..0c9b2efaa2e830c8096315aae89a0bc42229c1f0 100644
--- a/bilby/core/series.py
+++ b/bilby/core/series.py
@@ -35,7 +35,7 @@ class CoupledTimeAndFrequencySeries(object):
-------
array_like: The frequency array
"""
- if self._frequency_array_updated is False:
+ if not self._frequency_array_updated:
if self.sampling_frequency and self.duration:
self._frequency_array = utils.create_frequency_series(
sampling_frequency=self.sampling_frequency,
@@ -45,6 +45,7 @@ class CoupledTimeAndFrequencySeries(object):
'legitimate sampling_frequency ({}) or duration ({})'
.format(self.sampling_frequency, self.duration))
+ self._frequency_array_updated = True
return self._frequency_array
@frequency_array.setter
@@ -63,7 +64,7 @@ class CoupledTimeAndFrequencySeries(object):
array_like: The time array
"""
- if self._time_array_updated is False:
+ if not self._time_array_updated:
if self.sampling_frequency and self.duration:
self._time_array = utils.create_time_series(
sampling_frequency=self.sampling_frequency,
diff --git a/bilby/gw/calibration.py b/bilby/gw/calibration.py
index d95a12f002a391188ec5efcc3a9dbc44594b0e39..cf4fc7d9392e35f369edb31d2379b5430b65a058 100644
--- a/bilby/gw/calibration.py
+++ b/bilby/gw/calibration.py
@@ -2,7 +2,7 @@
"""
import numpy as np
-from scipy.interpolate import UnivariateSpline
+from scipy.interpolate import interp1d
class Recalibrate(object):
@@ -83,11 +83,12 @@ class CubicSpline(Recalibrate):
self.n_points = n_points
self.minimum_frequency = minimum_frequency
self.maximum_frequency = maximum_frequency
- self.__spline_points = np.logspace(np.log10(minimum_frequency), np.log10(maximum_frequency), n_points)
+ self._log_spline_points = np.linspace(
+ np.log10(minimum_frequency), np.log10(maximum_frequency), n_points)
@property
- def spline_points(self):
- return self.__spline_points
+ def log_spline_points(self):
+ return self._log_spline_points
def __repr__(self):
return self.__class__.__name__ + '(prefix=\'{}\', minimum_frequency={}, maximum_frequency={}, n_points={})'\
@@ -112,13 +113,17 @@ class CubicSpline(Recalibrate):
The factor to multiply the strain by.
"""
self.set_calibration_parameters(**params)
- amplitude_parameters = [self.params['amplitude_{}'.format(ii)] for ii in range(self.n_points)]
- amplitude_spline = UnivariateSpline(self.spline_points, amplitude_parameters)
- delta_amplitude = amplitude_spline(frequency_array)
-
- phase_parameters = [self.params['phase_{}'.format(ii)] for ii in range(self.n_points)]
- phase_spline = UnivariateSpline(self.spline_points, phase_parameters)
- delta_phase = phase_spline(frequency_array)
+ amplitude_parameters = [self.params['amplitude_{}'.format(ii)]
+ for ii in range(self.n_points)]
+ delta_amplitude = interp1d(
+ self.log_spline_points, amplitude_parameters, kind='cubic',
+ bounds_error=False, fill_value=0)(np.log10(frequency_array))
+
+ phase_parameters = [
+ self.params['phase_{}'.format(ii)] for ii in range(self.n_points)]
+ delta_phase = interp1d(
+ self.log_spline_points, phase_parameters, kind='cubic',
+ bounds_error=False, fill_value=0)(np.log10(frequency_array))
calibration_factor = (1 + delta_amplitude) * (2 + 1j * delta_phase) / (2 - 1j * delta_phase)
diff --git a/bilby/gw/likelihood.py b/bilby/gw/likelihood.py
index 7e913a50364a14d4ce3eda51c5fbe247600be43f..bbcb75ab770a891644436721e7320f12f3351c4a 100644
--- a/bilby/gw/likelihood.py
+++ b/bilby/gw/likelihood.py
@@ -70,6 +70,12 @@ class GravitationalWaveTransient(likelihood.Likelihood):
"""
+ _CalculatedSNRs = namedtuple('CalculatedSNRs',
+ ['d_inner_h',
+ 'optimal_snr_squared',
+ 'complex_matched_filter_snr',
+ 'd_inner_h_squared_tc_array'])
+
def __init__(self, interferometers, waveform_generator,
time_marginalization=False, distance_marginalization=False,
phase_marginalization=False, priors=None,
@@ -160,8 +166,6 @@ class GravitationalWaveTransient(likelihood.Likelihood):
d_inner_h = interferometer.inner_product(signal=signal)
optimal_snr_squared = interferometer.optimal_snr_squared(signal=signal)
complex_matched_filter_snr = d_inner_h / (optimal_snr_squared**0.5)
- CalculatedSNRs = namedtuple(
- 'CalculatedSNRs', ['d_inner_h', 'optimal_snr_squared', 'complex_matched_filter_snr', 'd_inner_h_squared_tc_array'])
if self.time_marginalization:
d_inner_h_squared_tc_array =\
@@ -172,7 +176,7 @@ class GravitationalWaveTransient(likelihood.Likelihood):
else:
d_inner_h_squared_tc_array = None
- return CalculatedSNRs(
+ return self._CalculatedSNRs(
d_inner_h=d_inner_h, optimal_snr_squared=optimal_snr_squared,
complex_matched_filter_snr=complex_matched_filter_snr,
d_inner_h_squared_tc_array=d_inner_h_squared_tc_array)
@@ -334,7 +338,7 @@ class GravitationalWaveTransient(likelihood.Likelihood):
dmax = self._distance_array[-1]
n = len(self._distance_array)
self._lookup_table_filename = (
- '.distance_marginalization_lookup_dmin{}_dmax{}_n{}.npz'
+ '.distance_marginalization_lookup.npz'
.format(dmin, dmax, n))
return self._lookup_table_filename
@@ -609,12 +613,9 @@ class ROQGravitationalWaveTransient(GravitationalWaveTransient):
h_plus_quadratic = f_plus * signal['quadratic']['plus']
h_cross_quadratic = f_cross * signal['quadratic']['cross']
- CalculatedSNRs = namedtuple(
- 'CalculatedSNRs', ['d_inner_h', 'optimal_snr_squared', 'complex_matched_filter_snr', 'd_inner_h_squared_tc_array'])
-
indices, in_bounds = self._closest_time_indices(ifo_time, self.time_samples)
if not in_bounds:
- return CalculatedSNRs(
+ return self._CalculatedSNRs(
d_inner_h=np.nan_to_num(-np.inf), optimal_snr_squared=0,
complex_matched_filter_snr=np.nan_to_num(-np.inf),
d_inner_h_squared_tc_array=None)
@@ -634,7 +635,7 @@ class ROQGravitationalWaveTransient(GravitationalWaveTransient):
complex_matched_filter_snr = d_inner_h / (optimal_snr_squared**0.5)
d_inner_h_squared_tc_array = None
- return CalculatedSNRs(
+ return self._CalculatedSNRs(
d_inner_h=d_inner_h, optimal_snr_squared=optimal_snr_squared,
complex_matched_filter_snr=complex_matched_filter_snr,
d_inner_h_squared_tc_array=d_inner_h_squared_tc_array)
diff --git a/bilby/gw/prior.py b/bilby/gw/prior.py
index 58790e41cff09904f4accfe656f3a01f7ead8e70..c751303d2c9270c8a48697db62236f70e95f3eb8 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 410fc927682662918f3d229640b564b672e4a86e..9641b77d906a058481bc7fb76512b8f9738540e8 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 ca49a506b37a1fd7d8e855029932e8c437335a4c..8be6dcd1eeea32d50eecd9598cd98c5be30b7b71 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 e79bd5baa754eb09aeb2a9de41b47adf82b9452c..6df97a1dd817948fa86d88410445bf89f11bb4fa 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 eff1a2f2a4923b520cf3e242ee1bc1545cb6c64f..09acbdd14d8df03ea842872efab4de803c14e611 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 0000000000000000000000000000000000000000..25ad2ca38576b83835854b025b36f3a73b6a1cb8
--- /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 0000000000000000000000000000000000000000..328fd98a89f051b634f99602a5cfba50b7a68dc6
--- /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/bilby/gw/utils.py b/bilby/gw/utils.py
index b0959727f202d07e6e7e69b31030df485ee481d9..3bb38eb7a896dd6cac4ca92aa310f7f61b11a570 100644
--- a/bilby/gw/utils.py
+++ b/bilby/gw/utils.py
@@ -443,7 +443,7 @@ def read_frame_file(file_name, start_time, end_time, channel=None, buffer_time=0
return None
-def get_gracedb(gracedb, outdir, duration, calibration, detectors, query_types=None):
+def get_gracedb(gracedb, outdir, duration, calibration, detectors, query_types=None, server=None):
candidate = gracedb_to_json(gracedb, outdir)
trigger_time = candidate['gpstime']
gps_start_time = trigger_time - duration
@@ -453,7 +453,7 @@ def get_gracedb(gracedb, outdir, duration, calibration, detectors, query_types=N
for i, det in enumerate(detectors):
output_cache_file = gw_data_find(
det, gps_start_time, duration, calibration,
- outdir=outdir, query_type=query_types[i])
+ outdir=outdir, query_type=query_types[i], server=server)
cache_files.append(output_cache_file)
return candidate, cache_files
@@ -476,7 +476,7 @@ def gracedb_to_json(gracedb, outdir=None):
logger.info('Initialise client and attempt to download')
try:
client = GraceDb()
- except FileNotFoundError:
+ except IOError:
raise ValueError(
'Failed to authenticate with gracedb: check your X509 '
'certificate is accessible and valid')
@@ -499,7 +499,7 @@ def gracedb_to_json(gracedb, outdir=None):
def gw_data_find(observatory, gps_start_time, duration, calibration,
- outdir='.', query_type=None):
+ outdir='.', query_type=None, server=None):
""" Builds a gw_data_find call and process output
Parameters
@@ -542,6 +542,13 @@ def gw_data_find(observatory, gps_start_time, duration, calibration,
output_cache_file = os.path.join(outdir, cache_file)
gps_end_time = gps_start_time + duration
+ if server is None:
+ server = os.environ.get("LIGO_DATAFIND_SERVER")
+ if server is None:
+ logger.warning('No data_find server found, defaulting to CIT server. '
+ 'To run on other clusters, pass the output of '
+ '`echo $LIGO_DATAFIND_SERVER`')
+ server = 'ldr.ldas.cit:80'
cl_list = ['gw_data_find']
cl_list.append('--observatory {}'.format(observatory_code))
@@ -549,6 +556,7 @@ def gw_data_find(observatory, gps_start_time, duration, calibration,
cl_list.append('--gps-end-time {}'.format(int(np.ceil(gps_end_time))))
cl_list.append('--type {}'.format(query_type))
cl_list.append('--output {}'.format(output_cache_file))
+ cl_list.append('--server {}'.format(server))
cl_list.append('--url-type file')
cl_list.append('--lal-cache')
cl = ' '.join(cl_list)
diff --git a/bilby/gw/waveform_generator.py b/bilby/gw/waveform_generator.py
index 00e24bde621c5b989ba3daca27d4f4ebf4ed2267..98d06914695af20473745e88d3916a0367832e2b 100644
--- a/bilby/gw/waveform_generator.py
+++ b/bilby/gw/waveform_generator.py
@@ -62,6 +62,7 @@ class WaveformGenerator(object):
self.waveform_arguments = dict()
if isinstance(parameters, dict):
self.parameters = parameters
+ self._cache = dict(parameters=None, waveform=None)
def __repr__(self):
if self.frequency_domain_source_model is not None:
@@ -147,6 +148,8 @@ class WaveformGenerator(object):
transformed_model_data_points, parameters):
if parameters is not None:
self.parameters = parameters
+ if self.parameters == self._cache['parameters']:
+ return self._cache['waveform']
if model is not None:
model_strain = self._strain_from_model(model_data_points, model)
elif transformed_model is not None:
@@ -154,6 +157,8 @@ class WaveformGenerator(object):
transformation_function)
else:
raise RuntimeError("No source model given")
+ self._cache['waveform'] = model_strain
+ self._cache['parameters'] = self.parameters.copy()
return model_strain
def _strain_from_model(self, model_data_points, model):
diff --git a/docs/conversion.txt b/docs/conversion.txt
index f130df063d75bbe1dd32888f737c874c48a209a0..d15671aa98f13e7eb2f54fb054ada4383d6d2281 100644
--- a/docs/conversion.txt
+++ b/docs/conversion.txt
@@ -7,5 +7,11 @@ E.g., sampling in chirp mass and mass ratio can be much more efficient than samp
We have many functions to do this.
+These are used in multiple places:
+- `PriorDict`s have a `conversion_function`, for the GW PriorDicts, these are from this module.
+- `WaveformGenerator`s can use a `parameter_conversion`, again these are from this module.
+- A `conversion_function` can be passed to `run_sampler`, this is done as a post-processing step.
+For CBCs either `generate_all_bbh_parameters` or `generate_all_bns_parameters` can be used.
+
.. automodule:: bilby.gw.conversion
:members:
\ No newline at end of file
diff --git a/docs/prior.txt b/docs/prior.txt
index 867a843516e4c4db7dff685cf8953e431f1f46e6..fe655557ff3ea611c5d4b7c7bfbc8b93f099e2fa 100644
--- a/docs/prior.txt
+++ b/docs/prior.txt
@@ -83,8 +83,9 @@ note that this list is incomplete.
:members:
:special-members:
+---------------------------
Multivariate Gaussian prior
-===========================
+---------------------------
We provide a prior class for correlated parameters in the form of a
`multivariate Gaussian distribution <https://en.wikipedia.org/wiki/Multivariate_normal_distribution>`_.
@@ -130,11 +131,85 @@ The detailed API information for the distribution and prior classes are below:
.. autoclass:: bilby.core.prior.MultivariateGaussian
:members:
+-----------------------
Defining your own prior
-=======================
+-----------------------
You can define your own by subclassing the :code:`bilby.prior.Prior` class.
.. autoclass:: bilby.core.prior.Prior
:members:
:special-members:
+
+-----------------
+Prior Constraints
+-----------------
+
+Added v. 0.4.3
+
+This allows cuts to be specified in the prior space.
+
+You can provide the `PriorDict` a `conversion_function` and a set of `Constraint` priors to remove parts of the prior space.
+
+*Note*: after doing this the prior probability will not be normalised.
+
+Simple example
+==============
+
+Sample from uniform distributions in two parameters x and y with the condition x >= y.
+
+First thing: define a function which generates z=x-y from x and y.
+
+.. code:: python
+
+ def convert_x_y_to_z(parameters):
+ """
+ Function to convert between sampled parameters and constraint parameter.
+
+ Parameters
+ ----------
+ parameters: dict
+ Dictionary containing sampled parameter values, 'x', 'y'.
+
+ Returns
+ -------
+ dict: Dictionary with constraint parameter 'z' added.
+ """
+ parameters['z'] = parameters['x'] - parameters['y']
+ return parameters
+
+Create our prior:
+
+.. code:: python
+
+ from bilby.core.prior import PriorDict, Uniform, Constraint
+
+ priors = PriorDict(conversion_function=convert_x_y_to_z)
+ priors['x'] = Uniform(minimum=0, maximum=10)
+ priors['y'] = Uniform(minimum=0, maximum=10)
+ priors['z'] = Constraint(minimum=0, maximum=10)
+
+Sample from this distribution and plot the samples.
+
+.. code:: python
+
+ import matplotlib.pyplot as plt
+
+ samples = priors.sample(1000000)
+ plt.hist2d(samples['x'], samples['y'], bins=100, cmap='Blues')
+ plt.xlabel('$x$')
+ plt.ylabel('$y$')
+ plt.tight_layout()
+ plt.show()
+ plt.close()
+
+------
+
+Gravitational wave priors
+=========================
+
+We provide default conversion functions for the BBH and BNS PriorDicts.
+
+For BBHs this generates all the BBH mass parameters so constraints can be placed on any mass parameters.
+
+For BNSs it also generates the tidal deformability parameters.
diff --git a/examples/injection_examples/custom_proposal_example.py b/examples/injection_examples/custom_proposal_example.py
new file mode 100644
index 0000000000000000000000000000000000000000..8480dd361a141b7407d7ab618028ad542d737cd1
--- /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/requirements.txt b/requirements.txt
index bb184b6a62790c07f5b232056a2d89b0291f398e..de58a6b16f36ea4614961dcdb0b54cce46167f8e 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -6,3 +6,4 @@ matplotlib>=2.0
scipy>=0.16
pandas
mock
+dill
diff --git a/setup.py b/setup.py
index 81575551bd6cfc9ae8a6682328739361c9141d8d..43e8190855ca0dbde139eb17cdd1f49817461822 100644
--- a/setup.py
+++ b/setup.py
@@ -57,7 +57,7 @@ def readfile(filename):
return filecontents
-VERSION = '0.4.3'
+VERSION = '0.4.4'
version_file = write_version_file(VERSION)
long_description = get_long_description()
@@ -79,6 +79,7 @@ setup(name='bilby',
'future',
'dynesty',
'corner',
+ 'dill',
'numpy>=1.9',
'matplotlib>=2.0',
'pandas',
diff --git a/test/gw_prior_test.py b/test/gw_prior_test.py
index 893d4f3de0ed3391df17b17603e9588d71d2f80e..485341dbfa4fc82895ca48b9afd7d438236cf8b1 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 e79bd5baa754eb09aeb2a9de41b47adf82b9452c..b9abc7a45e36e7f42e9a768128d9e5723643b82d 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 eff1a2f2a4923b520cf3e242ee1bc1545cb6c64f..09acbdd14d8df03ea842872efab4de803c14e611 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 ce3ab12041c2343744e44e1a5b6f1b7eafc10411..4734d80e2845e51d1b4877dda62f75579e915031 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):
@@ -163,7 +185,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'),
bilby.core.prior.MultivariateGaussian(mvg=mvg, name='testa', unit='unit'),
bilby.core.prior.MultivariateGaussian(mvg=mvg, name='testb', unit='unit'),
@@ -489,8 +510,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,
@@ -530,43 +552,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)
@@ -576,8 +601,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),
@@ -586,33 +613,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 0000000000000000000000000000000000000000..dc86444c455d019b8f44886d0ba6cd9c7be839e7
--- /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/result_test.py b/test/result_test.py
index 7648b9708948d9c293e4ae0a3df268d46b74d37b..a18a44ac6bb088253bfe014f99b01de87f33bddc 100644
--- a/test/result_test.py
+++ b/test/result_test.py
@@ -266,6 +266,17 @@ class TestResult(unittest.TestCase):
df = pd.read_csv(filename)
self.assertTrue(np.allclose(self.result.posterior.values, df.values))
+ def test_samples_to_posterior_simple(self):
+ self.result.posterior = None
+ x = [1, 2, 3]
+ y = [4, 6, 8]
+ self.result.samples = np.array([x, y]).T
+ self.result.samples_to_posterior()
+ self.assertTrue(all(self.result.posterior['x'] == x))
+ self.assertTrue(all(self.result.posterior['y'] == y))
+ self.assertTrue(np.all(
+ None == self.result.posterior.log_likelihood.values))
+
def test_samples_to_posterior(self):
self.result.posterior = None
x = [1, 2, 3]
diff --git a/test/sampler_test.py b/test/sampler_test.py
index f3063b4628e9d76b3808a557258ea8a3321cc7eb..6aa1d138b6df1275e81ca477f9290c973071f8ed 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']
diff --git a/test/waveform_generator_test.py b/test/waveform_generator_test.py
index 8eba859ac03f0a1d73ade2a74f0785aff1162953..15aa65e467f2ecccb379effab951417904bed6b7 100644
--- a/test/waveform_generator_test.py
+++ b/test/waveform_generator_test.py
@@ -263,6 +263,19 @@ class TestFrequencyDomainStrainMethod(unittest.TestCase):
self.assertListEqual(sorted(self.waveform_generator.parameters.keys()),
sorted(['amplitude', 'mu', 'sigma', 'ra', 'dec', 'geocent_time', 'psi']))
+ def test_caching_with_parameters(self):
+ original_waveform = self.waveform_generator.frequency_domain_strain(
+ parameters=self.simulation_parameters)
+ new_waveform = self.waveform_generator.frequency_domain_strain(
+ parameters=self.simulation_parameters)
+ self.assertDictEqual(original_waveform, new_waveform)
+
+ def test_caching_without_parameters(self):
+ original_waveform = self.waveform_generator.frequency_domain_strain(
+ parameters=self.simulation_parameters)
+ new_waveform = self.waveform_generator.frequency_domain_strain()
+ self.assertDictEqual(original_waveform, new_waveform)
+
class TestTimeDomainStrainMethod(unittest.TestCase):
@@ -354,6 +367,20 @@ class TestTimeDomainStrainMethod(unittest.TestCase):
self.assertListEqual(sorted(self.waveform_generator.parameters.keys()),
sorted(['amplitude', 'mu', 'sigma', 'ra', 'dec', 'geocent_time', 'psi']))
+ def test_caching_with_parameters(self):
+ original_waveform = self.waveform_generator.time_domain_strain(
+ parameters=self.simulation_parameters)
+ new_waveform = self.waveform_generator.time_domain_strain(
+ parameters=self.simulation_parameters)
+ self.assertDictEqual(original_waveform, new_waveform)
+
+ def test_caching_without_parameters(self):
+ original_waveform = self.waveform_generator.time_domain_strain(
+ parameters=self.simulation_parameters)
+ new_waveform = self.waveform_generator.time_domain_strain()
+ self.assertDictEqual(original_waveform, new_waveform)
+
+
if __name__ == '__main__':
unittest.main()