Figure 8: Residuals and Contrast Curve#
1. Imports#
[1]:
import json
import warnings
from pathlib import Path
from copy import deepcopy
import pandas as pd
import numpy as np
from scipy import interpolate
import matplotlib.pyplot as plt
from matplotlib.lines import Line2D
import matplotlib.gridspec as gridspec
from matplotlib.legend_handler import HandlerTuple
from mpl_toolkits.axes_grid1 import make_axes_locatable
import cmocean
color_map = cmocean.cm.ice
import seaborn as sns
from applefy.detections.contrast import Contrast
from applefy.detections.uncertainty import compute_detection_uncertainty
from applefy.utils.file_handling import load_adi_data
from applefy.utils import flux_ratio2mag, mag2flux_ratio
from applefy.utils.contrast_grid import compute_contrast_from_grid
from applefy.utils.photometry import AperturePhotometryMode
from applefy.statistics import TTest, fpf_2_gaussian_sigma, gaussian_sigma_2_fpf
from fours.detection_limits.applefy_wrapper import cADIDataReductionGPU, \
PCADataReductionGPU, FourSDataReduction
from fours.utils.data_handling import read_fours_root_dir
from fours.utils.setups import contrast_grid_setup_1
2. Load the dataset#
[2]:
dataset_name = "HD22049_303_199_C-0065_C_"
root_dir = Path(read_fours_root_dir())
json_file = root_dir / Path("30_data/" + dataset_name + ".json")
Data in the FOURS_ROOT_DIR found. Location: /fast/mbonse/s4
[3]:
experiment_root_dir = root_dir / Path("70_results/x1_fake_planet_experiments/" + dataset_name)
experiment_root_dir.mkdir(exist_ok=True)
[4]:
with open(json_file) as f:
parameter_config = json.load(f)
dit_psf_template = float(parameter_config["dit_psf"])
dit_science = float(parameter_config["dit_science"])
fwhm = float(parameter_config["fwhm"])
scaling_factor = float(parameter_config["nd_scaling"])
lambda_reg = float(parameter_config["lambda_reg"])
svd_approx = int(parameter_config["svd_approx"])
pixel_scale=0.02718
[5]:
dataset_file = root_dir / Path("30_data/" + dataset_name + ".hdf5")
[6]:
science_data, angles, raw_psf_template_data = load_adi_data(
dataset_file,
data_tag="object_stacked_05",
psf_template_tag="psf_template",
para_tag="header_object_stacked_05/PARANG")
psf_template = np.median(raw_psf_template_data, axis=0)
[7]:
# we want the image to show the innermost 1.2 arcsec
print(1.2 / pixel_scale * 2)
88.30022075055187
[8]:
# we cut the image to 91 x 91 pixel to be slightly larger than 1.2 arcsec
cut_off = int((science_data.shape[1] - 91) / 2)
science_data = science_data[:, cut_off:-cut_off, cut_off:-cut_off]
3. Create Contrast Instance#
[9]:
contrast_instance = Contrast(
science_sequence=science_data,
psf_template=psf_template,
parang_rad=angles,
psf_fwhm_radius=fwhm / 2,
dit_psf_template=dit_psf_template,
dit_science=dit_science,
scaling_factor=scaling_factor,
checkpoint_dir=experiment_root_dir)
[10]:
# get fake planet setup
flux_ratios, separations, num_fake_planets = contrast_grid_setup_1(fwhm)
[11]:
contrast_instance.design_fake_planet_experiments(
flux_ratios=flux_ratios,
num_planets=num_fake_planets,
separations=separations,
overwrite=True)
Overwriting existing config files.
4. Compute cADI limits (raw contrast)#
[12]:
cadi_algorithm_function = cADIDataReductionGPU("cpu")
[13]:
contrast_instance.run_fake_planet_experiments(
algorithm_function=cadi_algorithm_function,
num_parallel=1)
Running fake planet experiments...
100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 946/946 [00:06<00:00, 139.11it/s]
[DONE]
[14]:
# for cADI the more stable mode "F" for the planet signal seems to cause problems
contrast_instance.prepare_contrast_results(
photometry_mode_planet=AperturePhotometryMode("P"),
photometry_mode_noise=AperturePhotometryMode("P"))
[15]:
contrast_grids = contrast_instance.compute_contrast_grids(
statistical_test=TTest(),
num_cores=32,
num_rot_iter=10,
compute_snr_grid=True,
safety_margin=2.5,
pixel_scale=pixel_scale)
Computing contrast grid for cADI
Computing contrast grid with multiprocessing:
...........................................................................................................................................................................................................................................................................................................................[DONE]
[16]:
tmp_contrast_grid = contrast_grids["cADI"]
tmp_contrast_grid_fpf = tmp_contrast_grid.applymap(gaussian_sigma_2_fpf)
cADI_contrast_curve = compute_contrast_from_grid(
tmp_contrast_grid_fpf.fillna(1e-2),
gaussian_sigma_2_fpf(5))
5. Restore PCA and 4S results#
5.1 Restore PCA results#
[17]:
pca_numbers = np.concatenate(
[np.arange(0, 30, 2)[1:],
np.arange(30, 100, 10),
np.arange(100, 200, 20),
np.arange(200, 550, 50)])
pca_algorithm_function = PCADataReductionGPU(
approx_svd=8000,
pca_numbers=pca_numbers,
device="cpu",
work_dir=None,
special_name="stacked_05",
verbose=False)
[18]:
contrast_instance.run_fake_planet_experiments(
algorithm_function=pca_algorithm_function,
num_parallel=1)
Running fake planet experiments...
100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 946/946 [03:25<00:00, 4.60it/s]
[DONE]
5.2 Restore 4S results#
[19]:
for lambda_reg, special_name in [
(100, "lambda_000100"),
(1000, "lambda_001000"),
(10000, "lambda_010000"),
(100000, "lambda_100000")]:
fours_model = FourSDataReduction(
device="cpu",
lambda_reg=lambda_reg,
psf_fwhm=fwhm,
right_reason_mask_factor = 1.5,
special_name=special_name,
rotation_grid_down_sample=1,
logging_interval=50,
save_models=True,
train_num_epochs=500,
work_dir=None,
verbose=False)
# We have to keep track of the 4s results to not loose them
old_results = contrast_instance.results_dict
contrast_instance.run_fake_planet_experiments(
algorithm_function=fours_model,
num_parallel=1)
contrast_instance.results_dict.update(old_results)
Running fake planet experiments...
100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 946/946 [00:13<00:00, 68.63it/s]
[DONE]
Running fake planet experiments...
100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 946/946 [00:13<00:00, 69.50it/s]
[DONE]
Running fake planet experiments...
100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 946/946 [00:13<00:00, 70.22it/s]
[DONE]
Running fake planet experiments...
100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 946/946 [00:13<00:00, 72.68it/s]
[DONE]
4. Compute the contrast grids#
[20]:
# Use pixel values spaced by the FWHM
photometry_mode_planet = AperturePhotometryMode("F")
photometry_mode_noise = AperturePhotometryMode("P")
[21]:
contrast_instance.prepare_contrast_results(
photometry_mode_planet=photometry_mode_planet,
photometry_mode_noise=photometry_mode_noise)
[22]:
with warnings.catch_warnings():
warnings.simplefilter("ignore")
contrast_grids = contrast_instance.compute_contrast_grids(
statistical_test=TTest(),
num_cores=32,
num_rot_iter=10,
compute_snr_grid=True,
safety_margin=2.5,
pixel_scale=pixel_scale)
Computing contrast grid for s4_mean_lambda_100000
Computing contrast grid with multiprocessing:
...........................................................................................................................................................................................................................................................................................................................[DONE]
Computing contrast grid for s4_median_lambda_100000
Computing contrast grid with multiprocessing:
...........................................................................................................................................................................................................................................................................................................................[DONE]
Computing contrast grid for s4_mean_lambda_010000
Computing contrast grid with multiprocessing:
...........................................................................................................................................................................................................................................................................................................................[DONE]
Computing contrast grid for s4_median_lambda_010000
Computing contrast grid with multiprocessing:
...........................................................................................................................................................................................................................................................................................................................[DONE]
Computing contrast grid for s4_mean_lambda_001000
Computing contrast grid with multiprocessing:
...........................................................................................................................................................................................................................................................................................................................[DONE]
Computing contrast grid for s4_median_lambda_001000
Computing contrast grid with multiprocessing:
...........................................................................................................................................................................................................................................................................................................................[DONE]
Computing contrast grid for s4_mean_lambda_000100
Computing contrast grid with multiprocessing:
...........................................................................................................................................................................................................................................................................................................................[DONE]
Computing contrast grid for s4_median_lambda_000100
Computing contrast grid with multiprocessing:
...........................................................................................................................................................................................................................................................................................................................[DONE]
Computing contrast grid for stacked_05_PCA_002_components
Computing contrast grid with multiprocessing:
...........................................................................................................................................................................................................................................................................................................................[DONE]
Computing contrast grid for stacked_05_PCA_004_components
Computing contrast grid with multiprocessing:
...........................................................................................................................................................................................................................................................................................................................[DONE]
Computing contrast grid for stacked_05_PCA_006_components
Computing contrast grid with multiprocessing:
...........................................................................................................................................................................................................................................................................................................................[DONE]
Computing contrast grid for stacked_05_PCA_008_components
Computing contrast grid with multiprocessing:
...........................................................................................................................................................................................................................................................................................................................[DONE]
Computing contrast grid for stacked_05_PCA_010_components
Computing contrast grid with multiprocessing:
...........................................................................................................................................................................................................................................................................................................................[DONE]
Computing contrast grid for stacked_05_PCA_012_components
Computing contrast grid with multiprocessing:
...........................................................................................................................................................................................................................................................................................................................[DONE]
Computing contrast grid for stacked_05_PCA_014_components
Computing contrast grid with multiprocessing:
...........................................................................................................................................................................................................................................................................................................................[DONE]
Computing contrast grid for stacked_05_PCA_016_components
Computing contrast grid with multiprocessing:
...........................................................................................................................................................................................................................................................................................................................[DONE]
Computing contrast grid for stacked_05_PCA_018_components
Computing contrast grid with multiprocessing:
...........................................................................................................................................................................................................................................................................................................................[DONE]
Computing contrast grid for stacked_05_PCA_020_components
Computing contrast grid with multiprocessing:
...........................................................................................................................................................................................................................................................................................................................[DONE]
Computing contrast grid for stacked_05_PCA_022_components
Computing contrast grid with multiprocessing:
...........................................................................................................................................................................................................................................................................................................................[DONE]
Computing contrast grid for stacked_05_PCA_024_components
Computing contrast grid with multiprocessing:
...........................................................................................................................................................................................................................................................................................................................[DONE]
Computing contrast grid for stacked_05_PCA_026_components
Computing contrast grid with multiprocessing:
...........................................................................................................................................................................................................................................................................................................................[DONE]
Computing contrast grid for stacked_05_PCA_028_components
Computing contrast grid with multiprocessing:
...........................................................................................................................................................................................................................................................................................................................[DONE]
Computing contrast grid for stacked_05_PCA_030_components
Computing contrast grid with multiprocessing:
...........................................................................................................................................................................................................................................................................................................................[DONE]
Computing contrast grid for stacked_05_PCA_040_components
Computing contrast grid with multiprocessing:
...........................................................................................................................................................................................................................................................................................................................[DONE]
Computing contrast grid for stacked_05_PCA_050_components
Computing contrast grid with multiprocessing:
...........................................................................................................................................................................................................................................................................................................................[DONE]
Computing contrast grid for stacked_05_PCA_060_components
Computing contrast grid with multiprocessing:
...........................................................................................................................................................................................................................................................................................................................[DONE]
Computing contrast grid for stacked_05_PCA_070_components
Computing contrast grid with multiprocessing:
...........................................................................................................................................................................................................................................................................................................................[DONE]
Computing contrast grid for stacked_05_PCA_080_components
Computing contrast grid with multiprocessing:
...........................................................................................................................................................................................................................................................................................................................[DONE]
Computing contrast grid for stacked_05_PCA_090_components
Computing contrast grid with multiprocessing:
...........................................................................................................................................................................................................................................................................................................................[DONE]
Computing contrast grid for stacked_05_PCA_100_components
Computing contrast grid with multiprocessing:
...........................................................................................................................................................................................................................................................................................................................[DONE]
Computing contrast grid for stacked_05_PCA_120_components
Computing contrast grid with multiprocessing:
...........................................................................................................................................................................................................................................................................................................................[DONE]
Computing contrast grid for stacked_05_PCA_140_components
Computing contrast grid with multiprocessing:
...........................................................................................................................................................................................................................................................................................................................[DONE]
Computing contrast grid for stacked_05_PCA_160_components
Computing contrast grid with multiprocessing:
...........................................................................................................................................................................................................................................................................................................................[DONE]
Computing contrast grid for stacked_05_PCA_180_components
Computing contrast grid with multiprocessing:
...........................................................................................................................................................................................................................................................................................................................[DONE]
Computing contrast grid for stacked_05_PCA_200_components
Computing contrast grid with multiprocessing:
...........................................................................................................................................................................................................................................................................................................................[DONE]
Computing contrast grid for stacked_05_PCA_250_components
Computing contrast grid with multiprocessing:
...........................................................................................................................................................................................................................................................................................................................[DONE]
Computing contrast grid for stacked_05_PCA_300_components
Computing contrast grid with multiprocessing:
...........................................................................................................................................................................................................................................................................................................................[DONE]
Computing contrast grid for stacked_05_PCA_350_components
Computing contrast grid with multiprocessing:
...........................................................................................................................................................................................................................................................................................................................[DONE]
Computing contrast grid for stacked_05_PCA_400_components
Computing contrast grid with multiprocessing:
...........................................................................................................................................................................................................................................................................................................................[DONE]
Computing contrast grid for stacked_05_PCA_450_components
Computing contrast grid with multiprocessing:
...........................................................................................................................................................................................................................................................................................................................[DONE]
Computing contrast grid for stacked_05_PCA_500_components
Computing contrast grid with multiprocessing:
...........................................................................................................................................................................................................................................................................................................................[DONE]
6. Plot an example of a Contrast Grid#
[23]:
def plot_contrast_grid(
contrast_grid_axis,
colorbar_axis,
contrast_grid):
c_bar_kargs = dict(
orientation = "vertical",
label = r"SNR [$T_{obs}$]")
heat = sns.heatmap(
contrast_grid,
vmax=2, vmin=7,
annot=True,
cmap="YlGnBu",
ax=contrast_grid_axis,
cbar_ax=colorbar_axis,
cbar_kws=c_bar_kargs)
ylabels = ['{:.1f}'.format(float(x.get_text()))
for x in heat.get_yticklabels()]
_=heat.set_yticklabels(ylabels)
xlabels = ['{:.1f}'.format(float(x.get_text()))
for x in heat.get_xticklabels()]
_=heat.set_xticklabels(xlabels)
[24]:
tmp_contrast_grid = deepcopy(contrast_grids['stacked_05_PCA_010_components'])
tmp_contrast_grid.index = flux_ratio2mag(tmp_contrast_grid.index)
[25]:
fig = plt.figure(figsize=(8, 4))
gs0 = fig.add_gridspec(1, 1)
gs0.update(wspace=0.0, hspace=0.2)
gs1 = gridspec.GridSpecFromSubplotSpec(
1, 2, subplot_spec = gs0[0],
wspace=0.05, width_ratios=[1, 0.03])
# All axis we need
contrast_ax = fig.add_subplot(gs1[0])
colorbar_ax = fig.add_subplot(gs1[1])
# Plot the contrast grid
plot_contrast_grid(
contrast_grid_axis=contrast_ax,
colorbar_axis=colorbar_ax,
contrast_grid=tmp_contrast_grid)
colorbar_ax.yaxis.label.set_size(14)
contrast_ax.set_ylabel(
"Contrast - $c = f_p / f_*$ - [mag]", size=14)
contrast_ax.set_xlabel(
r"Separation [$\lambda /D$]", size=14)
contrast_ax.set_title(
"Contrast Grid " + dataset_name,
fontsize=16,
fontweight="bold",
y=1.03)
contrast_ax.tick_params(
axis='both',
which='major',
labelsize=12)
# Save the figure
fig.patch.set_facecolor('white')
7. Get the S/N = 5 limits from contrast grids#
7.1 Threshold the PCA result#
[26]:
pca_results = dict()
with warnings.catch_warnings():
warnings.simplefilter("ignore")
for tmp_pca_number in pca_numbers:
tmp_contrast_grid = contrast_grids["stacked_05_PCA_" + str(tmp_pca_number).zfill(3) + "_components"]
tmp_contrast_grid_fpf = tmp_contrast_grid.applymap(gaussian_sigma_2_fpf)
tmp_contrast_curve = compute_contrast_from_grid(
tmp_contrast_grid_fpf.fillna(1e-2),
gaussian_sigma_2_fpf(5))
pca_results[tmp_pca_number] = tmp_contrast_curve.applymap(flux_ratio2mag)
[27]:
pca_results = pd.concat(pca_results.values(), keys=pca_results.keys(), axis=1)
pca_results.index = np.round(pca_results.index, 2)
pca_results.columns = pca_results.columns.get_level_values(0)
7.2 Get the optimal PCA#
[28]:
overall_best_pca = pca_results.max(axis=1)
7.3 Threshold the 4S result#
[29]:
s4_lambdas = [100, 1000, 10000, 100000]
[30]:
s4_results = dict()
with warnings.catch_warnings():
warnings.simplefilter("ignore")
for tmp_lambda in s4_lambdas:
for tmp_combine in ["mean", "median"]:
tmp_contrast_grid = contrast_grids["s4_" + tmp_combine + "_lambda_" + str(tmp_lambda).zfill(6)]
tmp_contrast_grid_fpf = tmp_contrast_grid.applymap(gaussian_sigma_2_fpf)
tmp_contrast_curve = compute_contrast_from_grid(
tmp_contrast_grid_fpf.fillna(1e-2),
gaussian_sigma_2_fpf(5))
s4_results["s4_" + tmp_combine + "_lambda_" + str(tmp_lambda).zfill(6)] = tmp_contrast_curve.applymap(flux_ratio2mag)
[31]:
s4_results = pd.concat(s4_results.values(), keys=s4_results.keys(), axis=1)
s4_results.index = np.round(s4_results.index, 2)
s4_results.columns = s4_results.columns.get_level_values(0)
[32]:
s4_contrast = s4_results["s4_mean_lambda_000100"]
8. Plot the S/N = 5 curves#
[33]:
separations_FWHM = s4_results.index
separations_arcsec = s4_results.index * pixel_scale * fwhm
[34]:
# Find one color for each number of PCA components used.
colors = sns.color_palette("rocket_r",
n_colors=pca_results.shape[1])
colors
[34]:
[35]:
# Examples for later
example_1 = (1.5, 7.5, 0)
example_2 = (2.5, 9.5, 2)
example_3 = (11., 13., 1)
[36]:
# 1.) Create Plot Layout
fig = plt.figure(constrained_layout=False, figsize=(8, 7))
gs0 = fig.add_gridspec(2, 1, height_ratios=[1, 0.1], hspace=0.3)
axis_contrast_curvse = fig.add_subplot(gs0[0, 0])
# ------------------- Create the PCA number legend -----------
axis_legend_lines = fig.add_subplot(gs0[1, 0])
for idx, tmp_components in enumerate(
pca_results.columns.values):
if tmp_components < 30:
lw=3
elif tmp_components < 100:
lw=12
elif tmp_components < 200:
lw=24
else:
lw=57
axis_legend_lines.vlines(
[tmp_components],
ymin=0, ymax=1, color=colors[idx],
lw=lw)
axis_legend_lines.set_xlim(0, 500)
axis_legend_lines.set_ylim(0, 1)
axis_legend_lines.set_yticks([])
axis_legend_lines.tick_params(
axis='both', which='major', labelsize=14)
axis_legend_lines.set_xlabel("Number of PCA components", size=16)
# ---------------------- Create the PCA plots --------------------
i = 0 # color picker
for pca_number, tmp_contrast_mag in pca_results.items():
axis_contrast_curvse.plot(
separations_arcsec,
mag2flux_ratio(tmp_contrast_mag),
color = colors[i],
lw=2,
alpha=0.7)
i+=1
axis_contrast_curvse.set_yscale("log")
# ----------- Plot the best PCA contrast -------------------------
pca_plot = axis_contrast_curvse.plot(
separations_arcsec,
mag2flux_ratio(overall_best_pca),
color = "black",
lw=2.5,
label="optimal PCA",
ls="dashed")
# ----------- Plot the raw contrast -------------------------
pca_plot = axis_contrast_curvse.plot(
separations_arcsec[1:],
cADI_contrast_curve.values[1:],
color = "black",
lw=2.5,
label="cADI",
ls="dotted")
# ----------- Plot the S4 contrast --------------------------
s4_plot = axis_contrast_curvse.plot(
separations_arcsec,
mag2flux_ratio(s4_contrast),
color = "darkblue",
lw=2.5,
label="${4S}_{\lambda=100}$",
ls="-")
# ------------- Double axis and limits -----------------------
lim_mag_y = (14.5, 6)
lim_arcsec_x = (0.1, 1.2)
sep_lambda_arcse = interpolate.interp1d(
separations_arcsec,
separations_FWHM,
fill_value='extrapolate')
axis_contrast_curvse_mag = axis_contrast_curvse.twinx()
axis_contrast_curvse_mag.plot(
separations_arcsec,
s4_contrast,
alpha=0.)
axis_contrast_curvse_mag.invert_yaxis()
axis_contrast_curvse_lambda = axis_contrast_curvse.twiny()
axis_contrast_curvse_lambda.plot(
separations_FWHM,
mag2flux_ratio(s4_contrast),
alpha=0.)
axis_contrast_curvse_mag.grid(
which='major', ls="dotted", lw=1,
alpha=0.2, color="black", zorder=0)
axis_contrast_curvse_mag.set_ylim(*lim_mag_y)
axis_contrast_curvse.set_ylim(
mag2flux_ratio(lim_mag_y[0]),
mag2flux_ratio(lim_mag_y[1]))
axis_contrast_curvse.set_xlim(
*lim_arcsec_x)
axis_contrast_curvse_mag.set_xlim(
*lim_arcsec_x)
axis_contrast_curvse_lambda.set_xlim(
*sep_lambda_arcse(lim_arcsec_x))
# ------------- Examples -----------------------
ex1=axis_contrast_curvse.scatter(
[example_1[0] * fwhm * pixel_scale],
[mag2flux_ratio(example_1[1])],
s=150,
color="deepskyblue",
marker="^",
zorder=10)
ex2=axis_contrast_curvse.scatter(
[example_2[0] * fwhm * pixel_scale],
[mag2flux_ratio(example_2[1])],
s=150,
color="deepskyblue",
marker="s",
zorder=10)
ex3=axis_contrast_curvse.scatter(
[example_3[0] * fwhm * pixel_scale],
[mag2flux_ratio(example_3[1])],
s=150,
color="deepskyblue",
marker="*",
zorder=10)
# ----------- Labels and fontsizes --------------------------
axis_contrast_curvse.set_xlabel(
r"Separation [arcsec]", size=16)
axis_contrast_curvse_lambda.set_xlabel(
r"Separation [$\lambda/D$]", size=16)
axis_contrast_curvse.set_ylabel(
r"Planet-to-star flux ratio", size=16)
axis_contrast_curvse_mag.set_ylabel(
r"$\Delta$ Magnitude", size=16)
axis_contrast_curvse.tick_params(
axis='both', which='major', labelsize=14)
axis_contrast_curvse_lambda.tick_params(
axis='both', which='major', labelsize=14)
axis_contrast_curvse_mag.tick_params(
axis='both', which='major', labelsize=14)
axis_contrast_curvse_mag.set_title(
r"$S/N = 5$ Contrast Limits - #5 (HD 22049)",
fontsize=18, fontweight="bold", y=1.15)
# --------------------------- Legend -----------------------
handles, labels = axis_contrast_curvse.\
get_legend_handles_labels()
labels.append("Examples")
# Create a legend with grouped markers and line plot
legend_elements = [
handles[0],
handles[1],
handles[2],
(ex1, ex2, ex3)]
# Add legend
leg1 = fig.legend(legend_elements, labels,
bbox_to_anchor=(0.14, 0.28),
fontsize=14,
handler_map={tuple: HandlerTuple(ndivide=None, pad=1.7)},
loc='lower left', ncol=2)
_=plt.setp(leg1.get_title(),fontsize=14)
fig.patch.set_facecolor('white')
plt.savefig("./final_plots/05_contrast_limits.pdf", bbox_inches='tight')
[37]:
# the contrast improvement compared to cADI! [mag]
(s4_contrast.values[1:] - flux_ratio2mag(cADI_contrast_curve.values.flatten()[1:]))
[37]:
array([3.50935351, 3.51034351, 3.61590362, 3.66838367, 3.36663337,
3.39481339, 3.99031399, 3.71503372, 3.73565374, 3.64799365,
3.37367337, 2.76830277, 2.3965824 , 2.003802 ])
9. Plot example residuals#
[38]:
best_number_of_components = pca_results.T.idxmax().values
[39]:
magnitudes = flux_ratio2mag(contrast_grids["s4_mean_lambda_000100"].index)
example_1_idx_sep = np.where(separations_FWHM == example_1[0])[0][0]
example_1_idx_mag = np.where(magnitudes == example_1[1])[0][0]
example_1_pcas = best_number_of_components[example_1_idx_sep]
example_2_idx_sep = np.where(separations_FWHM == example_2[0])[0][0]
example_2_idx_mag = np.where(magnitudes == example_2[1])[0][0]
example_2_pcas = best_number_of_components[example_2_idx_sep]
example_3_idx_sep = np.where(separations_FWHM == example_3[0])[0][0]
example_3_idx_mag = np.where(magnitudes == example_3[1])[0][0]
example_3_pcas = best_number_of_components[example_3_idx_sep]
9.1 Get the residual images#
[40]:
s4_residuals = contrast_instance.contrast_results["s4_mean_lambda_000100"].residuals
s4_residual_1 = s4_residuals[example_1_idx_sep, example_1_idx_mag, example_1[2]]
s4_residual_2 = s4_residuals[example_2_idx_sep, example_2_idx_mag, example_2[2]]
s4_residual_3 = s4_residuals[example_3_idx_sep, example_3_idx_mag, example_3[2]]
[41]:
tmp_residuals = contrast_instance.contrast_results[
"stacked_05_PCA_" + str(example_1_pcas).zfill(3) +"_components"].residuals
pca_residual_1 = tmp_residuals[example_1_idx_sep, example_1_idx_mag, example_1[2]]
tmp_residuals = contrast_instance.contrast_results[
"stacked_05_PCA_" + str(example_2_pcas).zfill(3) +"_components"].residuals
pca_residual_2 = tmp_residuals[example_2_idx_sep, example_2_idx_mag, example_2[2]]
tmp_residuals = contrast_instance.contrast_results[
"stacked_05_PCA_" + str(example_3_pcas).zfill(3) +"_components"].residuals
pca_residual_3 = tmp_residuals[example_3_idx_sep, example_3_idx_mag, example_3[2]]
9.2 compute the S/N#
[42]:
def get_snr(
method_name,
example):
contrast_result = contrast_instance.contrast_results[method_name]
contrast_result.idx_table.index = np.round(contrast_result.idx_table.index, 2)
example_idx = contrast_result.idx_table[mag2flux_ratio(example[1])][example[0] * fwhm]
residual, position = contrast_result.planet_dict[example_idx][example[2]]
with warnings.catch_warnings():
warnings.simplefilter("ignore")
_, _, snr = compute_detection_uncertainty(
frame=residual,
planet_position=position,
statistical_test=TTest(),
psf_fwhm_radius=fwhm,
photometry_mode_planet=photometry_mode_planet,
photometry_mode_noise=photometry_mode_noise,
safety_margin=2.5,
num_rot_iter=20)
snr = snr[~np.isnan(snr)]
return np.median(snr)
[43]:
pca_example_1_snr = np.round(
get_snr("stacked_05_PCA_" + str(example_1_pcas).zfill(3) +"_components",
example_1), 2)
pca_example_2_snr = np.round(
get_snr("stacked_05_PCA_" + str(example_2_pcas).zfill(3) +"_components",
example_2), 2)
pca_example_3_snr = np.round(
get_snr("stacked_05_PCA_" + str(example_3_pcas).zfill(3) +"_components",
example_3), 2)
s4_example_1_snr = np.round(get_snr("s4_mean_lambda_000100", example_1), 2)
s4_example_2_snr = np.round(get_snr("s4_mean_lambda_000100", example_2), 2)
s4_example_3_snr = np.round(get_snr("s4_mean_lambda_000100", example_3), 2)
[44]:
def add_colorbar(axis_in, plot_in, left=True):
if left:
tick_side = "left"
else:
tick_side = "right"
divider = make_axes_locatable(axis_in)
cax = divider.new_horizontal(size="5%", pad=0.1, pack_start=left)
fig.add_axes(cax)
cbar = fig.colorbar(plot_in, cax=cax,
orientation="vertical")
cbar.ax.yaxis.set_ticks_position(tick_side)
axis_in.set_xticks([])
axis_in.set_yticks([])
cbar.ax.tick_params(labelsize=10)
[45]:
def annotate_pca(axis_in, num_pca, snr):
t2 = axis_in.annotate(
xytext=(5., 9.),
xy=(0, 0),
text=str(num_pca) + " components; S/N = " + str(snr),
color='w',
va='top',
ha='left',
fontsize=10,
bbox=dict(facecolor='black',
edgecolor='none',
boxstyle='round',
alpha=0.5))
def annotate_s4(axis_in, snr):
t2 = axis_in.annotate(
xytext=(5., 9.),
xy=(0, 0),
text="S/N = " + str(snr),
color='w',
va='top',
ha='left',
fontsize=10,
bbox=dict(facecolor='black',
edgecolor='none',
boxstyle='round',
alpha=0.5))
[46]:
fig = plt.figure(constrained_layout=False, figsize=(8, 11.2))
gs0 = fig.add_gridspec(3, 2,hspace=0.05, wspace=0.05) # height_ratios=[1, 1]
axis_pca_1 = fig.add_subplot(gs0[0, 0])
axis_pca_2 = fig.add_subplot(gs0[1, 0])
axis_pca_3 = fig.add_subplot(gs0[2, 0])
axis_s4_1 = fig.add_subplot(gs0[0, 1])
axis_s4_2 = fig.add_subplot(gs0[1, 1])
axis_s4_3 = fig.add_subplot(gs0[2, 1])
# PCA plots -----------------
plot1 = axis_pca_1.imshow(pca_residual_1, cmap=color_map, origin='lower')
add_colorbar(axis_pca_1, plot1, True)
annotate_pca(axis_pca_1, example_1_pcas, pca_example_1_snr)
plot2 = axis_pca_2.imshow(pca_residual_2, cmap=color_map, origin='lower')
add_colorbar(axis_pca_2, plot2, True)
annotate_pca(axis_pca_2, example_2_pcas, pca_example_2_snr)
plot3 = axis_pca_3.imshow(pca_residual_3, cmap=color_map, origin='lower')
add_colorbar(axis_pca_3, plot3, True)
annotate_pca(axis_pca_3, example_3_pcas, pca_example_3_snr)
# S4 plots -----------------
plot4 = axis_s4_1.imshow(s4_residual_1, cmap=color_map, origin='lower',
vmin=-0.04, vmax=0.05)
add_colorbar(axis_s4_1, plot4, False)
annotate_s4(axis_s4_1, s4_example_1_snr)
plot5 = axis_s4_2.imshow(s4_residual_2, cmap=color_map, origin='lower',
vmin=-0.04, vmax=0.05)
add_colorbar(axis_s4_2, plot5, False)
annotate_s4(axis_s4_2, s4_example_2_snr)
plot6 = axis_s4_3.imshow(s4_residual_3, cmap=color_map, origin='lower',
vmin=-0.04, vmax=0.05)
add_colorbar(axis_s4_3, plot6, False)
annotate_s4(axis_s4_3, s4_example_3_snr)
# Titles -----------------
axis_pca_1.set_title("Optimal PCA", fontsize=16, fontweight="bold", y=1.02)
axis_s4_1.set_title("{4S}_{\lambda=100}$", fontsize=16, fontweight="bold", y=1.02)
axis_pca_1.set_ylabel(
r"Example: $\blacktriangle$", fontsize=14,
fontweight="bold",labelpad=85)
axis_pca_1.text(
-35, 20,
"separation: " + str(example_1[0]) + " [$\lambda/D$] \n contrast:" + str(example_1[1]) + " [mag]",
fontsize=12, rotation=90)
axis_pca_2.set_ylabel(
r"Example: $\blacksquare$", fontsize=14,
fontweight="bold",labelpad=85)
axis_pca_2.text(
-35, 20,
"separation: " + str(example_2[0]) + " [$\lambda/D$] \n contrast:" + str(example_2[1]) + " [mag]",
fontsize=12, rotation=90)
axis_pca_3.set_ylabel(
r"Example: $\bigstar$", fontsize=14,
fontweight="bold",labelpad=85)
axis_pca_3.text(
-35, 20,
"separation: " + str(example_3[0]) + " [$\lambda/D$] \n contrast:" + str(example_3[1]) + " [mag]",
fontsize=12, rotation=90)
# Save the plot
fig.patch.set_facecolor('white')
plt.savefig("./final_plots/05_residuals_PCA_S4.pdf", bbox_inches='tight')
10. Create one large plot for both#
10.1 Prepare the data#
[47]:
separations_FWHM = s4_results.index
separations_arcsec = s4_results.index * pixel_scale * fwhm
[48]:
# Find one color for each number of PCA components used.
colors = sns.color_palette("rocket_r",
n_colors=pca_results.shape[1])
colors
[48]:
[49]:
# Examples for later
example_1 = (1.5, 7.5, 0)
example_2 = (2.5, 9.5, 2)
example_3 = (11., 13., 1)
[50]:
best_number_of_components = pca_results.T.idxmax().values
[51]:
magnitudes = flux_ratio2mag(contrast_grids["s4_mean_lambda_000100"].index)
example_1_idx_sep = np.where(separations_FWHM == example_1[0])[0][0]
example_1_idx_mag = np.where(magnitudes == example_1[1])[0][0]
example_1_pcas = best_number_of_components[example_1_idx_sep]
example_2_idx_sep = np.where(separations_FWHM == example_2[0])[0][0]
example_2_idx_mag = np.where(magnitudes == example_2[1])[0][0]
example_2_pcas = best_number_of_components[example_2_idx_sep]
example_3_idx_sep = np.where(separations_FWHM == example_3[0])[0][0]
example_3_idx_mag = np.where(magnitudes == example_3[1])[0][0]
example_3_pcas = best_number_of_components[example_3_idx_sep]
[52]:
# Get the residual images
s4_residuals = contrast_instance.contrast_results["s4_mean_lambda_000100"].residuals
s4_residual_1 = s4_residuals[example_1_idx_sep, example_1_idx_mag, example_1[2]]
s4_residual_2 = s4_residuals[example_2_idx_sep, example_2_idx_mag, example_2[2]]
s4_residual_3 = s4_residuals[example_3_idx_sep, example_3_idx_mag, example_3[2]]
[53]:
tmp_residuals = contrast_instance.contrast_results[
"stacked_05_PCA_" + str(example_1_pcas).zfill(3) +"_components"].residuals
pca_residual_1 = tmp_residuals[example_1_idx_sep, example_1_idx_mag, example_1[2]]
tmp_residuals = contrast_instance.contrast_results[
"stacked_05_PCA_" + str(example_2_pcas).zfill(3) +"_components"].residuals
pca_residual_2 = tmp_residuals[example_2_idx_sep, example_2_idx_mag, example_2[2]]
tmp_residuals = contrast_instance.contrast_results[
"stacked_05_PCA_" + str(example_3_pcas).zfill(3) +"_components"].residuals
pca_residual_3 = tmp_residuals[example_3_idx_sep, example_3_idx_mag, example_3[2]]
10.2 compute the S/N#
[54]:
def get_snr(
method_name,
example):
contrast_result = contrast_instance.contrast_results[method_name]
contrast_result.idx_table.index = np.round(contrast_result.idx_table.index, 2)
example_idx = contrast_result.idx_table[mag2flux_ratio(example[1])][example[0] * fwhm]
residual, position = contrast_result.planet_dict[example_idx][example[2]]
with warnings.catch_warnings():
warnings.simplefilter("ignore")
_, _, snr = compute_detection_uncertainty(
frame=residual,
planet_position=position,
statistical_test=TTest(),
psf_fwhm_radius=fwhm,
photometry_mode_planet=photometry_mode_planet,
photometry_mode_noise=photometry_mode_noise,
safety_margin=2.5,
num_rot_iter=20)
snr = snr[~np.isnan(snr)]
return np.median(snr)
[55]:
pca_example_1_snr = np.round(
get_snr("stacked_05_PCA_" + str(example_1_pcas).zfill(3) +"_components",
example_1), 2)
pca_example_2_snr = np.round(
get_snr("stacked_05_PCA_" + str(example_2_pcas).zfill(3) +"_components",
example_2), 2)
pca_example_3_snr = np.round(
get_snr("stacked_05_PCA_" + str(example_3_pcas).zfill(3) +"_components",
example_3), 2)
s4_example_1_snr = np.round(get_snr("s4_mean_lambda_000100", example_1), 2)
s4_example_2_snr = np.round(get_snr("s4_mean_lambda_000100", example_2), 2)
s4_example_3_snr = np.round(get_snr("s4_mean_lambda_000100", example_3), 2)
10.3 define plot helper functions#
[56]:
def remove_image_ticks(axis_in):
axis_in.set_xticks([])
axis_in.set_yticks([])
[57]:
def annotate_pca(axis_in, num_pca, snr):
#return
t2 = axis_in.annotate(
xytext=(5., 9.), #(5., 16.),
xy=(0, 0),
#text="S/N = " + str(snr) + "\n" + str(num_pca) + " components",
text="S/N = " + str(snr),
color='w',
va='top',
ha='left',
fontsize=8,
bbox=dict(facecolor='black',
edgecolor='none',
boxstyle='round',
alpha=0.5))
def annotate_s4(axis_in, snr):
#return
t2 = axis_in.annotate(
xytext=(5., 9.),
xy=(0, 0),
text="S/N = " + str(snr),
color='w',
va='top',
ha='left',
fontsize=8,
bbox=dict(facecolor='black',
edgecolor='none',
boxstyle='round',
alpha=0.5))
10.4 Create the plot#
[ ]:
# 1.) Create Plot Layout
fig = plt.figure(constrained_layout=False, figsize=(12, 5))
gs0 = fig.add_gridspec(
1, 2, width_ratios=[2, 1.2],
hspace=0.05,
wspace=0.18
)
# Contrast curve and color bar
gs1 = gridspec.GridSpecFromSubplotSpec(
4, 1, subplot_spec = gs0[1],
height_ratios=[0.05, 1, 0.1, 0.25], hspace=0.7)
axis_contrast_curvse = fig.add_subplot(gs1[1, 0])
axis_legend_lines = fig.add_subplot(gs1[2, 0])
# Example residuals
gs2 = gridspec.GridSpecFromSubplotSpec(
2, 3, subplot_spec = gs0[0],
hspace=-0.12, wspace=0.05)
axis_pca_1 = fig.add_subplot(gs2[0, 0])
axis_pca_2 = fig.add_subplot(gs2[0, 1])
axis_pca_3 = fig.add_subplot(gs2[0, 2])
axis_s4_1 = fig.add_subplot(gs2[1, 0])
axis_s4_2 = fig.add_subplot(gs2[1, 1])
axis_s4_3 = fig.add_subplot(gs2[1, 2])
# ------------------- Create the PCA number legend -----------
for idx, tmp_components in enumerate(
pca_results.columns.values):
if tmp_components < 30:
lw=3
elif tmp_components < 100:
lw=12
elif tmp_components < 200:
lw=24
else:
lw=57
axis_legend_lines.vlines(
[tmp_components],
ymin=0, ymax=1, color=colors[idx],
lw=lw)
axis_legend_lines.set_xlim(0, 500)
axis_legend_lines.set_ylim(0, 1)
axis_legend_lines.set_yticks([])
axis_legend_lines.tick_params(
axis='both', which='major', labelsize=8)
axis_legend_lines.set_xlabel("Number of PCA components", size=10)
# ---------------------- Create the PCA plots --------------------
i = 0 # color picker
for pca_number, tmp_contrast_mag in pca_results.items():
axis_contrast_curvse.plot(
separations_arcsec,
mag2flux_ratio(tmp_contrast_mag),
color = colors[i],
lw=1,
alpha=0.7)
i+=1
axis_contrast_curvse.set_yscale("log")
# ----------- Plot the best PCA contrast -------------------------
pca_plot = axis_contrast_curvse.plot(
separations_arcsec,
mag2flux_ratio(overall_best_pca),
color = "black",
lw=1.5,
label="Optimal PCA",
ls="dashed")
# ----------- Plot the raw contrast -------------------------
pca_plot = axis_contrast_curvse.plot(
separations_arcsec[1:],
cADI_contrast_curve.values[1:],
color = "black",
lw=1.5,
label="cADI",
ls="dotted")
# ----------- Plot the S4 contrast --------------------------
s4_plot = axis_contrast_curvse.plot(
separations_arcsec,
mag2flux_ratio(s4_contrast),
color = "darkblue",
lw=1.5,
label="${4S}_{\lambda=100}$",
ls="-")
# ------------- Double axis and limits -----------------------
lim_mag_y = (13.7, 6)
lim_arcsec_x = (0.1, 1.2)
sep_lambda_arcse = interpolate.interp1d(
separations_arcsec,
separations_FWHM,
fill_value='extrapolate')
axis_contrast_curvse_mag = axis_contrast_curvse.twinx()
axis_contrast_curvse_mag.plot(
separations_arcsec,
s4_contrast,
alpha=0.)
axis_contrast_curvse_mag.invert_yaxis()
axis_contrast_curvse_lambda = axis_contrast_curvse.twiny()
axis_contrast_curvse_lambda.plot(
separations_FWHM,
mag2flux_ratio(s4_contrast),
alpha=0.)
axis_contrast_curvse_mag.grid(
which='major', ls="dotted", lw=1,
alpha=0.2, color="black", zorder=0)
axis_contrast_curvse_mag.set_ylim(*lim_mag_y)
axis_contrast_curvse.set_ylim(
mag2flux_ratio(lim_mag_y[0]),
mag2flux_ratio(lim_mag_y[1]))
axis_contrast_curvse.set_xlim(
*lim_arcsec_x)
axis_contrast_curvse_mag.set_xlim(
*lim_arcsec_x)
axis_contrast_curvse_lambda.set_xlim(
*sep_lambda_arcse(lim_arcsec_x))
# ------------- Examples -----------------------
ex1=axis_contrast_curvse.scatter(
[example_1[0] * fwhm * pixel_scale],
[mag2flux_ratio(example_1[1])],
s=50,
color="deepskyblue",
marker="^",
zorder=10)
ex2=axis_contrast_curvse.scatter(
[example_2[0] * fwhm * pixel_scale],
[mag2flux_ratio(example_2[1])],
s=50,
color="deepskyblue",
marker="s",
zorder=10)
ex3=axis_contrast_curvse.scatter(
[example_3[0] * fwhm * pixel_scale],
[mag2flux_ratio(example_3[1])],
s=50,
color="deepskyblue",
marker="*",
zorder=10)
# ----------- Labels and fontsizes --------------------------
axis_contrast_curvse.set_xlabel(
r"Separation [arcsec]", size=10)
axis_contrast_curvse_lambda.set_xlabel(
r"Separation [$\lambda/D$]", size=10)
axis_contrast_curvse.set_ylabel(
r"Planet-to-star flux ratio", size=10)
axis_contrast_curvse_mag.set_ylabel(
r"$\Delta$ Magnitude", size=10)
axis_contrast_curvse.tick_params(
axis='both', which='major', labelsize=8)
axis_contrast_curvse_lambda.tick_params(
axis='both', which='major', labelsize=8)
axis_contrast_curvse_mag.tick_params(
axis='both', which='major', labelsize=8)
axis_contrast_curvse_mag.set_title(
r"$S/N = 5$ Contrast Limits - #5 (HD 22049)",
fontsize=12, fontweight="bold", y=1.21)
# --------------------------- Legend -----------------------
handles, labels = axis_contrast_curvse.\
get_legend_handles_labels()
labels.append("Examples")
# Create a legend with grouped markers and line plot
legend_elements = [
handles[0],
handles[1],
handles[2],
(ex1, ex2, ex3)]
# Add legend
leg1 = fig.legend(legend_elements, labels,
bbox_to_anchor=(0.585, 0.13),
fontsize=9,
handler_map={tuple: HandlerTuple(ndivide=None, pad=1.7)},
loc='lower left', ncol=4)
# PCA plots -----------------
plot1 = axis_pca_1.imshow(pca_residual_1, cmap=color_map, origin='lower')
remove_image_ticks(axis_pca_1)
annotate_pca(axis_pca_1, example_1_pcas, pca_example_1_snr)
plot2 = axis_pca_2.imshow(pca_residual_2, cmap=color_map, origin='lower')
remove_image_ticks(axis_pca_2)
annotate_pca(axis_pca_2, example_2_pcas, pca_example_2_snr)
plot3 = axis_pca_3.imshow(pca_residual_3, cmap=color_map, origin='lower')
remove_image_ticks(axis_pca_3)
annotate_pca(axis_pca_3, example_3_pcas, pca_example_3_snr)
# S4 plots -----------------
plot4 = axis_s4_1.imshow(s4_residual_1, cmap=color_map, origin='lower',
vmin=-0.04, vmax=0.05)
remove_image_ticks(axis_s4_1)
annotate_s4(axis_s4_1, s4_example_1_snr)
plot5 = axis_s4_2.imshow(s4_residual_2, cmap=color_map, origin='lower',
vmin=-0.04, vmax=0.05)
remove_image_ticks(axis_s4_2)
annotate_s4(axis_s4_2, s4_example_2_snr)
plot6 = axis_s4_3.imshow(s4_residual_3, cmap=color_map, origin='lower',
vmin=-0.04, vmax=0.05)
remove_image_ticks(axis_s4_3)
annotate_s4(axis_s4_3, s4_example_3_snr)
# Titles -----------------
axis_pca_1.set_ylabel("Optimal PCA", fontsize=12, labelpad=5)
axis_s4_1.set_ylabel("4S", fontsize=14, labelpad=5)
axis_pca_1.set_title(
r"Example: $\blacktriangle$", fontsize=12,
fontweight="bold", y=1.02)
axis_pca_2.set_title(
r"Example: $\blacksquare$", fontsize=12,
fontweight="bold", y=1.02)
axis_pca_3.set_title(
r"Example: $\bigstar$", fontsize=12,
fontweight="bold", y=1.02)
_=plt.setp(leg1.get_title(),fontsize=14)
fig.patch.set_facecolor('white')
plt.savefig("./final_plots/05_PCA_vs_S4_with_label.pdf", bbox_inches='tight')
