|
1
|
|
|
""" |
|
2
|
|
|
Module containing additional utility functions for selecting a preferred model. |
|
3
|
|
|
""" |
|
4
|
|
|
|
|
5
|
|
|
import numpy as np |
|
6
|
|
|
import yaml |
|
7
|
|
|
from astropy.table import QTable |
|
8
|
|
|
|
|
9
|
|
|
from asgardpy.config.generator import AsgardpyConfig, all_model_templates |
|
10
|
|
|
from asgardpy.stats.stats import check_model_preference_lrt |
|
11
|
|
|
|
|
12
|
|
|
__all__ = [ |
|
13
|
|
|
"fetch_all_analysis_fit_info", |
|
14
|
|
|
"get_model_config_files", |
|
15
|
|
|
"tabulate_best_fit_stats", |
|
16
|
|
|
"copy_target_config", |
|
17
|
|
|
"write_output_config_yaml", |
|
18
|
|
|
] |
|
19
|
|
|
|
|
20
|
|
|
|
|
21
|
|
|
def get_model_config_files(select_model_tags): |
|
22
|
|
|
"""From the default model templates, select some.""" |
|
23
|
|
|
|
|
24
|
|
|
all_tags, template_files = all_model_templates() |
|
25
|
|
|
|
|
26
|
|
|
spec_model_temp_files = [] |
|
27
|
|
|
for tag in select_model_tags: |
|
28
|
|
|
spec_model_temp_files.append(template_files[np.where(all_tags == tag)[0][0]]) |
|
29
|
|
|
|
|
30
|
|
|
spec_model_temp_files = np.array(spec_model_temp_files) |
|
31
|
|
|
|
|
32
|
|
|
return spec_model_temp_files |
|
33
|
|
|
|
|
34
|
|
|
|
|
35
|
|
|
def get_spec_params_indices(aa_config): |
|
36
|
|
|
""" |
|
37
|
|
|
For copying the spectral flux amplitude and flux normalization energy, |
|
38
|
|
|
from one config to another, find the correct parameter indices within a |
|
39
|
|
|
given config. |
|
40
|
|
|
""" |
|
41
|
|
|
par_names = [] |
|
42
|
|
|
for p in aa_config.config.target.components[0].spectral.parameters: |
|
43
|
|
|
par_names.append(p.name) |
|
44
|
|
|
par_names = np.array(par_names) |
|
45
|
|
|
|
|
46
|
|
|
if "amplitude" in par_names: |
|
47
|
|
|
amp_idx = np.where(par_names == "amplitude")[0][0] |
|
48
|
|
|
else: |
|
49
|
|
|
# For models without this parameter name has not yet been included or |
|
50
|
|
|
# checked with Asgardpy |
|
51
|
|
|
amp_idx = None |
|
52
|
|
|
|
|
53
|
|
|
if "reference" in par_names: |
|
54
|
|
|
eref_idx = np.where(par_names == "reference")[0][0] |
|
55
|
|
|
else: |
|
56
|
|
|
eref_idx = np.where(par_names == "ebreak")[0][0] |
|
57
|
|
|
|
|
58
|
|
|
return amp_idx, eref_idx |
|
59
|
|
|
|
|
60
|
|
|
|
|
61
|
|
|
def copy_target_config(aa_config_1, aa_config_2): |
|
62
|
|
|
"""From aa_config_1 update information in aa_config_2.""" |
|
63
|
|
|
|
|
64
|
|
|
amp_idx_1, eref_idx_1 = get_spec_params_indices(aa_config_1) |
|
65
|
|
|
amp_idx_2, eref_idx_2 = get_spec_params_indices(aa_config_2) |
|
66
|
|
|
|
|
67
|
|
|
# Have the same value of amplitude |
|
68
|
|
|
aa_config_2.config.target.components[0].spectral.parameters[amp_idx_2].value = ( |
|
69
|
|
|
aa_config_1.config.target.components[0].spectral.parameters[amp_idx_1].value |
|
70
|
|
|
) |
|
71
|
|
|
# Have the same value of reference/e_break energy |
|
72
|
|
|
aa_config_2.config.target.components[0].spectral.parameters[eref_idx_2].value = ( |
|
73
|
|
|
aa_config_1.config.target.components[0].spectral.parameters[eref_idx_1].value |
|
74
|
|
|
) |
|
75
|
|
|
# Have the same value of redshift value and EBL reference model |
|
76
|
|
|
aa_config_2.config.target.components[0].spectral.ebl_abs.redshift = aa_config_1.config.target.components[ |
|
77
|
|
|
0 |
|
78
|
|
|
].spectral.ebl_abs.redshift |
|
79
|
|
|
|
|
80
|
|
|
# Make sure the source names are the same |
|
81
|
|
|
aa_config_2.config.target.source_name = aa_config_1.config.target.source_name |
|
82
|
|
|
aa_config_2.config.target.components[0].name = aa_config_1.config.target.components[0].name |
|
83
|
|
|
|
|
84
|
|
|
return aa_config_2 |
|
85
|
|
|
|
|
86
|
|
|
|
|
87
|
|
|
def fetch_all_analysis_fit_info(main_analysis_list, spec_models_list): |
|
88
|
|
|
""" |
|
89
|
|
|
For a list of spectral models, with the AsgardpyAnalysis run till the fit |
|
90
|
|
|
step, get the relevant information for testing the model preference. |
|
91
|
|
|
""" |
|
92
|
|
|
fit_success_list = [] |
|
93
|
|
|
pref_over_pl_chi2_list = [] |
|
94
|
|
|
stat_list = [] |
|
95
|
|
|
dof_list = [] |
|
96
|
|
|
|
|
97
|
|
|
for tag in spec_models_list: |
|
98
|
|
|
dict_tag = main_analysis_list[tag]["Analysis"].instrument_spectral_info |
|
99
|
|
|
dict_pl = main_analysis_list["pl"]["Analysis"].instrument_spectral_info |
|
100
|
|
|
|
|
101
|
|
|
# Collect parameters for AIC check |
|
102
|
|
|
stat = dict_tag["best_fit_stat"] |
|
103
|
|
|
dof = dict_tag["DoF"] |
|
104
|
|
|
|
|
105
|
|
|
fit_success = main_analysis_list[tag]["Analysis"].fit_result.success |
|
106
|
|
|
|
|
107
|
|
|
fit_success_list.append(fit_success) |
|
108
|
|
|
stat_list.append(stat) |
|
109
|
|
|
dof_list.append(dof) |
|
110
|
|
|
|
|
111
|
|
|
# Checking the preference of a "nested" spectral model (observed), |
|
112
|
|
|
# over Power Law. |
|
113
|
|
|
if tag == "pl": |
|
114
|
|
|
main_analysis_list[tag]["Pref_over_pl_chi2"] = 0 |
|
115
|
|
|
main_analysis_list[tag]["Pref_over_pl_pval"] = 0 |
|
116
|
|
|
main_analysis_list[tag]["DoF_over_pl"] = 0 |
|
117
|
|
|
pref_over_pl_chi2_list.append(0) |
|
118
|
|
|
continue |
|
119
|
|
|
|
|
120
|
|
|
p_pl_x, g_pl_x, ndof_pl_x = check_model_preference_lrt( |
|
121
|
|
|
dict_pl["best_fit_stat"], |
|
122
|
|
|
dict_tag["best_fit_stat"], |
|
123
|
|
|
dict_pl["DoF"], |
|
124
|
|
|
dict_tag["DoF"], |
|
125
|
|
|
) |
|
126
|
|
|
|
|
127
|
|
|
main_analysis_list[tag]["Pref_over_pl_chi2"] = g_pl_x |
|
128
|
|
|
pref_over_pl_chi2_list.append(g_pl_x) |
|
129
|
|
|
main_analysis_list[tag]["Pref_over_pl_pval"] = p_pl_x |
|
130
|
|
|
main_analysis_list[tag]["DoF_over_pl"] = ndof_pl_x |
|
131
|
|
|
|
|
132
|
|
|
fit_success_list = np.array(fit_success_list) |
|
133
|
|
|
|
|
134
|
|
|
# Only select fit results that were successful for comparisons |
|
135
|
|
|
stat_list = np.array(stat_list)[fit_success_list] |
|
136
|
|
|
dof_list = np.array(dof_list)[fit_success_list] |
|
137
|
|
|
pref_over_pl_chi2_list = np.array(pref_over_pl_chi2_list)[fit_success_list] |
|
138
|
|
|
|
|
139
|
|
|
return fit_success_list, stat_list, dof_list, pref_over_pl_chi2_list |
|
140
|
|
|
|
|
141
|
|
|
|
|
142
|
|
|
def tabulate_best_fit_stats(spec_models_list, fit_success_list, main_analysis_list, list_rel_p): |
|
143
|
|
|
"""For a list of spectral models, tabulate the best fit information.""" |
|
144
|
|
|
|
|
145
|
|
|
fit_stats_table = [] |
|
146
|
|
|
|
|
147
|
|
|
for i, tag in enumerate(spec_models_list[fit_success_list]): |
|
148
|
|
|
info_ = main_analysis_list[tag]["Analysis"].instrument_spectral_info |
|
149
|
|
|
|
|
150
|
|
|
t = main_analysis_list[tag] |
|
151
|
|
|
|
|
152
|
|
|
ts_gof = round(info_["best_fit_stat"] - info_["max_fit_stat"], 3) |
|
153
|
|
|
t_fits = { |
|
154
|
|
|
"Spectral Model": tag.upper(), |
|
155
|
|
|
"TS of Best Fit": round(info_["best_fit_stat"], 3), |
|
156
|
|
|
"TS of Max Fit": round(info_["max_fit_stat"], 3), |
|
157
|
|
|
"TS of Goodness of Fit": ts_gof, |
|
158
|
|
|
"DoF of Fit": info_["DoF"], |
|
159
|
|
|
r"Significance ($\sigma$) of Goodness of Fit": round(info_["fit_chi2_sig"], 3), |
|
160
|
|
|
"p-value of Goodness of Fit": float(f"{info_['fit_pval']:.4g}"), |
|
161
|
|
|
"Pref over PL (chi2)": round(t["Pref_over_pl_chi2"], 3), |
|
162
|
|
|
"Pref over PL (p-value)": float(f"{t['Pref_over_pl_pval']:.4g}"), |
|
163
|
|
|
"Pref over PL (DoF)": t["DoF_over_pl"], |
|
164
|
|
|
"Relative p-value (AIC)": float(f"{list_rel_p[i]:.4g}"), |
|
165
|
|
|
} |
|
166
|
|
|
fit_stats_table.append(t_fits) |
|
167
|
|
|
stats_table = QTable(fit_stats_table) |
|
168
|
|
|
|
|
169
|
|
|
return stats_table |
|
170
|
|
|
|
|
171
|
|
|
|
|
172
|
|
|
def write_output_config_yaml(model_): |
|
173
|
|
|
"""With the selected spectral model, update a default config in yaml.""" |
|
174
|
|
|
|
|
175
|
|
|
spec_model = model_.spectral_model.model1.to_dict() |
|
176
|
|
|
|
|
177
|
|
|
temp_config = AsgardpyConfig() |
|
178
|
|
|
temp_config.target.components[0] = spec_model |
|
179
|
|
|
|
|
180
|
|
|
# Update with the spectral model info |
|
181
|
|
|
temp_ = temp_config.dict(exclude_defaults=True) |
|
182
|
|
|
|
|
183
|
|
|
# Remove some of the unnecessary keys |
|
184
|
|
|
temp_["target"].pop("models_file", None) |
|
185
|
|
|
temp_["target"]["components"][0]["spectral"].pop("ebl_abs", None) |
|
186
|
|
|
|
|
187
|
|
|
yaml_ = yaml.dump( |
|
188
|
|
|
temp_, |
|
189
|
|
|
sort_keys=False, |
|
190
|
|
|
indent=4, |
|
191
|
|
|
width=80, |
|
192
|
|
|
default_flow_style=None, |
|
193
|
|
|
) |
|
194
|
|
|
return yaml_ |
|
195
|
|
|
|
chaimain /
asgardpy
