|
1
|
|
|
import numpy as np |
|
2
|
|
|
|
|
3
|
|
|
from asgardpy.analysis import AsgardpyAnalysis |
|
4
|
|
|
from asgardpy.stats import ( |
|
5
|
|
|
check_model_preference_aic, |
|
6
|
|
|
check_model_preference_lrt, |
|
7
|
|
|
copy_target_config, |
|
8
|
|
|
fetch_all_analysis_fit_info, |
|
9
|
|
|
get_model_config_files, |
|
10
|
|
|
tabulate_best_fit_stats, |
|
11
|
|
|
) |
|
12
|
|
|
|
|
13
|
|
|
|
|
14
|
|
|
def test_preferred_model(base_config_1d): |
|
15
|
|
|
""" |
|
16
|
|
|
Testing the script code of checking the preferred spectral model. |
|
17
|
|
|
""" |
|
18
|
|
|
select_model_tags = ["lp", "bpl2", "ecpl", "pl", "eclp"] |
|
19
|
|
|
spec_model_temp_files = [] |
|
20
|
|
|
spec_model_temp_files = get_model_config_files(select_model_tags) |
|
21
|
|
|
|
|
22
|
|
|
main_analysis_list = {} |
|
23
|
|
|
spec_models_list = [] |
|
24
|
|
|
|
|
25
|
|
|
for temp in spec_model_temp_files: |
|
26
|
|
|
temp_model = AsgardpyAnalysis(base_config_1d) |
|
27
|
|
|
temp_model.config.fit_params.fit_range.min = "100 GeV" |
|
28
|
|
|
|
|
29
|
|
|
temp_model.config.target.models_file = temp |
|
30
|
|
|
|
|
31
|
|
|
temp_model_2 = AsgardpyAnalysis(temp_model.config) |
|
32
|
|
|
|
|
33
|
|
|
copy_target_config(temp_model, temp_model_2) |
|
34
|
|
|
|
|
35
|
|
|
spec_tag = temp.name.split(".")[0].split("_")[-1] |
|
36
|
|
|
spec_models_list.append(spec_tag) |
|
37
|
|
|
main_analysis_list[spec_tag] = {} |
|
38
|
|
|
|
|
39
|
|
|
main_analysis_list[spec_tag]["Analysis"] = temp_model_2 |
|
40
|
|
|
|
|
41
|
|
|
spec_models_list = np.array(spec_models_list) |
|
42
|
|
|
|
|
43
|
|
|
# Run Analysis Steps till Fit |
|
44
|
|
|
for tag in spec_models_list: |
|
45
|
|
|
main_analysis_list[tag]["Analysis"].run(["datasets-1d", "fit"]) |
|
46
|
|
|
|
|
47
|
|
|
fit_success_list, stat_list, dof_list, pref_over_pl_chi2_list = fetch_all_analysis_fit_info( |
|
48
|
|
|
main_analysis_list, spec_models_list |
|
49
|
|
|
) |
|
50
|
|
|
|
|
51
|
|
|
# If any spectral model has at least 5 sigmas preference over PL |
|
52
|
|
|
best_sp_idx_lrt = np.nonzero(pref_over_pl_chi2_list == np.nanmax(pref_over_pl_chi2_list))[0] |
|
53
|
|
|
for idx in best_sp_idx_lrt: |
|
54
|
|
|
if pref_over_pl_chi2_list[idx] > 5: |
|
55
|
|
|
lrt_best_model = spec_models_list[idx] |
|
56
|
|
|
|
|
57
|
|
|
list_rel_p = check_model_preference_aic(stat_list, dof_list) |
|
58
|
|
|
|
|
59
|
|
|
best_sp_idx_aic = np.nonzero(list_rel_p == np.nanmax(list_rel_p))[0] |
|
60
|
|
|
|
|
61
|
|
|
aic_best_model = select_model_tags[best_sp_idx_aic[0]] |
|
62
|
|
|
|
|
63
|
|
|
stats_table = tabulate_best_fit_stats(spec_models_list, fit_success_list, main_analysis_list, list_rel_p) |
|
64
|
|
|
|
|
65
|
|
|
assert lrt_best_model == "lp" |
|
66
|
|
|
assert aic_best_model == "lp" |
|
67
|
|
|
assert len(stats_table.colnames) == 11 |
|
68
|
|
|
|
|
69
|
|
|
# Check for bad comparisons, same dof |
|
70
|
|
|
p_val_0, g_sig_0, dof_0 = check_model_preference_lrt(4.4, 2.2, 2, 2) |
|
71
|
|
|
|
|
72
|
|
|
assert np.isnan(p_val_0) |
|
73
|
|
|
|
chaimain /
asgardpy
