Minimal initial commit of FreePACT code

.github/workflows/ci.yml

@@ -58,7 +58,7 @@ jobs:
             python-version: "3.12"
             install-method: pip
             extra-args: ["codecov"]
-            extras: tests,all
+            extras: tests,all,impact
 
           - name: Linux (3.13, mamba)
             os: ubuntu-24.04

docs/changes/2705.feature.rst

@@ -0,0 +1 @@
+Make template paths and likelihood parameters TelecopeParameters in ImPACT

pyproject.toml

@@ -17,54 +17,51 @@ license = "BSD-3-Clause"
 classifiers = [
     "Intended Audience :: Science/Research",
     "Programming Language :: Python :: 3",
+    "Programming Language :: Python :: 3.10",
+    "Programming Language :: Python :: 3.11",
     "Programming Language :: Python :: 3.12",
     "Programming Language :: Python :: 3.13",
-    "Programming Language :: Python :: 3.14",
     "Programming Language :: Python :: Implementation :: CPython",
     "Topic :: Scientific/Engineering :: Astronomy",
     "Development Status :: 3 - Alpha",
 ]
 
 dynamic = ["version"]
-requires-python = ">=3.12"
+requires-python = ">=3.10"
 
 dependencies = [
-    "astropy >=6.1,<8.0.0",
+    "astropy >=6.1,<8.0.0a0",
     "docutils",
     "joblib",
-    "numba >=0.59.0",
-    "numpy >=1.26,<3.0.0",
+    "numba >=0.57",
+    "numpy >=1.24,<3.0.0a0",
     "packaging",
     "psutil",
     "pyyaml >=5.1",
     "requests",
     "scikit-learn !=1.4.0", # 1.4.0 breaks with astropy tables, before and after works
-    "scipy ~=1.15",
+    "scipy ~=1.10",
     "tables ~=3.4",
     "tqdm >=4.32",
     "traitlets ~=5.6",
 ]
 
 [project.optional-dependencies]
 
-eventio = [
-    "eventio >=1.9.1,<3.0.0",
-]
-
 # all is with all optional *runtime* dependencies
 # use `dev` to get really all dependencies
 all = [
     "bokeh ~=3.0",
-    "ctapipe[eventio]",
+    "eventio >=1.9.1,<2.0.0a0",
     "iminuit >=2",
     "matplotlib ~=3.0",
-    "pyirf ~=0.13.0"
+    "pyirf ~=0.13.0",
 ]
 
 tests = [
     # at the moment, essentially all tests rely on test data from simtel
     # it doesn't make sense to skip all of these.
-    "ctapipe[eventio]",
+    "eventio >=1.9.1,<2.0.0a0",
     "h5py",
     "pandas",
     "pytest >= 7.0",
@@ -75,6 +72,9 @@ tests = [
     "astroquery",
 ]
 
+impact = ["tensorflow >=2.19.0",
+    "iminuit >=2"]
+
 docs = [
     "ctapipe[all]",
     "ffmpeg-python",
@@ -93,6 +93,7 @@ docs = [
     "sphinx_automodapi",
     "sphinxcontrib-bibtex",
     "sphinx-changelog",
+    "tomli; python_version < '3.11'",
 ]
 
 dev = [
@@ -161,19 +162,8 @@ filterwarnings = [
     "error::ctapipe.utils.deprecation.CTAPipeDeprecationWarning",
     "error::ctapipe.instrument.FromNameWarning",
     "error::ctapipe.core.traits.NoneDefaultNotAllowedWarning",
-    "error::tables.UnclosedFileWarning",
     "ignore:`np.MachAr` is deprecated:DeprecationWarning",
     "ignore::ctapipe.core.provenance.MissingReferenceMetadata",
-    # inside matplotlib, already fixed but occurs in "oldest-deps" tests
-    "ignore:'parseString':DeprecationWarning:matplotlib",
-    "ignore:'resetCache':DeprecationWarning:matplotlib",
-    "ignore:'oneOf':DeprecationWarning:matplotlib",
-    "ignore:'leaveWhiteSpace':DeprecationWarning:matplotlib",
-    "ignore:'setName':DeprecationWarning:matplotlib",
-    "ignore:'setParseAction':DeprecationWarning:matplotlib",
-    "ignore:'endQuoteChar':DeprecationWarning:matplotlib",
-    "ignore:'unquoteResults':DeprecationWarning:matplotlib",
-    "ignore:'parseAll':DeprecationWarning:pyparsing",
 ]
 norecursedirs = [
     ".git",
@@ -225,7 +215,7 @@ omit = ["src/ctapipe/_version.py"]
 
 
 [tool.ruff]
-target-version = "py312"
+target-version = "py310"
 line-length = 88
 
 [tool.ruff.lint]

src/ctapipe/reco/__init__.py

@@ -2,6 +2,7 @@
 
 # reconstructors must be imported before ShowerProcessor, so
 # they are available there
+from .freepact import FreePACTProtonReconstructor, FreePACTReconstructor
 from .hillas_intersection import HillasIntersection
 from .hillas_reconstructor import HillasReconstructor
 from .impact import ImPACTReconstructor
@@ -26,6 +27,8 @@
     "ShowerProcessor",
     "HillasReconstructor",
     "ImPACTReconstructor",
+    "FreePACTReconstructor",
+    "FreePACTProtonReconstructor",
     "HillasIntersection",
     "EnergyRegressor",
     "ParticleClassifier",

src/ctapipe/reco/freepact.py

@@ -0,0 +1,158 @@
+"""FreePACT Reconstructor for ctapipe
+This module implements the FreePACT reconstructor, which is a subclass of ImPACTReconstructor.
+It uses a neural network to predict th likelihood camera images based on the shower parameters."""
+
+import numpy as np
+import numpy.ma as ma
+
+from ctapipe.core import traits
+from ctapipe.core.telescope_component import TelescopeParameter
+from ctapipe.reco.impact import ImPACTReconstructor
+from ctapipe.utils.template_network_interpolator import FreePACTInterpolator
+
+from ..compat import COPY_IF_NEEDED
+from .impact_utilities import rotate_translate
+
+__all__ = ["FreePACTReconstructor"]
+
+
+class FreePACTReconstructor(ImPACTReconstructor):
+    """
+    FreePACT Reconstructor
+    This class implements the FreePACT reconstructor, which uses a neural network
+    to predict the likelihood of camera images based on shower parameters.
+    """
+
+    image_template_path = TelescopeParameter(
+        trait=traits.Path(exists=True, directory_ok=True, allow_none=False),
+        allow_none=False,
+        help=("Path to the image templates to be used in the reconstruction"),
+    ).tag(config=True)
+
+    def __init__(self, *args, **kwargs):
+        super().__init__(*args, **kwargs)
+
+    def set_up_templates(self):
+        """Set up the templates for the FreePACT reconstructor."""
+        template_sort_dict = {}
+
+        for tel_id in self.subarray.tel_ids:
+            if self.image_template_path.tel[tel_id] not in template_sort_dict.keys():
+                template_sort_dict[self.image_template_path.tel[tel_id]] = [tel_id]
+            else:
+                template_sort_dict[self.image_template_path.tel[tel_id]].append(tel_id)
+
+        for template_path, tel_ids in template_sort_dict.items():
+            net_interpolator = FreePACTInterpolator(template_path)
+
+            self.prediction[tuple(tel_ids)] = net_interpolator
+
+    def get_likelihood(
+        self,
+        source_x,
+        source_y,
+        core_x,
+        core_y,
+        energy,
+        x_max_scale,
+        goodness_of_fit=False,
+    ):
+        """Get the likelihood that the image predicted at the given test
+        position matches the camera image.
+
+        Parameters
+        ----------
+        source_x: float
+            Source position of shower in the nominal system (in deg)
+        source_y: float
+            Source position of shower in the nominal system (in deg)
+        core_x: float
+            Core position of shower in tilted telescope system (in m)
+        core_y: float
+            Core position of shower in tilted telescope system (in m)
+        energy: float
+            Shower energy (in TeV)
+        x_max_scale: float
+            Scaling factor applied to geometrically calculated Xmax
+
+        Returns
+        -------
+        float: Likelihood the model represents the camera image at this position
+
+        """
+        if np.isnan(source_x) or np.isnan(source_y):
+            return 1e8
+
+        # First we add units back onto everything.  Currently not
+        # handled very well, maybe in future we could just put
+        # everything in the correct units when loading in the class
+        # and ignore them from then on
+
+        zenith = self.zenith
+        azimuth = self.azimuth
+
+        # Geometrically calculate the depth of maximum given this test position
+        x_max_guess = self.get_shower_max(source_x, source_y, core_x, core_y, zenith)
+        x_max_guess *= x_max_scale
+
+        # Convert to binning of Xmax
+        x_max_diff = x_max_guess
+
+        # Calculate impact distance for all telescopes
+        impact = np.sqrt(
+            (self.tel_pos_x - core_x) ** 2 + (self.tel_pos_y - core_y) ** 2
+        )
+        # And the expected rotation angle
+        phi = np.arctan2(self.tel_pos_y - core_y, self.tel_pos_x - core_x)
+
+        # Rotate and translate all pixels such that they match the
+        # template orientation
+        # numba does not support masked arrays, work on underlying array and add mask back
+        pix_x_rot, pix_y_rot = rotate_translate(
+            self.pixel_y.data, self.pixel_x.data, source_y, source_x, -phi
+        )
+        pix_x_rot = np.ma.array(pix_x_rot, mask=self.pixel_x.mask, copy=COPY_IF_NEEDED)
+        pix_y_rot = np.ma.array(pix_y_rot, mask=self.pixel_y.mask, copy=COPY_IF_NEEDED)
+
+        # In the interpolator class we can gain speed advantages by using masked arrays
+        # so we need to make sure here everything is masked
+        likelihood = ma.zeros(self.image.shape)
+        likelihood.mask = ma.getmask(self.image)
+
+        for tel_ids, template in self.prediction.items():
+            template_mask = self.template_masks[tel_ids]
+            if np.any(template_mask):
+                likelihood[template_mask] = template(
+                    np.rad2deg(zenith),
+                    azimuth,
+                    energy * np.ones_like(impact[template_mask]),
+                    impact[template_mask],
+                    x_max_diff * np.ones_like(impact[template_mask]),
+                    np.rad2deg(pix_x_rot[template_mask]),
+                    np.rad2deg(pix_y_rot[template_mask]),
+                    self.image[template_mask],
+                )
+        if goodness_of_fit:
+            # return -2 * np.sum(likelihood[self.image>5]) / np.sum(self.image>5)
+            # print(ma.getmask(self.image))
+            return -2 * ma.sum(likelihood) / np.sum(~ma.getmask(self.image))
+
+        likelihood = ma.sum(likelihood) * -2
+
+        return likelihood
+
+
+class FreePACTProtonReconstructor(FreePACTReconstructor):
+    """FreePACT Proton Reconstructor
+    This class implements the FreePACT reconstructor for proton showers.
+    It uses a neural network to predict the likelihood of camera images based on shower parameters.
+    """
+
+    image_template_path = TelescopeParameter(
+        trait=traits.Path(exists=True, directory_ok=True, allow_none=False),
+        allow_none=False,
+        help=("Path to the image templates to be used in the reconstruction"),
+    ).tag(config=True)
+
+    def __init__(self, *args, **kwargs):
+        super().__init__(*args, **kwargs)