Merge pull request #124 from HealthyPear/upgrade-ctapipe_0.10.5

Support for ctapipe 0.10.5

.github/install.sh

@@ -7,7 +7,7 @@ conda update -q conda  # get latest conda version
 # Useful for debugging any issues with conda
 conda info -a
 
-sed -i -e "s/- python=.*/- python=$PYTHON_VERSION/g" environment.yml
+sed -i -e "s/- python=.*/- python=$PYTHON_VERSION/g" environment_development.yml
 conda install -c conda-forge mamba
-mamba env create -n protopipe --file environment.yml
+mamba env create -n protopipe --file environment_development.yml
 conda activate protopipe

docs/install/development.rst

@@ -4,11 +4,11 @@ Development version
 ===================
 
   1. `fork <https://help.github.com/en/articles/fork-a-repo>`__ the `repository <https://github.com/cta-observatory/protopipe>`_
-  2. create and enter a basic virtual environment (or use the ``environment.yaml`` file)
+  2. create a virtual environment (Anaconda users can use the ``environment_development.yml`` file)
   3. ``pip install -e '.[all]'``
 
   The ``all`` keyword will install all extra requirements,
-  which can be also installed separately using ``tests`` and ``docs``.
+  which can be also installed separately using ``tests`` and/or ``docs``.
 
 Next steps:
 

docs/install/release.rst

@@ -3,9 +3,13 @@
 Released version
 ================
 
-To install a released version >= ``0.4.0`` it is sufficient to install the
+To install the latest released version it is sufficient to install the
 package from ``PyPI`` with ``pip install protopipe``.
 
+If you prefer to work from an Anaconda virtual environment you can create it with,
+
+``conda env create -f environment_latest_release.yml``
+
 For previous releases,
 
 1. download the corresponding tarball stored `here <https://github.com/cta-observatory/protopipe/releases>`__

docs/pipeline/index.rst

@@ -106,4 +106,4 @@ Reference/API
 .. automodapi:: protopipe.pipeline
     :no-inheritance-diagram:
     :include-all-objects:
-    :skip: event_source
+    :skip: EventSource

environment.yml → environment_development.yml

@@ -1,12 +1,11 @@
-name: protopipe
+name: protopipe-dev
 channels:
   - defaults
-  - cta-observatory
 dependencies:
+  - python >=3.7,<3.9
   - pip
-  - ctapipe=0.9.1
+  - conda-forge::ctapipe=0.10.5
   - conda-forge::gammapy
-  - ctapipe-extra
   - astropy>=4.0.1
   - h5py=2
   - ipython
@@ -23,3 +22,4 @@ dependencies:
   - conda-forge::vitables
   - pip:
     - pyirf
+    - pytest-dependency

environment_latest_release.yml

@@ -0,0 +1,8 @@
+name: protopipe
+channels:
+  - defaults
+dependencies:
+  - python >=3.7,<3.9
+  - pip
+  - pip:
+    - protopipe

protopipe/pipeline/event_preparer.py

@@ -11,7 +11,7 @@
 from ctapipe.containers import ReconstructedShowerContainer
 from ctapipe.calib import CameraCalibrator
 from ctapipe.image.extractor import TwoPassWindowSum
-from ctapipe.image import leakage, number_of_islands, largest_island
+from ctapipe.image import leakage_parameters, number_of_islands, largest_island
 from ctapipe.utils.CutFlow import CutFlow
 from ctapipe.coordinates import GroundFrame, TelescopeFrame, CameraFrame
 
@@ -317,7 +317,7 @@ def prepare_event(self, source, return_stub=True, save_images=False, debug=False
                 )
                 print(
                     bcolors.BOLD
-                    + f"has triggered telescopes {event.r0.tels_with_data}"
+                    + f"has triggered telescopes {event.r1.tel.keys()}"
                     + bcolors.ENDC
                 )
                 ievt += 1
@@ -326,7 +326,7 @@ def prepare_event(self, source, return_stub=True, save_images=False, debug=False
 
             self.event_cutflow.count("noCuts")
 
-            if self.event_cutflow.cut("min2Tels trig", len(event.dl0.tels_with_data)):
+            if self.event_cutflow.cut("min2Tels trig", len(event.r1.tel.keys())):
                 if return_stub:
                     print(
                         bcolors.WARNING
@@ -361,7 +361,7 @@ def prepare_event(self, source, return_stub=True, save_images=False, debug=False
             hillas_dict = {}  # for discrimination
             leakage_dict = {}
             n_tels = {
-                "Triggered": len(event.dl0.tels_with_data),
+                "Triggered": len(event.r1.tel.keys()),
                 "LST_LST_LSTCam": 0,
                 "MST_MST_NectarCam": 0,
                 "MST_MST_FlashCam": 0,
@@ -378,16 +378,16 @@ def prepare_event(self, source, return_stub=True, save_images=False, debug=False
 
             good_for_reco = {}  # 1 = success, 0 = fail
 
-            # Compute impact parameter in tilt system
-            run_array_direction = event.mcheader.run_array_direction
-            az, alt = run_array_direction[0], run_array_direction[1]
+            # Array pointing in AltAz frame
+            az = event.pointing.array_azimuth
+            alt = event.pointing.array_altitude
 
             ground_frame = GroundFrame()
 
-            for tel_id in event.dl0.tels_with_data:
+            for tel_id in event.r1.tel.keys():
 
-                point_azimuth_dict[tel_id] = event.mc.tel[tel_id].azimuth_raw * u.rad
-                point_altitude_dict[tel_id] = event.mc.tel[tel_id].altitude_raw * u.rad
+                point_azimuth_dict[tel_id] = event.pointing.tel[tel_id].azimuth
+                point_altitude_dict[tel_id] = event.pointing.tel[tel_id].altitude
 
                 if debug:
                     print(
@@ -411,7 +411,7 @@ def prepare_event(self, source, return_stub=True, save_images=False, debug=False
                 if save_images is True:
                     # Save the simulated and reconstructed image of the event
                     dl1_phe_image[tel_id] = pmt_signal
-                    mc_phe_image[tel_id] = event.mc.tel[tel_id].true_image
+                    mc_phe_image[tel_id] = event.simulation.tel[tel_id].true_image
 
                 # We now ASSUME that the event will be good
                 good_for_reco[tel_id] = 1
@@ -429,7 +429,7 @@ def prepare_event(self, source, return_stub=True, save_images=False, debug=False
                     # The check on SIZE shouldn't be here, but for the moment
                     # I prefer to sacrifice elegancy...
                     if np.sum(image_biggest[mask_reco]) != 0.0:
-                        leakage_biggest = leakage(camera, image_biggest, mask_reco)
+                        leakage_biggest = leakage_parameters(camera, image_biggest, mask_reco)
                         leakages["leak1_reco"] = leakage_biggest["intensity_width_1"]
                         leakages["leak2_reco"] = leakage_biggest["intensity_width_2"]
                     else:
@@ -469,7 +469,7 @@ def prepare_event(self, source, return_stub=True, save_images=False, debug=False
                     # calculate the leakage (before filtering)
                     # this part is not well coded, but for the moment it works
                     if np.sum(image_extended[mask_extended]) != 0.0:
-                        leakage_extended = leakage(
+                        leakage_extended = leakage_parameters(
                             camera, image_extended, mask_extended
                         )
                         leakages["leak1"] = leakage_extended["intensity_width_1"]

protopipe/pipeline/utils.py

@@ -2,12 +2,13 @@
 import yaml
 import argparse
 import math
+import joblib
 
 import astropy.units as u
 import matplotlib.pyplot as plt
 import os.path as path
 
-from ctapipe.io import event_source
+from ctapipe.io import EventSource
 
 
 class bcolors:
@@ -237,7 +238,7 @@ def prod3b_array(fileName, site, array):
         This will feed both the estimators and the image cleaning.
 
     """
-    source = event_source(input_url=fileName, max_events=1)
+    source = EventSource(input_url=fileName, max_events=1)
 
     for event in source:  # get only first event
         pass
@@ -440,7 +441,7 @@ def CTAMARS_radii(camera_name):
         Name of the camera.
 
     Returns
-    ----------
+    -------
     average_camera_radii_deg : dict
         Dictionary containing the hard-coded values.
     """
@@ -480,3 +481,36 @@ def get_camera_names(inputPath=None):
     camera_names = [x.name for x in group._f_list_nodes()]
     h5file.close()
     return camera_names
+
+
+def load_models(path, cam_id_list):
+    """Load the pickled dictionary of model from disk
+    and fill the model dictionary.
+    
+    Parameters
+    ----------
+    path : string
+        The path where the pre-trained, pickled models are
+        stored. `path` is assumed to contain a `{cam_id}` keyword
+        to be replaced by each camera identifier in `cam_id_list`
+        (or at least a naked `{}`).
+    cam_id_list : list
+        List of camera identifiers like telescope ID or camera ID
+        and the assumed distinguishing feature in the filenames of
+        the various pickled regressors.
+    
+    Returns
+    -------
+    model_dict: dict
+        Dictionary with `cam_id` as keys and pickled models as values.
+    
+    """
+
+    model_dict = {}
+    for key in cam_id_list:
+            try:
+                model_dict[key] = joblib.load(path.format(cam_id=key))
+            except IndexError:
+                model_dict[key] = joblib.load(path.format(key))
+
+    return model_dict

protopipe/scripts/data_training.py

@@ -10,8 +10,7 @@
 import tables as tb
 
 from ctapipe.utils.CutFlow import CutFlow
-from ctapipe.io import event_source
-from ctapipe.reco.energy_regressor import EnergyRegressor
+from ctapipe.io import EventSource
 
 from protopipe.pipeline import EventPreparer
 from protopipe.pipeline.utils import (
@@ -21,6 +20,7 @@
     load_config,
     SignalHandler,
     bcolors,
+    load_models,
 )
 
 
@@ -116,7 +116,7 @@ def main():
             }
         )
 
-        regressor = EnergyRegressor.load(reg_file, cam_id_list=cams_and_foclens.keys())
+        regressors = load_models(reg_file, cam_id_list=cams_and_foclens.keys())
 
     # COLUMN DESCRIPTOR AS DICTIONARY
     # Column descriptor for the file containing output training data."""
@@ -207,7 +207,7 @@ def main():
 
         print("file: {} filename = {}".format(i, filename))
 
-        source = event_source(
+        source = EventSource(
             input_url=filename, allowed_tels=allowed_tels, max_events=args.max_events
         )
 
@@ -234,18 +234,15 @@ def main():
             source, save_images=args.save_images, debug=args.debug
         ):
 
-            # Angular quantities
-            run_array_direction = event.mcheader.run_array_direction
-
             if good_event:
 
                 xi = angular_separation(
-                    event.mc.az, event.mc.alt, reco_result.az, reco_result.alt
+                    event.simulation.shower.az, event.simulation.shower.alt, reco_result.az, reco_result.alt
                 )
 
                 offset = angular_separation(
-                    run_array_direction[0],  # az
-                    run_array_direction[1],  # alt
+                    event.pointing.array_azimuth,
+                    event.pointing.array_altitude,
                     reco_result.az,
                     reco_result.alt,
                 )
@@ -295,7 +292,7 @@ def main():
 
                     cam_id = source.subarray.tel[tel_id].camera.camera_name
                     moments = hillas_dict[tel_id]
-                    model = regressor.model_dict[cam_id]
+                    model = regressors[cam_id]
 
                     features_img = np.array(
                         [
@@ -377,7 +374,7 @@ def main():
                 outData[cam_id]["h_max"] = h_max.to("m").value
                 outData[cam_id]["err_est_pos"] = np.nan
                 outData[cam_id]["err_est_dir"] = np.nan
-                outData[cam_id]["true_energy"] = event.mc.energy.to("TeV").value
+                outData[cam_id]["true_energy"] = event.simulation.shower.energy.to("TeV").value
                 outData[cam_id]["hillas_x"] = moments.x.to("deg").value
                 outData[cam_id]["hillas_y"] = moments.y.to("deg").value
                 outData[cam_id]["hillas_phi"] = moments.phi.to("deg").value
@@ -392,13 +389,13 @@ def main():
                 outData[cam_id]["hillas_ellipticity"] = ellipticity.value
                 outData[cam_id]["clusters"] = n_cluster_dict[tel_id]
                 outData[cam_id]["n_tel_discri"] = n_tels["GOOD images"]
-                outData[cam_id]["mc_core_x"] = event.mc.core_x.to("m").value
-                outData[cam_id]["mc_core_y"] = event.mc.core_y.to("m").value
+                outData[cam_id]["mc_core_x"] = event.simulation.shower.core_x.to("m").value
+                outData[cam_id]["mc_core_y"] = event.simulation.shower.core_y.to("m").value
                 outData[cam_id]["reco_core_x"] = reco_core_x.to("m").value
                 outData[cam_id]["reco_core_y"] = reco_core_y.to("m").value
-                outData[cam_id]["mc_h_first_int"] = event.mc.h_first_int.to("m").value
+                outData[cam_id]["mc_h_first_int"] = event.simulation.shower.h_first_int.to("m").value
                 outData[cam_id]["offset"] = offset.to("deg").value
-                outData[cam_id]["mc_x_max"] = event.mc.x_max.value  # g / cm2
+                outData[cam_id]["mc_x_max"] = event.simulation.shower.x_max.value  # g / cm2
                 outData[cam_id]["alt"] = reco_result.alt.to("deg").value
                 outData[cam_id]["az"] = reco_result.az.to("deg").value
                 outData[cam_id]["reco_energy_tel"] = reco_energy_tel[tel_id]