StereoDispCombiner: Add Transformation to Nominalframe so divergent pointing is supported

[AkramAki]

This introduces support for divergent pointing for the StereoDispCombiner. Since the StereoDispCombiner PR is not finished yet this PR targets the corresponding branch instead of main.

[LukasBeiske]

@maxnoe two things:

• Does this need a changelog, since it's not getting merged into main?
• FYI: I had to approve the starting of the CI manually. I guess this is related to the new github security rules... Just in case that is not the intended behaviour .

[maxnoe]

Does this need a changelog, since it's not getting merged into main?

no

FYI: I had to approve the starting of the CI manually.

This is not new, any workflow of a first time outside contributor has to be approved first.

[LukasBeiske]

@maxnoe And just to be sure: We can merge this into the stereo_combiner branch at our own discretion once the tests pass, correct?

[maxnoe]

Sure, that's your PR branch

[ctao-sonarqube]

Quality Gate passed

Issues
0 New issues
0 Accepted issues

Measures
0 Security Hotspots
100.0% Coverage on New Code
0.0% Duplication on New Code

See analysis details on SonarQube

[Hckjs]

#3076 (comment)

see above. Why is it necessary to add the unit here? Shouldn't they already have a unit?

A couple of comments to @maxnoe comment:
I guess you have a point here, that the current unit handling is a bit messy, because in an earlier version of the stereo_combiner branch, the idea was to get rid of the units at the beginning and add them back later in the end because of some njit'ed functions later on an
hillas_fov_lon + signs * disp * np.cos(hillas_psi) in the beginning
See some old discussion here: #2731 (comment)

After some generalization i merged the calculation of the minium distances to one function used by the table-wise and event-wise processing (calc_combs_min_distances), which can handle astropy units. I did unfortunatly not simplify the unit handling accordingly...

The only function, which currently cannot handle astropy units is weighted_mean_std_ufunc. For the rest i think we can save some overhead by reducing the conversions/strip of units to a minimum. i'll make some suggestions in the following

/edit see #3076 (comment)

[maxnoe]

It's fine to strip units for performance, but then one needs to be careful when and how they are added back.

E.g. this is bad:

SkyCoord(alt=alt * u.deg, az=az*u.deg, frame=AltAz())

because it makes unnecessary copies of potentially large arrays.

These can be avoided by not using the naive * unit syntax, e.g. using the u.Quantity(array, unit, copy=False) approach or even passing unit= to SkyCoord:

In [1]: import numpy as np

In [2]: import astropy.units as u

In [3]: from astropy.coordinates import SkyCoord, AltAz

In [4]: frame = AltAz()

In [5]: az = np.random.normal(90, 5, size=10000)

In [6]: alt = np.random.normal(70, 2, size=10000)

In [7]: %timeit SkyCoord(alt=alt * u.deg, az=az * u.deg, frame=frame)
78 μs ± 812 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

In [8]: %timeit SkyCoord(alt=alt, az=az, unit=u.deg, copy=False, frame=frame)
60.6 μs ± 935 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

[Hckjs]

in telescope_event_handling.py you also have to keep the quantities in calc_fov_lon_lat, so just:

    hillas_fov_lon = tel_table["hillas_fov_lon"].quantity
    hillas_fov_lat = tel_table["hillas_fov_lat"].quantity
    hillas_psi = tel_table["hillas_psi"].quantity
    disp = tel_table[f"{prefix}_parameter"].quantity

[Hckjs]

I think it would be also nice to have a unit test for an event with at least two participating telescopes with different pointings.

[LukasBeiske]

I think it would be also nice to have a unit test for an event with at least two participating telescopes with different pointings.

If we just change the telescope pointings in the dummy table to not be identical for all telescope events, this should be covered, no?

[Hckjs]

I just checked the performances for single array events in the event-loop and it looks like its much faster calculating the fov_lon/lats and the minimum distances without units (also adding the strip and wrap back of the units...):

import numpy as np
import astropy.units as u

signs = np.array([-1, 1])

hillas_fov_lon = 0.5 * u.deg
disp = 0.8 * u.deg
hillas_psi = 40 * u.deg

%timeit hillas_fov_lon + signs * disp * np.cos(hillas_psi)
%timeit (u.Quantity(hillas_fov_lon.to_value(u.deg) + signs * disp.to_value(u.deg) * np.cos(hillas_psi.to_value(u.rad)), u.deg, copy=False,))

28.8 μs ± 344 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)
7.99 μs ± 40.9 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)

# Without units
%timeit calc_combs_min_distances(index_tel_combs, fov_lon, fov_lat, weights)
68 μs ± 194 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

# With units
%timeit calc_combs_min_distances(index_tel_combs, fov_lon_q, fov_lat_q, weights)
386 μs ± 1.31 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

Im sorry, didn't expect to be that much of a difference. I would just leave it as it is.

[Hckjs]

I think it would be also nice to have a unit test for an event with at least two participating telescopes with different pointings.

If we just change the telescope pointings in the dummy table to not be identical for all telescope events, this should be covered, no?

If you write a new unit test, adapting the dummy table just inside the test, i think yes. Otherwise i guess you also have to adapt all the expected reconstructed parameter aswell in the other tests