Repository: cta-observatory/ctapipe
Pull request #1993 (2022-07-04 13:03:32)

Add flake8 pre-commit hook. Change black hook to also format notebooks.
I like having my code checked by pyflakes/flake8, but I dislike getting told so by our CI.

Personally, I think it is better to have this checked locally, to prevent "fix pyflakes"-commits.

I chose flake8, because it wraps pyflakes among others

--

`black-jupyter` formats python files *and* notebooks


Patch:
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index 2e19564e..3995f024 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -1,5 +1,17 @@
 repos:
+  # https://pycqa.github.io/isort/docs/configuration/black_compatibility.html#integration-with-pre-commit
+  - repo: https://github.com/pycqa/isort
+    rev: 5.6.4
+    hooks:
+      - id: isort
+        args: ["--profile", "black", "--filter-files"]
   - repo: https://github.com/psf/black
-    rev: 22.3.0
+    rev: 22.6.0
+    hooks:
+      - id: black-jupyter
+  # https://black.readthedocs.io/en/stable/guides/using_black_with_other_tools.html?highlight=other%20tools#flake8
+  - repo: https://github.com/PyCQA/flake8
+    rev: 4.0.1
     hooks:
-      - id: black
\ No newline at end of file
+      - id: flake8
+        args: [--max-line-length=88, --extend-ignore=E203]
diff --git a/ctapipe/image/cleaning.py b/ctapipe/image/cleaning.py
index eeace737..78fdde8e 100644
--- a/ctapipe/image/cleaning.py
+++ b/ctapipe/image/cleaning.py
@@ -18,6 +18,8 @@ __all__ = [
     "mars_cleaning_1st_pass",
     "fact_image_cleaning",
     "apply_time_delta_cleaning",
+    "apply_time_average_cleaning",
+    "time_constrained_clean",
     "ImageCleaner",
     "TailcutsImageCleaner",
 ]
@@ -32,6 +34,7 @@ from ..core.traits import (
     IntTelescopeParameter,
     BoolTelescopeParameter,
 )
+from .morphology import number_of_islands, brightest_island
 
 
 def tailcuts_clean(
@@ -238,6 +241,54 @@ def apply_time_delta_cleaning(
     return pixels_to_keep
 
 
+def apply_time_average_cleaning(
+    geom, mask, image, arrival_times, picture_thresh, time_limit
+):
+    """
+    Extract all pixels that arrived within a given timeframe
+    with respect to the time average of the pixels on the main island.
+
+    In order to avoid removing signal pixels of large impact-parameter events,
+    the time limit for bright pixels is doubled.
+
+    Parameters
+    ----------
+    geom: `ctapipe.instrument.CameraGeometry`
+        Camera geometry information
+    mask: array, boolean
+        boolean mask of *clean* pixels before time_delta_cleaning
+    image: array
+        pixel values
+    arrival_times: array
+        pixel timing information
+    picture_thresh: float
+        threshold above which time limit is extended twice its value
+    time_limit: int or float
+        arrival time limit w.r.t. the average time of the core pixels
+
+    Returns
+    -------
+
+    A boolean mask of *clean* pixels.
+    """
+    mask = mask.copy()
+    if np.count_nonzero(mask) > 0:
+
+        # use main island (maximum charge) for time average calculation
+        num_islands, island_labels = number_of_islands(geom, mask)
+        mask_main = brightest_island(num_islands, island_labels, image)
+        time_ave = np.average(arrival_times[mask_main], weights=image[mask_main] ** 2)
+
+        time_diffs = np.abs(arrival_times[mask] - time_ave)
+        time_limit_pixwise = np.where(
+            image < (2 * picture_thresh), time_limit, time_limit * 2
+        )[mask]
+
+        mask[mask] &= time_diffs < time_limit_pixwise
+
+    return mask
+
+
 def fact_image_cleaning(
     geom,
     image,
@@ -323,6 +374,92 @@ def fact_image_cleaning(
     return pixels_to_keep
 
 
+def time_constrained_clean(
+    geom,
+    image,
+    arrival_times,
+    picture_thresh=7,
+    boundary_thresh=5,
+    time_limit_core=4.5,
+    time_limit_boundary=1.5,
+    min_number_picture_neighbors=1,
+):
+    """
+    time constrained cleaning by MAGIC
+
+    Cleaning contains the following steps:
+    - Find core pixels (containing more photons than a picture threshold)
+    - Remove pixels with less than N neighbors
+    - Remove pixels whose arrival times are within a time limit of the average time
+    - Find boundary pixels (containing more photons than a boundary threshold)
+    - Remove pixels with less than N neighbors arriving within a given timeframe
+
+    Parameters
+    ----------
+    geom: `ctapipe.instrument.CameraGeometry`
+        Camera geometry information
+    image: array
+        pixel values
+    arrival_times: array
+        pixel timing information
+    picture_threshold: float or array
+        threshold above which all pixels are retained
+    boundary_threshold: float or array
+        threshold above which pixels are retained if they have a neighbor
+        already above the picture_thresh
+    time_limit_core: int or float
+        arrival time limit of core pixels w.r.t the average time
+    time_limit_boundary: int or float
+        arrival time limit of boundary pixels w.r.t neighboring core pixels
+    min_number_neighbors: int
+        Threshold to determine if a pixel survives cleaning steps.
+        These steps include checks of neighbor arrival time and value
+
+    Returns
+    -------
+    A boolean mask of *clean* pixels.
+    """
+
+    # find core pixels that pass a picture threshold
+    pixels_above_picture = image >= picture_thresh
+
+    # require at least min_number_picture_neighbors
+    number_of_neighbors_above_picture = geom.neighbor_matrix_sparse.dot(
+        pixels_above_picture.view(np.byte)
+    )
+    pixels_in_picture = pixels_above_picture & (
+        number_of_neighbors_above_picture >= min_number_picture_neighbors
+    )
+
+    # keep core pixels whose arrival times are within a certain time limit of the average
+    mask_core = apply_time_average_cleaning(
+        geom, pixels_in_picture, image, arrival_times, picture_thresh, time_limit_core
+    )
+
+    # find boundary pixels that pass a boundary threshold
+    pixels_above_boundary = image >= boundary_thresh
+    pixels_with_picture_neighbors = geom.neighbor_matrix_sparse.dot(mask_core)
+    mask_boundary = (pixels_above_boundary & pixels_with_picture_neighbors) & np.invert(
+        mask_core
+    )
+
+    # keep boundary pixels whose arrival times are within a certain time limit of the neighboring core pixels
+    mask_boundary = mask_boundary.copy()
+
+    time_diffs = np.abs(arrival_times[mask_boundary, None] - arrival_times)
+    valid_neighbors = (
+        (time_diffs < time_limit_boundary)
+        & geom.neighbor_matrix[mask_boundary]
+        & mask_core
+    )
+    enough_neighbors = (
+        np.count_nonzero(valid_neighbors, axis=1) >= min_number_picture_neighbors
+    )
+    mask_boundary[mask_boundary] &= enough_neighbors
+
+    return mask_core | mask_boundary
+
+
 class ImageCleaner(TelescopeComponent):
     """
     Abstract class for all configurable Image Cleaning algorithms.   Use
@@ -430,7 +567,7 @@ class FACTImageCleaner(TailcutsImageCleaner):
     def __call__(
         self, tel_id: int, image: np.ndarray, arrival_times=None
     ) -> np.ndarray:
-        """ Apply FACT-style image cleaning. see ImageCleaner.__call__()"""
+        """Apply FACT-style image cleaning. see ImageCleaner.__call__()"""
         return fact_image_cleaning(
             geom=self.subarray.tel[tel_id].camera.geometry,
             image=image,
@@ -440,3 +577,36 @@ class FACTImageCleaner(TailcutsImageCleaner):
             min_number_neighbors=self.min_picture_neighbors.tel[tel_id],
             time_limit=self.time_limit_ns.tel[tel_id],
         )
+
+
+class TimeConstrainedImageCleaner(TailcutsImageCleaner):
+    """
+    MAGIC-like Image cleaner with timing information (See `ctapipe.image.time_constrained_clean`)
+    """
+
+    time_limit_core_ns = FloatTelescopeParameter(
+        default_value=4.5,
+        help="arrival time limit for neighboring " "core pixels, in ns",
+    ).tag(config=True)
+    time_limit_boundary_ns = FloatTelescopeParameter(
+        default_value=1.5,
+        help="arrival time limit for neighboring " "boundary pixels, in ns",
+    ).tag(config=True)
+
+    def __call__(
+        self, tel_id: int, image: np.ndarray, arrival_times=None
+    ) -> np.ndarray:
+        """
+        Apply MAGIC-like image cleaning with timing information. See `ImageCleaner.__call__()`
+        """
+
+        return time_constrained_clean(
+            self.subarray.tel[tel_id].camera.geometry,
+            image,
+            arrival_times=arrival_times,
+            picture_thresh=self.picture_threshold_pe.tel[tel_id],
+            boundary_thresh=self.boundary_threshold_pe.tel[tel_id],
+            min_number_picture_neighbors=self.min_picture_neighbors.tel[tel_id],
+            time_limit_core=self.time_limit_core_ns.tel[tel_id],
+            time_limit_boundary=self.time_limit_boundary_ns.tel[tel_id],
+        )
diff --git a/ctapipe/instrument/camera/geometry.py b/ctapipe/instrument/camera/geometry.py
index 09a5eb8d..95e986f0 100644
--- a/ctapipe/instrument/camera/geometry.py
+++ b/ctapipe/instrument/camera/geometry.py
@@ -105,8 +105,6 @@ class CameraGeometry:
     cam_rotation: overall camera rotation with units
     """
 
-    _geometry_cache = {}  # dictionary CameraGeometry instances for speed
-
     def __init__(
         self,
         camera_name,
@@ -171,7 +169,7 @@ class CameraGeometry:
             self.rotate(self.cam_rotation)
 
             # cache border pixel mask per instance
-        self.border_cache = {}
+        self._border_cache = {}
 
     def __eq__(self, other):
         if self.camera_name != other.camera_name:
@@ -884,8 +882,8 @@ class CameraGeometry:
         mask: array
             A boolean mask, True if pixel is in the border of the specified width
         """
-        if width in self.border_cache:
-            return self.border_cache[width]
+        if width in self._border_cache:
+            return self._border_cache[width]
 
         # filter annoying deprecation warning from within scipy
         # scipy still uses np.matrix in scipy.sparse, but we do not
@@ -901,7 +899,7 @@ class CameraGeometry:
                 n = self.neighbor_matrix
                 mask = (n & self.get_border_pixel_mask(width - 1)).any(axis=1)
 
-        self.border_cache[width] = mask
+        self._border_cache[width] = mask
         return mask
 
     def position_to_pix_index(self, x, y):
diff --git a/ctapipe/instrument/camera/readout.py b/ctapipe/instrument/camera/readout.py
index 0c92e16e..39a92aef 100644
--- a/ctapipe/instrument/camera/readout.py
+++ b/ctapipe/instrument/camera/readout.py
@@ -20,6 +20,13 @@ def parse_dotted_version(version):
 
 
 class CameraReadout:
+    __slots__ = (
+        "camera_name",
+        "sampling_rate",
+        "reference_pulse_shape",
+        "reference_pulse_sample_width",
+    )
+
     def __init__(
         self,
         camera_name,
@@ -67,9 +74,9 @@ class CameraReadout:
         return hash(
             (
                 self.camera_name,
-                self.sampling_rate.to_value(u.GHz),
+                round(self.sampling_rate.to_value(u.GHz), 3),
                 self.reference_pulse_shape.size,
-                self.reference_pulse_sample_width.to_value(u.ns),
+                round(self.reference_pulse_sample_width.to_value(u.ns), 2),
             )
         )
 
diff --git a/ctapipe/instrument/optics.py b/ctapipe/instrument/optics.py
index 70146162..8103f067 100644
--- a/ctapipe/instrument/optics.py
+++ b/ctapipe/instrument/optics.py
@@ -40,7 +40,15 @@ class OpticsDescription:
         if the units of one of the inputs are missing or incompatible
     """
 
-    @u.quantity_input(mirror_area=u.m ** 2, equivalent_focal_length=u.m)
+    __slots__ = (
+        "equivalent_focal_length",
+        "mirror_area",
+        "name",
+        "num_mirror_tiles",
+        "num_mirrors",
+    )
+
+    @u.quantity_input(mirror_area=u.m**2, equivalent_focal_length=u.m)
     def __init__(
         self,
         name,
@@ -60,8 +68,8 @@ class OpticsDescription:
         """Make this hashable, so it can be used as dict keys or in sets"""
         return hash(
             (
-                self.equivalent_focal_length.to_value(u.m),
-                self.mirror_area,
+                round(self.equivalent_focal_length.to_value(u.m), 3),
+                round(self.mirror_area.to_value(u.m**2), 3),
                 self.num_mirrors,
                 self.num_mirror_tiles,
             )
diff --git a/ctapipe/instrument/telescope.py b/ctapipe/instrument/telescope.py
index 69c48ef8..e351b67d 100644
--- a/ctapipe/instrument/telescope.py
+++ b/ctapipe/instrument/telescope.py
@@ -39,6 +39,13 @@ class TelescopeDescription:
        the camera associated with this telescope
     """
 
+    __slots__ = (
+        "name",
+        "type",
+        "optics",
+        "camera",
+    )
+
     def __init__(
         self,
         name: str,
diff --git a/ctapipe/io/astropy_helpers.py b/ctapipe/io/astropy_helpers.py
index bcd2687c..5e75e4fc 100644
--- a/ctapipe/io/astropy_helpers.py
+++ b/ctapipe/io/astropy_helpers.py
@@ -89,9 +89,15 @@ def read_table(
 
         astropy_table = table_cls(array, meta=meta, copy=False)
         for column, tr in transforms.items():
+            if column not in astropy_table.colnames:
+                continue
+
             astropy_table[column] = tr.inverse(astropy_table[column])
 
         for column, desc in descriptions.items():
+            if column not in astropy_table.colnames:
+                continue
+
             astropy_table[column].description = desc
 
         return astropy_table
diff --git a/ctapipe/io/hdf5eventsource.py b/ctapipe/io/hdf5eventsource.py
index fa4aaeab..81cd64a8 100644
--- a/ctapipe/io/hdf5eventsource.py
+++ b/ctapipe/io/hdf5eventsource.py
@@ -275,11 +275,17 @@ class HDF5EventSource(EventSource):
             }
 
         if DataLevel.DL1_IMAGES in self.datalevels:
+            ignore_columns = {"parameters"}
+
+            # if there are no parameters, there are no image_mask, avoids warnings
+            if DataLevel.DL1_PARAMETERS not in self.datalevels:
+                ignore_columns.add("image_mask")
+
             image_readers = {
                 tel.name: self.reader.read(
                     f"/dl1/event/telescope/images/{tel.name}",
                     DL1CameraContainer,
-                    ignore_columns={"parameters"},
+                    ignore_columns=ignore_columns,
                 )
                 for tel in self.file_.root.dl1.event.telescope.images
             }
diff --git a/docs/ctapipe_api/instrument/camera_readout.rst b/docs/ctapipe_api/instrument/camera_readout.rst
index cc06bb56..1f94bb5e 100644
--- a/docs/ctapipe_api/instrument/camera_readout.rst
+++ b/docs/ctapipe_api/instrument/camera_readout.rst
@@ -32,9 +32,9 @@ Input/Output
 ------------
 
 
-You can write out a `CameraReadout` by using the `CameraReadout.to_table()`
- method to turn it into an ``astropy.table.Table``, and then call its ``write()``
-  function.  Reading it back in can be done with ``CameraReadout.from_table()``
+You can write out a `CameraReadout` by using the ``CameraReadout.to_table()``
+ method to turn it into a `~astropy.table.Table`, and then call its ``write()``
+  function.  Reading it back in can be done with ``CameraReadout.from_table()``.
 
 .. code-block:: python
 


Test patch:
diff --git a/ctapipe/image/tests/test_cleaning.py b/ctapipe/image/tests/test_cleaning.py
index eb921f8b..d004f7bc 100644
--- a/ctapipe/image/tests/test_cleaning.py
+++ b/ctapipe/image/tests/test_cleaning.py
@@ -198,7 +198,7 @@ def test_tailcuts_clean_min_neighbors_1():
 
 
 def test_tailcuts_clean_min_neighbors_2():
-    """ requiring that picture pixels have at least two neighbors above
+    """requiring that picture pixels have at least two neighbors above
     picture_thresh"""
 
     # start with simple 3-pixel camera
@@ -331,3 +331,128 @@ def test_apply_time_delta_cleaning():
     test_mask = mask.copy()
     test_mask[[41, 157]] = 0
     assert (test_mask == td_mask).all()
+
+
+def test_time_constrained_clean():
+    geom = CameraGeometry.from_name("LSTCam")
+    charge = np.zeros(geom.n_pixels, dtype=np.float64)
+    peak_time = np.zeros(geom.n_pixels, dtype=np.float64)
+
+    # define signal pixels and their charges/timings (1 core pixel + 6 neighboring core pixels + 12 neighboring boundary pixels)
+    core_pixel = 100
+    core_neighbors = geom.neighbors[core_pixel]
+    boundary_pixels = np.setdiff1d(
+        np.array([geom.neighbors[core_neighbor] for core_neighbor in core_neighbors]),
+        np.append(core_neighbors, core_pixel),
+    )
+    charge[core_pixel], charge[core_neighbors], charge[boundary_pixels] = 15, 10, 6
+    peak_time[core_pixel], peak_time[core_neighbors], peak_time[boundary_pixels] = (
+        18,
+        20,
+        21,
+    )
+
+    # define initial cleaning parameters
+    picture_thresh, boundary_thresh = 8, 4
+    time_limit_core, time_limit_boundary = 4.5, 1.5
+    min_number_picture_neighbors = 1
+
+    mask_signal = charge > 0
+
+    # 1. basic test
+    mask_reco = cleaning.time_constrained_clean(
+        geom,
+        charge,
+        peak_time,
+        picture_thresh=picture_thresh,
+        boundary_thresh=boundary_thresh,
+        time_limit_core=time_limit_core,
+        time_limit_boundary=time_limit_boundary,
+        min_number_picture_neighbors=min_number_picture_neighbors,
+    )
+    test_mask = mask_signal.copy()
+    assert (test_mask == mask_reco).all()
+
+    # 2. increased min_number_picture_neighbors test (here 3)
+    min_number_picture_neighbors = 3
+    mask_reco = cleaning.time_constrained_clean(
+        geom,
+        charge,
+        peak_time,
+        picture_thresh=picture_thresh,
+        boundary_thresh=boundary_thresh,
+        time_limit_core=time_limit_core,
+        time_limit_boundary=time_limit_boundary,
+        min_number_picture_neighbors=min_number_picture_neighbors,
+    )
+    ## removed pixels : boundary pixels
+    test_mask = mask_signal.copy()
+    test_mask[boundary_pixels] = 0
+    assert (test_mask == mask_reco).all()
+
+    # 3. strict time_limit_boundary test (here 0.5)
+    min_number_picture_neighbors = 1
+    time_limit_boundary = 0.5
+    mask_reco = cleaning.time_constrained_clean(
+        geom,
+        charge,
+        peak_time,
+        picture_thresh=picture_thresh,
+        boundary_thresh=boundary_thresh,
+        time_limit_core=time_limit_core,
+        time_limit_boundary=time_limit_boundary,
+        min_number_picture_neighbors=min_number_picture_neighbors,
+    )
+    ## removed pixels : boundary pixels
+    test_mask = mask_signal.copy()
+    test_mask[boundary_pixels] = 0
+    assert (test_mask == mask_reco).all()
+
+    # 4. time_limit_core test (one of core_neighbors have peak time >5 slice away from the average)
+    time_limit_boundary = 1.5
+    noise_core_neighbor = core_neighbors[0]
+    peak_time[noise_core_neighbor] = 25
+    mask_reco = cleaning.time_constrained_clean(
+        geom,
+        charge,
+        peak_time,
+        picture_thresh=picture_thresh,
+        boundary_thresh=boundary_thresh,
+        time_limit_core=time_limit_core,
+        time_limit_boundary=time_limit_boundary,
+        min_number_picture_neighbors=min_number_picture_neighbors,
+    )
+    ## removed pixels : the noise core neighbor pixel + one neighboring boundary
+    test_mask = mask_signal.copy()
+    test_mask[noise_core_neighbor] = 0
+    noise_boundary = np.setdiff1d(
+        geom.neighbors[noise_core_neighbor],
+        np.array(
+            [geom.neighbors[core_neighbor] for core_neighbor in core_neighbors[1:]]
+        ),
+    )
+    test_mask[noise_boundary] = 0
+    assert (test_mask == mask_reco).all()
+
+    # 5. time_limit_core test for brighter pixels (one of core_neighbors have peak time >5 slice away from the average)
+    charge[core_pixel], charge[core_neighbors] = 30, 20
+    mask_reco = cleaning.time_constrained_clean(
+        geom,
+        charge,
+        peak_time,
+        picture_thresh=picture_thresh,
+        boundary_thresh=boundary_thresh,
+        time_limit_core=time_limit_core,
+        time_limit_boundary=time_limit_boundary,
+        min_number_picture_neighbors=min_number_picture_neighbors,
+    )
+    ## removed pixels : one neighboring boundary to the noise core pixel
+    test_mask = mask_signal.copy()
+    noise_boundary = np.setdiff1d(
+        geom.neighbors[noise_core_neighbor],
+        np.array(
+            [geom.neighbors[core_neighbor] for core_neighbor in core_neighbors[1:]]
+        ),
+    )
+    test_mask[noise_boundary] = 0
+    assert (test_mask == mask_reco).all()
diff --git a/ctapipe/io/tests/test_astropy_helpers.py b/ctapipe/io/tests/test_astropy_helpers.py
index a7504b3c..718caf09 100644
--- a/ctapipe/io/tests/test_astropy_helpers.py
+++ b/ctapipe/io/tests/test_astropy_helpers.py
@@ -188,3 +188,22 @@ def test_read_table_astropy(tmp_path):
     table.write(path, "/group/table", serialize_meta=True)
     read = read_table(path, "/group/table")
     assert_table_equal(table, read)
+
+
+def test_unknown_meta(tmp_path):
+    """Test read_table ignores unknown metadata
+
+    Regression test for https://github.com/cta-observatory/ctapipe/issues/1785
+    """
+    container = ReconstructedEnergyContainer()
+    filename = tmp_path / "test_astropy_table.h5"
+
+    with HDF5TableWriter(filename) as writer:
+        for energy in [np.nan, 100, np.nan, 50, -1.0] * u.TeV:
+            container.energy = energy
+            writer.write("events", container)
+
+        # unknown column matching ctapipe metadata
+        writer.h5file.root.events._v_attrs["foo_UNIT"] = "TeV"
+
+    read_table(filename, "/events", condition="energy > 0")