diff --git a/autofit/mapper/prior/abstract.py b/autofit/mapper/prior/abstract.py index d0a2721c2..0de07a3c1 100644 --- a/autofit/mapper/prior/abstract.py +++ b/autofit/mapper/prior/abstract.py @@ -37,6 +37,27 @@ class attribute. self.width_modifier = None + def log_normalisation(self, xp=np) -> float: + """ + The additive constant dropped from :meth:`log_prior_from_value`. + + By contract ``log_prior_from_value`` returns the log prior density *up to + an additive constant*: the value-independent normaliser is dropped because + it is irrelevant to posterior shape (it cancels in the Metropolis ratio, + and nested samplers use the unit-cube transform rather than the prior + density). This is harmless for sampling but means absolute ``log_prior`` + values are not comparable across prior types. + + For code that needs the *fully normalised* log density — e.g. evidence or + Bayes-factor arithmetic — the normalised value is:: + + log_prior_from_value(value) + log_normalisation() + + The default is ``0.0`` (constant unknown / already normalised). Priors with + a known closed-form normaliser override this. + """ + return 0.0 + @classmethod def tree_unflatten(cls, aux_data, children): """ diff --git a/autofit/mapper/prior/gaussian.py b/autofit/mapper/prior/gaussian.py index fd536c147..7ed2f3640 100644 --- a/autofit/mapper/prior/gaussian.py +++ b/autofit/mapper/prior/gaussian.py @@ -114,6 +114,12 @@ def parameter_string(self) -> str: """ return f"mean = {self.mean}, sigma = {self.sigma}" + def log_normalisation(self, xp=np) -> float: + """The constant ``-log(sigma) - 0.5*log(2*pi)`` dropped from the density-form + quadratic returned by ``NormalMessage.log_prior_from_value``. See + ``Prior.log_normalisation``.""" + return -xp.log(self.sigma) - 0.5 * xp.log(2.0 * np.pi) + def value_for(self, unit, xp=np): """ Map a unit value in [0, 1] to a physical value drawn from this Gaussian prior. diff --git a/autofit/mapper/prior/log_gaussian.py b/autofit/mapper/prior/log_gaussian.py index 8e83e80c1..d244faac8 100644 --- a/autofit/mapper/prior/log_gaussian.py +++ b/autofit/mapper/prior/log_gaussian.py @@ -95,8 +95,6 @@ def with_limits(cls, lower_limit: float, upper_limit: float) -> "LogGaussianPrio return cls( mean=(lower_limit + upper_limit) / 2, sigma=upper_limit - lower_limit, - lower_limit=lower_limit, - upper_limit=upper_limit, ) def _new_for_base_message(self, message): @@ -108,9 +106,7 @@ def _new_for_base_message(self, message): """ return LogGaussianPrior( *message.parameters, - lower_limit=self.lower_limit, - upper_limit=self.upper_limit, - id_=self.instance().id, + id_=self.id, ) def value_for(self, unit, xp=np): @@ -141,6 +137,12 @@ def value_for(self, unit, xp=np): def parameter_string(self) -> str: return f"mean = {self.mean}, sigma = {self.sigma}" + def log_normalisation(self, xp=np) -> float: + """The constant ``-log(sigma) - 0.5*log(2*pi)`` dropped from the Gaussian-in-log + density in ``log_prior_from_value`` (the value-dependent ``-log(value)`` + change-of-variables Jacobian is kept). See ``Prior.log_normalisation``.""" + return -xp.log(self.sigma) - 0.5 * xp.log(2.0 * np.pi) + def log_prior_from_value(self, value, xp=np): """ Compute the log prior density of a given physical value under this log-Gaussian prior. diff --git a/autofit/mapper/prior/log_uniform.py b/autofit/mapper/prior/log_uniform.py index 7d0725468..06bdcdebf 100644 --- a/autofit/mapper/prior/log_uniform.py +++ b/autofit/mapper/prior/log_uniform.py @@ -153,6 +153,12 @@ def log_prior_from_value(self, value, xp=np): in_bounds = (value >= self.lower_limit) & (value <= self.upper_limit) return xp.where(in_bounds, -xp.log(value), -xp.inf) + def log_normalisation(self, xp=np) -> float: + """The constant ``-log(log(upper / lower))`` dropped from + ``log_prior_from_value`` (which returns ``-log(value)``). See + ``Prior.log_normalisation``.""" + return -xp.log(xp.log(self.upper_limit / self.lower_limit)) + def value_for(self, unit, xp=np): """ Returns a physical value from an input unit value according to the limits of the log10 uniform prior. diff --git a/autofit/mapper/prior/truncated_gaussian.py b/autofit/mapper/prior/truncated_gaussian.py index 0728bde47..20a07ce9a 100644 --- a/autofit/mapper/prior/truncated_gaussian.py +++ b/autofit/mapper/prior/truncated_gaussian.py @@ -133,6 +133,21 @@ def parameter_string(self) -> str: f"upper_limit = {self.upper_limit}" ) + def log_normalisation(self, xp=np) -> float: + """The constant ``-log(sigma) - 0.5*log(2*pi) - log(Z)`` dropped from + ``TruncatedNormalMessage.log_prior_from_value``, where + ``Z = Phi((upper - mean)/sigma) - Phi((lower - mean)/sigma)`` is the + truncation mass. See ``Prior.log_normalisation``.""" + if xp.__name__.startswith("jax"): + import jax.scipy.stats as jstats + norm = jstats.norm + else: + from scipy.stats import norm + a = (self.lower_limit - self.mean) / self.sigma + b = (self.upper_limit - self.mean) / self.sigma + Z = norm.cdf(b) - norm.cdf(a) + return -xp.log(self.sigma) - 0.5 * xp.log(2.0 * np.pi) - xp.log(Z) + def value_for(self, unit, xp=np): """ Map a unit value in [0, 1] to a physical value drawn from this truncated Gaussian prior. diff --git a/autofit/mapper/prior/uniform.py b/autofit/mapper/prior/uniform.py index 282303dbc..e2daa5397 100644 --- a/autofit/mapper/prior/uniform.py +++ b/autofit/mapper/prior/uniform.py @@ -106,13 +106,21 @@ def logpdf(self, x): Parameters ---------- x - The value at which to evaluate the log PDF. + The value(s) at which to evaluate the log PDF. May be a scalar or a + numpy array. """ - # TODO: handle x as a numpy array - if x == self.lower_limit: - x += epsilon - elif x == self.upper_limit: - x -= epsilon + # Nudge values sitting exactly on a boundary inwards by epsilon, where the + # PDF is otherwise undefined. The scalar path is kept bit-identical to the + # historical behaviour; the array path vectorises the same snap with + # ``np.where`` (the previous scalar-only ``==`` comparison raised on arrays). + if np.ndim(x) == 0: + if x == self.lower_limit: + x += epsilon + elif x == self.upper_limit: + x -= epsilon + else: + x = np.where(x == self.lower_limit, x + epsilon, x) + x = np.where(x == self.upper_limit, x - epsilon, x) return self.message.logpdf(x) def dict(self) -> dict: @@ -178,6 +186,11 @@ def log_prior_from_value(self, value, xp=np): in_bounds = (value >= self.lower_limit) & (value <= self.upper_limit) return xp.where(in_bounds, xp.zeros_like(value), -xp.inf) + def log_normalisation(self, xp=np) -> float: + """The constant ``-log(upper - lower)`` dropped from ``log_prior_from_value`` + (which returns ``0.0``). See ``Prior.log_normalisation``.""" + return -xp.log(self.upper_limit - self.lower_limit) + @property def limits(self) -> Tuple[float, float]: """The (lower_limit, upper_limit) bounds of this uniform prior.""" diff --git a/autofit/messages/beta.py b/autofit/messages/beta.py index 9b9e6b919..b489d2251 100644 --- a/autofit/messages/beta.py +++ b/autofit/messages/beta.py @@ -74,9 +74,12 @@ def inv_beta_suffstats( b Estimated beta parameter(s) of the Beta distribution. - Warnings - -------- - Emits a RuntimeWarning if negative parameters are found, and clamps them to 0.5. + Raises + ------ + ValueError + If the Newton-Raphson projection produces a negative alpha or beta, which + is not a valid Beta distribution. Previously this was silently warned and + clamped into a local that was immediately overwritten (a no-op). """ _lnX, _ln1X = np.ravel(lnX), np.ravel(ln1X) @@ -94,12 +97,14 @@ def inv_beta_suffstats( ab += np.linalg.solve(jac, - f) if np.any(ab < 0): - warnings.warn( - "invalid negative parameters found for inv_beta_suffstats, " - "clampling value to 0.5", - RuntimeWarning + raise ValueError( + "inv_beta_suffstats produced negative Beta parameters " + f"(alpha, beta):\n\n{ab}\n\n" + "A negative alpha or beta is not a valid Beta distribution, so the " + "moment-matching projection has failed. The previous behaviour " + "silently warned and clamped into a local that was immediately " + "overwritten (a no-op), letting the invalid parameters escape." ) - b = np.clip(ab, 0.5, None) shape = np.shape(lnX) if shape: diff --git a/autofit/messages/fixed.py b/autofit/messages/fixed.py index 1280d5e95..6e672847b 100644 --- a/autofit/messages/fixed.py +++ b/autofit/messages/fixed.py @@ -54,12 +54,12 @@ def sample(self, n_samples: Optional[int] = None) -> np.ndarray: return self.value return np.array([self.value]) - logpdf_cache = {} - def logpdf(self, x: np.ndarray) -> np.ndarray: - if x.shape not in FixedMessage.logpdf_cache: - FixedMessage.logpdf_cache[x.shape] = np.zeros_like(x) - return FixedMessage.logpdf_cache[x.shape] + # A fixed message contributes zero log-density everywhere. Return a + # fresh zero array each call: the previous class-level ``logpdf_cache`` + # was an unbounded dict keyed on shape that also handed back an aliased + # mutable array, so mutating one result silently corrupted later calls. + return np.zeros_like(x) @cached_property def mean(self) -> np.ndarray: diff --git a/autofit/messages/gamma.py b/autofit/messages/gamma.py index 3677186f0..c744571ba 100644 --- a/autofit/messages/gamma.py +++ b/autofit/messages/gamma.py @@ -76,8 +76,13 @@ def sample(self, n_samples=None): def from_mode(cls, mode, covariance, **kwargs): m, V = cls._get_mean_variance(mode, covariance) - alpha = 1 + m ** 2 * V # match variance - beta = alpha / m # match mean + # Match the requested mean m and variance V, consistent with the Normal + # family. For a Gamma, mean = alpha / beta and variance = alpha / beta**2, + # which inverts to alpha = m**2 / V, beta = m / V. The previous + # ``alpha = 1 + m**2 * V`` had the variance the wrong way up (requesting + # variance 0.25 produced 2.0, and 4.0 produced 0.235). + alpha = m ** 2 / V + beta = m / V return cls(alpha, beta, **kwargs) def kl(self, dist): diff --git a/autofit/messages/truncated_normal.py b/autofit/messages/truncated_normal.py index 812f80384..9915e1bff 100644 --- a/autofit/messages/truncated_normal.py +++ b/autofit/messages/truncated_normal.py @@ -30,11 +30,20 @@ def log_partition(self, xp=np) -> float: """ Compute the log-partition function (normalizer) of the truncated Gaussian. - This is the log of the normalization constant Z of the truncated normal: + For the exponential-family interface this is the full cumulant, made of + two parts: - Z = Φ((b - μ)/σ) - Φ((a - μ)/σ) + 1. the untruncated Gaussian log-partition ``A(η) = μ²/(2σ²) + log σ`` + (identical to ``NormalMessage.log_partition``), and + 2. the truncation-mass correction ``log Z`` with + ``Z = Φ((b - μ)/σ) - Φ((a - μ)/σ)``, where Φ is the standard normal + CDF and ``[a, b]`` are the truncation bounds. - where Φ is the standard normal CDF and [a, b] are the truncation bounds. + Previously only ``log Z`` was returned, dropping the Gaussian term. That + made the generic exponential-family pdf integrate to + ``σ·exp(μ²/2σ²)`` (e.g. 2.27 for the unit case) instead of 1.0 — the path + the EP machinery consumes. Sampling and ``log_prior_from_value`` use a + separate, correct path and were unaffected. Returns ------- @@ -43,10 +52,14 @@ def log_partition(self, xp=np) -> float: """ from scipy.stats import norm + # Untruncated Gaussian log-partition — see NormalMessage.log_partition. + gaussian = (self.mean ** 2) / (2 * self.sigma ** 2) + xp.log(self.sigma) + a = (self.lower_limit - self.mean) / self.sigma b = (self.upper_limit - self.mean) / self.sigma Z = norm.cdf(b) - norm.cdf(a) - return xp.log(Z) if Z > 0 else -xp.inf + log_Z = xp.log(Z) if Z > 0 else -xp.inf + return gaussian + log_Z log_base_measure = -0.5 * np.log(2 * np.pi) @@ -472,9 +485,19 @@ def log_prior_from_value(self, value: float, xp=np) -> float: """ Compute the log prior probability of a given physical value under this truncated Gaussian prior. - This accounts for truncation by normalizing the Gaussian density over the - interval [lower_limit, upper_limit], returning -inf if the value lies outside - these limits. + Returns ``log p(value)`` in density form, up to an additive constant, and + ``-inf`` for values outside ``[lower_limit, upper_limit]``. + + The value-independent constants ``-log(sigma) - 0.5*log(2*pi)`` (the + Gaussian normaliser) and ``-log(Z)`` (the truncation mass) are dropped, so + this matches the constant-dropping convention already used by + ``NormalMessage.log_prior_from_value`` (and Uniform / LogUniform / LogGaussian). + Previously this method returned the *fully normalised* truncated density, + making it the odd one out — harmless to posterior shape (constants cancel + in the Metropolis ratio and nested samplers use the unit-cube transform), + but inconsistent for anyone reading absolute ``log_prior`` values or doing + evidence arithmetic. The dropped constant is recoverable via + ``TruncatedGaussianPrior.log_normalisation``. Parameters ---------- @@ -483,33 +506,18 @@ def log_prior_from_value(self, value: float, xp=np) -> float: Returns ------- - The log prior probability of the given value, or -inf if outside truncation bounds. + The log prior density at the given value up to an additive constant, or + -inf if outside the truncation bounds. """ - if xp.__name__.startswith("jax"): - import jax.scipy.stats as jstats - norm = jstats.norm - else: - from scipy.stats import norm - - # Normalization term (truncation) - a = (self.lower_limit - self.mean) / self.sigma - b = (self.upper_limit - self.mean) / self.sigma - Z = norm.cdf(b) - norm.cdf(a) - - # Log pdf + # Density-form quadratic, constants dropped (see docstring / NormalMessage). z = (value - self.mean) / self.sigma - log_pdf = ( - -0.5 * z ** 2 - - xp.log(self.sigma) - - 0.5 * xp.log(2.0 * xp.pi) - ) - log_trunc_pdf = log_pdf - xp.log(Z) + log_pdf = -0.5 * z ** 2 - # Truncation mask (must be xp.where for JAX) + # Truncation mask (must be xp.where for JAX). in_bounds = (self.lower_limit <= value) & (value <= self.upper_limit) - return xp.where(in_bounds, log_trunc_pdf, -xp.inf) + return xp.where(in_bounds, log_pdf, -xp.inf) def __str__(self): """ diff --git a/test_autofit/mapper/prior/test_priors_messages_fixes_1331.py b/test_autofit/mapper/prior/test_priors_messages_fixes_1331.py new file mode 100644 index 000000000..782c5cd23 --- /dev/null +++ b/test_autofit/mapper/prior/test_priors_messages_fixes_1331.py @@ -0,0 +1,133 @@ +"""Regression tests for the priors/messages correctness batch (PyAutoFit #1344, +decision hub #1331; per-bug detail #1330). Numpy-only — no JAX in library unit +tests. +""" +import numpy as np +import pytest +from scipy import integrate +from scipy.stats import norm + +import autofit as af +from autofit.messages.fixed import FixedMessage +from autofit.messages.gamma import GammaMessage +from autofit.messages.truncated_normal import TruncatedNormalMessage +from autofit.messages.beta import inv_beta_suffstats + + +# --- #01: LogGaussianPrior.with_limits / _new_for_base_message no longer crash --- + +def test__log_gaussian_with_limits_no_crash(): + # Previously raised TypeError: the ctor does not accept lower/upper_limit. + prior = af.LogGaussianPrior.with_limits(0.5, 1.5) + assert prior.mean == pytest.approx(1.0) # (0.5 + 1.5) / 2 + assert prior.sigma == pytest.approx(1.0) # 1.5 - 0.5 + # _new_for_base_message (called during projection with the base NormalMessage) + # was also broken — the invalid kwargs plus a stale self.instance().id call. + remade = prior._new_for_base_message(prior.message.base_message) + assert isinstance(remade, af.LogGaussianPrior) + assert remade.id == prior.id + + +# --- #02: UniformPrior.logpdf handles array input; scalar path bit-identical --- + +def test__uniform_logpdf_array_matches_scalar(): + prior = af.UniformPrior(lower_limit=0.0, upper_limit=2.0) + xs = np.array([0.0, 0.5, 1.0, 1.5, 2.0]) + arr = np.asarray(prior.logpdf(xs)) + assert arr.shape == xs.shape + for i, x in enumerate(xs): + # scalar path must equal the corresponding array element exactly + assert float(prior.logpdf(float(x))) == float(arr[i]) + + +# --- #04: TruncatedNormalMessage.log_partition includes the Gaussian term, so the +# generic exponential-family pdf integrates to 1.0 (was sigma*exp(mu^2/2sigma^2)) --- + +def test__truncated_message_log_partition_normalises_pdf(): + # mu=1, sigma=2 -> old error factor sigma*exp(mu^2/2sigma^2) = 2*exp(1/8) ~= 2.27. + m = TruncatedNormalMessage(mean=1.0, sigma=2.0, lower_limit=-20.0, upper_limit=20.0) + integral, _ = integrate.quad(lambda x: float(np.exp(m.logpdf(np.array(x)))), -20.0, 20.0) + assert integral == pytest.approx(1.0, rel=1e-4) + + # log_partition = untruncated Gaussian cumulant + log(truncation mass). + a = (m.lower_limit - m.mean) / m.sigma + b = (m.upper_limit - m.mean) / m.sigma + expected = (m.mean ** 2) / (2 * m.sigma ** 2) + np.log(m.sigma) + np.log(norm.cdf(b) - norm.cdf(a)) + assert float(m.log_partition()) == pytest.approx(expected, rel=1e-9) + + +# --- #10: FixedMessage.logpdf returns a fresh array (no aliased class-level cache) --- + +def test__fixed_message_logpdf_not_aliased(): + assert not hasattr(FixedMessage, "logpdf_cache") + m = FixedMessage(np.array([1.0, 2.0, 3.0])) + x = np.zeros(3) + first = m.logpdf(x) + first[0] = 999.0 # mutating one result must not corrupt the next + second = m.logpdf(x) + assert np.all(second == 0.0) + assert first is not second + + +# --- Decision 1: inv_beta_suffstats raises on negative projection (was a no-op clamp) --- + +def test__inv_beta_suffstats_raises_on_negative(monkeypatch): + # Force the Newton-Raphson step negative (per the #05 reproducer) and confirm + # the projection now raises instead of silently warning + no-op clamping. + from scipy.special import digamma + a0, b0 = np.array([2.0, 4.0, 3.0]), np.array([3.0, 1.5, 5.0]) + lnX = digamma(a0) - digamma(a0 + b0) + ln1X = digamma(b0) - digamma(a0 + b0) + monkeypatch.setattr(np.linalg, "solve", lambda J, f: np.full(np.shape(f), -1e3)) + with pytest.raises(ValueError): + inv_beta_suffstats(lnX, ln1X) + + +def test__inv_beta_suffstats_positive_projection_does_not_raise(monkeypatch): + # The guard only fires on a negative projection: a positive Newton step must + # pass through unchanged (this is the in-scope behaviour Decision 1 changes — + # the underlying solve shape-handling is pre-existing and untouched here). + from scipy.special import digamma + a0, b0 = np.array([2.0, 4.0, 3.0]), np.array([3.0, 1.5, 5.0]) + lnX = digamma(a0) - digamma(a0 + b0) + ln1X = digamma(b0) - digamma(a0 + b0) + monkeypatch.setattr( + np.linalg, "solve", lambda J, f: np.zeros(np.shape(f)) + ) # no Newton movement -> ab stays at its positive initial guess + a, b = inv_beta_suffstats(lnX, ln1X) + assert np.all(np.asarray(a) > 0) and np.all(np.asarray(b) > 0) + + +# --- Decision 3: GammaMessage.from_mode matches the requested mean and variance --- + +@pytest.mark.parametrize("mean, variance", [(2.0, 0.25), (1.0, 4.0), (5.0, 0.5)]) +def test__gamma_from_mode_matches_mean_variance(mean, variance): + g = GammaMessage.from_mode(mean, variance) + assert float(g.mean) == pytest.approx(mean, rel=1e-9) + assert float(g.variance) == pytest.approx(variance, rel=1e-9) + + +# --- Decision 4: log_normalisation recovers the fully normalised density from the +# constant-dropping log_prior_from_value, for every prior with a known normaliser --- + +def test__uniform_log_normalisation(): + prior = af.UniformPrior(lower_limit=0.0, upper_limit=2.0) + # normalised density of a uniform on [0, 2] is log(1/2) + recovered = prior.log_prior_from_value(1.0) + prior.log_normalisation() + assert float(recovered) == pytest.approx(np.log(0.5), rel=1e-9) + + +def test__log_uniform_log_normalisation(): + prior = af.LogUniformPrior(lower_limit=1.0, upper_limit=100.0) + value = 10.0 + recovered = prior.log_prior_from_value(value) + prior.log_normalisation() + expected = -np.log(value) - np.log(np.log(100.0 / 1.0)) + assert float(recovered) == pytest.approx(expected, rel=1e-9) + + +def test__gaussian_log_normalisation(): + prior = af.GaussianPrior(mean=1.0, sigma=2.0) + value = 1.5 + recovered = prior.log_prior_from_value(value) + prior.log_normalisation() + expected = norm.logpdf(value, loc=1.0, scale=2.0) + assert float(recovered) == pytest.approx(expected, rel=1e-9) diff --git a/test_autofit/mapper/prior/test_truncated_gaussian.py b/test_autofit/mapper/prior/test_truncated_gaussian.py index 78ac12dd4..ce05a3935 100644 --- a/test_autofit/mapper/prior/test_truncated_gaussian.py +++ b/test_autofit/mapper/prior/test_truncated_gaussian.py @@ -26,16 +26,35 @@ def test__values(truncated_gaussian, unit, value): assert truncated_gaussian.value_for(unit) == pytest.approx(value, rel=0.1) @pytest.mark.parametrize( - "unit, value", + "value, expected", [ - (0.01, -np.inf), - (1.0, 2.3026892553), - (2.0, -np.inf), + (0.01, -np.inf), # below lower_limit -> out of support + (1.0, 0.0), # at the mode (== mean): bare quadratic, constants dropped + (2.0, -np.inf), # above upper_limit -> out of support ], ) -def test__log_prior_from_value(truncated_gaussian, unit, value): - - assert truncated_gaussian.log_prior_from_value(unit) == pytest.approx(value, rel=0.1) +def test__log_prior_from_value(truncated_gaussian, value, expected): + # Drop-constants contract (#1331 Decision 4): log_prior_from_value returns the + # density-form quadratic -0.5*((x-mean)/sigma)**2 up to an additive constant, + # and -inf outside [lower_limit, upper_limit]. Previously this method was the + # odd one out, returning the fully normalised truncated density (2.30 at the + # mode); now it matches NormalMessage / Uniform / LogUniform / LogGaussian. + assert truncated_gaussian.log_prior_from_value(value) == pytest.approx(expected, abs=1e-9) + + +def test__log_normalisation_recovers_normalised_density(truncated_gaussian): + # The dropped constant is recoverable via log_normalisation: adding it back to + # log_prior_from_value must reproduce the fully normalised truncated-normal + # log-density (scipy.stats.truncnorm.logpdf). This is the evidence-user path. + mean, sigma, lower, upper = 1.0, 2.0, 0.95, 1.05 + a, b = (lower - mean) / sigma, (upper - mean) / sigma + for value in (0.96, 1.0, 1.04): + recovered = ( + truncated_gaussian.log_prior_from_value(value) + + truncated_gaussian.log_normalisation() + ) + expected = truncnorm.logpdf(value, a=a, b=b, loc=mean, scale=sigma) + assert float(recovered) == pytest.approx(float(expected), rel=1e-9) # --- Numerical equivalence: new direct-ndtr path vs the OLD scipy.stats.norm