3998833878

Build Type

github

Committed by GitHub

Commit Message

Merge 0a1ac6757 into 7f688e005

Pull Request Pull Request #246: ResultsHandler improvements

Run Details

17 of 17 new or added lines in 2 files covered. (100.0%)

4394 of 5752 relevant lines covered (76.39%)

2.29 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

45.09

/flarestack/core/unblinding.py

import logging
import sys
import os
import numpy as np
from flarestack.core.minimisation import (
    MinimisationHandler,
    read_mh_dict,
    FlareMinimisationHandler,
)
from flarestack.core.injector import MockUnblindedInjector, TrueUnblindedInjector
from flarestack.core.results import ResultsHandler, OverfluctuationError
from flarestack.core.time_pdf import TimePDF
from flarestack.shared import (
    name_pickle_output_dir,
    plot_output_dir,
    analysis_pickle_path,
    limit_output_path,
    unblinding_output_path,
)
import pickle
from flarestack.core.ts_distributions import (
    plot_background_ts_distribution,
    get_ts_fit_type,
)
import matplotlib.pyplot as plt
from flarestack.utils.catalogue_loader import load_catalogue

logger = logging.getLogger(__name__)


def confirm():
    print("Is this correct? (y/n)")

    x = ""

    while x not in ["y", "n"]:
        x = input("")

    if x == "n":
        logger.warning("Please check carefully before unblinding!")
        sys.exit()


def create_unblinder(
    unblind_dict,
    mock_unblind=True,
    full_plots=False,
    disable_warning=False,
    reunblind=True,
    **kwargs,
):
    """Dynamically create an Unblinder class that inherits corectly from the
    appropriate MinimisationHandler. The name of the parent is specified in
    the unblinder dictionary as 'mh_name'.

    :param unblind_dict: Dictionary containing Unblinding arguments
    :param mock_unblind: By default, the unblinder returns a fixed-seed
    background scramble. mock_unblind must be explicitly set to False for
    unblinding to occur.
    :param full_plots: Boolean, determines whether likelihood scans and
    limits be generated (can be computationally expensive)
    :param disable_warning: By default, the unblinder gives a warning if real
    data is unblinded. This can be disabled.
    :param reunblind: By default, the unblinder performs the analysis and
    save the ns, gamma and TS results in a pickle file. If set to False, the
    pickle file is loaded and unblinding not performed.
    :return: Instance of dynamically-generated Unblinder class
    """

    unblind_dict = read_mh_dict(unblind_dict)

    try:
        mh_name = unblind_dict["mh_name"]
    except KeyError:
        raise KeyError("No MinimisationHandler specified.")

    # Set up dynamic inheritance

    try:
        ParentMiminisationHandler = MinimisationHandler.subclasses[mh_name]
    except KeyError:
        raise KeyError("Parent class {} not found.".format(mh_name))

    # Defines custom Unblinder class

    class Unblinder(ParentMiminisationHandler):
        def __init__(self, unblind_dict, seed=None, scan_2d=False, **kwargs):

            self.unblind_dict = unblind_dict
            unblind_dict["unblind_bool"] = True
            unblind_dict["mock_unblind_bool"] = mock_unblind
            unblind_dict["inj_dict"] = {}

            # specifies whether all available background distributions should be merged
            # to improve statistics (see loop over subdirs in self.compare_to_background_TS()
            self.use_merged_background_ts_distributions = kwargs.get(
                "merge_background_ts_distr", True
            )
            self.background_median = np.nan

            # gives the medians for the separate background distributions
            # this is used when background TS distributions are not merged
            self.background_median_gamma = dict()

            if np.logical_and(not mock_unblind, not disable_warning):
                self.check_unblind()

            if mock_unblind is False:
                self.mock_unblind = False
            else:
                self.mock_unblind = True

            ParentMiminisationHandler.__init__(self, unblind_dict)

            try:
                if self.mock_unblind:
                    self.limit_path = limit_output_path(
                        self.unblind_dict["name"] + "mock_unblind/"
                    )
                    self.unblind_res_path = unblinding_output_path(
                        self.unblind_dict["name"] + "mock_unblind/"
                    )

                else:
                    self.limit_path = limit_output_path(
                        self.unblind_dict["name"] + "real_unblind/"
                    )
                    self.unblind_res_path = unblinding_output_path(
                        self.unblind_dict["name"] + "real_unblind/"
                    )

            except KeyError:
                self.limit_path = np.nan
                self.unblind_res_path = np.nan

            if self.name != " /":
                if self.mock_unblind:
                    self.name += "mock_unblind/"
                else:
                    self.name += "real_unblind/"

            self.plot_dir = plot_output_dir(self.name)

            try:
                os.makedirs(os.path.dirname(self.unblind_res_path))
            except OSError:
                pass

            # Avoid redoing unblinding

            if not reunblind:

                try:
                    logger.info(
                        "Not re-unblinding, loading results from {0}".format(
                            self.unblind_res_path
                        )
                    )

                    with open(self.unblind_res_path, "rb") as f:
                        self.res_dict = pickle.load(f)

                except FileNotFoundError:
                    logger.warning(
                        "No pickle files containing unblinding results found. "
                        "Re-unblinding now."
                    )

            # Otherwise unblind

            if not hasattr(self, "res_dict"):
                logger.info("Unblinding catalogue")

                # Minimise likelihood and produce likelihood scans
                self.res_dict = self.simulate_and_run(0, seed)

            logger.info(self.res_dict)

            # Quantify the TS value significance
            self.bkg_ts_array = None
            self.ts = np.array([self.res_dict["TS"]])[0]
            self.sigma = np.nan
            # the raw pre-trial p-value will be calculated using only the TS distributions without fitting them
            # This is only valid if you ran enough background trials!
            self.raw_pre_trial_pvalue = np.nan
            self.ts_type = get_ts_fit_type(unblind_dict)

            logger.info("Test Statistic of: {0}".format(self.ts))

            ub_res_dict = {
                "Parameters": self.res_dict["Parameters"],
                "TS": self.res_dict["TS"],
                "Flag": self.res_dict["Flag"],
            }

            logger.info(
                "Saving unblinding results to {0}".format(self.unblind_res_path)
            )

            with open(self.unblind_res_path, "wb") as f:
                pickle.dump(ub_res_dict, f)

            try:
                path = self.unblind_dict["background_ts"]
                self.pickle_dir = name_pickle_output_dir(path)
                self.output_file = self.plot_dir + "TS.pdf"
                self.compare_to_background_TS()
            except KeyError as e:
                logger.warning(
                    f"No Background TS Distribution specified. "
                    f"Cannot assess significance of TS value: {e}!"
                )

            if full_plots:

                self.calculate_upper_limits()
                if isinstance(self, FlareMinimisationHandler):
                    self.neutrino_lightcurve()
                else:
                    self.scan_likelihood(scan_2d=scan_2d)

        def add_injector(self, season, sources):
            if self.mock_unblind is False:
                return TrueUnblindedInjector(season, sources, **self.inj_dict)
            else:
                return MockUnblindedInjector(season, sources, **self.inj_dict)

        def calculate_upper_limits(self):

            try:
                ul_dir = os.path.join(self.plot_dir, "upper_limits/")

                try:
                    os.makedirs(ul_dir)
                except OSError:
                    pass

                flux_uls = []
                fluence_uls = []
                energy_flux = []
                x_axis = []

                # e_per_source_uls = []

                for subdir in sorted(os.listdir(self.pickle_dir)):
                    # this is a loop over different `gamma` values!

                    root = os.path.join(self.unblind_dict["background_ts"], subdir)

                    new_path = analysis_pickle_path(name=root)

                    if os.path.isfile(new_path):
                        logger.info(f"Opening file {new_path}")

                        with open(new_path, "rb") as f:
                            mh_dict = pickle.load(f)
                            e_pdf_dict = mh_dict["inj_dict"]["injection_energy_pdf"]

                    else:
                        logger.warning(f"Could not load {new_path}! Not a file!")
                        continue

                    self.unblind_dict["name"] = mh_dict["name"]

                    if "scale" not in self.unblind_dict:
                        self.unblind_dict["scale"] = 0.0

                    rh = ResultsHandler(self.unblind_dict)

                    logger.debug(f"In calculate_upper_limits, ResultsHandler is {rh}")

                    savepath = ul_dir + subdir + ".pdf"

                    if self.ts < self.bkg_median:
                        logging.warning(
                            f"Specified TS ({self.ts}) less than the background median ({self.bkg_median}). "
                            f"There was thus a background underfluctuation. "
                            f"Using the sensitivity for an upper limit."
                        )

                    bkg_median = (
                        self.background_median
                        if self.use_merged_background_ts_distributions
                        else self.background_median_gamma[subdir]
                    )
                    logger.debug(f"background median is {bkg_median}")
                    ul, err, extrapolated = rh.find_overfluctuations(
                        max(self.ts, bkg_median), savepath
                    )

                    flux_uls.append(ul)

                    # Calculate mean injection time per source

                    n_sources = float(len(self.sources))
                    seasons = mh_dict["dataset"]

                    if "injection_time_seconds" in self.unblind_dict:
                        inj_time = self.unblind_dict["injection_time_seconds"]
                        logger.debug(
                            f"using provided injection time of {inj_time} seconds"
                        )

                    else:
                        logger.debug("calculating injection time from TimePDF")
                        inj_time = 0.0
                        for (name, season) in seasons.items():
                            t_pdf = TimePDF.create(
                                mh_dict["inj_dict"]["injection_sig_time_pdf"],
                                season.get_time_pdf(),
                            )
                            for src in self.sources:
                                inj_time += t_pdf.raw_injection_time(src) / n_sources

                    astro_dict = rh.nu_astronomy(ul, e_pdf_dict)

                    key = "Energy Flux (GeV cm^{-2} s^{-1})"

                    fluence_uls.append(astro_dict[key] * inj_time)

                    # this functionality has been retired
                    # e_per_source_uls.append(
                    #    astro_dict["Mean Luminosity (erg/s)"] * inj_time
                    # )

                    energy_flux.append(astro_dict[key])
                    x_axis.append(float(subdir))

                plt.figure()
                ax = plt.subplot(111)
                plt.plot(x_axis, flux_uls, marker="o", label="upper limit")

                if self.mock_unblind:
                    ax.text(
                        0.2,
                        0.5,
                        "MOCK DATA",
                        color="grey",
                        alpha=0.5,
                        transform=ax.transAxes,
                    )

                ax.set_xlabel("Gamma")
                ax.set_ylabel("Flux")
                plt.yscale("log")
                plt.savefig(self.plot_dir + "upper_limit_flux.pdf")
                plt.close()

                plt.figure()
                ax1 = plt.subplot(111)
                # ax2 = ax1.twinx()

                ax1.plot(x_axis, fluence_uls, marker="o")
                # ax2.plot(x_axis, e_per_source_uls, marker="o")

                if self.mock_unblind:
                    ax1.text(
                        0.2,
                        0.5,
                        "MOCK DATA",
                        color="grey",
                        alpha=0.5,
                        transform=ax1.transAxes,
                    )

                ax1.grid(True, which="both")
                # ax2.grid(True, which="both")
                ax1.set_xlabel("Gamma")
                # ax2.set_xlabel("Gamma")
                ax1.set_ylabel(r"Total Fluence [GeV cm$^{-2}$ s$^{-1}$]")
                # ax2.set_ylabel(r"Isotropic-Equivalent Luminosity $L_{\nu}$ (erg/s)")
                ax1.set_yscale("log")
                # ax2.set_yscale("log")

                axes = [ax1]  # was: axes = [ax1, ax2]
                ordinates = [fluence_uls]  # was: [fluence_uls, e_per_source_uls]

                for k, ax in enumerate(axes):
                    y = ordinates[k]
                    ax.set_ylim(0.95 * min(y), 1.1 * max(y))

                plt.tight_layout()
                plt.savefig(self.plot_dir + "upper_limit_fluence.pdf")
                plt.close()

                try:
                    os.makedirs(os.path.dirname(self.limit_path))
                except OSError:
                    pass
                logger.info("Saving limits to {0}".format(self.limit_path))

                res_dict = {
                    "gamma": x_axis,
                    "flux": flux_uls,
                    "energy_flux": energy_flux,
                    "fluence": fluence_uls,
                    # "energy": e_per_source_uls,
                }

                with open(self.limit_path, "wb") as f:
                    pickle.dump(res_dict, f)

            except AttributeError as e:
                logger.warning(f"Unable to set limits. No TS distributions found.: {e}")

            except ValueError as e:
                raise OverfluctuationError(
                    "Insufficent pseudotrial values above and below threshold"
                )

        def compare_to_background_TS(self):
            logger.info(
                "Retrieving Background TS Distribution from {0}".format(self.pickle_dir)
            )

            try:

                ts_array = list()

                for subdir in os.listdir(self.pickle_dir):

                    merged_pkl = os.path.join(
                        os.path.join(self.pickle_dir, subdir), "merged/0.pkl"
                    )

                    logger.info("Loading {0}".format(merged_pkl))

                    try:
                        with open(merged_pkl, "rb") as mp:
                            merged_data = pickle.load(mp)
                    except (IOError, OSError) as e:
                        if len(ts_array) > 0:
                            logger.warning(e)
                        else:
                            raise OSError(e)

                    ts_array += list(merged_data["TS"])
                    self.background_median_gamma[subdir] = np.median(merged_data["TS"])

                ts_array = np.array(ts_array)
                self.bkg_ts_array = ts_array

                self.sigma = plot_background_ts_distribution(
                    ts_array, self.output_file, self.ts_type, self.ts, self.mock_unblind
                )

                self.background_median = np.median(ts_array)
                self.raw_pre_trial_pvalue = len(ts_array[ts_array > self.ts]) / len(
                    ts_array
                )

            except (IOError, OSError):

                try:
                    ts_array = list()

                    merged_pkl = os.path.join(self.pickle_dir, "merged/0.pkl")

                    logger.info("No subfolders found, loading {0}".format(merged_pkl))

                    with open(merged_pkl, "rb") as mp:
                        merged_data = pickle.load(mp)

                    ts_array += list(merged_data["TS"])

                    ts_array = np.array(ts_array)
                    self.bkg_ts_array = ts_array

                    self.sigma = plot_background_ts_distribution(
                        ts_array,
                        self.output_file,
                        self.ts_type,
                        self.ts,
                        self.mock_unblind,
                    )

                    self.raw_pre_trial_pvalue = len(ts_array[ts_array > self.ts]) / len(
                        ts_array
                    )
                    self.background_median = np.median(ts_array)

                except (IOError, OSError):
                    logger.warning("No Background TS Distribution found")
                    pass

        def check_unblind(self):
            print("\n")
            print("You are proposing to unblind data.")
            print("\n")
            confirm()
            print("\n")
            print("You are unblinding the following catalogue:")
            print("\n")
            print(self.unblind_dict["catalogue"])
            print("\n")
            confirm()
            print("\n")
            print("The catalogue has the following entries:")
            print("\n")

            cat = load_catalogue(self.unblind_dict["catalogue"])

            print(cat.dtype.names)
            print(len(cat), np.sum(cat["base_weight"]))
            print(cat)
            print("\n")
            confirm()
            print("\n")
            print("The following datasets will be used:")
            print("\n")
            for x in self.unblind_dict["dataset"].values():
                print(x.sample_name, x.season_name)
                print("\n")
                print(x.exp_path)
                print(x.pseudo_mc_path)
                print("\n")
            confirm()
            print("\n")
            print("The following LLH will be used:")
            print("\n")
            for (key, val) in self.unblind_dict["llh_dict"].items():
                print(key, val)
            print("\n")
            confirm()
            print("\n")
            print("Are you really REALLY sure about this?")
            print("You will unblind. This is your final warning.")
            print("\n")
            confirm()
            print("\n")
            print("OK, you asked for it...")
            print("\n")

        def dump_results(self, *args, **kwargs):
            logger.info(f"Using Unblinder, not dumping results!")

    return Unblinder(unblind_dict, **kwargs)

1	import logging	3✔
2	import sys	3✔
3	import os	3✔
4	import numpy as np	3✔
5	from flarestack.core.minimisation import (	3✔
6	MinimisationHandler,
7	read_mh_dict,
8	FlareMinimisationHandler,
9	)
10	from flarestack.core.injector import MockUnblindedInjector, TrueUnblindedInjector	3✔
11	from flarestack.core.results import ResultsHandler, OverfluctuationError	3✔
12	from flarestack.core.time_pdf import TimePDF	3✔
13	from flarestack.shared import (	3✔
14	name_pickle_output_dir,
15	plot_output_dir,
16	analysis_pickle_path,
17	limit_output_path,
18	unblinding_output_path,
19	)
20	import pickle	3✔
21	from flarestack.core.ts_distributions import (	3✔
22	plot_background_ts_distribution,
23	get_ts_fit_type,
24	)
25	import matplotlib.pyplot as plt	3✔
26	from flarestack.utils.catalogue_loader import load_catalogue	3✔
27
28	logger = logging.getLogger(__name__)	3✔
29
30
31	def confirm():	3✔
32	print("Is this correct? (y/n)")	×
33
34	x = ""	×
35
36	while x not in ["y", "n"]:	×
37	x = input("")	×
38
39	if x == "n":	×
40	logger.warning("Please check carefully before unblinding!")	×
41	sys.exit()	×
42
43
44	def create_unblinder(	3✔
45	unblind_dict,
46	mock_unblind=True,
47	full_plots=False,
48	disable_warning=False,
49	reunblind=True,
50	**kwargs,
51	):
52	"""Dynamically create an Unblinder class that inherits corectly from the
53	appropriate MinimisationHandler. The name of the parent is specified in
54	the unblinder dictionary as 'mh_name'.
55
56	:param unblind_dict: Dictionary containing Unblinding arguments
57	:param mock_unblind: By default, the unblinder returns a fixed-seed
58	background scramble. mock_unblind must be explicitly set to False for
59	unblinding to occur.
60	:param full_plots: Boolean, determines whether likelihood scans and
61	limits be generated (can be computationally expensive)
62	:param disable_warning: By default, the unblinder gives a warning if real
63	data is unblinded. This can be disabled.
64	:param reunblind: By default, the unblinder performs the analysis and
65	save the ns, gamma and TS results in a pickle file. If set to False, the
66	pickle file is loaded and unblinding not performed.
67	:return: Instance of dynamically-generated Unblinder class
68	"""
69
70	unblind_dict = read_mh_dict(unblind_dict)	3✔
71
72	try:	3✔
73	mh_name = unblind_dict["mh_name"]	3✔
74	except KeyError:	×
75	raise KeyError("No MinimisationHandler specified.")	×
76
77	# Set up dynamic inheritance
78
79	try:	3✔
80	ParentMiminisationHandler = MinimisationHandler.subclasses[mh_name]	3✔
81	except KeyError:	×
82	raise KeyError("Parent class {} not found.".format(mh_name))	×
83
84	# Defines custom Unblinder class
85
86	class Unblinder(ParentMiminisationHandler):	3✔
87	def __init__(self, unblind_dict, seed=None, scan_2d=False, **kwargs):	3✔
88
89	self.unblind_dict = unblind_dict	3✔
90	unblind_dict["unblind_bool"] = True	3✔
91	unblind_dict["mock_unblind_bool"] = mock_unblind	3✔
92	unblind_dict["inj_dict"] = {}	3✔
93
94	# specifies whether all available background distributions should be merged
95	# to improve statistics (see loop over subdirs in self.compare_to_background_TS()
96	self.use_merged_background_ts_distributions = kwargs.get(	3✔
97	"merge_background_ts_distr", True
98	)
99	self.background_median = np.nan	3✔
100
101	# gives the medians for the separate background distributions
102	# this is used when background TS distributions are not merged
103	self.background_median_gamma = dict()	3✔
104
105	if np.logical_and(not mock_unblind, not disable_warning):	3✔
106	self.check_unblind()	×
107
108	if mock_unblind is False:	3✔
109	self.mock_unblind = False	×
110	else:
111	self.mock_unblind = True	3✔
112
113	ParentMiminisationHandler.__init__(self, unblind_dict)	3✔
114
115	try:	3✔
116	if self.mock_unblind:	3✔
117	self.limit_path = limit_output_path(	3✔
118	self.unblind_dict["name"] + "mock_unblind/"
119	)
120	self.unblind_res_path = unblinding_output_path(	3✔
121	self.unblind_dict["name"] + "mock_unblind/"
122	)
123
124	else:
125	self.limit_path = limit_output_path(	×
126	self.unblind_dict["name"] + "real_unblind/"
127	)
128	self.unblind_res_path = unblinding_output_path(	×
129	self.unblind_dict["name"] + "real_unblind/"
130	)
131
132	except KeyError:	×
133	self.limit_path = np.nan	×
134	self.unblind_res_path = np.nan	×
135
136	if self.name != " /":	3✔
137	if self.mock_unblind:	3✔
138	self.name += "mock_unblind/"	3✔
139	else:
140	self.name += "real_unblind/"	×
141
142	self.plot_dir = plot_output_dir(self.name)	3✔
143
144	try:	3✔
145	os.makedirs(os.path.dirname(self.unblind_res_path))	3✔
146	except OSError:	3✔
147	pass	3✔
148
149	# Avoid redoing unblinding
150
151	if not reunblind:	3✔
152
153	try:	×
154	logger.info(	×
155	"Not re-unblinding, loading results from {0}".format(
156	self.unblind_res_path
157	)
158	)
159
160	with open(self.unblind_res_path, "rb") as f:	×
161	self.res_dict = pickle.load(f)	×
162
163	except FileNotFoundError:	×
164	logger.warning(	×
165	"No pickle files containing unblinding results found. "
166	"Re-unblinding now."
167	)
168
169	# Otherwise unblind
170
171	if not hasattr(self, "res_dict"):	3✔
172	logger.info("Unblinding catalogue")	3✔
173
174	# Minimise likelihood and produce likelihood scans
175	self.res_dict = self.simulate_and_run(0, seed)	3✔
176
177	logger.info(self.res_dict)	3✔
178
179	# Quantify the TS value significance
180	self.bkg_ts_array = None	3✔
181	self.ts = np.array([self.res_dict["TS"]])[0]	3✔
182	self.sigma = np.nan	3✔
183	# the raw pre-trial p-value will be calculated using only the TS distributions without fitting them
184	# This is only valid if you ran enough background trials!
185	self.raw_pre_trial_pvalue = np.nan	3✔
186	self.ts_type = get_ts_fit_type(unblind_dict)	3✔
187
188	logger.info("Test Statistic of: {0}".format(self.ts))	3✔
189
190	ub_res_dict = {	3✔
191	"Parameters": self.res_dict["Parameters"],
192	"TS": self.res_dict["TS"],
193	"Flag": self.res_dict["Flag"],
194	}
195
196	logger.info(	3✔
197	"Saving unblinding results to {0}".format(self.unblind_res_path)
198	)
199
200	with open(self.unblind_res_path, "wb") as f:	3✔
201	pickle.dump(ub_res_dict, f)	3✔
202
203	try:	3✔
204	path = self.unblind_dict["background_ts"]	3✔
205	self.pickle_dir = name_pickle_output_dir(path)	3✔
206	self.output_file = self.plot_dir + "TS.pdf"	3✔
207	self.compare_to_background_TS()	3✔
208	except KeyError as e:	3✔
209	logger.warning(	3✔
210	f"No Background TS Distribution specified. "
211	f"Cannot assess significance of TS value: {e}!"
212	)
213
214	if full_plots:	3✔
215
216	self.calculate_upper_limits()	3✔
217	if isinstance(self, FlareMinimisationHandler):	3✔
218	self.neutrino_lightcurve()	3✔
219	else:
220	self.scan_likelihood(scan_2d=scan_2d)	3✔
221
222	def add_injector(self, season, sources):	3✔
223	if self.mock_unblind is False:	3✔
224	return TrueUnblindedInjector(season, sources, **self.inj_dict)	×
225	else:
226	return MockUnblindedInjector(season, sources, **self.inj_dict)	3✔
227
228	def calculate_upper_limits(self):	3✔
229
230	try:	3✔
231	ul_dir = os.path.join(self.plot_dir, "upper_limits/")	3✔
232
233	try:	3✔
234	os.makedirs(ul_dir)	3✔
235	except OSError:	3✔
236	pass	3✔
237
238	flux_uls = []	3✔
239	fluence_uls = []	3✔
240	energy_flux = []	3✔
241	x_axis = []	3✔
242
243	# e_per_source_uls = []
244
245	for subdir in sorted(os.listdir(self.pickle_dir)):	3✔
246	# this is a loop over different `gamma` values!
247
248	root = os.path.join(self.unblind_dict["background_ts"], subdir)	3✔
249
250	new_path = analysis_pickle_path(name=root)	3✔
251
252	if os.path.isfile(new_path):	3✔
253	logger.info(f"Opening file {new_path}")	3✔
254
255	with open(new_path, "rb") as f:	3✔
256	mh_dict = pickle.load(f)	3✔
257	e_pdf_dict = mh_dict["inj_dict"]["injection_energy_pdf"]	3✔
258
259	else:
260	logger.warning(f"Could not load {new_path}! Not a file!")	×
261	continue	×
262
263	self.unblind_dict["name"] = mh_dict["name"]	3✔
264
265	if "scale" not in self.unblind_dict:	3✔
266	self.unblind_dict["scale"] = 0.0	×
267
268	rh = ResultsHandler(self.unblind_dict)	3✔
269
270	logger.debug(f"In calculate_upper_limits, ResultsHandler is {rh}")	3✔
271
272	savepath = ul_dir + subdir + ".pdf"	3✔
273
274	if self.ts < self.bkg_median:	3✔
275	logging.warning(	×
276	f"Specified TS ({self.ts}) less than the background median ({self.bkg_median}). "
277	f"There was thus a background underfluctuation. "
278	f"Using the sensitivity for an upper limit."
279	)
280
281	bkg_median = (	×
282	self.background_median
283	if self.use_merged_background_ts_distributions
284	else self.background_median_gamma[subdir]
285	)
286	logger.debug(f"background median is {bkg_median}")	×
287	ul, err, extrapolated = rh.find_overfluctuations(	×
288	max(self.ts, bkg_median), savepath
289	)
290
291	flux_uls.append(ul)	×
292
293	# Calculate mean injection time per source
294
295	n_sources = float(len(self.sources))	×
296	seasons = mh_dict["dataset"]	×
297
298	if "injection_time_seconds" in self.unblind_dict:	×
299	inj_time = self.unblind_dict["injection_time_seconds"]	×
300	logger.debug(	×
301	f"using provided injection time of {inj_time} seconds"
302	)
303
304	else:
305	logger.debug("calculating injection time from TimePDF")	×
306	inj_time = 0.0	×
307	for (name, season) in seasons.items():	×
308	t_pdf = TimePDF.create(	×
309	mh_dict["inj_dict"]["injection_sig_time_pdf"],
310	season.get_time_pdf(),
311	)
312	for src in self.sources:	×
313	inj_time += t_pdf.raw_injection_time(src) / n_sources	×
314
315	astro_dict = rh.nu_astronomy(ul, e_pdf_dict)	×
316
317	key = "Energy Flux (GeV cm^{-2} s^{-1})"	×
318
319	fluence_uls.append(astro_dict[key] * inj_time)	×
320
321	# this functionality has been retired
322	# e_per_source_uls.append(
323	# astro_dict["Mean Luminosity (erg/s)"] * inj_time
324	# )
325
326	energy_flux.append(astro_dict[key])	×
327	x_axis.append(float(subdir))	×
328
329	plt.figure()	×
330	ax = plt.subplot(111)	×
331	plt.plot(x_axis, flux_uls, marker="o", label="upper limit")	×
332
333	if self.mock_unblind:	×
334	ax.text(	×
335	0.2,
336	0.5,
337	"MOCK DATA",
338	color="grey",
339	alpha=0.5,
340	transform=ax.transAxes,
341	)
342
343	ax.set_xlabel("Gamma")	×
344	ax.set_ylabel("Flux")	×
345	plt.yscale("log")	×
346	plt.savefig(self.plot_dir + "upper_limit_flux.pdf")	×
347	plt.close()	×
348
349	plt.figure()	×
350	ax1 = plt.subplot(111)	×
351	# ax2 = ax1.twinx()
352
353	ax1.plot(x_axis, fluence_uls, marker="o")	×
354	# ax2.plot(x_axis, e_per_source_uls, marker="o")
355
356	if self.mock_unblind:	×
357	ax1.text(	×
358	0.2,
359	0.5,
360	"MOCK DATA",
361	color="grey",
362	alpha=0.5,
363	transform=ax1.transAxes,
364	)
365
366	ax1.grid(True, which="both")	×
367	# ax2.grid(True, which="both")
368	ax1.set_xlabel("Gamma")	×
369	# ax2.set_xlabel("Gamma")
370	ax1.set_ylabel(r"Total Fluence [GeV cm$^{-2}$ s$^{-1}$]")	×
371	# ax2.set_ylabel(r"Isotropic-Equivalent Luminosity $L_{\nu}$ (erg/s)")
372	ax1.set_yscale("log")	×
373	# ax2.set_yscale("log")
374
375	axes = [ax1] # was: axes = [ax1, ax2]	×
376	ordinates = [fluence_uls] # was: [fluence_uls, e_per_source_uls]	×
377
378	for k, ax in enumerate(axes):	×
379	y = ordinates[k]	×
380	ax.set_ylim(0.95 * min(y), 1.1 * max(y))	×
381
382	plt.tight_layout()	×
383	plt.savefig(self.plot_dir + "upper_limit_fluence.pdf")	×
384	plt.close()	×
385
386	try:	×
387	os.makedirs(os.path.dirname(self.limit_path))	×
388	except OSError:	×
389	pass	×
390	logger.info("Saving limits to {0}".format(self.limit_path))	×
391
392	res_dict = {	×
393	"gamma": x_axis,
394	"flux": flux_uls,
395	"energy_flux": energy_flux,
396	"fluence": fluence_uls,
397	# "energy": e_per_source_uls,
398	}
399
400	with open(self.limit_path, "wb") as f:	×
401	pickle.dump(res_dict, f)	×
402
403	except AttributeError as e:	3✔
404	logger.warning(f"Unable to set limits. No TS distributions found.: {e}")	3✔
405
406	except ValueError as e:	×
407	raise OverfluctuationError(	×
408	"Insufficent pseudotrial values above and below threshold"
409	)
410
411	def compare_to_background_TS(self):	3✔
412	logger.info(	3✔
413	"Retrieving Background TS Distribution from {0}".format(self.pickle_dir)
414	)
415
416	try:	3✔
417
418	ts_array = list()	3✔
419
420	for subdir in os.listdir(self.pickle_dir):	3✔
421
422	merged_pkl = os.path.join(	3✔
423	os.path.join(self.pickle_dir, subdir), "merged/0.pkl"
424	)
425
426	logger.info("Loading {0}".format(merged_pkl))	3✔
427
428	try:	3✔
429	with open(merged_pkl, "rb") as mp:	3✔
430	merged_data = pickle.load(mp)	3✔
431	except (IOError, OSError) as e:	×
432	if len(ts_array) > 0:	×
433	logger.warning(e)	×
434	else:
435	raise OSError(e)	×
436
437	ts_array += list(merged_data["TS"])	3✔
438	self.background_median_gamma[subdir] = np.median(merged_data["TS"])	3✔
439
440	ts_array = np.array(ts_array)	3✔
441	self.bkg_ts_array = ts_array	3✔
442
443	self.sigma = plot_background_ts_distribution(	3✔
444	ts_array, self.output_file, self.ts_type, self.ts, self.mock_unblind
445	)
446
447	self.background_median = np.median(ts_array)	3✔
448	self.raw_pre_trial_pvalue = len(ts_array[ts_array > self.ts]) / len(	3✔
449	ts_array
450	)
451
452	except (IOError, OSError):	×
453
454	try:	×
455	ts_array = list()	×
456
457	merged_pkl = os.path.join(self.pickle_dir, "merged/0.pkl")	×
458
459	logger.info("No subfolders found, loading {0}".format(merged_pkl))	×
460
461	with open(merged_pkl, "rb") as mp:	×
462	merged_data = pickle.load(mp)	×
463
464	ts_array += list(merged_data["TS"])	×
465
466	ts_array = np.array(ts_array)	×
467	self.bkg_ts_array = ts_array	×
468
469	self.sigma = plot_background_ts_distribution(	×
470	ts_array,
471	self.output_file,
472	self.ts_type,
473	self.ts,
474	self.mock_unblind,
475	)
476
477	self.raw_pre_trial_pvalue = len(ts_array[ts_array > self.ts]) / len(	×
478	ts_array
479	)
480	self.background_median = np.median(ts_array)	×
481
482	except (IOError, OSError):	×
483	logger.warning("No Background TS Distribution found")	×
484	pass	×
485
486	def check_unblind(self):	3✔
487	print("\n")	×
488	print("You are proposing to unblind data.")	×
489	print("\n")	×
490	confirm()	×
491	print("\n")	×
492	print("You are unblinding the following catalogue:")	×
493	print("\n")	×
494	print(self.unblind_dict["catalogue"])	×
495	print("\n")	×
496	confirm()	×
497	print("\n")	×
498	print("The catalogue has the following entries:")	×
499	print("\n")	×
500
501	cat = load_catalogue(self.unblind_dict["catalogue"])	×
502
503	print(cat.dtype.names)	×
504	print(len(cat), np.sum(cat["base_weight"]))	×
505	print(cat)	×
506	print("\n")	×
507	confirm()	×
508	print("\n")	×
509	print("The following datasets will be used:")	×
510	print("\n")	×
511	for x in self.unblind_dict["dataset"].values():	×
512	print(x.sample_name, x.season_name)	×
513	print("\n")	×
514	print(x.exp_path)	×
515	print(x.pseudo_mc_path)	×
516	print("\n")	×
517	confirm()	×
518	print("\n")	×
519	print("The following LLH will be used:")	×
520	print("\n")	×
521	for (key, val) in self.unblind_dict["llh_dict"].items():	×
522	print(key, val)	×
523	print("\n")	×
524	confirm()	×
525	print("\n")	×
526	print("Are you really REALLY sure about this?")	×
527	print("You will unblind. This is your final warning.")	×
528	print("\n")	×
529	confirm()	×
530	print("\n")	×
531	print("OK, you asked for it...")	×
532	print("\n")	×
533
534	def dump_results(self, args, *kwargs):	3✔
535	logger.info(f"Using Unblinder, not dumping results!")	3✔
536
537	return Unblinder(unblind_dict, **kwargs)	3✔

icecube / flarestack / 3998833878

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous