diff --git a/autoarray/fit/fit_dataset.py b/autoarray/fit/fit_dataset.py index 8ef8dfc64..3f392ddcf 100644 --- a/autoarray/fit/fit_dataset.py +++ b/autoarray/fit/fit_dataset.py @@ -325,11 +325,20 @@ def log_likelihood_with_regularization(self) -> float: @property def log_evidence(self) -> float: """ - Returns the log evidence of the inversion's fit to a dataset, where the log evidence includes a number of terms - which quantify the complexity of an inversion's reconstruction (see the `Inversion` module): + Returns the log Bayesian evidence of the inversion's fit to a dataset, which extends the log likelihood by + including penalty terms that quantify the complexity of the inversion's reconstruction: - Log Evidence = -0.5*[Chi_Squared_Term + Regularization_Term + Log(Covariance_Regularization_Term) - - Log(Regularization_Matrix_Term) + Noise_Term] + Log Evidence = -0.5 * [χ² + s^T H s + ln(det(F + H)) - ln(det(H)) + Σ ln(2π σ²)] + + where: + - χ² is the chi-squared goodness-of-fit term + - s^T H s is the regularization term (smoothness penalty on the reconstructed source pixels) + - ln(det(F + H)) penalizes overly complex reconstructions (log determinant of the curvature + regularization matrix) + - ln(det(H)) normalizes the regularization matrix complexity (log determinant of the regularization matrix) + - Σ ln(2π σ²) is the noise normalization term + + This is described in Warren & Dye 2003 (https://arxiv.org/pdf/astro-ph/0302587.pdf) and + Nightingale & Dye 2015 (https://arxiv.org/abs/1708.07377). Returns `None` if no inversion is present, in which case `log_likelihood` is used as the figure of merit. """ diff --git a/autoarray/fit/fit_interferometer.py b/autoarray/fit/fit_interferometer.py index e849b07c0..61051bdc2 100644 --- a/autoarray/fit/fit_interferometer.py +++ b/autoarray/fit/fit_interferometer.py @@ -136,13 +136,20 @@ def noise_normalization(self) -> float: @property def log_evidence(self) -> float: """ - Returns the log evidence of the inversion's fit to a dataset, where the log evidence includes a number of terms - which quantify the complexity of an inversion's reconstruction (see the `Inversion` module): + Returns the log Bayesian evidence of the inversion's fit to a dataset, which extends the log likelihood by + including penalty terms that quantify the complexity of the inversion's reconstruction: - Log Evidence = -0.5*[Chi_Squared_Term + Regularization_Term + Log(Covariance_Regularization_Term) - - Log(Regularization_Matrix_Term) + Noise_Term] + Log Evidence = -0.5 * [χ² + s^T H s + ln(det(F + H)) - ln(det(H)) + Σ ln(2π σ²)] - For interferometer fits the chi-squared uses the fast inversion chi-squared (`inversion.fast_chi_squared`). + where: + - χ² is the chi-squared goodness-of-fit term + - s^T H s is the regularization term (smoothness penalty on the reconstructed source pixels) + - ln(det(F + H)) penalizes overly complex reconstructions (log determinant of the curvature + regularization matrix) + - ln(det(H)) normalizes the regularization matrix complexity (log determinant of the regularization matrix) + - Σ ln(2π σ²) is the noise normalization term + + For interferometer fits the chi-squared uses `inversion.fast_chi_squared`, which avoids computing the full + residual visibilities by evaluating χ² directly from the reconstruction and inversion matrices. Returns `None` if no inversion is present, in which case `log_likelihood` is used as the figure of merit. """ diff --git a/autoarray/fit/fit_util.py b/autoarray/fit/fit_util.py index 01e9a621e..cf210214c 100644 --- a/autoarray/fit/fit_util.py +++ b/autoarray/fit/fit_util.py @@ -404,11 +404,12 @@ def log_evidence_from( noise_normalization: float, ) -> float: """ - Returns the log evidence of an inversion's fit to a dataset, where the log evidence includes a number of terms which - quantify the complexity of an inversion's reconstruction (see the `Inversion` module): + Returns the log Bayesian evidence of an inversion's fit to a dataset, which extends the log likelihood by + including penalty terms that quantify the complexity of the inversion's reconstruction: - Log Evidence = -0.5*[Chi_Squared_Term + Regularization_Term + Log(Covariance_Regularization_Term) - - Log(Regularization_Matrix_Term) + Noise_Term] + Log Evidence = -0.5 * [χ² + s^T H s + ln(det(F + H)) - ln(det(H)) + Σ ln(2π σ²)] + + This is described in Warren & Dye 2003 (https://arxiv.org/pdf/astro-ph/0302587.pdf). Parameters ---------- diff --git a/autoarray/inversion/inversion/abstract.py b/autoarray/inversion/inversion/abstract.py index b45b1984a..abda0ef07 100644 --- a/autoarray/inversion/inversion/abstract.py +++ b/autoarray/inversion/inversion/abstract.py @@ -55,15 +55,16 @@ def __init__( Parameters ---------- - data - The data of the dataset (e.g. the `image` of `Imaging` data) which may have been changed. - noise_map - The noise_map of the noise_mapset (e.g. the `noise_map` of `Imaging` noise_map) which may have been changed. + dataset + The dataset being reconstructed (e.g. an `Imaging` or `Interferometer` dataset, or a `DatasetInterface` + whose attributes like `data` and `noise_map` may have been modified before being passed in). linear_obj_list The list of linear objects (e.g. analytic functions, a mapper with a pixelized grid) which reconstruct the input dataset's data and whose values are solved for via the inversion. settings - Settings controlling how an inversion is fitted for example which linear algebra formalism is used. + Settings controlling how an inversion is fitted, for example which linear algebra formalism is used. + xp + The array module to use (`numpy` by default; pass `jax.numpy` for JAX support). """ self.dataset = dataset @@ -96,11 +97,8 @@ def has(self, cls: Type) -> bool: Parameters ---------- - dict_values - A class dictionary of values which instances of the input `cls` are checked to see if it has an instance - of that class. cls - The type of class that is checked if the object has an instance of. + The type of class whose presence is checked among the linear objects and regularizations in this inversion. """ return misc_util.has( values=self.linear_obj_list + self.regularization_list, cls=cls diff --git a/autoarray/inversion/inversion/dataset_interface.py b/autoarray/inversion/inversion/dataset_interface.py index 4356b7aff..ba731ec21 100644 --- a/autoarray/inversion/inversion/dataset_interface.py +++ b/autoarray/inversion/inversion/dataset_interface.py @@ -36,16 +36,10 @@ def __init__( An array describing the RMS standard deviation error in each pixel used for computing quantities like the chi-squared in a fit (in PyAutoGalaxy and PyAutoLens the recommended units are electrons per second). grids - The grids of (y,x) Cartesian coordinates that the image data is paired with, which are used for evaluting + The grids of (y,x) Cartesian coordinates that the image data is paired with, which are used for evaluating light profiles and calculations associated with a pixelization. - over_sampler - Performs over-sampling whereby the masked image pixels are split into sub-pixels, which are all - mapped via the mapper with sub-fractional values of flux. - border_relocator - The border relocator, which relocates coordinates outside the border of the source-plane data grid to its - edge. psf - Perform 2D convolution of the imaigng data's PSF when computing the operated mapping matrix. + Perform 2D convolution of the imaging data's PSF when computing the operated mapping matrix. transformer Performs a Fourier transform of the image-data from real-space to visibilities when computing the operated mapping matrix. diff --git a/autoarray/inversion/inversion/imaging/abstract.py b/autoarray/inversion/inversion/imaging/abstract.py index ae7dc8f44..03341fb14 100644 --- a/autoarray/inversion/inversion/imaging/abstract.py +++ b/autoarray/inversion/inversion/imaging/abstract.py @@ -50,15 +50,16 @@ def __init__( Parameters ---------- - data - The data of the dataset (e.g. the `image` of `Imaging` data) which may have been changed. - noise_map - The noise_map of the noise_mapset (e.g. the `noise_map` of `Imaging` noise_map) which may have been changed. + dataset + The imaging dataset being reconstructed (e.g. an `Imaging` dataset or a `DatasetInterface` whose + attributes like `data`, `noise_map`, and `psf` may have been modified before being passed in). linear_obj_list The list of linear objects (e.g. analytic functions, a mapper with a pixelized grid) which reconstruct the input dataset's data and whose values are solved for via the inversion. settings - Settings controlling how an inversion is fitted for example which linear algebra formalism is used. + Settings controlling how an inversion is fitted, for example which linear algebra formalism is used. + xp + The array module to use (`numpy` by default; pass `jax.numpy` for JAX support). """ super().__init__( diff --git a/autoarray/inversion/inversion/imaging/inversion_imaging_util.py b/autoarray/inversion/inversion/imaging/inversion_imaging_util.py index 776938805..789478b74 100644 --- a/autoarray/inversion/inversion/imaging/inversion_imaging_util.py +++ b/autoarray/inversion/inversion/imaging/inversion_imaging_util.py @@ -154,7 +154,7 @@ def data_linear_func_matrix_from( The operated values of each linear function divided by the noise-map squared, in a matrix of shape [data_pixels, total_fixed_linear_functions]. kernel_native - The 2D PSf kernel. + The 2D PSF kernel. Returns ------- diff --git a/autoarray/inversion/inversion/interferometer/abstract.py b/autoarray/inversion/inversion/interferometer/abstract.py index 7b37f2c1f..9368eaa4b 100644 --- a/autoarray/inversion/inversion/interferometer/abstract.py +++ b/autoarray/inversion/inversion/interferometer/abstract.py @@ -31,14 +31,16 @@ def __init__( Parameters ---------- - noise_map - The noise-map of the observed interferometer data which values are solved for. - transformer - The transformer which performs a non-uniform fast Fourier transform operations on the mapping matrix - with the interferometer data's transformer. + dataset + The interferometer dataset being reconstructed (e.g. an `Interferometer` dataset or a `DatasetInterface` + whose attributes like `data`, `noise_map`, and `transformer` may have been modified). linear_obj_list The linear objects used to reconstruct the data's observed values. If multiple linear objects are passed the simultaneous linear equations are combined and solved simultaneously. + settings + Settings controlling how an inversion is fitted, for example which linear algebra formalism is used. + xp + The array module to use (`numpy` by default; pass `jax.numpy` for JAX support). """ super().__init__( @@ -124,7 +126,19 @@ def mapped_reconstructed_data_dict( @property def fast_chi_squared(self): + """ + Returns the chi-squared of the interferometer inversion without needing to form the full residual visibilities. + + This is computed directly from the reconstruction and the matrices of the inversion: + + chi_squared = s^T F s - 2 s^T D + sum(d_r^2/sigma_r^2) + sum(d_i^2/sigma_i^2) + where `s` is the reconstruction vector, `F` is the curvature matrix, `D` is the data vector, + and `d_r`/`d_i` are the real/imaginary parts of the observed visibilities. + + This avoids computing the full mapped reconstructed visibilities and is faster than computing + `chi_squared` via the residual visibilities when many source pixels are used. + """ xp = self._xp chi_squared_term_1 = xp.linalg.multi_dot( diff --git a/autoarray/inversion/inversion/interferometer/inversion_interferometer_util.py b/autoarray/inversion/inversion/interferometer/inversion_interferometer_util.py index fef730bb3..e2c8b27df 100644 --- a/autoarray/inversion/inversion/interferometer/inversion_interferometer_util.py +++ b/autoarray/inversion/inversion/interferometer/inversion_interferometer_util.py @@ -20,10 +20,10 @@ def data_vector_via_transformed_mapping_matrix_from( ---------- transformed_mapping_matrix The matrix representing the transformed mappings between sub-grid pixels and pixelization pixels. - image - Flattened 1D array of the observed image the inversion is fitting. + visibilities + The complex 1D array of observed visibilities the inversion is fitting, with real and imaginary components. noise_map - Flattened 1D array of the noise-map used by the inversion during the fit. + The complex 1D array of the noise-map used by the inversion during the fit. """ # Extract components vis_real = visibilities.real diff --git a/autoarray/inversion/inversion/inversion_util.py b/autoarray/inversion/inversion/inversion_util.py index a44311992..7c552c23a 100644 --- a/autoarray/inversion/inversion/inversion_util.py +++ b/autoarray/inversion/inversion/inversion_util.py @@ -92,10 +92,20 @@ def curvature_matrix_via_mapping_matrix_from( Parameters ---------- mapping_matrix - The matrix representing the mappings (these could be blurred or transfomed) between sub-grid pixels and + The matrix representing the mappings (these could be blurred or transformed) between sub-grid pixels and pixelization pixels. noise_map Flattened 1D array of the noise-map used by the inversion during the fit. + add_to_curvature_diag + If `True`, adds a small numerical value to the diagonal entries of the curvature matrix at the indices + specified by `no_regularization_index_list`, to ensure the matrix is positive-definite. + no_regularization_index_list + A list of parameter indices which have no regularization applied. A small value is added to their diagonal + entries if `add_to_curvature_diag` is `True`. + settings + Settings controlling how the inversion is performed, for example whether mixed precision is used. + xp + The array module to use (`numpy` by default; pass `jax.numpy` for JAX support). """ # NumPy path: keep it simple + stable if xp is np: diff --git a/autoarray/inversion/linear_obj/linear_obj.py b/autoarray/inversion/linear_obj/linear_obj.py index a4cc450d4..ec11ea0e6 100644 --- a/autoarray/inversion/linear_obj/linear_obj.py +++ b/autoarray/inversion/linear_obj/linear_obj.py @@ -14,7 +14,7 @@ def __init__(self, regularization: Optional[AbstractRegularization], xp=np): The values of linear object are computed via a regularized linear matrix inversion, which infers a solution which best fits the data given its noise-map (by minimizing a chi-squared term) whilst accounting for smoothing - due to regularizaiton. + due to regularization. For example, in `PyAutoGalaxy` and `PyAutoLens` the light of galaxies is represented using `LightProfile` objects, which describe the surface brightness of a galaxy as a function. This function can either be assigned @@ -46,7 +46,7 @@ def params(self) -> int: For example for the following linear objects: - `AbstractLinearObjFuncList` the number of analytic functions. - - `Mapper` the number of piels in the mesh used to reconstruct the data. + - `Mapper` the number of pixels in the mesh used to reconstruct the data. Returns ------- diff --git a/autoarray/inversion/mappers/abstract.py b/autoarray/inversion/mappers/abstract.py index d6d832b89..f94530109 100644 --- a/autoarray/inversion/mappers/abstract.py +++ b/autoarray/inversion/mappers/abstract.py @@ -30,11 +30,11 @@ def __init__( the four grids are explained (`image_plane_data_grid`, `source_plane_data_grid`, `image_plane_mesh_grid`,`source_plane_mesh_grid`) - If you are unfamliar withe above objects, read through the docstrings of the `pixelization`, `mesh` and + If you are unfamiliar with the above objects, read through the docstrings of the `pixelization`, `mesh` and `image_mesh` packages. A `Mapper` determines the mappings between the masked data grid's pixels (`image_plane_data_grid` and - `source_plane_data_grid`) and the pxelization's pixels (`image_plane_mesh_grid` and `source_plane_mesh_grid`). + `source_plane_data_grid`) and the pixelization's pixels (`image_plane_mesh_grid` and `source_plane_mesh_grid`). The 1D Indexing of each grid is identical in the `data` and `source` frames (e.g. the transformation does not change the indexing, such that `source_plane_data_grid[0]` corresponds to the transformed value @@ -66,28 +66,25 @@ def __init__( pixelization's 6th pixel (index 5). The mapper allows us to create a mapping matrix, which is a matrix representing the mapping between every - unmasked data pixel annd the pixels of a pixelization. This matrix is the basis of performing an `Inversion`, + unmasked data pixel and the pixels of a pixelization. This matrix is the basis of performing an `Inversion`, which reconstructs the data using the `source_plane_mesh_grid`. Parameters ---------- - source_plane_data_grid - A 2D grid of (y,x) coordinates associated with the unmasked 2D data after it has been transformed to the - `source` reference frame. - source_plane_mesh_grid - The 2D grid of (y,x) centres of every pixelization pixel in the `source` frame. - image_plane_mesh_grid - The sparse set of (y,x) coordinates computed from the unmasked data in the `data` frame. This has a - transformation applied to it to create the `source_plane_mesh_grid`. - adapt_data - An image which is used to determine the `image_plane_mesh_grid` and therefore adapt the distribution of - pixels of the Delaunay grid to the data it discretizes. - mesh_weight_map - The weight map used to weight the creation of the rectangular mesh grid, which is used for the - `RectangularBrightness` mesh which adapts the size of its pixels to where the source is reconstructed. + interpolator + The interpolator which computes the mappings between image-plane data pixels and source-plane mesh pixels, + including the `source_plane_data_grid`, `source_plane_mesh_grid`, interpolation weights, and + `adapt_data`. regularization The regularization scheme which may be applied to this linear object in order to smooth its solution, which for a mapper smooths neighboring pixels on the mesh. + settings + Settings controlling how an inversion is fitted, for example which linear algebra formalism is used. + image_plane_mesh_grid + The sparse set of (y,x) coordinates computed from the unmasked data in the image frame. This is + transformed to create the `source_plane_mesh_grid` accessed via `self.source_plane_mesh_grid`. + xp + The array module to use (`numpy` by default; pass `jax.numpy` for JAX support). """ super().__init__(regularization=regularization, xp=xp)