Faster overlap calc (#335)

nuztf/base_scanner.py

@@ -734,30 +734,29 @@ def text_summary(self):
             ]
             detection_jds = [x["jd"] for x in detections]
             first_detection = detections[detection_jds.index(min(detection_jds))]
+
             latest = [
                 x
                 for x in res["prv_candidates"] + [res["candidate"]]
                 if "isdiffpos" in x.keys()
             ][-1]
-            try:
-                last_upper_limit = [
-                    x
-                    for x in res["prv_candidates"]
-                    if np.logical_and(
-                        "isdiffpos" in x.keys(), x["jd"] < first_detection["jd"]
-                    )
-                ][-1]
 
-                text += self.candidate_text(
-                    name,
-                    first_detection["jd"],
-                    last_upper_limit["diffmaglim"],
-                    last_upper_limit["jd"],
-                )
-
-            # No pre-detection upper limit
-            except IndexError:
-                text += self.candidate_text(name, first_detection["jd"], None, None)
+            upper_limits = [
+                x for x in res["prv_candidates"] if x["jd"] < first_detection["jd"]
+            ]
+            if len(upper_limits) > 0:
+                last_upper_limit_mag = upper_limits[-1]["diffmaglim"]
+                last_upper_limit_jd = upper_limits[-1]["jd"]
+            else:
+                last_upper_limit_mag = None
+                last_upper_limit_jd = None
+
+            text += self.candidate_text(
+                ztf_id=name,
+                first_detection=first_detection["jd"],
+                lul_lim=last_upper_limit_mag,
+                lul_jd=last_upper_limit_jd,
+            )
 
             ned_z, ned_dist = query_ned_for_z(
                 ra_deg=latest["ra"],
@@ -833,7 +832,7 @@ def __init__(self, data):
         ras = np.ones_like(data["field"]) * np.nan
         decs = np.ones_like(data["field"]) * np.nan
 
-        # Actually load up ra/dec
+        # Actually load up RA/Dec
 
         veto_fields = []
 
@@ -892,6 +891,18 @@ def __init__(self, data):
         self.logger.info(f"Most recent observation found is {obs_times[-1]}")
         self.logger.info("Unpacking observations")
 
+        # re-read the map and subsample it to nside=64 if nside is too large, as overlap calculation for big maps with large nside take forever
+        if self.nside > 256:
+            (
+                self.map_coords,
+                self.pixel_nos,
+                self.nside,
+                self.map_probs,
+                self.data,
+                self.pixel_area,
+                self.key,
+            ) = self.unpack_skymap(output_nside=256)
+
         pix_map = dict()
         pix_obs_times = dict()
 
@@ -1030,6 +1041,17 @@ def calculate_overlap_with_depot_observations(
 
         self.logger.info("Unpacking observations")
 
+        if self.nside > 256:
+            (
+                self.map_coords,
+                self.pixel_nos,
+                self.nside,
+                self.map_probs,
+                self.data,
+                self.pixel_area,
+                self.key,
+            ) = self.unpack_skymap(output_nside=256)
+
         pix_map = dict()
         pix_obs_times = dict()
 

nuztf/neutrino_scanner.py

@@ -211,9 +211,10 @@ def in_contour(self, ra_deg, dec_deg):
 
         return np.logical_and(in_ra, in_dec)
 
-    def unpack_skymap(self, skymap=None):
+    def unpack_skymap(self, skymap=None, output_nside: None | int = None):
         """ """
-        nside = 1024
+        output_nside = 1024
+
         map_coords = []
         pixel_nos = []
 
@@ -226,15 +227,15 @@ def unpack_skymap(self, skymap=None):
 
         nearish_pixels = list(
             hp.query_disc(
-                nside=nside,
+                nside=output_nside,
                 vec=hp.ang2vec(np.pi / 2.0 - center_dec, center_ra),
                 radius=rad,
                 nest=True,
             )
         )
 
         for i in tqdm(nearish_pixels):
-            ra, dec = self.extract_ra_dec(nside, i)
+            ra, dec = self.extract_ra_dec(output_nside, i)
             if self.in_contour(ra, dec):
                 map_coords.append((ra, dec))
                 pixel_nos.append(i)
@@ -244,9 +245,11 @@ def unpack_skymap(self, skymap=None):
 
         key = "PROB"
 
-        data = np.zeros(hp.nside2npix(nside), dtype=np.dtype([(key, float)]))
+        data = np.zeros(hp.nside2npix(output_nside), dtype=np.dtype([(key, float)]))
         data[np.array(pixel_nos)] = map_probs
 
-        pixel_area = hp.nside2pixarea(nside, degrees=True) * float(len(map_coords))
+        pixel_area = hp.nside2pixarea(output_nside, degrees=True) * float(
+            len(map_coords)
+        )
 
-        return map_coords, pixel_nos, nside, map_probs, data, pixel_area, key
+        return map_coords, pixel_nos, output_nside, map_probs, data, pixel_area, key

nuztf/skymap_scanner.py

@@ -41,6 +41,10 @@ def __init__(
         self.logger = logging.getLogger(__name__)
         self.prob_threshold = prob_threshold
         self.n_days = n_days
+        self.event = event
+        self.rev = rev
+        self.prob_threshold = prob_threshold
+        self.output_nside = output_nside
 
         if config:
             self.config = config
@@ -50,10 +54,10 @@ def __init__(
                 self.config = yaml.safe_load(f)
 
         self.skymap = Skymap(
-            event=event,
-            rev=rev,
-            prob_threshold=prob_threshold,
-            output_nside=output_nside,
+            event=self.event,
+            rev=self.rev,
+            prob_threshold=self.prob_threshold,
+            output_nside=self.output_nside,
         )
 
         self.t_min = Time(self.skymap.t_obs, format="isot", scale="utc")
@@ -284,22 +288,22 @@ def remove_variability_line():
         )
 
     def candidate_text(
-        self, name: str, first_detection: float, lul_lim: float, lul_jd: float
+        self, ztf_id: str, first_detection: float, lul_lim: float, lul_jd: float
     ):
         """ """
         try:
             text = (
                 "{0}, first detected {1:.1f} hours after merger, "
                 "was not detected {2:.1f} days prior to a depth of {3:.2f}. ".format(
-                    name,
+                    ztf_id,
                     24.0 * (first_detection - self.t_min.jd),
                     first_detection - lul_jd,
                     lul_lim,
                 )
             )
         except TypeError:
             text = (
-                f"{name} had upper limit problems. PLEASE FILL IN NUMBERS BY HAND!!! "
+                f"{ztf_id} had upper limit problems. PLEASE FILL IN NUMBERS BY HAND!!! "
             )
 
         return text
@@ -397,8 +401,15 @@ def filter_f_history(self, res: dict, t_max_jd=None):
 
         return True
 
-    def unpack_skymap(self):
+    def unpack_skymap(self, output_nside: None | int = None):
         """ """
+        if output_nside is not None:
+            self.skymap = Skymap(
+                event=self.event,
+                rev=self.rev,
+                prob_threshold=self.prob_threshold,
+                output_nside=output_nside,
+            )
 
         nside = hp.npix2nside(len(self.skymap.data[self.skymap.key]))
 

slackbot.py

@@ -45,6 +45,7 @@ def __init__(
         dl_results: bool = False,
         do_gcn: bool = False,
         time_window: int | None = None,
+        prob_threshold: float | None = None,
     ):
         self.channel = channel
         self.ts = ts
@@ -63,13 +64,20 @@ def __init__(
         else:
             self.time_window = time_window
 
+        if prob_threshold is None:
+            self.prob_threshold = 0.95
+
         self.scanner: NeutrinoScanner | SkymapScanner
 
         if self.event_type == "nu":
             self.scanner = NeutrinoScanner(self.name)
         elif self.event_type == "gw":
             try:
-                self.scanner = SkymapScanner(self.name)
+                self.scanner = SkymapScanner(
+                    self.name,
+                    n_days=self.time_window,
+                    prob_threshold=self.prob_threshold,
+                )
             except EventNotFound:
                 self.post(
                     f"The specified LIGO event, {self.name}, was not found on GraceDB. Please check that you entered the correct event name."

slackbot_server.py

@@ -38,6 +38,7 @@ def scan(
         event_type=event_type,
         do_gcn=do_gcn,
         time_window=time_window,
+        prob_threshold=prob_threshold,
     )
 
 
@@ -114,6 +115,7 @@ def message(payload):
                 return
 
             time_window = None
+            prob_threshold = None
             dl_results = True
 
             for i, parameter in enumerate(split_text):
@@ -133,6 +135,15 @@ def message(payload):
                             thread_ts=ts,
                         )
                         return
+                elif parameter in fuzzy_parameters(["prob", "probability", "p"]):
+                    try:
+                        prob_threshold = float(split_text[i + 1])
+                    except ValueError:
+                        wc.chat_postMessage(
+                            channel=channel_id,
+                            text="Error: --prob has to be a float.",
+                            thread_ts=ts,
+                        )
 
             do_scan = True
             display_help = False

tests/test_skymap_scanner.py

@@ -59,7 +59,7 @@ def test_gw_scan(self):
 
         fig, coverage_summary = scanner.plot_coverage()
 
-        true_coverage_summary = "In total, 88.6 % of the contour was observed at least once.\nThis estimate includes 0.0 % of the contour at a galactic latitude <10 deg.\nIn total, 73.8 % of the contour was observed at least twice. \nIn total, 73.8 % of the contour was observed at least twice, and excluding low galactic latitudes.\nThese estimates account for chip gaps."
+        true_coverage_summary = "In total, 88.0 % of the contour was observed at least once.\nThis estimate includes 0.0 % of the contour at a galactic latitude <10 deg.\nIn total, 72.1 % of the contour was observed at least twice. \nIn total, 72.1 % of the contour was observed at least twice, and excluding low galactic latitudes.\nThese estimates account for chip gaps."
 
         self.assertEqual(coverage_summary, true_coverage_summary)