diff --git a/autogalaxy/analysis/model_util.py b/autogalaxy/analysis/model_util.py index c61f0124..66db8131 100644 --- a/autogalaxy/analysis/model_util.py +++ b/autogalaxy/analysis/model_util.py @@ -17,6 +17,7 @@ def mge_model_from( ell_comps_uniform_width: float = 0.2, ell_comps_sigma : float = 0.3, use_spherical: bool = False, + sigma_min: float = 1e-4, ) -> af.Collection: """ Construct a Multi-Gaussian Expansion (MGE) for the lens or source galaxy light. @@ -81,6 +82,11 @@ def mge_model_from( use_spherical If True, use ``GaussianSph`` (no ell_comps). If False (default), use ``Gaussian`` with ellipticity. + sigma_min + The smallest Gaussian width (`sigma`) in arcseconds, which sets the lower end + of the log-spaced sigma values. Defaults to ``1e-4``. Increase it (e.g. to a + tenth of the pixel scale) to stop the basis wasting components on scales the + data cannot resolve. Returns ------- @@ -97,8 +103,22 @@ def mge_model_from( from autogalaxy.profiles.light.linear import Gaussian, GaussianSph from autogalaxy.profiles.basis import Basis - # The sigma values of the Gaussians will be fixed to values spanning 0.01 to the mask radius. - log10_sigma_list = np.linspace(-4, np.log10(mask_radius), total_gaussians) + if sigma_min <= 0.0: + raise ValueError( + f"mge_model_from requires sigma_min > 0.0, got {sigma_min}." + ) + + if sigma_min > mask_radius: + raise ValueError( + f"mge_model_from requires sigma_min <= mask_radius, got sigma_min=" + f"{sigma_min} and mask_radius={mask_radius}." + ) + + # The sigma values of the Gaussians are fixed to log-spaced values spanning + # `sigma_min` (default 0.0001") to the mask radius. + log10_sigma_list = np.linspace( + np.log10(sigma_min), np.log10(mask_radius), total_gaussians + ) if use_spherical: model_cls = GaussianSph @@ -179,13 +199,14 @@ def mge_point_model_from( pixel_scales: float, total_gaussians: int = 10, centre: Tuple[float, float] = (0.0, 0.0), + sigma_min: float = 0.01, ) -> 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. + are logarithmically spaced between ``sigma_min`` 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. @@ -200,6 +221,11 @@ def mge_point_model_from( centre (y, x) centre of the point source in arc-seconds. A ±0.1 arcsecond uniform prior is placed on each coordinate. + sigma_min + The smallest Gaussian width (`sigma`) in arcseconds, which sets the lower end + of the log-spaced sigma values. Defaults to ``0.01``. Increase it (e.g. to a + tenth of the pixel scale) to stop the basis wasting components on scales the + data cannot resolve. Returns ------- @@ -220,14 +246,19 @@ def mge_point_model_from( 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) + if sigma_min <= 0.0: + raise ValueError( + f"mge_point_model_from requires sigma_min > 0.0, got {sigma_min}." + ) - log10_sigma_list = np.linspace(min_log10_sigma, max_log10_sigma, total_gaussians) + # Sigma values are logarithmically spaced between `sigma_min` (default 0.01") + # and twice the pixel scale, with a floor to keep the upper end of the list + # at or above `sigma_min` when the pixel scale is very small. + max_sigma = max(2.0 * pixel_scales, sigma_min) + + log10_sigma_list = np.linspace( + np.log10(sigma_min), np.log10(max_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) diff --git a/test_autogalaxy/analysis/test_model_util.py b/test_autogalaxy/analysis/test_model_util.py index ed7e79a7..1fa11e52 100644 --- a/test_autogalaxy/analysis/test_model_util.py +++ b/test_autogalaxy/analysis/test_model_util.py @@ -105,6 +105,72 @@ def test__mge_model_from__total_gaussians_per_basis(): assert len(instance.profile_list) == 10 +def test__mge_model_from__sigma_min_default_spans_1e4_to_mask_radius(): + model = ag.model_util.mge_model_from(mask_radius=3.0, total_gaussians=5) + + instance = model.instance_from_prior_medians() + sigma_list = [profile.sigma for profile in instance.profile_list] + + assert sigma_list[0] == pytest.approx(1e-4, 1.0e-8) + assert sigma_list[-1] == pytest.approx(3.0, 1.0e-8) + + +def test__mge_model_from__sigma_min_input_sets_smallest_gaussian(): + model = ag.model_util.mge_model_from( + mask_radius=3.0, total_gaussians=5, sigma_min=0.01 + ) + + instance = model.instance_from_prior_medians() + sigma_list = [profile.sigma for profile in instance.profile_list] + + assert sigma_list[0] == pytest.approx(0.01, 1.0e-8) + assert sigma_list[-1] == pytest.approx(3.0, 1.0e-8) + assert sigma_list == pytest.approx( + list(10 ** np.linspace(np.log10(0.01), np.log10(3.0), 5)), 1.0e-8 + ) + + +def test__mge_model_from__sigma_min_invalid_raises(): + with pytest.raises(ValueError): + ag.model_util.mge_model_from( + mask_radius=3.0, total_gaussians=5, sigma_min=0.0 + ) + + with pytest.raises(ValueError): + ag.model_util.mge_model_from( + mask_radius=3.0, total_gaussians=5, sigma_min=4.0 + ) + + +def test__mge_model_from__default_sigma_list_is_bitwise_unchanged(): + """ + The default `sigma_min=1e-4` must reproduce the hardcoded `np.linspace(-4, ...)` + ladder that predates the `sigma_min` argument EXACTLY, not approximately. + + Every fixed `sigma` feeds the PyAutoFit identifier of a run, so drift here gives + existing fits a new `unique_id`, orphaning their output directories and silently + restarting them from scratch. The identifier quantizes floats at + `RESOLUTION = 1e-8` (see `autofit.mapper.identifier`), so it does not in fact move + for drift below that -- but exact equality is the stronger guarantee and costs + nothing, catching drift ~8 orders of magnitude earlier than the identifier does. + + `pytest.approx(rel=1e-8)` is deliberately NOT used: it only fails once the ladder + has moved by a relative ~1e-7, which is already past the point where the + identifier changes. + """ + for mask_radius, total_gaussians in [(3.0, 20), (3.5, 30), (7.5, 10), (1.0, 5)]: + model = ag.model_util.mge_model_from( + mask_radius=mask_radius, total_gaussians=total_gaussians + ) + + instance = model.instance_from_prior_medians() + sigma_list = [profile.sigma for profile in instance.profile_list] + + assert sigma_list == list( + 10 ** np.linspace(-4, np.log10(mask_radius), total_gaussians) + ) + + 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 @@ -136,6 +202,51 @@ def test__mge_point_model_from__sigma_values_span_correct_range(): assert gaussian_list[-1].sigma == pytest.approx(pixel_scales * 2.0, rel=1.0e-4) +def test__mge_point_model_from__default_sigma_list_is_bitwise_unchanged(): + """ + As for `mge_model_from`, the default `sigma_min=0.01` must reproduce the + hardcoded `min_log10_sigma = -2.0` ladder that predates the argument EXACTLY, + so the identifier of an existing point-source fit does not change. + """ + for pixel_scales, total_gaussians in [(0.1, 10), (0.05, 5), (0.2, 3), (0.001, 4)]: + model = ag.model_util.mge_point_model_from( + pixel_scales=pixel_scales, total_gaussians=total_gaussians + ) + + sigma_list = [gaussian.sigma for gaussian in model.profile_list] + + max_sigma = max(2.0 * pixel_scales, 10**-2.0) + + assert sigma_list == list( + 10 ** np.linspace(-2.0, np.log10(max_sigma), total_gaussians) + ) + + +def test__mge_point_model_from__sigma_min_input_sets_smallest_gaussian(): + total_gaussians = 5 + + model = ag.model_util.mge_point_model_from( + pixel_scales=0.1, total_gaussians=total_gaussians, sigma_min=0.01 / 10.0 + ) + + sigma_list = [gaussian.sigma for gaussian in model.profile_list] + + assert sigma_list[0] == pytest.approx(0.001, 1.0e-8) + assert sigma_list[-1] == pytest.approx(0.2, 1.0e-8) + assert sigma_list == pytest.approx( + list(10 ** np.linspace(np.log10(0.001), np.log10(0.2), total_gaussians)), + 1.0e-8, + ) + + +def test__mge_point_model_from__sigma_min_invalid_raises(): + with pytest.raises(ValueError): + ag.model_util.mge_point_model_from(pixel_scales=0.1, sigma_min=0.0) + + with pytest.raises(ValueError): + ag.model_util.mge_point_model_from(pixel_scales=0.1, sigma_min=-1.0) + + def test__mge_point_model_from__shared_centre_and_ell_comps(): """ All Gaussians must share exactly the same centre prior objects and ell_comps