Vectorize operations in HillasReconstructor

[maxnoe]

Timing script:

from ctapipe.io import EventSource
from ctapipe.reco import HillasReconstructor
import numpy as np
from timeit import repeat
from cProfile import Profile



source = EventSource("../data/gamma_20deg_180deg_run3___cta-prod6-paranal-2147m-Paranal-dark.dl1.h5", max_events=200)
events = list(source)

reco = HillasReconstructor(subarray=source.subarray)

code = '''
print("start")
for event in events:
    reco(event)
'''

with Profile() as profile:
    exec(code)

profile.dump_stats("hillas_reco.profile")


n_repeat = 5
number = 1
times = np.array(repeat(code, number=number, repeat=n_repeat, globals=globals()))
times /= number

On master:

11.73 ± 0.31 s

on this branch:

2.25 ± 0.02 s

Units tests are passing, but I will also process a larger DL2 run and check that results are the same.

Most time is not spend in astropy coordinate transforms and project_to_ground.

[codecov]

Codecov Report

Base: 92.49% // Head: 92.49% // Decreases project coverage by -0.00% ⚠️

Coverage data is based on head (cf7f20e) compared to base (4905dc1).
Patch coverage: 94.53% of modified lines in pull request are covered.

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #2036      +/-   ##
==========================================
- Coverage   92.49%   92.49%   -0.01%     
==========================================
  Files         197      197              
  Lines       16246    16276      +30     
==========================================
+ Hits        15027    15054      +27     
- Misses       1219     1222       +3     

Impacted Files	Coverage Δ
ctapipe/reco/tests/test_reconstruction_methods.py	100.00% <ø> (ø)
ctapipe/coordinates/nominal_frame.py	89.65% <66.66%> (ø)
ctapipe/reco/hillas_reconstructor.py	95.77% <94.17%> (-2.62%)	⬇️
ctapipe/reco/tests/test_HillasReconstructor.py	99.23% <98.00%> (-0.77%)	⬇️
ctapipe/coordinates/__init__.py	100.00% <100.00%> (ø)
ctapipe/coordinates/ground_frames.py	100.00% <100.00%> (ø)
ctapipe/coordinates/telescope_frame.py	90.00% <100.00%> (+0.34%)	⬆️
ctapipe/instrument/camera/geometry.py	90.47% <0.00%> (+0.56%)	⬆️

Help us with your feedback. Take ten seconds to tell us how you rate us. Have a feature suggestion? Share it here.

☔ View full report at Codecov.
📢 Do you have feedback about the report comment? Let us know in this issue.

[kosack]

as suggested in #1922 perhaps we should also drop the CameraFrame support in HillasReconstructor (either here or in a separate PR) to simplify this code even more and remove a few branch points.

[maxnoe]

@kosack If you look at the code, the branching for CameraFrame vs. TelescopeFrame is now very minimal, so I'd say not really a concern here.

I am investigating a few issues with astropy concerning performance and I have identified several spots where astropy can be much improved.

For now, I would keep it in the current version I think, we already get a very noticeable speedup.

[kosack]

Looks good!

Before merging, It might be a good idea to do a regression test on a largish number of events to be sure the results have not changed since before the refactoring. (if you have not done that already).

[maxnoe]

Will do

[kosack]

One minor annoyance that would be nice to fix in this is that the azimuths end up in a strange quadrant (negative). It would be nice for sources with azimuth=180.0 that the reconstructed position is centered there, not at -180.0. It's still correct, but may be confusing when we use that value in a ML algorithm without any transforms.

E.g.: can it return values more like the right plot?

[maxnoe]

Sure, I'll wrap the angle into [-180, 180)

[kosack]

Astropy.coord seems to wrap to (0,360], not sure which is better (either way the user has to do a wrap transform if you are one of the usual North or South pointings). At least we should make all Reconstructors consistent. I think HillasIntersection returns something in (0,360]. Not sure about ImPACT.

I only noticed this in #2073 when using multiple reconstructors

[maxnoe]

Ok, let's go with the astropy default then

[maxnoe]

I processed 100 gamma diffuse runs (~25k successfully reconstructed events) using 0.16.0 and this branch.

The results are not exactly identical, but they are very close for almost all events and the angular resolution is the same or even slightly better with the new code:

[kosack]

I guess the small differences are just due to the different order of operations. In any case as the results are generally better, I think it's ok.

[maxnoe]

Had to rebase, please re-review @kosack @nbiederbeck @LukasNickel

[LukasNickel]

AttributeError: 'Quantity' object has no 'deg' member

error in docs

[LukasNickel]

Not directly related to this PR, but is anyone else really disliking this code block and would rather have the cam radius in meter and degree be a property of the camera geometry?

[maxnoe]

We currently have the geometry strictly in one frame and I quite like this...

[maxnoe]

This should rather be replaced by dropping support for doing high-level things in camera frame

[LukasNickel]

For my understanding:
Why not directly assign to self._cam_radius_x above?
I am guessing you save calculations by performing them once per camery type above and only then looping over the actual telescopes down here?

[maxnoe]

Errors should be fixed now