diff --git a/docs/general/images/model_cookbook/acceptance_group_scale.png b/docs/general/images/model_cookbook/acceptance_group_scale.png new file mode 100644 index 0000000000..58ec36a974 Binary files /dev/null and b/docs/general/images/model_cookbook/acceptance_group_scale.png differ diff --git a/docs/general/images/model_cookbook/acceptance_mge_pixelized.png b/docs/general/images/model_cookbook/acceptance_mge_pixelized.png new file mode 100644 index 0000000000..15fe26a6d9 Binary files /dev/null and b/docs/general/images/model_cookbook/acceptance_mge_pixelized.png differ diff --git a/docs/general/images/model_cookbook/acceptance_simple_lens.png b/docs/general/images/model_cookbook/acceptance_simple_lens.png new file mode 100644 index 0000000000..2019ca24f4 Binary files /dev/null and b/docs/general/images/model_cookbook/acceptance_simple_lens.png differ diff --git a/docs/general/images/model_cookbook/complex_lens.png b/docs/general/images/model_cookbook/complex_lens.png new file mode 100644 index 0000000000..fe45b699cf Binary files /dev/null and b/docs/general/images/model_cookbook/complex_lens.png differ diff --git a/docs/general/images/model_cookbook/concise.png b/docs/general/images/model_cookbook/concise.png new file mode 100644 index 0000000000..2f91e267fe Binary files /dev/null and b/docs/general/images/model_cookbook/concise.png differ diff --git a/docs/general/images/model_cookbook/four_galaxies.png b/docs/general/images/model_cookbook/four_galaxies.png new file mode 100644 index 0000000000..cd11b58152 Binary files /dev/null and b/docs/general/images/model_cookbook/four_galaxies.png differ diff --git a/docs/general/images/model_cookbook/make_figures.py b/docs/general/images/model_cookbook/make_figures.py new file mode 100644 index 0000000000..e07825385c --- /dev/null +++ b/docs/general/images/model_cookbook/make_figures.py @@ -0,0 +1,317 @@ +""" +Regenerate the committed model-figure renders for the PyAutoLens docs. + +Run from the PyAutoLens repo root, with the worktree environment active:: + + python docs/general/images/model_cookbook/make_figures.py + +Every PNG is written at ``width=14.0`` inches and DPI 100, so none exceeds the +1400 px width budget of the ``model-figures`` epic. + +Two families of figure are written: + +* the **cookbook stages** -- one per ``print(model.info)`` stage of + ``autolens_workspace/scripts/guides/modeling/cookbook.py``, plus the + ``Redshift Free`` and ``Solved Parameters`` stages, into this directory and + embedded by ``docs/general/model_cookbook.md``; +* the **acceptance trio** of the epic (the simple lens, the MGE + pixelized + source and the group-scale model), imported from + ``test_autolens/model_figure/lens_models.py`` so the docs and the acceptance + tests render the identical models. + +``docs/overview/images/overview_3/`` additionally receives the three feature +figures embedded by ``docs/overview/overview_3_features.md``. +""" + +import sys +from pathlib import Path + +import matplotlib + +matplotlib.use("Agg") + +REPO_ROOT = Path(__file__).resolve().parents[4] +sys.path.insert(0, str(REPO_ROOT)) +sys.path.insert(0, str(REPO_ROOT / "test_autolens" / "model_figure")) + +import autofit as af +import autolens as al + +import lens_models + +OUTPUT_PATH = Path(__file__).resolve().parent +OVERVIEW_PATH = REPO_ROOT / "docs" / "overview" / "images" / "overview_3" +WIDTH = 14.0 + + +def write(model, name: str, path: Path = OUTPUT_PATH, **kwargs): + """Render one model to ``.png`` in ``path``.""" + path.mkdir(parents=True, exist_ok=True) + + af.ModelPlotter(model).figure( + path=str(path), + filename=name, + format="png", + width=WIDTH, + **kwargs, + ) + print(f"wrote {path.name}/{name}.png") + + +# --------------------------------------------------------------------------- +# the cookbook stages -- `scripts/guides/modeling/cookbook.py` +# --------------------------------------------------------------------------- + + +def simple_lens(): + """``__Simple Lens Model__`` -- cookbook.py:53-82.""" + bulge = af.Model(al.lp_linear.Sersic) + mass = af.Model(al.mp.Isothermal) + + lens = af.Model(al.Galaxy, redshift=0.5, bulge=bulge, mass=mass) + + source = af.Model(al.Galaxy, redshift=1.0, bulge=al.lp_linear.SersicCore) + + return af.Collection(galaxies=af.Collection(lens=lens, source=source)) + + +def complex_lens(): + """``__More Complex Lens Models__`` -- cookbook.py:89-117.""" + bulge = af.Model(al.lp_linear.Sersic) + disk = af.Model(al.lp_linear.Exponential) + + mass = af.Model(al.mp.Isothermal) + shear = af.Model(al.mp.ExternalShear) + + lens = af.Model( + al.Galaxy, redshift=0.5, bulge=bulge, disk=disk, mass=mass, shear=shear + ) + + bulge = af.Model(al.lp_linear.SersicCore) + disk = af.Model(al.lp_linear.ExponentialCore) + + source = af.Model(al.Galaxy, redshift=1.0, bulge=bulge, disk=disk) + + return af.Collection(galaxies=af.Collection(lens=lens, source=source)) + + +def four_galaxies(): + """``__More Complex Lens Models__``, four galaxies -- cookbook.py:123-166.""" + lens_0 = af.Model( + al.Galaxy, + redshift=0.5, + bulge=af.Model(al.lp_linear.Sersic), + mass=af.Model(al.mp.Isothermal), + ) + lens_1 = af.Model( + al.Galaxy, + redshift=0.5, + bulge=af.Model(al.lp_linear.Sersic), + mass=af.Model(al.mp.Isothermal), + ) + source_0 = af.Model( + al.Galaxy, redshift=1.0, bulge=af.Model(al.lp_linear.SersicCore) + ) + source_1 = af.Model( + al.Galaxy, redshift=1.0, bulge=af.Model(al.lp_linear.SersicCore) + ) + + return af.Collection( + galaxies=af.Collection( + lens_0=lens_0, lens_1=lens_1, source_0=source_0, source_1=source_1 + ), + ) + + +def concise(): + """``__Concise API__`` -- cookbook.py:183-200.""" + lens = af.Model( + al.Galaxy, + redshift=0.5, + bulge=al.lp_linear.Sersic, + disk=al.lp_linear.Sersic, + mass=al.mp.Isothermal, + shear=al.mp.ExternalShear, + ) + + source = af.Model( + al.Galaxy, + redshift=1.0, + bulge=al.lp_linear.SersicCore, + disk=al.lp_linear.ExponentialCore, + ) + + return af.Collection(galaxies=af.Collection(lens=lens, source=source)) + + +def prior_custom(): + """``__Prior Customization__`` -- cookbook.py:221-256.""" + bulge = af.Model(al.lp_linear.Sersic) + bulge.sersic_index = af.TruncatedGaussianPrior( + mean=4.0, sigma=1.0, lower_limit=1.0, upper_limit=8.0 + ) + + mass = af.Model(al.mp.Isothermal) + mass.centre.centre_0 = af.TruncatedGaussianPrior( + mean=0.0, sigma=0.1, lower_limit=-0.5, upper_limit=0.5 + ) + mass.centre.centre_1 = af.TruncatedGaussianPrior( + mean=0.0, sigma=0.1, lower_limit=-0.5, upper_limit=0.5 + ) + mass.einstein_radius = af.UniformPrior(lower_limit=0.0, upper_limit=8.0) + + lens = af.Model(al.Galaxy, redshift=0.5, bulge=bulge, mass=mass) + + bulge = af.Model(al.lp_linear.SersicCore) + + source = af.Model(al.Galaxy, redshift=1.0, bulge=bulge) + source.effective_radius = af.TruncatedGaussianPrior( + mean=0.1, sigma=0.05, lower_limit=0.0, upper_limit=1.0 + ) + + return af.Collection(galaxies=af.Collection(lens=lens, source=source)) + + +def model_custom(): + """``__Model Customization__`` -- cookbook.py:263-319.""" + bulge = af.Model(al.lp_linear.Sersic) + disk = af.Model(al.lp_linear.Exponential) + + bulge.centre = disk.centre + bulge.sersic_index = 4.0 + + mass = af.Model(al.mp.Isothermal) + mass.centre.centre_0 = bulge.centre.centre_0 + 0.1 + mass.centre.centre_1 = bulge.centre.centre_1 + 0.1 + + shear = af.Model(al.mp.ExternalShear) + + lens = af.Model( + al.Galaxy, redshift=0.5, bulge=bulge, disk=disk, mass=mass, shear=shear + ) + + bulge = af.Model(al.lp_linear.SersicCore) + disk = af.Model(al.lp_linear.ExponentialCore) + + source = af.Model(al.Galaxy, redshift=1.0, bulge=bulge, disk=disk) + + model = af.Collection(galaxies=af.Collection(lens=lens, source=source)) + + model.add_assertion( + model.galaxies.lens.bulge.effective_radius + > model.galaxies.lens.disk.effective_radius + ) + model.add_assertion(model.galaxies.lens.mass.einstein_radius < 3.0) + + return model + + +def redshift_free(): + """``__Redshift Free__`` -- cookbook.py:327-329, brought into a full model.""" + redshift = af.Model(al.Redshift) + redshift.redshift = af.UniformPrior(lower_limit=0.0, upper_limit=2.0) + + lens = af.Model(al.Galaxy, redshift=redshift, mass=al.mp.Isothermal) + + source = af.Model(al.Galaxy, redshift=1.0, bulge=al.lp_linear.SersicCore) + + return af.Collection(galaxies=af.Collection(lens=lens, source=source)) + + +def solved_parameters(): + """ + ``__Solved Parameters__`` -- the three solved encodings in one model. + + A linear ``Sersic`` bulge (solved ``intensity``), a pixelized source (solved + ``reconstruction``, and a ``Delaunay`` mesh whose ``areas_factor`` is + ``missing``) and a ``PointSolved`` point source (solved ``centre``). + """ + lens = af.Model( + al.Galaxy, + redshift=0.5, + bulge=af.Model(al.lp_linear.Sersic), + mass=af.Model(al.mp.Isothermal), + ) + + pixelization = af.Model( + al.Pixelization, + mesh=af.Model(al.mesh.Delaunay, pixels=500, zeroed_pixels=0), + regularization=af.Model(al.reg.ConstantSplit), + ) + + source = af.Model(al.Galaxy, redshift=1.0, pixelization=pixelization) + + point_galaxy = af.Model(al.Galaxy, redshift=1.0, point=af.Model(al.ps.PointSolved)) + + return af.Collection( + galaxies=af.Collection(lens=lens, source=source, point_source=point_galaxy) + ) + + +# --------------------------------------------------------------------------- +# the overview_3 feature figures +# --------------------------------------------------------------------------- + + +def mge_overview(): + """The MGE lens light model of ``features/multi_gaussian_expansion/modeling.py``.""" + lens = af.Model( + al.Galaxy, + redshift=0.5, + bulge=lens_models.mge_basis(), + mass=af.Model(al.mp.Isothermal), + shear=af.Model(al.mp.ExternalShear), + ) + + source = af.Model(al.Galaxy, redshift=1.0, bulge=al.lp_linear.SersicCore) + + return af.Collection(galaxies=af.Collection(lens=lens, source=source)) + + +def pixelization_overview(): + """The pixelized source model of ``features/pixelization/modeling.py``:274-289.""" + mass = af.Model(al.mp.PowerLaw) + shear = af.Model(al.mp.ExternalShear) + + lens = af.Model(al.Galaxy, redshift=0.5, mass=mass, shear=shear) + + mesh = af.Model(al.mesh.RectangularBilinearAdaptDensity, shape=(30, 30)) + regularization = af.Model(al.reg.Constant) + + pixelization = af.Model(al.Pixelization, mesh=mesh, regularization=regularization) + + source = af.Model(al.Galaxy, redshift=1.0, pixelization=pixelization) + + return af.Collection(galaxies=af.Collection(lens=lens, source=source)) + + +def point_overview(): + """A point-source lens model -- an ``Isothermal`` lens and an ``al.ps.Point``.""" + lens = af.Model(al.Galaxy, redshift=0.5, mass=af.Model(al.mp.Isothermal)) + + source = af.Model(al.Galaxy, redshift=1.0, point_0=af.Model(al.ps.Point)) + + return af.Collection(galaxies=af.Collection(lens=lens, source=source)) + + +def main(): + write(simple_lens(), "simple_lens") + write(complex_lens(), "complex_lens") + write(four_galaxies(), "four_galaxies") + write(concise(), "concise") + write(prior_custom(), "prior_custom") + write(model_custom(), "model_custom") + write(redshift_free(), "redshift_free") + write(solved_parameters(), "solved_parameters") + + write(lens_models.simple_lens(), "acceptance_simple_lens") + write(lens_models.mge_pixelized(), "acceptance_mge_pixelized") + write(lens_models.group_scale(), "acceptance_group_scale") + + write(mge_overview(), "mge_model", path=OVERVIEW_PATH) + write(pixelization_overview(), "pixelization_model", path=OVERVIEW_PATH) + write(point_overview(), "point_model", path=OVERVIEW_PATH) + + +if __name__ == "__main__": + main() diff --git a/docs/general/images/model_cookbook/model_custom.png b/docs/general/images/model_cookbook/model_custom.png new file mode 100644 index 0000000000..567df980be Binary files /dev/null and b/docs/general/images/model_cookbook/model_custom.png differ diff --git a/docs/general/images/model_cookbook/prior_custom.png b/docs/general/images/model_cookbook/prior_custom.png new file mode 100644 index 0000000000..51e1c7e140 Binary files /dev/null and b/docs/general/images/model_cookbook/prior_custom.png differ diff --git a/docs/general/images/model_cookbook/redshift_free.png b/docs/general/images/model_cookbook/redshift_free.png new file mode 100644 index 0000000000..c39ecc0ee0 Binary files /dev/null and b/docs/general/images/model_cookbook/redshift_free.png differ diff --git a/docs/general/images/model_cookbook/simple_lens.png b/docs/general/images/model_cookbook/simple_lens.png new file mode 100644 index 0000000000..69aaac4b60 Binary files /dev/null and b/docs/general/images/model_cookbook/simple_lens.png differ diff --git a/docs/general/images/model_cookbook/solved_parameters.png b/docs/general/images/model_cookbook/solved_parameters.png new file mode 100644 index 0000000000..3bdc63b233 Binary files /dev/null and b/docs/general/images/model_cookbook/solved_parameters.png differ diff --git a/docs/general/model_cookbook.md b/docs/general/model_cookbook.md index 232bf0ee53..ed04a951d6 100644 --- a/docs/general/model_cookbook.md +++ b/docs/general/model_cookbook.md @@ -10,13 +10,13 @@ readable code for different use-cases. ## Simple Lens Model -A simple lens model has a lens galaxy with a Sersic light profile, Isothermal mass profile and -source galaxy with a Sersic light profile: +A simple lens model has a lens galaxy with a linear Sersic light profile, Isothermal mass profile and +source galaxy with a linear Sersic light profile: ```python # Lens: -bulge = af.Model(al.lp.Sersic) +bulge = af.Model(al.lp_linear.Sersic) mass = af.Model(al.mp.Isothermal) lens = af.Model( @@ -28,9 +28,7 @@ lens = af.Model( # Source: -bulge = af.Model(al.lp.SersicCore) - -source = af.Model(al.Galaxy, redshift=1.0, bulge=bulge) +source = af.Model(al.Galaxy, redshift=1.0, bulge=al.lp_linear.SersicCore) # Overall Lens Model: @@ -39,9 +37,10 @@ model = af.Collection(galaxies=af.Collection(lens=lens, source=source)) The redshifts in the above model are used to determine which galaxy is the lens and which is the source. -The model `total_free_parameters` attribute which tells us the total number of free parameters (which are fitted for -via a non-linear search), which in this case is 19 (7 from the lens Sersic, 5 from the lens Isothermal and 7 from the -source Sersic). +The model `total_free_parameters` attribute tells us the total number of free parameters (which are fitted for +via a non-linear search), which in this case is 17 (6 from the lens `Sersic`, 5 from the lens `Isothermal` and 6 from +the source `SersicCore`). The `intensity` of each linear light profile is *not* one of them: it is solved for by the +inversion during every likelihood evaluation. ```python print(f"Model Total Free Parameters = {model.total_free_parameters}") @@ -56,39 +55,67 @@ print(model.info) This gives the following output: ```bash +Total Free Parameters = 17 + +model Collection (N=17) + galaxies Collection (N=17) + lens Galaxy (N=11) + bulge Sersic (N=6) + mass Isothermal (N=5) + source Galaxy (N=6) + bulge SersicCore (N=6) + galaxies lens - redshift 0.5 + redshift 0.5 bulge centre - centre_0 GaussianPrior, mean = 0.0, sigma = 0.3 - centre_1 GaussianPrior, mean = 0.0, sigma = 0.3 + centre_0 GaussianPrior [0], mean = 0.0, sigma = 0.3 + centre_1 GaussianPrior [1], mean = 0.0, sigma = 0.3 ell_comps - ell_comps_0 GaussianPrior, mean = 0.0, sigma = 0.5 - ell_comps_1 GaussianPrior, mean = 0.0, sigma = 0.5 - intensity LogUniformPrior, lower_limit = 1e-06, upper_limit = 1000000.0 - effective_radius UniformPrior, lower_limit = 0.0, upper_limit = 30.0 + ell_comps_0 TruncatedGaussianPrior [2], mean = 0.0, sigma = 0.3, lower_limit = -1.0, upper_limit = 1.0 + ell_comps_1 TruncatedGaussianPrior [3], mean = 0.0, sigma = 0.3, lower_limit = -1.0, upper_limit = 1.0 + effective_radius UniformPrior [4], lower_limit = 0.0, upper_limit = 30.0 + sersic_index UniformPrior [5], lower_limit = 0.8, upper_limit = 5.0 mass centre - centre_0 GaussianPrior, mean = 0.0, sigma = 0.1 - centre_1 GaussianPrior, mean = 0.0, sigma = 0.1 + centre_0 GaussianPrior [6], mean = 0.0, sigma = 0.1 + centre_1 GaussianPrior [7], mean = 0.0, sigma = 0.1 ell_comps - ell_comps_0 GaussianPrior, mean = 0.0, sigma = 0.3 - ell_comps_1 GaussianPrior, mean = 0.0, sigma = 0.3 - einstein_radius UniformPrior, lower_limit = 0.0, upper_limit = 8.0 + ell_comps_0 TruncatedGaussianPrior [8], mean = 0.0, sigma = 0.3, lower_limit = -1.0, upper_limit = 1.0 + ell_comps_1 TruncatedGaussianPrior [9], mean = 0.0, sigma = 0.3, lower_limit = -1.0, upper_limit = 1.0 + einstein_radius UniformPrior [10], lower_limit = 0.0, upper_limit = 8.0 source - redshift 1.0 - disk + redshift 1.0 + bulge centre - centre_0 GaussianPrior, mean = 0.0, sigma = 0.3 - centre_1 GaussianPrior, mean = 0.0, sigma = 0.3 + centre_0 GaussianPrior [11], mean = 0.0, sigma = 0.3 + centre_1 GaussianPrior [12], mean = 0.0, sigma = 0.3 ell_comps - ell_comps_0 GaussianPrior, mean = 0.0, sigma = 0.5 - ell_comps_1 GaussianPrior, mean = 0.0, sigma = 0.5 - intensity LogUniformPrior, lower_limit = 1e-06, upper_limit = 1000000.0 - effective_radius UniformPrior, lower_limit = 0.0, upper_limit = 30.0 + ell_comps_0 TruncatedGaussianPrior [13], mean = 0.0, sigma = 0.3, lower_limit = -1.0, upper_limit = 1.0 + ell_comps_1 TruncatedGaussianPrior [14], mean = 0.0, sigma = 0.3, lower_limit = -1.0, upper_limit = 1.0 + effective_radius UniformPrior [15], lower_limit = 0.0, upper_limit = 30.0 + sersic_index UniformPrior [16], lower_limit = 0.8, upper_limit = 5.0 + radius_break 0.025 + gamma 0.25 + alpha 3.0 +``` + +The same model can be drawn as a figure, which shows its structure at a glance: + +```python +af.ModelPlotter(model).figure() ``` +```{image} https://raw.githubusercontent.com/PyAutoLabs/PyAutoLens/main/docs/general/images/model_cookbook/simple_lens.png +:alt: The simple lens model, drawn as nested component cards with one pill per parameter. +:width: 600 +``` + +Every stage below prints `model.info` in the same way; from here on only the figure is shown, because the figure is +the quicker read of the two. The `Reading the Figure` section at the end of this cookbook explains what every element +of the figure means. + ## More Complex Lens Models The API above can be easily extended to compose lens models where each galaxy has multiple light or mass profiles: @@ -96,8 +123,8 @@ The API above can be easily extended to compose lens models where each galaxy ha ```python # Lens: -bulge = af.Model(al.lp.Sersic) -disk = af.Model(al.lp.Exponential) +bulge = af.Model(al.lp_linear.Sersic) +disk = af.Model(al.lp_linear.Exponential) mass = af.Model(al.mp.Isothermal) shear = af.Model(al.mp.ExternalShear) @@ -113,14 +140,21 @@ lens = af.Model( # Source: -bulge = af.Model(al.lp.SersicCore) -disk = af.Model(al.lp.Exponential) +bulge = af.Model(al.lp_linear.SersicCore) +disk = af.Model(al.lp_linear.ExponentialCore) source = af.Model(al.Galaxy, redshift=1.0, bulge=bulge, disk=disk) # Overall Lens Model: model = af.Collection(galaxies=af.Collection(lens=lens, source=source)) + +af.ModelPlotter(model).figure() +``` + +```{image} https://raw.githubusercontent.com/PyAutoLabs/PyAutoLens/main/docs/general/images/model_cookbook/complex_lens.png +:alt: A lens model whose lens and source galaxies each have a bulge and a disk, with an external shear. +:width: 600 ``` The use of the words `bulge`, `disk`, `mass` and `shear` above are arbitrary. They can be replaced with any name you @@ -129,7 +163,7 @@ like, e.g. `bulge_0`, `bulge_1`, `mass_0`, `mass_1`, and the model will still be The API can also be extended to compose lens models where there are multiple galaxies: ```python -bulge = af.Model(al.lp.Sersic) +bulge = af.Model(al.lp_linear.Sersic) mass = af.Model(al.mp.Isothermal) lens_0 = af.Model( @@ -139,7 +173,7 @@ lens_0 = af.Model( mass=mass, ) -bulge = af.Model(al.lp.Sersic) +bulge = af.Model(al.lp_linear.Sersic) mass = af.Model(al.mp.Isothermal) lens_1 = af.Model( @@ -151,13 +185,13 @@ lens_1 = af.Model( # Source 0: -bulge = af.Model(al.lp.SersicCore) +bulge = af.Model(al.lp_linear.SersicCore) source_0 = af.Model(al.Galaxy, redshift=1.0, bulge=bulge) # Source 1 : -bulge = af.Model(al.lp.SersicCore) +bulge = af.Model(al.lp_linear.SersicCore) source_1 = af.Model(al.Galaxy, redshift=1.0, bulge=bulge) @@ -171,10 +205,21 @@ model = af.Collection( source_1=source_1 ) ) + +af.ModelPlotter(model).figure() +``` + +```{image} https://raw.githubusercontent.com/PyAutoLabs/PyAutoLens/main/docs/general/images/model_cookbook/four_galaxies.png +:alt: A lens model with two lens galaxies and two source galaxies. +:width: 600 ``` +The two lens galaxies are identical in structure, and so are the two source galaxies, so the figure draws each pair +once inside a dashed plate badged `2 components` rather than drawing four cards. Their parameters are badged +`independent`: two separate priors with the same configuration, which is not the same thing as one shared prior. + The above lens model consists of only two planes (an image-plane and source-plane), but has four galaxies in total. -This is because the lens galaxies have the same redshift and the souece galaxies have the same redshift. +This is because the lens galaxies have the same redshift and the source galaxies have the same redshift. If we gave one of the lens galaxies a different redshift, it would be included in a third plane, and the model would perform multi-plane ray tracing when the model-fit is performed. @@ -191,8 +236,8 @@ removing the comments reading Lens / Source / Overall Lens Model to make the cod lens = af.Model( al.Galaxy, redshift=0.5, - bulge=al.lp.Sersic, - disk=al.lp.Sersic, + bulge=al.lp_linear.Sersic, + disk=al.lp_linear.Sersic, mass=al.mp.Isothermal, shear=al.mp.ExternalShear, ) @@ -200,27 +245,56 @@ lens = af.Model( source = af.Model( al.Galaxy, redshift=1.0, - bulge=al.lp.SersicCore, - disk=al.lp.Exponential + bulge=al.lp_linear.SersicCore, + disk=al.lp_linear.ExponentialCore, ) model = af.Collection(galaxies=af.Collection(lens=lens, source=source)) + +af.ModelPlotter(model).figure() ``` +```{image} https://raw.githubusercontent.com/PyAutoLabs/PyAutoLens/main/docs/general/images/model_cookbook/concise.png +:alt: The same lens and source model composed with the concise API. +:width: 600 +``` + +The concise API is a shorthand for writing a model, not a different model. Note that this stage gives the lens galaxy +a `Sersic` bulge *and* a `Sersic` disk, so the two collapse into a single dashed plate badged `2 components`: the +figure groups repeated sibling components rather than drawing two identical cards, and its `independent` badges say +that each of the two has its own priors rather than sharing one. + ## Prior Customization -We can customize the priors of the lens model component individual parameters as follows: +We can customize the priors of the lens model component individual parameters, using the following three types +of priors: + +`UniformPrior`: The values of a parameter are randomly drawn between a `lower_limit` and `upper_limit`. For example, +the effective radius of elliptical Sersic profiles typically assumes a uniform prior between 0.0" and 30.0". + +`LogUniformPrior`: Like a `UniformPrior` this randomly draws values between a `lower_limit` and `upper_limit`, but the +values are drawn from a distribution with base 10. This is used for the `intensity` of a light profile, as the +luminosity of galaxies follows a log10 distribution. + +`GaussianPrior`: The values of a parameter are randomly drawn from a Gaussian distribution with a `mean` and width +`sigma`. For example, the $y$ and $x$ centre values in a light profile typically assume a mean of 0.0" and a +sigma of 0.3", indicating that we expect the profile centre to be located near the centre of the image. ```python # Lens: -bulge = af.Model(al.lp.Sersic) -bulge.intensity = af.LogUniformPrior(lower_limit=1e-4, upper_limit=1e4) -bulge.sersic_index = af.GaussianPrior(mean=4.0, sigma=1.0, lower_limit=1.0, upper_limit=8.0) +bulge = af.Model(al.lp_linear.Sersic) +bulge.sersic_index = af.TruncatedGaussianPrior( + mean=4.0, sigma=1.0, lower_limit=1.0, upper_limit=8.0 +) mass = af.Model(al.mp.Isothermal) -mass.centre.centre_0 = af.GaussianPrior(mean=0.0, sigma=0.1, lower_limit=-0.5, upper_limit=0.5) -mass.centre.centre_1 = af.GaussianPrior(mean=0.0, sigma=0.1, lower_limit=-0.5, upper_limit=0.5) +mass.centre.centre_0 = af.TruncatedGaussianPrior( + mean=0.0, sigma=0.1, lower_limit=-0.5, upper_limit=0.5 +) +mass.centre.centre_1 = af.TruncatedGaussianPrior( + mean=0.0, sigma=0.1, lower_limit=-0.5, upper_limit=0.5 +) mass.einstein_radius = af.UniformPrior(lower_limit=0.0, upper_limit=8.0) lens = af.Model( @@ -232,16 +306,34 @@ lens = af.Model( # Source -bulge = af.Model(al.lp.Sersic) +bulge = af.Model(al.lp_linear.SersicCore) source = af.Model(al.Galaxy, redshift=1.0, bulge=bulge) -source.effective_radius = af.GaussianPrior(mean=0.1, sigma=0.05, lower_limit=0.0, upper_limit=1.0) +source.effective_radius = af.TruncatedGaussianPrior( + mean=0.1, sigma=0.05, lower_limit=0.0, upper_limit=1.0 +) # Overall Lens Model: model = af.Collection(galaxies=af.Collection(lens=lens, source=source)) + +af.ModelPlotter(model).figure() ``` +```{image} https://raw.githubusercontent.com/PyAutoLabs/PyAutoLens/main/docs/general/images/model_cookbook/prior_custom.png +:alt: A lens model whose sersic index, mass centre and einstein radius priors have been customized. +:width: 600 +``` + +Customizing a prior does not change a parameter's state — a parameter with a customized prior is still sampled — so +the map does not change when only the numbers do. Print `model.info`, or call +`af.ModelPlotter(model).figure(detail="priors")`, to read the prior on each parameter. + +The figure does show one thing the code above is easy to misread, though: the last line sets `effective_radius` on the +**source galaxy**, not on the source galaxy's `bulge`. A `Galaxy` accepts any attribute, so this adds a *new* free +`effective_radius` parameter to the `source · Galaxy` card, drawn above its `bulge` card, rather than overriding the +bulge's own prior. Writing `source.bulge.effective_radius = ...` is what customizes the bulge. + ## Model Customization We can customize the lens model parameters in a number of different ways, as shown below: @@ -249,8 +341,8 @@ We can customize the lens model parameters in a number of different ways, as sho ```python # Lens: -bulge = af.Model(al.lp.Sersic) -disk = af.Model(al.lp.Exponential) +bulge = af.Model(al.lp_linear.Sersic) +disk = af.Model(al.lp_linear.Exponential) # Parameter Pairing: Pair the centre of the bulge and disk together, reducing # the complexity of non-linear parameter space by N = 2 @@ -283,8 +375,8 @@ lens = af.Model( # Source: -bulge = af.Model(al.lp.SersicCore) -disk = af.Model(al.lp.Exponential) +bulge = af.Model(al.lp_linear.SersicCore) +disk = af.Model(al.lp_linear.ExponentialCore) source = af.Model(al.Galaxy, redshift=1.0, bulge=bulge, disk=disk) @@ -294,13 +386,30 @@ model = af.Collection(galaxies=af.Collection(lens=lens, source=source)) # Assert that the effective radius of the bulge is larger than that of the disk. # (Assertions can only be added at the end of model composition, after all components -# have been bright together in a `Collection`. +# have been brought together in a `Collection`. model.add_assertion(model.galaxies.lens.bulge.effective_radius > model.galaxies.lens.disk.effective_radius) # Assert that the Einstein Radius is below 3.0": model.add_assertion(model.galaxies.lens.mass.einstein_radius < 3.0) + +af.ModelPlotter(model).figure() +``` + +```{image} https://raw.githubusercontent.com/PyAutoLabs/PyAutoLens/main/docs/general/images/model_cookbook/model_custom.png +:alt: A customized lens model showing a paired centre, a fixed sersic index, offset relations and two assertions. +:width: 600 ``` +This is the stage where the figure earns its keep: the paired `centre` is drawn once on its owner badged +`shared x2`, with a blue link from the `disk` that reuses it (`-> bulge.centre`); the fixed `sersic_index` is a grey +pill; and each assertion is drawn as a compact dashed-orange label naming both of its operands, rather than as a line +traced across the figure. + +The two offset `centre` components of the mass profile are a `relation`, and a relation on a *scalar* parameter is +drawn with its defining expression on the pill (the galaxy cookbook's `effective_radius` offset shows this). A +relation on one component of a **tuple** parameter such as `centre` is not yet annotated on the tuple's single pill, +so the mass `centre` above is drawn as an ordinary sampled pill; `model.info` is the place to read it. + ## Redshift Free The redshift of a galaxy can be treated as a free parameter in the model-fit by using the following API: @@ -314,8 +423,23 @@ lens = af.Model( redshift=redshift, mass=al.mp.Isothermal ) + +source = af.Model(al.Galaxy, redshift=1.0, bulge=al.lp_linear.SersicCore) + +model = af.Collection(galaxies=af.Collection(lens=lens, source=source)) + +af.ModelPlotter(model).figure() +``` + +```{image} https://raw.githubusercontent.com/PyAutoLabs/PyAutoLens/main/docs/general/images/model_cookbook/redshift_free.png +:alt: A lens model whose lens galaxy redshift is a free parameter, drawn as an ordinary pill on the galaxy card. +:width: 600 ``` +A free redshift is an ordinary sampled parameter, so it is drawn as an ordinary pill on the galaxy's card. A *fixed* +redshift is not drawn as a pill at all: it is printed under the galaxy's header as `redshift = 1.0`, which is the one +number the figure shows on the map rather than in the legend. The source galaxy above shows this. + The model-fit will automatically enable multi-plane ray tracing and alter the ordering of the planes depending on the redshifts of the galaxies. @@ -324,9 +448,101 @@ impact the model-fit. You should therefore never make the redshifts free if you system. This is because lensing calculations can be defined in arc-second coordinates, which do not change as a function of redshift. -Redshifts should be made free when modeling three or more planes, as the mulit-plane ray-tracing calculations have an +Redshifts should be made free when modeling three or more planes, as the multi-plane ray-tracing calculations have an obvious dependence on the redshifts of the galaxies which could be inferred by the model-fit. +## Solved Parameters + +Some parameters of a lens model are neither sampled by the non-linear search nor fixed by the user: they are **solved +for during the fit**, at every likelihood evaluation. The model below contains all three of them: + +```python +# Lens: a linear light profile, whose `intensity` is solved for by the inversion. + +lens = af.Model( + al.Galaxy, + redshift=0.5, + bulge=af.Model(al.lp_linear.Sersic), + mass=af.Model(al.mp.Isothermal), +) + +# Source: a pixelization, whose source `reconstruction` is solved for by the inversion. + +pixelization = af.Model( + al.Pixelization, + mesh=af.Model(al.mesh.Delaunay, pixels=500, zeroed_pixels=0), + regularization=af.Model(al.reg.ConstantSplit), +) + +source = af.Model(al.Galaxy, redshift=1.0, pixelization=pixelization) + +# Point source: a `PointSolved`, whose `centre` is solved for analytically. + +point_source = af.Model(al.Galaxy, redshift=1.0, point=af.Model(al.ps.PointSolved)) + +model = af.Collection( + galaxies=af.Collection(lens=lens, source=source, point_source=point_source) +) + +af.ModelPlotter(model).figure() +``` + +```{image} https://raw.githubusercontent.com/PyAutoLabs/PyAutoLens/main/docs/general/images/model_cookbook/solved_parameters.png +:alt: A lens model showing the three solved parameters and one missing configuration value. +:width: 600 +``` + +Three parameters in the figure carry a dashed `solved` pill: + +- the `intensity` of every linear light profile (`al.lp_linear.*`, and every member of an `al.lp_basis.Basis` built + from them), which the inversion solves for by linear algebra. A `Basis` declares no solved amplitude of its own: + the solved intensities belong to its member profiles, one per member. +- the `reconstruction` of an `al.Pixelization`, which is the solved surface brightness of every source pixel. The only + sampled parameters of a pixelization are its regularization coefficients. +- the `centre` of an `al.ps.PointSolved`, which is solved for analytically by the `*Solved` point-source fit classes. + Without the annotation this component would look like an empty box; `al.ps.Point`, by contrast, samples its centre + and draws an ordinary free `centre` pill. + +**These three have no counterpart in `model.info` at all.** They are additional information the figure supplies, which +is why they are drawn dashed and named in the legend as *solved during fitting*, and why the figure's footer counts +them separately from the sampled parameters. + +A fourth pill in the figure reads `areas_factor · missing`. That is a different thing again: `areas_factor` is a +parameter of the `Delaunay` mesh for which no prior is configured, so the model cannot be fitted until one is supplied +(via `config/priors` or by setting it in the script). It is **unset configuration** — not solved, and not absent from +the model. Absence from the figure would read as absence from the model, so `missing` is a state of its own. + +The red pill above reflects PyAutoLens' own default configuration, which ships no prior for `areas_factor`. The +`autolens_workspace` supplies one in `config/priors/mesh/delaunay.yaml`, as a constant `0.5`, so the same code run from +the workspace draws an ordinary grey fixed pill instead: `missing` is what you see when your own configuration has a +gap. + +## Reading the Figure + +The figure is the **map**; `model.info` is the **legend**. The map shows the shape of the model — which components +contain which, which parameters are shared, fixed, related or solved — and the legend gives the numbers: the prior on +every parameter, and the exact value of every fixed one. + +The correspondence between them is an explicit contract: + +> Every displayed model element resolves to its corresponding path or grouped paths in `model.info`; every omission and +> every added annotation (`solved`, `missing`) is explicit. + +It is deliberately **not** a line-for-line correspondence, because the two group the model differently: + +> The figure partitions **by component**: the MGE model shows two `30 components` plates, split because each basis +> holds its own `ell_comps` pair. `model.info` groups **per parameter**: it prints a single `0 - 59` block for +> `centre` spanning both figure plates, two separate ellipticity blocks, and individual `sigma` blocks. + +Two further reading notes: + +- The background tints follow **nesting depth only**. They are decorative: they help you see which card sits inside + which, and they carry no information about the class family of a component. A light profile, a mass profile and a + pixelization at the same depth are tinted identically. +- A plate badge reads `30 components` and means repetition; a parameter badge reads `shared across group` and means + one prior reused. These are different statements, so they never share a symbol. A grey `fixed, varies by member` + pill means every member has its own fixed value, not one value common to all. + ## Available Model Components The light profiles, mass profiles and other components that can be used for lens modeling are given at the following @@ -404,9 +620,25 @@ for i, (centre, luminosity) in enumerate(zip(centres, luminosities)): members on scaling relations) but change the source strategy: many point-source multiple-image position datasets, each at its own redshift, fitted via ``AnalysisPoint``. +The full group-scale composition — one main lens galaxy with an external shear, a `SersicCore` source and one +extra galaxy per detected group member, each with a `SersicSph` light profile and a `dPIEMassSph` mass profile whose +`centre` is fixed to its measured position — is the script +`autolens_workspace/scripts/group/modeling.py`. Drawn as a figure, with eight extra galaxies, it looks like this: + +```{image} https://raw.githubusercontent.com/PyAutoLabs/PyAutoLens/main/docs/general/images/model_cookbook/acceptance_group_scale.png +:alt: A group-scale lens model, whose eight extra galaxies collapse into a single plate of eight components. +:width: 600 +``` + +The eight extra galaxies are identical in structure, so the figure draws them once inside a dashed plate badged +`8 components` rather than eight times. Their `centre` pills read `fixed, varies by member` — each galaxy has its own +fixed centre, not one centre common to all — and their `sigma` priors are badged `independent`, because eight separate +priors with the same configuration is not the same thing as one shared prior. + The following example notebooks show each regime's full model composition: + @@ -416,6 +648,17 @@ The following example notebooks show each regime's full model composition: Features such as the Multi Gaussian Expansion (MGE) and shapelets compose models consisting of 50 - 500+ light profiles. +An MGE lens light model of 2 x 30 linear Gaussians, fitted alongside a pixelized source, draws like this: + +```{image} https://raw.githubusercontent.com/PyAutoLabs/PyAutoLens/main/docs/general/images/model_cookbook/acceptance_mge_pixelized.png +:alt: An MGE lens light model of two bases of thirty linear Gaussians, with a pixelized Delaunay source. +:width: 600 +``` + +The 60 Gaussians collapse into two plates of `30 components`, split because each basis holds its own `ell_comps` +pair while all 60 share one `centre`. Every member's `intensity` is `solved`, the source's `reconstruction` is +`solved`, and the mesh's unconfigured `areas_factor` is `missing`. + The following example notebooks show how to compose and fit these models: diff --git a/docs/overview/images/overview_3/mge_model.png b/docs/overview/images/overview_3/mge_model.png new file mode 100644 index 0000000000..68e89668ff Binary files /dev/null and b/docs/overview/images/overview_3/mge_model.png differ diff --git a/docs/overview/images/overview_3/pixelization_model.png b/docs/overview/images/overview_3/pixelization_model.png new file mode 100644 index 0000000000..9cf58cd54c Binary files /dev/null and b/docs/overview/images/overview_3/pixelization_model.png differ diff --git a/docs/overview/images/overview_3/point_model.png b/docs/overview/images/overview_3/point_model.png new file mode 100644 index 0000000000..07f6a69298 Binary files /dev/null and b/docs/overview/images/overview_3/point_model.png differ diff --git a/docs/overview/overview_3_features.md b/docs/overview/overview_3_features.md index 3e03120086..ff219c80d5 100644 --- a/docs/overview/overview_3_features.md +++ b/docs/overview/overview_3_features.md @@ -43,6 +43,15 @@ galaxy: :width: 600 ``` +The lens model a pixelized source composes is drawn below: the source's only sampled parameter is its +regularization coefficient, while the `reconstruction` — the surface brightness of every source pixel — is solved +for during the fit. + +```{image} https://raw.githubusercontent.com/PyAutoLabs/PyAutoLens/main/docs/overview/images/overview_3/pixelization_model.png +:alt: A pixelized source lens model, whose source reconstruction is solved for during the fit. +:width: 600 +``` + A complete overview of pixelized source reconstructions can be found at `notebooks/overview/overview_5_pixelizations.ipynb`. @@ -85,6 +94,15 @@ higher resolutions to determine the source's exact locations in the image-plane: :width: 400 ``` +The lens model this composes is drawn below: the source galaxy carries a single `al.ps.Point` component whose +`centre` is a sampled two-dimensional parameter. (Using `al.ps.PointSolved` instead solves that centre analytically, +and the figure draws it as a dashed `solved` pill.) + +```{image} https://raw.githubusercontent.com/PyAutoLabs/PyAutoLens/main/docs/overview/images/overview_3/point_model.png +:alt: A point-source lens model, with an isothermal lens and a point source whose centre is sampled. +:width: 600 +``` + Note that the image positions above include the fifth central image of the strong lens, which is often not seen in strong lens imaging data. It is easy to disable this image in the point source modeling. @@ -133,6 +151,14 @@ An MGE decomposes the light of a galaxy into tens or hundreds of two dimensional In the image above, 30 Gaussians are shown, where their sizes go from below the pixel scale (in order to resolve point emission) to beyond the size of the galaxy (to capture its extended emission). +The lens model this composes is drawn below: the 60 Gaussians of its two bases collapse into two plates badged +`30 components`, whose `intensity` parameters are solved for by the inversion rather than sampled. + +```{image} https://raw.githubusercontent.com/PyAutoLabs/PyAutoLens/main/docs/overview/images/overview_3/mge_model.png +:alt: The MGE lens model, whose sixty Gaussians are drawn as two plates of thirty components. +:width: 600 +``` + An MGE is an extremely powerful way to model and subtract the light of the foreground lens galaxy in strong lens imaging, and makes it possible to model the stellar mass of the lens galaxy in a way that is tied to its light. diff --git a/test_autolens/model_figure/__init__.py b/test_autolens/model_figure/__init__.py new file mode 100644 index 0000000000..e69de29bb2 diff --git a/test_autolens/model_figure/lens_models.py b/test_autolens/model_figure/lens_models.py new file mode 100644 index 0000000000..4be78bf128 --- /dev/null +++ b/test_autolens/model_figure/lens_models.py @@ -0,0 +1,181 @@ +""" +The three lens acceptance models of the ``model-figures`` epic, built from the +**real** PyAutoLens classes. + +Phases 1 and 2 rendered these three models from *structural doubles* inside +``test_autofit`` (classes with matching constructor signatures, so PyAutoFit +could render them without importing PyAutoLens). Phase 3 replaces the doubles +with the real thing: every model below is the composition of a shipped +``autolens_workspace`` script, so the figure a reader sees in the docs is the +model the workspace composes. + +* :func:`simple_lens` -- ``scripts/imaging/modeling.py`` +* :func:`mge_pixelized` -- ``scripts/imaging/features/multi_gaussian_expansion/modeling.py`` + with the pixelized source of ``scripts/imaging/features/pixelization/modeling.py`` + (no shipped script has both, so the two are assembled here) +* :func:`group_scale` -- ``scripts/group/modeling.py`` with eight synthetic + extra-galaxy centres (the shipped dataset has one; the epic's acceptance case + is eight) + +The builders take no arguments and load no data, so they are importable from the +docs ``make_figures.py`` script as well as from the tests. +""" + +import numpy as np + +import autofit as af +import autolens as al + +__all__ = [ + "simple_lens", + "mge_pixelized", + "group_scale", + "EXTRA_GALAXY_CENTRES", +] + +#: The eight synthetic extra-galaxy centres of the group-scale acceptance case. +EXTRA_GALAXY_CENTRES = [ + (1.0, 1.0), + (-1.0, 1.0), + (1.0, -1.0), + (-1.0, -1.0), + (2.0, 0.0), + (-2.0, 0.0), + (0.0, 2.0), + (0.0, -2.0), +] + + +def simple_lens(): + """ + Acceptance case (a) -- the simple lens of ``scripts/imaging/modeling.py``. + + A lens galaxy at ``z=0.5`` with a ``Sersic`` bulge, an ``Isothermal`` mass + profile and an ``ExternalShear``, and a source galaxy at ``z=1.0`` with a + ``SersicCore`` bulge. + """ + bulge = af.Model(al.lp.Sersic) + mass = af.Model(al.mp.Isothermal) + shear = af.Model(al.mp.ExternalShear) + + lens = af.Model(al.Galaxy, redshift=0.5, bulge=bulge, mass=mass, shear=shear) + + bulge = af.Model(al.lp.SersicCore) + + source = af.Model(al.Galaxy, redshift=1.0, bulge=bulge) + + return af.Collection(galaxies=af.Collection(lens=lens, source=source)) + + +def mge_basis(total_gaussians: int = 30, gaussian_per_basis: int = 2): + """ + The MGE ``Basis`` of ``features/multi_gaussian_expansion/modeling.py``. + + ``gaussian_per_basis`` bases of ``total_gaussians`` linear Gaussians each. + All 60 Gaussians share one ``centre_0`` and one ``centre_1`` prior; the + Gaussians of each basis share that basis's own ``ell_comps`` pair (which is + what splits the figure into two plates); every ``sigma`` is fixed to its own + log-spaced value. + """ + mask_radius = 3.0 + pixel_scale = 0.1 + + log10_sigma_list = np.linspace( + np.log10(pixel_scale / 10.0), np.log10(mask_radius), total_gaussians + ) + + centre_0 = af.UniformPrior(lower_limit=-0.1, upper_limit=0.1) + centre_1 = af.UniformPrior(lower_limit=-0.1, upper_limit=0.1) + + bulge_gaussian_list = [] + + for _ in range(gaussian_per_basis): + gaussian_list = af.Collection( + af.Model(al.lp_linear.Gaussian) for _ in range(total_gaussians) + ) + + for i, gaussian in enumerate(gaussian_list): + gaussian.centre.centre_0 = centre_0 + gaussian.centre.centre_1 = centre_1 + gaussian.ell_comps = gaussian_list[0].ell_comps + gaussian.sigma = 10 ** log10_sigma_list[i] + + bulge_gaussian_list += gaussian_list + + return af.Model(al.lp_basis.Basis, profile_list=bulge_gaussian_list) + + +def mge_pixelized(): + """ + Acceptance case (b) -- MGE lens light, pixelized source. + + The lens is the 2 x 30 MGE of + ``features/multi_gaussian_expansion/modeling.py`` plus an ``Isothermal`` and + an ``ExternalShear``; the source is the ``Pixelization`` of + ``features/pixelization/modeling.py`` built on a ``Delaunay`` mesh, whose + ``areas_factor`` has no prior configured and is therefore ``missing``. + """ + lens = af.Model( + al.Galaxy, + redshift=0.5, + bulge=mge_basis(), + mass=af.Model(al.mp.Isothermal), + shear=af.Model(al.mp.ExternalShear), + ) + + pixelization = af.Model( + al.Pixelization, + mesh=af.Model(al.mesh.Delaunay, pixels=500, zeroed_pixels=0), + regularization=af.Model(al.reg.ConstantSplit), + ) + + source = af.Model(al.Galaxy, redshift=1.0, pixelization=pixelization) + + return af.Collection(galaxies=af.Collection(lens=lens, source=source)) + + +def group_scale(centres=None): + """ + Acceptance case (e) -- the group-scale model of ``scripts/group/modeling.py``. + + One main lens galaxy (``Sersic`` + ``Isothermal`` + ``ExternalShear``), a + ``SersicCore`` source, and eight extra galaxies whose ``dPIEMassSph`` + centres are fixed to their own measured positions -- a plate of eight whose + ``centre`` is ``fixed, varies by member`` and whose ``sigma`` priors are + independent. + """ + centres = EXTRA_GALAXY_CENTRES if centres is None else centres + + lens = af.Model( + al.Galaxy, + redshift=0.5, + bulge=af.Model(al.lp.Sersic), + mass=af.Model(al.mp.Isothermal), + shear=af.Model(al.mp.ExternalShear), + ) + + extra_galaxies_list = [] + + for centre in centres: + bulge = af.Model(al.lp.SersicSph) + + mass = af.Model(al.mp.dPIEMassSph) + mass.centre = centre + mass.sigma = af.UniformPrior(lower_limit=0.0, upper_limit=300.0) + mass.r_core = 0.0 + mass.r_cut = 10.0 + mass.redshift_object = 0.5 + mass.redshift_source = 1.0 + mass.H0 = 67.66 + mass.Om0 = 0.30966 + + extra_galaxies_list.append( + af.Model(al.Galaxy, redshift=0.5, bulge=bulge, mass=mass) + ) + + source = af.Model(al.Galaxy, redshift=1.0, bulge=af.Model(al.lp.SersicCore)) + + return af.Collection( + galaxies=af.Collection(lens=lens, source=source), + extra_galaxies=af.Collection(extra_galaxies_list), + ) diff --git a/test_autolens/model_figure/test_acceptance_lens.py b/test_autolens/model_figure/test_acceptance_lens.py new file mode 100644 index 0000000000..346268b608 --- /dev/null +++ b/test_autolens/model_figure/test_acceptance_lens.py @@ -0,0 +1,313 @@ +""" +The ``model-figures`` epic's three lens acceptance cases, over real classes. + +Phases 1 and 2 rendered these from structural doubles inside ``test_autofit``. +This module renders the real thing (:mod:`test_autolens.model_figure.lens_models`) +and asserts, in this order: + +1. the **structural invariants** -- what the figure must say about the model, + which is what the epic and its independent review are actually about; +2. the **counts**, pinned as measured, with the epic's own numbers recorded in a + ``# epic target:`` comment wherever the two differ (the epic's numbers were + measured on the doubles, whose priors differ from the real classes'); +3. the **rendered PNG** -- inside the width budget, and identical across two + independent builds of the same model. +""" + +import matplotlib + +matplotlib.use("Agg") + +import pytest + +import autofit as af + +from test_autolens.model_figure import lens_models + + +# --------------------------------------------------------------------------- +# helpers +# --------------------------------------------------------------------------- + + +def walk_nodes(node): + yield node + for child in node.children: + yield from walk_nodes(child) + + +def rows_of(spec): + rows = {} + for node in walk_nodes(spec.root): + for row in node.rows: + rows[".".join(row.path)] = row + return rows + + +def plates_of(spec): + return [node for node in walk_nodes(spec.root) if node.plate is not None] + + +def cards_of(presentation): + return list(presentation.walk()) + + +def pill_named(card, name): + for pill in card.pills: + if pill.text == name or pill.text.startswith(f"{name} ·"): + return pill + return None + + +def card_of(presentation, key): + for card in presentation.walk(): + if card.key == key: + return card + raise AssertionError( + f"no card keyed {key!r}: {[card.key for card in presentation.walk()]}" + ) + + +def png_size(path): + """``(width, height)`` in pixels, read from the PNG IHDR chunk.""" + import struct + + with open(path, "rb") as f: + header = f.read(24) + + assert header[:8] == b"\x89PNG\r\n\x1a\n" + + return struct.unpack(">II", header[16:24]) + + +# --------------------------------------------------------------------------- +# (a) the simple lens +# --------------------------------------------------------------------------- + + +def test__simple_lens__six_cards_in_model_info_order(): + model = lens_models.simple_lens() + + presentation = af.ModelPlotter(model).presentation() + + component_titles = [ + card.title for card in cards_of(presentation) if card.kind == "model" + ] + + assert component_titles == [ + "lens · Galaxy", + "bulge · Sersic", + "mass · Isothermal", + "shear · ExternalShear", + "source · Galaxy", + "bulge · SersicCore", + ] + + +def test__simple_lens__counts(): + model = lens_models.simple_lens() + + spec = af.GraphSpec.from_model(model) + + # 8 components = the six drawn cards plus the root `Collection` and the + # `galaxies` `Collection` frame, which are nodes of the tree too. + assert spec.counts["components"] == 8 + assert spec.counts["unique_sampled_scalars"] == model.prior_count == 21 + assert spec.counts["shared_priors"] == 0 + assert spec.counts["solved"] == 0 + assert spec.counts["missing"] == 0 + # `SersicCore`'s `radius_break`, `gamma` and `alpha`, and the two redshifts. + assert spec.counts["fixed_leaf_slots"] == 5 + + +# --------------------------------------------------------------------------- +# (b) the MGE lens light + pixelized source +# --------------------------------------------------------------------------- + + +@pytest.fixture(name="mge_spec", scope="module") +def make_mge_spec(): + return af.GraphSpec.from_model(lens_models.mge_pixelized()) + + +def test__mge__two_plates_of_thirty_split_on_ell_comps(mge_spec): + plates = [node for node in plates_of(mge_spec) if node.plate.count == 30] + + assert len(plates) == 2 + + # The two plates exist *because* each basis holds its own `ell_comps` pair: + # the shared `centre` priors are carried by all 60 members and so cannot + # discriminate, while the `ell_comps` priors differ between the bases. + ell_comps_ids = { + tuple( + component.prior_id + for row in plate.rows + if row.name == "ell_comps" + for component in row.components + ) + for plate in plates + } + + assert len(ell_comps_ids) == 2 + + +def test__mge__every_member_intensity_is_solved(mge_spec): + plates = [node for node in plates_of(mge_spec) if node.plate.count == 30] + + for plate in plates: + intensity = next(row for row in plate.rows if row.name == "intensity") + + assert intensity.sampling == "solved" + assert intensity.in_model_info is False + + +def test__mge__basis_itself_declares_no_solved_intensity(mge_spec): + basis = next(node for node in walk_nodes(mge_spec.root) if node.cls_name == "Basis") + + assert [row.name for row in basis.rows if row.name == "intensity"] == [] + + +def test__mge__shared_centre_is_badged_on_its_owner(): + model = lens_models.mge_pixelized() + + presentation = af.ModelPlotter(model).presentation() + + badges = { + pill.badge + for card in cards_of(presentation) + for pill in card.pills + if pill.text.startswith("centre") and pill.badge is not None + } + + assert "shared across group" in badges + + +def test__mge__areas_factor_is_missing_on_the_mesh(mge_spec): + row = rows_of(mge_spec)["galaxies.source.pixelization.mesh.areas_factor"] + + assert row.sampling == "missing" + assert row.prior_cls_name == "ConfigException" + + +def test__mge__reconstruction_is_solved_on_the_pixelization_card(): + model = lens_models.mge_pixelized() + + presentation = af.ModelPlotter(model).presentation() + + card = card_of(presentation, "galaxies/source/pixelization") + pill = pill_named(card, "reconstruction") + + assert pill is not None + assert pill.state == "solved" + assert pill.text == "reconstruction · solved" + + +def test__mge__counts(mge_spec): + model = lens_models.mge_pixelized() + + assert mge_spec.counts["unique_sampled_scalars"] == model.prior_count == 14 + # epic target: 16 unique sampled scalars. The epic measured the phase-1 + # structural double, whose source was a `Sersic` (7 priors); this model's + # source is pixelized, so its only sampled parameter is the `ConstantSplit` + # coefficient (1 prior). 13 lens priors + 1 = 14. + assert mge_spec.counts["components"] == 13 + assert mge_spec.counts["components_raw"] == 71 + assert mge_spec.counts["plates"] == 2 + # `centre_0`, `centre_1`, and each basis's own `ell_comps` pair. + assert mge_spec.counts["shared_priors"] == 6 + assert mge_spec.counts["fixed_leaf_slots"] == 64 + # 60 member intensities plus the source reconstruction. + assert mge_spec.counts["solved"] == 61 + # `Delaunay.areas_factor`, which has no prior configured. + assert mge_spec.counts["missing"] == 1 + + +# --------------------------------------------------------------------------- +# (e) the group-scale model +# --------------------------------------------------------------------------- + + +@pytest.fixture(name="group_spec", scope="module") +def make_group_spec(): + return af.GraphSpec.from_model(lens_models.group_scale()) + + +def test__group__one_plate_of_eight(group_spec): + plates = plates_of(group_spec) + + assert len(plates) == 1 + assert plates[0].plate.count == 8 + + +def test__group__member_centres_are_fixed_and_vary_by_member(): + model = lens_models.group_scale() + + presentation = af.ModelPlotter(model).presentation() + + texts = [pill.text for pill in presentation.pills()] + + assert "centre · fixed, varies by member" in texts + + +def test__group__member_sigma_priors_are_independent(): + model = lens_models.group_scale() + + presentation = af.ModelPlotter(model).presentation() + + sigma = [pill for pill in presentation.pills() if pill.text.startswith("sigma")] + + assert sigma + assert any(pill.badge == "independent" for pill in sigma) + + +def test__group__counts(group_spec): + model = lens_models.group_scale() + + assert group_spec.counts["unique_sampled_scalars"] == model.prior_count == 69 + assert group_spec.counts["components"] == 12 + assert group_spec.counts["components_raw"] == 33 + assert group_spec.counts["plates"] == 1 + assert group_spec.counts["shared_priors"] == 0 + # Eight extra galaxies x (fixed centre pair + r_core, r_cut, the two + # redshifts, H0, Om0) plus the main lens and source redshifts. + assert group_spec.counts["fixed_leaf_slots"] == 77 + assert group_spec.counts["solved"] == 0 + assert group_spec.counts["missing"] == 0 + + +# --------------------------------------------------------------------------- +# the renders +# --------------------------------------------------------------------------- + + +@pytest.mark.parametrize( + "builder, filename", + [ + (lens_models.simple_lens, "acceptance_simple_lens"), + (lens_models.mge_pixelized, "acceptance_mge_pixelized"), + (lens_models.group_scale, "acceptance_group_scale"), + ], +) +def test__renders_inside_the_width_budget(builder, filename, tmp_path): + af.ModelPlotter(builder()).figure( + path=str(tmp_path), filename=filename, format="png", width=14.0 + ) + + png = tmp_path / f"{filename}.png" + + assert png.exists() + + width, _ = png_size(png) + + assert width <= 1400 + + +@pytest.mark.parametrize( + "builder", + [lens_models.simple_lens, lens_models.mge_pixelized, lens_models.group_scale], +) +def test__layout_is_identical_across_two_builds_of_the_same_model(builder): + first = af.ModelPlotter(builder()).layout(width=14.0).to_dict() + second = af.ModelPlotter(builder()).layout(width=14.0).to_dict() + + assert first == second