Neighbormatrix fails to compute for geometry with a single pixel

[maxnoe]

Discussed in https://github.com/cta-observatory/ctapipe/discussions/2309

^{Originally posted by **Elisa-Visentin** April 19, 2023} 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 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.

[maxnoe]

Minimal example:

In [1]: from ctapipe.instrument import CameraGeometry

In [2]: cam = CameraGeometry.from_name("LSTCam")

In [3]: one_pixel = cam[[0]]

In [4]: one_pixel.neighbor_matrix
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[4], line 1
----> 1 one_pixel.neighbor_matrix

File ~/.local/conda/envs/cta-dev/lib/python3.9/site-packages/astropy/utils/decorators.py:841, in lazyproperty.__get__(self, obj, owner)
    839         val = obj_dict.get(self._key, _NotFound)
    840         if val is _NotFound:
--> 841             val = self.fget(obj)
    842             obj_dict[self._key] = val
    843 return val

File ~/Uni/CTA/ctapipe/ctapipe/instrument/camera/geometry.py:671, in CameraGeometry.neighbor_matrix(self)
    669 @lazyproperty
    670 def neighbor_matrix(self):
--> 671     return self.neighbor_matrix_sparse.A

File ~/.local/conda/envs/cta-dev/lib/python3.9/site-packages/astropy/utils/decorators.py:841, in lazyproperty.__get__(self, obj, owner)
    839         val = obj_dict.get(self._key, _NotFound)
    840         if val is _NotFound:
--> 841             val = self.fget(obj)
    842             obj_dict[self._key] = val
    843 return val

File ~/Uni/CTA/ctapipe/ctapipe/instrument/camera/geometry.py:682, in CameraGeometry.neighbor_matrix_sparse(self)
    680     return self._neighbors
    681 else:
--> 682     return self.calc_pixel_neighbors(diagonal=False)

File ~/Uni/CTA/ctapipe/ctapipe/instrument/camera/geometry.py:733, in CameraGeometry.calc_pixel_neighbors(self, diagonal)
    730 neighbors = neighbor_candidates[pixels, neigbor_index]
    731 data = np.ones(len(pixels), dtype=bool)
--> 733 neighbor_matrix = csr_matrix((data, (pixels, neighbors)))
    735 # filter annoying deprecation warning from within scipy
    736 # scipy still uses np.matrix in scipy.sparse, but we do not
    737 # explicitly use any feature of np.matrix, so we can ignore this here
    738 with warnings.catch_warnings():

File ~/.local/conda/envs/cta-dev/lib/python3.9/site-packages/scipy/sparse/_compressed.py:53, in _cs_matrix.__init__(self, arg1, shape, dtype, copy)
     49 else:
     50     if len(arg1) == 2:
     51         # (data, ij) format
     52         other = self.__class__(
---> 53             self._coo_container(arg1, shape=shape, dtype=dtype)
     54         )
     55         self._set_self(other)
     56     elif len(arg1) == 3:
     57         # (data, indices, indptr) format

File ~/.local/conda/envs/cta-dev/lib/python3.9/site-packages/scipy/sparse/_coo.py:148, in coo_matrix.__init__(self, arg1, shape, dtype, copy)
    146 if shape is None:
    147     if len(row) == 0 or len(col) == 0:
--> 148         raise ValueError('cannot infer dimensions from zero '
    149                          'sized index arrays')
    150     M = operator.index(np.max(row)) + 1
    151     N = operator.index(np.max(col)) + 1

ValueError: cannot infer dimensions from zero sized index arrays

[Elisa-Visentin]

Hi. I made the script print some cleaning information (maybe this could help you): Masked CameraGeometry (pixels x-coordinate):

[-0.47245608 -0.48190485 -0.49135362 -0.5008024  -0.51025117 -0.4913551
 -0.50080387 -0.51025265 -0.51970142 -0.52915019 -0.53859897 -0.54804774
 -0.5102543  -0.51970307 -0.52915185 -0.53860062 -0.54804939 -0.55749817
 -0.56694694 -0.57639571 -0.58584449 -0.52915344 -0.53860221 -0.54805098
 -0.55749976 -0.56694853 -0.57639731 -0.58584608 -0.59529485 -0.60474363
 -0.6141924  -0.62364117 -0.56695012 -0.5763989  -0.58584767 -0.59529645
 -0.60474522 -0.61419399 -0.62364277 -0.63309154 -0.64254031 -0.65198909
 -0.51970466 -0.60474675 -0.61419553 -0.6236443  -0.63309307 -0.64254185
 -0.65199062 -0.66143939 -0.67088817 -0.68033694 -0.51025589 -0.64254344
 -0.65199221 -0.66144099 -0.67088976 -0.68033853 -0.68978731 -0.69923608
 -0.70868485 -0.54805258 -0.59529798 -0.68034013 -0.6897889  -0.69923767
 -0.70868644 -0.71813522 -0.72758399 -0.63309467 -0.72758552 -0.75593185
 -0.78427817 -0.66144258 -0.70868798 -0.74648466 -0.7653838 ] m

77 pixels surviving the cleaning Image:

Script: magic-cta-pipe repository, lstchain branch, magicctapipe/scripts/lst1_magic/lst1_magic_mc_dl0_to_dl1.py (now 'fixed' by a try/except)

[maxnoe]

You see the error you posted with a matrix with 77 pixels? Are you sure?

I get the exact same error only for 1 pixel in the CameraGeometry.

[Elisa-Visentin]

Here a more complete log:

count   20    tel_id    1 [ 1.47291807e-01  9.81927157e-02  1.80024222e-01  1.30925131e-01
  8.18260377e-02  3.27269503e-02 -6.54712247e-02 -1.14570327e-01
 -1.63669400e-01  1.63657557e-01  1.14558466e-01  6.54593746e-02
  1.63602799e-02 -3.27388076e-02 -8.18378951e-02 -1.30936997e-01
 -1.80036070e-01  1.47290892e-01  9.81917990e-02  4.90927116e-02
 -6.37588515e-06 -4.91054633e-02 -9.82045654e-02 -1.47303638e-01
 -1.96402740e-01  3.27260299e-02 -1.63730575e-02 -6.54721450e-02
 -1.14571247e-01 -1.63670320e-01 -2.12769422e-01  1.63593742e-02
 -3.27397133e-02 -8.18388154e-02 -1.30937888e-01 -1.80036990e-01
 -1.47304559e-01 -1.96403661e-01] m
count   20    tel_id    1
count   21    tel_id    1 [-0.16366303 -0.21276211 -0.19639544 -0.18002879 -0.24549453 -0.22912788] m
count   21    tel_id    1
--> 22 event (event ID: 11209, telescope 1) could not survive the image cleaning. Skipping...
count   23    tel_id    1 [0.13092696 0.21275848 0.16365939 0.19639181 0.14729272] m
count   23    tel_id    1
--> 24 event (event ID: 11708, telescope 1) could not survive the image cleaning. Skipping...
count   25    tel_id    1 [-6.54621192e-02 -1.63630262e-02  3.27360612e-02  8.18351487e-02
  1.30934236e-01 -1.96392705e-01 -1.47293614e-01 -9.81945229e-02
 -4.90954317e-02  3.65937145e-06  4.91027542e-02  9.82018416e-02
 -3.27323316e-01 -2.78224228e-01 -2.29125141e-01 -1.80026053e-01
 -1.30926960e-01 -8.18278692e-02 -3.27287763e-02  1.63703112e-02
  6.54693986e-02 -4.58253922e-01 -4.09154849e-01 -3.60055747e-01
 -3.10956660e-01 -2.61857572e-01 -2.12758478e-01 -1.63659387e-01
 -1.14560295e-01 -6.54612043e-02 -1.63621132e-02  3.27369816e-02
 -4.90986354e-01 -4.41887252e-01 -3.92788179e-01 -3.43689077e-01
 -2.94589990e-01 -2.45490902e-01 -1.96391815e-01 -1.47292722e-01
 -9.81936306e-02 -4.90945376e-02 -5.07353024e-01 -5.23718774e-01
 -4.74619672e-01 -4.25520600e-01 -3.76421497e-01 -3.27322410e-01
 -2.78223323e-01 -2.29124228e-01 -1.80025137e-01 -1.30926046e-01
 -5.56452097e-01 -5.56451206e-01 -5.07352104e-01 -4.58253002e-01
 -4.09153929e-01 -3.60054827e-01 -3.10955740e-01 -2.61856652e-01
 -2.12757565e-01 -5.89184529e-01 -5.72817847e-01 -5.89183609e-01
 -5.40084536e-01 -4.90985434e-01 -4.41886361e-01 -3.92787259e-01
 -3.43688171e-01 -6.05550279e-01 -5.72816956e-01 -4.25519679e-01
 -2.78222402e-01 -6.87381813e-01 -6.71015102e-01 -7.03747563e-01
 -8.34678170e-01] m
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 534, in <module>
    main()
  File "/home/evisentin/Desktop/lstchain/magic-cta-pipe/magicctapipe/scripts/lst1_magic/lst1_magic_mc_dl0_to_dl1.py", line 525, 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 354, 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'

code:

      print('count  ',event.count,'   tel_id   ',tel_id,camera_geom_masked.pix_y)
      #print(event.count,tel_id, signal_pixels)
      #print(event.count, signal_pixels[signal_pixels == True], len(signal_pixels[signal_pixels == True]))
      if len(signal_pixels[signal_pixels == True])<2:
          logger.info(
              f"--> {event.count} event (event ID: {event.index.event_id}, "
              f"telescope {tel_id}) only one pixel surviving cleaning. Skipping..."
          )
          continue
      # Parametrize the image
      hillas_params = hillas_parameters(camera_geom_masked, image_masked)

      timing_params = timing_parameters(
          camera_geom_masked, image_masked, peak_time_masked, hillas_params
      )
      concentration_params = concentration_parameters(
          camera_geom_masked, image_masked, hillas_params
      )

      morphology_params = morphology_parameters(
          camera_geom_masked, signal_pixels
      )
      print('count  ',event.count, '   tel_id   ',tel_id)

It apparently crashes on the 25th event, which has this huge Masked CameraGeometry array. In case of good events, the event count is printed both before and after morphology_parameters() (as for event 23). This is the first event that makes the script crash... (Sorry, again)

[maxnoe]

@Elisa-Visentin I see that in the stage1 tool, we pass the full geometry to the morphology parameter calculation.

I need to check if that is required or just an accident, but could you try the same? Pass the full geometry, not the goemetry selected (only for the morphology parameters)

[Elisa-Visentin]

By passing the full geometry (not masked) it works! Thank you very much! And sorry for having disturbed you, I should have checked this... (coding isn't my thing)

[maxnoe]

@Elisa-Visentin Please stop apologizing :wink: You identified a valid issue and we need to fix that.

At the very least, this is not properly documented but it still might be a bug that we need to fix!

[Elisa-Visentin]

With a full geometry, even the warnings disappear; maybe this is due to scipy (csr_matrix). But adding a note in the docs could help (and solve the problem).

Please stop apologizing 😉

I must go back to exploit somebody else for questions/issues... :sweat_smile: