From 9e66b85e3ecffc4d7e7ba30bc129c8f47940b5a4 Mon Sep 17 00:00:00 2001 From: James Nightingale Date: Wed, 22 Mar 2023 14:10:34 +0000 Subject: [PATCH 1/5] added method to compute brightness zeroth regularization weights --- .../regularization/adaptive_brightness.py | 4 +- .../adaptive_brightness_split.py | 4 +- .../regularization/brightness_zeroth.py | 89 +++++++++++++++++++ .../regularization/regularization_util.py | 47 ++++++++-- .../inversion/inversion/test_factory.py | 4 +- 5 files changed, 136 insertions(+), 12 deletions(-) create mode 100644 autoarray/inversion/regularization/brightness_zeroth.py diff --git a/autoarray/inversion/regularization/adaptive_brightness.py b/autoarray/inversion/regularization/adaptive_brightness.py index ff8794a31..9c10e2c2e 100644 --- a/autoarray/inversion/regularization/adaptive_brightness.py +++ b/autoarray/inversion/regularization/adaptive_brightness.py @@ -25,13 +25,13 @@ def __init__( different regions of a pixelization's mesh require different levels of regularization (e.g., high smoothing where the no signal is present and less smoothing where it is, see (Nightingale, Dye and Massey 2018)). - Unlike the instance regularization_matrix scheme, neighboring pixels must now be regularized with one another + Unlike ``Constant`` regularization, neighboring pixels must now be regularized with one another in both directions (e.g. if pixel 0 regularizes pixel 1, pixel 1 must also regularize pixel 0). For example: B = [-1, 1] [0->1] [-1, -1] 1 now also regularizes 0 - For a instance regularization coefficient this would NOT produce a positive-definite matrix. However, for + For ``Constant`` regularization this would NOT produce a positive-definite matrix. However, for the weighted scheme, it does! The regularize weight_list change the B matrix as shown below - we simply multiply each pixel's effective diff --git a/autoarray/inversion/regularization/adaptive_brightness_split.py b/autoarray/inversion/regularization/adaptive_brightness_split.py index 06f529ab5..ef23025f8 100644 --- a/autoarray/inversion/regularization/adaptive_brightness_split.py +++ b/autoarray/inversion/regularization/adaptive_brightness_split.py @@ -31,13 +31,13 @@ def __init__( different regions of a pixelization's mesh require different levels of regularization (e.g., high smoothing where the no signal is present and less smoothing where it is, see (Nightingale, Dye and Massey 2018)). - Unlike the instance regularization_matrix scheme, neighboring pixels must now be regularized with one another + Unlike ``Constant`` regularization, neighboring pixels must now be regularized with one another in both directions (e.g. if pixel 0 regularizes pixel 1, pixel 1 must also regularize pixel 0). For example: B = [-1, 1] [0->1] [-1, -1] 1 now also regularizes 0 - For a instance regularization coefficient this would NOT produce a positive-definite matrix. However, for + For ``Constant`` regularization this would NOT produce a positive-definite matrix. However, for the weighted scheme, it does! The regularize weight_list change the B matrix as shown below - we simply multiply each pixel's effective diff --git a/autoarray/inversion/regularization/brightness_zeroth.py b/autoarray/inversion/regularization/brightness_zeroth.py new file mode 100644 index 000000000..65e22187b --- /dev/null +++ b/autoarray/inversion/regularization/brightness_zeroth.py @@ -0,0 +1,89 @@ +from __future__ import annotations +import numpy as np +from typing import TYPE_CHECKING + +if TYPE_CHECKING: + from autoarray.inversion.linear_obj.linear_obj import LinearObj + +from autoarray.inversion.regularization.abstract import AbstractRegularization + +from autoarray.inversion.regularization import regularization_util + + +class BrightnessZeroth(AbstractRegularization): + def __init__( + self, + coefficient: float = 1.0, + signal_scale: float = 1.0, + ): + """ + An adaptive regularization scheme which applies zeroth order regularization to pixels with low expected + signal values. + + For the weighted regularization scheme, each pixel is given an 'effective regularization weight', which is + controls the degree of zeroth order regularization applied to each pixel. The motivation of this is that + the exterior regions different regions of a pixelization's mesh ought to have a signal consistent with zero, + but may have a low level of non-zero signal when fitting the data. + + To implement this regularization, values on the diagonal of the regularization matrix are increased + according to the regularization weight_list of each pixel. + + Parameters + ---------- + coefficient + The regularization coefficient which controls the degree of zeroth order regularizaiton applied to + the inversion reconstruction, in regions of low signal. + signal_scale + A factor which controls how rapidly the smoothness of regularization varies from high signal regions to + low signal regions. + """ + + super().__init__() + + self.coefficient = coefficient + self.signal_scale = signal_scale + + def regularization_weights_from(self, linear_obj: LinearObj) -> np.ndarray: + """ + Returns the regularization weights of the ``BrightnessZeroth`` regularization scheme. + + The regularization weights define the level of regularization applied to each parameter in the linear object + (e.g. the ``pixels`` in a ``Mapper``). + + Parameters + ---------- + linear_obj + The linear object (e.g. a ``Mapper``) which uses these weights when performing regularization. + + Returns + ------- + The regularization weights. + """ + pixel_signals = linear_obj.pixel_signals_from(signal_scale=self.signal_scale) + + return regularization_util.adaptive_regularization_weights_from( + inner_coefficient=self.inner_coefficient, + outer_coefficient=self.outer_coefficient, + pixel_signals=pixel_signals, + ) + + def regularization_matrix_from(self, linear_obj: LinearObj) -> np.ndarray: + """ + Returns the regularization matrix of this regularization scheme. + + Parameters + ---------- + linear_obj + The linear object (e.g. a ``Mapper``) which uses this matrix to perform regularization. + + Returns + ------- + The regularization matrix. + """ + regularization_weights = self.regularization_weights_from(linear_obj=linear_obj) + + return regularization_util.weighted_regularization_matrix_from( + regularization_weights=regularization_weights, + neighbors=linear_obj.source_plane_mesh_grid.neighbors, + neighbors_sizes=linear_obj.source_plane_mesh_grid.neighbors.sizes, + ) diff --git a/autoarray/inversion/regularization/regularization_util.py b/autoarray/inversion/regularization/regularization_util.py index 797e9239c..ebf60fa21 100644 --- a/autoarray/inversion/regularization/regularization_util.py +++ b/autoarray/inversion/regularization/regularization_util.py @@ -135,13 +135,17 @@ def adaptive_regularization_weights_from( inner_coefficient: float, outer_coefficient: float, pixel_signals: np.ndarray ) -> np.ndarray: """ - Returns the regularization weight_list (the effective regularization coefficient of every pixel). They are computed - using the pixel-signal of each pixel. + Returns the regularization weights for the adaptive regularization scheme (e.g. ``AdaptiveBrightness``). + + The weights define the effective regularization coefficient of every mesh parameter (typically pixels + of a ``Mapper``). + + They are computed using an estimate of the expected signal in each pixel. Two regularization coefficients are used, corresponding to the: - 1) (pixel_signals) - pixels with a high pixel-signal (i.e. where the signal is located in the pixelization). - 2) (1.0 - pixel_signals) - pixels with a low pixel-signal (i.e. where the signal is not located in the pixelization). + 1) pixel_signals: pixels with a high pixel-signal (i.e. where the signal is located in the pixelization). + 2) 1.0 - pixel_signals: pixels with a low pixel-signal (i.e. where the signal is not located in the pixelization). Parameters ---------- @@ -158,7 +162,7 @@ def adaptive_regularization_weights_from( Returns ------- np.ndarray - The weight_list of the adaptive regularization scheme which act as the effective regularization coefficients of + The adaptive regularization weights which act as the effective regularization coefficients of every source pixel. """ return ( @@ -166,6 +170,39 @@ def adaptive_regularization_weights_from( ) ** 2.0 +def brightness_zeroth_regularization_weights_from( + coefficient: float, pixel_signals: np.ndarray +) -> np.ndarray: + """ + Returns the regularization weights for the brighness zeroth regularization scheme (e.g. ``BrightnessZeroth``). + + The weights define the level of zeroth order regularization applied to every mesh parameter (typically pixels + of a ``Mapper``). + + They are computed using an estimate of the expected signal in each pixel. + + The zeroth order regularization coefficients is applied in combination with 1.0 - pixel_signals, which are + the pixels with a low pixel-signal (i.e. where the signal is not located near the source being reconstructed in + the pixelization). + + Parameters + ---------- + coefficient + The level of zeroth order regularization applied to every mesh parameter (typically pixels of a ``Mapper``), + with the degree applied varying based on the ``pixel_signals``. + pixel_signals + The estimated signal in every pixelization pixel, used to change the regularization weighting of high signal + and low signal pixelizations. + + Returns + ------- + np.ndarray + The zeroth order regularization weights which act as the effective level of zeroth order regularization + applied to every mesh parameter. + """ + return coefficient * (1.0 - pixel_signals) + + @numba_util.jit() def weighted_regularization_matrix_from( regularization_weights: np.ndarray, diff --git a/test_autoarray/inversion/inversion/test_factory.py b/test_autoarray/inversion/inversion/test_factory.py index f554c1ece..42f9535eb 100644 --- a/test_autoarray/inversion/inversion/test_factory.py +++ b/test_autoarray/inversion/inversion/test_factory.py @@ -652,9 +652,7 @@ def test__inversion_imaging__positive_only_solver(masked_imaging_7x7_no_blur): inversion = aa.Inversion( dataset=masked_imaging_7x7_no_blur, linear_obj_list=[linear_obj], - settings=aa.SettingsInversion( - use_w_tilde=False, use_positive_only_solver=True - ), + settings=aa.SettingsInversion(use_w_tilde=False, use_positive_only_solver=True), ) assert isinstance(inversion.linear_obj_list[0], aa.m.MockLinearObjFuncList) From 7a6b461d77006c0312915129e1536b4d518f9769 Mon Sep 17 00:00:00 2001 From: James Nightingale Date: Wed, 22 Mar 2023 14:15:07 +0000 Subject: [PATCH 2/5] unit test on brightness_zeroth_regularization_weights_fro --- .../regularization/brightness_zeroth.py | 12 ++++----- .../regularization/regularization_util.py | 2 +- .../test_regularization_util.py | 27 +++++++++++++++++++ 3 files changed, 34 insertions(+), 7 deletions(-) diff --git a/autoarray/inversion/regularization/brightness_zeroth.py b/autoarray/inversion/regularization/brightness_zeroth.py index 65e22187b..f21a2eacf 100644 --- a/autoarray/inversion/regularization/brightness_zeroth.py +++ b/autoarray/inversion/regularization/brightness_zeroth.py @@ -47,8 +47,10 @@ def regularization_weights_from(self, linear_obj: LinearObj) -> np.ndarray: """ Returns the regularization weights of the ``BrightnessZeroth`` regularization scheme. - The regularization weights define the level of regularization applied to each parameter in the linear object - (e.g. the ``pixels`` in a ``Mapper``). + The weights define the level of zeroth order regularization applied to every mesh parameter (typically pixels + of a ``Mapper``). + + They are computed using an estimate of the expected signal in each pixel. Parameters ---------- @@ -61,10 +63,8 @@ def regularization_weights_from(self, linear_obj: LinearObj) -> np.ndarray: """ pixel_signals = linear_obj.pixel_signals_from(signal_scale=self.signal_scale) - return regularization_util.adaptive_regularization_weights_from( - inner_coefficient=self.inner_coefficient, - outer_coefficient=self.outer_coefficient, - pixel_signals=pixel_signals, + return regularization_util.zeroth_regularization_matrix_from( + coefficient=self.coefficient, pixel_signals=pixel_signals ) def regularization_matrix_from(self, linear_obj: LinearObj) -> np.ndarray: diff --git a/autoarray/inversion/regularization/regularization_util.py b/autoarray/inversion/regularization/regularization_util.py index ebf60fa21..0e8430995 100644 --- a/autoarray/inversion/regularization/regularization_util.py +++ b/autoarray/inversion/regularization/regularization_util.py @@ -174,7 +174,7 @@ def brightness_zeroth_regularization_weights_from( coefficient: float, pixel_signals: np.ndarray ) -> np.ndarray: """ - Returns the regularization weights for the brighness zeroth regularization scheme (e.g. ``BrightnessZeroth``). + Returns the regularization weights for the brightness zeroth regularization scheme (e.g. ``BrightnessZeroth``). The weights define the level of zeroth order regularization applied to every mesh parameter (typically pixels of a ``Mapper``). diff --git a/test_autoarray/inversion/regularizations/test_regularization_util.py b/test_autoarray/inversion/regularizations/test_regularization_util.py index f86915df0..dc3646ea0 100644 --- a/test_autoarray/inversion/regularizations/test_regularization_util.py +++ b/test_autoarray/inversion/regularizations/test_regularization_util.py @@ -201,6 +201,33 @@ def test__adaptive_regularization_weights_from(): assert (weight_list == np.array([0.75**2.0, 0.5**2.0, 0.25**2.0])).all() +def test__brightness_zeroth_regularization_weights_from(): + + pixel_signals = np.array([1.0, 1.0, 1.0]) + + weight_list = aa.util.regularization.brightness_zeroth_regularization_weights_from( + coefficient=1.0, pixel_signals=pixel_signals + ) + + assert (weight_list == np.array([0.0, 0.0, 0.0])).all() + + pixel_signals = np.array([0.25, 0.5, 0.75]) + + weight_list = aa.util.regularization.brightness_zeroth_regularization_weights_from( + coefficient=1.0, pixel_signals=pixel_signals + ) + + assert (weight_list == np.array([0.75, 0.5, 0.25])).all() + + pixel_signals = np.array([0.25, 0.5, 0.75]) + + weight_list = aa.util.regularization.brightness_zeroth_regularization_weights_from( + coefficient=2.0, pixel_signals=pixel_signals + ) + + assert (weight_list == np.array([1.5, 1.0, 0.5])).all() + + def test__weighted_regularization_matrix_from(): neighbors = np.array([[2], [3], [0], [1]]) From df9f5199ba8887ac1febcb85b809a355acfbbe94 Mon Sep 17 00:00:00 2001 From: James Nightingale Date: Wed, 22 Mar 2023 14:23:41 +0000 Subject: [PATCH 3/5] brightness_zeroth_regularization_matrix_from --- .../regularization/regularization_util.py | 46 +++++++++++++++++-- .../test_regularization_util.py | 27 +++++++++++ 2 files changed, 70 insertions(+), 3 deletions(-) diff --git a/autoarray/inversion/regularization/regularization_util.py b/autoarray/inversion/regularization/regularization_util.py index 0e8430995..82ec9530e 100644 --- a/autoarray/inversion/regularization/regularization_util.py +++ b/autoarray/inversion/regularization/regularization_util.py @@ -210,15 +210,23 @@ def weighted_regularization_matrix_from( neighbors_sizes: np.ndarray, ) -> np.ndarray: """ - From the pixel-neighbors, setup the regularization matrix using the weighted regularization scheme. + Returns the regularization matrix of the adaptive regularization scheme (e.g. ``AdaptiveBrightness``). + + This matrix is computed using the regularization weights of every mesh pixel, which are computed using the + function ``adaptive_regularization_weights_from``. These act as the effective regularization coefficients of + every mesh pixel. + + The regularization matrix is computed using the pixel-neighbors array, which is setup using the appropriate + neighbor calculation of the corresponding ``Mapper`` class. Parameters ---------- regularization_weights - The regularization_ weight of each pixel, which governs how much smoothing is applied to that individual pixel. + The regularization weight of each pixel, adaptively governing the degree of gradient regularization + applied to each inversion parameter (e.g. mesh pixels of a ``Mapper``). neighbors An array of length (total_pixels) which provides the index of all neighbors of every pixel in - the Voronoi grid (entries of -1 correspond to no neighbor). + the mesh grid (entries of -1 correspond to no neighbor). neighbors_sizes An array of length (total_pixels) which gives the number of neighbors of every pixel in the Voronoi grid. @@ -254,6 +262,38 @@ def weighted_regularization_matrix_from( return regularization_matrix +@numba_util.jit() +def brightness_zeroth_regularization_matrix_from( + regularization_weights: np.ndarray, +) -> np.ndarray: + """ + Returns the regularization matrix of the brightness zeroth regularization scheme (e.g. ``BrightnessZeroth``). + + Parameters + ---------- + regularization_weights + The regularization weight of each pixel, adaptively governing the degree of zeroth order regularization + applied to each inversion parameter (e.g. mesh pixels of a ``Mapper``). + + Returns + ------- + np.ndarray + The regularization matrix computed using an adaptive regularization scheme where the effective regularization + coefficient of every source pixel is different. + """ + + parameters = len(regularization_weights) + + regularization_matrix = np.zeros(shape=(parameters, parameters)) + + regularization_weight = regularization_weights**2.0 + + for i in range(parameters): + regularization_matrix[i, i] += regularization_weight[i] + + return regularization_matrix + + def reg_split_from( splitted_mappings: np.ndarray, splitted_sizes: np.ndarray, diff --git a/test_autoarray/inversion/regularizations/test_regularization_util.py b/test_autoarray/inversion/regularizations/test_regularization_util.py index dc3646ea0..26a6dd012 100644 --- a/test_autoarray/inversion/regularizations/test_regularization_util.py +++ b/test_autoarray/inversion/regularizations/test_regularization_util.py @@ -514,6 +514,33 @@ def test__weighted_regularization_matrix_from(): assert regularization_matrix == pytest.approx(test_regularization_matrix, 1.0e-4) +def test__brightness_zeroth_regularization_matrix_from(): + + regularization_weights = np.ones(shape=(3,)) + + regularization_matrix = ( + aa.util.regularization.brightness_zeroth_regularization_matrix_from( + regularization_weights=regularization_weights, + ) + ) + + assert regularization_matrix == pytest.approx( + np.array([[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]), 1.0e-4 + ) + + regularization_weights = np.array([1.0, 2.0, 3.0]) + + regularization_matrix = ( + aa.util.regularization.brightness_zeroth_regularization_matrix_from( + regularization_weights=regularization_weights, + ) + ) + + assert regularization_matrix == pytest.approx( + np.array([[1.0, 0.0, 0.0], [0.0, 4.0, 0.0], [0.0, 0.0, 9.0]]), 1.0e-4 + ) + + def test__constant_pixel_splitted_regularization_matrix(): splitted_mappings = np.array( [ From f6cea8ffaa739c0064153c479f9f8c850c8ca3fc Mon Sep 17 00:00:00 2001 From: James Nightingale Date: Wed, 22 Mar 2023 14:53:07 +0000 Subject: [PATCH 4/5] added AdaptiveBrightnessSplitZeroth --- autoarray/inversion/mock/mock_mapper.py | 11 +- .../pixelization/mappers/abstract.py | 20 ++- .../inversion/regularization/__init__.py | 2 + .../adaptive_brightness_split_zeroth.py | 125 ++++++++++++++++++ .../regularization/brightness_zeroth.py | 8 +- .../test_adaptive_brightness_split_zeroth.py | 41 ++++++ 6 files changed, 186 insertions(+), 21 deletions(-) create mode 100644 autoarray/inversion/regularization/adaptive_brightness_split_zeroth.py create mode 100644 test_autoarray/inversion/regularizations/test_adaptive_brightness_split_zeroth.py diff --git a/autoarray/inversion/mock/mock_mapper.py b/autoarray/inversion/mock/mock_mapper.py index 9e9c69e60..5cd60624b 100644 --- a/autoarray/inversion/mock/mock_mapper.py +++ b/autoarray/inversion/mock/mock_mapper.py @@ -14,6 +14,7 @@ def __init__( edge_pixel_list=None, regularization=None, pix_sub_weights=None, + pix_sub_weights_split_cross=None, mapping_matrix=None, pixel_signals=None, parameters=None, @@ -29,15 +30,11 @@ def __init__( super().__init__(mapper_grids=mapper_grids, regularization=regularization) self._edge_pixel_list = edge_pixel_list - self._pix_sub_weights = pix_sub_weights - + self._pix_sub_weights_split_cross = pix_sub_weights_split_cross self._mapping_matrix = mapping_matrix - self._parameters = parameters - self._pixel_signals = pixel_signals - self._interpolated_array = interpolated_array def pixel_signals_from(self, signal_scale): @@ -59,6 +56,10 @@ def edge_pixel_list(self): def pix_sub_weights(self): return self._pix_sub_weights + @property + def pix_sub_weights_split_cross(self): + return self._pix_sub_weights_split_cross + @property def mapping_matrix(self): return self._mapping_matrix diff --git a/autoarray/inversion/pixelization/mappers/abstract.py b/autoarray/inversion/pixelization/mappers/abstract.py index fa262c788..dd8d71b2e 100644 --- a/autoarray/inversion/pixelization/mappers/abstract.py +++ b/autoarray/inversion/pixelization/mappers/abstract.py @@ -401,15 +401,8 @@ def interpolated_array_from( class PixSubWeights: def __init__(self, mappings: np.ndarray, sizes: np.ndarray, weights: np.ndarray): """ - Packages the following three quantities of a mapper: - - - `mappings` (`pix_indexes_for_sub_slim_index`): the mapping of every data pixel (given its `sub_slim_index`) - to pixelization pixels (given their `pix_indexes`). - - - `sizes`(`pix_sizes_for_sub_slim_index`): the number of mappings of every data pixel to pixelization pixels. - - - `weights` (`pix_weights_for_sub_slim_index`): the interpolation weights of every data pixel's pixelization - pixel mapping + Packages the mappings, sizes and weights of every data pixel to pixelization pixels, which are computed + from associated ``Mapper`` properties.. The need to store separately the mappings and sizes is so that the `sizes` can be easy iterated over when perform calculations for efficiency. @@ -417,9 +410,14 @@ def __init__(self, mappings: np.ndarray, sizes: np.ndarray, weights: np.ndarray) Parameters ---------- mappings - The mappings of the masked data's sub-pixels to the pixelization's pixels. + The mapping of every data pixel, given its `sub_slim_index`, to its corresponding pixelization mesh + pixels, given their `pix_indexes` (corresponds to the ``Mapper`` + property ``pix_indexes_for_sub_slim_index``) sizes - The number of pixelization pixels each masked data sub-pixel maps too. + The number of mappings of every data pixel to pixelization mesh pixels (corresponds to the ``Mapper`` + property ``pix_sizes_for_sub_slim_index``). + weights + The interpolation weights of every data pixel's pixelization pixel mapping. """ self.mappings = mappings self.sizes = sizes diff --git a/autoarray/inversion/regularization/__init__.py b/autoarray/inversion/regularization/__init__.py index 852e193df..2aaa6ee26 100644 --- a/autoarray/inversion/regularization/__init__.py +++ b/autoarray/inversion/regularization/__init__.py @@ -5,3 +5,5 @@ from .constant_split import ConstantSplit from .adaptive_brightness import AdaptiveBrightness from .adaptive_brightness_split import AdaptiveBrightnessSplit +from .brightness_zeroth import BrightnessZeroth +from .adaptive_brightness_split_zeroth import AdaptiveBrightnessSplitZeroth diff --git a/autoarray/inversion/regularization/adaptive_brightness_split_zeroth.py b/autoarray/inversion/regularization/adaptive_brightness_split_zeroth.py new file mode 100644 index 000000000..00aa514e4 --- /dev/null +++ b/autoarray/inversion/regularization/adaptive_brightness_split_zeroth.py @@ -0,0 +1,125 @@ +from __future__ import annotations +import numpy as np +from typing import TYPE_CHECKING + +if TYPE_CHECKING: + from autoarray.inversion.linear_obj.linear_obj import LinearObj + +from autoarray.inversion.regularization.adaptive_brightness import AdaptiveBrightness +from autoarray.inversion.regularization.brightness_zeroth import BrightnessZeroth +from autoarray.inversion.regularization import regularization_util + + +class AdaptiveBrightnessSplitZeroth(AdaptiveBrightness): + def __init__( + self, + zeroth_coefficient: float = 1.0, + zeroth_signal_scale: float = 1.0, + inner_coefficient: float = 1.0, + outer_coefficient: float = 1.0, + signal_scale: float = 1.0, + ): + """ + An adaptive regularization scheme which splits every source pixel into a cross of four regularization points + (regularization is described in the `Regularization` class above) and interpolates to these points in order + to apply smoothing on the solution of an `Inversion`. + + The size of this cross is determined via the size of the source-pixel, for example if the source pixel is a + Voronoi pixel the area of the pixel is computed and the distance of each point of the cross is given by + the area times 0.5. + + For the weighted regularization scheme, each pixel is given an 'effective regularization weight', which is + applied when each set of pixel neighbors are regularized with one another. The motivation of this is that + different regions of a pixelization's mesh require different levels of regularization (e.g., high smoothing where the + no signal is present and less smoothing where it is, see (Nightingale, Dye and Massey 2018)). + + Unlike ``Constant`` regularization, neighboring pixels must now be regularized with one another + in both directions (e.g. if pixel 0 regularizes pixel 1, pixel 1 must also regularize pixel 0). For example: + + B = [-1, 1] [0->1] + [-1, -1] 1 now also regularizes 0 + + For ``Constant`` regularization this would NOT produce a positive-definite matrix. However, for + the weighted scheme, it does! + + The regularize weight_list change the B matrix as shown below - we simply multiply each pixel's effective + regularization weight by each row of B it has a -1 in, so: + + regularization_weights = [1, 2, 3, 4] + + B = [-1, 1, 0 ,0] # [0->1] + [0, -2, 2 ,0] # [1->2] + [0, 0, -3 ,3] # [2->3] + [4, 0, 0 ,-4] # [3->0] + + If our -1's werent down the diagonal this would look like: + + B = [4, 0, 0 ,-4] # [3->0] + [0, -2, 2 ,0] # [1->2] + [-1, 1, 0 ,0] # [0->1] + [0, 0, -3 ,3] # [2->3] This is valid! + + Parameters + ---------- + coefficients + The regularization coefficients which controls the degree of smoothing of the inversion reconstruction in + high and low signal regions of the reconstruction. + signal_scale + A factor which controls how rapidly the smoothness of regularization varies from high signal regions to + low signal regions. + """ + + self.zeroth_coefficient = zeroth_coefficient + self.zeroth_signal_scale = zeroth_signal_scale + + super().__init__( + inner_coefficient=inner_coefficient, + outer_coefficient=outer_coefficient, + signal_scale=signal_scale, + ) + + def regularization_matrix_from(self, linear_obj: LinearObj) -> np.ndarray: + """ + Returns the regularization matrix of this regularization scheme. + + Parameters + ---------- + linear_obj + The linear object (e.g. a ``Mapper``) which uses this matrix to perform regularization. + + Returns + ------- + The regularization matrix. + """ + regularization_weights = self.regularization_weights_from(linear_obj=linear_obj) + + pix_sub_weights_split_cross = linear_obj.pix_sub_weights_split_cross + + ( + splitted_mappings, + splitted_sizes, + splitted_weights, + ) = regularization_util.reg_split_from( + splitted_mappings=pix_sub_weights_split_cross.mappings, + splitted_sizes=pix_sub_weights_split_cross.sizes, + splitted_weights=pix_sub_weights_split_cross.weights, + ) + + regularization_matrix = ( + regularization_util.pixel_splitted_regularization_matrix_from( + regularization_weights=regularization_weights, + splitted_mappings=splitted_mappings, + splitted_sizes=splitted_sizes, + splitted_weights=splitted_weights, + ) + ) + + brightness_zeroth = BrightnessZeroth( + coefficient=self.zeroth_coefficient, signal_scale=self.zeroth_signal_scale + ) + + regularization_matrix_zeroth = brightness_zeroth.regularization_matrix_from( + linear_obj=linear_obj + ) + + return regularization_matrix + regularization_matrix_zeroth diff --git a/autoarray/inversion/regularization/brightness_zeroth.py b/autoarray/inversion/regularization/brightness_zeroth.py index f21a2eacf..fc24a495c 100644 --- a/autoarray/inversion/regularization/brightness_zeroth.py +++ b/autoarray/inversion/regularization/brightness_zeroth.py @@ -63,7 +63,7 @@ def regularization_weights_from(self, linear_obj: LinearObj) -> np.ndarray: """ pixel_signals = linear_obj.pixel_signals_from(signal_scale=self.signal_scale) - return regularization_util.zeroth_regularization_matrix_from( + return regularization_util.brightness_zeroth_regularization_weights_from( coefficient=self.coefficient, pixel_signals=pixel_signals ) @@ -82,8 +82,6 @@ def regularization_matrix_from(self, linear_obj: LinearObj) -> np.ndarray: """ regularization_weights = self.regularization_weights_from(linear_obj=linear_obj) - return regularization_util.weighted_regularization_matrix_from( - regularization_weights=regularization_weights, - neighbors=linear_obj.source_plane_mesh_grid.neighbors, - neighbors_sizes=linear_obj.source_plane_mesh_grid.neighbors.sizes, + return regularization_util.brightness_zeroth_regularization_matrix_from( + regularization_weights=regularization_weights ) diff --git a/test_autoarray/inversion/regularizations/test_adaptive_brightness_split_zeroth.py b/test_autoarray/inversion/regularizations/test_adaptive_brightness_split_zeroth.py new file mode 100644 index 000000000..3cef13adf --- /dev/null +++ b/test_autoarray/inversion/regularizations/test_adaptive_brightness_split_zeroth.py @@ -0,0 +1,41 @@ +import numpy as np +import pytest + +import autoarray as aa + + +def test__regularization_matrix_from(delaunay_mapper_9_3x3): + + reg = aa.reg.AdaptiveBrightnessSplit( + inner_coefficient=1.0, outer_coefficient=2.0, signal_scale=1.0 + ) + + regularization_matrix_adaptive = reg.regularization_matrix_from( + linear_obj=delaunay_mapper_9_3x3 + ) + + reg = aa.reg.BrightnessZeroth(coefficient=3.0, signal_scale=2.0) + + regularization_matrix_zeroth = reg.regularization_matrix_from( + linear_obj=delaunay_mapper_9_3x3 + ) + + regularization_matrix = ( + regularization_matrix_adaptive + regularization_matrix_zeroth + ) + + reg = aa.reg.AdaptiveBrightnessSplitZeroth( + inner_coefficient=1.0, + outer_coefficient=2.0, + signal_scale=1.0, + zeroth_coefficient=3.0, + zeroth_signal_scale=2.0, + ) + + regularization_matrix_both = reg.regularization_matrix_from( + linear_obj=delaunay_mapper_9_3x3 + ) + + print(regularization_matrix_both - regularization_matrix) + + assert (regularization_matrix_both == regularization_matrix).all() From 81138be66411dfc6c046e01e2017da22760b3cab Mon Sep 17 00:00:00 2001 From: James Nightingale Date: Wed, 22 Mar 2023 16:06:24 +0000 Subject: [PATCH 5/5] remove print statements in tests --- .../inversion/imaging/test_inversion_imaging_util.py | 3 --- .../inversion/pixelization/mappers/test_factory.py | 2 -- .../regularizations/test_adaptive_brightness_split_zeroth.py | 2 -- test_autoarray/plot/get_visuals/test_two_d.py | 2 -- test_autoarray/structures/arrays/test_uniform_2d.py | 3 --- test_autoarray/structures/mesh/test_delaunay.py | 5 ----- test_autoarray/structures/test_structure_decorators.py | 2 -- test_autoarray/structures/vectors/test_vectors_uniform.py | 1 - test_autoarray/test_preloads.py | 3 --- 9 files changed, 23 deletions(-) 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 b573caa1f..e81d91e94 100644 --- a/test_autoarray/inversion/inversion/imaging/test_inversion_imaging_util.py +++ b/test_autoarray/inversion/inversion/imaging/test_inversion_imaging_util.py @@ -247,9 +247,6 @@ def test__data_vector_via_w_tilde_data_two_methods_agree(): ) ) - # print(data_vector_via_w_tilde) - # print(data_vector) - assert data_vector_via_w_tilde == pytest.approx(data_vector, 1.0e-4) diff --git a/test_autoarray/inversion/pixelization/mappers/test_factory.py b/test_autoarray/inversion/pixelization/mappers/test_factory.py index 58b0ed929..a50db0cde 100644 --- a/test_autoarray/inversion/pixelization/mappers/test_factory.py +++ b/test_autoarray/inversion/pixelization/mappers/test_factory.py @@ -144,8 +144,6 @@ def test__delaunay_mapper(): assert (mapper.source_plane_mesh_grid == sparse_grid).all() assert mapper.source_plane_mesh_grid.origin == pytest.approx((0.0, 0.0), 1.0e-4) - print(mapper.mapping_matrix) - assert mapper.mapping_matrix == pytest.approx( np.array( [ diff --git a/test_autoarray/inversion/regularizations/test_adaptive_brightness_split_zeroth.py b/test_autoarray/inversion/regularizations/test_adaptive_brightness_split_zeroth.py index 3cef13adf..72d888d93 100644 --- a/test_autoarray/inversion/regularizations/test_adaptive_brightness_split_zeroth.py +++ b/test_autoarray/inversion/regularizations/test_adaptive_brightness_split_zeroth.py @@ -36,6 +36,4 @@ def test__regularization_matrix_from(delaunay_mapper_9_3x3): linear_obj=delaunay_mapper_9_3x3 ) - print(regularization_matrix_both - regularization_matrix) - assert (regularization_matrix_both == regularization_matrix).all() diff --git a/test_autoarray/plot/get_visuals/test_two_d.py b/test_autoarray/plot/get_visuals/test_two_d.py index 385b3edda..81a69bf9b 100644 --- a/test_autoarray/plot/get_visuals/test_two_d.py +++ b/test_autoarray/plot/get_visuals/test_two_d.py @@ -77,8 +77,6 @@ def test__via_mapper_for_data_from(voronoi_mapper_9_3x3): == voronoi_mapper_9_3x3.source_plane_data_grid.mask.derive_grid.border_sub_1.binned ).all() - print(visuals_2d.mesh_grid) - print(voronoi_mapper_9_3x3.source_plane_mesh_grid) assert ( visuals_2d_via.mesh_grid == voronoi_mapper_9_3x3.image_plane_mesh_grid ).all() diff --git a/test_autoarray/structures/arrays/test_uniform_2d.py b/test_autoarray/structures/arrays/test_uniform_2d.py index 7d44a785f..67e00b4a7 100644 --- a/test_autoarray/structures/arrays/test_uniform_2d.py +++ b/test_autoarray/structures/arrays/test_uniform_2d.py @@ -610,9 +610,6 @@ def test__extent_of_zoomed_array(): extent = arr_masked.extent_of_zoomed_array(buffer=1) - print(extent) - print(arr_masked.origin) - assert extent == pytest.approx(np.array([-4.0, 6.0, -2.0, 3.0]), 1.0e-4) diff --git a/test_autoarray/structures/mesh/test_delaunay.py b/test_autoarray/structures/mesh/test_delaunay.py index 96956e923..721104c43 100644 --- a/test_autoarray/structures/mesh/test_delaunay.py +++ b/test_autoarray/structures/mesh/test_delaunay.py @@ -61,11 +61,6 @@ def test__edge_pixel_list(): mesh = aa.Mesh2DDelaunay(values=grid) - print(mesh.delaunay.neighbors) - print(mesh.delaunay.vertex_neighbor_vertices) - print(mesh.neighbors) - print(mesh.edge_pixel_list) - assert mesh.edge_pixel_list == [0, 1, 2, 3, 5, 6, 7, 8] diff --git a/test_autoarray/structures/test_structure_decorators.py b/test_autoarray/structures/test_structure_decorators.py index 21a906bfa..aec2da7ce 100644 --- a/test_autoarray/structures/test_structure_decorators.py +++ b/test_autoarray/structures/test_structure_decorators.py @@ -723,8 +723,6 @@ def test__in_grid_2d_iterate__out_ndarray_yx_2d_list__values_use_iteration(): values_sub_3 = ndarray_2d_from(grid=grid_sub_3, profile=None) values_sub_3 = grid_sub_3.structure_2d_from(result=values_sub_3) - print(type(ndarray_yx_2d_list[0])) - assert isinstance(ndarray_yx_2d_list[0], aa.VectorYX2D) assert (ndarray_yx_2d_list[0][0] == values_sub_3.binned[0]).all() assert (ndarray_yx_2d_list[0][1] == values_sub_3.binned[1]).all() diff --git a/test_autoarray/structures/vectors/test_vectors_uniform.py b/test_autoarray/structures/vectors/test_vectors_uniform.py index c51181567..33992d8ac 100644 --- a/test_autoarray/structures/vectors/test_vectors_uniform.py +++ b/test_autoarray/structures/vectors/test_vectors_uniform.py @@ -56,7 +56,6 @@ def test__no_mask(): ).all() assert (vectors.binned.native == np.array([[[4.0, 5.0]]])).all() - print(vectors.grid) assert (vectors.binned.slim == np.array([[4.0, 5.0]])).all() assert ( vectors.grid.native diff --git a/test_autoarray/test_preloads.py b/test_autoarray/test_preloads.py index eae6315bd..1f54ef3ba 100644 --- a/test_autoarray/test_preloads.py +++ b/test_autoarray/test_preloads.py @@ -276,9 +276,6 @@ def test__set_operated_mapping_matrix_with_preloads(): assert (preloads.operated_mapping_matrix == operated_mapping_matrix_0).all() - print(preloads.curvature_matrix_preload) - print(curvature_matrix_preload) - assert ( preloads.curvature_matrix_preload == curvature_matrix_preload.astype("int") ).all()