Skip to content
Snippets Groups Projects
Commit 02d4d60e authored by Bruce Edelman's avatar Bruce Edelman
Browse files

Revert "generalized the JointPrior object structure from Matthew Pitkin's... 

Revert "generalized the JointPrior object structure from Matthew Pitkin's MutlivariateGaussian prior formalism. TODO: add in the joint MapPrior for HEALPix priors"

This reverts commit 85cae594.
parent 39d52b56
No related branches found
No related tags found
1 merge request!638Add kombine to dockerfile pip installs
Hide whitespace changes
Inline Side-by-side
bilby/core/prior.py
+ 321
− 372
View file @ 02d4d60e
......@@ -2520,280 +2520,7 @@ class FermiDirac(Prior):
return lnp
class JointPriorDist(object):
def __init__(self, names, bounds=None):
"""
A class defining JointPriorDist that will be overwritten with child
classes defining the joint prior distribtuions between given parameters,
Parameters
----------
names: list
A list of the parameter names in the JointPriorDist. The
listed parameters must have the same order that they appear in
the lists of statistical parameters that may be passed in child class
bounds: list
A list of bounds on each parameter. The defaults are for bounds at
+/- infinity.
"""
if not isinstance(names, list):
self.names = [names]
else:
self.names = names
self.num_vars = len(self.names)
# set the bounds for each parameter
if isinstance(bounds, list):
if len(bounds) != len(self):
raise ValueError("Wrong number of parameter bounds")
# check bounds
for bound in bounds:
if isinstance(bounds, (list, tuple, np.ndarray)):
if len(bound) != 2:
raise ValueError("Bounds must contain an upper and "
"lower value.")
else:
if bound[1] <= bound[0]:
raise ValueError("Bounds are not properly set")
else:
raise TypeError("Bound must be a list")
logger.warning("If using bounded ranges on the multivariate "
"Gaussian this will lead to biased posteriors "
"for nested sampling routines that require "
"a prior transform.")
else:
bounds = [(-np.inf, np.inf) for _ in self.names]
self.bounds = {name: val for name, val in zip(self.names, bounds)}
self._current_sample = {} # initialise empty sample
self._uncorrelated = None
self._current_lnprob = None
# a dictionary of the parameters as requested by the prior
self.requested_parameters = OrderedDict()
self.reset_request()
# a dictionary of the rescaled parameters
self.rescale_parameters = OrderedDict()
self.reset_rescale()
# a list of sampled parameters
self.reset_sampled()
def reset_sampled(self):
self.sampled_parameters = []
self.current_sample = {}
def filled_request(self):
"""
Check if all requested parameters have been filled.
"""
return not np.any([val is None for val in
self.requested_parameters.values()])
def reset_request(self):
"""
Reset the requested parameters to None.
"""
for name in self.names:
self.requested_parameters[name] = None
def filled_rescale(self):
"""
Check is all the rescaled parameters have been filled.
"""
return not np.any([val is None for val in
self.rescale_parameters.values()])
def reset_rescale(self):
"""
Reset the rescaled parameters to None.
"""
for name in self.names:
self.rescale_parameters[name] = None
def _get_instantiation_dict(self):
subclass_args = infer_args_from_method(self.__init__)
property_names = [p for p in dir(self.__class__)
if isinstance(getattr(self.__class__, p), property)]
dict_with_properties = self.__dict__.copy()
for key in property_names:
dict_with_properties[key] = getattr(self, key)
instantiation_dict = OrderedDict()
for key in subclass_args:
if isinstance(dict_with_properties[key], list):
value = np.asarray(dict_with_properties[key]).tolist()
else:
value = dict_with_properties[key]
instantiation_dict[key] = value
return instantiation_dict
def __len__(self):
return len(self.names)
def __repr__(self):
"""Overrides the special method __repr__.
Returns a representation of this instance that resembles how it is instantiated.
Works correctly for all child classes
Returns
-------
str: A string representation of this instance
"""
dist_name = self.__class__.__name__
instantiation_dict = self._get_instantiation_dict()
args = ', '.join(['{}={}'.format(key, repr(instantiation_dict[key]))
for key in instantiation_dict])
return "{}({})".format(dist_name, args)
def prob(self, samp):
"""
Get the probability of a sample. For bounded priors the
probability will not be properly normalised.
"""
return np.exp(self.ln_prob(samp))
def _check_samp(self, value):
samp = np.asarray(value)
if len(samp.shape) == 1:
samp = samp.reshape(1, self.num_vars)
if len(samp.shape) != 2:
raise ValueError("Array is the wrong shape")
elif samp.shape[1] != self.num_vars:
raise ValueError("Array is the wrong shape")
# check sample(s) is within bounds
outbounds = np.ones(samp.shape[0], dtype=np.bool)
for s, bound in zip(samp.T, self.bounds.values()):
outbounds = (s < bound[0]) | (s > bound[1])
if np.any(outbounds):
break
return samp, outbounds
def ln_prob(self, value):
"""
Get the log-probability of a sample. For bounded priors the
probability will not be properly normalised.
Parameters
----------
value: array_like
A 1d vector of the sample, or 2d array of sample values with shape
NxM, where N is the number of samples and M is the number of
parameters.
"""
samp, outbounds = self._check_samp(value)
lnprob = -np.inf * np.ones(samp.shape[0])
lnprob = self._ln_prob(samp, lnprob, outbounds)
if samp.shape[0] == 1:
return lnprob[0]
else:
return lnprob
def _ln_prob(self, samp, lnprob, outbounds):
'''
CHILD CLASS OVERWRITES THIS METHOD AND FILLS IN THIS PART OF ln_prob METHOD to samplethe lnprob for
the value of this sample
'''
return lnprob
def sample(self, size=1, **kwargs):
"""
Draw, and set, a sample from the Dist
Parameters
----------
size: int
number of samples to generate, defualts to 1
"""
if size is None:
size = 1
# samples drawn from unit variance uncorrelated multivariate Gaussian
samps = self._draw_samp(size=size, **kwargs)
for i, name in enumerate(self.names):
if size == 1:
self.current_sample[name] = samps[:, i].flatten()[0]
else:
self.current_sample[name] = samps[:, i].flatten()
def _draw_samp(self, size, **kwargs):
"""
Draw, and set, a sample from the joint dist (needs to be ovewritten by child class)
Parameters
----------
size: int
number of samples to generate, defualts to 1
"""
samps = np.zeros((size, len(self)))
"""
Here is where the subclass where overwrite sampling method
"""
return samps
def rescale(self, value, **kwargs):
"""
Rescale from a unit hypercube to multivariate Gaussian. Note that no
bounds are applied in the rescale function.
Parameters
----------
value: array
A 1d vector sample (one for each parameter) drawn from a uniform
distribution between 0 and 1, or a 2d NxM array of samples where
N is the number of samples and M is the number of parameters.
mode: int
Specify which mode to sample from. If not set then a mode is
chosen randomly based on its weight.
Returns
-------
array:
An vector sample drawn from the multivariate Gaussian
distribution.
"""
# pick a mode (with a probability given by their weights)
samp = np.asarray(value)
if len(samp.shape) == 1:
samp = samp.reshape(1, self.num_vars)
if len(samp.shape) != 2:
raise ValueError("Array is the wrong shape")
elif samp.shape[1] != self.num_vars:
raise ValueError("Array is the wrong shape")
samp = self._rescale(samp, **kwargs)
return np.squeeze(samp)
def _rescale(self, samp, **kwargs):
'''
needs to be overwritten
:param samp:
:param kwargs:
:return:
'''
return samp
class MultivariateGaussianDist(JointPriorDist):
class MultivariateGaussianDist(object):
def __init__(self, names, nmodes=1, mus=None, sigmas=None, corrcoefs=None,
covs=None, weights=None, bounds=None):
......@@ -2841,7 +2568,40 @@ class MultivariateGaussianDist(JointPriorDist):
A list of bounds on each parameter. The defaults are for bounds at
+/- infinity.
"""
super(MultivariateGaussianDist, self).__init__(names=names, bounds=bounds)
if not isinstance(names, list):
self.names = [names]
else:
self.names = names
self.num_vars = len(self.names) # the number of parameters
# set the bounds for each parameter
if isinstance(bounds, list):
if len(bounds) != len(self):
raise ValueError("Wrong number of parameter bounds")
# check bounds
for bound in bounds:
if isinstance(bounds, (list, tuple, np.ndarray)):
if len(bound) != 2:
raise ValueError("Bounds must contain an upper and "
"lower value.")
else:
if bound[1] <= bound[0]:
raise ValueError("Bounds are not properly set")
else:
raise TypeError("Bound must be a list")
logger.warning("If using bounded ranges on the multivariate "
"Gaussian this will lead to biased posteriors "
"for nested sampling routines that require "
"a prior transform.")
else:
bounds = [(-np.inf, np.inf) for _ in self.names]
# set bounds as dictionary
self.bounds = {name: val for name, val in zip(self.names, bounds)}
self.mus = []
self.covs = []
......@@ -2898,23 +2658,70 @@ class MultivariateGaussianDist(JointPriorDist):
if len(weights) != 1:
raise ValueError("Wrong number of weights given")
for val in [mus, sigmas, covs, corrcoefs, weights]:
if val is not None and not isinstance(val, list):
raise TypeError("Value must be a list")
else:
if val is not None and len(val) != nmodes:
raise ValueError("Wrong number of modes given")
for val in [mus, sigmas, covs, corrcoefs, weights]:
if val is not None and not isinstance(val, list):
raise TypeError("Value must be a list")
else:
if val is not None and len(val) != nmodes:
raise ValueError("Wrong number of modes given")
# add the modes
self.nmodes = 0
for i in range(nmodes):
mu = mus[i] if mus is not None else None
sigma = sigmas[i] if sigmas is not None else None
corrcoef = corrcoefs[i] if corrcoefs is not None else None
cov = covs[i] if covs is not None else None
weight = weights[i] if weights is not None else 1.
self.add_mode(mu, sigma, corrcoef, cov, weight)
# a dictionary of the parameters as requested by the prior
self.requested_parameters = OrderedDict()
self.reset_request()
# a dictionary of the rescaled parameters
self.rescale_parameters = OrderedDict()
self.reset_rescale()
# a list of sampled parameters
self.reset_sampled()
def reset_sampled(self):
self.sampled_parameters = []
self.current_sample = {}
def filled_request(self):
"""
Check if all requested parameters have been filled.
"""
# add the modes
self.nmodes = 0
for i in range(nmodes):
mu = mus[i] if mus is not None else None
sigma = sigmas[i] if sigmas is not None else None
corrcoef = corrcoefs[i] if corrcoefs is not None else None
cov = covs[i] if covs is not None else None
weight = weights[i] if weights is not None else 1.
return not np.any([val is None for val in
self.requested_parameters.values()])
self.add_mode(mu, sigma, corrcoef, cov, weight)
def reset_request(self):
"""
Reset the requested parameters to None.
"""
for name in self.names:
self.requested_parameters[name] = None
def filled_rescale(self):
"""
Check is all the rescaled parameters have been filled.
"""
return not np.any([val is None for val in
self.rescale_parameters.values()])
def reset_rescale(self):
"""
Reset the rescaled parameters to None.
"""
for name in self.names:
self.rescale_parameters[name] = None
def add_mode(self, mus=None, sigmas=None, corrcoef=None, cov=None,
weight=1.):
......@@ -3023,37 +2830,69 @@ class MultivariateGaussianDist(JointPriorDist):
self.mvn.append(scipy.stats.multivariate_normal(mean=self.mus[-1],
cov=self.covs[-1]))
def _rescale(self, samp, **kwargs):
try:
mode = kwargs['mode']
except KeyError:
mode = None
def rescale(self, value, mode=None):
"""
Rescale from a unit hypercube to multivariate Gaussian. Note that no
bounds are applied in the rescale function.
Parameters
----------
value: array
A 1d vector sample (one for each parameter) drawn from a uniform
distribution between 0 and 1, or a 2d NxM array of samples where
N is the number of samples and M is the number of parameters.
mode: int
Specify which mode to sample from. If not set then a mode is
chosen randomly based on its weight.
Returns
-------
array:
An vector sample drawn from the multivariate Gaussian
distribution.
"""
# pick a mode (with a probability given by their weights)
if mode is None:
if self.nmodes == 1:
mode = 0
else:
mode = np.argwhere(self.cumweights - np.random.rand() > 0)[0][0]
samp = np.asarray(value)
if len(samp.shape) == 1:
samp = samp.reshape(1, self.num_vars)
if len(samp.shape) != 2:
raise ValueError("Array is the wrong shape")
elif samp.shape[1] != self.num_vars:
raise ValueError("Array is the wrong shape")
# draw points from unit variance, uncorrelated Gaussian
samp = erfinv(2. * samp - 1) * 2. ** 0.5
# rotate and scale to the multivariate normal shape
samp = self.mus[mode] + self.sigmas[mode] * np.einsum('ij,kj->ik',
samp * self.sqeigvalues[mode],
self.eigvectors[mode])
return samp
def _draw_samp(self, size, **kwargs):
try:
mode = kwargs['mode']
except KeyError:
mode = None
return np.squeeze(samp)
if mode is None:
if self.nmodes == 1:
mode = 0
else:
mode = np.argwhere(self.cumweights - np.random.rand() > 0)[0][0]
def sample(self, size=1, mode=None):
"""
Draw, and set, a sample from the multivariate Gaussian.
Parameters
----------
mode: int
Specify which mode to sample from. If not set then a mode is
chosen randomly based on its weight.
"""
if size is None:
size = 1
# samples drawn from unit variance uncorrelated multivariate Gaussian
samps = np.zeros((size, len(self)))
for i in range(size):
inbound = False
......@@ -3074,9 +2913,42 @@ class MultivariateGaussianDist(JointPriorDist):
if not outbound:
inbound = True
return samps
for i, name in enumerate(self.names):
if size == 1:
self.current_sample[name] = samps[:, i].flatten()[0]
else:
self.current_sample[name] = samps[:, i].flatten()
def ln_prob(self, value):
"""
Get the log-probability of a sample. For bounded priors the
probability will not be properly normalised.
Parameters
----------
value: array_like
A 1d vector of the sample, or 2d array of sample values with shape
NxM, where N is the number of samples and M is the number of
parameters.
"""
samp = np.asarray(value)
if len(samp.shape) == 1:
samp = samp.reshape(1, self.num_vars)
if len(samp.shape) != 2:
raise ValueError("Array is the wrong shape")
elif samp.shape[1] != self.num_vars:
raise ValueError("Array is the wrong shape")
# check sample(s) is within bounds
outbounds = np.ones(samp.shape[0], dtype=np.bool)
for s, bound in zip(samp.T, self.bounds.values()):
outbounds = (s < bound[0]) | (s > bound[1])
if np.any(outbounds):
break
def _ln_prob(self, samp, lnprob, outbounds):
lnprob = -np.inf * np.ones(samp.shape[0])
for j in range(samp.shape[0]):
# loop over the modes and sum the probabilities
for i in range(self.nmodes):
......@@ -3084,7 +2956,55 @@ class MultivariateGaussianDist(JointPriorDist):
# set out-of-bounds values to -inf
lnprob[outbounds] = -np.inf
return lnprob
if samp.shape[0] == 1:
return lnprob[0]
else:
return lnprob
def prob(self, samp):
"""
Get the probability of a sample. For bounded priors the
probability will not be properly normalised.
"""
return np.exp(self.ln_prob(samp))
def _get_instantiation_dict(self):
subclass_args = infer_args_from_method(self.__init__)
property_names = [p for p in dir(self.__class__)
if isinstance(getattr(self.__class__, p), property)]
dict_with_properties = self.__dict__.copy()
for key in property_names:
dict_with_properties[key] = getattr(self, key)
instantiation_dict = OrderedDict()
for key in subclass_args:
if isinstance(dict_with_properties[key], list):
value = np.asarray(dict_with_properties[key]).tolist()
else:
value = dict_with_properties[key]
instantiation_dict[key] = value
return instantiation_dict
def __len__(self):
return len(self.names)
def __repr__(self):
"""Overrides the special method __repr__.
Returns a representation of this instance that resembles how it is instantiated.
Works correctly for all child classes
Returns
-------
str: A string representation of this instance
"""
dist_name = self.__class__.__name__
instantiation_dict = self._get_instantiation_dict()
args = ', '.join(['{}={}'.format(key, repr(instantiation_dict[key]))
for key in instantiation_dict])
return "{}({})".format(dist_name, args)
def __eq__(self, other):
if self.__class__ != other.__class__:
......@@ -3122,62 +3042,74 @@ class MultivariateNormalDist(MultivariateGaussianDist):
""" A synonym for the :class:`~bilby.core.prior.MultivariateGaussianDist` distribution."""
class JointPrior(Prior):
def __init__(self, dist, name=None, latex_label=None, unit=None):
if dist.__class__.__bases__ [0] != JointPriorDist:
raise TypeError("Must supply a JointPriorDist object instance to be shared by all joint params")
class MultivariateGaussian(Prior):
if name not in dist.names:
raise ValueError("'{}' is not a parameter in the JointPriorDist")
def __init__(self, mvg, name=None, latex_label=None, unit=None):
"""
A prior class for a multivariate Gaussian (mixture model) prior.
self.dist = dist
super(JointPrior, self).__init__(name=name, latex_label=latex_label, unit=unit,
minimum=dist.bounds[name][0],
maximum=dist.bounds[name][1])
Parameters
----------
mvg: MultivariateGaussianDist
A :class:`bilby.core.prior.MultivariateGaussianDist` object defining
the multivariate Gaussian distribution. This object is not copied,
as it needs to be shared across multiple priors, and as such its
contents will be altered by the prior.
name: str
See superclass
latex_label: str
See superclass
unit: str
See superclass
@property
def minimum(self):
return self._minimum
"""
@minimum.setter
def minimum(self, minimum):
self._minimum = minimum
self.dist.bounds[self.name] = (minimum, self.dist.bounds[self.name][1])
if not isinstance(mvg, MultivariateGaussianDist):
raise TypeError("Must supply a multivariate Gaussian object")
@property
def maximum(self):
return self._maximum
# check name is in the MultivariateGaussianDist class
if name not in mvg.names:
raise ValueError("'{}' is not a parameter in the multivariate "
"Gaussian")
self.mvg = mvg
@maximum.setter
def maximum(self, maximum):
self._maximum = maximum
self.dist.bounds[self.name] = (self.dist.bounds[self.name][0], maximum)
super(MultivariateGaussian, self).__init__(name=name, latex_label=latex_label, unit=unit,
minimum=mvg.bounds[name][0],
maximum=mvg.bounds[name][1])
def rescale(self, val, **kwargs):
def rescale(self, val, mode=None):
"""
Scale a unit hypercube sample to the prior.
Parameters
----------
mode: int
Specify which mode to sample from. If not set then a mode is
chosen randomly based on its weight.
"""
self.test_valid_for_rescaling(val)
Prior.test_valid_for_rescaling(val)
# add parameter value to multivariate Gaussian
self.dist.rescale_parameters[self.name] = val
self.mvg.rescale_parameters[self.name] = val
if self.dist.filled_rescale():
values = np.array(list(self.dist.rescale_parameters.values())).T
samples = self.dist.rescale(values, **kwargs)
self.dist.reset_rescale()
if self.mvg.filled_rescale():
values = np.array(list(self.mvg.rescale_parameters.values())).T
samples = self.mvg.rescale(values, mode=mode)
self.mvg.reset_rescale()
return samples
else:
return [] # return empty list
def sample(self, size=1, **kwargs):
def sample(self, size=1, mode=None):
"""
Draw a sample from the prior.
Parameters
----------
mode: int
Specify which mode to sample from. If not set then a mode is
chosen randomly based on its weight.
Returns
-------
......@@ -3185,25 +3117,41 @@ class JointPrior(Prior):
A sample from the prior paramter.
"""
if self.name in self.dist.sampled_parameters:
if self.name in self.mvg.sampled_parameters:
logger.warning("You have already drawn a sample from parameter "
"'{}'. The same sample will be "
"returned".format(self.name))
if len(self.dist.current_sample) == 0:
if len(self.mvg.current_sample) == 0:
# generate a sample
self.dist.sample(size=size, **kwargs)
self.mvg.sample(size=size, mode=mode)
sample = self.dist.current_sample[self.name]
sample = self.mvg.current_sample[self.name]
if self.name not in self.dist.sampled_parameters:
self.dist.sampled_parameters.append(self.name)
if self.name not in self.mvg.sampled_parameters:
self.mvg.sampled_parameters.append(self.name)
if len(self.dist.sampled_parameters) == len(self.dist):
if len(self.mvg.sampled_parameters) == len(self.mvg):
# reset samples
self.dist.reset_sampled()
self.mvg.reset_sampled()
return sample
def prob(self, val):
"""Return the prior probability of val
Parameters
----------
val: float
Returns
-------
float:
"""
return np.exp(self.ln_prob(val))
def ln_prob(self, val):
"""
Return the natural logarithm of the prior probability. Note that this
......@@ -3212,14 +3160,14 @@ class JointPrior(Prior):
"""
# add parameter value to multivariate Gaussian
self.dist.requested_parameters[self.name] = val
self.mvg.requested_parameters[self.name] = val
if self.dist.filled_request():
if self.mvg.filled_request():
# all required parameters have been set
values = list(self.dist.requested_parameters.values())
values = list(self.mvg.requested_parameters.values())
# check for the same number of values for each parameter
for i in range(len(self.dist) - 1):
for i in range(len(self.mvg) - 1):
if (isinstance(values[i], (list, np.ndarray)) or
isinstance(values[i + 1], (list, np.ndarray))):
if (isinstance(values[i], (list, np.ndarray)) and
......@@ -3231,10 +3179,11 @@ class JointPrior(Prior):
raise ValueError("Each parameter must have the same "
"number of requested values.")
lnp = self.dist.ln_prob(np.asarray(values).T)
lnp = self.mvg.ln_prob(np.asarray(values).T)
# reset the requested parameters
self.dist.reset_request()
self.mvg.reset_request()
return lnp
else:
# if not all parameters have been requested yet, just return 0
......@@ -3252,27 +3201,27 @@ class JointPrior(Prior):
else:
return np.zeros_like(val)
def prob(self, val):
"""Return the prior probability of val
Parameters
----------
val: float
@property
def minimum(self):
return self._minimum
Returns
-------
float:
@minimum.setter
def minimum(self, minimum):
self._minimum = minimum
"""
# update the bounds in the MultivariateGaussianDist
self.mvg.bounds[self.name] = (minimum, self.mvg.bounds[self.name][1])
return np.exp(self.ln_prob(val))
@property
def maximum(self):
return self._maximum
@maximum.setter
def maximum(self, maximum):
self._maximum = maximum
class MultivariateGaussian(JointPrior):
"""
A synonmy class for MultiVariateNormal / deprecated, now use JointPrior with the dist object controlling
the type of joint prior
"""
# update the bounds in the MultivariateGaussianDist
self.mvg.bounds[self.name] = (self.mvg.bounds[self.name][0], maximum)
class MultivariateNormal(MultivariateGaussian):
......
examples/core_examples/multivariate_gaussian_prior.py
+ 2
− 2
View file @ 02d4d60e
......@@ -50,8 +50,8 @@ mvg = bilby.core.prior.MultivariateGaussianDist(names, nmodes=2, mus=mus,
corrcoefs=corrcoefs,
sigmas=sigmas, weights=weights)
priors = dict()
priors['m'] = bilby.core.prior.JointPrior(mvg, 'm')
priors['c'] = bilby.core.prior.JointPrior(mvg, 'c')
priors['m'] = bilby.core.prior.MultivariateGaussian(mvg, 'm')
priors['c'] = bilby.core.prior.MultivariateGaussian(mvg, 'c')
result = bilby.run_sampler(
likelihood=likelihood, priors=priors, sampler='dynesty', nlive=4000,
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment