feat: PowerLawIntermediate — intermediate-axis (COOLEST) Einstein radius power-law

autogalaxy/config/priors/mass/total/power_law.yaml

@@ -104,3 +104,69 @@ PowerLawSph:
     limits:
       lower: 1.0
       upper: 3.0
+
+PowerLawIntermediate:
+  centre_0:
+    type: Gaussian
+    mean: 0.0
+    sigma: 0.1
+    width_modifier:
+      type: Absolute
+      value: 0.05
+    limits:
+      lower: -inf
+      upper: inf
+  centre_1:
+    type: Gaussian
+    mean: 0.0
+    sigma: 0.1
+    width_modifier:
+      type: Absolute
+      value: 0.05
+    limits:
+      lower: -inf
+      upper: inf
+  einstein_radius:
+    type: Uniform
+    lower_limit: 0.0
+    upper_limit: 8.0
+    width_modifier:
+      type: Relative
+      value: 0.25
+    limits:
+      lower: 0.0
+      upper: inf
+  ell_comps_0:
+    type: TruncatedGaussian
+    mean: 0.0
+    sigma: 0.3
+    lower_limit: -1.0
+    upper_limit: 1.0
+    width_modifier:
+      type: Absolute
+      value: 0.2
+    limits:
+      lower: -1.0
+      upper: 1.0
+  ell_comps_1:
+    type: TruncatedGaussian
+    mean: 0.0
+    sigma: 0.3
+    lower_limit: -1.0
+    upper_limit: 1.0
+    width_modifier:
+      type: Absolute
+      value: 0.2
+    limits:
+      lower: -1.0
+      upper: 1.0
+  slope:
+    type: Uniform
+    lower_limit: 1.5
+    upper_limit: 3.0
+    width_modifier:
+      type: Absolute
+      value: 0.2
+    limits:
+      lower: 1.0
+      upper: 3.0

autogalaxy/interop/coolest/mass.py

@@ -46,6 +46,7 @@
 from autogalaxy.profiles.mass.total.isothermal import Isothermal
 from autogalaxy.profiles.mass.total.isothermal import IsothermalSph
 from autogalaxy.profiles.mass.total.power_law import PowerLaw
+from autogalaxy.profiles.mass.total.power_law import PowerLawIntermediate
 from autogalaxy.profiles.mass.total.power_law import PowerLawSph
 
 
@@ -174,6 +175,40 @@ def _power_law_from(parameters: Dict) -> PowerLaw:
     )
 
 
+def _pemd_dict_from_intermediate(profile, **kwargs) -> Dict:
+    """
+    The ``PowerLawIntermediate`` Einstein radius is already the COOLEST
+    intermediate-axis theta_E, so the mapping is an identity on all parameters
+    (only the ellipticity / angle / centre conventions change).
+    """
+    q, phi = conventions.q_phi_from_ell_comps(ell_comps=profile.ell_comps)
+    center_x, center_y = conventions.center_from_centre(centre=profile.centre)
+    return {
+        "type": "PEMD",
+        "parameters": {
+            "gamma": float(profile.slope),
+            "theta_E": float(profile.einstein_radius),
+            "q": q,
+            "phi": phi,
+            "center_x": center_x,
+            "center_y": center_y,
+        },
+    }
+
+
+def _power_law_intermediate_from(parameters: Dict) -> PowerLawIntermediate:
+    return PowerLawIntermediate(
+        centre=conventions.centre_from_center(
+            center_x=parameters["center_x"], center_y=parameters["center_y"]
+        ),
+        ell_comps=conventions.ell_comps_from_q_phi(
+            q=parameters["q"], phi=parameters["phi"]
+        ),
+        einstein_radius=float(parameters["theta_E"]),
+        slope=float(parameters["gamma"]),
+    )
+
+
 def _nfw_dict_from(profile, sigma_crit: Optional[float] = None, **kwargs) -> Dict:
     """
     The PyAutoGalaxy NFW convergence is ``kappa(xi) = 2 kappa_s g(xi / r_s)``
@@ -281,6 +316,7 @@ def _mass_sheet_from(parameters: Dict) -> MassSheet:
     Isothermal: _sie_dict_from,
     IsothermalSph: _sie_dict_from,
     PowerLaw: _pemd_dict_from,
+    PowerLawIntermediate: _pemd_dict_from_intermediate,
     PowerLawSph: _pemd_dict_from,
     NFW: _nfw_dict_from,
     NFWSph: _nfw_dict_from,
@@ -325,7 +361,10 @@ def coolest_dict_from_mass(profile, sigma_crit: Optional[float] = None) -> Dict:
 
 
 def mass_profile_from(
-    profile_type: str, parameters: Dict, sigma_crit: Optional[float] = None
+    profile_type: str,
+    parameters: Dict,
+    sigma_crit: Optional[float] = None,
+    intermediate: bool = False,
 ):
     """
     Build a PyAutoGalaxy mass profile from a COOLEST profile type name and its
@@ -342,9 +381,16 @@ def mass_profile_from(
     sigma_crit
         The critical surface mass density of the lens configuration, required
         only for NFW profiles.
+    intermediate
+        If True, a PEMD profile is built as a ``PowerLawIntermediate`` (its
+        ``einstein_radius`` is the template's ``theta_E`` unchanged) instead of
+        the default ``PowerLaw`` (which converts theta_E to the PyAutoGalaxy
+        major-axis convention). Both give the identical mass distribution.
     """
     if profile_type == "NFW":
         return _nfw_from(parameters=parameters, sigma_crit=sigma_crit)
+    if profile_type == "PEMD" and intermediate:
+        return _power_law_intermediate_from(parameters=parameters)
     try:
         from_coolest = _FROM_COOLEST[profile_type]
     except KeyError:

autogalaxy/profiles/mass/__init__.py

@@ -17,6 +17,7 @@
     IsothermalCore,
     IsothermalCoreSph,
     PowerLaw,
+    PowerLawIntermediate,
     PowerLawSph,
     Isothermal,
     IsothermalSph,

autogalaxy/profiles/mass/total/__init__.py

@@ -8,7 +8,7 @@
 )
 from .isothermal import Isothermal, IsothermalSph
 from .isothermal_core import IsothermalCore, IsothermalCoreSph
-from .power_law import PowerLaw, PowerLawSph
+from .power_law import PowerLaw, PowerLawIntermediate, PowerLawSph
 from .power_law_broken import PowerLawBroken, PowerLawBrokenSph
 from .power_law_core import PowerLawCore, PowerLawCoreSph
 from .power_law_multipole import PowerLawMultipole

autogalaxy/profiles/mass/total/power_law.py

@@ -103,7 +103,7 @@ def deflections_yx_2d_from(self, grid: aa.type.Grid2DLike, xp=np, **kwargs):
         slope = self.slope - 1.0
         einstein_radius = (
             2.0 / (self.axis_ratio(xp) ** -0.5 + self.axis_ratio(xp) ** 0.5)
-        ) * self.einstein_radius
+        ) * self.einstein_radius_major_from(xp)
 
         factor = xp.divide(1.0 - self.axis_ratio(xp), 1.0 + self.axis_ratio(xp))
         b = xp.multiply(einstein_radius, xp.sqrt(self.axis_ratio(xp)))
@@ -168,6 +168,61 @@ def potential_func(u, y, x, axis_ratio, slope, core_radius, xp=np):
         )
 
 
+class PowerLawIntermediate(PowerLaw):
+    r"""Elliptical power-law mass profile with an intermediate-axis Einstein radius.
+
+    Identical mass distribution to :class:`PowerLaw`, but the ``einstein_radius``
+    parameter follows the *intermediate-axis* convention used by the COOLEST
+    standard, lenstronomy and herculens: radii are measured as
+    :math:`r = \sqrt{a b} = \sqrt{q} \, a` on the elliptical isodensity
+    contours, and the convergence is
+
+    .. math::
+
+        \kappa(r) = \frac{3 - \gamma}{2}
+                    \left(\frac{\theta_{\rm E}}{r}\right)^{\gamma - 1}.
+
+    The relation to the :class:`PowerLaw` ``einstein_radius`` (:math:`\theta_{\rm PL}`) is:
+
+    .. math::
+
+        \theta_{\rm E} = \sqrt{q} \left(\frac{2}{1 + q}\right)^{1 / (\gamma - 1)} \theta_{\rm PL},
+
+    which for the isothermal case (:math:`\gamma = 2`) reduces to
+    :math:`2 \sqrt{q} / (1 + q) \, \theta_{\rm PL}`.
+
+    Use this profile when parameter values must line up directly with other lens
+    modeling codes — its COOLEST ``PEMD`` mapping is an identity
+    (``autogalaxy.interop.coolest``) — at the cost of departing from the
+    Einstein-radius convention of the rest of **PyAutoGalaxy**.
+
+    Parameters
+    ----------
+    centre : (float, float)
+        (y, x) arc-second coordinates of the profile centre.
+    ell_comps : (float, float)
+        Ellipticity components (e1, e2) of the elliptical coordinate system.
+    einstein_radius : float
+        Einstein radius in arcseconds, in the intermediate-axis (COOLEST) convention.
+    slope : float
+        Logarithmic density slope :math:`\gamma`.
+    """
+
+    def einstein_radius_major_from(self, xp=np):
+        """
+        Convert the stored intermediate-axis ``einstein_radius`` to the major-axis
+        convention the power-law formulae are written in (the inverse of the
+        relation in the class docstring). Threaded through ``xp`` so the profile
+        is JAX-traceable like its parent.
+        """
+        axis_ratio = self.axis_ratio(xp)
+        return (
+            self.einstein_radius
+            / xp.sqrt(axis_ratio)
+            * ((1.0 + axis_ratio) / 2.0) ** (1.0 / (self.slope - 1.0))
+        )
+
+
 class PowerLawSph(PowerLaw):
     r"""Spherical power-law mass profile.
 

autogalaxy/profiles/mass/total/power_law_core.py

@@ -71,14 +71,25 @@ def __init__(
         self.slope = slope
         self.core_radius = core_radius
 
+    def einstein_radius_major_from(self, xp=np):
+        """
+        The Einstein radius in the convention the profile's convergence and deflection formulae are written in,
+        where the elliptical radius is the major-axis radius of the isodensity contour.
+
+        For this class and its standard subclasses this is the ``einstein_radius`` parameter itself; subclasses
+        which parameterise the Einstein radius in a different convention (e.g. ``PowerLawIntermediate``) override
+        this with the conversion, which every calculation then inherits.
+        """
+        return self.einstein_radius
+
     def einstein_radius_rescaled(self, xp=np):
         """
         Rescale the einstein radius by slope and axis_ratio, to reduce its degeneracy with other mass-profiles
         parameters.
         """
         return (
             (3 - self.slope) / (1 + self.axis_ratio(xp))
-        ) * self.einstein_radius ** (self.slope - 1)
+        ) * self.einstein_radius_major_from(xp) ** (self.slope - 1)
 
     @aa.over_sample
     @aa.decorators.to_array

test_autogalaxy/interop/test_coolest_mass.py

@@ -112,6 +112,51 @@ def test__power_law__round_trip_is_numerically_identical(axis_ratio, angle, slop
     )
 
 
+def test__power_law_intermediate__coolest_mapping_is_identity():
+    profile = ag.mp.PowerLawIntermediate(
+        centre=(0.1, -0.2),
+        ell_comps=ag.convert.ell_comps_from(axis_ratio=0.7, angle=45.0),
+        einstein_radius=1.2,
+        slope=2.1,
+    )
+
+    profile_dict = coolest.coolest_dict_from_mass(profile=profile)
+
+    assert profile_dict["type"] == "PEMD"
+    # No sqrt(q) factor — theta_E is the parameter itself.
+    assert profile_dict["parameters"]["theta_E"] == 1.2
+    assert profile_dict["parameters"]["gamma"] == pytest.approx(2.1)
+    assert profile_dict["parameters"]["q"] == pytest.approx(0.7)
+
+    # And the analytic COOLEST convergence with these parameters matches.
+    convergence = profile.convergence_2d_from(grid=aa.Grid2DIrregular(grid_yx()))
+    convergence_coolest = coolest_power_law_convergence(
+        parameters=profile_dict["parameters"], grid=grid_yx()
+    )
+    assert np.asarray(convergence) == pytest.approx(convergence_coolest, rel=1e-8)
+
+    profile_back = coolest.mass_profile_from(
+        profile_type="PEMD",
+        parameters=profile_dict["parameters"],
+        intermediate=True,
+    )
+
+    assert isinstance(profile_back, ag.mp.PowerLawIntermediate)
+    assert profile_back.einstein_radius == pytest.approx(1.2, rel=1e-12)
+
+    # Default import still builds the PyAutoGalaxy-convention PowerLaw, with
+    # the identical mass distribution.
+    profile_back_default = coolest.mass_profile_from(
+        profile_type="PEMD", parameters=profile_dict["parameters"]
+    )
+    assert type(profile_back_default) is ag.mp.PowerLaw
+
+    grid = aa.Grid2DIrregular(grid_yx())
+    assert np.asarray(
+        profile_back_default.convergence_2d_from(grid=grid)
+    ) == pytest.approx(np.asarray(profile.convergence_2d_from(grid=grid)), rel=1e-8)
+
+
 def test__isothermal_sph__round_trip():
     profile = ag.mp.IsothermalSph(centre=(0.1, 0.2), einstein_radius=0.8)