From 38dcceb4c764f3d1e642312499931481acc495bd Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Mon, 20 Jul 2026 13:51:12 +0100 Subject: [PATCH] Add cluster mass-parameterization guide scripts MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Two companion guides for cluster / group (multi-galaxy lens) mass parameterization, using the model-composition API with free/fixed made explicit and numerical checks throughout: - cluster/mass_parameterizations.py — a tour of the standard Lenstool-style parameterizations (dPIE sigma vs b0; Faber-Jackson scaling relations; mass-anchoring to the BCG; the two-parameter relation freeing r_cut_ref), each model self-contained (halo + individual galaxies + scaling galaxies). - cluster/mass_parameterizations_pyautolens.py — for Lenstool users: maps Lenstool Model 1 onto PyAutoLens's native MGL parameterization. Section 1 is the Lenstool model; Section 2 the standard PyAutoLens mapping (halo dPIE->NFW, galaxies dPIE->Isothermal, sigma->einstein_radius, and the main/extra/scaling coupling spectrum: free / luminosity-bounded / BGC-tied); Section 3 keeps the dPIE truncation on the scaling galaxies while main and extra galaxies stay Isothermal, with the cross-parameterization b0<-> einstein_radius linking. Notebooks regenerated; catalogue (llms-full.txt, workspace_index.json) updated. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01K5MpYSYekRWVtoUmExAgKB --- llms-full.txt | 5 +- .../cluster/mass_parameterizations.ipynb | 741 ++++++++++++++++++ .../mass_parameterizations_pyautolens.ipynb | 509 ++++++++++++ scripts/cluster/mass_parameterizations.py | 593 ++++++++++++++ .../mass_parameterizations_pyautolens.py | 393 ++++++++++ workspace_index.json | 23 + 6 files changed, 2263 insertions(+), 1 deletion(-) create mode 100644 notebooks/cluster/mass_parameterizations.ipynb create mode 100644 notebooks/cluster/mass_parameterizations_pyautolens.ipynb create mode 100644 scripts/cluster/mass_parameterizations.py create mode 100644 scripts/cluster/mass_parameterizations_pyautolens.py diff --git a/llms-full.txt b/llms-full.txt index ecfddf9c1..6b85fd975 100644 --- a/llms-full.txt +++ b/llms-full.txt @@ -1,4 +1,4 @@ -AUTO-GENERATED by PyAutoBuild — do not edit by hand; regenerate with generate.py. +AUTO-GENERATED by PyAutoHands — do not edit by hand; regenerate with generate.py. # PyAutoLens Workspace — Full Catalogue @@ -390,6 +390,9 @@ AUTO-GENERATED by PyAutoBuild — do not edit by hand; regenerate with generate. - [Lenstool Users: dPIE Parameterization Mapping (Scaling Galaxies)](scripts/cluster/lenstool/parameterization_mapping.py): **If you want to fit cluster member (scaling) galaxies without knowing the source or lens redshift, this script shows you how — and proves it recovers the standard Lenstool model exactly once the redshifts become known.** - [Log Likelihood Function: Cluster Point Source](scripts/cluster/likelihood_function.py): This script provides a step-by-step guide of the cluster point-source ``log_likelihood_function``, the figure-of-merit Nautilus optimises when fitting a cluster lens model to ``point_datasets.csv``. - Contents: Dataset, Truth Model, Tracer, Source-Plane Chi Squared: Concept, Multi-Plane Ray Tracing, Back-Traced Source-Plane Positions, Source-Plane Centroid, Residual Map, Magnifications at Positions, Chi Squared Map (Source), Per-Source Chi Squared, Total Chi Squared (Source), Noise Normalization (Source), Source-Plane Log Likelihood, Source-Plane Validation, Image-Plane Chi Squared: Concept, Point Solver Setup, Forward Solving Model Positions, Pairing Model to Observed, Image-Plane Residual Map, Chi Squared Map (Image), Per-Source / Total Chi Squared, Noise Normalization (Image), Image-Plane Log Likelihood, Image-Plane Validation, Source-Plane vs Image-Plane: When to Use Which, Wrap Up +- [Cluster & Group Mass Parameterizations](scripts/cluster/mass_parameterizations.py): **A plain-language tour of the ways to parameterize and couple the masses in a Lenstool-style cluster model — the cluster-scale halo, the BCG and other individually-freed galaxies, and the member ("scaling") population.** + - Contents: Model 1 — The Standard Lenstool Model, Model 2 — The Angular (b0) Parameterization, Model 3 — Mass Anchoring, Model 4 — The Two-Parameter Relation +- [Mass Parameterizations II — Mapping Lenstool onto PyAutoLens](scripts/cluster/mass_parameterizations_pyautolens.py): **For Lenstool users: how and why PyAutoLens re-expresses the standard Lenstool cluster / group model in its native parameterization for multi-galaxy lenses (MGLs), and how to keep the dPIE's truncation on the outskirts if you want it.** - [Modeling: Cluster](scripts/cluster/modeling.py): This script models the small multi-plane cluster lens simulated by ``cluster/simulator.py``: a Brightest Cluster Galaxy (BCG) plus one satellite member at the lens redshift ``z = 0.5``, 10 lower-mass cluster members modelled collectively via a luminosity-mass scaling relation, a standalone NFW host dark matter halo *not* tied to any individual galaxy, and 2 background sources at *different* redshifts (``z = 1.0`` and ``z = 2.0``) which the cluster lenses into multiple images. - Contents: Example, Simulation, Dataset, Model CSVs, Scaling Galaxies Table, Point Solver, Chi Squared, Cluster Components, Redshifts, Model, Scaling Relation, Name Pairing, Search, Live Visual Update, Analysis, Factor Graph, Run Times, Output Folder Layout, Result - [Simulator: Cluster](scripts/cluster/simulator.py): This script simulates an example strong lens on the 'cluster' scale: a small cluster consisting of 2 main lens galaxies (a brightest cluster galaxy + a single satellite), 10 lower-mass cluster member galaxies on a luminosity-mass scaling relation, a single host dark matter halo not tied to any individual galaxy, and 2 multiply-imaged background source galaxies sitting at *different* redshifts (``z = 1.0`` and ``z = 2.0``) — making this a genuine multi-plane lens. diff --git a/notebooks/cluster/mass_parameterizations.ipynb b/notebooks/cluster/mass_parameterizations.ipynb new file mode 100644 index 000000000..ebfbf7ff8 --- /dev/null +++ b/notebooks/cluster/mass_parameterizations.ipynb @@ -0,0 +1,741 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Cluster & Group Mass Parameterizations\n", + "======================================\n", + "\n", + "**A plain-language tour of the ways to parameterize and couple the masses in a Lenstool-style cluster\n", + "model \u2014 the cluster-scale halo, the BCG and other individually-freed galaxies, and the member\n", + "(\"scaling\") population.**\n", + "\n", + "There are several parameterizations floating around \u2014 velocity dispersion vs the internal ``b0``;\n", + "magnitudes vs luminosities; members tied to a scaling relation vs left free; the members normalized to\n", + "their own free parameter vs to the BCG's mass. They describe the *same physics*, so it is easy to lose\n", + "track of what actually differs. This script builds each option **with the model-composition API and\n", + "prints ``model.info``**, so for every model you can see at a glance *exactly which parameters are free\n", + "and which are fixed* \u2014 the thing that really distinguishes the approaches.\n", + "\n", + "Nothing here is fitted. Each section composes an ``af.Collection`` model and prints it (plus a few\n", + "assertions); the point is to *read* the models, not run them.\n", + "\n", + "__The axes of choice__\n", + "\n", + "Every model below is a point in the same space of decisions:\n", + "\n", + " - **What number carries the mass?** ``sigma`` (velocity dispersion, km/s \u2014 Lenstool's native\n", + " parameter) or the internal ``b0`` (arcsec, the dPIE lens strength). Inter-convertible; which you\n", + " make *free* is the choice.\n", + " - **How are members coupled?** A Faber-Jackson scaling relation ties the whole member population to a\n", + " shared normalization; individual galaxies are left free (only feasible for a few).\n", + " - **What anchors the relation?** A reference *magnitude* ``mag0`` (Lenstool's ``potfile`` keyword).\n", + " - **What is the member normalization?** Its own free parameter (``sigma_ref``), or the BCG's own mass\n", + " (mass anchoring \u2014 the member tier then costs zero parameters).\n", + " - **How many relation parameters?** Just the dispersion normalization, or the truncation\n", + " normalization ``r_cut_ref`` as well.\n", + "\n", + "__Contents__\n", + "\n", + "Each model is a **complete, self-contained** Lenstool-style cluster model \u2014 a cluster-scale halo, the\n", + "individually-modelled galaxies, and the member relation \u2014 so you can copy any one block wholesale\n", + "(there is deliberate redundancy between them). All use dPIE / Lenstool parameterizations only; mapping\n", + "these onto PyAutoLens-native profiles (isothermals, etc.) is a separate guide.\n", + "\n", + " - **Model 1 \u2014 The Standard Lenstool Model** \u2014 members on a velocity-dispersion Faber-Jackson relation\n", + " anchored to ``mag0``; BCG + a freed member modelled individually; a free elliptical halo. No ``b0``.\n", + " - **Model 2 \u2014 The Angular (b0) Parameterization** \u2014 the same model with ``sigma -> b0`` (redshift-free).\n", + " - **Model 3 \u2014 Mass Anchoring** \u2014 the members normalized to the BCG's own mass (member tier -> 0 free).\n", + " - **Model 4 \u2014 The Two-Parameter Relation** \u2014 free the truncation normalization ``r_cut_ref`` too.\n", + "\n", + "Freeing the Faber-Jackson *exponent* is a one-line change (make ``sigma_exponent`` a prior instead of\n", + "``0.25``) and is noted inline rather than given its own model. How Lenstool handles *unknown source\n", + "redshifts* (``z_m_limit``) is about the source, not the mass, and is out of scope here.\n", + "\n", + "__Reading model.info__\n", + "\n", + "For every model we print ``model.info``. A line showing a ``UniformPrior`` (or any prior) is a **free**\n", + "parameter the search would sample; a bare number is **fixed**. A member's ``sigma`` shown as a\n", + "``MultiplePrior`` is *tied* to a shared free parameter, not free in its own right." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Google Colab Setup__\n", + "\n", + "This cell sets up the environment when the notebook is run on Google Colab: it installs the\n", + "required PyAuto packages, clones the workspace (configuration files and example datasets) and\n", + "points the configuration at it. If you are running the notebook elsewhere (e.g. locally via\n", + "your own installation) it does nothing, and you can run it safely.\n", + "\n", + "Colab tip: model-fits run much faster on a GPU \u2014 enable one via \"Runtime\" -> \"Change runtime\n", + "type\" -> \"Hardware accelerator\" before running the notebook." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "try:\n", + " import google.colab\n", + " import subprocess\n", + " import sys\n", + "\n", + " subprocess.check_call(\n", + " [sys.executable, \"-m\", \"pip\", \"install\", \"autonerves\", \"--no-deps\"]\n", + " )\n", + "except ImportError:\n", + " pass\n", + "\n", + "from autonerves import setup_colab\n", + "\n", + "setup_colab.setup(\"autolens\")" + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "\n", + "import numpy as np\n", + "\n", + "import autofit as af\n", + "import autolens as al\n", + "import autogalaxy as ag" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Shared Setup__\n", + "\n", + "Redshifts and cosmology are the same as the ``cluster/simple`` example. The cosmology is the flat\n", + "``H0 = 67.66``, ``Om0 = 0.30966`` that the dPIE profile uses internally, so nothing is approximate.\n", + "For a real cluster, set these to the measured lens redshift and the run's cosmology.\n", + "\n", + "The photometry below is illustrative. In a real analysis it comes from a member catalogue (one\n", + "magnitude per galaxy in the band the scaling relation is calibrated in \u2014 e.g. F160W) plus each\n", + "galaxy's sky position. ``scaling_*`` are the fainter cluster members that go on the scaling relation;\n", + "``individual_*`` are the galaxies modelled on their own \u2014 the BCG (first) plus any bright member worth\n", + "freeing (e.g. a galaxy near a multiple image). The BCG's magnitude also anchors the relation (mag0)." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "redshift_lens = 0.5\n", + "redshift_source = 2.0\n", + "\n", + "H0 = 67.66\n", + "Om0 = 0.30966\n", + "\n", + "# Scaling-relation members (the \"potfile\" population): (y, x) arc-second position + apparent magnitude.\n", + "scaling_centres = [(5.5, -6.5), (-7.5, 3.0), (3.0, 13.0), (12.0, -5.0), (-6.5, 11.0)]\n", + "scaling_magnitudes = [19.2, 19.8, 20.4, 21.0, 21.6]\n", + "\n", + "# Individually-modelled galaxies: the BCG (first) plus any bright member we want to free \u2014 e.g. a\n", + "# galaxy sitting near a multiple image, which the scaling relation cannot describe accurately enough\n", + "# locally. In Lenstool these are just separate \"potentiel\" sections; the BCG is not structurally\n", + "# special, only the first of them. The BCG's magnitude also anchors the scaling relation (mag0).\n", + "individual_centres = [(0.0, 0.0), (8.5, 5.5)]\n", + "individual_magnitudes = [17.8, 18.9]\n" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Model 1 \u2014 The Standard Lenstool Model__\n", + "\n", + "This is the model a Lenstool ``.par`` file describes. It has the three components of a real cluster\n", + "model, and this script builds all three:\n", + "\n", + " - a **cluster-scale dark-matter halo** \u2014 an elliptical dPIE, fully free, not tied to any light: the\n", + " dominant, smooth large-scale mass;\n", + " - **individually-modelled galaxies** \u2014 the BCG and any bright member worth freeing, each its own free\n", + " dPIE (Lenstool's individual \"potentiel\" sections; the BCG is simply the first);\n", + " - the **member population** (\"potfile\") \u2014 masses set by a velocity-dispersion Faber-Jackson relation\n", + " anchored to a reference magnitude.\n", + "\n", + "dPIE profiles throughout, composed directly in ``sigma`` (km/s), ``r_core`` and ``r_cut`` (arcsec) \u2014\n", + "exactly the numbers a Lenstool results table quotes. No ``b0`` appears.\n", + "\n", + "**The reference magnitude ``mag0``.** The relation is anchored to a reference magnitude (the ``potfile``\n", + "keyword ``mag0``); ``sigma_ref`` is the velocity dispersion of a galaxy *at that magnitude*. A member of\n", + "magnitude ``m`` enters through its brightness *relative* to the reference,\n", + "``L/L_ref = 10 ** (0.4 * (mag0 - m))``, which folds into Faber-Jackson (``sigma ~ L^0.25``) to give a\n", + "clean magnitude form:\n", + "\n", + " sigma_i = sigma_ref * 10 ** (0.1 * (mag0 - m_i)) # 0.25 * 0.4 = 0.1\n", + " r_cut_i = r_cut_ref * 10 ** (0.2 * (mag0 - m_i)) # 0.5 * 0.4 = 0.2 (constant M/L)\n", + " r_core_i = r_core_ref * 10 ** (0.2 * (mag0 - m_i))\n", + "\n", + "We set ``mag0`` to the BCG's magnitude \u2014 a natural anchor. Note this uses only the BCG's *brightness*;\n", + "the BCG's *mass* is one of the individually-modelled galaxies and is never coupled to the members.\n", + "\n", + "**Free vs fixed.**\n", + "\n", + " - Halo: free ``ellipticity``, ``angle_pos``, ``sigma``, ``r_core`` (4). Fixed centre (often freed near\n", + " the BCG) and ``r_cut`` (unconstrained within the field).\n", + " - Individually-modelled galaxies (BCG + freed member): each free ``sigma`` + ``r_cut`` (2 per galaxy);\n", + " fixed centre and ``r_core``.\n", + " - Members: one free parameter for the whole population, ``sigma_ref``. Adding members adds nothing.\n", + " - The exponents, reference radii, ``mag0``, redshifts and cosmology are fixed.\n", + "\n", + "Total: 9 free parameters (halo 4 + 2 individual galaxies x 2 + members 1)." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "\n", + "# --- Cluster-scale dark-matter halo: elliptical dPIE, fully free, not tied to light ---\n", + "halo_mass = af.Model(al.mp.dPIEMass)\n", + "halo_mass.centre = (0.0, 0.0) # [FIXED] often freed near the BCG (set centre_0 / centre_1 to priors)\n", + "halo_mass.ellipticity = af.UniformPrior(lower_limit=0.0, upper_limit=0.7) # [FREE]\n", + "halo_mass.angle_pos = af.UniformPrior(lower_limit=0.0, upper_limit=180.0) # [FREE] degrees\n", + "halo_mass.sigma = af.UniformPrior(lower_limit=500.0, upper_limit=1500.0) # [FREE] cluster-scale dispersion (km/s)\n", + "halo_mass.r_core = af.UniformPrior(lower_limit=20.0, upper_limit=150.0) # [FREE] large halo core (arcsec)\n", + "halo_mass.r_cut = 1000.0 # [FIXED] truncation unconstrained within the field; fixed large\n", + "halo_mass.redshift_object = redshift_lens\n", + "halo_mass.redshift_source = redshift_source\n", + "halo_mass.H0 = H0\n", + "halo_mass.Om0 = Om0\n", + "cluster_halo = af.Model(al.Galaxy, redshift=redshift_lens, mass=halo_mass)\n", + "\n", + "# --- Individually-modelled galaxies: the BCG (first) + any freed member, each its own free dPIE ---\n", + "individual_sigma_priors = [(200.0, 600.0), (150.0, 450.0)] # per-galaxy free sigma range (km/s)\n", + "individual_rcut_priors = [(20.0, 200.0), (20.0, 150.0)] # per-galaxy free r_cut range (arcsec)\n", + "\n", + "individual_galaxies = []\n", + "for centre, (sigma_lo, sigma_hi), (rcut_lo, rcut_hi) in zip(\n", + " individual_centres, individual_sigma_priors, individual_rcut_priors\n", + "):\n", + " mass = af.Model(al.mp.dPIEMassSph)\n", + " mass.centre = centre # [FIXED] observed position\n", + " mass.sigma = af.UniformPrior(lower_limit=sigma_lo, upper_limit=sigma_hi) # [FREE] this galaxy's own dispersion\n", + " mass.r_core = 0.3 # [FIXED] arcsec\n", + " mass.r_cut = af.UniformPrior(lower_limit=rcut_lo, upper_limit=rcut_hi) # [FREE] this galaxy's own truncation\n", + " mass.redshift_object = redshift_lens\n", + " mass.redshift_source = redshift_source\n", + " mass.H0 = H0\n", + " mass.Om0 = Om0\n", + " individual_galaxies.append(af.Model(al.Galaxy, redshift=redshift_lens, mass=mass))\n", + "\n", + "# --- Members: one free sigma_ref, Faber-Jackson exponent 0.25 ---\n", + "sigma_ref = af.UniformPrior(lower_limit=100.0, upper_limit=400.0) # [FREE] km/s (dispersion at mag0)\n", + "mag0 = individual_magnitudes[0] # [FIXED] reference magnitude (Lenstool's mag0) = the BCG's brightness\n", + "sigma_exponent = 0.25 # [FIXED] Faber-Jackson (make this an af.UniformPrior to free it \u2014 a one-line change)\n", + "radius_exponent = 0.5 # [FIXED] constant mass-to-light: r_core, r_cut ~ L^0.5\n", + "r_core_ref = 0.15 # [FIXED] arcsec\n", + "r_cut_ref = 20.0 # [FIXED] arcsec\n", + "\n", + "scaling_galaxies = []\n", + "for centre, magnitude in zip(scaling_centres, scaling_magnitudes):\n", + " luminosity_ratio = 10.0 ** (0.4 * (mag0 - magnitude)) # L / L_ref from magnitudes\n", + "\n", + " mass = af.Model(al.mp.dPIEMassSph)\n", + " mass.centre = centre # [FIXED]\n", + " mass.sigma = sigma_ref * luminosity_ratio**sigma_exponent # tied to the one free sigma_ref\n", + " mass.r_core = r_core_ref * luminosity_ratio**radius_exponent\n", + " mass.r_cut = r_cut_ref * luminosity_ratio**radius_exponent\n", + " mass.redshift_object = redshift_lens\n", + " mass.redshift_source = redshift_source\n", + " mass.H0 = H0\n", + " mass.Om0 = Om0\n", + " scaling_galaxies.append(af.Model(al.Galaxy, redshift=redshift_lens, mass=mass))\n", + "\n", + "model_1 = af.Collection(\n", + " cluster_halo=cluster_halo,\n", + " individual_galaxies=af.Collection(individual_galaxies),\n", + " scaling_galaxies=af.Collection(scaling_galaxies),\n", + ")\n", + "\n", + "print(\"=\" * 80)\n", + "print(\"Model 1 \u2014 The Standard Lenstool Model\")\n", + "print(\"=\" * 80)\n", + "print(f\"Total free parameters: {model_1.prior_count} (halo 4 + 2 individual x 2 + members 1 = 9)\")\n", + "assert model_1.prior_count == 9\n", + "print(model_1.info)\n" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Model 2 \u2014 The Angular Parameterization: dPIEMassB0__\n", + "\n", + "**How.** Swap every dPIE for its ``B0`` twin and compose in angular parameters: ``ra`` (core),\n", + "``rs`` (truncation) and ``b0`` (lens strength) \u2014 all arcsec, with **no redshift and no cosmology**. The\n", + "radii map one-to-one (``r_core -> ra``, ``r_cut -> rs``); the mass normalization moves ``sigma``\n", + "(km/s) ``-> b0`` (arcsec). The scaling relation moves with it, and the exponent **doubles**:\n", + "\n", + " Model 1 (sigma): sigma_i = sigma_ref * (L/L_ref)^0.25\n", + " Model 2 (b0): b0_i = b0_ref * (L/L_ref)^0.50 # because b0 ~ sigma^2\n", + "\n", + "It is the *same* Faber-Jackson relation, written in the quantity the lensing actually constrains. The\n", + "halo becomes an elliptical ``dPIEMassB0`` (ellipticity via ``ell_comps``); individual galaxies and\n", + "members become ``dPIEMassB0Sph``.\n", + "\n", + "**Why.** The ``sigma`` parameterization needs the redshifts (``b0 = 6*648000*(sigma/c)^2*(D_LS/D_S)``),\n", + "so ``sigma`` only means something once ``z_lens`` / ``z_source`` are known. When they are *not*, the\n", + "images constrain only the angular ``b0``, and ``sigma`` is degenerate with ``D_LS/D_S``. Fitting ``b0``\n", + "directly fits the quantity the data determines, in arcsec you can reason about. Convert back once the\n", + "redshifts are measured: ``sigma = c * sqrt(b0 / K)``, ``K = 6*648000*(D_LS/D_S)``. The ``B0`` mass\n", + "carries no redshift; the galaxy ``redshift`` is only a ray-tracing plane label.\n", + "\n", + "**Free vs fixed.** Identical structure to Model 1 (9 free): halo free ``ell_comps`` (2) + ``ra`` + ``b0``;\n", + "individual galaxies free ``b0`` + ``rs``; members one free ``b0_ref``." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "\n", + "# --- Halo: elliptical dPIEMassB0 (ellipticity via ell_comps) ---\n", + "halo_mass = af.Model(al.mp.dPIEMassB0)\n", + "halo_mass.centre = (0.0, 0.0) # [FIXED]\n", + "halo_mass.ell_comps.ell_comps_0 = af.UniformPrior(lower_limit=-0.5, upper_limit=0.5) # [FREE]\n", + "halo_mass.ell_comps.ell_comps_1 = af.UniformPrior(lower_limit=-0.5, upper_limit=0.5) # [FREE]\n", + "halo_mass.ra = af.UniformPrior(lower_limit=20.0, upper_limit=150.0) # [FREE] core (arcsec)\n", + "halo_mass.rs = 1000.0 # [FIXED] truncation, fixed large\n", + "halo_mass.b0 = af.UniformPrior(lower_limit=5.0, upper_limit=50.0) # [FREE] cluster-scale lens strength (arcsec)\n", + "cluster_halo = af.Model(al.Galaxy, redshift=redshift_lens, mass=halo_mass)\n", + "\n", + "# --- Individually-modelled galaxies: free b0 + free rs (angular, redshift-free) ---\n", + "individual_b0_priors = [(0.0, 5.0), (0.0, 3.0)] # arcsec, angular lens-strength range per galaxy\n", + "individual_rs_priors = [(20.0, 200.0), (20.0, 150.0)] # arcsec, truncation range per galaxy\n", + "\n", + "individual_galaxies_b0 = []\n", + "for centre, (b0_lo, b0_hi), (rs_lo, rs_hi) in zip(\n", + " individual_centres, individual_b0_priors, individual_rs_priors\n", + "):\n", + " mass = af.Model(al.mp.dPIEMassB0Sph)\n", + " mass.centre = centre # [FIXED]\n", + " mass.b0 = af.UniformPrior(lower_limit=b0_lo, upper_limit=b0_hi) # [FREE] angular lens strength (arcsec)\n", + " mass.ra = 0.3 # [FIXED] = Model 1's r_core\n", + " mass.rs = af.UniformPrior(lower_limit=rs_lo, upper_limit=rs_hi) # [FREE] truncation (arcsec)\n", + " individual_galaxies_b0.append(af.Model(al.Galaxy, redshift=redshift_lens, mass=mass))\n", + "\n", + "# --- Members: one free b0_ref, exponent 0.5 (b0 ~ L^0.5) ---\n", + "b0_ref = af.UniformPrior(lower_limit=0.0, upper_limit=1.0) # [FREE] arcsec (lens strength at mag0)\n", + "b0_exponent = 0.5 # [FIXED] b0 ~ L^0.5 (Model 1 used sigma ~ L^0.25; b0 ~ sigma^2 doubles the exponent)\n", + "mag0 = individual_magnitudes[0] # [FIXED]\n", + "ra_ref = 0.15 # [FIXED] = Model 1's r_core_ref\n", + "rs_ref = 20.0 # [FIXED] = Model 1's r_cut_ref\n", + "\n", + "scaling_galaxies_b0 = []\n", + "for centre, magnitude in zip(scaling_centres, scaling_magnitudes):\n", + " luminosity_ratio = 10.0 ** (0.4 * (mag0 - magnitude))\n", + "\n", + " mass = af.Model(al.mp.dPIEMassB0Sph)\n", + " mass.centre = centre # [FIXED]\n", + " mass.b0 = b0_ref * luminosity_ratio**b0_exponent # tied to the one free b0_ref\n", + " mass.ra = ra_ref * luminosity_ratio**0.5\n", + " mass.rs = rs_ref * luminosity_ratio**0.5\n", + " scaling_galaxies_b0.append(af.Model(al.Galaxy, redshift=redshift_lens, mass=mass))\n", + "\n", + "model_2 = af.Collection(\n", + " cluster_halo=cluster_halo,\n", + " individual_galaxies=af.Collection(individual_galaxies_b0),\n", + " scaling_galaxies=af.Collection(scaling_galaxies_b0),\n", + ")\n", + "\n", + "print(\"\\n\" + \"=\" * 80)\n", + "print(\"Model 2 \u2014 The Angular (b0) Parameterization\")\n", + "print(\"=\" * 80)\n", + "print(f\"Total free parameters: {model_2.prior_count} (same structure as Model 1)\")\n", + "assert model_2.prior_count == 9" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Model 2 \u2014 Numerical Checks__\n", + "\n", + "Model 2 is Model 1 re-expressed, so it must reproduce Model 1's physics once the redshifts are\n", + "supplied. We relate the two with the distance factor ``K`` (redshifts + cosmology used *only* for the\n", + "comparison, never inside the b0 model), and check three things:\n", + "\n", + " (a) the Faber-Jackson relation is consistent \u2014 a member's ``b0`` from the ``b0``-relation equals\n", + " ``K (sigma/c)^2`` with ``sigma`` from the ``sigma``-relation (the exponent doubling is exact);\n", + " (b) the deflections are identical, profile-for-profile;\n", + " (c) a fitted ``b0_ref`` converts back to the Model 1 ``sigma_ref`` exactly." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "C_KM_S = 299792.458\n", + "cosmology = ag.cosmo.FlatLambdaCDM(H0=H0, Om0=Om0) # matches dPIEMass's internal cosmology\n", + "d_s = cosmology.angular_diameter_distance_to_earth_in_kpc_from(redshift=redshift_source)\n", + "d_ls = cosmology.angular_diameter_distance_between_redshifts_in_kpc_from(\n", + " redshift_0=redshift_lens, redshift_1=redshift_source\n", + ")\n", + "K = 6.0 * 648000.0 * (d_ls / d_s) # b0 = K * (sigma / c)^2\n", + "\n", + "grid = al.Grid2DIrregular([[3.0, 4.0], [-8.0, 2.0], [10.0, -6.0]])\n", + "\n", + "# A worked value: the b0_ref a redshift-free fit would return for sigma_ref = 250 km/s.\n", + "sigma_ref_value = 250.0\n", + "b0_ref_value = K * (sigma_ref_value / C_KM_S) ** 2\n", + "\n", + "# (a) Faber-Jackson consistency.\n", + "worst_relation = 0.0\n", + "for magnitude in scaling_magnitudes:\n", + " ratio = 10.0 ** (0.4 * (mag0 - magnitude))\n", + " worst_relation = max(\n", + " worst_relation,\n", + " abs(b0_ref_value * ratio**0.5 - K * (sigma_ref_value * ratio**0.25 / C_KM_S) ** 2),\n", + " )\n", + "print(f\"(a) member b0 [b0-relation vs sigma-relation]: max diff = {worst_relation:.2e}\")\n", + "assert worst_relation < 1e-12\n", + "\n", + "# (b) Deflections identical: dPIEMassB0Sph(b0) == dPIEMassSph(sigma) for each member.\n", + "worst_defl = 0.0\n", + "for centre, magnitude in zip(scaling_centres, scaling_magnitudes):\n", + " ratio = 10.0 ** (0.4 * (mag0 - magnitude))\n", + " mass_b0 = al.mp.dPIEMassB0Sph(\n", + " centre=centre, ra=ra_ref * ratio**0.5, rs=rs_ref * ratio**0.5, b0=b0_ref_value * ratio**0.5\n", + " )\n", + " mass_sigma = al.mp.dPIEMassSph(\n", + " centre=centre,\n", + " sigma=sigma_ref_value * ratio**0.25,\n", + " r_core=ra_ref * ratio**0.5,\n", + " r_cut=rs_ref * ratio**0.5,\n", + " redshift_object=redshift_lens,\n", + " redshift_source=redshift_source,\n", + " H0=H0,\n", + " Om0=Om0,\n", + " )\n", + " worst_defl = max(\n", + " worst_defl,\n", + " np.max(\n", + " np.abs(\n", + " np.asarray(mass_b0.deflections_yx_2d_from(grid=grid))\n", + " - np.asarray(mass_sigma.deflections_yx_2d_from(grid=grid))\n", + " )\n", + " ),\n", + " )\n", + "print(f\"(b) member deflections [b0 vs sigma param]: max diff = {worst_defl:.2e}\")\n", + "assert worst_defl < 1e-12\n", + "\n", + "# (c) Round-trip.\n", + "sigma_ref_recovered = C_KM_S * (b0_ref_value / K) ** 0.5\n", + "print(f\"(c) b0_ref = {b0_ref_value:.5f} arcsec -> sigma_ref = {sigma_ref_recovered:.2f} km/s (input {sigma_ref_value})\")\n", + "assert np.isclose(sigma_ref_recovered, sigma_ref_value, rtol=1e-9)\n", + "print(\"Model 2 checks passed.\")\n" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Model 3 \u2014 Mass Anchoring: the Members Normalized to the BCG__\n", + "\n", + "**How.** Take Model 1 and *delete* the members' free normalization ``sigma_ref``. In its place the\n", + "members are normalized to the BCG's OWN velocity dispersion \u2014 the same free ``sigma`` that sets the\n", + "BCG's mass:\n", + "\n", + " Model 1: sigma_i = sigma_ref * (L/L_BCG)^0.25 # sigma_ref FREE, independent of the BCG\n", + " Model 3: sigma_i = sigma_BCG * (L/L_BCG)^0.25 # sigma_BCG = the BCG's own free mass\n", + "\n", + "Because ``sigma_BCG`` is already free (it *is* the BCG's mass), the member population now adds **zero**\n", + "free parameters. This is the referee's \"anchor alpha to the BCG Einstein radius\" \u2014 ``sigma_BCG`` and\n", + "the BCG Einstein radius are the same quantity up to the fixed conversion.\n", + "\n", + "**Why.** The most economical model: one free normalization (the BCG's, tightly pinned by the central\n", + "images) describes everyone, with an interpretable prior. You save a parameter.\n", + "\n", + "**The cost \u2014 and why this is NOT the default (Model 1).** It *assumes the BCG lies on the member\n", + "relation*; real BCGs frequently deviate (merger products at the bottom of the potential well), which is\n", + "why high-precision models free the BCG separately (Model 1). And it hard-couples every member's mass to\n", + "the BCG's \u2014 move the fitted BCG mass and all members move in fixed proportion. A real, used variant\n", + "(simpler models, the referee's convention), but a deliberate coupling choice, not the standard.\n", + "\n", + "**Free vs fixed.** Halo as Model 1 (4). BCG: free ``sigma`` (its own mass *and* the member anchor) +\n", + "``r_cut``. Other freed galaxy: its own free ``sigma`` + ``r_cut``. Members: **zero** free. Total: 8." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "\n", + "# --- Halo: identical to Model 1 (the halo is unaffected by how the members are anchored) ---\n", + "halo_mass = af.Model(al.mp.dPIEMass)\n", + "halo_mass.centre = (0.0, 0.0) # [FIXED]\n", + "halo_mass.ellipticity = af.UniformPrior(lower_limit=0.0, upper_limit=0.7) # [FREE]\n", + "halo_mass.angle_pos = af.UniformPrior(lower_limit=0.0, upper_limit=180.0) # [FREE]\n", + "halo_mass.sigma = af.UniformPrior(lower_limit=500.0, upper_limit=1500.0) # [FREE]\n", + "halo_mass.r_core = af.UniformPrior(lower_limit=20.0, upper_limit=150.0) # [FREE]\n", + "halo_mass.r_cut = 1000.0 # [FIXED]\n", + "halo_mass.redshift_object = redshift_lens\n", + "halo_mass.redshift_source = redshift_source\n", + "halo_mass.H0 = H0\n", + "halo_mass.Om0 = Om0\n", + "cluster_halo = af.Model(al.Galaxy, redshift=redshift_lens, mass=halo_mass)\n", + "\n", + "# --- The BCG: the ANCHOR. Its free sigma sets its own mass AND normalizes every member. ---\n", + "bcg_mass = af.Model(al.mp.dPIEMassSph)\n", + "bcg_mass.centre = individual_centres[0] # [FIXED]\n", + "bcg_mass.sigma = af.UniformPrior(lower_limit=200.0, upper_limit=600.0) # [FREE] BCG mass + member anchor\n", + "bcg_mass.r_core = 0.3 # [FIXED]\n", + "bcg_mass.r_cut = af.UniformPrior(lower_limit=20.0, upper_limit=200.0) # [FREE] BCG's own truncation\n", + "bcg_mass.redshift_object = redshift_lens\n", + "bcg_mass.redshift_source = redshift_source\n", + "bcg_mass.H0 = H0\n", + "bcg_mass.Om0 = Om0\n", + "bcg = af.Model(al.Galaxy, redshift=redshift_lens, mass=bcg_mass)\n", + "\n", + "# --- Other individually-freed galaxies (still independent \u2014 e.g. a member on an arc) ---\n", + "individual_sigma_priors = [(200.0, 600.0), (150.0, 450.0)]\n", + "individual_rcut_priors = [(20.0, 200.0), (20.0, 150.0)]\n", + "\n", + "other_individual_galaxies = []\n", + "for centre, (sigma_lo, sigma_hi), (rcut_lo, rcut_hi) in zip(\n", + " individual_centres[1:], individual_sigma_priors[1:], individual_rcut_priors[1:]\n", + "):\n", + " mass = af.Model(al.mp.dPIEMassSph)\n", + " mass.centre = centre\n", + " mass.sigma = af.UniformPrior(lower_limit=sigma_lo, upper_limit=sigma_hi) # [FREE] independent of the BCG\n", + " mass.r_core = 0.3\n", + " mass.r_cut = af.UniformPrior(lower_limit=rcut_lo, upper_limit=rcut_hi) # [FREE]\n", + " mass.redshift_object = redshift_lens\n", + " mass.redshift_source = redshift_source\n", + " mass.H0 = H0\n", + " mass.Om0 = Om0\n", + " other_individual_galaxies.append(af.Model(al.Galaxy, redshift=redshift_lens, mass=mass))\n", + "\n", + "# --- Members: normalized to the BCG's OWN sigma. ZERO new free parameters. ---\n", + "mag0 = individual_magnitudes[0] # [FIXED] = BCG magnitude\n", + "r_core_ref = 0.15 # [FIXED]\n", + "r_cut_ref = 20.0 # [FIXED]\n", + "\n", + "mass_anchored_members = []\n", + "for centre, magnitude in zip(scaling_centres, scaling_magnitudes):\n", + " luminosity_ratio = 10.0 ** (0.4 * (mag0 - magnitude)) # L / L_BCG\n", + "\n", + " mass = af.Model(al.mp.dPIEMassSph)\n", + " mass.centre = centre\n", + " mass.sigma = bcg_mass.sigma * luminosity_ratio**0.25 # tied to the BCG's free sigma -> NO new parameter\n", + " mass.r_core = r_core_ref * luminosity_ratio**0.5\n", + " mass.r_cut = r_cut_ref * luminosity_ratio**0.5\n", + " mass.redshift_object = redshift_lens\n", + " mass.redshift_source = redshift_source\n", + " mass.H0 = H0\n", + " mass.Om0 = Om0\n", + " mass_anchored_members.append(af.Model(al.Galaxy, redshift=redshift_lens, mass=mass))\n", + "\n", + "model_3 = af.Collection(\n", + " cluster_halo=cluster_halo,\n", + " individual_galaxies=af.Collection([bcg] + other_individual_galaxies),\n", + " scaling_galaxies=af.Collection(mass_anchored_members),\n", + ")\n", + "\n", + "print(\"\\n\" + \"=\" * 80)\n", + "print(\"Model 3 \u2014 Mass Anchoring (members normalized to the BCG)\")\n", + "print(\"=\" * 80)\n", + "print(f\"Total free parameters: {model_3.prior_count} (Model 1 was 9; the member tier now adds 0)\")\n", + "assert model_3.prior_count == 8" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Model 3 \u2014 Numerical Checks__\n", + "\n", + " (a) the member collection adds ZERO free parameters (all 8 belong to the halo + individual galaxies);\n", + " (b) at any instance, each member's sigma is exactly ``sigma_BCG * (L/L_BCG)^0.25`` \u2014 a fixed fraction\n", + " of the BCG, so a change in the fitted BCG mass rescales them all together." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "# (a) The member tier adds zero free parameters.\n", + "without_members = af.Collection(\n", + " cluster_halo=cluster_halo,\n", + " individual_galaxies=af.Collection([bcg] + other_individual_galaxies),\n", + ")\n", + "added = model_3.prior_count - without_members.prior_count\n", + "print(f\"(a) free params: without members = {without_members.prior_count}, with members = {model_3.prior_count} (members added {added})\")\n", + "assert added == 0\n", + "\n", + "# (b) The members are driven by the BCG's sigma.\n", + "instance = model_3.instance_from_prior_medians()\n", + "sigma_bcg_instance = instance.individual_galaxies[0].mass.sigma\n", + "print(f\"(b) instance BCG sigma = {sigma_bcg_instance:.2f} km/s; members = sigma_BCG * (L/L_BCG)^0.25:\")\n", + "worst_coupling = 0.0\n", + "for i, magnitude in enumerate(scaling_magnitudes):\n", + " ratio = 10.0 ** (0.4 * (mag0 - magnitude))\n", + " expected = sigma_bcg_instance * ratio**0.25\n", + " actual = instance.scaling_galaxies[i].mass.sigma\n", + " worst_coupling = max(worst_coupling, abs(actual - expected))\n", + " print(f\" member {i} (mag {magnitude}): sigma = {actual:6.2f} km/s = {ratio**0.25:.3f} x sigma_BCG\")\n", + "assert worst_coupling < 1e-9\n", + "print(\"Model 3 checks passed.\")\n" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Model 4 \u2014 The Two-Parameter Relation__\n", + "\n", + "**How.** Model 1 fixes the truncation reference ``r_cut_ref`` and frees only ``sigma_ref``. Here we\n", + "free ``r_cut_ref`` as well, so the member population has **two** shared free parameters \u2014 a dispersion\n", + "normalization and a truncation normalization:\n", + "\n", + " sigma_i = sigma_ref * (L/L_ref)^0.25 # sigma_ref FREE (as Model 1)\n", + " r_cut_i = r_cut_ref * (L/L_ref)^0.5 # r_cut_ref now FREE too\n", + "\n", + "Every member's truncation now scales from a *free* reference instead of a fixed one. It is a one-symbol\n", + "change (``r_cut_ref`` becomes a prior) \u2014 but it is a real modelling decision, so it earns a section.\n", + "\n", + "**Why.** The reference truncation ``r_cut*`` is genuinely uncertain \u2014 it sets how much mass each member\n", + "carries at large radius, and the data can constrain it. Lenstool's ``potfile`` routinely optimizes both\n", + "``sigma*`` and ``r_cut*``. You free it when the members contribute enough lensing that their outer mass\n", + "matters; you fix it (Model 1) when they do not, to save a parameter and avoid a weakly-constrained\n", + "direction.\n", + "\n", + "**Free vs fixed.** As Model 1, but the member tier now has 2 free parameters (``sigma_ref`` +\n", + "``r_cut_ref``) instead of 1. Total: 10." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "\n", + "# --- Halo: as Model 1 ---\n", + "halo_mass = af.Model(al.mp.dPIEMass)\n", + "halo_mass.centre = (0.0, 0.0) # [FIXED]\n", + "halo_mass.ellipticity = af.UniformPrior(lower_limit=0.0, upper_limit=0.7) # [FREE]\n", + "halo_mass.angle_pos = af.UniformPrior(lower_limit=0.0, upper_limit=180.0) # [FREE]\n", + "halo_mass.sigma = af.UniformPrior(lower_limit=500.0, upper_limit=1500.0) # [FREE]\n", + "halo_mass.r_core = af.UniformPrior(lower_limit=20.0, upper_limit=150.0) # [FREE]\n", + "halo_mass.r_cut = 1000.0 # [FIXED]\n", + "halo_mass.redshift_object = redshift_lens\n", + "halo_mass.redshift_source = redshift_source\n", + "halo_mass.H0 = H0\n", + "halo_mass.Om0 = Om0\n", + "cluster_halo = af.Model(al.Galaxy, redshift=redshift_lens, mass=halo_mass)\n", + "\n", + "# --- Individually-modelled galaxies: as Model 1 ---\n", + "individual_sigma_priors = [(200.0, 600.0), (150.0, 450.0)]\n", + "individual_rcut_priors = [(20.0, 200.0), (20.0, 150.0)]\n", + "\n", + "individual_galaxies = []\n", + "for centre, (sigma_lo, sigma_hi), (rcut_lo, rcut_hi) in zip(\n", + " individual_centres, individual_sigma_priors, individual_rcut_priors\n", + "):\n", + " mass = af.Model(al.mp.dPIEMassSph)\n", + " mass.centre = centre\n", + " mass.sigma = af.UniformPrior(lower_limit=sigma_lo, upper_limit=sigma_hi) # [FREE]\n", + " mass.r_core = 0.3\n", + " mass.r_cut = af.UniformPrior(lower_limit=rcut_lo, upper_limit=rcut_hi) # [FREE]\n", + " mass.redshift_object = redshift_lens\n", + " mass.redshift_source = redshift_source\n", + " mass.H0 = H0\n", + " mass.Om0 = Om0\n", + " individual_galaxies.append(af.Model(al.Galaxy, redshift=redshift_lens, mass=mass))\n", + "\n", + "# --- Members: TWO free normalizations, sigma_ref and r_cut_ref ---\n", + "sigma_ref = af.UniformPrior(lower_limit=100.0, upper_limit=400.0) # [FREE] km/s\n", + "r_cut_ref = af.UniformPrior(lower_limit=5.0, upper_limit=50.0) # [FREE] arcsec <-- now a prior, not a constant\n", + "mag0 = individual_magnitudes[0] # [FIXED]\n", + "sigma_exponent = 0.25 # [FIXED]\n", + "radius_exponent = 0.5 # [FIXED]\n", + "r_core_ref = 0.15 # [FIXED] (the core reference is usually kept fixed even in the 2-parameter relation)\n", + "\n", + "scaling_galaxies = []\n", + "for centre, magnitude in zip(scaling_centres, scaling_magnitudes):\n", + " luminosity_ratio = 10.0 ** (0.4 * (mag0 - magnitude))\n", + "\n", + " mass = af.Model(al.mp.dPIEMassSph)\n", + " mass.centre = centre\n", + " mass.sigma = sigma_ref * luminosity_ratio**sigma_exponent # tied to free sigma_ref\n", + " mass.r_core = r_core_ref * luminosity_ratio**radius_exponent # fixed reference\n", + " mass.r_cut = r_cut_ref * luminosity_ratio**radius_exponent # tied to free r_cut_ref\n", + " mass.redshift_object = redshift_lens\n", + " mass.redshift_source = redshift_source\n", + " mass.H0 = H0\n", + " mass.Om0 = Om0\n", + " scaling_galaxies.append(af.Model(al.Galaxy, redshift=redshift_lens, mass=mass))\n", + "\n", + "model_4 = af.Collection(\n", + " cluster_halo=cluster_halo,\n", + " individual_galaxies=af.Collection(individual_galaxies),\n", + " scaling_galaxies=af.Collection(scaling_galaxies),\n", + ")\n", + "\n", + "print(\"\\n\" + \"=\" * 80)\n", + "print(\"Model 4 \u2014 The Two-Parameter Relation\")\n", + "print(\"=\" * 80)\n", + "print(f\"Total free parameters: {model_4.prior_count} (Model 1 was 9; the member tier is now 2, not 1)\")\n", + "assert model_4.prior_count == 10\n", + "\n", + "# The member tier now contributes 2 free parameters (sigma_ref + r_cut_ref) instead of 1.\n", + "without_members = af.Collection(\n", + " cluster_halo=cluster_halo, individual_galaxies=af.Collection(individual_galaxies)\n", + ")\n", + "added = model_4.prior_count - without_members.prior_count\n", + "print(f\"member tier free parameters: {added} (sigma_ref + r_cut_ref)\")\n", + "assert added == 2\n", + "print(\"Model 4 checks passed.\")\n" + ], + "outputs": [], + "execution_count": null + } + ], + "metadata": { + "anaconda-cloud": {}, + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.6.1" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} \ No newline at end of file diff --git a/notebooks/cluster/mass_parameterizations_pyautolens.ipynb b/notebooks/cluster/mass_parameterizations_pyautolens.ipynb new file mode 100644 index 000000000..84ca58735 --- /dev/null +++ b/notebooks/cluster/mass_parameterizations_pyautolens.ipynb @@ -0,0 +1,509 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Mass Parameterizations II \u2014 Mapping Lenstool onto PyAutoLens\n", + "============================================================\n", + "\n", + "**For Lenstool users: how and why PyAutoLens re-expresses the standard Lenstool cluster / group model\n", + "in its native parameterization for multi-galaxy lenses (MGLs), and how to keep the dPIE's truncation on\n", + "the outskirts if you want it.**\n", + "\n", + "The companion script ``cluster/mass_parameterizations.py`` builds the standard Lenstool model. This\n", + "script takes its **Model 1** as the starting point and maps it, component by component, onto what\n", + "PyAutoLens uses at galaxy / group scale. Both are supported; this is an *alternative*, not a\n", + "replacement, and the focus is the **mass**.\n", + "\n", + "Three sections:\n", + "\n", + " - **Section 1 \u2014 Lenstool Model 1** (dPIE throughout, velocity dispersions, a magnitude-anchored\n", + " Faber-Jackson relation), the starting point.\n", + " - **Section 2 \u2014 The standard PyAutoLens model**, mapping every component and justifying each change.\n", + " - **Section 3 \u2014 Keeping truncation**: dPIE for the scaling galaxies while the main and extra galaxies\n", + " stay Isothermal, and how the Faber-Jackson linking bridges the two parameterizations.\n", + "\n", + "Every model has four tiers, and this is the structure of a PyAutoLens MGL:\n", + "\n", + " Tier Role Coupling to the BGC\n", + " ---------------- -------------------------------------- --------------------------------\n", + " cluster halo the smooth large-scale mass none (its own free profile)\n", + " main galaxies the primary lens galaxies (incl. BGC) none \u2014 free\n", + " extra galaxies bright members worth freeing free, but luminosity-BOUNDED\n", + " scaling galaxies the faint outskirts (\"potfile\") TIED \u2014 mass set by the relation\n", + "\n", + "The luminosity that drives the relation is computed from each galaxy's fitted MGE light,\n", + "``L = sum_g 2*pi*sigma_g^2/q_g * intensity_g`` (the ``2*np.pi*g.sigma**2/g.axis_ratio()*g.intensity``\n", + "you have seen); only *ratios* matter, so we drive everything off illustrative magnitudes here\n", + "(``L/L_BGC = 10 ** (0.4 * (mag0 - m))``).\n", + "\n", + "Everything uses the model-composition API and prints ``model.info`` so free vs fixed is explicit.\n", + "\n", + "__What changes, at a glance__\n", + "\n", + " Component Lenstool (Section 1) PyAutoLens (Sections 2-3)\n", + " ---------------- ---------------------------- --------------------------------------\n", + " cluster halo elliptical dPIE (sigma) NFW (the CDM halo profile)\n", + " main galaxies dPIE (sigma, r_core, r_cut) Isothermal (einstein_radius)\n", + " extra galaxies dPIE, free Isothermal, free but luminosity-bounded\n", + " scaling galaxies dPIE, free sigma_ref Isothermal, TIED to the BGC's Einstein radius\n", + " mass parameter sigma (needs redshifts) einstein_radius (arcsec; redshift-free)\n", + " luminosity input magnitude computed from the fitted MGE light" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Google Colab Setup__\n", + "\n", + "This cell sets up the environment when the notebook is run on Google Colab: it installs the\n", + "required PyAuto packages, clones the workspace (configuration files and example datasets) and\n", + "points the configuration at it. If you are running the notebook elsewhere (e.g. locally via\n", + "your own installation) it does nothing, and you can run it safely.\n", + "\n", + "Colab tip: model-fits run much faster on a GPU \u2014 enable one via \"Runtime\" -> \"Change runtime\n", + "type\" -> \"Hardware accelerator\" before running the notebook." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "try:\n", + " import google.colab\n", + " import subprocess\n", + " import sys\n", + "\n", + " subprocess.check_call(\n", + " [sys.executable, \"-m\", \"pip\", \"install\", \"autonerves\", \"--no-deps\"]\n", + " )\n", + "except ImportError:\n", + " pass\n", + "\n", + "from autonerves import setup_colab\n", + "\n", + "setup_colab.setup(\"autolens\")" + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "\n", + "import numpy as np\n", + "\n", + "import autofit as af\n", + "import autolens as al\n", + "import autogalaxy as ag" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Setup__\n", + "\n", + "Redshifts / cosmology are needed by the *Lenstool* (sigma) side and again to *check* the\n", + "sigma <-> einstein_radius bridge \u2014 they do not appear in the PyAutoLens mass model." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "redshift_lens = 0.5\n", + "redshift_source = 2.0\n", + "H0 = 67.66\n", + "Om0 = 0.30966\n", + "\n", + "main_centres = [(0.0, 0.0), (12.0, 8.0)] # primary lens galaxies (the brightest is the BGC)\n", + "main_magnitudes = [17.8, 18.9]\n", + "\n", + "extra_centres = [(8.5, 5.5)] # a bright member worth freeing (e.g. near a multiple image)\n", + "extra_magnitudes = [18.5]\n", + "\n", + "scaling_centres = [(5.5, -6.5), (-7.5, 3.0), (3.0, 13.0)] # the faint outskirts\n", + "scaling_magnitudes = [19.2, 20.0, 21.0]\n", + "\n", + "bgc_index = int(np.argmin(main_magnitudes)) # brightest (smallest magnitude) main galaxy\n", + "mag0 = main_magnitudes[bgc_index] # the BGC magnitude anchors the luminosity ratios\n", + "\n", + "\n", + "def luminosity_ratio(magnitude):\n", + " \"\"\"L / L_BGC from a magnitude, relative to the BGC (mag0).\"\"\"\n", + " return 10.0 ** (0.4 * (mag0 - magnitude))\n" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "================================================================================\n", + "__Section 1 \u2014 Lenstool Model 1__\n", + "================================================================================\n", + "\n", + "The full standard Lenstool model, dPIE throughout (see ``cluster/mass_parameterizations.py`` Model 1).\n", + "The main and extra galaxies are both individually-modelled dPIEs (free ``sigma`` + ``r_cut``); the\n", + "distinction between them is a PyAutoLens refinement (Section 2) \u2014 in Lenstool they are both just\n", + "\"potentiel\" sections. The scaling galaxies use a magnitude-anchored Faber-Jackson relation with a free\n", + "normalization ``sigma_ref``: ``sigma_i = sigma_ref * (L/L_ref) ** 0.25``.\n", + "\n", + "Total: 11 free (halo 4 + 2 main x 2 + 1 extra x 2 + scaling 1)." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "\n", + "\n", + "def dpie_galaxy(centre, sigma_lo, sigma_hi, rcut_lo, rcut_hi):\n", + " \"\"\"An individually-modelled Lenstool galaxy: free dPIE dispersion + truncation.\"\"\"\n", + " mass = af.Model(al.mp.dPIEMassSph)\n", + " mass.centre = centre\n", + " mass.sigma = af.UniformPrior(lower_limit=sigma_lo, upper_limit=sigma_hi) # [FREE] km/s\n", + " mass.r_core = 0.3\n", + " mass.r_cut = af.UniformPrior(lower_limit=rcut_lo, upper_limit=rcut_hi) # [FREE] arcsec\n", + " mass.redshift_object = redshift_lens\n", + " mass.redshift_source = redshift_source\n", + " mass.H0 = H0\n", + " mass.Om0 = Om0\n", + " return af.Model(al.Galaxy, redshift=redshift_lens, mass=mass)\n", + "\n", + "\n", + "# Cluster halo \u2014 elliptical dPIE, fully free.\n", + "halo_mass = af.Model(al.mp.dPIEMass)\n", + "halo_mass.centre = (0.0, 0.0) # [FIXED]\n", + "halo_mass.ellipticity = af.UniformPrior(lower_limit=0.0, upper_limit=0.7) # [FREE]\n", + "halo_mass.angle_pos = af.UniformPrior(lower_limit=0.0, upper_limit=180.0) # [FREE]\n", + "halo_mass.sigma = af.UniformPrior(lower_limit=500.0, upper_limit=1500.0) # [FREE] km/s\n", + "halo_mass.r_core = af.UniformPrior(lower_limit=20.0, upper_limit=150.0) # [FREE] arcsec\n", + "halo_mass.r_cut = 1000.0 # [FIXED]\n", + "halo_mass.redshift_object = redshift_lens\n", + "halo_mass.redshift_source = redshift_source\n", + "halo_mass.H0 = H0\n", + "halo_mass.Om0 = Om0\n", + "cluster_halo_lt = af.Model(al.Galaxy, redshift=redshift_lens, mass=halo_mass)\n", + "\n", + "main_galaxies_lt = [dpie_galaxy(c, 100.0, 600.0, 20.0, 200.0) for c in main_centres]\n", + "extra_galaxies_lt = [dpie_galaxy(c, 100.0, 450.0, 20.0, 150.0) for c in extra_centres]\n", + "\n", + "# Scaling galaxies \u2014 Faber-Jackson on a FREE normalization sigma_ref (independent of the BGC).\n", + "sigma_ref = af.UniformPrior(lower_limit=100.0, upper_limit=400.0) # [FREE] km/s\n", + "r_core_ref, r_cut_ref = 0.15, 20.0 # [FIXED] arcsec\n", + "scaling_galaxies_lt = []\n", + "for centre, magnitude in zip(scaling_centres, scaling_magnitudes):\n", + " ratio = luminosity_ratio(magnitude)\n", + " mass = af.Model(al.mp.dPIEMassSph)\n", + " mass.centre = centre # [FIXED]\n", + " mass.sigma = sigma_ref * ratio**0.25 # tied to the one free sigma_ref\n", + " mass.r_core = r_core_ref * ratio**0.5\n", + " mass.r_cut = r_cut_ref * ratio**0.5\n", + " mass.redshift_object = redshift_lens\n", + " mass.redshift_source = redshift_source\n", + " mass.H0 = H0\n", + " mass.Om0 = Om0\n", + " scaling_galaxies_lt.append(af.Model(al.Galaxy, redshift=redshift_lens, mass=mass))\n", + "\n", + "model_lenstool = af.Collection(\n", + " cluster_halo=cluster_halo_lt,\n", + " main_galaxies=af.Collection(main_galaxies_lt),\n", + " extra_galaxies=af.Collection(extra_galaxies_lt),\n", + " scaling_galaxies=af.Collection(scaling_galaxies_lt),\n", + ")\n", + "\n", + "print(\"=\" * 80)\n", + "print(\"Section 1 \u2014 Lenstool Model 1 (dPIE throughout)\")\n", + "print(\"=\" * 80)\n", + "print(f\"Free parameters: {model_lenstool.prior_count} (halo 4 + 2 main x 2 + 1 extra x 2 + scaling 1 = 11)\")\n", + "assert model_lenstool.prior_count == 11\n" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "================================================================================\n", + "__Section 2 \u2014 The Standard PyAutoLens Model__\n", + "================================================================================\n", + "\n", + "We map Model 1 onto the PyAutoLens MGL parameterization. Component by component:\n", + "\n", + "**(1) The halo: dPIE -> NFW.** Lenstool models the smooth halo as a dPIE for analytic convenience;\n", + "PyAutoLens uses the **NFW**, the CDM-motivated halo profile. A genuine *profile change*, not a\n", + "reparameterization \u2014 you refit ``kappa_s`` / ``scale_radius`` rather than convert. (Keep a dPIE halo if\n", + "you want to reproduce a Lenstool halo; both are supported.)\n", + "\n", + "**(2) Galaxies: dPIE -> Isothermal (SIE).** One mass parameter, ``einstein_radius`` \u2014 the reduced\n", + "deflection, in arcsec, redshift-free, exactly what the images constrain. It equals the dPIE ``b0`` in\n", + "the SIS limit (checked below). No core, no truncation.\n", + "\n", + "**(3) The coupling spectrum.** The individually-modelled galaxies split into two tiers, and the scaling\n", + "tier is anchored to the brightest galaxy (the BGC):\n", + "\n", + " - **main galaxies** \u2014 free ``einstein_radius``, no coupling.\n", + " - **extra galaxies** \u2014 free ``einstein_radius``, but with a luminosity-informed *upper bound*\n", + " ``min(2 * (upper_einstein_radius / L_BGC^0.5) * L^0.5, 5.0)`` \u2014 the BGC-scaled Faber-Jackson\n", + " prediction (x2, capped). It prevents a runaway mass while keeping the galaxy free. This is the\n", + " pipeline's middle tier.\n", + " - **scaling galaxies** \u2014 *tied*: ``einstein_radius_i = einstein_radius_BGC * (L/L_BGC) ** 0.5``. Zero\n", + " free parameters; the members inherit the BGC's mass (mass-anchored coupling, companion Model 3).\n", + "\n", + "So the magnitude of Model 1 splits its two jobs: the luminosity *ratio* still enters (now from the\n", + "fitted light), but the *normalization* is no longer a free ``sigma_ref`` \u2014 it is the BGC's own\n", + "``einstein_radius``.\n", + "\n", + "Total: 7 free (NFW halo 4 + 2 main einstein_radii + 1 extra einstein_radius + scaling 0)." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "\n", + "# (1) Halo -> elliptical NFW.\n", + "halo_mass = af.Model(al.mp.NFW)\n", + "halo_mass.centre = (0.0, 0.0) # [FIXED]\n", + "halo_mass.ell_comps.ell_comps_0 = af.UniformPrior(lower_limit=-0.5, upper_limit=0.5) # [FREE]\n", + "halo_mass.ell_comps.ell_comps_1 = af.UniformPrior(lower_limit=-0.5, upper_limit=0.5) # [FREE]\n", + "halo_mass.kappa_s = af.UniformPrior(lower_limit=0.05, upper_limit=0.5) # [FREE]\n", + "halo_mass.scale_radius = af.UniformPrior(lower_limit=10.0, upper_limit=60.0) # [FREE] arcsec\n", + "cluster_halo_pa = af.Model(al.Galaxy, redshift=redshift_lens, mass=halo_mass)\n", + "\n", + "# (2,3a) Main galaxies -> SIE, free einstein_radius. The BGC is the brightest.\n", + "main_galaxies_pa = []\n", + "for centre in main_centres:\n", + " mass = af.Model(al.mp.IsothermalSph)\n", + " mass.centre = centre # [FIXED]\n", + " mass.einstein_radius = af.UniformPrior(lower_limit=0.0, upper_limit=3.0) # [FREE] arcsec\n", + " main_galaxies_pa.append(af.Model(al.Galaxy, redshift=redshift_lens, mass=mass))\n", + "\n", + "einstein_radius_bgc = main_galaxies_pa[bgc_index].mass.einstein_radius # the BGC anchor\n", + "upper_einstein_radius = 3.0 # fixed proxy for the BGC Einstein radius (the main-lens prior upper bound)\n", + "\n", + "# (3b) Extra galaxies -> SIE, free einstein_radius with a luminosity-informed upper bound.\n", + "extra_galaxies_pa = []\n", + "for centre, magnitude in zip(extra_centres, extra_magnitudes):\n", + " ratio = luminosity_ratio(magnitude) # L / L_BGC\n", + " upper = min(2.0 * (upper_einstein_radius / 1.0**0.5) * ratio**0.5, 5.0) # L_BGC = 1 by construction\n", + " mass = af.Model(al.mp.IsothermalSph)\n", + " mass.centre = centre # [FIXED]\n", + " mass.einstein_radius = af.UniformPrior(lower_limit=0.0, upper_limit=upper) # [FREE, BOUNDED]\n", + " extra_galaxies_pa.append(af.Model(al.Galaxy, redshift=redshift_lens, mass=mass))\n", + "\n", + "# (3c) Scaling galaxies -> SIE, einstein_radius tied to the BGC.\n", + "scaling_galaxies_pa = []\n", + "for centre, magnitude in zip(scaling_centres, scaling_magnitudes):\n", + " ratio = luminosity_ratio(magnitude)\n", + " mass = af.Model(al.mp.IsothermalSph)\n", + " mass.centre = centre # [FIXED]\n", + " mass.einstein_radius = einstein_radius_bgc * ratio**0.5 # tied to the BGC -> NO new parameter\n", + " scaling_galaxies_pa.append(af.Model(al.Galaxy, redshift=redshift_lens, mass=mass))\n", + "\n", + "model_pyautolens = af.Collection(\n", + " cluster_halo=cluster_halo_pa,\n", + " main_galaxies=af.Collection(main_galaxies_pa),\n", + " extra_galaxies=af.Collection(extra_galaxies_pa),\n", + " scaling_galaxies=af.Collection(scaling_galaxies_pa),\n", + ")\n", + "\n", + "print(\"\\n\" + \"=\" * 80)\n", + "print(\"Section 2 \u2014 Standard PyAutoLens model (NFW halo, Isothermal galaxies)\")\n", + "print(\"=\" * 80)\n", + "print(f\"Free parameters: {model_pyautolens.prior_count} (NFW 4 + 2 main + 1 extra + scaling 0 = 7)\")\n", + "assert model_pyautolens.prior_count == 7\n", + "\n", + "# --- Numerical bridge: the scaling-galaxy parameter change is exact ---\n", + "C_KM_S = 299792.458\n", + "cosmology = ag.cosmo.FlatLambdaCDM(H0=H0, Om0=Om0)\n", + "d_s = cosmology.angular_diameter_distance_to_earth_in_kpc_from(redshift=redshift_source)\n", + "d_ls = cosmology.angular_diameter_distance_between_redshifts_in_kpc_from(\n", + " redshift_0=redshift_lens, redshift_1=redshift_source\n", + ")\n", + "grid = al.Grid2DIrregular([[1.5, 2.0], [-1.0, 0.8], [3.0, -2.5]])\n", + "\n", + "fj = [\n", + " (4.0 * np.pi * (s / C_KM_S) ** 2 * (d_ls / d_s) * (648000.0 / np.pi)) / s**2\n", + " for s in [150.0, 250.0, 350.0]\n", + "]\n", + "print(f\"(a) einstein_radius / sigma^2 = {fj[0]:.3e} (constant => einstein_radius ~ sigma^2, Faber-Jackson carries over)\")\n", + "assert np.allclose(fj, fj[0], rtol=1e-12)\n", + "\n", + "theta = 1.4\n", + "diff = np.max(\n", + " np.abs(\n", + " np.asarray(al.mp.IsothermalSph(centre=(0.0, 0.0), einstein_radius=theta).deflections_yx_2d_from(grid=grid))\n", + " - np.asarray(al.mp.dPIEMassB0Sph(centre=(0.0, 0.0), ra=1e-6, rs=1e7, b0=theta).deflections_yx_2d_from(grid=grid))\n", + " )\n", + ")\n", + "print(f\"(b) IsothermalSph(einstein_radius={theta}) vs dPIE b0 (SIS limit): max defl diff = {diff:.2e}\")\n", + "assert diff < 1e-5\n", + "print(\"Section 2 checks passed.\")\n" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "================================================================================\n", + "__Section 3 \u2014 Keeping the Truncation (dPIE outskirts, Isothermal main + extra)__\n", + "================================================================================\n", + "\n", + "Section 2 dropped the members' truncation. Here is how to keep it, changing as little as possible: the\n", + "halo (NFW), the main galaxies and the extra galaxies (all Isothermal) are unchanged from Section 2 \u2014\n", + "**only the scaling galaxies become truncated dPIEs**.\n", + "\n", + "**The difference truncation makes (all at large radius).**\n", + "\n", + " R (arcsec) SIE deflection dPIE (rs=20) deflection\n", + " 1 1.40 1.36\n", + " 20 1.40 0.82\n", + " 300 1.40 0.09\n", + "\n", + "The SIE has ``rho ~ r^-2`` forever: deflection stays at ``einstein_radius``, mass diverges (``M ~ R``).\n", + "The dPIE falls as ``r^-4`` beyond ``rs``: deflection dies away, total mass is *finite*. Inside the arc\n", + "region (R ~ 1-5\") they agree to a few percent \u2014 which is why the SIE is fine where the constraints live.\n", + "\n", + "**When MGLs do not need it.** ``rs`` models tidal *stripping* of a galaxy's outer halo by the host\n", + "potential. In a rich cluster stripping is strong and members are packed, so untruncated SIEs would each\n", + "contribute unbounded, overlapping mass \u2014 truncation is needed. In an MGL that is *not* a rich group the\n", + "potential is shallower, stripping weaker, and members are subdominant perturbers with the arcs near the\n", + "main lens \u2014 so the untruncated SIE is a defensible simplification. Keep truncation if the members are\n", + "numerous / massive enough to contribute meaningfully at large radius.\n", + "\n", + "**How, staying closest to the parameterization.** There is no separate \"truncated isothermal\" class \u2014\n", + "the dPIE *is* it (with ``ra -> 0`` it is a Pseudo-Jaffe). So the scaling galaxies become\n", + "``dPIEMassB0Sph`` with a tiny core and a truncation ``rs``, keeping ``b0`` (the angular lens strength) as\n", + "the mass parameter. The main and extra galaxies stay Isothermal.\n", + "\n", + "**The linking across parameterizations.** Both the SIE ``einstein_radius`` and the dPIE ``b0`` are the\n", + "*angular lens strength* in arcsec, both ``~ sigma^2``, equal in the SIS limit. So a dPIE scaling\n", + "galaxy's ``b0`` anchors *directly* to the SIE BGC's ``einstein_radius`` with the identical relation \u2014 no\n", + "conversion:\n", + "\n", + " b0_i = einstein_radius_BGC * (L / L_BGC) ** 0.5, rs_i = rs_ref * (L / L_BGC) ** 0.5\n", + "\n", + "The truncation ``rs`` is a *separate* parameter (the outer fall-off); it does not enter the anchoring.\n", + "Caveat: with finite ``rs`` the actual central deflection is a few percent below ``b0`` (truncation\n", + "removes outer mass) \u2014 negligible for subdominant members; you anchor on ``b0`` regardless." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "\n", + "# Halo (NFW), main and extra galaxies (SIE) \u2014 identical to Section 2.\n", + "halo_mass = af.Model(al.mp.NFW)\n", + "halo_mass.centre = (0.0, 0.0)\n", + "halo_mass.ell_comps.ell_comps_0 = af.UniformPrior(lower_limit=-0.5, upper_limit=0.5) # [FREE]\n", + "halo_mass.ell_comps.ell_comps_1 = af.UniformPrior(lower_limit=-0.5, upper_limit=0.5) # [FREE]\n", + "halo_mass.kappa_s = af.UniformPrior(lower_limit=0.05, upper_limit=0.5) # [FREE]\n", + "halo_mass.scale_radius = af.UniformPrior(lower_limit=10.0, upper_limit=60.0) # [FREE]\n", + "cluster_halo_trunc = af.Model(al.Galaxy, redshift=redshift_lens, mass=halo_mass)\n", + "\n", + "main_galaxies_trunc = []\n", + "for centre in main_centres:\n", + " mass = af.Model(al.mp.IsothermalSph)\n", + " mass.centre = centre # [FIXED]\n", + " mass.einstein_radius = af.UniformPrior(lower_limit=0.0, upper_limit=3.0) # [FREE]\n", + " main_galaxies_trunc.append(af.Model(al.Galaxy, redshift=redshift_lens, mass=mass))\n", + "\n", + "einstein_radius_bgc = main_galaxies_trunc[bgc_index].mass.einstein_radius # SIE Einstein radius = anchor\n", + "\n", + "extra_galaxies_trunc = []\n", + "for centre, magnitude in zip(extra_centres, extra_magnitudes):\n", + " ratio = luminosity_ratio(magnitude)\n", + " upper = min(2.0 * (3.0 / 1.0**0.5) * ratio**0.5, 5.0)\n", + " mass = af.Model(al.mp.IsothermalSph)\n", + " mass.centre = centre # [FIXED]\n", + " mass.einstein_radius = af.UniformPrior(lower_limit=0.0, upper_limit=upper) # [FREE, BOUNDED]\n", + " extra_galaxies_trunc.append(af.Model(al.Galaxy, redshift=redshift_lens, mass=mass))\n", + "\n", + "# Scaling galaxies -> dPIE (truncated), b0 anchored to the SIE BGC's einstein_radius.\n", + "rs_ref = 30.0 # [FIXED] reference truncation (arcsec); scales ~ L^0.5 like Lenstool's r_cut\n", + "scaling_galaxies_trunc = []\n", + "for centre, magnitude in zip(scaling_centres, scaling_magnitudes):\n", + " ratio = luminosity_ratio(magnitude)\n", + " mass = af.Model(al.mp.dPIEMassB0Sph)\n", + " mass.centre = centre # [FIXED]\n", + " mass.ra = 0.01 # [FIXED] ~0: a pure truncated isothermal (no core)\n", + " mass.b0 = einstein_radius_bgc * ratio**0.5 # SIE einstein_radius drives dPIE b0 -> NO new parameter\n", + " mass.rs = rs_ref * ratio**0.5 # [FIXED] truncation, scaled with L^0.5\n", + " scaling_galaxies_trunc.append(af.Model(al.Galaxy, redshift=redshift_lens, mass=mass))\n", + "\n", + "model_truncated = af.Collection(\n", + " cluster_halo=cluster_halo_trunc,\n", + " main_galaxies=af.Collection(main_galaxies_trunc),\n", + " extra_galaxies=af.Collection(extra_galaxies_trunc),\n", + " scaling_galaxies=af.Collection(scaling_galaxies_trunc),\n", + ")\n", + "\n", + "print(\"\\n\" + \"=\" * 80)\n", + "print(\"Section 3 \u2014 Truncated scaling galaxies (dPIE outskirts, Isothermal main + extra)\")\n", + "print(\"=\" * 80)\n", + "print(f\"Free parameters: {model_truncated.prior_count} (same 7 as Section 2; only the scaling profile changed)\")\n", + "assert model_truncated.prior_count == 7\n", + "\n", + "# (a) The cross-parameterization linking: a scaling galaxy's b0 = einstein_radius_BGC * (L/L_BGC)^0.5.\n", + "instance = model_truncated.instance_from_prior_medians()\n", + "er_bgc = instance.main_galaxies[bgc_index].mass.einstein_radius\n", + "print(f\"(a) BGC einstein_radius = {er_bgc:.3f} arcsec (SIE); dPIE scaling galaxies b0 = einstein_radius_BGC * (L/L_BGC)^0.5:\")\n", + "worst = 0.0\n", + "for i, magnitude in enumerate(scaling_magnitudes):\n", + " ratio = luminosity_ratio(magnitude)\n", + " expected = er_bgc * ratio**0.5\n", + " actual = instance.scaling_galaxies[i].mass.b0\n", + " worst = max(worst, abs(actual - expected))\n", + " print(f\" scaling {i} (mag {magnitude}): b0 = {actual:.3f} arcsec = {ratio**0.5:.3f} x einstein_radius_BGC\")\n", + "assert worst < 1e-9\n", + "\n", + "# (b) The truncation bites: the deflection dies away at large radius (an SIE would stay flat).\n", + "member0 = instance.scaling_galaxies[0].mass\n", + "c = scaling_centres[0]\n", + "print(\"(b) scaling galaxy 0 deflection vs radius (truncated dPIE \u2014 dies away; an SIE would stay flat):\")\n", + "for R in [1.0, 10.0, 50.0, 200.0]:\n", + " g = al.Grid2DIrregular([[c[0], c[1] + R]])\n", + " a = float(np.max(np.abs(np.asarray(member0.deflections_yx_2d_from(grid=g)))))\n", + " print(f\" R={R:>4.0f}\\\": alpha = {a:.4f} arcsec\")\n", + "print(\"Section 3 checks passed.\")\n" + ], + "outputs": [], + "execution_count": null + } + ], + "metadata": { + "anaconda-cloud": {}, + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.6.1" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} \ No newline at end of file diff --git a/scripts/cluster/mass_parameterizations.py b/scripts/cluster/mass_parameterizations.py new file mode 100644 index 000000000..8c996e3ae --- /dev/null +++ b/scripts/cluster/mass_parameterizations.py @@ -0,0 +1,593 @@ +""" +Cluster & Group Mass Parameterizations +====================================== + +**A plain-language tour of the ways to parameterize and couple the masses in a Lenstool-style cluster +model — the cluster-scale halo, the BCG and other individually-freed galaxies, and the member +("scaling") population.** + +There are several parameterizations floating around — velocity dispersion vs the internal ``b0``; +magnitudes vs luminosities; members tied to a scaling relation vs left free; the members normalized to +their own free parameter vs to the BCG's mass. They describe the *same physics*, so it is easy to lose +track of what actually differs. This script builds each option **with the model-composition API and +prints ``model.info``**, so for every model you can see at a glance *exactly which parameters are free +and which are fixed* — the thing that really distinguishes the approaches. + +Nothing here is fitted. Each section composes an ``af.Collection`` model and prints it (plus a few +assertions); the point is to *read* the models, not run them. + +__The axes of choice__ + +Every model below is a point in the same space of decisions: + + - **What number carries the mass?** ``sigma`` (velocity dispersion, km/s — Lenstool's native + parameter) or the internal ``b0`` (arcsec, the dPIE lens strength). Inter-convertible; which you + make *free* is the choice. + - **How are members coupled?** A Faber-Jackson scaling relation ties the whole member population to a + shared normalization; individual galaxies are left free (only feasible for a few). + - **What anchors the relation?** A reference *magnitude* ``mag0`` (Lenstool's ``potfile`` keyword). + - **What is the member normalization?** Its own free parameter (``sigma_ref``), or the BCG's own mass + (mass anchoring — the member tier then costs zero parameters). + - **How many relation parameters?** Just the dispersion normalization, or the truncation + normalization ``r_cut_ref`` as well. + +__Contents__ + +Each model is a **complete, self-contained** Lenstool-style cluster model — a cluster-scale halo, the +individually-modelled galaxies, and the member relation — so you can copy any one block wholesale +(there is deliberate redundancy between them). All use dPIE / Lenstool parameterizations only; mapping +these onto PyAutoLens-native profiles (isothermals, etc.) is a separate guide. + + - **Model 1 — The Standard Lenstool Model** — members on a velocity-dispersion Faber-Jackson relation + anchored to ``mag0``; BCG + a freed member modelled individually; a free elliptical halo. No ``b0``. + - **Model 2 — The Angular (b0) Parameterization** — the same model with ``sigma -> b0`` (redshift-free). + - **Model 3 — Mass Anchoring** — the members normalized to the BCG's own mass (member tier -> 0 free). + - **Model 4 — The Two-Parameter Relation** — free the truncation normalization ``r_cut_ref`` too. + +Freeing the Faber-Jackson *exponent* is a one-line change (make ``sigma_exponent`` a prior instead of +``0.25``) and is noted inline rather than given its own model. How Lenstool handles *unknown source +redshifts* (``z_m_limit``) is about the source, not the mass, and is out of scope here. + +__Reading model.info__ + +For every model we print ``model.info``. A line showing a ``UniformPrior`` (or any prior) is a **free** +parameter the search would sample; a bare number is **fixed**. A member's ``sigma`` shown as a +``MultiplePrior`` is *tied* to a shared free parameter, not free in its own right. +""" + +import numpy as np + +import autofit as af +import autolens as al +import autogalaxy as ag + +""" +__Shared Setup__ + +Redshifts and cosmology are the same as the ``cluster/simple`` example. The cosmology is the flat +``H0 = 67.66``, ``Om0 = 0.30966`` that the dPIE profile uses internally, so nothing is approximate. +For a real cluster, set these to the measured lens redshift and the run's cosmology. + +The photometry below is illustrative. In a real analysis it comes from a member catalogue (one +magnitude per galaxy in the band the scaling relation is calibrated in — e.g. F160W) plus each +galaxy's sky position. ``scaling_*`` are the fainter cluster members that go on the scaling relation; +``individual_*`` are the galaxies modelled on their own — the BCG (first) plus any bright member worth +freeing (e.g. a galaxy near a multiple image). The BCG's magnitude also anchors the relation (mag0). +""" +redshift_lens = 0.5 +redshift_source = 2.0 + +H0 = 67.66 +Om0 = 0.30966 + +# Scaling-relation members (the "potfile" population): (y, x) arc-second position + apparent magnitude. +scaling_centres = [(5.5, -6.5), (-7.5, 3.0), (3.0, 13.0), (12.0, -5.0), (-6.5, 11.0)] +scaling_magnitudes = [19.2, 19.8, 20.4, 21.0, 21.6] + +# Individually-modelled galaxies: the BCG (first) plus any bright member we want to free — e.g. a +# galaxy sitting near a multiple image, which the scaling relation cannot describe accurately enough +# locally. In Lenstool these are just separate "potentiel" sections; the BCG is not structurally +# special, only the first of them. The BCG's magnitude also anchors the scaling relation (mag0). +individual_centres = [(0.0, 0.0), (8.5, 5.5)] +individual_magnitudes = [17.8, 18.9] + + +""" +__Model 1 — The Standard Lenstool Model__ + +This is the model a Lenstool ``.par`` file describes. It has the three components of a real cluster +model, and this script builds all three: + + - a **cluster-scale dark-matter halo** — an elliptical dPIE, fully free, not tied to any light: the + dominant, smooth large-scale mass; + - **individually-modelled galaxies** — the BCG and any bright member worth freeing, each its own free + dPIE (Lenstool's individual "potentiel" sections; the BCG is simply the first); + - the **member population** ("potfile") — masses set by a velocity-dispersion Faber-Jackson relation + anchored to a reference magnitude. + +dPIE profiles throughout, composed directly in ``sigma`` (km/s), ``r_core`` and ``r_cut`` (arcsec) — +exactly the numbers a Lenstool results table quotes. No ``b0`` appears. + +**The reference magnitude ``mag0``.** The relation is anchored to a reference magnitude (the ``potfile`` +keyword ``mag0``); ``sigma_ref`` is the velocity dispersion of a galaxy *at that magnitude*. A member of +magnitude ``m`` enters through its brightness *relative* to the reference, +``L/L_ref = 10 ** (0.4 * (mag0 - m))``, which folds into Faber-Jackson (``sigma ~ L^0.25``) to give a +clean magnitude form: + + sigma_i = sigma_ref * 10 ** (0.1 * (mag0 - m_i)) # 0.25 * 0.4 = 0.1 + r_cut_i = r_cut_ref * 10 ** (0.2 * (mag0 - m_i)) # 0.5 * 0.4 = 0.2 (constant M/L) + r_core_i = r_core_ref * 10 ** (0.2 * (mag0 - m_i)) + +We set ``mag0`` to the BCG's magnitude — a natural anchor. Note this uses only the BCG's *brightness*; +the BCG's *mass* is one of the individually-modelled galaxies and is never coupled to the members. + +**Free vs fixed.** + + - Halo: free ``ellipticity``, ``angle_pos``, ``sigma``, ``r_core`` (4). Fixed centre (often freed near + the BCG) and ``r_cut`` (unconstrained within the field). + - Individually-modelled galaxies (BCG + freed member): each free ``sigma`` + ``r_cut`` (2 per galaxy); + fixed centre and ``r_core``. + - Members: one free parameter for the whole population, ``sigma_ref``. Adding members adds nothing. + - The exponents, reference radii, ``mag0``, redshifts and cosmology are fixed. + +Total: 9 free parameters (halo 4 + 2 individual galaxies x 2 + members 1). +""" + +# --- Cluster-scale dark-matter halo: elliptical dPIE, fully free, not tied to light --- +halo_mass = af.Model(al.mp.dPIEMass) +halo_mass.centre = (0.0, 0.0) # [FIXED] often freed near the BCG (set centre_0 / centre_1 to priors) +halo_mass.ellipticity = af.UniformPrior(lower_limit=0.0, upper_limit=0.7) # [FREE] +halo_mass.angle_pos = af.UniformPrior(lower_limit=0.0, upper_limit=180.0) # [FREE] degrees +halo_mass.sigma = af.UniformPrior(lower_limit=500.0, upper_limit=1500.0) # [FREE] cluster-scale dispersion (km/s) +halo_mass.r_core = af.UniformPrior(lower_limit=20.0, upper_limit=150.0) # [FREE] large halo core (arcsec) +halo_mass.r_cut = 1000.0 # [FIXED] truncation unconstrained within the field; fixed large +halo_mass.redshift_object = redshift_lens +halo_mass.redshift_source = redshift_source +halo_mass.H0 = H0 +halo_mass.Om0 = Om0 +cluster_halo = af.Model(al.Galaxy, redshift=redshift_lens, mass=halo_mass) + +# --- Individually-modelled galaxies: the BCG (first) + any freed member, each its own free dPIE --- +individual_sigma_priors = [(200.0, 600.0), (150.0, 450.0)] # per-galaxy free sigma range (km/s) +individual_rcut_priors = [(20.0, 200.0), (20.0, 150.0)] # per-galaxy free r_cut range (arcsec) + +individual_galaxies = [] +for centre, (sigma_lo, sigma_hi), (rcut_lo, rcut_hi) in zip( + individual_centres, individual_sigma_priors, individual_rcut_priors +): + mass = af.Model(al.mp.dPIEMassSph) + mass.centre = centre # [FIXED] observed position + mass.sigma = af.UniformPrior(lower_limit=sigma_lo, upper_limit=sigma_hi) # [FREE] this galaxy's own dispersion + mass.r_core = 0.3 # [FIXED] arcsec + mass.r_cut = af.UniformPrior(lower_limit=rcut_lo, upper_limit=rcut_hi) # [FREE] this galaxy's own truncation + mass.redshift_object = redshift_lens + mass.redshift_source = redshift_source + mass.H0 = H0 + mass.Om0 = Om0 + individual_galaxies.append(af.Model(al.Galaxy, redshift=redshift_lens, mass=mass)) + +# --- Members: one free sigma_ref, Faber-Jackson exponent 0.25 --- +sigma_ref = af.UniformPrior(lower_limit=100.0, upper_limit=400.0) # [FREE] km/s (dispersion at mag0) +mag0 = individual_magnitudes[0] # [FIXED] reference magnitude (Lenstool's mag0) = the BCG's brightness +sigma_exponent = 0.25 # [FIXED] Faber-Jackson (make this an af.UniformPrior to free it — a one-line change) +radius_exponent = 0.5 # [FIXED] constant mass-to-light: r_core, r_cut ~ L^0.5 +r_core_ref = 0.15 # [FIXED] arcsec +r_cut_ref = 20.0 # [FIXED] arcsec + +scaling_galaxies = [] +for centre, magnitude in zip(scaling_centres, scaling_magnitudes): + luminosity_ratio = 10.0 ** (0.4 * (mag0 - magnitude)) # L / L_ref from magnitudes + + mass = af.Model(al.mp.dPIEMassSph) + mass.centre = centre # [FIXED] + mass.sigma = sigma_ref * luminosity_ratio**sigma_exponent # tied to the one free sigma_ref + mass.r_core = r_core_ref * luminosity_ratio**radius_exponent + mass.r_cut = r_cut_ref * luminosity_ratio**radius_exponent + mass.redshift_object = redshift_lens + mass.redshift_source = redshift_source + mass.H0 = H0 + mass.Om0 = Om0 + scaling_galaxies.append(af.Model(al.Galaxy, redshift=redshift_lens, mass=mass)) + +model_1 = af.Collection( + cluster_halo=cluster_halo, + individual_galaxies=af.Collection(individual_galaxies), + scaling_galaxies=af.Collection(scaling_galaxies), +) + +print("=" * 80) +print("Model 1 — The Standard Lenstool Model") +print("=" * 80) +print(f"Total free parameters: {model_1.prior_count} (halo 4 + 2 individual x 2 + members 1 = 9)") +assert model_1.prior_count == 9 +print(model_1.info) + + +""" +__Model 2 — The Angular Parameterization: dPIEMassB0__ + +**How.** Swap every dPIE for its ``B0`` twin and compose in angular parameters: ``ra`` (core), +``rs`` (truncation) and ``b0`` (lens strength) — all arcsec, with **no redshift and no cosmology**. The +radii map one-to-one (``r_core -> ra``, ``r_cut -> rs``); the mass normalization moves ``sigma`` +(km/s) ``-> b0`` (arcsec). The scaling relation moves with it, and the exponent **doubles**: + + Model 1 (sigma): sigma_i = sigma_ref * (L/L_ref)^0.25 + Model 2 (b0): b0_i = b0_ref * (L/L_ref)^0.50 # because b0 ~ sigma^2 + +It is the *same* Faber-Jackson relation, written in the quantity the lensing actually constrains. The +halo becomes an elliptical ``dPIEMassB0`` (ellipticity via ``ell_comps``); individual galaxies and +members become ``dPIEMassB0Sph``. + +**Why.** The ``sigma`` parameterization needs the redshifts (``b0 = 6*648000*(sigma/c)^2*(D_LS/D_S)``), +so ``sigma`` only means something once ``z_lens`` / ``z_source`` are known. When they are *not*, the +images constrain only the angular ``b0``, and ``sigma`` is degenerate with ``D_LS/D_S``. Fitting ``b0`` +directly fits the quantity the data determines, in arcsec you can reason about. Convert back once the +redshifts are measured: ``sigma = c * sqrt(b0 / K)``, ``K = 6*648000*(D_LS/D_S)``. The ``B0`` mass +carries no redshift; the galaxy ``redshift`` is only a ray-tracing plane label. + +**Free vs fixed.** Identical structure to Model 1 (9 free): halo free ``ell_comps`` (2) + ``ra`` + ``b0``; +individual galaxies free ``b0`` + ``rs``; members one free ``b0_ref``. +""" + +# --- Halo: elliptical dPIEMassB0 (ellipticity via ell_comps) --- +halo_mass = af.Model(al.mp.dPIEMassB0) +halo_mass.centre = (0.0, 0.0) # [FIXED] +halo_mass.ell_comps.ell_comps_0 = af.UniformPrior(lower_limit=-0.5, upper_limit=0.5) # [FREE] +halo_mass.ell_comps.ell_comps_1 = af.UniformPrior(lower_limit=-0.5, upper_limit=0.5) # [FREE] +halo_mass.ra = af.UniformPrior(lower_limit=20.0, upper_limit=150.0) # [FREE] core (arcsec) +halo_mass.rs = 1000.0 # [FIXED] truncation, fixed large +halo_mass.b0 = af.UniformPrior(lower_limit=5.0, upper_limit=50.0) # [FREE] cluster-scale lens strength (arcsec) +cluster_halo = af.Model(al.Galaxy, redshift=redshift_lens, mass=halo_mass) + +# --- Individually-modelled galaxies: free b0 + free rs (angular, redshift-free) --- +individual_b0_priors = [(0.0, 5.0), (0.0, 3.0)] # arcsec, angular lens-strength range per galaxy +individual_rs_priors = [(20.0, 200.0), (20.0, 150.0)] # arcsec, truncation range per galaxy + +individual_galaxies_b0 = [] +for centre, (b0_lo, b0_hi), (rs_lo, rs_hi) in zip( + individual_centres, individual_b0_priors, individual_rs_priors +): + mass = af.Model(al.mp.dPIEMassB0Sph) + mass.centre = centre # [FIXED] + mass.b0 = af.UniformPrior(lower_limit=b0_lo, upper_limit=b0_hi) # [FREE] angular lens strength (arcsec) + mass.ra = 0.3 # [FIXED] = Model 1's r_core + mass.rs = af.UniformPrior(lower_limit=rs_lo, upper_limit=rs_hi) # [FREE] truncation (arcsec) + individual_galaxies_b0.append(af.Model(al.Galaxy, redshift=redshift_lens, mass=mass)) + +# --- Members: one free b0_ref, exponent 0.5 (b0 ~ L^0.5) --- +b0_ref = af.UniformPrior(lower_limit=0.0, upper_limit=1.0) # [FREE] arcsec (lens strength at mag0) +b0_exponent = 0.5 # [FIXED] b0 ~ L^0.5 (Model 1 used sigma ~ L^0.25; b0 ~ sigma^2 doubles the exponent) +mag0 = individual_magnitudes[0] # [FIXED] +ra_ref = 0.15 # [FIXED] = Model 1's r_core_ref +rs_ref = 20.0 # [FIXED] = Model 1's r_cut_ref + +scaling_galaxies_b0 = [] +for centre, magnitude in zip(scaling_centres, scaling_magnitudes): + luminosity_ratio = 10.0 ** (0.4 * (mag0 - magnitude)) + + mass = af.Model(al.mp.dPIEMassB0Sph) + mass.centre = centre # [FIXED] + mass.b0 = b0_ref * luminosity_ratio**b0_exponent # tied to the one free b0_ref + mass.ra = ra_ref * luminosity_ratio**0.5 + mass.rs = rs_ref * luminosity_ratio**0.5 + scaling_galaxies_b0.append(af.Model(al.Galaxy, redshift=redshift_lens, mass=mass)) + +model_2 = af.Collection( + cluster_halo=cluster_halo, + individual_galaxies=af.Collection(individual_galaxies_b0), + scaling_galaxies=af.Collection(scaling_galaxies_b0), +) + +print("\n" + "=" * 80) +print("Model 2 — The Angular (b0) Parameterization") +print("=" * 80) +print(f"Total free parameters: {model_2.prior_count} (same structure as Model 1)") +assert model_2.prior_count == 9 + +""" +__Model 2 — Numerical Checks__ + +Model 2 is Model 1 re-expressed, so it must reproduce Model 1's physics once the redshifts are +supplied. We relate the two with the distance factor ``K`` (redshifts + cosmology used *only* for the +comparison, never inside the b0 model), and check three things: + + (a) the Faber-Jackson relation is consistent — a member's ``b0`` from the ``b0``-relation equals + ``K (sigma/c)^2`` with ``sigma`` from the ``sigma``-relation (the exponent doubling is exact); + (b) the deflections are identical, profile-for-profile; + (c) a fitted ``b0_ref`` converts back to the Model 1 ``sigma_ref`` exactly. +""" +C_KM_S = 299792.458 +cosmology = ag.cosmo.FlatLambdaCDM(H0=H0, Om0=Om0) # matches dPIEMass's internal cosmology +d_s = cosmology.angular_diameter_distance_to_earth_in_kpc_from(redshift=redshift_source) +d_ls = cosmology.angular_diameter_distance_between_redshifts_in_kpc_from( + redshift_0=redshift_lens, redshift_1=redshift_source +) +K = 6.0 * 648000.0 * (d_ls / d_s) # b0 = K * (sigma / c)^2 + +grid = al.Grid2DIrregular([[3.0, 4.0], [-8.0, 2.0], [10.0, -6.0]]) + +# A worked value: the b0_ref a redshift-free fit would return for sigma_ref = 250 km/s. +sigma_ref_value = 250.0 +b0_ref_value = K * (sigma_ref_value / C_KM_S) ** 2 + +# (a) Faber-Jackson consistency. +worst_relation = 0.0 +for magnitude in scaling_magnitudes: + ratio = 10.0 ** (0.4 * (mag0 - magnitude)) + worst_relation = max( + worst_relation, + abs(b0_ref_value * ratio**0.5 - K * (sigma_ref_value * ratio**0.25 / C_KM_S) ** 2), + ) +print(f"(a) member b0 [b0-relation vs sigma-relation]: max diff = {worst_relation:.2e}") +assert worst_relation < 1e-12 + +# (b) Deflections identical: dPIEMassB0Sph(b0) == dPIEMassSph(sigma) for each member. +worst_defl = 0.0 +for centre, magnitude in zip(scaling_centres, scaling_magnitudes): + ratio = 10.0 ** (0.4 * (mag0 - magnitude)) + mass_b0 = al.mp.dPIEMassB0Sph( + centre=centre, ra=ra_ref * ratio**0.5, rs=rs_ref * ratio**0.5, b0=b0_ref_value * ratio**0.5 + ) + mass_sigma = al.mp.dPIEMassSph( + centre=centre, + sigma=sigma_ref_value * ratio**0.25, + r_core=ra_ref * ratio**0.5, + r_cut=rs_ref * ratio**0.5, + redshift_object=redshift_lens, + redshift_source=redshift_source, + H0=H0, + Om0=Om0, + ) + worst_defl = max( + worst_defl, + np.max( + np.abs( + np.asarray(mass_b0.deflections_yx_2d_from(grid=grid)) + - np.asarray(mass_sigma.deflections_yx_2d_from(grid=grid)) + ) + ), + ) +print(f"(b) member deflections [b0 vs sigma param]: max diff = {worst_defl:.2e}") +assert worst_defl < 1e-12 + +# (c) Round-trip. +sigma_ref_recovered = C_KM_S * (b0_ref_value / K) ** 0.5 +print(f"(c) b0_ref = {b0_ref_value:.5f} arcsec -> sigma_ref = {sigma_ref_recovered:.2f} km/s (input {sigma_ref_value})") +assert np.isclose(sigma_ref_recovered, sigma_ref_value, rtol=1e-9) +print("Model 2 checks passed.") + + +""" +__Model 3 — Mass Anchoring: the Members Normalized to the BCG__ + +**How.** Take Model 1 and *delete* the members' free normalization ``sigma_ref``. In its place the +members are normalized to the BCG's OWN velocity dispersion — the same free ``sigma`` that sets the +BCG's mass: + + Model 1: sigma_i = sigma_ref * (L/L_BCG)^0.25 # sigma_ref FREE, independent of the BCG + Model 3: sigma_i = sigma_BCG * (L/L_BCG)^0.25 # sigma_BCG = the BCG's own free mass + +Because ``sigma_BCG`` is already free (it *is* the BCG's mass), the member population now adds **zero** +free parameters. This is the referee's "anchor alpha to the BCG Einstein radius" — ``sigma_BCG`` and +the BCG Einstein radius are the same quantity up to the fixed conversion. + +**Why.** The most economical model: one free normalization (the BCG's, tightly pinned by the central +images) describes everyone, with an interpretable prior. You save a parameter. + +**The cost — and why this is NOT the default (Model 1).** It *assumes the BCG lies on the member +relation*; real BCGs frequently deviate (merger products at the bottom of the potential well), which is +why high-precision models free the BCG separately (Model 1). And it hard-couples every member's mass to +the BCG's — move the fitted BCG mass and all members move in fixed proportion. A real, used variant +(simpler models, the referee's convention), but a deliberate coupling choice, not the standard. + +**Free vs fixed.** Halo as Model 1 (4). BCG: free ``sigma`` (its own mass *and* the member anchor) + +``r_cut``. Other freed galaxy: its own free ``sigma`` + ``r_cut``. Members: **zero** free. Total: 8. +""" + +# --- Halo: identical to Model 1 (the halo is unaffected by how the members are anchored) --- +halo_mass = af.Model(al.mp.dPIEMass) +halo_mass.centre = (0.0, 0.0) # [FIXED] +halo_mass.ellipticity = af.UniformPrior(lower_limit=0.0, upper_limit=0.7) # [FREE] +halo_mass.angle_pos = af.UniformPrior(lower_limit=0.0, upper_limit=180.0) # [FREE] +halo_mass.sigma = af.UniformPrior(lower_limit=500.0, upper_limit=1500.0) # [FREE] +halo_mass.r_core = af.UniformPrior(lower_limit=20.0, upper_limit=150.0) # [FREE] +halo_mass.r_cut = 1000.0 # [FIXED] +halo_mass.redshift_object = redshift_lens +halo_mass.redshift_source = redshift_source +halo_mass.H0 = H0 +halo_mass.Om0 = Om0 +cluster_halo = af.Model(al.Galaxy, redshift=redshift_lens, mass=halo_mass) + +# --- The BCG: the ANCHOR. Its free sigma sets its own mass AND normalizes every member. --- +bcg_mass = af.Model(al.mp.dPIEMassSph) +bcg_mass.centre = individual_centres[0] # [FIXED] +bcg_mass.sigma = af.UniformPrior(lower_limit=200.0, upper_limit=600.0) # [FREE] BCG mass + member anchor +bcg_mass.r_core = 0.3 # [FIXED] +bcg_mass.r_cut = af.UniformPrior(lower_limit=20.0, upper_limit=200.0) # [FREE] BCG's own truncation +bcg_mass.redshift_object = redshift_lens +bcg_mass.redshift_source = redshift_source +bcg_mass.H0 = H0 +bcg_mass.Om0 = Om0 +bcg = af.Model(al.Galaxy, redshift=redshift_lens, mass=bcg_mass) + +# --- Other individually-freed galaxies (still independent — e.g. a member on an arc) --- +individual_sigma_priors = [(200.0, 600.0), (150.0, 450.0)] +individual_rcut_priors = [(20.0, 200.0), (20.0, 150.0)] + +other_individual_galaxies = [] +for centre, (sigma_lo, sigma_hi), (rcut_lo, rcut_hi) in zip( + individual_centres[1:], individual_sigma_priors[1:], individual_rcut_priors[1:] +): + mass = af.Model(al.mp.dPIEMassSph) + mass.centre = centre + mass.sigma = af.UniformPrior(lower_limit=sigma_lo, upper_limit=sigma_hi) # [FREE] independent of the BCG + mass.r_core = 0.3 + mass.r_cut = af.UniformPrior(lower_limit=rcut_lo, upper_limit=rcut_hi) # [FREE] + mass.redshift_object = redshift_lens + mass.redshift_source = redshift_source + mass.H0 = H0 + mass.Om0 = Om0 + other_individual_galaxies.append(af.Model(al.Galaxy, redshift=redshift_lens, mass=mass)) + +# --- Members: normalized to the BCG's OWN sigma. ZERO new free parameters. --- +mag0 = individual_magnitudes[0] # [FIXED] = BCG magnitude +r_core_ref = 0.15 # [FIXED] +r_cut_ref = 20.0 # [FIXED] + +mass_anchored_members = [] +for centre, magnitude in zip(scaling_centres, scaling_magnitudes): + luminosity_ratio = 10.0 ** (0.4 * (mag0 - magnitude)) # L / L_BCG + + mass = af.Model(al.mp.dPIEMassSph) + mass.centre = centre + mass.sigma = bcg_mass.sigma * luminosity_ratio**0.25 # tied to the BCG's free sigma -> NO new parameter + mass.r_core = r_core_ref * luminosity_ratio**0.5 + mass.r_cut = r_cut_ref * luminosity_ratio**0.5 + mass.redshift_object = redshift_lens + mass.redshift_source = redshift_source + mass.H0 = H0 + mass.Om0 = Om0 + mass_anchored_members.append(af.Model(al.Galaxy, redshift=redshift_lens, mass=mass)) + +model_3 = af.Collection( + cluster_halo=cluster_halo, + individual_galaxies=af.Collection([bcg] + other_individual_galaxies), + scaling_galaxies=af.Collection(mass_anchored_members), +) + +print("\n" + "=" * 80) +print("Model 3 — Mass Anchoring (members normalized to the BCG)") +print("=" * 80) +print(f"Total free parameters: {model_3.prior_count} (Model 1 was 9; the member tier now adds 0)") +assert model_3.prior_count == 8 + +""" +__Model 3 — Numerical Checks__ + + (a) the member collection adds ZERO free parameters (all 8 belong to the halo + individual galaxies); + (b) at any instance, each member's sigma is exactly ``sigma_BCG * (L/L_BCG)^0.25`` — a fixed fraction + of the BCG, so a change in the fitted BCG mass rescales them all together. +""" +# (a) The member tier adds zero free parameters. +without_members = af.Collection( + cluster_halo=cluster_halo, + individual_galaxies=af.Collection([bcg] + other_individual_galaxies), +) +added = model_3.prior_count - without_members.prior_count +print(f"(a) free params: without members = {without_members.prior_count}, with members = {model_3.prior_count} (members added {added})") +assert added == 0 + +# (b) The members are driven by the BCG's sigma. +instance = model_3.instance_from_prior_medians() +sigma_bcg_instance = instance.individual_galaxies[0].mass.sigma +print(f"(b) instance BCG sigma = {sigma_bcg_instance:.2f} km/s; members = sigma_BCG * (L/L_BCG)^0.25:") +worst_coupling = 0.0 +for i, magnitude in enumerate(scaling_magnitudes): + ratio = 10.0 ** (0.4 * (mag0 - magnitude)) + expected = sigma_bcg_instance * ratio**0.25 + actual = instance.scaling_galaxies[i].mass.sigma + worst_coupling = max(worst_coupling, abs(actual - expected)) + print(f" member {i} (mag {magnitude}): sigma = {actual:6.2f} km/s = {ratio**0.25:.3f} x sigma_BCG") +assert worst_coupling < 1e-9 +print("Model 3 checks passed.") + + +""" +__Model 4 — The Two-Parameter Relation__ + +**How.** Model 1 fixes the truncation reference ``r_cut_ref`` and frees only ``sigma_ref``. Here we +free ``r_cut_ref`` as well, so the member population has **two** shared free parameters — a dispersion +normalization and a truncation normalization: + + sigma_i = sigma_ref * (L/L_ref)^0.25 # sigma_ref FREE (as Model 1) + r_cut_i = r_cut_ref * (L/L_ref)^0.5 # r_cut_ref now FREE too + +Every member's truncation now scales from a *free* reference instead of a fixed one. It is a one-symbol +change (``r_cut_ref`` becomes a prior) — but it is a real modelling decision, so it earns a section. + +**Why.** The reference truncation ``r_cut*`` is genuinely uncertain — it sets how much mass each member +carries at large radius, and the data can constrain it. Lenstool's ``potfile`` routinely optimizes both +``sigma*`` and ``r_cut*``. You free it when the members contribute enough lensing that their outer mass +matters; you fix it (Model 1) when they do not, to save a parameter and avoid a weakly-constrained +direction. + +**Free vs fixed.** As Model 1, but the member tier now has 2 free parameters (``sigma_ref`` + +``r_cut_ref``) instead of 1. Total: 10. +""" + +# --- Halo: as Model 1 --- +halo_mass = af.Model(al.mp.dPIEMass) +halo_mass.centre = (0.0, 0.0) # [FIXED] +halo_mass.ellipticity = af.UniformPrior(lower_limit=0.0, upper_limit=0.7) # [FREE] +halo_mass.angle_pos = af.UniformPrior(lower_limit=0.0, upper_limit=180.0) # [FREE] +halo_mass.sigma = af.UniformPrior(lower_limit=500.0, upper_limit=1500.0) # [FREE] +halo_mass.r_core = af.UniformPrior(lower_limit=20.0, upper_limit=150.0) # [FREE] +halo_mass.r_cut = 1000.0 # [FIXED] +halo_mass.redshift_object = redshift_lens +halo_mass.redshift_source = redshift_source +halo_mass.H0 = H0 +halo_mass.Om0 = Om0 +cluster_halo = af.Model(al.Galaxy, redshift=redshift_lens, mass=halo_mass) + +# --- Individually-modelled galaxies: as Model 1 --- +individual_sigma_priors = [(200.0, 600.0), (150.0, 450.0)] +individual_rcut_priors = [(20.0, 200.0), (20.0, 150.0)] + +individual_galaxies = [] +for centre, (sigma_lo, sigma_hi), (rcut_lo, rcut_hi) in zip( + individual_centres, individual_sigma_priors, individual_rcut_priors +): + mass = af.Model(al.mp.dPIEMassSph) + mass.centre = centre + mass.sigma = af.UniformPrior(lower_limit=sigma_lo, upper_limit=sigma_hi) # [FREE] + mass.r_core = 0.3 + mass.r_cut = af.UniformPrior(lower_limit=rcut_lo, upper_limit=rcut_hi) # [FREE] + mass.redshift_object = redshift_lens + mass.redshift_source = redshift_source + mass.H0 = H0 + mass.Om0 = Om0 + individual_galaxies.append(af.Model(al.Galaxy, redshift=redshift_lens, mass=mass)) + +# --- Members: TWO free normalizations, sigma_ref and r_cut_ref --- +sigma_ref = af.UniformPrior(lower_limit=100.0, upper_limit=400.0) # [FREE] km/s +r_cut_ref = af.UniformPrior(lower_limit=5.0, upper_limit=50.0) # [FREE] arcsec <-- now a prior, not a constant +mag0 = individual_magnitudes[0] # [FIXED] +sigma_exponent = 0.25 # [FIXED] +radius_exponent = 0.5 # [FIXED] +r_core_ref = 0.15 # [FIXED] (the core reference is usually kept fixed even in the 2-parameter relation) + +scaling_galaxies = [] +for centre, magnitude in zip(scaling_centres, scaling_magnitudes): + luminosity_ratio = 10.0 ** (0.4 * (mag0 - magnitude)) + + mass = af.Model(al.mp.dPIEMassSph) + mass.centre = centre + mass.sigma = sigma_ref * luminosity_ratio**sigma_exponent # tied to free sigma_ref + mass.r_core = r_core_ref * luminosity_ratio**radius_exponent # fixed reference + mass.r_cut = r_cut_ref * luminosity_ratio**radius_exponent # tied to free r_cut_ref + mass.redshift_object = redshift_lens + mass.redshift_source = redshift_source + mass.H0 = H0 + mass.Om0 = Om0 + scaling_galaxies.append(af.Model(al.Galaxy, redshift=redshift_lens, mass=mass)) + +model_4 = af.Collection( + cluster_halo=cluster_halo, + individual_galaxies=af.Collection(individual_galaxies), + scaling_galaxies=af.Collection(scaling_galaxies), +) + +print("\n" + "=" * 80) +print("Model 4 — The Two-Parameter Relation") +print("=" * 80) +print(f"Total free parameters: {model_4.prior_count} (Model 1 was 9; the member tier is now 2, not 1)") +assert model_4.prior_count == 10 + +# The member tier now contributes 2 free parameters (sigma_ref + r_cut_ref) instead of 1. +without_members = af.Collection( + cluster_halo=cluster_halo, individual_galaxies=af.Collection(individual_galaxies) +) +added = model_4.prior_count - without_members.prior_count +print(f"member tier free parameters: {added} (sigma_ref + r_cut_ref)") +assert added == 2 +print("Model 4 checks passed.") diff --git a/scripts/cluster/mass_parameterizations_pyautolens.py b/scripts/cluster/mass_parameterizations_pyautolens.py new file mode 100644 index 000000000..345d36f96 --- /dev/null +++ b/scripts/cluster/mass_parameterizations_pyautolens.py @@ -0,0 +1,393 @@ +""" +Mass Parameterizations II — Mapping Lenstool onto PyAutoLens +============================================================ + +**For Lenstool users: how and why PyAutoLens re-expresses the standard Lenstool cluster / group model +in its native parameterization for multi-galaxy lenses (MGLs), and how to keep the dPIE's truncation on +the outskirts if you want it.** + +The companion script ``cluster/mass_parameterizations.py`` builds the standard Lenstool model. This +script takes its **Model 1** as the starting point and maps it, component by component, onto what +PyAutoLens uses at galaxy / group scale. Both are supported; this is an *alternative*, not a +replacement, and the focus is the **mass**. + +Three sections: + + - **Section 1 — Lenstool Model 1** (dPIE throughout, velocity dispersions, a magnitude-anchored + Faber-Jackson relation), the starting point. + - **Section 2 — The standard PyAutoLens model**, mapping every component and justifying each change. + - **Section 3 — Keeping truncation**: dPIE for the scaling galaxies while the main and extra galaxies + stay Isothermal, and how the Faber-Jackson linking bridges the two parameterizations. + +Every model has four tiers, and this is the structure of a PyAutoLens MGL: + + Tier Role Coupling to the BGC + ---------------- -------------------------------------- -------------------------------- + cluster halo the smooth large-scale mass none (its own free profile) + main galaxies the primary lens galaxies (incl. BGC) none — free + extra galaxies bright members worth freeing free, but luminosity-BOUNDED + scaling galaxies the faint outskirts ("potfile") TIED — mass set by the relation + +The luminosity that drives the relation is computed from each galaxy's fitted MGE light, +``L = sum_g 2*pi*sigma_g^2/q_g * intensity_g`` (the ``2*np.pi*g.sigma**2/g.axis_ratio()*g.intensity`` +you have seen); only *ratios* matter, so we drive everything off illustrative magnitudes here +(``L/L_BGC = 10 ** (0.4 * (mag0 - m))``). + +Everything uses the model-composition API and prints ``model.info`` so free vs fixed is explicit. + +__What changes, at a glance__ + + Component Lenstool (Section 1) PyAutoLens (Sections 2-3) + ---------------- ---------------------------- -------------------------------------- + cluster halo elliptical dPIE (sigma) NFW (the CDM halo profile) + main galaxies dPIE (sigma, r_core, r_cut) Isothermal (einstein_radius) + extra galaxies dPIE, free Isothermal, free but luminosity-bounded + scaling galaxies dPIE, free sigma_ref Isothermal, TIED to the BGC's Einstein radius + mass parameter sigma (needs redshifts) einstein_radius (arcsec; redshift-free) + luminosity input magnitude computed from the fitted MGE light +""" + +import numpy as np + +import autofit as af +import autolens as al +import autogalaxy as ag + +""" +__Setup__ + +Redshifts / cosmology are needed by the *Lenstool* (sigma) side and again to *check* the +sigma <-> einstein_radius bridge — they do not appear in the PyAutoLens mass model. +""" +redshift_lens = 0.5 +redshift_source = 2.0 +H0 = 67.66 +Om0 = 0.30966 + +main_centres = [(0.0, 0.0), (12.0, 8.0)] # primary lens galaxies (the brightest is the BGC) +main_magnitudes = [17.8, 18.9] + +extra_centres = [(8.5, 5.5)] # a bright member worth freeing (e.g. near a multiple image) +extra_magnitudes = [18.5] + +scaling_centres = [(5.5, -6.5), (-7.5, 3.0), (3.0, 13.0)] # the faint outskirts +scaling_magnitudes = [19.2, 20.0, 21.0] + +bgc_index = int(np.argmin(main_magnitudes)) # brightest (smallest magnitude) main galaxy +mag0 = main_magnitudes[bgc_index] # the BGC magnitude anchors the luminosity ratios + + +def luminosity_ratio(magnitude): + """L / L_BGC from a magnitude, relative to the BGC (mag0).""" + return 10.0 ** (0.4 * (mag0 - magnitude)) + + +""" +================================================================================ +__Section 1 — Lenstool Model 1__ +================================================================================ + +The full standard Lenstool model, dPIE throughout (see ``cluster/mass_parameterizations.py`` Model 1). +The main and extra galaxies are both individually-modelled dPIEs (free ``sigma`` + ``r_cut``); the +distinction between them is a PyAutoLens refinement (Section 2) — in Lenstool they are both just +"potentiel" sections. The scaling galaxies use a magnitude-anchored Faber-Jackson relation with a free +normalization ``sigma_ref``: ``sigma_i = sigma_ref * (L/L_ref) ** 0.25``. + +Total: 11 free (halo 4 + 2 main x 2 + 1 extra x 2 + scaling 1). +""" + + +def dpie_galaxy(centre, sigma_lo, sigma_hi, rcut_lo, rcut_hi): + """An individually-modelled Lenstool galaxy: free dPIE dispersion + truncation.""" + mass = af.Model(al.mp.dPIEMassSph) + mass.centre = centre + mass.sigma = af.UniformPrior(lower_limit=sigma_lo, upper_limit=sigma_hi) # [FREE] km/s + mass.r_core = 0.3 + mass.r_cut = af.UniformPrior(lower_limit=rcut_lo, upper_limit=rcut_hi) # [FREE] arcsec + mass.redshift_object = redshift_lens + mass.redshift_source = redshift_source + mass.H0 = H0 + mass.Om0 = Om0 + return af.Model(al.Galaxy, redshift=redshift_lens, mass=mass) + + +# Cluster halo — elliptical dPIE, fully free. +halo_mass = af.Model(al.mp.dPIEMass) +halo_mass.centre = (0.0, 0.0) # [FIXED] +halo_mass.ellipticity = af.UniformPrior(lower_limit=0.0, upper_limit=0.7) # [FREE] +halo_mass.angle_pos = af.UniformPrior(lower_limit=0.0, upper_limit=180.0) # [FREE] +halo_mass.sigma = af.UniformPrior(lower_limit=500.0, upper_limit=1500.0) # [FREE] km/s +halo_mass.r_core = af.UniformPrior(lower_limit=20.0, upper_limit=150.0) # [FREE] arcsec +halo_mass.r_cut = 1000.0 # [FIXED] +halo_mass.redshift_object = redshift_lens +halo_mass.redshift_source = redshift_source +halo_mass.H0 = H0 +halo_mass.Om0 = Om0 +cluster_halo_lt = af.Model(al.Galaxy, redshift=redshift_lens, mass=halo_mass) + +main_galaxies_lt = [dpie_galaxy(c, 100.0, 600.0, 20.0, 200.0) for c in main_centres] +extra_galaxies_lt = [dpie_galaxy(c, 100.0, 450.0, 20.0, 150.0) for c in extra_centres] + +# Scaling galaxies — Faber-Jackson on a FREE normalization sigma_ref (independent of the BGC). +sigma_ref = af.UniformPrior(lower_limit=100.0, upper_limit=400.0) # [FREE] km/s +r_core_ref, r_cut_ref = 0.15, 20.0 # [FIXED] arcsec +scaling_galaxies_lt = [] +for centre, magnitude in zip(scaling_centres, scaling_magnitudes): + ratio = luminosity_ratio(magnitude) + mass = af.Model(al.mp.dPIEMassSph) + mass.centre = centre # [FIXED] + mass.sigma = sigma_ref * ratio**0.25 # tied to the one free sigma_ref + mass.r_core = r_core_ref * ratio**0.5 + mass.r_cut = r_cut_ref * ratio**0.5 + mass.redshift_object = redshift_lens + mass.redshift_source = redshift_source + mass.H0 = H0 + mass.Om0 = Om0 + scaling_galaxies_lt.append(af.Model(al.Galaxy, redshift=redshift_lens, mass=mass)) + +model_lenstool = af.Collection( + cluster_halo=cluster_halo_lt, + main_galaxies=af.Collection(main_galaxies_lt), + extra_galaxies=af.Collection(extra_galaxies_lt), + scaling_galaxies=af.Collection(scaling_galaxies_lt), +) + +print("=" * 80) +print("Section 1 — Lenstool Model 1 (dPIE throughout)") +print("=" * 80) +print(f"Free parameters: {model_lenstool.prior_count} (halo 4 + 2 main x 2 + 1 extra x 2 + scaling 1 = 11)") +assert model_lenstool.prior_count == 11 + + +""" +================================================================================ +__Section 2 — The Standard PyAutoLens Model__ +================================================================================ + +We map Model 1 onto the PyAutoLens MGL parameterization. Component by component: + +**(1) The halo: dPIE -> NFW.** Lenstool models the smooth halo as a dPIE for analytic convenience; +PyAutoLens uses the **NFW**, the CDM-motivated halo profile. A genuine *profile change*, not a +reparameterization — you refit ``kappa_s`` / ``scale_radius`` rather than convert. (Keep a dPIE halo if +you want to reproduce a Lenstool halo; both are supported.) + +**(2) Galaxies: dPIE -> Isothermal (SIE).** One mass parameter, ``einstein_radius`` — the reduced +deflection, in arcsec, redshift-free, exactly what the images constrain. It equals the dPIE ``b0`` in +the SIS limit (checked below). No core, no truncation. + +**(3) The coupling spectrum.** The individually-modelled galaxies split into two tiers, and the scaling +tier is anchored to the brightest galaxy (the BGC): + + - **main galaxies** — free ``einstein_radius``, no coupling. + - **extra galaxies** — free ``einstein_radius``, but with a luminosity-informed *upper bound* + ``min(2 * (upper_einstein_radius / L_BGC^0.5) * L^0.5, 5.0)`` — the BGC-scaled Faber-Jackson + prediction (x2, capped). It prevents a runaway mass while keeping the galaxy free. This is the + pipeline's middle tier. + - **scaling galaxies** — *tied*: ``einstein_radius_i = einstein_radius_BGC * (L/L_BGC) ** 0.5``. Zero + free parameters; the members inherit the BGC's mass (mass-anchored coupling, companion Model 3). + +So the magnitude of Model 1 splits its two jobs: the luminosity *ratio* still enters (now from the +fitted light), but the *normalization* is no longer a free ``sigma_ref`` — it is the BGC's own +``einstein_radius``. + +Total: 7 free (NFW halo 4 + 2 main einstein_radii + 1 extra einstein_radius + scaling 0). +""" + +# (1) Halo -> elliptical NFW. +halo_mass = af.Model(al.mp.NFW) +halo_mass.centre = (0.0, 0.0) # [FIXED] +halo_mass.ell_comps.ell_comps_0 = af.UniformPrior(lower_limit=-0.5, upper_limit=0.5) # [FREE] +halo_mass.ell_comps.ell_comps_1 = af.UniformPrior(lower_limit=-0.5, upper_limit=0.5) # [FREE] +halo_mass.kappa_s = af.UniformPrior(lower_limit=0.05, upper_limit=0.5) # [FREE] +halo_mass.scale_radius = af.UniformPrior(lower_limit=10.0, upper_limit=60.0) # [FREE] arcsec +cluster_halo_pa = af.Model(al.Galaxy, redshift=redshift_lens, mass=halo_mass) + +# (2,3a) Main galaxies -> SIE, free einstein_radius. The BGC is the brightest. +main_galaxies_pa = [] +for centre in main_centres: + mass = af.Model(al.mp.IsothermalSph) + mass.centre = centre # [FIXED] + mass.einstein_radius = af.UniformPrior(lower_limit=0.0, upper_limit=3.0) # [FREE] arcsec + main_galaxies_pa.append(af.Model(al.Galaxy, redshift=redshift_lens, mass=mass)) + +einstein_radius_bgc = main_galaxies_pa[bgc_index].mass.einstein_radius # the BGC anchor +upper_einstein_radius = 3.0 # fixed proxy for the BGC Einstein radius (the main-lens prior upper bound) + +# (3b) Extra galaxies -> SIE, free einstein_radius with a luminosity-informed upper bound. +extra_galaxies_pa = [] +for centre, magnitude in zip(extra_centres, extra_magnitudes): + ratio = luminosity_ratio(magnitude) # L / L_BGC + upper = min(2.0 * (upper_einstein_radius / 1.0**0.5) * ratio**0.5, 5.0) # L_BGC = 1 by construction + mass = af.Model(al.mp.IsothermalSph) + mass.centre = centre # [FIXED] + mass.einstein_radius = af.UniformPrior(lower_limit=0.0, upper_limit=upper) # [FREE, BOUNDED] + extra_galaxies_pa.append(af.Model(al.Galaxy, redshift=redshift_lens, mass=mass)) + +# (3c) Scaling galaxies -> SIE, einstein_radius tied to the BGC. +scaling_galaxies_pa = [] +for centre, magnitude in zip(scaling_centres, scaling_magnitudes): + ratio = luminosity_ratio(magnitude) + mass = af.Model(al.mp.IsothermalSph) + mass.centre = centre # [FIXED] + mass.einstein_radius = einstein_radius_bgc * ratio**0.5 # tied to the BGC -> NO new parameter + scaling_galaxies_pa.append(af.Model(al.Galaxy, redshift=redshift_lens, mass=mass)) + +model_pyautolens = af.Collection( + cluster_halo=cluster_halo_pa, + main_galaxies=af.Collection(main_galaxies_pa), + extra_galaxies=af.Collection(extra_galaxies_pa), + scaling_galaxies=af.Collection(scaling_galaxies_pa), +) + +print("\n" + "=" * 80) +print("Section 2 — Standard PyAutoLens model (NFW halo, Isothermal galaxies)") +print("=" * 80) +print(f"Free parameters: {model_pyautolens.prior_count} (NFW 4 + 2 main + 1 extra + scaling 0 = 7)") +assert model_pyautolens.prior_count == 7 + +# --- Numerical bridge: the scaling-galaxy parameter change is exact --- +C_KM_S = 299792.458 +cosmology = ag.cosmo.FlatLambdaCDM(H0=H0, Om0=Om0) +d_s = cosmology.angular_diameter_distance_to_earth_in_kpc_from(redshift=redshift_source) +d_ls = cosmology.angular_diameter_distance_between_redshifts_in_kpc_from( + redshift_0=redshift_lens, redshift_1=redshift_source +) +grid = al.Grid2DIrregular([[1.5, 2.0], [-1.0, 0.8], [3.0, -2.5]]) + +fj = [ + (4.0 * np.pi * (s / C_KM_S) ** 2 * (d_ls / d_s) * (648000.0 / np.pi)) / s**2 + for s in [150.0, 250.0, 350.0] +] +print(f"(a) einstein_radius / sigma^2 = {fj[0]:.3e} (constant => einstein_radius ~ sigma^2, Faber-Jackson carries over)") +assert np.allclose(fj, fj[0], rtol=1e-12) + +theta = 1.4 +diff = np.max( + np.abs( + np.asarray(al.mp.IsothermalSph(centre=(0.0, 0.0), einstein_radius=theta).deflections_yx_2d_from(grid=grid)) + - np.asarray(al.mp.dPIEMassB0Sph(centre=(0.0, 0.0), ra=1e-6, rs=1e7, b0=theta).deflections_yx_2d_from(grid=grid)) + ) +) +print(f"(b) IsothermalSph(einstein_radius={theta}) vs dPIE b0 (SIS limit): max defl diff = {diff:.2e}") +assert diff < 1e-5 +print("Section 2 checks passed.") + + +""" +================================================================================ +__Section 3 — Keeping the Truncation (dPIE outskirts, Isothermal main + extra)__ +================================================================================ + +Section 2 dropped the members' truncation. Here is how to keep it, changing as little as possible: the +halo (NFW), the main galaxies and the extra galaxies (all Isothermal) are unchanged from Section 2 — +**only the scaling galaxies become truncated dPIEs**. + +**The difference truncation makes (all at large radius).** + + R (arcsec) SIE deflection dPIE (rs=20) deflection + 1 1.40 1.36 + 20 1.40 0.82 + 300 1.40 0.09 + +The SIE has ``rho ~ r^-2`` forever: deflection stays at ``einstein_radius``, mass diverges (``M ~ R``). +The dPIE falls as ``r^-4`` beyond ``rs``: deflection dies away, total mass is *finite*. Inside the arc +region (R ~ 1-5") they agree to a few percent — which is why the SIE is fine where the constraints live. + +**When MGLs do not need it.** ``rs`` models tidal *stripping* of a galaxy's outer halo by the host +potential. In a rich cluster stripping is strong and members are packed, so untruncated SIEs would each +contribute unbounded, overlapping mass — truncation is needed. In an MGL that is *not* a rich group the +potential is shallower, stripping weaker, and members are subdominant perturbers with the arcs near the +main lens — so the untruncated SIE is a defensible simplification. Keep truncation if the members are +numerous / massive enough to contribute meaningfully at large radius. + +**How, staying closest to the parameterization.** There is no separate "truncated isothermal" class — +the dPIE *is* it (with ``ra -> 0`` it is a Pseudo-Jaffe). So the scaling galaxies become +``dPIEMassB0Sph`` with a tiny core and a truncation ``rs``, keeping ``b0`` (the angular lens strength) as +the mass parameter. The main and extra galaxies stay Isothermal. + +**The linking across parameterizations.** Both the SIE ``einstein_radius`` and the dPIE ``b0`` are the +*angular lens strength* in arcsec, both ``~ sigma^2``, equal in the SIS limit. So a dPIE scaling +galaxy's ``b0`` anchors *directly* to the SIE BGC's ``einstein_radius`` with the identical relation — no +conversion: + + b0_i = einstein_radius_BGC * (L / L_BGC) ** 0.5, rs_i = rs_ref * (L / L_BGC) ** 0.5 + +The truncation ``rs`` is a *separate* parameter (the outer fall-off); it does not enter the anchoring. +Caveat: with finite ``rs`` the actual central deflection is a few percent below ``b0`` (truncation +removes outer mass) — negligible for subdominant members; you anchor on ``b0`` regardless. +""" + +# Halo (NFW), main and extra galaxies (SIE) — identical to Section 2. +halo_mass = af.Model(al.mp.NFW) +halo_mass.centre = (0.0, 0.0) +halo_mass.ell_comps.ell_comps_0 = af.UniformPrior(lower_limit=-0.5, upper_limit=0.5) # [FREE] +halo_mass.ell_comps.ell_comps_1 = af.UniformPrior(lower_limit=-0.5, upper_limit=0.5) # [FREE] +halo_mass.kappa_s = af.UniformPrior(lower_limit=0.05, upper_limit=0.5) # [FREE] +halo_mass.scale_radius = af.UniformPrior(lower_limit=10.0, upper_limit=60.0) # [FREE] +cluster_halo_trunc = af.Model(al.Galaxy, redshift=redshift_lens, mass=halo_mass) + +main_galaxies_trunc = [] +for centre in main_centres: + mass = af.Model(al.mp.IsothermalSph) + mass.centre = centre # [FIXED] + mass.einstein_radius = af.UniformPrior(lower_limit=0.0, upper_limit=3.0) # [FREE] + main_galaxies_trunc.append(af.Model(al.Galaxy, redshift=redshift_lens, mass=mass)) + +einstein_radius_bgc = main_galaxies_trunc[bgc_index].mass.einstein_radius # SIE Einstein radius = anchor + +extra_galaxies_trunc = [] +for centre, magnitude in zip(extra_centres, extra_magnitudes): + ratio = luminosity_ratio(magnitude) + upper = min(2.0 * (3.0 / 1.0**0.5) * ratio**0.5, 5.0) + mass = af.Model(al.mp.IsothermalSph) + mass.centre = centre # [FIXED] + mass.einstein_radius = af.UniformPrior(lower_limit=0.0, upper_limit=upper) # [FREE, BOUNDED] + extra_galaxies_trunc.append(af.Model(al.Galaxy, redshift=redshift_lens, mass=mass)) + +# Scaling galaxies -> dPIE (truncated), b0 anchored to the SIE BGC's einstein_radius. +rs_ref = 30.0 # [FIXED] reference truncation (arcsec); scales ~ L^0.5 like Lenstool's r_cut +scaling_galaxies_trunc = [] +for centre, magnitude in zip(scaling_centres, scaling_magnitudes): + ratio = luminosity_ratio(magnitude) + mass = af.Model(al.mp.dPIEMassB0Sph) + mass.centre = centre # [FIXED] + mass.ra = 0.01 # [FIXED] ~0: a pure truncated isothermal (no core) + mass.b0 = einstein_radius_bgc * ratio**0.5 # SIE einstein_radius drives dPIE b0 -> NO new parameter + mass.rs = rs_ref * ratio**0.5 # [FIXED] truncation, scaled with L^0.5 + scaling_galaxies_trunc.append(af.Model(al.Galaxy, redshift=redshift_lens, mass=mass)) + +model_truncated = af.Collection( + cluster_halo=cluster_halo_trunc, + main_galaxies=af.Collection(main_galaxies_trunc), + extra_galaxies=af.Collection(extra_galaxies_trunc), + scaling_galaxies=af.Collection(scaling_galaxies_trunc), +) + +print("\n" + "=" * 80) +print("Section 3 — Truncated scaling galaxies (dPIE outskirts, Isothermal main + extra)") +print("=" * 80) +print(f"Free parameters: {model_truncated.prior_count} (same 7 as Section 2; only the scaling profile changed)") +assert model_truncated.prior_count == 7 + +# (a) The cross-parameterization linking: a scaling galaxy's b0 = einstein_radius_BGC * (L/L_BGC)^0.5. +instance = model_truncated.instance_from_prior_medians() +er_bgc = instance.main_galaxies[bgc_index].mass.einstein_radius +print(f"(a) BGC einstein_radius = {er_bgc:.3f} arcsec (SIE); dPIE scaling galaxies b0 = einstein_radius_BGC * (L/L_BGC)^0.5:") +worst = 0.0 +for i, magnitude in enumerate(scaling_magnitudes): + ratio = luminosity_ratio(magnitude) + expected = er_bgc * ratio**0.5 + actual = instance.scaling_galaxies[i].mass.b0 + worst = max(worst, abs(actual - expected)) + print(f" scaling {i} (mag {magnitude}): b0 = {actual:.3f} arcsec = {ratio**0.5:.3f} x einstein_radius_BGC") +assert worst < 1e-9 + +# (b) The truncation bites: the deflection dies away at large radius (an SIE would stay flat). +member0 = instance.scaling_galaxies[0].mass +c = scaling_centres[0] +print("(b) scaling galaxy 0 deflection vs radius (truncated dPIE — dies away; an SIE would stay flat):") +for R in [1.0, 10.0, 50.0, 200.0]: + g = al.Grid2DIrregular([[c[0], c[1] + R]]) + a = float(np.max(np.abs(np.asarray(member0.deflections_yx_2d_from(grid=g))))) + print(f" R={R:>4.0f}\": alpha = {a:.4f} arcsec") +print("Section 3 checks passed.") diff --git a/workspace_index.json b/workspace_index.json index f8e4fe130..722fa3700 100644 --- a/workspace_index.json +++ b/workspace_index.json @@ -119,6 +119,29 @@ "summary": "This script provides a step-by-step guide of the cluster point-source ``log_likelihood_function``, the figure-of-merit Nautilus optimises when fitting a cluster lens model to ``point_datasets.csv``.", "title": "Log Likelihood Function: Cluster Point Source" }, + { + "contents": [ + "Model 1 \u2014 The Standard Lenstool Model", + "Model 2 \u2014 The Angular (b0) Parameterization", + "Model 3 \u2014 Mass Anchoring", + "Model 4 \u2014 The Two-Parameter Relation" + ], + "cross_refs": [], + "notebook": "notebooks/cluster/mass_parameterizations.ipynb", + "path": "scripts/cluster/mass_parameterizations.py", + "summary": "**A plain-language tour of the ways to parameterize and couple the masses in a Lenstool-style cluster model \u2014 the cluster-scale halo, the BCG and other individually-freed galaxies, and the member (\"scaling\") population.**", + "title": "Cluster & Group Mass Parameterizations" + }, + { + "contents": [], + "cross_refs": [ + "cluster/mass_parameterizations.py" + ], + "notebook": "notebooks/cluster/mass_parameterizations_pyautolens.ipynb", + "path": "scripts/cluster/mass_parameterizations_pyautolens.py", + "summary": "**For Lenstool users: how and why PyAutoLens re-expresses the standard Lenstool cluster / group model in its native parameterization for multi-galaxy lenses (MGLs), and how to keep the dPIE's truncation on the outskirts if you want it.**", + "title": "Mass Parameterizations II \u2014 Mapping Lenstool onto PyAutoLens" + }, { "contents": [ "Example",