diff --git a/.script_sizes.json b/.script_sizes.json index 06e478e3c..00e0a1669 100644 --- a/.script_sizes.json +++ b/.script_sizes.json @@ -1,17 +1,17 @@ { "scripts/cluster/__init__.py": 0, - "scripts/cluster/csv_api.py": 18602, + "scripts/cluster/csv_api.py": 18648, "scripts/cluster/lenstool/__init__.py": 0, - "scripts/cluster/lenstool/data.py": 15172, - "scripts/cluster/lenstool/modeling.py": 21932, - "scripts/cluster/lenstool/parameterization_mapping.py": 24021, - "scripts/cluster/likelihood_function.py": 37956, - "scripts/cluster/mass_parameterizations.py": 34990, - "scripts/cluster/mass_parameterizations_pyautolens.py": 20824, - "scripts/cluster/modeling.py": 32864, - "scripts/cluster/plot.py": 4443, + "scripts/cluster/lenstool/data.py": 15103, + "scripts/cluster/lenstool/modeling.py": 22297, + "scripts/cluster/lenstool/parameterization_mapping.py": 23956, + "scripts/cluster/likelihood_function.py": 39137, + "scripts/cluster/mass_parameterizations.py": 35001, + "scripts/cluster/mass_parameterizations_pyautolens.py": 17701, + "scripts/cluster/modeling.py": 32813, + "scripts/cluster/plot.py": 7381, "scripts/cluster/simulator.py": 33421, - "scripts/cluster/start_here.py": 23785, + "scripts/cluster/start_here.py": 24027, "scripts/group/__init__.py": 0, "scripts/group/data_preparation/__init__.py": 0, "scripts/group/data_preparation/start_here.py": 19804, @@ -86,20 +86,20 @@ "scripts/group/simulator.py": 13175, "scripts/group/slam.py": 36772, "scripts/group/source_science.py": 12692, - "scripts/group/start_here.py": 23475, + "scripts/group/start_here.py": 22849, "scripts/guides/__init__.py": 0, "scripts/guides/advanced/__init__.py": 0, "scripts/guides/advanced/add_a_profile.py": 28656, "scripts/guides/advanced/custom_analysis.py": 25045, - "scripts/guides/advanced/multi_plane.py": 34456, + "scripts/guides/advanced/multi_plane.py": 23386, "scripts/guides/advanced/over_sampling.py": 23958, "scripts/guides/advanced/over_sampling_chaining.py": 10482, "scripts/guides/coolest_interop.py": 5640, - "scripts/guides/data_structures.py": 18294, - "scripts/guides/galaxies.py": 10636, + "scripts/guides/data_structures.py": 15004, + "scripts/guides/galaxies.py": 8231, "scripts/guides/hpc/__init__.py": 0, "scripts/guides/hpc/example_cpu_and_gpu.py": 14065, - "scripts/guides/lens_calc.py": 25664, + "scripts/guides/lens_calc.py": 18780, "scripts/guides/modeling/__init__.py": 0, "scripts/guides/modeling/advanced/__init__.py": 0, "scripts/guides/modeling/advanced/expectation_propagation.py": 14753, @@ -107,7 +107,7 @@ "scripts/guides/modeling/advanced/hierarchical.py": 13049, "scripts/guides/modeling/bug_fix.py": 5769, "scripts/guides/modeling/chaining.py": 18024, - "scripts/guides/modeling/cookbook.py": 14384, + "scripts/guides/modeling/cookbook.py": 14400, "scripts/guides/modeling/customize.py": 11887, "scripts/guides/modeling/searches.py": 20117, "scripts/guides/modeling/slam_start_here.py": 22642, @@ -116,7 +116,7 @@ "scripts/guides/plot/searches.py": 23865, "scripts/guides/plot/start_here.py": 10152, "scripts/guides/plot/visuals.py": 7782, - "scripts/guides/point_source_pairing.py": 10038, + "scripts/guides/point_source_pairing.py": 17021, "scripts/guides/profiles/__init__.py": 0, "scripts/guides/profiles/light.py": 27532, "scripts/guides/profiles/light_and_mass_profiles.py": 27758, @@ -124,10 +124,10 @@ "scripts/guides/results/__init__.py": 0, "scripts/guides/results/_quick_fit.py": 5639, "scripts/guides/results/aggregator/__init__.py": 0, - "scripts/guides/results/aggregator/data_fitting.py": 7674, + "scripts/guides/results/aggregator/data_fitting.py": 7682, "scripts/guides/results/aggregator/galaxies_fits.py": 9258, "scripts/guides/results/aggregator/interferometer.py": 2355, - "scripts/guides/results/aggregator/models.py": 10146, + "scripts/guides/results/aggregator/models.py": 10154, "scripts/guides/results/aggregator/queries.py": 5694, "scripts/guides/results/aggregator/samples.py": 22417, "scripts/guides/results/aggregator/samples_via_aggregator.py": 22018, @@ -139,12 +139,12 @@ "scripts/guides/results/workflow/csv_make.py": 11839, "scripts/guides/results/workflow/fits_make.py": 13021, "scripts/guides/results/workflow/png_make.py": 11639, - "scripts/guides/tracer.py": 23319, + "scripts/guides/tracer.py": 20777, "scripts/guides/units/__init__.py": 0, "scripts/guides/units/cosmology.py": 11036, "scripts/guides/units/flux.py": 13965, "scripts/guides/units/mass_to_light_ratio_units.py": 5251, - "scripts/guides/using_jax.py": 9936, + "scripts/guides/using_jax.py": 22319, "scripts/imaging/__init__.py": 0, "scripts/imaging/data_preparation/__init__.py": 0, "scripts/imaging/data_preparation/examples/__init__.py": 0, @@ -243,13 +243,13 @@ "scripts/imaging/features/scaling_relation/likelihood_function.py": 12298, "scripts/imaging/features/scaling_relation/modeling.py": 17450, "scripts/imaging/features/scaling_relation/simulator.py": 11212, - "scripts/imaging/features/scaling_relation/slam.py": 31120, + "scripts/imaging/features/scaling_relation/slam.py": 32254, "scripts/imaging/features/simulator_manual_signal_to_noise_ratio.py": 7624, - "scripts/imaging/fit.py": 16848, + "scripts/imaging/fit.py": 16863, "scripts/imaging/likelihood_function.py": 21241, "scripts/imaging/modeling.py": 33768, - "scripts/imaging/plot.py": 6596, - "scripts/imaging/simulator.py": 19684, + "scripts/imaging/plot.py": 8544, + "scripts/imaging/simulator.py": 19656, "scripts/imaging/simulator_sample.py": 8120, "scripts/imaging/source_science.py": 10283, "scripts/imaging/start_here.py": 25845, @@ -276,9 +276,9 @@ "scripts/interferometer/features/datacube/delaunay.py": 9601, "scripts/interferometer/features/datacube/likelihood_function.py": 29864, "scripts/interferometer/features/datacube/modeling.py": 13043, - "scripts/interferometer/features/datacube/modeling_parametric.py": 9663, + "scripts/interferometer/features/datacube/modeling_parametric.py": 9679, "scripts/interferometer/features/datacube/simulator.py": 16873, - "scripts/interferometer/features/datacube/start_here.py": 15401, + "scripts/interferometer/features/datacube/start_here.py": 15409, "scripts/interferometer/features/extra_galaxies/__init__.py": 0, "scripts/interferometer/features/extra_galaxies/modeling.py": 13067, "scripts/interferometer/features/extra_galaxies/simulator.py": 7881, @@ -310,55 +310,73 @@ "scripts/interferometer/fit.py": 13211, "scripts/interferometer/likelihood_function.py": 20288, "scripts/interferometer/modeling.py": 28829, - "scripts/interferometer/plot.py": 6401, - "scripts/interferometer/simulator.py": 14522, + "scripts/interferometer/plot.py": 8568, + "scripts/interferometer/simulator.py": 14494, "scripts/interferometer/source_science.py": 11051, "scripts/interferometer/start_here.py": 23464, "scripts/multi_dataset/__init__.py": 0, "scripts/multi_dataset/features/__init__.py": 0, "scripts/multi_dataset/features/dataset_offsets/__init__.py": 0, - "scripts/multi_dataset/features/dataset_offsets/modeling.py": 11027, - "scripts/multi_dataset/features/dataset_offsets/simulator.py": 9222, + "scripts/multi_dataset/features/dataset_offsets/modeling.py": 11059, + "scripts/multi_dataset/features/dataset_offsets/simulator.py": 9238, "scripts/multi_dataset/features/imaging_and_interferometer/__init__.py": 0, - "scripts/multi_dataset/features/imaging_and_interferometer/modeling.py": 10259, - "scripts/multi_dataset/features/imaging_and_interferometer/simulator.py": 6659, + "scripts/multi_dataset/features/imaging_and_interferometer/modeling.py": 10291, + "scripts/multi_dataset/features/imaging_and_interferometer/simulator.py": 6683, "scripts/multi_dataset/features/imaging_and_point_source/__init__.py": 0, - "scripts/multi_dataset/features/imaging_and_point_source/modeling.py": 12218, + "scripts/multi_dataset/features/imaging_and_point_source/modeling.py": 12242, "scripts/multi_dataset/features/one_by_one/__init__.py": 0, - "scripts/multi_dataset/features/one_by_one/modeling.py": 10959, + "scripts/multi_dataset/features/one_by_one/modeling.py": 11007, "scripts/multi_dataset/features/pixelization/__init__.py": 0, - "scripts/multi_dataset/features/pixelization/modeling.py": 8423, - "scripts/multi_dataset/features/pixelization/simulator.py": 7287, + "scripts/multi_dataset/features/pixelization/modeling.py": 8447, + "scripts/multi_dataset/features/pixelization/simulator.py": 7303, "scripts/multi_dataset/features/same_wavelength/__init__.py": 0, - "scripts/multi_dataset/features/same_wavelength/modeling.py": 9785, - "scripts/multi_dataset/features/same_wavelength/simulator.py": 6388, + "scripts/multi_dataset/features/same_wavelength/modeling.py": 9809, + "scripts/multi_dataset/features/same_wavelength/simulator.py": 6390, "scripts/multi_dataset/features/slam/__init__.py": 0, - "scripts/multi_dataset/features/slam/independent.py": 25407, - "scripts/multi_dataset/features/slam/simultaneous.py": 23950, + "scripts/multi_dataset/features/slam/independent.py": 25463, + "scripts/multi_dataset/features/slam/simultaneous.py": 23982, "scripts/multi_dataset/features/wavelength_dependence/__init__.py": 0, - "scripts/multi_dataset/features/wavelength_dependence/modeling.py": 10825, - "scripts/multi_dataset/features/wavelength_dependence/simulator.py": 9037, - "scripts/multi_dataset/modeling.py": 17322, - "scripts/multi_dataset/plot.py": 7315, - "scripts/multi_dataset/simulator.py": 9479, - "scripts/multi_dataset/start_here.py": 19867, + "scripts/multi_dataset/features/wavelength_dependence/modeling.py": 10881, + "scripts/multi_dataset/features/wavelength_dependence/simulator.py": 9061, + "scripts/multi_dataset/modeling.py": 17346, + "scripts/multi_dataset/plot.py": 7347, + "scripts/multi_dataset/simulator.py": 9503, + "scripts/multi_dataset/start_here.py": 19875, "scripts/multi_galaxy/__init__.py": 0, "scripts/multi_galaxy/features/__init__.py": 0, "scripts/multi_galaxy/features/extra_galaxies/__init__.py": 0, "scripts/multi_galaxy/features/extra_galaxies/modeling.py": 15650, "scripts/multi_galaxy/features/extra_galaxies/simulator.py": 10059, + "scripts/multi_galaxy/features/extra_galaxies/slam.py": 22490, + "scripts/multi_galaxy/features/linear_light_profiles/__init__.py": 0, + "scripts/multi_galaxy/features/linear_light_profiles/fit.py": 11477, + "scripts/multi_galaxy/features/linear_light_profiles/likelihood_function.py": 13947, + "scripts/multi_galaxy/features/linear_light_profiles/modeling.py": 14887, + "scripts/multi_galaxy/features/linear_light_profiles/slam.py": 19308, + "scripts/multi_galaxy/features/multi_gaussian_expansion/__init__.py": 0, + "scripts/multi_galaxy/features/multi_gaussian_expansion/fit.py": 10772, + "scripts/multi_galaxy/features/multi_gaussian_expansion/likelihood_function.py": 12672, + "scripts/multi_galaxy/features/multi_gaussian_expansion/modeling.py": 14885, + "scripts/multi_galaxy/features/multi_gaussian_expansion/simulator.py": 10154, + "scripts/multi_galaxy/features/multi_gaussian_expansion/slam.py": 17869, + "scripts/multi_galaxy/features/multi_gaussian_expansion/source_science.py": 9833, + "scripts/multi_galaxy/features/no_lens_light/__init__.py": 0, + "scripts/multi_galaxy/features/no_lens_light/modeling.py": 15086, + "scripts/multi_galaxy/features/no_lens_light/simulator.py": 12776, + "scripts/multi_galaxy/features/no_lens_light/slam.py": 17293, "scripts/multi_galaxy/features/scaling_relation/__init__.py": 0, - "scripts/multi_galaxy/features/scaling_relation/fit.py": 10065, - "scripts/multi_galaxy/features/scaling_relation/likelihood_function.py": 9523, - "scripts/multi_galaxy/features/scaling_relation/modeling.py": 12918, - "scripts/multi_galaxy/features/scaling_relation/simulator.py": 9705, - "scripts/multi_galaxy/features/scaling_relation/slam.py": 32383, + "scripts/multi_galaxy/features/scaling_relation/fit.py": 10209, + "scripts/multi_galaxy/features/scaling_relation/likelihood_function.py": 9769, + "scripts/multi_galaxy/features/scaling_relation/modeling.py": 13509, + "scripts/multi_galaxy/features/scaling_relation/simulator.py": 9887, + "scripts/multi_galaxy/features/scaling_relation/slam.py": 35494, "scripts/multi_galaxy/fit.py": 26174, "scripts/multi_galaxy/likelihood_function.py": 28830, - "scripts/multi_galaxy/modeling.py": 46981, + "scripts/multi_galaxy/modeling.py": 47036, "scripts/multi_galaxy/plot.py": 7858, "scripts/multi_galaxy/simulator.py": 19058, "scripts/multi_galaxy/simulator_sample.py": 12601, + "scripts/multi_galaxy/slam.py": 26113, "scripts/multi_galaxy/source_science.py": 16739, "scripts/multi_galaxy/start_here.py": 30033, "scripts/point_source/__init__.py": 0, @@ -369,33 +387,33 @@ "scripts/point_source/features/extra_galaxies/__init__.py": 0, "scripts/point_source/features/extra_galaxies/modeling.py": 19537, "scripts/point_source/features/extra_galaxies/simulator.py": 14376, - "scripts/point_source/features/fluxes.py": 7398, + "scripts/point_source/features/fluxes.py": 8825, "scripts/point_source/features/multiple_sources/__init__.py": 0, - "scripts/point_source/features/multiple_sources/modeling.py": 13008, - "scripts/point_source/features/multiple_sources/simulator.py": 10262, + "scripts/point_source/features/multiple_sources/modeling.py": 13048, + "scripts/point_source/features/multiple_sources/simulator.py": 10270, "scripts/point_source/features/scaling_relation/__init__.py": 0, - "scripts/point_source/features/scaling_relation/fit.py": 8974, - "scripts/point_source/features/scaling_relation/likelihood_function.py": 8916, + "scripts/point_source/features/scaling_relation/fit.py": 9530, + "scripts/point_source/features/scaling_relation/likelihood_function.py": 8936, "scripts/point_source/features/scaling_relation/modeling.py": 13488, "scripts/point_source/features/scaling_relation/simulator.py": 11704, - "scripts/point_source/features/time_delays.py": 8936, - "scripts/point_source/fit.py": 28476, - "scripts/point_source/modeling.py": 24900, - "scripts/point_source/plot.py": 4530, + "scripts/point_source/features/time_delays.py": 10230, + "scripts/point_source/fit.py": 31250, + "scripts/point_source/modeling.py": 25219, + "scripts/point_source/plot.py": 4918, "scripts/point_source/simulator.py": 22113, "scripts/point_source/simulator_sample.py": 10427, - "scripts/point_source/start_here.py": 27156, + "scripts/point_source/start_here.py": 27172, "scripts/weak/__init__.py": 0, "scripts/weak/features/__init__.py": 0, "scripts/weak/features/strong_lensing/__init__.py": 0, "scripts/weak/features/strong_lensing/a2744.py": 11339, "scripts/weak/features/strong_lensing/fit.py": 5204, - "scripts/weak/features/strong_lensing/modeling.py": 8256, + "scripts/weak/features/strong_lensing/modeling.py": 8272, "scripts/weak/features/strong_lensing/simulator.py": 6461, "scripts/weak/fit.py": 7595, "scripts/weak/likelihood_function.py": 10693, "scripts/weak/modeling.py": 12975, - "scripts/weak/plot.py": 4741, + "scripts/weak/plot.py": 5749, "scripts/weak/real_data/__init__.py": 0, "scripts/weak/real_data/a2744.py": 11906, "scripts/weak/simulator.py": 8318, diff --git a/llms-full.txt b/llms-full.txt index dcf686689..7793f8e5d 100644 --- a/llms-full.txt +++ b/llms-full.txt @@ -410,7 +410,7 @@ AUTO-GENERATED by PyAutoHands — do not edit by hand; regenerate with generate. - [Modeling: Time Delays](scripts/point_source/features/time_delays.py): A measurable quantity of a point source is its time delay—the time it takes for light to travel from the source to the observer for each multiple image of the point source (e.g., the quasar images). This is often expressed as the relative time delay between each image and the image with the shortest time delay, which is often referred to as the "reference image." - Contents: Model, Dataset, Point Solver, Search, Analysis, Run Times, Result, Cosmology - [Guide: Point Sources](scripts/point_source/fit.py): -------------------- - - Contents: Lensed Point Source, Point Source, Point Solver, Number of Solutions, Solving the Lens Equation, Triangle Tracing, Dataset, Name Pairing, Fitting, Chi Squared, Fluxes, Flux Point Dataset, Flux Fitting, Time Delays, Time Delay Fitting, New User Wrap Up, Shape Solver + - Contents: Lensed Point Source, Point Source, Point Solver, Number of Solutions, Solving the Lens Equation, Triangle Tracing, Dataset, Name Pairing, Fitting, Chi Squared, Solved Source Centre, Fluxes, Flux Point Dataset, Flux Fitting, Time Delays, Time Delay Fitting, New User Wrap Up, Shape Solver - [Modeling: Start Here](scripts/point_source/modeling.py): This script is the starting point for lens modeling of point-source lens datasets, for example the multiple image positions of a lensed quasar. - Contents: Not Using Light Profiles, Model, Dataset, Point Solver, Model Composition, Name Pairing, Coordinates, Search, Unique Identifier, Live Visual Update, Chi Squared, Analysis, JAX, VRAM Use, Run Times, Output Folder Layout, Result, Results, Modeling Customization - [Plots: Point Source](scripts/point_source/plot.py): This example shows how to plot a `PointDataset` and a `FitPointDataset` fit. diff --git a/notebooks/cluster/lenstool/modeling.ipynb b/notebooks/cluster/lenstool/modeling.ipynb index 2cb127955..157b26b02 100644 --- a/notebooks/cluster/lenstool/modeling.ipynb +++ b/notebooks/cluster/lenstool/modeling.ipynb @@ -37,6 +37,9 @@ " ``arcs.dat`` ``point_datasets.csv`` \u2192 ``al.PointDataset`` list\n", " sigposArcsec ``positions_noise`` column (same chi-squared)\n", " source-plane optimization ``al.FitPositionsSource`` (this script)\n", + " (``al.FitPositionsSourceSolved`` solves the source centre\n", + " analytically; its ``weighting = \"magnification\"`` option\n", + " matches Lenstool's scalar \u00b5\u00b2 convention \u2014 see guides)\n", " image-plane optimization ``al.FitPositionsImagePair*`` (heavier; see guides)\n", " ``best.par`` the max-likelihood instance of the PyAutoLens fit\n", "\n", diff --git a/notebooks/cluster/likelihood_function.ipynb b/notebooks/cluster/likelihood_function.ipynb index 20f70012a..ed1bc40fa 100644 --- a/notebooks/cluster/likelihood_function.ipynb +++ b/notebooks/cluster/likelihood_function.ipynb @@ -21,6 +21,10 @@ " observed position, and measure the image-plane residuals. More intuitive (residuals in\n", " arc-seconds), but slower and carries pairing pathologies the source-plane variant avoids.\n", "\n", + "Each flavour also has a ``*Solved`` sibling (``FitPositionsSourceSolved`` and the image-plane\n", + "solved variants) which solves the source-plane centre analytically instead of sampling it \u2014 see\n", + "the ``Source-Plane Centroid`` section below and ``guides/point_source_pairing.py``.\n", + "\n", "We walk through both, end to end, with the actual library formulae. The standard cluster model is\n", "assumed: all lens-plane galaxies at ``z = 0.5``, two background sources at *different* redshifts\n", "(``z = 1.0`` and ``z = 2.0``). Multi-plane ray tracing therefore applies and we explain how the\n", @@ -431,20 +435,23 @@ "\n", "The \"reference point\" against which we measure the source-plane scatter has two options:\n", "\n", - " 1. **Truth ``Point`` centre.** Each source carries a ``Point`` profile in the model whose\n", + " 1. **Free ``Point`` centre.** Each source carries a ``Point`` profile in the model whose\n", " ``centre`` is a free parameter (or fixed to the truth here). At a model fit, ``Point.centre``\n", - " *is* the source-plane (y, x) the multiple images should converge to.\n", - "\n", - " 2. **Barycenter of back-traced positions.** Pretend you don't know the truth centre. Compute the\n", - " centroid of the back-traced positions \u2014 at the right model the centroid sits where the source\n", - " actually is, and the residuals are scatter around it.\n", - "\n", - "Option (1) is what ``al.FitPositionsSource(profile=point_profile)`` uses. Option (2) is what\n", - "``al.FitPositionsSource(profile=None)`` uses (the default during model fits, because at search\n", - "time the model doesn't yet know the truth).\n", - "\n", - "For this walkthrough we use option (1) so the residuals have a well-defined physical meaning\n", - "(distance from truth, not from a derived centroid)." + " *is* the source-plane (y, x) the multiple images should converge to, and the sampler explores\n", + " it as two non-linear parameters per source.\n", + "\n", + " 2. **Analytically solved centre.** Give each source the parameter-free ``al.ps.PointSolved``\n", + " profile and fit with ``al.FitPositionsSourceSolved``: the centre is computed in closed form as\n", + " the precision-weighted mean of the back-traced positions (Lombardi 2024, arXiv:2406.15280,\n", + " \u00a75.1), with a tensor weighting derived from the lensing Jacobian, and the likelihood is\n", + " analytically marginalized over it. This removes 2 free parameters per source from the search \u2014\n", + " at cluster scale, with many sources, a substantial dimensionality reduction. See\n", + " ``guides/point_source_pairing.py`` for the full solved-variant matrix.\n", + "\n", + "Option (1) is what ``al.FitPositionsSource(profile=point_profile)`` uses and is what this\n", + "walkthrough demonstrates, so the residuals have a well-defined physical meaning (distance from\n", + "truth). Option (2) is its solved sibling ``FitPositionsSourceSolved``, whose reference point is a\n", + "weighted barycenter of the back-traced positions rather than a sampled parameter." ] }, { @@ -1096,13 +1103,18 @@ "| **Best for** | Fast Nautilus fits; JAX-jit'd parameter estimation | Final residual visualisation; cases where pairing is unambiguous |\n", "\n", "For most cluster fits the source-plane chi\u00b2 is the right default: it's faster, JAX-compatible,\n", - "and doesn't suffer pairing pathologies. The image-plane chi\u00b2 is most useful for diagnostic\n", - "visualisation of where the model's predicted images sit relative to the observed ones, and for\n", - "cases where the source-plane chi\u00b2 has bias issues that need cross-checking.\n", - "\n", - "The cluster modelling script at ``scripts/cluster/modeling.py`` uses ``AnalysisPoint`` which\n", - "selects the chi\u00b2 flavour via its constructor; consult the ``AnalysisPoint`` docstring for the\n", - "current default.\n", + "and doesn't suffer pairing pathologies. Its solved sibling ``FitPositionsSourceSolved`` (paired\n", + "with ``al.ps.PointSolved`` sources) sharpens this further \u2014 the source centres drop out of the\n", + "non-linear space entirely (2 parameters per source) at timing-noise-level extra cost per call,\n", + "and its tensor weighting is a better error model than the scalar ``\u03bc\u00b2`` used above. The\n", + "recommended cluster workflow is therefore: **search with ``FitPositionsSourceSolved``, validate\n", + "with the image-plane chi\u00b2** on the max-likelihood model. The image-plane chi\u00b2 remains the\n", + "diagnostic tool \u2014 visualising where predicted images sit relative to observed ones, and\n", + "cross-checking any source-plane bias.\n", + "\n", + "The cluster modelling script at ``scripts/cluster/modeling.py`` uses ``AnalysisPoint``, which\n", + "selects the chi\u00b2 flavour via its ``fit_positions_cls`` input (default\n", + "``FitPositionsImagePairRepeat``).\n", "\n", "__Wrap Up__\n", "\n", diff --git a/notebooks/cluster/modeling.ipynb b/notebooks/cluster/modeling.ipynb index 7f5771e64..2121369cd 100644 --- a/notebooks/cluster/modeling.ipynb +++ b/notebooks/cluster/modeling.ipynb @@ -669,7 +669,12 @@ "no matching dataset, **PyAutoLens** raises an error.\n", "\n", "In multi-source cluster lenses, this name pairing is what ensures every source's positions are fitted by\n", - "the correct model component." + "the correct model component.\n", + "\n", + "A centre-free alternative exists for the sources: the parameter-free ``al.ps.PointSolved`` paired with\n", + "``fit_positions_cls=al.FitPositionsSourceSolved`` solves each source centre analytically, removing 2 free\n", + "parameters per source \u2014 the recommended search-stage configuration at cluster scale (validate image-plane;\n", + "see ``guides/point_source_pairing.py``)." ] }, { diff --git a/notebooks/cluster/start_here.ipynb b/notebooks/cluster/start_here.ipynb index add6689cd..e19d5b84f 100644 --- a/notebooks/cluster/start_here.ipynb +++ b/notebooks/cluster/start_here.ipynb @@ -381,7 +381,11 @@ "\n", "The model is composed below in four blocks: main-tier loop, host halo, source-tier loop, scaling-tier\n", "loop (defining the shared ``sigma_ref`` normalization once outside the loop). The four blocks are then\n", - "bundled into a single ``af.Collection`` model that the analysis will receive." + "bundled into a single ``af.Collection`` model that the analysis will receive.\n", + "\n", + "Each source's ``Point`` centre is sampled as 2 free parameters. A centre-free alternative\n", + "(``al.ps.PointSolved`` + the ``*Solved`` fit classes) solves the centres analytically instead \u2014 see\n", + "``guides/point_source_pairing.py``." ] }, { diff --git a/notebooks/guides/point_source_pairing.ipynb b/notebooks/guides/point_source_pairing.ipynb index d55cbde51..1fcece0f9 100644 --- a/notebooks/guides/point_source_pairing.ipynb +++ b/notebooks/guides/point_source_pairing.ipynb @@ -15,10 +15,11 @@ "rather than explicit choices.\n", "\n", "This guide documents PyAutoLens's choices: the three image-plane pairing schemes, the\n", - "over/under-prediction policies, the solver settings that interact with them at cluster scale, and\n", - "the source-plane vs image-plane chi-squared trade-off. It is the reference the cluster examples\n", - "(``scripts/cluster/``, including the Lenstool walkthrough in ``scripts/cluster/lenstool/``) point\n", - "at for likelihood choices.\n", + "over/under-prediction policies, the solver settings that interact with them at cluster scale, the\n", + "source-plane vs image-plane chi-squared trade-off, and the second axis of the likelihood-option\n", + "matrix \u2014 whether the source-plane centre is a free model parameter or analytically solved\n", + "(``al.ps.PointSolved``). It is the reference the cluster examples (``scripts/cluster/``, including\n", + "the Lenstool walkthrough in ``scripts/cluster/lenstool/``) point at for likelihood choices.\n", "\n", "__The two failure modes__\n", "\n", @@ -107,7 +108,96 @@ "The pragmatic workflow at cluster scale: **search with the source-plane chi-squared, validate with\n", "the image-plane chi-squared** \u2014 run the image-plane fit (and inspect ``n_unmatched_model_positions``\n", "plus the per-system image counts) on the max-likelihood model before publishing, exactly as the\n", - "Lenstool-users example does.\n", + "Lenstool-users example does. The solved-centre variants below sharpen this advice further.\n", + "\n", + "__The second axis: free vs analytically-solved source centre__\n", + "\n", + "Every fit above anchors its prediction to a source-plane centre \u2014 the ``centre`` of the model's\n", + "``al.ps.Point`` (or ``al.ps.PointFlux``), sampled as two non-linear parameters per point source.\n", + "That centre can instead be **solved analytically** from the observed positions and the current mass\n", + "model, dropping out of the non-linear parameter space entirely. The model component for this is\n", + "``al.ps.PointSolved``, which is parameter-free, and each fit class gains a ``*Solved`` sibling:\n", + "\n", + "| Free-centre fit | Solved-centre sibling |\n", + "|---|---|\n", + "| ``FitPositionsSource`` | ``FitPositionsSourceSolved`` |\n", + "| ``FitPositionsImagePairRepeat`` | ``FitPositionsImagePairRepeatSolved`` |\n", + "| ``FitPositionsImagePairAll`` | ``FitPositionsImagePairAllSolved`` |\n", + "| ``FitPositionsImagePair`` | \u2014 (deliberately none; see below) |\n", + "| ``FitFluxes`` | ``FitFluxesSolved`` |\n", + "| ``FitTimeDelays`` | ``FitTimeDelaysSolved`` |\n", + "\n", + "The Hungarian ``FitPositionsImagePair`` has no solved sibling: its linear-sum-assignment pairing is\n", + "numpy-only (not JAX-differentiable) and the repeat/all-pairs solved variants supersede it.\n", + "\n", + "The payoff is dimensionality: each point source loses 2 free parameters (and ``FitFluxesSolved``\n", + "drops the ``flux`` parameter too, so ``PointSolved`` alone covers positions + fluxes + time-delay\n", + "datasets). A 10-source cluster fit loses 20 parameters \u2014 at cluster scale, where sources are many\n", + "and each adds little individual constraint on the mass model, this is where the gain is largest.\n", + "\n", + "Mixing the two conventions is an error in both directions, raised loudly as\n", + "``PointProfileMismatchException``: a ``*Solved`` fit given a centre-bearing profile would sample two\n", + "parameters the analytic solve silently ignores, and a free-centre fit given ``PointSolved`` has no\n", + "centre to read.\n", + "\n", + "__The solved source-plane chi-squared (FitPositionsSourceSolved)__\n", + "\n", + "The solved source-plane fit follows Lombardi (2024, arXiv:2406.15280, \u00a75.1): Taylor-expanding the\n", + "lens equation around each observed image position makes the back-traced source position linear in\n", + "the source centre, so the optimal centre has a closed form. Each back-traced position ``\u03b2\u0302\u1d62`` is\n", + "weighted by its precision tensor ``W\u1d62 = A\u1d62\u207b\u1d40 \u0398\u1d62 A\u1d62\u207b\u00b9`` (``A = \u2202\u03b2/\u2202\u03b8`` is the lensing Jacobian,\n", + "``\u0398\u1d62 = \u03c3\u1d62\u207b\u00b2 I`` the image-plane precision), and the solved centre is the precision-weighted mean\n", + "\n", + " \u03b2* = (\u03a3\u1d62 W\u1d62)\u207b\u00b9 \u03a3\u1d62 W\u1d62 \u03b2\u0302\u1d62\n", + "\n", + "with the likelihood analytically marginalized over the centre (a flat prior; the fit's\n", + "``marginalization_term`` property is the resulting log-determinant contribution).\n", + "\n", + "The tensor weighting is also a better error model than the scalar ``\u00b5\u00b2/\u03c3\u00b2`` weighting of\n", + "``FitPositionsSource``. The eigenvalues of ``W\u1d62`` are ``\u03bb\u00b2/\u03c3\u1d62\u00b2`` with ``\u03bb`` the linear stretch\n", + "along each eigendirection: isotropically each stretch is ``\u221a\u00b5``, while near a critical curve the\n", + "tangential stretch is ``\u2248 \u00b5`` \u2014 the scalar ``\u00b5\u00b2`` convention coincides with the tensor only in that\n", + "near-critical tangential limit, so it over-weights images everywhere else. The\n", + "``weighting`` class attribute selects the convention: ``\"jacobian\"`` (default, the tensor) or\n", + "``\"magnification\"`` (the scalar, retained for comparisons with the traditional Lenstool-style\n", + "convention).\n", + "\n", + "__Solved image-plane variants__\n", + "\n", + "``FitPositionsImagePairRepeatSolved`` and ``FitPositionsImagePairAllSolved`` reuse the same solved\n", + "``\u03b2*`` to drive the forward lens-equation solve, then apply their scheme's image-plane pairing\n", + "chi-squared unchanged. These are **not** from Lombardi (2024) \u2014 the paper keeps the centre free in\n", + "its image-plane likelihoods. They are a PyAutoLens extension in the spirit of glafic's\n", + "source-position optimization (Oguri 2010, PASJ 62, 1017), which likewise eliminates the source\n", + "position from the sampled space of an image-plane chi-squared.\n", + "\n", + "__Solved fluxes and time delays__\n", + "\n", + "``FitFluxesSolved`` solves the source flux the same way (following Lombardi 2024 \u00a76.1, ported to\n", + "flux space to match PyAutoLens's flux-space Gaussian noise maps): ``F* = \u03a3\u1d62 \u00b5\u1d62 f\u0302\u1d62/\u03c3\u1d62\u00b2 / \u03a3\u1d62 \u00b5\u1d62\u00b2/\u03c3\u1d62\u00b2``\n", + "\u2014 magnification-first, mirroring ``FitFluxes.model_data = |\u00b5\u1d62|\u00b7F`` \u2014 plus its marginalization term.\n", + "``FitTimeDelaysSolved`` replaces the reference-image min-subtraction of ``FitTimeDelays`` with a\n", + "precision-weighted analytic reference time ``T*``. Both are selected via the ``fit_flux_cls`` /\n", + "``fit_time_delays_cls`` inputs of ``FitPointDataset`` / ``AnalysisPoint``, mirroring\n", + "``fit_positions_cls``.\n", + "\n", + "__Missing-image penalty__\n", + "\n", + "``FitPositionsImagePairAll``'s mixture normalization already implements the principled\n", + "over-prediction Occam factor (the ``1/P^I`` term of Lombardi 2024's mixture likelihood).\n", + "``FitPositionsImagePairRepeat``'s ``unmatched_model_policy`` heuristics stay as they are:\n", + "best-match pairing is not a normalized mixture, so no principled combinatorial term applies to it.\n", + "\n", + "__Choosing at cluster scale__\n", + "\n", + "Profiling on the standard cluster model (see the likelihood-breakdown scripts referenced above)\n", + "puts the analytic ``\u03b2*`` solve at timing-noise-level overhead per likelihood call \u2014 net +3% on the\n", + "image-plane likelihood, +9\u201330% on sub-0.1 s eager source-plane totals \u2014 while removing 2 free\n", + "parameters per point source from the non-linear space. The recommendation for cluster fits is\n", + "therefore to sharpen the workflow above: **search with ``FitPositionsSourceSolved`` (with\n", + "``al.ps.PointSolved`` sources), validate with the image-plane chi-squared** on the max-likelihood\n", + "model. Reserve free-centre ``FitPositionsSource`` for direct comparisons with codes that sample the\n", + "source position (its ``weighting = \"magnification\"`` scalar convention matches Lenstool's).\n", "\n", "__Demonstration__\n", "\n", @@ -273,6 +363,37 @@ "outputs": [], "execution_count": null }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Case 4 \u2014 solved source centre__\n", + "\n", + "The solved variants swap the source's ``al.ps.Point`` (free centre) for the parameter-free\n", + "``al.ps.PointSolved``. The source-plane fit needs no solver at all: it back-traces the observed\n", + "positions and solves the centre analytically. The fit exposes the solved centre via\n", + "``source_plane_coordinate`` and the analytic-marginalization contribution via\n", + "``marginalization_term``." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "source_solved = al.Galaxy(redshift=1.0, point_0=al.ps.PointSolved())\n", + "tracer_solved = al.Tracer(galaxies=[lens, source_solved])\n", + "\n", + "fit = al.FitPositionsSourceSolved(\n", + " name=\"point_0\", data=data, noise_map=noise_map, tracer=tracer_solved, solver=None\n", + ")\n", + "print(\"Case 4 \u2014 solved source centre (source-plane fit, no free centre parameters):\")\n", + "print(f\" solved centre beta* = {tuple(round(float(c), 4) for c in fit.source_plane_coordinate)}\")\n", + "print(f\" chi_squared = {float(fit.chi_squared):.4f}\")\n", + "print(f\" marginalization_term = {float(fit.marginalization_term):.4f}\")" + ], + "outputs": [], + "execution_count": null + }, { "cell_type": "markdown", "metadata": {}, diff --git a/notebooks/point_source/features/fluxes.ipynb b/notebooks/point_source/features/fluxes.ipynb index 589acda62..4d1dfdac0 100644 --- a/notebooks/point_source/features/fluxes.ipynb +++ b/notebooks/point_source/features/fluxes.ipynb @@ -356,6 +356,49 @@ "outputs": [], "execution_count": null }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Analytic Flux (Solved)__\n", + "\n", + "An alternative to sampling `flux` is to solve it analytically: `al.FitFluxesSolved` computes the\n", + "magnification-weighted best-fit source flux from the data in closed form, so the model needs no\n", + "`flux` parameter at all. The model component is then the parameter-free `al.ps.PointSolved` (which\n", + "also solves the source centre analytically \u2014 solved fits require it, and mixing solved fits with\n", + "`Point`/`PointFlux` raises an error), reducing the model above from N=8 to N=5.\n", + "\n", + "The composition is shown below (not fitted here \u2014 the model-fit in this script demonstrates the\n", + "free-`flux` convention; see `guides/point_source_pairing.py` for when to prefer each):\n", + "\n", + "The microlensing caveat at the top of this script applies equally to the solved flux: solving the\n", + "source flux analytically does not make image fluxes any less affected by microlensing." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "lens_solved = af.Model(al.Galaxy, redshift=0.5, mass=al.mp.Isothermal)\n", + "\n", + "source_solved = af.Model(al.Galaxy, redshift=1.0, point_0=af.Model(al.ps.PointSolved))\n", + "\n", + "model_solved = af.Collection(\n", + " galaxies=af.Collection(lens=lens_solved, source=source_solved)\n", + ")\n", + "\n", + "analysis_solved = al.AnalysisPoint(\n", + " dataset=dataset,\n", + " solver=solver,\n", + " fit_positions_cls=al.FitPositionsSourceSolved, # Solved fits pair with `PointSolved`.\n", + " fit_flux_cls=al.FitFluxesSolved, # Flux solved analytically, no free `flux` parameter.\n", + ")\n", + "\n", + "print(model_solved.info)" + ], + "outputs": [], + "execution_count": null + }, { "cell_type": "markdown", "metadata": {}, diff --git a/notebooks/point_source/features/time_delays.ipynb b/notebooks/point_source/features/time_delays.ipynb index 8d4c0b0ab..680d47af6 100644 --- a/notebooks/point_source/features/time_delays.ipynb +++ b/notebooks/point_source/features/time_delays.ipynb @@ -356,6 +356,47 @@ "outputs": [], "execution_count": null }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Analytic Reference Time (Solved)__\n", + "\n", + "The default `FitTimeDelays` compares *relative* delays by subtracting the minimum delay from both\n", + "the data and the model (the \"reference image\" convention). An alternative is\n", + "`al.FitTimeDelaysSolved`, which instead solves a precision-weighted analytic reference time from\n", + "the data in closed form and analytically marginalizes over it \u2014 a smooth (and JAX-differentiable)\n", + "alternative to the min-subtraction.\n", + "\n", + "Solved fit classes pair with the parameter-free `al.ps.PointSolved` model component (which also\n", + "solves the source centre, dropping its 2 free parameters; mixing solved fits with `Point` raises an\n", + "error). The composition below is shown for reference and not fitted in this script:" + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "lens_solved = af.Model(al.Galaxy, redshift=0.5, mass=al.mp.Isothermal)\n", + "\n", + "source_solved = af.Model(al.Galaxy, redshift=1.0, point_0=af.Model(al.ps.PointSolved))\n", + "\n", + "model_solved = af.Collection(\n", + " galaxies=af.Collection(lens=lens_solved, source=source_solved)\n", + ")\n", + "\n", + "analysis_solved = al.AnalysisPoint(\n", + " dataset=dataset,\n", + " solver=solver,\n", + " fit_positions_cls=al.FitPositionsSourceSolved, # Solved fits pair with `PointSolved`.\n", + " fit_time_delays_cls=al.FitTimeDelaysSolved, # Analytic reference time instead of min-subtraction.\n", + ")\n", + "\n", + "print(model_solved.info)" + ], + "outputs": [], + "execution_count": null + }, { "cell_type": "markdown", "metadata": {}, diff --git a/notebooks/point_source/fit.ipynb b/notebooks/point_source/fit.ipynb index a7a99f595..4dd42549b 100644 --- a/notebooks/point_source/fit.ipynb +++ b/notebooks/point_source/fit.ipynb @@ -37,6 +37,7 @@ "- **Name Pairing:** The names of the point-source datasets have an even more important role, the names are used to pair.\n", "- **Fitting:** Fit the lens model to the dataset and inspect the results.\n", "- **Chi Squared:** For point-source modeling, there are many different ways to define the likelihood function, broadly.\n", + "- **Solved Source Centre:** The `*Solved` fit variants solve the source-plane centre analytically instead of sampling it.\n", "- **Fluxes:** Another measurable quantity of a point source is its flux\u2014the total amount of light received from.\n", "- **Flux Point Dataset:** The fluxes are not input a `PointDataset` object, alongside the image-plane coordinates of the.\n", "- **Flux Fitting:** Above, we used a `FitPointDataset` to fit the positions of the point source in the image-plane.\n", @@ -602,6 +603,61 @@ "cell_type": "markdown", "metadata": {}, "source": [ + "A note on defaults: `FitPointDataset` itself defaults to `FitPositionsImagePair` (Hungarian pairing,\n", + "no repeats), whereas `AnalysisPoint` \u2014 the object used in the model-fitting examples \u2014 defaults to\n", + "`FitPositionsImagePairRepeat`. The examples above pass `fit_positions_cls` explicitly so there is no\n", + "ambiguity about which chi-squared is being used.\n", + "\n", + "__Solved Source Centre__\n", + "\n", + "Every fit above reads the source-plane centre from the model's `Point` profile, whose `centre` is a\n", + "free parameter during model-fitting. Each fit class also has a `*Solved` variant which instead solves\n", + "the centre analytically from the observed positions and the mass model, removing those 2 free\n", + "parameters per point source from the non-linear search: `FitPositionsSourceSolved` (a\n", + "precision-weighted source-plane centre, following Lombardi 2024, arXiv:2406.15280) and\n", + "`FitPositionsImagePairRepeatSolved` / `FitPositionsImagePairAllSolved` (the same solved centre\n", + "driving the image-plane forward solve).\n", + "\n", + "The solved variants pair with the parameter-free `al.ps.PointSolved` model component instead\n", + "of `Point` \u2014 mixing the two conventions (a `Solved` fit with a `Point`, or vice versa) raises\n", + "a `PointProfileMismatchException`, so a model cannot silently sample centre parameters the fit\n", + "ignores." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "point_source_solved = al.ps.PointSolved()\n", + "\n", + "source_galaxy_solved = al.Galaxy(redshift=1.0, point_0=point_source_solved)\n", + "\n", + "tracer_solved = al.Tracer(galaxies=[lens_galaxy, source_galaxy_solved])\n", + "\n", + "fit = al.FitPointDataset(\n", + " dataset=dataset,\n", + " tracer=tracer_solved,\n", + " solver=solver,\n", + " fit_positions_cls=al.FitPositionsSourceSolved, # Solved-centre source-plane chi-squared\n", + ")\n", + "\n", + "print(\"Analytically Solved Source-Plane Centre:\")\n", + "print(fit.positions.source_plane_coordinate)\n", + "\n", + "print(\"Log Likelihood with Solved Centre:\")\n", + "print(fit.positions.log_likelihood)" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The full matrix of pairing schemes \u00d7 centre treatments \u2014 including the solved flux and time-delay\n", + "variants used below, the tensor-vs-scalar magnification weighting, and when each option is\n", + "preferable \u2014 is documented in `guides/point_source_pairing.py`.\n", + "\n", "__Fluxes__\n", "\n", "Another measurable quantity of a point source is its flux\u2014the total amount of light received from each multiple image \n", @@ -765,6 +821,10 @@ "cell_type": "markdown", "metadata": {}, "source": [ + "An analytic-flux variant also exists: `fit_flux_cls=al.FitFluxesSolved` (paired with\n", + "`al.ps.PointSolved`) solves the source flux from the data and magnifications instead of sampling\n", + "a `flux` parameter \u2014 see `point_source/features/fluxes.py` and `guides/point_source_pairing.py`.\n", + "\n", "__Time Delays__\n", "\n", "Another measurable quantity of a point source is its time delay\u2014the time it takes for light to travel from the\n", @@ -905,6 +965,10 @@ "cell_type": "markdown", "metadata": {}, "source": [ + "As with fluxes, an analytic variant exists: `fit_time_delays_cls=al.FitTimeDelaysSolved` replaces\n", + "the reference-image subtraction with a precision-weighted analytic reference time \u2014 see\n", + "`point_source/features/time_delays.py` and `guides/point_source_pairing.py`.\n", + "\n", "__New User Wrap Up__\n", "\n", "The `point_source` package of the `autolens_workspace` contains numerous example scripts for performing point source\n", diff --git a/notebooks/point_source/modeling.ipynb b/notebooks/point_source/modeling.ipynb index c5ede80fe..8e3c7d926 100644 --- a/notebooks/point_source/modeling.ipynb +++ b/notebooks/point_source/modeling.ipynb @@ -324,10 +324,14 @@ "The API is fairly self explanatory and is straight forward to extend, for example adding more light profiles\n", "to the lens and source or using a different mass profile.\n", "\n", - "A full description of model composition is provided by the model cookbook: \n", + "A full description of model composition is provided by the model cookbook:\n", "\n", "https://pyautolens.readthedocs.io/en/latest/general/model_cookbook.html\n", "\n", + "A centre-free alternative exists: `al.ps.PointSolved` has no free parameters \u2014 the source centre is\n", + "solved analytically by the paired `*Solved` fit classes (e.g. `fit_positions_cls=al.FitPositionsSourceSolved`),\n", + "removing 2 parameters per point source. See `guides/point_source_pairing.py` for the full option matrix.\n", + "\n", "__Name Pairing__\n", "\n", "Every point-source dataset in the `PointDataset` has a name, which in this example was `point_0`. This `name` pairs \n", diff --git a/scripts/cluster/lenstool/modeling.py b/scripts/cluster/lenstool/modeling.py index a666680a8..3761ab6b6 100644 --- a/scripts/cluster/lenstool/modeling.py +++ b/scripts/cluster/lenstool/modeling.py @@ -32,6 +32,9 @@ ``arcs.dat`` ``point_datasets.csv`` → ``al.PointDataset`` list sigposArcsec ``positions_noise`` column (same chi-squared) source-plane optimization ``al.FitPositionsSource`` (this script) + (``al.FitPositionsSourceSolved`` solves the source centre + analytically; its ``weighting = "magnification"`` option + matches Lenstool's scalar µ² convention — see guides) image-plane optimization ``al.FitPositionsImagePair*`` (heavier; see guides) ``best.par`` the max-likelihood instance of the PyAutoLens fit diff --git a/scripts/cluster/likelihood_function.py b/scripts/cluster/likelihood_function.py index 1466e1581..46552ff46 100644 --- a/scripts/cluster/likelihood_function.py +++ b/scripts/cluster/likelihood_function.py @@ -16,6 +16,10 @@ observed position, and measure the image-plane residuals. More intuitive (residuals in arc-seconds), but slower and carries pairing pathologies the source-plane variant avoids. +Each flavour also has a ``*Solved`` sibling (``FitPositionsSourceSolved`` and the image-plane +solved variants) which solves the source-plane centre analytically instead of sampling it — see +the ``Source-Plane Centroid`` section below and ``guides/point_source_pairing.py``. + We walk through both, end to end, with the actual library formulae. The standard cluster model is assumed: all lens-plane galaxies at ``z = 0.5``, two background sources at *different* redshifts (``z = 1.0`` and ``z = 2.0``). Multi-plane ray tracing therefore applies and we explain how the @@ -314,20 +318,23 @@ The "reference point" against which we measure the source-plane scatter has two options: - 1. **Truth ``Point`` centre.** Each source carries a ``Point`` profile in the model whose + 1. **Free ``Point`` centre.** Each source carries a ``Point`` profile in the model whose ``centre`` is a free parameter (or fixed to the truth here). At a model fit, ``Point.centre`` - *is* the source-plane (y, x) the multiple images should converge to. - - 2. **Barycenter of back-traced positions.** Pretend you don't know the truth centre. Compute the - centroid of the back-traced positions — at the right model the centroid sits where the source - actually is, and the residuals are scatter around it. - -Option (1) is what ``al.FitPositionsSource(profile=point_profile)`` uses. Option (2) is what -``al.FitPositionsSource(profile=None)`` uses (the default during model fits, because at search -time the model doesn't yet know the truth). - -For this walkthrough we use option (1) so the residuals have a well-defined physical meaning -(distance from truth, not from a derived centroid). + *is* the source-plane (y, x) the multiple images should converge to, and the sampler explores + it as two non-linear parameters per source. + + 2. **Analytically solved centre.** Give each source the parameter-free ``al.ps.PointSolved`` + profile and fit with ``al.FitPositionsSourceSolved``: the centre is computed in closed form as + the precision-weighted mean of the back-traced positions (Lombardi 2024, arXiv:2406.15280, + §5.1), with a tensor weighting derived from the lensing Jacobian, and the likelihood is + analytically marginalized over it. This removes 2 free parameters per source from the search — + at cluster scale, with many sources, a substantial dimensionality reduction. See + ``guides/point_source_pairing.py`` for the full solved-variant matrix. + +Option (1) is what ``al.FitPositionsSource(profile=point_profile)`` uses and is what this +walkthrough demonstrates, so the residuals have a well-defined physical meaning (distance from +truth). Option (2) is its solved sibling ``FitPositionsSourceSolved``, whose reference point is a +weighted barycenter of the back-traced positions rather than a sampled parameter. """ source_plane_centroids = [] for i, dataset in enumerate(dataset_list): @@ -799,13 +806,18 @@ def _pair_hungarian(model_positions, observed_positions): | **Best for** | Fast Nautilus fits; JAX-jit'd parameter estimation | Final residual visualisation; cases where pairing is unambiguous | For most cluster fits the source-plane chi² is the right default: it's faster, JAX-compatible, -and doesn't suffer pairing pathologies. The image-plane chi² is most useful for diagnostic -visualisation of where the model's predicted images sit relative to the observed ones, and for -cases where the source-plane chi² has bias issues that need cross-checking. - -The cluster modelling script at ``scripts/cluster/modeling.py`` uses ``AnalysisPoint`` which -selects the chi² flavour via its constructor; consult the ``AnalysisPoint`` docstring for the -current default. +and doesn't suffer pairing pathologies. Its solved sibling ``FitPositionsSourceSolved`` (paired +with ``al.ps.PointSolved`` sources) sharpens this further — the source centres drop out of the +non-linear space entirely (2 parameters per source) at timing-noise-level extra cost per call, +and its tensor weighting is a better error model than the scalar ``μ²`` used above. The +recommended cluster workflow is therefore: **search with ``FitPositionsSourceSolved``, validate +with the image-plane chi²** on the max-likelihood model. The image-plane chi² remains the +diagnostic tool — visualising where predicted images sit relative to observed ones, and +cross-checking any source-plane bias. + +The cluster modelling script at ``scripts/cluster/modeling.py`` uses ``AnalysisPoint``, which +selects the chi² flavour via its ``fit_positions_cls`` input (default +``FitPositionsImagePairRepeat``). __Wrap Up__ diff --git a/scripts/cluster/modeling.py b/scripts/cluster/modeling.py index ae0ebdda3..b3826c32b 100644 --- a/scripts/cluster/modeling.py +++ b/scripts/cluster/modeling.py @@ -462,6 +462,11 @@ In multi-source cluster lenses, this name pairing is what ensures every source's positions are fitted by the correct model component. + +A centre-free alternative exists for the sources: the parameter-free ``al.ps.PointSolved`` paired with +``fit_positions_cls=al.FitPositionsSourceSolved`` solves each source centre analytically, removing 2 free +parameters per source — the recommended search-stage configuration at cluster scale (validate image-plane; +see ``guides/point_source_pairing.py``). """ print(model) diff --git a/scripts/cluster/start_here.py b/scripts/cluster/start_here.py index 5c1d92007..8cdf9ea2c 100644 --- a/scripts/cluster/start_here.py +++ b/scripts/cluster/start_here.py @@ -300,6 +300,10 @@ The model is composed below in four blocks: main-tier loop, host halo, source-tier loop, scaling-tier loop (defining the shared ``sigma_ref`` normalization once outside the loop). The four blocks are then bundled into a single ``af.Collection`` model that the analysis will receive. + +Each source's ``Point`` centre is sampled as 2 free parameters. A centre-free alternative +(``al.ps.PointSolved`` + the ``*Solved`` fit classes) solves the centres analytically instead — see +``guides/point_source_pairing.py``. """ redshift_lens = 0.308 source_redshifts = [dataset.redshift for dataset in dataset_list] diff --git a/scripts/guides/point_source_pairing.py b/scripts/guides/point_source_pairing.py index 943f0c3cd..d2e3e2d96 100644 --- a/scripts/guides/point_source_pairing.py +++ b/scripts/guides/point_source_pairing.py @@ -10,10 +10,11 @@ rather than explicit choices. This guide documents PyAutoLens's choices: the three image-plane pairing schemes, the -over/under-prediction policies, the solver settings that interact with them at cluster scale, and -the source-plane vs image-plane chi-squared trade-off. It is the reference the cluster examples -(``scripts/cluster/``, including the Lenstool walkthrough in ``scripts/cluster/lenstool/``) point -at for likelihood choices. +over/under-prediction policies, the solver settings that interact with them at cluster scale, the +source-plane vs image-plane chi-squared trade-off, and the second axis of the likelihood-option +matrix — whether the source-plane centre is a free model parameter or analytically solved +(``al.ps.PointSolved``). It is the reference the cluster examples (``scripts/cluster/``, including +the Lenstool walkthrough in ``scripts/cluster/lenstool/``) point at for likelihood choices. __The two failure modes__ @@ -102,7 +103,96 @@ class FitStrict(al.FitPositionsImagePairRepeat): The pragmatic workflow at cluster scale: **search with the source-plane chi-squared, validate with the image-plane chi-squared** — run the image-plane fit (and inspect ``n_unmatched_model_positions`` plus the per-system image counts) on the max-likelihood model before publishing, exactly as the -Lenstool-users example does. +Lenstool-users example does. The solved-centre variants below sharpen this advice further. + +__The second axis: free vs analytically-solved source centre__ + +Every fit above anchors its prediction to a source-plane centre — the ``centre`` of the model's +``al.ps.Point`` (or ``al.ps.PointFlux``), sampled as two non-linear parameters per point source. +That centre can instead be **solved analytically** from the observed positions and the current mass +model, dropping out of the non-linear parameter space entirely. The model component for this is +``al.ps.PointSolved``, which is parameter-free, and each fit class gains a ``*Solved`` sibling: + +| Free-centre fit | Solved-centre sibling | +|---|---| +| ``FitPositionsSource`` | ``FitPositionsSourceSolved`` | +| ``FitPositionsImagePairRepeat`` | ``FitPositionsImagePairRepeatSolved`` | +| ``FitPositionsImagePairAll`` | ``FitPositionsImagePairAllSolved`` | +| ``FitPositionsImagePair`` | — (deliberately none; see below) | +| ``FitFluxes`` | ``FitFluxesSolved`` | +| ``FitTimeDelays`` | ``FitTimeDelaysSolved`` | + +The Hungarian ``FitPositionsImagePair`` has no solved sibling: its linear-sum-assignment pairing is +numpy-only (not JAX-differentiable) and the repeat/all-pairs solved variants supersede it. + +The payoff is dimensionality: each point source loses 2 free parameters (and ``FitFluxesSolved`` +drops the ``flux`` parameter too, so ``PointSolved`` alone covers positions + fluxes + time-delay +datasets). A 10-source cluster fit loses 20 parameters — at cluster scale, where sources are many +and each adds little individual constraint on the mass model, this is where the gain is largest. + +Mixing the two conventions is an error in both directions, raised loudly as +``PointProfileMismatchException``: a ``*Solved`` fit given a centre-bearing profile would sample two +parameters the analytic solve silently ignores, and a free-centre fit given ``PointSolved`` has no +centre to read. + +__The solved source-plane chi-squared (FitPositionsSourceSolved)__ + +The solved source-plane fit follows Lombardi (2024, arXiv:2406.15280, §5.1): Taylor-expanding the +lens equation around each observed image position makes the back-traced source position linear in +the source centre, so the optimal centre has a closed form. Each back-traced position ``β̂ᵢ`` is +weighted by its precision tensor ``Wᵢ = Aᵢ⁻ᵀ Θᵢ Aᵢ⁻¹`` (``A = ∂β/∂θ`` is the lensing Jacobian, +``Θᵢ = σᵢ⁻² I`` the image-plane precision), and the solved centre is the precision-weighted mean + + β* = (Σᵢ Wᵢ)⁻¹ Σᵢ Wᵢ β̂ᵢ + +with the likelihood analytically marginalized over the centre (a flat prior; the fit's +``marginalization_term`` property is the resulting log-determinant contribution). + +The tensor weighting is also a better error model than the scalar ``µ²/σ²`` weighting of +``FitPositionsSource``. The eigenvalues of ``Wᵢ`` are ``λ²/σᵢ²`` with ``λ`` the linear stretch +along each eigendirection: isotropically each stretch is ``√µ``, while near a critical curve the +tangential stretch is ``≈ µ`` — the scalar ``µ²`` convention coincides with the tensor only in that +near-critical tangential limit, so it over-weights images everywhere else. The +``weighting`` class attribute selects the convention: ``"jacobian"`` (default, the tensor) or +``"magnification"`` (the scalar, retained for comparisons with the traditional Lenstool-style +convention). + +__Solved image-plane variants__ + +``FitPositionsImagePairRepeatSolved`` and ``FitPositionsImagePairAllSolved`` reuse the same solved +``β*`` to drive the forward lens-equation solve, then apply their scheme's image-plane pairing +chi-squared unchanged. These are **not** from Lombardi (2024) — the paper keeps the centre free in +its image-plane likelihoods. They are a PyAutoLens extension in the spirit of glafic's +source-position optimization (Oguri 2010, PASJ 62, 1017), which likewise eliminates the source +position from the sampled space of an image-plane chi-squared. + +__Solved fluxes and time delays__ + +``FitFluxesSolved`` solves the source flux the same way (following Lombardi 2024 §6.1, ported to +flux space to match PyAutoLens's flux-space Gaussian noise maps): ``F* = Σᵢ µᵢ f̂ᵢ/σᵢ² / Σᵢ µᵢ²/σᵢ²`` +— magnification-first, mirroring ``FitFluxes.model_data = |µᵢ|·F`` — plus its marginalization term. +``FitTimeDelaysSolved`` replaces the reference-image min-subtraction of ``FitTimeDelays`` with a +precision-weighted analytic reference time ``T*``. Both are selected via the ``fit_flux_cls`` / +``fit_time_delays_cls`` inputs of ``FitPointDataset`` / ``AnalysisPoint``, mirroring +``fit_positions_cls``. + +__Missing-image penalty__ + +``FitPositionsImagePairAll``'s mixture normalization already implements the principled +over-prediction Occam factor (the ``1/P^I`` term of Lombardi 2024's mixture likelihood). +``FitPositionsImagePairRepeat``'s ``unmatched_model_policy`` heuristics stay as they are: +best-match pairing is not a normalized mixture, so no principled combinatorial term applies to it. + +__Choosing at cluster scale__ + +Profiling on the standard cluster model (see the likelihood-breakdown scripts referenced above) +puts the analytic ``β*`` solve at timing-noise-level overhead per likelihood call — net +3% on the +image-plane likelihood, +9–30% on sub-0.1 s eager source-plane totals — while removing 2 free +parameters per point source from the non-linear space. The recommendation for cluster fits is +therefore to sharpen the workflow above: **search with ``FitPositionsSourceSolved`` (with +``al.ps.PointSolved`` sources), validate with the image-plane chi-squared** on the max-likelihood +model. Reserve free-centre ``FitPositionsSource`` for direct comparisons with codes that sample the +source position (its ``weighting = "magnification"`` scalar convention matches Lenstool's). __Demonstration__ @@ -178,6 +268,26 @@ class FitStrict(al.FitPositionsImagePairRepeat): print(f" residuals = {[round(float(r), 3) for r in np.asarray(fit.residual_map)]}") print(f" chi_squared = {float(fit.chi_squared):.4f}") +""" +__Case 4 — solved source centre__ + +The solved variants swap the source's ``al.ps.Point`` (free centre) for the parameter-free +``al.ps.PointSolved``. The source-plane fit needs no solver at all: it back-traces the observed +positions and solves the centre analytically. The fit exposes the solved centre via +``source_plane_coordinate`` and the analytic-marginalization contribution via +``marginalization_term``. +""" +source_solved = al.Galaxy(redshift=1.0, point_0=al.ps.PointSolved()) +tracer_solved = al.Tracer(galaxies=[lens, source_solved]) + +fit = al.FitPositionsSourceSolved( + name="point_0", data=data, noise_map=noise_map, tracer=tracer_solved, solver=None +) +print("Case 4 — solved source centre (source-plane fit, no free centre parameters):") +print(f" solved centre beta* = {tuple(round(float(c), 4) for c in fit.source_plane_coordinate)}") +print(f" chi_squared = {float(fit.chi_squared):.4f}") +print(f" marginalization_term = {float(fit.marginalization_term):.4f}") + """ For the production-scale picture — real solver, multi-plane cluster tracer, timings — see ``scripts/cluster/likelihood_function.py`` and the profiling breakdowns referenced above. diff --git a/scripts/point_source/features/fluxes.py b/scripts/point_source/features/fluxes.py index c68ba1a7a..d0df04cc7 100644 --- a/scripts/point_source/features/fluxes.py +++ b/scripts/point_source/features/fluxes.py @@ -184,6 +184,38 @@ fit_positions_cls=al.FitPositionsImagePairRepeat, # Image-plane chi-squared with repeat image pairs. ) +""" +__Analytic Flux (Solved)__ + +An alternative to sampling `flux` is to solve it analytically: `al.FitFluxesSolved` computes the +magnification-weighted best-fit source flux from the data in closed form, so the model needs no +`flux` parameter at all. The model component is then the parameter-free `al.ps.PointSolved` (which +also solves the source centre analytically — solved fits require it, and mixing solved fits with +`Point`/`PointFlux` raises an error), reducing the model above from N=8 to N=5. + +The composition is shown below (not fitted here — the model-fit in this script demonstrates the +free-`flux` convention; see `guides/point_source_pairing.py` for when to prefer each): + +The microlensing caveat at the top of this script applies equally to the solved flux: solving the +source flux analytically does not make image fluxes any less affected by microlensing. +""" +lens_solved = af.Model(al.Galaxy, redshift=0.5, mass=al.mp.Isothermal) + +source_solved = af.Model(al.Galaxy, redshift=1.0, point_0=af.Model(al.ps.PointSolved)) + +model_solved = af.Collection( + galaxies=af.Collection(lens=lens_solved, source=source_solved) +) + +analysis_solved = al.AnalysisPoint( + dataset=dataset, + solver=solver, + fit_positions_cls=al.FitPositionsSourceSolved, # Solved fits pair with `PointSolved`. + fit_flux_cls=al.FitFluxesSolved, # Flux solved analytically, no free `flux` parameter. +) + +print(model_solved.info) + """ __Run Times__ diff --git a/scripts/point_source/features/time_delays.py b/scripts/point_source/features/time_delays.py index ad3b4d303..6dbf546a5 100644 --- a/scripts/point_source/features/time_delays.py +++ b/scripts/point_source/features/time_delays.py @@ -184,6 +184,36 @@ fit_positions_cls=al.FitPositionsImagePairRepeat, # Image-plane chi-squared with repeat image pairs. ) +""" +__Analytic Reference Time (Solved)__ + +The default `FitTimeDelays` compares *relative* delays by subtracting the minimum delay from both +the data and the model (the "reference image" convention). An alternative is +`al.FitTimeDelaysSolved`, which instead solves a precision-weighted analytic reference time from +the data in closed form and analytically marginalizes over it — a smooth (and JAX-differentiable) +alternative to the min-subtraction. + +Solved fit classes pair with the parameter-free `al.ps.PointSolved` model component (which also +solves the source centre, dropping its 2 free parameters; mixing solved fits with `Point` raises an +error). The composition below is shown for reference and not fitted in this script: +""" +lens_solved = af.Model(al.Galaxy, redshift=0.5, mass=al.mp.Isothermal) + +source_solved = af.Model(al.Galaxy, redshift=1.0, point_0=af.Model(al.ps.PointSolved)) + +model_solved = af.Collection( + galaxies=af.Collection(lens=lens_solved, source=source_solved) +) + +analysis_solved = al.AnalysisPoint( + dataset=dataset, + solver=solver, + fit_positions_cls=al.FitPositionsSourceSolved, # Solved fits pair with `PointSolved`. + fit_time_delays_cls=al.FitTimeDelaysSolved, # Analytic reference time instead of min-subtraction. +) + +print(model_solved.info) + """ __Run Times__ diff --git a/scripts/point_source/fit.py b/scripts/point_source/fit.py index da8ff3ead..070f18f94 100644 --- a/scripts/point_source/fit.py +++ b/scripts/point_source/fit.py @@ -32,6 +32,7 @@ - **Name Pairing:** The names of the point-source datasets have an even more important role, the names are used to pair. - **Fitting:** Fit the lens model to the dataset and inspect the results. - **Chi Squared:** For point-source modeling, there are many different ways to define the likelihood function, broadly. +- **Solved Source Centre:** The `*Solved` fit variants solve the source-plane centre analytically instead of sampling it. - **Fluxes:** Another measurable quantity of a point source is its flux—the total amount of light received from. - **Flux Point Dataset:** The fluxes are not input a `PointDataset` object, alongside the image-plane coordinates of the. - **Flux Fitting:** Above, we used a `FitPointDataset` to fit the positions of the point source in the image-plane. @@ -406,6 +407,50 @@ print(fit.positions.log_likelihood) """ +A note on defaults: `FitPointDataset` itself defaults to `FitPositionsImagePair` (Hungarian pairing, +no repeats), whereas `AnalysisPoint` — the object used in the model-fitting examples — defaults to +`FitPositionsImagePairRepeat`. The examples above pass `fit_positions_cls` explicitly so there is no +ambiguity about which chi-squared is being used. + +__Solved Source Centre__ + +Every fit above reads the source-plane centre from the model's `Point` profile, whose `centre` is a +free parameter during model-fitting. Each fit class also has a `*Solved` variant which instead solves +the centre analytically from the observed positions and the mass model, removing those 2 free +parameters per point source from the non-linear search: `FitPositionsSourceSolved` (a +precision-weighted source-plane centre, following Lombardi 2024, arXiv:2406.15280) and +`FitPositionsImagePairRepeatSolved` / `FitPositionsImagePairAllSolved` (the same solved centre +driving the image-plane forward solve). + +The solved variants pair with the parameter-free `al.ps.PointSolved` model component instead +of `Point` — mixing the two conventions (a `Solved` fit with a `Point`, or vice versa) raises +a `PointProfileMismatchException`, so a model cannot silently sample centre parameters the fit +ignores. +""" +point_source_solved = al.ps.PointSolved() + +source_galaxy_solved = al.Galaxy(redshift=1.0, point_0=point_source_solved) + +tracer_solved = al.Tracer(galaxies=[lens_galaxy, source_galaxy_solved]) + +fit = al.FitPointDataset( + dataset=dataset, + tracer=tracer_solved, + solver=solver, + fit_positions_cls=al.FitPositionsSourceSolved, # Solved-centre source-plane chi-squared +) + +print("Analytically Solved Source-Plane Centre:") +print(fit.positions.source_plane_coordinate) + +print("Log Likelihood with Solved Centre:") +print(fit.positions.log_likelihood) + +""" +The full matrix of pairing schemes × centre treatments — including the solved flux and time-delay +variants used below, the tensor-vs-scalar magnification weighting, and when each option is +preferable — is documented in `guides/point_source_pairing.py`. + __Fluxes__ Another measurable quantity of a point source is its flux—the total amount of light received from each multiple image @@ -503,6 +548,10 @@ print(fit.flux.log_likelihood) """ +An analytic-flux variant also exists: `fit_flux_cls=al.FitFluxesSolved` (paired with +`al.ps.PointSolved`) solves the source flux from the data and magnifications instead of sampling +a `flux` parameter — see `point_source/features/fluxes.py` and `guides/point_source_pairing.py`. + __Time Delays__ Another measurable quantity of a point source is its time delay—the time it takes for light to travel from the @@ -588,6 +637,10 @@ print(fit.time_delays.log_likelihood) """ +As with fluxes, an analytic variant exists: `fit_time_delays_cls=al.FitTimeDelaysSolved` replaces +the reference-image subtraction with a precision-weighted analytic reference time — see +`point_source/features/time_delays.py` and `guides/point_source_pairing.py`. + __New User Wrap Up__ The `point_source` package of the `autolens_workspace` contains numerous example scripts for performing point source diff --git a/scripts/point_source/modeling.py b/scripts/point_source/modeling.py index c879489bf..55c29dcc4 100644 --- a/scripts/point_source/modeling.py +++ b/scripts/point_source/modeling.py @@ -182,10 +182,14 @@ The API is fairly self explanatory and is straight forward to extend, for example adding more light profiles to the lens and source or using a different mass profile. -A full description of model composition is provided by the model cookbook: +A full description of model composition is provided by the model cookbook: https://pyautolens.readthedocs.io/en/latest/general/model_cookbook.html +A centre-free alternative exists: `al.ps.PointSolved` has no free parameters — the source centre is +solved analytically by the paired `*Solved` fit classes (e.g. `fit_positions_cls=al.FitPositionsSourceSolved`), +removing 2 parameters per point source. See `guides/point_source_pairing.py` for the full option matrix. + __Name Pairing__ Every point-source dataset in the `PointDataset` has a name, which in this example was `point_0`. This `name` pairs diff --git a/workspace_index.json b/workspace_index.json index a69372de2..b5b0d885d 100644 --- a/workspace_index.json +++ b/workspace_index.json @@ -108,6 +108,7 @@ "autolens_workspace/scripts/guides/point_source_pairing.py", "autolens_workspace_test/scripts/cluster/likelihood_sanity.py", "csv_api.py", + "guides/point_source_pairing.py", "modeling.py", "scripts/cluster/csv_api.py", "scripts/cluster/modeling.py", @@ -172,6 +173,7 @@ "/cluster/simulator.py", "autolens_workspace_test/scripts/cluster/visualization.py", "cluster/simulator.py", + "guides/point_source_pairing.py", "point_source/start_here.py", "scripts/cluster/csv_api.py", "scripts/group/features/scaling_relation/modeling.py", @@ -268,6 +270,7 @@ "cluster/modeling.py", "csv_api.py", "group/start_here.ipynb", + "guides/point_source_pairing.py", "imaging/start_here.ipynb", "multi_galaxy/start_here.ipynb", "prep.py", @@ -6942,6 +6945,7 @@ "Result" ], "cross_refs": [ + "guides/point_source_pairing.py", "point_source/start_here.ipynb", "start.here.py", "start_here.ipynb", @@ -7129,6 +7133,7 @@ "Name Pairing", "Fitting", "Chi Squared", + "Solved Source Centre", "Fluxes", "Flux Point Dataset", "Flux Fitting", @@ -7138,7 +7143,10 @@ "Shape Solver" ], "cross_refs": [ + "guides/point_source_pairing.py", "modeling.py", + "point_source/features/fluxes.py", + "point_source/features/time_delays.py", "scripts/point_source/simulator.py", "simulator.py", "start_here.py" @@ -7172,6 +7180,7 @@ ], "cross_refs": [ "/point_source/start_here.py", + "guides/point_source_pairing.py", "start_here.py" ], "notebook": "notebooks/point_source/modeling.ipynb",