Use HillasParametersContainer only

ctapipe/image/hillas.py

@@ -4,8 +4,6 @@
 Hillas-style moment-based shower image parametrization.
 """
 
-from collections import namedtuple
-
 import astropy.units as u
 import numpy as np
 from astropy.coordinates import Angle
@@ -16,7 +14,6 @@
 
 
 __all__ = [
-    'MomentParameters',
     'hillas_parameters',
     'hillas_parameters_1',
     'hillas_parameters_2',
@@ -26,11 +23,6 @@
     'HillasParameterizationError',
 ]
 
-MomentParameters = namedtuple(
-    "MomentParameters",
-    "size,cen_x,cen_y,length,width,r,phi,psi,miss,skewness,kurtosis"
-)
-
 
 class HillasParameterizationError(RuntimeError):
     pass
@@ -55,7 +47,7 @@ def hillas_parameters_1(geom: CameraGeometry, image):
 
     Returns
     -------
-    hillas_parameters : `MomentParameters`
+    hillas_parameters : `HillasParametersContainer`
     """
     unit = Quantity(geom.pix_x).unit
     pix_x = Quantity(np.asanyarray(geom.pix_x, dtype=np.float64)).value
@@ -118,7 +110,7 @@ def hillas_parameters_1(geom: CameraGeometry, image):
     width = 0. if width_2 < 0. else np.sqrt(width_2)
     length = 0. if length_2 < 0. else np.sqrt(length_2)
 
-    miss = np.abs(b / np.sqrt(1 + a * a))
+    # miss = np.abs(b / np.sqrt(1 + a * a))
     r = np.sqrt(mean_x * mean_x + mean_y * mean_y)
     phi = np.arctan2(mean_y, mean_x)
 
@@ -152,17 +144,16 @@ def hillas_parameters_1(geom: CameraGeometry, image):
     # azwidth_2 = m_qq - m_q * m_q
     # azwidth = np.sqrt(azwidth_2)
 
-    return MomentParameters(size=size,
-                            cen_x=mean_x * unit,
-                            cen_y=mean_y * unit,
-                            length=length * unit,
-                            width=width * unit,
-                            r=r * unit,
-                            phi=Angle(phi * u.rad),
-                            psi=Angle(delta * u.rad),
-                            miss=miss * unit,
-                            skewness=skewness,
-                            kurtosis=kurtosis)
+    return HillasParametersContainer(
+        x=mean_x * unit, y=mean_y * unit,
+        r=r * unit, phi=Angle(phi * u.rad),
+        intensity=size,
+        length=length * unit,
+        width=width * unit,
+        psi=Angle(delta * u.rad),
+        skewness=skewness,
+        kurtosis=kurtosis
+    )
 
 
 def hillas_parameters_2(geom: CameraGeometry, image):
@@ -182,7 +173,7 @@ def hillas_parameters_2(geom: CameraGeometry, image):
 
     Returns
     -------
-    hillas_parameters : `MomentParameters`
+    hillas_parameters : `HillasParametersContainer`
     """
 
     if type(geom.pix_x) == Quantity:
@@ -208,12 +199,6 @@ def hillas_parameters_2(geom: CameraGeometry, image):
                                             "Cannot calculate image parameters."
                                             "Exiting...")))
 
-    pixdata = np.row_stack([pix_x,
-                            pix_y,
-                            pix_x * pix_x,
-                            pix_x * pix_y,
-                            pix_y * pix_y])
-
     # Compute image moments (done in a bit faster way, but putting all
     # into one 2D array, where each row will be summed to calculate a
     # moment) However, this doesn't avoid a temporary created for the
@@ -311,12 +296,15 @@ def hillas_parameters_2(geom: CameraGeometry, image):
                 vy4 * spsi2 * spsi2)
         kurtosis = kurt / (length * length * length * length)
 
-    return MomentParameters(size=size, cen_x=xm * unit, cen_y=ym * unit,
-                            length=length * unit, width=width * unit,
-                            r=rr * unit,
-                            phi=Angle(phi), psi=Angle(psi),
-                            miss=miss * unit, skewness=skewness,
-                            kurtosis=kurtosis)
+    return HillasParametersContainer(
+        x=xm * unit, y=ym * unit,
+        r=rr * unit, phi=Angle(phi),
+        intensity=size,
+        length=length * unit, width=width * unit,
+        psi=Angle(psi),
+        skewness=skewness,
+        kurtosis=kurtosis,
+    )
 
 
 def hillas_parameters_3(geom: CameraGeometry, image):
@@ -335,7 +323,7 @@ def hillas_parameters_3(geom: CameraGeometry, image):
 
     Returns
     -------
-    hillas_parameters : `MomentParameters`
+    hillas_parameters : `HillasParametersContainer`
     """
 
     if type(geom.pix_x) == Quantity:
@@ -483,15 +471,19 @@ def hillas_parameters_3(geom: CameraGeometry, image):
 
     # Code to de-interface with historical code
 
-    return MomentParameters(size=size, cen_x=m_x * unit, cen_y=m_y * unit,
-                            length=length * unit, width=width * unit,
-                            r=r * unit, phi=Angle(phi * u.rad),
-                            psi=Angle(psi * u.rad),
-                            miss=miss * unit,
-                            skewness=skewness, kurtosis=kurtosis)
+    return HillasParametersContainer(
+        x=m_x * unit, y=m_y * unit,
+        r=r * unit, phi=Angle(phi * u.rad),
+        intensity=size,
+        length=length * unit,
+        width=width * unit,
+        psi=Angle(psi * u.rad),
+        skewness=skewness,
+        kurtosis=kurtosis,
+    )
 
 
-def hillas_parameters_4(geom: CameraGeometry, image, container=False):
+def hillas_parameters_4(geom: CameraGeometry, image):
     """Compute Hillas parameters for a given shower image.
 
     As for hillas_parameters_3 (old Whipple Fortran code), but more Pythonized
@@ -511,16 +503,9 @@ def hillas_parameters_4(geom: CameraGeometry, image, container=False):
 
     Returns
     -------
-    MomentParameters:
-        tuple of hillas parameters
     HillasParametersContainer:
         container of hillas parametesr
     """
-
-    if not isinstance(geom, CameraGeometry):
-        raise ValueError("Hillas Parameters API has changed: hillas_parameters("
-                         "geom, image). Please update your code")
-
     unit = geom.pix_x.unit
 
     # MP: Actually, I don't know why we need to strip the units... shouldn'
@@ -641,23 +626,16 @@ def hillas_parameters_4(geom: CameraGeometry, image, container=False):
                 vy4 * spsi2 * spsi2)
         kurtosis = kurt / (length * length * length * length)
 
-    if container:
-        return HillasParametersContainer(x=m_x * unit, y=m_y * unit, r=r * unit,
-                                         phi=Angle(phi * u.rad),
-                                         intensity=size,
-                                         length=length * unit,
-                                         width=width * unit,
-                                         psi=Angle(psi * u.rad),
-                                         skewness=skewness,
-                                         kurtosis=kurtosis)
-    else:
-        return MomentParameters(size=size, cen_x=m_x * unit, cen_y=m_y * unit,
-                                length=length * unit, width=width * unit,
-                                r=r * unit,
-                                phi=Angle(phi * u.rad),
-                                psi=Angle(psi * u.rad),
-                                miss=miss * unit,
-                                skewness=skewness, kurtosis=kurtosis)
+    return HillasParametersContainer(
+        x=m_x * unit, y=m_y * unit,
+        r=r * unit, phi=Angle(phi * u.rad),
+        intensity=size,
+        length=length * unit,
+        width=width * unit,
+        psi=Angle(psi * u.rad),
+        skewness=skewness,
+        kurtosis=kurtosis,
+    )
 
 
 def hillas_parameters_5(geom: CameraGeometry, image):

ctapipe/image/tests/test_hillas.py

@@ -81,13 +81,11 @@ def test_hillas():
     for psi_angle in ['30d', '120d', '-30d', '-120d']:
 
         geom, image = create_sample_image_zeros(psi_angle)
-        results = {}
+        results = {
+            'v{}'.format(i): method(geom, image)
+            for i, method in enumerate(methods, start=1)
+        }
 
-        results['v1'] = hillas_parameters_1(geom, image)
-        results['v2'] = hillas_parameters_2(geom, image)
-        results['v3'] = hillas_parameters_3(geom, image)
-        results['v4'] = hillas_parameters_4(geom, image)
-        results['v5'] = hillas_parameters_5(geom, image)
         # compare each method's output
         for aa in results:
             for bb in results:
@@ -133,12 +131,13 @@ def test_hillas_failure():
 
 
 def test_hillas_api_change():
-    with pytest.raises(ValueError):
+    with pytest.raises(TypeError):
         hillas_parameters_4(arange(10), arange(10), arange(10))
 
 
 def test_hillas_container():
     geom, image = create_sample_image_zeros(psi='0d')
 
-    params = hillas_parameters_4(geom, image, container=True)
-    assert isinstance(params, HillasParametersContainer)
+    for method in methods:
+        params = method(geom, image)
+        assert isinstance(params, HillasParametersContainer)

ctapipe/plotting/array.py

@@ -24,8 +24,8 @@ def __init__(self, hillas_parameters, draw_axes=False, ax=None, **kwargs):
         """
         self.axes = ax if ax is not None else plt.gca()
 
-        self.cen_x = [i.cen_x.to(u.deg).value for i in hillas_parameters.values()]
-        self.cen_y = [i.cen_y.to(u.deg).value for i in hillas_parameters.values()]
+        self.cen_x = [i.x.to(u.deg).value for i in hillas_parameters.values()]
+        self.cen_y = [i.y.to(u.deg).value for i in hillas_parameters.values()]
 
         self.centre = (0, 0)
         self.array = ArrayDisplay(telx=np.asarray(self.cen_x),

ctapipe/reco/HillasReconstructor.py

@@ -87,7 +87,7 @@ class are set to np.nan
         -----------
         hillas_dict : python dictionary
             dictionary with telescope IDs as key and
-            MomentParameters instances as values
+            HillasParametersContainer instances as values
         inst : ctapipe.io.InstrumentContainer
             instrumental description
         pointing_alt:
@@ -136,7 +136,7 @@ class are set to np.nan
         result.core_uncert = err_est_pos
 
         result.tel_ids = [h for h in hillas_dict.keys()]
-        result.average_size = np.mean([h.size for h in hillas_dict.values()])
+        result.average_size = np.mean([h.intensity for h in hillas_dict.values()])
         result.is_valid = True
 
         result.alt_uncert = err_est_dir
@@ -174,8 +174,8 @@ def inititialize_hillas_planes(
 
         self.hillas_planes = {}
         for tel_id, moments in hillas_dict.items():
-            p2_x = moments.cen_x + 0.1 * u.m * np.cos(moments.psi)
-            p2_y = moments.cen_y + 0.1 * u.m * np.sin(moments.psi)
+            p2_x = moments.x + 0.1 * u.m * np.cos(moments.psi)
+            p2_y = moments.y + 0.1 * u.m * np.sin(moments.psi)
             focal_length = subarray.tel[tel_id].optics.equivalent_focal_length
 
             pointing = SkyCoord(
@@ -194,7 +194,7 @@ def inititialize_hillas_planes(
                 pointing_direction=pointing
             )
 
-            cog_coord = SkyCoord(x=moments.cen_x, y=moments.cen_y, frame=cf)
+            cog_coord = SkyCoord(x=moments.x, y=moments.y, frame=cf)
             cog_coord = cog_coord.transform_to(hf)
 
             p2_coord = SkyCoord(x=p2_x, y=p2_y, frame=cf)
@@ -204,7 +204,7 @@ def inititialize_hillas_planes(
                 p1=cog_coord,
                 p2=p2_coord,
                 telescope_position=subarray.positions[tel_id],
-                weight=moments.size * (moments.length / moments.width),
+                weight=moments.intensity * (moments.length / moments.width),
             )
             self.hillas_planes[tel_id] = circle
 
@@ -382,7 +382,7 @@ def estimate_h_max(self, hillas_dict, subarray, pointing_alt, pointing_az):
                 pointing_direction=pointing
             )
 
-            cog_coord = SkyCoord(x=moments.cen_x, y=moments.cen_y, frame=cf)
+            cog_coord = SkyCoord(x=moments.x, y=moments.y, frame=cf)
             cog_coord = cog_coord.transform_to(hf)
 
             cog_direction = spherical_to_cartesian(1, cog_coord.alt, cog_coord.az)

ctapipe/reco/ImPACT.py

@@ -28,7 +28,7 @@
 
 def guess_shower_depth(energy):
     """
-    Simple estimation of depth of shower max based on the expected gamma-ray elongation 
+    Simple estimation of depth of shower max based on the expected gamma-ray elongation
     rate.
 
     Parameters
@@ -192,9 +192,9 @@ def get_brightest_mean(self):
         tel_num = 0
         for tel in self.hillas:
 
-            weight = self.hillas[tel].size
-            weighted_x = self.hillas[tel].cen_x.to(u.rad).value * weight
-            weighted_y = self.hillas[tel].cen_y.to(u.rad).value * weight
+            weight = self.hillas[tel].intensity
+            weighted_x = self.hillas[tel].x.to(u.rad).value * weight
+            weighted_y = self.hillas[tel].y.to(u.rad).value * weight
 
             ppx = np.sum(weighted_x) / np.sum(weight)
             ppy = np.sum(weighted_y) / np.sum(weight)
@@ -595,7 +595,7 @@ def predict(self, shower_seed, energy_seed):
         Returns
         -------
         ReconstructedShowerContainer, ReconstructedEnergyContainer:
-        Reconstructed ImPACT shower geometry and energy        
+        Reconstructed ImPACT shower geometry and energy
         """
 
         horizon_seed = HorizonFrame(az=shower_seed.az, alt=shower_seed.alt)
@@ -751,14 +751,14 @@ def minimise(self, params, step, limits, minimiser_name="minuit"):
     def draw_nominal_surface(self, shower_seed, energy_seed, bins=30,
                              nominal_range=2.5 * u.deg):
         """
-        Simple reconstruction for evaluating the likelihood in a grid across the 
-        nominal system, fixing all values but the source position of the gamma rays. 
+        Simple reconstruction for evaluating the likelihood in a grid across the
+        nominal system, fixing all values but the source position of the gamma rays.
         Useful for checking the reconstruction performance of the algorithm
 
         Parameters
         ----------
         shower_seed: ReconstructedShowerContainer
-            Best fit ImPACT shower geometry 
+            Best fit ImPACT shower geometry
         energy_seed: ReconstructedEnergyContainer
             Best fit ImPACT energy
         bins: int
@@ -768,8 +768,8 @@ def draw_nominal_surface(self, shower_seed, energy_seed, bins=30,
 
         Returns
         -------
-        ndarray, ndarray, ndarray: 
-        Bin centres in X and Y coordinates and the values of the likelihood at each 
+        ndarray, ndarray, ndarray:
+        Bin centres in X and Y coordinates and the values of the likelihood at each
         position
         """
         horizon_seed = HorizonFrame(az=shower_seed.az, alt=shower_seed.alt)
@@ -825,7 +825,7 @@ def draw_tilted_surface(self, shower_seed, energy_seed,
         Parameters
         ----------
         shower_seed: ReconstructedShowerContainer
-            Best fit ImPACT shower geometry 
+            Best fit ImPACT shower geometry
         energy_seed: ReconstructedEnergyContainer
             Best fit ImPACT energy
         bins: int
@@ -835,7 +835,7 @@ def draw_tilted_surface(self, shower_seed, energy_seed,
 
         Returns
         -------
-        ndarray, ndarray, ndarray: 
+        ndarray, ndarray, ndarray:
             Bin centres in X and Y coordinates and the values of the likelihood
             at each position
         """