diff --git a/autogalaxy/config/priors/light/linear/point_source.yaml b/autogalaxy/config/priors/light/linear/point_source.yaml new file mode 100644 index 00000000..f6c3cb22 --- /dev/null +++ b/autogalaxy/config/priors/light/linear/point_source.yaml @@ -0,0 +1,21 @@ +PointSource: + centre_0: + type: Gaussian + mean: 0.0 + sigma: 0.3 + width_modifier: + type: Absolute + value: 0.05 + limits: + lower: -inf + upper: inf + centre_1: + type: Gaussian + mean: 0.0 + sigma: 0.3 + width_modifier: + type: Absolute + value: 0.05 + limits: + lower: -inf + upper: inf diff --git a/autogalaxy/config/priors/light/standard/point_source.yaml b/autogalaxy/config/priors/light/standard/point_source.yaml new file mode 100644 index 00000000..40523758 --- /dev/null +++ b/autogalaxy/config/priors/light/standard/point_source.yaml @@ -0,0 +1,31 @@ +PointSource: + centre_0: + type: Gaussian + mean: 0.0 + sigma: 0.3 + width_modifier: + type: Absolute + value: 0.05 + limits: + lower: -inf + upper: inf + centre_1: + type: Gaussian + mean: 0.0 + sigma: 0.3 + width_modifier: + type: Absolute + value: 0.05 + limits: + lower: -inf + upper: inf + intensity: + type: LogUniform + lower_limit: 1.0e-06 + upper_limit: 1000000.0 + width_modifier: + type: Relative + value: 0.5 + limits: + lower: 0.0 + upper: inf diff --git a/autogalaxy/galaxy/galaxies.py b/autogalaxy/galaxy/galaxies.py index d972bd51..b8607f01 100644 --- a/autogalaxy/galaxy/galaxies.py +++ b/autogalaxy/galaxy/galaxies.py @@ -10,6 +10,7 @@ In a typical modeling workflow, a list of fitted galaxies is always wrapped in a `Galaxies` object, which is then passed to a `Fit*` class (e.g. `FitImaging`) for comparison against the observed data. """ + import numpy as np from typing import Dict, List, Optional, Tuple, Type, Union @@ -109,6 +110,26 @@ def image_2d_list_from( for galaxy in self ] + def image_2d_list_unbinned_from( + self, grid: aa.Grid2D, xp=np, operated_only: Optional[bool] = None + ) -> List[np.ndarray]: + return [ + galaxy.image_2d_unbinned_from(grid=grid, xp=xp, operated_only=operated_only) + for galaxy in self + ] + + def image_2d_unbinned_from( + self, grid: aa.Grid2D, xp=np, operated_only: Optional[bool] = None + ) -> np.ndarray: + image_2d_list = self.image_2d_list_unbinned_from( + grid=grid, xp=xp, operated_only=operated_only + ) + + if image_2d_list: + return sum(image_2d_list) + + return xp.zeros((grid.over_sampled.shape[0],)) + @aa.decorators.to_array def image_2d_from( self, grid: aa.type.Grid2DLike, xp=np, operated_only: Optional[bool] = None @@ -141,7 +162,7 @@ def galaxy_image_2d_dict_from( grid: aa.type.Grid2DLike, xp=np, operated_only: Optional[bool] = None, - ) -> {Galaxy: np.ndarray}: + ) -> Dict[Galaxy, np.ndarray]: """ Returns a dictionary associating every `Galaxy` object with its corresponding 2D image, using the instance of each galaxy as the dictionary keys. @@ -175,6 +196,17 @@ def galaxy_image_2d_dict_from( return galaxy_image_2d_dict + def galaxy_image_2d_dict_unbinned_from( + self, + grid: aa.Grid2D, + xp=np, + operated_only: Optional[bool] = None, + ) -> Dict[Galaxy, np.ndarray]: + image_2d_list = self.image_2d_list_unbinned_from( + grid=grid, xp=xp, operated_only=operated_only + ) + return {galaxy: image_2d_list[index] for index, galaxy in enumerate(self)} + @aa.decorators.to_vector_yx def deflections_yx_2d_from( self, grid: aa.type.Grid2DLike, xp=np, **kwargs diff --git a/autogalaxy/galaxy/galaxy.py b/autogalaxy/galaxy/galaxy.py index 74934688..43a92521 100644 --- a/autogalaxy/galaxy/galaxy.py +++ b/autogalaxy/galaxy/galaxy.py @@ -10,6 +10,7 @@ The `Galaxies` class (in `galaxies.py`) wraps a list of `Galaxy` objects and provides the same aggregate interface over the whole ensemble. """ + from typing import Dict, List, Optional, Type, Union import numpy as np @@ -210,6 +211,31 @@ def image_2d_list_from( ) ] + def image_2d_list_unbinned_from( + self, grid: aa.Grid2D, xp=np, operated_only: Optional[bool] = None + ) -> List[np.ndarray]: + """Return each non-linear profile image before over-sampling binning.""" + return [ + light_profile.image_2d_unbinned_from( + grid=grid, xp=xp, operated_only=operated_only + ) + for light_profile in self.cls_list_from( + cls=LightProfile, cls_filtered=LightProfileLinear + ) + ] + + def image_2d_unbinned_from( + self, grid: aa.Grid2D, xp=np, operated_only: Optional[bool] = None + ) -> np.ndarray: + image_2d_list = self.image_2d_list_unbinned_from( + grid=grid, xp=xp, operated_only=operated_only + ) + + if image_2d_list: + return sum(image_2d_list) + + return xp.zeros((grid.over_sampled.shape[0],)) + @aa.decorators.to_array def image_2d_from( self, diff --git a/autogalaxy/operate/image.py b/autogalaxy/operate/image.py index 9c8f3cda..834dccda 100644 --- a/autogalaxy/operate/image.py +++ b/autogalaxy/operate/image.py @@ -14,6 +14,7 @@ (`LightProfile`, `Galaxy`, `Galaxies`, `Tracer`) to expose a consistent API for blurring and Fourier transforming images. """ + from __future__ import annotations import numpy as np from typing import TYPE_CHECKING, Dict, List, Optional @@ -56,6 +57,18 @@ def image_2d_from( """ raise NotImplementedError + def image_2d_unbinned_from( + self, grid: aa.Grid2D, xp=np, operated_only: Optional[bool] = None + ): + """Evaluate on every over-sampled coordinate without binning. + + Subclasses that need the parent grid's mask or over-sampling metadata + override this method. The fallback preserves the historical behaviour. + """ + return self.image_2d_from( + grid=grid.over_sampled, xp=xp, operated_only=operated_only + ) + def has(self, cls) -> bool: """ Returns `True` if any attribute of this object is an instance of `cls`, else `False`. @@ -129,14 +142,12 @@ def _psf_evaluation_grids_from(grid, blurring_grid, psf): For a regular PSF the input grids are returned unchanged (evaluation is binned to image resolution and the mask travels on the arrays). For an - oversampled PSF (`convolve_over_sample_size > 1`) the over-sampled - coordinates are returned — `grid.over_sampled` is a `Grid2DIrregular` in - per-pixel sub-block order, which the `over_sample` decorator passes through - unbinned and which is the oversampled Convolver's input format — along with - the image mask, which those coordinate arrays cannot carry themselves. + oversampled PSF (`convolve_over_sample_size > 1`) the parent grids are also + returned, along with the image mask. The caller requests unbinned values + separately so discrete profiles retain access to the parent-pixel geometry. """ if psf.convolve_over_sample_size > 1: - return grid.over_sampled, blurring_grid.over_sampled, grid.mask + return grid, blurring_grid, grid.mask return grid, blurring_grid, None @@ -208,12 +219,20 @@ def blurred_image_2d_from( ) ) - image_2d_not_operated = self.image_2d_from( - grid=evaluation_grid, xp=xp, operated_only=False - ) - blurring_image_2d_not_operated = self.image_2d_from( - grid=evaluation_blurring_grid, xp=xp, operated_only=False - ) + if convolution_mask is not None: + image_2d_not_operated = self.image_2d_unbinned_from( + grid=evaluation_grid, xp=xp, operated_only=False + ) + blurring_image_2d_not_operated = self.image_2d_unbinned_from( + grid=evaluation_blurring_grid, xp=xp, operated_only=False + ) + else: + image_2d_not_operated = self.image_2d_from( + grid=evaluation_grid, xp=xp, operated_only=False + ) + blurring_image_2d_not_operated = self.image_2d_from( + grid=evaluation_blurring_grid, xp=xp, operated_only=False + ) blurred_image_2d = self._convolved_from_evaluations( image_2d=image_2d_not_operated, @@ -314,8 +333,8 @@ def convolved_padded_image_2d_from(self, grid, psf: aa.Convolver, xp=np): mask=padded_mask, over_sample_size=over_sample_size ) - image_over_sampled = self.image_2d_from( - grid=padded_grid.over_sampled, xp=xp, operated_only=False + image_over_sampled = self.image_2d_unbinned_from( + grid=padded_grid, xp=xp, operated_only=False ) convolved = psf.convolved_image_from( @@ -433,6 +452,13 @@ class OperateImageList(OperateImage): def image_2d_list_from(self, grid: aa.Grid2D, operated_only: Optional[bool] = None): raise NotImplementedError + def image_2d_list_unbinned_from( + self, grid: aa.Grid2D, operated_only: Optional[bool] = None, xp=np + ): + return self.image_2d_list_from( + grid=grid.over_sampled, operated_only=operated_only, xp=xp + ) + def blurred_image_2d_list_from( self, grid: aa.Grid2D, @@ -468,12 +494,20 @@ def blurred_image_2d_list_from( ) ) - image_2d_not_operated_list = self.image_2d_list_from( - grid=evaluation_grid, operated_only=False - ) - blurring_image_2d_not_operated_list = self.image_2d_list_from( - grid=evaluation_blurring_grid, operated_only=False - ) + if convolution_mask is not None: + image_2d_not_operated_list = self.image_2d_list_unbinned_from( + grid=evaluation_grid, operated_only=False + ) + blurring_image_2d_not_operated_list = self.image_2d_list_unbinned_from( + grid=evaluation_blurring_grid, operated_only=False + ) + else: + image_2d_not_operated_list = self.image_2d_list_from( + grid=evaluation_grid, operated_only=False + ) + blurring_image_2d_not_operated_list = self.image_2d_list_from( + grid=evaluation_blurring_grid, operated_only=False + ) blurred_image_2d_list = [] @@ -597,6 +631,13 @@ def galaxy_image_2d_dict_from( ) -> Dict[Galaxy, aa.Array2D]: raise NotImplementedError + def galaxy_image_2d_dict_unbinned_from( + self, grid: aa.Grid2D, xp=np, operated_only: Optional[bool] = None + ): + return self.galaxy_image_2d_dict_from( + grid=grid.over_sampled, xp=xp, operated_only=operated_only + ) + def galaxy_blurred_image_2d_dict_from( self, grid, psf, blurring_grid, xp=np ) -> Dict[Galaxy, aa.Array2D]: @@ -627,13 +668,22 @@ def galaxy_blurred_image_2d_dict_from( ) ) - galaxy_image_2d_not_operated_dict = self.galaxy_image_2d_dict_from( - grid=evaluation_grid, operated_only=False, xp=xp - ) - - galaxy_blurring_image_2d_not_operated_dict = self.galaxy_image_2d_dict_from( - grid=evaluation_blurring_grid, operated_only=False, xp=xp - ) + if convolution_mask is not None: + galaxy_image_2d_not_operated_dict = self.galaxy_image_2d_dict_unbinned_from( + grid=evaluation_grid, operated_only=False, xp=xp + ) + galaxy_blurring_image_2d_not_operated_dict = ( + self.galaxy_image_2d_dict_unbinned_from( + grid=evaluation_blurring_grid, operated_only=False, xp=xp + ) + ) + else: + galaxy_image_2d_not_operated_dict = self.galaxy_image_2d_dict_from( + grid=evaluation_grid, operated_only=False, xp=xp + ) + galaxy_blurring_image_2d_not_operated_dict = self.galaxy_image_2d_dict_from( + grid=evaluation_blurring_grid, operated_only=False, xp=xp + ) galaxy_image_2d_operated_dict = self.galaxy_image_2d_dict_from( grid=grid, operated_only=True, xp=xp diff --git a/autogalaxy/profiles/basis.py b/autogalaxy/profiles/basis.py index 1debe8da..f54f184e 100644 --- a/autogalaxy/profiles/basis.py +++ b/autogalaxy/profiles/basis.py @@ -9,6 +9,7 @@ linear inversion (a single matrix solve), making the inference highly efficient regardless of how many basis components are included. """ + import numpy as np from typing import Dict, List, Optional, Union @@ -114,11 +115,17 @@ def image_2d_from( The image of the light profiles in the basis summed together. """ return sum( - self.image_2d_list_from(grid=grid, xp=xp, operated_only=operated_only) + self.image_2d_list_from( + grid=grid, xp=xp, operated_only=operated_only, **kwargs + ) ) def image_2d_list_from( - self, grid: aa.type.Grid2DLike, xp=np, operated_only: Optional[bool] = None + self, + grid: aa.type.Grid2DLike, + xp=np, + operated_only: Optional[bool] = None, + **kwargs, ) -> List[aa.Array2D]: """ Returns each image of each light profiles in the basis as a list, from a 2D grid of Cartesian (y,x) coordinates. @@ -141,18 +148,23 @@ def image_2d_list_from( ------- The image of the light profiles in the basis summed together. """ - return [ - ( - light_profile.image_2d_from( - grid=grid, xp=xp, operated_only=operated_only + image_2d_list = [] + + for light_profile in self.light_profile_list: + if not isinstance(light_profile, lp_linear.LightProfileLinear): + image_2d_list.append( + light_profile.image_2d_from( + grid=grid, xp=xp, operated_only=operated_only, **kwargs + ) ) - if not isinstance(light_profile, lp_linear.LightProfileLinear) - else aa.Array2D( - values=xp.zeros((grid.shape[0],)), mask=grid.mask + elif kwargs.get("binned", True) is False: + image_2d_list.append(xp.zeros((grid.over_sampled.shape[0],))) + else: + image_2d_list.append( + aa.Array2D(values=xp.zeros((grid.shape[0],)), mask=grid.mask) ) - ) - for light_profile in self.light_profile_list - ] + + return image_2d_list def convergence_2d_from( self, grid: aa.type.Grid2DLike, xp=np, **kwargs diff --git a/autogalaxy/profiles/light/abstract.py b/autogalaxy/profiles/light/abstract.py index 480fc405..396fbe78 100644 --- a/autogalaxy/profiles/light/abstract.py +++ b/autogalaxy/profiles/light/abstract.py @@ -8,6 +8,7 @@ The `LightProfile` class is the root of the light profile hierarchy. All concrete profiles (e.g. `Sersic`, `Exponential`, `Gaussian`) inherit from it and implement `image_2d_from` and `image_2d_via_radii_from`. """ + import numpy as np from typing import Optional, Tuple @@ -84,6 +85,22 @@ def image_2d_from( """ raise NotImplementedError() + def image_2d_unbinned_from( + self, + grid: aa.Grid2D, + xp=np, + operated_only: Optional[bool] = None, + ) -> np.ndarray: + """Evaluate the profile on every over-sampled coordinate without binning. + + Keeping the parent :class:`~autoarray.Grid2D` available is important for + discrete profiles, which need its mask and per-pixel over-sampling + metadata to normalize their image correctly. + """ + return self.image_2d_from( + grid=grid, xp=xp, operated_only=operated_only, binned=False + ) + def image_2d_via_radii_from(self, grid_radii: np.ndarray, xp=np) -> np.ndarray: """ Returns the light profile's 2D image from a 1D grid of coordinates which are the radial distance of each diff --git a/autogalaxy/profiles/light/linear/__init__.py b/autogalaxy/profiles/light/linear/__init__.py index 1a6b3fdb..c07c1f08 100644 --- a/autogalaxy/profiles/light/linear/__init__.py +++ b/autogalaxy/profiles/light/linear/__init__.py @@ -1,5 +1,6 @@ from .abstract import LightProfile, LightProfileLinear, LightProfileLinearObjFuncList from .gaussian import Gaussian, GaussianSph +from .point_source import PointSource from .moffat import Moffat, MoffatSph from .sersic import Sersic, SersicSph from .exponential import Exponential, ExponentialSph diff --git a/autogalaxy/profiles/light/linear/abstract.py b/autogalaxy/profiles/light/linear/abstract.py index 0aed16f6..d287f8ce 100644 --- a/autogalaxy/profiles/light/linear/abstract.py +++ b/autogalaxy/profiles/light/linear/abstract.py @@ -12,6 +12,7 @@ The `LightProfileLinearObjFuncList` subclass additionally supports regularization, allowing the solved intensities to be penalized by a smoothness prior. """ + import inspect import itertools import numpy as np @@ -238,14 +239,12 @@ def __init__( """ for light_profile in light_profile_list: if not isinstance(light_profile, LightProfileLinear): - raise exc.ProfileException( - """ + raise exc.ProfileException(""" A light profile that is not a LightProfileLinear object has been input into the LightProfileLinearObjFuncList object. Only children of the LightProfileLinear class can be used in a linear inversion. - """ - ) + """) super().__init__( grid=grid, regularization=regularization, settings=settings, xp=xp @@ -341,12 +340,11 @@ def operated_mapping_matrix_override(self) -> Optional[np.ndarray]: return self.mapping_matrix if self.psf.convolve_over_sample_size > 1: - # Evaluate each profile on the over-sampled coordinates (per-pixel - # sub-block order, unbinned — the oversampled Convolver's input format) - # so convolution runs at the fine resolution, mirroring - # OperateImage.blurred_image_2d_from. - evaluation_grid = self.grid.over_sampled - evaluation_blurring_grid = self.blurring_grid.over_sampled + # Retain the parent Grid2D while requesting unbinned values. A + # discrete profile needs the parent-pixel geometry to place its + # flux, whereas a bare Grid2DIrregular only contains coordinates. + evaluation_grid = self.grid + evaluation_blurring_grid = self.blurring_grid convolution_mask = self.grid.mask else: evaluation_grid = self.grid @@ -358,11 +356,20 @@ def operated_mapping_matrix_override(self) -> Optional[np.ndarray]: from autogalaxy.operate.image import OperateImage for pixel, light_profile in enumerate(self.light_profile_list): - image_2d = light_profile.image_2d_from(grid=evaluation_grid, xp=self._xp) - - blurring_image_2d = light_profile.image_2d_from( - grid=evaluation_blurring_grid, xp=self._xp - ) + if convolution_mask is not None: + image_2d = light_profile.image_2d_unbinned_from( + grid=evaluation_grid, xp=self._xp + ) + blurring_image_2d = light_profile.image_2d_unbinned_from( + grid=evaluation_blurring_grid, xp=self._xp + ) + else: + image_2d = light_profile.image_2d_from( + grid=evaluation_grid, xp=self._xp + ) + blurring_image_2d = light_profile.image_2d_from( + grid=evaluation_blurring_grid, xp=self._xp + ) if convolution_mask is not None: image_2d = OperateImage._binned_for_convolution( diff --git a/autogalaxy/profiles/light/linear/point_source.py b/autogalaxy/profiles/light/linear/point_source.py new file mode 100644 index 00000000..f48c2497 --- /dev/null +++ b/autogalaxy/profiles/light/linear/point_source.py @@ -0,0 +1,11 @@ +from typing import Tuple + +from autogalaxy.profiles.light.linear.abstract import LightProfileLinear +from autogalaxy.profiles.light import standard as lp + + +class PointSource(lp.PointSource, LightProfileLinear): + """A point source whose total flux is solved for by linear inversion.""" + + def __init__(self, centre: Tuple[float, float] = (0.0, 0.0)): + super().__init__(centre=centre, intensity=1.0) diff --git a/autogalaxy/profiles/light/standard/__init__.py b/autogalaxy/profiles/light/standard/__init__.py index 2b96031e..57dba46c 100644 --- a/autogalaxy/profiles/light/standard/__init__.py +++ b/autogalaxy/profiles/light/standard/__init__.py @@ -1,4 +1,5 @@ from .gaussian import Gaussian, GaussianSph +from .point_source import PointSource from .moffat import Moffat, MoffatSph from .sersic import ( Sersic, diff --git a/autogalaxy/profiles/light/standard/point_source.py b/autogalaxy/profiles/light/standard/point_source.py new file mode 100644 index 00000000..edce91e8 --- /dev/null +++ b/autogalaxy/profiles/light/standard/point_source.py @@ -0,0 +1,108 @@ +import numpy as np +from typing import Optional, Tuple + +import autoarray as aa + +from autogalaxy.profiles.light.abstract import LightProfile + + +class PointSource(LightProfile): + """A discrete point source of light evaluated on an image-plane grid. + + ``intensity`` is the total flux assigned to the detector pixel containing + ``centre``. When the grid is over-sampled, the flux is placed at the nearest + sub-pixel and scaled before mean-binning so the total remains invariant. + + This profile describes direct image-plane emission, such as an unresolved + star or AGN. Multiple lensed point images require the lens equation and are + instead modelled with ``ag.ps.PointFlux`` and PyAutoLens's ``PointSolver``. + """ + + def __init__( + self, + centre: Tuple[float, float] = (0.0, 0.0), + intensity: float = 0.1, + ): + super().__init__(centre=centre, ell_comps=(0.0, 0.0), intensity=intensity) + + def _over_sampled_values_from(self, grid: aa.Grid2D, xp=np) -> np.ndarray: + if grid.shape[0] == 0: + return xp.zeros((0,)) + + centre = xp.asarray(self.centre) + pixel_centres = xp.asarray(grid.array) + + pixel_distances = xp.sum(xp.square(pixel_centres - centre), axis=1) + nearest_pixel_index = xp.argmin(pixel_distances) + nearest_pixel_centre = pixel_centres[nearest_pixel_index] + + pixel_scales = xp.asarray(grid.mask.pixel_scales) + inside_pixel = xp.all( + xp.abs(centre - nearest_pixel_centre) <= 0.5 * pixel_scales + ) + + over_sampled_grid = xp.asarray(grid.over_sampled.array) + sub_to_pixel = xp.asarray(grid.over_sampler.slim_for_sub_slim) + in_nearest_pixel = sub_to_pixel == nearest_pixel_index + + sub_distances = xp.sum(xp.square(over_sampled_grid - centre), axis=1) + nearest_sub_index = xp.argmin(xp.where(in_nearest_pixel, sub_distances, xp.inf)) + + local_sub_size = xp.asarray(grid.over_sample_size.array)[nearest_pixel_index] + amplitude = self._intensity * xp.square(local_sub_size) + + return xp.where( + xp.arange(over_sampled_grid.shape[0]) == nearest_sub_index, + xp.where(inside_pixel, amplitude, 0.0), + 0.0, + ) + + def image_2d_from( + self, + grid: aa.type.Grid2DLike, + xp=np, + operated_only: Optional[bool] = None, + binned: bool = True, + **kwargs, + ): + """Return the point-source image on ``grid``. + + A uniform ``Grid2D`` provides the pixel geometry required for + flux-conserving binning. For an irregular grid, the result is the + corresponding discrete sample-space delta at the closest coordinate. + """ + if isinstance(grid, aa.Grid2D): + values = self._over_sampled_values_from(grid=grid, xp=xp) + + if operated_only is True: + values = xp.zeros_like(values) + + if not binned: + return values + + return grid.over_sampler.binned_array_2d_from(array=values, xp=xp) + + values = xp.asarray(grid.array if hasattr(grid, "array") else grid) + + if values.shape[0] == 0: + result = xp.zeros((0,)) + else: + distances = xp.sum(xp.square(values - xp.asarray(self.centre)), axis=1) + nearest_index = xp.argmin(distances) + result = xp.where( + xp.arange(values.shape[0]) == nearest_index, self._intensity, 0.0 + ) + + if operated_only is True: + result = xp.zeros_like(result) + + if isinstance(grid, aa.Grid2DIrregular): + return aa.ArrayIrregular(values=result) + + return result + + def image_2d_via_radii_from(self, grid_radii: np.ndarray, xp=np) -> np.ndarray: + raise NotImplementedError( + "PointSource is a discrete image-plane profile and cannot be " + "evaluated from radial coordinates." + ) diff --git a/test_autogalaxy/profiles/light/linear/test_abstract.py b/test_autogalaxy/profiles/light/linear/test_abstract.py index 62f72920..ee6e0d29 100644 --- a/test_autogalaxy/profiles/light/linear/test_abstract.py +++ b/test_autogalaxy/profiles/light/linear/test_abstract.py @@ -1,6 +1,7 @@ import numpy as np import pytest +import autoarray as aa import autogalaxy as ag from autogalaxy.profiles.light.linear import LightProfileLinear @@ -100,6 +101,58 @@ def test__lp_instance_from__returns_instance_with_correct_intensity(): assert lp_non_linear.intensity == 3.0 +def test__point_source_lp_instance_from__returns_standard_point_source(): + lp_linear = ag.lp_linear.PointSource(centre=(1.0, 2.0)) + + lp_non_linear = lp_linear.lp_instance_from( + linear_light_profile_intensity_dict={lp_linear: 3.0} + ) + + assert type(lp_non_linear) is ag.lp.PointSource + assert lp_non_linear.centre == (1.0, 2.0) + assert lp_non_linear.intensity == 3.0 + + +def test__point_source_operated_mapping_matrix__oversampled_psf_conserves_flux(): + mask = ag.Mask2D.all_false(shape_native=(11, 11), pixel_scales=1.0) + + over_sample_size = 2 + kernel = aa.Array2D.no_mask( + values=np.array( + [ + [0.0, 0.0, 0.0], + [0.0, 1.0, 0.0], + [0.0, 0.0, 0.0], + ] + ), + pixel_scales=1.0 / over_sample_size, + ) + psf = aa.Convolver(kernel=kernel, convolve_over_sample_size=over_sample_size) + + grid = ag.Grid2D.from_mask(mask=mask, over_sample_size=over_sample_size) + blurring_mask = mask.derive_mask.blurring_from( + kernel_shape_native=psf.kernel_shape_image_resolution, + allow_padding=True, + ) + blurring_grid = ag.Grid2D.from_mask( + mask=blurring_mask, over_sample_size=over_sample_size + ) + + lp = ag.lp_linear.PointSource(centre=(0.3, 0.3)) + func_list = LightProfileLinearObjFuncList( + grid=grid, + blurring_grid=blurring_grid, + psf=psf, + light_profile_list=[lp], + regularization=None, + ) + + operated_mapping_matrix = func_list.operated_mapping_matrix_override + + assert operated_mapping_matrix.shape == (mask.pixels_in_mask, 1) + assert np.sum(operated_mapping_matrix[:, 0]) == pytest.approx(1.0) + + def test__pytree_token_is_int_and_unique(): lp_0 = ag.lp_linear.Sersic() lp_1 = ag.lp_linear.Sersic() diff --git a/test_autogalaxy/profiles/light/standard/test_point_source.py b/test_autogalaxy/profiles/light/standard/test_point_source.py new file mode 100644 index 00000000..3b2b77ab --- /dev/null +++ b/test_autogalaxy/profiles/light/standard/test_point_source.py @@ -0,0 +1,94 @@ +import numpy as np +import pytest + +import autoarray as aa +import autogalaxy as ag + + +@pytest.mark.parametrize("over_sample_size", [1, 2, 4]) +def test__image_places_total_flux_in_pixel_containing_centre(over_sample_size): + grid = aa.Grid2D.uniform( + shape_native=(3, 3), + pixel_scales=1.0, + over_sample_size=over_sample_size, + ) + + image = ag.lp.PointSource(centre=(0.2, -0.2), intensity=3.0).image_2d_from( + grid=grid + ) + + assert image.native == pytest.approx( + np.array([[0.0, 0.0, 0.0], [0.0, 3.0, 0.0], [0.0, 0.0, 0.0]]) + ) + assert np.sum(image) == pytest.approx(3.0) + + +def test__image_is_zero_when_centre_is_outside_unmasked_grid(): + grid = aa.Grid2D.uniform(shape_native=(3, 3), pixel_scales=1.0) + + image = ag.lp.PointSource(centre=(2.0, 2.0), intensity=3.0).image_2d_from(grid=grid) + + assert image.native == pytest.approx(np.zeros((3, 3))) + + +def test__irregular_grid_uses_nearest_sample_as_discrete_delta(): + grid = aa.Grid2DIrregular(values=[(1.0, 1.0), (0.1, -0.2), (-1.0, -1.0)]) + + image = ag.lp.PointSource(centre=(0.0, 0.0), intensity=3.0).image_2d_from(grid=grid) + + assert image == pytest.approx([0.0, 3.0, 0.0]) + + +def test__oversampled_psf_convolution_conserves_flux_and_resolves_sub_pixel_shift(): + mask = aa.Mask2D.all_false(shape_native=(11, 11), pixel_scales=1.0) + + over_sample_size = 2 + kernel_size = 9 + coordinates = (np.arange(kernel_size) - (kernel_size - 1) / 2.0) / over_sample_size + yy, xx = np.meshgrid(-coordinates, coordinates, indexing="ij") + kernel = np.exp(-0.5 * (yy**2 + xx**2) / 0.8**2) + kernel = aa.Array2D.no_mask( + values=kernel / kernel.sum(), pixel_scales=1.0 / over_sample_size + ) + psf = aa.Convolver(kernel=kernel, convolve_over_sample_size=over_sample_size) + + grid = aa.Grid2D.from_mask(mask=mask, over_sample_size=over_sample_size) + blurring_mask = mask.derive_mask.blurring_from( + kernel_shape_native=psf.kernel_shape_image_resolution, + allow_padding=True, + ) + blurring_grid = aa.Grid2D.from_mask( + mask=blurring_mask, over_sample_size=over_sample_size + ) + + negative = ag.lp.PointSource( + centre=(-0.3, -0.3), intensity=3.0 + ).blurred_image_2d_from(grid=grid, blurring_grid=blurring_grid, psf=psf) + positive = ag.lp.PointSource( + centre=(0.3, 0.3), intensity=3.0 + ).blurred_image_2d_from(grid=grid, blurring_grid=blurring_grid, psf=psf) + + assert np.sum(negative) == pytest.approx(3.0) + assert np.sum(positive) == pytest.approx(3.0) + + negative_centroid = np.sum(np.asarray(negative)[:, None] * grid.array, axis=0) / 3.0 + positive_centroid = np.sum(np.asarray(positive)[:, None] * grid.array, axis=0) / 3.0 + + assert negative_centroid == pytest.approx((-0.25, -0.25), abs=0.002) + assert positive_centroid == pytest.approx((0.25, 0.25), abs=0.002) + + galaxy = ag.Galaxy( + redshift=0.5, + point_source=ag.lp.PointSource(centre=(0.3, 0.3), intensity=3.0), + ) + galaxy_image = galaxy.blurred_image_2d_from( + grid=grid, blurring_grid=blurring_grid, psf=psf + ) + + galaxies = ag.Galaxies(galaxies=[galaxy]) + galaxy_image_dict = galaxies.galaxy_blurred_image_2d_dict_from( + grid=grid, blurring_grid=blurring_grid, psf=psf + ) + + assert np.asarray(galaxy_image) == pytest.approx(np.asarray(positive)) + assert np.asarray(galaxy_image_dict[galaxy]) == pytest.approx(np.asarray(positive))