Making a corner plot of free parameters

Making a corner plot of free parameters#

This notebook shows an example of making a corner plot of free model parameters from two COOLEST files/models. A chain file must be contained in the directory of each COOLEST file. In this case we plot the mass model and external shear parameters (they must match between the two COOLEST files, otherwise an exception is raised).

Note: this notebook is aiming to reproduce fig. 12 from Vernardos & Koopmans (2022).

authors: @gvernard, @aymgal

created on: 20/07/23

last update: 06/11/23

%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, ParametersPlotter
from coolest.api import util

def print_info(coolest_object):
    print("Lensing entities:", [type(le).__name__ for le in coolest_object.lensing_entities])

Load a few COOLEST files#

We are loading two COOLEST files, the ones corresponding to fig. 12 in Vernardos & Koopmans (2022). We add a quick check to ensure that each file has a chain file associated to it in the ‘meta’ field.

model_1_path = 'database/vkl_paper_3.3_dummy/fff_gauss_gauss_n3/VKL_coolest/coolest_vkl'
coolest_1 = util.get_coolest_object(model_1_path, verbose=True, check_external_files=False)
dir_1 = os.path.dirname(model_1_path)
print_info(coolest_1)

model_2_path = 'database/vkl_paper_3.4_dummy/fff_exp_gauss_n3/VKL_coolest/coolest_vkl'
coolest_2 = util.get_coolest_object(model_2_path, verbose=True, check_external_files=False)
dir_2 = os.path.dirname(model_2_path)
print_info(coolest_2)

truth_path = 'database/corner_truth/coolest_vkl'
truth = util.get_coolest_object(truth_path, verbose=True, check_external_files=False)
print_info(truth)

truth_path2 = 'database/corner_truth_2/coolest_vkl'
truth2 = util.get_coolest_object(truth_path2, verbose=True, check_external_files=False)
print_info(truth2)

Template file '/Users/aymgal/Science/packages/my_packages/coolest/docs/notebooks/database/vkl_paper_3.3_dummy/fff_gauss_gauss_n3/VKL_coolest/coolest_vkl_pyAPI.json' not found, now trying to read '/Users/aymgal/Science/packages/my_packages/coolest/docs/notebooks/database/vkl_paper_3.3_dummy/fff_gauss_gauss_n3/VKL_coolest/coolest_vkl.json'.
Lensing entities: ['MassField', 'Galaxy', 'Galaxy', 'Galaxy']
Template file '/Users/aymgal/Science/packages/my_packages/coolest/docs/notebooks/database/vkl_paper_3.4_dummy/fff_exp_gauss_n3/VKL_coolest/coolest_vkl_pyAPI.json' not found, now trying to read '/Users/aymgal/Science/packages/my_packages/coolest/docs/notebooks/database/vkl_paper_3.4_dummy/fff_exp_gauss_n3/VKL_coolest/coolest_vkl.json'.
Lensing entities: ['MassField', 'Galaxy', 'Galaxy', 'Galaxy']
Template file '/Users/aymgal/Science/packages/my_packages/coolest/docs/notebooks/database/corner_truth/coolest_vkl_pyAPI.json' not found, now trying to read '/Users/aymgal/Science/packages/my_packages/coolest/docs/notebooks/database/corner_truth/coolest_vkl.json'.
Lensing entities: ['MassField', 'Galaxy', 'Galaxy', 'Galaxy']
Template file '/Users/aymgal/Science/packages/my_packages/coolest/docs/notebooks/database/corner_truth_2/coolest_vkl_pyAPI.json' not found, now trying to read '/Users/aymgal/Science/packages/my_packages/coolest/docs/notebooks/database/corner_truth_2/coolest_vkl.json'.
Lensing entities: ['MassField', 'Galaxy', 'Galaxy', 'Galaxy']

Get the list of free parameters from one of the models#

tmp_free_pars = coolest_1.lensing_entities.get_parameter_ids()
free_pars = tmp_free_pars[:-2] # Remove the last parameters that refer to the light of the source and the perturbations
#print("Removed parameter(s): ",tmp_free_pars[-2:])

# Re-order parameters
reorder = [2,3,4,5,6,0,1]
pars = [free_pars[i] for i in reorder]
free_pars = pars

pprint(free_pars)

['1-galaxy-mass-0-SIE-theta_E',
 '1-galaxy-mass-0-SIE-q',
 '1-galaxy-mass-0-SIE-phi',
 '1-galaxy-mass-0-SIE-center_x',
 '1-galaxy-mass-0-SIE-center_y',
 '0-massfield-mass-0-ExternalShear-gamma_ext',
 '0-massfield-mass-0-ExternalShear-phi_ext']

Set up the plotter#

colors = ['#7FB6F5', '#E03424']

param_plotter = ParametersPlotter(
    free_pars, [coolest_1, coolest_2], 
    coolest_directories=[dir_1, dir_2], 
    coolest_names=["Smooth source", "Complex source"], 
    ref_coolest_objects=[truth, truth2],
    colors=colors,
)


# initialize the GetDist plots
settings = {
    "ignore_rows": 0.0,
    "fine_bins_2D": 800,
    "smooth_scale_2D": 0.5,
    "mult_bias_correction_order": 5
}
param_plotter.init_getdist(settings_mcsamples=settings)

Make the corner plot#

g = param_plotter.plot_triangle_getdist(filled_contours=True, subplot_size=2)

WARNING:root:outlier fraction 0.001990049751243781 
WARNING:root:outlier fraction 0.0008496176720475786 
WARNING:root:fine_bins not large enough to well sample smoothing scale - 1-galaxy-mass-0-SIE-theta_E
WARNING:root:fine_bins not large enough to well sample smoothing scale - 1-galaxy-mass-0-SIE-center_x
WARNING:root:fine_bins_2D not large enough for optimal density: 1-galaxy-mass-0-SIE-theta_E, 1-galaxy-mass-0-SIE-center_x
WARNING:root:fine_bins_2D not large enough for optimal density: 1-galaxy-mass-0-SIE-theta_E, 1-galaxy-mass-0-SIE-center_y
WARNING:root:fine_bins_2D not large enough for optimal density: 1-galaxy-mass-0-SIE-center_x, 1-galaxy-mass-0-SIE-center_y
../_images/1f1560e86d9da8aeb551f1cc9a7f4375446f288130ececdf853d1b0f1d0e75c6.png
Contents