Add mge_point_model_from to model_util for compact point-source model…

autogalaxy/analysis/model_util.py

@@ -129,6 +129,81 @@ def mge_model_from(
     )
 
 
+def mge_point_model_from(
+    pixel_scales: float,
+    total_gaussians: int = 10,
+    centre: Tuple[float, float] = (0.0, 0.0),
+) -> af.Model:
+    """
+    Construct a Multi-Gaussian Expansion (MGE) model for a compact or unresolved
+    point-like component (e.g. a nuclear starburst, AGN, or unresolved bulge).
+
+    The model is composed of ``total_gaussians`` linear Gaussians whose sigma values
+    are logarithmically spaced between 0.01 arcseconds and twice the pixel scale.
+    All Gaussians share the same centre and ellipticity components, keeping the
+    parameter count low while capturing a realistic PSF-convolved point source.
+
+    Parameters
+    ----------
+    pixel_scales
+        The pixel scale of the image in arcseconds per pixel.  The maximum Gaussian
+        width is set to ``2 * pixel_scales`` so that the model is compact relative to
+        the resolution of the data.
+    total_gaussians
+        Number of Gaussian components in the basis.
+    centre
+        (y, x) centre of the point source in arc-seconds.  A ±0.1 arcsecond uniform
+        prior is placed on each coordinate.
+
+    Returns
+    -------
+    af.Model
+        An ``autofit.Model`` wrapping a ``Basis`` of linear Gaussians.
+    """
+
+    from autogalaxy.profiles.light.linear import Gaussian
+    from autogalaxy.profiles.basis import Basis
+
+    if total_gaussians < 1:
+        raise ValueError(
+            f"mge_point_model_from requires total_gaussians >= 1, got {total_gaussians}."
+        )
+
+    if pixel_scales <= 0:
+        raise ValueError(
+            f"mge_point_model_from requires pixel_scales > 0, got {pixel_scales}."
+        )
+
+    # Sigma values are logarithmically spaced between 0.01 arcsec (10**-2)
+    # and twice the pixel scale, with a floor to avoid taking log10 of
+    # very small or non-positive values.
+    min_log10_sigma = -2.0  # corresponds to 0.01 arcsec
+    max_sigma = max(2.0 * pixel_scales, 10 ** min_log10_sigma)
+    max_log10_sigma = np.log10(max_sigma)
+
+    log10_sigma_list = np.linspace(
+        min_log10_sigma, max_log10_sigma, total_gaussians
+    )
+    centre_0 = af.UniformPrior(
+        lower_limit=centre[0] - 0.1, upper_limit=centre[0] + 0.1
+    )
+    centre_1 = af.UniformPrior(
+        lower_limit=centre[1] - 0.1, upper_limit=centre[1] + 0.1
+    )
+
+    gaussian_list = af.Collection(
+        af.Model(Gaussian) for _ in range(total_gaussians)
+    )
+
+    for i, gaussian in enumerate(gaussian_list):
+        gaussian.centre.centre_0 = centre_0
+        gaussian.centre.centre_1 = centre_1
+        gaussian.ell_comps = gaussian_list[0].ell_comps
+        gaussian.sigma = 10 ** log10_sigma_list[i]
+
+    return af.Model(Basis, profile_list=gaussian_list)
+
+
 def simulator_start_here_model_from():
 
     from autogalaxy.profiles.light.snr import Sersic

test_autogalaxy/analysis/test_model_util.py

@@ -0,0 +1,81 @@
+import pytest
+import numpy as np
+
+import autofit as af
+import autogalaxy as ag
+
+
+def test__mge_point_model_from__returns_basis_model_with_correct_gaussians():
+    """
+    mge_point_model_from should return an af.Model wrapping a Basis whose
+    profile_list contains the requested number of linear Gaussian components.
+    """
+    model = ag.model_util.mge_point_model_from(pixel_scales=0.1, total_gaussians=5)
+
+    instance = model.instance_from_prior_medians()
+
+    assert isinstance(instance, ag.lp_basis.Basis)
+    assert len(instance.profile_list) == 5
+
+
+def test__mge_point_model_from__sigma_values_span_correct_range():
+    """
+    Sigma values should run from 10^-2 = 0.01 arcseconds up to 2 * pixel_scales,
+    logarithmically spaced.
+    """
+    pixel_scales = 0.1
+    total_gaussians = 10
+
+    model = ag.model_util.mge_point_model_from(
+        pixel_scales=pixel_scales, total_gaussians=total_gaussians
+    )
+
+    gaussian_list = list(model.profile_list)
+
+    assert gaussian_list[0].sigma == pytest.approx(0.01, rel=1.0e-4)
+    assert gaussian_list[-1].sigma == pytest.approx(pixel_scales * 2.0, rel=1.0e-4)
+
+
+def test__mge_point_model_from__shared_centre_and_ell_comps():
+    """
+    All Gaussians must share exactly the same centre prior objects and ell_comps
+    prior objects so the model has only 4 free parameters total.
+    """
+    model = ag.model_util.mge_point_model_from(pixel_scales=0.1, total_gaussians=5)
+
+    gaussian_list = list(model.profile_list)
+
+    # Centres are all the same prior objects
+    for gaussian in gaussian_list[1:]:
+        assert gaussian.centre.centre_0 is gaussian_list[0].centre.centre_0
+        assert gaussian.centre.centre_1 is gaussian_list[0].centre.centre_1
+
+    # Ell_comps are all the same prior objects
+    for gaussian in gaussian_list[1:]:
+        assert gaussian.ell_comps is gaussian_list[0].ell_comps
+
+    # Only 4 free parameters: centre_0, centre_1, ell_comps_0, ell_comps_1
+    assert model.prior_count == 4
+
+
+def test__mge_point_model_from__centre_prior_bounds():
+    """
+    When a custom centre is supplied the UniformPrior limits shift by ±0.1
+    arcseconds around that centre.
+    """
+    centre = (0.3, -0.2)
+    model = ag.model_util.mge_point_model_from(
+        pixel_scales=0.1, total_gaussians=3, centre=centre
+    )
+
+    gaussian_list = list(model.profile_list)
+    centre_0_prior = gaussian_list[0].centre.centre_0
+    centre_1_prior = gaussian_list[0].centre.centre_1
+
+    assert isinstance(centre_0_prior, af.UniformPrior)
+    assert centre_0_prior.lower_limit == pytest.approx(centre[0] - 0.1, rel=1.0e-6)
+    assert centre_0_prior.upper_limit == pytest.approx(centre[0] + 0.1, rel=1.0e-6)
+
+    assert isinstance(centre_1_prior, af.UniformPrior)
+    assert centre_1_prior.lower_limit == pytest.approx(centre[1] - 0.1, rel=1.0e-6)
+    assert centre_1_prior.upper_limit == pytest.approx(centre[1] + 0.1, rel=1.0e-6)