Add additional cuts and update the module to operate on higher resolution

fink_science/rubin/hostless_detection/config.json

@@ -0,0 +1,5 @@
+{
+  "minimum_number_of_alerts": 3,
+  "cutout_timeframe": 30,
+  "cutout_magnitude": 20
+}

fink_science/rubin/hostless_detection/config_base.json

@@ -0,0 +1,12 @@
+{
+  "image_shape": [30, 30],
+  "sigma_clipping_kwargs": {
+    "sigma": 3,
+    "maxiters": 10
+  },
+  "hostless_detection_with_clipping": {
+    "crop_radius": 15,
+    "max_number_of_pixels_clipped": 5,
+    "min_number_of_pixels_clipped": 3
+  }
+}

fink_science/rubin/hostless_detection/pipeline_utils.py

@@ -15,10 +15,17 @@
 """Implementation of the paper: ELEPHANT: ExtragaLactic alErt Pipeline for Hostless AstroNomical Transients https://arxiv.org/abs/2404.18165"""
 
 import io
+from typing import Dict
 
 from astropy.io import fits
 import numpy as np
 
+from fink_science.ztf.hostless_detection.pipeline_utils import (
+    apply_sigma_clipping,
+    _check_hostless_conditions,
+    crop_center_patch,
+)
+
 
 def read_cutout_stamp(fits_bytes: bytes) -> np.ndarray:
     """
@@ -32,3 +39,69 @@ def read_cutout_stamp(fits_bytes: bytes) -> np.ndarray:
     fits_buffer = io.BytesIO(fits_bytes)
     with fits.open(fits_buffer) as hdulist:
         return hdulist[0].data
+
+
+def is_outliers_in_template(clipped_image, number_of_pixels: int = 20):
+    """Checks if is a big host pixels by applying sigma clipping
+
+    The threshold value is decided based on the samples we have seen.
+    But should be updated eventually
+    """
+    num_template_pixels_masked = np.ma.count_masked(clipped_image)
+    print(num_template_pixels_masked)
+    if num_template_pixels_masked > number_of_pixels:
+        return True
+    return False
+
+
+def run_hostless_detection_with_clipped_data(
+    science_stamp: np.ndarray,
+    template_stamp: np.ndarray,
+    configs: Dict,
+    check_outliers_in_template: bool = True,
+) -> bool:
+    """Detects potential hostless candidates
+
+    Notes
+    -----
+    We use sigma clipped stamp images by
+     cropping an image patch from the center of the image.
+    If pixels are rejected in scientific image but not in corresponding
+     template image, such candidates are flagged as potential hostless
+
+    Parameters
+    ----------
+    science_stamp
+       science image
+    template_stamp
+        template image
+    configs
+        detection configs with detection threshold
+    """
+    sigma_clipping_config = configs["sigma_clipping_kwargs"]
+
+    science_clipped = apply_sigma_clipping(science_stamp, sigma_clipping_config)
+    template_clipped = apply_sigma_clipping(template_stamp, sigma_clipping_config)
+
+    if check_outliers_in_template:
+        if is_outliers_in_template(template_clipped):
+            return False
+    detection_config = configs["hostless_detection_with_clipping"]
+
+    is_hostless_candidate = _check_hostless_conditions(
+        science_clipped, template_clipped, detection_config
+    )
+    if is_hostless_candidate:
+        return is_hostless_candidate
+
+    # Check again at half resolution
+    crop_radius = int(detection_config["crop_radius"] / 2)
+    science_stamp = crop_center_patch(science_stamp, crop_radius)
+    template_stamp = crop_center_patch(template_stamp, crop_radius)
+
+    science_clipped = apply_sigma_clipping(science_stamp, sigma_clipping_config)
+    template_clipped = apply_sigma_clipping(template_stamp, sigma_clipping_config)
+    is_hostless_candidate = _check_hostless_conditions(
+        science_clipped, template_clipped, detection_config
+    )
+    return is_hostless_candidate

fink_science/rubin/hostless_detection/processor.py

@@ -32,10 +32,10 @@
 
 # Share configuration with ZTF
 CONFIGS_BASE = load_json(
-    "{}/ztf/hostless_detection/config_base.json".format(os.path.dirname(__file__))
+    "{}/rubin/hostless_detection/config_base.json".format(os.path.dirname(__file__))
 )
 CONFIGS = load_json(
-    "{}/ztf/hostless_detection/config.json".format(os.path.dirname(__file__))
+    "{}/rubin/hostless_detection/config.json".format(os.path.dirname(__file__))
 )
 CONFIGS.update(CONFIGS_BASE)
 
@@ -52,10 +52,16 @@ def run_potential_hostless(
     ssObjectId: pd.Series,
     nDiaSources: pd.Series,
     psfFlux: pd.Series,
+    templateFlux: pd.Series,
+    templateFluxErr: pd.Series,
     midpointMjdTai: pd.Series,
     firstDiaSourceMjdTaiFink: pd.Series,
     simbad_otype: pd.Series,
     gaiadr3_DR3Name: pd.Series,
+    mangrove_2MASS_name: pd.Series,
+    mangrove_HyperLEDA_name: pd.Series,
+    legacydr8_zphot: pd.Series,
+    spicy_class: pd.Series,
 ) -> pd.Series:
     """Runs potential hostless candidate detection for Rubin without any filtering
 
@@ -71,6 +77,10 @@ def run_potential_hostless(
         Number of previous detections
     psfFlux: pd.Series
         Flux in the difference image
+    templateFlux: pd.Series,
+        Flux in template image
+    templateFluxErr: pd.Series
+        Flux error in template image
     midpointMjdTai: pd.Series
         Emission date for the alert
     firstDiaSourceMjdTaiFink: pd.Series
@@ -79,6 +89,15 @@ def run_potential_hostless(
         SIMBAD type. NaN if not existing
     gaiadr3_DR3Name: pd.Series
         Name in Gaia DR3. NaN if not existing
+    mangrove_2MASS_name: pd.Series
+        Name in mangrove 2MASS. NaN if not existing
+    mangrove_HyperLEDA_name: pd.Series
+        Name in mangrove HyperLEDA. NaN if not existing
+    legacydr8_zphot: pd.Series
+        Photo-z estimate from Legacy Surveys DR8 South Photometric Redshifts catalog.
+         NaN if not existing
+    spicy_class
+        closest source from SPICY catalog. Nan if not existing
 
     Notes
     -----
@@ -114,12 +133,18 @@ def run_potential_hostless(
     ...         F.lit(None),
     ...         F.lit(3),
     ...         F.lit(100000000),
+    ...         F.lit(100000000),
+    ...         F.lit(100),
     ...         df["diaSource.midpointMjdTai"],
     ...         df["diaSource.midpointMjdTai"],
     ...         F.lit(None),
+    ...         F.lit(None),
+    ...         F.lit(None),
+    ...         F.lit(None),
+    ...         F.lit(None),
     ...         F.lit(None),))
     >>> df.filter(df.elephant_kstest_no_cuts.kstest_science.isNotNull()).count()
-    3
+    0
 
     # SSO cuts
     >>> df = df.withColumn('elephant_kstest_sso_cuts',
@@ -129,12 +154,18 @@ def run_potential_hostless(
     ...         df["ssSource.ssObjectId"],
     ...         F.lit(3),
     ...         F.lit(100000000),
+    ...         F.lit(100000000),
+    ...         F.lit(100),
     ...         df["diaSource.midpointMjdTai"],
     ...         df["diaSource.midpointMjdTai"],
     ...         F.lit(None),
+    ...         F.lit(None),
+    ...         F.lit(None),
+    ...         F.lit(None),
+    ...         F.lit(None),
     ...         F.lit(None),))
     >>> df.filter(df.elephant_kstest_sso_cuts.kstest_science.isNotNull()).count()
-    1
+    0
 
     # All cuts
     >>> df = df.withColumn('elephant_kstest',
@@ -144,27 +175,69 @@ def run_potential_hostless(
     ...         df["ssSource.ssObjectId"],
     ...         df["diaObject.nDiaSources"],
     ...         df["diaSource.psfFlux"],
+    ...         df["diaSource.templateFlux"],
+    ...         df["diaSource.templateFluxErr"],
     ...         df["diaSource.midpointMjdTai"],
     ...         F.array_min("prvDiaSources.midpointMjdTai"),
     ...         F.lit("star"),
-    ...         F.lit("DR3 toto"),))
+    ...         F.lit("DR3 toto"),
+    ...         F.lit("galaxy"),
+    ...         F.lit("galaxy"),
+    ...         F.lit("galaxy"),
+    ...         F.lit("galaxy"),))
     >>> df.filter(df.elephant_kstest.kstest_science.isNotNull()).count()
     0
     """
-    f_min_point = nDiaSources >= 2  # N + 1
+    f_min_point = nDiaSources >= CONFIGS["minimum_number_of_alerts"] - 1  # N + 1
     # FIXME: put the conversion formula in fink-utils
-    f_bright = (-2.5 * np.log10(psfFlux) + ZP_NJY) < 20
+    f_bright = (-2.5 * np.log10(psfFlux) + ZP_NJY) < CONFIGS["cutout_magnitude"]
     f_not_in_simbad = simbad_otype.apply(lambda val: val in BAD_VALUES or pd.isna(val))
     f_not_in_gaia = gaiadr3_DR3Name.apply(lambda val: val in BAD_VALUES or pd.isna(val))
+    f_not_in_2mass = mangrove_2MASS_name.apply(
+        lambda val: val in BAD_VALUES or pd.isna(val)
+    )
+    f_not_in_hyperlda = mangrove_HyperLEDA_name.apply(
+        lambda val: val in BAD_VALUES or pd.isna(val)
+    )
+    f_not_in_legacydr8 = legacydr8_zphot.apply(
+        lambda val: val in BAD_VALUES or pd.isna(val)
+    )
+    f_not_in_spicy = spicy_class.apply(lambda val: val in BAD_VALUES or pd.isna(val))
+
     f_not_sso = ssObjectId.apply(lambda val: val in BAD_VALUES or pd.isna(val))
-    f_early = (midpointMjdTai - firstDiaSourceMjdTaiFink) < 30
+    f_early = (midpointMjdTai - firstDiaSourceMjdTaiFink) < CONFIGS["cutout_timeframe"]
 
+    # Cuts on psfFlux and templateFlux
+    PSF_FLUX_THRESHOLD = 0
+    TEMPLATE_FLUX_THRESHOLD = [300, 2000]
+    psfFlux_cut = psfFlux > PSF_FLUX_THRESHOLD
+    templateFlux_cut = (templateFlux > TEMPLATE_FLUX_THRESHOLD[0]) & (
+        templateFlux < TEMPLATE_FLUX_THRESHOLD[1]
+    )
+    # SNR cut > 5
+    template_snr_cut = (templateFlux / templateFluxErr) > 5
     good_candidate = (
-        f_min_point & f_bright & f_not_in_simbad & f_not_in_gaia & f_not_sso & f_early
+        f_min_point
+        & f_bright
+        & f_not_in_simbad
+        & f_not_in_gaia
+        & f_not_sso
+        & f_early
+        & psfFlux_cut
+        & templateFlux_cut
+        & template_snr_cut
+        & f_not_in_2mass
+        & f_not_in_hyperlda
+        & f_not_in_legacydr8
+        & f_not_in_spicy
     )
-
-    default_result = {"kstest_science": None, "kstest_template": None}
+    # Process full image
+    default_result = {
+        "kstest_science": None,
+        "kstest_template": None,
+    }
     kstest_results = []
+
     hostless_science_class = HostLessExtragalacticRubin(CONFIGS_BASE)
     for index in range(cutoutScience.shape[0]):
         if good_candidate[index]:

fink_science/rubin/hostless_detection/run_pipeline.py

@@ -14,14 +14,17 @@
 # limitations under the License.
 """Implementation of the paper: ELEPHANT: ExtragaLactic alErt Pipeline for Hostless AstroNomical Transients https://arxiv.org/abs/2404.18165"""
 
-from typing import Tuple
+from typing import Union, Any
 
 from fink_science.ztf.hostless_detection.pipeline_utils import (
-    run_hostless_detection_with_clipped_data,
     run_powerspectrum_analysis,
+    crop_center_patch,
 )
 from fink_science.ztf.hostless_detection.run_pipeline import HostLessExtragalactic
-from fink_science.rubin.hostless_detection.pipeline_utils import read_cutout_stamp
+from fink_science.rubin.hostless_detection.pipeline_utils import (
+    read_cutout_stamp,
+    run_hostless_detection_with_clipped_data,
+)
 
 
 class HostLessExtragalacticRubin(HostLessExtragalactic):
@@ -34,7 +37,7 @@ class HostLessExtragalacticRubin(HostLessExtragalactic):
 
     def process_candidate_fink_rubin(
         self, science_stamp: bytes, template_stamp: bytes
-    ) -> Tuple[float, float]:
+    ) -> Union[tuple[None, None], tuple[Any, Any]]:
         """
         Processes each candidate
 
@@ -46,27 +49,37 @@ def process_candidate_fink_rubin(
            template stamp data
         """
         science_stamp = read_cutout_stamp(science_stamp)
-
         template_stamp = read_cutout_stamp(template_stamp)
+        crop_radius = self.configs["hostless_detection_with_clipping"]["crop_radius"]
+        science_stamp = crop_center_patch(science_stamp, crop_radius)
+        template_stamp = crop_center_patch(template_stamp, crop_radius)
 
         if science_stamp.shape != template_stamp.shape:
             return None, None
 
-        science_stamp_clipped, template_stamp_clipped = self._run_sigma_clipping(
-            science_stamp, template_stamp
-        )
+        if (science_stamp.shape[0] < self._image_shape[0]) or (
+            template_stamp.shape[0] < self._image_shape[0]
+        ):
+            return None, None
+
+        science_stamp = crop_center_patch(science_stamp, crop_radius)
+        template_stamp = crop_center_patch(template_stamp, crop_radius)
+
         is_hostless_candidate = run_hostless_detection_with_clipped_data(
-            science_stamp_clipped, template_stamp_clipped, self.configs
+            science_stamp, template_stamp, self.configs
         )
         if is_hostless_candidate:
+            science_stamp_clipped, template_stamp_clipped = self._run_sigma_clipping(
+                science_stamp, template_stamp
+            )
             power_spectrum_results = run_powerspectrum_analysis(
                 science_stamp,
                 template_stamp,
                 science_stamp_clipped.mask.astype(int),
                 template_stamp_clipped.mask.astype(int),
-                self._image_shape,
+                crop_radius=crop_radius,
             )
             return power_spectrum_results[
-                "kstest_SCIENCE_15_statistic"
-            ], power_spectrum_results["kstest_TEMPLATE_15_statistic"]
+                "kstest_SCIENCE_statistic"
+            ], power_spectrum_results["kstest_TEMPLATE_statistic"]
         return None, None

fink_science/ztf/hostless_detection/pipeline_utils.py

@@ -150,12 +150,14 @@ def run_hostless_detection_with_clipped_data(
 
     science_clipped = apply_sigma_clipping(science_stamp, sigma_clipping_config)
     template_clipped = apply_sigma_clipping(template_stamp, sigma_clipping_config)
+
     detection_config = configs["hostless_detection_with_clipping"]
     is_hostless_candidate = _check_hostless_conditions(
         science_clipped, template_clipped, detection_config
     )
     if is_hostless_candidate:
         return is_hostless_candidate
+
     science_stamp = crop_center_patch(science_stamp, detection_config["crop_radius"])
     template_stamp = crop_center_patch(template_stamp, detection_config["crop_radius"])
     science_clipped = apply_sigma_clipping(science_stamp, sigma_clipping_config)
@@ -195,8 +197,8 @@ def run_powerspectrum_analysis(
     template_image: np.ndarray,
     science_mask: np.ndarray,
     template_mask: np.ndarray,
-    image_size: List,
     number_of_iterations: int = 200,
+    crop_radius: int = 15,
 ) -> Dict:
     """Runs powerspectrum analysis
 
@@ -232,5 +234,6 @@ def run_powerspectrum_analysis(
         template_data,
         number_of_iterations=number_of_iterations,
         metric="kstest",
+        cutout_size=crop_radius * 2,
     )
     return kstest_results_dict