diff --git a/slsim/Deflectors/DeflectorPopulation/galaxy_deflectors.py b/slsim/Deflectors/DeflectorPopulation/galaxy_deflectors.py index a29bc691f..e50ad21af 100644 --- a/slsim/Deflectors/DeflectorPopulation/galaxy_deflectors.py +++ b/slsim/Deflectors/DeflectorPopulation/galaxy_deflectors.py @@ -3,6 +3,7 @@ vel_disp_abundance_matching, ) from slsim.Deflectors.DeflectorPopulation.deflectors_base import DeflectorsBase +from slsim.Util.color_gradient import attach_foreground_deflector_color_gradient from astropy.table import vstack @@ -26,6 +27,8 @@ def __init__( catalog_type="skypy", mass_type="EPL", light_type="single_sersic", + foreground_color_gradient=None, + foreground_component_weights=(0.4, 0.6), ): """ :param red_galaxy_list: list of dictionary with elliptical galaxy @@ -50,6 +53,13 @@ def __init__( :type mass_type: string :param light_type: type of Source() model class for the light distribution :type light_type: string + :param foreground_color_gradient: Optional configuration for a + band-dependent two-component deflector light profile. When provided, + the light is handled by the standard ``Source``/``DoubleSersic`` path. + :type foreground_color_gradient: dict or None + :param foreground_component_weights: Reference-band flux weights of the + two Sersic components. + :type foreground_component_weights: tuple or list :param catalog_type: type of the catalog. If user is using deflector catalog other than generated from skypy pipeline, we require them to provide angular size of the galaxy in arcsec and specify catalog_type as None. Otherwise, by @@ -57,6 +67,27 @@ def __init__( pipeline. :type catalog_type: str. "skypy" or None. """ + if foreground_color_gradient is not None: + if light_type not in ("single_sersic", "double_sersic"): + raise ValueError( + "foreground_color_gradient requires a single_sersic or " + "double_sersic light_type." + ) + light_type = "double_sersic" + red_galaxy_list = red_galaxy_list.copy() + attach_foreground_deflector_color_gradient( + red_galaxy_list, + foreground_color_gradient, + component_weights=foreground_component_weights, + ) + if blue_galaxy_list is not None: + blue_galaxy_list = blue_galaxy_list.copy() + attach_foreground_deflector_color_gradient( + blue_galaxy_list, + foreground_color_gradient, + component_weights=foreground_component_weights, + ) + red_column_names = red_galaxy_list.colnames if "galaxy_type" not in red_column_names: red_galaxy_list["galaxy_type"] = "red" diff --git a/slsim/Deflectors/deflector_util.py b/slsim/Deflectors/deflector_util.py index e37d6b67b..3c4681138 100644 --- a/slsim/Deflectors/deflector_util.py +++ b/slsim/Deflectors/deflector_util.py @@ -121,7 +121,35 @@ def light2mass( ) else: # scale light to mass ellipticity - e1_light, e2_light = kwargs_source["e1"], kwargs_source["e2"] + if "e1" in kwargs_source and "e2" in kwargs_source: + e1_light, e2_light = kwargs_source["e1"], kwargs_source["e2"] + elif all( + key in kwargs_source + for key in ("e1_0", "e2_0", "e1_1", "e2_1", "w0", "w1") + ): + # Use the total reference-band light shape rather than an + # arbitrarily selected DoubleSersic component. + weights = np.asarray( + [kwargs_source["w0"], kwargs_source["w1"]], dtype=float + ) + if np.any(weights < 0) or np.sum(weights) <= 0: + raise ValueError( + "DoubleSersic weights 'w0' and 'w1' must be non-negative " + "and have a positive sum when deriving mass ellipticity." + ) + weights /= np.sum(weights) + e1_light = np.dot( + weights, [kwargs_source["e1_0"], kwargs_source["e1_1"]] + ) + e2_light = np.dot( + weights, [kwargs_source["e2_0"], kwargs_source["e2_1"]] + ) + else: + raise ValueError( + "Cannot derive mass ellipticity from the light model. Provide " + "either ('e1', 'e2') or all DoubleSersic fields " + "('e1_0', 'e2_0', 'e1_1', 'e2_1', 'w0', 'w1')." + ) e1_mass, e2_mass = ( light2mass_e_scaling * e1_light, light2mass_e_scaling * e2_light, diff --git a/slsim/ImageSimulation/image_quality_lenstronomy.py b/slsim/ImageSimulation/image_quality_lenstronomy.py index 35bd8546c..11bdeb134 100644 --- a/slsim/ImageSimulation/image_quality_lenstronomy.py +++ b/slsim/ImageSimulation/image_quality_lenstronomy.py @@ -1,7 +1,9 @@ +import numpy as np +import speclite.filters + from lenstronomy.SimulationAPI.ObservationConfig.LSST import LSST from lenstronomy.SimulationAPI.ObservationConfig.Roman import Roman from lenstronomy.SimulationAPI.ObservationConfig.Euclid import Euclid -import speclite.filters _OBSERVATORY_REGISTRY = {} @@ -10,6 +12,10 @@ LSST_BAND_LIST = ["u", "g", "r", "i", "z", "y"] EUCLID_BAND_LIST = ["VIS", "Y", "J", "H"] +# Ancillary bandpasses that are used by catalog sources but are not imaging +# observatories registered below. +_ADDITIONAL_BAND_EFFECTIVE_WAVELENGTH_MICRON = {"F814W": 0.805} + def check_speclite_name(band): """Checks if the raw band name is a valid speclite filter. @@ -33,6 +39,7 @@ def register_observatory( bands: list, speclite_fmt=check_speclite_name, sncosmo_fmt=None, + effective_wavelengths=None, ): """Register a new observatory to integrate it with image simulation tools. @@ -56,6 +63,10 @@ def register_observatory( :param sncosmo_fmt: A callable function that takes a ``band`` string and returns the corresponding sncosmo bandpass name. Set to ``None`` to use the raw band name as the sncosmo bandpass name. :type sncosmo_fmt: callable, optional + :param effective_wavelengths: Optional mapping from registered band names to + throughput-weighted effective wavelengths in microns. Use this only + when the responses are not available through ``speclite_fmt``. + :type effective_wavelengths: dict or None Given below is a simple example of how to define a custom observatory and register it using this function. A sophisticated example demonstrating full image simulation capabilities can be found at https://github.com/timedilatesme/MidEx-sims/blob/main/v1/lagn_sims.ipynb @@ -106,6 +117,7 @@ def kwargs_single_band(self): "bands": list(bands), "speclite_fmt": speclite_fmt, "sncosmo_fmt": sncosmo_fmt, + "effective_wavelengths": dict(effective_wavelengths or {}), } @@ -123,6 +135,16 @@ def kwargs_single_band(self): bands=ROMAN_BAND_LIST, speclite_fmt=lambda band: f"Roman-{band}", sncosmo_fmt=lambda band: f"{band}", + effective_wavelengths={ + "F062": 0.620, + "F087": 0.870, + "F106": 1.060, + "F129": 1.290, + "F146": 1.460, + "F158": 1.580, + "F184": 1.840, + "F213": 2.130, + }, ) register_observatory( name="Euclid", @@ -256,3 +278,54 @@ def get_all_supported_bands(): for info in _OBSERVATORY_REGISTRY.values(): all_bands.extend(info["bands"]) return all_bands + + +def get_band_effective_wavelength(band): + """Return the throughput-weighted effective wavelength of a band. + + The registered speclite response is used when available. Roman and + HST bands that are not shipped by speclite use explicitly configured + fallback values. Band-list position is deliberately not used as a + wavelength proxy. + + :param band: Imaging band name. + :type band: str + :return: Effective wavelength in microns. + :rtype: float + :raises ValueError: if the band is not registered. + """ + if band in _ADDITIONAL_BAND_EFFECTIVE_WAVELENGTH_MICRON: + return _ADDITIONAL_BAND_EFFECTIVE_WAVELENGTH_MICRON[band] + + obs_name = get_observatory(band) + configured_wavelengths = _OBSERVATORY_REGISTRY[obs_name]["effective_wavelengths"] + if band in configured_wavelengths: + return float(configured_wavelengths[band]) + + filter_name = get_speclite_filtername(band) + if filter_name is None: + raise ValueError( + f"Band '{band}' has neither an explicitly configured effective " + f"wavelength nor a registered speclite filter response." + ) + response = speclite.filters.load_filter(filter_name) + return float(response.effective_wavelength.to("micron").value) + + +def get_band_central_wavelength(band): + """Backward-compatible alias for :func:`get_band_effective_wavelength`.""" + return get_band_effective_wavelength(band) + + +def get_band_log_wavelength_ratio(band, reference_band): + """Return ``log(lambda_band / lambda_reference)``. + + Both wavelengths are throughput-weighted effective wavelengths. This + is the chromatic coordinate used by the local power-law SED + approximation. + """ + wavelength = get_band_effective_wavelength(band) + reference_wavelength = get_band_effective_wavelength(reference_band) + if wavelength <= 0 or reference_wavelength <= 0: + raise ValueError("Band effective wavelengths must be positive.") + return float(np.log(wavelength / reference_wavelength)) diff --git a/slsim/Lenses/lens.py b/slsim/Lenses/lens.py index cca7eabda..fd6e1d6b9 100644 --- a/slsim/Lenses/lens.py +++ b/slsim/Lenses/lens.py @@ -1388,8 +1388,8 @@ def lenstronomy_kwargs( field_galaxies_lens_model_list, kwargs_field_galaxies = ( self.field_galaxy_light_model_lenstronomy(band=band) ) - lens_light_model_list += field_galaxies_lens_model_list - kwargs_lens_light += kwargs_field_galaxies + lens_light_model_list.extend(field_galaxies_lens_model_list) + kwargs_lens_light.extend(kwargs_field_galaxies) kwargs_model = { "lens_light_model_list": lens_light_model_list, diff --git a/slsim/Sources/SourcePopulation/galaxies.py b/slsim/Sources/SourcePopulation/galaxies.py index 706097e4d..6b1221e06 100644 --- a/slsim/Sources/SourcePopulation/galaxies.py +++ b/slsim/Sources/SourcePopulation/galaxies.py @@ -381,78 +381,16 @@ def convert_catalog_to_source( kwargs_source["n_sersic"] = float(galaxy["n_sersic"]) if extended_source_type == "double_sersic": - if "e1_0" not in colnames or "e2_0" not in colnames: - if "ellipticity0" in colnames: - ellipticity0 = galaxy["ellipticity0"] - elif "a0" in colnames and "b0" in colnames: - axis_ratio_0 = axis_ratio(a=galaxy["a0"], b=galaxy["b0"]) - ellipticity0 = eccentricity(q=axis_ratio_0) - else: - raise ValueError( - "ellipticity or semi-major and semi-minor axis are missing for" - " the first light profile in galaxy_list columns" - ) - - e1_0, e2_0 = galaxy_projected_eccentricity( - float(ellipticity0), rotation_angle=phi_rot + kwargs_source.update( + _double_sersic_source_kwargs( + galaxy=galaxy, + colnames=colnames, + phi_rot=phi_rot, + size_model=size_model, + catalog_type=catalog_type, + cosmo=cosmo, ) - kwargs_source["e1_0"] = e1_0 - kwargs_source["e2_0"] = e2_0 - - else: - kwargs_source["e1_0"] = galaxy["e1_0"] - kwargs_source["e2_0"] = galaxy["e2_0"] - - if "e1_1" not in colnames or "e2_1" not in colnames: - if "ellipticity1" in colnames: - ellipticity1 = galaxy["ellipticity1"] - elif "a1" in colnames and "b1" in colnames: - axis_ratio_1 = axis_ratio(a=galaxy["a1"], b=galaxy["b1"]) - ellipticity1 = eccentricity(q=axis_ratio_1) - else: - raise ValueError( - "ellipticity or semi-major and semi-minor axis are missing for" - " the second light profile in galaxy_list columns" - ) - - e1_1, e2_1 = galaxy_projected_eccentricity( - float(ellipticity1), rotation_angle=phi_rot - ) - kwargs_source["e1_1"] = e1_1 - kwargs_source["e2_1"] = e2_1 - else: - kwargs_source["e1_1"] = galaxy["e1_1"] - kwargs_source["e2_1"] = galaxy["e2_1"] - if "angular_size_0" not in colnames or "angular_size_1" not in colnames: - if "a0" in colnames and "b0" in colnames: - kwargs_source["angular_size_0"] = average_angular_size( - a=galaxy["a0"], b=galaxy["b0"] - ) - else: - raise ValueError( - "semi-major and semi-minor axis are missing for the first light" - " profile in galaxy_list columns %s" % colnames - ) - if "a1" in colnames and "b1" in colnames: - kwargs_source["angular_size_1"] = average_angular_size( - a=galaxy["a1"], b=galaxy["b1"] - ) - else: - raise ValueError( - "semi-major and semi-minor axis are missing for the second" - " light profile in galaxy_list columns" - ) - else: - kwargs_source["angular_size_0"] = galaxy["angular_size_0"] - kwargs_source["angular_size_1"] = galaxy["angular_size_1"] - if "n_sersic_0" not in colnames or "n_sersic_1" not in colnames: - kwargs_source["n_sersic_0"] = 1 - kwargs_source["n_sersic_1"] = 4 - else: - kwargs_source["n_sersic_0"] = galaxy["n_sersic_0"] - kwargs_source["n_sersic_1"] = galaxy["n_sersic_1"] - kwargs_source["w0"] = galaxy["w0"] - kwargs_source["w1"] = galaxy["w1"] + ) if "vel_disp" in colnames: kwargs_source["vel_disp"] = float(galaxy["vel_disp"]) if "stellar_mass" in colnames: @@ -464,6 +402,164 @@ def convert_catalog_to_source( return kwargs_source +def _double_sersic_source_kwargs( + galaxy, + colnames, + phi_rot, + size_model=None, + catalog_type=None, + cosmo=None, +): + """Build the DoubleSersic and colour-gradient Source dictionary. + + The input ``galaxy`` dictionary may provide component-specific values or + single-component catalog values that are shared by both components: + + .. code-block:: python + + galaxy = { + # Required reference-band flux fractions + "w0": 0.4, + "w1": 0.6, + # Shape: component-specific, shared Cartesian, or shared scalar + "e1_0": 0.1, + "e2_0": 0.0, + "e1_1": 0.1, + "e2_1": 0.0, + # Size and Sersic-index component values + "angular_size_0": 0.3, + "angular_size_1": 0.9, + "n_sersic_0": 1.0, + "n_sersic_1": 4.0, + # Optional chromatic configuration + "color_gradient": { + "component_spectral_slopes": [0.5, -0.5], + "reference_band": "i", + "component_radius_factors": [0.5, 1.5], + "component_sersic_indices": [1.0, 4.0], + "min_weight": 1e-3, + }, + } + + When component sizes are absent, ``angular_size`` is multiplied by + ``color_gradient['component_radius_factors']``. When component Sersic + indices are absent, ``n_sersic`` and + ``color_gradient['component_sersic_indices']`` provide the defaults. + + :param galaxy: Catalog row or dictionary containing light-profile fields. + :param colnames: Available keys/column names in ``galaxy``. + :param phi_rot: Position angle in radians used for scalar ellipticities. + :param size_model: Galaxy size model passed to :func:`_galaxy_size`. + :param catalog_type: Optional catalog convention identifier. + :param cosmo: Astropy cosmology used to convert physical/angular size. + :return: Dictionary containing ``e1_0``, ``e2_0``, ``e1_1``, ``e2_1``, + ``angular_size_0``, ``angular_size_1``, ``n_sersic_0``, + ``n_sersic_1``, ``w0``, ``w1``, and optional ``color_gradient``. + """ + kwargs_double_sersic = {} + color_gradient = galaxy["color_gradient"] if "color_gradient" in colnames else {} + if not isinstance(color_gradient, dict): + raise ValueError( + "galaxy['color_gradient'] must be a dictionary when constructing " + f"a DoubleSersic source; received {color_gradient!r}." + ) + + for component in (0, 1): + e1_key, e2_key = f"e1_{component}", f"e2_{component}" + if e1_key in colnames and e2_key in colnames: + e1, e2 = galaxy[e1_key], galaxy[e2_key] + elif "e1" in colnames and "e2" in colnames: + e1, e2 = galaxy["e1"], galaxy["e2"] + elif "e1_light" in colnames and "e2_light" in colnames: + e1, e2 = galaxy["e1_light"], galaxy["e2_light"] + else: + ellipticity_key = f"ellipticity{component}" + a_key, b_key = f"a{component}", f"b{component}" + if ellipticity_key in colnames: + ellipticity = galaxy[ellipticity_key] + elif a_key in colnames and b_key in colnames: + ellipticity = eccentricity( + q=axis_ratio(a=galaxy[a_key], b=galaxy[b_key]) + ) + elif "ellipticity" in colnames or "e" in colnames: + ellipticity = ( + galaxy["ellipticity"] if "ellipticity" in colnames else galaxy["e"] + ) + else: + raise ValueError( + f"Cannot determine ellipticity for DoubleSersic component " + f"{component}; available galaxy_list columns are {colnames}." + ) + e1, e2 = galaxy_projected_eccentricity( + float(ellipticity), rotation_angle=phi_rot + ) + kwargs_double_sersic[e1_key] = e1 + kwargs_double_sersic[e2_key] = e2 + + if "angular_size_0" in colnames and "angular_size_1" in colnames: + kwargs_double_sersic["angular_size_0"] = galaxy["angular_size_0"] + kwargs_double_sersic["angular_size_1"] = galaxy["angular_size_1"] + elif all(key in colnames for key in ("a0", "b0", "a1", "b1")): + kwargs_double_sersic["angular_size_0"] = average_angular_size( + a=galaxy["a0"], b=galaxy["b0"] + ) + kwargs_double_sersic["angular_size_1"] = average_angular_size( + a=galaxy["a1"], b=galaxy["b1"] + ) + elif "angular_size" in colnames or catalog_type is not None: + angular_size, _ = _galaxy_size( + galaxy, + size_model=size_model, + catalog_type=catalog_type, + cosmo=cosmo, + ) + radius_factors = color_gradient.get("component_radius_factors", (0.5, 1.5)) + if len(radius_factors) != 2: + raise ValueError( + "color_gradient['component_radius_factors'] must contain two values." + ) + kwargs_double_sersic["angular_size_0"] = angular_size * float(radius_factors[0]) + kwargs_double_sersic["angular_size_1"] = angular_size * float(radius_factors[1]) + else: + raise ValueError( + "Cannot determine DoubleSersic component sizes. Provide one of: " + "(1) both 'angular_size_0' and 'angular_size_1'; " + "(2) all four axis fields 'a0', 'b0', 'a1', and 'b1'; " + "(3) a shared 'angular_size'; or " + "(4) a catalog_type from which the shared size can be inferred. " + f"Available galaxy_list columns are {colnames}." + ) + + if "n_sersic_0" in colnames and "n_sersic_1" in colnames: + indices = (galaxy["n_sersic_0"], galaxy["n_sersic_1"]) + else: + if "n_sersic" in colnames: + n_sersic = float(galaxy["n_sersic"]) + elif "galaxy_type" in colnames and galaxy["galaxy_type"] == "red": + n_sersic = 4.0 + else: + n_sersic = 1.0 + indices = color_gradient.get("component_sersic_indices", (n_sersic, n_sersic)) + if len(indices) != 2: + raise ValueError( + "color_gradient['component_sersic_indices'] must contain two values." + ) + kwargs_double_sersic["n_sersic_0"] = float(indices[0]) + kwargs_double_sersic["n_sersic_1"] = float(indices[1]) + + missing_weights = [key for key in ("w0", "w1") if key not in colnames] + if missing_weights: + raise ValueError( + "DoubleSersic source dictionary is missing reference-band flux " + f"weights {missing_weights}; provide both 'w0' and 'w1'." + ) + kwargs_double_sersic["w0"] = galaxy["w0"] + kwargs_double_sersic["w1"] = galaxy["w1"] + if "color_gradient" in colnames: + kwargs_double_sersic["color_gradient"] = color_gradient + return kwargs_double_sersic + + def down_sample_to_dc2(galaxy_pop, sky_area): """Downsamples given galaxy pop above redshift 1.5 to DC2 galaxy population. diff --git a/slsim/Sources/SourceTypes/catalog_source.py b/slsim/Sources/SourceTypes/catalog_source.py index c56ab8292..773779a7e 100644 --- a/slsim/Sources/SourceTypes/catalog_source.py +++ b/slsim/Sources/SourceTypes/catalog_source.py @@ -1,7 +1,9 @@ +from slsim.Sources.SourceTypes.double_sersic import DoubleSersic from slsim.Sources.SourceTypes.single_sersic import SingleSersic from slsim.Sources.SourceTypes.source_base import SourceBase from slsim.Sources.SourceCatalogues.CosmosWebCatalog import galaxy_match as CosmosWeb from slsim.Sources.SourceCatalogues.HSTCosmosCatalog import galaxy_match as HSTCosmos +from slsim.Util.color_gradient import radial_color_gradient_image from lenstronomy.Util.param_util import ellipticity2phi_q CATALOG_TYPES = ["HST_COSMOS, COSMOS_WEB"] @@ -26,6 +28,9 @@ def __init__( max_scale=1, match_n_sersic=False, sersic_fallback=False, + band_dependent_color_gradient=False, + color_gradient=None, + fallback_double_sersic_kwargs=None, **source_dict, ): """ @@ -54,6 +59,18 @@ def __init__( :type match_n_sersic: bool :param sersic_fallback: If the matching process returns no matches, then fall back on a single sersic profile. :type sersic_fallback: bool + :param band_dependent_color_gradient: If True, apply an opt-in radial + colour-gradient transfer to matched HST_COSMOS images. Failed matches + fall back to a DoubleSersic model with the same ``color_gradient``. + :type band_dependent_color_gradient: bool + :param color_gradient: Dictionary containing colour-gradient settings. + Matched HST_COSMOS images use ``grad_color`` (mag/dex) with + ``reference_band`` defaulting to ``F814W``. DoubleSersic fallback uses + ``component_spectral_slopes``. + :type color_gradient: dict or None + :param fallback_double_sersic_kwargs: Optional overrides for the + DoubleSersic parameters used after a failed HST_COSMOS match. + :type fallback_double_sersic_kwargs: dict or None """ super().__init__(extended_source=True, point_source=False, **source_dict) self.name = "GAL" @@ -65,6 +82,9 @@ def __init__( self._max_scale = max_scale self._match_n_sersic = match_n_sersic self._sersic_fallback = sersic_fallback + self._band_dependent_color_gradient = band_dependent_color_gradient + self._color_gradient = color_gradient + self._fallback_double_sersic_kwargs = fallback_double_sersic_kwargs self.source_dict = source_dict # Process catalog and store as class attribute @@ -93,6 +113,27 @@ def __init__( f"Catalog_type {catalog_type} not supported. Currently only {CATALOG_TYPES} are supported." ) + if self._band_dependent_color_gradient: + if catalog_type != "HST_COSMOS": + raise ValueError( + "band_dependent_color_gradient is currently supported only " + "for catalog_type='HST_COSMOS'; received " + f"catalog_type={catalog_type!r}." + ) + if not isinstance(self._color_gradient, dict): + raise ValueError( + "color_gradient must be a dictionary when " + "band_dependent_color_gradient is enabled; received " + f"{self._color_gradient!r} (type " + f"{type(self._color_gradient).__name__})." + ) + if self._fallback_double_sersic_kwargs is not None and not isinstance( + self._fallback_double_sersic_kwargs, dict + ): + raise ValueError( + "fallback_double_sersic_kwargs must be a dictionary or None." + ) + self._catalog_type = catalog_type self._catalog_path = catalog_path @@ -146,8 +187,10 @@ def kwargs_extended_light(self, band=None): match_n_sersic=self._match_n_sersic, ) ) - # If the matching failed, fall back on a regular sersic profile + # If matching fails, the optional chromatic mode uses DoubleSersic. if self._image_list is None: + if self._band_dependent_color_gradient: + return self._double_sersic_fallback().kwargs_extended_light(band=band) if self._sersic_fallback: if not hasattr(self, "single_sersic"): self.single_sersic = SingleSersic( @@ -158,11 +201,10 @@ def kwargs_extended_light(self, band=None): **self.source_dict, ) return self.single_sersic.kwargs_extended_light(band=band) - else: - raise ValueError( - "No valid matches found! Try reducing the desired angular size or increasing max_scale." - "Alternatively, enable sersic_fallback to use a single sersic whenever the matching fails." - ) + raise ValueError( + "No valid matches found! Try reducing the desired angular size or increasing max_scale." + "Alternatively, enable sersic_fallback to use a single sersic whenever the matching fails." + ) if band is None: mag_source = 1 @@ -170,7 +212,7 @@ def kwargs_extended_light(self, band=None): mag_source = self.extended_source_magnitude(band=band) center_source = self.extended_source_position - image = self._select_image_from_band(band) + image = self._image_for_band(band) light_model_list = ["INTERPOL"] kwargs_extended_source = [ @@ -185,6 +227,51 @@ def kwargs_extended_light(self, band=None): ] return light_model_list, kwargs_extended_source + def _image_for_band(self, band): + """Return the catalog image, optionally with HST chromatic + morphology.""" + if self._band_dependent_color_gradient: + image = self._image_list[0] + else: + image = self._select_image_from_band(band) + + if not self._band_dependent_color_gradient or band is None: + return image + + return radial_color_gradient_image( + image=image, + band=band, + color_gradient=self._color_gradient, + angular_size=self.angular_size, + pixel_scale=self._scale, + default_reference="F814W", + ) + + def _double_sersic_fallback(self): + """Build the chromatic fallback model after a failed HST match.""" + if hasattr(self, "double_sersic"): + return self.double_sersic + + fallback_kwargs = { + "angular_size_0": 0.5 * self.angular_size, + "angular_size_1": self.angular_size, + "n_sersic_0": 4.0, + "n_sersic_1": 1.0, + "w0": 0.4, + "w1": 0.6, + "e1_0": self._e1, + "e2_0": self._e2, + "e1_1": self._e1, + "e2_1": self._e2, + } + fallback_kwargs.update(self._fallback_double_sersic_kwargs or {}) + fallback_kwargs["color_gradient"] = self._color_gradient + self.double_sersic = DoubleSersic( + **fallback_kwargs, + **self.source_dict, + ) + return self.double_sersic + def _select_image_from_band(self, band): """Selects an image based off of the input band. Only relevant for source catalogs that provide images for multiple bands. diff --git a/slsim/Sources/SourceTypes/double_sersic.py b/slsim/Sources/SourceTypes/double_sersic.py index d297c2343..0e4921558 100644 --- a/slsim/Sources/SourceTypes/double_sersic.py +++ b/slsim/Sources/SourceTypes/double_sersic.py @@ -2,6 +2,7 @@ from slsim.Sources.SourceTypes.source_base import SourceBase from slsim.Util.param_util import ellipticity_slsim_to_lenstronomy from slsim.Util.param_util import surface_brightness_reff +from slsim.Util.color_gradient import component_weights_for_band class DoubleSersic(SourceBase): @@ -19,6 +20,7 @@ def __init__( e2_0=0, e1_1=0, e2_1=0, + color_gradient=None, **source_dict ): """ @@ -33,6 +35,10 @@ def __init__( :param e2_1: eccentricity component of second Sersic :param w0: flux weight of first Sersic component :param w1: flux weight of second Sersic component, if =None, will be set w1 = 1 - w0, otherwise it has to match. + :param color_gradient: Optional dictionary defining a band-dependent + two-component colour gradient with lightweight SED slopes. Supported + keys are ``component_spectral_slopes``, ``reference_band``, and + ``min_weight``. Components with larger spectral slopes are redder. :param source_dict: dictionary for SourceBase() option (see documentation) :type source_dict: dict or astropy.table.Table @@ -50,6 +56,7 @@ def __init__( w1 = 1 - w0 assert np.isclose(w0 + w1, 1, rtol=1e-3) self._w1 = w1 + self._color_gradient = color_gradient self._light_model_list = [ "SERSIC_ELLIPSE", @@ -118,8 +125,9 @@ def kwargs_extended_light(self, band=None): center_source = self.extended_source_position # compute magnitude for each sersic component based on weight flux = 10 ** (-mag_source / 2.5) - mag_source0 = -2.5 * np.log10(self._w0 * flux) - mag_source1 = -2.5 * np.log10(self._w1 * flux) + w0, w1 = self._weights_for_band(band) + mag_source0 = -2.5 * np.log10(w0 * flux) + mag_source1 = -2.5 * np.log10(w1 * flux) # convert from slsim to lenstronomy convention. e1_light_source_1_lenstronomy, e2_light_source_1_lenstronomy = ( ellipticity_slsim_to_lenstronomy( @@ -154,7 +162,17 @@ def kwargs_extended_light(self, band=None): "center_y": center_source[1], }, ] - return self._light_model_list, kwargs_extended_source + return self._light_model_list.copy(), kwargs_extended_source + + def _weights_for_band(self, band): + """Return Sersic component weights for an imaging band.""" + return component_weights_for_band( + base_weights=(self._w0, self._w1), + band=band, + color_gradient=self._color_gradient, + source_dict=self.source_dict, + default_reference="i", + ) def _shape_light_model(self): """ diff --git a/slsim/Util/color_gradient.py b/slsim/Util/color_gradient.py new file mode 100644 index 000000000..5f927725d --- /dev/null +++ b/slsim/Util/color_gradient.py @@ -0,0 +1,212 @@ +import numpy as np + + +def default_reference_band(source_dict, default="i"): + """Select the available magnitude band nearest a preferred reference band. + + Bands are inferred from ``mag_`` keys in ``source_dict``. The band + closest to ``default`` in logarithmic effective-wavelength space is + returned, matching the model's power-law wavelength dependence. If no + magnitude band is available, ``default`` is returned unchanged. + + :param source_dict: Source parameters that may contain magnitude entries + such as ``mag_g``, ``mag_i``, or ``mag_F158``. + :type source_dict: dict or mapping + :param default: Preferred reference band and fallback when no magnitude + bands are available. + :type default: str + :return: Available band with the smallest absolute + ``log(lambda_band / lambda_default)``, or ``default`` when none exists. + :rtype: str + """ + from slsim.ImageSimulation.image_quality_lenstronomy import ( + get_band_log_wavelength_ratio, + ) + + available_bands = [ + key.replace("mag_", "", 1) + for key in source_dict + if isinstance(key, str) and key.startswith("mag_") + ] + if not available_bands: + return default + + positions = [ + abs(get_band_log_wavelength_ratio(band=band, reference_band=default)) + for band in available_bands + ] + return available_bands[int(np.argmin(positions))] + + +def component_weights_for_band( + base_weights, + band, + color_gradient=None, + source_dict=None, + default_reference="i", +): + """Return band-dependent component weights from local SED slopes. + + This is a lightweight, effective-wavelength approximation to chromatic + light components. Instead of integrating a full stellar-population SED + through each bandpass, each component is assigned a local power-law SED, + ``S_k(lambda) proportional lambda**alpha_k``, evaluated at the central + wavelength of the requested band. The reference-band component weights are + then reweighted as + + ``w_k(b) = w_k(ref) * (lambda_b/lambda_ref)**alpha_k / normalization``. + + The approximation is intended to introduce controlled colour gradients in + analytic multi-component light profiles; it is not a replacement for a + stellar population synthesis model. + + Ref: + Hogg et al. 2002, "The K correction", astro-ph/0210394: + broadband fluxes are formally filter-response weighted SED integrals; + this function uses the corresponding effective-wavelength limit. + Conroy 2013, ARA&A, 51, 393: + review of full stellar-population SED modelling, useful context for + what is intentionally omitted by this lightweight approximation. + La Barbera et al. 2005, MNRAS, 358, 1116; La Barbera & de Carvalho 2009, + ApJ, 699, L76: + observational motivation for radial colour gradients in galaxies. + """ + weights = np.asarray(base_weights, dtype=float) + weights = weights / np.sum(weights) + + if band is None or color_gradient is None: + return tuple(float(weight) for weight in weights) + if not isinstance(color_gradient, dict): + raise ValueError("color_gradient must be a dictionary or None.") + + slopes = color_gradient.get( + "component_spectral_slopes", color_gradient.get("sed_slopes") + ) + if slopes is None: + return tuple(float(weight) for weight in weights) + + slopes = np.asarray(slopes, dtype=float) + if slopes.shape != weights.shape: + raise ValueError( + "color_gradient['component_spectral_slopes'] must match the " + "number of components." + ) + + reference_band = color_gradient.get("reference_band") + if reference_band is None: + reference_band = default_reference_band( + source_dict or {}, default=default_reference + ) + + min_weight = float(color_gradient.get("min_weight", 1e-4)) + if not 0 <= min_weight < 1 / len(weights): + raise ValueError( + "color_gradient['min_weight'] must be in [0, 1 / n_components)." + ) + + from slsim.ImageSimulation.image_quality_lenstronomy import ( + get_band_effective_wavelength, + ) + + wavelength = get_band_effective_wavelength(band) + reference_wavelength = get_band_effective_wavelength(reference_band) + if reference_wavelength <= 0: + raise ValueError("The reference band wavelength must be positive.") + + sed_factors = (wavelength / reference_wavelength) ** slopes + sed_weights = weights * sed_factors + sed_weights = sed_weights / np.sum(sed_weights) + sed_weights = np.clip(sed_weights, min_weight, 1 - min_weight) + sed_weights = sed_weights / np.sum(sed_weights) + return tuple(float(weight) for weight in sed_weights) + + +def attach_foreground_deflector_color_gradient( + galaxy_table, + color_gradient, + component_weights=(0.4, 0.6), +): + """Attach opt-in foreground colour-gradient columns to a deflector table. + + The resulting columns are consumed by the standard + ``Source``/``DoubleSersic`` light-model path to split foreground light into + two chromatic Sersic components. + The operation is in-place and returns ``galaxy_table`` for convenience. + + :param galaxy_table: galaxy/deflector table to annotate + :param color_gradient: dictionary with ``component_spectral_slopes`` and + optional foreground component settings + :param component_weights: two reference-band flux weights for the Sersic + components + :return: annotated galaxy table + """ + if color_gradient is None: + return galaxy_table + if not isinstance(color_gradient, dict): + raise ValueError("color_gradient must be a dictionary or None.") + + weights = np.asarray(component_weights, dtype=float) + if weights.shape != (2,): + raise ValueError("component_weights must contain two values.") + if np.any(weights < 0) or np.sum(weights) <= 0: + raise ValueError("component_weights must be non-negative with positive sum.") + weights = weights / np.sum(weights) + + galaxy_table["color_gradient"] = [ + dict(color_gradient) for _ in range(len(galaxy_table)) + ] + galaxy_table["w0"] = np.full(len(galaxy_table), weights[0]) + galaxy_table["w1"] = np.full(len(galaxy_table), weights[1]) + return galaxy_table + + +def radial_color_gradient_image( + image, + band, + color_gradient, + angular_size, + pixel_scale, + default_reference="F814W", +): + """Apply a d(color)/dlog10(r) gradient to an image and preserve flux. + + See https://arxiv.org/pdf/1006.4056 for details. + """ + if band is None or color_gradient is None: + return image + if not isinstance(color_gradient, dict): + raise ValueError("color_gradient must be a dictionary or None.") + + grad_color = float( + color_gradient.get("grad_color", color_gradient.get("gradient", 0.0)) + ) + if grad_color == 0: + return image + + reference_band = color_gradient.get("reference_band") or default_reference + from slsim.ImageSimulation.image_quality_lenstronomy import ( + get_band_log_wavelength_ratio, + ) + + band_offset = get_band_log_wavelength_ratio( + band=band, reference_band=reference_band + ) + if band_offset == 0: + return image + + image = np.asarray(image, dtype=float) + y_grid, x_grid = np.indices(image.shape, dtype=float) + center_x = (image.shape[1] - 1) / 2 + center_y = (image.shape[0] - 1) / 2 + radius = np.hypot(x_grid - center_x, y_grid - center_y) + half_light_radius_pixels = angular_size / pixel_scale + radius = np.maximum(radius, 0.5) + radius_ratio = radius / max(half_light_radius_pixels, 0.5) + + delta_mag = band_offset * grad_color * np.log10(radius_ratio) + chromatic_image = image * 10 ** (-0.4 * delta_mag) + original_flux = np.sum(image) + chromatic_flux = np.sum(chromatic_image) + if chromatic_flux != 0: + chromatic_image *= original_flux / chromatic_flux + return chromatic_image diff --git a/tests/test_Deflectors/test_DeflectorPopulation/test_galaxy_deflectors.py b/tests/test_Deflectors/test_DeflectorPopulation/test_galaxy_deflectors.py index 9477229e4..29ba00ff1 100644 --- a/tests/test_Deflectors/test_DeflectorPopulation/test_galaxy_deflectors.py +++ b/tests/test_Deflectors/test_DeflectorPopulation/test_galaxy_deflectors.py @@ -7,6 +7,7 @@ ) from slsim.Pipelines.skypy_pipeline import SkyPyPipeline from astropy.units import Quantity +from astropy.table import Table import numpy as np import pytest import copy @@ -171,5 +172,80 @@ def test_elliptical_galaxies(): assert galaxy_class1.draw_deflector().mass_properties["gamma_pl"] == 2.05 +def test_galaxy_deflectors_foreground_color_gradient(): + red_galaxies = foreground_test_galaxy_table() + blue_galaxies = foreground_test_galaxy_table() + kwargs_deflector_cut = {} + kwargs_mass2light = {} + cosmo = FlatLambdaCDM(H0=70, Om0=0.3) + sky_area = Quantity(value=0.05, unit="deg2") + foreground_color_gradient = { + "component_spectral_slopes": [2.0, -1.0], + "reference_band": "i", + } + + galaxy_class = GalaxyDeflectors( + red_galaxies, + blue_galaxy_list=blue_galaxies, + kwargs_cut=kwargs_deflector_cut, + kwargs_mass2light=kwargs_mass2light, + cosmo=cosmo, + sky_area=sky_area, + catalog_type=None, + foreground_color_gradient=foreground_color_gradient, + foreground_component_weights=(0.4, 0.6), + ) + assert "color_gradient" not in red_galaxies.colnames + assert "color_gradient" not in blue_galaxies.colnames + deflector = galaxy_class.draw_deflector() + model_list, kwargs_i = deflector.light_model_lenstronomy(band="i") + _, kwargs_g = deflector.light_model_lenstronomy(band="g") + _, kwargs_y = deflector.light_model_lenstronomy(band="y") + + assert model_list == ["SERSIC_ELLIPSE", "SERSIC_ELLIPSE"] + assert kwargs_i[0]["R_sersic"] < kwargs_i[1]["R_sersic"] + assert component_flux_fraction(kwargs_i) == pytest.approx(0.4) + assert component_flux_fraction(kwargs_y) > component_flux_fraction(kwargs_g) + + +def test_foreground_color_gradient_rejects_unsupported_light_type(): + with pytest.raises(ValueError, match="requires a single_sersic or double_sersic"): + GalaxyDeflectors( + foreground_test_galaxy_table(), + kwargs_mass2light={}, + cosmo=FlatLambdaCDM(H0=70, Om0=0.3), + sky_area=Quantity(value=0.05, unit="deg2"), + catalog_type=None, + light_type="catalog_source", + foreground_color_gradient={"component_spectral_slopes": [2.0, -1.0]}, + ) + + +def foreground_test_galaxy_table(): + return Table( + { + "z": [0.2, 0.3], + "stellar_mass": [10**11, 2 * 10**11], + "angular_size": [0.7, 0.8], + "ellipticity": [0.2, 0.25], + "mag_i": [19.0, 20.0], + "mag_g": [20.0, 21.0], + "mag_y": [18.0, 19.0], + "e1_light": [0.1, 0.1], + "e2_light": [0.0, 0.0], + "e1_mass": [0.1, 0.1], + "e2_mass": [0.0, 0.0], + "n_sersic": [4.0, 4.0], + "vel_disp": [200.0, 210.0], + } + ) + + +def component_flux_fraction(kwargs_light): + flux0 = 10 ** (-kwargs_light[0]["magnitude"] / 2.5) + flux1 = 10 ** (-kwargs_light[1]["magnitude"] / 2.5) + return flux0 / (flux0 + flux1) + + if __name__ == "__main__": pytest.main() diff --git a/tests/test_Deflectors/test_MassTypes/test_epl_sersic.py b/tests/test_Deflectors/test_MassTypes/test_epl_sersic.py index e16fd963e..6f68575d8 100644 --- a/tests/test_Deflectors/test_MassTypes/test_epl_sersic.py +++ b/tests/test_Deflectors/test_MassTypes/test_epl_sersic.py @@ -1,5 +1,4 @@ import pytest - from slsim.Deflectors.MassTypes.epl import EPL from slsim.Sources.source import Source from astropy.cosmology import FlatLambdaCDM diff --git a/tests/test_Deflectors/test_deflector.py b/tests/test_Deflectors/test_deflector.py index f190f7753..c994cca48 100644 --- a/tests/test_Deflectors/test_deflector.py +++ b/tests/test_Deflectors/test_deflector.py @@ -152,6 +152,8 @@ def test_light_model_lenstronomy(self): light_model, kwargs_lens_light = self.deflector.light_model_lenstronomy( band=band ) + assert isinstance(light_model, list) + assert isinstance(kwargs_lens_light, list) assert light_model[0] == "SERSIC_ELLIPSE" assert kwargs_lens_light[0]["R_sersic"] == 0.5 diff --git a/tests/test_Deflectors/test_deflector_util.py b/tests/test_Deflectors/test_deflector_util.py index 48aed7cb3..86374094a 100644 --- a/tests/test_Deflectors/test_deflector_util.py +++ b/tests/test_Deflectors/test_deflector_util.py @@ -70,3 +70,46 @@ def test_light2mass(): kwargs_source, mass_type="EPL", halo_dict={"theta_E": 1, "gamma_pl": 2} ) npt.assert_almost_equal(kwargs_mass_epl["theta_E"], 1, decimal=5) + + +def test_light2mass_uses_effective_double_sersic_ellipticity(): + kwargs_mass = deflector_util.light2mass( + { + "e1_0": 0.2, + "e2_0": -0.1, + "e1_1": -0.1, + "e2_1": 0.3, + "w0": 0.25, + "w1": 0.75, + }, + mass_type="EPL", + halo_dict={"theta_E": 1, "gamma_pl": 2}, + light2mass_e_scaling=1, + light2mass_e_scatter=0, + ) + + npt.assert_almost_equal(kwargs_mass["e1"], -0.025) + npt.assert_almost_equal(kwargs_mass["e2"], 0.2) + + +def test_light2mass_rejects_invalid_double_sersic_light_shapes(): + common = { + "e1_0": 0.2, + "e2_0": -0.1, + "e1_1": -0.1, + "e2_1": 0.3, + } + + with npt.assert_raises_regex(ValueError, "must be non-negative"): + deflector_util.light2mass( + {**common, "w0": -0.25, "w1": 0.75}, + mass_type="EPL", + halo_dict={"theta_E": 1, "gamma_pl": 2}, + ) + + with npt.assert_raises_regex(ValueError, "Cannot derive mass ellipticity"): + deflector_util.light2mass( + {"e1_0": 0.2, "e2_0": -0.1}, + mass_type="EPL", + halo_dict={"theta_E": 1, "gamma_pl": 2}, + ) diff --git a/tests/test_ImageSimulation/test_image_quality_lenstronomy.py b/tests/test_ImageSimulation/test_image_quality_lenstronomy.py index cb161679d..363f118a6 100644 --- a/tests/test_ImageSimulation/test_image_quality_lenstronomy.py +++ b/tests/test_ImageSimulation/test_image_quality_lenstronomy.py @@ -1,4 +1,5 @@ import pytest +import numpy as np from slsim.ImageSimulation.image_quality_lenstronomy import ( check_speclite_name, get_speclite_filtername, @@ -8,6 +9,9 @@ get_observatory, register_observatory, get_all_supported_bands, + get_band_central_wavelength, + get_band_effective_wavelength, + get_band_log_wavelength_ratio, ) @@ -248,5 +252,50 @@ def test_get_all_supported_bands_contains_defaults(): assert band in all_bands +def test_default_band_wavelength_ordering(): + assert get_band_effective_wavelength("g") == pytest.approx( + get_band_central_wavelength("g") + ) + assert get_band_central_wavelength("g") < get_band_central_wavelength("i") + assert get_band_central_wavelength("F106") < get_band_central_wavelength("F184") + assert get_band_central_wavelength("VIS") < get_band_central_wavelength("H") + assert get_band_log_wavelength_ratio("F184", reference_band="g") > 0 + + +def test_band_wavelength_helpers_cover_hst_and_custom_registry_fallback(): + assert get_band_central_wavelength("F814W") == pytest.approx(0.805) + assert get_band_log_wavelength_ratio("F814W", reference_band="F814W") == 0 + + register_observatory( + name="WavelengthTestObs", + observatory_class=DummyObservatory, + bands=["bessell-B", "bessell-V", "bessell-R"], + ) + assert get_band_central_wavelength("bessell-B") < get_band_central_wavelength( + "bessell-R" + ) + assert get_band_log_wavelength_ratio( + "bessell-R", reference_band="bessell-B" + ) == pytest.approx( + np.log( + get_band_central_wavelength("bessell-R") + / get_band_central_wavelength("bessell-B") + ) + ) + + with pytest.raises(ValueError, match="not recognised"): + get_band_central_wavelength("UnknownBand") + + register_observatory( + name="NonPositiveWavelengthTestObs", + observatory_class=DummyObservatory, + bands=["NPW-band", "positive-band"], + speclite_fmt=None, + effective_wavelengths={"NPW-band": 0.0, "positive-band": 1.0}, + ) + with pytest.raises(ValueError, match="effective wavelengths must be positive"): + get_band_log_wavelength_ratio("NPW-band", reference_band="positive-band") + + if __name__ == "__main__": pytest.main() diff --git a/tests/test_Lenses/test_lens.py b/tests/test_Lenses/test_lens.py index b000454db..bc57692bd 100644 --- a/tests/test_Lenses/test_lens.py +++ b/tests/test_Lenses/test_lens.py @@ -854,6 +854,18 @@ def test_field_galaxies_integration(self): kwargs_params_no_fg["kwargs_lens_light"] ) + # Repeated calls in different bands must not accumulate field-galaxy + # models in the deflector's internal light-model list. + kwargs_model_fg_r, kwargs_params_fg_r = lens_with_fg.lenstronomy_kwargs( + band="r" + ) + assert len(kwargs_model_fg_r["lens_light_model_list"]) == len( + kwargs_model_fg["lens_light_model_list"] + ) + assert len(kwargs_params_fg_r["kwargs_lens_light"]) == len( + kwargs_model_fg_r["lens_light_model_list"] + ) + # Test C: Ensure field_galaxy_light_model_lenstronomy correctly handles None (fallback) empty_model_list, empty_kwargs_list = ( self.gg_lens.field_galaxy_light_model_lenstronomy(band="i") diff --git a/tests/test_Sources/test_SourcePopulation/test_galaxies.py b/tests/test_Sources/test_SourcePopulation/test_galaxies.py index 78411e844..038e39307 100644 --- a/tests/test_Sources/test_SourcePopulation/test_galaxies.py +++ b/tests/test_Sources/test_SourcePopulation/test_galaxies.py @@ -402,6 +402,126 @@ def test_galaxy_projected_eccentricity(): assert e2 == 0 +def test_double_sersic_catalog_single_component_fallbacks(): + """Single-component catalog fields can seed both Sersic components.""" + common = { + "z": 0.5, + "w0": 0.4, + "w1": 0.6, + "angular_size_0": 0.2, + "angular_size_1": 0.6, + "n_sersic_0": 1.0, + "n_sersic_1": 4.0, + } + + cartesian = convert_catalog_to_source( + {**common, "e1": 0.1, "e2": -0.2}, + extended_source_type="double_sersic", + catalog_type=None, + ) + assert cartesian["e1_1"] == pytest.approx(0.1) + assert cartesian["e2_1"] == pytest.approx(-0.2) + + projected = convert_catalog_to_source( + {**common, "ellipticity": 0.2, "a_rot": 0.0}, + extended_source_type="double_sersic", + catalog_type=None, + ) + assert projected["e1_1"] == pytest.approx(projected["e1_0"]) + assert projected["e2_1"] == pytest.approx(projected["e2_0"]) + + +def test_double_sersic_catalog_component_defaults_and_validation(): + common = { + "z": 0.5, + "w0": 0.4, + "w1": 0.6, + "e1": 0.1, + "e2": 0.0, + } + + red_galaxy = convert_catalog_to_source( + { + **common, + "angular_size_0": 0.2, + "angular_size_1": 0.6, + "galaxy_type": "red", + }, + extended_source_type="double_sersic", + catalog_type=None, + ) + assert red_galaxy["n_sersic_0"] == pytest.approx(4.0) + assert red_galaxy["n_sersic_1"] == pytest.approx(4.0) + + with pytest.raises(ValueError, match="component_radius_factors"): + convert_catalog_to_source( + { + **common, + "angular_size": 0.5, + "n_sersic_0": 1.0, + "n_sersic_1": 4.0, + "color_gradient": {"component_radius_factors": [1.0]}, + }, + extended_source_type="double_sersic", + catalog_type=None, + cosmo=FlatLambdaCDM(H0=70, Om0=0.3), + ) + + with pytest.raises( + ValueError, match="Cannot determine DoubleSersic component sizes" + ): + convert_catalog_to_source( + { + **common, + "n_sersic_0": 1.0, + "n_sersic_1": 4.0, + }, + extended_source_type="double_sersic", + catalog_type=None, + ) + + with pytest.raises(ValueError, match="component_sersic_indices"): + convert_catalog_to_source( + { + **common, + "angular_size_0": 0.2, + "angular_size_1": 0.6, + "color_gradient": {"component_sersic_indices": [1.0]}, + }, + extended_source_type="double_sersic", + catalog_type=None, + ) + + with pytest.raises(ValueError, match="missing reference-band flux weights"): + convert_catalog_to_source( + { + "z": 0.5, + "e1": 0.1, + "e2": 0.0, + "angular_size_0": 0.2, + "angular_size_1": 0.6, + "n_sersic_0": 1.0, + "n_sersic_1": 4.0, + }, + extended_source_type="double_sersic", + catalog_type=None, + ) + + with pytest.raises(ValueError, match="color_gradient.*must be a dictionary"): + convert_catalog_to_source( + { + **common, + "angular_size_0": 0.2, + "angular_size_1": 0.6, + "n_sersic_0": 1.0, + "n_sersic_1": 4.0, + "color_gradient": "invalid", + }, + extended_source_type="double_sersic", + catalog_type=None, + ) + + def test_down_sample_to_dc2(): galaxy_pop = Table( { diff --git a/tests/test_Sources/test_SourceTypes/test_catalog_source.py b/tests/test_Sources/test_SourceTypes/test_catalog_source.py index 3f806e970..3cfb849b7 100644 --- a/tests/test_Sources/test_SourceTypes/test_catalog_source.py +++ b/tests/test_Sources/test_SourceTypes/test_catalog_source.py @@ -10,6 +10,7 @@ from slsim.Pipelines import SkyPyPipeline from slsim.Sources.SourcePopulation.galaxies import Galaxies from slsim.Sources.SourceTypes.single_sersic import SingleSersic +from slsim.Sources.SourceTypes.double_sersic import DoubleSersic from slsim.Sources.SourceTypes.catalog_source import CatalogSource from slsim.Sources.source import Source from slsim.Deflectors.deflector import Deflector @@ -134,6 +135,108 @@ def test_select_image_from_band(self): band="wrong", ) + def test_hst_band_dependent_color_gradient(self): + source_dict = dict(self.source1.source_dict) + source = CatalogSource( + angular_size=self.source1.angular_size, + e1=self.source1.ellipticity[0], + e2=self.source1.ellipticity[1], + n_sersic=0.8, + cosmo=self.source1._cosmo, + catalog_path=hst_cosmos_path, + catalog_type="HST_COSMOS", + band_dependent_color_gradient=True, + color_gradient={ + "grad_color": -0.3, + "reference_band": "F814W", + }, + **source_dict, + ) + _, reference_kwargs = source.kwargs_extended_light(band="i") + reference_image = source._image_for_band(band=None) + + assert not np.allclose(reference_kwargs[0]["image"], reference_image) + np.testing.assert_allclose( + np.sum(reference_kwargs[0]["image"]), np.sum(reference_image) + ) + np.testing.assert_allclose(reference_image, source._image_list[0]) + + source._color_gradient["grad_color"] = 0.0 + np.testing.assert_allclose( + source._image_for_band(band="i"), source._image_for_band(band=None) + ) + + def test_hst_chromatic_double_sersic_fallback(self): + source_dict = { + "z": 0.5, + "mag_g": 20.3, + "mag_i": 20.3, + "mag_y": 20.3, + "n_sersic": 0.8, + "angular_size": 1.3, + "e1": 0.09697001616620306, + "e2": 0.040998265256000574, + "center_x": 0.0, + "center_y": 0.0, + } + source = CatalogSource( + cosmo=FlatLambdaCDM(H0=70, Om0=0.3), + catalog_path=hst_cosmos_path, + catalog_type="HST_COSMOS", + max_scale=0.1, + band_dependent_color_gradient=True, + color_gradient={ + "component_spectral_slopes": [2.0, -1.0], + "reference_band": "i", + }, + **source_dict, + ) + source_model, kwargs_light = source.kwargs_extended_light(band="y") + _, kwargs_light_blue = source.kwargs_extended_light(band="g") + + assert source_model == ["SERSIC_ELLIPSE", "SERSIC_ELLIPSE"] + assert isinstance(source.double_sersic, DoubleSersic) + assert len(kwargs_light) == 2 + + flux_y = 10 ** (-np.array([item["magnitude"] for item in kwargs_light]) / 2.5) + flux_g = 10 ** ( + -np.array([item["magnitude"] for item in kwargs_light_blue]) / 2.5 + ) + assert flux_y[0] / np.sum(flux_y) > flux_g[0] / np.sum(flux_g) + + def test_chromatic_catalog_source_validation(self): + source_dict = dict(self.source1.source_dict) + common_kwargs = { + "angular_size": self.source1.angular_size, + "e1": self.source1.ellipticity[0], + "e2": self.source1.ellipticity[1], + "n_sersic": 0.8, + "cosmo": self.source1._cosmo, + "catalog_path": hst_cosmos_path, + "band_dependent_color_gradient": True, + } + + cosmos_web_kwargs = dict(common_kwargs) + cosmos_web_kwargs["catalog_path"] = cosmos_web_path + with pytest.raises(ValueError, match="received catalog_type='COSMOS_WEB'"): + CatalogSource( + catalog_type="COSMOS_WEB", + **cosmos_web_kwargs, + **source_dict, + ) + + with pytest.raises(ValueError, match=r"received None \(type NoneType\)"): + CatalogSource(catalog_type="HST_COSMOS", **common_kwargs, **source_dict) + + with pytest.raises(ValueError, match="fallback_double_sersic_kwargs"): + CatalogSource( + catalog_type="HST_COSMOS", + color_gradient={"grad_color": -0.1}, + fallback_double_sersic_kwargs="invalid", + **common_kwargs, + **source_dict, + ) + def test_redshift(self): assert self.source1.redshift == 3.5 diff --git a/tests/test_Sources/test_SourceTypes/test_double_sersic.py b/tests/test_Sources/test_SourceTypes/test_double_sersic.py index 30d0c6b46..3e1f9de1a 100644 --- a/tests/test_Sources/test_SourceTypes/test_double_sersic.py +++ b/tests/test_Sources/test_SourceTypes/test_double_sersic.py @@ -1,5 +1,7 @@ from slsim.Sources.SourceTypes.double_sersic import DoubleSersic +from slsim.Util.color_gradient import default_reference_band from slsim.Util.param_util import ellipticity_slsim_to_lenstronomy +import numpy as np import pytest from numpy import testing as npt @@ -78,10 +80,103 @@ def test_kwargs_extended_source_light(self): assert source_model[0] == "SERSIC_ELLIPSE" assert source_model[1] == "SERSIC_ELLIPSE" + def test_kwargs_extended_light_returns_independent_model_list(self): + source_model, _ = self.source.kwargs_extended_light(band="i") + source_model.append("SERSIC") + + assert self.source._light_model_list == [ + "SERSIC_ELLIPSE", + "SERSIC_ELLIPSE", + ] + def test_surface_brightness_reff(self): result = self.source.surface_brightness_reff(band="i") npt.assert_almost_equal(result, 21.313, decimal=3) + def test_band_dependent_color_gradient(self): + source_dict = dict(self.source_dict) + source_dict.update( + { + "mag_g": 23, + "mag_y": 23, + "color_gradient": { + "component_spectral_slopes": [2.0, -1.0], + "reference_band": "i", + }, + } + ) + source = DoubleSersic(**source_dict) + + _, kwargs_g = source.kwargs_extended_light(band="g") + _, kwargs_y = source.kwargs_extended_light(band="y") + + flux_g0 = 10 ** (-kwargs_g[0]["magnitude"] / 2.5) + flux_g1 = 10 ** (-kwargs_g[1]["magnitude"] / 2.5) + flux_y0 = 10 ** (-kwargs_y[0]["magnitude"] / 2.5) + flux_y1 = 10 ** (-kwargs_y[1]["magnitude"] / 2.5) + + assert flux_y0 / (flux_y0 + flux_y1) > flux_g0 / (flux_g0 + flux_g1) + + def test_color_gradient_disabled_or_missing_slopes_uses_base_weights(self): + assert self.source._weights_for_band("i") == (0.4, 0.6) + + source_dict = dict(self.source_dict) + source_dict["color_gradient"] = {} + no_gradient_source = DoubleSersic(**source_dict) + assert no_gradient_source._weights_for_band("i") == (0.4, 0.6) + + def test_default_reference_band_and_weight_validation(self): + source_dict = dict(self.source_dict) + source_dict.update( + { + "mag_g": 23, + "mag_y": 23, + "color_gradient": {"component_spectral_slopes": [1.0, 0.0]}, + } + ) + source = DoubleSersic(**source_dict) + + assert default_reference_band(source.source_dict) == "i" + assert source._weights_for_band("i") == (0.4, 0.6) + + source_without_magnitudes = dict(self.source_dict) + source_without_magnitudes.pop("mag_i") + source_without_magnitudes["color_gradient"] = { + "component_spectral_slopes": [1.0, 0.0] + } + no_magnitude_source = DoubleSersic(**source_without_magnitudes) + assert default_reference_band(no_magnitude_source.source_dict) == "i" + + source_dict["color_gradient"] = "invalid" + with pytest.raises(ValueError, match="must be a dictionary"): + DoubleSersic(**source_dict)._weights_for_band("i") + + source_dict["color_gradient"] = { + "component_spectral_slopes": [1.0], + } + with pytest.raises(ValueError, match="must match the number of components"): + DoubleSersic(**source_dict)._weights_for_band("i") + + source_dict["color_gradient"] = { + "component_spectral_slopes": [1.0, 0.0], + "min_weight": 0.5, + } + with pytest.raises(ValueError, match="must be in"): + DoubleSersic(**source_dict)._weights_for_band("i") + + def test_color_gradient_clips_component_weight(self): + source_dict = dict(self.source_dict) + source_dict["color_gradient"] = { + "component_spectral_slopes": [100.0, 0.0], + "reference_band": "i", + "min_weight": 0.2, + } + source = DoubleSersic(**source_dict) + + w0, w1 = source._weights_for_band("F213") + assert np.isclose(w0, 0.8) + assert np.isclose(w0 + w1, 1.0) + if __name__ == "__main__": pytest.main() diff --git a/tests/test_Util/test_color_gradient.py b/tests/test_Util/test_color_gradient.py new file mode 100644 index 000000000..305c5a964 --- /dev/null +++ b/tests/test_Util/test_color_gradient.py @@ -0,0 +1,164 @@ +import numpy as np +import pytest + +from slsim.Util.color_gradient import ( + attach_foreground_deflector_color_gradient, + component_weights_for_band, + default_reference_band, + radial_color_gradient_image, +) +from slsim.ImageSimulation.image_quality_lenstronomy import register_observatory +from astropy.table import Table + + +class DummyObservatory: + def __init__(self, band, **kwargs): + self.band = band + + +def test_component_weights_for_band_uses_power_law_sed_slopes(): + color_gradient = { + "component_spectral_slopes": [2.0, -1.0], + "reference_band": "i", + } + + w_g = component_weights_for_band((0.4, 0.6), "g", color_gradient) + w_i = component_weights_for_band((0.4, 0.6), "i", color_gradient) + w_y = component_weights_for_band((0.4, 0.6), "y", color_gradient) + + assert w_i == pytest.approx((0.4, 0.6)) + assert w_y[0] > w_g[0] + assert np.sum(w_y) == pytest.approx(1.0) + + +def test_attach_foreground_deflector_color_gradient_adds_table_columns(): + table = Table({"mag_i": [20.0, 21.0]}) + color_gradient = { + "component_spectral_slopes": [2.0, -1.0], + "reference_band": "i", + } + + result = attach_foreground_deflector_color_gradient( + table, color_gradient, component_weights=(2, 3) + ) + + assert result is table + assert table["color_gradient"][0] == color_gradient + assert table["w0"][0] == pytest.approx(0.4) + assert table["w1"][0] == pytest.approx(0.6) + + unchanged = Table({"mag_i": [20.0]}) + assert attach_foreground_deflector_color_gradient(unchanged, None) is unchanged + + with pytest.raises(ValueError, match="must be a dictionary"): + attach_foreground_deflector_color_gradient(table, "bad") + with pytest.raises(ValueError, match="two values"): + attach_foreground_deflector_color_gradient(table, color_gradient, (1, 2, 3)) + with pytest.raises(ValueError, match="positive sum"): + attach_foreground_deflector_color_gradient(table, color_gradient, (0, 0)) + + +def test_component_weights_validation_and_default_reference(): + source_dict = {"mag_g": 22, "mag_i": 22, "mag_y": 22} + assert default_reference_band(source_dict) == "i" + + assert component_weights_for_band((2, 3), "i") == pytest.approx((0.4, 0.6)) + assert component_weights_for_band((0.4, 0.6), "i", {}) == pytest.approx((0.4, 0.6)) + + with pytest.raises(ValueError, match="must be a dictionary"): + component_weights_for_band((0.4, 0.6), "i", "bad") + + with pytest.raises(ValueError, match="must match the number of components"): + component_weights_for_band( + (0.4, 0.6), + "i", + {"component_spectral_slopes": [1.0]}, + ) + + with pytest.raises(ValueError, match="must be in"): + component_weights_for_band( + (0.4, 0.6), + "i", + {"component_spectral_slopes": [1.0, 0.0], "min_weight": 0.5}, + ) + + register_observatory( + name="ZeroWavelengthReferenceObs", + observatory_class=DummyObservatory, + bands=["ZW1", "ZW2"], + ) + with pytest.raises(ValueError, match="neither an explicitly configured"): + component_weights_for_band( + (0.4, 0.6), + "i", + {"component_spectral_slopes": [1.0, 0.0], "reference_band": "ZW1"}, + ) + + register_observatory( + name="NonPositiveReferenceWavelengthObs", + observatory_class=DummyObservatory, + bands=["NPW0", "NPW1"], + speclite_fmt=None, + effective_wavelengths={"NPW0": 0.0, "NPW1": 1.0}, + ) + with pytest.raises(ValueError, match="reference band wavelength must be positive"): + component_weights_for_band( + (0.4, 0.6), + "NPW1", + { + "component_spectral_slopes": [1.0, 0.0], + "reference_band": "NPW0", + }, + ) + + +def test_radial_color_gradient_image_preserves_flux_and_reference_band(): + image = np.ones((9, 9)) + color_gradient = {"grad_color": -0.4, "reference_band": "F814W"} + + assert ( + radial_color_gradient_image( + image=image, + band=None, + color_gradient=color_gradient, + angular_size=0.3, + pixel_scale=0.03, + ) + is image + ) + with pytest.raises(ValueError, match="must be a dictionary"): + radial_color_gradient_image( + image=image, + band="i", + color_gradient="bad", + angular_size=0.3, + pixel_scale=0.03, + ) + + reference = radial_color_gradient_image( + image=image, + band="F814W", + color_gradient=color_gradient, + angular_size=0.3, + pixel_scale=0.03, + ) + np.testing.assert_allclose(reference, image) + + chromatic = radial_color_gradient_image( + image=image, + band="i", + color_gradient=color_gradient, + angular_size=0.3, + pixel_scale=0.03, + ) + assert not np.allclose(chromatic, image) + assert np.sum(chromatic) == pytest.approx(np.sum(image)) + + no_gradient = radial_color_gradient_image( + image=image, + band="i", + color_gradient={"grad_color": 0.0}, + angular_size=0.3, + pixel_scale=0.03, + ) + np.testing.assert_allclose(no_gradient, image)