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
7 changes: 6 additions & 1 deletion .github/scripts/smoke_install.sh
Original file line number Diff line number Diff line change
Expand Up @@ -7,4 +7,9 @@
set -e

pip install ./PyAutoNerves ./PyAutoFit
pip install numba
# The gradient tutorials use BlackJAXNUTS. BlackJAX is a PyAutoFit optional
# dependency, so the base install above does not provide it.
pip install numba "blackjax>=1.6.2"
# tutorial_6_gradients.py also runs Nautilus, another PyAutoFit optional
# dependency; pin it as PyAutoFit's `optional` extra does.
pip install "nautilus-sampler==1.0.5"
2 changes: 1 addition & 1 deletion config/build/no_run.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -10,4 +10,4 @@
# every mega-run surfaces them with a loud warning banner. Investigate the
# failure, fix the underlying bug, and remove the NEEDS_FIX marker.

- tutorial_5_expectation_propagation # NEEDS_FIX 2026-08-05 - LinearRegressionAnalysis.log_likelihood_function returns a constant -1, ignoring `instance`, so every initial sample has an identical figure of merit and the per-factor search cannot initialise. The tutorial only "passes" today because CI runs at PYAUTO_TEST_MODE=2, which bypasses the sampler; at real sampling the linear_regression factor never completes a single update. Fix needs a real likelihood for the m/c regression over fwhm_list (currently computed and unused), which is tutorial authoring, not a mechanical repair. See PyAutoFit#1454.
[]
4 changes: 4 additions & 0 deletions config/build/profile_smoke.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -31,6 +31,10 @@ defaults:
# tutorials. Never leave the `overrides:` key with an empty value — a bare
# `overrides:` loads as None; write `overrides: []` if every entry is removed.)
overrides:
# Exercise the EP factor likelihoods with real sampling: mode 2 concealed
# the former constant regression likelihood and unshared centre.
- pattern: "tutorial_5_expectation_propagation"
set: {PYAUTO_TEST_MODE: "0"}
# The gradient tutorials teach JAX autodiff searches (MultiStartAdam, BlackJAXNUTS),
# which refuse to run with use_jax=False, so JAX stays on for them.
- pattern: "tutorial_6_gradients"
Expand Down
2 changes: 1 addition & 1 deletion llms-full.txt
Original file line number Diff line number Diff line change
Expand Up @@ -39,7 +39,7 @@ AUTO-GENERATED by PyAutoHands — do not edit by hand; regenerate with generate.
- [Tutorial 4: Hierarchical](scripts/chapter_advanced/tutorial_4_hierachical_models.py): In the previous tutorial, we fitted a graphical model with the aim of determining an estimate of shared parameters, the `centre`'s of a dataset of 1D Gaussians. We did this by fitting all datasets simultaneously. When there are shared parameters in a model, this is a powerful and effective tool, but things may not always be so simple.
- Contents: Example Source Code (`af.ex`), Dataset, Analysis, Model Individual Factors, Analysis Factors, Model, Factor Graph, Search, Result, Comparison to One-by-One Fits, Benefits of Graphical Model, Wrap Up
- [Tutorial 5: Expectation Propagation](scripts/chapter_advanced/tutorial_5_expectation_propagation.py): In the previous tutorials, we fitted graphical models to dataset comprising many noisy 1D Gaussians. These had a shared and global value of their `centre`, or assumed their centres were hierarchically drawn from a parent Gaussian distribution. This provides the basis of composing and fitting complex graphical models to large datasets.
- Contents: Example Source Code (`af.ex`), Dataset, Analysis, Model, Analysis Factors, Factor Graph, Expectation Propagation, Cyclic Fitting, Result, Output, Results
- Contents: Example Source Code (`af.ex`), Dataset, Analysis, Model, Analysis Factors, A Width Relation, Factor Graph, Expectation Propagation, Cyclic Fitting, Result, Output, Results
- [Tutorial Optional: Hierarchical Expectation Propagation (EP)](scripts/chapter_advanced/tutorial_optional_hierarchical_ep.py): This optional tutorial gives an example of fitting a hierarchical model using EP.
- Contents: Example Source Code (`af.ex`), Dataset, Analysis, Model Individual Factors, Analysis Factors, Model, Factor Graph, Model Fit
- [Tutorial Optional: Hierarchical Individual](scripts/chapter_advanced/tutorial_optional_hierarchical_individual.py): In tutorial 4, we fit a hierarchical model using a graphical model, whereby all datasets are fitted simultaneously and the hierarchical parameters are fitted for simultaneously with the model parameters of each 1D Gaussian in each dataset.
Expand Down
87 changes: 75 additions & 12 deletions notebooks/chapter_advanced/tutorial_5_expectation_propagation.ipynb
Original file line number Diff line number Diff line change
Expand Up @@ -32,6 +32,7 @@
"- **Analysis**: Create Analysis objects for each dataset.\n",
"- **Model**: Set up the model with a shared centre prior across all datasets.\n",
"- **Analysis Factors**: Create Analysis Factors with individual searches for each dataset.\n",
"- **A Width Relation**: Connect the inferred widths with an explicitly assumed linear relation.\n",
"- **Factor Graph**: Compose the factor graph for the EP framework.\n",
"- **Expectation Propagation**: Explain the EP message passing algorithm.\n",
"- **Cyclic Fitting**: Describe the iterative EP convergence process.\n",
Expand Down Expand Up @@ -245,7 +246,7 @@
"for model_index in range(len(data_list)):\n",
" gaussian = af.Model(af.ex.Gaussian)\n",
"\n",
" gaussian.centre = af.GaussianPrior(mean=50.0, sigma=30.0)\n",
" gaussian.centre = centre_shared_prior\n",
" gaussian.normalization = af.TruncatedGaussianPrior(\n",
" mean=3.0, sigma=5.0, lower_limit=0.0\n",
" )\n",
Expand Down Expand Up @@ -297,23 +298,56 @@
" prior_model=model, analysis=analysis, optimiser=search, name=dataset_name\n",
" )\n",
"\n",
" analysis_factor_list.append(analysis_factor)\n",
" analysis_factor_list.append(analysis_factor)\n"
],
"outputs": [],
"execution_count": null
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"__A Width Relation__\n",
"\n",
"To show how another factor can connect to the same datasets, we also model their full widths at half maximum\n",
"(FWHMs). For this teaching example, assume the widths scatter around a line `m * x + c`, where `x` is the\n",
"dataset index (0, 1 or 2). The scatter is fixed at 5 pixels. This is an additional modelling assumption, not\n",
"a second set of measured widths: the widths remain parameters inferred from the noisy profiles above.\n",
"\n",
"The regression factor evaluates a normalized Gaussian density for each width about that line. It therefore\n",
"passes information about the widths between datasets, while the shared centre connects the profile factors\n",
"directly. The simulated widths are similar, so a slope near zero is plausible; we do not fix it to zero."
]
},
{
"cell_type": "code",
"metadata": {},
"source": [
"\n",
"\n",
"class LinearRegressionAnalysis(af.Analysis):\n",
" def __init__(self, xvalues, scatter):\n",
" self.xvalues = xvalues\n",
" self.scatter = scatter\n",
"\n",
" def log_likelihood_function(self, instance, xp=np):\n",
" return -1\n",
" widths = xp.asarray(list(instance.fwhm))\n",
" prediction = instance.m * xp.asarray(self.xvalues) + instance.c\n",
" residuals = (widths - prediction) / self.scatter\n",
" return -0.5 * xp.sum(\n",
" residuals**2 + xp.log(2.0 * xp.pi * self.scatter**2)\n",
" )\n",
"\n",
"\n",
"fwhm_list = [2 * np.sqrt(2 * np.log(2)) * model.gaussian.sigma for model in model_list]\n",
"\n",
"linear_regression_factor = af.AnalysisFactor(\n",
" prior_model=af.Collection(\n",
" m=af.GaussianPrior(mean=0.0, sigma=1.0),\n",
" c=af.GaussianPrior(mean=0.0, sigma=1.0),\n",
" *fwhm_list,\n",
" m=af.GaussianPrior(mean=0.0, sigma=5.0),\n",
" c=af.GaussianPrior(mean=2 * np.sqrt(2 * np.log(2)) * 10.0, sigma=10.0),\n",
" fwhm=af.Collection(*fwhm_list),\n",
" ),\n",
" analysis=LinearRegressionAnalysis(),\n",
" analysis=LinearRegressionAnalysis(xvalues=np.arange(total_datasets), scatter=5.0),\n",
" optimiser=search,\n",
" name=\"linear_regression\",\n",
")"
Expand Down Expand Up @@ -398,7 +432,7 @@
"\n",
"2) Once the model-fit is complete, inspect the model for parameters that are shared with other nodes on the factor\n",
"graph. In this example, the `centre` of the `Gaussian` fitted to the first dataset is global, and therefore connects\n",
"to the other nodes on the factor graph (the `AnalysisFactor`'s) of the second and first `Gaussian` datasets.\n",
"to the other nodes on the factor graph (the `AnalysisFactor`'s) of the second and third `Gaussian` datasets.\n",
"\n",
"3) The EP framework now creates a 'message' that is to be passed to the connecting nodes on the factor graph. This\n",
"message informs them of the results of the model-fit, so they can update their priors on the `Gaussian`'s centre \n",
Expand All @@ -410,7 +444,7 @@
"\n",
"__Cyclic Fitting__\n",
"\n",
"After every `AnalysisFactor` has been fitted (e.g. after each fit to each of the 5 datasets in this example), we have a \n",
"After every `AnalysisFactor` has been fitted (the three datasets and the width relation in this example), we have a\n",
"new estimate of the shared parameter `centre`. This updates our priors on the shared parameter `centre`, which needs \n",
"to be reflected in each model-fit we perform on each `AnalysisFactor`. \n",
"\n",
Expand Down Expand Up @@ -447,6 +481,35 @@
"outputs": [],
"execution_count": null
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Five sweeps demonstrate message passing, but do not guarantee convergence. Before interpreting the result,\n",
"check that every data or regression factor contributed at least one update. A run that exits normally can\n",
"still contain a stale factor; the guard below makes that failure explicit, including in automated smoke tests.\n",
"An update can also reject an individual marginal projection, which the state figure below marks separately."
]
},
{
"cell_type": "code",
"metadata": {},
"source": [
"factor_updates = {\n",
" factor.name: len(history.updates)\n",
" for factor, history in factor_graph_result.ep_history.items()\n",
"}\n",
"stale_factors = [\n",
" name\n",
" for name in [*dataset_name_list, \"linear_regression\"]\n",
" if factor_updates.get(name, 0) == 0\n",
"]\n",
"if stale_factors:\n",
" raise RuntimeError(f\"EP factors never updated: {', '.join(stale_factors)}\")"
],
"outputs": [],
"execution_count": null
},
{
"cell_type": "markdown",
"metadata": {},
Expand Down Expand Up @@ -480,9 +543,9 @@
"to read it: `box = factor`, `pill = variable`, `dashed frame = repeated structure`. Every `AnalysisFactor` is a box, every\n",
"parameter it touches is a pill, and an edge between them means \"this factor has something to say about this\n",
"variable\". That is the picture the message passing above lives on: a message leaves a box, arrives at a pill, and is\n",
"read by every other box attached to it. The three datasets sit inside a dashed plate badged `3 datasets`, drawn once\n",
"as an aggregate `AnalysisFactor` box, with the individual `dataset_0`, `dataset_1` and `dataset_2` boxes expanded\n",
"alongside it because their states differ from one another.\n",
"read by every other box attached to it. Follow the shared `centre` pill to see the connection between the three\n",
"dataset factors. The width relation connects their `sigma` parameters through the FWHM expressions. Repeated\n",
"structure may be grouped into a plate; individual factors are expanded when their run states differ.\n",
"\n",
"Onto that structure the `ep_history` paints what the sweeps actually did. Every box carries a status:\n",
"\n",
Expand Down
64 changes: 53 additions & 11 deletions scripts/chapter_advanced/tutorial_5_expectation_propagation.py
Original file line number Diff line number Diff line change
Expand Up @@ -27,6 +27,7 @@
- **Analysis**: Create Analysis objects for each dataset.
- **Model**: Set up the model with a shared centre prior across all datasets.
- **Analysis Factors**: Create Analysis Factors with individual searches for each dataset.
- **A Width Relation**: Connect the inferred widths with an explicitly assumed linear relation.
- **Factor Graph**: Compose the factor graph for the EP framework.
- **Expectation Propagation**: Explain the EP message passing algorithm.
- **Cyclic Fitting**: Describe the iterative EP convergence process.
Expand Down Expand Up @@ -142,7 +143,7 @@
for model_index in range(len(data_list)):
gaussian = af.Model(af.ex.Gaussian)

gaussian.centre = af.GaussianPrior(mean=50.0, sigma=30.0)
gaussian.centre = centre_shared_prior
gaussian.normalization = af.TruncatedGaussianPrior(
mean=3.0, sigma=5.0, lower_limit=0.0
)
Expand Down Expand Up @@ -186,20 +187,43 @@
analysis_factor_list.append(analysis_factor)


"""
__A Width Relation__

To show how another factor can connect to the same datasets, we also model their full widths at half maximum
(FWHMs). For this teaching example, assume the widths scatter around a line `m * x + c`, where `x` is the
dataset index (0, 1 or 2). The scatter is fixed at 5 pixels. This is an additional modelling assumption, not
a second set of measured widths: the widths remain parameters inferred from the noisy profiles above.

The regression factor evaluates a normalized Gaussian density for each width about that line. It therefore
passes information about the widths between datasets, while the shared centre connects the profile factors
directly. The simulated widths are similar, so a slope near zero is plausible; we do not fix it to zero.
"""


class LinearRegressionAnalysis(af.Analysis):
def __init__(self, xvalues, scatter):
self.xvalues = xvalues
self.scatter = scatter

def log_likelihood_function(self, instance, xp=np):
return -1
widths = xp.asarray(list(instance.fwhm))
prediction = instance.m * xp.asarray(self.xvalues) + instance.c
residuals = (widths - prediction) / self.scatter
return -0.5 * xp.sum(
residuals**2 + xp.log(2.0 * xp.pi * self.scatter**2)
)


fwhm_list = [2 * np.sqrt(2 * np.log(2)) * model.gaussian.sigma for model in model_list]

linear_regression_factor = af.AnalysisFactor(
prior_model=af.Collection(
m=af.GaussianPrior(mean=0.0, sigma=1.0),
c=af.GaussianPrior(mean=0.0, sigma=1.0),
*fwhm_list,
m=af.GaussianPrior(mean=0.0, sigma=5.0),
c=af.GaussianPrior(mean=2 * np.sqrt(2 * np.log(2)) * 10.0, sigma=10.0),
fwhm=af.Collection(*fwhm_list),
),
analysis=LinearRegressionAnalysis(),
analysis=LinearRegressionAnalysis(xvalues=np.arange(total_datasets), scatter=5.0),
optimiser=search,
name="linear_regression",
)
Expand Down Expand Up @@ -245,7 +269,7 @@ def log_likelihood_function(self, instance, xp=np):

2) Once the model-fit is complete, inspect the model for parameters that are shared with other nodes on the factor
graph. In this example, the `centre` of the `Gaussian` fitted to the first dataset is global, and therefore connects
to the other nodes on the factor graph (the `AnalysisFactor`'s) of the second and first `Gaussian` datasets.
to the other nodes on the factor graph (the `AnalysisFactor`'s) of the second and third `Gaussian` datasets.

3) The EP framework now creates a 'message' that is to be passed to the connecting nodes on the factor graph. This
message informs them of the results of the model-fit, so they can update their priors on the `Gaussian`'s centre
Expand All @@ -257,7 +281,7 @@ def log_likelihood_function(self, instance, xp=np):

__Cyclic Fitting__

After every `AnalysisFactor` has been fitted (e.g. after each fit to each of the 5 datasets in this example), we have a
After every `AnalysisFactor` has been fitted (the three datasets and the width relation in this example), we have a
new estimate of the shared parameter `centre`. This updates our priors on the shared parameter `centre`, which needs
to be reflected in each model-fit we perform on each `AnalysisFactor`.

Expand Down Expand Up @@ -286,6 +310,24 @@ def log_likelihood_function(self, instance, xp=np):
visualise_interval=1,
)

"""
Five sweeps demonstrate message passing, but do not guarantee convergence. Before interpreting the result,
check that every data or regression factor contributed at least one update. A run that exits normally can
still contain a stale factor; the guard below makes that failure explicit, including in automated smoke tests.
An update can also reject an individual marginal projection, which the state figure below marks separately.
"""
factor_updates = {
factor.name: len(history.updates)
for factor, history in factor_graph_result.ep_history.items()
}
stale_factors = [
name
for name in [*dataset_name_list, "linear_regression"]
if factor_updates.get(name, 0) == 0
]
if stale_factors:
raise RuntimeError(f"EP factors never updated: {', '.join(stale_factors)}")

"""
__Result__

Expand All @@ -305,9 +347,9 @@ def log_likelihood_function(self, instance, xp=np):
to read it: `box = factor`, `pill = variable`, `dashed frame = repeated structure`. Every `AnalysisFactor` is a box, every
parameter it touches is a pill, and an edge between them means "this factor has something to say about this
variable". That is the picture the message passing above lives on: a message leaves a box, arrives at a pill, and is
read by every other box attached to it. The three datasets sit inside a dashed plate badged `3 datasets`, drawn once
as an aggregate `AnalysisFactor` box, with the individual `dataset_0`, `dataset_1` and `dataset_2` boxes expanded
alongside it because their states differ from one another.
read by every other box attached to it. Follow the shared `centre` pill to see the connection between the three
dataset factors. The width relation connects their `sigma` parameters through the FWHM expressions. Repeated
structure may be grouped into a plate; individual factors are expanded when their run states differ.

Onto that structure the `ep_history` paints what the sweeps actually did. Every box carries a status:

Expand Down
1 change: 1 addition & 0 deletions workspace_index.json
Original file line number Diff line number Diff line change
Expand Up @@ -308,6 +308,7 @@
"Analysis",
"Model",
"Analysis Factors",
"A Width Relation",
"Factor Graph",
"Expectation Propagation",
"Cyclic Fitting",
Expand Down
Loading