Repository: cta-observatory/ctapipe
Pull request #2604 (2024-07-26 11:19:28)

OutlierDetector components
This PR added outlier detection approaches. The sub classes can be applied to any astropy column of shape (n_entries, n_channels, n_pix) returning an outlier mask of the same shape. The code utilizing this components later need to handle multiple instances and therefore combine the returning outlier masks.

–––
Related: #2542

ctao-dpps-sonarqube[bot]: ![Passed](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/checks/QualityGateBadge/passed.svg?sanitize=true)

# Analysis Details
## 0 Issues
- ![Bug](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/common/bug.svg?sanitize=true) 0 Bugs
- ![Vulnerability](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/common/vulnerability.svg?sanitize=true) 0 Vulnerabilities
- ![Code Smell](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/common/code_smell.svg?sanitize=true) 0 Code Smells

## Coverage and Duplications
- ![Coverage](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/checks/CoverageChart/90.svg?sanitize=true) 98.60% Coverage (94.20% Estimated after merge)
- ![Duplications](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/checks/Duplications/3.svg?sanitize=true) 0.00% Duplicated Code (0.60% Estimated after merge)

**Project ID:** cta-observatory_ctapipe_AY52EYhuvuGcMFidNyUs

[View in SonarQube](https://sonar-cta-dpps.zeuthen.desy.de/dashboard?id=cta-observatory_ctapipe_AY52EYhuvuGcMFidNyUs&pullRequest=2604)


ctao-dpps-sonarqube[bot]: ![Passed](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/checks/QualityGateBadge/passed.svg?sanitize=true)

# Analysis Details
## 0 Issues
- ![Bug](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/common/bug.svg?sanitize=true) 0 Bugs
- ![Vulnerability](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/common/vulnerability.svg?sanitize=true) 0 Vulnerabilities
- ![Code Smell](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/common/code_smell.svg?sanitize=true) 0 Code Smells

## Coverage and Duplications
- ![Coverage](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/checks/CoverageChart/90.svg?sanitize=true) 98.60% Coverage (93.60% Estimated after merge)
- ![Duplications](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/checks/Duplications/3.svg?sanitize=true) 0.00% Duplicated Code (0.60% Estimated after merge)

**Project ID:** cta-observatory_ctapipe_AY52EYhuvuGcMFidNyUs

[View in SonarQube](https://sonar-cta-dpps.zeuthen.desy.de/dashboard?id=cta-observatory_ctapipe_AY52EYhuvuGcMFidNyUs&pullRequest=2604)


TjarkMiener: This Pr needs to be rebased and "outlier files" have to moved to the `monitoring` module once #2554 is merged.

Patch:
diff --git a/docs/api-reference/monitoring/index.rst b/docs/api-reference/monitoring/index.rst
index 9d1fbc38..51268d35 100644
--- a/docs/api-reference/monitoring/index.rst
+++ b/docs/api-reference/monitoring/index.rst
@@ -10,7 +10,7 @@ Monitoring data are time-series used to monitor the status or quality of hardwar
 
 This module provides some code to help to generate monitoring data from processed event data, particularly for the purposes of calibration and data quality assessment.
 
-Currently, only code related to :ref:`stats_aggregator` is implemented here.
+Currently, only code related to :ref:`stats_aggregator` and :ref:`outlier_detector` is implemented here.
 
 
 Submodules
@@ -21,6 +21,7 @@ Submodules
   :glob:
 
   aggregator
+  outlier
 
 
 Reference/API
diff --git a/docs/api-reference/monitoring/outlier.rst b/docs/api-reference/monitoring/outlier.rst
new file mode 100644
index 00000000..e4d02af3
--- /dev/null
+++ b/docs/api-reference/monitoring/outlier.rst
@@ -0,0 +1,11 @@
+.. _outlier_detector:
+
+****************
+Outlier Detector
+****************
+
+
+Reference/API
+=============
+
+.. automodapi:: ctapipe.monitoring.outlier
diff --git a/docs/changes/2604.feature.rst b/docs/changes/2604.feature.rst
new file mode 100644
index 00000000..751bd591
--- /dev/null
+++ b/docs/changes/2604.feature.rst
@@ -0,0 +1,1 @@
+Add outlier detection components to identify faulty pixels.
diff --git a/src/ctapipe/monitoring/outlier.py b/src/ctapipe/monitoring/outlier.py
new file mode 100644
index 00000000..71bcfd36
--- /dev/null
+++ b/src/ctapipe/monitoring/outlier.py
@@ -0,0 +1,138 @@
+"""
+Outlier detection algorithms to identify faulty pixels
+"""
+
+__all__ = [
+    "OutlierDetector",
+    "RangeOutlierDetector",
+    "MedianOutlierDetector",
+    "StdOutlierDetector",
+]
+
+from abc import abstractmethod
+
+import numpy as np
+
+from ctapipe.core import TelescopeComponent
+from ctapipe.core.traits import Float, List
+
+
+class OutlierDetector(TelescopeComponent):
+    """
+    Base class for outlier detection algorithms.
+    """
+
+    @abstractmethod
+    def __call__(self, column) -> bool:
+        """
+        Detect outliers in the provided column.
+
+        This function should be implemented by subclasses to define the specific
+        outlier detection approach. The function examines the statistics in the
+        given column of the table and returns a boolean mask indicating which
+        entries are considered as outliers.
+
+        Parameters
+        ----------
+        column : astropy.table.Column
+            column with chunk-wise aggregated statistic values (mean, median, or std)
+            of shape (n_entries, n_channels, n_pix)
+
+        Returns
+        -------
+        boolean mask
+            mask of outliers of shape (n_entries, n_channels, n_pix)
+        """
+        pass
+
+
+class RangeOutlierDetector(OutlierDetector):
+    """
+    Detect outliers based on a valid range.
+
+    The clipping interval to set the thresholds for detecting outliers corresponds to
+    a configurable range of valid statistic values.
+    """
+
+    validity_range = List(
+        trait=Float(),
+        default_value=[1.0, 2.0],
+        help=(
+            "Defines the range of acceptable values (lower, upper) in units of image values. "
+            "Values outside this range will be flagged as outliers."
+        ),
+        minlen=2,
+        maxlen=2,
+    ).tag(config=True)
+
+    def __call__(self, column):
+        # Remove outliers is statistical values out a given range
+        outliers = np.logical_or(
+            column < self.validity_range[0],
+            column > self.validity_range[1],
+        )
+        return outliers
+
+
+class MedianOutlierDetector(OutlierDetector):
+    """
+    Detect outliers based on the deviation from the camera median.
+
+    The clipping interval to set the thresholds for detecting outliers is computed by multiplying
+    the configurable factors and the camera median of the statistic values.
+    """
+
+    median_range_factors = List(
+        trait=Float(),
+        default_value=[-1.0, 1.0],
+        help=(
+            "Defines the range of acceptable values (lower, upper) as fractional deviations from the mean. "
+            "Values outside of this range will be flagged as outliers."
+        ),
+        minlen=2,
+        maxlen=2,
+    ).tag(config=True)
+
+    def __call__(self, column):
+        # Camera median
+        camera_median = np.ma.median(column, axis=2)
+        # Detect outliers based on the deviation of the median distribution
+        deviation = column - camera_median[:, :, np.newaxis]
+        outliers = np.logical_or(
+            deviation < self.median_range_factors[0] * camera_median[:, :, np.newaxis],
+            deviation > self.median_range_factors[1] * camera_median[:, :, np.newaxis],
+        )
+        return outliers
+
+
+class StdOutlierDetector(OutlierDetector):
+    """
+    Detect outliers based on the deviation from the camera standard deviation.
+
+    The clipping interval to set the thresholds for detecting outliers is computed by multiplying
+    the configurable factors and the camera standard deviation of the statistic values.
+    """
+
+    std_range_factors = List(
+        trait=Float(),
+        default_value=[-1.0, 1.0],
+        help=(
+            "Defines the range of acceptable values (lower, upper) in units of standard deviations. "
+            "Values outside of this range will be flagged as outliers."
+        ),
+        minlen=2,
+        maxlen=2,
+    ).tag(config=True)
+
+    def __call__(self, column):
+        # Camera median
+        camera_median = np.ma.median(column, axis=2)
+        # Camera std
+        camera_std = np.ma.std(column, axis=2)
+        # Detect outliers based on the deviation of the standard deviation distribution
+        deviation = column - camera_median[:, :, np.newaxis]
+        outliers = np.logical_or(
+            deviation < self.std_range_factors[0] * camera_std[:, :, np.newaxis],
+            deviation > self.std_range_factors[1] * camera_std[:, :, np.newaxis],
+        )
+        return outliers


Test patch:
diff --git a/src/ctapipe/monitoring/tests/test_outlier.py b/src/ctapipe/monitoring/tests/test_outlier.py
new file mode 100644
index 00000000..da7d7619
--- /dev/null
+++ b/src/ctapipe/monitoring/tests/test_outlier.py
@@ -0,0 +1,104 @@
+"""
+Tests for OutlierDetector and related functions
+"""
+
+import numpy as np
+from astropy.table import Table
+
+from ctapipe.monitoring.outlier import (
+    MedianOutlierDetector,
+    RangeOutlierDetector,
+    StdOutlierDetector,
+)
+
+
+def test_range_detection(example_subarray):
+    """test the RangeOutlierDetector"""
+
+    # Create dummy data for testing
+    rng = np.random.default_rng(0)
+    # Distribution mimics the mean values of peak times of flat-field events
+    mean = rng.normal(18.0, 0.4, size=(50, 2, 1855))
+    # Insert outliers
+    mean[12, 0, 120] = 3.0
+    mean[42, 1, 67] = 35.0
+    # Create astropy table
+    table = Table([mean], names=("mean",))
+    # Initialize the outlier detector based on the range of valid values
+    # In this test, the interval [15, 25] corresponds to the range (in waveform samples)
+    # of accepted mean (or median) values of peak times of flat-field events
+    detector = RangeOutlierDetector(
+        subarray=example_subarray, validity_range=[15.0, 25.0]
+    )
+    # Detect outliers
+    outliers = detector(table["mean"])
+    # Construct the expected outliers
+    expected_outliers = np.zeros((50, 2, 1855), dtype=bool)
+    expected_outliers[12, 0, 120] = True
+    expected_outliers[42, 1, 67] = True
+    # Check if outliers where detected correctly
+    np.testing.assert_array_equal(outliers, expected_outliers)
+
+
+def test_median_detection(example_subarray):
+    """test the MedianOutlierDetector"""
+
+    # Create dummy data for testing
+    rng = np.random.default_rng(0)
+    # Distribution mimics the median values of charge images of flat-field events
+    median = rng.normal(77.0, 0.6, size=(50, 2, 1855))
+    # Insert outliers
+    median[12, 0, 120] = 1.2
+    median[42, 1, 67] = 21045.1
+    # Create astropy table
+    table = Table([median], names=("median",))
+    # Initialize the outlier detector based on the deviation from the camera median
+    # In this test, the interval [-0.9, 8] corresponds to multiplication factors
+    # typical used for the median values of charge images of flat-field events
+    detector = MedianOutlierDetector(
+        subarray=example_subarray, median_range_factors=[-0.9, 8.0]
+    )
+    # Detect outliers
+    outliers = detector(table["median"])
+    # Construct the expected outliers
+    expected_outliers = np.zeros((50, 2, 1855), dtype=bool)
+    expected_outliers[12, 0, 120] = True
+    expected_outliers[42, 1, 67] = True
+    # Check if outliers where detected correctly
+    np.testing.assert_array_equal(outliers, expected_outliers)
+
+
+def test_std_detection(example_subarray):
+    """test the StdOutlierDetector"""
+
+    # Create dummy data for testing
+    rng = np.random.default_rng(0)
+    # Distribution mimics the std values of charge images of flat-field events
+    ff_std = rng.normal(10.0, 2.0, size=(50, 2, 1855))
+    # Distribution mimics the median values of charge images of pedestal events
+    ped_median = rng.normal(2.0, 1.5, size=(50, 2, 1855))
+    # Insert outliers
+    ff_std[12, 0, 120] = 45.5
+    ped_median[42, 1, 67] = 77.2
+    # Create astropy table
+    ff_table = Table([ff_std], names=("std",))
+    ped_table = Table([ped_median], names=("median",))
+
+    # Initialize the outlier detector based on the deviation from the camera standard deviation
+    # In this test, the interval [-15, 15] corresponds to multiplication factors
+    # typical used for the std values of charge images of flat-field events
+    # and median (and std) values of charge images of pedestal events
+    detector = StdOutlierDetector(
+        subarray=example_subarray, std_range_factors=[-15.0, 15.0]
+    )
+    ff_outliers = detector(ff_table["std"])
+    ped_outliers = detector(ped_table["median"])
+    # Construct the expected outliers
+    ff_expected_outliers = np.zeros((50, 2, 1855), dtype=bool)
+    ff_expected_outliers[12, 0, 120] = True
+    ped_expected_outliers = np.zeros((50, 2, 1855), dtype=bool)
+    ped_expected_outliers[42, 1, 67] = True
+
+    # Check if outliers where detected correctly
+    np.testing.assert_array_equal(ff_outliers, ff_expected_outliers)
+    np.testing.assert_array_equal(ped_outliers, ped_expected_outliers)