diff --git a/autofit/mapper/prior/abstract.py b/autofit/mapper/prior/abstract.py index 6ab3061d0..b3048021f 100644 --- a/autofit/mapper/prior/abstract.py +++ b/autofit/mapper/prior/abstract.py @@ -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 @@ -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. diff --git a/autofit/mapper/prior/log_gaussian.py b/autofit/mapper/prior/log_gaussian.py index d244faac8..a97328802 100644 --- a/autofit/mapper/prior/log_gaussian.py +++ b/autofit/mapper/prior/log_gaussian.py @@ -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, @@ -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, diff --git a/autofit/non_linear/clipper.py b/autofit/non_linear/clipper.py index 2503d85e1..752fefb5c 100644 --- a/autofit/non_linear/clipper.py +++ b/autofit/non_linear/clipper.py @@ -99,8 +99,6 @@ import numpy as np -from autofit.mapper.prior.log_gaussian import LogGaussianPrior - logger = logging.getLogger(__name__) @@ -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 ---------- @@ -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. @@ -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) diff --git a/test_autofit/mapper/prior/test_log_gaussian.py b/test_autofit/mapper/prior/test_log_gaussian.py index 6bef0d064..596fa7b59 100644 --- a/test_autofit/mapper/prior/test_log_gaussian.py +++ b/test_autofit/mapper/prior/test_log_gaussian.py @@ -1,5 +1,7 @@ import pickle +from copy import copy +import numpy as np import pytest import autofit as af @@ -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 diff --git a/test_autofit/mapper/prior/test_prior_properties.py b/test_autofit/mapper/prior/test_prior_properties.py index 5ca42a9e8..bca1c322e 100644 --- a/test_autofit/mapper/prior/test_prior_properties.py +++ b/test_autofit/mapper/prior/test_prior_properties.py @@ -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) diff --git a/test_autofit/non_linear/test_clipper.py b/test_autofit/non_linear/test_clipper.py index 0e5a89c2e..61d9133ed 100644 --- a/test_autofit/non_linear/test_clipper.py +++ b/test_autofit/non_linear/test_clipper.py @@ -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(