Repository: cta-observatory/ctapipe
Pull request #1941 (2022-06-20 15:47:11)

Combining telescope-level dl2 predictions
Working on ways to combine the individual predictions to array-wide predictions.
Requires #1938 and could be improved with #1888 and #1935 merged.

This adds components, that get a table of monoscopic predictions, finds groups of unique `(obs_id, event_id)` pairs and
performs some operations on them. Currently this supports the mean and median, but everything available as `numpy.ufunc` is straigthforward to add.
The same is true for measurement functions in `scipy.ndimage`, which is probably even easier, but I already had the mean implemented using `np.add`.

My biggest issue right now is the handling of nans or, more generally, selecting rows to use for the calculation.
I did not find a solution I am really happy with, but it probably comes down to it being a usecase for quality queries on tables.
Maybe one could use the normal query in a loop to hide the issue from the user. I did not time it, but given that we spend most time in lower datalevels anyway, it might not be too bad...

review-notebook-app[bot]: Check out this pull request on&nbsp; <a href="https://app.reviewnb.com/cta-observatory/ctapipe/pull/1941"><img align="absmiddle"  alt="ReviewNB" height="28" class="BotMessageButtonImage" src="https://raw.githubusercontent.com/ReviewNB/support/master/images/button_reviewnb.png"/></a> 

 See visual diffs & provide feedback on Jupyter Notebooks. 

---

 <i>Powered by <a href='https://www.reviewnb.com/?utm_source=gh'>ReviewNB</a></i>
kosack: This is intended to operate on an existing DL2 file I guess.  It would be nice to also support event-wise prediction, as then we can include this in the list of Reconstructors that are run during `ctapipe-process` in the `ShowerProcessor`.  

Remember that the final goal is to have DL2 files with many reconstructions in them, so that in the DL2→DL3 step we can select which reconstruction to used on an event-by-event basis . 
maxnoe: @kosack yes, this is the idea. Have `__call__(event:ArrayEvent)` for process and `predict(table: Table)` for the `ctapipe-ml-apply` tool

Patch:
diff --git a/ctapipe/io/simteleventsource.py b/ctapipe/io/simteleventsource.py
index 8eb52425..8e064f66 100644
--- a/ctapipe/io/simteleventsource.py
+++ b/ctapipe/io/simteleventsource.py
@@ -20,6 +20,7 @@ from ..containers import (
     SimulatedShowerContainer,
     TelescopePointingContainer,
     TelescopeTriggerContainer,
+    DL2Container,
 )
 from ..coordinates import CameraFrame
 from ..core.traits import (
@@ -411,6 +412,7 @@ class SimTelEventSource(EventSource):
             data.r1.tel.clear()
             data.dl0.tel.clear()
             data.dl1.tel.clear()
+            data.dl2 = DL2Container()
             data.pointing.tel.clear()
             data.simulation.tel.clear()
             data.trigger.tel.clear()
diff --git a/ctapipe/io/tableloader.py b/ctapipe/io/tableloader.py
index f3ddc8cb..e2754fd9 100644
--- a/ctapipe/io/tableloader.py
+++ b/ctapipe/io/tableloader.py
@@ -173,7 +173,7 @@ class TableLoader(Component):
 
         self.instrument_table = None
         if self.load_instrument:
-            self.instrument_table = self.subarray.to_table('joined')
+            self.instrument_table = self.subarray.to_table("joined")
 
     def close(self):
         """Close the underlying hdf5 file"""
@@ -233,9 +233,12 @@ class TableLoader(Component):
 
                         # add the algorithm as prefix to distinguish multiple
                         # algorithms predicting the same quantities
-                        _add_column_prefix(
-                            dl2, prefix=algorithm, ignore=SUBARRAY_EVENT_KEYS
-                        )
+                        if group_name == "geometry":
+                            _add_column_prefix(
+                                dl2,
+                                prefix=algorithm,
+                                ignore=SUBARRAY_EVENT_KEYS,
+                            )
                         table = _join_subarray_events(table, dl2)
         return table
 
@@ -280,9 +283,12 @@ class TableLoader(Component):
 
                         # add the algorithm as prefix to distinguish multiple
                         # algorithms predicting the same quantities
-                        _add_column_prefix(
-                            dl2, prefix=algorithm, ignore=TELESCOPE_EVENT_KEYS
-                        )
+                        if group_name == "geometry":
+                            _add_column_prefix(
+                                dl2,
+                                prefix=algorithm,
+                                ignore=TELESCOPE_EVENT_KEYS,
+                            )
                         table = _join_telescope_events(table, dl2)
 
         if self.load_true_images:
@@ -315,7 +321,13 @@ class TableLoader(Component):
     def _join_subarray_info(self, table):
         subarray_events = self.read_subarray_events()
         table = join_allow_empty(
-            table, subarray_events, keys=SUBARRAY_EVENT_KEYS, join_type="left"
+            table,
+            subarray_events,
+            keys=SUBARRAY_EVENT_KEYS,
+            join_type="left",
+            # add suffix mono on duplicated columns, avoid underscore for stereo
+            table_names=["_mono", ""],
+            uniq_col_name="{col_name}{table_name}",
         )
         return table
 
diff --git a/ctapipe/ml/__init__.py b/ctapipe/ml/__init__.py
index 6a8beba7..5299750d 100644
--- a/ctapipe/ml/__init__.py
+++ b/ctapipe/ml/__init__.py
@@ -6,6 +6,7 @@ from .apply import (
     EnergyRegressor,
     ParticleIdClassifier,
 )
+from .stereo_combination import StereoCombiner, StereoMeanCombiner
 
 __all__ = [
     "ClassificationReconstructor",
@@ -14,4 +15,6 @@ __all__ = [
     "Reconstructor",
     "RegressionReconstructor",
     "TrainEnergyRegressor",
+    "StereoMeanCombiner",
+    "StereoCombiner",
 ]
diff --git a/ctapipe/ml/apply.py b/ctapipe/ml/apply.py
index b0efaabe..1d4487b5 100644
--- a/ctapipe/ml/apply.py
+++ b/ctapipe/ml/apply.py
@@ -15,11 +15,11 @@ from .sklearn import Classifier, Regressor, Model
 
 
 __all__ = [
-    'Reconstructor',
-    'ClassificationReconstructor',
-    'RegressionReconstructor',
-    'EnergyRegressor',
-    'ParticleIdClassifier',
+    "Reconstructor",
+    "ClassificationReconstructor",
+    "RegressionReconstructor",
+    "EnergyRegressor",
+    "ParticleIdClassifier",
 ]
 
 
@@ -34,7 +34,7 @@ class Reconstructor(Component):
     def __init__(self, subarray, *args, **kwargs):
         super().__init__(*args, **kwargs)
         self.subarray = subarray
-        self.instrument_table = self.subarray.to_table('joined')
+        self.instrument_table = self.subarray.to_table("joined")
 
     @classmethod
     def read(cls, path, subarray, *args, **kwargs):
@@ -82,12 +82,12 @@ class Reconstructor(Component):
             for algorithm, container in containers.items():
                 prefix = container.prefix
                 if prefix:
-                    container.prefix = f'{algorithm}_{prefix}'
+                    container.prefix = f"{algorithm}_{prefix}"
                 else:
                     container.prefix = algorithm
 
                 features.update(container.as_dict(add_prefix=True))
-                container.prefix = ''
+                container.prefix = ""
 
         features.update(self.instrument_table.loc[tel_id])
         return Table({k: [v] for k, v in features.items()})
@@ -95,12 +95,14 @@ class Reconstructor(Component):
 
 class RegressionReconstructor(Reconstructor):
     """Base class for sklearn regressors."""
+
     model_cls = Regressor
     model = Instance(model_cls).tag(config=True)
 
 
 class ClassificationReconstructor(Reconstructor):
     """Base class for sklearn regressors."""
+
     model_cls = Classifier
     model = Instance(model_cls).tag(config=True)
 
@@ -113,8 +115,8 @@ class EnergyRegressor(RegressionReconstructor):
             features = self._collect_features(event, tel_id)
             prediction, valid = self.model.predict(features)
             container = ReconstructedEnergyContainer(
-                energy=prediction,
-                is_valid=valid,
+                energy=prediction[0],
+                is_valid=valid[0],
             )
             event.dl2.tel[tel_id].energy[self.model.model_cls] = container
 
@@ -128,7 +130,7 @@ class ParticleIdClassifier(ClassificationReconstructor):
             prediction, valid = self.model.predict(features)
 
             container = ParticleClassificationContainer(
-                prediction=prediction,
-                is_valid=valid,
+                prediction=prediction[0],
+                is_valid=valid[0],
             )
             event.dl2.tel[tel_id].classification[self.model.model_cls] = container
diff --git a/ctapipe/ml/stereo_combination.py b/ctapipe/ml/stereo_combination.py
new file mode 100644
index 00000000..6d12c0fb
--- /dev/null
+++ b/ctapipe/ml/stereo_combination.py
@@ -0,0 +1,238 @@
+from abc import abstractmethod
+import numpy as np
+from astropy.table import Table
+import astropy.units as u
+from ctapipe.core import Component, Container
+from ctapipe.core.traits import CaselessStrEnum, Unicode
+from ..containers import (
+    ArrayEventContainer,
+    ParticleClassificationContainer,
+    ReconstructedEnergyContainer,
+)
+
+
+def _calculate_ufunc_of_telescope_values(tel_data, n_array_events, indices, ufunc):
+    combined_values = np.zeros(n_array_events)
+    ufunc.at(combined_values, indices, tel_data)
+    return combined_values
+
+
+def _weighted_mean_ufunc(tel_values, weights, n_array_events, indices):
+    # avoid numerical problems by very large or small weights
+    weights = weights / weights.max()
+    sum_prediction = _calculate_ufunc_of_telescope_values(
+        tel_values * weights,
+        n_array_events,
+        indices,
+        np.add,
+    )
+    sum_of_weights = _calculate_ufunc_of_telescope_values(
+        weights, n_array_events, indices, np.add
+    )
+    mean = np.full(n_array_events, np.nan)
+    valid = sum_of_weights > 0
+    mean[valid] = sum_prediction[valid] / sum_of_weights[valid]
+    return mean
+
+
+class StereoCombiner(Component):
+    # TODO: Add quality query (after #1888)
+    algorithm = Unicode().tag(config=True)
+    combine_property = CaselessStrEnum(["energy", "classification", "direction"]).tag(
+        config=True
+    )
+
+    @abstractmethod
+    def __call__(self, event: ArrayEventContainer) -> None:
+        """
+        Fill event container with stereo predictions
+        """
+
+    @abstractmethod
+    def predict(self, mono_predictions: Table) -> Table:
+        """
+        Constructs stereo predictions from a table of
+        telescope events.
+        """
+
+
+class StereoMeanCombiner(StereoCombiner):
+    """
+    Calculates array-wide (stereo) predictions as the mean of
+    the reconstruction on telescope-level with an optional weighting.
+    """
+
+    weights = CaselessStrEnum(
+        ["none", "intensity", "konrad"],
+        default_value="none",
+    ).tag(config=True)
+
+    def _calculate_weights(self, data):
+        """"""
+
+        if isinstance(data, Container):
+            if self.weights == "intensity":
+                return data.hillas.intensity
+
+            if self.weights == "konrad":
+                return data.hillas.intensity * data.hillas.length / data.hillas.width
+
+            return 1
+
+        if isinstance(data, Table):
+            if self.weights == "intensity":
+                return data["hillas_intensity"]
+
+            if self.weights == "konrad":
+                return (
+                    data["hillas_intensity"]
+                    * data["hillas_length"]
+                    / data["hillas_width"]
+                )
+
+            return np.ones(len(data))
+
+        raise TypeError(
+            "Dl1 data needs to be provided in the form of a container or astropy.table.Table"
+        )
+
+    def _combine_energy(self, event):
+        ids = []
+        values = []
+        weights = []
+
+        for tel_id, dl2 in event.dl2.tel.items():
+            mono = dl2.energy[self.algorithm]
+            if mono.is_valid:
+                values.append(mono.energy.to_value(u.TeV))
+                if tel_id not in event.dl1.tel:
+                    raise ValueError("No parameters for weighting available")
+                weights.append(
+                    self._calculate_weights(event.dl1.tel[tel_id].parameters)
+                )
+                ids.append(tel_id)
+
+        if len(values) > 0:
+            weights = np.array(weights, dtype=np.float64)
+            weights /= weights.max()
+            mean = np.average(values, weights=weights)
+            std = np.sqrt(np.cov(values, aweights=weights))
+            valid = True
+        else:
+            mean = std = np.nan
+            valid = False
+
+        event.dl2.stereo.energy[self.algorithm] = ReconstructedEnergyContainer(
+            energy=u.Quantity(mean, u.TeV, copy=False),
+            energy_uncert=u.Quantity(std, u.TeV, copy=False),
+            tel_ids=ids,
+            is_valid=valid,
+        )
+        event.dl2.stereo.energy[self.algorithm].prefix = self.algorithm
+
+    def _combine_classification(self, event):
+        ids = []
+        values = []
+        weights = []
+
+        for tel_id, dl2 in event.dl2.tel.items():
+            mono = dl2.classification[self.algorithm]
+            if mono.is_valid:
+                values.append(mono.prediction)
+                dl1 = event.dl1.tel[tel_id].parameters
+                weights.append(self._calculate_weights(dl1) if dl1 else 1)
+                ids.append(tel_id)
+
+        if len(values) > 0:
+            mean = np.average(values, weights=weights)
+            valid = True
+        else:
+            mean = np.nan
+            valid = False
+
+        container = ParticleClassificationContainer(
+            prediction=mean, tel_ids=ids, is_valid=valid
+        )
+        container.prefix = self.algorithm
+        event.dl2.stereo.classification[self.algorithm] = container
+
+    def __call__(self, event: ArrayEventContainer) -> None:
+        """
+        Calculate the mean prediction for a single array event.
+        """
+        if self.combine_property == "energy":
+            self._combine_energy(event)
+        elif self.combine_property == "classification":
+            self._combine_classification(event)
+        else:
+            raise NotImplementedError(f"Cannot combine {self.combine_property}")
+
+    def predict(self, mono_predictions: Table) -> Table:
+        """
+        Calculates the (array-)event-wise mean.
+        Telescope events, that are nan, get discarded.
+        This means you might end up with less events if
+        all telescope predictions of a shower are invalid.
+        """
+
+        prefix = self.algorithm
+        # TODO: Integrate table quality query once its done
+        valid = mono_predictions[f"{prefix}_is_valid"]
+        valid_predictions = mono_predictions[valid]
+
+        array_events, indices, multiplicity = np.unique(
+            mono_predictions[["obs_id", "event_id"]],
+            return_inverse=True,
+            return_counts=True,
+        )
+        stereo_table = Table(array_events)
+        n_array_events = len(array_events)
+        weights = self._calculate_weights(valid_predictions)
+
+        if self.combine_property == "classification":
+            mono_predictions = valid_predictions[f"{prefix}_prediction"]
+            stereo_predictions = _weighted_mean_ufunc(
+                mono_predictions, weights, n_array_events, indices[valid]
+            )
+            stereo_table[f"{prefix}_prediction"] = stereo_predictions
+            stereo_table[f"{prefix}_is_valid"] = np.isfinite(stereo_predictions)
+            stereo_table[f"{prefix}_goodness_of_fit"] = np.nan
+
+        elif self.combine_property == "energy":
+            mono_energies = valid_predictions[f"{prefix}_energy"].quantity.to_value(
+                u.TeV
+            )
+            stereo_energy = _weighted_mean_ufunc(
+                mono_energies,
+                weights,
+                n_array_events,
+                indices[valid],
+            )
+
+            stereo_table[f"{prefix}_energy"] = u.Quantity(
+                stereo_energy, u.TeV, copy=False
+            )
+
+            variance = _weighted_mean_ufunc(
+                (mono_energies - np.repeat(stereo_energy, multiplicity)[valid]) ** 2,
+                weights,
+                n_array_events,
+                indices[valid],
+            )
+            stereo_table[f"{prefix}_energy_uncert"] = u.Quantity(
+                np.sqrt(variance), u.TeV, copy=False
+            )
+            stereo_table[f"{prefix}_is_valid"] = np.isfinite(stereo_energy)
+            stereo_table[f"{prefix}_goodness_of_fit"] = np.nan
+
+        else:
+            raise NotImplementedError()
+
+        tel_ids = [[] for _ in range(n_array_events)]
+
+        for index, tel_id in zip(indices[valid], valid_predictions["tel_id"]):
+            tel_ids[index].append(tel_id)
+
+        k = f"{prefix}_tel_ids"
+        stereo_table[k] = tel_ids
+        return stereo_table
diff --git a/ctapipe/ml/tools/apply_energy_regressor.py b/ctapipe/ml/tools/apply_energy_regressor.py
index 80c28ab3..1a801124 100644
--- a/ctapipe/ml/tools/apply_energy_regressor.py
+++ b/ctapipe/ml/tools/apply_energy_regressor.py
@@ -1,12 +1,15 @@
-from astropy.table import Table
+from astropy.table.operations import vstack
 import tables
 from ctapipe.core.tool import Tool
-from ctapipe.core.traits import Bool, Path, flag
+from ctapipe.core.traits import Bool, Path, flag, create_class_enum_trait
 from ctapipe.io import TableLoader, write_table
 from tqdm.auto import tqdm
+from ctapipe.io.tableio import TelListToMaskTransform
+import numpy as np
 
 from ..sklearn import Regressor
 from ..apply import EnergyRegressor
+from ..stereo_combination import StereoCombiner
 
 
 class ApplyEnergyRegressor(Tool):
@@ -19,10 +22,15 @@ class ApplyEnergyRegressor(Tool):
         directory_ok=False,
         exists=True,
     ).tag(config=True)
+
     model_path = Path(
         default_value=None, allow_none=False, exists=True, directory_ok=False
     ).tag(config=True)
 
+    stereo_combiner_type = create_class_enum_trait(
+        base_class=StereoCombiner, default_value="StereoMeanCombiner"
+    ).tag(config=True)
+
     aliases = {
         ("i", "input"): "ApplyEnergyRegressor.input_url",
         ("m", "model"): "ApplyEnergyRegressor.model_path",
@@ -44,11 +52,11 @@ class ApplyEnergyRegressor(Tool):
     classes = [
         TableLoader,
         Regressor,
+        StereoCombiner,
     ]
 
     def setup(self):
         """"""
-        self.estimator = EnergyRegressor.read(self.model_path, parent=self)
         self.h5file = tables.open_file(self.input_url, mode="r+")
         self.loader = TableLoader(
             parent=self,
@@ -59,30 +67,59 @@ class ApplyEnergyRegressor(Tool):
             load_simulated=True,
             load_instrument=True,
         )
+        self.estimator = EnergyRegressor.read(
+            self.model_path,
+            self.loader.subarray,
+            parent=self,
+        )
+        self.combine = StereoCombiner.from_name(
+            self.stereo_combiner_type,
+            combine_property="energy",
+            algorithm=self.estimator.model.model_cls,
+            parent=self,
+        )
 
     def start(self):
         self.log.info("Applying model")
 
+        tables = []
         for tel_id in tqdm(self.loader.subarray.tel):
             table = self.loader.read_telescope_events([tel_id])
+            if len(table) == 0:
+                continue
 
             prediction, valid = self.estimator.predict(table)
-            table = Table(
-                {
-                    "obs_id": table["obs_id"],
-                    "event_id": table["event_id"],
-                    "tel_id": table["tel_id"],
-                    "reconstructed_energy_energy": prediction,
-                    "reconstructed_energy_is_valid": valid,
-                }
-            )
+            prefix = self.estimator.model.model_cls
+
+            energy_col = f"{prefix}_energy"
+            valid_col = f"{prefix}_is_valid"
+            table[energy_col] = prediction
+            table[valid_col] = valid
+
             write_table(
-                table,
+                table[["obs_id", "event_id", "tel_id", energy_col, valid_col]],
                 self.loader.input_url,
-                f"/dl2/event/telescope/energy/{self.estimator.model.model_cls}/tel_{tel_id:03d}",
+                f"/dl2/event/telescope/energy/{prefix}/tel_{tel_id:03d}",
                 mode="a",
                 overwrite=self.overwrite,
             )
+            tables.append(table)
+
+        mono_predictions = vstack(tables)
+        stereo_predictions = self.combine.predict(mono_predictions)
+        trafo = TelListToMaskTransform(self.loader.subarray)
+        for c in filter(
+            lambda c: c.name.endswith("tel_ids"), stereo_predictions.columns.values()
+        ):
+            stereo_predictions[c.name] = np.array([trafo(r) for r in c])
+
+        write_table(
+            stereo_predictions,
+            self.loader.input_url,
+            f"/dl2/event/subarray/energy/{self.estimator.model.model_cls}",
+            mode="a",
+            overwrite=self.overwrite,
+        )
 
     def finish(self):
         self.h5file.close()
diff --git a/ctapipe/tools/process.py b/ctapipe/tools/process.py
index 1b562b6b..534fe63d 100644
--- a/ctapipe/tools/process.py
+++ b/ctapipe/tools/process.py
@@ -8,7 +8,14 @@ from tqdm.auto import tqdm
 
 from ..calib import CameraCalibrator, GainSelector
 from ..core import QualityQuery, Tool
-from ..core.traits import Bool, classes_with_traits, flag, Path
+from ..core.traits import (
+    Bool,
+    Dict,
+    List,
+    classes_with_traits,
+    flag,
+    Path,
+)
 from ..image import ImageCleaner, ImageModifier, ImageProcessor
 from ..image.extractor import ImageExtractor
 from ..io import DataLevel, DataWriter, EventSource, SimTelEventSource, write_table
@@ -16,7 +23,7 @@ from ..io.datawriter import DATA_MODEL_VERSION
 from ..reco import ShowerProcessor
 from ..utils import EventTypeFilter
 from ..io import metadata
-from ..ml import EnergyRegressor
+from ..ml import EnergyRegressor, StereoCombiner
 
 COMPATIBLE_DATALEVELS = [
     DataLevel.R1,
@@ -59,8 +66,11 @@ class ProcessorTool(Tool):
     ).tag(config=True)
 
     energy_regressor_path = Path(
-        default_value=None, allow_none=True, exists=True, directory_ok=False,
-        help='Path to a trained energy regression model (see ctapipe-ml-train-energy-regressor)'
+        default_value=None,
+        allow_none=True,
+        exists=True,
+        directory_ok=False,
+        help="Path to a trained energy regression model (see ctapipe-ml-train-energy-regressor)",
     ).tag(config=True)
 
     force_recompute_dl2 = Bool(
@@ -68,6 +78,8 @@ class ProcessorTool(Tool):
         default_value=False,
     ).tag(config=True)
 
+    stereo_combiner_configs = List(Dict()).tag(config=True)
+
     aliases = {
         ("i", "input"): "EventSource.input_url",
         ("o", "output"): "DataWriter.output_path",
@@ -193,6 +205,14 @@ class ProcessorTool(Tool):
                 parent=self,
             )
 
+        self.stereo_combiners = []
+        for stereo_combiner in self.stereo_combiner_configs:
+            cfg = stereo_combiner.copy()
+            name = cfg.pop("type", "StereoMeanCombiner")
+            self.stereo_combiners.append(
+                StereoCombiner.from_name(name, **cfg, parent=self)
+            )
+
         # warn if max_events prevents writing the histograms
         if (
             isinstance(self.event_source, SimTelEventSource)
@@ -297,6 +317,9 @@ class ProcessorTool(Tool):
                 if self.energy_regressor is not None:
                     self.energy_regressor(event)
 
+            for combiner in self.stereo_combiners:
+                combiner(event)
+
             self.write(event)
 
     def finish(self):


Test patch:
diff --git a/ctapipe/ml/tests/test_stereo_combination.py b/ctapipe/ml/tests/test_stereo_combination.py
new file mode 100644
index 00000000..64ea6c2d
--- /dev/null
+++ b/ctapipe/ml/tests/test_stereo_combination.py
@@ -0,0 +1,161 @@
+import pytest
+from astropy.table import Table
+import astropy.units as u
+import numpy as np
+from numpy.testing import assert_array_equal
+from ctapipe.containers import (
+    HillasParametersContainer,
+    ImageParametersContainer,
+    ReconstructedEnergyContainer,
+    ParticleClassificationContainer,
+    ArrayEventContainer,
+    ReconstructedContainer,
+)
+from ctapipe.ml.stereo_combination import StereoMeanCombiner
+
+
+@pytest.fixture(scope="module")
+def mono_table():
+    """
+    Dummy table of telescope events with a
+    prediction and weights.
+    """
+    return Table(
+        {
+            "obs_id": [1, 1, 1, 1, 1, 2],
+            "event_id": [1, 1, 1, 2, 2, 1],
+            "tel_id": [1, 2, 3, 5, 7, 1],
+            "hillas_intensity": [1, 2, 0, 1, 5, 9],
+            "hillas_width": [0.1, 0.2, 0.1, 0.1, 0.2, 0.1] * u.deg,
+            "hillas_length": 3 * ([0.1, 0.2, 0.1, 0.1, 0.2, 0.1] * u.deg),
+            "dummy_energy": [1, 10, 4, 0.5, 0.7, 1] * u.TeV,
+            "dummy_is_valid": [
+                True,
+                True,
+                True,
+                True,
+                False,
+                False,
+            ],
+            "classifier_prediction": [1, 0, 0.5, 0, 0.6, 1],
+            "classifier_is_valid": [
+                True,
+                True,
+                False,
+                True,
+                True,
+                True,
+            ],
+        }
+    )
+
+
+@pytest.mark.parametrize("weights", ["konrad", "intensity", "none"])
+def test_predict_mean_energy(weights, mono_table):
+    combine = StereoMeanCombiner(
+        algorithm="dummy",
+        combine_property="energy",
+        weights=weights,
+    )
+    stereo = combine.predict(mono_table)
+    assert stereo.colnames == [
+        "obs_id",
+        "event_id",
+        "dummy_energy",
+        "dummy_energy_uncert",
+        "dummy_is_valid",
+        "dummy_goodness_of_fit",
+        "dummy_tel_ids",
+    ]
+    assert_array_equal(stereo["obs_id"], np.array([1, 1, 2]))
+    assert_array_equal(stereo["event_id"], np.array([1, 2, 1]))
+    if weights == "intensity":
+        assert_array_equal(stereo["dummy_energy"], [7, 0.5, np.nan] * u.TeV)
+    elif weights == "none":
+        assert_array_equal(stereo["dummy_energy"], [5, 0.5, np.nan] * u.TeV)
+
+    assert_array_equal(stereo["dummy_tel_ids"][0], np.array([1, 2, 3]))
+    assert_array_equal(stereo["dummy_tel_ids"][1], 5)
+
+
+def test_predict_mean_classification(mono_table):
+    combine = StereoMeanCombiner(
+        algorithm="classifier",
+        combine_property="classification",
+    )
+    stereo = combine.predict(mono_table)
+    assert stereo.colnames == [
+        "obs_id",
+        "event_id",
+        "classifier_prediction",
+        "classifier_is_valid",
+        "classifier_goodness_of_fit",
+        "classifier_tel_ids",
+    ]
+    assert_array_equal(stereo["obs_id"], np.array([1, 1, 2]))
+    assert_array_equal(stereo["event_id"], np.array([1, 2, 1]))
+    assert_array_equal(
+        stereo["classifier_prediction"],
+        [0.5, 0.3, 1],
+    )
+    tel_ids = stereo["classifier_tel_ids"]
+    assert_array_equal(tel_ids[0], [1, 2])
+    assert_array_equal(tel_ids[1], [5, 7])
+    assert_array_equal(tel_ids[2], [1])
+
+
+@pytest.mark.parametrize("weights", ["konrad", "intensity", "none"])
+def test_mean_prediction_single_event(weights):
+    event = ArrayEventContainer()
+
+    for tel_id, intensity in zip((25, 125, 130), (100, 200, 400)):
+        event.dl1.tel[tel_id].parameters = ImageParametersContainer(
+            hillas=HillasParametersContainer(
+                intensity=intensity,
+                width=0.1 * u.deg,
+                length=0.3 * u.deg,
+            )
+        )
+
+    event.dl2.tel[25] = ReconstructedContainer(
+        energy={
+            "dummy": ReconstructedEnergyContainer(energy=10 * u.GeV, is_valid=True)
+        },
+        classification={
+            "dummy": ParticleClassificationContainer(prediction=1.0, is_valid=True)
+        },
+    )
+    event.dl2.tel[125] = ReconstructedContainer(
+        energy={
+            "dummy": ReconstructedEnergyContainer(energy=20 * u.GeV, is_valid=True)
+        },
+        classification={
+            "dummy": ParticleClassificationContainer(prediction=0.0, is_valid=True)
+        },
+    )
+    event.dl2.tel[130] = ReconstructedContainer(
+        energy={
+            "dummy": ReconstructedEnergyContainer(energy=0.04 * u.TeV, is_valid=True)
+        },
+        classification={
+            "dummy": ParticleClassificationContainer(prediction=0.8, is_valid=True)
+        },
+    )
+
+    combine_energy = StereoMeanCombiner(
+        algorithm="dummy",
+        combine_property="energy",
+        weights=weights,
+    )
+    combine_classification = StereoMeanCombiner(
+        algorithm="dummy",
+        combine_property="classification",
+        weights=weights,
+    )
+    combine_energy(event)
+    combine_classification(event)
+    if weights == "none":
+        assert u.isclose(event.dl2.stereo.energy["dummy"].energy, (70 / 3) * u.GeV)
+    elif weights == "intensity":
+        assert u.isclose(event.dl2.stereo.energy["dummy"].energy, 30 * u.GeV)
+    assert event.dl2.stereo.classification["dummy"].prediction == 0.6
diff --git a/ctapipe/ml/tools/tests/test_apply.py b/ctapipe/ml/tools/tests/test_apply.py
new file mode 100644
index 00000000..63c51c59
--- /dev/null
+++ b/ctapipe/ml/tools/tests/test_apply.py
@@ -0,0 +1,40 @@
+from ctapipe.core import run_tool
+from ctapipe.io import TableLoader
+import shutil
+
+
+def test_apply_energy_regressor(
+    energy_regressor_path, dl2_shower_geometry_file, tmp_path
+):
+    from ctapipe.ml.tools.apply_energy_regressor import ApplyEnergyRegressor
+
+    input_path = tmp_path / dl2_shower_geometry_file.name
+
+    # create copy to not mutate common test file
+    shutil.copy2(dl2_shower_geometry_file, input_path)
+
+    ret = run_tool(
+        ApplyEnergyRegressor(),
+        argv=[
+            f"--input={input_path}",
+            f"--model={energy_regressor_path}",
+            "--ApplyEnergyRegressor.StereoMeanCombiner.weights=konrad",
+        ],
+    )
+    assert ret == 0
+
+    loader = TableLoader(input_path, load_dl2=True)
+    events = loader.read_subarray_events()
+    assert "ExtraTreesRegressor_energy" in events.colnames
+    assert "ExtraTreesRegressor_energy_uncert" in events.colnames
+    assert "ExtraTreesRegressor_is_valid" in events.colnames
+    assert "ExtraTreesRegressor_tel_ids" in events.colnames
+
+    events = loader.read_telescope_events()
+    assert "ExtraTreesRegressor_energy" in events.colnames
+    assert "ExtraTreesRegressor_energy_uncert" in events.colnames
+    assert "ExtraTreesRegressor_is_valid" in events.colnames
+    assert "ExtraTreesRegressor_tel_ids" in events.colnames
+
+    assert "ExtraTreesRegressor_energy_mono" in events.colnames
+    assert "ExtraTreesRegressor_is_valid_mono" in events.colnames
diff --git a/ctapipe/ml/tools/tests/test_process.py b/ctapipe/ml/tools/tests/test_process.py
index d5f630c3..82d845aa 100644
--- a/ctapipe/ml/tools/tests/test_process.py
+++ b/ctapipe/ml/tools/tests/test_process.py
@@ -1,5 +1,6 @@
 from ctapipe.core import run_tool
 from ctapipe.io import read_table
+import json
 
 
 def test_process_apply_energy(tmp_path, energy_regressor_path):
@@ -7,16 +8,34 @@ def test_process_apply_energy(tmp_path, energy_regressor_path):
 
     output = tmp_path / "gamma_prod5.dl2_energy.h5"
 
+    config_path = tmp_path / "config.json"
+
+    config = {
+        "ProcessorTool": {
+            "stereo_combiner_configs": [
+                {
+                    "type": "StereoMeanCombiner",
+                    "combine_property": "energy",
+                    "algorithm": "ExtraTreesRegressor",
+                    "weights": "konrad",
+                }
+            ]
+        }
+    }
+
+    with config_path.open("w") as f:
+        json.dump(config, f)
+
     argv = [
         "--input=dataset://gamma_prod5.simtel.zst",
         f"--output={output}",
         "--write-images",
         "--write-stereo-shower",
         "--write-mono-shower",
-        f"--energy-regressor={energy_regressor_path}"
-
+        f"--energy-regressor={energy_regressor_path}",
+        f"--config={config_path}",
     ]
     assert run_tool(ProcessorTool(), argv=argv, cwd=tmp_path) == 0
 
-    print(read_table(output, '/dl2/event/telescope/energy/ExtraTreesRegressor/tel_004'))
-
+    print(read_table(output, "/dl2/event/telescope/energy/ExtraTreesRegressor/tel_004"))
+    print(read_table(output, "/dl2/event/subarray/energy/ExtraTreesRegressor"))