diff --git a/.github/scripts/smoke_install.sh b/.github/scripts/smoke_install.sh index 47ed9a6..e236953 100755 --- a/.github/scripts/smoke_install.sh +++ b/.github/scripts/smoke_install.sh @@ -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" diff --git a/config/build/no_run.yaml b/config/build/no_run.yaml index ab92eef..e273e2a 100644 --- a/config/build/no_run.yaml +++ b/config/build/no_run.yaml @@ -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. +[] diff --git a/config/build/profile_smoke.yaml b/config/build/profile_smoke.yaml index 52a10d7..ce0f149 100644 --- a/config/build/profile_smoke.yaml +++ b/config/build/profile_smoke.yaml @@ -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" diff --git a/llms-full.txt b/llms-full.txt index 9218edd..553a9fe 100644 --- a/llms-full.txt +++ b/llms-full.txt @@ -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. diff --git a/notebooks/chapter_advanced/tutorial_5_expectation_propagation.ipynb b/notebooks/chapter_advanced/tutorial_5_expectation_propagation.ipynb index 3081c09..0b3f2e6 100644 --- a/notebooks/chapter_advanced/tutorial_5_expectation_propagation.ipynb +++ b/notebooks/chapter_advanced/tutorial_5_expectation_propagation.ipynb @@ -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", @@ -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", @@ -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", ")" @@ -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", @@ -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", @@ -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": {}, @@ -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", diff --git a/scripts/chapter_advanced/tutorial_5_expectation_propagation.py b/scripts/chapter_advanced/tutorial_5_expectation_propagation.py index 4d263ed..3d55c55 100644 --- a/scripts/chapter_advanced/tutorial_5_expectation_propagation.py +++ b/scripts/chapter_advanced/tutorial_5_expectation_propagation.py @@ -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. @@ -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 ) @@ -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", ) @@ -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 @@ -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`. @@ -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__ @@ -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: diff --git a/workspace_index.json b/workspace_index.json index 3bcadab..af7ab5e 100644 --- a/workspace_index.json +++ b/workspace_index.json @@ -308,6 +308,7 @@ "Analysis", "Model", "Analysis Factors", + "A Width Relation", "Factor Graph", "Expectation Propagation", "Cyclic Fitting",