Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
26 changes: 23 additions & 3 deletions autofit/mapper/prior/abstract.py
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,20 @@ class Prior(Variable, ABC, ArithmeticMixin):

_ids = itertools.count()

#: Whether the support *excludes* the corresponding limit itself, i.e. whether
#: the bound is ``value > limit`` rather than ``value >= limit``.
#:
#: Most priors are inclusive at their limits: ``UniformPrior.log_prior_from_value``
#: is finite exactly on the bound. ``LogGaussianPrior`` is not — its support is the
#: open ``(0, inf)`` — and a consumer that clips onto an exclusive bound lands on a
#: point of zero density. These flags let such a consumer tell the two apart without
#: a per-type ``isinstance`` switch (see ``non_linear.clipper.ClipperPriorBox``).
#:
#: Class attributes, deliberately: they are a property of the prior *family*, never
#: of an instance, so they stay out of ``__dict__`` and out of the identifier.
lower_limit_strict = False
upper_limit_strict = False

def __init__(self, message, id_=None):
"""
An object used to mappers a unit value to an attribute value for a specific
Expand Down Expand Up @@ -347,10 +361,16 @@ def name_of_class(cls) -> str:
def limits(self) -> Tuple[float, float]:
"""The (lower, upper) bounds of this prior.

Returns (-inf, inf) by default. Subclasses with finite bounds
(e.g. UniformPrior) override this.
Derived from ``lower_limit`` / ``upper_limit`` rather than stated
separately, so the two cannot disagree. They used to: this property
returned a hardcoded ``(-inf, inf)`` for every prior that did not
override it, which was right for ``GaussianPrior`` and wrong for
``LogGaussianPrior``, whose support is ``(0, inf)``.

Use ``lower_limit_strict`` / ``upper_limit_strict`` to tell whether the
support includes the bounds returned here.
"""
return (float("-inf"), float("inf"))
return (float(self.lower_limit), float(self.upper_limit))

def gaussian_prior_model_for_arguments(self, arguments):
"""Look up this prior in an arguments dict and return the mapped value.
Expand Down
13 changes: 13 additions & 0 deletions autofit/mapper/prior/log_gaussian.py
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,10 @@ class LogGaussianPrior(Prior):
__identifier_fields__ = ("mean", "sigma")
__database_args__ = ("mean", "sigma", "id_")

#: The support is the *open* interval ``(0, inf)``: ``log_prior_from_value``
#: returns ``-inf`` at ``0`` itself, so a consumer must stay strictly above it.
lower_limit_strict = True

def __init__(
self,
mean: float,
Expand Down Expand Up @@ -51,6 +55,15 @@ def __init__(
self.mean = mean
self.sigma = sigma

# Declared on the prior, not on the message below. `TransformedMessage`
# defaults its limits to +/-inf and derives its `_support` separately, so
# a prior that left them to delegation reported (-inf, inf) for a strictly
# positive parameter -- the bug this shadowing fixes. Keeping the message's
# own limits untouched keeps the EP/Laplace machinery, which reads them,
# behaving exactly as before.
self.lower_limit = 0.0
self.upper_limit = float("inf")

message = TransformedMessage(
NormalMessage(mean, sigma),
log_transform,
Expand Down
45 changes: 19 additions & 26 deletions autofit/non_linear/clipper.py
Original file line number Diff line number Diff line change
Expand Up @@ -99,8 +99,6 @@

import numpy as np

from autofit.mapper.prior.log_gaussian import LogGaussianPrior

logger = logging.getLogger(__name__)


Expand Down Expand Up @@ -225,7 +223,9 @@ class ClipperPriorBox(AbstractClipper):
``(0, inf)``) — inset by an *absolute* ``strict_epsilon``. A relative margin
is identically zero here, there being no finite width, so it would clip
exactly onto ``0.0``, where the support is strict and ``log_prior`` is
``-inf``. Only an absolute nudge lands strictly inside.
``-inf``. Only an absolute nudge lands strictly inside. Which bounds are
exclusive is read off the prior's ``lower_limit_strict`` /
``upper_limit_strict``, never inferred from its type.

Parameters
----------
Expand All @@ -249,23 +249,20 @@ def _limits_from_model(self, model):
The raw ``(lower, upper, lower_strict, upper_strict)`` arrays for ``model``,
in physical parameter order.

Limits are read off ``prior.lower_limit`` / ``prior.upper_limit``, which
resolves for **every** prior type without a type switch: ``Prior.__getattr__``
delegates to the prior's message, and ``AbstractMessage`` defaults both to
``±inf``. ``UniformPrior`` and ``LogUniformPrior`` shadow them with their
own attributes, ``TruncatedGaussianPrior`` picks up real limits from
``TruncatedNormalMessage``, and ``GaussianPrior`` correctly falls through
to ``±inf``.

``LogGaussianPrior`` is the one prior that read gets *wrong*, and silently.
Its message is a ``TransformedMessage``, which defaults its limits to
``±inf`` and is never passed any — yet
``LogGaussianPrior.log_prior_from_value`` returns ``-inf`` for
``value <= 0``. Left uncorrected, the clipper would report that coordinate
as unbounded and fail to protect precisely the mechanism it exists to fix,
so its ``(0, inf)`` support is declared here. Declaring it on the prior
itself is the cleaner fix, but it would change a class the EP machinery and
the nested samplers also read, so it is deliberately kept local.
Limits are read off ``prior.lower_limit`` / ``prior.upper_limit``, and
strictness off ``prior.lower_limit_strict`` / ``prior.upper_limit_strict``.
Both resolve for **every** prior type without a type switch: ``UniformPrior``
and ``LogUniformPrior`` shadow the limits with their own attributes,
``LogGaussianPrior`` declares its ``(0, inf)`` support and flags the lower
bound strict, ``TruncatedGaussianPrior`` picks up real limits from
``TruncatedNormalMessage``, and ``GaussianPrior`` falls through
``Prior.__getattr__`` to its message's ``±inf``.

This clipper used to declare ``LogGaussianPrior``'s support itself, because
that prior reported ``(-inf, inf)`` while ``log_prior_from_value`` returned
``-inf`` for ``value <= 0``. That was a workaround for a defect in the prior,
fixed in PyAutoFit#1526: any consumer of ``lower_limit`` — not just this one —
was being told a strictly positive parameter could go negative.

The ``strict`` flags mark bounds the support *excludes* (``value > limit``
rather than ``value >= limit``); only those need the absolute inset.
Expand All @@ -275,12 +272,8 @@ def _limits_from_model(self, model):
for prior in model.priors_ordered_by_id:
low = float(prior.lower_limit)
high = float(prior.upper_limit)
low_strict = False
high_strict = False

if isinstance(prior, LogGaussianPrior):
low = 0.0
low_strict = True
low_strict = bool(prior.lower_limit_strict)
high_strict = bool(prior.upper_limit_strict)

lower.append(low)
upper.append(high)
Expand Down
128 changes: 128 additions & 0 deletions test_autofit/mapper/prior/test_log_gaussian.py
Original file line number Diff line number Diff line change
@@ -1,5 +1,7 @@
import pickle
from copy import copy

import numpy as np
import pytest

import autofit as af
Expand Down Expand Up @@ -30,3 +32,129 @@ def test_pickle(log_gaussian):

def test_identifier(log_gaussian):
Identifier(log_gaussian)


# === PyAutoFit#1526: the prior declares its own (0, inf) support ===
#
# The support was always (0, inf) -- ``log_prior_from_value`` returns -inf for
# ``value <= 0`` -- but the prior reported (-inf, inf), because ``Prior.__getattr__``
# delegated to a ``TransformedMessage`` whose limits default to +/-inf and were never
# set. ``ClipperPriorBox`` worked around it with an ``isinstance`` switch; every other
# consumer of ``lower_limit`` was simply told the wrong thing.
#
# The pinned values below are the pre-change ones, measured on the commit before the
# fix. They are the guarantee that declaring the support moved *what the prior says*
# and nothing about *what it computes*.


def test__reports_its_own_support(log_gaussian):
assert log_gaussian.lower_limit == 0.0
assert log_gaussian.upper_limit == float("inf")
assert log_gaussian.limits == (0.0, float("inf"))


def test__lower_bound_is_strict__upper_is_not(log_gaussian):
"""
The support is the *open* (0, inf): ``log_prior_from_value(0.0)`` is -inf, so a
consumer that clips onto the reported bound must know to stay strictly above it.
"""
assert log_gaussian.lower_limit_strict is True
assert log_gaussian.upper_limit_strict is False

assert log_gaussian.log_prior_from_value(log_gaussian.lower_limit) == -np.inf


def test__other_prior_families_keep_their_limits():
"""
``Prior.limits`` now derives from ``lower_limit``/``upper_limit`` rather than
returning a hardcoded (-inf, inf). GaussianPrior must still be unbounded.
"""
assert af.GaussianPrior(mean=0.0, sigma=1.0).limits == (-np.inf, np.inf)
assert af.UniformPrior(lower_limit=0.0, upper_limit=2.0).limits == (0.0, 2.0)
assert af.LogUniformPrior(lower_limit=0.01, upper_limit=100.0).limits == (
0.01,
100.0,
)

for cls in (af.GaussianPrior, af.UniformPrior, af.LogUniformPrior):
assert cls.lower_limit_strict is False
assert cls.upper_limit_strict is False


@pytest.mark.parametrize(
"value, expected",
[
(-3.0, -np.inf),
(-1e-09, -np.inf),
(0.0, -np.inf),
(1e-12, -204.83588563011642),
(0.001, -8.892033838618836),
(0.1, 0.14164835190717184),
(0.5, 0.3396055360729164),
(1.0, -0.04733727810650888),
(2.0, -0.7185718164978876),
(10.0, -3.3735407249713125),
(1000.0, -19.437601069254285),
],
)
def test__log_prior_from_value_is_unchanged(value, expected):
"""
The density was always correct; only the reported limits were wrong. Pinned
against the pre-change values either side of zero.
"""
prior = af.LogGaussianPrior(mean=0.4, sigma=1.3)
assert prior.log_prior_from_value(value) == pytest.approx(expected, rel=1e-12)


@pytest.mark.parametrize(
"unit, expected",
[
(1e-09, 0.0006129978595719644),
(0.0001, 0.011858368820420245),
(0.01, 0.07249394511581292),
(0.1, 0.2819523947584061),
(0.25, 0.6207439038545902),
(0.5, 1.4918246976412703),
(0.75, 3.5852803622761034),
(0.9, 7.893321602745905),
(0.99, 30.69968015862632),
(0.999999999, 3630.585196715403),
],
)
def test__unit_cube_mapping_is_unchanged(unit, expected):
"""
The nested samplers work in unit-cube coordinates and map through the prior. The
limits live on the prior while the mapping lives on the message stack, which this
change does not touch -- so no stored nested-sampling result shifts.
"""
prior = af.LogGaussianPrior(mean=0.4, sigma=1.3)
assert prior.value_for(unit) == pytest.approx(expected, rel=1e-12)
assert prior.unit_value_for(prior.value_for(unit)) == pytest.approx(unit, rel=1e-6)


def test__identifier_is_unchanged():
"""
If the declared limits fed the identifier, every existing output directory would
re-key and its stored results would be orphaned. ``__identifier_fields__`` is
("mean", "sigma"), and the limits are instance/class attributes outside it, so the
hash is untouched. Pinned to the pre-change value.
"""
prior = af.LogGaussianPrior(mean=0.4, sigma=1.3)
assert str(Identifier(prior)) == "34cb61ade6bafa6050229e8b6b390235"


def test__declared_support_survives_copy_pickle_and_projection(log_gaussian):
"""
The limits are derived in ``__init__`` rather than carried as parameters, which is
what makes them survive every path that rebuilds the prior -- including the JAX
pytree round-trip, where only (mean, sigma, id) are flattened.
"""
assert copy(log_gaussian).lower_limit == 0.0
assert pickle.loads(pickle.dumps(log_gaussian)).lower_limit == 0.0

samples = np.exp(np.random.default_rng(0).normal(1.0, 2.0, 500))
projected = log_gaussian.project(samples, np.zeros(500))
assert projected.lower_limit == 0.0

rebuilt = af.LogGaussianPrior.tree_unflatten((), log_gaussian.tree_flatten()[0])
assert rebuilt.lower_limit == 0.0
80 changes: 80 additions & 0 deletions test_autofit/mapper/prior/test_prior_properties.py
Original file line number Diff line number Diff line change
Expand Up @@ -271,3 +271,83 @@ def test__from_mode_matches_mean_and_variance(cls, mean, variance):
message = cls.from_mode(np.asarray(mean), np.asarray(variance))
assert float(message.mean) == pytest.approx(mean, rel=1e-6)
assert float(message.variance) == pytest.approx(variance, rel=1e-6)


# === P6: the reported support matches the actual support ===


def interior_probes(prior):
"""
Points that lie strictly inside the prior's *reported* ``limits``.

Every one of them must have finite log prior. If a prior reports a box wider
than its true support, some point in here falls in the gap and the density
there is -inf -- which is precisely the shape of PyAutoFit#1526, where
``LogGaussianPrior`` reported (-inf, inf) for a support of (0, inf).
"""
lo, hi = prior.limits
probes = []

if np.isfinite(lo) and np.isfinite(hi):
width = hi - lo
probes += [lo + 0.25 * width, lo + 0.5 * width, lo + 0.75 * width]
probes += [lo + 1e-9 * width, hi - 1e-9 * width]
elif np.isfinite(lo):
probes += [lo + 1e-9, lo + 1.0, lo + 1e3]
elif np.isfinite(hi):
probes += [hi - 1e-9, hi - 1.0, hi - 1e3]
else:
probes += [-1e6, -1.0, 0.0, 1.0, 1e3]

return probes


@pytest.mark.parametrize("prior", all_priors(), ids=prior_id)
def test__log_prior_is_finite_everywhere_inside_the_reported_limits(prior):
"""
The general form of PyAutoFit#1526. A prior that reports a bound it does not
actually have hands every consumer a licence to evaluate where the density is
zero -- and for a strictly positive parameter, ``log(0)`` or a division by it
is the failure that follows.
"""
for value in interior_probes(prior):
assert np.isfinite(
prior.log_prior_from_value(value)
), f"{prior} reports limits {prior.limits} but log_prior({value}) is not finite"


@pytest.mark.parametrize("prior", all_priors(), ids=prior_id)
def test__log_prior_is_minus_inf_outside_the_reported_limits(prior):
"""
The converse: the reported box must not be *narrower* than the support either.
"""
lo, hi = prior.limits

if np.isfinite(lo):
assert prior.log_prior_from_value(lo - 1.0) == -np.inf
if np.isfinite(hi):
assert prior.log_prior_from_value(hi + 1.0) == -np.inf


@pytest.mark.parametrize("prior", all_priors(), ids=prior_id)
def test__strictness_flags_agree_with_the_density_at_the_bounds(prior):
"""
A bound flagged strict must have zero density on it; a bound not flagged strict
must have finite density on it. This is what lets a consumer clip onto a bound
without a per-type switch -- ``ClipperPriorBox`` insets only the strict ones.
"""
lo, hi = prior.limits

if np.isfinite(lo):
at_lower = prior.log_prior_from_value(lo)
if prior.lower_limit_strict:
assert at_lower == -np.inf
else:
assert np.isfinite(at_lower)

if np.isfinite(hi):
at_upper = prior.log_prior_from_value(hi)
if prior.upper_limit_strict:
assert at_upper == -np.inf
else:
assert np.isfinite(at_upper)
18 changes: 13 additions & 5 deletions test_autofit/non_linear/test_clipper.py
Original file line number Diff line number Diff line change
Expand Up @@ -69,16 +69,24 @@ def test__gaussian_prior__passes_through_as_infinite(self):
assert not np.isnan(lower).any()
assert not np.isnan(upper).any()

def test__log_gaussian_prior__lower_bound_is_declared_by_the_clipper(self):
def test__log_gaussian_prior__lower_bound_is_declared_by_the_prior(self):
"""
LogGaussianPrior reports (-inf, inf) because its TransformedMessage is never
given limits, yet ``log_prior_from_value`` is -inf for value <= 0. The
clipper declares the real (0, inf) support, strictly, so the projected value
LogGaussianPrior declares its own (0, inf) support and flags the lower bound
strict, so the clipper reads it like any other prior and the projected value
lands *above* zero rather than on it.

This test used to assert ``prior.lower_limit == -np.inf`` and the clipper
supplied the real bound itself via an ``isinstance`` switch. That was a
workaround for PyAutoFit#1526: the prior was telling every consumer, not just
this one, that a strictly positive parameter could go negative. The bounds
assertions below are unchanged by the fix — that equivalence is the point.
"""
prior = af.LogGaussianPrior(mean=0.0, sigma=1.0)

assert prior.lower_limit == -np.inf
assert prior.lower_limit == 0.0
assert prior.upper_limit == np.inf
assert prior.lower_limit_strict is True
assert prior.upper_limit_strict is False

model = _model(alpha=prior)
lower, upper = ClipperPriorBox(strict_epsilon=1.0e-12).bounds_from_model(
Expand Down
Loading