Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
11 changes: 6 additions & 5 deletions autoarray/inversion/mock/mock_mapper.py
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand All @@ -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):
Expand All @@ -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
Expand Down
20 changes: 9 additions & 11 deletions autoarray/inversion/pixelization/mappers/abstract.py
Original file line number Diff line number Diff line change
Expand Up @@ -401,25 +401,23 @@ 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.

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
Expand Down
2 changes: 2 additions & 0 deletions autoarray/inversion/regularization/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
4 changes: 2 additions & 2 deletions autoarray/inversion/regularization/adaptive_brightness.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
125 changes: 125 additions & 0 deletions autoarray/inversion/regularization/adaptive_brightness_split_zeroth.py
Original file line number Diff line number Diff line change
@@ -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
87 changes: 87 additions & 0 deletions autoarray/inversion/regularization/brightness_zeroth.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,87 @@
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 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
----------
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.brightness_zeroth_regularization_weights_from(
coefficient=self.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.brightness_zeroth_regularization_matrix_from(
regularization_weights=regularization_weights
)
Loading