Side-to-side model comparison plots

Side-to-side model comparison plots#

This notebook shows an example of lens model comparison side-to-side featuring very different source modeling strategies, using the Plotter classes. It also shows how to use the Analysis class to compute effective quantities in a model-independent way.

Note: this notebook was used to produce the figure in the JOSS paper.

authors: @aymgal

created on: 27/04/23

last update: 16/02/26

%config InlineBackend.figure_format = 'retina'

import os
import time
import glob
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import Normalize, LogNorm, TwoSlopeNorm
from pprint import pprint

from coolest.api.analysis import Analysis
from coolest.api.plotting import ModelPlotter, MultiModelPlotter
from coolest.api import util
from coolest.api import plot_util as plut

def print_info(coolest_object):
    source_index = 2  # index of the source galaxy in the list of `lensing entities`
    print("Lensing entities:", [type(le).__name__ for le in coolest_object.lensing_entities])
    print("Source light model:", [type(m).__name__ for m in coolest_object.lensing_entities[source_index].light_model])

Load the lens models stored in COOLEST format#

Code n°1, based on an analytical source model#

code_1_path = 'database/code_1/coolest_code_1'
coolest_1 = util.get_coolest_object(code_1_path, verbose=True)
print_info(coolest_1)

Lensing entities: ['MassField', 'Galaxy', 'Galaxy']
Source light model: ['Sersic', 'Shapelets']

/Users/aymgal/Science/packages/my_packages/coolest/coolest/template/json.py:497: UserWarning: Warning: the entity 'an external shear' is a mass field, hence it is likely a lensing mass component. Setting 'lensed' to False.
  warnings.warn(f"Warning: the entity '{entity.name}' is a mass field, hence it is likely a lensing mass component. Setting 'lensed' to False.")
/Users/aymgal/Science/packages/my_packages/coolest/coolest/template/json.py:491: UserWarning: Warning: the galaxy 'a lens galaxy' has a mass model, hence it is likely a lens galaxy. Setting 'lensed' to False.
  warnings.warn(f"Warning: the galaxy '{entity.name}' has a mass model, hence it is likely a lens galaxy. Setting 'lensed' to False.")
/Users/aymgal/Science/packages/my_packages/coolest/coolest/template/json.py:494: UserWarning: Warning: the galaxy 'an analytical source galaxy' does not have a mass model, hence it is likely a source galaxy. Setting 'lensed' to True.
  warnings.warn(f"Warning: the galaxy '{entity.name}' does not have a mass model, hence it is likely a source galaxy. Setting 'lensed' to True.")

Code n°2, based on pixelated source model on a regular grid#

code_2_path = 'database/code_2/coolest_code_2'
coolest_2 = util.get_coolest_object(code_2_path, verbose=True)
print_info(coolest_2)

Template file '/Users/aymgal/Science/packages/my_packages/coolest/docs/notebooks/database/code_2/coolest_code_2_pyAPI.json' not found, now trying to read '/Users/aymgal/Science/packages/my_packages/coolest/docs/notebooks/database/code_2/coolest_code_2.json'.
Lensing entities: ['MassField', 'Galaxy', 'Galaxy']
Source light model: ['PixelatedRegularGrid']

/Users/aymgal/Science/packages/my_packages/coolest/coolest/template/json.py:480: UserWarning: Warning: the entity 'extshear' is a mass field, hence it is likely a lensing mass component. Setting 'lensed' to False.
  warnings.warn(f"Warning: the entity '{entity['name']}' is a mass field, hence it is likely a lensing mass component. Setting 'lensed' to False.")
/Users/aymgal/Science/packages/my_packages/coolest/coolest/template/json.py:474: UserWarning: Warning: the galaxy 'lens' has a mass model, hence it is likely a lens galaxy. Setting 'lensed' to False.
  warnings.warn(f"Warning: the galaxy '{entity['name']}' has a mass model, hence it is likely a lens galaxy. Setting 'lensed' to False.")
/Users/aymgal/Science/packages/my_packages/coolest/coolest/template/json.py:477: UserWarning: Warning: the galaxy 'source' does not have a mass model, hence it is likely a source galaxy. Setting 'lensed' to True.
  warnings.warn(f"Warning: the galaxy '{entity['name']}' does not have a mass model, hence it is likely a source galaxy. Setting 'lensed' to True.")

Code n°3, based on pixelated source model on an irregular grid#

code_3_path = 'database/code_3/coolest_code_3'
coolest_3 = util.get_coolest_object(code_3_path, verbose=True)
print_info(coolest_3)

Template file '/Users/aymgal/Science/packages/my_packages/coolest/docs/notebooks/database/code_3/coolest_code_3_pyAPI.json' not found, now trying to read '/Users/aymgal/Science/packages/my_packages/coolest/docs/notebooks/database/code_3/coolest_code_3.json'.
Lensing entities: ['MassField', 'Galaxy', 'Galaxy']
Source light model: ['IrregularGrid']

Setup a custom coordinates#

so that we can evaluate light profiles consistently on that specific coordinates grid.

coord_orig = util.get_coordinates(coolest_1)
x_orig, y_orig = coord_orig.pixel_coordinates
print(coord_orig.plt_extent)

coord_src = coord_orig.create_new_coordinates(pixel_scale_factor=0.1, grid_shape=(1.42, 1.42))
x_src, y_src = coord_src.pixel_coordinates
print(coord_src.plt_extent)

[np.float64(-3.4000000000000004), np.float64(3.3999999999999995), np.float64(-3.4000000000000004), np.float64(3.3999999999999995)]
[np.float64(-0.7100000000000005), np.float64(0.7099999999999999), np.float64(-0.7100000000000006), np.float64(0.7099999999999997)]

Compute effective quantities#

analysis_1 = Analysis(coolest_1, code_1_path, supersampling=5)
analysis_2 = Analysis(coolest_2, code_2_path, supersampling=5)
analysis_3 = Analysis(coolest_3, code_3_path, supersampling=5)

Effective radius of the source surface brigtness#

r_eff_src_1 = analysis_1.effective_radius_light(center=(0, 0), coordinates=coord_src, 
                                                outer_radius=1., entity_selection=[2])
print(r_eff_src_1)

r_eff_src_2 = analysis_2.effective_radius_light(center=(0, 0), coordinates=coord_src, 
                                                outer_radius=1., entity_selection=[2])
print(r_eff_src_2)

r_eff_src_3 = analysis_3.effective_radius_light(center=(0, 0), coordinates=coord_src, 
                                                outer_radius=1., entity_selection=[2])
print(r_eff_src_3)

WARNING:root:Outer limit of integration exceeds FoV; effective radius may not be accurate.
WARNING:root:Outer limit of integration exceeds FoV; effective radius may not be accurate.
WARNING:root:Outer limit of integration exceeds FoV; effective radius may not be accurate.

4076171874999994
36914062499999944
3695312499999994

Effective Einstein radius#

theta_E_1 = analysis_1.effective_einstein_radius(entity_selection=[0, 1])
print(theta_E_1)

theta_E_2 = analysis_2.effective_einstein_radius(entity_selection=[0, 1])
print(theta_E_2)

theta_E_3 = analysis_3.effective_einstein_radius(entity_selection=[0, 1])
print(theta_E_3)

INFO:root:Picked center from profile 'PEMD'
INFO:root:Picked center from profile 'PEMD'
INFO:root:Picked center from profile 'PEMD'

5625000000000007
5750000000000006
5625000000000007

Generate the comparison figure#

Instantiate the plotter engines#

#x_src, y_src = np.meshgrid(np.linspace(-1, 1, 20), np.linspace(-1, 1, 20))

# initialize the Plotters
mplotter = MultiModelPlotter(
    [coolest_1, coolest_2, coolest_3],
    coolest_directories=[
        os.path.dirname(code_1_path), 
        os.path.dirname(code_2_path), 
        os.path.dirname(code_3_path)
    ]
)
splotter_2 = ModelPlotter(coolest_2, coolest_directory=os.path.dirname(code_2_path))
splotter_3 = ModelPlotter(coolest_3, coolest_directory=os.path.dirname(code_3_path))

Select what to plot on which axes, and show final plot#

# Normalizing the scale across the data and the three models
plotters = mplotter.plotter_list
plotters.insert(0, mplotter.plotter_list[0])
vmin, vmax = plut.normalize_across_images(
    plotters, [0, 1, 1, 1], 
    auto_selection=True,
    supersampling=5, 
    convolved=True, 
    super_convolution=True
)

norm = Normalize(vmin, vmax) # LogNorm(2e-3, 5e-2)

# create the figure
fig, axes = plt.subplots(2, 4, figsize=(14, 5.5))

# multi-plotter panels
mplotter.plot_data_image(
    [axes[0, 0], None, None],
    titles=["Observed data", None, None],
    norm=norm
)
mplotter.plot_model_image(
    [axes[0, 1], axes[0, 2], axes[0, 3]],
    titles=["Image model 1", "Image model 2", "Image model 3"],
    supersampling=5, convolved=True,
    kwargs_source=dict(entity_selection=[[2], [2], [2]]),
    kwargs_lens_mass=dict(entity_selection=[[0, 1], [0, 1], [0, 1]]),
    norm=norm
)
axes[0, 1].text(0.05, 0.05, r'$\theta_{\rm E}$ = '+f'{theta_E_1:.2f}"', color='white', fontsize=12, alpha=0.8, 
                va='bottom', ha='left', transform=axes[0, 1].transAxes)
axes[0, 2].text(0.05, 0.05, r'$\theta_{\rm E}$ = '+f'{theta_E_2:.2f}"', color='white', fontsize=12, alpha=0.8, 
                va='bottom', ha='left', transform=axes[0, 2].transAxes)
axes[0, 3].text(0.05, 0.05, r'$\theta_{\rm E}$ = '+f'{theta_E_3:.2f}"', color='white', fontsize=12, alpha=0.8, 
                va='bottom', ha='left', transform=axes[0, 3].transAxes)
mplotter.plot_model_residuals(
    [axes[1, 0], None, None],
    titles=["Normalized residuals", None, None],
    supersampling=5, add_chi2_label=True, chi2_fontsize=12,
    kwargs_source=dict(entity_selection=[[2], [2], [2]]),
    kwargs_lens_mass=dict(entity_selection=[[0, 1], [0, 1], [0, 1]]),
)
res = mplotter.plot_surface_brightness(
    [axes[1, 1], None, None], 
    titles=[
        "Source model 1", 
        # "Model 2\nUnlensed source (interpolated)", 
        None,
        None, #"Model 3: source model", 
    ],
    kwargs_light=dict(entity_selection=[[2], [2], [2]]),
    norm=norm,
    neg_values_as_bad=False,
    coordinates=coord_src,
)
axes[1, 1].text(0.05, 0.05, r'$\theta_{\rm eff}$ = '+f'{r_eff_src_1:.2f}"', color='white', fontsize=12, alpha=0.8, 
                va='bottom', ha='left', transform=axes[1, 1].transAxes)

# single-plotter panels
splotter_2.plot_surface_brightness(
    axes[1, 2], 
    title="Source model 2",
    kwargs_light=dict(entity_selection=[2]),
    norm=norm,
    neg_values_as_bad=False,
)
axes[1, 2].text(0.05, 0.05, r'$\theta_{\rm eff}$ = '+f'{r_eff_src_2:.2f}"', color='white', fontsize=12, alpha=0.8, 
                va='bottom', ha='left', transform=axes[1, 2].transAxes)
splotter_3.plot_surface_brightness(
    axes[1, 3], 
    title="Source model 3",
    kwargs_light=dict(entity_selection=[2]),
    norm=norm,
    neg_values_as_bad=False,
)
axes[1, 3].text(0.05, 0.05, r'$\theta_{\rm eff}$ = '+f'{r_eff_src_3:.2f}"', color='white', fontsize=12, alpha=0.8, 
                va='bottom', ha='left', transform=axes[1, 3].transAxes)


axes[1, 0].set_xlabel(r"$x$ (arcsec)")
axes[1, 0].set_ylabel(r"$y$ (arcsec)")

fig.tight_layout()
# fig.subplots_adjust(
#     wspace=0.02,
#     #hspace=0.15,
# )
plt.show()

Auto-selected the following entities and profiles:
- source: {'entity_selection': [2], 'profile_selection': 'all'}
- lens mass: {'entity_selection': [0, 1], 'profile_selection': 'all'}
- lens light: {'entity_selection': [], 'profile_selection': 'all'}
Auto-selected the following entities and profiles:
- source: {'entity_selection': [2], 'profile_selection': 'all'}
- lens mass: {'entity_selection': [0, 1], 'profile_selection': 'all'}
- lens light: {'entity_selection': [], 'profile_selection': 'all'}

WARNING:root:Found NaN values in image prior to convolution; they have been replaced by zeros.

Auto-selected the following entities and profiles:
- source: {'entity_selection': [2], 'profile_selection': 'all'}
- lens mass: {'entity_selection': [0, 1], 'profile_selection': 'all'}
- lens light: {'entity_selection': [], 'profile_selection': 'all'}
Auto-selected the following entities and profiles:
- source: {'entity_selection': [2], 'profile_selection': 'all'}
- lens mass: {'entity_selection': [0, 1], 'profile_selection': 'all'}
- lens light: {'entity_selection': [], 'profile_selection': 'all'}

../_images/3190db8bc723f3310d372bf693fbf739695b31ea4e550a5d59c7b977628ee225.png

# uncomment to save figure for the JOSS paper
# fig.savefig("../../joss/coolest_plot_example.png", dpi=150, bbox_inches='tight')