Residual Gallery - Fake Planets#
1. Imports#
[1]:
import json
from pathlib import Path
import numpy as np
import cmocean
color_map = cmocean.cm.ice
from applefy.detections.contrast import Contrast
from applefy.utils.file_handling import load_adi_data
from applefy.utils import mag2flux_ratio
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]:
from bokeh.plotting import figure, show, output_file, save
from bokeh.models import ColumnDataSource, CustomJS, Select, Label, Div, Legend, Span
from bokeh.layouts import column, row
from bokeh.models import ColorBar, LinearColorMapper, LogColorMapper
from bokeh.models import Slider, RangeSlider
from bokeh.core.validation import silence
from bokeh.core.validation.warnings import MISSING_RENDERERS
_=silence(MISSING_RENDERERS, True)
[3]:
# Start the bokeh for jupyter notebooks
from bokeh.io import output_notebook
output_notebook()
2. Load the data#
[4]:
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
[5]:
experiment_root_dir = root_dir / Path("70_results/x1_fake_planet_experiments/" + dataset_name)
[6]:
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
[7]:
dataset_file = root_dir / Path("30_data/" + dataset_name + ".hdf5")
[8]:
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)
[9]:
# we want the image to show the innermost 1.2 arcsec
print(1.2 / pixel_scale * 2)
88.30022075055187
[10]:
# 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. Restore the Residuals from the fake planet experiments#
3.1. Create the Contrast Grid#
[11]:
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)
[12]:
# get fake planet setup
flux_ratios, separations, num_fake_planets = contrast_grid_setup_1(fwhm)
[13]:
contrast_instance.design_fake_planet_experiments(
flux_ratios=flux_ratios,
num_planets=num_fake_planets,
separations=separations,
overwrite=True)
Overwriting existing config files.
3.2 Restore the PCA results#
[27]:
# We have more residuals. But to save space on the webpage we only show some PCA residuals
pca_numbers = np.concatenate(
[np.arange(2, 30, 4)[1:],
np.arange(30, 60, 10),
np.arange(60, 200, 20)[:-2],
np.arange(200, 600, 100)])
pca_numbers
[27]:
array([ 6, 10, 14, 18, 22, 26, 30, 40, 50, 60, 80, 100, 120,
140, 200, 300, 400, 500])
[28]:
pca_algorithm_function = PCADataReductionGPU(
approx_svd=8000,
pca_numbers=pca_numbers,
device="cpu",
work_dir=None,
special_name="stacked_05",
verbose=False)
[29]:
contrast_instance.run_fake_planet_experiments(
algorithm_function=pca_algorithm_function,
num_parallel=1)
Running fake planet experiments...
100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 946/946 [00:48<00:00, 19.66it/s]
[DONE]
3.3 Restore the 4S results#
[34]:
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:01<00:00, 488.62it/s]
[DONE]
Running fake planet experiments...
100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 946/946 [00:01<00:00, 483.34it/s]
[DONE]
Running fake planet experiments...
100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 946/946 [00:02<00:00, 456.28it/s]
[DONE]
Running fake planet experiments...
100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 946/946 [00:02<00:00, 442.95it/s]
[DONE]
4. Sort find the right residuals#
[35]:
# Examples for later
example_1 = (1.5, 7.5, "a")
example_2 = (2.5, 9.5, "c")
example_3 = (11., 13., "b")
[36]:
def find_residual(
method_tag,
example):
for tmp_result in contrast_instance.results_dict[method_tag]:
tmp_config, tmp_residual = tmp_result
if tmp_config["type"] != "TP estimation":
continue
if np.round(tmp_config["separation"] / 3.6, 2) != example[0]:
continue
if tmp_config["flux_ratio"] != mag2flux_ratio(example[1]):
continue
if tmp_config["exp_id"][-1] != example[2]:
continue
return tmp_residual
4.1 PCA#
[37]:
pca_start_comp = 50
all_residuals_pca = dict()
for tmp_component in pca_numbers:
tmp_name = "stacked_05_PCA_" + str(tmp_component).zfill(3) + "_components"
all_residuals_pca[str(tmp_component) + "_e1"] = find_residual(tmp_name, example_1)
all_residuals_pca[str(tmp_component) + "_e2"] = find_residual(tmp_name, example_2)
all_residuals_pca[str(tmp_component) + "_e3"] = find_residual(tmp_name, example_3)
residual_selected_pca = {"residual" : [all_residuals_pca[str(pca_start_comp) + "_e2"]]}
start_min = np.min(all_residuals_pca[str(pca_start_comp) + "_e2"]) / 1.8
start_max = np.max(all_residuals_pca[str(pca_start_comp) + "_e2"]) * 1.5
4.2 4S#
[38]:
s4_start_comp = 1000
lambda_regs = [100, 1000, 10000, 100000]
all_residuals_4s = dict()
for tmp_lambda in lambda_regs:
tmp_name = "s4_mean_lambda_" + str(tmp_lambda).zfill(6)
all_residuals_4s[str(tmp_lambda) + "_e1"] = find_residual(tmp_name, example_1)
all_residuals_4s[str(tmp_lambda) + "_e2"] = find_residual(tmp_name, example_2)
all_residuals_4s[str(tmp_lambda) + "_e3"] = find_residual(tmp_name, example_3)
residual_selected_4s = {"residual" : [all_residuals_4s[str(s4_start_comp) + "_e2"]]}
5. Interactive Plot#
[39]:
residual_size=all_residuals_pca['100_e1'].shape[0]
fig_size=350
[43]:
# 1.) Setup all the data needed for the plot ------------------------------
source_all_residuals_pca = ColumnDataSource(data=all_residuals_pca)
source_selected_residual_pca = ColumnDataSource(data=residual_selected_pca)
source_all_residuals_4s = ColumnDataSource(data=all_residuals_4s)
source_selected_residual_4s = ColumnDataSource(data=residual_selected_4s)
# 2. Residual Plot PCA ----------------------------------------------------------
plot_residual_pca = figure(
width=fig_size, height=fig_size,
toolbar_location='left',
tools="pan,wheel_zoom,reset,save",
active_scroll="wheel_zoom",
title="Fake Planet (PCA)")
plot_residual_pca.x_range.range_padding = 0
plot_residual_pca.y_range.range_padding = 0
plot_residual_pca.xaxis.visible = False
plot_residual_pca.yaxis.visible = False
plot_residual_pca.title.text_font_size = '12pt'
color_mapper_pca = LinearColorMapper(
palette="Blues256",
low=start_min, high=start_max)
image = plot_residual_pca.image(
image="residual",
color_mapper=color_mapper_pca,
x=0, y=0,
dw=residual_size,
dh=residual_size,
source=source_selected_residual_pca)
# 2. Residual Plot 4S ----------------------------------------------------------
plot_residual_4s = figure(
width=fig_size, height=fig_size,
toolbar_location='right',
tools="pan,wheel_zoom,reset,save",
active_scroll="wheel_zoom",
title="Fake Planet (4S)")
plot_residual_4s.x_range.range_padding = 0
plot_residual_4s.y_range.range_padding = 0
plot_residual_4s.xaxis.visible = False
plot_residual_4s.yaxis.visible = False
plot_residual_4s.title.text_font_size = '12pt'
color_mapper_4s = LinearColorMapper(
palette="Blues256",
low=-0.03, high=0.1)
image = plot_residual_4s.image(
image="residual",
color_mapper=color_mapper_4s,
x=0, y=0,
dw=residual_size,
dh=residual_size,
source=source_selected_residual_4s)
example_1 = (1.5, 7.5, "a")
example_2 = (2.5, 9.5, "c")
example_3 = (11., 13., "b")
# 3.) Example selection ------------------------------------------------
example_select = Select(
title="Example",
width=fig_size * 2,
value=2,
options=[
(1, r"Separation: 1.5 lambda /D; Contrast: 7.5 [mag]"),
(2, r"Separation: 2.5 lambda /D; Contrast: 9.5 [mag]"),
(3, r"Separation: 11 lambda /D; Contrast: 13 [mag]")])
# 4.) Widgets PCA ------------------------------------------------------------
color_slider_pca = RangeSlider(
width=fig_size,
start=-20, end=20, step=.05,
value=(start_min, start_max),
title="color range")
color_slider_pca.js_link("value", color_mapper_pca, "low", attr_selector=0)
color_slider_pca.js_link("value", color_mapper_pca, "high", attr_selector=1)
# PCA components
pca_setup_select = Select(
title="PCA components",
width=fig_size,
value=str(pca_start_comp),
options=[str(i) for i in pca_numbers])
# 4.) Widgets 4S ------------------------------------------------------------
color_slider_4s = RangeSlider(
width=fig_size,
start=-0.3, end=0.3, step=.005,
value=(-0.03, 0.1),
title="color range")
color_slider_4s.js_link("value", color_mapper_4s, "low", attr_selector=0)
color_slider_4s.js_link("value", color_mapper_4s, "high", attr_selector=1)
# 4S setup
s4_setup_select = Select(
title="Lambda",
width=fig_size,
value=str(s4_start_comp),
options=[str(i) for i in lambda_regs])
# 5.) Interaction ------------------------------------------------------------
# setups
code_widget = """
// 1.) Update the PCA residual
Array.prototype.max = function() {
return Math.max.apply(null, this);
};
Array.prototype.min = function() {
return Math.min.apply(null, this);
};
const current_id = components.value + "_e" + example.value;
source_selected_residual_pca.data = {"residual": [source_all_residuals_pca.data[current_id],]};
let min_val = source_all_residuals_pca.data[current_id].min() / 1.8;
let max_val = source_all_residuals_pca.data[current_id].max() * 1.5;
color_slider.value = [min_val, max_val];
color_mapper.low = color_slider.value[0]
color_mapper.high = color_slider.value[1]
"""
callback_pca = CustomJS(
args=dict(
source_selected_residual_pca = source_selected_residual_pca,
source_all_residuals_pca = source_all_residuals_pca,
components = pca_setup_select,
color_slider = color_slider_pca,
color_mapper = color_mapper_pca,
example = example_select
),
code=code_widget)
pca_setup_select.js_on_change('value', callback_pca)
example_select.js_on_change('value', callback_pca)
# 4S
code_widget = """
// 1.) Update the 4S residual
const current_id = setup.value + "_e" + example.value;
source_selected_residual_4s.data = {"residual": [source_all_residuals_4s.data[current_id],]};
"""
callback_4s = CustomJS(
args=dict(
source_selected_residual_4s = source_selected_residual_4s,
source_all_residuals_4s = source_all_residuals_4s,
setup = s4_setup_select,
example = example_select
),
code=code_widget)
s4_setup_select.js_on_change('value', callback_4s)
example_select.js_on_change('value', callback_4s)
# 6.) Create the Layout
part_1 = row(
plot_residual_pca,
plot_residual_4s)
part_2 = row(
color_slider_pca,
color_slider_4s)
part_3 = row(
pca_setup_select,
s4_setup_select)
final_plot = column(
example_select,
part_1,
part_2,
part_3)
#output_file("../../_static/gallery_fake_planets.html")
#save(final_plot)
show(final_plot)