diff --git a/autoarray/__init__.py b/autoarray/__init__.py index dd2cd4476..7a32c6477 100644 --- a/autoarray/__init__.py +++ b/autoarray/__init__.py @@ -14,7 +14,6 @@ from .dataset.interferometer.dataset import Interferometer from .dataset.interferometer.simulator import SimulatorInterferometer from .dataset.interferometer.w_tilde import WTildeInterferometer -from .dataset.over_sampling import OverSamplingDataset from .dataset.dataset_model import DatasetModel from .fit.fit_dataset import AbstractFit from .fit.fit_dataset import FitDataset @@ -68,14 +67,8 @@ from .structures.arrays.irregular import ArrayIrregular from .structures.grids.uniform_1d import Grid1D from .structures.grids.uniform_2d import Grid2D -from .operators.over_sampling.decorator import perform_over_sampling_from -from .operators.over_sampling.grid_oversampled import Grid2DOverSampled -from .operators.over_sampling.uniform import OverSamplingUniform -from .operators.over_sampling.iterate import OverSamplingIterate -from .operators.over_sampling.uniform import OverSamplerUniform -from .operators.over_sampling.iterate import OverSamplerIterate +from .operators.over_sampling.over_sampler import OverSampler from .structures.grids.irregular_2d import Grid2DIrregular -from .structures.grids.irregular_2d import Grid2DIrregularUniform from .structures.mesh.rectangular_2d import Mesh2DRectangular from .structures.mesh.voronoi_2d import Mesh2DVoronoi from .structures.mesh.delaunay_2d import Mesh2DDelaunay @@ -98,4 +91,4 @@ conf.instance.register(__file__) -__version__ = "2024.11.6.1" +__version__ = "2025.1.18.7" diff --git a/autoarray/config/visualize/mat_wrap_2d.yaml b/autoarray/config/visualize/mat_wrap_2d.yaml index 6ca78d967..3cd6f4b4e 100644 --- a/autoarray/config/visualize/mat_wrap_2d.yaml +++ b/autoarray/config/visualize/mat_wrap_2d.yaml @@ -146,6 +146,15 @@ VectorYXQuiver: # wrapper for `plt.quiver()`: customize (y,x) vectors appea linewidth: 5 pivot: middle units: xy +DelaunayDrawer: # wrapper for `plt.fill()`: customize the appearance of Delaunay mesh's. + figure: + alpha: 0.7 + edgecolor: k + linewidth: 0.0 + subplot: + alpha: 0.7 + edgecolor: k + linewidth: 0.0 VoronoiDrawer: # wrapper for `plt.fill()`: customize the appearance of Voronoi mesh's. figure: alpha: 0.7 diff --git a/autoarray/config/visualize/plots.yaml b/autoarray/config/visualize/plots.yaml index ce0137511..aa2664506 100644 --- a/autoarray/config/visualize/plots.yaml +++ b/autoarray/config/visualize/plots.yaml @@ -5,20 +5,11 @@ dataset: # Settings for plots of all datasets (e.g. ImagingPlotter, InterferometerPlotter). subplot_dataset: true # Plot subplot containing all dataset quantities (e.g. the data, noise-map, etc.)? - data: false # Plot the individual data of every dataset? - noise_map: false # Plot the individual noise-map of every dataset? - signal_to_noise_map: false # Plot the individual signal-to-noise-map of every dataset? - over_sampling: false # Plot the over-sampling sub-size, used to evaluate light profiles, of every dataset? - over_sampling_non_uniform: false # Plot the over-sampling sub-size, used to evaluate non uniform grids, of every dataset? - over_sampling_pixelization: false # Plot the over-sampling sub-size, used to evaluate pixelizations, of every dataset? imaging: # Settings for plots of imaging datasets (e.g. ImagingPlotter) psf: false fit: # Settings for plots of all fits (e.g. FitImagingPlotter, FitInterferometerPlotter). subplot_fit: true # Plot subplot of all fit quantities for any dataset (e.g. the model data, residual-map, etc.)? subplot_fit_log10: true # Plot subplot of all fit quantities for any dataset using log10 color maps (e.g. the model data, residual-map, etc.)? - all_at_end_png: true # Plot all individual plots listed below as .png (even if False)? - all_at_end_fits: true # Plot all individual plots listed below as .fits (even if False)? - all_at_end_pdf: false # Plot all individual plots listed below as publication-quality .pdf (even if False)? data: false # Plot individual plots of the data? noise_map: false # Plot individual plots of the noise-map? signal_to_noise_map: false # Plot individual plots of the signal-to-noise-map? @@ -31,33 +22,14 @@ fit_imaging: {} # Settings for plots of fits to imagi inversion: # Settings for plots of inversions (e.g. InversionPlotter). subplot_inversion: true # Plot subplot of all quantities in each inversion (e.g. reconstrucuted image, reconstruction)? subplot_mappings: true # Plot subplot of the image-to-source pixels mappings of each pixelization? - all_at_end_png: true # Plot all individual plots listed below as .png (even if False)? - all_at_end_fits: true # Plot all individual plots listed below as .fits (even if False)? - all_at_end_pdf: false # Plot all individual plots listed below as publication-quality .pdf (even if False)? data_subtracted: false # Plot individual plots of the data with the other inversion linear objects subtracted? - errors: false # Plot image of the errors of every mesh-pixel reconstructed value? + reconstruction_noise_map: false # Plot image of the noise of every mesh-pixel reconstructed value? sub_pixels_per_image_pixels: false # Plot the number of sub pixels per masked data pixels? mesh_pixels_per_image_pixels: false # Plot the number of image-plane mesh pixels per masked data pixels? image_pixels_per_mesh_pixels: false # Plot the number of image pixels in each pixel of the mesh? reconstructed_image: false # Plot image of the reconstructed data (e.g. in the image-plane)? reconstruction: false # Plot the reconstructed inversion (e.g. the pixelization's mesh in the source-plane)? regularization_weights: false # Plot the effective regularization weight of every inversion mesh pixel? -interferometer: # Settings for plots of interferometer datasets (e.g. InterferometerPlotter). - amplitudes_vs_uv_distances: false - phases_vs_uv_distances: false - uv_wavelengths: false - dirty_image: false - dirty_noise_map: false - dirty_signal_to_noise_map: false fit_interferometer: # Settings for plots of fits to interferometer datasets (e.g. FitInterferometerPlotter). subplot_fit_dirty_images: false # Plot subplot of the dirty-images of all interferometer datasets? - subplot_fit_real_space: false # Plot subplot of the real-space images of all interferometer datasets? - amplitudes_vs_uv_distances: false - phases_vs_uv_distances: false - uv_wavelengths: false - dirty_image: false - dirty_noise_map: false - dirty_signal_to_noise_map: false - dirty_residual_map: false - dirty_normalized_residual_map: false - dirty_chi_squared_map: false \ No newline at end of file + subplot_fit_real_space: false # Plot subplot of the real-space images of all interferometer datasets? \ No newline at end of file diff --git a/autoarray/dataset/abstract/dataset.py b/autoarray/dataset/abstract/dataset.py index 62bf602fb..726211e7f 100644 --- a/autoarray/dataset/abstract/dataset.py +++ b/autoarray/dataset/abstract/dataset.py @@ -4,7 +4,8 @@ import warnings from typing import Optional, Union -from autoarray.dataset.over_sampling import OverSamplingDataset +from autoconf import cached_property + from autoarray.dataset.grids import GridsDataset from autoarray import exc @@ -12,7 +13,8 @@ from autoarray.mask.mask_2d import Mask2D from autoarray.structures.abstract_structure import Structure from autoarray.structures.arrays.uniform_2d import Array2D -from autoconf import cached_property + +from autoarray.operators.over_sampling import over_sample_util logger = logging.getLogger(__name__) @@ -24,7 +26,8 @@ def __init__( data: Structure, noise_map: Structure, noise_covariance_matrix: Optional[np.ndarray] = None, - over_sampling: Optional[OverSamplingDataset] = OverSamplingDataset(), + over_sample_size_lp: Union[int, Array2D] = 4, + over_sample_size_pixelization: Union[int, Array2D] = 4, ): """ An abstract dataset, containing the image data, noise-map, PSF and associated quantities for calculations @@ -45,6 +48,32 @@ def __init__( over sampling calculations built in which approximate the 2D line integral of these calculations within a pixel. This is explained in more detail in the `GridsDataset` class. + **Over Sampling** + + If a grid is uniform and the centre of each point on the grid is the centre of a 2D pixel, evaluating + the value of a function on the grid requires a 2D line integral to compute it precisely. This can be + computationally expensive and difficult to implement. + + Over sampling is a numerical technique where the function is evaluated on a sub-grid within each grid pixel + which is higher resolution than the grid itself. This approximates more closely the value of the function + within a 2D line intergral of the values in the square pixel that the grid is centred. + + For example, in PyAutoGalaxy and PyAutoLens the light profiles and galaxies are evaluated in order to determine + how much light falls in each pixel. This uses over sampling and therefore a higher resolution grid than the + image data to ensure the calculation is accurate. + + This class controls how over sampling is performed for 2 different types of grids: + + - `lp`: A grids of (y,x) coordinates which aligns with the centre of every image pixel of the image data + and is used to evaluate light profiles for model-fititng. + + - `pixelization`: A grid of (y,x) coordinates which again align with the centre of every image pixel of + the image data. This grid is used specifically for pixelizations computed via the `inversion` module, which + can benefit from using different oversampling schemes than the normal grid. + + Different calculations typically benefit from different over sampling, which this class enables + the customization of. + Parameters ---------- data @@ -56,10 +85,13 @@ def __init__( noise_covariance_matrix A noise-map covariance matrix representing the covariance between noise in every `data` value, which can be used via a bespoke fit to account for correlated noise in the data. - over_sampling - The over sampling schemes which divide the grids into sub grids of smaller pixels within their host image - pixels when using the grid to evaluate a function (e.g. images) to better approximate the 2D line integral - This class controls over sampling for all the different grids (e.g. `grid`, `grids.pixelization). + over_sample_size_lp + The over sampling scheme size, which divides the grid into a sub grid of smaller pixels when computing + values (e.g. images) from the grid to approximate the 2D line integral of the amount of light that falls + into each pixel. + over_sample_size_pixelization + How over sampling is performed for the grid which is associated with a pixelization, which is therefore + passed into the calculations performed in the `inversion` module. """ self.data = data @@ -93,15 +125,28 @@ def __init__( self.noise_map = noise_map - self.over_sampling = over_sampling + self.over_sample_size_lp = ( + over_sample_util.over_sample_size_convert_to_array_2d_from( + over_sample_size=over_sample_size_lp, mask=self.mask + ) + ) + self.over_sample_size_pixelization = ( + over_sample_util.over_sample_size_convert_to_array_2d_from( + over_sample_size=over_sample_size_pixelization, mask=self.mask + ) + ) @property def grid(self): - return self.grids.uniform + return self.grids.lp @cached_property def grids(self): - return GridsDataset(mask=self.data.mask, over_sampling=self.over_sampling) + return GridsDataset( + mask=self.data.mask, + over_sample_size_lp=self.over_sample_size_lp, + over_sample_size_pixelization=self.over_sample_size_pixelization, + ) @property def shape_native(self): @@ -161,4 +206,16 @@ def trimmed_after_convolution_from(self, kernel_shape) -> "AbstractDataset": kernel_shape=kernel_shape ) + dataset.over_sample_size_lp = ( + dataset.over_sample_size_lp.trimmed_after_convolution_from( + kernel_shape=kernel_shape + ) + ) + + dataset.over_sample_size_pixelization = ( + dataset.over_sample_size_pixelization.trimmed_after_convolution_from( + kernel_shape=kernel_shape + ) + ) + return dataset diff --git a/autoarray/dataset/grids.py b/autoarray/dataset/grids.py index 52390ae33..c460bf820 100644 --- a/autoarray/dataset/grids.py +++ b/autoarray/dataset/grids.py @@ -1,12 +1,11 @@ from typing import Optional, Union -from autoarray.dataset.over_sampling import OverSamplingDataset from autoarray.mask.mask_2d import Mask2D +from autoarray.structures.arrays.uniform_2d import Array2D from autoarray.structures.arrays.kernel_2d import Kernel2D from autoarray.structures.grids.uniform_1d import Grid1D from autoarray.structures.grids.uniform_2d import Grid2D -from autoarray.operators.over_sampling.uniform import OverSamplingUniform from autoarray.inversion.pixelization.border_relocator import BorderRelocator from autoconf import cached_property @@ -15,7 +14,8 @@ class GridsDataset: def __init__( self, mask: Mask2D, - over_sampling: OverSamplingDataset, + over_sample_size_lp: Union[int, Array2D], + over_sample_size_pixelization: Union[int, Array2D], psf: Optional[Kernel2D] = None, ): """ @@ -28,11 +28,6 @@ def __init__( which is used for most normal calculations (e.g. evaluating the amount of light that falls in an pixel from a light profile). - - `non_uniform`: A grid of (y,x) coordinates which aligns with the centre of every image pixel of the image - data, but where their values are going to be deflected to become non-uniform such that the adaptive over - sampling scheme used for the main grid does not apply. This is used to compute over sampled light profiles of - lensed sources in PyAutoLens. - - `pixelization`: A grid of (y,x) coordinates which again align with the centre of every image pixel of the image data. This grid is used specifically for pixelizations computed via the `inversion` module, which can benefit from using different oversampling schemes than the normal grid. @@ -49,16 +44,24 @@ def __init__( Parameters ---------- - mask - over_sampling + over_sample_size_lp + The over sampling scheme size, which divides the grid into a sub grid of smaller pixels when computing + values (e.g. images) from the grid to approximate the 2D line integral of the amount of light that falls + into each pixel. + over_sample_size_pixelization + How over sampling is performed for the grid which is associated with a pixelization, which is therefore + passed into the calculations performed in the `inversion` module. psf + The Point Spread Function kernel of the image which accounts for diffraction due to the telescope optics + via 2D convolution. """ self.mask = mask - self.over_sampling = over_sampling + self.over_sample_size_lp = over_sample_size_lp + self.over_sample_size_pixelization = over_sample_size_pixelization self.psf = psf @cached_property - def uniform(self) -> Union[Grid1D, Grid2D]: + def lp(self) -> Union[Grid1D, Grid2D]: """ Returns the grid of (y,x) Cartesian coordinates at the centre of every pixel in the masked data, which is used to perform most normal calculations (e.g. evaluating the amount of light that falls in an pixel from a light @@ -70,38 +73,9 @@ def uniform(self) -> Union[Grid1D, Grid2D]: ------- The (y,x) coordinates of every pixel in the data. """ - - return Grid2D.from_mask( - mask=self.mask, - over_sampling=self.over_sampling.uniform, - ) - - @cached_property - def non_uniform(self) -> Optional[Union[Grid1D, Grid2D]]: - """ - Returns the grid of (y,x) Cartesian coordinates at the centre of every pixel in the masked data, but - with a different over sampling scheme designed for - - where - their values are going to be deflected to become non-uniform such that the adaptive over sampling scheme used - for the main grid does not apply. - - This is used to compute over sampled light profiles of lensed sources in PyAutoLens. - - - This grid is computed based on the mask, in particular its pixel-scale and sub-grid size. - - Returns - ------- - The (y,x) coordinates of every pixel in the data. - """ - - if self.over_sampling.non_uniform is None: - return None - return Grid2D.from_mask( mask=self.mask, - over_sampling=self.over_sampling.non_uniform, + over_sample_size=self.over_sample_size_lp, ) @cached_property @@ -120,15 +94,9 @@ def pixelization(self) -> Grid2D: ------- The (y,x) coordinates of every pixel in the data, used for pixelization / inversion calculations. """ - - over_sampling = self.over_sampling.pixelization - - if over_sampling is None: - over_sampling = OverSamplingUniform(sub_size=4) - return Grid2D.from_mask( mask=self.mask, - over_sampling=over_sampling, + over_sample_size=self.over_sample_size_pixelization, ) @cached_property @@ -151,36 +119,26 @@ def blurring(self) -> Optional[Grid2D]: if self.psf is None: return None - return self.uniform.blurring_grid_via_kernel_shape_from( + return self.lp.blurring_grid_via_kernel_shape_from( kernel_shape_native=self.psf.shape_native, ) - @cached_property - def over_sampler_non_uniform(self): - return self.non_uniform.over_sampling.over_sampler_from(mask=self.mask) - - @cached_property - def over_sampler_pixelization(self): - return self.pixelization.over_sampling.over_sampler_from(mask=self.mask) - @cached_property def border_relocator(self) -> BorderRelocator: return BorderRelocator( - mask=self.mask, sub_size=self.pixelization.over_sampling.sub_size + mask=self.mask, sub_size=self.over_sample_size_pixelization ) class GridsInterface: def __init__( self, - uniform=None, - non_uniform=None, + lp=None, pixelization=None, blurring=None, border_relocator=None, ): - self.uniform = uniform - self.non_uniform = non_uniform + self.lp = lp self.pixelization = pixelization self.blurring = blurring self.border_relocator = border_relocator diff --git a/autoarray/dataset/imaging/dataset.py b/autoarray/dataset/imaging/dataset.py index 85bc3faac..78449ee47 100644 --- a/autoarray/dataset/imaging/dataset.py +++ b/autoarray/dataset/imaging/dataset.py @@ -8,7 +8,6 @@ from autoarray.dataset.abstract.dataset import AbstractDataset from autoarray.dataset.grids import GridsDataset from autoarray.dataset.imaging.w_tilde import WTildeImaging -from autoarray.dataset.over_sampling import OverSamplingDataset from autoarray.structures.arrays.uniform_2d import Array2D from autoarray.operators.convolver import Convolver from autoarray.structures.arrays.kernel_2d import Kernel2D @@ -16,6 +15,7 @@ from autoarray import type as ty from autoarray import exc +from autoarray.operators.over_sampling import over_sample_util from autoarray.inversion.inversion.imaging import inversion_imaging_util logger = logging.getLogger(__name__) @@ -28,7 +28,8 @@ def __init__( noise_map: Optional[Array2D] = None, psf: Optional[Kernel2D] = None, noise_covariance_matrix: Optional[np.ndarray] = None, - over_sampling: Optional[OverSamplingDataset] = OverSamplingDataset(), + over_sample_size_lp: Union[int, Array2D] = 4, + over_sample_size_pixelization: Union[int, Array2D] = 4, pad_for_convolver: bool = False, use_normalized_psf: Optional[bool] = True, check_noise_map: bool = True, @@ -69,10 +70,13 @@ def __init__( noise_covariance_matrix A noise-map covariance matrix representing the covariance between noise in every `data` value, which can be used via a bespoke fit to account for correlated noise in the data. - over_sampling - The over sampling schemes which divide the grids into sub grids of smaller pixels within their host image - pixels when using the grid to evaluate a function (e.g. images) to better approximate the 2D line integral - This class controls over sampling for all the different grids (e.g. `grid`, `grids.pixelization). + over_sample_size_lp + The over sampling scheme size, which divides the grid into a sub grid of smaller pixels when computing + values (e.g. images) from the grid to approximate the 2D line integral of the amount of light that falls + into each pixel. + over_sample_size_pixelization + How over sampling is performed for the grid which is associated with a pixelization, which is therefore + passed into the calculations performed in the `inversion` module. pad_for_convolver The PSF convolution may extend beyond the edges of the image mask, which can lead to edge effects in the convolved image. If `True`, the image and noise-map are padded to ensure the PSF convolution does not @@ -94,6 +98,28 @@ def __init__( kernel_shape_native=psf.shape_native ) except exc.MaskException: + over_sample_size_lp = ( + over_sample_util.over_sample_size_convert_to_array_2d_from( + over_sample_size=over_sample_size_lp, mask=data.mask + ) + ) + over_sample_size_lp = ( + over_sample_size_lp.padded_before_convolution_from( + kernel_shape=psf.shape_native, mask_pad_value=1 + ) + ) + print(over_sample_size_lp.shape_native) + over_sample_size_pixelization = ( + over_sample_util.over_sample_size_convert_to_array_2d_from( + over_sample_size=over_sample_size_pixelization, mask=data.mask + ) + ) + over_sample_size_pixelization = ( + over_sample_size_pixelization.padded_before_convolution_from( + kernel_shape=psf.shape_native, mask_pad_value=1 + ) + ) + data = data.padded_before_convolution_from( kernel_shape=psf.shape_native, mask_pad_value=1 ) @@ -114,7 +140,8 @@ def __init__( data=data, noise_map=noise_map, noise_covariance_matrix=noise_covariance_matrix, - over_sampling=over_sampling, + over_sample_size_lp=over_sample_size_lp, + over_sample_size_pixelization=over_sample_size_pixelization, ) self.use_normalized_psf = use_normalized_psf @@ -143,7 +170,10 @@ def __init__( @cached_property def grids(self): return GridsDataset( - mask=self.data.mask, over_sampling=self.over_sampling, psf=self.psf + 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, ) @cached_property @@ -214,7 +244,8 @@ def from_fits( psf_hdu: int = 0, noise_covariance_matrix: Optional[np.ndarray] = None, check_noise_map: bool = True, - over_sampling: Optional[OverSamplingDataset] = OverSamplingDataset(), + over_sample_size_lp: Union[int, Array2D] = 4, + over_sample_size_pixelization: Union[int, Array2D] = 4, ) -> "Imaging": """ Load an imaging dataset from multiple .fits file. @@ -253,10 +284,13 @@ def from_fits( A noise-map covariance matrix representing the covariance between noise in every `data` value. check_noise_map If True, the noise-map is checked to ensure all values are above zero. - over_sampling - The over sampling schemes which divide the grids into sub grids of smaller pixels within their host image - pixels when using the grid to evaluate a function (e.g. images) to better approximate the 2D line integral - This class controls over sampling for all the different grids (e.g. `grid`, `grids.pixelization). + over_sample_size_lp + The over sampling scheme size, which divides the grid into a sub grid of smaller pixels when computing + values (e.g. images) from the grid to approximate the 2D line integral of the amount of light that falls + into each pixel. + over_sample_size_pixelization + How over sampling is performed for the grid which is associated with a pixelization, which is therefore + passed into the calculations performed in the `inversion` module. """ data = Array2D.from_fits( @@ -284,7 +318,8 @@ def from_fits( psf=psf, noise_covariance_matrix=noise_covariance_matrix, check_noise_map=check_noise_map, - over_sampling=over_sampling, + over_sample_size_lp=over_sample_size_lp, + over_sample_size_pixelization=over_sample_size_pixelization, ) def apply_mask(self, mask: Mask2D) -> "Imaging": @@ -323,12 +358,18 @@ def apply_mask(self, mask: Mask2D) -> "Imaging": else: noise_covariance_matrix = None + over_sample_size_lp = Array2D(values=self.over_sample_size_lp.native, mask=mask) + over_sample_size_pixelization = Array2D( + values=self.over_sample_size_pixelization.native, mask=mask + ) + dataset = Imaging( data=data, noise_map=noise_map, psf=self.psf, noise_covariance_matrix=noise_covariance_matrix, - over_sampling=self.over_sampling, + over_sample_size_lp=over_sample_size_lp, + over_sample_size_pixelization=over_sample_size_pixelization, pad_for_convolver=True, ) @@ -340,11 +381,27 @@ def apply_mask(self, mask: Mask2D) -> "Imaging": return dataset - def apply_noise_scaling(self, mask: Mask2D, noise_value: float = 1e8) -> "Imaging": + def apply_noise_scaling( + self, + mask: Mask2D, + noise_value: float = 1e8, + signal_to_noise_value: Optional[float] = None, + should_zero_data: bool = True, + ) -> "Imaging": """ - Apply a mask to the imaging dataset using noise scaling, whereby the mask increases noise-map values to be - extremely large such that they are never included in the likelihood calculation, but it does - not remove the image data values, which are set to zero. + Apply a mask to the imaging dataset using noise scaling, whereby the maskmay zero the data and increase + noise-map values to change how they enter the likelihood calculation. + + Given this data region is masked, it is likely thr data itself should not be included and therefore + the masked data values are set to zero. This can be disabled by setting `should_zero_data=False`. + + Two forms of scaling are supported depending on whether the `signal_to_noise_value` is input: + + - `noise_value`: The noise-map values in the masked region are set to this value, typically a very large value, + such that they are never included in the likelihood calculation. + + - `signal_to_noise_value`: The noise-map values in the masked region are set to values such that they give + this signal-to-noise ratio. This overwrites the `noise_value` parameter. For certain modeling tasks, the mask defines regions of the data that are used to calculate the likelihood. For example, all data points in a mask may be used to create a pixel-grid, which is used in the likelihood. @@ -358,32 +415,63 @@ def apply_noise_scaling(self, mask: Mask2D, noise_value: float = 1e8) -> "Imagin ---------- mask The 2D mask that is applied to the image and noise-map, to scale the noise-map values to large values. + noise_value + The value that the noise-map values are set to in the masked region where noise scaling is applied. + signal_to_noise_value + The noise-map values are instead set to values such that they give this signal-to-noise_maps ratio. + This overwrites the noise_value parameter. + should_zero_data + If True, the data values in the masked region are set to zero. """ - data = np.where(np.invert(mask), 0.0, self.data.native) - data = Array2D.no_mask( + + if signal_to_noise_value is None: + noise_map = self.noise_map.native + noise_map[mask == False] = noise_value + else: + noise_map = np.where( + mask == False, + np.median(self.data.native[mask.derive_mask.edge == False]) + / signal_to_noise_value, + self.noise_map.native, + ) + + if should_zero_data: + data = np.where(np.invert(mask), 0.0, self.data.native) + else: + data = self.data.native + + data_unmasked = Array2D.no_mask( values=data, shape_native=self.data.shape_native, pixel_scales=self.data.pixel_scales, ) - noise_map = self.noise_map.native - noise_map[mask == False] = noise_value - noise_map = Array2D.no_mask( + noise_map_unmasked = Array2D.no_mask( values=noise_map, - shape_native=self.data.shape_native, - pixel_scales=self.data.pixel_scales, + shape_native=self.noise_map.shape_native, + pixel_scales=self.noise_map.pixel_scales, ) + data = Array2D(values=data, mask=self.data.mask) + + noise_map = Array2D(values=noise_map, mask=self.data.mask) + dataset = Imaging( data=data, noise_map=noise_map, psf=self.psf, noise_covariance_matrix=self.noise_covariance_matrix, - over_sampling=self.over_sampling, + over_sample_size_lp=self.over_sample_size_lp, + over_sample_size_pixelization=self.over_sample_size_pixelization, pad_for_convolver=False, - check_noise_map=False + check_noise_map=False, ) + if self.unmasked is not None: + dataset.unmasked = self.unmasked + dataset.unmasked.data = data_unmasked + dataset.unmasked.noise_map = noise_map_unmasked + logger.info( f"IMAGING - Data noise scaling applied, a total of {mask.pixels_in_mask} pixels were scaled to large noise values." ) @@ -392,7 +480,8 @@ def apply_noise_scaling(self, mask: Mask2D, noise_value: float = 1e8) -> "Imagin def apply_over_sampling( self, - over_sampling: Optional[OverSamplingDataset] = OverSamplingDataset(), + over_sample_size_lp: Union[int, Array2D] = None, + over_sample_size_pixelization: Union[int, Array2D] = None, ) -> "AbstractDataset": """ Apply new over sampling objects to the grid and grid pixelization of the dataset. @@ -404,37 +493,26 @@ def apply_over_sampling( This function resets the cached properties so that the new over sampling is used in the grid and grid pixelization. - The `default_galaxy_mode` parameter is used to set up default over sampling for galaxy light profiles in - the project PyAutoGalaxy. This sets up the over sampling such that there is high over sampling in the centre - of the mask, where the galaxy is located, and lower over sampling in the outer regions of the mask. It - does this based on the pixel scale, which gives a good estimate of how large the central region - requiring over sampling is. - Parameters ---------- - over_sampling - The over sampling schemes which divide the grids into sub grids of smaller pixels within their host image - pixels when using the grid to evaluate a function (e.g. images) to better approximate the 2D line integral - This class controls over sampling for all the different grids (e.g. `grid`, `grids.pixelization). + over_sample_size_lp + The over sampling scheme size, which divides the grid into a sub grid of smaller pixels when computing + values (e.g. images) from the grid to approximate the 2D line integral of the amount of light that falls + into each pixel. + over_sample_size_pixelization + How over sampling is performed for the grid which is associated with a pixelization, which is therefore + passed into the calculations performed in the `inversion` module. """ - uniform = over_sampling.uniform or self.over_sampling.uniform - non_uniform = over_sampling.non_uniform or self.over_sampling.non_uniform - pixelization = over_sampling.pixelization or self.over_sampling.pixelization - - over_sampling = OverSamplingDataset( - uniform=uniform, - non_uniform=non_uniform, - pixelization=pixelization, - ) - return Imaging( data=self.data, noise_map=self.noise_map, psf=self.psf, - over_sampling=over_sampling, + over_sample_size_lp=over_sample_size_lp or self.over_sample_size_lp, + over_sample_size_pixelization=over_sample_size_pixelization + or self.over_sample_size_pixelization, pad_for_convolver=False, - check_noise_map=False + check_noise_map=False, ) def output_to_fits( diff --git a/autoarray/dataset/imaging/simulator.py b/autoarray/dataset/imaging/simulator.py index 35babbc07..449a3d4b6 100644 --- a/autoarray/dataset/imaging/simulator.py +++ b/autoarray/dataset/imaging/simulator.py @@ -1,5 +1,6 @@ import logging import numpy as np +from typing import Optional, Union from autoarray.dataset.imaging.dataset import Imaging from autoarray.structures.arrays.uniform_2d import Array2D @@ -20,28 +21,45 @@ def __init__( subtract_background_sky: bool = True, psf: Kernel2D = None, normalize_psf: bool = True, - add_poisson_noise: bool = True, + add_poisson_noise_to_data: bool = True, + include_poisson_noise_in_noise_map: bool = True, noise_if_add_noise_false: float = 0.1, noise_seed: int = -1, ): """ - Simulations observations of imaging data, including simulation of the image, noise-map, PSF, etc. as - an `Imaging` object. + Simulates observations of `Imaging` data, including simulating the image, noise, blurring due to the telescope + optics via the Point Spread Function (PSF) and the background sky. The simulation of an `Imaging` dataset uses the following steps: - 1) Receive as input the raw image which is simulated via the steps below. - 2) Convolve the image with the Point Spread Function of the simulated dataset. + 1) Receive as input a raw image of what the data looks like before any simulaiton process is applied. + 2) Include dirrection due to the telescope optics by convolve the image with an input Point Spread + Function (PSF). 3) Use input values of the background sky level in every pixel of the image to add the background sky to the PSF convolved image. 4) Add Poisson noise to the image, which represents noise due to whether photons hits the CCD and are converted - to photo-electrons which are succcessfully detected by the CCD and converted to counts. + to photo-electrons which are successfully detected by the CCD and converted to counts. 5) Subtract the background sky from the image, so that the returned simulated dataset is background sky subtracted. The inputs of the `SimulatorImaging` object can toggle these steps on and off, for example if `psf=None` the PSF convolution step is omitted. + The inputs `add_poisson_noise_to_data` and `include_poisson_noise_in_noise_map` provide flexibility in + handling Poisson noise during simulations: + + - `add_poisson_noise_to_data` determines whether Poisson noise is added to the simulated data. + - `include_poisson_noise_in_noise_map` decides if Poisson noise is included in the `noise_map` of the dataset. + + These options allow you to create a dataset without adding Poisson noise to the data while still using + a `noise_map` that represents the noise levels expected if Poisson noise had been included. + + This approach is particularly useful for diagnostic purposes. By simulating noise-free data but retaining + realistic noise levels in the `noise_map`, you can test whether the model-fitting process accurately recovers + the input parameters. Unlike noisy data, where random Poisson noise introduces offsets in each simulation, + this method ensures the posterior peaks precisely around the input values. This is a powerful way to identify + potential systematic biases in a model-fitting analysis. + Parameters ---------- exposure_time @@ -54,8 +72,11 @@ def __init__( An array describing the PSF kernel of the image. normalize_psf If `True`, the PSF kernel is normalized so all values sum to 1.0. - add_poisson_noise + add_poisson_noise_to_data Whether Poisson noise corresponding to photon count statistics on the imaging observation is added. + include_poisson_noise_in_noise_map + Whether Poisson noise is included in the noise-map of the dataset. If `False` the noise-map is filled with + the value `noise_if_add_noise_false`. noise_if_add_noise_false If noise is not added to the simulated dataset a `noise_map` must still be returned. This value gives the value of noise assigned to every pixel in the noise-map. @@ -73,11 +94,14 @@ def __init__( self.exposure_time = exposure_time self.background_sky_level = background_sky_level self.subtract_background_sky = subtract_background_sky - self.add_poisson_noise = add_poisson_noise + self.add_poisson_noise_to_data = add_poisson_noise_to_data + self.include_poisson_noise_in_noise_map = include_poisson_noise_in_noise_map self.noise_if_add_noise_false = noise_if_add_noise_false self.noise_seed = noise_seed - def via_image_from(self, image: Array2D) -> Imaging: + def via_image_from( + self, image: Array2D, over_sample_size: Optional[Union[int, np.ndarray]] = None + ) -> Imaging: """ Simulate an `Imaging` dataset from an input image. @@ -106,15 +130,18 @@ def via_image_from(self, image: Array2D) -> Imaging: image = image + background_sky_map - if self.add_poisson_noise is True: - image = preprocess.data_eps_with_poisson_noise_added( - data_eps=image, - exposure_time_map=exposure_time_map, - seed=self.noise_seed, - ) + image_with_poisson_noise = preprocess.data_eps_with_poisson_noise_added( + data_eps=image, + exposure_time_map=exposure_time_map, + seed=self.noise_seed, + ) + if self.add_poisson_noise_to_data: + image = image_with_poisson_noise + + if self.include_poisson_noise_in_noise_map: noise_map = preprocess.noise_map_via_data_eps_and_exposure_time_map_from( - data_eps=image, exposure_time_map=exposure_time_map + data_eps=image_with_poisson_noise, exposure_time_map=exposure_time_map ) else: @@ -139,4 +166,13 @@ def via_image_from(self, image: Array2D) -> Imaging: image = Array2D(values=image, mask=mask) - return Imaging(data=image, psf=self.psf, noise_map=noise_map) + dataset = Imaging( + data=image, psf=self.psf, noise_map=noise_map, check_noise_map=False + ) + + if over_sample_size is not None: + dataset = dataset.apply_over_sampling( + over_sample_size_lp=over_sample_size.native + ) + + return dataset diff --git a/autoarray/dataset/interferometer/dataset.py b/autoarray/dataset/interferometer/dataset.py index 58023719f..f69753dd0 100644 --- a/autoarray/dataset/interferometer/dataset.py +++ b/autoarray/dataset/interferometer/dataset.py @@ -7,8 +7,8 @@ from autoarray.dataset.abstract.dataset import AbstractDataset from autoarray.dataset.interferometer.w_tilde import WTildeInterferometer from autoarray.dataset.grids import GridsDataset -from autoarray.dataset.over_sampling import OverSamplingDataset from autoarray.operators.transformer import TransformerNUFFT + from autoarray.structures.visibilities import Visibilities from autoarray.structures.visibilities import VisibilitiesNoiseMap @@ -25,7 +25,6 @@ def __init__( uv_wavelengths: np.ndarray, real_space_mask, transformer_class=TransformerNUFFT, - over_sampling: Optional[OverSamplingDataset] = OverSamplingDataset(), ): """ An interferometer dataset, containing the visibilities data, noise-map, real-space msk, Fourier transformer and @@ -68,10 +67,6 @@ def __init__( noise_covariance_matrix A noise-map covariance matrix representing the covariance between noise in every `data` value, which can be used via a bespoke fit to account for correlated noise in the data. - over_sampling - The over sampling schemes which divide the grids into sub grids of smaller pixels within their host image - pixels when using the grid to evaluate a function (e.g. images) to better approximate the 2D line integral - This class controls over sampling for all the different grids (e.g. `grid`, `grids.pixelization). transformer_class The class of the Fourier Transform which maps images from real space to Fourier space visibilities and the uv-plane. @@ -81,7 +76,8 @@ def __init__( super().__init__( data=data, noise_map=noise_map, - over_sampling=over_sampling, + over_sample_size_lp=1, + over_sample_size_pixelization=1, ) self.uv_wavelengths = uv_wavelengths @@ -92,53 +88,10 @@ def __init__( @cached_property def grids(self): - return GridsDataset(mask=self.real_space_mask, over_sampling=self.over_sampling) - - def apply_over_sampling( - self, - over_sampling: Optional[OverSamplingDataset] = OverSamplingDataset(), - ) -> "Interferometer": - """ - Apply new over sampling objects to the grid and grid pixelization of the dataset. - - This method is used to change the over sampling of the grid and grid pixelization, for example when the - user wishes to perform over sampling with a higher sub grid size or with an iterative over sampling strategy. - - The `grid` and grids.pixelization` are cached properties which after use are stored in memory for efficiency. - This function resets the cached properties so that the new over sampling is used in the grid and grid - pixelization. - - The `default_galaxy_mode` parameter is used to set up default over sampling for galaxy light profiles in - the project PyAutoGalaxy. This sets up the over sampling such that there is high over sampling in the centre - of the mask, where the galaxy is located, and lower over sampling in the outer regions of the mask. It - does this based on the pixel scale, which gives a good estimate of how large the central region - requiring over sampling is. - - Parameters - ---------- - over_sampling - The over sampling schemes which divide the grids into sub grids of smaller pixels within their host image - pixels when using the grid to evaluate a function (e.g. images) to better approximate the 2D line integral - This class controls over sampling for all the different grids (e.g. `grid`, `grids.pixelization). - """ - - uniform = over_sampling.uniform or self.over_sampling.uniform - non_uniform = over_sampling.non_uniform or self.over_sampling.non_uniform - pixelization = over_sampling.pixelization or self.over_sampling.pixelization - - over_sampling = OverSamplingDataset( - uniform=uniform, - non_uniform=non_uniform, - pixelization=pixelization, - ) - - return Interferometer( - data=self.data, - noise_map=self.noise_map, - uv_wavelengths=self.uv_wavelengths, - real_space_mask=self.real_space_mask, - transformer_class=self.transformer.__class__, - over_sampling=over_sampling, + return GridsDataset( + mask=self.real_space_mask, + over_sample_size_lp=self.over_sample_size_lp, + over_sample_size_pixelization=self.over_sample_size_pixelization, ) @classmethod @@ -152,7 +105,6 @@ def from_fits( noise_map_hdu=0, uv_wavelengths_hdu=0, transformer_class=TransformerNUFFT, - over_sampling: Optional[OverSamplingDataset] = OverSamplingDataset(), ): """ Factory for loading the interferometer data_type from .fits files, as well as computing properties like the @@ -178,7 +130,6 @@ def from_fits( noise_map=noise_map, uv_wavelengths=uv_wavelengths, transformer_class=transformer_class, - over_sampling=over_sampling, ) @cached_property diff --git a/autoarray/dataset/interferometer/simulator.py b/autoarray/dataset/interferometer/simulator.py index 8a6a38719..1f8ea2210 100644 --- a/autoarray/dataset/interferometer/simulator.py +++ b/autoarray/dataset/interferometer/simulator.py @@ -56,7 +56,7 @@ def via_image_from(self, image): An array representing the effective exposure time of each pixel. psf: PSF An array describing the PSF the simulated image is blurred with. - add_poisson_noise: Bool + add_poisson_noise_to_data: Bool If `True` poisson noise_maps is simulated and added to the image, based on the total counts in each image pixel noise_seed: int diff --git a/autoarray/dataset/mock/mock_dataset.py b/autoarray/dataset/mock/mock_dataset.py index 29ece502d..c425c6a27 100644 --- a/autoarray/dataset/mock/mock_dataset.py +++ b/autoarray/dataset/mock/mock_dataset.py @@ -6,8 +6,7 @@ class MockDataset: def __init__(self, grids=None, psf=None, mask=None): if grids is None: self.grids = GridsInterface( - uniform=None, - non_uniform=None, + lp=None, pixelization=Grid2D.no_mask(values=[[[1.0, 1.0]]], pixel_scales=1.0), blurring=None, border_relocator=None, diff --git a/autoarray/dataset/over_sampling.py b/autoarray/dataset/over_sampling.py deleted file mode 100644 index e707e3251..000000000 --- a/autoarray/dataset/over_sampling.py +++ /dev/null @@ -1,62 +0,0 @@ -import logging -from typing import Optional - -from autoarray.operators.over_sampling.abstract import AbstractOverSampling - -logger = logging.getLogger(__name__) - - -class OverSamplingDataset: - def __init__( - self, - uniform: Optional[AbstractOverSampling] = None, - non_uniform: Optional[AbstractOverSampling] = None, - pixelization: Optional[AbstractOverSampling] = None, - ): - """ - Customizes how over sampling calculations are performed using the grids of the data. - - If a grid is uniform and the centre of each point on the grid is the centre of a 2D pixel, evaluating - the value of a function on the grid requires a 2D line integral to compute it precisely. This can be - computationally expensive and difficult to implement. - - Over sampling is a numerical technique where the function is evaluated on a sub-grid within each grid pixel - which is higher resolution than the grid itself. This approximates more closely the value of the function - within a 2D line intergral of the values in the square pixel that the grid is centred. - - For example, in PyAutoGalaxy and PyAutoLens the light profiles and galaxies are evaluated in order to determine - how much light falls in each pixel. This uses over sampling and therefore a higher resolution grid than the - image data to ensure the calculation is accurate. - - This class controls how over sampling is performed for 3 different types of grids: - - `grid`: A grids of (y,x) coordinates which aligns with the centre of every image pixel of the image data. - - - `grids.pixelization`: A grid of (y,x) coordinates which again align with the centre of every image pixel of - the image data. This grid is used specifically for pixelizations computed via the `inversion` module, which - can benefit from using different oversampling schemes than the normal grid. - - - `grid_non_uniform`: A grid of (y,x) coordinates which are mapped from the image pixel centres but have had - their values deflected to become non-uniform. This is used to compute over sampled light profiles of lensed - sources in PyAutoLens. - - Different calculations typically benefit from different over sampling, which this class enables - the customization of. - - Parameters - ---------- - uniform - The over sampling scheme, which divides the grid into a sub grid of smaller pixels when computing values - (e.g. images) from the grid so as to approximate the 2D line integral of the amount of light that falls - into each pixel. - pixelization - How over sampling is performed for the grid which is associated with a pixelization, which is therefore - passed into the calculations performed in the `inversion` module. - non_uniform - The over sampling scheme when the grid input into a function is not a uniform grid. This is used - by **PyAutoLens** when the grid has been deflected and ray-traced and therefore some of the default - over sampling schemes are not appropriate. - """ - self.uniform = uniform - self.pixelization = pixelization - self.non_uniform = non_uniform diff --git a/autoarray/dataset/plot/imaging_plotters.py b/autoarray/dataset/plot/imaging_plotters.py index a173675ba..8fed6f601 100644 --- a/autoarray/dataset/plot/imaging_plotters.py +++ b/autoarray/dataset/plot/imaging_plotters.py @@ -61,9 +61,8 @@ def figures_2d( noise_map: bool = False, psf: bool = False, signal_to_noise_map: bool = False, - over_sampling: bool = False, - over_sampling_non_uniform: bool = False, - over_sampling_pixelization: bool = False, + over_sample_size_lp: bool = False, + over_sample_size_pixelization: bool = False, title_str: Optional[str] = None, ): """ @@ -82,10 +81,10 @@ def figures_2d( Whether to make a 2D plot (via `imshow`) of the psf. signal_to_noise_map Whether to make a 2D plot (via `imshow`) of the signal-to-noise map. - over_sampling + over_sample_size_lp Whether to make a 2D plot (via `imshow`) of the Over Sampling for input light profiles. If adaptive sub size is used, the sub size grid for a centre of (0.0, 0.0) is used. - over_sampling_pixelization + over_sample_size_pixelization Whether to make a 2D plot (via `imshow`) of the Over Sampling for pixelizations. """ @@ -125,56 +124,24 @@ def figures_2d( ), ) - if over_sampling: - if self.dataset.over_sampling.uniform is None: - from autoarray.operators.over_sampling.uniform import ( - OverSamplingUniform, - ) - - over_sampling = OverSamplingUniform.from_adaptive_scheme( - grid=self.dataset.grids.uniform, - name="PlotExample", - centre=(0.0, 0.0), - ) - title = title_str or f"Over Sampling (Adaptive)" - - else: - over_sampling = self.dataset.over_sampling.uniform - title = title_str or f"Over Sampling" - - over_sampler = over_sampling.over_sampler_from( - mask=self.dataset.mask, - ) - + if over_sample_size_lp: self.mat_plot_2d.plot_array( - array=over_sampler.sub_size, + array=self.dataset.grids.over_sample_size_lp, visuals_2d=self.get_visuals_2d(), auto_labels=AutoLabels( - title=title, - filename="over_sampling", + title=title_str or f"Over Sample Size (Light Profiles)", + filename="over_sample_size_lp", cb_unit="", ), ) - if over_sampling_non_uniform: - if self.dataset.over_sampling.non_uniform is not None: - self.mat_plot_2d.plot_array( - array=self.dataset.grids.over_sampler_non_uniform.sub_size, - visuals_2d=self.get_visuals_2d(), - auto_labels=AutoLabels( - title=title_str or f"Over Sampling (Non Uniform)", - filename="over_sampling_non_uniform", - cb_unit="", - ), - ) - - if over_sampling_pixelization: + if over_sample_size_pixelization: self.mat_plot_2d.plot_array( - array=self.dataset.grids.over_sampler_pixelization.sub_size, + array=self.dataset.grids.over_sample_size_pixelization, visuals_2d=self.get_visuals_2d(), auto_labels=AutoLabels( - title=title_str or f"Over Sampling (Pixelization)", - filename="over_sampling_pixelization", + title=title_str or f"Over Sample Size (Pixelization)", + filename="over_sample_size_pixelization", cb_unit="", ), ) @@ -251,9 +218,8 @@ def subplot_dataset(self): self.figures_2d(signal_to_noise_map=True) - self.figures_2d(over_sampling=True) - self.figures_2d(over_sampling_non_uniform=True) - self.figures_2d(over_sampling_pixelization=True) + self.figures_2d(over_sample_size_lp=True) + self.figures_2d(over_sample_size_pixelization=True) self.mat_plot_2d.output.subplot_to_figure(auto_filename="subplot_dataset") self.close_subplot_figure() diff --git a/autoarray/fit/fit_dataset.py b/autoarray/fit/fit_dataset.py index 1ace1795f..c2f498949 100644 --- a/autoarray/fit/fit_dataset.py +++ b/autoarray/fit/fit_dataset.py @@ -161,8 +161,7 @@ def grids(self) -> GridsInterface: if offset[0] == 0.0 and offset[1] == 0.0: return GridsInterface( - uniform=self.dataset.grids.uniform, - non_uniform=self.dataset.grids.non_uniform, + lp=self.dataset.grids.lp, pixelization=self.dataset.grids.pixelization, blurring=self.dataset.grids.blurring, border_relocator=self.dataset.grids.border_relocator, @@ -174,11 +173,8 @@ def subtracted_from(grid, offset): return grid.subtracted_from(offset=offset) - uniform = subtracted_from( - grid=self.dataset.grids.uniform, offset=self.dataset_model.grid_offset - ) - non_uniform = subtracted_from( - grid=self.dataset.grids.non_uniform, offset=self.dataset_model.grid_offset + lp = subtracted_from( + grid=self.dataset.grids.lp, offset=self.dataset_model.grid_offset ) pixelization = subtracted_from( grid=self.dataset.grids.pixelization, offset=self.dataset_model.grid_offset @@ -188,8 +184,7 @@ def subtracted_from(grid, offset): ) return GridsInterface( - uniform=uniform, - non_uniform=non_uniform, + lp=lp, pixelization=pixelization, blurring=blurring, border_relocator=self.dataset.grids.border_relocator, diff --git a/autoarray/fit/mock/mock_fit_imaging.py b/autoarray/fit/mock/mock_fit_imaging.py index 81df04a03..181d3d56e 100644 --- a/autoarray/fit/mock/mock_fit_imaging.py +++ b/autoarray/fit/mock/mock_fit_imaging.py @@ -8,7 +8,7 @@ class MockFitImaging(FitImaging): def __init__( self, - dataset=MockDataset(), + dataset: Optional[MockDataset] = None, dataset_model: Optional[DatasetModel] = None, use_mask_in_fit: bool = False, noise_map=None, @@ -18,7 +18,7 @@ def __init__( run_time_dict: Optional[Dict] = None, ): super().__init__( - dataset=dataset, + dataset=dataset or MockDataset(), dataset_model=dataset_model, use_mask_in_fit=use_mask_in_fit, run_time_dict=run_time_dict, diff --git a/autoarray/fit/mock/mock_fit_interferometer.py b/autoarray/fit/mock/mock_fit_interferometer.py index 40378b531..8d00dad9c 100644 --- a/autoarray/fit/mock/mock_fit_interferometer.py +++ b/autoarray/fit/mock/mock_fit_interferometer.py @@ -1,3 +1,5 @@ +from typing import Optional + from autoarray.dataset.mock.mock_dataset import MockDataset from autoarray.fit.fit_interferometer import FitInterferometer @@ -5,13 +7,15 @@ class MockFitInterferometer(FitInterferometer): def __init__( self, - dataset=MockDataset(), + dataset: Optional[MockDataset] = None, use_mask_in_fit: bool = False, model_data=None, inversion=None, noise_map=None, ): - super().__init__(dataset=dataset, use_mask_in_fit=use_mask_in_fit) + super().__init__( + dataset=dataset or MockDataset(), use_mask_in_fit=use_mask_in_fit + ) self._model_data = model_data self._inversion = inversion diff --git a/autoarray/fixtures.py b/autoarray/fixtures.py index 15cec61d4..acfd277df 100644 --- a/autoarray/fixtures.py +++ b/autoarray/fixtures.py @@ -89,6 +89,14 @@ def make_grid_2d_7x7(): return aa.Grid2D.from_mask(mask=make_mask_2d_7x7()) +def make_grid_2d_sub_1_7x7(): + return aa.Grid2D.from_mask(mask=make_mask_2d_7x7(), over_sample_size=1) + + +def make_grid_2d_sub_2_7x7(): + return aa.Grid2D.from_mask(mask=make_mask_2d_7x7(), over_sample_size=2) + + def make_grid_2d_7x7_simple(): grid_2d_7x7 = make_grid_2d_7x7() grid_2d_7x7[0] = np.array([1.0, 1.0]) @@ -152,9 +160,7 @@ def make_imaging_7x7(): data=make_image_7x7(), psf=make_psf_3x3(), noise_map=make_noise_map_7x7(), - over_sampling=aa.OverSamplingDataset( - uniform=aa.OverSamplingUniform(sub_size=1) - ), + over_sample_size_lp=1, ) @@ -163,9 +169,7 @@ def make_imaging_7x7_sub_2(): data=make_image_7x7(), psf=make_psf_3x3(), noise_map=make_noise_map_7x7(), - over_sampling=aa.OverSamplingDataset( - uniform=aa.OverSamplingUniform(sub_size=2) - ), + over_sample_size_lp=2, ) @@ -174,9 +178,7 @@ def make_imaging_covariance_7x7(): data=make_image_7x7(), psf=make_psf_3x3(), noise_covariance_matrix=make_noise_covariance_matrix_7x7(), - over_sampling=aa.OverSamplingDataset( - uniform=aa.OverSamplingUniform(sub_size=1) - ), + over_sample_size_lp=1, ) @@ -185,9 +187,7 @@ def make_imaging_7x7_no_blur(): data=make_image_7x7(), psf=make_psf_3x3_no_blur(), noise_map=make_noise_map_7x7(), - over_sampling=aa.OverSamplingDataset( - uniform=aa.OverSamplingUniform(sub_size=1) - ), + over_sample_size_lp=1, ) @@ -196,9 +196,7 @@ def make_imaging_7x7_no_blur_sub_2(): data=make_image_7x7(), psf=make_psf_3x3_no_blur(), noise_map=make_noise_map_7x7(), - over_sampling=aa.OverSamplingDataset( - uniform=aa.OverSamplingUniform(sub_size=2) - ), + over_sample_size_lp=2, ) @@ -235,9 +233,6 @@ def make_interferometer_7(): uv_wavelengths=make_uv_wavelengths_7x2(), real_space_mask=make_mask_2d_7x7(), transformer_class=aa.TransformerDFT, - over_sampling=aa.OverSamplingDataset( - pixelization=aa.OverSamplingUniform(sub_size=1) - ), ) @@ -248,9 +243,6 @@ def make_interferometer_7_no_fft(): uv_wavelengths=make_uv_wavelengths_7x2_no_fft(), real_space_mask=make_mask_2d_7x7(), transformer_class=aa.TransformerDFT, - over_sampling=aa.OverSamplingDataset( - pixelization=aa.OverSamplingUniform(sub_size=1) - ), ) @@ -261,9 +253,6 @@ def make_interferometer_7_grid(): uv_wavelengths=make_uv_wavelengths_7x2(), real_space_mask=make_mask_2d_7x7(), transformer_class=aa.TransformerDFT, - over_sampling=aa.OverSamplingDataset( - pixelization=aa.OverSamplingUniform(sub_size=1) - ), ) @@ -370,7 +359,7 @@ def make_regularization_matern_kernel(): def make_rectangular_mesh_grid_3x3(): return aa.Mesh2DRectangular.overlay_grid( - grid=make_over_sampler_2d_7x7().over_sampled_grid, shape_native=(3, 3) + grid=make_grid_2d_sub_2_7x7().over_sampled, shape_native=(3, 3) ) @@ -417,7 +406,7 @@ def make_voronoi_mesh_grid_9(): def make_over_sampler_2d_7x7(): - return aa.OverSamplerUniform(mask=make_mask_2d_7x7(), sub_size=2) + return aa.OverSampler(mask=make_mask_2d_7x7(), sub_size=2) def make_border_relocator_2d_7x7(): @@ -429,7 +418,7 @@ def make_border_relocator_2d_7x7(): def make_rectangular_mapper_7x7_3x3(): mapper_grids = aa.MapperGrids( mask=make_mask_2d_7x7(), - source_plane_data_grid=make_over_sampler_2d_7x7().over_sampled_grid, + source_plane_data_grid=make_grid_2d_sub_2_7x7(), source_plane_mesh_grid=make_rectangular_mesh_grid_3x3(), image_plane_mesh_grid=None, adapt_data=aa.Array2D.ones(shape_native=(3, 3), pixel_scales=0.1), @@ -437,7 +426,6 @@ def make_rectangular_mapper_7x7_3x3(): return aa.MapperRectangular( mapper_grids=mapper_grids, - over_sampler=make_over_sampler_2d_7x7(), border_relocator=make_border_relocator_2d_7x7(), regularization=make_regularization_constant(), ) @@ -446,7 +434,7 @@ def make_rectangular_mapper_7x7_3x3(): def make_delaunay_mapper_9_3x3(): mapper_grids = aa.MapperGrids( mask=make_mask_2d_7x7(), - source_plane_data_grid=make_over_sampler_2d_7x7().over_sampled_grid, + source_plane_data_grid=make_grid_2d_sub_2_7x7(), source_plane_mesh_grid=make_delaunay_mesh_grid_9(), image_plane_mesh_grid=aa.Grid2D.uniform(shape_native=(3, 3), pixel_scales=0.1), adapt_data=aa.Array2D.ones(shape_native=(3, 3), pixel_scales=0.1), @@ -454,7 +442,6 @@ def make_delaunay_mapper_9_3x3(): return aa.MapperDelaunay( mapper_grids=mapper_grids, - over_sampler=make_over_sampler_2d_7x7(), border_relocator=make_border_relocator_2d_7x7(), regularization=make_regularization_constant(), ) @@ -463,7 +450,7 @@ def make_delaunay_mapper_9_3x3(): def make_voronoi_mapper_9_3x3(): mapper_grids = aa.MapperGrids( mask=make_mask_2d_7x7(), - source_plane_data_grid=make_over_sampler_2d_7x7().over_sampled_grid, + source_plane_data_grid=make_grid_2d_sub_2_7x7(), source_plane_mesh_grid=make_voronoi_mesh_grid_9(), image_plane_mesh_grid=aa.Grid2D.uniform(shape_native=(3, 3), pixel_scales=0.1), adapt_data=aa.Array2D.ones(shape_native=(3, 3), pixel_scales=0.1), @@ -471,7 +458,6 @@ def make_voronoi_mapper_9_3x3(): return aa.MapperVoronoi( mapper_grids=mapper_grids, - over_sampler=make_over_sampler_2d_7x7(), border_relocator=make_border_relocator_2d_7x7(), regularization=make_regularization_constant(), ) diff --git a/autoarray/inversion/inversion/abstract.py b/autoarray/inversion/inversion/abstract.py index 1bca77c00..757a24ef0 100644 --- a/autoarray/inversion/inversion/abstract.py +++ b/autoarray/inversion/inversion/abstract.py @@ -739,16 +739,55 @@ def log_det_regularization_matrix_term(self) -> float: raise exc.InversionException() from e @property - def errors_with_covariance(self) -> np.ndarray: - return np.linalg.inv(self.curvature_reg_matrix) + def reconstruction_noise_map_with_covariance(self) -> np.ndarray: + """ + Returns the noise-map of the reconstruction as a two dimension matrix which accounts for the covariance + of the noise between pixels. + + The diagonal of this matrix is the noise-map of the reconstruction, which can be used for analysing the + reconstruction with noise properties that are representative of the fit and therefore should be used + for any scientific analysis (e.g. source reconstructions of strong lenses). + + This noise-map is defined as the RMS standard deviation of the noise in every pixel of the reconstruction. + This definition is identical to the `noise_map` attributes of dataset objects. + + It is computed as the square root of the inverse of the curvature matrix with regularization, which is the + same matrix used to solve for the reconstruction via the linear inversion. + + Returns + ------- + The noise-map of the reconstruction as a two dimension matrix which accounts for the covariance of the noise + between pixels. + """ + return np.sqrt(np.linalg.inv(self.curvature_reg_matrix)) @property - def errors(self): - return np.diagonal(self.errors_with_covariance) + def reconstruction_noise_map(self): + """ + Returns the noise-map of the reconstruction as a one dimensional ndarray, which does not account for the + covariance of the noise between pixels. + + This matrix is representative of the noise properties of the fit and should be used for any scientific + analysis (e.g. source reconstructions of strong lenses). + + The noise-map of the reconstruction is the RMS standard deviation of the noise in every pixel of the + reconstruction. This definition is identical to the `noise_map` attributes of dataset objects. + + It is computed as the square root of the diagonal of the `reconstruction_noise_map_with_covariance` matrix, + which is the same matrix used to solve for the reconstruction via the linear inversion. + + Returns + ------- + The noise-map of the reconstruction as a one dimensional ndarray, which does not account for the covariance + of the noise between pixels. + """ + return np.diagonal(self.reconstruction_noise_map_with_covariance) @property - def errors_dict(self) -> Dict[LinearObj, np.ndarray]: - return self.source_quantity_dict_from(source_quantity=self.errors) + def reconstruction_noise_map_dict(self) -> Dict[LinearObj, np.ndarray]: + return self.source_quantity_dict_from( + source_quantity=self.reconstruction_noise_map + ) def regularization_weights_from(self, index: int) -> np.ndarray: linear_obj = self.linear_obj_list[index] diff --git a/autoarray/inversion/inversion/dataset_interface.py b/autoarray/inversion/inversion/dataset_interface.py index 660253552..abc780411 100644 --- a/autoarray/inversion/inversion/dataset_interface.py +++ b/autoarray/inversion/inversion/dataset_interface.py @@ -66,4 +66,4 @@ def __init__( @property def mask(self): - return self.grids.uniform.mask + return self.grids.lp.mask diff --git a/autoarray/inversion/inversion/interferometer/lop.py b/autoarray/inversion/inversion/interferometer/lop.py index 2c8f23f77..fdd6b8adf 100644 --- a/autoarray/inversion/inversion/interferometer/lop.py +++ b/autoarray/inversion/inversion/interferometer/lop.py @@ -142,5 +142,5 @@ def log_det_curvature_reg_matrix_term(self): ) @property - def errors(self): + def reconstruction_noise_map(self): return None diff --git a/autoarray/inversion/inversion/mapper_valued.py b/autoarray/inversion/inversion/mapper_valued.py index 273753886..175046b9d 100644 --- a/autoarray/inversion/inversion/mapper_valued.py +++ b/autoarray/inversion/inversion/mapper_valued.py @@ -29,7 +29,7 @@ def __init__(self, mapper, values, mesh_pixel_mask: Optional[np.ndarray] = None) The `Mapper` object which pairs with the values, for example a `MapperVoronoi` object. values The values of each pixel of the mapper, which could be the `reconstruction` values of an `Inversion`, - but alternatively could be other quantities such as the errors on these values. + but alternatively could be other quantities such as the noise-map of these values. mesh_pixel_mask The mask of pixels that are omitted from the reconstruction when computing the image, for example to remove pixels with low signal-to-noise so they do not impact the magnification calculation. @@ -259,7 +259,7 @@ def magnification_via_mesh_from( def magnification_via_interpolation_from( self, - shape_native: Tuple[int, int] = (401, 401), + shape_native: Tuple[int, int] = (201, 201), extent: Optional[Tuple[float, float, float, float]] = None, ) -> float: """ diff --git a/autoarray/inversion/linear_obj/func_list.py b/autoarray/inversion/linear_obj/func_list.py index 547a9427b..ce72bdcbe 100644 --- a/autoarray/inversion/linear_obj/func_list.py +++ b/autoarray/inversion/linear_obj/func_list.py @@ -96,7 +96,7 @@ 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_sampling.sub_size) + 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 diff --git a/autoarray/inversion/mock/mock_inversion.py b/autoarray/inversion/mock/mock_inversion.py index b221163e3..1fb125861 100644 --- a/autoarray/inversion/mock/mock_inversion.py +++ b/autoarray/inversion/mock/mock_inversion.py @@ -24,8 +24,8 @@ def __init__( reconstruction_dict: List[np.ndarray] = None, mapped_reconstructed_data_dict=None, mapped_reconstructed_image_dict=None, - errors: np.ndarray = None, - errors_dict: List[np.ndarray] = None, + reconstruction_noise_map: np.ndarray = None, + reconstruction_noise_map_dict: List[np.ndarray] = None, regularization_term=None, log_det_curvature_reg_matrix_term=None, log_det_regularization_matrix_term=None, @@ -56,8 +56,8 @@ def __init__( self._mapped_reconstructed_data_dict = mapped_reconstructed_data_dict self._mapped_reconstructed_image_dict = mapped_reconstructed_image_dict - self._errors = errors - self._errors_dict = errors_dict + self._reconstruction_noise_map = reconstruction_noise_map + self._reconstruction_noise_map_dict = reconstruction_noise_map_dict self._regularization_term = regularization_term self._log_det_curvature_reg_matrix_term = log_det_curvature_reg_matrix_term @@ -168,16 +168,16 @@ def mapped_reconstructed_image_dict(self): return self._mapped_reconstructed_image_dict @property - def errors(self): - if self._errors is None: - return super().errors - return self._errors + def reconstruction_noise_map(self): + if self._reconstruction_noise_map is None: + return super().reconstruction_noise_map + return self._reconstruction_noise_map @property - def errors_dict(self): - if self._errors_dict is None: - return super().errors_dict - return self._errors_dict + def reconstruction_noise_map_dict(self): + if self._reconstruction_noise_map_dict is None: + return super().reconstruction_noise_map_dict + return self._reconstruction_noise_map_dict @property def regularization_term(self): diff --git a/autoarray/inversion/mock/mock_mapper.py b/autoarray/inversion/mock/mock_mapper.py index 6d98aa174..84f904250 100644 --- a/autoarray/inversion/mock/mock_mapper.py +++ b/autoarray/inversion/mock/mock_mapper.py @@ -32,11 +32,11 @@ def __init__( super().__init__( mapper_grids=mapper_grids, - over_sampler=over_sampler, border_relocator=border_relocator, regularization=regularization, ) + self._over_sampler = over_sampler self._edge_pixel_list = edge_pixel_list self._pix_sub_weights = pix_sub_weights self._pix_sub_weights_split_cross = pix_sub_weights_split_cross @@ -56,6 +56,12 @@ def params(self): return super().params return self._parameters + @property + def over_sampler(self): + if self._over_sampler is None: + return super().over_sampler + return self._over_sampler + @property def edge_pixel_list(self): return self._edge_pixel_list diff --git a/autoarray/inversion/pixelization/border_relocator.py b/autoarray/inversion/pixelization/border_relocator.py index 8af010a7f..73a8d0d28 100644 --- a/autoarray/inversion/pixelization/border_relocator.py +++ b/autoarray/inversion/pixelization/border_relocator.py @@ -132,12 +132,47 @@ class BorderRelocator: def __init__(self, mask: Mask2D, sub_size: Union[int, Array2D]): self.mask = mask - if isinstance(sub_size, int): - sub_size = Array2D( - values=np.full(fill_value=sub_size, shape=mask.shape_slim), mask=mask + self.sub_size = over_sample_util.over_sample_size_convert_to_array_2d_from( + over_sample_size=sub_size, mask=mask + ) + + @cached_property + def border_slim(self): + """ + Returns the 1D ``slim`` indexes of border pixels in the ``Mask2D``, representing all unmasked + sub-pixels (given by ``False``) which neighbor any masked value (give by ``True``) and which are on the + extreme exterior of the mask. + + The indexes are the extended below to form the ``sub_border_slim`` which is illustrated above. + + This quantity is too complicated to write-out in a docstring, and it is recommended you print it in + Python code to understand it if anything is unclear. + + Examples + -------- + + .. code-block:: python + + import autoarray as aa + + mask_2d = aa.Mask2D( + mask=[[True, True, True, True, True, True, True, True, True], + [True, False, False, False, False, False, False, False, True], + [True, False, True, True, True, True, True, False, True], + [True, False, True, False, False, False, True, False, True], + [True, False, True, False, True, False, True, False, True], + [True, False, True, False, False, False, True, False, True], + [True, False, True, True, True, True, True, False, True], + [True, False, False, False, False, False, False, False, True], + [True, True, True, True, True, True, True, True, True]] + pixel_scales=1.0, ) - self.sub_size = sub_size + derive_indexes_2d = aa.DeriveIndexes2D(mask=mask_2d) + + print(derive_indexes_2d.border_slim) + """ + return self.mask.derive_indexes.border_slim @cached_property def sub_border_slim(self) -> np.ndarray: @@ -176,18 +211,9 @@ def sub_border_slim(self) -> np.ndarray: print(derive_indexes_2d.sub_border_slim) """ return sub_border_pixel_slim_indexes_from( - mask_2d=np.array(self.mask), sub_size=self.sub_size + mask_2d=np.array(self.mask), sub_size=np.array(self.sub_size).astype("int") ).astype("int") - @property - def sub_grid(self): - return over_sample_util.grid_2d_slim_over_sampled_via_mask_from( - mask_2d=np.array(self.mask), - pixel_scales=self.mask.pixel_scales, - sub_size=np.array(self.sub_size), - origin=self.mask.origin, - ) - @cached_property def border_grid(self) -> np.ndarray: """ @@ -206,9 +232,16 @@ def sub_border_grid(self) -> np.ndarray: This is NOT all sub-pixels which are in mask pixels at the mask's border, but specifically the sub-pixels within these border pixels which are at the extreme edge of the border. """ - return self.sub_grid[self.sub_border_slim] + sub_grid = over_sample_util.grid_2d_slim_over_sampled_via_mask_from( + mask_2d=np.array(self.mask), + pixel_scales=self.mask.pixel_scales, + sub_size=np.array(self.sub_size).astype("int"), + origin=self.mask.origin, + ) - def relocated_grid_from(self, grid: Grid2DIrregular) -> Grid2DIrregular: + return sub_grid[self.sub_border_slim] + + def relocated_grid_from(self, grid: Grid2D) -> Grid2D: """ Relocate the coordinates of a grid to the border of this grid 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*). @@ -235,11 +268,21 @@ def relocated_grid_from(self, grid: Grid2DIrregular) -> Grid2DIrregular: if len(self.sub_border_grid) == 0: return grid - return Grid2DIrregular( - values=grid_2d_util.relocated_grid_via_jit_from( - grid=np.array(grid), - border_grid=np.array(grid[self.sub_border_slim]), - ), + values = grid_2d_util.relocated_grid_via_jit_from( + grid=np.array(grid), + border_grid=np.array(grid[self.border_slim]), + ) + + over_sampled = grid_2d_util.relocated_grid_via_jit_from( + grid=np.array(grid.over_sampled), + border_grid=np.array(grid.over_sampled[self.sub_border_slim]), + ) + + return Grid2D( + values=values, + mask=grid.mask, + over_sample_size=self.sub_size, + over_sampled=over_sampled, ) def relocated_mesh_grid_from( diff --git a/autoarray/inversion/pixelization/image_mesh/hilbert.py b/autoarray/inversion/pixelization/image_mesh/hilbert.py index 0b7cc1bb4..964457ded 100644 --- a/autoarray/inversion/pixelization/image_mesh/hilbert.py +++ b/autoarray/inversion/pixelization/image_mesh/hilbert.py @@ -9,7 +9,7 @@ from autoarray.inversion.pixelization.image_mesh.abstract_weighted import ( AbstractImageMeshWeighted, ) -from autoarray.operators.over_sampling.uniform import OverSamplerUniform +from autoarray.operators.over_sampling.over_sampler import OverSampler from autoarray.inversion.inversion.settings import SettingsInversion from autoarray.structures.grids.irregular_2d import Grid2DIrregular @@ -87,42 +87,6 @@ def generate2d(x, y, ax, ay, bx, by): ) -def super_resolution_grid_from(img_2d, mask, mask_radius, pixel_scales, sub_scale=11): - """ - This function will create a higher resolution grid for the img_2d. The new grid and its - interpolated values will be used to generate a sparse image grid. - - img_2d: the hyper image in 2d (e.g. hyper_source_model_image.native) - mask: the mask used for the fitting. - mask_radius: the circular mask radius. Currently, the code only works with a circular mask. - sub_scale: oversampling scale for each image pixel. - """ - - shape_nnn = np.shape(mask)[0] - - grid = Grid2D.uniform( - shape_native=(shape_nnn, shape_nnn), - pixel_scales=pixel_scales, - ) - - new_mask = Mask2D.circular( - shape_native=(shape_nnn, shape_nnn), - pixel_scales=pixel_scales, - centre=mask.origin, - radius=mask_radius, - ) - - over_sampler = OverSamplerUniform(mask=new_mask, sub_size=sub_scale) - - new_grid = over_sampler.over_sampled_grid - - new_img = griddata( - points=grid, values=img_2d.ravel(), xi=new_grid, fill_value=0.0, method="linear" - ) - - return new_img, new_grid - - def grid_hilbert_order_from(length, mask_radius): """ This function will create a grid in the Hilbert space-filling curve order. diff --git a/autoarray/inversion/pixelization/mappers/abstract.py b/autoarray/inversion/pixelization/mappers/abstract.py index 6d6107585..fd0ec1038 100644 --- a/autoarray/inversion/pixelization/mappers/abstract.py +++ b/autoarray/inversion/pixelization/mappers/abstract.py @@ -11,7 +11,6 @@ from autoarray.inversion.pixelization.border_relocator import BorderRelocator from autoarray.inversion.pixelization.mappers.mapper_grids import MapperGrids from autoarray.inversion.regularization.abstract import AbstractRegularization -from autoarray.operators.over_sampling.abstract import AbstractOverSampler from autoarray.structures.arrays.uniform_2d import Array2D from autoarray.structures.grids.uniform_2d import Grid2D from autoarray.structures.mesh.abstract_2d import Abstract2DMesh @@ -26,7 +25,6 @@ def __init__( self, mapper_grids: MapperGrids, regularization: Optional[AbstractRegularization], - over_sampler: AbstractOverSampler, border_relocator: BorderRelocator, run_time_dict: Optional[Dict] = None, ): @@ -82,9 +80,6 @@ def __init__( 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. - 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. @@ -94,7 +89,6 @@ def __init__( super().__init__(regularization=regularization, run_time_dict=run_time_dict) - self.over_sampler = over_sampler self.border_relocator = border_relocator self.mapper_grids = mapper_grids @@ -118,6 +112,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 over_sampler(self): + return self.mapper_grids.source_plane_data_grid.over_sampler + @property def edge_pixel_list(self) -> List[int]: return self.source_plane_mesh_grid.edge_pixel_list @@ -262,7 +260,7 @@ def unique_mappings(self) -> UniqueMappings: pix_sizes_for_sub_slim_index=self.pix_sizes_for_sub_slim_index, pix_weights_for_sub_slim_index=self.pix_weights_for_sub_slim_index, pix_pixels=self.params, - sub_size=np.array(self.over_sampler.sub_size), + sub_size=np.array(self.over_sampler.sub_size).astype("int"), ) return UniqueMappings( @@ -394,12 +392,16 @@ def mapped_to_source_from(self, array: Array2D) -> np.ndarray: ) def extent_from( - self, values: np.ndarray = None, zoom_to_brightest: bool = True + self, + values: np.ndarray = None, + zoom_to_brightest: bool = True, + zoom_percent: Optional[float] = None, ) -> Tuple[float, float, float, float]: if zoom_to_brightest and values is not None: - zoom_percent = conf.instance["visualize"]["general"]["zoom"][ - "inversion_percent" - ] + if zoom_percent is None: + zoom_percent = conf.instance["visualize"]["general"]["zoom"][ + "inversion_percent" + ] fractional_value = np.max(values) * zoom_percent fractional_bool = values > fractional_value diff --git a/autoarray/inversion/pixelization/mappers/delaunay.py b/autoarray/inversion/pixelization/mappers/delaunay.py index c35c4d9b7..737247b0b 100644 --- a/autoarray/inversion/pixelization/mappers/delaunay.py +++ b/autoarray/inversion/pixelization/mappers/delaunay.py @@ -112,12 +112,12 @@ def pix_sub_weights(self) -> PixSubWeights: delaunay = self.delaunay simplex_index_for_sub_slim_index = delaunay.find_simplex( - self.source_plane_data_grid + self.source_plane_data_grid.over_sampled ) pix_indexes_for_simplex_index = delaunay.simplices mappings, sizes = mapper_util.pix_indexes_for_sub_slim_index_delaunay_from( - source_plane_data_grid=np.array(self.source_plane_data_grid), + source_plane_data_grid=np.array(self.source_plane_data_grid.over_sampled), simplex_index_for_sub_slim_index=simplex_index_for_sub_slim_index, pix_indexes_for_simplex_index=pix_indexes_for_simplex_index, delaunay_points=delaunay.points, @@ -127,7 +127,7 @@ def pix_sub_weights(self) -> PixSubWeights: sizes = sizes.astype("int") weights = mapper_util.pixel_weights_delaunay_from( - source_plane_data_grid=np.array(self.source_plane_data_grid), + source_plane_data_grid=np.array(self.source_plane_data_grid.over_sampled), source_plane_mesh_grid=np.array(self.source_plane_mesh_grid), slim_index_for_sub_slim_index=self.slim_index_for_sub_slim_index, pix_indexes_for_sub_slim_index=mappings, diff --git a/autoarray/inversion/pixelization/mappers/factory.py b/autoarray/inversion/pixelization/mappers/factory.py index 8cdfa721e..310133202 100644 --- a/autoarray/inversion/pixelization/mappers/factory.py +++ b/autoarray/inversion/pixelization/mappers/factory.py @@ -2,7 +2,6 @@ from autoarray.inversion.pixelization.mappers.mapper_grids import MapperGrids from autoarray.inversion.pixelization.border_relocator import BorderRelocator -from autoarray.operators.over_sampling.abstract import AbstractOverSampler from autoarray.inversion.regularization.abstract import AbstractRegularization from autoarray.structures.mesh.rectangular_2d import Mesh2DRectangular from autoarray.structures.mesh.delaunay_2d import Mesh2DDelaunay @@ -12,7 +11,6 @@ def mapper_from( mapper_grids: MapperGrids, regularization: Optional[AbstractRegularization], - over_sampler: AbstractOverSampler, border_relocator: Optional[BorderRelocator] = None, run_time_dict: Optional[Dict] = None, ): @@ -50,7 +48,6 @@ def mapper_from( if isinstance(mapper_grids.source_plane_mesh_grid, Mesh2DRectangular): return MapperRectangular( mapper_grids=mapper_grids, - over_sampler=over_sampler, border_relocator=border_relocator, regularization=regularization, run_time_dict=run_time_dict, @@ -58,7 +55,6 @@ def mapper_from( elif isinstance(mapper_grids.source_plane_mesh_grid, Mesh2DDelaunay): return MapperDelaunay( mapper_grids=mapper_grids, - over_sampler=over_sampler, border_relocator=border_relocator, regularization=regularization, run_time_dict=run_time_dict, @@ -66,7 +62,6 @@ def mapper_from( elif isinstance(mapper_grids.source_plane_mesh_grid, Mesh2DVoronoi): return MapperVoronoi( mapper_grids=mapper_grids, - over_sampler=over_sampler, border_relocator=border_relocator, regularization=regularization, run_time_dict=run_time_dict, diff --git a/autoarray/inversion/pixelization/mappers/rectangular.py b/autoarray/inversion/pixelization/mappers/rectangular.py index 815e19227..9a78c1b8a 100644 --- a/autoarray/inversion/pixelization/mappers/rectangular.py +++ b/autoarray/inversion/pixelization/mappers/rectangular.py @@ -100,7 +100,7 @@ def pix_sub_weights(self) -> PixSubWeights: are equal to 1.0. """ mappings = geometry_util.grid_pixel_indexes_2d_slim_from( - grid_scaled_2d_slim=np.array(self.source_plane_data_grid), + grid_scaled_2d_slim=np.array(self.source_plane_data_grid.over_sampled), shape_native=self.source_plane_mesh_grid.shape_native, pixel_scales=self.source_plane_mesh_grid.pixel_scales, origin=self.source_plane_mesh_grid.origin, @@ -109,7 +109,9 @@ def pix_sub_weights(self) -> PixSubWeights: mappings = mappings.reshape((len(mappings), 1)) return PixSubWeights( - mappings=mappings.reshape((len(mappings), 1)), + mappings=mappings, sizes=np.ones(len(mappings), dtype="int"), - weights=np.ones((len(self.source_plane_data_grid), 1), dtype="int"), + weights=np.ones( + (len(self.source_plane_data_grid.over_sampled), 1), dtype="int" + ), ) diff --git a/autoarray/inversion/pixelization/mappers/voronoi.py b/autoarray/inversion/pixelization/mappers/voronoi.py index 766a9bcbf..c8e54cbf3 100644 --- a/autoarray/inversion/pixelization/mappers/voronoi.py +++ b/autoarray/inversion/pixelization/mappers/voronoi.py @@ -130,7 +130,8 @@ def pix_sub_weights(self) -> PixSubWeights: """ mappings, sizes, weights = mapper_util.pix_size_weights_voronoi_nn_from( - grid=self.source_plane_data_grid, mesh_grid=self.source_plane_mesh_grid + grid=self.source_plane_data_grid.over_sampled, + mesh_grid=self.source_plane_mesh_grid, ) mappings = mappings.astype("int") diff --git a/autoarray/inversion/pixelization/mesh/rectangular.py b/autoarray/inversion/pixelization/mesh/rectangular.py index 514e5b6d4..b1b5a7017 100644 --- a/autoarray/inversion/pixelization/mesh/rectangular.py +++ b/autoarray/inversion/pixelization/mesh/rectangular.py @@ -104,6 +104,7 @@ def mapper_grids_from( border_relocator=border_relocator, source_plane_data_grid=source_plane_data_grid, ) + mesh_grid = self.mesh_grid_from(source_plane_data_grid=relocated_grid) return MapperGrids( @@ -135,7 +136,7 @@ def mesh_grid_from( by overlaying the `source_plane_data_grid` with the rectangular pixelization. """ return Mesh2DRectangular.overlay_grid( - shape_native=self.shape, grid=source_plane_data_grid + shape_native=self.shape, grid=source_plane_data_grid.over_sampled ) @property diff --git a/autoarray/inversion/pixelization/mesh/triangulation.py b/autoarray/inversion/pixelization/mesh/triangulation.py index 17b7536f8..a72236e50 100644 --- a/autoarray/inversion/pixelization/mesh/triangulation.py +++ b/autoarray/inversion/pixelization/mesh/triangulation.py @@ -65,28 +65,28 @@ def mapper_grids_from( self.run_time_dict = run_time_dict - source_plane_data_grid = self.relocated_grid_from( + relocated_grid = self.relocated_grid_from( border_relocator=border_relocator, source_plane_data_grid=source_plane_data_grid, ) - relocated_source_plane_mesh_grid = self.relocated_mesh_grid_from( + relocated_mesh_grid = self.relocated_mesh_grid_from( border_relocator=border_relocator, - source_plane_data_grid=source_plane_data_grid, + source_plane_data_grid=relocated_grid.over_sampled, source_plane_mesh_grid=source_plane_mesh_grid, ) try: source_plane_mesh_grid = self.mesh_grid_from( - source_plane_data_grid=source_plane_data_grid, - source_plane_mesh_grid=relocated_source_plane_mesh_grid, + source_plane_data_grid=relocated_grid.over_sampled, + source_plane_mesh_grid=relocated_mesh_grid, ) except ValueError as e: raise e return MapperGrids( mask=mask, - source_plane_data_grid=source_plane_data_grid, + source_plane_data_grid=relocated_grid, source_plane_mesh_grid=source_plane_mesh_grid, image_plane_mesh_grid=image_plane_mesh_grid, adapt_data=adapt_data, diff --git a/autoarray/inversion/pixelization/mesh/voronoi.py b/autoarray/inversion/pixelization/mesh/voronoi.py index ef76e7c1a..dc8d9310f 100644 --- a/autoarray/inversion/pixelization/mesh/voronoi.py +++ b/autoarray/inversion/pixelization/mesh/voronoi.py @@ -55,7 +55,6 @@ def mesh_grid_from( settings Settings controlling the pixelization for example if a border is used to relocate its exterior coordinates. """ - return Mesh2DVoronoi( values=source_plane_mesh_grid, ) diff --git a/autoarray/inversion/plot/inversion_plotters.py b/autoarray/inversion/plot/inversion_plotters.py index f3d58b635..98cede938 100644 --- a/autoarray/inversion/plot/inversion_plotters.py +++ b/autoarray/inversion/plot/inversion_plotters.py @@ -121,7 +121,7 @@ def figures_2d_of_pixelization( data_subtracted: bool = False, reconstructed_image: bool = False, reconstruction: bool = False, - errors: bool = False, + reconstruction_noise_map: bool = False, signal_to_noise_map: bool = False, regularization_weights: bool = False, sub_pixels_per_image_pixels: bool = False, @@ -145,8 +145,8 @@ def figures_2d_of_pixelization( Whether to make a 2D plot (via `imshow`) of the mapper's reconstructed image data. reconstruction Whether to make a 2D plot (via `imshow` or `fill`) of the mapper's source-plane reconstruction. - errors - Whether to make a 2D plot (via `imshow` or `fill`) of the mapper's source-plane errors. + reconstruction_noise_map + Whether to make a 2D plot (via `imshow` or `fill`) of the mapper's source-plane noise-map. signal_to_noise_map Whether to make a 2D plot (via `imshow` or `fill`) of the mapper's source-plane signal-to-noise-map. sub_pixels_per_image_pixels @@ -184,7 +184,7 @@ def figures_2d_of_pixelization( self.mat_plot_2d.plot_array( array=array, visuals_2d=self.get_visuals_2d_for_data(), - grid_indexes=mapper_plotter.mapper.over_sampler.over_sampled_grid, + grid_indexes=mapper_plotter.mapper.over_sampler.uniform_over_sampled, auto_labels=AutoLabels( title="Data Subtracted", filename="data_subtracted" ), @@ -200,7 +200,7 @@ def figures_2d_of_pixelization( self.mat_plot_2d.plot_array( array=array, visuals_2d=self.get_visuals_2d_for_data(), - grid_indexes=mapper_plotter.mapper.over_sampler.over_sampled_grid, + grid_indexes=mapper_plotter.mapper.over_sampler.uniform_over_sampled, auto_labels=AutoLabels( title="Reconstructed Image", filename="reconstructed_image" ), @@ -235,11 +235,15 @@ def figures_2d_of_pixelization( if vmax_custom: self.mat_plot_2d.cmap.kwargs["vmax"] = None - if errors: + if reconstruction_noise_map: try: mapper_plotter.plot_source_from( - pixel_values=self.inversion.errors_dict[mapper_plotter.mapper], - auto_labels=AutoLabels(title="Errors", filename="errors"), + pixel_values=self.inversion.reconstruction_noise_map_dict[ + mapper_plotter.mapper + ], + auto_labels=AutoLabels( + title="Noise Map", filename="reconstruction_noise_map" + ), ) except TypeError: @@ -249,7 +253,9 @@ def figures_2d_of_pixelization( try: signal_to_noise_values = ( self.inversion.reconstruction_dict[mapper_plotter.mapper] - / self.inversion.errors_dict[mapper_plotter.mapper] + / self.inversion.reconstruction_noise_map_dict[ + mapper_plotter.mapper + ] ) mapper_plotter.plot_source_from( @@ -387,9 +393,11 @@ def subplot_of_mapper( ) self.set_title(label=None) - self.set_title(label="Errors (Unzoomed)") + self.set_title(label="Noise-Map (Unzoomed)") self.figures_2d_of_pixelization( - pixelization_index=mapper_index, errors=True, zoom_to_brightest=False + pixelization_index=mapper_index, + reconstruction_noise_map=True, + zoom_to_brightest=False, ) self.set_title(label="Regularization Weights (Unzoomed)") diff --git a/autoarray/inversion/plot/mapper_plotters.py b/autoarray/inversion/plot/mapper_plotters.py index d6fa791d6..9fd6608d7 100644 --- a/autoarray/inversion/plot/mapper_plotters.py +++ b/autoarray/inversion/plot/mapper_plotters.py @@ -80,7 +80,7 @@ def figure_2d( interpolate_to_uniform=interpolate_to_uniform, pixel_values=solution_vector, auto_labels=AutoLabels( - title="Pixelization Mesh (Image-Plane)", filename="mapper" + title="Pixelization Mesh (Source-Plane)", filename="mapper" ), ) @@ -117,7 +117,7 @@ def subplot_image_and_mapper( ) self.mat_plot_2d.index_scatter.scatter_grid_indexes( - grid=self.mapper.over_sampler.over_sampled_grid, + grid=self.mapper.over_sampler.uniform_over_sampled, indexes=indexes, ) diff --git a/autoarray/mask/abstract_mask.py b/autoarray/mask/abstract_mask.py index fabee9c96..bc80d6a1b 100644 --- a/autoarray/mask/abstract_mask.py +++ b/autoarray/mask/abstract_mask.py @@ -76,24 +76,6 @@ def pixel_scale(self) -> float: For a mask with dimensions two or above check that are pixel scales are the same, and if so return this single value as a float. """ - def exception_message(): - raise exc.MaskException( - "Cannot return a pixel_scale for a grid where each dimension has a " - "different pixel scale (e.g. pixel_scales[0] != pixel_scales[1])" - ) - - for pixel_scale in self.pixel_scales: - cond = abs(pixel_scale - self.pixel_scales[0]) > 1.0e-8 - if use_jax: - jax.lax.cond( - cond, - lambda _: jax.debug.callback(exception_message), - lambda _: None, - None - ) - elif cond: - exception_message() - return self.pixel_scales[0] @property diff --git a/autoarray/mock.py b/autoarray/mock.py index a9344b516..78a857e7e 100644 --- a/autoarray/mock.py +++ b/autoarray/mock.py @@ -21,4 +21,3 @@ from autoarray.structures.mock.mock_decorators import MockGridRadialMinimum from autoarray.structures.mock.mock_decorators import MockGrid1DLikeObj from autoarray.structures.mock.mock_decorators import MockGrid2DLikeObj -from autoarray.structures.mock.mock_decorators import MockGridLikeIteratorObj diff --git a/autoarray/operators/over_sampling/abstract.py b/autoarray/operators/over_sampling/abstract.py deleted file mode 100644 index a29dc1177..000000000 --- a/autoarray/operators/over_sampling/abstract.py +++ /dev/null @@ -1,18 +0,0 @@ -from autoarray.numpy_wrapper import register_pytree_node_class - -from autoarray.mask.mask_2d import Mask2D - - -class AbstractOverSampling: - def over_sampler_from(self, mask: Mask2D) -> "AbstractOverSampler": - raise NotImplementedError() - - -@register_pytree_node_class -class AbstractOverSampler: - def tree_flatten(self): - return (self.mask,), () - - @classmethod - def tree_unflatten(cls, aux_data, children): - return cls(mask=children[0]) diff --git a/autoarray/operators/over_sampling/decorator.py b/autoarray/operators/over_sampling/decorator.py index 584baf35d..c028ee31a 100644 --- a/autoarray/operators/over_sampling/decorator.py +++ b/autoarray/operators/over_sampling/decorator.py @@ -4,9 +4,6 @@ from typing import List, Union -from autoarray.operators.over_sampling.grid_oversampled import Grid2DOverSampled -from autoarray.operators.over_sampling.uniform import OverSamplingUniform - from autoarray.structures.arrays.irregular import ArrayIrregular from autoarray.structures.arrays.uniform_1d import Array1D from autoarray.structures.arrays.uniform_2d import Array2D @@ -15,27 +12,6 @@ from autoarray.structures.grids.uniform_2d import Grid2D -def perform_over_sampling_from(grid, **kwargs): - if kwargs.get("over_sampling_being_performed"): - return False - - perform_over_sampling = False - - if isinstance(grid, Grid2D): - if grid.over_sampling is not None: - perform_over_sampling = True - - if isinstance(grid.over_sampling, OverSamplingUniform): - try: - if grid.over_sampling.sub_size == 1: - perform_over_sampling = False - except ValueError: - if sum(grid.over_sampling.sub_size) == grid.mask.pixels_in_mask: - perform_over_sampling = False - - return perform_over_sampling - - def over_sample(func): """ Homogenize the inputs and outputs of functions that take 1D or 2D grids of coordinates and return a 1D ndarray @@ -72,33 +48,14 @@ def wrapper( The function values evaluated on the grid with the same structure as the input grid_like object. """ - if isinstance(grid, Grid2DOverSampled): - result = func(obj, grid.grid, *args, **kwargs) - - return grid.over_sampler.binned_array_2d_from(array=result) - - if isinstance(grid, Grid2D): - if grid.over_sampling is None: - if grid.is_uniform: - over_sampling = OverSamplingUniform.from_adaptive_scheme( - grid=grid, - name=obj.__class__.__name__, - centre=obj.centre, - ) - - grid = Grid2D( - values=grid, mask=grid.mask, over_sampling=over_sampling - ) - - perform_over_sampling = perform_over_sampling_from(grid=grid, kwargs=kwargs) - - if not perform_over_sampling: + if isinstance(grid, Grid2DIrregular) or isinstance(grid, Grid1D): return func(obj=obj, grid=grid, *args, **kwargs) - kwargs["over_sampling_being_performed"] = True + if obj is not None: + values = func(obj, grid.over_sampled, *args, **kwargs) + else: + values = func(grid.over_sampled, *args, **kwargs) - return grid.over_sampler.array_via_func_from( - func=func, obj=obj, *args, **kwargs - ) + return grid.over_sampler.binned_array_2d_from(array=values) return wrapper diff --git a/autoarray/operators/over_sampling/grid_oversampled.py b/autoarray/operators/over_sampling/grid_oversampled.py deleted file mode 100644 index e3b0e6667..000000000 --- a/autoarray/operators/over_sampling/grid_oversampled.py +++ /dev/null @@ -1,9 +0,0 @@ -class Grid2DOverSampled: - def __init__(self, grid, over_sampler, pixels_in_mask): - self.grid = grid - self.over_sampler = over_sampler - self.pixels_in_mask = pixels_in_mask - - @property - def mask(self): - return self.grid.mask diff --git a/autoarray/operators/over_sampling/iterate.py b/autoarray/operators/over_sampling/iterate.py deleted file mode 100644 index 57540bfc6..000000000 --- a/autoarray/operators/over_sampling/iterate.py +++ /dev/null @@ -1,295 +0,0 @@ -import numpy as np -from typing import Callable, List, Optional - -from autoarray import numba_util -from autoarray.mask.mask_2d import Mask2D -from autoarray.operators.over_sampling.abstract import AbstractOverSampling -from autoarray.operators.over_sampling.abstract import AbstractOverSampler -from autoarray.operators.over_sampling.uniform import OverSamplerUniform -from autoarray.structures.arrays.uniform_2d import Array2D - - -class OverSamplingIterate(AbstractOverSampling): - def __init__( - self, - fractional_accuracy: float = 0.9999, - relative_accuracy: Optional[float] = None, - sub_steps: List[int] = None, - ): - """ - Over samples grid calculations using an iterative sub-grid that increases the sampling until a threshold - accuracy is met. - - When a 2D grid of (y,x) coordinates is input into a function, the result is evaluated at every coordinate - on the grid. When the grid is paired to a 2D image (e.g. an `Array2D`) the solution needs to approximate - the 2D integral of that function in each pixel. Over sample objects define how this over-sampling is performed. - - This object iteratively recomputes the analytic function at increasing sub-grid resolutions until an input - fractional accuracy is reached. The sub-grid is increase in each pixel, therefore it will gradually better - approximate the 2D integral after each iteration. - - Iteration is performed on a per pixel basis, such that the sub-grid size will stop at lower values - in pixels where the fractional accuracy is met quickly. It will only go to high values where high sampling is - required to meet the accuracy. This ensures the function is evaluated accurately in a computationally - efficient manner. - - Parameters - ---------- - fractional_accuracy - The fractional accuracy the function evaluated must meet to be accepted, where this accuracy is the ratio - of the value at a higher sub size to the value computed using the previous sub_size. The fractional - accuracy does not depend on the units or magnitude of the function being evaluated. - relative_accuracy - The relative accuracy the function evaluted must meet to be accepted, where this value is the absolute - difference of the values computed using the higher sub size and lower sub size grids. The relative - accuracy depends on the units / magnitude of the function being evaluated. - sub_steps - The sub-size values used to iteratively evaluated the function at high levels of sub-gridding. If None, - they are setup as the default values [2, 4, 8, 16]. - """ - - if sub_steps is None: - sub_steps = [2, 4, 8, 16] - - self.fractional_accuracy = fractional_accuracy - self.relative_accuracy = relative_accuracy - self.sub_steps = sub_steps - - def over_sampler_from(self, mask: Mask2D) -> "OverSamplerIterate": - return OverSamplerIterate( - mask=mask, - sub_steps=self.sub_steps, - fractional_accuracy=self.fractional_accuracy, - relative_accuracy=self.relative_accuracy, - ) - - -@numba_util.jit() -def threshold_mask_via_arrays_jit_from( - fractional_accuracy_threshold: float, - relative_accuracy_threshold: Optional[float], - threshold_mask: np.ndarray, - array_higher_sub_2d: np.ndarray, - array_lower_sub_2d: np.ndarray, - array_higher_mask: np.ndarray, -) -> np.ndarray: - """ - Jitted function to determine the fractional mask, which is a mask where: - - - ``True`` entries signify the function has been evaluated in that pixel to desired accuracy and - therefore does not need to be iteratively computed at higher levels of sub-gridding. - - - ``False`` entries signify the function has not been evaluated in that pixel to desired fractional accuracy and - therefore must be iterative computed at higher levels of sub-gridding to meet this accuracy. - """ - - if fractional_accuracy_threshold is not None: - for y in range(threshold_mask.shape[0]): - for x in range(threshold_mask.shape[1]): - if not array_higher_mask[y, x]: - if array_lower_sub_2d[y, x] > 0: - fractional_accuracy = ( - array_lower_sub_2d[y, x] / array_higher_sub_2d[y, x] - ) - - if fractional_accuracy > 1.0: - fractional_accuracy = 1.0 / fractional_accuracy - - else: - fractional_accuracy = 0.0 - - if fractional_accuracy < fractional_accuracy_threshold: - threshold_mask[y, x] = False - - if relative_accuracy_threshold is not None: - for y in range(threshold_mask.shape[0]): - for x in range(threshold_mask.shape[1]): - if not array_higher_mask[y, x]: - if ( - abs(array_lower_sub_2d[y, x] - array_higher_sub_2d[y, x]) - > relative_accuracy_threshold - ): - threshold_mask[y, x] = False - - return threshold_mask - - -@numba_util.jit() -def iterated_array_jit_from( - iterated_array: np.ndarray, - threshold_mask_higher_sub: np.ndarray, - threshold_mask_lower_sub: np.ndarray, - array_higher_sub_2d: np.ndarray, -) -> np.ndarray: - """ - Create the iterated array from a result array that is computed at a higher sub size leel than the previous grid. - - The iterated array is only updated for pixels where the fractional accuracy is met. - """ - - for y in range(iterated_array.shape[0]): - for x in range(iterated_array.shape[1]): - if threshold_mask_higher_sub[y, x] and not threshold_mask_lower_sub[y, x]: - iterated_array[y, x] = array_higher_sub_2d[y, x] - - return iterated_array - - -class OverSamplerIterate(AbstractOverSampler): - def __init__( - self, - mask: Mask2D, - fractional_accuracy: float = 0.9999, - relative_accuracy: Optional[float] = None, - sub_steps: List[int] = None, - ): - self.mask = mask - self.fractional_accuracy = fractional_accuracy - self.relative_accuracy = relative_accuracy - self.sub_steps = sub_steps - - def array_at_sub_size_from( - self, func: Callable, cls, mask: Mask2D, sub_size, *args, **kwargs - ) -> Array2D: - over_sample_uniform = OverSamplerUniform(mask=mask, sub_size=sub_size) - - over_sampled_grid = over_sample_uniform.over_sampled_grid - - array_higher_sub = func(cls, over_sampled_grid, *args, **kwargs) - - return over_sample_uniform.binned_array_2d_from(array=array_higher_sub).native - - def threshold_mask_from( - self, array_lower_sub_2d: Array2D, array_higher_sub_2d: Array2D - ) -> Mask2D: - """ - Returns a fractional mask from a result array, where the fractional mask describes whether the evaluated - value in the result array is within the ``OverSamplingIterate``'s specified fractional accuracy. The fractional mask thus - determines whether a pixel on the grid needs to be reevaluated at a higher level of sub-gridding to meet the - specified fractional accuracy. If it must be re-evaluated, the fractional masks's entry is ``False``. - - The fractional mask is computed by comparing the results evaluated at one level of sub-gridding to another - at a higher level of sub-griding. Thus, the sub-grid size in chosen on a per-pixel basis until the function - is evaluated at the specified fractional accuracy. - - Parameters - ---------- - array_lower_sub_2d - The results computed by a function using a lower sub-grid size - array_higher_sub_2d - The results computed by a function using a higher sub-grid size. - """ - - threshold_mask = Mask2D.all_false( - shape_native=array_lower_sub_2d.shape_native, - pixel_scales=array_lower_sub_2d.pixel_scales, - invert=True, - ) - - threshold_mask = threshold_mask_via_arrays_jit_from( - fractional_accuracy_threshold=self.fractional_accuracy, - relative_accuracy_threshold=self.relative_accuracy, - threshold_mask=np.array(threshold_mask), - array_higher_sub_2d=np.array(array_higher_sub_2d), - array_lower_sub_2d=np.array(array_lower_sub_2d), - array_higher_mask=np.array(array_higher_sub_2d.mask), - ) - - return Mask2D( - mask=threshold_mask, - pixel_scales=array_lower_sub_2d.pixel_scales, - origin=array_lower_sub_2d.origin, - ) - - def array_via_func_from( - self, func: Callable, obj: object, *args, **kwargs - ) -> Array2D: - """ - Iterate over a function that returns an array of values until the it meets a specified fractional accuracy. - The function returns a result on a pixel-grid where evaluating it on more points on a higher resolution - sub-grid followed by binning lead to a more precise evaluation of the function. The function is assumed to - belong to a class, which is input into tthe method. - - The function is first called for a sub-grid size of 1 and a higher resolution grid. The ratio of values give - the fractional accuracy of each function evaluation. Pixels which do not meet the fractional accuracy are - iteratively revaluated on higher resolution sub-grids. This is repeated until all pixels meet the fractional - accuracy or the highest sub-size specified in the *sub_steps* attribute is computed. - - If the function return all zeros, the iteration is terminated early given that all levels of sub-gridding will - return zeros. This occurs when a function is missing optional objects that contribute to the calculation. - - An example use case of this function is when a "image_2d_from" methods in **PyAutoGalaxy**'s - ``LightProfile`` module is comomputed, which by evaluating the function on a higher resolution sub-grids sample - the analytic light profile at more points and thus more precisely. - - Iterate over a function that returns an array or grid of values until the it meets a specified fractional - accuracy. The function returns a result on a pixel-grid where evaluating it on more points on a higher - resolution sub-grid followed by binning lead to a more precise evaluation of the function. - - A full description of the iteration method can be found in the functions *array_via_func_from* and - *iterated_grid_from*. This function computes the result on a grid with a sub-size of 1, and uses its - shape to call the correct function. - - Parameters - ---------- - func : func - The function which is iterated over to compute a more precise evaluation. - obj : cls - The class the function belongs to. - grid_lower_sub_2d - The results computed by the function using a lower sub-grid size - """ - - unmasked_grid = self.mask.derive_grid.unmasked - - array_sub_1 = func(obj, unmasked_grid, *args, **kwargs) - - array_sub_1 = Array2D(values=array_sub_1, mask=self.mask).native - - if not np.any(array_sub_1): - return array_sub_1.slim - - iterated_array = np.zeros(shape=self.mask.shape_native) - - threshold_mask_lower_sub = self.mask - - for sub_size in self.sub_steps[:-1]: - array_higher_sub = self.array_at_sub_size_from( - func=func, cls=obj, mask=threshold_mask_lower_sub, sub_size=sub_size - ) - - try: - threshold_mask_higher_sub = self.threshold_mask_from( - array_lower_sub_2d=array_sub_1, - array_higher_sub_2d=array_higher_sub, - ) - - iterated_array = iterated_array_jit_from( - iterated_array=iterated_array, - threshold_mask_higher_sub=np.array(threshold_mask_higher_sub), - threshold_mask_lower_sub=np.array(threshold_mask_lower_sub), - array_higher_sub_2d=np.array(array_higher_sub), - ) - - except ZeroDivisionError: - return Array2D(values=iterated_array, mask=self.mask) - - if threshold_mask_higher_sub.is_all_true: - return Array2D(values=iterated_array, mask=self.mask) - - array_sub_1 = array_higher_sub - threshold_mask_lower_sub = threshold_mask_higher_sub - threshold_mask_higher_sub.pixel_scales = self.mask.pixel_scales - - array_higher_sub = self.array_at_sub_size_from( - func=func, - cls=obj, - mask=threshold_mask_lower_sub, - sub_size=self.sub_steps[-1], - *args, - **kwargs - ) - - iterated_array_2d = iterated_array + array_higher_sub - - return Array2D(values=iterated_array_2d, mask=self.mask) diff --git a/autoarray/operators/over_sampling/over_sample_util.py b/autoarray/operators/over_sampling/over_sample_util.py index 0a8d43c43..b0df135a8 100644 --- a/autoarray/operators/over_sampling/over_sample_util.py +++ b/autoarray/operators/over_sampling/over_sample_util.py @@ -1,15 +1,57 @@ +from __future__ import annotations +import numpy as np +from typing import TYPE_CHECKING, Union, List, Tuple from autoarray.numpy_wrapper import np, register_pytree_node_class, use_jax, jit from typing import List, Tuple -from autoarray.mask.mask_2d import Mask2D +from autoarray.structures.arrays.uniform_2d import Array2D + +if TYPE_CHECKING: + from autoarray.structures.grids.uniform_2d import Grid2D from autoarray.geometry import geometry_util +from autoarray.mask.mask_2d import Mask2D + from autoarray import numba_util from autoarray.mask import mask_2d_util + from autoarray import type as ty +def over_sample_size_convert_to_array_2d_from( + over_sample_size: Union[int, np.ndarray], mask: Union[np.ndarray, Mask2D] +): + """ + Returns the over sample size as an `Array2D` object, for example converting it from a single integer. + + The interface allows a user to specify the `over_sample_size` as either: + + - A single integer, whereby over sampling is performed to this degree for every pixel. + - An ndarray with the same number of entries as the mask, to enable adaptive over sampling. + + This function converts these input structures to an `Array2D` which is used internally in the source code + to perform computations. + + Parameters + ---------- + over_sample_size + The over sampling scheme size, which divides the grid into a sub grid of smaller pixels when computing + values (e.g. images) from the grid to approximate the 2D line integral of the amount of light that falls + into each pixel. + + Returns + ------- + + """ + if isinstance(over_sample_size, int): + over_sample_size = np.full( + fill_value=over_sample_size, shape=mask.pixels_in_mask + ).astype("int") + + return Array2D(values=over_sample_size, mask=mask) + + @numba_util.jit() def total_sub_pixels_2d_from(sub_size: np.ndarray) -> int: """ @@ -342,7 +384,7 @@ def grid_2d_slim_over_sampled_via_mask_from( """ For a sub-grid, every unmasked pixel of its 2D mask with shape (total_y_pixels, total_x_pixels) is divided into a finer uniform grid of shape (total_y_pixels*sub_size, total_x_pixels*sub_size). This routine computes the (y,x) - scaled coordinates a the centre of every sub-pixel defined by this 2D mask array. + scaled coordinates at the centre of every sub-pixel defined by this 2D mask array. The sub-grid is returned on an array of shape (total_unmasked_pixels*sub_size**2, 2). y coordinates are stored in the 0 index of the second dimension, x coordinates in the 1 index. Masked coordinates are therefore @@ -517,3 +559,114 @@ def binned_array_2d_from( index += 1 return binned_array_2d_slim + + +def over_sample_size_via_radial_bins_from( + grid: Grid2D, + sub_size_list: List[int], + radial_list: List[float], + centre_list: List[Tuple] = None, +) -> Array2D: + """ + Returns an adaptive sub-grid size based on the radial distance of every pixel from the centre of the mask. + + The adaptive sub-grid size is computed as follows: + + 1) Compute the radial distance of every pixel in the mask from the centre of the mask (or input centres). + 2) For every pixel, determine the sub-grid size based on the radial distance of that pixel. For example, if + the first entry in `radial_list` is 0.5 and the first entry in `sub_size_list` 8, all pixels with a radial + distance less than 0.5 will have a sub-grid size of 8x8. + + This scheme can produce high sub-size values towards the centre of the mask, where the galaxy is brightest and + has the most rapidly changing light profile which requires a high sub-grid size to resolve accurately. + + If the data has multiple galaxies, the `centre_list` can be used to define the centre of each galaxy + and therefore increase the sub-grid size based on the light profile of each individual galaxy. + + Parameters + ---------- + mask + The mask defining the 2D region where the over-sampled grid is computed. + sub_size_list + The sub-grid size for every radial bin. + radial_list + The radial distance defining each bin, which are refeneced based on the previous entry. For example, if + the first entry is 0.5, the second 1.0 and the third 1.5, the adaptive sub-grid size will be between 0.5 + and 1.0 for the first sub-grid size, between 1.0 and 1.5 for the second sub-grid size, etc. + centre_list + A list of centres for each galaxy whose centres require higher sub-grid sizes. + + Returns + ------- + A uniform over-sampling object with an adaptive sub-grid size based on the radial distance of every pixel from + the centre of the mask. + """ + + if centre_list is None: + centre_list = [grid.mask.mask_centre] + + sub_size = np.zeros(grid.shape_slim) + + for centre in centre_list: + radial_grid = grid.distances_to_coordinate_from(coordinate=centre) + + sub_size_of_centre = sub_size_radial_bins_from( + radial_grid=np.array(radial_grid), + sub_size_list=np.array(sub_size_list), + radial_list=np.array(radial_list), + ) + + sub_size = np.where(sub_size_of_centre > sub_size, sub_size_of_centre, sub_size) + + return Array2D(values=sub_size, mask=grid.mask) + + +def over_sample_size_via_adapt_from( + data: Array2D, + noise_map: Array2D, + signal_to_noise_cut: float = 5.0, + sub_size_lower: int = 2, + sub_size_upper: int = 4, +) -> Array2D: + """ + Returns an adaptive sub-grid size based on the signal-to-noise of the data. + + The adaptive sub-grid size is computed as follows: + + 1) The signal-to-noise of every pixel is computed as the data divided by the noise-map. + 2) For all pixels with signal-to-noise above the signal-to-noise cut, the sub-grid size is set to the upper + value. For all other pixels, the sub-grid size is set to the lower value. + + This scheme can produce low sub-size values over entire datasets if the data has a low signal-to-noise. However, + just because the data has a low signal-to-noise does not mean that the sub-grid size should be low. + + To mitigate this, the signal-to-noise cut is set to the maximum signal-to-noise of the data divided by 2.0 if + it this value is below the signal-to-noise cut. + + Parameters + ---------- + data + The data which is to be fitted via a calculation using this over-sampling sub-grid. + noise_map + The noise-map of the data. + signal_to_noise_cut + The signal-to-noise cut which defines whether the sub-grid size is the upper or lower value. + sub_size_lower + The sub-grid size for pixels with signal-to-noise below the signal-to-noise cut. + sub_size_upper + The sub-grid size for pixels with signal-to-noise above the signal-to-noise cut. + + Returns + ------- + The adaptive sub-grid sizes. + """ + signal_to_noise = data / noise_map + + if np.max(signal_to_noise) < (2.0 * signal_to_noise_cut): + signal_to_noise_cut = np.max(signal_to_noise) / 2.0 + + sub_size = np.where( + signal_to_noise > signal_to_noise_cut, sub_size_upper, sub_size_lower + ) + + return Array2D(values=sub_size, mask=data.mask) diff --git a/autoarray/operators/over_sampling/uniform.py b/autoarray/operators/over_sampling/over_sampler.py similarity index 56% rename from autoarray/operators/over_sampling/uniform.py rename to autoarray/operators/over_sampling/over_sampler.py index b6c8d71a3..d7187a33d 100644 --- a/autoarray/operators/over_sampling/uniform.py +++ b/autoarray/operators/over_sampling/over_sampler.py @@ -5,23 +5,17 @@ from autoconf import cached_property from autoarray.mask.mask_2d import Mask2D -from autoarray.operators.over_sampling.abstract import AbstractOverSampling -from autoarray.operators.over_sampling.abstract import AbstractOverSampler from autoarray.structures.arrays.uniform_2d import Array2D -from autoarray.structures.grids.irregular_2d import Grid2DIrregular -from autoarray.structures.grids.uniform_2d import Grid2D -from autoarray import exc from autoarray.operators.over_sampling import over_sample_util from autofit.jax_wrapper import register_pytree_node_class - @register_pytree_node_class -class OverSamplingUniform(AbstractOverSampling): - def __init__(self, sub_size: Union[int, Array2D]): +class OverSampler: + def __init__(self, mask: Mask2D, sub_size: Union[int, Array2D]): """ - Over samples grid calculations using a uniform sub-grid. + Over samples grid calculations using a uniform sub-grid. When a 2D grid of (y,x) coordinates is input into a function, the result is evaluated at every coordinate on the grid. When the grid is paired to a 2D image (e.g. an `Array2D`) the solution needs to approximate @@ -134,231 +128,29 @@ def __init__(self, sub_size: Union[int, Array2D]): The over sampling class has functions dedicated to mapping between the sub-grid and pixel-grid, for example `sub_mask_native_for_sub_mask_slim` and `slim_for_sub_slim`. - Parameters - ---------- - sub_size - The size (sub_size x sub_size) of each unmasked pixels sub-grid. - """ - - self.sub_size = sub_size - - @classmethod - def from_radial_bins( - cls, - grid: Grid2D, - sub_size_list: List[int], - radial_list: List[float], - centre_list: List[Tuple] = None, - ) -> "OverSamplingUniform": - """ - Returns an adaptive sub-grid size based on the radial distance of every pixel from the centre of the mask. - - The adaptive sub-grid size is computed as follows: - - 1) Compute the radial distance of every pixel in the mask from the centre of the mask (or input centres). - 2) For every pixel, determine the sub-grid size based on the radial distance of that pixel. For example, if - the first entry in `radial_list` is 0.5 and the first entry in `sub_size_list` 8, all pixels with a radial - distance less than 0.5 will have a sub-grid size of 8x8. - - This scheme can produce high sub-size values towards the centre of the mask, where the galaxy is brightest and - has the most rapidly changing light profile which requires a high sub-grid size to resolve accurately. - - If the data has multiple galaxies, the `centre_list` can be used to define the centre of each galaxy - and therefore increase the sub-grid size based on the light profile of each individual galaxy. + The class `OverSampling` is used for the high level API, whereby this is where users input their + preferred over-sampling configuration. This class, `OverSampler`, contains the functionality + which actually performs the over-sampling calculations, but is hidden from the user. Parameters ---------- mask The mask defining the 2D region where the over-sampled grid is computed. - sub_size_list - The sub-grid size for every radial bin. - radial_list - The radial distance defining each bin, which are refeneced based on the previous entry. For example, if - the first entry is 0.5, the second 1.0 and the third 1.5, the adaptive sub-grid size will be between 0.5 - and 1.0 for the first sub-grid size, between 1.0 and 1.5 for the second sub-grid size, etc. - centre_list - A list of centres for each galaxy whose centres require higher sub-grid sizes. - - Returns - ------- - A uniform over-sampling object with an adaptive sub-grid size based on the radial distance of every pixel from - the centre of the mask. - """ - - if centre_list is None: - centre_list = [grid.mask.mask_centre] - - sub_size = np.zeros(grid.shape_slim) - - for centre in centre_list: - radial_grid = grid.distances_to_coordinate_from(coordinate=centre).array - - sub_size_of_centre = over_sample_util.sub_size_radial_bins_from( - radial_grid=np.array(radial_grid), - sub_size_list=np.array(sub_size_list), - radial_list=np.array(radial_list), - ) - - sub_size = np.where( - sub_size_of_centre > sub_size, sub_size_of_centre, sub_size - ) - - sub_size = Array2D(values=sub_size, mask=grid.mask) - - return cls(sub_size=sub_size) - - @classmethod - def from_adaptive_scheme( - cls, grid: Grid2D, name: str, centre: Tuple[float, float] - ) -> "OverSamplingUniform": - """ - Returns a 2D grid whose over sampling is adaptive, placing a high number of sub-pixels in the regions of the - grid closest to the centre input (y,x) coordinates. - - This adaptive over sampling is primarily used in PyAutoGalaxy, to evaluate the image of a light profile - (e.g. a Sersic function) with high levels of sub gridding in its centre and lower levels of sub gridding - further away from the centre (saving computational time). - - The `autogalaxy_workspace` and `autolens_workspace` packages have guides called `over_sampling.ipynb` - which describe over sampling in more detail. - - The inputs `name` and `centre` typically correspond to the name of the light profile (e.g. `Sersic`) and - its `centre`, so that the adaptive over sampling parameters for that light profile are loaded from config - files and used to setup the grid. - - Parameters - ---------- - name - The name of the light profile, which is used to load the adaptive over sampling parameters from config files. - centre - The (y,x) centre of the adaptive over sampled grid, around which the highest sub-pixel resolution is placed. - - Returns - ------- - A new Grid2D with adaptive over sampling. - - """ - - if not grid.is_uniform: - raise exc.GridException( - "You cannot make an adaptive over-sampled grid from a non-uniform grid." - ) - - sub_size_list = conf.instance["grids"]["over_sampling"]["sub_size_list"][name] - radial_factor_list = conf.instance["grids"]["over_sampling"][ - "radial_factor_list" - ][name] - - centre = grid.geometry.scaled_coordinate_2d_to_scaled_at_pixel_centre_from( - scaled_coordinate_2d=centre - ) - - return OverSamplingUniform.from_radial_bins( - grid=grid, - sub_size_list=sub_size_list, - radial_list=[ - min(grid.pixel_scales) * radial_factor - for radial_factor in radial_factor_list - ], - centre_list=[centre], - ) - - @classmethod - def from_adapt( - cls, - data: Array2D, - noise_map: Array2D, - signal_to_noise_cut: float = 5.0, - sub_size_lower: int = 2, - sub_size_upper: int = 4, - ): - """ - Returns an adaptive sub-grid size based on the signal-to-noise of the data. - - The adaptive sub-grid size is computed as follows: - - 1) The signal-to-noise of every pixel is computed as the data divided by the noise-map. - 2) For all pixels with signal-to-noise above the signal-to-noise cut, the sub-grid size is set to the upper - value. For all other pixels, the sub-grid size is set to the lower value. - - This scheme can produce low sub-size values over entire datasets if the data has a low signal-to-noise. However, - just because the data has a low signal-to-noise does not mean that the sub-grid size should be low. - - To mitigate this, the signal-to-noise cut is set to the maximum signal-to-noise of the data divided by 2.0 if - it this value is below the signal-to-noise cut. - - Parameters - ---------- - data - The data which is to be fitted via a calculation using this over-sampling sub-grid. - noise_map - The noise-map of the data. - signal_to_noise_cut - The signal-to-noise cut which defines whether the sub-grid size is the upper or lower value. - sub_size_lower - The sub-grid size for pixels with signal-to-noise below the signal-to-noise cut. - sub_size_upper - The sub-grid size for pixels with signal-to-noise above the signal-to-noise cut. - - Returns - ------- - The adaptive sub-grid sizes. + sub_size + The size (sub_size x sub_size) of each unmasked pixels sub-grid. """ - signal_to_noise = data / noise_map - - if np.max(signal_to_noise) < (2.0 * signal_to_noise_cut): - signal_to_noise_cut = np.max(signal_to_noise) / 2.0 + self.mask = mask - sub_size = np.where( - signal_to_noise > signal_to_noise_cut, sub_size_upper, sub_size_lower + self.sub_size = over_sample_util.over_sample_size_convert_to_array_2d_from( + over_sample_size=sub_size, mask=mask ) - sub_size = Array2D(values=sub_size, mask=data.mask) - - return cls(sub_size=sub_size) - - def over_sampler_from(self, mask: Mask2D) -> "OverSamplerUniform": - return OverSamplerUniform( - mask=mask, - sub_size=self.sub_size, - ) - def tree_flatten(self): - children = (self.sub_size,) - aux_data = None - return (children, aux_data) + return (self.mask,), () @classmethod def tree_unflatten(cls, aux_data, children): - return cls(*children) - - -class OverSamplerUniform(AbstractOverSampler): - def __init__(self, mask: Mask2D, sub_size: Union[int, Array2D]): - """ - Over samples grid calculations using a uniform sub-grid. - - See the class `OverSamplingUniform` for a description of how the over-sampling works in full. - - The class `OverSamplingUniform` is used for the high level API, whereby this is where users input their - preferred over-sampling configuration. This class, `OverSamplerUniform`, contains the functionality - which actually performs the over-sampling calculations, but is hidden from the user. - - Parameters - ---------- - mask - The mask defining the 2D region where the over-sampled grid is computed. - sub_size - The size (sub_size x sub_size) of each unmasked pixels sub-grid. - """ - self.mask = mask - - if isinstance(sub_size, int): - sub_size = Array2D( - values=np.full(fill_value=sub_size, shape=mask.shape_slim), mask=mask - ) - - self.sub_size = sub_size + return cls(mask=children[0]) @property def sub_total(self): @@ -404,17 +196,6 @@ def sub_pixel_areas(self) -> np.ndarray: return sub_pixel_areas - @cached_property - def over_sampled_grid(self) -> Grid2DIrregular: - grid = over_sample_util.grid_2d_slim_over_sampled_via_mask_from( - mask_2d=np.array(self.mask), - pixel_scales=self.mask.pixel_scales, - sub_size=np.array(self.sub_size.array).astype("int"), - origin=self.mask.origin, - ) - - return Grid2DIrregular(values=grid) - def binned_array_2d_from(self, array: Array2D) -> "Array2D": """ Convenience method to access the binned-up array in its 1D representation, which is a Grid2D stored as an @@ -439,9 +220,9 @@ def binned_array_2d_from(self, array: Array2D) -> "Array2D": pass binned_array_2d = over_sample_util.binned_array_2d_from( - array_2d=np.array(array.array), + array_2d=np.array(array), mask_2d=np.array(self.mask), - sub_size=np.array(self.sub_size.array).astype("int"), + sub_size=np.array(self.sub_size).astype("int"), ) return Array2D( @@ -449,16 +230,6 @@ def binned_array_2d_from(self, array: Array2D) -> "Array2D": mask=self.mask, ) - def array_via_func_from(self, func, obj, *args, **kwargs): - over_sampled_grid = self.over_sampled_grid - - if obj is not None: - values = func(obj, over_sampled_grid, *args, **kwargs) - else: - values = func(over_sampled_grid, *args, **kwargs) - - return self.binned_array_2d_from(array=values) - @cached_property def sub_mask_native_for_sub_mask_slim(self) -> np.ndarray: """ @@ -569,3 +340,29 @@ def slim_for_sub_slim(self) -> np.ndarray: return over_sample_util.slim_index_for_sub_slim_index_via_mask_2d_from( mask_2d=np.array(self.mask), sub_size=np.array(self.sub_size) ).astype("int") + + @property + def uniform_over_sampled(self): + """ + For a sub-grid, every unmasked pixel of its 2D mask with shape (total_y_pixels, total_x_pixels) is divided into + a finer uniform grid of shape (total_y_pixels*sub_size, total_x_pixels*sub_size). This routine computes the (y,x) + scaled coordinates at the centre of every sub-pixel defined by this 2D mask array. + + The sub-grid is returned on an array of shape (total_unmasked_pixels*sub_size**2, 2). y coordinates are + stored in the 0 index of the second dimension, x coordinates in the 1 index. Masked coordinates are therefore + removed and not included in the slimmed grid. + + Grid2D are defined from the top-left corner, where the first unmasked sub-pixel corresponds to index 0. + Sub-pixels that are part of the same mask array pixel are indexed next to one another, such that the second + sub-pixel in the first pixel has index 1, its next sub-pixel has index 2, and so forth. + """ + from autoarray.structures.grids.irregular_2d import Grid2DIrregular + + grid = over_sample_util.grid_2d_slim_over_sampled_via_mask_from( + mask_2d=np.array(self.mask), + pixel_scales=self.mask.pixel_scales, + sub_size=np.array(self.sub_size).astype("int"), + origin=self.mask.origin, + ) + + return Grid2DIrregular(values=grid) diff --git a/autoarray/operators/transformer.py b/autoarray/operators/transformer.py index 89e48ecf2..acbe6bb73 100644 --- a/autoarray/operators/transformer.py +++ b/autoarray/operators/transformer.py @@ -47,7 +47,7 @@ def pylops_exception(): "\n--------------------\n" "You are attempting to perform interferometer analysis.\n\n" "However, the optional library PyLops (https://github.com/PyLops/pylops) is not installed.\n\n" - "Install it via the command `pip install pylops==1.18.3`.\n\n" + "Install it via the command `pip install pylops==2.3.1`.\n\n" "----------------------" ) diff --git a/autoarray/plot/__init__.py b/autoarray/plot/__init__.py index c8553060a..71d7abb45 100644 --- a/autoarray/plot/__init__.py +++ b/autoarray/plot/__init__.py @@ -30,6 +30,7 @@ from autoarray.plot.wrap.two_d.interpolated_reconstruction import ( InterpolatedReconstruction, ) +from autoarray.plot.wrap.two_d.delaunay_drawer import DelaunayDrawer from autoarray.plot.wrap.two_d.voronoi_drawer import VoronoiDrawer from autoarray.plot.wrap.two_d.origin_scatter import OriginScatter from autoarray.plot.wrap.two_d.mask_scatter import MaskScatter diff --git a/autoarray/plot/abstract_plotters.py b/autoarray/plot/abstract_plotters.py index d179533d9..07db07291 100644 --- a/autoarray/plot/abstract_plotters.py +++ b/autoarray/plot/abstract_plotters.py @@ -37,7 +37,6 @@ def __init__( self.subplot_figsize = None - def set_title(self, label): if self.mat_plot_1d is not None: self.mat_plot_1d.title.manual_label = label diff --git a/autoarray/plot/get_visuals/two_d.py b/autoarray/plot/get_visuals/two_d.py index 816a9ed89..c2b99a173 100644 --- a/autoarray/plot/get_visuals/two_d.py +++ b/autoarray/plot/get_visuals/two_d.py @@ -183,11 +183,13 @@ def via_mapper_for_source_from(self, mapper: MapperRectangular) -> Visuals2D: ) grid = self.get( - "grid", mapper.source_plane_data_grid, "mapper_source_plane_data_grid" + "grid", + mapper.source_plane_data_grid.over_sampled, + "mapper_source_plane_data_grid", ) try: - border_grid = mapper.mapper_grids.source_plane_data_grid[ + border_grid = mapper.mapper_grids.source_plane_data_grid.over_sampled[ mapper.border_relocator.sub_border_slim ] border = self.get("border", border_grid) diff --git a/autoarray/plot/mat_plot/two_d.py b/autoarray/plot/mat_plot/two_d.py index b527fd7ad..9fcec3657 100644 --- a/autoarray/plot/mat_plot/two_d.py +++ b/autoarray/plot/mat_plot/two_d.py @@ -48,6 +48,7 @@ def __init__( vector_yx_quiver: Optional[w2d.VectorYXQuiver] = None, patch_overlay: Optional[w2d.PatchOverlay] = None, interpolated_reconstruction: Optional[w2d.InterpolatedReconstruction] = None, + delaunay_drawer: Optional[w2d.DelaunayDrawer] = None, voronoi_drawer: Optional[w2d.VoronoiDrawer] = None, origin_scatter: Optional[w2d.OriginScatter] = None, mask_scatter: Optional[w2d.MaskScatter] = None, @@ -128,6 +129,8 @@ def __init__( Plots a `VectorField` object using the matplotlib function `plt.quiver`. patch_overlay Overlays matplotlib `patches.Patch` objects over the figure, such as an `Ellipse`. + delaunay_drawer + Draws a colored Delaunay mesh of pixels using `plt.tripcolor`. voronoi_drawer Draws a colored Voronoi mesh of pixels using `plt.fill`. interpolated_reconstruction @@ -188,6 +191,7 @@ def __init__( interpolated_reconstruction or w2d.InterpolatedReconstruction(is_default=True) ) + self.delaunay_drawer = delaunay_drawer or w2d.DelaunayDrawer(is_default=True) self.voronoi_drawer = voronoi_drawer or w2d.VoronoiDrawer(is_default=True) self.origin_scatter = origin_scatter or w2d.OriginScatter(is_default=True) @@ -417,7 +421,7 @@ def plot_grid( ax = self.setup_subplot() if plot_over_sampled_grid: - grid_plot = grid.over_sampler.over_sampled_grid + grid_plot = grid.over_sampled else: grid_plot = grid @@ -606,7 +610,7 @@ def _plot_rectangular_mapper( visuals_2d.plot_via_plotter( plotter=self, - grid_indexes=mapper.source_plane_data_grid, + grid_indexes=mapper.source_plane_data_grid.over_sampled, mapper=mapper, geometry=mapper.mapper_grids.mask.geometry, ) @@ -639,6 +643,8 @@ def _plot_delaunay_mapper( self.axis.set(extent=extent, grid=mapper.source_plane_mesh_grid) + plt.gca().set_aspect(aspect_inv) + self.tickparams.set() self.yticks.set(min_value=extent[2], max_value=extent[3], units=self.units) self.xticks.set(min_value=extent[0], max_value=extent[1], units=self.units) @@ -653,17 +659,34 @@ def _plot_delaunay_mapper( else: [annotate.set() for annotate in self.annotate] - interpolation_array = self.interpolated_reconstruction.imshow_reconstruction( - mapper=mapper, - pixel_values=pixel_values, - units=self.units, - cmap=self.cmap, - colorbar=self.colorbar, - colorbar_tickparams=self.colorbar_tickparams, - aspect=aspect_inv, - ax=ax, - use_log10=self.use_log10, - ) + interpolation_array = None + + if interpolate_to_uniform: + interpolation_array = ( + self.interpolated_reconstruction.imshow_reconstruction( + mapper=mapper, + pixel_values=pixel_values, + units=self.units, + cmap=self.cmap, + colorbar=self.colorbar, + colorbar_tickparams=self.colorbar_tickparams, + aspect=aspect_inv, + ax=ax, + use_log10=self.use_log10, + ) + ) + + else: + self.delaunay_drawer.draw_delaunay_pixels( + mapper=mapper, + pixel_values=pixel_values, + units=self.units, + cmap=self.cmap, + colorbar=self.colorbar, + colorbar_tickparams=self.colorbar_tickparams, + ax=ax, + use_log10=self.use_log10, + ) self.title.set(auto_title=auto_labels.title) self.ylabel.set() @@ -671,7 +694,7 @@ def _plot_delaunay_mapper( visuals_2d.plot_via_plotter( plotter=self, - grid_indexes=mapper.source_plane_data_grid, + grid_indexes=mapper.source_plane_data_grid.over_sampled, mapper=mapper, geometry=mapper.mapper_grids.mask.geometry, ) @@ -754,7 +777,7 @@ def _plot_voronoi_mapper( visuals_2d.plot_via_plotter( plotter=self, - grid_indexes=mapper.source_plane_data_grid, + grid_indexes=mapper.source_plane_data_grid.over_sampled, mapper=mapper, geometry=mapper.mapper_grids.mask.geometry, ) diff --git a/autoarray/plot/multi_plotters.py b/autoarray/plot/multi_plotters.py index 01ed8364f..5c2c5071e 100644 --- a/autoarray/plot/multi_plotters.py +++ b/autoarray/plot/multi_plotters.py @@ -1,3 +1,5 @@ +import os +from pathlib import Path from typing import List, Optional, Tuple from autoarray.plot.wrap.base.ticks import YTicks @@ -104,29 +106,6 @@ def plot_via_func(self, plotter, figure_name: str, func_name: str, kwargs): else: func(**{**{figure_name: True}, **kwargs}) - def output_subplot(self, filename_suffix: str = ""): - """ - Outplot the subplot to a file after all figures have been plotted on the subplot. - - The multi-plotter requires its own output function to ensure that the subplot is output to a file, which - this provides. - - Parameters - ---------- - filename_suffix - The suffix of the filename that the subplot is output to. - """ - - if self.plotter_list[0].mat_plot_1d is not None: - self.plotter_list[0].mat_plot_1d.output.subplot_to_figure( - auto_filename=f"subplot_{filename_suffix}" - ) - if self.plotter_list[0].mat_plot_2d is not None: - self.plotter_list[0].mat_plot_2d.output.subplot_to_figure( - auto_filename=f"subplot_{filename_suffix}" - ) - self.plotter_list[0].close_subplot_figure() - def subplot_of_figure( self, func_name: str, figure_name: str, filename_suffix: str = "", **kwargs ): @@ -266,6 +245,99 @@ def subplot_of_multi_yx_1d(self, filename_suffix="", **kwargs): ) self.plotter_list[0].plotter_list[0].close_subplot_figure() + def output_subplot(self, filename_suffix: str = ""): + """ + Outplot the subplot to a file after all figures have been plotted on the subplot. + + The multi-plotter requires its own output function to ensure that the subplot is output to a file, which + this provides. + + Parameters + ---------- + filename_suffix + The suffix of the filename that the subplot is output to. + """ + + plotter = self.plotter_list[0] + + if plotter.mat_plot_1d is not None: + plotter.mat_plot_1d.output.subplot_to_figure( + auto_filename=f"subplot_{filename_suffix}" + ) + if plotter.mat_plot_2d is not None: + plotter.mat_plot_2d.output.subplot_to_figure( + auto_filename=f"subplot_{filename_suffix}" + ) + plotter.close_subplot_figure() + + def output_to_fits( + self, + func_name_list: List[str], + figure_name_list: List[str], + filename: str, + tag_list: Optional[List[str]] = None, + remove_fits_first: bool = False, + **kwargs, + ): + """ + Outputs a list of figures of the plotter objects in the `plotter_list` to a single .fits file. + + This function takes as input lists of function names and figure names and then calls them via + the `plotter_list` with an interface that outputs each to a .fits file. + + For example, if you have multiple `ImagingPlotter` objects and want to output the `data` and `noise_map` of + each to a single .fits files, you would input: + + - `func_name_list=['figures_2d', 'figures_2d']` and + - `figure_name_list=['data', 'noise_map']`. + + The implementation of this code is hacky, with it using a specific interface in the `Output` object + which sets the format to `fits_multi` to call a function which outputs the .fits files. A major visualuzation + refactor is required to make this more elegant. + + Parameters + ---------- + func_name_list + The list of function names that are called to plot the figures on the subplot. + figure_name_list + The list of figure names that are plotted on the subplot. + filenane + The filename that the .fits file is output to. + tag_list + The list of tags that are used to set the `EXTNAME` of each hdu of the .fits file. + remove_fits_first + If the .fits file already exists, it is removed before the new .fits file is output, else it is updated + with the figure going into the next hdu. + kwargs + Any additional keyword arguments that are passed to the function that plots the figure on the subplot. + """ + + output_path = self.plotter_list[0].mat_plot_2d.output.output_path_from( + format="fits_multi" + ) + output_fits_file = Path(output_path)/ f"{filename}.fits" + + if remove_fits_first: + output_fits_file.unlink(missing_ok=True) + + for i, plotter in enumerate(self.plotter_list): + plotter.mat_plot_2d.output._format = "fits_multi" + + plotter.set_filename(filename=f"{filename}") + + for j, (func_name, figure_name) in enumerate( + zip(func_name_list, figure_name_list) + ): + if tag_list is not None: + plotter.mat_plot_2d.output._tag_fits_multi = tag_list[j] + + self.plot_via_func( + plotter=plotter, + figure_name=figure_name, + func_name=func_name, + kwargs=kwargs, + ) + class MultiYX1DPlotter: def __init__( diff --git a/autoarray/plot/visuals/two_d.py b/autoarray/plot/visuals/two_d.py index 9a5629c91..d1da534ff 100644 --- a/autoarray/plot/visuals/two_d.py +++ b/autoarray/plot/visuals/two_d.py @@ -100,6 +100,6 @@ def plot_via_plotter(self, plotter, grid_indexes=None, mapper=None, geometry=Non else: plotter.index_scatter.scatter_grid_indexes( - grid=mapper.source_plane_data_grid, + grid=mapper.source_plane_data_grid.over_sampled, indexes=indexes, ) diff --git a/autoarray/plot/wrap/base/figure.py b/autoarray/plot/wrap/base/figure.py index fb25da2e4..238264fd9 100644 --- a/autoarray/plot/wrap/base/figure.py +++ b/autoarray/plot/wrap/base/figure.py @@ -1,4 +1,5 @@ from enum import Enum +import gc import matplotlib.pyplot as plt from typing import Union, Tuple @@ -92,5 +93,5 @@ def close(self): """ Wraps the Matplotlib method 'plt.close' for closing a figure. """ - if plt.fignum_exists(num=1): - plt.close() + plt.close() + gc.collect() diff --git a/autoarray/plot/wrap/base/output.py b/autoarray/plot/wrap/base/output.py index 6646e8c21..945f91abd 100644 --- a/autoarray/plot/wrap/base/output.py +++ b/autoarray/plot/wrap/base/output.py @@ -61,6 +61,7 @@ def __init__( self.format_folder = format_folder self.bypass = bypass self.bbox_inches = bbox_inches + self._tag_fits_multi = None self.kwargs = kwargs @@ -105,6 +106,7 @@ def savefig(self, filename: str, output_path: str, format: str): plt.savefig( path.join(output_path, f"{filename}.{format}"), bbox_inches=self.bbox_inches, + pad_inches=0, ) except ValueError as e: logger.info( @@ -151,6 +153,17 @@ def to_figure( file_path=path.join(output_path, f"{filename}.fits"), overwrite=True, ) + elif format == "fits_multi": + if structure is not None: + from autoarray.structures.arrays.array_2d_util import ( + update_fits_file, + ) + + update_fits_file( + arr=structure.native, + file_path=path.join(output_path, f"{filename}.fits"), + tag=self._tag_fits_multi, + ) def subplot_to_figure(self, auto_filename: Optional[str] = None): """ diff --git a/autoarray/plot/wrap/two_d/__init__.py b/autoarray/plot/wrap/two_d/__init__.py index 3ca56f6e3..5eb85eeab 100644 --- a/autoarray/plot/wrap/two_d/__init__.py +++ b/autoarray/plot/wrap/two_d/__init__.py @@ -8,6 +8,7 @@ from .interpolated_reconstruction import ( InterpolatedReconstruction, ) +from .delaunay_drawer import DelaunayDrawer from .voronoi_drawer import VoronoiDrawer from .origin_scatter import OriginScatter from .mask_scatter import MaskScatter diff --git a/autoarray/plot/wrap/two_d/delaunay_drawer.py b/autoarray/plot/wrap/two_d/delaunay_drawer.py new file mode 100644 index 000000000..5d168213b --- /dev/null +++ b/autoarray/plot/wrap/two_d/delaunay_drawer.py @@ -0,0 +1,116 @@ +import matplotlib.pyplot as plt +import numpy as np +from typing import Optional + +from autoarray.plot.wrap.two_d.abstract import AbstractMatWrap2D +from autoarray.plot.wrap.base.units import Units +from autoarray.inversion.pixelization.mappers.voronoi import MapperVoronoi + +from autoarray.plot.wrap import base as wb + + +def facecolors_from(values, simplices): + facecolors = np.zeros(shape=simplices.shape[0]) + for i in range(simplices.shape[0]): + facecolors[i] = np.sum(1.0 / 3.0 * values[simplices[i, :]]) + + return facecolors + + +class DelaunayDrawer(AbstractMatWrap2D): + """ + Draws Voronoi pixels from a `MapperVoronoi` object (see `inversions.mapper`). This includes both drawing + each Voronoi cell and coloring it according to a color value. + + The mapper contains the grid of (y,x) coordinate where the centre of each Voronoi cell is plotted. + + This object wraps methods described in below: + + https://matplotlib.org/3.3.2/api/_as_gen/matplotlib.pyplot.fill.html + """ + + def draw_delaunay_pixels( + self, + mapper: MapperVoronoi, + pixel_values: Optional[np.ndarray], + units: Units, + cmap: Optional[wb.Cmap], + colorbar: Optional[wb.Colorbar], + colorbar_tickparams: Optional[wb.ColorbarTickParams] = None, + ax=None, + use_log10: bool = False, + ): + """ + Draws the Voronoi pixels of the input `mapper` using its `mesh_grid` which contains the (y,x) + coordinate of the centre of every Voronoi cell. This uses the method `plt.fill`. + + Parameters + ---------- + mapper + A mapper object which contains the Voronoi mesh. + pixel_values + An array used to compute the color values that every Voronoi cell is plotted using. + cmap + The colormap used to plot each Voronoi cell. + colorbar + The `Colorbar` object in `mat_base` used to set the colorbar of the figure the Voronoi mesh is plotted on. + colorbar_tickparams + The `ColorbarTickParams` object in `mat_base` used to set the tick labels of the colorbar. + ax + The matplotlib axis the Voronoi mesh is plotted on. + use_log10 + If `True`, the colorbar is plotted using a log10 scale. + """ + + if pixel_values is None: + raise ValueError( + "pixel_values input to DelaunayPlotter are None and thus cannot be plotted." + ) + + if ax is None: + ax = plt.gca() + + source_pixelization_grid = mapper.mapper_grids.source_plane_mesh_grid + + simplices = mapper.delaunay.simplices + + facecolors = facecolors_from(values=pixel_values, simplices=simplices) + + norm = cmap.norm_from(array=pixel_values, use_log10=use_log10) + + if use_log10: + pixel_values[pixel_values < 1e-4] = 1e-4 + pixel_values = np.log10(pixel_values) + + vmin = cmap.vmin_from(array=pixel_values, use_log10=use_log10) + vmax = cmap.vmax_from(array=pixel_values, use_log10=use_log10) + + color_values = np.where(pixel_values > vmax, vmax, pixel_values) + color_values = np.where(pixel_values < vmin, vmin, color_values) + + cmap = plt.get_cmap(cmap.cmap) + + if colorbar is not None: + cb = colorbar.set_with_color_values( + units=units, + norm=norm, + cmap=cmap, + color_values=color_values, + ax=ax, + use_log10=use_log10, + ) + + if cb is not None and colorbar_tickparams is not None: + colorbar_tickparams.set(cb=cb) + + ax.tripcolor( + source_pixelization_grid[:, 1], + source_pixelization_grid[:, 0], + simplices, + facecolors=facecolors, + edgecolors="None", + cmap=cmap, + vmin=vmin, + vmax=vmax, + **self.config_dict + ) diff --git a/autoarray/plot/wrap/two_d/grid_plot.py b/autoarray/plot/wrap/two_d/grid_plot.py index 123f3194a..a727bec30 100644 --- a/autoarray/plot/wrap/two_d/grid_plot.py +++ b/autoarray/plot/wrap/two_d/grid_plot.py @@ -47,11 +47,15 @@ def plot_rectangular_grid_lines( ys = np.linspace(extent[2], extent[3], shape_native[1] + 1) xs = np.linspace(extent[0], extent[1], shape_native[0] + 1) + config_dict = self.config_dict + config_dict.pop("c") + config_dict["c"] = "k" + # grid lines for x in xs: - plt.plot([x, x], [ys[0], ys[-1]], **self.config_dict) + plt.plot([x, x], [ys[0], ys[-1]], **config_dict) for y in ys: - plt.plot([xs[0], xs[-1]], [y, y], **self.config_dict) + plt.plot([xs[0], xs[-1]], [y, y], **config_dict) def plot_grid(self, grid: Union[np.ndarray, Grid2D]): """ diff --git a/autoarray/plot/wrap/two_d/voronoi_drawer.py b/autoarray/plot/wrap/two_d/voronoi_drawer.py index 9e189413b..aec6661cc 100644 --- a/autoarray/plot/wrap/two_d/voronoi_drawer.py +++ b/autoarray/plot/wrap/two_d/voronoi_drawer.py @@ -81,7 +81,7 @@ def draw_voronoi_pixels( cmap = plt.get_cmap(cmap.cmap) - if colorbar is not None: + if colorbar is not None and colorbar is not False: cb = colorbar.set_with_color_values( units=units, norm=norm, diff --git a/autoarray/preloads.py b/autoarray/preloads.py new file mode 100644 index 000000000..6808f0f6c --- /dev/null +++ b/autoarray/preloads.py @@ -0,0 +1,568 @@ +import logging +import numpy as np +import os +from typing import List + +from autoconf import conf + +from autoarray.inversion.linear_obj.func_list import AbstractLinearObjFuncList +from autoarray.inversion.pixelization.mappers.abstract import AbstractMapper + +from autoarray import exc +from autoarray.inversion.inversion.imaging import inversion_imaging_util + +logger = logging.getLogger(__name__) + +logger.setLevel(level="INFO") + + +class Preloads: + def __init__( + self, + w_tilde=None, + use_w_tilde=None, + image_plane_mesh_grid_pg_list=None, + relocated_grid=None, + mapper_list=None, + operated_mapping_matrix=None, + linear_func_operated_mapping_matrix_dict=None, + data_linear_func_matrix_dict=None, + mapper_operated_mapping_matrix_dict=None, + curvature_matrix=None, + data_vector_mapper=None, + curvature_matrix_mapper_diag=None, + regularization_matrix=None, + log_det_regularization_matrix_term=None, + traced_mesh_grids_list_of_planes=None, + image_plane_mesh_grid_list=None, + ): + self.w_tilde = w_tilde + self.use_w_tilde = use_w_tilde + + self.image_plane_mesh_grid_pg_list = image_plane_mesh_grid_pg_list + self.relocated_grid = relocated_grid + self.mapper_list = mapper_list + self.operated_mapping_matrix = operated_mapping_matrix + self.linear_func_operated_mapping_matrix_dict = ( + linear_func_operated_mapping_matrix_dict + ) + self.data_linear_func_matrix_dict = data_linear_func_matrix_dict + self.mapper_operated_mapping_matrix_dict = mapper_operated_mapping_matrix_dict + self.curvature_matrix = curvature_matrix + self.data_vector_mapper = data_vector_mapper + self.curvature_matrix_mapper_diag = curvature_matrix_mapper_diag + self.regularization_matrix = regularization_matrix + self.log_det_regularization_matrix_term = log_det_regularization_matrix_term + + self.traced_mesh_grids_list_of_planes = traced_mesh_grids_list_of_planes + self.image_plane_mesh_grid_list = image_plane_mesh_grid_list + + @property + def check_threshold(self): + return conf.instance["general"]["test"]["preloads_check_threshold"] + + def set_w_tilde_imaging(self, fit_0, fit_1): + """ + The w-tilde linear algebra formalism speeds up inversions by computing beforehand quantities that enable + efficiently construction of the curvature matrix. These quantities can only be used if the noise-map is + fixed, therefore this function preloads these w-tilde quantities if the noise-map does not change. + + This function compares the noise map of two fit's corresponding to two model instances, and preloads wtilde + if the noise maps of both fits are the same. + + The preload is typically used through search chaining pipelines, as it is uncommon for the noise map to be + scaled during the model-fit (although it is common for a fixed but scaled noise map to be used). + + Parameters + ---------- + fit_0 + The first fit corresponding to a model with a specific set of unit-values. + fit_1 + The second fit corresponding to a model with a different set of unit-values. + """ + + self.w_tilde = None + self.use_w_tilde = False + + if fit_0.inversion is None: + return + + if not fit_0.inversion.has(cls=AbstractMapper): + return + + if np.max(abs(fit_0.noise_map - fit_1.noise_map)) < 1e-8: + logger.info("PRELOADS - Computing W-Tilde... May take a moment.") + + from autoarray.dataset.imaging.w_tilde import WTildeImaging + + ( + preload, + indexes, + lengths, + ) = inversion_imaging_util.w_tilde_curvature_preload_imaging_from( + noise_map_native=np.array(fit_0.noise_map.native), + kernel_native=np.array(fit_0.dataset.psf.native), + native_index_for_slim_index=np.array( + fit_0.dataset.mask.derive_indexes.native_for_slim + ), + ) + + self.w_tilde = WTildeImaging( + curvature_preload=preload, + indexes=indexes.astype("int"), + lengths=lengths.astype("int"), + noise_map_value=fit_0.noise_map[0], + ) + + self.use_w_tilde = True + + logger.info("PRELOADS - W-Tilde preloaded for this model-fit.") + + def set_relocated_grid(self, fit_0, fit_1): + """ + If the `MassProfile`'s in a model are fixed their traced grids (which may have had coordinates relocated at + the border) does not change during the model=fit and can therefore be preloaded. + + This function compares the relocated grids of the mappers of two fit corresponding to two model instances, and + preloads the grid if the grids of both fits are the same. + + The preload is typically used in adapt searches, where the mass model is fixed and the parameters are + varied. + + Parameters + ---------- + fit_0 + The first fit corresponding to a model with a specific set of unit-values. + fit_1 + The second fit corresponding to a model with a different set of unit-values. + """ + + self.relocated_grid = None + + if fit_0.inversion is None: + return + + if ( + fit_0.inversion.total(cls=AbstractMapper) > 1 + or fit_0.inversion.total(cls=AbstractMapper) == 0 + ): + return + + mapper_0 = fit_0.inversion.cls_list_from(cls=AbstractMapper)[0] + mapper_1 = fit_1.inversion.cls_list_from(cls=AbstractMapper)[0] + + if ( + mapper_0.source_plane_data_grid.shape[0] + == mapper_1.source_plane_data_grid.shape[0] + ): + if ( + np.max( + abs( + mapper_0.source_plane_data_grid + - mapper_1.source_plane_data_grid + ) + ) + < 1.0e-8 + ): + self.relocated_grid = mapper_0.source_plane_data_grid + + logger.info( + "PRELOADS - Relocated grid of pxielization preloaded for this model-fit." + ) + + def set_mapper_list(self, fit_0, fit_1): + """ + If the `MassProfile`'s and `Mesh`'s in a model are fixed, the mapping of image-pixels to the + source-pixels does not change during the model-fit and the list of `Mapper`'s containing this information can + be preloaded. This includes preloading the `mapping_matrix`. + + This function compares the mapping matrix of two fit's corresponding to two model instances, and preloads the + list of mappers if the mapping matrix of both fits are the same. + + The preload is typically used in searches where only light profiles vary (e.g. when only the lens's light is + being fitted for). + + Parameters + ---------- + fit_0 + The first fit corresponding to a model with a specific set of unit-values. + fit_1 + The second fit corresponding to a model with a different set of unit-values. + """ + + self.mapper_list = None + + if fit_0.inversion is None: + return + + if fit_0.inversion.total(cls=AbstractMapper) == 0: + return + + from autoarray.inversion.inversion.interferometer.lop import ( + InversionInterferometerMappingPyLops, + ) + + if isinstance(fit_0.inversion, InversionInterferometerMappingPyLops): + return + + inversion_0 = fit_0.inversion + inversion_1 = fit_1.inversion + + if inversion_0.mapping_matrix.shape[1] == inversion_1.mapping_matrix.shape[1]: + if np.allclose(inversion_0.mapping_matrix, inversion_1.mapping_matrix): + self.mapper_list = inversion_0.cls_list_from(cls=AbstractMapper) + + logger.info( + "PRELOADS - Mappers of planes preloaded for this model-fit." + ) + + def set_operated_mapping_matrix_with_preloads(self, fit_0, fit_1): + """ + If the `MassProfile`'s and `Mesh`'s in a model are fixed, the mapping of image-pixels to the + source-pixels does not change during the model-fit and matrices used to perform the linear algebra in an + inversion can be preloaded, which help efficiently construct the curvature matrix. + + This function compares the operated mapping matrix of two fit's corresponding to two model instances, and + preloads the mapper if the mapping matrix of both fits are the same. + + The preload is typically used in searches where only light profiles vary (e.g. when only the lens's light is + being fitted for). + + Parameters + ---------- + fit_0 + The first fit corresponding to a model with a specific set of unit-values. + fit_1 + The second fit corresponding to a model with a different set of unit-values. + """ + + self.operated_mapping_matrix = None + + from autoarray.inversion.inversion.interferometer.lop import ( + InversionInterferometerMappingPyLops, + ) + + if isinstance(fit_0.inversion, InversionInterferometerMappingPyLops): + return + + inversion_0 = fit_0.inversion + inversion_1 = fit_1.inversion + + if inversion_0 is None: + return + + if ( + inversion_0.operated_mapping_matrix.shape[1] + == inversion_1.operated_mapping_matrix.shape[1] + ): + if ( + np.max( + abs( + inversion_0.operated_mapping_matrix + - inversion_1.operated_mapping_matrix + ) + ) + < 1e-8 + ): + self.operated_mapping_matrix = inversion_0.operated_mapping_matrix + + logger.info( + "PRELOADS - Inversion linear algebra quantities preloaded for this model-fit." + ) + + def set_linear_func_inversion_dicts(self, fit_0, fit_1): + """ + If the `MassProfile`'s and `Mesh`'s in a model are fixed, the mapping of image-pixels to the + source-pixels does not change during the model-fit and matrices used to perform the linear algebra in an + inversion can be preloaded, which help efficiently construct the curvature matrix. + + This function compares the operated mapping matrix of two fit's corresponding to two model instances, and + preloads the mapper if the mapping matrix of both fits are the same. + + The preload is typically used in searches where only light profiles vary (e.g. when only the lens's light is + being fitted for). + + Parameters + ---------- + fit_0 + The first fit corresponding to a model with a specific set of unit-values. + fit_1 + The second fit corresponding to a model with a different set of unit-values. + """ + + from autoarray.inversion.pixelization.pixelization import Pixelization + + self.linear_func_operated_mapping_matrix_dict = None + + inversion_0 = fit_0.inversion + inversion_1 = fit_1.inversion + + if inversion_0 is None: + return + + if not inversion_0.has(cls=AbstractMapper): + return + + if not inversion_0.has(cls=AbstractLinearObjFuncList): + return + + try: + inversion_0.linear_func_operated_mapping_matrix_dict + except NotImplementedError: + return + + if not hasattr(inversion_0, "linear_func_operated_mapping_matrix_dict"): + return + + should_preload = False + + for operated_mapping_matrix_0, operated_mapping_matrix_1 in zip( + inversion_0.linear_func_operated_mapping_matrix_dict.values(), + inversion_1.linear_func_operated_mapping_matrix_dict.values(), + ): + if ( + np.max(abs(operated_mapping_matrix_0 - operated_mapping_matrix_1)) + < 1e-8 + ): + should_preload = True + else: + should_preload = False + break + + if should_preload: + self.linear_func_operated_mapping_matrix_dict = ( + inversion_0.linear_func_operated_mapping_matrix_dict + ) + self.data_linear_func_matrix_dict = inversion_0.data_linear_func_matrix_dict + + logger.info( + "PRELOADS - Inversion linear light profile operated mapping matrix / data linear func matrix preloaded for this model-fit." + ) + + def set_curvature_matrix(self, fit_0, fit_1): + """ + If the `MassProfile`'s and `Mesh`'s in a model are fixed, the mapping of image-pixels to the + source-pixels does not change during the model-fit and therefore its associated curvature matrix is also + fixed, meaning the curvature matrix preloaded. + + If linear ``LightProfiles``'s are included, the regions of the curvature matrix associatd with these + objects vary but the diagonals of the mapper do not change. In this case, the `curvature_matrix_mapper_diag` + is preloaded. + + This function compares the curvature matrix of two fit's corresponding to two model instances, and preloads + this value if it is the same for both fits. + + The preload is typically used in **PyAutoGalaxy** inversions using a `Rectangular` pixelization. + + Parameters + ---------- + fit_0 + The first fit corresponding to a model with a specific set of unit-values. + fit_1 + The second fit corresponding to a model with a different set of unit-values. + """ + + self.curvature_matrix = None + self.data_vector_mapper = None + self.curvature_matrix_mapper_diag = None + self.mapper_operated_mapping_matrix_dict = None + + inversion_0 = fit_0.inversion + inversion_1 = fit_1.inversion + + if inversion_0 is None: + return + + try: + inversion_0._curvature_matrix_mapper_diag + except NotImplementedError: + return + + if inversion_0.curvature_matrix.shape == inversion_1.curvature_matrix.shape: + if ( + np.max(abs(inversion_0.curvature_matrix - inversion_1.curvature_matrix)) + < 1e-8 + ): + self.curvature_matrix = inversion_0.curvature_matrix + + logger.info( + "PRELOADS - Inversion Curvature Matrix preloaded for this model-fit." + ) + + return + + if inversion_0._curvature_matrix_mapper_diag is not None: + if ( + np.max( + abs( + inversion_0._curvature_matrix_mapper_diag + - inversion_1._curvature_matrix_mapper_diag + ) + ) + < 1e-8 + ): + self.mapper_operated_mapping_matrix_dict = ( + inversion_0.mapper_operated_mapping_matrix_dict + ) + self.data_vector_mapper = inversion_0._data_vector_mapper + self.curvature_matrix_mapper_diag = ( + inversion_0._curvature_matrix_mapper_diag + ) + + logger.info( + "PRELOADS - Inversion Curvature Matrix Mapper Diag preloaded for this model-fit." + ) + + def set_regularization_matrix_and_term(self, fit_0, fit_1): + """ + If the `MassProfile`'s and `Mesh`'s in a model are fixed, the mapping of image-pixels to the + source-pixels does not change during the model-fit and therefore its associated regularization matrices are + also fixed, meaning the log determinant of the regularization matrix term of the Bayesian evidence can be + preloaded. + + This function compares the value of the log determinant of the regularization matrix of two fit's corresponding + to two model instances, and preloads this value if it is the same for both fits. + + The preload is typically used in searches where only light profiles vary (e.g. when only the lens's light is + being fitted for). + + Parameters + ---------- + fit_0 + The first fit corresponding to a model with a specific set of unit-values. + fit_1 + The second fit corresponding to a model with a different set of unit-values. + """ + self.regularization_matrix = None + self.log_det_regularization_matrix_term = None + + inversion_0 = fit_0.inversion + inversion_1 = fit_1.inversion + + if inversion_0 is None: + return + + if inversion_0.total(cls=AbstractMapper) == 0: + return + + if ( + abs( + inversion_0.log_det_regularization_matrix_term + - inversion_1.log_det_regularization_matrix_term + ) + < 1e-8 + ): + self.regularization_matrix = inversion_0.regularization_matrix + self.log_det_regularization_matrix_term = ( + inversion_0.log_det_regularization_matrix_term + ) + + logger.info( + "PRELOADS - Inversion Log Det Regularization Matrix Term preloaded for this model-fit." + ) + + def check_via_fit(self, fit): + import copy + + settings_inversion = copy.deepcopy(fit.settings_inversion) + + fit_with_preloads = fit.refit_with_new_preloads( + preloads=self, settings_inversion=settings_inversion + ) + + fit_without_preloads = fit.refit_with_new_preloads( + preloads=self.__class__(use_w_tilde=False), + settings_inversion=settings_inversion, + ) + + if os.environ.get("PYAUTOFIT_TEST_MODE") == "1": + return + + try: + if ( + abs( + fit_with_preloads.figure_of_merit + - fit_without_preloads.figure_of_merit + ) + > self.check_threshold + ): + raise exc.PreloadsException( + f""" + The log likelihood of fits using and not using preloads are not consistent by a value larger than + the preloads check threshold of {self.check_threshold}, indicating preloading has gone wrong. + + The likelihood values are: + + With Preloads: {fit_with_preloads.figure_of_merit} + Without Preloads: {fit_without_preloads.figure_of_merit} + + Double check that the model-fit is set up correctly and that the preloads are being used correctly. + + This exception can be turned off by setting the general.yaml -> test -> check_preloads to False + in the config files. However, care should be taken when doing this. + """ + ) + + except exc.InversionException: + data_vector_difference = np.max( + np.abs( + fit_with_preloads.inversion.data_vector + - fit_without_preloads.inversion.data_vector + ) + ) + + if data_vector_difference > 1.0e-4: + raise exc.PreloadsException( + f""" + The data vectors of fits using and not using preloads are not consistent, indicating + preloading has gone wrong. + + The maximum value a data vector absolute value difference is: {data_vector_difference} + """ + ) + + curvature_reg_matrix_difference = np.max( + np.abs( + fit_with_preloads.inversion.curvature_reg_matrix + - fit_without_preloads.inversion.curvature_reg_matrix + ) + ) + + if curvature_reg_matrix_difference > 1.0e-4: + raise exc.PreloadsException( + f""" + The curvature matrices of fits using and not using preloads are not consistent, indicating + preloading has gone wrong. + + The maximum value of a curvature matrix absolute value difference is: {curvature_reg_matrix_difference} + """ + ) + + @property + def info(self) -> List[str]: + """ + The information on what has or has not been preloaded, which is written to the file `preloads.summary`. + + Returns + ------- + A list of strings containing statements on what has or has not been preloaded. + """ + line = [f"W Tilde = {self.w_tilde is not None}\n"] + line += [f"Relocated Grid = {self.relocated_grid is not None}\n"] + line += [f"Mapper = {self.mapper_list is not None}\n"] + line += [ + f"Blurred Mapping Matrix = {self.operated_mapping_matrix is not None}\n" + ] + line += [ + f"Inversion Linear Func (Linear Light Profile) Dicts = {self.linear_func_operated_mapping_matrix_dict is not None}\n" + ] + line += [f"Curvature Matrix = {self.curvature_matrix is not None}\n"] + line += [ + f"Curvature Matrix Mapper Diag = {self.curvature_matrix_mapper_diag is not None}\n" + ] + line += [f"Regularization Matrix = {self.regularization_matrix is not None}\n"] + line += [ + f"Log Det Regularization Matrix Term = {self.log_det_regularization_matrix_term is not None}\n" + ] + + return line diff --git a/autoarray/structures/arrays/array_2d_util.py b/autoarray/structures/arrays/array_2d_util.py index 53e0b1890..f147baf41 100644 --- a/autoarray/structures/arrays/array_2d_util.py +++ b/autoarray/structures/arrays/array_2d_util.py @@ -913,3 +913,61 @@ def header_obj_from(file_path: Union[Path, str], hdu: int) -> Dict: hdu_list = fits.open(file_path) return hdu_list[hdu].header + + +def update_fits_file( + arr: np.ndarray, + file_path: str, + tag: Optional[str] = None, + header: Optional[fits.Header] = None, +): + """ + Update a .fits file with a new array. + + This function is used by the `fits_multi` output interface so that a single .fits file with groups of data + in hdu's can be created. + + It may receive a `tag` which is used to set the `EXTNAME` of the HDU in the .fits file and therefore is the name + of the hdu seen by the user when they open it with DS9 or other .fits software. + + A header may also be provided, which by default has the pixel scales of the array added to it. + + Parameters + ---------- + arr + The array that is written to the .fits file. + file_path + The full path of the file that is output, including the file name and ``.fits`` extension. + tag + The `EXTNAME` of the HDU in the .fits file. + header + The header of the .fits file that the array is written to, which if blank will still contain the pixel scales + of the array. + """ + + if header is None: + header = fits.Header() + + try: + y, x = map(str, arr.pixel_scales) + header["PIXSCAY"] = y + header["PIXSCAX"] = x + except AttributeError: + pass + + if conf.instance["general"]["fits"]["flip_for_ds9"]: + arr = np.flipud(arr) + + if os.path.exists(file_path): + with fits.open(file_path, mode="update") as hdul: + hdul.append(fits.ImageHDU(arr, header)) + if tag is not None: + hdul[-1].header["EXTNAME"] = tag.upper() + hdul.flush() + + else: + hdu = fits.PrimaryHDU(arr, header) + if tag is not None: + hdu.header["EXTNAME"] = tag.upper() + hdul = fits.HDUList([hdu]) + hdul.writeto(file_path, overwrite=True) diff --git a/autoarray/structures/arrays/uniform_2d.py b/autoarray/structures/arrays/uniform_2d.py index b7af699ba..823558d80 100644 --- a/autoarray/structures/arrays/uniform_2d.py +++ b/autoarray/structures/arrays/uniform_2d.py @@ -358,6 +358,103 @@ def binned_across_columns(self) -> Array1D: # binned_array = (self.native*np.invert(self.mask)).sum(axis=1)/np.invert(self.mask).sum(axis=1) return Array1D.no_mask(values=binned_array, pixel_scales=self.pixel_scale) + def brightest_coordinate_in_region_from( + self, region: Optional[Tuple[float, float, float, float]] + ) -> Tuple[float, float]: + """ + Returns the brightest pixel in the array inside an input region of form (y0, y1, x0, x1) where + the region is in scaled coordinates. + + For example, if the input region is `region=(-0.15, 0.25, 0.35, 0.55)` the code finds all pixels inside of + this region in scaled units, finds the brightest pixel of those pixels and then returns the scaled + coordinate location of that brightest pixel. + + The centre of the brightest pixel is returned, the function `brightest_sub_pixel_coordinate_in_region_from` + performs a sub pixel calculation to return the brightest sub pixel coordinate. + + Parameters + ---------- + region + The (y0, y1, x0, x1) region in scaled coordinates over which the brightest coordinate is located. + + Returns + ------- + The coordinates of the brightest pixel in scaled units (converted from pixels units). + """ + + y1, y0 = self.geometry.pixel_coordinates_2d_from( + scaled_coordinates_2d=( + region[0] - self.pixel_scales[0] / 2.0, + region[2] + self.pixel_scales[0] / 2.0, + ) + ) + x0, x1 = self.geometry.pixel_coordinates_2d_from( + scaled_coordinates_2d=( + region[1] - self.pixel_scales[1] / 2.0, + region[3] + self.pixel_scales[1] / 2.0, + ) + ) + + extracted_region = self.native[y0:y1, x0:x1] + + brightest_pixel_value = np.max(extracted_region) + extracted_pixels = np.argwhere(extracted_region == brightest_pixel_value)[0] + pixel_coordinates_2d = (y0 + extracted_pixels[0], x0 + extracted_pixels[1]) + + return self.geometry.scaled_coordinates_2d_from( + pixel_coordinates_2d=pixel_coordinates_2d + ) + + def brightest_sub_pixel_coordinate_in_region_from( + self, region: Optional[Tuple[float, float, float, float]], box_size: int = 1 + ) -> Tuple[float, float]: + """ + Returns the brightest sub-pixel in the array inside an input region of form (y0, y1, x0, x1) where + the region is in scaled coordinates. + + For example, if the input region is `region=(-0.15, 0.25, 0.35, 0.55)` the code finds all pixels inside of + this region in scaled units, finds the brightest pixel of those pixels, and then on a 3x3 grid surrounding + this pixel determines the locaiton of the brightest sub pixel using a weighted centre calculation. + + The centre of the brightest pixel is returned, the function `brightest_sub_pixel_coordinate_in_region_from` + performs a sub pixel calculation to return the brightest sub pixel coordinate. + + Parameters + ---------- + region + The (y0, y1, x0, x1) region in scaled coordinates over which the brightest coordinate is located. + + Returns + ------- + The coordinates of the brightest pixel in scaled units (converted from pixels units). + """ + brightest_pixel = self.brightest_coordinate_in_region_from(region=region) + + y, x = self.geometry.pixel_coordinates_2d_from( + scaled_coordinates_2d=brightest_pixel + ) + region_start_y, region_end_y = max(0, y - box_size), min( + self.shape_native[0], y + box_size + 1 + ) + region_start_x, region_end_x = max(0, x - box_size), min( + self.shape_native[1], x + box_size + 1 + ) + region = self.native[region_start_y:region_end_y, region_start_x:region_end_x] + + y_indices, x_indices = np.meshgrid( + range(region_start_y, region_end_y), + range(region_start_x, region_end_x), + indexing="ij", + ) + + weights = region.flatten() + subpixel_y = np.sum(weights * y_indices.flatten()) / np.sum(weights) + subpixel_x = np.sum(weights * x_indices.flatten()) / np.sum(weights) + + return self.geometry.scaled_coordinates_2d_from( + pixel_coordinates_2d=(subpixel_y, subpixel_x) + ) + def zoomed_around_mask(self, buffer: int = 1) -> "Array2D": """ Extract the 2D region of an array corresponding to the rectangle encompassing all unmasked values. diff --git a/autoarray/structures/decorators/abstract.py b/autoarray/structures/decorators/abstract.py index a0b850367..c9e5fca87 100644 --- a/autoarray/structures/decorators/abstract.py +++ b/autoarray/structures/decorators/abstract.py @@ -1,8 +1,9 @@ from typing import List, Union +import numpy as np + from autoarray.mask.mask_1d import Mask1D from autoarray.mask.mask_2d import Mask2D -from autoarray.operators.over_sampling.abstract import AbstractOverSampling from autoarray.structures.grids.uniform_1d import Grid1D from autoarray.structures.grids.irregular_2d import Grid2DIrregular from autoarray.structures.grids.uniform_2d import Grid2D @@ -59,8 +60,8 @@ def mask(self) -> Union[Mask1D, Mask2D]: return self.grid.mask @property - def over_sampling(self) -> AbstractOverSampling: - return self.grid.over_sampling + def over_sample_size(self) -> np.ndarray: + return self.grid.over_sample_size def via_grid_2d(self, result): raise NotImplementedError diff --git a/autoarray/structures/decorators/to_grid.py b/autoarray/structures/decorators/to_grid.py index 94d91e15e..144137c69 100644 --- a/autoarray/structures/decorators/to_grid.py +++ b/autoarray/structures/decorators/to_grid.py @@ -22,18 +22,31 @@ def via_grid_2d(self, result) -> Union[Grid2D, List[Grid2D]]: The input result (e.g. of a decorated function) that is converted to a Grid2D or list of Grid2D objects. """ if not isinstance(result, list): + try: + over_sampled = result.over_sampled + except AttributeError: + over_sampled = None + return Grid2D( values=result, mask=self.mask, - over_sampling=self.over_sampling, + over_sample_size=self.over_sample_size, + over_sampled=over_sampled, ) + + try: + grid_over_sampled_list = [res.over_sampled for res in result] + except AttributeError: + grid_over_sampled_list = [None] * len(result) + return [ Grid2D( values=res, mask=self.mask, - over_sampling=self.over_sampling, + over_sample_size=self.over_sample_size, + over_sampled=over_sampled, ) - for res in result + for res, over_sampled in zip(result, grid_over_sampled_list) ] def via_grid_2d_irr(self, result) -> Union[Grid2DIrregular, List[Grid2DIrregular]]: diff --git a/autoarray/structures/grids/irregular_2d.py b/autoarray/structures/grids/irregular_2d.py index 520031153..68ad0dce4 100644 --- a/autoarray/structures/grids/irregular_2d.py +++ b/autoarray/structures/grids/irregular_2d.py @@ -7,7 +7,6 @@ from autoarray.mask.mask_2d import Mask2D from autoarray.structures.arrays.irregular import ArrayIrregular -from autoarray import exc from autoarray.structures.grids import grid_2d_util from autoarray.geometry import geometry_util @@ -271,152 +270,3 @@ def grid_of_closest_from(self, grid_pair: "Grid2DIrregular") -> "Grid2DIrregular grid_of_closest[i, :] = self[np.argmin(radial_distances), :] return Grid2DIrregular(values=grid_of_closest) - - -class Grid2DIrregularUniform(Grid2DIrregular): - def __init__( - self, - values: np.ndarray, - shape_native: Optional[Tuple[float, float]] = None, - pixel_scales: Optional[Tuple[float, float]] = None, - ): - """ - A collection of (y,x) coordinates which is structured as follows: - - :: - - [[x0, x1], [x0, x1]] - - The grid object does not store the coordinates as a list of tuples, but instead a 2D NumPy array of - shape [total_coordinates, 2]. They are stored as a NumPy array so the coordinates can be used efficiently for - calculations. - - The coordinates input to this function can have any of the following forms: - - :: - - [(y0,x0), (y1,x1)] - - In all cases, they will be converted to a list of tuples followed by a 2D NumPy array. - - Print methods are overidden so a user always "sees" the coordinates as the list structure. - - Like the `Grid2D` structure, `Grid2DIrregularUniform` lie on a uniform grid corresponding to values that - originate from a uniform grid. This contrasts the `Grid2DIrregular` class above. However, although this class - stores the pixel-scale and 2D shape of this grid, it does not store the mask that a `Grid2D` does that enables - the coordinates to be mapped from 1D to 2D. This is for calculations that utilize the 2d information of the - grid but do not want the memory overheads associated with the 2D mask. - - Parameters - ---------- - values - A collection of (y,x) coordinates that. - """ - - if len(values) == 0: - super().__init__(values=values) - return - - if isinstance(values[0], float): - values = [values] - - if isinstance(values[0], tuple): - values = [values] - elif isinstance(values[0], (np.ndarray, AbstractNDArray)): - if len(values[0].shape) == 1: - values = [values] - elif isinstance(values[0], list) and isinstance(values[0][0], (float)): - values = [values] - - coordinates_arr = np.concatenate([np.array(i) for i in values]) - - self._internal_list = values - - pixel_scales = geometry_util.convert_pixel_scales_2d(pixel_scales=pixel_scales) - - self.shape_native = shape_native - self.pixel_scales = pixel_scales - - super().__init__(coordinates_arr) - - def __array_finalize__(self, obj): - if hasattr(obj, "_internal_list"): - self._internal_list = obj._internal_list - - if hasattr(obj, "shape_native"): - self.shape_native = obj.shape_native - - if hasattr(obj, "pixel_scales"): - self.pixel_scales = obj.pixel_scales - - @property - def pixel_scale(self) -> float: - if self.pixel_scales[0] == self.pixel_scales[1]: - return self.pixel_scales[0] - else: - logger.warning( - f""" - The `Grid2DIrregular` has pixel scales of {self.pixel_scales}, which are not the same in both - dimensions. This means that the pixel scale of the grid is not a single value and may cause - issues with calculations that assume a uniform pixel scale. - """ - ) - - @classmethod - def from_grid_sparse_uniform_upscale( - cls, - grid_sparse_uniform: np.ndarray, - upscale_factor: int, - pixel_scales, - shape_native=None, - ) -> "Grid2DIrregularUniform": - pixel_scales = geometry_util.convert_pixel_scales_2d(pixel_scales=pixel_scales) - - grid_upscaled_1d = grid_2d_util.grid_2d_slim_upscaled_from( - grid_slim=np.array(grid_sparse_uniform), - upscale_factor=upscale_factor, - pixel_scales=pixel_scales, - ) - - pixel_scales = ( - pixel_scales[0] / upscale_factor, - pixel_scales[1] / upscale_factor, - ) - - return Grid2DIrregularUniform( - values=grid_upscaled_1d, - pixel_scales=pixel_scales, - shape_native=shape_native, - ) - - def grid_from(self, grid_slim: np.ndarray) -> "Grid2DIrregularUniform": - """ - Create a `Grid2DIrregularUniform` object from a 2D NumPy array of values of shape [total_coordinates, 2]. The - `Grid2DIrregularUniform` are structured following this *GridIrregular2D* instance. - """ - return Grid2DIrregularUniform( - values=grid_slim, - pixel_scales=self.pixel_scales, - shape_native=self.shape_native, - ) - - def grid_2d_via_deflection_grid_from( - self, deflection_grid: np.ndarray - ) -> "Grid2DIrregularUniform": - """ - Returns a new Grid2DIrregular from this grid coordinates, where the (y,x) coordinates of this grid have a - grid of (y,x) values, termed the deflection grid, subtracted from them to determine the new grid of (y,x) - values. - - This is used by PyAutoLens to perform grid ray-tracing. - - Parameters - ---------- - deflection_grid - The grid of (y,x) coordinates which is subtracted from this grid. - """ - return Grid2DIrregularUniform( - values=self - deflection_grid, - pixel_scales=self.pixel_scales, - shape_native=self.shape_native, - ) diff --git a/autoarray/structures/grids/uniform_2d.py b/autoarray/structures/grids/uniform_2d.py index b1c147af2..23e02f97f 100644 --- a/autoarray/structures/grids/uniform_2d.py +++ b/autoarray/structures/grids/uniform_2d.py @@ -1,16 +1,12 @@ from __future__ import annotations from autoarray.numpy_wrapper import np, use_jax from pathlib import Path -from typing import TYPE_CHECKING, List, Optional, Tuple, Union - -if TYPE_CHECKING: - from autoarray.operators.over_sampling.abstract import AbstractOverSampling +from typing import List, Optional, Tuple, Union from autoconf import conf from autoconf import cached_property from autoarray.mask.mask_2d import Mask2D -from autoarray.operators.over_sampling.abstract import AbstractOverSampler from autoarray.structures.abstract_structure import Structure from autoarray.structures.arrays.uniform_2d import Array2D from autoarray.structures.grids.irregular_2d import Grid2DIrregular @@ -18,7 +14,9 @@ from autoarray.structures.arrays import array_2d_util from autoarray.structures.grids import grid_2d_util from autoarray.geometry import geometry_util +from autoarray.operators.over_sampling import over_sample_util +from autoarray import exc from autoarray import type as ty @@ -28,8 +26,8 @@ def __init__( values: Union[np.ndarray, List], mask: Mask2D, store_native: bool = False, - over_sampling: Optional[AbstractOverSampling] = None, - over_sampling_non_uniform: Optional[AbstractOverSampling] = None, + over_sample_size: Union[int, Array2D] = 4, + over_sampled: Optional[Grid2D] = None, *args, **kwargs, ): @@ -127,10 +125,6 @@ def __init__( - grid[3,4,0] = 1.5 - grid[3,4,1] = -0.5 - - - - **Grid2D Mapping:** Every set of (y,x) coordinates in a pixel of the grid maps to an unmasked pixel in the mask. For a uniform @@ -155,10 +149,14 @@ def __init__( store_native If True, the ndarray is stored in its native format [total_y_pixels, total_x_pixels, 2]. This avoids mapping large data arrays to and from the slim / native formats, which can be a computational bottleneck. - over_sampling - The over sampling scheme, which divides the grid into a sub grid of smaller pixels when computing values - (e.g. images) from the grid so as to approximate the 2D line integral of the amount of light that falls + over_sample_size + The over sampling scheme size, which divides the grid into a sub grid of smaller pixels when computing + values (e.g. images) from the grid to approximate the 2D line integral of the amount of light that falls into each pixel. + over_sampled + The over sampled grid of (y,x) coordinates, which can be passed in manually because if the grid is + not uniform (e.g. due to gravitational lensing) is cannot be computed internally in this function. If the + over sampled grid is not passed in it is computed assuming uniformity. """ values = grid_2d_util.convert_grid_2d( grid_2d=values, @@ -172,8 +170,26 @@ def __init__( grid_2d_util.check_grid_2d(grid_2d=values) - self.over_sampling = over_sampling - self.over_sampling_non_uniform = over_sampling_non_uniform + over_sample_size = over_sample_util.over_sample_size_convert_to_array_2d_from( + over_sample_size=over_sample_size, mask=mask + ) + + from autoarray.operators.over_sampling.over_sampler import OverSampler + + self.over_sampler = OverSampler(sub_size=over_sample_size, mask=mask) + + if over_sampled is None: + self.over_sampled = ( + over_sample_util.grid_2d_slim_over_sampled_via_mask_from( + mask_2d=np.array(self.mask), + pixel_scales=self.mask.pixel_scales, + sub_size=np.array(self.over_sampler.sub_size).astype("int"), + origin=self.mask.origin, + ) + ) + + else: + self.over_sampled = over_sampled @classmethod def no_mask( @@ -182,7 +198,7 @@ def no_mask( pixel_scales: ty.PixelScales, shape_native: Tuple[int, int] = None, origin: Tuple[float, float] = (0.0, 0.0), - over_sampling: Optional[AbstractOverSampling] = None, + over_sample_size: Union[int, Array2D] = 4, ) -> "Grid2D": """ Create a Grid2D (see *Grid2D.__new__*) by inputting the grid coordinates in 1D or 2D, automatically @@ -229,7 +245,7 @@ def no_mask( return Grid2D( values=values, mask=mask, - over_sampling=over_sampling, + over_sample_size=over_sample_size, ) @classmethod @@ -240,7 +256,7 @@ def from_yx_1d( shape_native: Tuple[int, int], pixel_scales: ty.PixelScales, origin: Tuple[float, float] = (0.0, 0.0), - over_sampling: Optional[AbstractOverSampling] = None, + over_sample_size: Union[int, Array2D] = 4, ) -> "Grid2D": """ Create a Grid2D (see *Grid2D.__new__*) by inputting the grid coordinates as 1D y and x values. @@ -304,7 +320,7 @@ def from_yx_1d( shape_native=shape_native, pixel_scales=pixel_scales, origin=origin, - over_sampling=over_sampling, + over_sample_size=over_sample_size, ) @classmethod @@ -314,7 +330,7 @@ def from_yx_2d( x: Union[np.ndarray, List], pixel_scales: ty.PixelScales, origin: Tuple[float, float] = (0.0, 0.0), - over_sampling: Optional[AbstractOverSampling] = None, + over_sample_size: Union[int, Array2D] = 4, ) -> "Grid2D": """ Create a Grid2D (see *Grid2D.__new__*) by inputting the grid coordinates as 2D y and x values. @@ -359,7 +375,7 @@ def from_yx_2d( values=np.stack((y, x), axis=-1), pixel_scales=pixel_scales, origin=origin, - over_sampling=over_sampling, + over_sample_size=over_sample_size, ) @classmethod @@ -367,7 +383,7 @@ def from_extent( cls, extent: Tuple[float, float, float, float], shape_native: Tuple[int, int], - over_sampling: Optional[AbstractOverSampling] = None, + over_sample_size: Union[int, Array2D] = 4, ) -> "Grid2D": """ Create a Grid2D (see *Grid2D.__new__*) by inputting the extent of the (y,x) grid coordinates as an input @@ -416,7 +432,7 @@ def from_extent( return Grid2D.no_mask( values=grid_2d, pixel_scales=pixel_scales, - over_sampling=over_sampling, + over_sample_size=over_sample_size, ) @classmethod @@ -425,7 +441,7 @@ def uniform( shape_native: Tuple[int, int], pixel_scales: ty.PixelScales, origin: Tuple[float, float] = (0.0, 0.0), - over_sampling: Optional[AbstractOverSampling] = None, + over_sample_size: Union[int, Array2D] = 4, ) -> "Grid2D": """ Create a `Grid2D` (see *Grid2D.__new__*) as a uniform grid of (y,x) values given an input `shape_native` and @@ -454,7 +470,7 @@ def uniform( shape_native=shape_native, pixel_scales=pixel_scales, origin=origin, - over_sampling=over_sampling, + over_sample_size=over_sample_size, ) @classmethod @@ -463,7 +479,7 @@ def bounding_box( bounding_box: np.ndarray, shape_native: Tuple[int, int], buffer_around_corners: bool = False, - over_sampling: Optional[AbstractOverSampling] = None, + over_sample_size: Union[int, Array2D] = 4, ) -> "Grid2D": """ Create a Grid2D (see *Grid2D.__new__*) from an input bounding box with coordinates [y_min, y_max, x_min, x_max], @@ -506,14 +522,14 @@ def bounding_box( shape_native=shape_native, pixel_scales=pixel_scales, origin=origin, - over_sampling=over_sampling, + over_sample_size=over_sample_size, ) @classmethod def from_mask( cls, mask: Mask2D, - over_sampling: Optional[AbstractOverSampling] = None, + over_sample_size: Union[int, Array2D] = 4, ) -> "Grid2D": """ Create a Grid2D (see *Grid2D.__new__*) from a mask, where only unmasked pixels are included in the grid (if the @@ -536,7 +552,7 @@ def from_mask( return Grid2D( values=grid_1d, mask=mask, - over_sampling=over_sampling, + over_sample_size=over_sample_size, ) @classmethod @@ -545,7 +561,7 @@ def from_fits( file_path: Union[Path, str], pixel_scales: ty.PixelScales, origin: Tuple[float, float] = (0.0, 0.0), - over_sampling: Optional[AbstractOverSampling] = None, + over_sample_size: Union[int, Array2D] = 4, ) -> "Grid2D": """ Create a Grid2D (see *Grid2D.__new__*) from a mask, where only unmasked pixels are included in the grid (if the @@ -565,7 +581,7 @@ def from_fits( values=grid_2d, pixel_scales=pixel_scales, origin=origin, - over_sampling=over_sampling, + over_sample_size=over_sample_size, ) @classmethod @@ -573,7 +589,7 @@ def blurring_grid_from( cls, mask: Mask2D, kernel_shape_native: Tuple[int, int], - over_sampling: Optional[AbstractOverSampling] = None, + over_sample_size: Union[int, Array2D] = 4, ) -> "Grid2D": """ Setup a blurring-grid from a mask, where a blurring grid consists of all pixels that are masked (and @@ -661,7 +677,7 @@ def blurring_grid_from( return cls.from_mask( mask=blurring_mask, - over_sampling=over_sampling, + over_sample_size=over_sample_size, ) def subtracted_from(self, offset: Tuple[(float, float), np.ndarray]) -> "Grid2D": @@ -677,9 +693,13 @@ def subtracted_from(self, offset: Tuple[(float, float), np.ndarray]) -> "Grid2D" return Grid2D( values=self - np.array(offset), mask=mask, - over_sampling=self.over_sampling, + over_sample_size=self.over_sample_size, ) + @property + def over_sample_size(self): + return self.over_sampler.sub_size + @property def slim(self) -> "Grid2D": """ @@ -692,7 +712,7 @@ def slim(self) -> "Grid2D": return Grid2D( values=self, mask=self.mask, - over_sampling=self.over_sampling, + over_sample_size=self.over_sample_size, ) @property @@ -709,14 +729,10 @@ def native(self) -> "Grid2D": return Grid2D( values=self, mask=self.mask, - over_sampling=self.over_sampling, + over_sample_size=self.over_sample_size, store_native=True, ) - @cached_property - def over_sampler(self) -> AbstractOverSampler: - return self.over_sampling.over_sampler_from(mask=self.mask) - @property def flipped(self) -> "Grid2D": """ @@ -751,7 +767,7 @@ def grid_2d_via_deflection_grid_from(self, deflection_grid: "Grid2D") -> "Grid2D return Grid2D( values=self - deflection_grid, mask=self.mask, - over_sampling=self.over_sampling, + over_sample_size=self.over_sample_size, ) def blurring_grid_via_kernel_shape_from( @@ -768,12 +784,10 @@ def blurring_grid_via_kernel_shape_from( The 2D shape of the kernel which convolves signal from masked pixels to unmasked pixels. """ - from autoarray.operators.over_sampling.uniform import OverSamplingUniform - return Grid2D.blurring_grid_from( mask=self.mask, kernel_shape_native=kernel_shape_native, - over_sampling=OverSamplingUniform(sub_size=1), + over_sample_size=1, ) def grid_with_coordinates_within_distance_removed_from( @@ -810,7 +824,7 @@ def grid_with_coordinates_within_distance_removed_from( return Grid2D.from_mask( mask=mask, - over_sampling=self.over_sampling, + over_sample_size=self.over_sample_size.apply_mask(mask=mask), ) def squared_distances_to_coordinate_from( @@ -1065,6 +1079,8 @@ def padded_grid_from(self, kernel_shape_native: Tuple[int, int]) -> "Grid2D": kernel_shape_native The 2D shape of the kernel which convolves signal from masked pixels to unmasked pixels. """ + if kernel_shape_native[0] % 2 == 0 or kernel_shape_native[1] % 2 == 0: + raise exc.KernelException("Kernel2D Kernel2D must be odd") shape = self.mask.shape @@ -1078,11 +1094,22 @@ def padded_grid_from(self, kernel_shape_native: Tuple[int, int]) -> "Grid2D": pixel_scales=self.mask.pixel_scales, ) - return Grid2D.from_mask( - mask=padded_mask, - over_sampling=self.over_sampling, + pad_width = ( + (padded_shape[0] - shape[0]) // 2, + (padded_shape[1] - shape[1]) // 2, ) + over_sample_size = np.pad( + self.over_sample_size.native, + pad_width, + mode="constant", + constant_values=1, + ) + + over_sample_size[over_sample_size == 0] = 1 + + return Grid2D.from_mask(mask=padded_mask, over_sample_size=over_sample_size) + @cached_property def is_uniform(self) -> bool: """ @@ -1105,3 +1132,32 @@ def is_uniform(self) -> bool: return False return True + + def apply_over_sampling( + self, over_sample_size: Union[int, np.ndarray] + ) -> "AbstractDataset": + """ + Apply new over sampling to the grid. + + This method is used to change the over sampling of the grid, for example when the user wishes to perform + over sampling with a higher sub grid size. + + Parameters + ---------- + over_sample_size + The over sampling scheme size, which divides the grid into a sub grid of smaller pixels when computing + values (e.g. images) from the grid to approximate the 2D line integral of the amount of light that falls + into each pixel. + """ + if not self.is_uniform: + raise exc.GridException( + """ + Cannot apply over sampling to a Grid2D which is not uniform. + """ + ) + + return Grid2D( + values=self, + mask=self.mask, + over_sample_size=over_sample_size, + ) diff --git a/autoarray/structures/mock/mock_decorators.py b/autoarray/structures/mock/mock_decorators.py index 366a78670..876b456d7 100644 --- a/autoarray/structures/mock/mock_decorators.py +++ b/autoarray/structures/mock/mock_decorators.py @@ -81,144 +81,11 @@ def ndarray_2d_yx_from(profile, grid, *args, **kwargs): return 2.0 * grid -class MockGridLikeIteratorObj: - def __init__(self): - self.centre = (0.0, 0.0) - - @property - def sersic_constant(self): - return ( - (2 * 2.0) - - (1.0 / 3.0) - + (4.0 / (405.0 * 2.0)) - + (46.0 / (25515.0 * 2.0**2)) - + (131.0 / (1148175.0 * 2.0**3)) - - (2194697.0 / (30690717750.0 * 2.0**4)) - ) - - def radial_grid_from(self, grid): - return np.sqrt(np.add(np.square(grid[:, 0]), np.square(grid[:, 1]))) - - def angle_to_profile_grid_from(self, grid_angles): - """The angle between each (y,x) coordinate on the grid and the profile, in radians. - - Parameters - ---------- - grid_angles - The angle theta counter-clockwise from the positive x-axis to each coordinate in radians. - """ - return np.cos(grid_angles), np.sin(grid_angles) - - def _cartesian_grid_via_radial_from(self, grid, radius): - """ - Convert a grid of (y,x) coordinates with their specified circular radii to their original (y,x) Cartesian - coordinates. - - Parameters - ---------- - grid - The (y, x) coordinates in the reference frame of the profile. - radius - The circular radius of each coordinate from the profile center. - """ - grid_angles = np.arctan2(grid[:, 0], grid[:, 1]) - cos_theta, sin_theta = self.angle_to_profile_grid_from(grid_angles=grid_angles) - return np.multiply(radius[:, None], np.vstack((sin_theta, cos_theta)).T) - - @over_sample - @decorators.to_array - def ndarray_1d_from(self, grid, *args, **kwargs) -> np.ndarray: - """ - Mock function mimicking the behaviour of a class function which given an input 1D grid, returns a 1D ndarray - of shape [total_masked_grid_pixels]. - - Such functions are common in **PyAutoGalaxy** for light and mass profile objects. - """ - grid_radii = self.radial_grid_from(grid=grid) - return np.exp( - np.multiply( - -self.sersic_constant, - np.add(np.power(np.divide(grid_radii, 0.2), 1.0 / 2.0), -1), - ) - ) - - @decorators.to_grid - def ndarray_2d_from(self, grid, *args, **kwargs): - """ - Mock function mimicking the behaviour of a class function which given an input grid, returns a 2D ndarray - of shape [total_masked_grid_pixels, 2]. - - Such functions are common in **PyAutoGalaxy** for light and mass profile objects. - """ - return self._cartesian_grid_via_radial_from( - grid=grid, radius=np.full(grid.shape[0], 2.0) - ) - - @decorators.to_vector_yx - def ndarray_yx_2d_from(self, grid, *args, **kwargs): - """ - Mock function mimicking the behaviour of a class function which given an input grid, returns a 2D ndarray - of shape [total_masked_grid_pixels] which represents a vector field. - - Such functions are common in **PyAutoGalaxy** for light and mass profile objects. - """ - return self._cartesian_grid_via_radial_from( - grid=grid, radius=np.full(grid.shape[0], 2.0) - ) - - @decorators.to_array - def ndarray_1d_list_from(self, grid, *args, **kwargs): - """ - Mock function mimicking the behaviour of a class function which given an input 1D grid, returns a list of 1D - ndarrays of shape [total_masked_grid_pixels]. - - Such functions are common in **PyAutoGalaxy** for light and mass profile objects. - """ - grid_radii = self.radial_grid_from(grid=grid) - return [ - np.exp( - np.multiply( - -self.sersic_constant, - np.add(np.power(np.divide(grid_radii, 0.2), 1.0 / 2.0), -1), - ) - ) - ] - - @decorators.to_grid - def ndarray_2d_list_from(self, grid, *args, **kwargs): - """ - Mock function mimicking the behaviour of a class function which given an input grid, returns a 2D list of - ndarrays of shape [total_masked_grid_pixels, 2]. - - Such functions are common in **PyAutoGalaxy** for light and mass profile objects. - """ - return [ - self._cartesian_grid_via_radial_from( - grid=grid, radius=np.full(grid.shape[0], 2.0) - ) - ] - - @decorators.to_vector_yx - def ndarray_yx_2d_list_from(self, grid, *args, **kwargs): - """ - Mock function mimicking the behaviour of a class function which given an input grid, returns a list of 2D - ndarrays of shape [total_masked_grid_pixels] which represents a vector field. - - Such functions are common in **PyAutoGalaxy** for light and mass profile objects. - """ - return [ - self._cartesian_grid_via_radial_from( - grid=grid, radius=np.full(grid.shape[0], 2.0) - ) - ] - - class MockGrid1DLikeObj: def __init__(self, centre=(0.0, 0.0), angle=0.0): self.centre = centre self.angle = angle - @over_sample @decorators.project_grid def ndarray_1d_from(self, grid, *args, **kwargs): return np.ones(shape=grid.shape[0]) diff --git a/autoarray/structures/triangles/abstract.py b/autoarray/structures/triangles/abstract.py index 919b5909b..880eea2f7 100644 --- a/autoarray/structures/triangles/abstract.py +++ b/autoarray/structures/triangles/abstract.py @@ -9,81 +9,25 @@ class AbstractTriangles(ABC): - def __init__( - self, - indices, - vertices, - **kwargs, - ): - """ - Represents a set of triangles in efficient NumPy arrays. - - Parameters - ---------- - indices - The indices of the vertices of the triangles. This is a 2D array where each row is a triangle - with the three indices of the vertices. - vertices - The vertices of the triangles. - """ - self.indices = indices - self.vertices = vertices - def __len__(self): return len(self.triangles) @property + @abstractmethod def area(self) -> float: """ The total area covered by the triangles. """ - triangles = self.triangles - return 0.5 * np.nansum( - np.abs( - (triangles[:, 0, 0] * (triangles[:, 1, 1] - triangles[:, 2, 1])) - + (triangles[:, 1, 0] * (triangles[:, 2, 1] - triangles[:, 0, 1])) - + (triangles[:, 2, 0] * (triangles[:, 0, 1] - triangles[:, 1, 1])) - ) - ) @property @abstractmethod - def numpy(self): + def indices(self): pass - def _up_sample_triangle(self): - triangles = self.triangles - - m01 = (triangles[:, 0] + triangles[:, 1]) / 2 - m12 = (triangles[:, 1] + triangles[:, 2]) / 2 - m20 = (triangles[:, 2] + triangles[:, 0]) / 2 - - return self.numpy.concatenate( - [ - self.numpy.stack([triangles[:, 1], m12, m01], axis=1), - self.numpy.stack([triangles[:, 2], m20, m12], axis=1), - self.numpy.stack([m01, m12, m20], axis=1), - self.numpy.stack([triangles[:, 0], m01, m20], axis=1), - ], - axis=0, - ) - - def _neighborhood_triangles(self): - triangles = self.triangles - - new_v0 = triangles[:, 1] + triangles[:, 2] - triangles[:, 0] - new_v1 = triangles[:, 0] + triangles[:, 2] - triangles[:, 1] - new_v2 = triangles[:, 0] + triangles[:, 1] - triangles[:, 2] - - return self.numpy.concatenate( - [ - self.numpy.stack([new_v0, triangles[:, 1], triangles[:, 2]], axis=1), - self.numpy.stack([triangles[:, 0], new_v1, triangles[:, 2]], axis=1), - self.numpy.stack([triangles[:, 0], triangles[:, 1], new_v2], axis=1), - triangles, - ], - axis=0, - ) + @property + @abstractmethod + def vertices(self): + pass def __str__(self): return f"{self.__class__.__name__} with {len(self.indices)} triangles" @@ -97,6 +41,7 @@ def triangles(self): pass @classmethod + @abstractmethod def for_limits_and_scale( cls, y_min: float, @@ -106,69 +51,7 @@ def for_limits_and_scale( scale: float, **kwargs, ) -> "AbstractTriangles": - height = scale * HEIGHT_FACTOR - - vertices = [] - indices = [] - vertex_dict = {} - - def add_vertex(v): - if v not in vertex_dict: - vertex_dict[v] = len(vertices) - vertices.append(v) - return vertex_dict[v] - - rows = [] - for row_y in np.arange(y_min, y_max + height, height): - row = [] - offset = (len(rows) % 2) * scale / 2 - for col_x in np.arange(x_min - offset, x_max + scale, scale): - row.append((row_y, col_x)) - rows.append(row) - - for i in range(len(rows) - 1): - row = rows[i] - next_row = rows[i + 1] - for j in range(len(row)): - if i % 2 == 0 and j < len(next_row) - 1: - t1 = [ - add_vertex(row[j]), - add_vertex(next_row[j]), - add_vertex(next_row[j + 1]), - ] - if j < len(row) - 1: - t2 = [ - add_vertex(row[j]), - add_vertex(row[j + 1]), - add_vertex(next_row[j + 1]), - ] - indices.append(t2) - elif i % 2 == 1 and j < len(next_row) - 1: - t1 = [ - add_vertex(row[j]), - add_vertex(next_row[j]), - add_vertex(row[j + 1]), - ] - indices.append(t1) - if j < len(next_row) - 1: - t2 = [ - add_vertex(next_row[j]), - add_vertex(next_row[j + 1]), - add_vertex(row[j + 1]), - ] - indices.append(t2) - else: - continue - indices.append(t1) - - vertices = np.array(vertices) - indices = np.array(indices) - - return cls( - indices=indices, - vertices=vertices, - **kwargs, - ) + pass @classmethod def for_grid( @@ -251,7 +134,7 @@ def containing_indices(self, shape: Shape) -> np.ndarray: Returns ------- - The triangles that intersect the shape. + The indices of triangles that intersect the shape. """ @abstractmethod diff --git a/autoarray/structures/triangles/array/__init__.py b/autoarray/structures/triangles/array/__init__.py new file mode 100644 index 000000000..0fade4b81 --- /dev/null +++ b/autoarray/structures/triangles/array/__init__.py @@ -0,0 +1,6 @@ +from .array import ArrayTriangles + +try: + from .jax_array import ArrayTriangles as JAXArrayTriangles +except ImportError: + pass diff --git a/autoarray/structures/triangles/array/abstract_array.py b/autoarray/structures/triangles/array/abstract_array.py new file mode 100644 index 000000000..d0f8620ee --- /dev/null +++ b/autoarray/structures/triangles/array/abstract_array.py @@ -0,0 +1,211 @@ +from abc import abstractmethod + +import numpy as np + +from autoarray import Grid2D, AbstractTriangles +from autoarray.structures.triangles.abstract import HEIGHT_FACTOR + + +class AbstractArrayTriangles(AbstractTriangles): + def __init__( + self, + indices, + vertices, + **kwargs, + ): + """ + Represents a set of triangles in efficient NumPy arrays. + + Parameters + ---------- + indices + The indices of the vertices of the triangles. This is a 2D array where each row is a triangle + with the three indices of the vertices. + vertices + The vertices of the triangles. + """ + self._indices = indices + self._vertices = vertices + + @property + def indices(self): + return self._indices + + @property + def vertices(self): + return self._vertices + + def __len__(self): + return len(self.triangles) + + @property + def area(self) -> float: + """ + The total area covered by the triangles. + """ + triangles = self.triangles + return ( + 0.5 + * np.abs( + (triangles[:, 0, 0] * (triangles[:, 1, 1] - triangles[:, 2, 1])) + + (triangles[:, 1, 0] * (triangles[:, 2, 1] - triangles[:, 0, 1])) + + (triangles[:, 2, 0] * (triangles[:, 0, 1] - triangles[:, 1, 1])) + ).sum() + ) + + @property + @abstractmethod + def numpy(self): + pass + + def _up_sample_triangle(self): + triangles = self.triangles + + m01 = (triangles[:, 0] + triangles[:, 1]) / 2 + m12 = (triangles[:, 1] + triangles[:, 2]) / 2 + m20 = (triangles[:, 2] + triangles[:, 0]) / 2 + + return self.numpy.concatenate( + [ + self.numpy.stack([triangles[:, 1], m12, m01], axis=1), + self.numpy.stack([triangles[:, 2], m20, m12], axis=1), + self.numpy.stack([m01, m12, m20], axis=1), + self.numpy.stack([triangles[:, 0], m01, m20], axis=1), + ], + axis=0, + ) + + def _neighborhood_triangles(self): + triangles = self.triangles + + new_v0 = triangles[:, 1] + triangles[:, 2] - triangles[:, 0] + new_v1 = triangles[:, 0] + triangles[:, 2] - triangles[:, 1] + new_v2 = triangles[:, 0] + triangles[:, 1] - triangles[:, 2] + + return self.numpy.concatenate( + [ + self.numpy.stack([new_v0, triangles[:, 1], triangles[:, 2]], axis=1), + self.numpy.stack([triangles[:, 0], new_v1, triangles[:, 2]], axis=1), + self.numpy.stack([triangles[:, 0], triangles[:, 1], new_v2], axis=1), + triangles, + ], + axis=0, + ) + + def __str__(self): + return f"{self.__class__.__name__} with {len(self.indices)} triangles" + + def __repr__(self): + return str(self) + + @classmethod + def for_limits_and_scale( + cls, + y_min: float, + y_max: float, + x_min: float, + x_max: float, + scale: float, + **kwargs, + ) -> "AbstractTriangles": + height = scale * HEIGHT_FACTOR + + vertices = [] + indices = [] + vertex_dict = {} + + def add_vertex(v): + if v not in vertex_dict: + vertex_dict[v] = len(vertices) + vertices.append(v) + return vertex_dict[v] + + rows = [] + for row_y in np.arange(y_min, y_max + height, height): + row = [] + offset = (len(rows) % 2) * scale / 2 + for col_x in np.arange(x_min - offset, x_max + scale, scale): + row.append((row_y, col_x)) + rows.append(row) + + for i in range(len(rows) - 1): + row = rows[i] + next_row = rows[i + 1] + for j in range(len(row)): + if i % 2 == 0 and j < len(next_row) - 1: + t1 = [ + add_vertex(row[j]), + add_vertex(next_row[j]), + add_vertex(next_row[j + 1]), + ] + if j < len(row) - 1: + t2 = [ + add_vertex(row[j]), + add_vertex(row[j + 1]), + add_vertex(next_row[j + 1]), + ] + indices.append(t2) + elif i % 2 == 1 and j < len(next_row) - 1: + t1 = [ + add_vertex(row[j]), + add_vertex(next_row[j]), + add_vertex(row[j + 1]), + ] + indices.append(t1) + if j < len(next_row) - 1: + t2 = [ + add_vertex(next_row[j]), + add_vertex(next_row[j + 1]), + add_vertex(row[j + 1]), + ] + indices.append(t2) + else: + continue + indices.append(t1) + + vertices = np.array(vertices) + indices = np.array(indices) + + return cls( + indices=indices, + vertices=vertices, + **kwargs, + ) + + @classmethod + def for_grid( + cls, + grid: Grid2D, + **kwargs, + ) -> "AbstractTriangles": + """ + Create a grid of equilateral triangles from a regular grid. + + Parameters + ---------- + grid + The regular grid to convert to a grid of triangles. + + Returns + ------- + The grid of triangles. + """ + + scale = grid.pixel_scale + + y = grid[:, 0] + x = grid[:, 1] + + y_min = y.min() + y_max = y.max() + x_min = x.min() + x_max = x.max() + + return cls.for_limits_and_scale( + y_min, + y_max, + x_min, + x_max, + scale, + **kwargs, + ) diff --git a/autoarray/structures/triangles/array.py b/autoarray/structures/triangles/array/array.py similarity index 95% rename from autoarray/structures/triangles/array.py rename to autoarray/structures/triangles/array/array.py index 4bf03f9ef..06bb5dc89 100644 --- a/autoarray/structures/triangles/array.py +++ b/autoarray/structures/triangles/array/array.py @@ -1,10 +1,12 @@ +from abc import ABC + import numpy as np -from autoarray.structures.triangles.abstract import AbstractTriangles +from autoarray.structures.triangles.array.abstract_array import AbstractArrayTriangles from autoarray.structures.triangles.shape import Shape -class ArrayTriangles(AbstractTriangles): +class ArrayTriangles(AbstractArrayTriangles, ABC): @property def triangles(self): return self.vertices[self.indices] diff --git a/autoarray/structures/triangles/jax_array.py b/autoarray/structures/triangles/array/jax_array.py similarity index 98% rename from autoarray/structures/triangles/jax_array.py rename to autoarray/structures/triangles/array/jax_array.py index 9142330ff..23b9ad3b5 100644 --- a/autoarray/structures/triangles/jax_array.py +++ b/autoarray/structures/triangles/array/jax_array.py @@ -2,13 +2,14 @@ from jax.tree_util import register_pytree_node_class from autoarray.structures.triangles.abstract import AbstractTriangles +from autoarray.structures.triangles.array.abstract_array import AbstractArrayTriangles from autoarray.structures.triangles.shape import Shape -MAX_CONTAINING_SIZE = 10 +MAX_CONTAINING_SIZE = 15 @register_pytree_node_class -class ArrayTriangles(AbstractTriangles): +class ArrayTriangles(AbstractArrayTriangles): def __init__( self, indices, diff --git a/autoarray/structures/triangles/coordinate_array/__init__.py b/autoarray/structures/triangles/coordinate_array/__init__.py new file mode 100644 index 000000000..f70bc8a9a --- /dev/null +++ b/autoarray/structures/triangles/coordinate_array/__init__.py @@ -0,0 +1,8 @@ +from .coordinate_array import CoordinateArrayTriangles + +try: + from .jax_coordinate_array import ( + CoordinateArrayTriangles as JAXCoordinateArrayTriangles, + ) +except ImportError: + pass diff --git a/autoarray/structures/triangles/abstract_coordinate_array.py b/autoarray/structures/triangles/coordinate_array/abstract_coordinate_array.py similarity index 74% rename from autoarray/structures/triangles/abstract_coordinate_array.py rename to autoarray/structures/triangles/coordinate_array/abstract_coordinate_array.py index 0a5a36b26..5c0fe799c 100644 --- a/autoarray/structures/triangles/abstract_coordinate_array.py +++ b/autoarray/structures/triangles/coordinate_array/abstract_coordinate_array.py @@ -1,4 +1,4 @@ -from abc import ABC, abstractmethod +from abc import abstractmethod, ABC import numpy as np @@ -6,11 +6,11 @@ from autoconf import cached_property -class AbstractCoordinateArray(ABC): +class AbstractCoordinateArray(AbstractTriangles, ABC): def __init__( self, coordinates: np.ndarray, - side_length: float, + side_length: float = 1.0, x_offset: float = 0.0, y_offset: float = 0.0, flipped: bool = False, @@ -33,33 +33,53 @@ def __init__( self.side_length = side_length self.flipped = flipped - self.scaling_factors = np.array( + self.scaling_factors = self.numpy.array( [0.5 * side_length, HEIGHT_FACTOR * side_length] ) self.x_offset = x_offset self.y_offset = y_offset + @property + @abstractmethod + def numpy(self): + pass + + @cached_property + def vertex_coordinates(self) -> np.ndarray: + """ + The vertices of the triangles as an Nx3x2 array. + """ + coordinates = self.coordinates + return self.numpy.concatenate( + [ + coordinates + self.flip_array * np.array([0, 1], dtype=np.int32), + coordinates + self.flip_array * np.array([1, -1], dtype=np.int32), + coordinates + self.flip_array * np.array([-1, -1], dtype=np.int32), + ], + dtype=np.int32, + ) + @cached_property def triangles(self) -> np.ndarray: """ The vertices of the triangles as an Nx3x2 array. """ centres = self.centres - return np.stack( + return self.numpy.stack( ( centres + self.flip_array - * np.array( + * self.numpy.array( [0.0, 0.5 * self.side_length * HEIGHT_FACTOR], ), centres + self.flip_array - * np.array( + * self.numpy.array( [0.5 * self.side_length, -0.5 * self.side_length * HEIGHT_FACTOR] ), centres + self.flip_array - * np.array( + * self.numpy.array( [-0.5 * self.side_length, -0.5 * self.side_length * HEIGHT_FACTOR] ), ), @@ -71,7 +91,7 @@ def centres(self) -> np.ndarray: """ The centres of the triangles. """ - return self.scaling_factors * self.coordinates + np.array( + return self.scaling_factors * self.coordinates + self.numpy.array( [self.x_offset, self.y_offset] ) @@ -130,42 +150,13 @@ def with_vertices(self, vertices: np.ndarray) -> AbstractTriangles: The new set of triangles with the new vertices. """ - @classmethod - def for_limits_and_scale( - cls, - x_min: float, - x_max: float, - y_min: float, - y_max: float, - scale: float = 1.0, - **_, - ): - coordinates = [] - x = x_min - while x < x_max: - y = y_min - while y < y_max: - coordinates.append([x, y]) - y += scale - x += scale - - x_mean = (x_min + x_max) / 2 - y_mean = (y_min + y_max) / 2 - - return cls( - coordinates=np.array(coordinates), - side_length=scale, - x_offset=x_mean, - y_offset=y_mean, - ) - @property def means(self): - return np.mean(self.triangles, axis=1) + return self.numpy.mean(self.triangles, axis=1) @property def area(self): return (3**0.5 / 4 * self.side_length**2) * len(self) def __len__(self): - return np.count_nonzero(~np.isnan(self.coordinates).any(axis=1)) + return self.numpy.count_nonzero(~self.numpy.isnan(self.coordinates).any(axis=1)) diff --git a/autoarray/structures/triangles/coordinate_array.py b/autoarray/structures/triangles/coordinate_array/coordinate_array.py similarity index 72% rename from autoarray/structures/triangles/coordinate_array.py rename to autoarray/structures/triangles/coordinate_array/coordinate_array.py index e45f3eed4..997c8ab7f 100644 --- a/autoarray/structures/triangles/coordinate_array.py +++ b/autoarray/structures/triangles/coordinate_array/coordinate_array.py @@ -1,10 +1,11 @@ import numpy as np from autoarray.structures.triangles.abstract import HEIGHT_FACTOR -from autoarray.structures.triangles.abstract_coordinate_array import ( +from autoarray.structures.triangles.coordinate_array.abstract_coordinate_array import ( AbstractCoordinateArray, ) from autoarray.structures.triangles.array import ArrayTriangles +from autoarray.structures.triangles.shape import Shape from autoconf import cached_property @@ -14,16 +15,50 @@ def flip_array(self) -> np.ndarray: """ An array of 1s and -1s to flip the triangles. """ - array = np.ones(self.coordinates.shape[0]) + array = np.ones( + self.coordinates.shape[0], + dtype=np.int32, + ) array[self.flip_mask] = -1 return array[:, np.newaxis] + @property + def numpy(self): + return np + + @classmethod + def for_limits_and_scale( + cls, + x_min: float, + x_max: float, + y_min: float, + y_max: float, + scale: float = 1.0, + **_, + ): + x_shift = int(2 * x_min / scale) + y_shift = int(y_min / (HEIGHT_FACTOR * scale)) + + coordinates = [] + + for x in range(x_shift, int(2 * x_max / scale) + 1): + for y in range(y_shift - 1, int(y_max / (HEIGHT_FACTOR * scale)) + 2): + coordinates.append([x, y]) + + return cls( + coordinates=np.array(coordinates, dtype=np.int32), + side_length=scale, + ) + def up_sample(self) -> "CoordinateArrayTriangles": """ Up-sample the triangles by adding a new vertex at the midpoint of each edge. """ - new_coordinates = np.zeros((4 * self.coordinates.shape[0], 2)) + new_coordinates = np.zeros( + (4 * self.coordinates.shape[0], 2), + dtype=np.int32, + ) n_normal = 4 * np.sum(~self.flip_mask) new_coordinates[:n_normal] = np.vstack( @@ -57,7 +92,10 @@ def neighborhood(self) -> "CoordinateArrayTriangles": Ensures that the new triangles are unique. """ - new_coordinates = np.zeros((4 * self.coordinates.shape[0], 2)) + new_coordinates = np.zeros( + (4 * self.coordinates.shape[0], 2), + dtype=np.int32, + ) n_normal = 4 * np.sum(~self.flip_mask) new_coordinates[:n_normal] = np.vstack( @@ -88,7 +126,9 @@ def neighborhood(self) -> "CoordinateArrayTriangles": def _vertices_and_indices(self): flat_triangles = self.triangles.reshape(-1, 2) vertices, inverse_indices = np.unique( - flat_triangles, axis=0, return_inverse=True + flat_triangles, + axis=0, + return_inverse=True, ) indices = inverse_indices.reshape(-1, 3) return vertices, indices @@ -131,3 +171,18 @@ def for_indexes(self, indexes: np.ndarray) -> "CoordinateArrayTriangles": x_offset=self.x_offset, flipped=self.flipped, ) + + def containing_indices(self, shape: Shape) -> np.ndarray: + """ + Find the triangles that insect with a given shape. + + Parameters + ---------- + shape + The shape + + Returns + ------- + The indices of triangles that intersect the shape. + """ + return self.with_vertices(self.vertices).containing_indices(shape) diff --git a/autoarray/structures/triangles/jax_coordinate_array.py b/autoarray/structures/triangles/coordinate_array/jax_coordinate_array.py similarity index 73% rename from autoarray/structures/triangles/jax_coordinate_array.py rename to autoarray/structures/triangles/coordinate_array/jax_coordinate_array.py index 30aec4eb4..e80facd1f 100644 --- a/autoarray/structures/triangles/jax_coordinate_array.py +++ b/autoarray/structures/triangles/coordinate_array/jax_coordinate_array.py @@ -1,16 +1,47 @@ from jax import numpy as np +import jax from autoarray.structures.triangles.abstract import HEIGHT_FACTOR -from autoarray.structures.triangles.abstract_coordinate_array import ( +from autoarray.structures.triangles.coordinate_array.abstract_coordinate_array import ( AbstractCoordinateArray, ) -from autoarray.structures.triangles.jax_array import ArrayTriangles +from autoarray.structures.triangles.array.jax_array import ArrayTriangles from autoarray.numpy_wrapper import register_pytree_node_class from autoconf import cached_property +jax.config.update("jax_enable_x64", True) + @register_pytree_node_class class CoordinateArrayTriangles(AbstractCoordinateArray): + @property + def numpy(self): + return jax.numpy + + @classmethod + def for_limits_and_scale( + cls, + x_min: float, + x_max: float, + y_min: float, + y_max: float, + scale: float = 1.0, + **_, + ): + x_shift = int(2 * x_min / scale) + y_shift = int(y_min / (HEIGHT_FACTOR * scale)) + + coordinates = [] + + for x in range(x_shift, int(2 * x_max / scale) + 1): + for y in range(y_shift - 1, int(y_max / (HEIGHT_FACTOR * scale)) + 2): + coordinates.append([x, y]) + + return cls( + coordinates=np.array(coordinates), + side_length=scale, + ) + def tree_flatten(self): return ( self.coordinates, @@ -43,9 +74,10 @@ def centres(self) -> np.ndarray: """ The centres of the triangles. """ - return self.scaling_factors * self.coordinates + np.array( + centres = self.scaling_factors * self.coordinates + np.array( [self.x_offset, self.y_offset] ) + return centres @cached_property def flip_mask(self) -> np.ndarray: @@ -64,7 +96,7 @@ def flip_array(self) -> np.ndarray: """ An array of 1s and -1s to flip the triangles. """ - array = np.where(self.flip_mask, -1.0, 1.0) + array = np.where(self.flip_mask, -1, 1) return array[:, None] def __iter__(self): @@ -83,8 +115,8 @@ def up_sample(self) -> "CoordinateArrayTriangles": shift0 = np.zeros((n, 2)) shift3 = np.tile(np.array([0, 1]), (n, 1)) - shift1 = np.stack([np.ones(n), np.where(flip_mask, 1.0, 0.0)], axis=1) - shift2 = np.stack([-np.ones(n), np.where(flip_mask, 1.0, 0.0)], axis=1) + shift1 = np.stack([np.ones(n), np.where(flip_mask, 1, 0)], axis=1) + shift2 = np.stack([-np.ones(n), np.where(flip_mask, 1, 0)], axis=1) shifts = np.stack([shift0, shift1, shift2, shift3], axis=1) coordinates_expanded = coordinates[:, None, :] @@ -103,7 +135,7 @@ def neighborhood(self) -> "CoordinateArrayTriangles": """ Create a new set of triangles that are the neighborhood of the current triangles. - Ensures that the new triangles are unique. + Ensures that the new triangles are unique and adjusts the mask accordingly. """ coordinates = self.coordinates flip_mask = self.flip_mask @@ -111,7 +143,6 @@ def neighborhood(self) -> "CoordinateArrayTriangles": shift0 = np.zeros((coordinates.shape[0], 2)) shift1 = np.tile(np.array([1, 0]), (coordinates.shape[0], 1)) shift2 = np.tile(np.array([-1, 0]), (coordinates.shape[0], 1)) - shift3 = np.where( flip_mask[:, None], np.tile(np.array([0, 1]), (coordinates.shape[0], 1)), @@ -122,16 +153,19 @@ def neighborhood(self) -> "CoordinateArrayTriangles": coordinates_expanded = coordinates[:, None, :] new_coordinates = coordinates_expanded + shifts - new_coordinates = new_coordinates.reshape(-1, 2) + expected_size = 4 * coordinates.shape[0] + unique_coords, indices = np.unique( + new_coordinates, + axis=0, + size=expected_size, + fill_value=np.nan, + return_index=True, + ) + return CoordinateArrayTriangles( - coordinates=np.unique( - new_coordinates, - axis=0, - size=4 * self.coordinates.shape[0], - fill_value=np.nan, - ), + coordinates=unique_coords, side_length=self.side_length, flipped=self.flipped, y_offset=self.y_offset, @@ -146,7 +180,13 @@ def _vertices_and_indices(self): axis=0, return_inverse=True, size=3 * self.coordinates.shape[0], + equal_nan=True, + fill_value=np.nan, ) + + nan_mask = np.isnan(vertices).any(axis=1) + inverse_indices = np.where(nan_mask[inverse_indices], -1, inverse_indices) + indices = inverse_indices.reshape(-1, 3) return vertices, indices @@ -183,7 +223,7 @@ def for_indexes(self, indexes: np.ndarray) -> "CoordinateArrayTriangles": """ mask = indexes == -1 safe_indexes = np.where(mask, 0, indexes) - coordinates = self.coordinates[safe_indexes] + coordinates = np.take(self.coordinates, safe_indexes, axis=0) coordinates = np.where(mask[:, None], np.nan, coordinates) return CoordinateArrayTriangles( @@ -193,3 +233,6 @@ def for_indexes(self, indexes: np.ndarray) -> "CoordinateArrayTriangles": x_offset=self.x_offset, flipped=self.flipped, ) + + def containing_indices(self, shape: np.ndarray) -> np.ndarray: + raise NotImplementedError("JAX ArrayTriangles are used for this method.") diff --git a/test_autoarray/config/grids.yaml b/test_autoarray/config/grids.yaml index 33f7c08c3..61c268a27 100644 --- a/test_autoarray/config/grids.yaml +++ b/test_autoarray/config/grids.yaml @@ -1,8 +1,3 @@ -interpolate: - ndarray_1d_from_grid: - MockGridLikeIteratorObj: true - ndarray_2d_from_grid: - MockGridLikeIteratorObj: true # Certain light and mass profile calculations become ill defined at (0.0, 0.0) or close to this value. This can lead # to numerical issues in the calculation of the profile, for example a np.nan may arise, crashing the code. @@ -14,41 +9,4 @@ interpolate: radial_minimum: radial_minimum: - MockGridRadialMinimum: 2.5 - - -# Over sampling is an important numerical technique, whereby light profiles images are evaluated on a higher resolution -# grid than the image data to ensure the calculation is accurate. - -# By default, a user does not specify the over sampling factor, and a default over sampling scheme is used for each -# profile. This scheme first goes to the centre of the profile, and computes circles with certain radial values -# (e.g. radii). It then assigns an over sampling `sub_size` to each circle, where the central circles have the highest -# over sampling factor and the outer circles have the lowest. - -# The size of the circles that are appropriate for determining the over sampling factor are dependent on the resolution -# of the grid. For a high resolution grid (e.g. low pixel scale), a smaller circle central circle is necessary to -# over sample the profile accurately. The config file below therefore specifies the "radial factors" used for -# automatically determining the over sampling factors for each profile, which is the factor the pixel scale is multiplied -# by to determine the circle size. - -# The config entry below defines the default over sampling factor for each profile, where: - -# radial_factor_list: The factors that are multiplied by the pixel scale to determine the circle size that is used. -# sub_size_list: The over sampling factor that is used for each circle size. - -# For the default entries below, oversampling of degree 32 x 32 is used within a circle of radius 3.01 x pixel scale, -# 4 x 4 within a circle of radius 10.01 x pixel scale and 2 x 2 for all pixels outside of this radius. - -# For unit tests, we disable this feature by setting the over sampling factors to 1.0 and the sub sizes to 1. - -over_sampling: - radial_factor_list: - MockGrid1DLikeObj: [1.0] - MockGrid2DLikeObj: [1.0] - MockGridLikeIteratorObj: [1.0] - MockGridRadialMinimum: [1.0] - sub_size_list: - MockGrid1DLikeObj: [1, 1] - MockGrid2DLikeObj: [1, 1] - MockGridLikeIteratorObj: [1, 1] - MockGridRadialMinimum: [1, 1] \ No newline at end of file + MockGridRadialMinimum: 2.5 \ No newline at end of file diff --git a/test_autoarray/conftest.py b/test_autoarray/conftest.py index 0d5437d7d..1dbe19e19 100644 --- a/test_autoarray/conftest.py +++ b/test_autoarray/conftest.py @@ -84,6 +84,11 @@ def make_grid_2d_7x7(): return fixtures.make_grid_2d_7x7() +@pytest.fixture(name="grid_2d_sub_1_7x7") +def make_grid_2d_sub_1_7x7(): + return fixtures.make_grid_2d_sub_1_7x7() + + @pytest.fixture(name="grid_2d_irregular_7x7") def make_grid_2d_irregular_7x7(): return fixtures.make_grid_2d_irregular_7x7() diff --git a/test_autoarray/dataset/abstract/test_dataset.py b/test_autoarray/dataset/abstract/test_dataset.py index 42cd37b2c..dd99c1201 100644 --- a/test_autoarray/dataset/abstract/test_dataset.py +++ b/test_autoarray/dataset/abstract/test_dataset.py @@ -65,25 +65,12 @@ def test__grid__uses_mask_and_settings( masked_imaging_7x7 = ds.AbstractDataset( data=masked_image_7x7, noise_map=masked_noise_map_7x7, - over_sampling=aa.OverSamplingDataset( - uniform=aa.OverSamplingUniform(sub_size=2) - ), + over_sample_size_lp=2, ) - assert isinstance(masked_imaging_7x7.grids.uniform, aa.Grid2D) - assert (masked_imaging_7x7.grids.uniform == grid_2d_7x7).all() - assert (masked_imaging_7x7.grids.uniform.slim == grid_2d_7x7).all() - - masked_imaging_7x7 = ds.AbstractDataset( - data=masked_image_7x7, - noise_map=masked_noise_map_7x7, - over_sampling=aa.OverSamplingDataset(uniform=aa.OverSamplingIterate()), - ) - - assert isinstance( - masked_imaging_7x7.grids.uniform.over_sampling, aa.OverSamplingIterate - ) - assert (masked_imaging_7x7.grids.uniform == grid_2d_7x7).all() + assert isinstance(masked_imaging_7x7.grids.lp, aa.Grid2D) + assert (masked_imaging_7x7.grids.lp == grid_2d_7x7).all() + assert (masked_imaging_7x7.grids.lp.slim == grid_2d_7x7).all() def test__grids_pixelization__uses_mask_and_settings( @@ -99,40 +86,32 @@ def test__grids_pixelization__uses_mask_and_settings( masked_imaging_7x7 = ds.AbstractDataset( data=masked_image_7x7, noise_map=masked_noise_map_7x7, - over_sampling=aa.OverSamplingDataset( - pixelization=aa.OverSamplingIterate(sub_steps=[2, 4]) - ), ) - assert masked_imaging_7x7.grids.pixelization.over_sampling.sub_steps == [2, 4] assert (masked_imaging_7x7.grids.pixelization == grid_2d_7x7).all() assert (masked_imaging_7x7.grids.pixelization.slim == grid_2d_7x7).all() masked_imaging_7x7 = ds.AbstractDataset( data=masked_image_7x7, noise_map=masked_noise_map_7x7, - over_sampling=aa.OverSamplingDataset( - uniform=aa.OverSamplingUniform(sub_size=2), - pixelization=aa.OverSamplingUniform(sub_size=4), - ), + over_sample_size_lp=2, + over_sample_size_pixelization=4, ) assert isinstance(masked_imaging_7x7.grids.pixelization, aa.Grid2D) - assert masked_imaging_7x7.grids.pixelization.over_sampling.sub_size == 4 + assert masked_imaging_7x7.grids.over_sample_size_pixelization[0] == 4 def test__grid_settings__sub_size(image_7x7, noise_map_7x7): dataset_7x7 = ds.AbstractDataset( data=image_7x7, noise_map=noise_map_7x7, - over_sampling=aa.OverSamplingDataset( - uniform=aa.OverSamplingUniform(sub_size=2), - pixelization=aa.OverSamplingUniform(sub_size=4), - ), + over_sample_size_lp=2, + over_sample_size_pixelization=4, ) - assert dataset_7x7.grids.uniform.over_sampling.sub_size == 2 - assert dataset_7x7.grids.pixelization.over_sampling.sub_size == 4 + assert dataset_7x7.grids.over_sample_size_lp[0] == 2 + assert dataset_7x7.grids.over_sample_size_pixelization[0] == 4 def test__new_imaging_with_arrays_trimmed_via_kernel_shape(): @@ -157,33 +136,29 @@ def test__apply_over_sampling(image_7x7, noise_map_7x7): dataset_7x7 = aa.Imaging( data=image_7x7, noise_map=noise_map_7x7, - over_sampling=aa.OverSamplingDataset( - uniform=aa.OverSamplingUniform(sub_size=2), - pixelization=aa.OverSamplingUniform(sub_size=2), - ), + over_sample_size_lp=2, + over_sample_size_pixelization=2, ) # The grid and grid_pixelizaiton are a cached_property which needs to be reset, # Which the code below tests. - grid_sub_2 = dataset_7x7.grids.uniform + grid_sub_2 = dataset_7x7.grids.lp grids_pixelization_sub_2 = dataset_7x7.grids.pixelization - dataset_7x7.grids.__dict__["uniform"][0][0] = 100.0 + dataset_7x7.grids.__dict__["lp"][0][0] = 100.0 dataset_7x7.grids.__dict__["pixelization"][0][0] = 100.0 - assert dataset_7x7.grids.uniform[0][0] == pytest.approx(100.0, 1.0e-4) + assert dataset_7x7.grids.lp[0][0] == pytest.approx(100.0, 1.0e-4) assert dataset_7x7.grids.pixelization[0][0] == pytest.approx(100.0, 1.0e-4) dataset_7x7 = dataset_7x7.apply_over_sampling( - over_sampling=aa.OverSamplingDataset( - uniform=aa.OverSamplingUniform(sub_size=4), - pixelization=aa.OverSamplingUniform(sub_size=4), - ) + over_sample_size_lp=2, + over_sample_size_pixelization=4, ) - assert dataset_7x7.over_sampling.uniform.sub_size == 4 - assert dataset_7x7.over_sampling.pixelization.sub_size == 4 + assert dataset_7x7.over_sample_size_lp.slim[0] == 2 + assert dataset_7x7.over_sample_size_pixelization.slim[0] == 4 - assert dataset_7x7.grids.uniform[0][0] == pytest.approx(3.0, 1.0e-4) + assert dataset_7x7.grids.lp[0][0] == pytest.approx(3.0, 1.0e-4) assert dataset_7x7.grids.pixelization[0][0] == pytest.approx(3.0, 1.0e-4) diff --git a/test_autoarray/dataset/imaging/test_dataset.py b/test_autoarray/dataset/imaging/test_dataset.py index 46f6adf4d..6025733f5 100644 --- a/test_autoarray/dataset/imaging/test_dataset.py +++ b/test_autoarray/dataset/imaging/test_dataset.py @@ -1,3 +1,4 @@ +import copy import os from os import path @@ -158,6 +159,21 @@ def test__apply_noise_scaling(imaging_7x7, mask_2d_7x7): assert masked_imaging_7x7.noise_map.native[4, 4] == 1e5 +def test__apply_noise_scaling__use_signal_to_noise_value(imaging_7x7, mask_2d_7x7): + imaging_7x7 = copy.copy(imaging_7x7) + + imaging_7x7.data[24] = 2.0 + + masked_imaging_7x7 = imaging_7x7.apply_noise_scaling( + mask=mask_2d_7x7, signal_to_noise_value=0.1, should_zero_data=False + ) + + assert masked_imaging_7x7.data.native[3, 4] == 1.0 + assert masked_imaging_7x7.noise_map.native[3, 4] == 10.0 + assert masked_imaging_7x7.data.native[3, 3] == 2.0 + assert masked_imaging_7x7.noise_map.native[3, 3] == 10.0 + + def test__apply_mask__noise_covariance_matrix(): image = aa.Array2D.ones(shape_native=(2, 2), pixel_scales=(1.0, 1.0)) diff --git a/test_autoarray/dataset/imaging/test_simulator.py b/test_autoarray/dataset/imaging/test_simulator.py index 9d71cff5d..3cc4182d3 100644 --- a/test_autoarray/dataset/imaging/test_simulator.py +++ b/test_autoarray/dataset/imaging/test_simulator.py @@ -13,7 +13,11 @@ def make_array_2d_7x7(): def test__via_image_from__all_features_off(image_central_delta_3x3): - simulator = aa.SimulatorImaging(exposure_time=1.0, add_poisson_noise=False) + simulator = aa.SimulatorImaging( + exposure_time=1.0, + add_poisson_noise_to_data=False, + include_poisson_noise_in_noise_map=False, + ) dataset = simulator.via_image_from(image=image_central_delta_3x3) @@ -24,10 +28,37 @@ def test__via_image_from__all_features_off(image_central_delta_3x3): assert dataset.pixel_scales == (0.1, 0.1) -def test__via_image_from__noise_off__noise_map_is_noise_value(image_central_delta_3x3): +def test__via_image_from__psf_off__include_poisson_noise_in_noise_map( + image_central_delta_3x3, +): + image = image_central_delta_3x3 + 1.0 + + simulator = aa.SimulatorImaging( + exposure_time=20.0, + add_poisson_noise_to_data=True, + include_poisson_noise_in_noise_map=True, + noise_seed=1, + ) + + dataset = simulator.via_image_from(image=image) + + assert dataset.data.native == pytest.approx( + np.array([[1.05, 1.3, 1.25], [1.05, 2.1, 1.2], [1.05, 1.3, 1.15]]), 1e-2 + ) + + assert dataset.noise_map.native == pytest.approx( + np.array([[0.229, 0.255, 0.25], [0.229, 0.324, 0.245], [0.229, 0.255, 0.240]]), + 1e-2, + ) + + +def test__via_image_from__psf_off__noise_off_value_is_noise_value( + image_central_delta_3x3, +): simulator = aa.SimulatorImaging( exposure_time=1.0, - add_poisson_noise=False, + add_poisson_noise_to_data=False, + include_poisson_noise_in_noise_map=False, noise_if_add_noise_false=0.2, noise_seed=1, ) @@ -42,48 +73,62 @@ def test__via_image_from__noise_off__noise_map_is_noise_value(image_central_delt assert np.allclose(dataset.noise_map.native, 0.2 * np.ones((3, 3))) -def test__via_image_from__psf_blurs_image_with_edge_trimming(image_central_delta_3x3): - psf = aa.Kernel2D.no_mask( - values=np.array([[0.0, 1.0, 0.0], [1.0, 2.0, 1.0], [0.0, 1.0, 0.0]]), - pixel_scales=1.0, - ) - +def test__via_image_from__psf_off__background_sky_on(image_central_delta_3x3): simulator = aa.SimulatorImaging( - exposure_time=1.0, psf=psf, add_poisson_noise=False, normalize_psf=False + exposure_time=1.0, + background_sky_level=16.0, + add_poisson_noise_to_data=True, + noise_seed=1, ) dataset = simulator.via_image_from(image=image_central_delta_3x3) assert ( dataset.data.native - == np.array([[0.0, 1.0, 0.0], [1.0, 2.0, 1.0], [0.0, 1.0, 0.0]]) + == np.array([[1.0, 5.0, 4.0], [1.0, 2.0, 1.0], [5.0, 2.0, 7.0]]) ).all() + assert dataset.noise_map.native[0, 0] == pytest.approx(4.12310, 1.0e-4) + -def test__via_image_from__setup_with_noise(image_central_delta_3x3): - image = image_central_delta_3x3 + 1.0 +def test__via_image_from__psf_on__psf_blurs_image_with_edge_trimming( + image_central_delta_3x3, +): + psf = aa.Kernel2D.no_mask( + values=np.array([[0.0, 1.0, 0.0], [1.0, 2.0, 1.0], [0.0, 1.0, 0.0]]), + pixel_scales=1.0, + ) simulator = aa.SimulatorImaging( - exposure_time=20.0, add_poisson_noise=True, noise_seed=1 + exposure_time=1.0, + psf=psf, + add_poisson_noise_to_data=False, + include_poisson_noise_in_noise_map=False, + normalize_psf=False, ) - dataset = simulator.via_image_from(image=image) + dataset = simulator.via_image_from(image=image_central_delta_3x3) - assert dataset.data.native == pytest.approx( - np.array([[1.05, 1.3, 1.25], [1.05, 2.1, 1.2], [1.05, 1.3, 1.15]]), 1e-2 - ) + assert ( + dataset.data.native + == np.array([[0.0, 1.0, 0.0], [1.0, 2.0, 1.0], [0.0, 1.0, 0.0]]) + ).all() - assert dataset.noise_map.native == pytest.approx( - np.array([[0.229, 0.255, 0.25], [0.229, 0.324, 0.245], [0.229, 0.255, 0.240]]), - 1e-2, - ) +def test__via_image_from__psf_on__disable_poisson_noise_in_data( + image_central_delta_3x3, +): + psf = aa.Kernel2D.no_mask( + values=np.array([[0.0, 1.0, 0.0], [1.0, 2.0, 1.0], [0.0, 1.0, 0.0]]), + pixel_scales=1.0, + ) -def test__via_image_from__background_sky_on(image_central_delta_3x3): simulator = aa.SimulatorImaging( - exposure_time=1.0, - background_sky_level=16.0, - add_poisson_noise=True, + exposure_time=20.0, + psf=psf, + normalize_psf=False, + add_poisson_noise_to_data=False, + include_poisson_noise_in_noise_map=True, noise_seed=1, ) @@ -91,13 +136,18 @@ def test__via_image_from__background_sky_on(image_central_delta_3x3): assert ( dataset.data.native - == np.array([[1.0, 5.0, 4.0], [1.0, 2.0, 1.0], [5.0, 2.0, 7.0]]) + == np.array([[0.0, 1.0, 0.0], [1.0, 2.0, 1.0], [0.0, 1.0, 0.0]]) ).all() - assert dataset.noise_map.native[0, 0] == pytest.approx(4.12310, 1.0e-4) + assert dataset.noise_map.native == pytest.approx( + np.array( + [[0.0, 0.22912, 0.0], [0.25495, 0.34278, 0.22912], [0.0, 0.229128, 0.0]] + ), + 1e-2, + ) -def test__via_image_from__psf_and_noise_both_on(image_central_delta_3x3): +def test__via_image_from__psf_on__psf_and_noise_both_on(image_central_delta_3x3): image = image_central_delta_3x3 + 1.0 psf = aa.Kernel2D.no_mask( @@ -108,7 +158,7 @@ def test__via_image_from__psf_and_noise_both_on(image_central_delta_3x3): simulator = aa.SimulatorImaging( exposure_time=20.0, psf=psf, - add_poisson_noise=True, + add_poisson_noise_to_data=True, noise_seed=1, normalize_psf=False, ) diff --git a/test_autoarray/dataset/plot/test_imaging_plotters.py b/test_autoarray/dataset/plot/test_imaging_plotters.py index ebec9be47..d6c6673aa 100644 --- a/test_autoarray/dataset/plot/test_imaging_plotters.py +++ b/test_autoarray/dataset/plot/test_imaging_plotters.py @@ -19,12 +19,6 @@ def test__individual_attributes_are_output( ): visuals = aplt.Visuals2D(mask=mask_2d_7x7, positions=grid_2d_irregular_7x7_list) - imaging_7x7 = imaging_7x7.apply_over_sampling( - over_sampling=aa.OverSamplingDataset( - non_uniform=aa.OverSamplingUniform(sub_size=1) - ) - ) - dataset_plotter = aplt.ImagingPlotter( dataset=imaging_7x7, visuals_2d=visuals, @@ -36,18 +30,16 @@ def test__individual_attributes_are_output( noise_map=True, psf=True, signal_to_noise_map=True, - over_sampling=True, - over_sampling_non_uniform=True, - over_sampling_pixelization=True, + over_sample_size_lp=True, + over_sample_size_pixelization=True, ) assert path.join(plot_path, "data.png") in plot_patch.paths assert path.join(plot_path, "noise_map.png") in plot_patch.paths assert path.join(plot_path, "psf.png") in plot_patch.paths assert path.join(plot_path, "signal_to_noise_map.png") in plot_patch.paths - assert path.join(plot_path, "over_sampling.png") in plot_patch.paths - assert path.join(plot_path, "over_sampling_non_uniform.png") in plot_patch.paths - assert path.join(plot_path, "over_sampling_pixelization.png") in plot_patch.paths + assert path.join(plot_path, "over_sample_size_lp.png") in plot_patch.paths + assert path.join(plot_path, "over_sample_size_pixelization.png") in plot_patch.paths plot_patch.paths = [] diff --git a/test_autoarray/fit/test_fit_dataset.py b/test_autoarray/fit/test_fit_dataset.py index 2ddb13ad6..1f8ec6380 100644 --- a/test_autoarray/fit/test_fit_dataset.py +++ b/test_autoarray/fit/test_fit_dataset.py @@ -74,5 +74,5 @@ def test__grid_offset_via_data_model(masked_imaging_7x7, model_image_7x7): assert fit.dataset_model.grid_offset == (1.0, 2.0) - assert fit.grids.uniform[0] == pytest.approx((0.0, -3.0), 1.0e-4) + assert fit.grids.lp[0] == pytest.approx((0.0, -3.0), 1.0e-4) assert fit.grids.pixelization[0] == pytest.approx((0.0, -3.0), 1.0e-4) diff --git a/test_autoarray/inversion/inversion/imaging/test_inversion_imaging_util.py b/test_autoarray/inversion/inversion/imaging/test_inversion_imaging_util.py index d2875b2bd..711135ebc 100644 --- a/test_autoarray/inversion/inversion/imaging/test_inversion_imaging_util.py +++ b/test_autoarray/inversion/inversion/imaging/test_inversion_imaging_util.py @@ -188,20 +188,16 @@ def test__data_vector_via_w_tilde_data_two_methods_agree(): # TODO : Use pytest.parameterize for sub_size in range(1, 3): - over_sampler = aa.OverSamplerUniform( - mask=mask, - sub_size=sub_size, - ) - - grid = over_sampler.over_sampled_grid + grid = aa.Grid2D.from_mask(mask=mask, over_sample_size=sub_size) mapper_grids = pixelization.mapper_grids_from( - mask=mask, border_relocator=None, source_plane_data_grid=grid + mask=mask, + border_relocator=None, + source_plane_data_grid=grid, ) mapper = aa.Mapper( mapper_grids=mapper_grids, - over_sampler=over_sampler, regularization=None, ) @@ -236,7 +232,7 @@ def test__data_vector_via_w_tilde_data_two_methods_agree(): pix_sizes_for_sub_slim_index=mapper.pix_sizes_for_sub_slim_index, pix_weights_for_sub_slim_index=mapper.pix_weights_for_sub_slim_index, pix_pixels=mapper.params, - sub_size=np.array(over_sampler.sub_size), + sub_size=np.array(grid.over_sample_size), ) data_vector_via_w_tilde = ( @@ -272,14 +268,7 @@ def test__curvature_matrix_via_w_tilde_two_methods_agree(): source_plane_data_grid=mask.derive_grid.unmasked, ) - over_sampler = aa.OverSamplerUniform( - mask=mask, - sub_size=1, - ) - - mapper = aa.Mapper( - mapper_grids=mapper_grids, over_sampler=over_sampler, regularization=None - ) + mapper = aa.Mapper(mapper_grids=mapper_grids, regularization=None) mapping_matrix = mapper.mapping_matrix @@ -319,20 +308,16 @@ def test__curvature_matrix_via_w_tilde_preload_two_methods_agree(): pixelization = aa.mesh.Rectangular(shape=(20, 20)) for sub_size in range(1, 2, 3): - over_sampler = aa.OverSamplerUniform( - mask=mask, - sub_size=sub_size, - ) - - grid = over_sampler.over_sampled_grid + grid = aa.Grid2D.from_mask(mask=mask, over_sample_size=sub_size) mapper_grids = pixelization.mapper_grids_from( - mask=mask, border_relocator=None, source_plane_data_grid=grid + mask=mask, + border_relocator=None, + source_plane_data_grid=grid, ) mapper = aa.Mapper( mapper_grids=mapper_grids, - over_sampler=over_sampler, regularization=None, ) @@ -358,7 +343,7 @@ def test__curvature_matrix_via_w_tilde_preload_two_methods_agree(): pix_sizes_for_sub_slim_index=mapper.pix_sizes_for_sub_slim_index, pix_weights_for_sub_slim_index=mapper.pix_weights_for_sub_slim_index, pix_pixels=mapper.params, - sub_size=np.array(over_sampler.sub_size), + sub_size=np.array(grid.over_sample_size), ) curvature_matrix_via_w_tilde = aa.util.inversion_imaging.curvature_matrix_via_w_tilde_curvature_preload_imaging_from( diff --git a/test_autoarray/inversion/inversion/test_abstract.py b/test_autoarray/inversion/inversion/test_abstract.py index 149e73b91..5ce2ad473 100644 --- a/test_autoarray/inversion/inversion/test_abstract.py +++ b/test_autoarray/inversion/inversion/test_abstract.py @@ -114,7 +114,7 @@ def test__curvature_matrix__via_w_tilde__identical_to_mapping(): pixel_scales=2.0, ) - grid = aa.Grid2D.from_mask(mask=mask) + grid = aa.Grid2D.from_mask(mask=mask, over_sample_size=1) mesh_0 = aa.mesh.Rectangular(shape=(3, 3)) mesh_1 = aa.mesh.Rectangular(shape=(4, 4)) @@ -135,14 +135,8 @@ def test__curvature_matrix__via_w_tilde__identical_to_mapping(): reg = aa.reg.Constant(coefficient=1.0) - over_sampler = aa.OverSamplerUniform(mask=mask, sub_size=1) - - mapper_0 = aa.Mapper( - mapper_grids=mapper_grids_0, over_sampler=over_sampler, regularization=reg - ) - mapper_1 = aa.Mapper( - mapper_grids=mapper_grids_1, over_sampler=over_sampler, regularization=reg - ) + mapper_0 = aa.Mapper(mapper_grids=mapper_grids_0, regularization=reg) + mapper_1 = aa.Mapper(mapper_grids=mapper_grids_1, regularization=reg) image = aa.Array2D.no_mask(values=np.random.random((7, 7)), pixel_scales=1.0) noise_map = aa.Array2D.no_mask(values=np.random.random((7, 7)), pixel_scales=1.0) @@ -184,7 +178,7 @@ def test__curvature_matrix_via_w_tilde__includes_source_interpolation__identical pixel_scales=2.0, ) - grid = aa.Grid2D.from_mask(mask=mask) + grid = aa.Grid2D.from_mask(mask=mask, over_sample_size=1) mesh_0 = aa.mesh.Delaunay() mesh_1 = aa.mesh.Delaunay() @@ -216,14 +210,8 @@ def test__curvature_matrix_via_w_tilde__includes_source_interpolation__identical reg = aa.reg.Constant(coefficient=1.0) - over_sampler = aa.OverSamplerUniform(mask=mask, sub_size=1) - - mapper_0 = aa.Mapper( - mapper_grids=mapper_grids_0, over_sampler=over_sampler, regularization=reg - ) - mapper_1 = aa.Mapper( - mapper_grids=mapper_grids_1, over_sampler=over_sampler, regularization=reg - ) + mapper_0 = aa.Mapper(mapper_grids=mapper_grids_0, regularization=reg) + mapper_1 = aa.Mapper(mapper_grids=mapper_grids_1, regularization=reg) image = aa.Array2D.no_mask(values=np.random.random((7, 7)), pixel_scales=1.0) noise_map = aa.Array2D.no_mask(values=np.random.random((7, 7)), pixel_scales=1.0) @@ -565,12 +553,14 @@ def test__determinant_of_positive_definite_matrix_via_cholesky(): ) -def test__errors_and_errors_with_covariance(): +def test__reconstruction_noise_map(): curvature_reg_matrix = np.array([[1.0, 1.0, 1.0], [1.0, 2.0, 1.0], [1.0, 1.0, 3.0]]) inversion = aa.m.MockInversion(curvature_reg_matrix=curvature_reg_matrix) - assert inversion.errors_with_covariance == pytest.approx( - np.array([[2.5, -1.0, -0.5], [-1.0, 1.0, 0.0], [-0.5, 0.0, 0.5]]), 1.0e-2 + assert inversion.reconstruction_noise_map_with_covariance[0, 0] == pytest.approx( + np.sqrt(2.5), 1.0e-2 + ) + assert inversion.reconstruction_noise_map == pytest.approx( + np.sqrt(np.array([2.5, 1.0, 0.5])), 1.0e-3 ) - assert inversion.errors == pytest.approx(np.array([2.5, 1.0, 0.5]), 1.0e-3) diff --git a/test_autoarray/inversion/pixelization/mappers/test_abstract.py b/test_autoarray/inversion/pixelization/mappers/test_abstract.py index 0adbe4911..7217dc79a 100644 --- a/test_autoarray/inversion/pixelization/mappers/test_abstract.py +++ b/test_autoarray/inversion/pixelization/mappers/test_abstract.py @@ -140,7 +140,7 @@ def test__adaptive_pixel_signals_from___matches_util(grid_2d_7x7, image_7x7): ) pix_weights_for_sub_slim_index = np.ones((9, 1), dtype="int") - over_sampler = aa.OverSamplerUniform(mask=grid_2d_7x7.mask, sub_size=1) + over_sampler = aa.OverSampler(mask=grid_2d_7x7.mask, sub_size=1) mapper = aa.m.MockMapper( source_plane_data_grid=grid_2d_7x7, @@ -180,9 +180,7 @@ def test__interpolated_array_from(grid_2d_7x7): source_plane_mesh_grid=mesh_grid, ) - mapper = aa.Mapper( - mapper_grids=mapper_grids, over_sampler=None, regularization=None - ) + mapper = aa.Mapper(mapper_grids=mapper_grids, regularization=None) interpolated_array_via_mapper = mapper.interpolated_array_from( values=np.array([1.0, 2.0, 3.0, 4.0, 5.0, 6.0]), @@ -204,6 +202,7 @@ def test__mapped_to_source_from(grid_2d_7x7): values=[[0.1, 0.1], [1.1, 0.6], [2.1, 0.1], [0.4, 1.1], [1.1, 7.1], [2.1, 1.1]], shape_native=(3, 2), pixel_scales=1.0, + over_sample_size=1, ) mesh_grid = aa.Mesh2DDelaunay(values=mesh_grid) @@ -214,11 +213,7 @@ def test__mapped_to_source_from(grid_2d_7x7): source_plane_mesh_grid=mesh_grid, ) - over_sampler = aa.OverSamplerUniform(mask=grid_2d_7x7.mask, sub_size=1) - - mapper = aa.Mapper( - mapper_grids=mapper_grids, over_sampler=over_sampler, regularization=None - ) + mapper = aa.Mapper(mapper_grids=mapper_grids, regularization=None) array_slim = aa.Array2D.no_mask( [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0], diff --git a/test_autoarray/inversion/pixelization/mappers/test_delaunay.py b/test_autoarray/inversion/pixelization/mappers/test_delaunay.py index 970e8555a..33b03fbb4 100644 --- a/test_autoarray/inversion/pixelization/mappers/test_delaunay.py +++ b/test_autoarray/inversion/pixelization/mappers/test_delaunay.py @@ -2,26 +2,23 @@ import autoarray as aa -def test__pix_indexes_for_sub_slim_index__matches_util(grid_2d_7x7): +def test__pix_indexes_for_sub_slim_index__matches_util(grid_2d_sub_1_7x7): mesh_grid = aa.Grid2D.no_mask( values=[[0.1, 0.1], [1.1, 0.6], [2.1, 0.1], [0.4, 1.1], [1.1, 7.1], [2.1, 1.1]], shape_native=(3, 2), pixel_scales=1.0, + over_sample_size=1, ) mesh_grid = aa.Mesh2DDelaunay(values=mesh_grid) mapper_grids = aa.MapperGrids( - mask=grid_2d_7x7.mask, - source_plane_data_grid=grid_2d_7x7, + mask=grid_2d_sub_1_7x7.mask, + source_plane_data_grid=grid_2d_sub_1_7x7, source_plane_mesh_grid=mesh_grid, ) - over_sampler = aa.OverSamplerUniform(mask=grid_2d_7x7.mask, sub_size=1) - - mapper = aa.Mapper( - mapper_grids=mapper_grids, over_sampler=over_sampler, regularization=None - ) + mapper = aa.Mapper(mapper_grids=mapper_grids, regularization=None) simplex_index_for_sub_slim_index = mapper.delaunay.find_simplex( mapper.source_plane_data_grid diff --git a/test_autoarray/inversion/pixelization/mappers/test_factory.py b/test_autoarray/inversion/pixelization/mappers/test_factory.py index 2fb227716..c08bca937 100644 --- a/test_autoarray/inversion/pixelization/mappers/test_factory.py +++ b/test_autoarray/inversion/pixelization/mappers/test_factory.py @@ -20,23 +20,20 @@ def test__rectangular_mapper(): ) # Slightly manipulate input grid so sub gridding is evidence in first source pixel. - over_sampler = aa.OverSamplerUniform(mask=mask, sub_size=2) - over_sampled_grid = over_sampler.over_sampled_grid - over_sampled_grid[0, 0] = -2.0 - over_sampled_grid[0, 1] = 2.0 + grid = aa.Grid2D.from_mask(mask=mask, over_sample_size=2) + grid.over_sampled[0, 0] = -2.0 + grid.over_sampled[0, 1] = 2.0 mesh = aa.mesh.Rectangular(shape=(3, 3)) mapper_grids = mesh.mapper_grids_from( mask=mask, border_relocator=None, - source_plane_data_grid=over_sampled_grid, + source_plane_data_grid=grid, source_plane_mesh_grid=None, ) - mapper = aa.Mapper( - mapper_grids=mapper_grids, over_sampler=over_sampler, regularization=None - ) + mapper = aa.Mapper(mapper_grids=mapper_grids, regularization=None) assert isinstance(mapper, aa.MapperRectangular) assert mapper.image_plane_mesh_grid == None @@ -73,10 +70,10 @@ def test__delaunay_mapper(): ) # Slightly manipulate input grid so sub gridding is evidence in first source pixel. - over_sampler = aa.OverSamplerUniform(mask=mask, sub_size=2) - over_sampled_grid = over_sampler.over_sampled_grid - over_sampled_grid[0, 0] = -2.0 - over_sampled_grid[0, 1] = 2.0 + grid = aa.Grid2D.from_mask(mask=mask, over_sample_size=2) + + grid.over_sampled[0, 0] = -2.0 + grid.over_sampled[0, 1] = 2.0 mesh = aa.mesh.Delaunay() image_mesh = aa.image_mesh.Overlay(shape=(3, 3)) @@ -87,13 +84,11 @@ def test__delaunay_mapper(): mapper_grids = mesh.mapper_grids_from( mask=mask, border_relocator=None, - source_plane_data_grid=over_sampled_grid, + source_plane_data_grid=grid, source_plane_mesh_grid=image_plane_mesh_grid, ) - mapper = aa.Mapper( - mapper_grids=mapper_grids, over_sampler=over_sampler, regularization=None - ) + mapper = aa.Mapper(mapper_grids=mapper_grids, regularization=None) assert isinstance(mapper, aa.MapperDelaunay) assert (mapper.source_plane_mesh_grid == image_plane_mesh_grid).all() @@ -131,10 +126,10 @@ def test__voronoi_mapper(): ) # Slightly manipulate input grid so sub gridding is evidence in first source pixel. - over_sampler = aa.OverSamplerUniform(mask=mask, sub_size=2) - over_sampled_grid = over_sampler.over_sampled_grid - over_sampled_grid[0, 0] = -2.0 - over_sampled_grid[0, 1] = 2.0 + grid = aa.Grid2D.from_mask(mask=mask, over_sample_size=2) + + grid.over_sampled[0, 0] = -2.0 + grid.over_sampled[0, 1] = 2.0 mesh = aa.mesh.Voronoi() image_mesh = aa.image_mesh.Overlay(shape=(3, 3)) @@ -145,13 +140,11 @@ def test__voronoi_mapper(): mapper_grids = mesh.mapper_grids_from( mask=mask, border_relocator=None, - source_plane_data_grid=over_sampled_grid, + source_plane_data_grid=grid, source_plane_mesh_grid=image_plane_mesh_grid, ) - mapper = aa.Mapper( - mapper_grids=mapper_grids, over_sampler=over_sampler, regularization=None - ) + mapper = aa.Mapper(mapper_grids=mapper_grids, regularization=None) assert (mapper.source_plane_mesh_grid == image_plane_mesh_grid).all() assert mapper.source_plane_mesh_grid.origin == pytest.approx((0.0, 0.0), 1.0e-4) diff --git a/test_autoarray/inversion/pixelization/mappers/test_rectangular.py b/test_autoarray/inversion/pixelization/mappers/test_rectangular.py index 9e948a106..026e230df 100644 --- a/test_autoarray/inversion/pixelization/mappers/test_rectangular.py +++ b/test_autoarray/inversion/pixelization/mappers/test_rectangular.py @@ -18,24 +18,23 @@ def test__pix_indexes_for_sub_slim_index__matches_util(): ], pixel_scales=1.0, shape_native=(3, 3), + over_sample_size=1, ) - mesh_grid = aa.Mesh2DRectangular.overlay_grid(shape_native=(3, 3), grid=grid) + mesh_grid = aa.Mesh2DRectangular.overlay_grid( + shape_native=(3, 3), grid=grid.over_sampled + ) mapper_grids = aa.MapperGrids( mask=grid.mask, source_plane_data_grid=grid, source_plane_mesh_grid=mesh_grid ) - over_sampler = aa.OverSamplerUniform(mask=grid.mask, sub_size=1) - - mapper = aa.Mapper( - mapper_grids=mapper_grids, over_sampler=over_sampler, regularization=None - ) + mapper = aa.Mapper(mapper_grids=mapper_grids, regularization=None) pix_indexes_for_sub_slim_index_util = np.array( [ aa.util.geometry.grid_pixel_indexes_2d_slim_from( - grid_scaled_2d_slim=np.array(grid), + grid_scaled_2d_slim=np.array(grid.over_sampled), shape_native=mesh_grid.shape_native, pixel_scales=mesh_grid.pixel_scales, origin=mesh_grid.origin, @@ -48,21 +47,19 @@ def test__pix_indexes_for_sub_slim_index__matches_util(): ).all() -def test__pixel_signals_from__matches_util(grid_2d_7x7, image_7x7): - mesh_grid = aa.Mesh2DRectangular.overlay_grid(shape_native=(3, 3), grid=grid_2d_7x7) - - over_sampler = aa.OverSamplerUniform(mask=grid_2d_7x7.mask, sub_size=1) +def test__pixel_signals_from__matches_util(grid_2d_sub_1_7x7, image_7x7): + mesh_grid = aa.Mesh2DRectangular.overlay_grid( + shape_native=(3, 3), grid=grid_2d_sub_1_7x7.over_sampled + ) mapper_grids = aa.MapperGrids( - mask=grid_2d_7x7.mask, - source_plane_data_grid=grid_2d_7x7, + mask=grid_2d_sub_1_7x7.mask, + source_plane_data_grid=grid_2d_sub_1_7x7, source_plane_mesh_grid=mesh_grid, adapt_data=image_7x7, ) - mapper = aa.Mapper( - mapper_grids=mapper_grids, over_sampler=over_sampler, regularization=None - ) + mapper = aa.Mapper(mapper_grids=mapper_grids, regularization=None) pixel_signals = mapper.pixel_signals_from(signal_scale=2.0) @@ -72,7 +69,7 @@ def test__pixel_signals_from__matches_util(grid_2d_7x7, image_7x7): pix_indexes_for_sub_slim_index=mapper.pix_indexes_for_sub_slim_index, pix_size_for_sub_slim_index=mapper.pix_sizes_for_sub_slim_index, pixel_weights=mapper.pix_weights_for_sub_slim_index, - slim_index_for_sub_slim_index=over_sampler.slim_for_sub_slim, + slim_index_for_sub_slim_index=grid_2d_sub_1_7x7.over_sampler.slim_for_sub_slim, adapt_data=np.array(image_7x7), ) diff --git a/test_autoarray/inversion/pixelization/mappers/test_voronoi.py b/test_autoarray/inversion/pixelization/mappers/test_voronoi.py index f6a735237..57e971aae 100644 --- a/test_autoarray/inversion/pixelization/mappers/test_voronoi.py +++ b/test_autoarray/inversion/pixelization/mappers/test_voronoi.py @@ -3,26 +3,25 @@ import autoarray as aa -def test__pix_indexes_for_sub_slim_index__matches_util(grid_2d_7x7): +def test__pix_indexes_for_sub_slim_index__matches_util(grid_2d_sub_1_7x7): source_plane_mesh_grid = aa.Grid2D.no_mask( values=[[0.1, 0.1], [1.1, 0.6], [2.1, 0.1], [0.4, 1.1], [1.1, 7.1], [2.1, 1.1]], shape_native=(3, 2), pixel_scales=1.0, + over_sample_size=1, ) source_plane_mesh_grid = aa.Mesh2DVoronoi( values=source_plane_mesh_grid, ) - over_sampler = aa.OverSamplerUniform(mask=grid_2d_7x7.mask, sub_size=1) - source_plane_mesh_grid = aa.Mesh2DVoronoi( values=source_plane_mesh_grid, ) mapper_grids = aa.MapperGrids( - mask=grid_2d_7x7.mask, - source_plane_data_grid=grid_2d_7x7, + mask=grid_2d_sub_1_7x7.mask, + source_plane_data_grid=grid_2d_sub_1_7x7, source_plane_mesh_grid=source_plane_mesh_grid, ) pytest.importorskip( @@ -30,16 +29,14 @@ def test__pix_indexes_for_sub_slim_index__matches_util(grid_2d_7x7): reason="Voronoi C library not installed, see util.nn README.md", ) - mapper = aa.Mapper( - mapper_grids=mapper_grids, over_sampler=over_sampler, regularization=None - ) + mapper = aa.Mapper(mapper_grids=mapper_grids, regularization=None) ( pix_indexes_for_sub_slim_index_util, sizes, weights, ) = aa.util.mapper.pix_size_weights_voronoi_nn_from( - grid=grid_2d_7x7, mesh_grid=source_plane_mesh_grid + grid=grid_2d_sub_1_7x7, mesh_grid=source_plane_mesh_grid ) assert ( diff --git a/test_autoarray/inversion/pixelization/mesh/test_abstract.py b/test_autoarray/inversion/pixelization/mesh/test_abstract.py index ca431fe8b..c00b68178 100644 --- a/test_autoarray/inversion/pixelization/mesh/test_abstract.py +++ b/test_autoarray/inversion/pixelization/mesh/test_abstract.py @@ -29,9 +29,7 @@ def test__grid_is_relocated_via_border(grid_2d_7x7): source_plane_mesh_grid=image_mesh, ) - mapper = aa.Mapper( - mapper_grids=mapper_grids, over_sampler=None, regularization=None - ) + mapper = aa.Mapper(mapper_grids=mapper_grids, regularization=None) assert grid[8, 0] != mapper.source_plane_data_grid[8, 0] assert mapper.source_plane_data_grid[8, 0] < 5.0 @@ -48,9 +46,7 @@ def test__grid_is_relocated_via_border(grid_2d_7x7): source_plane_mesh_grid=image_mesh, ) - mapper = aa.Mapper( - mapper_grids=mapper_grids, over_sampler=None, regularization=None - ) + mapper = aa.Mapper(mapper_grids=mapper_grids, regularization=None) assert isinstance(mapper, aa.MapperDelaunay) assert image_mesh[0, 0] != mapper.source_plane_mesh_grid[0, 0] @@ -67,9 +63,7 @@ def test__grid_is_relocated_via_border(grid_2d_7x7): source_plane_mesh_grid=image_mesh, ) - mapper = aa.Mapper( - mapper_grids=mapper_grids, over_sampler=None, regularization=None - ) + mapper = aa.Mapper(mapper_grids=mapper_grids, regularization=None) assert isinstance(mapper, aa.MapperDelaunay) assert image_mesh[0, 0] != mapper.source_plane_mesh_grid[0, 0] diff --git a/test_autoarray/inversion/pixelization/test_border_relocator.py b/test_autoarray/inversion/pixelization/test_border_relocator.py index eb6cd3301..db1902692 100644 --- a/test_autoarray/inversion/pixelization/test_border_relocator.py +++ b/test_autoarray/inversion/pixelization/test_border_relocator.py @@ -325,11 +325,10 @@ def test__relocated_grid_from__inside_border_no_relocations(): shape_native=(30, 30), radius=1.0, pixel_scales=(0.1, 0.1) ) - over_sampling = aa.OverSamplerUniform( - mask=mask, sub_size=np.array(mask.pixels_in_mask * [2]) + grid = aa.Grid2D.from_mask( + mask=mask, over_sample_size=np.array(mask.pixels_in_mask * [2]) ) - grid = over_sampling.over_sampled_grid - grid[1, :] = [0.1, 0.1] + grid.over_sampled[1, :] = [0.1, 0.1] border_relocator = aa.BorderRelocator( mask=mask, sub_size=np.array(mask.pixels_in_mask * [2]) @@ -337,7 +336,7 @@ def test__relocated_grid_from__inside_border_no_relocations(): relocated_grid = border_relocator.relocated_grid_from(grid=grid) - assert (relocated_grid[1] == np.array([0.1, 0.1])).all() + assert (relocated_grid.over_sampled[1] == np.array([0.1, 0.1])).all() def test__relocated_grid_from__outside_border_includes_relocations(): @@ -345,11 +344,10 @@ def test__relocated_grid_from__outside_border_includes_relocations(): shape_native=(30, 30), radius=1.0, pixel_scales=(0.1, 0.1) ) - over_sampling = aa.OverSamplerUniform( - mask=mask, sub_size=np.array(mask.pixels_in_mask * [2]) + grid = aa.Grid2D.from_mask( + mask=mask, over_sample_size=np.array(mask.pixels_in_mask * [2]) ) - grid = over_sampling.over_sampled_grid - grid[1, :] = [10.1, 0.1] + grid.over_sampled[1, :] = [10.1, 0.1] border_relocator = aa.BorderRelocator( mask=mask, sub_size=np.array(mask.pixels_in_mask * [2]) @@ -357,7 +355,9 @@ def test__relocated_grid_from__outside_border_includes_relocations(): relocated_grid = border_relocator.relocated_grid_from(grid=grid) - assert relocated_grid[1] == pytest.approx([0.97783243, 0.00968151], 1e-4) + assert relocated_grid.over_sampled[1] == pytest.approx( + [0.97783243, 0.00968151], 1e-4 + ) def test__relocated_grid_from__positive_origin_included_in_relocate(): @@ -368,16 +368,11 @@ def test__relocated_grid_from__positive_origin_included_in_relocate(): centre=(1.0, 1.0), ) - over_sampling = aa.OverSamplerUniform( - mask=mask, sub_size=np.array(mask.pixels_in_mask * [2]) - ) - grid = over_sampling.over_sampled_grid - grid[1, :] = [11.1, 1.0] + grid = aa.Grid2D.from_mask(mask=mask, over_sample_size=2) + grid.over_sampled[1, :] = [11.1, 1.0] - border_relocator = aa.BorderRelocator( - mask=mask, sub_size=np.array(mask.pixels_in_mask * [2]) - ) + border_relocator = aa.BorderRelocator(mask=mask, sub_size=grid.over_sample_size) relocated_grid = border_relocator.relocated_grid_from(grid=grid) - assert relocated_grid[1] == pytest.approx([1.97783243, 1.0], 1e-4) + assert relocated_grid.over_sampled[1] == pytest.approx([1.97783243, 1.0], 1e-4) diff --git a/test_autoarray/inversion/plot/test_inversion_plotters.py b/test_autoarray/inversion/plot/test_inversion_plotters.py index 11dee829a..735365ecd 100644 --- a/test_autoarray/inversion/plot/test_inversion_plotters.py +++ b/test_autoarray/inversion/plot/test_inversion_plotters.py @@ -37,13 +37,13 @@ def test__individual_attributes_are_output_for_all_mappers( pixelization_index=0, reconstructed_image=True, reconstruction=True, - errors=True, + reconstruction_noise_map=True, regularization_weights=True, ) assert path.join(plot_path, "reconstructed_image.png") in plot_patch.paths assert path.join(plot_path, "reconstruction.png") in plot_patch.paths - assert path.join(plot_path, "errors.png") in plot_patch.paths + assert path.join(plot_path, "reconstruction_noise_map.png") in plot_patch.paths assert path.join(plot_path, "regularization_weights.png") in plot_patch.paths pytest.importorskip( @@ -66,7 +66,7 @@ def test__individual_attributes_are_output_for_all_mappers( sub_pixels_per_image_pixels=True, mesh_pixels_per_image_pixels=True, image_pixels_per_mesh_pixel=True, - errors=True, + reconstruction_noise_map=True, signal_to_noise_map=True, regularization_weights=True, ) @@ -76,7 +76,7 @@ def test__individual_attributes_are_output_for_all_mappers( assert path.join(plot_path, "sub_pixels_per_image_pixels.png") in plot_patch.paths assert path.join(plot_path, "mesh_pixels_per_image_pixels.png") in plot_patch.paths assert path.join(plot_path, "image_pixels_per_mesh_pixel.png") in plot_patch.paths - assert path.join(plot_path, "errors.png") in plot_patch.paths + assert path.join(plot_path, "reconstruction_noise_map.png") in plot_patch.paths assert path.join(plot_path, "signal_to_noise_map.png") in plot_patch.paths assert path.join(plot_path, "regularization_weights.png") in plot_patch.paths @@ -85,12 +85,12 @@ def test__individual_attributes_are_output_for_all_mappers( inversion_plotter.figures_2d_of_pixelization( pixelization_index=0, reconstructed_image=True, - errors=True, + reconstruction_noise_map=True, ) assert path.join(plot_path, "reconstructed_image.png") in plot_patch.paths assert path.join(plot_path, "reconstruction.png") not in plot_patch.paths - assert path.join(plot_path, "errors.png") in plot_patch.paths + assert path.join(plot_path, "reconstruction_noise_map.png") in plot_patch.paths def test__inversion_subplot_of_mapper__is_output_for_all_inversions( diff --git a/test_autoarray/inversion/plot/test_mapper_plotters.py b/test_autoarray/inversion/plot/test_mapper_plotters.py index 85ac07f89..3dabe5514 100644 --- a/test_autoarray/inversion/plot/test_mapper_plotters.py +++ b/test_autoarray/inversion/plot/test_mapper_plotters.py @@ -67,11 +67,15 @@ def test__get_2d__via_mapper_for_source_from(rectangular_mapper_7x7_3x3): assert mapper_plotter.visuals_2d.origin == None assert get_2d.origin.in_list == [(0.0, 0.0)] - assert (get_2d.grid == rectangular_mapper_7x7_3x3.source_plane_data_grid).all() + assert ( + get_2d.grid == rectangular_mapper_7x7_3x3.source_plane_data_grid.over_sampled + ).all() assert (get_2d.mesh_grid == rectangular_mapper_7x7_3x3.source_plane_mesh_grid).all() - border_grid = rectangular_mapper_7x7_3x3.mapper_grids.source_plane_data_grid[ - rectangular_mapper_7x7_3x3.border_relocator.sub_border_slim - ] + border_grid = ( + rectangular_mapper_7x7_3x3.mapper_grids.source_plane_data_grid.over_sampled[ + rectangular_mapper_7x7_3x3.border_relocator.sub_border_slim + ] + ) assert (get_2d.border == border_grid).all() include = aplt.Include2D( diff --git a/test_autoarray/inversion/test_linear_obj.py b/test_autoarray/inversion/test_linear_obj.py index ab5b447c8..c54204350 100644 --- a/test_autoarray/inversion/test_linear_obj.py +++ b/test_autoarray/inversion/test_linear_obj.py @@ -44,11 +44,7 @@ def test__data_to_pix_unique_from(): mask = aa.Mask2D.all_false(shape_native=(1, 2), pixel_scales=0.1) - over_sampling = aa.OverSamplerUniform(mask=mask, sub_size=2) - - grid = aa.Grid2D.uniform( - shape_native=(1, 2), pixel_scales=0.1, over_sampling=over_sampling - ) + grid = aa.Grid2D.uniform(shape_native=(1, 2), pixel_scales=0.1, over_sample_size=2) linear_obj = aa.AbstractLinearObjFuncList(grid=grid, regularization=None) diff --git a/test_autoarray/operators/over_sample/test_decorator.py b/test_autoarray/operators/over_sample/test_decorator.py index 73fde84a7..47dae1e88 100644 --- a/test_autoarray/operators/over_sample/test_decorator.py +++ b/test_autoarray/operators/over_sample/test_decorator.py @@ -20,15 +20,13 @@ def test__in_grid_2d__over_sample_uniform__out_ndarray_1d(): pixel_scales=(1.0, 1.0), ) - over_sampling = aa.OverSamplingUniform(sub_size=2) + grid_2d = aa.Grid2D.from_mask(mask=mask, over_sample_size=2) - grid_2d = aa.Grid2D.from_mask(mask=mask, over_sampling=over_sampling) - - obj = aa.m.MockGridLikeIteratorObj() + obj = aa.m.MockGrid2DLikeObj() ndarray_1d = obj.ndarray_1d_from(grid=grid_2d) - over_sample_uniform = aa.OverSamplerUniform(mask=mask, sub_size=2) + over_sample_uniform = aa.OverSampler(mask=mask, sub_size=2) mask_sub_2 = aa.util.over_sample.oversample_mask_2d_from( mask=np.array(mask), sub_size=2 @@ -36,7 +34,7 @@ def test__in_grid_2d__over_sample_uniform__out_ndarray_1d(): mask_sub_2 = aa.Mask2D(mask=mask_sub_2, pixel_scales=(0.5, 0.5)) - grid = aa.Grid2D(values=over_sample_uniform.over_sampled_grid, mask=mask_sub_2) + grid = aa.Grid2D(values=grid_2d.over_sampled, mask=mask_sub_2) ndarray_1d_via_grid = obj.ndarray_1d_from(grid=grid) @@ -48,76 +46,3 @@ def test__in_grid_2d__over_sample_uniform__out_ndarray_1d(): assert isinstance(ndarray_1d, aa.Array2D) assert (ndarray_1d == ndarray_1d_via_grid).all() - - -def test__in_grid_2d_over_sample_iterate__out_ndarray_1d__values_use_iteration(): - mask = aa.Mask2D( - mask=[ - [True, True, True, True, True], - [True, False, False, False, True], - [True, False, False, False, True], - [True, False, False, False, True], - [True, True, True, True, True], - ], - pixel_scales=(1.0, 1.0), - origin=(0.001, 0.001), - ) - - over_sampling = aa.OverSamplingIterate(fractional_accuracy=1.0, sub_steps=[2, 3]) - - grid_2d = aa.Grid2D.from_mask(mask=mask, over_sampling=over_sampling) - - obj = aa.m.MockGridLikeIteratorObj() - - ndarray_1d = obj.ndarray_1d_from(grid=grid_2d) - - over_sample_uniform = aa.OverSamplerUniform(mask=mask, sub_size=3) - - values_sub_3 = over_sample_uniform.array_via_func_from( - func=ndarray_1d_from, obj=object - ) - - assert ndarray_1d == pytest.approx(values_sub_3, 1.0e-4) - - grid_2d = aa.Grid2D.from_mask( - mask=mask, - over_sampling=aa.OverSamplingIterate( - fractional_accuracy=0.000001, sub_steps=[2, 4, 8, 16, 32] - ), - ) - - obj = aa.m.MockGridLikeIteratorObj() - - ndarray_1d = obj.ndarray_1d_from(grid=grid_2d) - - over_sample_uniform = aa.OverSamplerUniform(mask=mask, sub_size=2) - - values_sub_2 = over_sample_uniform.array_via_func_from( - func=ndarray_1d_from, obj=object - ) - - assert ndarray_1d == pytest.approx(values_sub_2, 1.0e-4) - - grid_2d = aa.Grid2D.from_mask( - mask=mask, - over_sampling=aa.OverSamplingIterate(fractional_accuracy=0.5, sub_steps=[2, 4]), - ) - - iterate_obj = aa.m.MockGridLikeIteratorObj() - - ndarray_1d = iterate_obj.ndarray_1d_from(grid=grid_2d) - - over_sample_uniform = aa.OverSamplerUniform(mask=mask, sub_size=2) - values_sub_2 = over_sample_uniform.array_via_func_from( - func=ndarray_1d_from, obj=object - ) - over_sample_uniform = aa.OverSamplerUniform(mask=mask, sub_size=4) - values_sub_4 = over_sample_uniform.array_via_func_from( - func=ndarray_1d_from, obj=object - ) - - assert ndarray_1d.native[1, 1] == values_sub_2.native[1, 1] - assert ndarray_1d.native[2, 2] != values_sub_2.native[2, 2] - - assert ndarray_1d.native[1, 1] != values_sub_4.native[1, 1] - assert ndarray_1d.native[2, 2] == values_sub_4.native[2, 2] diff --git a/test_autoarray/operators/over_sample/test_iterate.py b/test_autoarray/operators/over_sample/test_iterate.py deleted file mode 100644 index 7330eaecc..000000000 --- a/test_autoarray/operators/over_sample/test_iterate.py +++ /dev/null @@ -1,230 +0,0 @@ -import numpy as np - -import autoarray as aa - -from autoarray.structures.mock.mock_decorators import ( - ndarray_1d_from, - ndarray_1d_zeros_from, - ndarray_2d_from, -) - - -def test__threshold_mask_from(): - mask = aa.Mask2D( - mask=[ - [True, True, True, True], - [True, False, False, True], - [True, False, False, True], - [True, True, True, True], - ], - pixel_scales=(1.0, 1.0), - ) - - arr = aa.Array2D( - values=[ - [0.0, 0.0, 0.0, 0.0], - [0.0, 1.0, 1.0, 0.0], - [0.0, 1.0, 1.0, 0.0], - [0.0, 0.0, 0.0, 0.0], - ], - mask=mask, - ) - - over_sampling = aa.OverSamplerIterate(mask=mask, fractional_accuracy=0.9999) - - threshold_mask = over_sampling.threshold_mask_from( - array_lower_sub_2d=arr.native, array_higher_sub_2d=arr.native - ) - - assert ( - threshold_mask - == np.array( - [ - [True, True, True, True], - [True, True, True, True], - [True, True, True, True], - [True, True, True, True], - ] - ) - ).all() - - mask_lower_sub = aa.Mask2D( - mask=[ - [True, True, True, True], - [True, False, False, True], - [True, False, False, True], - [True, True, True, True], - ], - pixel_scales=(1.0, 1.0), - ) - - mask_higher_sub = aa.Mask2D( - mask=[ - [True, True, True, True], - [True, False, True, True], - [True, False, False, True], - [True, True, True, True], - ], - pixel_scales=(1.0, 1.0), - ) - - over_sampling = aa.OverSamplerIterate(mask=mask, fractional_accuracy=0.5) - - array_lower_sub = aa.Array2D( - [ - [0.0, 0.0, 0.0, 0.0], - [0.0, 2.0, 2.0, 0.0], - [0.0, 2.0, 2.0, 0.0], - [0.0, 0.0, 0.0, 0.0], - ], - mask=mask_lower_sub, - ) - - array_higher_sub = aa.Array2D( - [ - [0.0, 0.0, 0.0, 0.0], - [0.0, 5.0, 5.0, 0.0], - [0.0, 5.0, 5.0, 0.0], - [0.0, 0.0, 0.0, 0.0], - ], - mask=mask_higher_sub, - ) - - threshold_mask = over_sampling.threshold_mask_from( - array_lower_sub_2d=array_lower_sub.native, - array_higher_sub_2d=array_higher_sub.native, - ) - - assert ( - threshold_mask - == np.array( - [ - [True, True, True, True], - [True, False, True, True], - [True, False, False, True], - [True, True, True, True], - ] - ) - ).all() - - -def test__array_via_func_from__extreme_fractional_accuracies_uses_last_or_first_sub(): - mask = aa.Mask2D( - mask=[ - [True, True, True, True, True], - [True, False, False, False, True], - [True, False, False, False, True], - [True, False, False, False, True], - [True, True, True, True, True], - ], - pixel_scales=(1.0, 1.0), - origin=(0.001, 0.001), - ) - - over_sampling = aa.OverSamplerIterate( - mask=mask, fractional_accuracy=1.0, sub_steps=[2, 3] - ) - - values = over_sampling.array_via_func_from( - func=ndarray_1d_from, - obj=None, - ) - - over_sample_uniform = aa.OverSamplerUniform(mask=mask, sub_size=3) - - values_sub_3 = over_sample_uniform.array_via_func_from( - func=ndarray_1d_from, obj=object - ) - - assert (values == values_sub_3).all() - - # This test ensures that if the fractional accuracy is met on the last sub_size jump (e.g. 2 doesnt meet it, - # but 3 does) that the sub_size of 3 is used. There was a bug where the mask was not updated correctly and the - # iterated array double counted the values. - - values = over_sampling.array_via_func_from( - func=ndarray_1d_from, - obj=None, - ) - - assert (values == values_sub_3).all() - - over_sampling = aa.OverSamplerIterate( - mask=mask, fractional_accuracy=0.000001, sub_steps=[2, 4, 8, 16, 32] - ) - - values = over_sampling.array_via_func_from( - func=ndarray_1d_from, - obj=None, - ) - - over_sample_uniform = aa.OverSamplerUniform(mask=mask, sub_size=2) - - values_sub_2 = over_sample_uniform.array_via_func_from( - func=ndarray_1d_from, obj=object - ) - - assert (values == values_sub_2).all() - - -def test__array_via_func_from__check_values_computed_to_fractional_accuracy(): - mask = aa.Mask2D( - mask=[ - [True, True, True, True, True], - [True, False, False, False, True], - [True, False, False, False, True], - [True, False, False, False, True], - [True, True, True, True, True], - ], - pixel_scales=(1.0, 1.0), - origin=(0.001, 0.001), - ) - - over_sampling = aa.OverSamplerIterate( - mask=mask, fractional_accuracy=0.5, sub_steps=[2, 4] - ) - - values = over_sampling.array_via_func_from( - func=ndarray_1d_from, - obj=None, - ) - - over_sample_uniform = aa.OverSamplerUniform(mask=mask, sub_size=2) - values_sub_2 = over_sample_uniform.array_via_func_from( - func=ndarray_1d_from, obj=object - ) - over_sample_uniform = aa.OverSamplerUniform(mask=mask, sub_size=4) - values_sub_4 = over_sample_uniform.array_via_func_from( - func=ndarray_1d_from, obj=object - ) - - assert values.native[1, 1] == values_sub_2.native[1, 1] - assert values.native[2, 2] != values_sub_2.native[2, 2] - - assert values.native[1, 1] != values_sub_4.native[1, 1] - assert values.native[2, 2] == values_sub_4.native[2, 2] - - -def test__array_via_func_from__func_returns_all_zeros__iteration_terminated(): - mask = aa.Mask2D( - mask=[ - [True, True, True, True, True], - [True, False, False, False, True], - [True, False, False, False, True], - [True, False, False, False, True], - [True, True, True, True, True], - ], - pixel_scales=(1.0, 1.0), - origin=(0.001, 0.001), - ) - - over_sampling = aa.OverSamplerIterate( - mask=mask, fractional_accuracy=1.0, sub_steps=[2, 3] - ) - - values = over_sampling.array_via_func_from( - func=ndarray_1d_zeros_from, - obj=None, - ) - - assert (values == np.zeros((9,))).all() diff --git a/test_autoarray/operators/over_sample/test_over_sample_util.py b/test_autoarray/operators/over_sample/test_over_sample_util.py index 017900c5b..c0552d8a0 100644 --- a/test_autoarray/operators/over_sample/test_over_sample_util.py +++ b/test_autoarray/operators/over_sample/test_over_sample_util.py @@ -367,3 +367,81 @@ def test__grid_2d_slim_over_sampled_via_mask_from(): assert grid[0:4] == pytest.approx( np.array([[1.75, -0.5], [1.75, 2.5], [0.25, -0.5], [0.25, 2.5]]), 1e-4 ) + + +def test__from_manual_adapt_radial_bin(): + mask = aa.Mask2D.circular(shape_native=(5, 5), pixel_scales=2.0, radius=3.0) + + grid = aa.Grid2D.from_mask(mask=mask) + + sub_size = aa.util.over_sample.over_sample_size_via_radial_bins_from( + grid=grid, sub_size_list=[8, 4, 2], radial_list=[1.5, 2.5] + ) + assert sub_size.native == pytest.approx( + np.array( + [ + [0, 0, 0, 0, 0], + [0, 2, 4, 2, 0], + [0, 4, 8, 4, 0], + [0, 2, 4, 2, 0], + [0, 0, 0, 0, 0], + ] + ), + 1.0e-4, + ) + + +def test__from_manual_adapt_radial_bin__centre_list_input(): + mask = aa.Mask2D.circular(shape_native=(5, 5), pixel_scales=2.0, radius=3.0) + + grid = aa.Grid2D.from_mask(mask=mask) + + sub_size = aa.util.over_sample.over_sample_size_via_radial_bins_from( + grid=grid, + sub_size_list=[8, 4, 2], + radial_list=[1.5, 2.5], + centre_list=[(0.0, -2.0), (0.0, 2.0)], + ) + + assert sub_size.native == pytest.approx( + np.array( + [ + [0, 0, 0, 0, 0], + [0, 4, 2, 4, 0], + [0, 8, 4, 8, 0], + [0, 4, 2, 4, 0], + [0, 0, 0, 0, 0], + ] + ), + 1.0e-4, + ) + + +def test__from_adapt(): + mask = aa.Mask2D( + mask=[[True, True, True], [True, False, False], [True, True, False]], + pixel_scales=1.0, + ) + + data = aa.Array2D(values=[1.0, 2.0, 3.0], mask=mask) + noise_map = aa.Array2D(values=[1.0, 2.0, 1.0], mask=mask) + + sub_size = aa.util.over_sample.over_sample_size_via_adapt_from( + data=data, + noise_map=noise_map, + signal_to_noise_cut=1.5, + sub_size_lower=2, + sub_size_upper=4, + ) + + assert sub_size == pytest.approx([2, 2, 4], 1.0e-4) + + sub_size = aa.util.over_sample.over_sample_size_via_adapt_from( + data=data, + noise_map=noise_map, + signal_to_noise_cut=0.5, + sub_size_lower=2, + sub_size_upper=4, + ) + + assert sub_size == pytest.approx([4, 4, 4], 1.0e-4) diff --git a/test_autoarray/operators/over_sample/test_uniform.py b/test_autoarray/operators/over_sample/test_over_sampler.py similarity index 51% rename from test_autoarray/operators/over_sample/test_uniform.py rename to test_autoarray/operators/over_sample/test_over_sampler.py index 06d7cc67e..7da24d79a 100644 --- a/test_autoarray/operators/over_sample/test_uniform.py +++ b/test_autoarray/operators/over_sample/test_over_sampler.py @@ -30,7 +30,7 @@ def test__from_sub_size_int(): pixel_scales=1.0, ) - over_sampling = aa.OverSamplerUniform(mask=mask, sub_size=2) + over_sampling = aa.OverSampler(mask=mask, sub_size=2) assert over_sampling.sub_size.slim == pytest.approx([2, 2, 2], 1.0e-4) assert over_sampling.sub_size.native == pytest.approx( @@ -44,127 +44,33 @@ def test__sub_pixel_areas(): pixel_scales=1.0, ) - over_sampling = aa.OverSamplerUniform(mask=mask, sub_size=np.array([1, 2])) + over_sampling = aa.OverSampler(mask=mask, sub_size=np.array([1, 2])) areas = over_sampling.sub_pixel_areas assert areas == pytest.approx([1.0, 0.25, 0.25, 0.25, 0.25], 1.0e-4) -def test__from_manual_adapt_radial_bin(): - mask = aa.Mask2D.circular(shape_native=(5, 5), pixel_scales=2.0, radius=3.0) - - grid = aa.Grid2D.from_mask(mask=mask) - - over_sampling = aa.OverSamplingUniform.from_radial_bins( - grid=grid, sub_size_list=[8, 4, 2], radial_list=[1.5, 2.5] - ) - assert over_sampling.sub_size.native == pytest.approx( - np.array( - [ - [0, 0, 0, 0, 0], - [0, 2, 4, 2, 0], - [0, 4, 8, 4, 0], - [0, 2, 4, 2, 0], - [0, 0, 0, 0, 0], - ] - ), - 1.0e-4, - ) - - -def test__from_manual_adapt_radial_bin__centre_list_input(): - mask = aa.Mask2D.circular(shape_native=(5, 5), pixel_scales=2.0, radius=3.0) - - grid = aa.Grid2D.from_mask(mask=mask) - - over_sampling = aa.OverSamplingUniform.from_radial_bins( - grid=grid, - sub_size_list=[8, 4, 2], - radial_list=[1.5, 2.5], - centre_list=[(0.0, -2.0), (0.0, 2.0)], - ) - - assert over_sampling.sub_size.native == pytest.approx( - np.array( - [ - [0, 0, 0, 0, 0], - [0, 4, 2, 4, 0], - [0, 8, 4, 8, 0], - [0, 4, 2, 4, 0], - [0, 0, 0, 0, 0], - ] - ), - 1.0e-4, - ) - - -def test__from_adapt(): - mask = aa.Mask2D( - mask=[[True, True, True], [True, False, False], [True, True, False]], - pixel_scales=1.0, - ) - - data = aa.Array2D(values=[1.0, 2.0, 3.0], mask=mask) - noise_map = aa.Array2D(values=[1.0, 2.0, 1.0], mask=mask) - - over_sampling = aa.OverSamplingUniform.from_adapt( - data=data, - noise_map=noise_map, - signal_to_noise_cut=1.5, - sub_size_lower=2, - sub_size_upper=4, - ) - - assert over_sampling.sub_size == pytest.approx([2, 2, 4], 1.0e-4) - - over_sampling = aa.OverSamplingUniform.from_adapt( - data=data, - noise_map=noise_map, - signal_to_noise_cut=0.5, - sub_size_lower=2, - sub_size_upper=4, - ) - - assert over_sampling.sub_size == pytest.approx([4, 4, 4], 1.0e-4) - - def test__sub_fraction(): mask = aa.Mask2D( mask=[[False, False], [True, True]], pixel_scales=1.0, ) - over_sampling = aa.OverSamplerUniform( + over_sampling = aa.OverSampler( mask=mask, sub_size=aa.Array2D(values=[1, 2], mask=mask) ) assert over_sampling.sub_fraction.slim == pytest.approx([1.0, 0.25], 1.0e-4) -def test__over_sampled_grid(): - mask = aa.Mask2D( - mask=[[False, False], [True, True]], - pixel_scales=1.0, - ) - - over_sampling = aa.OverSamplerUniform( - mask=mask, sub_size=aa.Array2D(values=[1, 2], mask=mask) - ) - - assert over_sampling.over_sampled_grid.native == pytest.approx( - np.array([[0.5, -0.5], [0.75, 0.25], [0.75, 0.75], [0.25, 0.25], [0.25, 0.75]]), - 1.0e-4, - ) - - def test__binned_array_2d_from(): mask = aa.Mask2D( mask=[[False, False], [True, True]], pixel_scales=1.0, ) - over_sampling = aa.OverSamplerUniform( + over_sampling = aa.OverSampler( mask=mask, sub_size=aa.Array2D(values=[1, 2], mask=mask) ) @@ -182,7 +88,7 @@ def test__sub_mask_index_for_sub_mask_1d_index(): sub_size=2, ) - over_sampling = aa.OverSamplerUniform(mask=mask, sub_size=2) + over_sampling = aa.OverSampler(mask=mask, sub_size=2) sub_mask_index_for_sub_mask_1d_index = ( aa.util.over_sample.native_sub_index_for_slim_sub_index_2d_from( @@ -202,7 +108,7 @@ def test__slim_index_for_sub_slim_index(): sub_size=2, ) - over_sampling = aa.OverSamplerUniform(mask=mask, sub_size=2) + over_sampling = aa.OverSampler(mask=mask, sub_size=2) slim_index_for_sub_slim_index_util = ( aa.util.over_sample.slim_index_for_sub_slim_index_via_mask_2d_from( diff --git a/test_autoarray/plot/get_visuals/test_two_d.py b/test_autoarray/plot/get_visuals/test_two_d.py index c74a2ff0e..5494bcf22 100644 --- a/test_autoarray/plot/get_visuals/test_two_d.py +++ b/test_autoarray/plot/get_visuals/test_two_d.py @@ -107,11 +107,14 @@ def test__via_mapper_for_source_from(rectangular_mapper_7x7_3x3): assert visuals_2d.origin == (1.0, 1.0) assert ( - visuals_2d_via.grid == rectangular_mapper_7x7_3x3.source_plane_data_grid + visuals_2d_via.grid + == rectangular_mapper_7x7_3x3.source_plane_data_grid.over_sampled ).all() - border_grid = rectangular_mapper_7x7_3x3.mapper_grids.source_plane_data_grid[ - rectangular_mapper_7x7_3x3.border_relocator.sub_border_slim - ] + border_grid = ( + rectangular_mapper_7x7_3x3.mapper_grids.source_plane_data_grid.over_sampled[ + rectangular_mapper_7x7_3x3.border_relocator.sub_border_slim + ] + ) assert (visuals_2d_via.border == border_grid).all() assert ( visuals_2d_via.mesh_grid == rectangular_mapper_7x7_3x3.source_plane_mesh_grid diff --git a/test_autoarray/plot/wrap/two_d/test_delaunay_drawer.py b/test_autoarray/plot/wrap/two_d/test_delaunay_drawer.py new file mode 100644 index 000000000..f7f9e4521 --- /dev/null +++ b/test_autoarray/plot/wrap/two_d/test_delaunay_drawer.py @@ -0,0 +1,26 @@ +import autoarray.plot as aplt + +import numpy as np + + +def test__draws_delaunay_pixels_for_sensible_input(delaunay_mapper_9_3x3): + delaunay_drawer = aplt.DelaunayDrawer(linewidth=0.5, edgecolor="r", alpha=1.0) + + delaunay_drawer.draw_delaunay_pixels( + mapper=delaunay_mapper_9_3x3, + pixel_values=np.ones(9), + units=aplt.Units(), + cmap=aplt.Cmap(), + colorbar=None, + ) + + values = np.ones(9) + values[0] = 0.0 + + delaunay_drawer.draw_delaunay_pixels( + mapper=delaunay_mapper_9_3x3, + pixel_values=values, + units=aplt.Units(), + cmap=aplt.Cmap(), + colorbar=aplt.Colorbar(fraction=0.1, pad=0.05), + ) diff --git a/test_autoarray/structures/arrays/test_uniform_2d.py b/test_autoarray/structures/arrays/test_uniform_2d.py index efe3b2136..040fcfae9 100644 --- a/test_autoarray/structures/arrays/test_uniform_2d.py +++ b/test_autoarray/structures/arrays/test_uniform_2d.py @@ -565,6 +565,98 @@ def test__binned_across_columns(): assert (array.binned_across_columns == np.array([1.5, 6.0, 9.0])).all() +def test__brightest_pixel_in_region_from(): + mask = aa.Mask2D.all_false(shape_native=(4, 4), pixel_scales=0.1) + array_2d = aa.Array2D( + values=[ + [1.0, 2.0, 3.0, 4.0], + [6.0, 7.0, 8.0, 9.0], + [11.0, 12.0, 13.0, 15.0], + [16.0, 17.0, 18.0, 20.0], + ], + mask=mask, + ) + + brightest_coordinate = array_2d.brightest_coordinate_in_region_from( + region=(-0.26, 0.06, -0.06, 0.06) + ) + + assert brightest_coordinate == pytest.approx((-0.15, 0.05), 1.0e-4) + + brightest_coordinate = array_2d.brightest_coordinate_in_region_from( + region=(-0.051, 0.151, -0.151, 0.149) + ) + + assert brightest_coordinate == pytest.approx((-0.05, 0.05), 1.0e-4) + + mask = aa.Mask2D.all_false(shape_native=(5, 5), pixel_scales=0.1) + array_2d = aa.Array2D( + values=[ + [1.0, 2.0, 3.0, 4.0, 5.0], + [6.0, 7.0, 8.0, 9.0, 10.0], + [11.0, 12.0, 13.0, 14.0, 15.0], + [16.0, 17.0, 18.0, 19.0, 20.0], + [21.0, 22.0, 23.0, 24.0, 25.0], + ], + mask=mask, + ) + + brightest_coordinate = array_2d.brightest_coordinate_in_region_from( + region=(-0.15, 0.15, -0.15, 0.15) + ) + + assert brightest_coordinate == pytest.approx((-0.1, 0.1), 1.0e-4) + + brightest_coordinate = array_2d.brightest_coordinate_in_region_from( + region=(-0.25, 0.15, -0.15, 0.15) + ) + + assert brightest_coordinate == pytest.approx((-0.2, 0.1), 1.0e-4) + + brightest_coordinate = array_2d.brightest_coordinate_in_region_from( + region=(-0.15, 0.15, -0.15, 0.25) + ) + + assert brightest_coordinate == pytest.approx((-0.1, 0.2), 1.0e-4) + + +def test__brightest_sub_pixel_in_region_from(): + mask = aa.Mask2D.all_false(shape_native=(4, 4), pixel_scales=0.1) + array_2d = aa.Array2D( + values=[ + [1.0, 2.0, 3.0, 4.0], + [6.0, 7.0, 8.0, 9.0], + [11.0, 12.0, 13.0, 15.0], + [16.0, 17.0, 18.0, 20.0], + ], + mask=mask, + ) + + brightest_coordinate = array_2d.brightest_sub_pixel_coordinate_in_region_from( + region=(-0.26, 0.06, -0.06, 0.06) + ) + + assert brightest_coordinate == pytest.approx((-0.1078947, 0.056315), 1.0e-4) + + mask = aa.Mask2D.all_false(shape_native=(5, 5), pixel_scales=0.1) + array_2d = aa.Array2D( + values=[ + [1.0, 2.0, 3.0, 4.0, 5.0], + [6.0, 7.0, 8.0, 9.0, 10.0], + [11.0, 12.0, 13.0, 14.0, 15.0], + [16.0, 17.0, 18.0, 19.0, 20.0], + [21.0, 22.0, 23.0, 24.0, 25.0], + ], + mask=mask, + ) + + brightest_coordinate = array_2d.brightest_sub_pixel_coordinate_in_region_from( + region=(-0.15, 0.15, -0.15, 0.15) + ) + + assert brightest_coordinate == pytest.approx((-0.11754, 0.103508), 1.0e-4) + + def test__header__modified_julian_date(): header_sci_obj = {"DATE-OBS": "2000-01-01", "TIME-OBS": "00:00:00"} diff --git a/test_autoarray/structures/grids/test_irregular_2d.py b/test_autoarray/structures/grids/test_irregular_2d.py index 8039624cf..245848224 100644 --- a/test_autoarray/structures/grids/test_irregular_2d.py +++ b/test_autoarray/structures/grids/test_irregular_2d.py @@ -112,88 +112,3 @@ def test__grid_of_closest_from(): assert ( grid_of_closest == np.array([[0.0, 0.0], [0.0, 0.0], [0.0, 1.0], [0.0, 0.0]]) ).all() - - -def test__uniform__from_grid_sparse_uniform_upscale(): - grid_sparse_uniform = aa.Grid2DIrregularUniform( - values=[[(1.0, 1.0), (1.0, 3.0)]], pixel_scales=2.0 - ) - - grid_upscale = aa.Grid2DIrregularUniform.from_grid_sparse_uniform_upscale( - grid_sparse_uniform=grid_sparse_uniform, upscale_factor=2, pixel_scales=2.0 - ) - - assert ( - grid_upscale - == np.array( - [ - [1.5, 0.5], - [1.5, 1.5], - [0.5, 0.5], - [0.5, 1.5], - [1.5, 2.5], - [1.5, 3.5], - [0.5, 2.5], - [0.5, 3.5], - ] - ) - ).all() - - grid_sparse_uniform = aa.Grid2DIrregularUniform( - values=[[(1.0, 1.0), (1.0, 3.0), (1.0, 5.0), (3.0, 3.0)]], pixel_scales=2.0 - ) - - grid_upscale = aa.Grid2DIrregularUniform.from_grid_sparse_uniform_upscale( - grid_sparse_uniform=grid_sparse_uniform, upscale_factor=4, pixel_scales=2.0 - ) - - grid_upscale_util = aa.util.grid_2d.grid_2d_slim_upscaled_from( - grid_slim=np.array(grid_sparse_uniform), - upscale_factor=4, - pixel_scales=(2.0, 2.0), - ) - - assert (grid_upscale == grid_upscale_util).all() - - -def test__uniform__grid_2d_via_deflection_grid_from(): - grid = aa.Grid2DIrregularUniform( - values=[(1.0, 1.0), (2.0, 2.0)], pixel_scales=(1.0, 1.0), shape_native=(3, 3) - ) - - grid = grid.grid_2d_via_deflection_grid_from( - deflection_grid=np.array([[1.0, 0.0], [1.0, 1.0]]) - ) - - assert type(grid) == aa.Grid2DIrregularUniform - assert grid.in_list == [(0.0, 1.0), (1.0, 1.0)] - assert grid.pixel_scales == (1.0, 1.0) - assert grid.shape_native == (3, 3) - - -def test__uniform__grid_from(): - grid = aa.Grid2DIrregularUniform( - values=[(1.0, 1.0), (2.0, 2.0)], pixel_scales=(1.0, 1.0), shape_native=(3, 3) - ) - - grid_from_1d = grid.grid_from(grid_slim=np.array([[1.0, 1.0], [2.0, 2.0]])) - - assert type(grid_from_1d) == aa.Grid2DIrregularUniform - assert grid_from_1d.in_list == [(1.0, 1.0), (2.0, 2.0)] - assert grid.pixel_scales == (1.0, 1.0) - assert grid.shape_native == (3, 3) - - grid = aa.Grid2DIrregularUniform( - values=[(1.0, 1.0), (2.0, 2.0), (3.0, 3.0)], - pixel_scales=(1.0, 3.0), - shape_native=(5, 5), - ) - - grid_from_1d = grid.grid_from( - grid_slim=np.array([[1.0, 1.0], [2.0, 2.0], [3.0, 3.0]]) - ) - - assert type(grid_from_1d) == aa.Grid2DIrregularUniform - assert grid_from_1d.in_list == [(1.0, 1.0), (2.0, 2.0), (3.0, 3.0)] - assert grid.pixel_scales == (1.0, 3.0) - assert grid.shape_native == (5, 5) diff --git a/test_autoarray/structures/grids/test_uniform_2d.py b/test_autoarray/structures/grids/test_uniform_2d.py index c81f701ea..4084908ba 100644 --- a/test_autoarray/structures/grids/test_uniform_2d.py +++ b/test_autoarray/structures/grids/test_uniform_2d.py @@ -665,17 +665,6 @@ def test__padded_grid_from(): assert padded_grid.shape == (42, 2) assert (padded_grid == padded_grid_util).all() - grid_2d = aa.Grid2D.uniform(shape_native=(5, 5), pixel_scales=8.0) - - padded_grid = grid_2d.padded_grid_from(kernel_shape_native=(2, 5)) - - padded_grid_util = aa.util.grid_2d.grid_2d_slim_via_mask_from( - mask_2d=np.full((6, 9), False), pixel_scales=(8.0, 8.0) - ) - - assert padded_grid.shape == (54, 2) - assert (padded_grid == padded_grid_util).all() - def test__squared_distances_to_coordinate_from(): mask = aa.Mask2D( @@ -851,3 +840,11 @@ def test__is_uniform(): ) assert grid_2d.is_uniform == False + + +def test__apply_over_sampling(): + grid = aa.Grid2D.uniform(shape_native=(2, 2), pixel_scales=1.0, over_sample_size=1) + + grid = grid.apply_over_sampling(over_sample_size=2) + + assert grid.over_sampled.shape[0] == 16 diff --git a/test_autoarray/structures/triangles/coordinate/__init__.py b/test_autoarray/structures/triangles/coordinate/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/test_autoarray/structures/triangles/coordinate/conftest.py b/test_autoarray/structures/triangles/coordinate/conftest.py new file mode 100644 index 000000000..302b565f7 --- /dev/null +++ b/test_autoarray/structures/triangles/coordinate/conftest.py @@ -0,0 +1,21 @@ +import pytest + +import numpy as np + +from autoarray.structures.triangles.coordinate_array import CoordinateArrayTriangles + + +@pytest.fixture +def one_triangle(): + return CoordinateArrayTriangles( + coordinates=np.array([[0, 0]]), + side_length=1.0, + ) + + +@pytest.fixture +def two_triangles(): + return CoordinateArrayTriangles( + coordinates=np.array([[0, 0], [1, 0]]), + side_length=1.0, + ) diff --git a/test_autoarray/structures/triangles/test_coordinate_implementation.py b/test_autoarray/structures/triangles/coordinate/test_coordinate_implementation.py similarity index 87% rename from test_autoarray/structures/triangles/test_coordinate_implementation.py rename to test_autoarray/structures/triangles/coordinate/test_coordinate_implementation.py index d7c7a3388..545d16da0 100644 --- a/test_autoarray/structures/triangles/test_coordinate_implementation.py +++ b/test_autoarray/structures/triangles/coordinate/test_coordinate_implementation.py @@ -4,14 +4,7 @@ from autoarray.structures.triangles.abstract import HEIGHT_FACTOR from autoarray.structures.triangles.coordinate_array import CoordinateArrayTriangles - - -@pytest.fixture -def two_triangles(): - return CoordinateArrayTriangles( - coordinates=np.array([[0, 0], [1, 0]]), - side_length=1.0, - ) +from autoarray.structures.triangles.shape import Point def test_two(two_triangles): @@ -33,14 +26,6 @@ def test_two(two_triangles): ) -@pytest.fixture -def one_triangle(): - return CoordinateArrayTriangles( - coordinates=np.array([[0, 0]]), - side_length=1.0, - ) - - def test_trivial_triangles(one_triangle): assert one_triangle.flip_array == np.array([1]) assert np.all(one_triangle.centres == np.array([[0, 0]])) @@ -257,15 +242,49 @@ def test_for_indexes(two_triangles): ) -def test_for_limits_and_scale(): +def test_means(one_triangle): + assert np.all(one_triangle.means == [[0.0, -0.14433756729740643]]) + + +@pytest.mark.parametrize( + "x, y", + [ + (0.0, 0.0), + (-0.5, -HEIGHT_FACTOR / 2), + (0.5, -HEIGHT_FACTOR / 2), + (0.0, HEIGHT_FACTOR / 2), + ], +) +def test_containment(one_triangle, x, y): + assert one_triangle.containing_indices(Point(x, y)) == [0] + + +def test_triangles_touch(): + triangles = CoordinateArrayTriangles( + np.array([[0, 0], [2, 0]]), + ) + + assert max(triangles.triangles[0][:, 0]) == min(triangles.triangles[1][:, 0]) + + triangles = CoordinateArrayTriangles( + np.array([[0, 0], [0, 1]]), + ) + assert max(triangles.triangles[0][:, 1]) == min(triangles.triangles[1][:, 1]) + + +def test_from_grid_regression(): triangles = CoordinateArrayTriangles.for_limits_and_scale( - x_min=-1.0, - x_max=1.0, - y_min=-1.0, - y_max=1.0, + x_min=-4.75, + x_max=4.75, + y_min=-4.75, + y_max=4.75, + scale=0.5, ) - assert triangles.triangles.shape == (4, 3, 2) + x = triangles.vertices[:, 0] + assert min(x) <= -4.75 + assert max(x) >= 4.75 -def test_means(one_triangle): - assert np.all(one_triangle.means == [[0.0, -0.14433756729740643]]) + y = triangles.vertices[:, 1] + assert min(y) <= -4.75 + assert max(y) >= 4.75 diff --git a/test_autoarray/structures/triangles/test_coordinate_jax.py b/test_autoarray/structures/triangles/coordinate/test_coordinate_jax.py similarity index 87% rename from test_autoarray/structures/triangles/test_coordinate_jax.py rename to test_autoarray/structures/triangles/coordinate/test_coordinate_jax.py index b22701371..1f37a1c90 100644 --- a/test_autoarray/structures/triangles/test_coordinate_jax.py +++ b/test_autoarray/structures/triangles/coordinate/test_coordinate_jax.py @@ -9,7 +9,7 @@ import jax jax.config.update("jax_log_compiles", True) - from autoarray.structures.triangles.jax_coordinate_array import ( + from autoarray.structures.triangles.coordinate_array.jax_coordinate_array import ( CoordinateArrayTriangles, ) except ImportError: @@ -34,21 +34,21 @@ def full_routine(triangles): return up_sampled.for_indexes(indexes) -def test_full_routine(one_triangle, compare_with_nans): - result = full_routine(one_triangle) - - assert compare_with_nans( - result.triangles, - np.array( - [ - [ - [0.0, 0.4330126941204071], - [0.25, 0.0], - [-0.25, 0.0], - ] - ] - ), - ) +# def test_full_routine(one_triangle, compare_with_nans): +# result = full_routine(one_triangle) +# +# assert compare_with_nans( +# result.triangles, +# np.array( +# [ +# [ +# [0.0, 0.4330126941204071], +# [0.25, 0.0], +# [-0.25, 0.0], +# ] +# ] +# ), +# ) def test_neighborhood(one_triangle): diff --git a/test_autoarray/structures/triangles/coordinate/test_vertex_coordinates.py b/test_autoarray/structures/triangles/coordinate/test_vertex_coordinates.py new file mode 100644 index 000000000..88b7910c7 --- /dev/null +++ b/test_autoarray/structures/triangles/coordinate/test_vertex_coordinates.py @@ -0,0 +1,6 @@ +def test_vertex_coordinates(one_triangle): + print(one_triangle.vertex_coordinates) + + +def test_vertex_coordinates_two_triangles(two_triangles): + print(two_triangles.vertex_coordinates) diff --git a/test_autoarray/structures/triangles/test_jax.py b/test_autoarray/structures/triangles/test_jax.py index c934bada0..def239849 100644 --- a/test_autoarray/structures/triangles/test_jax.py +++ b/test_autoarray/structures/triangles/test_jax.py @@ -5,7 +5,7 @@ import jax jax.config.update("jax_log_compiles", True) - from autoarray.structures.triangles.jax_array import ArrayTriangles + from autoarray.structures.triangles.array.jax_array import ArrayTriangles except ImportError: import numpy as np from autoarray.structures.triangles.array import ArrayTriangles diff --git a/test_autoarray/structures/triangles/test_nan_triangles.py b/test_autoarray/structures/triangles/test_nan_triangles.py index 1c9fee576..725cf5257 100644 --- a/test_autoarray/structures/triangles/test_nan_triangles.py +++ b/test_autoarray/structures/triangles/test_nan_triangles.py @@ -2,7 +2,7 @@ try: from jax import numpy as np - from autoarray.structures.triangles.jax_array import ArrayTriangles + from autoarray.structures.triangles.array.jax_array import ArrayTriangles except ImportError: import numpy as np from autoarray.structures.triangles.array import ArrayTriangles