diff --git a/autofit/mapper/prior/log_uniform.py b/autofit/mapper/prior/log_uniform.py index 06bdcdebf..8ae4e7ed4 100644 --- a/autofit/mapper/prior/log_uniform.py +++ b/autofit/mapper/prior/log_uniform.py @@ -123,22 +123,18 @@ def log_prior_from_value(self, value, xp=np): used by ``UniformPrior.log_prior_from_value`` which drops ``-log(b - a)`` to return ``0.0``. - Non-positive ``value`` (``value <= 0``) returns ``-inf``. Emcee's stretch - move proposes physical values that can leave the support and go - non-positive; ``-log`` of a non-positive value is ``NaN``, which propagates - into the summed figure-of-merit and crashes the search with - ``ValueError: Probability function returned NaN``. Returning ``-inf`` - (zero density -> rejected move) keeps the figure-of-merit finite. The - "double where" pattern (a safe surrogate inside the ``log``) ensures no + Outside ``[lower_limit, upper_limit]`` the density is zero, so ``-inf`` is + returned — on both the NumPy and JAX paths. This bound test is the only + support enforcement in the search fitness path (see + ``UniformPrior.log_prior_from_value`` and PyAutoFit#1489): for MCMC an + out-of-box proposal becomes a rejected move. It also keeps the + figure-of-merit finite where Emcee's stretch move proposes a non-positive + value: ``-log`` of a non-positive value is ``NaN``, which would otherwise + crash the search with ``ValueError: Probability function returned NaN``. + The "double where" pattern (a safe surrogate inside the ``log``) ensures no ``log`` of a non-positive value is evaluated, avoiding NumPy ``RuntimeWarning``s. - The NumPy path is otherwise unnormalised and unbounded: for any positive - ``value`` it returns ``-log(value)`` regardless of ``[lower_limit, - upper_limit]`` (dropping the normalisation constant, matching - ``UniformPrior``'s convention of returning ``0.0``). The JAX path - additionally returns ``-inf`` outside ``[lower_limit, upper_limit]``. - Parameters ---------- value @@ -147,10 +143,9 @@ def log_prior_from_value(self, value, xp=np): xp Array-module to dispatch on (``numpy`` or ``jax.numpy``). Default ``numpy``. """ - if xp is np: - positive = value > 0.0 - return xp.where(positive, -xp.log(xp.where(positive, value, 1.0)), -xp.inf) in_bounds = (value >= self.lower_limit) & (value <= self.upper_limit) + if xp is np: + return xp.where(in_bounds, -xp.log(xp.where(in_bounds, value, 1.0)), -xp.inf) return xp.where(in_bounds, -xp.log(value), -xp.inf) def log_normalisation(self, xp=np) -> float: diff --git a/autofit/mapper/prior/uniform.py b/autofit/mapper/prior/uniform.py index e2daa5397..4390cdf43 100644 --- a/autofit/mapper/prior/uniform.py +++ b/autofit/mapper/prior/uniform.py @@ -179,10 +179,20 @@ def log_prior_from_value(self, value, xp=np): This is used by certain non-linear searches (e.g. Emcee) in the log likelihood function evaluation. - For a UniformPrior this is always zero, provided the value is between the lower and upper limit. + For a UniformPrior this is zero inside ``[lower_limit, upper_limit]`` (the constant + ``-log(upper - lower)`` is dropped, see ``log_normalisation``) and ``-inf`` outside it, + on both the NumPy and JAX paths. + + The bound test here is the only support enforcement in the search fitness path: since + the removal of ``assert_within_limits`` / ``PriorLimitException`` (shipped 2025.10.16.1), + ``instance_from_vector`` accepts out-of-support vectors without raising, so a + ``log_prior_from_value`` that skipped the test left the box entirely unenforced for + the NumPy-path searches that form a log posterior (Emcee, Zeus, Drawer, LBFGS) — + see PyAutoFit#1489. The ``-inf`` is what makes an out-of-box MCMC proposal a + rejected move. """ - if xp is np: - return 0.0 + if xp is np and np.ndim(value) == 0: + return 0.0 if self.lower_limit <= value <= self.upper_limit else -np.inf in_bounds = (value >= self.lower_limit) & (value <= self.upper_limit) return xp.where(in_bounds, xp.zeros_like(value), -xp.inf) diff --git a/autofit/non_linear/clipper.py b/autofit/non_linear/clipper.py index ceae70aae..eb11ea97b 100644 --- a/autofit/non_linear/clipper.py +++ b/autofit/non_linear/clipper.py @@ -22,6 +22,15 @@ the walker simply stays put. **Rejection is the restoring mechanism the gradient methods lack**, and a ``Clipper`` supplies it. +(Historical correction, PyAutoFit#1489: the rejection sentence above was written +while it was FALSE on the NumPy path — between release 2025.10.16.1, which removed +``assert_within_limits`` / ``PriorLimitException``, and the strict +``log_prior_from_value`` bounds restored for #1489, a NumPy-path ``UniformPrior`` +returned ``0.0`` outside its box, so there was no ``-inf`` for the MCMC samplers +to reject and walkers escaped the declared support. It was true before +2025.10.16.1, via exception-driven resampling, and is true again now via the +strict priors. The gradient-method reasoning below was unaffected.) + Two consumers, one source of truth ---------------------------------- diff --git a/autofit/non_linear/fitness.py b/autofit/non_linear/fitness.py index 2604d2569..e2367bef0 100644 --- a/autofit/non_linear/fitness.py +++ b/autofit/non_linear/fitness.py @@ -312,7 +312,13 @@ def log_likelihood_from(self, figure_of_merit, parameters): if not self.fom_is_log_likelihood: log_prior_list = np.array(self.model.log_prior_list_from_vector(vector=parameters, xp=np)) - log_likelihood = log_likelihood - np.sum(log_prior_list) + log_prior_sum = np.sum(log_prior_list) + # A non-finite prior sum marks a rejected out-of-support point (the strict + # priors return -inf outside their bounds, PyAutoFit#1489). The figure of + # merit there is -inf too, and subtracting -inf from -inf is NaN — keep the + # -inf so history / quick-update bookkeeping stays comparable. + if np.isfinite(log_prior_sum): + log_likelihood = log_likelihood - log_prior_sum return log_likelihood diff --git a/autofit/non_linear/search/mle/bfgs/search.py b/autofit/non_linear/search/mle/bfgs/search.py index 27fdc718a..f588663b4 100644 --- a/autofit/non_linear/search/mle/bfgs/search.py +++ b/autofit/non_linear/search/mle/bfgs/search.py @@ -244,11 +244,12 @@ def _fit( # ``minimize(bounds=None)`` is the same call as omitting it, so the # default path is unchanged. # - # Only the JAX branch is actually *exposed* to the prior-support - # problem: ``UniformPrior.log_prior_from_value`` returns ``0.0`` - # unconditionally on the NumPy path, with no bound test, so there - # is no hard wall to fall off there. Bounds are still passed on - # both branches — they are correct on both, and having them + # Both branches are exposed to the prior-support problem: + # ``UniformPrior.log_prior_from_value`` returns ``-inf`` outside + # its bounds on the NumPy and JAX paths alike (PyAutoFit#1489 + # restored the NumPy-side wall), so an unclipped step out of the + # box makes the objective non-finite on either branch. Bounds are + # passed on both — they are correct on both, and having them # diverge by branch would be a trap of its own. if analysis._use_jax: diff --git a/test_autofit/mapper/prior/test_prior.py b/test_autofit/mapper/prior/test_prior.py index 64f9702c1..1adada76b 100644 --- a/test_autofit/mapper/prior/test_prior.py +++ b/test_autofit/mapper/prior/test_prior.py @@ -175,6 +175,16 @@ def test__log_prior_from_value(self): assert log_prior == 0.0 + # Outside the declared box the density is zero: -inf, so the searches + # that form a log posterior reject the sample (PyAutoFit#1489). + log_prior = gaussian_simple.log_prior_from_value(value=71.0) + + assert log_prior == float("-inf") + + log_prior = gaussian_simple.log_prior_from_value(value=-41.0) + + assert log_prior == float("-inf") + class TestLogUniformPrior: def test__simple_assumptions(self): @@ -192,10 +202,12 @@ def test__non_zero_lower_limit(self): assert log_uniform_half.value_for(0.5) == pytest.approx(0.70710678118, 1.0e-4) def test__log_prior_from_value(self): - # LogUniformPrior log-density: -log(value), dropping the normalisation - # constant -log(log(upper / lower)). Consistent with UniformPrior's - # convention of returning 0.0 (dropping -log(b - a)). - log_uniform = af.LogUniformPrior(lower_limit=1e-8, upper_limit=1.0) + # LogUniformPrior log-density inside the support: -log(value), dropping + # the normalisation constant -log(log(upper / lower)). Consistent with + # UniformPrior's convention of returning 0.0 (dropping -log(b - a)). + # Outside [lower_limit, upper_limit] the density is zero: -inf + # (PyAutoFit#1489). + log_uniform = af.LogUniformPrior(lower_limit=1e-8, upper_limit=10.0) assert log_uniform.log_prior_from_value(value=1.0) == 0.0 assert log_uniform.log_prior_from_value(value=2.0) == pytest.approx( @@ -204,10 +216,11 @@ def test__log_prior_from_value(self): assert log_uniform.log_prior_from_value(value=4.0) == pytest.approx( -np.log(4.0), 1.0e-12 ) + assert log_uniform.log_prior_from_value(value=11.0) == float("-inf") - # The normalisation constant being dropped means the returned values - # do NOT depend on the (lower_limit, upper_limit) pair — only on `value`. - log_uniform = af.LogUniformPrior(lower_limit=50.0, upper_limit=100.0) + # Inside the support the dropped constant means the returned value does + # NOT depend on the (lower_limit, upper_limit) pair — only on `value`. + log_uniform = af.LogUniformPrior(lower_limit=0.5, upper_limit=100.0) assert log_uniform.log_prior_from_value(value=1.0) == 0.0 assert log_uniform.log_prior_from_value(value=2.0) == pytest.approx( @@ -216,6 +229,7 @@ def test__log_prior_from_value(self): assert log_uniform.log_prior_from_value(value=4.0) == pytest.approx( -np.log(4.0), 1.0e-12 ) + assert log_uniform.log_prior_from_value(value=0.4) == float("-inf") def test__log_prior_from_value__non_positive_returns_neg_inf(self): # Regression (PyAutoHeart #27 / release run 28784914443): Emcee's stretch @@ -232,14 +246,13 @@ def test__log_prior_from_value__non_positive_returns_neg_inf(self): assert log_uniform.log_prior_from_value(value=-1.0) == float("-inf") assert log_uniform.log_prior_from_value(value=0.0) == float("-inf") - # Positive values are unchanged: the NumPy path stays unnormalised and - # unbounded, returning -log(value) regardless of (lower_limit, upper_limit). + # In-support values are unchanged, returning the unnormalised -log(value); + # values above upper_limit are now -inf like the rest of the excluded + # support (PyAutoFit#1489). assert log_uniform.log_prior_from_value(value=10.0) == pytest.approx( -np.log(10.0), 1.0e-12 ) - assert log_uniform.log_prior_from_value(value=1.0e4) == pytest.approx( - -np.log(1.0e4), 1.0e-12 - ) + assert log_uniform.log_prior_from_value(value=1.0e4) == float("-inf") def test__lower_limit_zero_or_below_raises_error(self): with pytest.raises(exc.PriorException): diff --git a/test_autofit/mapper/prior/test_prior_bounds_1489.py b/test_autofit/mapper/prior/test_prior_bounds_1489.py new file mode 100644 index 000000000..ee6251754 --- /dev/null +++ b/test_autofit/mapper/prior/test_prior_bounds_1489.py @@ -0,0 +1,199 @@ +""" +Regression tests for PyAutoFit#1489: prior bounds must be enforced in the search +objective on the NumPy path. + +Between release 2025.10.16.1 (which removed ``assert_within_limits`` / +``PriorLimitException``) and this fix, ``UniformPrior.log_prior_from_value`` +returned ``0.0`` unconditionally on the NumPy path and ``LogUniformPrior`` +returned a finite ``-log(value)`` outside ``[lower_limit, upper_limit]``, so +nothing in the fitness path penalised an out-of-support parameter for the +searches that form a log posterior (Emcee, Zeus, Drawer, LBFGS). An Emcee walker +on a poorly-constrained parameter then escaped its declared box without bound. + +These tests pin the restored contract: ``-inf`` outside the support on the NumPy +path, matching the JAX path; a ``-inf`` figure of merit for an out-of-box +vector; a NaN-free figure-of-merit inversion; and end-to-end Emcee containment. +""" +import warnings + +import numpy as np +import pytest + +import autofit as af +from autofit.non_linear.fitness import Fitness + + +class TestUniformPriorBounds: + def test__scalar__zero_inside__neg_inf_outside(self): + prior = af.UniformPrior(lower_limit=-0.1, upper_limit=0.1) + + assert prior.log_prior_from_value(value=0.0) == 0.0 + assert prior.log_prior_from_value(value=0.05) == 0.0 + assert prior.log_prior_from_value(value=0.11) == -np.inf + assert prior.log_prior_from_value(value=-0.11) == -np.inf + assert prior.log_prior_from_value(value=5.0) == -np.inf + + def test__limits_are_inclusive(self): + prior = af.UniformPrior(lower_limit=-0.1, upper_limit=0.1) + + assert prior.log_prior_from_value(value=-0.1) == 0.0 + assert prior.log_prior_from_value(value=0.1) == 0.0 + + def test__array_input(self): + prior = af.UniformPrior(lower_limit=-0.1, upper_limit=0.1) + + log_prior = prior.log_prior_from_value(value=np.array([0.0, 0.1, 0.2, -5.0])) + + assert log_prior == pytest.approx(np.array([0.0, 0.0, -np.inf, -np.inf])) + + +class TestLogUniformPriorBounds: + def test__scalar__log_density_inside__neg_inf_outside(self): + prior = af.LogUniformPrior(lower_limit=1e-3, upper_limit=1e3) + + assert prior.log_prior_from_value(value=10.0) == pytest.approx( + -np.log(10.0), 1.0e-12 + ) + # Above the upper limit was previously a finite -log(value) — the gap + # that left the box unenforced from above. + assert prior.log_prior_from_value(value=1.0e4) == -np.inf + # Between zero and the lower limit was previously finite too. + assert prior.log_prior_from_value(value=1.0e-4) == -np.inf + + def test__non_positive_returns_neg_inf_without_warning(self): + prior = af.LogUniformPrior(lower_limit=1e-3, upper_limit=1e3) + + with warnings.catch_warnings(): + warnings.simplefilter("error", category=RuntimeWarning) + assert prior.log_prior_from_value(value=0.0) == -np.inf + assert prior.log_prior_from_value(value=-1.0) == -np.inf + + def test__array_input(self): + prior = af.LogUniformPrior(lower_limit=1e-3, upper_limit=1e3) + + with warnings.catch_warnings(): + warnings.simplefilter("error", category=RuntimeWarning) + log_prior = prior.log_prior_from_value( + value=np.array([10.0, 1.0e4, 1.0e-4, -1.0]) + ) + + assert log_prior == pytest.approx( + np.array([-np.log(10.0), -np.inf, -np.inf, -np.inf]) + ) + + +class GaussianWithOffset: + def __init__(self, centre=50.0, normalization=1.0, sigma=5.0, offset=0.0): + self.centre = centre + self.normalization = normalization + self.sigma = sigma + self.offset = offset # ignored by the likelihood: deliberately unconstrained + + def model_data_from(self, xvalues): + transformed = xvalues - self.centre + return ( + self.normalization + / (self.sigma * np.sqrt(2.0 * np.pi)) + * np.exp(-0.5 * (transformed / self.sigma) ** 2) + ) + + +class AnalysisIgnoringOffset(af.Analysis): + def __init__(self, data, noise_map): + super().__init__() + self.data = data + self.noise_map = noise_map + self.xvalues = np.arange(data.shape[0], dtype=float) + + def log_likelihood_function(self, instance): + model_data = instance.model_data_from(self.xvalues) + residual_map = self.data - model_data + chi_squared = float(np.sum((residual_map / self.noise_map) ** 2.0)) + noise_normalization = float(np.sum(np.log(2 * np.pi * self.noise_map**2.0))) + return -0.5 * (chi_squared + noise_normalization) + + +def _make_model_and_analysis(): + rng = np.random.default_rng(1) + xvalues = np.arange(100, dtype=float) + truth = GaussianWithOffset(centre=50.0, normalization=25.0, sigma=5.0) + noise = 0.1 + data = truth.model_data_from(xvalues) + rng.normal(0.0, noise, size=xvalues.shape) + noise_map = np.full(xvalues.shape, noise) + + model = af.Model(GaussianWithOffset) + model.centre = af.UniformPrior(lower_limit=0.0, upper_limit=100.0) + model.normalization = af.UniformPrior(lower_limit=1e-2, upper_limit=1e2) + model.sigma = af.UniformPrior(lower_limit=0.1, upper_limit=25.0) + model.offset = af.UniformPrior(lower_limit=-0.1, upper_limit=0.1) + + return model, AnalysisIgnoringOffset(data=data, noise_map=noise_map) + + +class TestFitnessObjective: + def test__out_of_box_vector__log_posterior_is_neg_inf(self): + model, analysis = _make_model_and_analysis() + + fitness = Fitness( + model=model, + analysis=analysis, + paths=None, + fom_is_log_likelihood=False, + ) + + in_bounds = [50.0, 25.0, 5.0, 0.05] + out_of_box = [50.0, 25.0, 5.0, 5.0] + + assert np.isfinite(fitness(parameters=in_bounds)) + assert fitness(parameters=out_of_box) == -np.inf + + def test__log_likelihood_from__out_of_box_is_not_nan(self): + # -inf figure of merit minus -inf log prior is NaN unless guarded; the + # inversion feeds history / quick-update / resume bookkeeping, which + # compare against it. + model, analysis = _make_model_and_analysis() + + fitness = Fitness( + model=model, + analysis=analysis, + paths=None, + fom_is_log_likelihood=False, + ) + + out_of_box = [50.0, 25.0, 5.0, 5.0] + figure_of_merit = fitness(parameters=out_of_box) + + log_likelihood = fitness.log_likelihood_from( + figure_of_merit=figure_of_merit, parameters=out_of_box + ) + + assert not np.isnan(log_likelihood) + assert log_likelihood == -np.inf + + +class TestEmceeContainment: + def test__unconstrained_parameter_stays_inside_declared_box(self): + """ + End-to-end: with the strict NumPy-path priors every out-of-box proposal + has a -inf log posterior, so Emcee rejects it and — because walkers are + initialised inside the box — containment is guaranteed, independent of + the random seed. Before the fix, the unconstrained ``offset`` escaped + ``UniformPrior(-0.1, 0.1)`` in 100% of accepted samples (PyAutoFit#1489). + """ + model, analysis = _make_model_and_analysis() + + search = af.Emcee( + name="test_uniform_bounds_containment_1489", + nwalkers=10, + nsteps=300, + iterations_per_full_update=1_000_000, + number_of_cores=1, + ) + + result = search.fit(model=model, analysis=analysis) + + offset_samples = np.asarray(result.samples.parameter_lists)[:, 3] + + assert offset_samples.size > 0 + assert np.all(offset_samples >= -0.1) + assert np.all(offset_samples <= 0.1) diff --git a/test_autofit/non_linear/test_clipper.py b/test_autofit/non_linear/test_clipper.py index 8404fb7f3..850936036 100644 --- a/test_autofit/non_linear/test_clipper.py +++ b/test_autofit/non_linear/test_clipper.py @@ -368,8 +368,13 @@ def model(): assert clipped_centre <= 20.0 assert not np.isnan(clipped_centre) - # The unconstrained fit is free to chase the truth past the prior edge. - assert unbounded_centre > 20.0 + # Since PyAutoFit#1489 the strict NumPy-path priors make the objective + # -inf outside the box, so even the clipper-less fit can no longer chase + # the truth past the prior edge — it hits an infinite wall instead of + # escaping. (This test previously asserted `unbounded_centre > 20.0`, + # pinning the pre-#1489 escape as the foil for the clipped run.) + assert unbounded_centre <= 20.0 + assert not np.isnan(unbounded_centre) def test__non_bound_supporting_method_raises_rather_than_being_ignored(self): """