From 31100e86dac1c26f0ba81c607fb40fb06deef2ef Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Thu, 16 Jul 2026 10:30:43 +0100 Subject: [PATCH] Add autogalaxy.interop.coolest: COOLEST-standard profile parameter converters Bidirectional, dependency-free (dict-in/dict-out) conversion between PyAutoGalaxy light/mass profiles and COOLEST conventions: ellipticity (e1,e2) <-> (q, phi East-of-North), (y,x) <-> (x,y) centres, intermediate-axis radii, and the power-law Einstein-radius rescaling theta_cool = sqrt(q) (2/(1+q))^(1/(gamma-1)) theta_ag. Covers Sersic(+Sph), Isothermal(+Sph)->SIE, PowerLaw(+Sph)->PEMD, NFW(+Sph), ExternalShear and MassSheet->ConvergenceSheet; unsupported profiles raise ProfileException. Exposed as lazy ag.interop. Part of PyAutoLabs/PyAutoLens#612. Co-Authored-By: Claude Fable 5 --- autogalaxy/__init__.py | 6 + autogalaxy/interop/__init__.py | 1 + autogalaxy/interop/coolest/__init__.py | 37 ++ autogalaxy/interop/coolest/conventions.py | 129 +++++++ autogalaxy/interop/coolest/light.py | 111 ++++++ autogalaxy/interop/coolest/mass.py | 355 ++++++++++++++++++ test_autogalaxy/interop/__init__.py | 0 .../interop/test_coolest_conventions.py | 48 +++ test_autogalaxy/interop/test_coolest_light.py | 134 +++++++ test_autogalaxy/interop/test_coolest_mass.py | 225 +++++++++++ 10 files changed, 1046 insertions(+) create mode 100644 autogalaxy/interop/__init__.py create mode 100644 autogalaxy/interop/coolest/__init__.py create mode 100644 autogalaxy/interop/coolest/conventions.py create mode 100644 autogalaxy/interop/coolest/light.py create mode 100644 autogalaxy/interop/coolest/mass.py create mode 100644 test_autogalaxy/interop/__init__.py create mode 100644 test_autogalaxy/interop/test_coolest_conventions.py create mode 100644 test_autogalaxy/interop/test_coolest_light.py create mode 100644 test_autogalaxy/interop/test_coolest_mass.py diff --git a/autogalaxy/__init__.py b/autogalaxy/__init__.py index 01742c891..b69de46f0 100644 --- a/autogalaxy/__init__.py +++ b/autogalaxy/__init__.py @@ -141,4 +141,10 @@ def __getattr__(name): globals()["plot"] = plot return plot + if name == "interop": + import importlib + + interop = importlib.import_module(f"{__name__}.interop") + globals()["interop"] = interop + return interop raise AttributeError(f"module {__name__!r} has no attribute {name!r}") diff --git a/autogalaxy/interop/__init__.py b/autogalaxy/interop/__init__.py new file mode 100644 index 000000000..9f5956460 --- /dev/null +++ b/autogalaxy/interop/__init__.py @@ -0,0 +1 @@ +from autogalaxy.interop import coolest diff --git a/autogalaxy/interop/coolest/__init__.py b/autogalaxy/interop/coolest/__init__.py new file mode 100644 index 000000000..e27302e1c --- /dev/null +++ b/autogalaxy/interop/coolest/__init__.py @@ -0,0 +1,37 @@ +""" +Conversion between PyAutoGalaxy profiles and the COOLEST standard +(https://github.com/aymgal/COOLEST). + +This package converts profile *parameterisations* only — every function is +dict-in / dict-out on plain floats, so ``autogalaxy`` gains no dependency on +the ``coolest`` package. Reading and writing COOLEST JSON template files is +performed one layer up, in ``autolens.interop.coolest``, which consumes the +converters defined here. + +COOLEST conventions (see https://coolest.readthedocs.io -> Conventions): + +- Cartesian coordinates (x to the right, y up), positions in arcseconds. +- Position angles ``phi`` measured counter-clockwise from the positive y axis + ("East-of-North"), in degrees, in the interval (-90, +90]. +- Ellipticity described by the axis ratio ``q = b / a`` and ``phi``. +- Characteristic radii (Einstein radius, effective radius, scale radius) are + defined along the *intermediate axis* of the ellipse, r = sqrt(a * b). + +PyAutoGalaxy conventions: + +- (y, x) coordinates in arcseconds. +- Position angles measured counter-clockwise from the positive x axis. +- Ellipticity described by the elliptical components ``ell_comps`` (e1, e2). +- Radius conventions are profile specific (e.g. the ``PowerLaw`` Einstein + radius follows an average-axis convention, the ``Sersic`` effective radius + is already an intermediate-axis radius) — see the converters in + ``mass.py`` / ``light.py`` for the exact factor each profile requires. +""" + +from autogalaxy.interop.coolest import conventions +from autogalaxy.interop.coolest import light +from autogalaxy.interop.coolest import mass +from autogalaxy.interop.coolest.light import coolest_dict_from_light +from autogalaxy.interop.coolest.light import light_profile_from +from autogalaxy.interop.coolest.mass import coolest_dict_from_mass +from autogalaxy.interop.coolest.mass import mass_profile_from diff --git a/autogalaxy/interop/coolest/conventions.py b/autogalaxy/interop/coolest/conventions.py new file mode 100644 index 000000000..d3a4b1a47 --- /dev/null +++ b/autogalaxy/interop/coolest/conventions.py @@ -0,0 +1,129 @@ +from typing import Tuple + +import numpy as np + +from autogalaxy import convert + + +def phi_coolest_from(angle: float) -> float: + """ + Convert a PyAutoGalaxy position angle to a COOLEST position angle. + + PyAutoGalaxy position angles are measured counter-clockwise from the + positive x axis; COOLEST position angles are measured counter-clockwise + from the positive y axis ("East-of-North") and lie in (-90, +90]. + + Parameters + ---------- + angle + Position angle in degrees, counter-clockwise from the positive x axis. + + Returns + ------- + Position angle in degrees, counter-clockwise from the positive y axis, + normalized to (-90, +90]. + """ + phi = float(angle) - 90.0 + while phi <= -90.0: + phi += 180.0 + while phi > 90.0: + phi -= 180.0 + return phi + + +def angle_from_phi_coolest(phi: float) -> float: + """ + Convert a COOLEST position angle (counter-clockwise from the positive y + axis) to a PyAutoGalaxy position angle (counter-clockwise from the + positive x axis). + + Parameters + ---------- + phi + Position angle in degrees, counter-clockwise from the positive y axis. + """ + return float(phi) + 90.0 + + +def q_phi_from_ell_comps(ell_comps: Tuple[float, float]) -> Tuple[float, float]: + """ + Convert PyAutoGalaxy elliptical components (e1, e2) to the COOLEST axis + ratio ``q`` and position angle ``phi`` (East-of-North, degrees). + + Parameters + ---------- + ell_comps + The elliptical components (e1, e2) of a light or mass profile. + """ + axis_ratio, angle = convert.axis_ratio_and_angle_from(ell_comps=ell_comps) + return float(axis_ratio), phi_coolest_from(angle=float(angle)) + + +def ell_comps_from_q_phi(q: float, phi: float) -> Tuple[float, float]: + """ + Convert a COOLEST axis ratio ``q`` and position angle ``phi`` + (East-of-North, degrees) to PyAutoGalaxy elliptical components (e1, e2). + + Parameters + ---------- + q + The axis ratio (b/a) of the ellipse. + phi + Position angle in degrees, counter-clockwise from the positive y axis. + """ + ell_comps = convert.ell_comps_from( + axis_ratio=q, angle=angle_from_phi_coolest(phi=phi) + ) + return float(ell_comps[0]), float(ell_comps[1]) + + +def center_from_centre(centre: Tuple[float, float]) -> Tuple[float, float]: + """ + Convert a PyAutoGalaxy (y, x) profile centre to a COOLEST + (center_x, center_y) position. Both are arcseconds with x increasing to + the right and y increasing upwards, so only the ordering changes. + + Parameters + ---------- + centre + The (y, x) arc-second coordinates of the profile centre. + """ + return float(centre[1]), float(centre[0]) + + +def centre_from_center(center_x: float, center_y: float) -> Tuple[float, float]: + """ + Convert a COOLEST (center_x, center_y) position to a PyAutoGalaxy (y, x) + profile centre. + """ + return float(center_y), float(center_x) + + +def radius_intermediate_from(radius_major: float, q: float) -> float: + """ + Convert a major-axis radius to the COOLEST intermediate-axis radius + r = sqrt(a * b) = sqrt(q) * a of the same elliptical contour. + + Parameters + ---------- + radius_major + The radius measured along the ellipse's major axis. + q + The axis ratio (b/a) of the ellipse. + """ + return float(radius_major) * np.sqrt(q) + + +def radius_major_from(radius_intermediate: float, q: float) -> float: + """ + Convert a COOLEST intermediate-axis radius r = sqrt(a * b) to the + major-axis radius of the same elliptical contour. + + Parameters + ---------- + radius_intermediate + The intermediate-axis radius sqrt(a * b) of the elliptical contour. + q + The axis ratio (b/a) of the ellipse. + """ + return float(radius_intermediate) / np.sqrt(q) diff --git a/autogalaxy/interop/coolest/light.py b/autogalaxy/interop/coolest/light.py new file mode 100644 index 000000000..464d04836 --- /dev/null +++ b/autogalaxy/interop/coolest/light.py @@ -0,0 +1,111 @@ +""" +Bidirectional parameter conversion between PyAutoGalaxy light profiles and the +COOLEST standard's light profiles. + +The PyAutoGalaxy ``Sersic`` evaluates its image on *eccentric* radii +``sqrt(q) * sqrt(x^2 + (y/q)^2)`` (``geometry_profiles.py``), which is exactly +the COOLEST intermediate-axis radius ``sqrt(q x^2 + y^2 / q)``, and its +``intensity`` is the amplitude at ``effective_radius``. The Sersic mapping is +therefore an identity on (I_eff, theta_eff, n) — only the centre ordering and +the ellipticity / position-angle conventions change. +""" + +from typing import Callable, Dict + +from autogalaxy import exc +from autogalaxy.interop.coolest import conventions +from autogalaxy.profiles.light.standard.sersic import Sersic +from autogalaxy.profiles.light.standard.sersic import SersicSph + + +def _sersic_dict_from(profile, **kwargs) -> Dict: + 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": "Sersic", + "parameters": { + "I_eff": float(profile.intensity), + "theta_eff": float(profile.effective_radius), + "n": float(profile.sersic_index), + "q": q, + "phi": phi, + "center_x": center_x, + "center_y": center_y, + }, + } + + +def _sersic_from(parameters: Dict) -> Sersic: + q = parameters["q"] + centre = conventions.centre_from_center( + center_x=parameters["center_x"], center_y=parameters["center_y"] + ) + if q == 1.0: + return SersicSph( + centre=centre, + intensity=float(parameters["I_eff"]), + effective_radius=float(parameters["theta_eff"]), + sersic_index=float(parameters["n"]), + ) + return Sersic( + centre=centre, + ell_comps=conventions.ell_comps_from_q_phi(q=q, phi=parameters["phi"]), + intensity=float(parameters["I_eff"]), + effective_radius=float(parameters["theta_eff"]), + sersic_index=float(parameters["n"]), + ) + + +_TO_COOLEST: Dict[type, Callable] = { + Sersic: _sersic_dict_from, + SersicSph: _sersic_dict_from, +} + +_FROM_COOLEST: Dict[str, Callable] = { + "Sersic": _sersic_from, +} + + +def coolest_dict_from_light(profile) -> Dict: + """ + Convert a PyAutoGalaxy light profile to its COOLEST representation, a dict + ``{"type": , "parameters": {: }}`` with + all parameters in COOLEST conventions. + + Parameters + ---------- + profile + The light profile to convert. Supported: ``Sersic`` / ``SersicSph``. + """ + try: + to_coolest = _TO_COOLEST[type(profile)] + except KeyError: + raise exc.ProfileException( + f"The light profile {type(profile).__name__} has no COOLEST " + f"converter. Supported profiles: " + f"{sorted(cls.__name__ for cls in _TO_COOLEST)}." + ) + return to_coolest(profile) + + +def light_profile_from(profile_type: str, parameters: Dict): + """ + Build a PyAutoGalaxy light profile from a COOLEST profile type name and + its parameter dict (in COOLEST conventions). + + Parameters + ---------- + profile_type + The COOLEST profile name, e.g. "Sersic". + parameters + The COOLEST parameters of the profile, e.g. point-estimate values read + from a COOLEST template file. + """ + try: + from_coolest = _FROM_COOLEST[profile_type] + except KeyError: + raise exc.ProfileException( + f"The COOLEST light profile '{profile_type}' has no PyAutoGalaxy " + f"converter. Supported: {sorted(_FROM_COOLEST)}." + ) + return from_coolest(parameters) diff --git a/autogalaxy/interop/coolest/mass.py b/autogalaxy/interop/coolest/mass.py new file mode 100644 index 000000000..1bc554d00 --- /dev/null +++ b/autogalaxy/interop/coolest/mass.py @@ -0,0 +1,355 @@ +""" +Bidirectional parameter conversion between PyAutoGalaxy mass profiles and the +COOLEST standard's mass profiles. + +Every converter is dict-in / dict-out on plain floats — no ``coolest`` +import — see ``autogalaxy.interop.coolest.__init__`` for the convention +summary. Each mapping is defined by equating the profile's *convergence* +in the two parameterisations at every point, so a profile exported to +COOLEST and re-imported is numerically identical. + +Einstein radius conventions +--------------------------- + +The PyAutoGalaxy power-law convergence is (``power_law_core.py``): + + kappa(xi) = (3 - gamma) / (1 + q) * (theta_ag / xi)^(gamma - 1) + +where ``xi = sqrt(x^2 + (y/q)^2)`` is the *major-axis* radius of the +elliptical contour through (x, y) in the rotated profile frame. + +The COOLEST PEMD convergence (following lenstronomy / Tessore & Metcalf 2015) +is: + + kappa(r) = (3 - gamma) / 2 * (theta_cool / r)^(gamma - 1) + +where ``r = sqrt(q x^2 + y^2 / q) = sqrt(q) * xi`` is the *intermediate-axis* +radius. Equating the two at every point gives: + + theta_cool = sqrt(q) * (2 / (1 + q))^(1 / (gamma - 1)) * theta_ag + +which for the isothermal case (gamma = 2) reduces to the familiar +``theta_cool = 2 sqrt(q) / (1 + q) * theta_ag``. +""" + +from typing import Callable, Dict, Optional + +import numpy as np + +from autogalaxy import convert +from autogalaxy import exc +from autogalaxy.interop.coolest import conventions +from autogalaxy.profiles.mass.dark.nfw import NFW +from autogalaxy.profiles.mass.dark.nfw import NFWSph +from autogalaxy.profiles.mass.sheets.external_shear import ExternalShear +from autogalaxy.profiles.mass.sheets.mass_sheet import MassSheet +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 PowerLawSph + + +def einstein_radius_coolest_from( + einstein_radius: float, axis_ratio: float, slope: float = 2.0 +) -> float: + """ + Convert a PyAutoGalaxy power-law Einstein radius to the COOLEST + (intermediate-axis) Einstein radius — see the module docstring for the + derivation. + + Parameters + ---------- + einstein_radius + The PyAutoGalaxy ``einstein_radius`` parameter of the profile. + axis_ratio + The axis ratio q = b/a of the profile. + slope + The logarithmic density slope gamma of the power-law. + """ + return ( + einstein_radius + * np.sqrt(axis_ratio) + * (2.0 / (1.0 + axis_ratio)) ** (1.0 / (slope - 1.0)) + ) + + +def einstein_radius_ag_from( + theta_E: float, axis_ratio: float, slope: float = 2.0 +) -> float: + """ + Convert a COOLEST (intermediate-axis) Einstein radius ``theta_E`` to the + PyAutoGalaxy power-law ``einstein_radius`` parameter — the inverse of + ``einstein_radius_coolest_from``. + """ + return ( + theta_E + / np.sqrt(axis_ratio) + * ((1.0 + axis_ratio) / 2.0) ** (1.0 / (slope - 1.0)) + ) + + +def _sie_dict_from(profile, **kwargs) -> Dict: + 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": "SIE", + "parameters": { + "theta_E": float( + einstein_radius_coolest_from( + einstein_radius=profile.einstein_radius, + axis_ratio=q, + slope=2.0, + ) + ), + "q": q, + "phi": phi, + "center_x": center_x, + "center_y": center_y, + }, + } + + +def _isothermal_from(parameters: Dict) -> Isothermal: + q = parameters["q"] + centre = conventions.centre_from_center( + center_x=parameters["center_x"], center_y=parameters["center_y"] + ) + einstein_radius = float( + einstein_radius_ag_from(theta_E=parameters["theta_E"], axis_ratio=q, slope=2.0) + ) + # An exactly-round COOLEST profile maps to the spherical class — the + # elliptical Isothermal clips its axis ratio to 0.99999 for the stability + # of its analytic deflections, so it is not numerically exact at q = 1. + if q == 1.0: + return IsothermalSph(centre=centre, einstein_radius=einstein_radius) + return Isothermal( + centre=centre, + ell_comps=conventions.ell_comps_from_q_phi(q=q, phi=parameters["phi"]), + einstein_radius=einstein_radius, + ) + + +def _pemd_dict_from(profile, **kwargs) -> Dict: + 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( + einstein_radius_coolest_from( + einstein_radius=profile.einstein_radius, + axis_ratio=q, + slope=profile.slope, + ) + ), + "q": q, + "phi": phi, + "center_x": center_x, + "center_y": center_y, + }, + } + + +def _power_law_from(parameters: Dict) -> PowerLaw: + q = parameters["q"] + slope = parameters["gamma"] + centre = conventions.centre_from_center( + center_x=parameters["center_x"], center_y=parameters["center_y"] + ) + einstein_radius = float( + einstein_radius_ag_from( + theta_E=parameters["theta_E"], axis_ratio=q, slope=slope + ) + ) + if q == 1.0: + return PowerLawSph( + centre=centre, einstein_radius=einstein_radius, slope=slope + ) + return PowerLaw( + centre=centre, + ell_comps=conventions.ell_comps_from_q_phi(q=q, phi=parameters["phi"]), + einstein_radius=einstein_radius, + slope=slope, + ) + + +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)`` + with ``xi`` the major-axis elliptical radius and ``kappa_s`` dimensionless; + COOLEST parameterises the NFW by its scale radius ``r_s`` (intermediate + axis) and characteristic density ``rho_c``, with + ``kappa_s = rho_c * r_s / Sigma_crit``. The conversion therefore requires + the critical surface mass density ``Sigma_crit`` of the lens + configuration, which depends on redshifts and cosmology and must be + supplied by the caller (PyAutoLens computes it from the lens model). + + ``rho_c`` is returned in units of ``[sigma_crit] / arcsec`` — supply + ``sigma_crit`` in mass per arcsec**2 to obtain mass per arcsec**3. + """ + if sigma_crit is None: + raise exc.ProfileException( + "Converting an NFW profile to COOLEST requires sigma_crit, the " + "critical surface mass density of the lens configuration, because " + "COOLEST parameterises the NFW by a physical density rho_c whereas " + "the PyAutoGalaxy NFW kappa_s is dimensionless. Compute sigma_crit " + "from the lens and source redshifts and a cosmology (see " + "autolens.interop.coolest, which does this automatically)." + ) + q, phi = conventions.q_phi_from_ell_comps(ell_comps=profile.ell_comps) + center_x, center_y = conventions.center_from_centre(centre=profile.centre) + r_s = conventions.radius_intermediate_from( + radius_major=profile.scale_radius, q=q + ) + return { + "type": "NFW", + "parameters": { + "r_s": r_s, + "rho_c": float(profile.kappa_s * sigma_crit / r_s), + "q": q, + "phi": phi, + "center_x": center_x, + "center_y": center_y, + }, + } + + +def _nfw_from(parameters: Dict, sigma_crit: Optional[float] = None) -> NFW: + if sigma_crit is None: + raise exc.ProfileException( + "Converting a COOLEST NFW profile to PyAutoGalaxy requires " + "sigma_crit, the critical surface mass density of the lens " + "configuration — see autogalaxy.interop.coolest.mass._nfw_dict_from." + ) + q = parameters["q"] + r_s = parameters["r_s"] + centre = conventions.centre_from_center( + center_x=parameters["center_x"], center_y=parameters["center_y"] + ) + kappa_s = float(parameters["rho_c"] * r_s / sigma_crit) + scale_radius = float(conventions.radius_major_from(radius_intermediate=r_s, q=q)) + if q == 1.0: + return NFWSph(centre=centre, kappa_s=kappa_s, scale_radius=scale_radius) + return NFW( + centre=centre, + ell_comps=conventions.ell_comps_from_q_phi(q=q, phi=parameters["phi"]), + kappa_s=kappa_s, + scale_radius=scale_radius, + ) + + +def _external_shear_dict_from(profile, **kwargs) -> Dict: + magnitude, angle = convert.shear_magnitude_and_angle_from( + gamma_1=profile.gamma_1, gamma_2=profile.gamma_2 + ) + return { + "type": "ExternalShear", + "parameters": { + "gamma_ext": float(magnitude), + "phi_ext": conventions.phi_coolest_from(angle=float(angle)), + }, + } + + +def _external_shear_from(parameters: Dict) -> ExternalShear: + gamma_1, gamma_2 = convert.shear_gamma_1_2_from( + magnitude=parameters["gamma_ext"], + angle=conventions.angle_from_phi_coolest(phi=parameters["phi_ext"]), + ) + return ExternalShear(gamma_1=float(gamma_1), gamma_2=float(gamma_2)) + + +def _convergence_sheet_dict_from(profile, **kwargs) -> Dict: + if tuple(profile.centre) != (0.0, 0.0): + raise exc.ProfileException( + "The COOLEST ConvergenceSheet is defined with its origin fixed at " + f"(0, 0), but this MassSheet has centre={profile.centre}. Only " + "sheets centred on the origin can be converted." + ) + return { + "type": "ConvergenceSheet", + "parameters": {"kappa_s": float(profile.kappa)}, + } + + +def _mass_sheet_from(parameters: Dict) -> MassSheet: + return MassSheet(centre=(0.0, 0.0), kappa=float(parameters["kappa_s"])) + + +_TO_COOLEST: Dict[type, Callable] = { + Isothermal: _sie_dict_from, + IsothermalSph: _sie_dict_from, + PowerLaw: _pemd_dict_from, + PowerLawSph: _pemd_dict_from, + NFW: _nfw_dict_from, + NFWSph: _nfw_dict_from, + ExternalShear: _external_shear_dict_from, + MassSheet: _convergence_sheet_dict_from, +} + +_FROM_COOLEST: Dict[str, Callable] = { + "SIE": _isothermal_from, + "PEMD": _power_law_from, + "ExternalShear": _external_shear_from, + "ConvergenceSheet": _mass_sheet_from, +} + + +def coolest_dict_from_mass(profile, sigma_crit: Optional[float] = None) -> Dict: + """ + Convert a PyAutoGalaxy mass profile to its COOLEST representation, a dict + ``{"type": , "parameters": {: }}`` with + all parameters in COOLEST conventions. + + Parameters + ---------- + profile + The mass profile to convert. Supported: ``Isothermal`` / + ``IsothermalSph`` (-> SIE), ``PowerLaw`` / ``PowerLawSph`` (-> PEMD), + ``NFW`` / ``NFWSph`` (-> NFW), ``ExternalShear``, ``MassSheet`` + (-> ConvergenceSheet). + sigma_crit + The critical surface mass density of the lens configuration, required + only for NFW profiles (COOLEST uses a physical density normalization). + """ + try: + to_coolest = _TO_COOLEST[type(profile)] + except KeyError: + raise exc.ProfileException( + f"The mass profile {type(profile).__name__} has no COOLEST " + f"converter. Supported profiles: " + f"{sorted(cls.__name__ for cls in _TO_COOLEST)}." + ) + return to_coolest(profile, sigma_crit=sigma_crit) + + +def mass_profile_from( + profile_type: str, parameters: Dict, sigma_crit: Optional[float] = None +): + """ + Build a PyAutoGalaxy mass profile from a COOLEST profile type name and its + parameter dict (in COOLEST conventions). + + Parameters + ---------- + profile_type + The COOLEST profile name, e.g. "SIE", "PEMD", "NFW", "ExternalShear", + "ConvergenceSheet". + parameters + The COOLEST parameters of the profile, e.g. point-estimate values read + from a COOLEST template file. + sigma_crit + The critical surface mass density of the lens configuration, required + only for NFW profiles. + """ + if profile_type == "NFW": + return _nfw_from(parameters=parameters, sigma_crit=sigma_crit) + try: + from_coolest = _FROM_COOLEST[profile_type] + except KeyError: + raise exc.ProfileException( + f"The COOLEST mass profile '{profile_type}' has no PyAutoGalaxy " + f"converter. Supported: {sorted(_FROM_COOLEST) + ['NFW']}." + ) + return from_coolest(parameters) diff --git a/test_autogalaxy/interop/__init__.py b/test_autogalaxy/interop/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/test_autogalaxy/interop/test_coolest_conventions.py b/test_autogalaxy/interop/test_coolest_conventions.py new file mode 100644 index 000000000..efd0e5660 --- /dev/null +++ b/test_autogalaxy/interop/test_coolest_conventions.py @@ -0,0 +1,48 @@ +import numpy as np +import pytest + +import autogalaxy as ag +from autogalaxy.interop.coolest import conventions + + +def test__phi_coolest_from__shifts_and_normalizes(): + assert conventions.phi_coolest_from(angle=90.0) == pytest.approx(0.0) + assert conventions.phi_coolest_from(angle=135.0) == pytest.approx(45.0) + assert conventions.phi_coolest_from(angle=45.0) == pytest.approx(-45.0) + # angle=0 gives phi=-90, which normalizes to the (-90, 90] boundary value. + assert conventions.phi_coolest_from(angle=0.0) == pytest.approx(90.0) + + +def test__phi_round_trip__same_ellipse_orientation(): + for angle in [0.0, 30.0, 90.0, 120.0]: + phi = conventions.phi_coolest_from(angle=angle) + angle_back = conventions.angle_from_phi_coolest(phi=phi) + # Orientations are degenerate under 180 deg rotations. + assert (angle_back - angle) % 180.0 == pytest.approx(0.0, abs=1e-10) + + +def test__q_phi_and_ell_comps__round_trip(): + ell_comps = ag.convert.ell_comps_from(axis_ratio=0.6, angle=40.0) + + q, phi = conventions.q_phi_from_ell_comps(ell_comps=ell_comps) + + assert q == pytest.approx(0.6) + assert phi == pytest.approx(-50.0) + + ell_comps_back = conventions.ell_comps_from_q_phi(q=q, phi=phi) + + assert ell_comps_back == pytest.approx(ell_comps) + + +def test__centre_conversion__swaps_ordering(): + assert conventions.center_from_centre(centre=(1.0, 2.0)) == (2.0, 1.0) + assert conventions.centre_from_center(center_x=2.0, center_y=1.0) == (1.0, 2.0) + + +def test__radius_conversion__intermediate_axis(): + assert conventions.radius_intermediate_from( + radius_major=2.0, q=0.25 + ) == pytest.approx(1.0) + assert conventions.radius_major_from( + radius_intermediate=1.0, q=0.25 + ) == pytest.approx(2.0) diff --git a/test_autogalaxy/interop/test_coolest_light.py b/test_autogalaxy/interop/test_coolest_light.py new file mode 100644 index 000000000..aa29122e8 --- /dev/null +++ b/test_autogalaxy/interop/test_coolest_light.py @@ -0,0 +1,134 @@ +import numpy as np +import pytest + +import autoarray as aa +import autogalaxy as ag +from autogalaxy import exc +from autogalaxy.interop import coolest + + +def grid_yx(): + return np.array( + [ + [1.0, 0.5], + [-0.3, 0.7], + [0.2, -1.1], + [0.8, 0.9], + ] + ) + + +def coolest_sersic_image(parameters, sersic_constant, grid): + """ + The COOLEST Sersic surface brightness evaluated directly from + COOLEST-convention parameters: + + I(r) = I_eff * exp(-b_n * ((r / theta_eff)^(1/n) - 1)), + r = sqrt(q x'^2 + y'^2 / q), + + with (x', y') the coordinates rotated into the profile frame whose major + axis lies at position angle phi counter-clockwise from the +y axis. + """ + x = grid[:, 1] - parameters["center_x"] + y = grid[:, 0] - parameters["center_y"] + + alpha = np.radians(parameters["phi"] + 90.0) + x_prime = np.cos(alpha) * x + np.sin(alpha) * y + y_prime = -np.sin(alpha) * x + np.cos(alpha) * y + + q = parameters["q"] + r = np.sqrt(q * x_prime**2 + y_prime**2 / q) + + return parameters["I_eff"] * np.exp( + -sersic_constant + * ((r / parameters["theta_eff"]) ** (1.0 / parameters["n"]) - 1.0) + ) + + +def test__sersic__image_matches_coolest_analytic_form(): + profile = ag.lp.Sersic( + centre=(0.1, -0.2), + ell_comps=ag.convert.ell_comps_from(axis_ratio=0.6, angle=40.0), + intensity=2.0, + effective_radius=0.8, + sersic_index=3.0, + ) + + profile_dict = coolest.coolest_dict_from_light(profile=profile) + + assert profile_dict["type"] == "Sersic" + assert profile_dict["parameters"]["I_eff"] == pytest.approx(2.0) + assert profile_dict["parameters"]["theta_eff"] == pytest.approx(0.8) + assert profile_dict["parameters"]["n"] == pytest.approx(3.0) + assert profile_dict["parameters"]["q"] == pytest.approx(0.6) + assert profile_dict["parameters"]["phi"] == pytest.approx(-50.0) + + image = profile.image_2d_from(grid=aa.Grid2DIrregular(grid_yx())) + image_coolest = coolest_sersic_image( + parameters=profile_dict["parameters"], + sersic_constant=profile.sersic_constant, + grid=grid_yx(), + ) + + assert np.asarray(image) == pytest.approx(image_coolest, rel=1e-8) + + +def test__sersic__round_trip_is_numerically_identical(): + profile = ag.lp.Sersic( + centre=(-0.05, 0.15), + ell_comps=ag.convert.ell_comps_from(axis_ratio=0.75, angle=100.0), + intensity=1.3, + effective_radius=1.1, + sersic_index=2.2, + ) + + profile_dict = coolest.coolest_dict_from_light(profile=profile) + profile_back = coolest.light_profile_from( + profile_type=profile_dict["type"], parameters=profile_dict["parameters"] + ) + + grid = aa.Grid2DIrregular(grid_yx()) + + assert np.asarray(profile_back.image_2d_from(grid=grid)) == pytest.approx( + np.asarray(profile.image_2d_from(grid=grid)), rel=1e-8 + ) + + +def test__sersic_sph__round_trip(): + profile = ag.lp.SersicSph( + centre=(0.1, 0.2), intensity=0.5, effective_radius=0.7, sersic_index=1.5 + ) + + profile_dict = coolest.coolest_dict_from_light(profile=profile) + + assert profile_dict["parameters"]["q"] == pytest.approx(1.0) + + profile_back = coolest.light_profile_from( + profile_type=profile_dict["type"], parameters=profile_dict["parameters"] + ) + + grid = aa.Grid2DIrregular(grid_yx()) + + assert np.asarray(profile_back.image_2d_from(grid=grid)) == pytest.approx( + np.asarray(profile.image_2d_from(grid=grid)), rel=1e-8 + ) + + +def test__major_axis_along_y__coolest_phi_is_zero(): + profile = ag.lp.Sersic( + ell_comps=ag.convert.ell_comps_from(axis_ratio=0.5, angle=90.0) + ) + + profile_dict = coolest.coolest_dict_from_light(profile=profile) + + assert profile_dict["parameters"]["phi"] == pytest.approx(0.0) + + +def test__unsupported_profiles__raise_with_named_error(): + # Subclasses of Sersic (e.g. Exponential) must not silently convert as if + # they were plain Sersic profiles. + with pytest.raises(exc.ProfileException): + coolest.coolest_dict_from_light(profile=ag.lp.Exponential()) + + with pytest.raises(exc.ProfileException): + coolest.light_profile_from(profile_type="Shapelets", parameters={}) diff --git a/test_autogalaxy/interop/test_coolest_mass.py b/test_autogalaxy/interop/test_coolest_mass.py new file mode 100644 index 000000000..416abb0fb --- /dev/null +++ b/test_autogalaxy/interop/test_coolest_mass.py @@ -0,0 +1,225 @@ +import numpy as np +import pytest + +import autoarray as aa +import autogalaxy as ag +from autogalaxy import exc +from autogalaxy.interop import coolest + + +def grid_yx(): + return np.array( + [ + [1.0, 0.5], + [-0.3, 0.7], + [0.2, -1.1], + [0.8, 0.9], + [-1.2, -0.4], + ] + ) + + +def coolest_power_law_convergence(parameters, grid): + """ + The COOLEST / lenstronomy PEMD convergence evaluated directly from + COOLEST-convention parameters: + + kappa(r) = (3 - gamma) / 2 * (theta_E / r)^(gamma - 1), + r = sqrt(q x'^2 + y'^2 / q), + + with (x', y') the coordinates rotated into the profile frame whose major + axis lies at position angle phi counter-clockwise from the +y axis. + """ + x = grid[:, 1] - parameters["center_x"] + y = grid[:, 0] - parameters["center_y"] + + alpha = np.radians(parameters["phi"] + 90.0) + x_prime = np.cos(alpha) * x + np.sin(alpha) * y + y_prime = -np.sin(alpha) * x + np.cos(alpha) * y + + q = parameters["q"] + gamma = parameters.get("gamma", 2.0) + r = np.sqrt(q * x_prime**2 + y_prime**2 / q) + + return (3.0 - gamma) / 2.0 * (parameters["theta_E"] / r) ** (gamma - 1.0) + + +def test__power_law__convergence_matches_coolest_analytic_form(): + profile = ag.mp.PowerLaw( + centre=(0.1, -0.2), + ell_comps=ag.convert.ell_comps_from(axis_ratio=0.6, angle=40.0), + einstein_radius=1.4, + slope=2.3, + ) + + profile_dict = coolest.coolest_dict_from_mass(profile=profile) + + assert profile_dict["type"] == "PEMD" + + 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) + + +def test__isothermal__convergence_matches_coolest_analytic_form(): + profile = ag.mp.Isothermal( + centre=(-0.05, 0.15), + ell_comps=ag.convert.ell_comps_from(axis_ratio=0.7, angle=110.0), + einstein_radius=1.1, + ) + + profile_dict = coolest.coolest_dict_from_mass(profile=profile) + + assert profile_dict["type"] == "SIE" + assert profile_dict["parameters"]["theta_E"] == pytest.approx( + 2.0 * np.sqrt(0.7) / (1.0 + 0.7) * 1.1 + ) + + 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) + + +@pytest.mark.parametrize("axis_ratio,angle,slope", [(0.6, 40.0, 2.3), (0.9, -20.0, 1.7)]) +def test__power_law__round_trip_is_numerically_identical(axis_ratio, angle, slope): + profile = ag.mp.PowerLaw( + centre=(0.1, -0.2), + ell_comps=ag.convert.ell_comps_from(axis_ratio=axis_ratio, angle=angle), + einstein_radius=1.4, + slope=slope, + ) + + profile_dict = coolest.coolest_dict_from_mass(profile=profile) + profile_back = coolest.mass_profile_from( + profile_type=profile_dict["type"], parameters=profile_dict["parameters"] + ) + + grid = aa.Grid2DIrregular(grid_yx()) + + assert np.asarray(profile_back.convergence_2d_from(grid=grid)) == pytest.approx( + np.asarray(profile.convergence_2d_from(grid=grid)), rel=1e-8 + ) + assert np.asarray( + profile_back.deflections_yx_2d_from(grid=grid) + ) == pytest.approx( + np.asarray(profile.deflections_yx_2d_from(grid=grid)), rel=1e-6 + ) + + +def test__isothermal_sph__round_trip(): + profile = ag.mp.IsothermalSph(centre=(0.1, 0.2), einstein_radius=0.8) + + profile_dict = coolest.coolest_dict_from_mass(profile=profile) + + assert profile_dict["parameters"]["q"] == pytest.approx(1.0) + assert profile_dict["parameters"]["theta_E"] == pytest.approx(0.8) + + profile_back = coolest.mass_profile_from( + profile_type=profile_dict["type"], parameters=profile_dict["parameters"] + ) + + grid = aa.Grid2DIrregular(grid_yx()) + + assert np.asarray(profile_back.convergence_2d_from(grid=grid)) == pytest.approx( + np.asarray(profile.convergence_2d_from(grid=grid)), rel=1e-8 + ) + + +def test__external_shear__coolest_angle_is_east_of_north(): + # gamma_1 > 0 only: shear aligned with the x axis (ag angle 0) = COOLEST + # phi_ext of 90 (the (-90, 90] representative of -90). + profile_dict = coolest.coolest_dict_from_mass( + profile=ag.mp.ExternalShear(gamma_1=0.05, gamma_2=0.0) + ) + assert profile_dict["type"] == "ExternalShear" + assert profile_dict["parameters"]["gamma_ext"] == pytest.approx(0.05) + assert abs(profile_dict["parameters"]["phi_ext"]) == pytest.approx(90.0) + + # gamma_2 > 0 only: ag angle 45 = COOLEST phi_ext -45. + profile_dict = coolest.coolest_dict_from_mass( + profile=ag.mp.ExternalShear(gamma_1=0.0, gamma_2=0.05) + ) + assert profile_dict["parameters"]["phi_ext"] == pytest.approx(-45.0) + + +def test__external_shear__round_trip(): + profile = ag.mp.ExternalShear(gamma_1=0.03, gamma_2=-0.04) + + profile_dict = coolest.coolest_dict_from_mass(profile=profile) + profile_back = coolest.mass_profile_from( + profile_type=profile_dict["type"], parameters=profile_dict["parameters"] + ) + + assert profile_back.gamma_1 == pytest.approx(0.03) + assert profile_back.gamma_2 == pytest.approx(-0.04) + + +def test__nfw__round_trip_requires_sigma_crit_and_is_identical(): + profile = ag.mp.NFW( + centre=(0.05, -0.1), + ell_comps=ag.convert.ell_comps_from(axis_ratio=0.8, angle=60.0), + kappa_s=0.12, + scale_radius=8.0, + ) + + with pytest.raises(exc.ProfileException): + coolest.coolest_dict_from_mass(profile=profile) + + sigma_crit = 2.5e9 + + profile_dict = coolest.coolest_dict_from_mass( + profile=profile, sigma_crit=sigma_crit + ) + + assert profile_dict["type"] == "NFW" + assert profile_dict["parameters"]["r_s"] == pytest.approx(np.sqrt(0.8) * 8.0) + assert profile_dict["parameters"]["rho_c"] == pytest.approx( + 0.12 * sigma_crit / (np.sqrt(0.8) * 8.0) + ) + + profile_back = coolest.mass_profile_from( + profile_type="NFW", + parameters=profile_dict["parameters"], + sigma_crit=sigma_crit, + ) + + assert profile_back.kappa_s == pytest.approx(0.12) + assert profile_back.scale_radius == pytest.approx(8.0) + + grid = aa.Grid2DIrregular(grid_yx()) + + assert np.asarray(profile_back.convergence_2d_from(grid=grid)) == pytest.approx( + np.asarray(profile.convergence_2d_from(grid=grid)), rel=1e-8 + ) + + +def test__mass_sheet__round_trip_and_off_centre_raises(): + profile_dict = coolest.coolest_dict_from_mass(profile=ag.mp.MassSheet(kappa=0.05)) + + assert profile_dict["type"] == "ConvergenceSheet" + assert profile_dict["parameters"]["kappa_s"] == pytest.approx(0.05) + + profile_back = coolest.mass_profile_from( + profile_type="ConvergenceSheet", parameters=profile_dict["parameters"] + ) + + assert profile_back.kappa == pytest.approx(0.05) + + with pytest.raises(exc.ProfileException): + coolest.coolest_dict_from_mass( + profile=ag.mp.MassSheet(centre=(0.1, 0.0), kappa=0.05) + ) + + +def test__unsupported_profiles__raise_with_named_error(): + with pytest.raises(exc.ProfileException): + coolest.coolest_dict_from_mass(profile=ag.mp.gNFW()) + + with pytest.raises(exc.ProfileException): + coolest.mass_profile_from(profile_type="Chameleon", parameters={})