Feature/numba speed up

autoarray/dataset/grids.py

@@ -17,6 +17,7 @@ def __init__(
         over_sample_size_lp: Union[int, Array2D],
         over_sample_size_pixelization: Union[int, Array2D],
         psf: Optional[Kernel2D] = None,
+        use_w_tilde: bool = False,
     ):
         """
         Contains grids of (y,x) Cartesian coordinates at the centre of every pixel in the dataset's image and
@@ -64,13 +65,11 @@ def __init__(
             mask=self.mask,
             over_sample_size=self.over_sample_size_lp,
         )
-        self.lp.over_sampled
 
         self.pixelization = Grid2D.from_mask(
             mask=self.mask,
             over_sample_size=self.over_sample_size_pixelization,
         )
-        self.pixelization.over_sampled
 
         if self.psf is None:
             self.blurring = None
@@ -79,12 +78,13 @@ def __init__(
                 self.blurring = self.lp.blurring_grid_via_kernel_shape_from(
                     kernel_shape_native=self.psf.shape_native,
                 )
-                self.blurring.over_sampled
             except exc.MaskException:
                 self.blurring = None
 
         self.border_relocator = BorderRelocator(
-            mask=self.mask, sub_size=self.over_sample_size_pixelization
+            mask=self.mask,
+            sub_size=self.over_sample_size_pixelization,
+            use_w_tilde=use_w_tilde,
         )
 
 

autoarray/dataset/imaging/dataset.py

@@ -191,11 +191,14 @@ def __init__(
             if psf.mask.shape[0] % 2 == 0 or psf.mask.shape[1] % 2 == 0:
                 raise exc.KernelException("Kernel2D Kernel2D must be odd")
 
+        use_w_tilde = True if w_tilde is not None else False
+
         self.grids = GridsDataset(
             mask=self.data.mask,
             over_sample_size_lp=self.over_sample_size_lp,
             over_sample_size_pixelization=self.over_sample_size_pixelization,
             psf=self.psf,
+            use_w_tilde=use_w_tilde,
         )
 
         self.w_tilde = w_tilde

autoarray/dataset/interferometer/dataset.py

@@ -100,10 +100,13 @@ def __init__(
 
         self.dft_preload_transform = dft_preload_transform
 
+        use_w_tilde = True if w_tilde is not None else False
+
         self.grids = GridsDataset(
             mask=self.real_space_mask,
             over_sample_size_lp=self.over_sample_size_lp,
             over_sample_size_pixelization=self.over_sample_size_pixelization,
+            use_w_tilde=use_w_tilde,
         )
 
         self.w_tilde = w_tilde

autoarray/inversion/inversion/abstract.py

@@ -293,7 +293,7 @@ def mapping_matrix(self) -> np.ndarray:
     def operated_mapping_matrix_list(self) -> np.ndarray:
         raise NotImplementedError
 
-    @property
+    @cached_property
     def operated_mapping_matrix(self) -> np.ndarray:
         """
         The `operated_mapping_matrix` of a linear object describes the mappings between the observed data's values and
@@ -314,7 +314,7 @@ def data_vector(self) -> np.ndarray:
     def curvature_matrix(self) -> np.ndarray:
         raise NotImplementedError
 
-    @property
+    @cached_property
     def regularization_matrix(self) -> Optional[np.ndarray]:
         """
         The regularization matrix H is used to impose smoothness on our inversion's reconstruction. This enters the
@@ -346,7 +346,7 @@ def regularization_matrix(self) -> Optional[np.ndarray]:
             *[linear_obj.regularization_matrix for linear_obj in self.linear_obj_list]
         )
 
-    @property
+    @cached_property
     def regularization_matrix_reduced(self) -> Optional[np.ndarray]:
         """
         The regularization matrix H is used to impose smoothness on our inversion's reconstruction. This enters the
@@ -360,7 +360,6 @@ def regularization_matrix_reduced(self) -> Optional[np.ndarray]:
         The scipy function `block_diag` has an overhead associated with it and if there is only one mapper and
         regularization it is bypassed.
         """
-
         if self.all_linear_obj_have_regularization:
             return self.regularization_matrix
 
@@ -380,7 +379,6 @@ def curvature_reg_matrix(self) -> np.ndarray:
         to ensure if we access it after computing the `curvature_reg_matrix` it is correctly recalculated in a new
         array of memory.
         """
-
         if not self.has(cls=AbstractRegularization):
             return self.curvature_matrix
 
@@ -400,7 +398,6 @@ def curvature_reg_matrix_reduced(self) -> Optional[np.ndarray]:
         The scipy function `block_diag` has an overhead associated with it and if there is only one mapper and
         regularization it is bypassed.
         """
-
         if self.all_linear_obj_have_regularization:
             return self.curvature_reg_matrix
 
@@ -426,7 +423,6 @@ def reconstruction(self) -> np.ndarray:
             ZTZ := np.dot(Z.T, Z)
             ZTx := np.dot(Z.T, x)
         """
-
         if self.settings.use_positive_only_solver:
 
             if (
@@ -481,15 +477,14 @@ def reconstruction(self) -> np.ndarray:
             xp=self._xp,
         )
 
-    @property
+    @cached_property
     def reconstruction_reduced(self) -> np.ndarray:
         """
         Solve the linear system [F + reg_coeff*H] S = D -> S = [F + reg_coeff*H]^-1 D given by equation (12)
         of https://arxiv.org/pdf/astro-ph/0302587.pdf
 
         S is the vector of reconstructed inversion values.
         """
-
         if self.all_linear_obj_have_regularization:
             return self.reconstruction
 

autoarray/inversion/inversion/imaging/inversion_imaging_numba_util.py

@@ -857,3 +857,71 @@ def mapped_reconstructed_data_via_image_to_pix_unique_from(
             )
 
     return mapped_reconstructed_data
+
+
+@numba_util.jit()
+def relocated_grid_via_jit_from(grid, border_grid):
+    """
+    Relocate the coordinates of a grid to its border if they are outside the border, where the border is
+    defined as all pixels at the edge of the grid's mask (see *mask._border_1d_indexes*).
+
+    This is performed as follows:
+
+    1: Use the mean value of the grid's y and x coordinates to determine the origin of the grid.
+    2: Compute the radial distance of every grid coordinate from the origin.
+    3: For every coordinate, find its nearest pixel in the border.
+    4: Determine if it is outside the border, by comparing its radial distance from the origin to its paired
+    border pixel's radial distance.
+    5: If its radial distance is larger, use the ratio of radial distances to move the coordinate to the
+    border (if its inside the border, do nothing).
+
+    The method can be used on uniform or irregular grids, however for irregular grids the border of the
+    'image-plane' mask is used to define border pixels.
+
+    Parameters
+    ----------
+    grid
+        The grid (uniform or irregular) whose pixels are to be relocated to the border edge if outside it.
+    border_grid : Grid2D
+        The grid of border (y,x) coordinates.
+    """
+
+    grid_relocated = np.zeros(grid.shape)
+    grid_relocated[:, :] = grid[:, :]
+
+    border_origin = np.zeros(2)
+    border_origin[0] = np.mean(border_grid[:, 0])
+    border_origin[1] = np.mean(border_grid[:, 1])
+    border_grid_radii = np.sqrt(
+        np.add(
+            np.square(np.subtract(border_grid[:, 0], border_origin[0])),
+            np.square(np.subtract(border_grid[:, 1], border_origin[1])),
+        )
+    )
+    border_min_radii = np.min(border_grid_radii)
+
+    grid_radii = np.sqrt(
+        np.add(
+            np.square(np.subtract(grid[:, 0], border_origin[0])),
+            np.square(np.subtract(grid[:, 1], border_origin[1])),
+        )
+    )
+
+    for pixel_index in range(grid.shape[0]):
+        if grid_radii[pixel_index] > border_min_radii:
+            closest_pixel_index = np.argmin(
+                np.square(grid[pixel_index, 0] - border_grid[:, 0])
+                + np.square(grid[pixel_index, 1] - border_grid[:, 1])
+            )
+
+            move_factor = (
+                border_grid_radii[closest_pixel_index] / grid_radii[pixel_index]
+            )
+
+            if move_factor < 1.0:
+                grid_relocated[pixel_index, :] = (
+                    move_factor * (grid[pixel_index, :] - border_origin[:])
+                    + border_origin[:]
+                )
+
+    return grid_relocated

autoarray/inversion/inversion/imaging/mapping.py

@@ -1,6 +1,8 @@
 import numpy as np
 from typing import Dict, List, Optional, Union
 
+from autoconf import cached_property
+
 from autoarray.dataset.imaging.dataset import Imaging
 from autoarray.inversion.inversion.dataset_interface import DatasetInterface
 from autoarray.inversion.inversion.imaging.abstract import AbstractInversionImaging
@@ -89,7 +91,7 @@ def _data_vector_mapper(self) -> np.ndarray:
 
         return data_vector
 
-    @property
+    @cached_property
     def data_vector(self) -> np.ndarray:
         """
         The `data_vector` is a 1D vector whose values are solved for by the simultaneous linear equations constructed
@@ -156,7 +158,7 @@ def _curvature_matrix_mapper_diag(self) -> Optional[np.ndarray]:
 
         return curvature_matrix
 
-    @property
+    @cached_property
     def curvature_matrix(self):
         """
         The `curvature_matrix` is a 2D matrix which uses the mappings between the data and the linear objects to

autoarray/inversion/inversion/imaging/w_tilde.py

@@ -1,6 +1,8 @@
 import numpy as np
 from typing import Dict, List, Optional, Union
 
+from autoconf import cached_property
+
 from autoarray.dataset.imaging.dataset import Imaging
 from autoarray.dataset.imaging.w_tilde import WTildeImaging
 from autoarray.inversion.inversion.dataset_interface import DatasetInterface
@@ -64,9 +66,8 @@ def __init__(
 
         self.w_tilde = dataset.w_tilde
 
-    @property
+    @cached_property
     def w_tilde_data(self):
-
         return inversion_imaging_numba_util.w_tilde_data_imaging_from(
             image_native=np.array(self.data.native.array),
             noise_map_native=self.noise_map.native.array,
@@ -110,7 +111,7 @@ def _data_vector_mapper(self) -> np.ndarray:
 
         return data_vector
 
-    @property
+    @cached_property
     def data_vector(self) -> np.ndarray:
         """
         Returns the `data_vector`, a 1D vector whose values are solved for by the simultaneous linear equations
@@ -210,7 +211,7 @@ def _data_vector_func_list_and_mapper(self) -> np.ndarray:
 
         return data_vector
 
-    @property
+    @cached_property
     def curvature_matrix(self) -> np.ndarray:
         """
         Returns the `curvature_matrix`, a 2D matrix which uses the mappings between the data and the linear objects to
@@ -231,7 +232,6 @@ def curvature_matrix(self) -> np.ndarray:
         to ensure if we access it after computing the `curvature_reg_matrix` it is correctly recalculated in a new
         array of memory.
         """
-
         if self.has(cls=AbstractLinearObjFuncList):
             curvature_matrix = self._curvature_matrix_func_list_and_mapper
         elif self.total(cls=AbstractMapper) == 1:

autoarray/inversion/linear_obj/func_list.py

@@ -1,6 +1,8 @@
 import numpy as np
 from typing import Optional
 
+from autoconf import cached_property
+
 from autoarray.inversion.linear_obj.linear_obj import LinearObj
 from autoarray.inversion.linear_obj.neighbors import Neighbors
 from autoarray.inversion.linear_obj.unique_mappings import UniqueMappings
@@ -44,7 +46,7 @@ def __init__(
 
         self.grid = grid
 
-    @property
+    @cached_property
     def neighbors(self) -> Neighbors:
         """
         An object describing how the different parameters in the linear object neighbor one another, which is used
@@ -76,7 +78,7 @@ def neighbors(self) -> Neighbors:
             arr=neighbors.astype("int"), sizes=neighbors_sizes.astype("int")
         )
 
-    @property
+    @cached_property
     def unique_mappings(self) -> UniqueMappings:
         """
         Returns the unique mappings of every unmasked data pixel's (e.g. `grid_slim`) sub-pixels (e.g. `grid_sub_slim`)
@@ -88,7 +90,6 @@ def unique_mappings(self) -> UniqueMappings:
 
         For a `LinearObjFuncList` every data pixel's group of sub-pixels maps directly to the linear function.
         """
-
         sub_size = np.max(self.grid.over_sample_size)
 
         # TODO : This shape slim is prob unreliable and needs to be divided by sub_size**2