|
| 1 | +from __future__ import annotations |
| 2 | +import numpy as np |
| 3 | +from typing import TYPE_CHECKING |
| 4 | + |
| 5 | +from autoarray.inversion.regularization.matern_kernel import MaternKernel |
| 6 | + |
| 7 | +if TYPE_CHECKING: |
| 8 | + from autoarray.inversion.linear_obj.linear_obj import LinearObj |
| 9 | + |
| 10 | +from autoarray.inversion.regularization.matern_kernel import matern_kernel |
| 11 | + |
| 12 | + |
| 13 | +def matern_cov_matrix_from( |
| 14 | + scale: float, |
| 15 | + nu: float, |
| 16 | + pixel_points, |
| 17 | + weights=None, |
| 18 | + xp=np, |
| 19 | +): |
| 20 | + """ |
| 21 | + Construct the regularization covariance matrix (N x N) using a Matérn kernel, |
| 22 | + optionally modulated by per-pixel weights. |
| 23 | +
|
| 24 | + If `weights` is provided (shape [N]), the covariance is: |
| 25 | + C_ij = K(d_ij; scale, nu) * w_i * w_j |
| 26 | + with a small diagonal jitter added for numerical stability. |
| 27 | +
|
| 28 | + Parameters |
| 29 | + ---------- |
| 30 | + scale |
| 31 | + Typical correlation length of the Matérn kernel. |
| 32 | + nu |
| 33 | + Smoothness parameter of the Matérn kernel. |
| 34 | + pixel_points |
| 35 | + Array-like of shape [N, 2] with (y, x) coordinates (or any 2D coords; only distances matter). |
| 36 | + weights |
| 37 | + Optional array-like of shape [N]. If None, treated as all ones. |
| 38 | + xp |
| 39 | + Backend (numpy or jax.numpy). |
| 40 | +
|
| 41 | + Returns |
| 42 | + ------- |
| 43 | + covariance_matrix |
| 44 | + Array of shape [N, N]. |
| 45 | + """ |
| 46 | + |
| 47 | + # -------------------------------- |
| 48 | + # Pairwise distances (broadcasted) |
| 49 | + # -------------------------------- |
| 50 | + diff = pixel_points[:, None, :] - pixel_points[None, :, :] # (N, N, 2) |
| 51 | + d_ij = xp.sqrt(diff[..., 0] ** 2 + diff[..., 1] ** 2) # (N, N) |
| 52 | + |
| 53 | + # -------------------------------- |
| 54 | + # Base Matérn covariance |
| 55 | + # -------------------------------- |
| 56 | + covariance_matrix = matern_kernel(d_ij, l=scale, v=nu, xp=xp) # (N, N) |
| 57 | + |
| 58 | + # -------------------------------- |
| 59 | + # Apply weights: C_ij *= w_i * w_j |
| 60 | + # (broadcasted outer product, JAX-safe) |
| 61 | + # -------------------------------- |
| 62 | + if weights is not None: |
| 63 | + w = xp.asarray(weights) |
| 64 | + # Ensure shape (N,) -> outer product (N,1)*(1,N) -> (N,N) |
| 65 | + covariance_matrix = covariance_matrix * (w[:, None] * w[None, :]) |
| 66 | + |
| 67 | + # -------------------------------- |
| 68 | + # Add diagonal jitter (JAX-safe) |
| 69 | + # -------------------------------- |
| 70 | + pixels = pixel_points.shape[0] |
| 71 | + covariance_matrix = covariance_matrix + 1e-8 * xp.eye(pixels) |
| 72 | + |
| 73 | + return covariance_matrix |
| 74 | + |
| 75 | + |
| 76 | +class MaternAdaptiveBrightnessKernel(MaternKernel): |
| 77 | + def __init__( |
| 78 | + self, |
| 79 | + coefficient: float = 1.0, |
| 80 | + scale: float = 1.0, |
| 81 | + nu: float = 0.5, |
| 82 | + rho: float = 1.0, |
| 83 | + ): |
| 84 | + """ |
| 85 | + Regularization which uses a Matern smoothing kernel to regularize the solution with regularization weights |
| 86 | + that adapt to the brightness of the source being reconstructed. |
| 87 | +
|
| 88 | + For this regularization scheme, every pixel is regularized with every other pixel. This contrasts many other |
| 89 | + schemes, where regularization is based on neighboring (e.g. do the pixels share a Delaunay edge?) or computing |
| 90 | + derivatives around the center of the pixel (where nearby pixels are regularization locally in similar ways). |
| 91 | +
|
| 92 | + This makes the regularization matrix fully dense and therefore may change the run times of the solution. |
| 93 | + It also leads to more overall smoothing which can lead to more stable linear inversions. |
| 94 | +
|
| 95 | + For the weighted regularization scheme, each pixel is given an 'effective regularization weight', which is |
| 96 | + applied when each set of pixel neighbors are regularized with one another. The motivation of this is that |
| 97 | + different regions of a pixelization's mesh require different levels of regularization (e.g., high smoothing where the |
| 98 | + no signal is present and less smoothing where it is, see (Nightingale, Dye and Massey 2018)). |
| 99 | +
|
| 100 | + This scheme is not used by Vernardos et al. (2022): https://arxiv.org/abs/2202.09378, but it follows |
| 101 | + a similar approach. |
| 102 | +
|
| 103 | + A full description of regularization and this matrix can be found in the parent `AbstractRegularization` class. |
| 104 | +
|
| 105 | + Parameters |
| 106 | + ---------- |
| 107 | + coefficient |
| 108 | + The regularization coefficient which controls the degree of smooth of the inversion reconstruction. |
| 109 | + scale |
| 110 | + The typical scale (correlation length) of the Matérn regularization kernel. |
| 111 | + nu |
| 112 | + Controls the smoothness (differentiability) of the Matérn kernel; ``nu=0.5`` corresponds to an |
| 113 | + exponential (Ornstein–Uhlenbeck) kernel, while a Gaussian covariance is obtained in the limit |
| 114 | + as ``nu`` approaches infinity. |
| 115 | + rho |
| 116 | + Controls how strongly the kernel weights adapt to pixel brightness. Larger values make bright pixels |
| 117 | + receive significantly higher weights (and faint pixels lower weights), while smaller values produce a |
| 118 | + more uniform weighting. Typical values are of order unity (e.g. 0.5–2.0). |
| 119 | + """ |
| 120 | + super().__init__(coefficient=coefficient, scale=scale, nu=nu) |
| 121 | + self.rho = rho |
| 122 | + |
| 123 | + def covariance_kernel_weights_from( |
| 124 | + self, linear_obj: LinearObj, xp=np |
| 125 | + ) -> np.ndarray: |
| 126 | + """ |
| 127 | + Returns per-pixel kernel weights that adapt to the reconstructed pixel brightness. |
| 128 | + """ |
| 129 | + # Assumes linear_obj.pixel_signals_from is xp-aware elsewhere in the codebase. |
| 130 | + pixel_signals = linear_obj.pixel_signals_from(signal_scale=1.0, xp=xp) |
| 131 | + |
| 132 | + max_signal = xp.max(pixel_signals) |
| 133 | + max_signal = xp.maximum(max_signal, 1e-8) # avoid divide-by-zero (JAX-safe) |
| 134 | + |
| 135 | + return xp.exp(-self.rho * (1.0 - pixel_signals / max_signal)) |
| 136 | + |
| 137 | + def regularization_matrix_from(self, linear_obj: LinearObj, xp=np) -> np.ndarray: |
| 138 | + kernel_weights = self.covariance_kernel_weights_from( |
| 139 | + linear_obj=linear_obj, xp=xp |
| 140 | + ) |
| 141 | + |
| 142 | + # Follow the xp pattern used in the Matérn kernel module (often `.array` for grids). |
| 143 | + pixel_points = linear_obj.source_plane_mesh_grid.array |
| 144 | + |
| 145 | + covariance_matrix = matern_cov_matrix_from( |
| 146 | + scale=self.scale, |
| 147 | + pixel_points=pixel_points, |
| 148 | + nu=self.nu, |
| 149 | + weights=kernel_weights, |
| 150 | + xp=xp, |
| 151 | + ) |
| 152 | + |
| 153 | + return self.coefficient * xp.linalg.inv(covariance_matrix) |
| 154 | + |
| 155 | + def regularization_weights_from(self, linear_obj: LinearObj, xp=np) -> np.ndarray: |
| 156 | + """ |
| 157 | + Returns the regularization weights of this regularization scheme. |
| 158 | + """ |
| 159 | + return 1.0 / self.covariance_kernel_weights_from(linear_obj=linear_obj, xp=xp) |
0 commit comments