lenstronomy.ImSim.SourceReconstruction package

Submodules

lenstronomy.ImSim.SourceReconstruction.pixelated_source_reconstruction module

class PixelatedSourceReconstruction(data_class, psf_class, lens_model_class, source_pixel_grid_class)

Bases: object

This class defines useful functions for pixelated source plane reconstruction in gravitational lensing. It handles initialization of data, lens model, and PSF kernel, and provides methods for generating the M and b matrices (definition see ‘placeholder for Nan Zhang’s paper) for source reconstruction based on different likelihood methods.

All sparse matrices (sp) in this class are represented as a list of lists: [[y_coord, x_coord, pixel_value], …], where [y_coord, x_coord] are the pixel coordinates and pixel_value is the corresponding pixel’s value.

__init__(data_class, psf_class, lens_model_class, source_pixel_grid_class)

Initializes the PixelatedSourceReconstruction class. This sets up the necessary data, PSF, lens model, and source grid for subsequent source reconstruction matrix generation.

Parameters:

• data_class – ImageData() class instance (for the observed image data)

• psf_class – PSF() class instance (for the observed image data)

• lens_model_class – LensModel class instance

• source_pixel_grid_class – PixelGrid() class instance (defining the source plane grid)

Raises:

ValueError –

• If the source pixel grid has rotational components or non-uniform pixel widths.

• If the PSF kernel size is improperly sized for interferometric likelihood methods.

generate_M_b(kwargs_lens, verbose=False, show_progress=True)

Generates the M matrix and the b vector for source reconstruction based on the selected likelihood method. Definitions of the M and b is given by (placeholder for Nan Zhang’s paper).

Parameters:

• kwargs_lens – List of keyword arguments for the lens_model_class.

• verbose – If True, print progress messages during matrix generation steps. Defaults to False.

• show_progress – If True, show progress bar of generating the M matrix. Defaults to True.

Returns:

(M, b) tuple, where M is the matrix and b is the vector.

generate_M_b_diagonal_likelihood(kwargs_lens, verbose=False, show_progress=True)

Generates M and b matrices assuming spatially uncorrelated noise with uniform RMS across the image. This approach is typically used for CCD image data.

This method performs lensing, convolution, and then computes M and b.

Parameters:

• kwargs_lens – List of keyword arguments for the lens_model_class.

• verbose – If True, print progress messages during matrix generation steps. Defaults to False.

• show_progress – If True, show progress bar of generating the M matrix. Defaults to True.

Returns:

(M, b) tuple, where M is the matrix and b is the vector.

generate_M_b_interferometry_natwt_likelihood(kwargs_lens, verbose=False, show_progress=True)

Generates the M and b matrices for interferometric data with natural weighting.

This method integrates the convolution step implicitly via specialized sparse product functions.

Parameters:

• kwargs_lens – List of keyword arguments for the lens_model_class.

• verbose – If True, print progress messages during matrix generation steps. Defaults to False.

• show_progress – If True, show progress bar of generating the M matrix. Defaults to True.

Returns:

(M, b) tuple, where M is the matrix and b is the vector.

lens_pixel_source_of_a_rectangular_region(kwargs_lens)

Maps image plane pixels to source plane pixels within a specified rectangular source grid, considering lensing deflections and applying bilinear interpolation.

Parameters:

kwargs_lens – List of keyword arguments for the lens_model_class.

Returns:

A list of lists. Each element in the outer list corresponds to a single source pixel within the defined source grid (ordered by their linear index). Each inner list contains [y_coord, x_coord, weight] tuples, indicating that the source pixel at this index contributes with weight to the image plane pixel at (y_coord, x_coord) when lensed.

Return type:

list

lens_an_image_by_rayshooting(kwargs_lens, source_image)

Lenses a pixelated source plane image to the image plane using ray-shooting and bilinear interpolation. The imput image should have the same dimension and coordinates defined by source_pixel_grid_class.

This method works by iterating through each pixel in the image plane, ray-shooting back to the source plane to find the corresponding source coordinate, and then interpolating the flux from the input source image at that coordinate. This provides an approximate lensed image, as flux between exact source pixels is derived via interpolation. Note that the primary beam will NOT be applied on the lensed image in this function.

Parameters:

• kwargs_lens – List of keyword arguments for the lens_model_class.

• image (numpy.ndarray) – 2D NumPy array representing the pixelated source plane image. Expected to have dimensions defined by source_pixel_grid_class.

Returns:

2D NumPy array representing the lensed image in the image plane.

Return type:

numpy.ndarray

Raises:

ValueError – If the input image dimensions do not match the dimensions of the defined source pixel grid (self._ny_source, self._nx_source).

sparse_to_array(sparse)

Converts a sparse image representation (list of [y_coord, x_coord, value] tuples) to a 2D NumPy array.

Parameters:

sparse – A list representing non-zero elements of the sparse image.

Returns:

A 2D NumPy array representing the full image.

Return type:

numpy.ndarray

sum_sparse_elementwise_product(sparse, ordinary)

Computes the element-wise sum of products between a sparse matrix and a dense 2D NumPy array image.

Parameters:

• sparse – Sparse matrix representation (list of [y_coord, x_coord, value] tuples).

• ordinary – A 2D NumPy array (dense matrix).

Returns:

The sum of the element-wise products.

Return type:

float

sparse_convolve_and_dot_product(sp1, sp2, kernel=None)

Computes the convolution product of two sparse matrices using a given kernel. Equivalent to (sp1 * kernel) . (sp2) where ‘*’ is convolution and ‘.’ is dot product.

Parameters:

• sp1 – First sparse matrix representation.

• sp2 – Second sparse matrix representation.

• kernel – The 2D PSF kernel (NumPy array). Assumed to be square with odd dimensions, with its center at the central pixel. If None, self._kernel is used

Returns:

The result of the convolution product.

Return type:

float

sparse_convolution(sp, kernel=None)

Performs convolution of a sparse matrix with a given kernel.

Parameters:

• sp – Sparse matrix representation.

• kernel – The 2D PSF kernel (NumPy array). Assumed to be square with odd dimensions, with its center at the central pixel. If None, self._kernel is used

Returns:

A 2D NumPy array representing the convolved image.

Return type:

numpy.ndarray

lenstronomy.ImSim.SourceReconstruction.regularization_matrix_pixel module

This module provides functions to construct regularization matrices for pixelated source plane reconstructions in gravitational lensing. These matrices are used to constrain the pixel amplitudes in linear inversion problems, promoting desired properties such as smoothness.

Currently supported regularization types include: - Zeroth-order (L2 norm on pixel amplitudes) - Gradient (L2 norm on spatial gradients, promoting smoothness) - Curvature (L2 norm on second-order derivatives, promoting smoother variations)

pixelated_regularization_matrix(xlen, ylen, regularization_type)

Constructs the regularization matrix for a rectangular pixelated source region.

The regularization term for pixel amplitudes \(\mathbf{a}\) (flattened into a 1D vector) is generally expressed as \(\mathbf{a}^T U \mathbf{a}\), where \(U\) is the regularization matrix. This function returns the matrix \(U\).

The pixels are aligned in a 1D vector a by flattening the 2D source grid in a row-major order. Specifically, if a pixel has 2D indices \((i_x, i_y)\) (where \(i_x\) is the column index from 0 to xlen-1 and \(i_y\) is the row index from 0 to ylen-1), its 1D index \(i\) is given by \(i = i_y \cdot xlen + i_x\).

Parameters:

• xlen – int, Number of pixels in the x-direction (horizontal dimension of the source grid).

• ylen – int, Number of pixels in the y-direction (vertical dimension of the source grid).

• regularization_type – str, Type of regularization to apply. Supported options are ‘zeroth_order’, ‘gradient’, ‘curvature’.

Returns:

numpy.ndarray, The regularization matrix \(U\) with shape (xlen * ylen, xlen * ylen).

Raises:

• TypeError – If xlen or ylen are not integers.

• ValueError – If an unsupported regularization_type is provided.

lenstronomy.ImSim.SourceReconstruction.solve_regularization_strength module

d_log_evi_d_lambda(l: float, U: ndarray, M: ndarray, b: ndarray) → float

Computes the derivative of the logarithm of the Bayesian evidence with respect to the regularization strength (lambda, l).

This function calculates the derivative as: d(ln(Evidence))/d(lambda) ~ N_s/lambda - tr[(M+lambda*U)^-1 * U] - b^T * (M+lambda*U)^-1 * U * (M+lambda*U)^-1 * b

Where: - N_s: Number of source pixels (U.shape[0]) - lambda: The regularization strength - U: The regularization matrix - M: The M matrix - b: The b vector

Parameters:

• l – The current value of the regularization strength (lambda).

• U – The regularization matrix (numpy.ndarray).

• M – The M matrix (numpy.ndarray).

• b – The b vector (numpy.ndarray).

Returns:

The computed derivative value (float).

solve_optimal_lambda(derivative_function: Callable[[float, ndarray, ndarray, ndarray], float], U: ndarray, M: ndarray, b: ndarray, initial_lower_bound: float, initial_upper_bound: float, tolerance: float = 1e-07, max_iterations: int = 20, check_initial_bounds: bool = True) → float

Finds the optimal regularization strength (lambda) by solving for the root of the log-evidence derivative using a bisection method.

The optimal lambda is typically the value where the derivative of the log-evidence is zero. This function assumes that the derivative d(ln(Evidence))/d(lambda) is monotonically decreasing and crosses zero within the specified bounds.

Parameters:

• derivative_function – A callable function that computes the derivative d(ln(Evidence))/d(lambda). It must accept (regularization_strength, data_matrix, regularization_matrix, data_vector) as its arguments.

• U – The regularization matrix (numpy.ndarray).

• M – The M matrix (numpy.ndarray).

• b – The b vector (numpy.ndarray).

• initial_lower_bound – The lower bound for the search range of lambda. It is expected that derivative_function(initial_lower_bound, …) > 0.

• initial_upper_bound – The upper bound for the search range of lambda. It is expected that derivative_function(initial_upper_bound, …) < 0.

• tolerance – float, The desired absolute tolerance for the lambda value. The search stops when the width of the search interval is less than this value. Defaults to 1e-7.

• max_iterations – int, The maximum number of bisection iterations to perform. Defaults to 20.

• check_initial_bounds – bool, If True, perform checks to ensure that initial_lower_bound < initial_upper_bound and that the derivative function returns the expected signs at the boundaries (positive at lower, negative at upper). Setting this to False can speed up repeated calls if the bounds are guaranteed to be valid, but disables critical error checking. Defaults to True.

Returns:

float, The optimized regularization strength (lambda) that maximizes the log-evidence.

Raises:

ValueError – If check_initial_bounds is True and initial_lower_bound is not strictly less than initial_upper_bound, or if the derivative function does not yield the expected signs at the initial bounds (i.e., the root is not bracketed).

Module contents