r"""Doublet-DALI likelihood and posterior sampling.
At ``order = 1`` the log-likelihood is exactly quadratic, so the Fisher posterior is
Gaussian and :math:`\Sigma = F^{-1}` describes it completely -- every marginal is
an ellipse available in closed form.
The doublet adds cubic and quartic terms,
.. math::
\log L = -\tfrac{1}{2}F_{ab}\Delta^a\Delta^b
- \tfrac{1}{2}G_{abc}\Delta^a\Delta^b\Delta^c
- \tfrac{1}{8}H_{abcd}\Delta^a\Delta^b\Delta^c\Delta^d,
which has neither a closed-form normalisation nor closed-form marginals.
The contractions are written as flattened matrix products rather than
:func:`numpy.einsum` so they reach BLAS: measured at 1.9 million evaluations per
second versus 0.03 million for the equivalent ``einsum`` call.
"""
import logging
import numpy
from bilby.core.likelihood import Likelihood
from bilby.core.prior import Constraint, PriorDict
# Hard physical ranges applied when ``enforce_physicality`` is on. Sampling a
# DALI posterior is sampling a truncated Taylor expansion, which knows nothing
# about these and will happily wander past them.
[docs]
PHYSICAL_RANGES = {
"symmetric_mass_ratio": (0.0, 0.25),
"mass_ratio": (0.0, 1.0),
"chi_1": (-1.0, 1.0),
"chi_2": (-1.0, 1.0),
"a_1": (0.0, 1.0),
"a_2": (0.0, 1.0),
"chirp_mass": (0.0, None),
"total_mass": (0.0, None),
"mass_1": (0.0, None),
"mass_2": (0.0, None),
"luminosity_distance": (0.0, None),
"theta_jn": (0.0, numpy.pi),
"tilt_1": (0.0, numpy.pi),
"tilt_2": (0.0, numpy.pi),
"dec": (-numpy.pi / 2, numpy.pi / 2),
"lambda_1": (0.0, None),
"lambda_2": (0.0, None),
}
[docs]
class DALILikelihood(Likelihood):
"""Doublet-DALI log-likelihood as a bilby likelihood.
Examples
--------
>>> likelihood = DALILikelihood(fisher)
>>> result = sample_posterior(likelihood, priors, outdir="out", label="src0")
Attributes
----------
names : list of str
Parameter names, in tensor order.
fiducial : dict
The point the expansion is about; samples are offsets from it.
"""
def __init__(self, fisher, order=2):
"""
Parameters
----------
fisher : GWForge.fisher.matrix.FisherMatrix
A matrix built with a matching ``order``.
order : int
``1`` reduces this to the Gaussian Fisher likelihood, which is
useful as a cross-check that the sampler reproduces ``F**-1``.
``2`` is doublet-DALI.
Raises
------
ValueError
If ``order = 2`` but the Fisher matrix carries no DALI tensors.
"""
[docs]
self.names = list(fisher.names)
[docs]
self.fiducial = {name: fisher.parameters[name] for name in self.names}
[docs]
self.order = int(order)
count = len(self.names)
self._fisher = numpy.asarray(fisher.matrix, dtype=float)
if self.order >= 2:
if fisher.doublet is None or fisher.quartic is None:
raise ValueError(
"order = 2 needs the DALI tensors; build the FisherMatrix "
"with order=2."
)
self._doublet = numpy.asarray(fisher.doublet, dtype=float).reshape(
count * count, count
)
self._quartic = numpy.asarray(fisher.quartic, dtype=float).reshape(
count * count, count * count
)
else:
self._doublet = self._quartic = None
super().__init__(parameters={name: None for name in self.names})
[docs]
def log_likelihood(self, parameters=None):
"""Log-likelihood at ``parameters`` (or at ``self.parameters``).
Parameters
----------
parameters : dict or None
Absolute parameter values, not offsets. Accepting them as an
argument is the bilby 3 interface; bilby 2 still sets
``self.parameters`` instead, and both work here.
Returns
-------
float
"""
source = self.parameters if parameters is None else parameters
offset = numpy.fromiter(
(source[name] - self.fiducial[name] for name in self.names),
dtype=float,
count=len(self.names),
)
return float(self.log_likelihood_from_offsets(offset[numpy.newaxis, :])[0])
[docs]
def log_likelihood_from_offsets(self, offsets):
"""Vectorised log-likelihood for a stack of parameter offsets.
Parameters
----------
offsets : numpy.ndarray
``(n_points, n_parameters)`` displacements from the fiducial point.
Returns
-------
numpy.ndarray
``(n_points,)`` log-likelihood values.
"""
offsets = numpy.atleast_2d(numpy.asarray(offsets, dtype=float))
value = -0.5 * numpy.einsum(
"wa,ab,wb->w", offsets, self._fisher, offsets, optimize=True
)
if self.order >= 2:
outer = (offsets[:, :, None] * offsets[:, None, :]).reshape(
len(offsets), -1
)
value = value - 0.5 * ((outer @ self._doublet) * offsets).sum(axis=1)
value = value - 0.125 * ((outer @ self._quartic) * outer).sum(axis=1)
return value
[docs]
def noise_log_likelihood(self):
"""Zero: the expansion is about the maximum, where the residual vanishes."""
return 0.0
# Warn when a prior spans fewer than this many Fisher sigmas either side of the
# fiducial value. A prior narrower than the likelihood silently dominates the
# posterior, and the resulting "DALI" width is really just the prior width --
# an easy mistake to make when the prior file is written by hand.
# Will use it for lower SNR signals :)
[docs]
NARROW_PRIOR_SIGMAS = 3.0
[docs]
def build_priors(prior_file, names, fiducial, enforce_physicality=False, sigmas=None):
"""Load and validate the sampling priors.
Parameters
----------
prior_file : str
Path to a bilby ``.prior`` file.
names : sequence of str
Fisher parameter names that must all be covered.
fiducial : dict
Fiducial values, checked to lie inside the prior support -- a prior that
excludes the point the expansion is about would be sampling nothing.
enforce_physicality : bool
Also reject samples outside :data:`PHYSICAL_RANGES`.
sigmas : dict or None
Fisher one-sigma widths. When given, priors narrower than
:data:`NARROW_PRIOR_SIGMAS` sigmas are warned about.
Returns
-------
bilby.core.prior.PriorDict
Raises
------
KeyError
If the file does not define every parameter. bilby's own prior files use
``mass_ratio``, ``a_1`` and friends, which do not line up with the
default Fisher parameters, so this is a likely mistake and is worth
failing loudly for rather than silently sampling the wrong space.
ValueError
If a fiducial value falls outside its prior.
"""
priors = PriorDict(filename=prior_file)
missing = [
name
for name in names
if name not in priors or isinstance(priors.get(name), Constraint)
]
if missing:
raise KeyError(
"The prior file {} does not define a sampling prior for: {}. The "
"Fisher parameters are {}.".format(
prior_file, ", ".join(sorted(missing)), ", ".join(names)
)
)
for name in list(priors):
if name not in names:
del priors[name]
outside = [
name
for name in names
if not numpy.isfinite(priors[name].ln_prob(fiducial[name]))
]
if outside:
raise ValueError(
"The fiducial values for {} lie outside the prior in {}. The DALI "
"expansion is about those values, so the prior has to contain "
"them.".format(", ".join(sorted(outside)), prior_file)
)
if enforce_physicality:
_add_physical_constraints(priors, names)
if sigmas:
_warn_about_narrow_priors(priors, names, fiducial, sigmas)
return priors
def _warn_about_narrow_priors(priors, names, fiducial, sigmas):
"""Flag priors too tight for the posterior they are supposed to bound."""
for name in names:
sigma = sigmas.get(name)
prior = priors[name]
lower = getattr(prior, "minimum", None)
upper = getattr(prior, "maximum", None)
if not sigma or lower is None or upper is None:
continue
room = min(fiducial[name] - lower, upper - fiducial[name]) / sigma
if room < NARROW_PRIOR_SIGMAS:
logging.warning(
"The prior on %s reaches only %.1f sigma from the fiducial "
"value, so it will truncate the posterior rather than bound it "
"-- the sampled width will mostly reflect the prior. Widen it "
"to at least %.0f sigma (%.3e) unless the truncation is "
"intended.",
name,
room,
NARROW_PRIOR_SIGMAS,
NARROW_PRIOR_SIGMAS * sigma,
)
def _add_physical_constraints(priors, names):
"""Clip each prior to its physical range, in place."""
for name in names:
limits = PHYSICAL_RANGES.get(name)
if limits is None:
continue
lower, upper = limits
prior = priors[name]
if lower is not None and getattr(prior, "minimum", None) is not None:
prior.minimum = max(prior.minimum, lower)
if upper is not None and getattr(prior, "maximum", None) is not None:
prior.maximum = min(prior.maximum, upper)
if (
getattr(prior, "minimum", None) is not None
and getattr(prior, "maximum", None) is not None
and prior.minimum >= prior.maximum
):
raise ValueError(
"Enforcing physicality leaves {} with an empty prior range "
"[{}, {}].".format(name, prior.minimum, prior.maximum)
)
# If I don't stick to bilby interface this can be sped up by vectorisation.
# But I am being lazy and this is fast enough for now.
[docs]
def sample_posterior(
likelihood,
priors,
outdir,
label,
sampler="emcee",
nwalkers=200,
nsteps=2000,
nburn=500,
**kwargs,
):
"""Sample a DALI posterior through bilby.
Parameters
----------
likelihood : DALILikelihood
priors : bilby.core.prior.PriorDict
outdir : str
Output directory for the sampler.
label : str
Run label, used for the output filenames.
sampler : str
Any bilby sampler; ``emcee`` by default.
nwalkers, nsteps, nburn : int
emcee settings, ignored by nested samplers.
**kwargs
Passed straight to :func:`bilby.run_sampler`.
Returns
-------
bilby.core.result.Result
"""
from bilby import run_sampler
settings = dict(outdir=outdir, label=label, sampler=sampler, injection_parameters=likelihood.fiducial)
if sampler == "emcee":
settings.update(nwalkers=nwalkers, iterations=nsteps, nburn=nburn)
settings.update(kwargs)
logging.info("Sampling the order-%d DALI posterior with %s", likelihood.order, sampler)
return run_sampler(likelihood=likelihood, priors=priors, **settings)