check_preferred_model.get_best_preferred_model_aic() - Code Metrics - Inspection of "Merge pull request #204 from chaimain/cleanup" - chaimain/asgardpy - Measure and Improve Code Quality continuously with Scrutinizer

Completed

Push — main ( 2edd10...a5707d )

by Chaitanya

created 2025-03-21 20:49 UTC

get_best_preferred_model_aic() A

↳ Parent: check_preferred_model

Complexity

Conditions

Size

Total Lines	15
Code Lines	8

Duplication

Lines	0
Ratio	0 %

Importance

Changes

Metric	Value
cc	3
eloc	8
nop	6
dl	0
loc	15
rs	10
c	0
b	0
f	0

import argparse
import logging
from pathlib import Path

import numpy as np

from asgardpy.analysis import AsgardpyAnalysis
from asgardpy.config import AsgardpyConfig, write_asgardpy_model_to_file
from asgardpy.stats import (
    check_model_preference_aic,
    copy_target_config,
    fetch_all_analysis_fit_info,
    get_model_config_files,
    tabulate_best_fit_stats,
)

log = logging.getLogger(__name__)

parser = argparse.ArgumentParser(description="Get preferred best-fit spectral model")

parser.add_argument(
    "--config",
    "-c",
    help="Path to the config file",
)

parser.add_argument("--ebl-scale-factor", "-e", help="Value of EBL Norm Scale Factor", default=1.0, type=float)

parser.add_argument(
    "--ebl-model-name",
    "-m",
    help="Name of EBL model as used by Gammapy",
    default="dominguez",
    type=str,
)

parser.add_argument(
    "--write-config",
    help="Boolean to write the best-fit model into a separate file.",
    default=True,
    type=bool,
)


def fetch_all_analysis_objects(main_config, spec_model_temp_files, ebl_scale_factor, ebl_model_name):
    """For a list of spectral models, initiate AsgardpyAnalysis objects."""
    main_analysis_list = {}
    spec_models_list = []

    for temp in spec_model_temp_files:
        temp_model = AsgardpyAnalysis(main_config)
        temp_model.config.target.models_file = temp

        temp_model_2 = AsgardpyAnalysis(temp_model.config)

        copy_target_config(temp_model, temp_model_2)

        if ebl_scale_factor != 1.0:
            temp_model_2.config.target.components[0].spectral.ebl_abs.alpha_norm = ebl_scale_factor

        if ebl_model_name != "dominguez":
            temp_model_2.config.target.components[0].spectral.ebl_abs.reference = ebl_model_name.replace("_", "-")
        else:
            temp_model_2.config.target.components[
                0
            ].spectral.ebl_abs.reference = temp_model.config.target.components[0].spectral.ebl_abs.reference

        spec_tag = temp.name.split(".")[0].split("_")[-1]
        spec_models_list.append(spec_tag)
        main_analysis_list[spec_tag] = {}

        main_analysis_list[spec_tag]["Analysis"] = temp_model_2

    spec_models_list = np.array(spec_models_list)

    return main_analysis_list, spec_models_list


def get_best_preferred_model_lrt(best_sp_idx_lrt, pref_over_pl_chi2_list, spec_models_list, PL_idx, log):
    """
    From a list of a given spectral model's preference over PL model as per LRT,
    get the index of the best spectral model and write appropriate logs.
    """
    for idx in best_sp_idx_lrt:
        if pref_over_pl_chi2_list[idx] > 5:
            sp_idx_lrt = idx
            log.info("Best preferred spectral model over PL is %s", spec_models_list[idx])
        else:
            sp_idx_lrt = PL_idx
            log.info("No other model preferred over PL")
    return sp_idx_lrt, log



def get_best_preferred_model_aic(best_sp_idx_aic, list_rel_p, spec_models_list, fit_success_list, PL_idx, log):
    """
    From a list of a given spectral model's relative p-value from a list of
    spectral models, as per AIC, get the index of the best spectral model and
    write appropriate logs.
    """
    for idx in best_sp_idx_aic:
        if list_rel_p[idx] > 0.95:
            sp_idx_aic = idx
            log.info("Best preferred spectral model is %s", spec_models_list[fit_success_list][idx])
        else:
            sp_idx_aic = PL_idx
            log.info("No other model preferred, hence PL is selected")

    return sp_idx_aic, log



def main():
    args = parser.parse_args()

    main_config = AsgardpyConfig.read(args.config)
    config_path = Path(args.config)
    config_path_file_name = config_path.name.split(".")[0]
    target_source_name = main_config.target.source_name

    steps_list = []
    for s in main_config.general.steps:
        if s != "flux-points":
            steps_list.append(s)
    log.info("Target source is: %s", target_source_name)


    spec_model_temp_files = get_model_config_files(["lp", "bpl", "ecpl", "pl", "eclp", "sbpl"])

    main_analysis_list, spec_models_list = fetch_all_analysis_objects(
        main_config, spec_model_temp_files, args.ebl_scale_factor, args.ebl_model_name
    )

    # Run Analysis Steps till Fit
    PL_idx = 0

    for i, tag in enumerate(spec_models_list):
        log.info("Spectral model being tested: %s", tag)

        main_analysis_list[tag]["Analysis"].run(steps_list)

        if tag == "pl":
            PL_idx = i

    fit_success_list, stat_list, dof_list, pref_over_pl_chi2_list = fetch_all_analysis_fit_info(
        main_analysis_list, spec_models_list
    )

    # If any spectral model has at least 5 sigmas preference over PL
    best_sp_idx_lrt = np.nonzero(pref_over_pl_chi2_list == np.nanmax(pref_over_pl_chi2_list))[0]
    sp_idx_lrt, log = get_best_preferred_model_lrt(
        best_sp_idx_lrt,
        pref_over_pl_chi2_list,
        spec_models_list,
        PL_idx,
        log,
    )

    list_rel_p = check_model_preference_aic(stat_list, dof_list)

    best_sp_idx_aic = np.nonzero(list_rel_p == np.nanmax(list_rel_p))[0]

    sp_idx_aic, log = get_best_preferred_model_aic(
        best_sp_idx_aic,
        list_rel_p,
        spec_models_list,
        fit_success_list,
        PL_idx,
        log,
    )

    stats_table = tabulate_best_fit_stats(spec_models_list, fit_success_list, main_analysis_list, list_rel_p)

    stats_table.meta["Target source name"] = target_source_name
    stats_table.meta["EBL model"] = args.ebl_model_name
    stats_table.meta["EBL scale factor"] = args.ebl_scale_factor

    file_name = f"{config_path_file_name}_{args.ebl_model_name}_{args.ebl_scale_factor}_fit_stats.ecsv"
    stats_table.write(
        main_config.general.outdir / file_name,
        format="ascii.ecsv",
        overwrite=True,
    )

    if args.write_config:
        log.info("Write the spectral model")

        for idx, name in zip([sp_idx_lrt, sp_idx_aic], ["lrt", "aic"], strict=False):
            tag = spec_models_list[fit_success_list][idx]

            path = config_path.parent / f"{config_path_file_name}_model_most_pref_{name}.yaml"

            write_asgardpy_model_to_file(
                gammapy_model=main_analysis_list[tag]["Analysis"].final_model[0], output_file=path
            )


if __name__ == "__main__":
    main()


1			import argparse
2			import logging
3			from pathlib import Path
4
5			import numpy as np
6
7			from asgardpy.analysis import AsgardpyAnalysis
8			from asgardpy.config import AsgardpyConfig, write_asgardpy_model_to_file
9			from asgardpy.stats import (
10			check_model_preference_aic,
11			copy_target_config,
12			fetch_all_analysis_fit_info,
13			get_model_config_files,
14			tabulate_best_fit_stats,
15			)
16
17			log = logging.getLogger(__name__)
18
19			parser = argparse.ArgumentParser(description="Get preferred best-fit spectral model")
20
21			parser.add_argument(
22			"--config",
23			"-c",
24			help="Path to the config file",
25			)
26
27			parser.add_argument("--ebl-scale-factor", "-e", help="Value of EBL Norm Scale Factor", default=1.0, type=float)
28
29			parser.add_argument(
30			"--ebl-model-name",
31			"-m",
32			help="Name of EBL model as used by Gammapy",
33			default="dominguez",
34			type=str,
35			)
36
37			parser.add_argument(
38			"--write-config",
39			help="Boolean to write the best-fit model into a separate file.",
40			default=True,
41			type=bool,
42			)
43
44
45			def fetch_all_analysis_objects(main_config, spec_model_temp_files, ebl_scale_factor, ebl_model_name):
46			"""For a list of spectral models, initiate AsgardpyAnalysis objects."""
47			main_analysis_list = {}
48			spec_models_list = []
49
50			for temp in spec_model_temp_files:
51			temp_model = AsgardpyAnalysis(main_config)
52			temp_model.config.target.models_file = temp
53
54			temp_model_2 = AsgardpyAnalysis(temp_model.config)
55
56			copy_target_config(temp_model, temp_model_2)
57
58			if ebl_scale_factor != 1.0:
59			temp_model_2.config.target.components[0].spectral.ebl_abs.alpha_norm = ebl_scale_factor
60
61			if ebl_model_name != "dominguez":
62			temp_model_2.config.target.components[0].spectral.ebl_abs.reference = ebl_model_name.replace("_", "-")
63			else:
64			temp_model_2.config.target.components[
65			0
66			].spectral.ebl_abs.reference = temp_model.config.target.components[0].spectral.ebl_abs.reference
67
68			spec_tag = temp.name.split(".")[0].split("_")[-1]
69			spec_models_list.append(spec_tag)
70			main_analysis_list[spec_tag] = {}
71
72			main_analysis_list[spec_tag]["Analysis"] = temp_model_2
73
74			spec_models_list = np.array(spec_models_list)
75
76			return main_analysis_list, spec_models_list
77
78
79			def get_best_preferred_model_lrt(best_sp_idx_lrt, pref_over_pl_chi2_list, spec_models_list, PL_idx, log):
80			"""
81			From a list of a given spectral model's preference over PL model as per LRT,
82			get the index of the best spectral model and write appropriate logs.
83			"""
84			for idx in best_sp_idx_lrt:
85			if pref_over_pl_chi2_list[idx] > 5:
86			sp_idx_lrt = idx
87			log.info("Best preferred spectral model over PL is %s", spec_models_list[idx])
88			else:
89			sp_idx_lrt = PL_idx
90			log.info("No other model preferred over PL")
91			return sp_idx_lrt, log
			0 ignored issues – show introduced 2024-01-06 22:55 UTC by Report Bug Copy Issue Report The variable `sp_idx_lrt` does not seem to be defined in case the `for` loop on line `84` is not entered. Are you sure this can never be the case? Loading history...
92
93
94			def get_best_preferred_model_aic(best_sp_idx_aic, list_rel_p, spec_models_list, fit_success_list, PL_idx, log):
95			"""
96			From a list of a given spectral model's relative p-value from a list of
97			spectral models, as per AIC, get the index of the best spectral model and
98			write appropriate logs.
99			"""
100			for idx in best_sp_idx_aic:
101			if list_rel_p[idx] > 0.95:
102			sp_idx_aic = idx
103			log.info("Best preferred spectral model is %s", spec_models_list[fit_success_list][idx])
104			else:
105			sp_idx_aic = PL_idx
106			log.info("No other model preferred, hence PL is selected")
107
108			return sp_idx_aic, log
			0 ignored issues – show introduced 2024-01-06 22:55 UTC by Report Bug Copy Issue Report The variable `sp_idx_aic` does not seem to be defined in case the `for` loop on line `100` is not entered. Are you sure this can never be the case? Loading history...
109
110
111			def main():
112			args = parser.parse_args()
113
114			main_config = AsgardpyConfig.read(args.config)
115			config_path = Path(args.config)
116			config_path_file_name = config_path.name.split(".")[0]
117			target_source_name = main_config.target.source_name
118
119			steps_list = []
120			for s in main_config.general.steps:
121			if s != "flux-points":
122			steps_list.append(s)
123			log.info("Target source is: %s", target_source_name)
			0 ignored issues – show Comprehensibility Best Practice introduced 2025-03-21 20:30 UTC by Report Bug Copy Issue Report The variable `log` does not seem to be defined. Loading history...
124
125			spec_model_temp_files = get_model_config_files(["lp", "bpl", "ecpl", "pl", "eclp", "sbpl"])
126
127			main_analysis_list, spec_models_list = fetch_all_analysis_objects(
128			main_config, spec_model_temp_files, args.ebl_scale_factor, args.ebl_model_name
129			)
130
131			# Run Analysis Steps till Fit
132			PL_idx = 0
133
134			for i, tag in enumerate(spec_models_list):
135			log.info("Spectral model being tested: %s", tag)
136
137			main_analysis_list[tag]["Analysis"].run(steps_list)
138
139			if tag == "pl":
140			PL_idx = i
141
142			fit_success_list, stat_list, dof_list, pref_over_pl_chi2_list = fetch_all_analysis_fit_info(
143			main_analysis_list, spec_models_list
144			)
145
146			# If any spectral model has at least 5 sigmas preference over PL
147			best_sp_idx_lrt = np.nonzero(pref_over_pl_chi2_list == np.nanmax(pref_over_pl_chi2_list))[0]
148			sp_idx_lrt, log = get_best_preferred_model_lrt(
149			best_sp_idx_lrt,
150			pref_over_pl_chi2_list,
151			spec_models_list,
152			PL_idx,
153			log,
154			)
155
156			list_rel_p = check_model_preference_aic(stat_list, dof_list)
157
158			best_sp_idx_aic = np.nonzero(list_rel_p == np.nanmax(list_rel_p))[0]
159
160			sp_idx_aic, log = get_best_preferred_model_aic(
161			best_sp_idx_aic,
162			list_rel_p,
163			spec_models_list,
164			fit_success_list,
165			PL_idx,
166			log,
167			)
168
169			stats_table = tabulate_best_fit_stats(spec_models_list, fit_success_list, main_analysis_list, list_rel_p)
170
171			stats_table.meta["Target source name"] = target_source_name
172			stats_table.meta["EBL model"] = args.ebl_model_name
173			stats_table.meta["EBL scale factor"] = args.ebl_scale_factor
174
175			file_name = f"{config_path_file_name}_{args.ebl_model_name}_{args.ebl_scale_factor}_fit_stats.ecsv"
176			stats_table.write(
177			main_config.general.outdir / file_name,
178			format="ascii.ecsv",
179			overwrite=True,
180			)
181
182			if args.write_config:
183			log.info("Write the spectral model")
184
185			for idx, name in zip([sp_idx_lrt, sp_idx_aic], ["lrt", "aic"], strict=False):
186			tag = spec_models_list[fit_success_list][idx]
187
188			path = config_path.parent / f"{config_path_file_name}_model_most_pref_{name}.yaml"
189
190			write_asgardpy_model_to_file(
191			gammapy_model=main_analysis_list[tag]["Analysis"].final_model[0], output_file=path
192			)
193
194
195			if __name__ == "__main__":
196			main()
197

chaimain / asgardpy

Push — main ( 2edd10...a5707d )

get_best_preferred_model_aic() A

Complexity

Size

Duplication

Importance

Duplication Side-by-Side

Filter issues like