Feature/mapper edge

autoarray/config/general.yaml

@@ -1,6 +1,5 @@
 inversion:
   reconstruction_vmax_factor: 0.5   # Plots of an Inversion's reconstruction use the reconstructed data's bright value multiplied by this factor.
-  check_solution_default: true      # When an Inversion checks solutions in case they are numerically unstable by default or not (e.g. InversionSettings.check_solution).
 numba:
   use_numba: true
   cache: false

autoarray/config/visualize/general.yaml

@@ -4,6 +4,5 @@ general:
   zoom_around_mask: true            # If True, plots of data structures with a mask automatically zoom in the masked region.
 inversion:
   reconstruction_vmax_factor: 0.5   # Plots of an Inversion's reconstruction use the reconstructed data's bright value multiplied by this factor.
-  check_solution_default: true      # When an Inversion checks solutions in case they are numerically unstable by default or not (e.g. InversionSettings.check_solution).
 units:
   in_kpc: false                     # If True, plots that are normally in arc-seconds are instead plotted using kiloparsecs.

autoarray/inversion/inversion/abstract.py

@@ -1,3 +1,5 @@
+import copy
+
 import numpy as np
 from scipy.linalg import block_diag
 from scipy.sparse import csc_matrix
@@ -223,6 +225,33 @@ def all_linear_obj_have_regularization(self) -> bool:
             list(filter(None, self.regularization_list))
         )
 
+    @property
+    def mapper_edge_pixel_list(self) -> List[int]:
+        """
+        Returns the edge pixels of all mappers in the inversion.
+
+        This uses the `edge_pixel_list` property of the `Mesh` of the `Mapper` class, and updates their values to
+        correspond to the indexing of the overall inversion's `curvature_matrix`.
+
+        This is used to regulareze the edge pixels of the inversion's `reconstruction` or remove them from the
+        inversion procedure entirely (e.g. make these values of these edge pixels zero).
+
+        Returns
+        -------
+        A list of the edge pixels of all mappers in the inversion, where the values are updated to correspond to the
+        indexing of the overall inversion's `curvature_matrix`.
+        """
+        mapper_edge_pixel_list = []
+
+        param_range_list = self.param_range_list_from(cls=LinearObj)
+
+        for param_range, linear_obj in zip(param_range_list, self.linear_obj_list):
+            if isinstance(linear_obj, AbstractMapper):
+                for edge_pixel in linear_obj.edge_pixel_list:
+                    mapper_edge_pixel_list.append(edge_pixel + param_range[0])
+
+        return mapper_edge_pixel_list
+
     @property
     def total_regularizations(self) -> int:
         return sum(
@@ -323,6 +352,9 @@ def regularization_matrix(self) -> Optional[np.ndarray]:
         For multiple mappers, the regularization matrix is computed as the block diagonal of each individual mapper.
         The scipy function `block_diag` has an overhead associated with it and if there is only one mapper and
         regularization it is bypassed.
+
+        If the `settings.force_edge_pixels_to_zeros` is `True`, the edge pixels of each mapper in the inversion
+        are regularized so high their value is forced to zero.
         """
         if self.preloads.regularization_matrix is not None:
             return self.preloads.regularization_matrix
@@ -375,8 +407,6 @@ def curvature_reg_matrix(self) -> np.ndarray:
         if not self.has(cls=AbstractRegularization):
             return self.curvature_matrix
 
-        # TODO : Done for speed, instead check if there is one regularization
-
         if len(self.regularization_list) == 1:
 
             curvature_matrix = self.curvature_matrix
@@ -388,6 +418,19 @@ def curvature_reg_matrix(self) -> np.ndarray:
 
         return np.add(self.curvature_matrix, self.regularization_matrix)
 
+    @cached_property
+    def curvature_reg_matrix_solver(self):
+
+        if self.settings.force_edge_pixels_to_zeros:
+
+            curvature_reg_matrix_solver = copy.copy(self.curvature_reg_matrix)
+
+            curvature_reg_matrix_solver[:, self.mapper_edge_pixel_list] = 0.0
+
+            return curvature_reg_matrix_solver
+
+        return self.curvature_reg_matrix
+
     @cached_property
     @profile_func
     def curvature_reg_matrix_reduced(self) -> np.ndarray:
@@ -411,36 +454,6 @@ def curvature_reg_matrix_reduced(self) -> np.ndarray:
 
         return curvature_reg_matrix
 
-    @cached_property
-    @profile_func
-    def curvature_reg_matrix_cholesky(self) -> np.ndarray:
-        """
-        Performs a Cholesky decomposition of the `curvature_reg_matrix`, the result of which is used to solve the
-        linear system of equations of the `Inversion`.
-
-        The method `np.linalg.solve` is faster to do this, but the Cholesky decomposition is used later in the code
-        to speed up the calculation of `log_det_curvature_reg_matrix_term`.
-        """
-        try:
-            return np.linalg.cholesky(self.curvature_reg_matrix)
-        except np.linalg.LinAlgError:
-            raise exc.InversionException()
-
-    @cached_property
-    @profile_func
-    def curvature_reg_matrix_reduced_cholesky(self) -> np.ndarray:
-        """
-        Performs a Cholesky decomposition of the `curvature_reg_matrix`, the result of which is used to solve the
-        linear system of equations of the `Inversion`.
-
-        The method `np.linalg.solve` is faster to do this, but the Cholesky decomposition is used later in the code
-        to speed up the calculation of `log_det_curvature_reg_matrix_term`.
-        """
-        try:
-            return np.linalg.cholesky(self.curvature_reg_matrix_reduced)
-        except np.linalg.LinAlgError:
-            raise exc.InversionException()
-
     @cached_property
     @profile_func
     def reconstruction(self) -> np.ndarray:
@@ -453,12 +466,13 @@ def reconstruction(self) -> np.ndarray:
         if not self.settings.use_positive_only_solver:
             return inversion_util.reconstruction_positive_negative_from(
                 data_vector=self.data_vector,
-                curvature_reg_matrix_cholesky=self.curvature_reg_matrix_cholesky,
+                curvature_reg_matrix=self.curvature_reg_matrix_solver,
                 settings=self.settings,
             )
+
         return inversion_util.reconstruction_positive_only_from(
             data_vector=self.data_vector,
-            curvature_reg_matrix=self.curvature_reg_matrix,
+            curvature_reg_matrix=self.curvature_reg_matrix_solver,
             settings=self.settings,
         )
 
@@ -610,7 +624,12 @@ def log_det_curvature_reg_matrix_term(self) -> float:
         if not self.has(cls=AbstractRegularization):
             return 0.0
 
-        return 2.0 * np.sum(np.log(np.diag(self.curvature_reg_matrix_reduced_cholesky)))
+        try:
+            return 2.0 * np.sum(
+                np.log(np.diag(np.linalg.cholesky(self.curvature_reg_matrix_reduced)))
+            )
+        except np.linalg.LinAlgError:
+            raise exc.InversionException()
 
     @cached_property
     @profile_func

autoarray/inversion/inversion/factory.py

@@ -227,7 +227,7 @@ def inversion_imaging_unpacked_from(
     An `Inversion` whose type is determined by the input `dataset` and `settings`.
     """
 
-    if all(
+    if any(
         isinstance(linear_obj, AbstractLinearObjFuncList)
         for linear_obj in linear_obj_list
     ):

autoarray/inversion/inversion/inversion_util.py

@@ -257,7 +257,7 @@ def mapped_reconstructed_data_via_mapping_matrix_from(
 
 def reconstruction_positive_negative_from(
     data_vector: np.ndarray,
-    curvature_reg_matrix_cholesky: np.ndarray,
+    curvature_reg_matrix: np.ndarray,
     settings: SettingsInversion = SettingsInversion(),
 ):
     """
@@ -278,8 +278,8 @@ def reconstruction_positive_negative_from(
     ----------
     data_vector
         The `data_vector` D which is solved for.
-    curvature_reg_matrix_cholesky
-        The cholesky decomposition of the sum of the curvature and regularization matrices.
+    curvature_reg_matrix
+        The sum of the curvature and regularization matrices.
     settings
         Controls the settings of the inversion, for this function where the solution is checked to not be all
         the same values.
@@ -290,14 +290,10 @@ def reconstruction_positive_negative_from(
         The curvature_matrix plus regularization matrix, overwriting the curvature_matrix in memory.
     """
     try:
-        reconstruction = cho_solve((curvature_reg_matrix_cholesky, True), data_vector)
+        reconstruction = np.linalg.solve(curvature_reg_matrix, data_vector)
     except np.linalg.LinAlgError as e:
         raise exc.InversionException() from e
 
-    reconstruction_check_solution(
-        data_vector=data_vector, reconstruction=reconstruction, settings=settings
-    )
-
     return reconstruction
 
 
@@ -324,7 +320,7 @@ def reconstruction_positive_only_from(
     ----------
     data_vector
         The `data_vector` D which is solved for.
-    curvature_reg_matrix_cholesky
+    curvature_reg_matrix
         The sum of the curvature and regularization matrices.
     settings
         Controls the settings of the inversion, for this function where the solution is checked to not be all
@@ -347,48 +343,9 @@ def reconstruction_positive_only_from(
     except (RuntimeError, np.linalg.LinAlgError) as e:
         raise exc.InversionException() from e
 
-    reconstruction_check_solution(
-        data_vector=data_vector, reconstruction=reconstruction, settings=settings
-    )
-
     return reconstruction
 
 
-def reconstruction_check_solution(
-    data_vector: np.ndarray,
-    reconstruction: np.ndarray,
-    settings: SettingsInversion = SettingsInversion(),
-):
-    """
-    Check that the solution after reconstructing an inversion is not one where all reconstructed values go to the same
-    value.
-
-    These solutions occur when the linear system is numerically unstable, and can lead to high `log_evidence` values
-    due to numerically issues not currently fully understood.
-
-    Parameters
-    ----------
-    data_vector
-        The `data_vector` D which is solved for.
-    reconstruction
-        The solution of the linear system whose solution is checked.
-    settings
-        Controls the settings of the inversion, for this function where the solution is checked to not be all
-        the same values.
-    """
-    if settings.check_solution and len(reconstruction) > 1:
-        if np.isclose(a=reconstruction[0], b=reconstruction[1], atol=1e-6).all():
-            raise exc.InversionException()
-
-    index = int(data_vector.shape[0] / 2)
-
-    if settings.check_solution and len(reconstruction) > 1:
-        if np.isclose(
-            a=reconstruction[index], b=reconstruction[index + 1], atol=1e-6
-        ).all():
-            raise exc.InversionException()
-
-
 def preconditioner_matrix_via_mapping_matrix_from(
     mapping_matrix: np.ndarray,
     regularization_matrix: np.ndarray,

autoarray/inversion/inversion/settings.py

@@ -12,8 +12,8 @@ def __init__(
         use_w_tilde: bool = True,
         use_positive_only_solver: bool = False,
         positive_only_maxiter: int = 5000,
+        force_edge_pixels_to_zeros: bool = False,
         no_regularization_add_to_curvature_diag: bool = True,
-        check_solution: Optional[bool] = None,
         use_w_tilde_numpy: bool = False,
         use_source_loop: bool = False,
         use_linear_operators: bool = False,
@@ -40,10 +40,6 @@ def __init__(
         no_regularization_add_to_curvature_diag
             When True, if a linear object in the inversion has no regularization, values of 1.0e-8 are added to the
             diagonal of its `curvature_matrix` to stablelize the linear algebra solver.
-        check_solution
-            If True, the `reconstruction` of the inversion is checked to ensure that no two source pixels have
-            numerically identical values, which indicates a spurious solution where the linear algebra solver
-            reconstructs every value as an identical value (for reasons that are currently not understood).
         use_w_tilde_numpy
             If True, the curvature_matrix is computed via numpy matrix multiplication (as opposed to numba functions
             which exploit sparsity to do the calculation normally in a more efficient way).
@@ -64,21 +60,11 @@ def __init__(
         self.use_positive_only_solver = use_positive_only_solver
         self.positive_only_maxiter = positive_only_maxiter
         self.use_linear_operators = use_linear_operators
+        self.force_edge_pixels_to_zeros = force_edge_pixels_to_zeros
         self.no_regularization_add_to_curvature_diag = (
             no_regularization_add_to_curvature_diag
         )
         self.tolerance = tolerance
         self.maxiter = maxiter
         self.use_w_tilde_numpy = use_w_tilde_numpy
         self.use_source_loop = use_source_loop
-
-        self._check_solution = check_solution
-
-    @property
-    def check_solution(self):
-
-        if self._check_solution is None:
-
-            return conf.instance["general"]["inversion"]["check_solution_default"]
-
-        return self._check_solution

autoarray/inversion/mock/mock_mapper.py

@@ -11,6 +11,7 @@ def __init__(
         source_plane_data_grid=None,
         source_plane_mesh_grid=None,
         hyper_data=None,
+        edge_pixel_list=None,
         regularization=None,
         pix_sub_weights=None,
         mapping_matrix=None,
@@ -27,6 +28,8 @@ def __init__(
 
         super().__init__(mapper_grids=mapper_grids, regularization=regularization)
 
+        self._edge_pixel_list = edge_pixel_list
+
         self._pix_sub_weights = pix_sub_weights
 
         self._mapping_matrix = mapping_matrix
@@ -48,6 +51,10 @@ def params(self):
             return super().params
         return self._parameters
 
+    @property
+    def edge_pixel_list(self):
+        return self._edge_pixel_list
+
     @property
     def pix_sub_weights(self):
         return self._pix_sub_weights

autoarray/inversion/pixelization/mappers/abstract.py

@@ -104,6 +104,10 @@ def source_plane_mesh_grid(self) -> Abstract2DMesh:
     def image_plane_mesh_grid(self) -> Grid2D:
         return self.mapper_grids.image_plane_mesh_grid
 
+    @property
+    def edge_pixel_list(self) -> List[int]:
+        return self.source_plane_mesh_grid.edge_pixel_list
+
     @property
     def hyper_data(self) -> np.ndarray:
         return self.mapper_grids.hyper_data

autoarray/inversion/pixelization/mappers/mapper_util.py

@@ -152,7 +152,7 @@ def pix_indexes_for_sub_slim_index_delaunay_from(
     delaunay_points,
 ) -> Tuple[np.ndarray, np.ndarray]:
     """
-    The indexes mappings between the sub pixels and Voronoi pixelization pixels.
+    The indexes mappings between the sub pixels and Voronoi mesh pixels.
     For Delaunay tessellation, most sub pixels should have contribution of 3 pixelization pixels. However,
     for those ones not belonging to any triangle, we link its value to its closest point.
 

autoarray/inversion/pixelization/mappers/voronoi.py

@@ -136,10 +136,10 @@ def pix_sub_weights(self) -> PixSubWeights:
         neighbor of the Voronoi mesh's third (index 2) pixel.
 
         - `pix_indexes_for_sub_slim_index[0, 1] = 5`: The data's first (index 0) sub-pixel also maps to the natural
-        neighbor of the Voronoi pixelization's sixth (index 5) pixel.
+        neighbor of the Voronoi mesh's sixth (index 5) pixel.
 
         - `pix_indexes_for_sub_slim_index[0, 2] = 8`: The data's first (index 0) sub-pixel also maps to the natural
-        neighbor of the Voronoi pixelization's ninth (index 8) pixel.
+        neighbor of the Voronoi mesh's ninth (index 8) pixel.
 
         The interpolation weights of these multiple mappings are stored in the array `pix_weights_for_sub_slim_index`.
         """
@@ -222,7 +222,7 @@ def pix_sub_weights(self) -> PixSubWeights:
         where each data pixel maps to 3 Delaunay triangles with interpolation weights). The weights of multiple mappings
         are stored in the array `pix_weights_for_sub_slim_index`.
 
-        For this Voronoi pixelization each data sub-pixel maps to a single pixelization pixel, thus the second
+        For this Voronoi mesh each data sub-pixel maps to a single pixelization pixel, thus the second
         dimension of the array `pix_indexes_for_sub_slim_index` 1 and all entries in `pix_weights_for_sub_slim_index`
         are equal to 1.0.