Repository: cta-observatory/ctapipe
Pull request #2014 (2022-07-14 14:27:36)

ExpressionEngine for QualityQuery and FeatureGenerator
Breakout the logic of the expression handling/compiling to be used from the QualityQuery (current) and a FeatureGenerator (new, for the ML functionality).


Patch:
diff --git a/ctapipe/core/expression_engine.py b/ctapipe/core/expression_engine.py
new file mode 100644
index 00000000..e459e57f
--- /dev/null
+++ b/ctapipe/core/expression_engine.py
@@ -0,0 +1,51 @@
+"""
+Expression Engine
+"""
+import astropy.units as u  # for use in selection functions
+import numpy as np  # for use in selection functions
+
+from .component import Component
+from .traits import List, Tuple, Unicode
+
+# the following are what are allowed to be used
+# in selection functions (passed to eval())
+ALLOWED_GLOBALS = {"u": u, "np": np}  # astropy units  # numpy
+
+for func in ("sin", "cos", "tan", "arctan2", "log", "log10", "exp", "sqrt"):
+    ALLOWED_GLOBALS[func] = getattr(np, func)
+
+
+class ExpressionError(TypeError):
+    """Signal a problem with a user-defined selection criteria function"""
+
+
+class ExpressionEngine(Component):
+    """
+    Compile expressions on init, evaluate on call.
+    """
+
+    expressions = List(
+        Tuple(Unicode(), Unicode()),
+        help="List of 2-Tuples of Strings: ('name', 'expression').",
+    ).tag(config=True)
+
+    def __init__(self, config=None, parent=None, **kwargs):
+        super().__init__(config=config, parent=parent, **kwargs)
+
+        self.compiled = []
+        for name, expression in self.expressions:
+            try:
+                self.compiled.append(compile(expression, __name__, mode="eval"))
+            except Exception:
+                raise ExpressionError(
+                    f"Error compiling expression '{expression}' for {name}"
+                )
+
+    def __call__(self, locals):
+        for expr in self.compiled:
+            try:
+                yield eval(expr, ALLOWED_GLOBALS, locals)
+            except NameError as err:
+                raise ExpressionError(f"Error evaluating expression '{expr}': {err}")
+            except Exception as err:
+                raise ExpressionError(f"Error evaluating expression '{expr}'") from err
diff --git a/ctapipe/core/feature_generator.py b/ctapipe/core/feature_generator.py
new file mode 100644
index 00000000..f1ceab18
--- /dev/null
+++ b/ctapipe/core/feature_generator.py
@@ -0,0 +1,48 @@
+"""
+Generate Features.
+"""
+from .component import Component
+from .expression_engine import ExpressionEngine
+from .traits import List, Tuple, Unicode
+
+
+class FeatureGeneratorException(TypeError):
+    """Signal a problem with a user-defined selection criteria function"""
+
+
+class FeatureGenerator(Component):
+    """
+    Generate features for astropy.table.Table.
+
+    Raises Exceptions in two cases:
+    1. If a feature already exists in the table
+    2. If a feature cannot be built with the given expression
+    """
+
+    features = List(
+        Tuple(Unicode(), Unicode()),
+        help=(
+            "List of 2-Tuples of Strings: ('new_feature_name', 'expression to generate feature'). "
+            "You can use `numpy` as `np` and `astropy.units` as `u`. "
+            "Several math functions are usable without the `np`-prefix. "
+            "Use `feature.quantity.to_value(unit)` to create features without units."
+        ),
+    ).tag(config=True)
+
+    def __init__(self, config=None, parent=None, **kwargs):
+        super().__init__(config=config, parent=parent, **kwargs)
+        self.engine = ExpressionEngine(parent=self, expressions=self.features)
+        self._feature_names = [name for name, _ in self.features]
+
+    def __call__(self, table):
+        for result, name in zip(self.engine(table), self._feature_names):
+            if name in table.colnames:
+                raise FeatureGeneratorException(f"{name} is already a column of table.")
+            try:
+                table[name] = result
+            except Exception as err:
+                raise err
+        return table
+
+    def __len__(self):
+        return len(self.features)
diff --git a/ctapipe/core/qualityquery.py b/ctapipe/core/qualityquery.py
index fe3c3458..84b82c6b 100644
--- a/ctapipe/core/qualityquery.py
+++ b/ctapipe/core/qualityquery.py
@@ -4,18 +4,11 @@ Data Quality selection
 
 __all__ = ["QualityQuery", "QualityCriteriaError"]
 
-import astropy.units as u  # for use in selection functions
 import numpy as np  # for use in selection functions
 
 from .component import Component
-from .traits import List
-
-# the following are what are allowed to be used
-# in selection functions (passed to eval())
-ALLOWED_GLOBALS = {"u": u, "np": np}  # astropy units  # numpy
-
-for func in ("sin", "cos", "tan", "arctan2", "log", "log10", "exp", "sqrt"):
-    ALLOWED_GLOBALS[func] = getattr(np, func)
+from .expression_engine import ExpressionEngine
+from .traits import List, Tuple, Unicode
 
 
 class QualityCriteriaError(TypeError):
@@ -32,12 +25,13 @@ class QualityQuery(Component):
     """
 
     quality_criteria = List(
+        Tuple(Unicode(), Unicode()),
         help=(
             "list of tuples of ('<description', 'expression string') to accept "
             "(select) a given data value.  E.g. `[('mycut', 'x > 3'),]. "
             "You may use `numpy` as `np` and `astropy.units` as `u`, but no other"
             " modules."
-        )
+        ),
     ).tag(config=True)
 
     def __init__(self, config=None, parent=None, **kwargs):
@@ -46,21 +40,19 @@ class QualityQuery(Component):
         # add a selection to count all entries and make it the first one
         self.criteria_names = [n for (n, _) in self.quality_criteria]
         self.expressions = [e for (_, e) in self.quality_criteria]
-        self._compiled_expressions = []
-        for name, e in self.quality_criteria:
-            if "lambda" in e:
+
+        self.engine = ExpressionEngine(
+            parent=self,
+            expressions=self.quality_criteria,
+        )
+        for _, expr in self.quality_criteria:
+            if "lambda" in expr:
                 raise ValueError(
                     "As of ctapipe 0.16, do not give lambda expressions"
                     " to QualityQuery. Directly give the expression."
                     " E.g. instead of `lambda p: p.hillas.width.value > 0`"
                     " use `parameters.hillas.width.value > 0`"
                 )
-            try:
-                self._compiled_expressions.append(compile(e, __name__, mode="eval"))
-            except Exception:
-                raise QualityCriteriaError(
-                    f"Error compiling expression '{e}' for check {name}"
-                )
 
         # arrays for recording overall statistics, add one for total count
         n = len(self.quality_criteria) + 1
@@ -121,7 +113,9 @@ class QualityQuery(Component):
         # add 1 for total
         result = np.ones(len(self.quality_criteria) + 1, dtype=bool)
 
-        _evaluate_expression(self._compiled_expressions, result, kwargs)
+        for i, res in enumerate(self.engine(kwargs), start=1):
+            result[i] = res
+
         self._counts += result.astype(int)
         self._cumulative_counts += result.cumprod()
         return result[1:]  # strip off TOTAL criterion, since redundant
@@ -144,16 +138,9 @@ class QualityQuery(Component):
         n_criteria = len(self.quality_criteria) + 1
         result = np.ones((n_criteria, len(table)), dtype=bool)
 
-        _evaluate_expression(self._compiled_expressions, result, table)
+        for i, res in enumerate(self.engine(table), start=1):
+            result[i] = res
 
         self._counts += np.count_nonzero(result, axis=1)
         self._cumulative_counts += np.count_nonzero(np.cumprod(result, axis=0), axis=1)
         return np.all(result, axis=0)
-
-
-def _evaluate_expression(expressions, result, locals):
-    try:
-        for i, expression in enumerate(expressions, start=1):
-            result[i] = eval(expression, ALLOWED_GLOBALS, locals)
-    except Exception:
-        raise QualityCriteriaError("Error evaluating quality query")
diff --git a/ctapipe/core/traits.py b/ctapipe/core/traits.py
index c336e44c..793951bc 100644
--- a/ctapipe/core/traits.py
+++ b/ctapipe/core/traits.py
@@ -41,6 +41,7 @@ __all__ = [
     "Long",
     "Set",
     "TraitError",
+    "Tuple",
     "Unicode",
     "flag",
     "observe",
@@ -68,6 +69,7 @@ CaselessStrEnum = traitlets.CaselessStrEnum
 UseEnum = traitlets.UseEnum
 TraitError = traitlets.TraitError
 TraitType = traitlets.TraitType
+Tuple = traitlets.Tuple
 observe = traitlets.observe
 flag = traitlets.config.boolean_flag
 
diff --git a/ctapipe/io/tableloader.py b/ctapipe/io/tableloader.py
index f52e1758..5d724cf1 100644
--- a/ctapipe/io/tableloader.py
+++ b/ctapipe/io/tableloader.py
@@ -28,8 +28,6 @@ SHOWER_TABLE = "/simulation/event/subarray/shower"
 TRUE_IMAGES_GROUP = "/simulation/event/telescope/images"
 TRUE_PARAMETERS_GROUP = "/simulation/event/telescope/parameters"
 TRUE_IMPACT_GROUP = "/simulation/event/telescope/impact"
-SIMULATION_CONFIG_TABLE = "/configuration/simulation/run"
-SHOWER_DISTRIBUTION_TABLE = "/simulation/service/shower_distribution"
 
 DL2_SUBARRAY_GROUP = "/dl2/event/subarray"
 DL2_TELESCOPE_GROUP = "/dl2/event/telescope"
@@ -301,18 +299,6 @@ class TableLoader(Component):
         if "__index__" not in table.colnames:
             table["__index__"] = np.arange(len(table))
 
-    def read_simulation_configuration(self):
-        """
-        Read the simulation configuration table
-        """
-        return read_table(self.h5file, SIMULATION_CONFIG_TABLE)
-
-    def read_shower_distribution(self):
-        """
-        Read the simulated shower distribution histograms
-        """
-        return read_table(self.h5file, SHOWER_DISTRIBUTION_TABLE)
-
     def read_subarray_events(self, start=None, stop=None, keep_order=True):
         """Read subarray-based event information.
 
diff --git a/docs/conf.py b/docs/conf.py
index 904334eb..2d507da0 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -83,6 +83,7 @@ def setup(app):
         "Long",
         "Set",
         "TraitError",
+        "Tuple",
         "Unicode",
         "UseEnum",
     ]


Test patch:
diff --git a/ctapipe/core/tests/test_feature_generator.py b/ctapipe/core/tests/test_feature_generator.py
new file mode 100644
index 00000000..3ccb9ddd
--- /dev/null
+++ b/ctapipe/core/tests/test_feature_generator.py
@@ -0,0 +1,69 @@
+""" Tests of Selectors """
+import numpy as np
+import pytest
+from astropy.table import Table
+
+from ctapipe.core.expression_engine import ExpressionError
+from ctapipe.core.feature_generator import FeatureGenerator, FeatureGeneratorException
+
+
+def test_generate_features():
+    """Test if generating features works."""
+    expressions = [
+        ("log_intensity", "log10(intensity)"),
+        ("area", "length * width"),
+        ("eccentricity", "sqrt(1 - width ** 2 / length ** 2)"),
+    ]
+
+    generator = FeatureGenerator(features=expressions)
+
+    table = Table({"intensity": [1, 10, 100], "length": [2, 4, 8], "width": [1, 2, 4]})
+    log_intensity = [0, 1, 2]
+    area = [2, 8, 32]
+    eccentricity = np.sqrt(0.75)
+
+    table = generator(table)
+
+    assert "log_intensity" in table.colnames
+    assert "area" in table.colnames
+    assert "eccentricity" in table.colnames
+
+    assert np.all(table["log_intensity"] == log_intensity)
+    assert np.all(table["area"] == area)
+    assert np.all(table["eccentricity"] == eccentricity)
+
+
+def test_existing_feature():
+    """If the feature already exists, fail"""
+    expressions = [("foo", "bar")]
+    generator = FeatureGenerator(features=expressions)
+    table = Table({"foo": [1], "bar": [1]})
+
+    with pytest.raises(FeatureGeneratorException):
+        generator(table)
+
+
+def test_missing_colname():
+    """If the column to create a feature misses, fail"""
+    expressions = [("foo", "bar")]
+    generator = FeatureGenerator(features=expressions)
+    table = Table({"baz": [1]})
+
+    with pytest.raises(ExpressionError):
+        generator(table)
+
+
+def test_to_unit():
+    """Test chaning the unit of a feature"""
+    from astropy import units as u
+
+    expressions = [
+        ("length_meter", "length.to(u.m)"),
+        ("log_length_meter", "log10(length.quantity.to_value(u.m))"),
+    ]
+    generator = FeatureGenerator(features=expressions)
+    table = Table({"length": [1 * u.km]})
+
+    table = generator(table)
+    assert table["length_meter"] == 1000
+    assert table["log_length_meter"] == 3
diff --git a/ctapipe/core/tests/test_qualityquery.py b/ctapipe/core/tests/test_qualityquery.py
index 85c631c0..602e2e26 100644
--- a/ctapipe/core/tests/test_qualityquery.py
+++ b/ctapipe/core/tests/test_qualityquery.py
@@ -1,9 +1,10 @@
 """ Tests of Selectors """
-import pytest
 import numpy as np
+import pytest
 from astropy.table import Table
 
-from ctapipe.core.qualityquery import QualityQuery, QualityCriteriaError
+from ctapipe.core.expression_engine import ExpressionError
+from ctapipe.core.qualityquery import QualityQuery
 from ctapipe.core.traits import List
 
 
@@ -34,7 +35,7 @@ def test_selector():
     assert (criteria3 == [True, False, False]).all()
 
     criteria4 = query(x=8)  # pass all
-    assert (criteria4 == True).all()
+    assert np.all(criteria4)
 
     tab = query.to_table()
     html = query._repr_html_()
@@ -73,7 +74,7 @@ def test_invalid_input():
         ).tag(config=True)
 
     query = ExampleQualityQuery()
-    with pytest.raises(QualityCriteriaError):
+    with pytest.raises(ExpressionError):
         query(y=5)
 
 
@@ -87,19 +88,19 @@ def test_bad_selector():
             ("smallish", "x < 10"),
         ]
     )
-    with pytest.raises(QualityCriteriaError):
+    with pytest.raises(ExpressionError):
         query(x=5)
 
     # ensure we can't run arbitrary code.
     # try to construct something that is not in the
     # ALLOWED_GLOBALS list, but which is imported in selector.py
     # and see if it works in a function
-    with pytest.raises(QualityCriteriaError):
+    with pytest.raises(ExpressionError):
         query = QualityQuery(quality_criteria=[("dangerous", "Component()")])
         query(x=10)
 
     # test we only support expressions, not statements
-    with pytest.raises(QualityCriteriaError):
+    with pytest.raises(ExpressionError):
         query = QualityQuery(
             quality_criteria=[("dangerous", "import numpy; np.array([])")]
         )
@@ -132,3 +133,15 @@ def test_printing():
     )
 
     assert isinstance(str(query), str)
+
+
+def test_setup():
+    """Test that tuples can only have length 2"""
+    from traitlets import TraitError
+
+    # 2-tuple works
+    QualityQuery(quality_criteria=[("1", "2")])
+
+    # 3-tuple fails
+    with pytest.raises(TraitError):
+        QualityQuery(quality_criteria=[("1", "2", "3")])
diff --git a/ctapipe/io/tests/test_table_loader.py b/ctapipe/io/tests/test_table_loader.py
index 727593ef..fb0e32ec 100644
--- a/ctapipe/io/tests/test_table_loader.py
+++ b/ctapipe/io/tests/test_table_loader.py
@@ -1,10 +1,8 @@
-import astropy.units as u
-import numpy as np
 import pytest
 import tables
+import numpy as np
 from astropy.table import Table
 
-from ctapipe.instrument.subarray import SubarrayDescription
 from ctapipe.io.astropy_helpers import read_table
 
 
@@ -65,6 +63,9 @@ def test_check_order():
     )
 
 
+from ctapipe.instrument.subarray import SubarrayDescription
+
+
 @pytest.fixture(params=["by_type", "by_id"])
 def test_file(request, dl1_file, dl1_by_type_file):
     """Test dl1 files in both structures"""
@@ -342,24 +343,3 @@ def test_chunked(dl2_shower_geometry_file):
             assert n_events_by_type == len(tel_events)
 
     assert n_read == n_events
-
-
-def test_read_simulation_config(dl2_merged_file):
-    from ctapipe.io import TableLoader
-
-    with TableLoader(dl2_merged_file) as table_loader:
-        runs = table_loader.read_simulation_configuration()
-        assert len(runs) == 2
-        assert np.all(runs["obs_id"] == [4, 1])
-        assert u.allclose(runs["energy_range_min"].quantity, [0.004, 0.003] * u.TeV)
-
-
-def test_read_shower_distributions(dl2_merged_file):
-    from ctapipe.io import TableLoader
-
-    with TableLoader(dl2_merged_file) as table_loader:
-        histograms = table_loader.read_shower_distribution()
-        assert len(histograms) == 2
-        assert np.all(histograms["obs_id"] == [4, 1])
-        assert np.all(histograms["num_entries"] == [2000, 1000])
-        assert np.all(histograms["histogram"].sum(axis=(1, 2)) == [2000, 1000])