A time clustering algorithm for image cleaning

.github/workflows/ci.yml

@@ -75,28 +75,28 @@ jobs:
           key: ctapipe-test-data
 
       - name: Prepare mamba installation
-        if: matrix.install-method == 'mamba'
+        if: matrix.install-method == 'mamba' &&  contains(github.event.pull_request.labels.*.name, 'documentation-only') == false
         env:
           PYTHON_VERSION: ${{ matrix.python-version }}
         run: |
           # setup correct python version
           sed -i -e "s/- python=.*/- python=$PYTHON_VERSION/g" environment.yml
 
       - name: mamba setup
-        if: matrix.install-method == 'mamba'
+        if: matrix.install-method == 'mamba' && contains(github.event.pull_request.labels.*.name, 'documentation-only') == false
         uses: mamba-org/setup-micromamba@v1
         with:
           environment-file: environment.yml
           cache-downloads: true
 
       - name: Python setup
-        if: matrix.install-method == 'pip'
+        if: matrix.install-method == 'pip' && contains(github.event.pull_request.labels.*.name, 'documentation-only') == false 
         uses: actions/setup-python@v4
         with:
           python-version: ${{ matrix.python-version }}
           check-latest: true
 
-      - if: ${{ matrix.install-method == 'pip' && runner.os == 'macOS' }}
+      - if: ${{ matrix.install-method == 'pip' && runner.os == 'macOS' }} && contains(github.event.pull_request.labels.*.name, 'documentation-only') == false
         name: Fix Python PATH on macOS
         # See https://github.com/actions/setup-python/issues/132 and
         # https://github.com/actions/setup-python/issues/132#issuecomment-779406058
@@ -107,6 +107,7 @@ jobs:
           tee -a ~/.bash_profile <<<'export PATH="$pythonLocation/bin:$PATH"'
 
       - name: Install dependencies
+        if: contains(github.event.pull_request.labels.*.name, 'documentation-only') == false
         run: |
           python --version
           pip install pytest-cov restructuredtext-lint pytest-xdist 'coverage!=6.3.0'
@@ -115,14 +116,17 @@ jobs:
           pip freeze
 
       - name: Static codechecks
+        if: contains(github.event.pull_request.labels.*.name, 'documentation-only') == false
         run: |
           restructuredtext-lint README.rst
 
       - name: ctapipe-info
+        if: contains(github.event.pull_request.labels.*.name, 'documentation-only') == false
         run: |
           ctapipe-info --all
 
       - name: Tests
+        if: contains(github.event.pull_request.labels.*.name, 'documentation-only') == false
         run: |
           cd $(mktemp -d)
           pytest -n auto --dist loadscope \

README.rst

@@ -71,11 +71,11 @@ or via::
 
   pip install ctapipe
 
-**Note**: to install a specific version of ctapipe take look at the documentation `here <https://ctapipe.readthedocs.org/en/latest/getting_started_users/>`__.
+**Note**: to install a specific version of ctapipe take look at the documentation `here <https://ctapipe.readthedocs.io/en/latest/user-guide/index.html>`__.
 
 **Note**: ``mamba`` is a C++ reimplementation of conda and can be found `here <https://github.com/mamba-org/mamba>`__.
 
 Note this is *pre-alpha* software and is not yet stable enough for end-users (expect large API changes until the first stable 1.0 release).
 
 Developers should follow the development install instructions found in the
-`documentation <https://ctapipe.readthedocs.org/en/latest/getting_started/>`__.
+`documentation <https://ctapipe.readthedocs.io/en/latest/developer-guide/getting-started.html>`__.

ctapipe/containers.py

@@ -772,6 +772,7 @@ class SimulationConfigContainer(Container):
     Configuration parameters of the simulation
     """
 
+    run_number = Field(np.int32(-1), description="Original sim_telarray run number")
     corsika_version = Field(nan, description="CORSIKA version * 1000")
     simtel_version = Field(nan, description="sim_telarray version * 1000")
     energy_range_min = Field(
@@ -920,12 +921,17 @@ class ReconstructedGeometryContainer(Container):
         "reconstructed core position uncertainty along tilted frame Y axis",
         unit=u.m,
     )
-    h_max = Field(nan * u.m, "reconstructed height of the shower maximum", unit=u.m)
+    h_max = Field(
+        nan * u.m,
+        "reconstructed vertical height above sea level of the shower maximum",
+        unit=u.m,
+    )
     h_max_uncert = Field(nan * u.m, "uncertainty of h_max", unit=u.m)
+
     is_valid = Field(
         False,
         (
-            "direction validity flag. True if the shower direction"
+            "Geometry validity flag. True if the shower geometry"
             "was properly reconstructed by the algorithm"
         ),
     )

ctapipe/core/telescope_component.py

@@ -217,8 +217,10 @@ def __getitem__(self, tel: Optional[Union[int, str]]):
             try:
                 return self._value_for_tel_id[tel]
             except KeyError:
+                if tel not in self._subarray.tel:
+                    raise KeyError(f"No telescope with id {tel} in subarray")
                 raise KeyError(
-                    f"TelescopeParameterLookup: no "
+                    "TelescopeParameterLookup: no "
                     f"parameter value was set for telescope with tel_id="
                     f"{tel}. Please set it explicitly, "
                     f"or by telescope type or '*'."

ctapipe/core/tests/test_telescope_component.py

@@ -342,3 +342,16 @@ class Foo(TelescopeComponent):
 
     f = Foo(mock_subarray, bar=[("type", "*", 1), ("id", 1, None)])
     assert f.bar.tel[1] is None
+
+
+def test_telescope_parameter_nonexistent_telescope(mock_subarray):
+    class Foo(TelescopeComponent):
+        bar = IntTelescopeParameter(
+            default_value=None,
+            allow_none=True,
+        ).tag(config=True)
+
+    foo = Foo(subarray=mock_subarray)
+
+    with pytest.raises(KeyError, match="No telescope with id 0"):
+        foo.bar.tel[0]

ctapipe/image/cleaning.py

@@ -15,6 +15,7 @@
 __all__ = [
     "tailcuts_clean",
     "dilate",
+    "time_clustering",
     "mars_cleaning_1st_pass",
     "fact_image_cleaning",
     "apply_time_delta_cleaning",
@@ -27,6 +28,7 @@
 from abc import abstractmethod
 
 import numpy as np
+from sklearn.cluster import DBSCAN
 
 from ..core import TelescopeComponent
 from ..core.traits import (
@@ -111,6 +113,71 @@ def tailcuts_clean(
         )
 
 
+def time_clustering(
+    geom,
+    image,
+    time,
+    minpts=5,
+    eps=1.0,
+    time_scale_ns=4.0,
+    space_scale_m=0.25,
+    hard_cut_pe=4,
+):
+    """
+    Clean an image by selecting pixels which pass a time clustering algorithm using DBSCAN.
+    Previously used for HESS [timecleaning]_.
+
+    As a neighbor-based image extractor algorithm can lead to biases in the time reconstruction of noise pixels,
+    specially those next to the shower, a cut in the minimum signal image is applied.
+
+    DBSCAN runs with the reconstructed times and pixel positions after scaling. Scaling is needed because eps
+    is not dimension dependent. If scaling is performed properly, eps can be set to 1. DBSCAN returns the
+    cluster IDs of each point. Pixels associated to cluster ID -1 are classified as noise.
+
+    Parameters
+    ----------
+    geom: `ctapipe.instrument.CameraGeometry`
+        Camera geometry information
+    image: array
+        pixel charge information
+    time: array
+        pixel timing information
+    minpts: int
+        Minimum number of points to consider a cluster
+    eps: float
+        Minimum distance in dbscan
+    time_scale_ns: float
+        Time scale in ns
+    space_scale: float
+        Space scale in m
+    hard_cut_pe: float
+        Hard cut in the number of signal pe
+
+    Returns
+    -------
+    A boolean mask of *clean* pixels.
+    """
+    precut_mask = image > hard_cut_pe
+
+    arr = np.zeros(len(image), dtype=float)
+    arr[~precut_mask] = -1
+
+    pix_x = geom.pix_x.value[precut_mask] / space_scale_m
+    pix_y = geom.pix_y.value[precut_mask] / space_scale_m
+
+    X = np.column_stack((time[precut_mask] / time_scale_ns, pix_x, pix_y))
+
+    labels = DBSCAN(eps=eps, min_samples=minpts).fit_predict(X)
+
+    # no_clusters = len(np.unique(labels))-1  # Could be used for gh separation
+
+    y = np.array(arr[(arr == 0)])
+    y[(labels == -1)] = -1
+    arr[arr == 0] = y
+    mask = arr == 0  # we keep these events
+    return mask
+
+
 def mars_cleaning_1st_pass(
     geom,
     image,
@@ -532,6 +599,43 @@ def __call__(
         )
 
 
+class TimeCleaner(ImageCleaner):
+    """
+    Clean images using the time clustering cleaning method
+    """
+
+    space_scale_m = FloatTelescopeParameter(
+        default_value=0.25, help="Pixel space scaling parameter in m"
+    ).tag(config=True)
+    time_scale_ns = FloatTelescopeParameter(
+        default_value=4.0, help="Time scale parameter in ns"
+    ).tag(config=True)
+    minpts = IntTelescopeParameter(
+        default_value=5, help="minimum number of points to form a cluster"
+    ).tag(config=True)
+    eps = FloatTelescopeParameter(
+        default_value=1.0, help="minimum distance in DBSCAN"
+    ).tag(config=True)
+    hard_cut_pe = FloatTelescopeParameter(
+        default_value=2.5, help="Hard cut in the number of pe"
+    ).tag(config=True)
+
+    def __call__(
+        self, tel_id: int, image: np.ndarray, arrival_times=None
+    ) -> np.ndarray:
+        """Apply HESS image cleaning. see ImageCleaner.__call__()"""
+        return time_clustering(
+            geom=self.subarray.tel[tel_id].camera.geometry,
+            image=image,
+            time=arrival_times,
+            eps=self.eps.tel[tel_id],
+            space_scale_m=self.space_scale_m.tel[tel_id],
+            time_scale_ns=self.time_scale_ns.tel[tel_id],
+            minpts=self.minpts.tel[tel_id],
+            hard_cut_pe=self.hard_cut_pe.tel[tel_id],
+        )
+
+
 class MARSImageCleaner(TailcutsImageCleaner):
     """
     1st-pass MARS-like Image cleaner (See `ctapipe.image.mars_cleaning_1st_pass`)

ctapipe/image/tests/test_cleaning.py

@@ -334,6 +334,73 @@ def test_apply_time_delta_cleaning(prod3_lst):
     assert (test_mask == td_mask).all()
 
 
+def test_time_cleaning():
+    geom = CameraGeometry.from_name("LSTCam")
+    charge = np.zeros(geom.n_pixels, dtype=np.float64)
+    peak_time = np.zeros(geom.n_pixels, dtype=np.float64)
+
+    core_pixel = 10
+    core_neighbors = geom.neighbors[core_pixel]
+
+    charge[core_pixel], charge[core_neighbors] = 15, 5
+    peak_time[core_pixel], peak_time[core_neighbors] = (0.0, 0.0)
+
+    mask = cleaning.time_clustering(
+        geom,
+        charge,
+        peak_time,
+        eps=1.0,
+        space_scale_m=0.25,
+        time_scale_ns=4.0,
+        hard_cut_pe=10,
+        minpts=5,
+    )
+
+    assert np.count_nonzero(mask) == 0
+
+    mask = cleaning.time_clustering(
+        geom,
+        charge,
+        peak_time,
+        eps=1.0,
+        space_scale_m=0.25,
+        time_scale_ns=4.0,
+        hard_cut_pe=1,
+        minpts=5,
+    )
+
+    assert np.count_nonzero(mask) == 7
+
+    mask = cleaning.time_clustering(
+        geom,
+        charge,
+        peak_time,
+        eps=1.0,
+        space_scale_m=0.25,
+        time_scale_ns=4.0,
+        hard_cut_pe=1,
+        minpts=8,
+    )
+
+    assert np.count_nonzero(mask) == 0
+
+    charge[core_pixel], charge[core_neighbors] = 15, 5
+    peak_time[core_pixel], peak_time[core_neighbors] = (0.0, 30.0)
+
+    mask = cleaning.time_clustering(
+        geom,
+        charge,
+        peak_time,
+        eps=1.0,
+        space_scale_m=0.25,
+        time_scale_ns=4.0,
+        hard_cut_pe=1,
+        minpts=5,
+    )
+
+    assert np.count_nonzero(mask) == 6
+
+
 def test_time_constrained_clean():
     geom = CameraGeometry.from_name("LSTCam")
     charge = np.zeros(geom.n_pixels, dtype=np.float64)

ctapipe/instrument/camera/geometry.py

@@ -800,7 +800,7 @@ def pixel_moment_matrix(self):
         x = self.pix_x.value
         y = self.pix_y.value
 
-        return np.row_stack(
+        return np.vstack(
             [
                 x,
                 y,
@@ -954,12 +954,13 @@ def position_to_pix_index(self, x, y):
         pix_indices: Pixel index or array of pixel indices. Returns -1 if position falls
                     outside camera
         """
-
         if not self._all_pixel_areas_equal:
             logger.warning(
                 " Method not implemented for cameras with varying pixel sizes"
             )
         unit = x.unit
+        scalar = x.ndim == 0
+
         points_searched = np.dstack([x.to_value(unit), y.to_value(unit)])
         circum_rad = self._pixel_circumradius[0].to_value(unit)
         kdtree = self._kdtree
@@ -969,8 +970,9 @@ def position_to_pix_index(self, x, y):
         del dist
         pix_indices = pix_indices.flatten()
 
+        invalid = np.iinfo(pix_indices.dtype).min
         # 1. Mark all points outside pixel circumeference as lying outside camera
-        pix_indices[pix_indices == self.n_pixels] = -1
+        pix_indices[pix_indices == self.n_pixels] = invalid
 
         # 2. Accurate check for the remaing cases (within circumference, but still outside
         # camera). It is first checked if any border pixel numbers are returned.
@@ -1006,17 +1008,9 @@ def position_to_pix_index(self, x, y):
                 )
                 del dist_check
                 if index_check != insidepix_index:
-                    pix_indices[index] = -1
-
-        # print warning:
-        for index in np.where(pix_indices == -1)[0]:
-            logger.warning(
-                " Coordinate ({} m, {} m) lies outside camera".format(
-                    points_searched[0][index, 0], points_searched[0][index, 1]
-                )
-            )
+                    pix_indices[index] = invalid
 
-        return pix_indices if len(pix_indices) > 1 else pix_indices[0]
+        return np.squeeze(pix_indices) if scalar else pix_indices
 
     @staticmethod
     def simtel_shape_to_type(pixel_shape):