Update ImPACT code to work again

docs/changes/2305.feature.rst

@@ -0,0 +1,8 @@
+Large updates to the Image Pixel-wise fit for Atmospheric Cherenkov Telescopes reconstruction method (https://doi.org/10.48550/arXiv.1403.2993)
+
+* ImPACT - General code clean up and optimisation. Now updated to work similarly to other reconstructors using the standardised interface, such that it can be used ctapipe-process. Significant improvements to tests too
+* ImPACT_utilities - Created new file to hold general usage functions, numba used in some areas for speedup
+* template_network_interpolator - Now works with templates with different zenith and azimuth angles
+* unstructured_interpolator - Significant speed improvements
+* pixel_likelihood - Constants added back to neg_log_likelihood_approx, these are quite important to obtaining a well normalised goodness of fit.
+* hillas_intersection - Fixed bug in core position being incorrectly calculated, fixed tests too

src/ctapipe/image/muon/intensity_fitter.py

@@ -403,7 +403,7 @@ def negative_log_likelihood(
         # scale prediction by optical efficiency of the telescope
         prediction *= optical_efficiency_muon
 
-        return neg_log_likelihood_approx(image, prediction, spe_width, pedestal)
+        return neg_log_likelihood_approx(image, prediction, spe_width, pedestal).sum()
 
     return negative_log_likelihood
 

src/ctapipe/image/pixel_likelihood.py

@@ -89,15 +89,13 @@ def neg_log_likelihood_approx(image, prediction, spe_width, pedestal):
 
     Returns
     -------
-    float
+    ndarray
     """
-    theta = pedestal**2 + prediction * (1 + spe_width**2)
 
-    neg_log_l = 0.5 * (
-        np.log(2 * np.pi * theta + EPSILON) + (image - prediction) ** 2 / theta
-    )
+    theta = 2 * (pedestal**2 + prediction * (1 + spe_width**2))
+    neg_log_l = np.log(np.pi * theta) / 2.0 + (image - prediction) ** 2 / theta
 
-    return np.sum(neg_log_l)
+    return neg_log_l
 
 
 def neg_log_likelihood_numeric(
@@ -123,14 +121,14 @@ def neg_log_likelihood_numeric(
 
     Returns
     -------
-    float
+    ndarray
     """
 
     epsilon = np.finfo(np.float64).eps
 
     prediction = prediction + epsilon
 
-    likelihood = epsilon
+    likelihood = np.full_like(prediction, epsilon, dtype=np.float64)
 
     n_signal = np.arange(poisson(np.max(prediction)).ppf(confidence) + 1)
 
@@ -147,7 +145,7 @@ def neg_log_likelihood_numeric(
         )
         likelihood += _l
 
-    return -np.sum(np.log(likelihood + EPSILON))
+    return -np.log(likelihood)
 
 
 def neg_log_likelihood(image, prediction, spe_width, pedestal, prediction_safety=20.0):
@@ -173,27 +171,33 @@ def neg_log_likelihood(image, prediction, spe_width, pedestal, prediction_safety
 
     Returns
     -------
-    float
+    ndarray
     """
 
     approx_mask = prediction > prediction_safety
 
-    neg_log_l = 0
+    neg_log_l = np.zeros_like(image, dtype=np.float64)
     if np.any(approx_mask):
-        neg_log_l += neg_log_likelihood_approx(
-            image[approx_mask], prediction[approx_mask], spe_width, pedestal
+        neg_log_l[approx_mask] += neg_log_likelihood_approx(
+            image[approx_mask],
+            prediction[approx_mask],
+            spe_width[approx_mask],
+            pedestal[approx_mask],
         )
 
     if not np.all(approx_mask):
-        neg_log_l += neg_log_likelihood_numeric(
-            image[~approx_mask], prediction[~approx_mask], spe_width, pedestal
+        neg_log_l[~approx_mask] += neg_log_likelihood_numeric(
+            image[~approx_mask],
+            prediction[~approx_mask],
+            spe_width[~approx_mask],
+            pedestal[~approx_mask],
         )
 
     return neg_log_l
 
 
 def mean_poisson_likelihood_gaussian(prediction, spe_width, pedestal):
-    """Calculation of the mean likelihood for a give expectation
+    """Calculation of the mean of twice the negative log likelihood for a give expectation
     value of pixel intensity in the gaussian approximation.
     This is useful in the calculation of the goodness of fit.
 
@@ -208,12 +212,12 @@ def mean_poisson_likelihood_gaussian(prediction, spe_width, pedestal):
 
     Returns
     -------
-    float
+    ndarray
     """
     theta = pedestal**2 + prediction * (1 + spe_width**2)
     mean_log_likelihood = 1 + np.log(2 * np.pi) + np.log(theta + EPSILON)
 
-    return np.sum(mean_log_likelihood)
+    return mean_log_likelihood
 
 
 def _integral_poisson_likelihood_full(image, prediction, spe_width, ped):
@@ -229,7 +233,7 @@ def _integral_poisson_likelihood_full(image, prediction, spe_width, ped):
 
 def mean_poisson_likelihood_full(prediction, spe_width, ped):
     """
-    Calculation of the mean  likelihood for a give expectation value
+    Calculation of the mean of twice the negative log likelihood for a give expectation value
     of pixel intensity using the full numerical integration.
     This is useful in the calculation of the goodness of fit.
     This numerical integration is very slow and really doesn't
@@ -246,21 +250,21 @@ def mean_poisson_likelihood_full(prediction, spe_width, ped):
 
     Returns
     -------
-    float
+    ndarray
     """
 
     if len(spe_width) == 1:
-        spe_width = np.full_like(prediction, spe_width)
+        spe_width = np.full_like(prediction, spe_width, dtype=np.float64)
 
     if len(ped) == 1:
-        ped = np.full_like(prediction, ped)
+        ped = np.full_like(prediction, ped, dtype=np.float64)
 
-    mean_like = 0
+    mean_like = np.zeros_like(prediction, dtype=np.float64)
 
     width = ped**2 + prediction * spe_width**2
     width = np.sqrt(width)
 
-    for pred, w, spe, p in zip(prediction, width, spe_width, ped):
+    for i, (pred, w, spe, p) in enumerate(zip(prediction, width, spe_width, ped)):
         lower_integration_bound = pred - 10 * w
         upper_integration_bound = pred + 10 * w
 
@@ -272,7 +276,7 @@ def mean_poisson_likelihood_full(prediction, spe_width, ped):
             epsrel=0.05,
         )
 
-        mean_like += integral
+        mean_like[i] = integral
 
     return mean_like
 
@@ -294,10 +298,10 @@ def chi_squared(image, prediction, pedestal, error_factor=2.9):
 
     Returns
     -------
-    float
+    ndarray
     """
 
     chi_square = (image - prediction) ** 2 / (pedestal + 0.5 * (image - prediction))
     chi_square *= 1.0 / error_factor
 
-    return np.sum(chi_square)
+    return chi_square

src/ctapipe/image/tests/test_pixel_likelihood.py

@@ -20,7 +20,7 @@ def test_chi_squared():
     chi = chi_squared(image, prediction, ped)
     bad_chi = chi_squared(image, bad_prediction, ped)
 
-    assert chi < bad_chi
+    assert np.sum(chi) < np.sum(bad_chi)
 
 
 def test_mean_poisson_likelihoood_gaussian():
@@ -30,7 +30,7 @@ def test_mean_poisson_likelihoood_gaussian():
     small_mean_likelihood = mean_poisson_likelihood_gaussian(prediction, spe, 0)
     large_mean_likelihood = mean_poisson_likelihood_gaussian(prediction, spe, 1)
 
-    assert small_mean_likelihood < large_mean_likelihood
+    assert np.all(small_mean_likelihood < large_mean_likelihood)
 
     # Test that the mean likelihood of abunch of samples drawn from the gaussian
     # behind the approximate log likelihood is indeed the precalculated mean
@@ -48,7 +48,7 @@ def test_mean_poisson_likelihoood_gaussian():
     )
 
     rel_diff = (
-        2 * neg_log_likelihood_approx(normal_samples, prediction[0], spe, ped) / 100000
+        np.mean(2 * neg_log_likelihood_approx(normal_samples, prediction[0], spe, ped))
         - mean_likelihood
     ) / mean_likelihood
 
@@ -63,7 +63,7 @@ def test_mean_poisson_likelihood_full():
     small_mean_likelihood = mean_poisson_likelihood_full(prediction, spe, [0.1])
     large_mean_likelihood = mean_poisson_likelihood_full(prediction, spe, [1])
 
-    assert small_mean_likelihood < large_mean_likelihood
+    assert np.all(small_mean_likelihood < large_mean_likelihood)
 
 
 def test_full_likelihood():
@@ -72,30 +72,30 @@ def test_full_likelihood():
     signal cases. Check that full calculation and the gaussian approx become
     equal at high signal.
     """
-    spe = 0.5  # Single photo-electron width
-    pedestal = 1  # width of the pedestal distribution
+    spe = 0.5 * np.ones(3)  # Single photo-electron width
+    pedestal = np.ones(3)  # width of the pedestal distribution
 
     image_small = np.array([0, 1, 2])
     expectation_small = np.array([1, 1, 1])
 
     full_like_small = neg_log_likelihood(image_small, expectation_small, spe, pedestal)
-    exp_diff = full_like_small - np.sum(
-        np.asarray([1.37815294, 1.31084662, 1.69627197])
-    )
+    exp_rel_diff = (
+        full_like_small - np.asarray([1.37815294, 1.31084662, 1.69627197])
+    ) / full_like_small
 
     # Check against known values
-    assert np.abs(exp_diff / np.sum(full_like_small)) < 1e-4
+    assert np.all(np.abs(exp_rel_diff) < 3e-4)
 
     image_large = np.array([40, 50, 60])
     expectation_large = np.array([50, 50, 50])
 
     full_like_large = neg_log_likelihood(image_large, expectation_large, spe, pedestal)
     # Check against known values
-    exp_diff = full_like_large - np.sum(
-        np.asarray([3.72744569, 2.99652694, 3.83113004])
-    )
+    exp_rel_diff = (
+        full_like_large - np.asarray([3.78183004, 2.99452694, 3.78183004])
+    ) / full_like_large
 
-    assert np.abs(exp_diff / np.sum(full_like_large)) < 3e-4
+    assert np.all(np.abs(exp_rel_diff) < 3e-5)
 
     gaus_like_large = neg_log_likelihood_approx(
         image_large, expectation_large, spe, pedestal

src/ctapipe/reco/hillas_intersection.py

@@ -32,6 +32,7 @@
 from .reconstructor import (
     HillasGeometryReconstructor,
     InvalidWidthException,
+    ReconstructionProperty,
     TooFewTelescopesException,
 )
 
@@ -90,19 +91,17 @@ class HillasIntersection(HillasGeometryReconstructor):
     Note: only input from CameraFrame is currently supported
     """
 
-    atmosphere_profile_name = traits.CaselessStrEnum(
-        ["paranal"], default_value="paranal", help="name of atmosphere profile to use"
-    ).tag(config=True)
-
     weighting = traits.CaselessStrEnum(
         ["Konrad", "hess"], default_value="Konrad", help="Weighting Method name"
     ).tag(config=True)
 
-    def __init__(self, subarray, **kwargs):
+    property = ReconstructionProperty.GEOMETRY
+
+    def __init__(self, subarray, atmosphere_profile=None, **kwargs):
         """
         Weighting must be a function similar to the weight_konrad already implemented
         """
-        super().__init__(subarray, **kwargs)
+        super().__init__(subarray, atmosphere_profile, **kwargs)
 
         # We need a conversion function from height above ground to depth of maximum
         # To do this we need the conversion table from CORSIKA
@@ -423,8 +422,8 @@ def reconstruct_tilted(self, hillas_parameters, tel_x, tel_y):
         hillas2 = np.transpose(hillas2)
 
         # Perform intersection
-        crossing_x, crossing_y = self.intersect_lines(
-            tel_x[:, 0], tel_y[:, 0], hillas1[0], tel_x[:, 1], tel_y[:, 1], hillas2[0]
+        crossing_y, crossing_x = self.intersect_lines(
+            tel_y[:, 0], tel_x[:, 0], hillas1[0], tel_y[:, 1], tel_x[:, 1], hillas2[0]
         )
 
         # Weight by chosen method

src/ctapipe/reco/hillas_reconstructor.py

@@ -107,8 +107,12 @@ class that reconstructs the direction of an atmospheric shower
 
     """
 
-    def __init__(self, subarray: SubarrayDescription, **kwargs):
-        super().__init__(subarray=subarray, **kwargs)
+    def __init__(
+        self, subarray: SubarrayDescription, atmosphere_profile=None, **kwargs
+    ):
+        super().__init__(
+            subarray=subarray, atmosphere_profile=atmosphere_profile, **kwargs
+        )
         _cam_radius_m = {
             cam: cam.geometry.guess_radius().to_value(u.m)
             for cam in subarray.camera_types