From 09be29c0a6b5a8b290e6ab22c1c79ab6a325157d Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Fri, 13 Feb 2026 10:47:01 +0000 Subject: [PATCH 01/11] added adaptive matern kernel --- .../matern_adaptive_brightness_kernel.py | 52 +++++++++++++++++++ 1 file changed, 52 insertions(+) create mode 100644 autoarray/inversion/regularization/matern_adaptive_brightness_kernel.py diff --git a/autoarray/inversion/regularization/matern_adaptive_brightness_kernel.py b/autoarray/inversion/regularization/matern_adaptive_brightness_kernel.py new file mode 100644 index 000000000..96d99c4aa --- /dev/null +++ b/autoarray/inversion/regularization/matern_adaptive_brightness_kernel.py @@ -0,0 +1,52 @@ +from autoarray.inversion.regularization.abstract import AbstractRegularization + +from autoarray.inversion.regularization.matern_kernel import matern_cov_matrix_from + + +class AdaptiveBrightnessMatern(AbstractRegularization): + def __init__( + self, + coefficient: float = 1.0, + scale: float = 1.0, + nu: float = 0.5, + rho: float = 1.0, + ): + super().__init__(coefficient=coefficient, scale=scale, rho=rho) + self.nu = nu + + def covariance_kernel_weights_from(self, linear_obj: LinearObj) -> np.ndarray: + pixel_signals = linear_obj.pixel_signals_from(signal_scale=1.0) + return np.exp(-self.rho * (1 - pixel_signals / pixel_signals.max())) + + def regularization_matrix_from(self, linear_obj: LinearObj) -> np.ndarray: + kernel_weights = self.covariance_kernel_weights_from(linear_obj=linear_obj) + + covariance_matrix = matern_cov_matrix_from( + scale=self.scale, + pixel_points=linear_obj.source_plane_mesh_grid, + nu=self.nu, + weights=kernel_weights, + ) + + return self.coefficient * np.linalg.inv(covariance_matrix) + + def regularization_weights_from(self, linear_obj: LinearObj) -> np.ndarray: + """ + Returns the regularization weights of this regularization scheme. + + The regularization weights define the level of regularization applied to each parameter in the linear object + (e.g. the ``pixels`` in a ``Mapper``). + + For standard regularization (e.g. ``Constant``) are weights are equal, however for adaptive schemes + (e.g. ``AdaptiveBrightness``) they vary to adapt to the data being reconstructed. + + Parameters + ---------- + linear_obj + The linear object (e.g. a ``Mapper``) which uses these weights when performing regularization. + + Returns + ------- + The regularization weights. + """ + return 1.0/self.covariance_kernel_weights_from(linear_obj=linear_obj) #meaningless, but consistent with other regularization schemes \ No newline at end of file From ac9805b528fc7e6a9a788bf0ea58b5636024e035 Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Fri, 13 Feb 2026 11:03:25 +0000 Subject: [PATCH 02/11] unit tests added with updates for JAX and correct functionality --- .../inversion/regularization/__init__.py | 1 + .../matern_adaptive_brightness_kernel.py | 164 ++++++++++++++---- autoarray/structures/arrays/kernel_2d.py | 2 +- .../test_matern_adaptive_brightness_kernel.py | 55 ++++++ 4 files changed, 187 insertions(+), 35 deletions(-) create mode 100644 test_autoarray/inversion/regularizations/test_matern_adaptive_brightness_kernel.py diff --git a/autoarray/inversion/regularization/__init__.py b/autoarray/inversion/regularization/__init__.py index e34d07b6b..c5d696062 100644 --- a/autoarray/inversion/regularization/__init__.py +++ b/autoarray/inversion/regularization/__init__.py @@ -10,3 +10,4 @@ from .gaussian_kernel import GaussianKernel from .exponential_kernel import ExponentialKernel from .matern_kernel import MaternKernel +from .matern_adaptive_brightness_kernel import MaternAdaptiveBrightnessKernel diff --git a/autoarray/inversion/regularization/matern_adaptive_brightness_kernel.py b/autoarray/inversion/regularization/matern_adaptive_brightness_kernel.py index 96d99c4aa..b629a10a0 100644 --- a/autoarray/inversion/regularization/matern_adaptive_brightness_kernel.py +++ b/autoarray/inversion/regularization/matern_adaptive_brightness_kernel.py @@ -1,52 +1,148 @@ -from autoarray.inversion.regularization.abstract import AbstractRegularization +from __future__ import annotations +import numpy as np +from typing import TYPE_CHECKING -from autoarray.inversion.regularization.matern_kernel import matern_cov_matrix_from +from autoarray.inversion.regularization.matern_kernel import MaternKernel +if TYPE_CHECKING: + from autoarray.inversion.linear_obj.linear_obj import LinearObj -class AdaptiveBrightnessMatern(AbstractRegularization): +from autoarray.inversion.regularization.matern_kernel import matern_kernel + +def matern_cov_matrix_from( + scale: float, + nu: float, + pixel_points, + weights=None, + xp=np, +): + """ + Construct the regularization covariance matrix (N x N) using a Matérn kernel, + optionally modulated by per-pixel weights. + + If `weights` is provided (shape [N]), the covariance is: + C_ij = K(d_ij; scale, nu) * w_i * w_j + with a small diagonal jitter added for numerical stability. + + Parameters + ---------- + scale + Typical correlation length of the Matérn kernel. + nu + Smoothness parameter of the Matérn kernel. + pixel_points + Array-like of shape [N, 2] with (y, x) coordinates (or any 2D coords; only distances matter). + weights + Optional array-like of shape [N]. If None, treated as all ones. + xp + Backend (numpy or jax.numpy). + + Returns + ------- + covariance_matrix + Array of shape [N, N]. + """ + + # -------------------------------- + # Pairwise distances (broadcasted) + # -------------------------------- + diff = pixel_points[:, None, :] - pixel_points[None, :, :] # (N, N, 2) + d_ij = xp.sqrt(diff[..., 0] ** 2 + diff[..., 1] ** 2) # (N, N) + + # -------------------------------- + # Base Matérn covariance + # -------------------------------- + covariance_matrix = matern_kernel(d_ij, l=scale, v=nu, xp=xp) # (N, N) + + # -------------------------------- + # Apply weights: C_ij *= w_i * w_j + # (broadcasted outer product, JAX-safe) + # -------------------------------- + if weights is not None: + w = xp.asarray(weights) + # Ensure shape (N,) -> outer product (N,1)*(1,N) -> (N,N) + covariance_matrix = covariance_matrix * (w[:, None] * w[None, :]) + + # -------------------------------- + # Add diagonal jitter (JAX-safe) + # -------------------------------- + pixels = pixel_points.shape[0] + covariance_matrix = covariance_matrix + 1e-8 * xp.eye(pixels) + + return covariance_matrix + + +class MaternAdaptiveBrightnessKernel(MaternKernel): def __init__( - self, - coefficient: float = 1.0, - scale: float = 1.0, - nu: float = 0.5, - rho: float = 1.0, + self, + coefficient: float = 1.0, + scale: float = 1.0, + nu: float = 0.5, + rho: float = 1.0, ): - super().__init__(coefficient=coefficient, scale=scale, rho=rho) - self.nu = nu + """ + Regularization which uses a Matern smoothing kernel to regularize the solution with regularization weights + that adapt to the brightness of the source being reconstructed. + + For this regularization scheme, every pixel is regularized with every other pixel. This contrasts many other + schemes, where regularization is based on neighboring (e.g. do the pixels share a Delaunay edge?) or computing + derivates around the center of the pixel (where nearby pixels are regularization locally in similar ways). + + This makes the regularization matrix fully dense and therefore maybe change the run times of the solution. + It also leads to more overall smoothing which can lead to more stable linear inversions. + + For the weighted regularization scheme, each pixel is given an 'effective regularization weight', which is + applied when each set of pixel neighbors are regularized with one another. The motivation of this is that + different regions of a pixelization's mesh require different levels of regularization (e.g., high smoothing where the + no signal is present and less smoothing where it is, see (Nightingale, Dye and Massey 2018)). + + This scheme is not used by Vernardos et al. (2022): https://arxiv.org/abs/2202.09378, but it follows + a similar approach. - def covariance_kernel_weights_from(self, linear_obj: LinearObj) -> np.ndarray: - pixel_signals = linear_obj.pixel_signals_from(signal_scale=1.0) - return np.exp(-self.rho * (1 - pixel_signals / pixel_signals.max())) + A full description of regularization and this matrix can be found in the parent `AbstractRegularization` class. - def regularization_matrix_from(self, linear_obj: LinearObj) -> np.ndarray: - kernel_weights = self.covariance_kernel_weights_from(linear_obj=linear_obj) + Parameters + ---------- + coefficient + The regularization coefficient which controls the degree of smooth of the inversion reconstruction. + scale + The typical scale of the exponential regularization pattern. + nu + Controls the derivative of the regularization pattern (`nu=0.5` is a Gaussian). + """ + super().__init__(coefficient=coefficient, scale=scale, nu=nu) + self.rho = rho + + def covariance_kernel_weights_from(self, linear_obj: LinearObj, xp=np) -> np.ndarray: + """ + Returns per-pixel kernel weights that adapt to the reconstructed pixel brightness. + """ + # Assumes linear_obj.pixel_signals_from is xp-aware elsewhere in the codebase. + pixel_signals = linear_obj.pixel_signals_from(signal_scale=1.0, xp=xp) + + max_signal = xp.max(pixel_signals) + max_signal = xp.maximum(max_signal, 1e-8) # avoid divide-by-zero (JAX-safe) + + return xp.exp(-self.rho * (1.0 - pixel_signals / max_signal)) + + def regularization_matrix_from(self, linear_obj: LinearObj, xp=np) -> np.ndarray: + kernel_weights = self.covariance_kernel_weights_from(linear_obj=linear_obj, xp=xp) + + # Follow the xp pattern used in the Matérn kernel module (often `.array` for grids). + pixel_points = linear_obj.source_plane_mesh_grid covariance_matrix = matern_cov_matrix_from( scale=self.scale, - pixel_points=linear_obj.source_plane_mesh_grid, + pixel_points=pixel_points, nu=self.nu, weights=kernel_weights, + xp=xp, ) - return self.coefficient * np.linalg.inv(covariance_matrix) + return self.coefficient * xp.linalg.inv(covariance_matrix) - def regularization_weights_from(self, linear_obj: LinearObj) -> np.ndarray: + def regularization_weights_from(self, linear_obj: LinearObj, xp=np) -> np.ndarray: """ Returns the regularization weights of this regularization scheme. - - The regularization weights define the level of regularization applied to each parameter in the linear object - (e.g. the ``pixels`` in a ``Mapper``). - - For standard regularization (e.g. ``Constant``) are weights are equal, however for adaptive schemes - (e.g. ``AdaptiveBrightness``) they vary to adapt to the data being reconstructed. - - Parameters - ---------- - linear_obj - The linear object (e.g. a ``Mapper``) which uses these weights when performing regularization. - - Returns - ------- - The regularization weights. """ - return 1.0/self.covariance_kernel_weights_from(linear_obj=linear_obj) #meaningless, but consistent with other regularization schemes \ No newline at end of file + return 1.0 / self.covariance_kernel_weights_from(linear_obj=linear_obj, xp=xp) \ No newline at end of file diff --git a/autoarray/structures/arrays/kernel_2d.py b/autoarray/structures/arrays/kernel_2d.py index 51c37975e..11aeaf927 100644 --- a/autoarray/structures/arrays/kernel_2d.py +++ b/autoarray/structures/arrays/kernel_2d.py @@ -612,7 +612,7 @@ def convolved_image_from( image, blurring_image, jax_method="direct", - use_mixed_precision : bool = False, + use_mixed_precision: bool = False, xp=np, ): """ diff --git a/test_autoarray/inversion/regularizations/test_matern_adaptive_brightness_kernel.py b/test_autoarray/inversion/regularizations/test_matern_adaptive_brightness_kernel.py new file mode 100644 index 000000000..4a8554d83 --- /dev/null +++ b/test_autoarray/inversion/regularizations/test_matern_adaptive_brightness_kernel.py @@ -0,0 +1,55 @@ +import pytest + +import autoarray as aa +import numpy as np + +np.set_printoptions(threshold=np.inf) + + +def test__regularization_matrix(): + + reg = aa.reg.MaternAdaptiveBrightnessKernel(coefficient=1.0, scale=2.0, nu=2.0, rho=1.0) + + neighbors = np.array( + [ + [1, 4, -1, -1], + [2, 4, 0, -1], + [3, 4, 5, 1], + [5, 2, -1, -1], + [5, 0, 1, 2], + [2, 3, 4, -1], + ] + ) + + neighbors_sizes = np.array([2, 3, 4, 2, 4, 3]) + pixel_signals = np.array([1.0, 2.0, 3.0, 4.0, 5.0, 6.0]) + + mesh_grid = aa.m.MockMeshGrid(neighbors=neighbors, neighbors_sizes=neighbors_sizes) + + source_plane_mesh_grid = aa.Grid2D.no_mask( + values=[[0.1, 0.1], [1.1, 0.6], [2.1, 0.1], [0.4, 1.1], [1.1, 7.1], [2.1, 1.1]], + shape_native=(3, 2), + pixel_scales=1.0, + ) + + mapper = aa.m.MockMapper( + source_plane_mesh_grid=source_plane_mesh_grid, pixel_signals=pixel_signals + ) + + regularization_matrix = reg.regularization_matrix_from(linear_obj=mapper) + + assert regularization_matrix[0, 0] == pytest.approx(18.7439565009, 1.0e-4) + assert regularization_matrix[0, 1] == pytest.approx(-8.786547368, 1.0e-4) + + reg = aa.reg.MaternAdaptiveBrightnessKernel(coefficient=1.5, scale=2.5, nu=2.5, rho=1.5) + + pixel_signals = np.array([1.0, 2.0, 3.0, 4.0, 5.0, 6.0]) + + mapper = aa.m.MockMapper( + source_plane_mesh_grid=source_plane_mesh_grid, pixel_signals=pixel_signals + ) + + regularization_matrix = reg.regularization_matrix_from(linear_obj=mapper) + + assert regularization_matrix[0, 0] == pytest.approx(121.0190770, 1.0e-4) + assert regularization_matrix[0, 1] == pytest.approx(-66.9580331, 1.0e-4) \ No newline at end of file From e366b497b5b354ed1d4cc417d908c1247841ec12 Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Fri, 13 Feb 2026 11:29:40 +0000 Subject: [PATCH 03/11] end to end matern tests and skip areas added --- .../matern_adaptive_brightness_kernel.py | 15 ++- .../inversion/regularization/matern_kernel.py | 20 +++- autoarray/preloads.py | 2 + autoarray/structures/mesh/delaunay_2d.py | 112 ++++++++++++++++-- .../test_matern_adaptive_brightness_kernel.py | 10 +- 5 files changed, 135 insertions(+), 24 deletions(-) diff --git a/autoarray/inversion/regularization/matern_adaptive_brightness_kernel.py b/autoarray/inversion/regularization/matern_adaptive_brightness_kernel.py index b629a10a0..952a34e3b 100644 --- a/autoarray/inversion/regularization/matern_adaptive_brightness_kernel.py +++ b/autoarray/inversion/regularization/matern_adaptive_brightness_kernel.py @@ -9,6 +9,7 @@ from autoarray.inversion.regularization.matern_kernel import matern_kernel + def matern_cov_matrix_from( scale: float, nu: float, @@ -47,7 +48,7 @@ def matern_cov_matrix_from( # Pairwise distances (broadcasted) # -------------------------------- diff = pixel_points[:, None, :] - pixel_points[None, :, :] # (N, N, 2) - d_ij = xp.sqrt(diff[..., 0] ** 2 + diff[..., 1] ** 2) # (N, N) + d_ij = xp.sqrt(diff[..., 0] ** 2 + diff[..., 1] ** 2) # (N, N) # -------------------------------- # Base Matérn covariance @@ -113,7 +114,9 @@ def __init__( super().__init__(coefficient=coefficient, scale=scale, nu=nu) self.rho = rho - def covariance_kernel_weights_from(self, linear_obj: LinearObj, xp=np) -> np.ndarray: + def covariance_kernel_weights_from( + self, linear_obj: LinearObj, xp=np + ) -> np.ndarray: """ Returns per-pixel kernel weights that adapt to the reconstructed pixel brightness. """ @@ -126,10 +129,12 @@ def covariance_kernel_weights_from(self, linear_obj: LinearObj, xp=np) -> np.nda return xp.exp(-self.rho * (1.0 - pixel_signals / max_signal)) def regularization_matrix_from(self, linear_obj: LinearObj, xp=np) -> np.ndarray: - kernel_weights = self.covariance_kernel_weights_from(linear_obj=linear_obj, xp=xp) + kernel_weights = self.covariance_kernel_weights_from( + linear_obj=linear_obj, xp=xp + ) # Follow the xp pattern used in the Matérn kernel module (often `.array` for grids). - pixel_points = linear_obj.source_plane_mesh_grid + pixel_points = linear_obj.source_plane_mesh_grid.array covariance_matrix = matern_cov_matrix_from( scale=self.scale, @@ -145,4 +150,4 @@ def regularization_weights_from(self, linear_obj: LinearObj, xp=np) -> np.ndarra """ Returns the regularization weights of this regularization scheme. """ - return 1.0 / self.covariance_kernel_weights_from(linear_obj=linear_obj, xp=xp) \ No newline at end of file + return 1.0 / self.covariance_kernel_weights_from(linear_obj=linear_obj, xp=xp) diff --git a/autoarray/inversion/regularization/matern_kernel.py b/autoarray/inversion/regularization/matern_kernel.py index c80c5d233..3e38757ec 100644 --- a/autoarray/inversion/regularization/matern_kernel.py +++ b/autoarray/inversion/regularization/matern_kernel.py @@ -44,6 +44,20 @@ def kv_xp(v, z, xp=np): ) +def gamma_xp(x, xp=np): + """ + XP-compatible Gamma(x). + """ + if xp is np: + import scipy.special as sc + + return sc.gamma(x) + else: + import jax.scipy.special as jsp + + return jsp.gamma(x) + + def matern_kernel(r, l: float = 1.0, v: float = 0.5, xp=np): """ XP-compatible Matérn kernel. @@ -55,7 +69,7 @@ def matern_kernel(r, l: float = 1.0, v: float = 0.5, xp=np): z = xp.sqrt(2.0 * v) * r / l - part1 = 2.0 ** (1.0 - v) / math.gamma(v) # scalar constant + part1 = 2.0 ** (1.0 - v) / gamma_xp(v, xp) # scalar constant part2 = z**v part3 = kv_xp(v, z, xp) @@ -141,8 +155,8 @@ def __init__(self, coefficient: float = 1.0, scale: float = 1.0, nu: float = 0.5 """ self.coefficient = coefficient - self.scale = float(scale) - self.nu = float(nu) + self.scale = scale + self.nu = nu super().__init__() def regularization_weights_from(self, linear_obj: LinearObj, xp=np) -> np.ndarray: diff --git a/autoarray/preloads.py b/autoarray/preloads.py index ba9e4e4dd..f9c858e9b 100644 --- a/autoarray/preloads.py +++ b/autoarray/preloads.py @@ -26,6 +26,7 @@ def __init__( linear_light_profile_blurred_mapping_matrix=None, use_voronoi_areas: bool = True, areas_factor: float = 0.5, + skip_areas: bool = False, ): """ Stores preloaded arrays and matrices used during pixelized linear inversions, improving both performance @@ -123,3 +124,4 @@ def __init__( self.use_voronoi_areas = use_voronoi_areas self.areas_factor = areas_factor + self.skip_areas = skip_areas diff --git a/autoarray/structures/mesh/delaunay_2d.py b/autoarray/structures/mesh/delaunay_2d.py index 74fb18112..dc550d3b4 100644 --- a/autoarray/structures/mesh/delaunay_2d.py +++ b/autoarray/structures/mesh/delaunay_2d.py @@ -339,6 +339,65 @@ def pix_indexes_for_sub_slim_index_delaunay_from( return out +def scipy_delaunay_matern(points_np, query_points_np): + """ + Minimal SciPy Delaunay callback for Matérn regularization. + + Returns only what’s needed for mapping: + - points (tri.points) + - simplices_padded + - mappings (pix indexes for each query point) + """ + + max_simplices = 2 * points_np.shape[0] + + # --- Delaunay mesh --- + tri = Delaunay(points_np) + + points = tri.points.astype(points_np.dtype) + simplices = tri.simplices.astype(np.int32) + + # --- Pad simplices to fixed shape for JAX --- + simplices_padded = -np.ones((max_simplices, 3), dtype=np.int32) + simplices_padded[: simplices.shape[0]] = simplices + + # --- find_simplex for query points --- + simplex_idx = tri.find_simplex(query_points_np).astype(np.int32) # (Q,) + + mappings = pix_indexes_for_sub_slim_index_delaunay_from( + source_plane_data_grid=query_points_np, + simplex_index_for_sub_slim_index=simplex_idx, + pix_indexes_for_simplex_index=simplices, + delaunay_points=points_np, + ) + + return points, simplices_padded, mappings + + +def jax_delaunay_matern(points, query_points): + """ + JAX wrapper using pure_callback to run SciPy Delaunay on CPU, + returning only the minimal outputs needed for Matérn usage. + """ + import jax + import jax.numpy as jnp + + N = points.shape[0] + Q = query_points.shape[0] + max_simplices = 2 * N + + points_shape = jax.ShapeDtypeStruct((N, 2), points.dtype) + simplices_padded_shape = jax.ShapeDtypeStruct((max_simplices, 3), jnp.int32) + mappings_shape = jax.ShapeDtypeStruct((Q, 3), jnp.int32) + + return jax.pure_callback( + lambda pts, qpts: scipy_delaunay_matern(np.asarray(pts), np.asarray(qpts)), + (points_shape, simplices_padded_shape, mappings_shape), + points, + query_points, + ) + + class DelaunayInterface: def __init__( @@ -466,33 +525,60 @@ def delaunay(self) -> "scipy.spatial.Delaunay": use_voronoi_areas = self.preloads.use_voronoi_areas areas_factor = self.preloads.areas_factor + skip_areas = self.preloads.skip_areas else: use_voronoi_areas = True areas_factor = 0.5 + skip_areas = False - if self._xp.__name__.startswith("jax"): + if not skip_areas: - import jax.numpy as jnp + if self._xp.__name__.startswith("jax"): - points, simplices, mappings, split_points, splitted_mappings = jax_delaunay( - points=self.mesh_grid_xy, - query_points=self._source_plane_data_grid_over_sampled, - use_voronoi_areas=use_voronoi_areas, - areas_factor=areas_factor, - ) + import jax.numpy as jnp + + points, simplices, mappings, split_points, splitted_mappings = ( + jax_delaunay( + points=self.mesh_grid_xy, + query_points=self._source_plane_data_grid_over_sampled, + use_voronoi_areas=use_voronoi_areas, + areas_factor=areas_factor, + ) + ) + + else: + + points, simplices, mappings, split_points, splitted_mappings = ( + scipy_delaunay( + points_np=self.mesh_grid_xy, + query_points_np=self._source_plane_data_grid_over_sampled, + use_voronoi_areas=use_voronoi_areas, + areas_factor=areas_factor, + ) + ) else: - points, simplices, mappings, split_points, splitted_mappings = ( - scipy_delaunay( + if self._xp.__name__.startswith("jax"): + + import jax.numpy as jnp + + points, simplices, mappings = jax_delaunay_matern( + points=self.mesh_grid_xy, + query_points=self._source_plane_data_grid_over_sampled, + ) + + else: + + points, simplices, mappings = scipy_delaunay_matern( points_np=self.mesh_grid_xy, query_points_np=self._source_plane_data_grid_over_sampled, - use_voronoi_areas=use_voronoi_areas, - areas_factor=areas_factor, ) - ) + + split_points = None + splitted_mappings = None return DelaunayInterface( points=points, diff --git a/test_autoarray/inversion/regularizations/test_matern_adaptive_brightness_kernel.py b/test_autoarray/inversion/regularizations/test_matern_adaptive_brightness_kernel.py index 4a8554d83..c0316376c 100644 --- a/test_autoarray/inversion/regularizations/test_matern_adaptive_brightness_kernel.py +++ b/test_autoarray/inversion/regularizations/test_matern_adaptive_brightness_kernel.py @@ -8,7 +8,9 @@ def test__regularization_matrix(): - reg = aa.reg.MaternAdaptiveBrightnessKernel(coefficient=1.0, scale=2.0, nu=2.0, rho=1.0) + reg = aa.reg.MaternAdaptiveBrightnessKernel( + coefficient=1.0, scale=2.0, nu=2.0, rho=1.0 + ) neighbors = np.array( [ @@ -41,7 +43,9 @@ def test__regularization_matrix(): assert regularization_matrix[0, 0] == pytest.approx(18.7439565009, 1.0e-4) assert regularization_matrix[0, 1] == pytest.approx(-8.786547368, 1.0e-4) - reg = aa.reg.MaternAdaptiveBrightnessKernel(coefficient=1.5, scale=2.5, nu=2.5, rho=1.5) + reg = aa.reg.MaternAdaptiveBrightnessKernel( + coefficient=1.5, scale=2.5, nu=2.5, rho=1.5 + ) pixel_signals = np.array([1.0, 2.0, 3.0, 4.0, 5.0, 6.0]) @@ -52,4 +56,4 @@ def test__regularization_matrix(): regularization_matrix = reg.regularization_matrix_from(linear_obj=mapper) assert regularization_matrix[0, 0] == pytest.approx(121.0190770, 1.0e-4) - assert regularization_matrix[0, 1] == pytest.approx(-66.9580331, 1.0e-4) \ No newline at end of file + assert regularization_matrix[0, 1] == pytest.approx(-66.9580331, 1.0e-4) From f72440413869211c1fcfd0393dacaa8d63d43b03 Mon Sep 17 00:00:00 2001 From: James Nightingale Date: Fri, 13 Feb 2026 12:56:00 +0000 Subject: [PATCH 04/11] Update autoarray/inversion/regularization/matern_adaptive_brightness_kernel.py Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- .../regularization/matern_adaptive_brightness_kernel.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/autoarray/inversion/regularization/matern_adaptive_brightness_kernel.py b/autoarray/inversion/regularization/matern_adaptive_brightness_kernel.py index 952a34e3b..1eb8dd3a7 100644 --- a/autoarray/inversion/regularization/matern_adaptive_brightness_kernel.py +++ b/autoarray/inversion/regularization/matern_adaptive_brightness_kernel.py @@ -110,6 +110,10 @@ def __init__( The typical scale of the exponential regularization pattern. nu Controls the derivative of the regularization pattern (`nu=0.5` is a Gaussian). + rho + Controls how strongly the kernel weights adapt to pixel brightness. Larger values make bright pixels + receive significantly higher weights (and faint pixels lower weights), while smaller values produce a + more uniform weighting. Typical values are of order unity (e.g. 0.5–2.0). """ super().__init__(coefficient=coefficient, scale=scale, nu=nu) self.rho = rho From 9184f2860954b9b677b03d478f6083e5d3b402fc Mon Sep 17 00:00:00 2001 From: James Nightingale Date: Fri, 13 Feb 2026 12:56:17 +0000 Subject: [PATCH 05/11] Update autoarray/inversion/regularization/matern_adaptive_brightness_kernel.py Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- .../regularization/matern_adaptive_brightness_kernel.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/autoarray/inversion/regularization/matern_adaptive_brightness_kernel.py b/autoarray/inversion/regularization/matern_adaptive_brightness_kernel.py index 1eb8dd3a7..d14c11541 100644 --- a/autoarray/inversion/regularization/matern_adaptive_brightness_kernel.py +++ b/autoarray/inversion/regularization/matern_adaptive_brightness_kernel.py @@ -107,9 +107,11 @@ def __init__( coefficient The regularization coefficient which controls the degree of smooth of the inversion reconstruction. scale - The typical scale of the exponential regularization pattern. + The typical scale (correlation length) of the Matérn regularization kernel. nu - Controls the derivative of the regularization pattern (`nu=0.5` is a Gaussian). + Controls the smoothness (differentiability) of the Matérn kernel; ``nu=0.5`` corresponds to an + exponential (Ornstein–Uhlenbeck) kernel, while a Gaussian covariance is obtained in the limit + as ``nu`` approaches infinity. rho Controls how strongly the kernel weights adapt to pixel brightness. Larger values make bright pixels receive significantly higher weights (and faint pixels lower weights), while smaller values produce a From a14f39a7a30dd8f14584213e4bd24381c591455a Mon Sep 17 00:00:00 2001 From: James Nightingale Date: Fri, 13 Feb 2026 12:56:27 +0000 Subject: [PATCH 06/11] Update autoarray/inversion/regularization/matern_adaptive_brightness_kernel.py Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- .../regularization/matern_adaptive_brightness_kernel.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/autoarray/inversion/regularization/matern_adaptive_brightness_kernel.py b/autoarray/inversion/regularization/matern_adaptive_brightness_kernel.py index d14c11541..95a6d2f16 100644 --- a/autoarray/inversion/regularization/matern_adaptive_brightness_kernel.py +++ b/autoarray/inversion/regularization/matern_adaptive_brightness_kernel.py @@ -87,7 +87,7 @@ def __init__( For this regularization scheme, every pixel is regularized with every other pixel. This contrasts many other schemes, where regularization is based on neighboring (e.g. do the pixels share a Delaunay edge?) or computing - derivates around the center of the pixel (where nearby pixels are regularization locally in similar ways). + derivatives around the center of the pixel (where nearby pixels are regularization locally in similar ways). This makes the regularization matrix fully dense and therefore maybe change the run times of the solution. It also leads to more overall smoothing which can lead to more stable linear inversions. From 244caf092904832ab057d34f2bccd969618f4538 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Fri, 13 Feb 2026 12:56:49 +0000 Subject: [PATCH 07/11] Initial plan From 12de13be4e0e070719ccbe3402c326812604faad Mon Sep 17 00:00:00 2001 From: James Nightingale Date: Fri, 13 Feb 2026 12:56:52 +0000 Subject: [PATCH 08/11] Update test_autoarray/inversion/regularizations/test_matern_adaptive_brightness_kernel.py Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- .../test_matern_adaptive_brightness_kernel.py | 14 -------------- 1 file changed, 14 deletions(-) diff --git a/test_autoarray/inversion/regularizations/test_matern_adaptive_brightness_kernel.py b/test_autoarray/inversion/regularizations/test_matern_adaptive_brightness_kernel.py index c0316376c..de514975f 100644 --- a/test_autoarray/inversion/regularizations/test_matern_adaptive_brightness_kernel.py +++ b/test_autoarray/inversion/regularizations/test_matern_adaptive_brightness_kernel.py @@ -12,22 +12,8 @@ def test__regularization_matrix(): coefficient=1.0, scale=2.0, nu=2.0, rho=1.0 ) - neighbors = np.array( - [ - [1, 4, -1, -1], - [2, 4, 0, -1], - [3, 4, 5, 1], - [5, 2, -1, -1], - [5, 0, 1, 2], - [2, 3, 4, -1], - ] - ) - - neighbors_sizes = np.array([2, 3, 4, 2, 4, 3]) pixel_signals = np.array([1.0, 2.0, 3.0, 4.0, 5.0, 6.0]) - mesh_grid = aa.m.MockMeshGrid(neighbors=neighbors, neighbors_sizes=neighbors_sizes) - source_plane_mesh_grid = aa.Grid2D.no_mask( values=[[0.1, 0.1], [1.1, 0.6], [2.1, 0.1], [0.4, 1.1], [1.1, 7.1], [2.1, 1.1]], shape_native=(3, 2), From 29a2d893275e6c25eded1b85c2985304c756df23 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Fri, 13 Feb 2026 13:00:00 +0000 Subject: [PATCH 09/11] Add documentation for skip_areas, use_voronoi_areas, and areas_factor parameters Co-authored-by: Jammy2211 <23455639+Jammy2211@users.noreply.github.com> --- autoarray/preloads.py | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/autoarray/preloads.py b/autoarray/preloads.py index f9c858e9b..1b4e00d28 100644 --- a/autoarray/preloads.py +++ b/autoarray/preloads.py @@ -82,6 +82,16 @@ def __init__( inversion, with the other component being the pixelization's pixels. These are fixed when the lens light is fixed to the maximum likelihood solution, allowing the blurred mapping matrix to be preloaded, but the intensity values will still be solved for during the inversion. + use_voronoi_areas + Whether to use Voronoi areas during Delaunay triangulation. When True, computes areas for each Voronoi + region which can be used in certain regularization schemes. Default is True. + areas_factor + Factor used to scale the Voronoi areas during split point computation. Default is 0.5. + skip_areas + Whether to skip Voronoi area calculations and split point computations during Delaunay triangulation. + When True, the Delaunay interface returns only the minimal set of outputs (points, simplices, mappings) + without computing split_points or splitted_mappings. This optimization is useful for regularization + schemes like Matérn kernels that don't require area-based calculations. Default is False. """ self.mapper_indices = None self.source_pixel_zeroed_indices = None From 917b78cef51722103b0f33becff32b00e0036d2a Mon Sep 17 00:00:00 2001 From: James Nightingale Date: Fri, 13 Feb 2026 13:30:43 +0000 Subject: [PATCH 10/11] Update autoarray/inversion/regularization/matern_adaptive_brightness_kernel.py Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- .../regularization/matern_adaptive_brightness_kernel.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/autoarray/inversion/regularization/matern_adaptive_brightness_kernel.py b/autoarray/inversion/regularization/matern_adaptive_brightness_kernel.py index 95a6d2f16..6442d0343 100644 --- a/autoarray/inversion/regularization/matern_adaptive_brightness_kernel.py +++ b/autoarray/inversion/regularization/matern_adaptive_brightness_kernel.py @@ -89,7 +89,7 @@ def __init__( schemes, where regularization is based on neighboring (e.g. do the pixels share a Delaunay edge?) or computing derivatives around the center of the pixel (where nearby pixels are regularization locally in similar ways). - This makes the regularization matrix fully dense and therefore maybe change the run times of the solution. + This makes the regularization matrix fully dense and therefore may change the run times of the solution. It also leads to more overall smoothing which can lead to more stable linear inversions. For the weighted regularization scheme, each pixel is given an 'effective regularization weight', which is From 95cf5b74e5aa12c0a4664800bf6b598dc1f3b4a8 Mon Sep 17 00:00:00 2001 From: James Nightingale Date: Fri, 13 Feb 2026 13:30:57 +0000 Subject: [PATCH 11/11] Update autoarray/structures/mesh/delaunay_2d.py Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- autoarray/structures/mesh/delaunay_2d.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/autoarray/structures/mesh/delaunay_2d.py b/autoarray/structures/mesh/delaunay_2d.py index dc550d3b4..961944b75 100644 --- a/autoarray/structures/mesh/delaunay_2d.py +++ b/autoarray/structures/mesh/delaunay_2d.py @@ -346,7 +346,9 @@ def scipy_delaunay_matern(points_np, query_points_np): Returns only what’s needed for mapping: - points (tri.points) - simplices_padded - - mappings (pix indexes for each query point) + - mappings: integer array of pixel indices for each query point, + typically of shape (Q, 3), where each row gives the indices of the + Delaunay mesh vertices ("pixels") associated with that query point. """ max_simplices = 2 * points_np.shape[0]