Repository: cta-observatory/ctapipe (python)

Issue:
Neighbormatrix fails to compute for geometry with a single pixel
### Discussed in https://github.com/cta-observatory/ctapipe/discussions/2309

<div type='discussions-op-text'>

<sup>Originally posted by **Elisa-Visentin** April 19, 2023</sup>
Hi,
Sorry for the banal question but... While computing morphology parameters (`morphology_parameters(...)`) on a set of images with `ctapipe 0.19.0`, some of the events make the script crash with the following error lines
```
Traceback (most recent call last):
  File "/home/evisentin/Desktop/lstchain/magic-cta-pipe/magicctapipe/scripts/lst1_magic/lst1_magic_mc_dl0_to_dl1.py", line 522, in <module>
    main()
  File "/home/evisentin/Desktop/lstchain/magic-cta-pipe/magicctapipe/scripts/lst1_magic/lst1_magic_mc_dl0_to_dl1.py", line 513, in main
    mc_dl0_to_dl1(args.input_file, args.output_dir, config)
  File "/home/evisentin/Desktop/lstchain/magic-cta-pipe/magicctapipe/scripts/lst1_magic/lst1_magic_mc_dl0_to_dl1.py", line 343, in mc_dl0_to_dl1
    morphology_params = morphology_parameters(
  File "/home/evisentin/mambaforge/envs/lstchain/lib/python3.9/site-packages/ctapipe/image/morphology.py", line 194, in morphology_parameters
    n_islands, island_labels = number_of_islands(geom=geom, mask=image_mask)
  File "/home/evisentin/mambaforge/envs/lstchain/lib/python3.9/site-packages/ctapipe/image/morphology.py", line 74, in number_of_islands
    neighbors = geom.neighbor_matrix_sparse
  File "/home/evisentin/mambaforge/envs/lstchain/lib/python3.9/site-packages/astropy/utils/decorators.py", line 841, in __get__
    val = self.fget(obj)
  File "/home/evisentin/mambaforge/envs/lstchain/lib/python3.9/site-packages/ctapipe/instrument/camera/geometry.py", line 682, in neighbor_matrix_sparse
    return self.calc_pixel_neighbors(diagonal=False)
  File "/home/evisentin/mambaforge/envs/lstchain/lib/python3.9/site-packages/ctapipe/instrument/camera/geometry.py", line 740, in calc_pixel_neighbors
    if (neighbor_matrix.T != neighbor_matrix).sum() > 0:
AttributeError: 'bool' object has no attribute 'sum'

```
whereas other events (same script, same input file) are well reconstructed or only return a warning 
```
/home/evisentin/mambaforge/envs/lstchain/lib/python3.9/site-packages/ctapipe/instrument/camera/geometry.py:741: UserWarning: Neighbor matrix is not symmetric. Is camera geometry irregular?
  warnings.warn(

```
which needs this `if` to be executed.
Maybe some events need a lot of memory to be processed? Or are there some 'anomalous' images on which the `morphology_parameters` algorithm fails? Other issues (dependencies...)?

Thanks.

</div>

Pull request description:
Fix neighbor matrix for no pixels, fixes #2316

Patch:
diff --git a/.mailmap b/.mailmap
index d431d96f..466d1714 100644
--- a/.mailmap
+++ b/.mailmap
@@ -88,6 +88,3 @@ Orel Gueta <orel.gueta@desy.de> orelgueta <22728856+orelgueta@users.noreply.gith
 Rune Michael Dominik <rune.dominik@tu-dortmund.de>
 
 Henning Ptaszyk <henning.ptaszyk@udo.edu>
-
-
-Tomas Bylund <tomas.bylund@cea.fr> Tomas Bylund <tobychev@live.com>
diff --git a/ctapipe/core/tool.py b/ctapipe/core/tool.py
index e39dcaff..56de387a 100644
--- a/ctapipe/core/tool.py
+++ b/ctapipe/core/tool.py
@@ -405,14 +405,11 @@ class Tool(Application):
 
         with self._exit_stack:
             try:
+                self.initialize(argv)
                 self.log.info(f"Starting: {self.name}")
                 Provenance().start_activity(self.name)
-
-                self.initialize(argv)
-
                 self.setup()
                 self.is_setup = True
-
                 self.log.debug(f"CONFIG: {self.get_current_config()}")
                 Provenance().add_config(self.get_current_config())
 
diff --git a/ctapipe/image/morphology.py b/ctapipe/image/morphology.py
index 97ca48ae..834efbdc 100644
--- a/ctapipe/image/morphology.py
+++ b/ctapipe/image/morphology.py
@@ -57,20 +57,26 @@ def number_of_islands(geom, mask):
 
     Parameters
     ----------
-    geom: `~ctapipe.instrument.CameraGeometry`
-        Camera geometry information
-    mask: ndarray
-        input mask (array of booleans)
+    geom : `~ctapipe.instrument.CameraGeometry`
+        Camera geometry information, needs to be the full camrea geometry,
+        not already masked with ``mask``
+    mask : ndarray
+        input mask (array of booleans) of pixels surviving the cleaning
 
     Returns
     -------
-    n_islands: int
+    n_islands : int
         Total number of clusters
-    island_labels: ndarray
+    island_labels : ndarray
         Contains cluster membership of each pixel.
         Dimension equals input geometry.
         Entries range from 0 (not in the pixel mask) to n_islands.
     """
+    if geom.n_pixels != len(mask):
+        raise ValueError(
+            "CameraGeometry has less pixels than mask"
+            ", number_of_islands needs the full CameraGeometry"
+        )
     neighbors = geom.neighbor_matrix_sparse
     n_islands, island_labels = _n_islands_sparse_indices(
         neighbors.indices, neighbors.indptr, mask
diff --git a/ctapipe/instrument/camera/geometry.py b/ctapipe/instrument/camera/geometry.py
index 079fd45a..042f945d 100644
--- a/ctapipe/instrument/camera/geometry.py
+++ b/ctapipe/instrument/camera/geometry.py
@@ -691,6 +691,9 @@ class CameraGeometry:
         diagonal: bool
             If rectangular geometry, also add diagonal neighbors
         """
+        if self.n_pixels <= 1:
+            return csr_matrix(np.ones((self.n_pixels, self.n_pixels), dtype=bool))
+
         # assume circle pixels are also on a hex grid
         if self.pix_type in (PixelShape.HEXAGON, PixelShape.CIRCLE):
             max_neighbors = 6
diff --git a/docs/changes/2312.bugfix.rst b/docs/changes/2312.bugfix.rst
deleted file mode 100644
index cdcf85d1..00000000
--- a/docs/changes/2312.bugfix.rst
+++ /dev/null
@@ -1,1 +0,0 @@
-Fix for config files not being included as inputs in provenance log.
diff --git a/docs/changes/2317.bugfix.rst b/docs/changes/2317.bugfix.rst
new file mode 100644
index 00000000..b6c24d72
--- /dev/null
+++ b/docs/changes/2317.bugfix.rst
@@ -0,0 +1,1 @@
+Fix calculation of the neighbor matrix of ``CameraGeometry`` for empty and single-pixel geometries.


Test patch:
diff --git a/ctapipe/core/tests/test_tool.py b/ctapipe/core/tests/test_tool.py
index a65ecf83..8db62479 100644
--- a/ctapipe/core/tests/test_tool.py
+++ b/ctapipe/core/tests/test_tool.py
@@ -426,38 +426,3 @@ def test_exit_stack():
     assert run_tool(tool, raises=False) == 1
     assert tool.manager.enter_called
     assert tool.manager.exit_called
-
-
-def test_activity(tmp_path):
-    """check that the config is correctly in the provenance"""
-
-    class MyTool(Tool):
-        name = "test_prov_log"
-        description = "test"
-        userparam = Float(5.0, help="parameter").tag(config=True)
-
-    tool = MyTool()
-
-    config_path = tmp_path / "config.json"
-    config_path.write_text(json.dumps({"MyTool": {"userparam": 10.0}}))
-    provenance_path = tmp_path / "provenance.json"
-
-    run_tool(
-        tool,
-        [
-            "--config",
-            str(config_path),
-            f"--provenance-log={provenance_path}",
-        ],
-    )
-
-    activities = json.loads(tool.provenance_log.read_text())
-    # provlog contains all activities from all tests, last one is the tool we just ran
-    provlog = activities[-1]
-    assert provlog["activity_name"] == MyTool.name
-
-    # test config file is in inputs, regression test for #2313
-    inputs = provlog["input"]
-    assert len(inputs) == 1
-    assert inputs[0]["role"] == "Tool Configuration"
-    assert inputs[0]["url"] == str(config_path)
diff --git a/ctapipe/image/tests/test_morphology.py b/ctapipe/image/tests/test_morphology.py
index cabab0eb..2b3904ad 100644
--- a/ctapipe/image/tests/test_morphology.py
+++ b/ctapipe/image/tests/test_morphology.py
@@ -1,4 +1,5 @@
 import numpy as np
+import pytest
 from numpy.testing import assert_allclose
 
 from ctapipe.instrument import CameraGeometry
@@ -124,3 +125,19 @@ def test_brightest_island():
     no_island_mask = brightest_island(0, island_labels, image)
     assert len(no_island_mask) == 3
     assert np.all(no_island_mask == False)
+
+
+def test_number_of_islands_masked(prod3_lst):
+    from ctapipe.image import number_of_islands
+
+    geom = prod3_lst.camera.geometry
+
+    # create 18 triggered pixels grouped to 5 clusters
+    mask = np.zeros(geom.n_pixels).astype("bool")
+    triggered_pixels = np.array(
+        [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 37, 38, 111, 222]
+    )
+    mask[triggered_pixels] = True
+
+    with pytest.raises(ValueError, match="needs the full CameraGeometry"):
+        number_of_islands(geom[mask], mask)
diff --git a/ctapipe/instrument/camera/tests/test_geometry.py b/ctapipe/instrument/camera/tests/test_geometry.py
index 536da776..1745319a 100644
--- a/ctapipe/instrument/camera/tests/test_geometry.py
+++ b/ctapipe/instrument/camera/tests/test_geometry.py
@@ -372,3 +372,19 @@ def test_guess_radius(prod5_lst, prod5_sst):
 
     prod5_chec = prod5_sst.camera.geometry
     assert u.isclose(prod5_chec.guess_radius(), 0.16 * u.m, rtol=0.05)
+
+
+def test_single_pixel(prod5_lst):
+    """Regression test for #2316"""
+    single_pixel = prod5_lst.camera.geometry[[0]]
+
+    assert single_pixel.neighbor_matrix.shape == (1, 1)
+    assert single_pixel.neighbor_matrix[0, 0]
+
+
+def test_empty(prod5_lst):
+    geometry = prod5_lst.camera.geometry
+    mask = np.zeros(len(geometry), dtype=bool)
+    empty = geometry[mask]
+
+    assert empty.neighbor_matrix.shape == (0, 0)