From 996d2c7b2c69e9348d57e74372dd282bd16c2550 Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Thu, 9 Jul 2026 14:56:48 +0100 Subject: [PATCH] Add real-data A2744 weak lensing example (weak series step 7b) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit scripts/weak/real_data/a2744.py — the first PyAutoLens weak-lensing analysis of real sky data: downloads the public pyRRG A2744 shape catalogue (pinned commit, provenance documented; not redistributed), projects RA/Dec to the tangent plane, applies standard quality cuts (6585 -> 382 usable shapes), builds a reduced-shear WeakDataset via the new from_arrays loader, shows the model-independent Kaiser-Squires map, and fits a spherical NFW with Nautilus. Framed honestly: training-data-grade depth gives a marginal (~1.5 sigma) tangential signal with the correct radial concentration and a null B-mode — the example demonstrates the real-data workflow, which transfers unchanged to deep survey catalogues via WeakDataset.from_fits/from_csv. Step 7b of the weak-lensing series (PyAutoLabs/PyAutoLens#588). Co-Authored-By: Claude Fable 5 --- llms-full.txt | 2 + scripts/weak/real_data/__init__.py | 0 scripts/weak/real_data/a2744.py | 278 +++++++++++++++++++++++++++++ workspace_index.json | 16 ++ 4 files changed, 296 insertions(+) create mode 100644 scripts/weak/real_data/__init__.py create mode 100644 scripts/weak/real_data/a2744.py diff --git a/llms-full.txt b/llms-full.txt index 4dca8cef3..f79e17a50 100644 --- a/llms-full.txt +++ b/llms-full.txt @@ -442,6 +442,8 @@ AUTO-GENERATED by PyAutoBuild — do not edit by hand; regenerate with generate. - Contents: Dataset, Lens Galaxy, Shear Field Evaluation, Residuals, Chi Squared, Noise Normalization Term, Calculate The Log Likelihood, Fit, Analysis, Wrap Up - [Modeling: Weak Lensing](scripts/weak/modeling.py): This script fits a lens mass model to a weak gravitational lensing shear catalogue using a non-linear search. It is the weak-lensing analogue of `scripts/imaging/modeling.py`: where that script infers a lens model from a 2D image of a lensed source, this script infers one from the (gamma_2, gamma_1) shear measurements of background source galaxies — a `WeakDataset` produced by `scripts/weak/simulator.py`. - Contents: Dataset, Model, Search, Analysis, Run Times, Model-Fit, Output Folder Layout, Result +- [Real Data: Weak Lensing of Abell 2744](scripts/weak/real_data/a2744.py): This script fits a dark-matter halo mass model to a **real weak-lensing shape catalogue** of the merging galaxy cluster Abell 2744 ("Pandora's Cluster", z = 0.308) — the first PyAutoLens weak-lensing analysis of real sky data, and the capstone of the weak-lensing example series. + - Contents: Catalogue Download, Catalogue Load & Projection, Quality Cuts, Weak Dataset, Mass Map, Model & Search, Result - [Simulator: Weak Lensing](scripts/weak/simulator.py): This script simulates a weak gravitational lensing shear catalogue. Unlike the imaging simulator (which produces a 2D image of the lensed source) the weak-lensing simulator produces a *catalogue* of (gamma_2, gamma_1) shear measurements at the (y, x) positions of a population of background source galaxies. - Contents: Model, Dataset Paths, Ray Tracing, Source Positions, Simulator, Output, Visualize diff --git a/scripts/weak/real_data/__init__.py b/scripts/weak/real_data/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/scripts/weak/real_data/a2744.py b/scripts/weak/real_data/a2744.py new file mode 100644 index 000000000..4910b7c50 --- /dev/null +++ b/scripts/weak/real_data/a2744.py @@ -0,0 +1,278 @@ +""" +Real Data: Weak Lensing of Abell 2744 +===================================== + +This script fits a dark-matter halo mass model to a **real weak-lensing shape catalogue** of the merging +galaxy cluster Abell 2744 ("Pandora's Cluster", z = 0.308) — the first PyAutoLens weak-lensing analysis of +real sky data, and the capstone of the weak-lensing example series. + +__The Data & Its Provenance__ + +The catalogue is the Abell 2744 galaxy shape catalogue shipped with the public pyRRG weak-lensing shape +measurement code (https://github.com/davidharvey1986/pyRRG, `jwst` branch), whose JWST application to this +cluster is described in Harvey & Massey 2024 (MNRAS 529, 802, arXiv:2401.16478). It contains 6,585 sources +with SExtractor photometry, RRG shape moments and per-galaxy ellipticity measurements `(e1, e2)` with +uncertainties. + +**Honest framing:** the file ships inside pyRRG's star/galaxy-classification training data, and we apply our +own quality cuts below rather than the exact selections of the published analysis — so this example is a +*real-data demonstration* whose results we sanity-check against published Abell 2744 analyses, not a +metrology-grade reproduction of any single paper. That is the right expectation for a first weak-lensing +look at any new catalogue. + +__What this script shows__ + +- Downloading and loading a real shape catalogue (FITS binary table) and converting RA/Dec to the + tangent-plane arc-second coordinates PyAutoLens uses. +- The standard weak-lensing quality cuts every real catalogue needs. +- Building a `WeakDataset` flagged `is_reduced` — real surveys measure galaxy ellipticities, i.e. the + *reduced* shear g = gamma / (1 - kappa), and `FitWeak` computes the matching model quantity. +- A model-independent Kaiser-Squires mass map — Abell 2744 is a famous merger, and the map's structure is + the first sanity check. +- Fitting a spherical NFW dark-matter halo with a Nautilus search and reading the tangential-shear profile + in which cluster weak-lensing results are usually shown. + +__Contents__ + +- **Catalogue Download:** Fetch the public catalogue (cached on disk after the first run). +- **Catalogue Load & Projection:** RA/Dec to tangent-plane arc-seconds about the cluster centre. +- **Quality Cuts:** The standard selections that turn raw shapes into a usable shear sample. +- **Weak Dataset:** Build the reduced-shear `WeakDataset`. +- **Mass Map:** Model-independent Kaiser-Squires reconstruction. +- **Model & Search:** Spherical NFW halo, Nautilus. +- **Result:** Tangential-shear profile, and how the numbers compare to the literature. +""" + +from autoconf import jax_wrapper # Sets JAX environment before other imports + +# from autoconf import setup_notebook; setup_notebook() + +from pathlib import Path + +import numpy as np + +import autofit as af +import autolens as al +import autolens.plot as aplt + +""" +__Catalogue Download__ + +The catalogue is fetched once from the public pyRRG repository (pinned to a specific commit for +reproducibility) and cached in the dataset folder. We do not redistribute the file with the workspace — +provenance stays with the pyRRG project. +""" +dataset_path = Path("dataset") / "weak" / "a2744_pyrrg" +catalogue_path = dataset_path / "abell2744_galaxies.fits" + +CATALOGUE_URL = ( + "https://raw.githubusercontent.com/davidharvey1986/pyRRG/" + "0ccc29fb4513137da61b1afb632ca492093bd609/" + "trainStarGalClass/TrainingData/abell2744_galaxies.fits" +) + +if not catalogue_path.exists(): + import urllib.request + + dataset_path.mkdir(parents=True, exist_ok=True) + print(f"Downloading A2744 catalogue from pyRRG (one-off, ~3 MB) ...") + urllib.request.urlretrieve(CATALOGUE_URL, catalogue_path) + +""" +__Catalogue Load & Projection__ + +The table stores sky positions as RA/Dec in degrees. PyAutoLens works in tangent-plane arc-second offsets +`(y, x)` about a chosen centre, so we project about the cluster core (the catalogue's density peak, +consistent with the BCG region used by published analyses): + + - `x = (RA - RA0) * cos(Dec0) * 3600` (arc-seconds East) + - `y = (Dec - Dec0) * 3600` (arc-seconds North) + +A note on conventions: whether East points left or right on the sky is a *parity* choice that rotates or +mirrors the shear components' frame. The tangential shear — the quantity our fit constrains — is invariant +under this mirror (only the B-mode cross component flips sign), so the halo-profile fit below is robust to +it. Precision studies of shear *systematics* must track the convention carefully. +""" +from astropy.io import fits as astropy_fits + +with astropy_fits.open(catalogue_path) as hdul: + table = hdul[1].data + +ra = np.asarray(table["ra"], dtype=float) +dec = np.asarray(table["dec"], dtype=float) +e1 = np.asarray(table["e1"], dtype=float) +e2 = np.asarray(table["e2"], dtype=float) +e1_err = np.asarray(table["e1_err"], dtype=float) +e2_err = np.asarray(table["e2_err"], dtype=float) + +ra_centre, dec_centre = 3.5875, -30.3972 # A2744 core (J2000 degrees) + +x = (ra - ra_centre) * np.cos(np.deg2rad(dec_centre)) * 3600.0 +y = (dec - dec_centre) * 3600.0 +radii = np.sqrt(x**2.0 + y**2.0) + +print(f"catalogue sources : {len(ra)}") + +""" +__Quality Cuts__ + +Raw shape catalogues always contain unusable measurements — blends, noise detections, objects whose moments +diverged. The cuts below are the standard minimum for any weak-lensing sample: + + - finite, physical ellipticities: |e1|, |e2| < 1 (a galaxy ellipticity cannot exceed 1; the raw catalogue + contains outliers far beyond this from failed moment measurements). + - measured uncertainties in a sane range: 0 < e_err < 0.4 per component. + - a radial window 10" < r < 130": inside ~10" we are in the strong-lensing core where cluster members + dominate and the weak-lensing (linear shear) approximation is worst; ~130" is the edge of this + catalogue's contiguous coverage. + +On this catalogue the cuts are severe: only ~1,600 of the 6,585 sources have measured shapes at all, and +the physical-ellipticity cut trims those to ~400 — this is training-data-grade depth, an order of magnitude +shallower than the selections behind published A2744 analyses. Keep that in mind when reading the results. +""" +finite = np.isfinite(e1) & np.isfinite(e2) & np.isfinite(e1_err) & np.isfinite(e2_err) +physical = (np.abs(e1) < 1.0) & (np.abs(e2) < 1.0) +well_measured = (e1_err > 0.0) & (e1_err < 0.4) & (e2_err > 0.0) & (e2_err < 0.4) +radial = (radii > 10.0) & (radii < 130.0) + +use = finite & physical & well_measured & radial + +print(f"after quality cuts : {use.sum()}") + +""" +__Weak Dataset__ + +The per-galaxy noise combines the intrinsic shape dispersion (each galaxy has a random unlensed ellipticity; +sigma_int ~ 0.25 per component is the standard value) with the measurement uncertainty, in quadrature. + +`from_arrays` builds the `WeakDataset`; `is_reduced=True` (the loader default) records that these are +measured ellipticities — reduced shear — so `FitWeak` will compare them against the model's +g = gamma / (1 - kappa), not the bare shear. +""" +sigma_int = 0.25 + +noise = np.sqrt(sigma_int**2.0 + 0.5 * (e1_err[use] ** 2.0 + e2_err[use] ** 2.0)) + +dataset = al.WeakDataset.from_arrays( + positions=np.stack([y[use], x[use]], axis=1), + gamma_1=e1[use], + gamma_2=e2[use], + noise_map=list(noise), + is_reduced=True, + name="a2744_pyrrg", +) + +print(dataset.info) + +aplt.subplot_weak_dataset(dataset=dataset, output_path=dataset_path, output_format="png") + +""" +__Mass Map__ + +Before any model is fitted, the Kaiser-Squires inversion gives a model-independent mass map. Abell 2744 is +one of the most disturbed clusters known — published lensing maps (Merten et al. 2011; Medezinski et al. +2016; Harvey & Massey 2024) show multiple substructures around the main core from an ongoing merger — so we +should *not* expect a clean single peak, and the structure in this map is the first indication the +catalogue's shear signal is real. +""" +aplt.plot_convergence_map( + shear_yx=dataset.shear_yx, + shape_native=(30, 30), + smoothing_sigma_pixels=1.5, + output_path=dataset_path, + output_format="png", +) + +""" +__Model & Search__ + +The model is a spherical NFW dark-matter halo — the standard first-order description of a cluster halo and +deliberately simple for a merging system (the published analyses use multiple halos; a single NFW measures +the dominant mass concentration). + +The halo centre gets Gaussian priors of width 10" about the projected cluster core, and the fit assumes a +single effective source plane at z = 1.0 behind the z = 0.308 cluster (the catalogue provides no per-galaxy +redshifts; this is the standard effective-depth approximation and its choice rescales the inferred halo +normalisation). +""" +mass = af.Model(al.mp.NFWSph) +mass.centre.centre_0 = af.GaussianPrior(mean=0.0, sigma=10.0) +mass.centre.centre_1 = af.GaussianPrior(mean=0.0, sigma=10.0) + +lens = af.Model(al.Galaxy, redshift=0.308, mass=mass) + +source = af.Model(al.Galaxy, redshift=1.0) + +model = af.Collection(galaxies=af.Collection(lens=lens, source=source)) + +print(model.info) + +analysis = al.AnalysisWeak(dataset=dataset) + +search = af.Nautilus( + path_prefix=Path("weak") / "real_data", + name="a2744_nfw", + unique_tag="a2744_pyrrg", + n_live=100, + iterations_per_quick_update=10000, +) + +print( + """ + The non-linear search has begun running — a few hundred galaxies and a 4-parameter model, + so expect minutes, not hours. + """ +) + +result = search.fit(model=model, analysis=analysis) + +print("The search has finished run - you may now continue the notebook.") + +""" +__Result__ + +The tangential-shear profile is the standard presentation of a cluster weak-lensing measurement: binned +data (with the cross-component B-mode null test) against the maximum-likelihood NFW curve. + +What this shallow sample can and cannot show — read the numbers with the sample size in mind: + + - With only ~400 usable shapes, the overall tangential-shear detection is *marginal* (a weighted mean + gamma_t of ~0.02 at ~1.5 sigma), though it behaves exactly as a real lensing signal should: it + concentrates at small radii (gamma_t is several times larger inside 40" than outside) and the + cross-component B-mode is consistent with zero. + - The NFW posterior is correspondingly broad: the halo normalisation and scale radius are each uncertain + at the tens-of-percent-to-factors level, and the centre is only loosely pinned near the projected core. + That *is* the honest result of fitting ~400 galaxies around one cluster — published A2744 analyses + (Medezinski et al. 2016 quote a virial mass ~2 x 10^15 solar masses; Harvey & Massey 2024 map the + merger's substructure) rest on samples an order of magnitude deeper with survey-grade calibration. + - What the example therefore demonstrates is the *workflow* on real sky data — download, projection, + quality cuts, reduced-shear dataset, model-independent map, likelihood fit, profile — which transfers + unchanged to a deep catalogue via `WeakDataset.from_fits` / `from_csv`. +""" +print(result.info) + +aplt.subplot_fit_weak( + fit=result.max_log_likelihood_fit, output_path=dataset_path, output_format="png" +) + +aplt.plot_shear_profile( + result.max_log_likelihood_fit, + centre=(0.0, 0.0), + bins=8, + output_path=dataset_path, + output_format="png", +) + +aplt.corner_anesthetic(samples=result.samples) + +""" +__Wrap Up__ + +The weak-lensing series is now complete end to end: simulator, visualization, fit, modeling, likelihood +guide, combined strong+weak analysis — and real sky data. From here: + + - Replace the single NFW with a multi-halo model (`scripts/cluster`) to chase A2744's substructures. + - Combine this shear catalogue with the cluster's strong-lensing constraints, exactly as in + `scripts/weak/features/strong_lensing` — the hybrid approach of Niemiec et al. 2020. + - Swap in your own survey's catalogue via `WeakDataset.from_fits` / `from_csv`. +""" diff --git a/workspace_index.json b/workspace_index.json index a1b989587..797a00451 100644 --- a/workspace_index.json +++ b/workspace_index.json @@ -5916,6 +5916,22 @@ "summary": "This script fits a lens mass model to a weak gravitational lensing shear catalogue using a non-linear search. It is the weak-lensing analogue of `scripts/imaging/modeling.py`: where that script infers a lens model from a 2D image of a lensed source, this script infers one from the (gamma_2, gamma_1) shear measurements of background source galaxies \u2014 a `WeakDataset` produced by `scripts/weak/simulator.py`.", "title": "Modeling: Weak Lensing" }, + { + "contents": [ + "Catalogue Download", + "Catalogue Load & Projection", + "Quality Cuts", + "Weak Dataset", + "Mass Map", + "Model & Search", + "Result" + ], + "cross_refs": [], + "notebook": null, + "path": "scripts/weak/real_data/a2744.py", + "summary": "This script fits a dark-matter halo mass model to a **real weak-lensing shape catalogue** of the merging galaxy cluster Abell 2744 (\"Pandora's Cluster\", z = 0.308) \u2014 the first PyAutoLens weak-lensing analysis of real sky data, and the capstone of the weak-lensing example series.", + "title": "Real Data: Weak Lensing of Abell 2744" + }, { "contents": [ "Model",