From 07072db45ce135a9e2c9edf7d1b4e019deff9045 Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Wed, 8 Jul 2026 17:27:23 +0100 Subject: [PATCH 1/2] docs: formal EP specification for autofit.graphical (README + statistical docstrings) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - autofit/graphical/README.md (new): the mathematics as implemented — factor-graph model, exponential-family mean field, the EP step (cavity / tilted / moment-matching projection / damped division as natural-parameter EMA), the three factor-optimiser paths, convergence criterion + KL direction contract, evidence decomposition, deterministic-variable mechanisms, parallel-EP semantics; each equation anchored to its implementing class/method. - MeanField.update_factor_mean_field: damped-update equations, delta semantics, invalid-projection fallback. - MeanField.kl: KL direction contract (m.kl(other) == KL(m||other)). - FactorApproximation: fix damped-update docstring error ((q_a^f)^(1-d) -> (q_a)^(1-d)). - composed_transform.py: module docstring stating the transform composition-order convention (innermost-first storage; reversal asymmetry) with a worked UniformPrior trace; direction notes on _transform/_inverse_transform; validity condition on TransformedMessage.kl. - AbstractMessage.project: moment-matching projection equations. EP framework review phase 2 (#1333); audit findings context in #1332. Co-Authored-By: Claude Fable 5 --- autofit/graphical/README.md | 228 +++++++++++++++++++++++++ autofit/graphical/mean_field.py | 54 +++++- autofit/messages/abstract.py | 23 ++- autofit/messages/composed_transform.py | 68 +++++++- 4 files changed, 362 insertions(+), 11 deletions(-) create mode 100644 autofit/graphical/README.md diff --git a/autofit/graphical/README.md b/autofit/graphical/README.md new file mode 100644 index 000000000..c833b1a8f --- /dev/null +++ b/autofit/graphical/README.md @@ -0,0 +1,228 @@ +# `autofit.graphical` — the mathematics, exactly as implemented + +This package fits **factor graphs** — joint models over many datasets +and/or shared parameters — using **expectation propagation (EP)**. This +document states the statistical machinery in formal equations, with each +equation anchored to the class or method that implements it, so that the +code can be verified against the math line-by-line (by a human or an AI +agent). It documents what the code *does*; design intent beyond that is +marked explicitly. + +Notation: `a, b` index factors; `i, j` index variables; `η` are natural +parameters; `T(x)` sufficient statistics; `A(η)` the log-partition +function; `h(x)` the base measure. + +--- + +## 1. The model + +A factor graph over variables `x = (x₁, …, x_n)`: + + p(x) ∝ ∏ₐ fₐ(xₐ) (1) + +where `xₐ` is the subset of variables factor `a` touches. Factors are +`Factor` objects (`factor_graphs/factor.py`) wrapping arbitrary +callables; a `FactorGraph` (`factor_graphs/graph.py`) is their product. +Priors participate as factors of their own (`declarative/factor/prior.py`, +`MeanField.from_priors`): a prior is one more term in Eq. (1). + +In the declarative interface (`declarative/`), each dataset contributes +an `AnalysisFactor` whose value is `Analysis.log_likelihood_function`, +and `FactorGraphModel` assembles Eq. (1) from those plus the prior +factors. + +## 2. The approximating family + +The approximation is fully factorised over both factors and variables +("mean field"): + + q(x) = ∏ₐ qₐ(x), qₐ(x) = ∏ᵢ qₐᵢ(xᵢ) (2) + +- `qₐᵢ` — one **message** per (factor, variable): an exponential-family + distribution (`autofit/messages/`), + + qₐᵢ(x) = h(x) exp( ηₐᵢ · T(x) − A(ηₐᵢ) ) (3) + +- `MeanField` (`mean_field.py`) is one factor's `{Variable → message}` + dictionary, i.e. one `qₐ`. +- `EPMeanField` (`expectation_propagation/ep_mean_field.py`) is the full + `{Factor → MeanField}` map, i.e. `q(x)`. + +Because messages are exponential-family, products and quotients are +sums and differences of natural parameters +(`MessageInterface.sum_natural_parameters`, message `__mul__` / +`__truediv__`; powers scale them, `AbstractMessage.__pow__`): + + ∏ₖ qₖ(x) ∝ h(x) exp( (Σₖ ηₖ) · T(x) ) (4) + +The global approximation to the posterior of `xᵢ` is the product of +every factor's message on it: `q(xᵢ) ∝ ∏ₐ qₐᵢ(xᵢ)` +(`EPMeanField.mean_field`). + +Non-Gaussian priors (Uniform, LogUniform, LogGaussian) are represented +as a base `NormalMessage` composed with deterministic transforms +(`TransformedMessage`, `messages/composed_transform.py`); the message +algebra operates on the base-space natural parameters and the transform +stack maps to physical space. See the module docstring of +`composed_transform.py` for the composition-order convention. + +## 3. One EP update + +For each factor `a` in turn (`EPOptimiser.run`, +`expectation_propagation/optimiser.py`): + +### 3.1 Cavity — `EPMeanField.factor_approximation` + +Everything except factor `a`: + + q^{\a}(x) = ∏_{b ≠ a} q_b(x) η_cav,i = Σ_{b≠a} η_{bi} (5) + +The code builds this as `MeanField({v: 1.0}).prod(*other_factor_fields)` +and packages `(factor, cavity_dist, factor_dist, model_dist)` into a +`FactorApproximation` (`mean_field.py`), where +`model_dist = factor_dist × cavity_dist` is the current full +approximation restricted to factor `a`'s variables. + +### 3.2 Tilted distribution — the factor optimiser + +The exact factor times the cavity: + + p̂ₐ(x) = fₐ(x) q^{\a}(x) / Ẑₐ , Ẑₐ = ∫ fₐ(x) q^{\a}(x) dx (6) + +`FactorApproximation.__call__` evaluates `log fₐ + log q^{\a}`. The +tilted distribution is then *fitted* by the factor's optimiser: + +- **Sampling path** (`AbstractSearch.optimise`, + `non_linear/search/abstract_search.py`): the cavity messages are + installed as the model's priors (`prior.with_message`), a non-linear + search (Dynesty, Nautilus, Emcee, …) samples Eq. (6), and the result + is projected per §3.3. +- **Laplace path** (`LaplaceOptimiser`, `graphical/laplace/`): + quasi-Newton ascent on `log p̂ₐ` (`newton.py`: BFGS/SR1 with Wolfe + line search, `line_search.py`), then a Gaussian at the optimum with + covariance from the (approximate) Hessian inverse + (`MeanField.from_mode_covariance` via `from_mode` on each message). +- **Exact path** (`ExactFactorFit`, + `expectation_propagation/factor_optimiser.py`): if the factor is + itself a message of the same family as the cavity + (`Factor.has_exact_projection`), the tilted distribution is conjugate + and the update is the closed-form natural-parameter sum — no sampler. + `EPOptimiser.from_meanfield` auto-selects this path. + +### 3.3 Projection (moment matching) — `AbstractMessage.project` + +Find the family member closest to the tilted distribution in inclusive +KL: + + q* = argmin_q KL( p̂ₐ ‖ q ) (7) + +For an exponential family this is **moment matching**: + + E_{q*}[ T(x) ] = E_{p̂ₐ}[ T(x) ] (8) + +Implemented as importance-weighted sample moments over the search's +weighted samples `{x_s, log w_s}`: + + E[T] ≈ Σ_s w̃_s T(x_s) / Σ_s w̃_s , w̃_s = exp(log w_s − max log w) (9) + +(`AbstractMessage.project`; the max-log-weight shift is the standard +numerical stabilisation, and the mean weight supplies the projection's +`log_norm`). `from_sufficient_statistics` then inverts Eq. (8) to +natural parameters per family. On the Laplace path the "projection" is +the Gaussian mode/covariance construction instead. + +### 3.4 Factor update with damping — `MeanField.update_factor_mean_field` + +Divide out the cavity and damp with `δ ∈ (0, 1]`: + + qₐ^new = (q*)^δ (qₐ^old)^{1−δ} / (q^{\a})^δ (10) + +equivalently, an exponential moving average on natural parameters: + + ηₐ ← (1 − δ) ηₐ + δ ( η_{q*} − η_cav ) (11) + +`δ` is supplied by an `ApproxUpdater` (`optimiser.py`): `SimplerUpdater` +(one fixed value), `FactorUpdater` (per factor), or `DynamicUpdater` +(per variable, `δᵢ ∝ min_count / count(i)` — variables shared by more +factors update more slowly). `δ` may therefore be a scalar or a +per-variable `MeanField` of scalars. + +**Invalid-projection fallback**: if the division produces an invalid +message (e.g. negative variance — possible because Eq. (10)'s +subtraction of natural parameters is not closed in the family), +`update_factor_mean_field` reverts the invalid parameters to the +previous message per-parameter (`update_invalid`) and flags +`StatusFlag.BAD_PROJECTION`. + +## 4. Convergence — `EPHistory` (`expectation_propagation/history.py`) + +After each factor update the history records the new `EPMeanField`. +Termination when either: + + KL( q_t ‖ q_{t−1} ) = Σᵢ KL( q_t(xᵢ) ‖ q_{t−1}(xᵢ) ) < kl_tol (12) + +(`FactorHistory.kl_divergence`, summed per-variable via `MeanField.kl`; +default `kl_tol = 1e-1`), or the log-evidence change drops below +`evidence_tol`, or a user callback fires. + +**KL direction contract**: `m.kl(other)` means `KL(m ‖ other)`. (As of +the 2026-07 audit, `NormalMessage` and `TruncatedNormalMessage` satisfy +this; `GammaMessage` and `BetaMessage` compute the reverse direction, +and `TruncatedNormalMessage` uses the untruncated formula — tracked in +issue #1332, findings F2/F6.) + +## 5. Evidence — `EPMeanField.log_evidence` + +The EP approximation to the model evidence factorises as: + + Zᵢ = ∫ ∏ₐ qₐᵢ(xᵢ) dxᵢ (per variable) (13) + Zₐ = Ẑₐ / ∏_{i ∈ a} Zᵢ (per factor) (14) + Z = ∏ᵢ Zᵢ ∏ₐ Zₐ (15) + +`Zᵢ` is computed in closed form from log-partitions +(`MessageInterface.log_normalisation`): + + log ∫ ∏ₖ qₖ = Σₖ (log hₖ − A(ηₖ)) − ( log h − A(Σₖ ηₖ) ) (16) + +and the per-factor `Ẑₐ` is carried on `MeanField.log_norm` by the +projection. (Audit note: as of 2026-07 the `log_norm` bookkeeping is +broken in three places — issue #1332 finding F7 — so Eq. (15) should +not be used for model comparison until those fixes land; the +decomposition itself is correct.) + +## 6. Deterministic variables + +Three composition mechanisms exist; they are **not** interchangeable +and reconciling them is an open design item (EP review Phase 5): + +1. **Graph-level deterministic variables**: `Factor(..., factor_out=v)` + declares outputs computed by the factor + (`FactorValue.deterministic_values`); the Laplace path maintains a + separate curvature estimate for them + (`quasi_deterministic_update`, `laplace/newton.py`), and their + approximate distributions live in the cavity + (`FactorApproximation.deterministic_dist`). +2. **Compound prior arithmetic** (`mapper/prior/arithmetic/`): + deterministic relations expressed on priors at model-composition + time (e.g. `prior_a * matrix + prior_b`). +3. **Free shared variables**: share a prior across factors and encode + the relation inside the likelihood. + +## 7. Parallel EP — `ParallelEPOptimiser` + +All factor approximations for a sweep are built from the **same** +mean field, factor optimisations run in a process pool, and the updates +are applied sequentially afterwards. This is standard "parallel EP" +semantics: cavities within a sweep are stale relative to serial EP, so +serial and parallel runs converge along different trajectories (to the +same fixed points when EP converges). + +## 8. Reading + +- T. Minka (2001), *Expectation Propagation for Approximate Bayesian + Inference* — the algorithm of §3. +- A. Vehtari et al. (2020), *Expectation propagation as a way of life* + (JMLR 21(17)) — the data-partitioned framing this package follows + (one factor per dataset, Eq. (1)). +- M. Seeger (2005), *Expectation propagation for exponential families* + — the natural-parameter algebra of §2–§3. diff --git a/autofit/graphical/mean_field.py b/autofit/graphical/mean_field.py index afd5071a3..d27bc91e8 100755 --- a/autofit/graphical/mean_field.py +++ b/autofit/graphical/mean_field.py @@ -388,6 +388,17 @@ def sample(self, n_samples=None): return VariableData({v: dist.sample(n_samples) for v, dist in self.items()}) def kl(self, mean_field: "MeanField") -> np.ndarray: + """ + The KL divergence between two mean fields, summed over variables: + + KL(self || mean_field) = Σᵢ KL( self[xᵢ] || mean_field[xᵢ] ) + + which is exact because both mean fields factorise over variables. + + Direction contract: ``message.kl(other)`` means ``KL(message || other)`` + for every message family. ``EPHistory`` uses this sum between + successive global mean fields as its convergence criterion. + """ return sum(np.sum(dist.kl(mean_field[k])) for k, dist in self.items()) __hash__ = Factor.__hash__ @@ -405,7 +416,44 @@ def update_factor_mean_field( delta: Delta = 1.0, status: Status = Status(), ) -> Tuple["MeanField", Status]: + """ + The EP factor update: given ``self`` = the newly-fitted model + distribution q* (the moment-matched projection of the tilted + distribution), divide out the cavity to recover the factor's new + message, damped by ``delta``: + + q_a_new = (q*)^δ · (q_a_old)^(1−δ) / (cavity)^δ + + which on natural parameters is an exponential moving average, + + η_a ← (1 − δ)·η_a + δ·(η_q* − η_cavity) + ``delta`` may be a scalar or a per-variable ``MeanField`` of scalars + (see ``DynamicUpdater``). At ``delta == 1`` this reduces to the + undamped ``q* / cavity``. + + The natural-parameter subtraction is not closed in the family (it can + produce e.g. a negative variance). If the result is invalid, the + invalid parameters are reverted per-parameter to ``last_dist`` via + ``update_invalid`` and the status is flagged ``BAD_PROJECTION``. + + Parameters + ---------- + cavity_dist + The cavity distribution: the product of every other factor's + messages on this factor's variables. + last_dist + The factor's previous message distribution ``q_a_old`` (required + when ``delta < 1`` and as the fallback for invalid projections). + delta + Damping coefficient in (0, 1]; scalar or per-variable MeanField. + status + Status carried through from the factor optimisation. + + Returns + ------- + The factor's new message distribution and the updated status. + """ success, messages, _, flag = status updated = False try: @@ -505,9 +553,9 @@ class FactorApproximation(AbstractNode): returns q⁺ᵃ(x₀), {xₐ: dq⁺ᵃ(x₀)/dxₐ} project_mean_field(mean_field, delta=1., status=Status()) - for qᶠ = mean_field, finds qₐ such that qᶠₐ * q⁻ᵃ = qᶠ - delta controls how much to change from the original factor qₐ - so qʳₐ = (qᶠₐ)ᵟ * (qᶠₐ)¹⁻ᵟ + for qᶠ = mean_field, finds qᶠₐ such that qᶠₐ * q⁻ᵃ = qᶠ + delta controls how much to change from the previous factor + distribution qₐ, so qʳₐ = (qᶠₐ)ᵟ * (qₐ)¹⁻ᵟ returns qʳₐ, status """ diff --git a/autofit/messages/abstract.py b/autofit/messages/abstract.py index e6a9fd860..de2b097ed 100644 --- a/autofit/messages/abstract.py +++ b/autofit/messages/abstract.py @@ -266,9 +266,26 @@ def factor(self, x): def project( cls, samples: np.ndarray, log_weight_list: Optional[np.ndarray] = None, **kwargs ) -> "AbstractMessage": - """Calculates the sufficient statistics of a set of samples - and returns the distribution with the appropriate parameters - that match the sufficient statistics + """ + Moment-matching projection: the exponential-family member closest + (in inclusive KL) to the weighted sample distribution. + + For an exponential family, argmin_q KL(p || q) is the member whose + expected sufficient statistics match p's: + + E_q[T(x)] = E_p[T(x)] + + Estimated by importance-weighted sample moments over weighted + samples {x_s, log w_s} (e.g. from a nested sampler run on the EP + tilted distribution): + + E[T] ≈ Σ_s w̃_s T(x_s) / Σ_s w̃_s , w̃_s = exp(log w_s − max log w) + + The max-log-weight shift is the standard numerical stabilisation; + the mean unshifted weight supplies the projection's ``log_norm`` + (the tilted distribution's normalisation estimate). + ``from_sufficient_statistics`` then inverts the moment equations to + natural parameters for the concrete family. """ # if weight_list aren't passed then equally weight all samples diff --git a/autofit/messages/composed_transform.py b/autofit/messages/composed_transform.py index 84792b062..4d42dcebf 100644 --- a/autofit/messages/composed_transform.py +++ b/autofit/messages/composed_transform.py @@ -1,3 +1,45 @@ +""" +Messages composed with deterministic transforms. + +A ``TransformedMessage`` represents a distribution over a *physical* +variable as a base exponential-family message over a *base* variable plus +a stack of invertible transforms. For example:: + + UniformPrior(0, 2) = TransformedMessage( + UniformNormalMessage, # base: N(0, 1) + Φ + LinearShiftTransform(shift=0, scale=2), # base → physical + ) + +Composition-order convention (the single most important fact in this +module): + +- ``transforms`` are stored **innermost-first**: going base → physical + they apply in tuple order (first element first), so the outermost + transform — the last function applied — is the *last* tuple element. +- ``_transform`` maps **physical → base** and therefore applies the + transforms in *reverse* tuple order, unwinding the outermost transform + (the last element) first. +- ``_inverse_transform`` maps **base → physical** and applies + ``inv_transform`` in *forward* tuple order, rebuilding the composition + from the inside out. + +Worked example for ``UniformPrior(0, 2).value_for(0.5)`` (base → +physical, forward order): + + 1. ``NormalMessage(0, 1).value_for(0.5)`` → 0.0 + 2. ``phi_transform.inv_transform(0.0)`` = Φ(0) → 0.5 + 3. ``LinearShiftTransform(0, 2).inv_transform(0.5)`` = 0.5·2 + 0 → 1.0 + +Evaluating ``logpdf(x)`` at a physical ``x`` runs the same three steps in +reverse (physical → base) via ``_transform``, accumulating the +log-Jacobian ``log_det`` of each transform for the change of variables. + +A new transform added to the stack must therefore implement +``transform`` (physical → base), ``inv_transform`` (base → physical) and +``log_det`` (the log-Jacobian of ``transform``); a ``LinearTransform``'s +stored operator is the Jacobian of the physical → base direction, which +is why ``LinearShiftTransform(scale=s)`` stores ``1/s``. +""" import functools import numpy as np import warnings @@ -155,6 +197,15 @@ def project( ) def kl(self, dist): + """ + KL divergence computed between the base messages. + + Valid because KL is invariant under a common invertible change of + variables — which requires ``dist`` to share this message's + transform stack (true inside EP, where both messages describe the + same variable). For messages with different transforms the result + is meaningless; no check is performed here. + """ return self.base_message.kl(dist.base_message) def natural_parameters(self, xp=np) -> np.ndarray: @@ -166,8 +217,13 @@ def sample(self, n_samples: Optional[int] = None): def _transform(self, x: float) -> float: """ - Transform some value in the space of the transformed message to - the space of the underlying message. + Map a value from *physical* space (the transformed message) to + *base* space (the underlying message). + + Applies the transforms in REVERSE tuple order: the stack is stored + innermost-first, so going physical -> base must unwind the outermost + transform (the last tuple element) first (see the module docstring + for the worked UniformPrior example). For example, a UniformPrior with limits 10 and 20 could be passed a value 15. If the underlying message is a NormalMessage with a @@ -188,10 +244,12 @@ def _transform(self, x: float) -> float: def _inverse_transform(self, x: float) -> float: """ - Transform some value in the space of the base message to a value in - the space of the transformed message. + Map a value from *base* space (the underlying message) to *physical* + space (the transformed message). Inverts ``_transform``. - Inverts transform (above) + Applies each transform's ``inv_transform`` in FORWARD tuple order, + rebuilding the composition from the inside out — the asymmetry with + ``_transform`` is deliberate and load-bearing (module docstring). """ for _transform in self.transforms: x = _transform.inv_transform(x) From 39eaf55ac93924f1d1e0343bbd67037d4dcbc1f8 Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Wed, 8 Jul 2026 19:59:20 +0100 Subject: [PATCH 2/2] =?UTF-8?q?docs/test:=20the=20EP=20lowering=20contract?= =?UTF-8?q?=20=E2=80=94=20seam=20documentation,=20rule,=20and=20tests?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - autofit/graphical/README.md: new §8 'The lowering contract (declarative → graph)' — every declarative concept and what it lowers to, with invariants (compound priors add no variables; factor_out not reachable by design; PriorFactor exact-hooks gap noted → #1337/#1338). - AGENTS.md: the EP seam rule — inner-layer capability lands in the same PR as its declarative expression or an explicit 'not exposed' entry, plus a seam test. - test_autofit/graphical/test_declarative_deterministic.py: replaced the commented-out PR #1153 illustration with 4 working seam tests pinning the explicit compound-prior lowering (no new variables, graph shape, in-factor realisation, shared-variable connectivity). The model. sugar remains deliberately reverted (be6411755). EP review seam recommendations (#1330 wrap-up discussion); Phase 5 recommendation A item 2 delivered independently of the #1336 decision. Co-Authored-By: Claude Fable 5 --- AGENTS.md | 9 + autofit/graphical/README.md | 30 ++- .../test_declarative_deterministic.py | 175 +++++++++++------- 3 files changed, 150 insertions(+), 64 deletions(-) diff --git a/AGENTS.md b/AGENTS.md index 692f75fdf..e44f5d7a8 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -72,6 +72,15 @@ nautilus, NSS; MLE: LBFGS/BFGS/drawer), `mapper/` (model + priors), - Import direction: autoconf only — never `autoarray` / `autogalaxy` / `autolens`. +- **The EP seam rule**: `autofit/graphical` is two layers — the inner + factor-graph/message engine and the `declarative/` user layer. A new + statistical capability in the inner layer must land **in the same PR** + with its declarative expression *or* an explicit "not exposed" row in + the lowering-contract table (`autofit/graphical/README.md` §8), plus a + seam test where behaviour crosses the boundary + (`test_autofit/graphical/test_declarative_deterministic.py` is the + pattern). Capabilities that exist below but are silently absent above + are the seam's known failure mode (see PyAutoFit#1336/#1337). - The `[nss]` extra (for `af.NSS`) needs a pinned **handley-lab/blackjax** fork installed manually *after* the extras step; that fork conflicts with the mainline `blackjax` pinned in `[optional]`. Do not naively combine or bump diff --git a/autofit/graphical/README.md b/autofit/graphical/README.md index c833b1a8f..57ced1a88 100644 --- a/autofit/graphical/README.md +++ b/autofit/graphical/README.md @@ -217,7 +217,35 @@ semantics: cavities within a sweep are stale relative to serial EP, so serial and parallel runs converge along different trajectories (to the same fixed points when EP converges). -## 8. Reading +## 8. The lowering contract (declarative → graph) + +`autofit.graphical` is two layers: the **inner layer** above (factor +graphs, messages, EP updates — this document's §1–§7) and the +**declarative layer** (`declarative/`: `FactorGraphModel`, +`AnalysisFactor`, `HierarchicalFactor`) that scientists actually use. +This section is the seam contract: what every declarative concept +*lowers to* on the graph, and which inner-layer capabilities survive +the translation. When either layer changes, this table is the thing to +re-verify (see the seam tests in +`test_autofit/graphical/test_declarative_deterministic.py`). + +| Declarative concept | Lowers to | Notes / invariants | +|---|---|---| +| `FactorGraphModel(*factors)` | `FactorGraph` = the product Eq. (1) | one graph; `mean_field_approximation()` builds the `EPMeanField` | +| `AnalysisFactor(prior_model, analysis, optimiser)` | one `Factor` whose value is `analysis.log_likelihood_function` on the instance built from its variables | carries its own tilted-fit optimiser (§3.2) | +| each free `Prior` | one graph `Variable` **and** one `PriorFactor` | priors are ordinary factors (§1); `PriorFactor` currently wraps the message's bound `factor` method, which strips the exact-update hooks — the conjugate update of §3.2 is *not* auto-selected declaratively (tracked: #1337 / plan #1338 WP1) | +| the *same prior object* assigned to several models | one shared `Variable` connecting those factors | this is how information flows between datasets | +| compound prior (`prior_a * x + prior_b`, `mapper/prior/arithmetic/`) | **no graph variable** — the arithmetic is evaluated at instance-creation inside every factor that references it; only its component priors are variables | the relation is enforced *exactly* inside each tilted fit (no extra approximation — cf. §6); consequently the compound quantity has no message, no marginal, no evidence contribution of its own. (A `model.` sugar for building these was deliberately reverted in `be6411755`.) | +| `HierarchicalFactor` | one `Factor` per drawn variable (plus the distribution's parameter variables) | deliberate dimensionality choice | +| — (no declarative expression) | `Factor(..., factor_out=v)` graph-level deterministic variables (§6.1) | **not reachable** from the declarative layer, by design as of the 2026-07 review (Phase 5, #1336): it trades the exact in-factor relation for a factorised q(v) with messages | + +Contract for contributors: a new statistical capability lands in the +inner layer **together with** its row in this table — either a +declarative expression, or an explicit "not exposed" entry with the +reason. Capabilities that exist below but are silently absent above +(the `PriorFactor` exact-hooks row) are the seam's known failure mode. + +## 9. Reading - T. Minka (2001), *Expectation Propagation for Approximate Bayesian Inference* — the algorithm of §3. diff --git a/test_autofit/graphical/test_declarative_deterministic.py b/test_autofit/graphical/test_declarative_deterministic.py index ce7b00e5b..63eb4f397 100644 --- a/test_autofit/graphical/test_declarative_deterministic.py +++ b/test_autofit/graphical/test_declarative_deterministic.py @@ -1,67 +1,116 @@ +""" +Seam tests for the declarative → graph lowering contract +(`autofit/graphical/README.md` §8), focused on deterministic quantities +expressed as compound priors. + +History: PR #1153 introduced this pattern with an illustration test that +was committed commented-out; the accompanying `model.` sugar +was later deliberately reverted (`be6411755`, "undo annoying change"). +The *explicit* compound-prior pattern below is the supported one — these +tests pin its lowering behaviour so the seam cannot silently drift again +(EP review Phase 5, PyAutoFit #1336). +""" +import pytest + import autofit as af from autofit.mock import MockAnalysis +FWHM_FACTOR = 2.354820045 + + +@pytest.fixture(name="factor_graph_pieces") +def make_factor_graph_pieces(): + model_1 = af.Model(af.ex.Gaussian) + model_2 = af.Model(af.ex.Gaussian) + + # Deterministic quantities built explicitly from other models' priors. + deterministic_model = af.Collection( + model_1.sigma * FWHM_FACTOR, + model_2.sigma * FWHM_FACTOR, + ) + + deterministic_factor = af.AnalysisFactor( + prior_model=deterministic_model, analysis=MockAnalysis() + ) + factor_graph = af.FactorGraphModel( + af.AnalysisFactor(prior_model=model_1, analysis=MockAnalysis()), + af.AnalysisFactor(prior_model=model_2, analysis=MockAnalysis()), + deterministic_factor, + ) + return model_1, model_2, deterministic_model, deterministic_factor, factor_graph + + +def test_compound_lowering_adds_no_variables(factor_graph_pieces): + """ + A compound prior lowers to *no* graph variable: only its component + priors are variables, and they are shared with the source models + (README §8 — the relation is enforced exactly inside each factor). + """ + model_1, model_2, deterministic_model, _, factor_graph = factor_graph_pieces + + source_priors = set(model_1.priors) | set(model_2.priors) + compound_priors = set(deterministic_model.priors) + + assert len(compound_priors) == 2 + assert compound_priors <= source_priors + + +def test_graph_shape(factor_graph_pieces): + """ + 3 AnalysisFactors + one PriorFactor per free parameter (2 models x 3 + parameters = 6). The deterministic AnalysisFactor contributes no new + PriorFactors because its priors are the shared components. + """ + *_, factor_graph = factor_graph_pieces + approx = factor_graph.mean_field_approximation() + + names = [factor.name for factor in approx.factor_graph.factors] + n_analysis = sum("AnalysisFactor" in name for name in names) + n_prior = sum("PriorFactor" in name for name in names) + + assert n_analysis == 3 + assert n_prior == 6 + + # Every variable in the mean field carries a message. + mean_field = approx.mean_field + assert len(mean_field) == 6 + + +def test_compound_realises_inside_factor(factor_graph_pieces): + """ + The deterministic relation is evaluated at instance-creation inside + the factor: realising the deterministic model at the prior medians + gives exactly sigma_median * FWHM_FACTOR. + """ + model_1, _, deterministic_model, _, _ = factor_graph_pieces + + instance = deterministic_model.instance_from_prior_medians() + sigma_median = model_1.sigma.value_for(0.5) + + assert instance[0] == pytest.approx(sigma_median * FWHM_FACTOR) + + +def test_deterministic_factor_connects_through_shared_variables( + factor_graph_pieces, +): + """ + The deterministic AnalysisFactor's graph variables are precisely the + shared component priors — this is the channel through which its + likelihood constrains the source models under EP. + """ + ( + model_1, + model_2, + deterministic_model, + deterministic_factor, + factor_graph, + ) = factor_graph_pieces + approx = factor_graph.mean_field_approximation() + + # Identity-based lookup: AnalysisFactor names come from a global + # counter, so name matching is test-order-dependent. + assert deterministic_factor in approx.factor_graph.factors -def test(): - - pass - - # model_1 = af.Model(af.ex.Gaussian) - # analysis_factor_1 = af.AnalysisFactor( - # prior_model=model_1, - # analysis=MockAnalysis(), - # ) - # - # model_2 = af.Model(af.ex.Gaussian) - # analysis_factor_2 = af.AnalysisFactor( - # prior_model=model_2, - # analysis=MockAnalysis(), - # ) - # - # model_3 = af.Collection( - # model_1.fwhm, - # model_2.fwhm, - # ) - # analysis_factor_3 = af.AnalysisFactor( - # prior_model=model_3, - # analysis=MockAnalysis(), - # ) - # - # factor_graph = af.FactorGraphModel( - # analysis_factor_1, - # analysis_factor_2, - # analysis_factor_3, - # ) - -# assert ( -# factor_graph.info -# == """PriorFactors -# -# PriorFactor0 (AnalysisFactor1.sigma, AnalysisFactor2.1.self) UniformPrior [5], lower_limit = 0.0, upper_limit = 1.0 -# PriorFactor1 (AnalysisFactor1.normalization) UniformPrior [4], lower_limit = 0.0, upper_limit = 1.0 -# PriorFactor2 (AnalysisFactor1.centre) UniformPrior [3], lower_limit = 0.0, upper_limit = 1.0 -# PriorFactor3 (AnalysisFactor0.sigma, AnalysisFactor2.0.self) UniformPrior [2], lower_limit = 0.0, upper_limit = 1.0 -# PriorFactor4 (AnalysisFactor0.normalization) UniformPrior [1], lower_limit = 0.0, upper_limit = 1.0 -# PriorFactor5 (AnalysisFactor0.centre) UniformPrior [0], lower_limit = 0.0, upper_limit = 1.0 -# -# AnalysisFactors -# -# AnalysisFactor0 -# -# centre (PriorFactor5) UniformPrior [0], lower_limit = 0.0, upper_limit = 1.0 -# normalization (PriorFactor4) UniformPrior [1], lower_limit = 0.0, upper_limit = 1.0 -# sigma (AnalysisFactor2.0.self, PriorFactor3) UniformPrior [2], lower_limit = 0.0, upper_limit = 1.0 -# -# AnalysisFactor1 -# -# centre (PriorFactor2) UniformPrior [3], lower_limit = 0.0, upper_limit = 1.0 -# normalization (PriorFactor1) UniformPrior [4], lower_limit = 0.0, upper_limit = 1.0 -# sigma (AnalysisFactor2.1.self, PriorFactor0) UniformPrior [5], lower_limit = 0.0, upper_limit = 1.0 -# -# AnalysisFactor2 -# -# 0 -# self (AnalysisFactor0.sigma, PriorFactor3) UniformPrior [2], lower_limit = 0.0, upper_limit = 1.0 -# 1 -# self (AnalysisFactor1.sigma, PriorFactor0) UniformPrior [5], lower_limit = 0.0, upper_limit = 1.0""" -# ) + factor_variables = set(deterministic_factor.variables) + assert factor_variables == set(deterministic_model.priors) + assert factor_variables <= set(model_1.priors) | set(model_2.priors)