From b508ffa81331f61b83f82fbe6cbb9748af292a71 Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Thu, 9 Jul 2026 11:17:15 +0100 Subject: [PATCH] Add tangential/cross shear profiles + Kaiser-Squires convergence map - autolens/weak/plot/shear_profile_plots.py: shear_tangential_cross_from / shear_profile_from (binned gamma_t with standard errors, gamma_x as the B-mode null test) + plot_shear_profile (WeakDataset or FitWeak, model overlay for fits); uses only the public .ellipticities/.phis accessors - autolens/weak/plot/convergence_plots.py: convergence_via_kaiser_squires_from (irregular->grid binning, FFT KS inversion, E/B-mode maps, dependency-free Gaussian smoothing) + plot_convergence_map - aplt re-exports; SIS-analytic unit tests (gamma_t = theta_E/2r, pure tangential null, KS peak at lens centre, E>>B) Step 6 of the weak-lensing series (PyAutoLabs/PyAutoLens#581). Co-Authored-By: Claude Fable 5 --- autolens/plot/__init__.py | 6 + autolens/weak/plot/convergence_plots.py | 218 ++++++++++++++++++++++ autolens/weak/plot/shear_profile_plots.py | 216 +++++++++++++++++++++ test_autolens/weak/test_convergence.py | 94 ++++++++++ test_autolens/weak/test_shear_profile.py | 86 +++++++++ 5 files changed, 620 insertions(+) create mode 100644 autolens/weak/plot/convergence_plots.py create mode 100644 autolens/weak/plot/shear_profile_plots.py create mode 100644 test_autolens/weak/test_convergence.py create mode 100644 test_autolens/weak/test_shear_profile.py diff --git a/autolens/plot/__init__.py b/autolens/plot/__init__.py index f67604d9f..c309ad091 100644 --- a/autolens/plot/__init__.py +++ b/autolens/plot/__init__.py @@ -85,6 +85,12 @@ plot_chi_squared_map, subplot_fit_weak, ) +from autolens.weak.plot.shear_profile_plots import ( + plot_shear_profile, +) +from autolens.weak.plot.convergence_plots import ( + plot_convergence_map, +) from autolens.lens.plot.subhalo_plots import ( subplot_detection_imaging, diff --git a/autolens/weak/plot/convergence_plots.py b/autolens/weak/plot/convergence_plots.py new file mode 100644 index 000000000..1a24375fa --- /dev/null +++ b/autolens/weak/plot/convergence_plots.py @@ -0,0 +1,218 @@ +""" +Kaiser-Squires convergence-map reconstruction from a shear catalogue. + +Kaiser & Squires (1993) showed the convergence :math:`\\kappa` and shear +:math:`\\gamma` are related algebraically in Fourier space: + +.. math:: + \\hat{\\kappa}(\\mathbf{k}) = D^*(\\mathbf{k}) \\, \\hat{\\gamma}(\\mathbf{k}), + \\qquad + D(\\mathbf{k}) = \\frac{k_x^2 - k_y^2 + 2 i k_x k_y}{k_x^2 + k_y^2} + +so a mass map can be reconstructed directly from the measured shear field with two +FFTs — no mass model assumed. This is the classic "dark matter map" technique used +for merging clusters (e.g. the Bullet cluster) and survey mass maps. + +Because the catalogue's galaxy positions are irregular, the shears are first binned +onto a regular grid (plain per-cell mean, empty cells zero) over the catalogue's +extent, then optionally smoothed with a Gaussian kernel — standard practice, since +the inversion is only defined on a grid and raw per-cell shears are shape-noise +dominated. The imaginary part of the reconstruction (the B-mode map) carries no +lensing signal and is returned alongside the E-mode map as a systematics check. + +The reconstruction inherits Kaiser-Squires' well-known caveats: the +:math:`\\mathbf{k} = 0` mode is unconstrained (the mass-sheet degeneracy — maps are +zero-mean by construction) and FFT periodicity causes edge artefacts on small +fields. For quantitative masses, fit a mass model (``scripts/weak/modeling.py``); +the map is a visualization and model-independent cross-check. +""" +from typing import Optional, Tuple + +import numpy as np + +import autoarray as aa + +from autoarray.plot.utils import save_figure, subplots + +from autolens.weak.plot.weak_dataset_plots import _positions_yx +from autolens.weak.plot.shear_profile_plots import _gamma_1_2_from + + +def convergence_via_kaiser_squires_from( + shear_yx, + shape_native: Tuple[int, int] = (50, 50), + smoothing_sigma_pixels: float = 1.0, + extent: Optional[Tuple[float, float, float, float]] = None, +) -> Tuple[aa.Array2D, aa.Array2D]: + """ + Reconstruct E-mode (convergence) and B-mode maps from a shear field. + + Parameters + ---------- + shear_yx + The shear field (e.g. ``dataset.shear_yx``). + shape_native + The ``(rows, cols)`` shape of the regular grid the shears are binned onto. + smoothing_sigma_pixels + The sigma (in pixels) of the Gaussian kernel applied to the binned shear + maps before inversion; ``0.0`` disables smoothing. + extent + The ``(x_min, x_max, y_min, y_max)`` field extent; defaults to the + catalogue's bounding box (with a small buffer, matching + ``WeakDataset.extent_from``). + + Returns + ------- + An ``(e_mode, b_mode)`` pair of ``aa.Array2D`` maps: the convergence + reconstruction and its B-mode systematics check. + """ + positions = _positions_yx(shear_yx) + gamma_1, gamma_2 = _gamma_1_2_from(shear_yx) + + if extent is None: + buffer = 0.1 + x_min, x_max = positions[:, 1].min() - buffer, positions[:, 1].max() + buffer + y_min, y_max = positions[:, 0].min() - buffer, positions[:, 0].max() + buffer + else: + x_min, x_max, y_min, y_max = extent + + rows, cols = shape_native + + # Bin the irregular shears onto the grid: plain mean per cell, zero where empty. + y_edges = np.linspace(y_min, y_max, rows + 1) + x_edges = np.linspace(x_min, x_max, cols + 1) + + counts, _, _ = np.histogram2d( + positions[:, 0], positions[:, 1], bins=[y_edges, x_edges] + ) + sum_1, _, _ = np.histogram2d( + positions[:, 0], positions[:, 1], bins=[y_edges, x_edges], weights=gamma_1 + ) + sum_2, _, _ = np.histogram2d( + positions[:, 0], positions[:, 1], bins=[y_edges, x_edges], weights=gamma_2 + ) + + with np.errstate(invalid="ignore"): + grid_1 = np.where(counts > 0, sum_1 / np.maximum(counts, 1), 0.0) + grid_2 = np.where(counts > 0, sum_2 / np.maximum(counts, 1), 0.0) + + if smoothing_sigma_pixels > 0.0: + kernel = _gaussian_kernel_2d(sigma=smoothing_sigma_pixels) + grid_1 = _convolve_2d_same(grid_1, kernel) + grid_2 = _convolve_2d_same(grid_2, kernel) + + # Kaiser-Squires inversion: kappa_hat = D*(k) gamma_hat, D = (kx^2 - ky^2 + 2i kx ky) / k^2. + ky = np.fft.fftfreq(rows)[:, None] + kx = np.fft.fftfreq(cols)[None, :] + k_squared = kx**2.0 + ky**2.0 + k_squared[0, 0] = 1.0 # k = 0 mode is unconstrained (mass-sheet degeneracy); set below. + + d_conj = ((kx**2.0 - ky**2.0) - 2.0j * kx * ky) / k_squared + + gamma_hat = np.fft.fft2(grid_1) + 1.0j * np.fft.fft2(grid_2) + kappa_hat = d_conj * gamma_hat + kappa_hat[0, 0] = 0.0 + + kappa = np.fft.ifft2(kappa_hat) + + pixel_scales = ((y_max - y_min) / rows, (x_max - x_min) / cols) + + # Row 0 of the binned grids is y_min; autoarray native arrays put y_max at row 0. + e_mode = aa.Array2D.no_mask( + values=np.flipud(kappa.real), pixel_scales=pixel_scales + ) + b_mode = aa.Array2D.no_mask( + values=np.flipud(kappa.imag), pixel_scales=pixel_scales + ) + + return e_mode, b_mode + + +def _gaussian_kernel_2d(sigma: float) -> np.ndarray: + """A normalised 2D Gaussian kernel truncated at 3 sigma.""" + half_width = max(int(np.ceil(3.0 * sigma)), 1) + x = np.arange(-half_width, half_width + 1) + kernel_1d = np.exp(-0.5 * (x / sigma) ** 2.0) + kernel = np.outer(kernel_1d, kernel_1d) + return kernel / kernel.sum() + + +def _convolve_2d_same(image: np.ndarray, kernel: np.ndarray) -> np.ndarray: + """Same-size 2D convolution via zero-padded FFT (avoids a scipy dependency).""" + kh, kw = kernel.shape + ih, iw = image.shape + padded = np.zeros((ih + kh - 1, iw + kw - 1)) + padded[:ih, :iw] = image + kernel_padded = np.zeros_like(padded) + kernel_padded[:kh, :kw] = kernel + result = np.fft.ifft2(np.fft.fft2(padded) * np.fft.fft2(kernel_padded)).real + h0, w0 = (kh - 1) // 2, (kw - 1) // 2 + return result[h0 : h0 + ih, w0 : w0 + iw] + + +def plot_convergence_map( + shear_yx, + shape_native: Tuple[int, int] = (50, 50), + smoothing_sigma_pixels: float = 1.0, + show_positions: bool = True, + ax=None, + title: str = "Kaiser-Squires Convergence", + output_path: Optional[str] = None, + output_filename: str = "convergence_map", + output_format: Optional[str] = None, +): + """ + Plot the Kaiser-Squires E-mode convergence reconstruction of a shear field. + + Parameters + ---------- + shear_yx + The shear field (e.g. ``dataset.shear_yx``). + shape_native + The ``(rows, cols)`` reconstruction grid shape. + smoothing_sigma_pixels + Gaussian smoothing applied to the binned shears before inversion. + show_positions + Whether to overlay the catalogue's galaxy positions on the map. + ax + An existing matplotlib axes to draw on; a new figure is created if ``None``. + """ + e_mode, _ = convergence_via_kaiser_squires_from( + shear_yx=shear_yx, + shape_native=shape_native, + smoothing_sigma_pixels=smoothing_sigma_pixels, + ) + + positions = _positions_yx(shear_yx) + extent = [ + positions[:, 1].min() - 0.1, + positions[:, 1].max() + 0.1, + positions[:, 0].min() - 0.1, + positions[:, 0].max() + 0.1, + ] + + fig = None + if ax is None: + fig, ax = subplots(1, 1) + + image = ax.imshow( + e_mode.native, + origin="upper", + extent=extent, + cmap="magma", + ) + if show_positions: + ax.scatter( + positions[:, 1], positions[:, 0], s=2, c="cyan", alpha=0.5, linewidths=0 + ) + ax.set_xlabel("x (arcsec)") + ax.set_ylabel("y (arcsec)") + ax.set_title(title) + if fig is not None: + fig.colorbar(image, ax=ax, label=r"$\kappa$ (E-mode)") + save_figure( + fig, + path=output_path, + filename=output_filename, + format=output_format, + ) diff --git a/autolens/weak/plot/shear_profile_plots.py b/autolens/weak/plot/shear_profile_plots.py new file mode 100644 index 000000000..0a7e597d2 --- /dev/null +++ b/autolens/weak/plot/shear_profile_plots.py @@ -0,0 +1,216 @@ +""" +Tangential / cross shear radial profiles for a ``WeakDataset`` or ``FitWeak``. + +The azimuthally averaged tangential shear profile :math:`\\gamma_t(r)` is the standard +observable of cluster weak lensing (e.g. Oguri et al. 2012's SGAS analysis; Medezinski +et al. 2016's A2744 analysis): a foreground mass distribution stretches background +galaxies tangentially around it, so binning the tangential component of the shear in +radius traces the projected mass profile. The cross component :math:`\\gamma_x` (the +45-degree rotated component) is not produced by gravitational lensing at leading order, +so its binned profile doubles as the standard B-mode null test — a systematic that +leaks power into :math:`\\gamma_x` would contaminate :math:`\\gamma_t` too. + +Component conventions follow the rest of ``autolens/weak``: shears are accessed +exclusively through the public ``.ellipticities`` (:math:`|\\gamma|`) and ``.phis`` +(position angle in **degrees**) accessors, from which +:math:`\\gamma_1 = |\\gamma| \\cos 2\\phi`, :math:`\\gamma_2 = |\\gamma| \\sin 2\\phi` +are reconstructed — the raw ``[gamma_2, gamma_1]`` storage is never indexed directly. + +With :math:`\\gamma = \\gamma_1 + i\\gamma_2` and :math:`\\varphi` the position angle of +a galaxy about the profile centre: + +.. math:: + \\gamma_t = -\\mathrm{Re}[\\gamma e^{-2i\\varphi}], \\qquad + \\gamma_x = -\\mathrm{Im}[\\gamma e^{-2i\\varphi}] + +so a tangentially aligned shear field has :math:`\\gamma_t > 0`. +""" +from typing import Optional, Tuple + +import numpy as np + +from autoarray.plot.utils import save_figure, subplots + +from autolens.weak.plot.weak_dataset_plots import _positions_yx + + +def _gamma_1_2_from(shear_yx) -> Tuple[np.ndarray, np.ndarray]: + """Reconstruct ``(gamma_1, gamma_2)`` from the public accessors.""" + magnitudes = np.asarray(shear_yx.ellipticities) + phis_rad = np.deg2rad(np.asarray(shear_yx.phis)) + return magnitudes * np.cos(2.0 * phis_rad), magnitudes * np.sin(2.0 * phis_rad) + + +def shear_tangential_cross_from( + shear_yx, + centre: Tuple[float, float] = (0.0, 0.0), +) -> Tuple[np.ndarray, np.ndarray, np.ndarray]: + """ + The per-galaxy tangential and cross shear components about ``centre``. + + Parameters + ---------- + shear_yx + The shear field (e.g. ``dataset.shear_yx`` or ``fit.model_shear``). + centre + The ``(y, x)`` centre about which tangential / cross components are defined, + e.g. the lens-galaxy centre. + + Returns + ------- + A ``(gamma_t, gamma_x, radii)`` tuple of ``(N,)`` arrays, where ``radii`` are the + galaxy distances from ``centre``. + """ + positions = _positions_yx(shear_yx) + dy = positions[:, 0] - centre[0] + dx = positions[:, 1] - centre[1] + + radii = np.sqrt(dy**2.0 + dx**2.0) + varphi = np.arctan2(dy, dx) + + gamma_1, gamma_2 = _gamma_1_2_from(shear_yx) + + cos_2varphi = np.cos(2.0 * varphi) + sin_2varphi = np.sin(2.0 * varphi) + + gamma_t = -(gamma_1 * cos_2varphi + gamma_2 * sin_2varphi) + gamma_x = -(gamma_2 * cos_2varphi - gamma_1 * sin_2varphi) + + return gamma_t, gamma_x, radii + + +def shear_profile_from( + shear_yx, + centre: Tuple[float, float] = (0.0, 0.0), + bins: int = 10, + noise_map=None, +) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray]: + """ + Azimuthally averaged tangential / cross shear profiles in radial bins. + + Bin edges are linearly spaced between the innermost and outermost galaxy radius. + Each bin's value is the plain mean of its galaxies' components and its error the + standard error on that mean — when ``noise_map`` is given, the per-galaxy + measurement noise replaces the empirical scatter in bins with fewer than two + galaxies (where a standard error is undefined). + + Returns + ------- + A ``(bin_radii, gamma_t, gamma_t_err, gamma_x, gamma_x_err)`` tuple of ``(bins,)`` + arrays; bins containing no galaxies hold ``np.nan``. + """ + gamma_t, gamma_x, radii = shear_tangential_cross_from( + shear_yx=shear_yx, centre=centre + ) + + edges = np.linspace(radii.min(), radii.max(), bins + 1) + bin_radii = 0.5 * (edges[:-1] + edges[1:]) + + profile_t = np.full(bins, np.nan) + error_t = np.full(bins, np.nan) + profile_x = np.full(bins, np.nan) + error_x = np.full(bins, np.nan) + + noise = np.asarray(noise_map) if noise_map is not None else None + + for i in range(bins): + # Right-inclusive final bin so the outermost galaxy is not dropped. + if i == bins - 1: + in_bin = (radii >= edges[i]) & (radii <= edges[i + 1]) + else: + in_bin = (radii >= edges[i]) & (radii < edges[i + 1]) + n = in_bin.sum() + if n == 0: + continue + profile_t[i] = gamma_t[in_bin].mean() + profile_x[i] = gamma_x[in_bin].mean() + if n > 1: + error_t[i] = gamma_t[in_bin].std(ddof=1) / np.sqrt(n) + error_x[i] = gamma_x[in_bin].std(ddof=1) / np.sqrt(n) + elif noise is not None: + error_t[i] = error_x[i] = noise[in_bin][0] + + return bin_radii, profile_t, error_t, profile_x, error_x + + +def plot_shear_profile( + obj, + centre: Tuple[float, float] = (0.0, 0.0), + bins: int = 10, + ax=None, + title: str = "Shear Profile", + output_path: Optional[str] = None, + output_filename: str = "shear_profile", + output_format: Optional[str] = None, +): + """ + Plot binned tangential and cross shear profiles of a ``WeakDataset`` or ``FitWeak``. + + The tangential profile is drawn with error bars; the cross profile (offset points) + should scatter around zero — the B-mode null test. When ``obj`` is a ``FitWeak`` + the model shear's tangential profile is overlaid as a line, so data and model can + be compared in the space where cluster weak-lensing results are usually shown. + + Parameters + ---------- + obj + A ``WeakDataset`` or ``FitWeak``. + centre + The ``(y, x)`` centre about which the profiles are computed. + bins + The number of linearly spaced radial bins. + ax + An existing matplotlib axes to draw on; a new figure is created if ``None``. + """ + dataset = obj.dataset if hasattr(obj, "model_shear") else obj + + bin_radii, gamma_t, gamma_t_err, gamma_x, gamma_x_err = shear_profile_from( + shear_yx=dataset.shear_yx, + centre=centre, + bins=bins, + noise_map=dataset.noise_map, + ) + + fig = None + if ax is None: + fig, ax = subplots(1, 1) + + ax.errorbar( + bin_radii, + gamma_t, + yerr=gamma_t_err, + fmt="o", + color="k", + capsize=3, + label=r"$\gamma_t$ (data)", + ) + ax.errorbar( + bin_radii, + gamma_x, + yerr=gamma_x_err, + fmt="s", + color="tab:blue", + alpha=0.6, + capsize=3, + label=r"$\gamma_\times$ (B-mode null test)", + ) + + if hasattr(obj, "model_shear"): + _, model_t, _, _, _ = shear_profile_from( + shear_yx=obj.model_shear, centre=centre, bins=bins + ) + ax.plot(bin_radii, model_t, color="r", alpha=0.8, label=r"$\gamma_t$ (model)") + + ax.axhline(0.0, color="0.6", lw=0.8, zorder=0) + ax.set_xlabel("Radius from centre (arcsec)") + ax.set_ylabel("Shear") + ax.set_title(title) + ax.legend() + + if fig is not None: + save_figure( + fig, + path=output_path, + filename=output_filename, + format=output_format, + ) diff --git a/test_autolens/weak/test_convergence.py b/test_autolens/weak/test_convergence.py new file mode 100644 index 000000000..34d133465 --- /dev/null +++ b/test_autolens/weak/test_convergence.py @@ -0,0 +1,94 @@ +import numpy as np + +import autolens as al + +from autolens.weak.plot.convergence_plots import convergence_via_kaiser_squires_from + + +def _sis_dataset(n_galaxies=2000, einstein_radius=1.6, seed=4): + lens = al.Galaxy( + redshift=0.5, + mass=al.mp.IsothermalSph(centre=(0.0, 0.0), einstein_radius=einstein_radius), + ) + tracer = al.Tracer(galaxies=[lens, al.Galaxy(redshift=1.0)]) + return al.SimulatorShearYX( + noise_sigma=0.0, seed=seed + ).via_tracer_random_positions_from( + tracer=tracer, n_galaxies=n_galaxies, grid_extent=3.0 + ) + + +def test__kaiser_squires__peak_at_lens_centre_and_weak_b_mode(): + dataset = _sis_dataset() + + e_mode, b_mode = convergence_via_kaiser_squires_from( + shear_yx=dataset.shear_yx, + shape_native=(32, 32), + smoothing_sigma_pixels=1.0, + ) + + e = np.asarray(e_mode.native) + b = np.asarray(b_mode.native) + + # The reconstruction must peak at the lens centre — the central pixels of the map. + peak_row, peak_col = np.unravel_index(np.argmax(e), e.shape) + centre = (e.shape[0] - 1) / 2.0 + assert abs(peak_row - centre) <= 2.0 + assert abs(peak_col - centre) <= 2.0 + + # E-mode signal dominates the B-mode systematics map for a pure lensing field. + inner = np.s_[8:24, 8:24] # interior region, away from FFT edge artefacts + assert np.abs(e[inner]).max() > 5.0 * np.abs(b[inner]).max() + + +def test__kaiser_squires__map_shapes_and_zero_mean(): + dataset = _sis_dataset(n_galaxies=500) + + e_mode, b_mode = convergence_via_kaiser_squires_from( + shear_yx=dataset.shear_yx, shape_native=(20, 24), smoothing_sigma_pixels=0.0 + ) + + assert e_mode.shape_native == (20, 24) + assert b_mode.shape_native == (20, 24) + + # The k = 0 mode is zeroed (mass-sheet degeneracy), so maps are zero-mean. + assert abs(np.asarray(e_mode.native).mean()) < 1e-12 + + +def test__kaiser_squires__elliptical_lens_orientation_parity(): + """ + The SIS tests are rotation-invariant, so they cannot catch an axis-parity (mirror) + error in the inversion. An elliptical lens breaks that degeneracy: the KS map's + major axis must align with the lens mass's major axis (45 deg), not its mirror. + """ + lens = al.Galaxy( + redshift=0.5, + mass=al.mp.Isothermal( + centre=(0.0, 0.0), + einstein_radius=1.6, + ell_comps=al.convert.ell_comps_from(axis_ratio=0.5, angle=45.0), + ), + ) + tracer = al.Tracer(galaxies=[lens, al.Galaxy(redshift=1.0)]) + dataset = al.SimulatorShearYX( + noise_sigma=0.0, seed=1 + ).via_tracer_random_positions_from(tracer=tracer, n_galaxies=5000, grid_extent=3.0) + + e_mode, _ = convergence_via_kaiser_squires_from( + shear_yx=dataset.shear_yx, shape_native=(64, 64), smoothing_sigma_pixels=1.0 + ) + e = np.asarray(e_mode.native) + + ys, xs = np.mgrid[0:64, 0:64] + y_coord = 3.0 - (ys + 0.5) * (6.0 / 64) + x_coord = -3.0 + (xs + 0.5) * (6.0 / 64) + weights = np.clip(np.where(e > np.percentile(e, 95), e, 0.0), 0.0, None) + x_bar = (weights * x_coord).sum() / weights.sum() + y_bar = (weights * y_coord).sum() / weights.sum() + q_xx = (weights * (x_coord - x_bar) ** 2.0).sum() / weights.sum() + q_yy = (weights * (y_coord - y_bar) ** 2.0).sum() / weights.sum() + q_xy = (weights * (x_coord - x_bar) * (y_coord - y_bar)).sum() / weights.sum() + + angle = 0.5 * np.degrees(np.arctan2(2.0 * q_xy, q_xx - q_yy)) + + assert abs(angle - 45.0) < 10.0 diff --git a/test_autolens/weak/test_shear_profile.py b/test_autolens/weak/test_shear_profile.py new file mode 100644 index 000000000..ea8f88528 --- /dev/null +++ b/test_autolens/weak/test_shear_profile.py @@ -0,0 +1,86 @@ +import numpy as np + +import autoarray as aa +import autolens as al + +from autolens.weak.plot.shear_profile_plots import ( + shear_profile_from, + shear_tangential_cross_from, +) + + +def _sis_dataset(n_positions=400, noise_sigma=0.0, einstein_radius=1.6, seed=2): + lens = al.Galaxy( + redshift=0.5, + mass=al.mp.IsothermalSph(centre=(0.0, 0.0), einstein_radius=einstein_radius), + ) + source = al.Galaxy(redshift=1.0) + tracer = al.Tracer(galaxies=[lens, source]) + simulator = al.SimulatorShearYX(noise_sigma=noise_sigma, seed=seed) + return simulator.via_tracer_random_positions_from( + tracer=tracer, n_galaxies=n_positions, grid_extent=3.0 + ) + + +def test__tangential_shear__sis_is_pure_tangential(): + """ + An SIS shear field is purely tangential about the lens centre, so the tangential + component must equal the shear magnitude and the cross component must vanish — + this validates the sign and angle conventions end to end against the simulator. + """ + dataset = _sis_dataset() + + gamma_t, gamma_x, radii = shear_tangential_cross_from( + shear_yx=dataset.shear_yx, centre=(0.0, 0.0) + ) + + np.testing.assert_allclose( + gamma_t, np.asarray(dataset.shear_yx.ellipticities), atol=1e-6 + ) + np.testing.assert_allclose(gamma_x, 0.0, atol=1e-6) + assert (radii > 0.0).all() + + +def test__tangential_profile__matches_sis_analytic(): + """Binned gamma_t must track the SIS analytic profile theta_E / 2r, noise-free.""" + einstein_radius = 1.6 + dataset = _sis_dataset(einstein_radius=einstein_radius) + + bin_radii, gamma_t, _, gamma_x, _ = shear_profile_from( + shear_yx=dataset.shear_yx, centre=(0.0, 0.0), bins=8 + ) + + valid = ~np.isnan(gamma_t) + assert valid.sum() >= 6 + + analytic = einstein_radius / (2.0 * bin_radii[valid]) + + # Binned means of 1/r inside finite bins deviate from the bin-centre value; + # 15% tolerance comfortably covers that discretisation while catching sign or + # convention errors (which are order-unity). + np.testing.assert_allclose(gamma_t[valid], analytic, rtol=0.15) + np.testing.assert_allclose(gamma_x[valid], 0.0, atol=1e-6) + + +def test__profile__off_centre_recentres(): + """Recomputing about the true (offset) lens centre must stay pure-tangential.""" + lens = al.Galaxy( + redshift=0.5, + mass=al.mp.IsothermalSph(centre=(0.5, -0.3), einstein_radius=1.0), + ) + tracer = al.Tracer(galaxies=[lens, al.Galaxy(redshift=1.0)]) + dataset = al.SimulatorShearYX(noise_sigma=0.0, seed=3).via_tracer_random_positions_from( + tracer=tracer, n_galaxies=200, grid_extent=3.0 + ) + + _, gamma_x_true_centre, _ = shear_tangential_cross_from( + shear_yx=dataset.shear_yx, centre=(0.5, -0.3) + ) + _, gamma_x_wrong_centre, _ = shear_tangential_cross_from( + shear_yx=dataset.shear_yx, centre=(0.0, 0.0) + ) + + # One random galaxy lands ~on the centre, where the hessian finite-difference + # error grows like 1/r — hence the looser tolerance than the on-centre test. + np.testing.assert_allclose(gamma_x_true_centre, 0.0, atol=1e-5) + assert np.abs(gamma_x_wrong_centre).max() > 1e-3