From 5bdb71870300622448c427424a0d3c44818429ec Mon Sep 17 00:00:00 2001 From: Christoph Toennis Date: Mon, 22 Jul 2024 12:04:13 +0200 Subject: [PATCH 01/23] Overhaul to resolve issues --- src/ctapipe/io/interpolating.py | 247 ++++++++++++++++++++++ src/ctapipe/io/tests/test_interpolator.py | 214 +++++++++++++++++++ 2 files changed, 461 insertions(+) create mode 100644 src/ctapipe/io/interpolating.py create mode 100644 src/ctapipe/io/tests/test_interpolator.py diff --git a/src/ctapipe/io/interpolating.py b/src/ctapipe/io/interpolating.py new file mode 100644 index 00000000000..98157dee68c --- /dev/null +++ b/src/ctapipe/io/interpolating.py @@ -0,0 +1,247 @@ +from typing import Any + +import astropy.units as u +import numpy as np +import tables +from astropy.time import Time +from scipy.interpolate import interp1d + +from ctapipe.core import Component, traits + +from .astropy_helpers import read_table + + +class StepInterpolator: + + """ + Step function Interpolator for the gain and pedestals + + Parameters + ---------- + values : None | np.array + Numpy array of the data that is to be interpolated. + The first dimension needs to be an index over time + times : None | np.array + Time values over which data are to be interpolated + need to be sorted and have same length as first dimension of values + """ + + def __init__( + self, + times, + values, + bounds_error=True, + fill_value="extrapolate", + assume_sorted=True, + copy=False, + ): + self.values = values + self.times = times + self.bounds_error = bounds_error + self.fill_value = fill_value + + def __call__(self, point): + if point < self.times[0]: + if self.bounds_error: + raise ValueError("below the interpolation range") + + if self.fill_value == "extrapolate": + return self.values[0] + + else: + a = np.empty(self.values[0].shape) + a[:] = np.nan + return a + + else: + i = np.searchsorted(self.times, point, side="left") + return self.values[i - 1] + + +class Interpolator(Component): + """ + Interpolate pointing and calibration parameters from a monitoring table to a given timestamp. + + Parameters + ---------- + h5file : None | tables.File + A open hdf5 file with read access. + table_location: | str + location where the monitoring data is expected to be stored in that file + interpolation_method: | str + method of interpolation to be used + """ + + bounds_error = traits.Bool( + default_value=True, + help="If true, raises an exception when trying to extrapolate out of the given table", + ).tag(config=True) + + extrapolate = traits.Bool( + help="If bounds_error is False, this flag will specify whether values outside" + "the available values are filled with nan (False) or extrapolated (True).", + default_value=False, + ).tag(config=True) + + def __init__( + self, h5file=None, table_location=None, interpolation_method=None, **kwargs + ): + super().__init__(**kwargs) + + if h5file is not None and not isinstance(h5file, tables.File): + raise TypeError("h5file must be a tables.File") + self.h5file = h5file + + if table_location is not None and not isinstance(table_location, str): + raise TypeError("table_location must be a string") + self.table_location = table_location + + self.interp_options: dict[str, Any] = dict(assume_sorted=True, copy=False) + if self.bounds_error: + self.interp_options["bounds_error"] = True + elif self.extrapolate: + self.interp_options["bounds_error"] = False + self.interp_options["fill_value"] = "extrapolate" + else: + self.interp_options["bounds_error"] = False + self.interp_options["fill_value"] = np.nan + + self._alt_interpolators = {} + self._az_interpolators = {} + self._ped_interpolators = {} + self._gain_interpolators = {} + + def add_table(self, tel_id, input_table): + """ + Add a table to this interpolator + + Parameters + ---------- + tel_id : int + Telescope id + pointing_table : astropy.table.Table + Table of pointing values, expected columns + are ``time`` as ``Time`` column, ``azimuth`` and ``altitude`` + as quantity columns. + """ + + if "gain" in set(input_table.colnames): + self.parameter_type = "gain" + # here i want to detect if the table contains pointing, gain or pedestal information + # can this be a loop? Should this be a function? + missing = {"time", "gain"} - set(input_table.colnames) + if len(missing) > 0: + raise ValueError(f"Table is missing required column(s): {missing}") + + elif "pedestal" in set(input_table.colnames): + self.parameter_type = "pedestal" + missing = {"time", "pedestal"} - set(input_table.colnames) + if len(missing) > 0: + raise ValueError(f"Table is missing required column(s): {missing}") + + elif "azimuth" in set(input_table.colnames): + self.parameter_type = "pointing" + missing = {"time", "azimuth", "altitude"} - set(input_table.colnames) + if len(missing) > 0: + raise ValueError(f"Table is missing required column(s): {missing}") + + for col in ("azimuth", "altitude"): + unit = input_table[col].unit + if unit is None or not u.rad.is_equivalent(unit): + raise ValueError(f"{col} must have units compatible with 'rad'") + + if not isinstance(input_table["time"], Time): + raise TypeError( + "'time' column of pointing table must be astropy.time.Time" + ) + # sort first, so it's not done twice for each interpolator + input_table.sort("time") + if self.parameter_type == "pointing": + az = input_table["azimuth"].quantity.to_value(u.rad) + # prepare azimuth for interpolation by "unwrapping": i.e. turning + # [359, 1] into [359, 361]. This assumes that if we get values like + # [359, 1] the telescope moved 2 degrees through 0, not 358 degrees + # the other way around. This should be true for all telescopes given + # the sampling speed of pointing values and their maximum movement speed. + # No telescope can turn more than 180° in 2 seconds. + az = np.unwrap(az) + alt = input_table["altitude"].quantity.to_value(u.rad) + mjd = input_table["time"].tai.mjd + self._az_interpolators[tel_id] = interp1d(mjd, az, **self.interp_options) + self._alt_interpolators[tel_id] = interp1d(mjd, alt, **self.interp_options) + + elif self.parameter_type == "pedestal": + time = input_table["time"] + + ped = input_table["pedestal"] + + self._ped_interpolators[tel_id] = StepInterpolator( + time, ped, **self.interp_options + ) + + elif self.parameter_type == "gain": + time = input_table["time"] + + gain = input_table["gain"] + + self._gain_interpolators[tel_id] = StepInterpolator( + time, gain, **self.interp_options + ) + + def _read_parameter_table(self, tel_id, table_location=None): + if table_location is None: + table_location = f"/dl0/monitoring/telescope/pointing/tel_{tel_id:03d}" + + pointing_table = read_table( + self.h5file, + table_location, + ) + self.add_table(tel_id, pointing_table) + + def __call__(self, tel_id, time): + """ + Interpolate alt/az for given time and tel_id. + + Parameters + ---------- + tel_id : int + telescope id + time : astropy.time.Time + time for which to interpolate the pointing + + Returns + ------- + altitude : astropy.units.Quantity[deg] + interpolated altitude angle + azimuth : astropy.units.Quantity[deg] + interpolated azimuth angle + """ + + if not any( + [ + tel_id in x + for x in ( + self._az_interpolators, + self._ped_interpolators, + self._gain_interpolators, + ) + ] + ): + if self.h5file is not None: + self._read_parameter_table(tel_id) + else: + raise KeyError(f"No table available for tel_id {tel_id}") + + if self.parameter_type == "pointing": + mjd = time.tai.mjd + az = u.Quantity(self._az_interpolators[tel_id](mjd), u.rad, copy=False) + alt = u.Quantity(self._alt_interpolators[tel_id](mjd), u.rad, copy=False) + return alt, az + + elif self.parameter_type == "pedestal": + ped = self._ped_interpolators[tel_id](time) + return ped + + elif self.parameter_type == "gain": + gain = self._gain_interpolators[tel_id](time) + return gain diff --git a/src/ctapipe/io/tests/test_interpolator.py b/src/ctapipe/io/tests/test_interpolator.py new file mode 100644 index 00000000000..128b3bc2e04 --- /dev/null +++ b/src/ctapipe/io/tests/test_interpolator.py @@ -0,0 +1,214 @@ +import astropy.units as u +import numpy as np +import pytest +import tables +from astropy.table import Table +from astropy.time import Time + + +def test_simple(): + """Test interpolator""" + from ctapipe.io.interpolating import Interpolator + + t0 = Time("2022-01-01T00:00:00") + + table = Table( + { + "time": t0 + np.arange(0.0, 10.1, 2.0) * u.s, + "azimuth": np.linspace(0.0, 10.0, 6) * u.deg, + "altitude": np.linspace(70.0, 60.0, 6) * u.deg, + }, + ) + + table_pedestal = Table( + { + "time": np.arange(0.0, 10.1, 2.0), + "pedestal": np.reshape(np.random.normal(4.0, 1.0, 1850 * 6), (6, 1850)), + }, + ) + + table_gain = Table( + { + "time": np.arange(0.0, 10.1, 2.0), + "gain": np.reshape(np.random.normal(1.0, 0.2, 1850 * 6), (6, 1850)), + }, + ) + + interpolator = Interpolator() + interpolator_pedestal = Interpolator() + interpolator_gain = Interpolator() + interpolator.add_table(1, table) + interpolator_pedestal.add_table(1, table_pedestal) + interpolator_gain.add_table(1, table_gain) + + alt, az = interpolator(tel_id=1, time=t0 + 1 * u.s) + assert u.isclose(alt, 69 * u.deg) + assert u.isclose(az, 1 * u.deg) + + pedestal = interpolator_pedestal(tel_id=1, time=1.0) + assert all(pedestal == table_pedestal["pedestal"][0]) + + gain = interpolator_gain(tel_id=1, time=1.0) + assert all(gain == table_gain["gain"][0]) + + with pytest.raises(KeyError): + interpolator(tel_id=2, time=t0 + 1 * u.s) + with pytest.raises(KeyError): + interpolator_pedestal(tel_id=2, time=1.0) + with pytest.raises(KeyError): + interpolator_gain(tel_id=2, time=1.0) + + +def test_azimuth_switchover(): + """Test pointing interpolation""" + from ctapipe.io.interpolating import Interpolator + + t0 = Time("2022-01-01T00:00:00") + + table = Table( + { + "time": t0 + [0, 1, 2] * u.s, + "azimuth": [359, 1, 3] * u.deg, + "altitude": [60, 61, 62] * u.deg, + }, + ) + + interpolator = Interpolator() + interpolator.add_table(1, table) + + alt, az = interpolator(tel_id=1, time=t0 + 0.5 * u.s) + assert u.isclose(az, 360 * u.deg) + assert u.isclose(alt, 60.5 * u.deg) + + +def test_invalid_input(): + """Test invalid pointing tables raise nice errors""" + from ctapipe.io.interpolating import Interpolator + + wrong_time = Table( + { + "time": [1, 2, 3] * u.s, + "azimuth": [1, 2, 3] * u.deg, + "altitude": [1, 2, 3] * u.deg, + } + ) + + interpolator = Interpolator() + with pytest.raises(TypeError, match="astropy.time.Time"): + interpolator.add_table(1, wrong_time) + + wrong_unit = Table( + { + "time": Time(1.7e9 + np.arange(3), format="unix"), + "azimuth": [1, 2, 3] * u.m, + "altitude": [1, 2, 3] * u.deg, + } + ) + with pytest.raises(ValueError, match="compatible with 'rad'"): + interpolator.add_table(1, wrong_unit) + + wrong_unit = Table( + { + "time": Time(1.7e9 + np.arange(3), format="unix"), + "azimuth": [1, 2, 3] * u.deg, + "altitude": [1, 2, 3], + } + ) + with pytest.raises(ValueError, match="compatible with 'rad'"): + interpolator.add_table(1, wrong_unit) + + +def test_hdf5(tmp_path): + from ctapipe.io import write_table + from ctapipe.io.interpolating import Interpolator + + t0 = Time("2022-01-01T00:00:00") + + table = Table( + { + "time": t0 + np.arange(0.0, 10.1, 2.0) * u.s, + "azimuth": np.linspace(0.0, 10.0, 6) * u.deg, + "altitude": np.linspace(70.0, 60.0, 6) * u.deg, + }, + ) + + path = tmp_path / "pointing.h5" + write_table(table, path, "/dl0/monitoring/telescope/pointing/tel_001") + with tables.open_file(path) as h5file: + interpolator = Interpolator(h5file) + alt, az = interpolator(tel_id=1, time=t0 + 1 * u.s) + assert u.isclose(alt, 69 * u.deg) + assert u.isclose(az, 1 * u.deg) + + +def test_bounds(): + """Test invalid pointing tables raise nice errors""" + from ctapipe.io.interpolating import Interpolator + + t0 = Time("2022-01-01T00:00:00") + + table_pointing = Table( + { + "time": t0 + np.arange(0.0, 10.1, 2.0) * u.s, + "azimuth": np.linspace(0.0, 10.0, 6) * u.deg, + "altitude": np.linspace(70.0, 60.0, 6) * u.deg, + }, + ) + + table_pedestal = Table( + { + "time": np.arange(0.0, 10.1, 2.0), + "pedestal": np.reshape(np.random.normal(4.0, 1.0, 1850 * 6), (6, 1850)), + }, + ) + + table_gain = Table( + { + "time": np.arange(0.0, 10.1, 2.0), + "gain": np.reshape(np.random.normal(1.0, 0.2, 1850 * 6), (6, 1850)), + }, + ) + + interpolator_pointing = Interpolator() + interpolator_pedestal = Interpolator() + interpolator_gain = Interpolator() + interpolator_pointing.add_table(1, table_pointing) + interpolator_pedestal.add_table(1, table_pedestal) + interpolator_gain.add_table(1, table_gain) + + with pytest.raises(ValueError, match="below the interpolation range"): + interpolator_pointing(tel_id=1, time=t0 - 0.1 * u.s) + + with pytest.raises(ValueError, match="below the interpolation range"): + interpolator_pedestal(tel_id=1, time=-0.1) + + with pytest.raises(ValueError, match="below the interpolation range"): + interpolator_gain(tel_id=1, time=-0.1) + + with pytest.raises(ValueError, match="above the interpolation range"): + interpolator_pointing(tel_id=1, time=t0 + 10.2 * u.s) + + interpolator_pointing = Interpolator(bounds_error=False) + interpolator_pedestal = Interpolator(bounds_error=False) + interpolator_gain = Interpolator(bounds_error=False) + interpolator_pointing.add_table(1, table_pointing) + interpolator_pedestal.add_table(1, table_pedestal) + interpolator_gain.add_table(1, table_gain) + + for dt in (-0.1, 10.1) * u.s: + alt, az = interpolator_pointing(tel_id=1, time=t0 + dt) + assert np.isnan(alt.value) + assert np.isnan(az.value) + + assert all(np.isnan(interpolator_pedestal(tel_id=1, time=-0.1))) + assert all(np.isnan(interpolator_gain(tel_id=1, time=-0.1))) + + interpolator = Interpolator(bounds_error=False, extrapolate=True) + interpolator.add_table(1, table_pointing) + alt, az = interpolator(tel_id=1, time=t0 - 1 * u.s) + assert u.isclose(alt, 71 * u.deg) + assert u.isclose(az, -1 * u.deg) + + alt, az = interpolator(tel_id=1, time=t0 + 11 * u.s) + assert u.isclose(alt, 59 * u.deg) + assert u.isclose(az, 11 * u.deg) From 28d6869032164cbd3b8d36c946ca3223c10416ff Mon Sep 17 00:00:00 2001 From: Christoph Toennis Date: Mon, 22 Jul 2024 14:12:20 +0200 Subject: [PATCH 02/23] I switched the pointing interpolator to the generalized on in the event source and table loader --- src/ctapipe/io/hdf5eventsource.py | 4 ++-- src/ctapipe/io/tableloader.py | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/ctapipe/io/hdf5eventsource.py b/src/ctapipe/io/hdf5eventsource.py index afbe2153d2c..e498da56af7 100644 --- a/src/ctapipe/io/hdf5eventsource.py +++ b/src/ctapipe/io/hdf5eventsource.py @@ -52,7 +52,7 @@ from .datalevels import DataLevel from .eventsource import EventSource from .hdf5tableio import HDF5TableReader -from .pointing import PointingInterpolator +from .interpolating import Interpolator from .tableloader import DL2_SUBARRAY_GROUP, DL2_TELESCOPE_GROUP, POINTING_GROUP __all__ = ["HDF5EventSource"] @@ -585,7 +585,7 @@ def _generator(self): pointing_interpolator = None if POINTING_GROUP in self.file_.root: - pointing_interpolator = PointingInterpolator( + pointing_interpolator = Interpolator( h5file=self.file_, parent=self, ) diff --git a/src/ctapipe/io/tableloader.py b/src/ctapipe/io/tableloader.py index 36a954fb7de..033bfc8e53e 100644 --- a/src/ctapipe/io/tableloader.py +++ b/src/ctapipe/io/tableloader.py @@ -14,7 +14,7 @@ from ..core import Component, Provenance, traits from ..instrument import FocalLengthKind, SubarrayDescription from .astropy_helpers import join_allow_empty, read_table -from .pointing import PointingInterpolator +from .interpolating import Interpolator __all__ = ["TableLoader"] @@ -248,7 +248,7 @@ def __init__(self, input_url=None, h5file=None, **kwargs): Provenance().add_input_file(self.input_url, role="Event data") - self._pointing_interpolator = PointingInterpolator( + self._pointing_interpolator = Interpolator( h5file=self.h5file, parent=self, ) From 2886a7d2a4aaa1b6d05df929da809af098527efb Mon Sep 17 00:00:00 2001 From: Christoph Toennis Date: Mon, 22 Jul 2024 15:42:21 +0200 Subject: [PATCH 03/23] I simplifid the interpolator instances inside the class --- src/ctapipe/io/interpolating.py | 47 +++++++++++++-------------------- 1 file changed, 18 insertions(+), 29 deletions(-) diff --git a/src/ctapipe/io/interpolating.py b/src/ctapipe/io/interpolating.py index 98157dee68c..8a0d1043914 100644 --- a/src/ctapipe/io/interpolating.py +++ b/src/ctapipe/io/interpolating.py @@ -106,10 +106,8 @@ def __init__( self.interp_options["bounds_error"] = False self.interp_options["fill_value"] = np.nan - self._alt_interpolators = {} - self._az_interpolators = {} - self._ped_interpolators = {} - self._gain_interpolators = {} + self._interpolators = {} + self._secondary_interpolators = {} def add_table(self, tel_id, input_table): """ @@ -167,15 +165,17 @@ def add_table(self, tel_id, input_table): az = np.unwrap(az) alt = input_table["altitude"].quantity.to_value(u.rad) mjd = input_table["time"].tai.mjd - self._az_interpolators[tel_id] = interp1d(mjd, az, **self.interp_options) - self._alt_interpolators[tel_id] = interp1d(mjd, alt, **self.interp_options) + self._interpolators[tel_id] = interp1d(mjd, az, **self.interp_options) + self._secondary_interpolators[tel_id] = interp1d( + mjd, alt, **self.interp_options + ) elif self.parameter_type == "pedestal": time = input_table["time"] ped = input_table["pedestal"] - self._ped_interpolators[tel_id] = StepInterpolator( + self._interpolators[tel_id] = StepInterpolator( time, ped, **self.interp_options ) @@ -184,7 +184,7 @@ def add_table(self, tel_id, input_table): gain = input_table["gain"] - self._gain_interpolators[tel_id] = StepInterpolator( + self._interpolators[tel_id] = StepInterpolator( time, gain, **self.interp_options ) @@ -192,11 +192,11 @@ def _read_parameter_table(self, tel_id, table_location=None): if table_location is None: table_location = f"/dl0/monitoring/telescope/pointing/tel_{tel_id:03d}" - pointing_table = read_table( + input_table = read_table( self.h5file, table_location, ) - self.add_table(tel_id, pointing_table) + self.add_table(tel_id, input_table) def __call__(self, tel_id, time): """ @@ -217,16 +217,7 @@ def __call__(self, tel_id, time): interpolated azimuth angle """ - if not any( - [ - tel_id in x - for x in ( - self._az_interpolators, - self._ped_interpolators, - self._gain_interpolators, - ) - ] - ): + if tel_id not in self._interpolators: if self.h5file is not None: self._read_parameter_table(tel_id) else: @@ -234,14 +225,12 @@ def __call__(self, tel_id, time): if self.parameter_type == "pointing": mjd = time.tai.mjd - az = u.Quantity(self._az_interpolators[tel_id](mjd), u.rad, copy=False) - alt = u.Quantity(self._alt_interpolators[tel_id](mjd), u.rad, copy=False) + az = u.Quantity(self._interpolators[tel_id](mjd), u.rad, copy=False) + alt = u.Quantity( + self._secondary_interpolators[tel_id](mjd), u.rad, copy=False + ) return alt, az - elif self.parameter_type == "pedestal": - ped = self._ped_interpolators[tel_id](time) - return ped - - elif self.parameter_type == "gain": - gain = self._gain_interpolators[tel_id](time) - return gain + elif self.parameter_type in ("pedestal", "gain"): + cal = self._interpolators[tel_id](time) + return cal From c4dd62d115896a295213689090490b3318ea841b Mon Sep 17 00:00:00 2001 From: Christoph Toennis Date: Tue, 23 Jul 2024 10:43:08 +0200 Subject: [PATCH 04/23] Updated docstrings and variable names --- src/ctapipe/io/interpolating.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/ctapipe/io/interpolating.py b/src/ctapipe/io/interpolating.py index 8a0d1043914..a9e2d980fa2 100644 --- a/src/ctapipe/io/interpolating.py +++ b/src/ctapipe/io/interpolating.py @@ -117,10 +117,12 @@ def add_table(self, tel_id, input_table): ---------- tel_id : int Telescope id - pointing_table : astropy.table.Table + input_table : astropy.table.Table Table of pointing values, expected columns are ``time`` as ``Time`` column, ``azimuth`` and ``altitude`` - as quantity columns. + as quantity columns for pointing and pointing correction data. + For pedestal data the quantity column is expected as + ``pedestal`` and for gain data as ``gain``. """ if "gain" in set(input_table.colnames): @@ -188,8 +190,8 @@ def add_table(self, tel_id, input_table): time, gain, **self.interp_options ) - def _read_parameter_table(self, tel_id, table_location=None): - if table_location is None: + def _read_parameter_table(self, tel_id, table_location="pointing"): + if table_location == "pointing": table_location = f"/dl0/monitoring/telescope/pointing/tel_{tel_id:03d}" input_table = read_table( From 6d185a1cea4b09c031213b58ac1ab20d33d7f2b5 Mon Sep 17 00:00:00 2001 From: Christoph Toennis Date: Tue, 23 Jul 2024 11:35:30 +0200 Subject: [PATCH 05/23] Added notes --- docs/changes/2600.feature.rst | 2 ++ 1 file changed, 2 insertions(+) create mode 100644 docs/changes/2600.feature.rst diff --git a/docs/changes/2600.feature.rst b/docs/changes/2600.feature.rst new file mode 100644 index 00000000000..aead92bc84b --- /dev/null +++ b/docs/changes/2600.feature.rst @@ -0,0 +1,2 @@ +Add Interpolator class to generalize the PointingInterpolator in the io collection. +Add new interpolator class called StepInterpolator to be used for calibration data. From 9579dbc3791fb64c22ff7bc06e2f7d6e11d2f79a Mon Sep 17 00:00:00 2001 From: Christoph Toennis Date: Tue, 23 Jul 2024 13:50:36 +0200 Subject: [PATCH 06/23] I reorganized a bit of code to simplify things --- src/ctapipe/io/interpolating.py | 70 +++++++++++------------ src/ctapipe/io/tests/test_interpolator.py | 8 ++- 2 files changed, 37 insertions(+), 41 deletions(-) diff --git a/src/ctapipe/io/interpolating.py b/src/ctapipe/io/interpolating.py index a9e2d980fa2..0ce7d54b817 100644 --- a/src/ctapipe/io/interpolating.py +++ b/src/ctapipe/io/interpolating.py @@ -125,35 +125,8 @@ def add_table(self, tel_id, input_table): ``pedestal`` and for gain data as ``gain``. """ - if "gain" in set(input_table.colnames): - self.parameter_type = "gain" - # here i want to detect if the table contains pointing, gain or pedestal information - # can this be a loop? Should this be a function? - missing = {"time", "gain"} - set(input_table.colnames) - if len(missing) > 0: - raise ValueError(f"Table is missing required column(s): {missing}") - - elif "pedestal" in set(input_table.colnames): - self.parameter_type = "pedestal" - missing = {"time", "pedestal"} - set(input_table.colnames) - if len(missing) > 0: - raise ValueError(f"Table is missing required column(s): {missing}") - - elif "azimuth" in set(input_table.colnames): - self.parameter_type = "pointing" - missing = {"time", "azimuth", "altitude"} - set(input_table.colnames) - if len(missing) > 0: - raise ValueError(f"Table is missing required column(s): {missing}") + self._identify_data(input_table) - for col in ("azimuth", "altitude"): - unit = input_table[col].unit - if unit is None or not u.rad.is_equivalent(unit): - raise ValueError(f"{col} must have units compatible with 'rad'") - - if not isinstance(input_table["time"], Time): - raise TypeError( - "'time' column of pointing table must be astropy.time.Time" - ) # sort first, so it's not done twice for each interpolator input_table.sort("time") if self.parameter_type == "pointing": @@ -172,23 +145,44 @@ def add_table(self, tel_id, input_table): mjd, alt, **self.interp_options ) - elif self.parameter_type == "pedestal": + elif self.parameter_type in ("pedestal", "gain"): time = input_table["time"] - - ped = input_table["pedestal"] + if self.parameter_type == "gain": + cal = input_table["gain"] + else: + cal = input_table["pedestal"] self._interpolators[tel_id] = StepInterpolator( - time, ped, **self.interp_options + time, cal, **self.interp_options ) - elif self.parameter_type == "gain": - time = input_table["time"] + def _identify_data(self, input_table): + missing = {} + + if "gain" in set(input_table.colnames): + self.parameter_type = "gain" + # here i want to detect if the table contains pointing, gain or pedestal information + # can this be a loop? Should this be a function? + missing = {"time", "gain"} - set(input_table.colnames) - gain = input_table["gain"] + elif "pedestal" in set(input_table.colnames): + self.parameter_type = "pedestal" + missing = {"time", "pedestal"} - set(input_table.colnames) - self._interpolators[tel_id] = StepInterpolator( - time, gain, **self.interp_options - ) + elif "azimuth" in set(input_table.colnames): + self.parameter_type = "pointing" + missing = {"time", "azimuth", "altitude"} - set(input_table.colnames) + for col in ("azimuth", "altitude"): + unit = input_table[col].unit + if unit is None or not u.rad.is_equivalent(unit): + raise ValueError(f"{col} must have units compatible with 'rad'") + + if not isinstance(input_table["time"], Time): + raise TypeError( + "'time' column of pointing table must be astropy.time.Time" + ) + if len(missing) > 0: + raise ValueError(f"Table is missing required column(s): {missing}") def _read_parameter_table(self, tel_id, table_location="pointing"): if table_location == "pointing": diff --git a/src/ctapipe/io/tests/test_interpolator.py b/src/ctapipe/io/tests/test_interpolator.py index 128b3bc2e04..46dcffc8017 100644 --- a/src/ctapipe/io/tests/test_interpolator.py +++ b/src/ctapipe/io/tests/test_interpolator.py @@ -176,13 +176,15 @@ def test_bounds(): interpolator_pedestal.add_table(1, table_pedestal) interpolator_gain.add_table(1, table_gain) - with pytest.raises(ValueError, match="below the interpolation range"): + error_message = "below the interpolation range" + + with pytest.raises(ValueError, match=error_message): interpolator_pointing(tel_id=1, time=t0 - 0.1 * u.s) - with pytest.raises(ValueError, match="below the interpolation range"): + with pytest.raises(ValueError, match=error_message): interpolator_pedestal(tel_id=1, time=-0.1) - with pytest.raises(ValueError, match="below the interpolation range"): + with pytest.raises(ValueError, match=error_message): interpolator_gain(tel_id=1, time=-0.1) with pytest.raises(ValueError, match="above the interpolation range"): From bcd9ea9f5de1b2e07c25090da3a4a4935ed4af6e Mon Sep 17 00:00:00 2001 From: Christoph Toennis Date: Wed, 24 Jul 2024 10:21:23 +0200 Subject: [PATCH 07/23] I split the interpolator into a parent and 2 child classes --- src/ctapipe/io/hdf5eventsource.py | 4 +- src/ctapipe/io/interpolating.py | 194 ++++++++++++---------- src/ctapipe/io/tableloader.py | 4 +- src/ctapipe/io/tests/test_interpolator.py | 80 +++------ 4 files changed, 131 insertions(+), 151 deletions(-) diff --git a/src/ctapipe/io/hdf5eventsource.py b/src/ctapipe/io/hdf5eventsource.py index e498da56af7..e7bcb1a5073 100644 --- a/src/ctapipe/io/hdf5eventsource.py +++ b/src/ctapipe/io/hdf5eventsource.py @@ -52,7 +52,7 @@ from .datalevels import DataLevel from .eventsource import EventSource from .hdf5tableio import HDF5TableReader -from .interpolating import Interpolator +from .interpolating import PointingInterpolator from .tableloader import DL2_SUBARRAY_GROUP, DL2_TELESCOPE_GROUP, POINTING_GROUP __all__ = ["HDF5EventSource"] @@ -585,7 +585,7 @@ def _generator(self): pointing_interpolator = None if POINTING_GROUP in self.file_.root: - pointing_interpolator = Interpolator( + pointing_interpolator = PointingInterpolator( h5file=self.file_, parent=self, ) diff --git a/src/ctapipe/io/interpolating.py b/src/ctapipe/io/interpolating.py index 0ce7d54b817..4b764b5bfe9 100644 --- a/src/ctapipe/io/interpolating.py +++ b/src/ctapipe/io/interpolating.py @@ -11,7 +11,7 @@ from .astropy_helpers import read_table -class StepInterpolator: +class StepFunction: """ Step function Interpolator for the gain and pedestals @@ -60,16 +60,12 @@ def __call__(self, point): class Interpolator(Component): """ - Interpolate pointing and calibration parameters from a monitoring table to a given timestamp. + Interpolator parent class. Parameters ---------- h5file : None | tables.File A open hdf5 file with read access. - table_location: | str - location where the monitoring data is expected to be stored in that file - interpolation_method: | str - method of interpolation to be used """ bounds_error = traits.Bool( @@ -83,19 +79,13 @@ class Interpolator(Component): default_value=False, ).tag(config=True) - def __init__( - self, h5file=None, table_location=None, interpolation_method=None, **kwargs - ): + def __init__(self, h5file=None, **kwargs): super().__init__(**kwargs) if h5file is not None and not isinstance(h5file, tables.File): raise TypeError("h5file must be a tables.File") self.h5file = h5file - if table_location is not None and not isinstance(table_location, str): - raise TypeError("table_location must be a string") - self.table_location = table_location - self.interp_options: dict[str, Any] = dict(assume_sorted=True, copy=False) if self.bounds_error: self.interp_options["bounds_error"] = True @@ -109,6 +99,48 @@ def __init__( self._interpolators = {} self._secondary_interpolators = {} + def add_table(self, tel_id, input_table): + # sort first, so it's not done twice for each interpolator + input_table.sort("time") + + def __call__(self, tel_id, time): + pass + + def _read_parameter_table(self, tel_id): + pass + + +class PointingInterpolator(Interpolator): + def __call__(self, tel_id, time): + """ + Interpolate alt/az for given time and tel_id. + + Parameters + ---------- + tel_id : int + telescope id + time : astropy.time.Time + time for which to interpolate the pointing + + Returns + ------- + altitude : astropy.units.Quantity[deg] + interpolated altitude angle + azimuth : astropy.units.Quantity[deg] + interpolated azimuth angle + """ + + if tel_id not in self._interpolators: + if self.h5file is not None: + self._read_parameter_table(tel_id) + else: + raise KeyError(f"No table available for tel_id {tel_id}") + + mjd = time.tai.mjd + az = u.Quantity(self._interpolators[tel_id](mjd), u.rad, copy=False) + alt = u.Quantity(self._secondary_interpolators[tel_id](mjd), u.rad, copy=False) + return alt, az + def add_table(self, tel_id, input_table): """ Add a table to this interpolator @@ -121,82 +153,46 @@ def add_table(self, tel_id, input_table): Table of pointing values, expected columns are ``time`` as ``Time`` column, ``azimuth`` and ``altitude`` as quantity columns for pointing and pointing correction data. - For pedestal data the quantity column is expected as - ``pedestal`` and for gain data as ``gain``. """ - self._identify_data(input_table) - - # sort first, so it's not done twice for each interpolator - input_table.sort("time") - if self.parameter_type == "pointing": - az = input_table["azimuth"].quantity.to_value(u.rad) - # prepare azimuth for interpolation by "unwrapping": i.e. turning - # [359, 1] into [359, 361]. This assumes that if we get values like - # [359, 1] the telescope moved 2 degrees through 0, not 358 degrees - # the other way around. This should be true for all telescopes given - # the sampling speed of pointing values and their maximum movement speed. - # No telescope can turn more than 180° in 2 seconds. - az = np.unwrap(az) - alt = input_table["altitude"].quantity.to_value(u.rad) - mjd = input_table["time"].tai.mjd - self._interpolators[tel_id] = interp1d(mjd, az, **self.interp_options) - self._secondary_interpolators[tel_id] = interp1d( - mjd, alt, **self.interp_options - ) - - elif self.parameter_type in ("pedestal", "gain"): - time = input_table["time"] - if self.parameter_type == "gain": - cal = input_table["gain"] - else: - cal = input_table["pedestal"] - - self._interpolators[tel_id] = StepInterpolator( - time, cal, **self.interp_options - ) - - def _identify_data(self, input_table): - missing = {} - - if "gain" in set(input_table.colnames): - self.parameter_type = "gain" - # here i want to detect if the table contains pointing, gain or pedestal information - # can this be a loop? Should this be a function? - missing = {"time", "gain"} - set(input_table.colnames) - - elif "pedestal" in set(input_table.colnames): - self.parameter_type = "pedestal" - missing = {"time", "pedestal"} - set(input_table.colnames) - - elif "azimuth" in set(input_table.colnames): - self.parameter_type = "pointing" - missing = {"time", "azimuth", "altitude"} - set(input_table.colnames) - for col in ("azimuth", "altitude"): - unit = input_table[col].unit - if unit is None or not u.rad.is_equivalent(unit): - raise ValueError(f"{col} must have units compatible with 'rad'") - - if not isinstance(input_table["time"], Time): - raise TypeError( - "'time' column of pointing table must be astropy.time.Time" - ) + missing = {"time", "azimuth", "altitude"} - set(input_table.colnames) if len(missing) > 0: raise ValueError(f"Table is missing required column(s): {missing}") + for col in ("azimuth", "altitude"): + unit = input_table[col].unit + if unit is None or not u.rad.is_equivalent(unit): + raise ValueError(f"{col} must have units compatible with 'rad'") + + if not isinstance(input_table["time"], Time): + raise TypeError("'time' column of pointing table must be astropy.time.Time") + + az = input_table["azimuth"].quantity.to_value(u.rad) + # prepare azimuth for interpolation by "unwrapping": i.e. turning + # [359, 1] into [359, 361]. This assumes that if we get values like + # [359, 1] the telescope moved 2 degrees through 0, not 358 degrees + # the other way around. This should be true for all telescopes given + # the sampling speed of pointing values and their maximum movement speed. + # No telescope can turn more than 180° in 2 seconds. + az = np.unwrap(az) + alt = input_table["altitude"].quantity.to_value(u.rad) + mjd = input_table["time"].tai.mjd + self._interpolators[tel_id] = interp1d(mjd, az, **self.interp_options) + self._secondary_interpolators[tel_id] = interp1d( + mjd, alt, **self.interp_options + ) - def _read_parameter_table(self, tel_id, table_location="pointing"): - if table_location == "pointing": - table_location = f"/dl0/monitoring/telescope/pointing/tel_{tel_id:03d}" - + def _read_parameter_table(self, tel_id): input_table = read_table( self.h5file, - table_location, + f"/dl0/monitoring/telescope/pointing/tel_{tel_id:03d}", ) self.add_table(tel_id, input_table) + +class CalibrationInterpolator(Interpolator): def __call__(self, tel_id, time): """ - Interpolate alt/az for given time and tel_id. + Interpolate pedestal or gain for a given time and tel_id. Parameters ---------- @@ -207,10 +203,8 @@ def __call__(self, tel_id, time): Returns ------- - altitude : astropy.units.Quantity[deg] - interpolated altitude angle - azimuth : astropy.units.Quantity[deg] - interpolated azimuth angle + cal : array [float] + interpolated calibration quantity """ if tel_id not in self._interpolators: @@ -219,14 +213,30 @@ def __call__(self, tel_id, time): else: raise KeyError(f"No table available for tel_id {tel_id}") - if self.parameter_type == "pointing": - mjd = time.tai.mjd - az = u.Quantity(self._interpolators[tel_id](mjd), u.rad, copy=False) - alt = u.Quantity( - self._secondary_interpolators[tel_id](mjd), u.rad, copy=False - ) - return alt, az - - elif self.parameter_type in ("pedestal", "gain"): - cal = self._interpolators[tel_id](time) - return cal + cal = self._interpolators[tel_id](time) + return cal + + def add_table(self, tel_id, input_table, par_name): + """ + Add a table to this interpolator + + Parameters + ---------- + tel_id : int + Telescope id + input_table : astropy.table.Table + Table of pointing values, expected columns + are ``time`` as ``Time`` column. The calibration + parameter column is given through the variable ``par_name`` + par_name : str + Name of the parameter that is to be interpolated + """ + + missing = {"time", par_name} - set(input_table.colnames) + if len(missing) > 0: + raise ValueError(f"Table is missing required column(s): {missing}") + + input_table.sort("time") + time = input_table["time"] + cal = input_table[par_name] + self._interpolators[tel_id] = StepFunction(time, cal, **self.interp_options) diff --git a/src/ctapipe/io/tableloader.py b/src/ctapipe/io/tableloader.py index 033bfc8e53e..771058ffad7 100644 --- a/src/ctapipe/io/tableloader.py +++ b/src/ctapipe/io/tableloader.py @@ -14,7 +14,7 @@ from ..core import Component, Provenance, traits from ..instrument import FocalLengthKind, SubarrayDescription from .astropy_helpers import join_allow_empty, read_table -from .interpolating import Interpolator +from .interpolating import PointingInterpolator __all__ = ["TableLoader"] @@ -248,7 +248,7 @@ def __init__(self, input_url=None, h5file=None, **kwargs): Provenance().add_input_file(self.input_url, role="Event data") - self._pointing_interpolator = Interpolator( + self._pointing_interpolator = PointingInterpolator( h5file=self.h5file, parent=self, ) diff --git a/src/ctapipe/io/tests/test_interpolator.py b/src/ctapipe/io/tests/test_interpolator.py index 46dcffc8017..b014cd51d98 100644 --- a/src/ctapipe/io/tests/test_interpolator.py +++ b/src/ctapipe/io/tests/test_interpolator.py @@ -8,7 +8,7 @@ def test_simple(): """Test interpolator""" - from ctapipe.io.interpolating import Interpolator + from ctapipe.io.interpolating import CalibrationInterpolator, PointingInterpolator t0 = Time("2022-01-01T00:00:00") @@ -20,48 +20,33 @@ def test_simple(): }, ) - table_pedestal = Table( + table_cal = Table( { "time": np.arange(0.0, 10.1, 2.0), "pedestal": np.reshape(np.random.normal(4.0, 1.0, 1850 * 6), (6, 1850)), }, ) - table_gain = Table( - { - "time": np.arange(0.0, 10.1, 2.0), - "gain": np.reshape(np.random.normal(1.0, 0.2, 1850 * 6), (6, 1850)), - }, - ) - - interpolator = Interpolator() - interpolator_pedestal = Interpolator() - interpolator_gain = Interpolator() + interpolator = PointingInterpolator() + interpolator_cal = CalibrationInterpolator() interpolator.add_table(1, table) - interpolator_pedestal.add_table(1, table_pedestal) - interpolator_gain.add_table(1, table_gain) + interpolator_cal.add_table(1, table_cal, "pedestal") alt, az = interpolator(tel_id=1, time=t0 + 1 * u.s) assert u.isclose(alt, 69 * u.deg) assert u.isclose(az, 1 * u.deg) - pedestal = interpolator_pedestal(tel_id=1, time=1.0) - assert all(pedestal == table_pedestal["pedestal"][0]) - - gain = interpolator_gain(tel_id=1, time=1.0) - assert all(gain == table_gain["gain"][0]) - + pedestal = interpolator_cal(tel_id=1, time=1.0) + assert all(pedestal == table_cal["pedestal"][0]) with pytest.raises(KeyError): interpolator(tel_id=2, time=t0 + 1 * u.s) with pytest.raises(KeyError): - interpolator_pedestal(tel_id=2, time=1.0) - with pytest.raises(KeyError): - interpolator_gain(tel_id=2, time=1.0) + interpolator_cal(tel_id=2, time=1.0) def test_azimuth_switchover(): """Test pointing interpolation""" - from ctapipe.io.interpolating import Interpolator + from ctapipe.io.interpolating import PointingInterpolator t0 = Time("2022-01-01T00:00:00") @@ -73,7 +58,7 @@ def test_azimuth_switchover(): }, ) - interpolator = Interpolator() + interpolator = PointingInterpolator() interpolator.add_table(1, table) alt, az = interpolator(tel_id=1, time=t0 + 0.5 * u.s) @@ -83,7 +68,7 @@ def test_azimuth_switchover(): def test_invalid_input(): """Test invalid pointing tables raise nice errors""" - from ctapipe.io.interpolating import Interpolator + from ctapipe.io.interpolating import PointingInterpolator wrong_time = Table( { @@ -93,7 +78,7 @@ def test_invalid_input(): } ) - interpolator = Interpolator() + interpolator = PointingInterpolator() with pytest.raises(TypeError, match="astropy.time.Time"): interpolator.add_table(1, wrong_time) @@ -120,7 +105,7 @@ def test_invalid_input(): def test_hdf5(tmp_path): from ctapipe.io import write_table - from ctapipe.io.interpolating import Interpolator + from ctapipe.io.interpolating import PointingInterpolator t0 = Time("2022-01-01T00:00:00") @@ -135,7 +120,7 @@ def test_hdf5(tmp_path): path = tmp_path / "pointing.h5" write_table(table, path, "/dl0/monitoring/telescope/pointing/tel_001") with tables.open_file(path) as h5file: - interpolator = Interpolator(h5file) + interpolator = PointingInterpolator(h5file) alt, az = interpolator(tel_id=1, time=t0 + 1 * u.s) assert u.isclose(alt, 69 * u.deg) assert u.isclose(az, 1 * u.deg) @@ -143,7 +128,7 @@ def test_hdf5(tmp_path): def test_bounds(): """Test invalid pointing tables raise nice errors""" - from ctapipe.io.interpolating import Interpolator + from ctapipe.io.interpolating import CalibrationInterpolator, PointingInterpolator t0 = Time("2022-01-01T00:00:00") @@ -155,26 +140,17 @@ def test_bounds(): }, ) - table_pedestal = Table( + table_cal = Table( { "time": np.arange(0.0, 10.1, 2.0), "pedestal": np.reshape(np.random.normal(4.0, 1.0, 1850 * 6), (6, 1850)), }, ) - table_gain = Table( - { - "time": np.arange(0.0, 10.1, 2.0), - "gain": np.reshape(np.random.normal(1.0, 0.2, 1850 * 6), (6, 1850)), - }, - ) - - interpolator_pointing = Interpolator() - interpolator_pedestal = Interpolator() - interpolator_gain = Interpolator() + interpolator_pointing = PointingInterpolator() + interpolator_cal = CalibrationInterpolator() interpolator_pointing.add_table(1, table_pointing) - interpolator_pedestal.add_table(1, table_pedestal) - interpolator_gain.add_table(1, table_gain) + interpolator_cal.add_table(1, table_cal, "pedestal") error_message = "below the interpolation range" @@ -182,30 +158,24 @@ def test_bounds(): interpolator_pointing(tel_id=1, time=t0 - 0.1 * u.s) with pytest.raises(ValueError, match=error_message): - interpolator_pedestal(tel_id=1, time=-0.1) - - with pytest.raises(ValueError, match=error_message): - interpolator_gain(tel_id=1, time=-0.1) + interpolator_cal(tel_id=1, time=-0.1) with pytest.raises(ValueError, match="above the interpolation range"): interpolator_pointing(tel_id=1, time=t0 + 10.2 * u.s) - interpolator_pointing = Interpolator(bounds_error=False) - interpolator_pedestal = Interpolator(bounds_error=False) - interpolator_gain = Interpolator(bounds_error=False) + interpolator_pointing = PointingInterpolator(bounds_error=False) + interpolator_cal = CalibrationInterpolator(bounds_error=False) interpolator_pointing.add_table(1, table_pointing) - interpolator_pedestal.add_table(1, table_pedestal) - interpolator_gain.add_table(1, table_gain) + interpolator_cal.add_table(1, table_cal, "pedestal") for dt in (-0.1, 10.1) * u.s: alt, az = interpolator_pointing(tel_id=1, time=t0 + dt) assert np.isnan(alt.value) assert np.isnan(az.value) - assert all(np.isnan(interpolator_pedestal(tel_id=1, time=-0.1))) - assert all(np.isnan(interpolator_gain(tel_id=1, time=-0.1))) + assert all(np.isnan(interpolator_cal(tel_id=1, time=-0.1))) - interpolator = Interpolator(bounds_error=False, extrapolate=True) + interpolator = PointingInterpolator(bounds_error=False, extrapolate=True) interpolator.add_table(1, table_pointing) alt, az = interpolator(tel_id=1, time=t0 - 1 * u.s) assert u.isclose(alt, 71 * u.deg) From 1036d7312575fe052e798edaa5c8167079bbe123 Mon Sep 17 00:00:00 2001 From: Christoph Toennis Date: Wed, 24 Jul 2024 11:31:27 +0200 Subject: [PATCH 08/23] I fixed moved some functions into the parent class --- src/ctapipe/io/interpolating.py | 36 +++++++++++++---------- src/ctapipe/io/tests/test_interpolator.py | 10 ++----- 2 files changed, 22 insertions(+), 24 deletions(-) diff --git a/src/ctapipe/io/interpolating.py b/src/ctapipe/io/interpolating.py index 4b764b5bfe9..2608d5d5756 100644 --- a/src/ctapipe/io/interpolating.py +++ b/src/ctapipe/io/interpolating.py @@ -1,3 +1,4 @@ +from abc import abstractmethod from typing import Any import astropy.units as u @@ -99,15 +100,21 @@ def __init__(self, h5file=None, **kwargs): self._interpolators = {} self._secondary_interpolators = {} + @abstractmethod def add_table(self, tel_id, input_table): - # sort first, so it's not done twice for each interpolator - input_table.sort("time") + pass def __call__(self, tel_id, time): - pass + if tel_id not in self._interpolators: + if self.h5file is not None: + self._read_parameter_table(tel_id) + else: + raise KeyError(f"No table available for tel_id {tel_id}") def _read_parameter_table(self, tel_id): - pass + pass # this method is called when you try to interpolate, but no table has been added. + # It will try adding a table from the hdf5 file specified in __init__ + # Each type of interpolator will need to get the data from a different point in the file class PointingInterpolator(Interpolator): @@ -130,11 +137,7 @@ def __call__(self, tel_id, time): interpolated azimuth angle """ - if tel_id not in self._interpolators: - if self.h5file is not None: - self._read_parameter_table(tel_id) - else: - raise KeyError(f"No table available for tel_id {tel_id}") + super().__call__(tel_id, time) mjd = time.tai.mjd az = u.Quantity(self._interpolators[tel_id](mjd), u.rad, copy=False) @@ -166,6 +169,8 @@ def add_table(self, tel_id, input_table): if not isinstance(input_table["time"], Time): raise TypeError("'time' column of pointing table must be astropy.time.Time") + input_table.sort("time") + az = input_table["azimuth"].quantity.to_value(u.rad) # prepare azimuth for interpolation by "unwrapping": i.e. turning # [359, 1] into [359, 361]. This assumes that if we get values like @@ -199,7 +204,7 @@ def __call__(self, tel_id, time): tel_id : int telescope id time : astropy.time.Time - time for which to interpolate the pointing + time for which to interpolate the calibration data Returns ------- @@ -207,16 +212,12 @@ def __call__(self, tel_id, time): interpolated calibration quantity """ - if tel_id not in self._interpolators: - if self.h5file is not None: - self._read_parameter_table(tel_id) - else: - raise KeyError(f"No table available for tel_id {tel_id}") + super().__call__(tel_id, time) cal = self._interpolators[tel_id](time) return cal - def add_table(self, tel_id, input_table, par_name): + def add_table(self, tel_id, input_table, par_name="pedestal"): """ Add a table to this interpolator @@ -230,6 +231,9 @@ def add_table(self, tel_id, input_table, par_name): parameter column is given through the variable ``par_name`` par_name : str Name of the parameter that is to be interpolated + ``pedestal`` is used for pedestals, ``gain`` for gain + can also be the name of statistical parameters to + interpolate the content of StatisticsContainers """ missing = {"time", par_name} - set(input_table.colnames) diff --git a/src/ctapipe/io/tests/test_interpolator.py b/src/ctapipe/io/tests/test_interpolator.py index b014cd51d98..2d0777ef44b 100644 --- a/src/ctapipe/io/tests/test_interpolator.py +++ b/src/ctapipe/io/tests/test_interpolator.py @@ -5,13 +5,13 @@ from astropy.table import Table from astropy.time import Time +t0 = Time("2022-01-01T00:00:00") + def test_simple(): """Test interpolator""" from ctapipe.io.interpolating import CalibrationInterpolator, PointingInterpolator - t0 = Time("2022-01-01T00:00:00") - table = Table( { "time": t0 + np.arange(0.0, 10.1, 2.0) * u.s, @@ -48,8 +48,6 @@ def test_azimuth_switchover(): """Test pointing interpolation""" from ctapipe.io.interpolating import PointingInterpolator - t0 = Time("2022-01-01T00:00:00") - table = Table( { "time": t0 + [0, 1, 2] * u.s, @@ -107,8 +105,6 @@ def test_hdf5(tmp_path): from ctapipe.io import write_table from ctapipe.io.interpolating import PointingInterpolator - t0 = Time("2022-01-01T00:00:00") - table = Table( { "time": t0 + np.arange(0.0, 10.1, 2.0) * u.s, @@ -130,8 +126,6 @@ def test_bounds(): """Test invalid pointing tables raise nice errors""" from ctapipe.io.interpolating import CalibrationInterpolator, PointingInterpolator - t0 = Time("2022-01-01T00:00:00") - table_pointing = Table( { "time": t0 + np.arange(0.0, 10.1, 2.0) * u.s, From 942f9fd4d7b3c07c1b7575b963551b53ee3f1aca Mon Sep 17 00:00:00 2001 From: Christoph Toennis Date: Wed, 24 Jul 2024 11:58:08 +0200 Subject: [PATCH 09/23] I remove the old pointing interpolator and related tests --- src/ctapipe/io/pointing.py | 137 ----------------------- src/ctapipe/io/tests/test_pointing.py | 154 -------------------------- 2 files changed, 291 deletions(-) delete mode 100644 src/ctapipe/io/pointing.py delete mode 100644 src/ctapipe/io/tests/test_pointing.py diff --git a/src/ctapipe/io/pointing.py b/src/ctapipe/io/pointing.py deleted file mode 100644 index 412c0510a98..00000000000 --- a/src/ctapipe/io/pointing.py +++ /dev/null @@ -1,137 +0,0 @@ -from typing import Any - -import astropy.units as u -import numpy as np -import tables -from astropy.time import Time -from scipy.interpolate import interp1d - -from ctapipe.core import Component, traits - -from .astropy_helpers import read_table - - -class PointingInterpolator(Component): - """ - Interpolate pointing from a monitoring table to a given timestamp. - - Parameters - ---------- - h5file : None | tables.File - A open hdf5 file with read access. - The monitoring table is expected to be stored in that file at - ``/dl0/monitoring/telescope/pointing/tel_{tel_id:03d}`` - - If not given, monitoring tables can be added via `PointingInterpolator.add_table`. - """ - - bounds_error = traits.Bool( - default_value=True, - help="If true, raises an exception when trying to extrapolate out of the given table", - ).tag(config=True) - - extrapolate = traits.Bool( - help="If bounds_error is False, this flag will specify whether values outside" - "the available values are filled with nan (False) or extrapolated (True).", - default_value=False, - ).tag(config=True) - - def __init__(self, h5file=None, **kwargs): - super().__init__(**kwargs) - - if h5file is not None and not isinstance(h5file, tables.File): - raise TypeError("h5file must be a tables.File") - self.h5file = h5file - - self.interp_options: dict[str, Any] = dict(assume_sorted=True, copy=False) - if self.bounds_error: - self.interp_options["bounds_error"] = True - elif self.extrapolate: - self.interp_options["bounds_error"] = False - self.interp_options["fill_value"] = "extrapolate" - else: - self.interp_options["bounds_error"] = False - self.interp_options["fill_value"] = np.nan - - self._alt_interpolators = {} - self._az_interpolators = {} - - def add_table(self, tel_id, pointing_table): - """ - Add a table to this interpolator - - Parameters - ---------- - tel_id : int - Telescope id - pointing_table : astropy.table.Table - Table of pointing values, expected columns - are ``time`` as ``Time`` column, ``azimuth`` and ``altitude`` - as quantity columns. - """ - missing = {"time", "azimuth", "altitude"} - set(pointing_table.colnames) - if len(missing) > 0: - raise ValueError(f"Table is missing required column(s): {missing}") - - if not isinstance(pointing_table["time"], Time): - raise TypeError("'time' column of pointing table must be astropy.time.Time") - - for col in ("azimuth", "altitude"): - unit = pointing_table[col].unit - if unit is None or not u.rad.is_equivalent(unit): - raise ValueError(f"{col} must have units compatible with 'rad'") - - # sort first, so it's not done twice for each interpolator - pointing_table.sort("time") - # interpolate in mjd TAI. Float64 mjd is precise enough for pointing - # and TAI is contiguous, so no issues with leap seconds. - mjd = pointing_table["time"].tai.mjd - - az = pointing_table["azimuth"].quantity.to_value(u.rad) - # prepare azimuth for interpolation by "unwrapping": i.e. turning - # [359, 1] into [359, 361]. This assumes that if we get values like - # [359, 1] the telescope moved 2 degrees through 0, not 358 degrees - # the other way around. This should be true for all telescopes given - # the sampling speed of pointing values and their maximum movement speed. - # No telescope can turn more than 180° in 2 seconds. - az = np.unwrap(az) - alt = pointing_table["altitude"].quantity.to_value(u.rad) - - self._az_interpolators[tel_id] = interp1d(mjd, az, **self.interp_options) - self._alt_interpolators[tel_id] = interp1d(mjd, alt, **self.interp_options) - - def _read_pointing_table(self, tel_id): - pointing_table = read_table( - self.h5file, - f"/dl0/monitoring/telescope/pointing/tel_{tel_id:03d}", - ) - self.add_table(tel_id, pointing_table) - - def __call__(self, tel_id, time): - """ - Interpolate alt/az for given time and tel_id. - - Parameters - ---------- - tel_id : int - telescope id - time : astropy.time.Time - time for which to interpolate the pointing - - Returns - ------- - altitude : astropy.units.Quantity[deg] - interpolated altitude angle - azimuth : astropy.units.Quantity[deg] - interpolated azimuth angle - """ - if tel_id not in self._az_interpolators: - if self.h5file is not None: - self._read_pointing_table(tel_id) - else: - raise KeyError(f"No pointing table available for tel_id {tel_id}") - - mjd = time.tai.mjd - az = u.Quantity(self._az_interpolators[tel_id](mjd), u.rad, copy=False) - alt = u.Quantity(self._alt_interpolators[tel_id](mjd), u.rad, copy=False) - return alt, az diff --git a/src/ctapipe/io/tests/test_pointing.py b/src/ctapipe/io/tests/test_pointing.py deleted file mode 100644 index 10913053c18..00000000000 --- a/src/ctapipe/io/tests/test_pointing.py +++ /dev/null @@ -1,154 +0,0 @@ -import astropy.units as u -import numpy as np -import pytest -import tables -from astropy.table import Table -from astropy.time import Time - - -def test_simple(): - """Test pointing interpolation""" - from ctapipe.io.pointing import PointingInterpolator - - t0 = Time("2022-01-01T00:00:00") - - table = Table( - { - "time": t0 + np.arange(0.0, 10.1, 2.0) * u.s, - "azimuth": np.linspace(0.0, 10.0, 6) * u.deg, - "altitude": np.linspace(70.0, 60.0, 6) * u.deg, - }, - ) - - interpolator = PointingInterpolator() - interpolator.add_table(1, table) - - alt, az = interpolator(tel_id=1, time=t0 + 1 * u.s) - assert u.isclose(alt, 69 * u.deg) - assert u.isclose(az, 1 * u.deg) - - with pytest.raises(KeyError): - interpolator(tel_id=2, time=t0 + 1 * u.s) - - -def test_azimuth_switchover(): - """Test pointing interpolation""" - from ctapipe.io.pointing import PointingInterpolator - - t0 = Time("2022-01-01T00:00:00") - - table = Table( - { - "time": t0 + [0, 1, 2] * u.s, - "azimuth": [359, 1, 3] * u.deg, - "altitude": [60, 61, 62] * u.deg, - }, - ) - - interpolator = PointingInterpolator() - interpolator.add_table(1, table) - - alt, az = interpolator(tel_id=1, time=t0 + 0.5 * u.s) - assert u.isclose(az, 360 * u.deg) - assert u.isclose(alt, 60.5 * u.deg) - - -def test_invalid_input(): - """Test invalid pointing tables raise nice errors""" - from ctapipe.io.pointing import PointingInterpolator - - wrong_time = Table( - { - "time": [1, 2, 3] * u.s, - "azimuth": [1, 2, 3] * u.deg, - "altitude": [1, 2, 3] * u.deg, - } - ) - - interpolator = PointingInterpolator() - with pytest.raises(TypeError, match="astropy.time.Time"): - interpolator.add_table(1, wrong_time) - - wrong_unit = Table( - { - "time": Time(1.7e9 + np.arange(3), format="unix"), - "azimuth": [1, 2, 3] * u.m, - "altitude": [1, 2, 3] * u.deg, - } - ) - with pytest.raises(ValueError, match="compatible with 'rad'"): - interpolator.add_table(1, wrong_unit) - - wrong_unit = Table( - { - "time": Time(1.7e9 + np.arange(3), format="unix"), - "azimuth": [1, 2, 3] * u.deg, - "altitude": [1, 2, 3], - } - ) - with pytest.raises(ValueError, match="compatible with 'rad'"): - interpolator.add_table(1, wrong_unit) - - -def test_hdf5(tmp_path): - from ctapipe.io import write_table - from ctapipe.io.pointing import PointingInterpolator - - t0 = Time("2022-01-01T00:00:00") - - table = Table( - { - "time": t0 + np.arange(0.0, 10.1, 2.0) * u.s, - "azimuth": np.linspace(0.0, 10.0, 6) * u.deg, - "altitude": np.linspace(70.0, 60.0, 6) * u.deg, - }, - ) - - path = tmp_path / "pointing.h5" - write_table(table, path, "/dl0/monitoring/telescope/pointing/tel_001") - with tables.open_file(path) as h5file: - interpolator = PointingInterpolator(h5file) - alt, az = interpolator(tel_id=1, time=t0 + 1 * u.s) - assert u.isclose(alt, 69 * u.deg) - assert u.isclose(az, 1 * u.deg) - - -def test_bounds(): - """Test invalid pointing tables raise nice errors""" - from ctapipe.io.pointing import PointingInterpolator - - t0 = Time("2022-01-01T00:00:00") - - table = Table( - { - "time": t0 + np.arange(0.0, 10.1, 2.0) * u.s, - "azimuth": np.linspace(0.0, 10.0, 6) * u.deg, - "altitude": np.linspace(70.0, 60.0, 6) * u.deg, - }, - ) - - interpolator = PointingInterpolator() - interpolator.add_table(1, table) - - with pytest.raises(ValueError, match="below the interpolation range"): - interpolator(tel_id=1, time=t0 - 0.1 * u.s) - - with pytest.raises(ValueError, match="above the interpolation range"): - interpolator(tel_id=1, time=t0 + 10.2 * u.s) - - interpolator = PointingInterpolator(bounds_error=False) - interpolator.add_table(1, table) - for dt in (-0.1, 10.1) * u.s: - alt, az = interpolator(tel_id=1, time=t0 + dt) - assert np.isnan(alt.value) - assert np.isnan(az.value) - - interpolator = PointingInterpolator(bounds_error=False, extrapolate=True) - interpolator.add_table(1, table) - alt, az = interpolator(tel_id=1, time=t0 - 1 * u.s) - assert u.isclose(alt, 71 * u.deg) - assert u.isclose(az, -1 * u.deg) - - alt, az = interpolator(tel_id=1, time=t0 + 11 * u.s) - assert u.isclose(alt, 59 * u.deg) - assert u.isclose(az, 11 * u.deg) From c4e6d5370ed13f61eb13eb5a5e0bb91bf27de324 Mon Sep 17 00:00:00 2001 From: Christoph Toennis Date: Wed, 24 Jul 2024 13:00:11 +0200 Subject: [PATCH 10/23] I changed the text to remove duplicate code --- src/ctapipe/io/interpolating.py | 11 +++- src/ctapipe/io/tests/test_interpolator.py | 63 +++++++---------------- 2 files changed, 28 insertions(+), 46 deletions(-) diff --git a/src/ctapipe/io/interpolating.py b/src/ctapipe/io/interpolating.py index 2608d5d5756..afabc513e54 100644 --- a/src/ctapipe/io/interpolating.py +++ b/src/ctapipe/io/interpolating.py @@ -107,14 +107,14 @@ def add_table(self, tel_id, input_table): def __call__(self, tel_id, time): if tel_id not in self._interpolators: if self.h5file is not None: - self._read_parameter_table(tel_id) + self._read_parameter_table(tel_id) # might need to be removed else: raise KeyError(f"No table available for tel_id {tel_id}") + @abstractmethod def _read_parameter_table(self, tel_id): pass # this method is called when you try to interpolate, but no table has been added. # It will try adding a table from the hdf5 file specified in __init__ - # Each type of interpolator will need to get the data from a different point in the file class PointingInterpolator(Interpolator): @@ -244,3 +244,10 @@ def add_table(self, tel_id, input_table, par_name="pedestal"): time = input_table["time"] cal = input_table[par_name] self._interpolators[tel_id] = StepFunction(time, cal, **self.interp_options) + + def _read_parameter_table(self, tel_id): + input_table = read_table( + self.h5file, + f"/dl0/calibration/tel_{tel_id:03d}", # needs to be updated once we determine where calibration data is put, might remove method from base class + ) + self.add_table(tel_id, input_table) diff --git a/src/ctapipe/io/tests/test_interpolator.py b/src/ctapipe/io/tests/test_interpolator.py index 2d0777ef44b..1fa9ca52972 100644 --- a/src/ctapipe/io/tests/test_interpolator.py +++ b/src/ctapipe/io/tests/test_interpolator.py @@ -5,48 +5,13 @@ from astropy.table import Table from astropy.time import Time -t0 = Time("2022-01-01T00:00:00") - - -def test_simple(): - """Test interpolator""" - from ctapipe.io.interpolating import CalibrationInterpolator, PointingInterpolator - - table = Table( - { - "time": t0 + np.arange(0.0, 10.1, 2.0) * u.s, - "azimuth": np.linspace(0.0, 10.0, 6) * u.deg, - "altitude": np.linspace(70.0, 60.0, 6) * u.deg, - }, - ) - - table_cal = Table( - { - "time": np.arange(0.0, 10.1, 2.0), - "pedestal": np.reshape(np.random.normal(4.0, 1.0, 1850 * 6), (6, 1850)), - }, - ) - - interpolator = PointingInterpolator() - interpolator_cal = CalibrationInterpolator() - interpolator.add_table(1, table) - interpolator_cal.add_table(1, table_cal, "pedestal") +from ctapipe.io.interpolating import CalibrationInterpolator, PointingInterpolator - alt, az = interpolator(tel_id=1, time=t0 + 1 * u.s) - assert u.isclose(alt, 69 * u.deg) - assert u.isclose(az, 1 * u.deg) - - pedestal = interpolator_cal(tel_id=1, time=1.0) - assert all(pedestal == table_cal["pedestal"][0]) - with pytest.raises(KeyError): - interpolator(tel_id=2, time=t0 + 1 * u.s) - with pytest.raises(KeyError): - interpolator_cal(tel_id=2, time=1.0) +t0 = Time("2022-01-01T00:00:00") def test_azimuth_switchover(): """Test pointing interpolation""" - from ctapipe.io.interpolating import PointingInterpolator table = Table( { @@ -66,7 +31,6 @@ def test_azimuth_switchover(): def test_invalid_input(): """Test invalid pointing tables raise nice errors""" - from ctapipe.io.interpolating import PointingInterpolator wrong_time = Table( { @@ -102,8 +66,8 @@ def test_invalid_input(): def test_hdf5(tmp_path): + """Test writing interpolated data to file""" from ctapipe.io import write_table - from ctapipe.io.interpolating import PointingInterpolator table = Table( { @@ -124,7 +88,7 @@ def test_hdf5(tmp_path): def test_bounds(): """Test invalid pointing tables raise nice errors""" - from ctapipe.io.interpolating import CalibrationInterpolator, PointingInterpolator + from ctapipe.io.interpolating import PointingInterpolator table_pointing = Table( { @@ -157,6 +121,17 @@ def test_bounds(): with pytest.raises(ValueError, match="above the interpolation range"): interpolator_pointing(tel_id=1, time=t0 + 10.2 * u.s) + alt, az = interpolator_pointing(tel_id=1, time=t0 + 1 * u.s) + assert u.isclose(alt, 69 * u.deg) + assert u.isclose(az, 1 * u.deg) + + pedestal = interpolator_cal(tel_id=1, time=1.0) + assert all(pedestal == table_cal["pedestal"][0]) + with pytest.raises(KeyError): + interpolator_pointing(tel_id=2, time=t0 + 1 * u.s) + with pytest.raises(KeyError): + interpolator_cal(tel_id=2, time=1.0) + interpolator_pointing = PointingInterpolator(bounds_error=False) interpolator_cal = CalibrationInterpolator(bounds_error=False) interpolator_pointing.add_table(1, table_pointing) @@ -169,12 +144,12 @@ def test_bounds(): assert all(np.isnan(interpolator_cal(tel_id=1, time=-0.1))) - interpolator = PointingInterpolator(bounds_error=False, extrapolate=True) - interpolator.add_table(1, table_pointing) - alt, az = interpolator(tel_id=1, time=t0 - 1 * u.s) + interpolator_pointing = PointingInterpolator(bounds_error=False, extrapolate=True) + interpolator_pointing.add_table(1, table_pointing) + alt, az = interpolator_pointing(tel_id=1, time=t0 - 1 * u.s) assert u.isclose(alt, 71 * u.deg) assert u.isclose(az, -1 * u.deg) - alt, az = interpolator(tel_id=1, time=t0 + 11 * u.s) + alt, az = interpolator_pointing(tel_id=1, time=t0 + 11 * u.s) assert u.isclose(alt, 59 * u.deg) assert u.isclose(az, 11 * u.deg) From dce4ef24a05e26823c239ad313cf9d6f7a7b66ab Mon Sep 17 00:00:00 2001 From: Christoph Toennis Date: Thu, 25 Jul 2024 10:26:20 +0200 Subject: [PATCH 11/23] Implementing recent comments and change requests --- src/ctapipe/io/__init__.py | 3 ++ src/ctapipe/io/hdf5eventsource.py | 2 +- .../io/{interpolating.py => interpolation.py} | 54 +++++++++---------- src/ctapipe/io/tableloader.py | 2 +- src/ctapipe/io/tests/test_interpolator.py | 3 +- 5 files changed, 32 insertions(+), 32 deletions(-) rename src/ctapipe/io/{interpolating.py => interpolation.py} (85%) diff --git a/src/ctapipe/io/__init__.py b/src/ctapipe/io/__init__.py index afe96f39430..f6e7a0b12ba 100644 --- a/src/ctapipe/io/__init__.py +++ b/src/ctapipe/io/__init__.py @@ -18,6 +18,7 @@ from .datawriter import DATA_MODEL_VERSION, DataWriter +from .interpolation import PointingInterpolator, CalibrationInterpolator __all__ = [ "HDF5TableWriter", @@ -36,4 +37,6 @@ "DataWriter", "DATA_MODEL_VERSION", "get_hdf5_datalevels", + "PointingInterpolator", + "CalibrationInterpolator", ] diff --git a/src/ctapipe/io/hdf5eventsource.py b/src/ctapipe/io/hdf5eventsource.py index e7bcb1a5073..eb91b216c27 100644 --- a/src/ctapipe/io/hdf5eventsource.py +++ b/src/ctapipe/io/hdf5eventsource.py @@ -52,7 +52,7 @@ from .datalevels import DataLevel from .eventsource import EventSource from .hdf5tableio import HDF5TableReader -from .interpolating import PointingInterpolator +from .interpolation import PointingInterpolator from .tableloader import DL2_SUBARRAY_GROUP, DL2_TELESCOPE_GROUP, POINTING_GROUP __all__ = ["HDF5EventSource"] diff --git a/src/ctapipe/io/interpolating.py b/src/ctapipe/io/interpolation.py similarity index 85% rename from src/ctapipe/io/interpolating.py rename to src/ctapipe/io/interpolation.py index afabc513e54..cc4a091fc79 100644 --- a/src/ctapipe/io/interpolating.py +++ b/src/ctapipe/io/interpolation.py @@ -80,6 +80,8 @@ class Interpolator(Component): default_value=False, ).tag(config=True) + table_location = None + def __init__(self, h5file=None, **kwargs): super().__init__(**kwargs) @@ -98,26 +100,30 @@ def __init__(self, h5file=None, **kwargs): self.interp_options["fill_value"] = np.nan self._interpolators = {} - self._secondary_interpolators = {} @abstractmethod def add_table(self, tel_id, input_table): pass - def __call__(self, tel_id, time): + def _check_interpolators(self, tel_id): if tel_id not in self._interpolators: if self.h5file is not None: self._read_parameter_table(tel_id) # might need to be removed else: raise KeyError(f"No table available for tel_id {tel_id}") - @abstractmethod def _read_parameter_table(self, tel_id): - pass # this method is called when you try to interpolate, but no table has been added. - # It will try adding a table from the hdf5 file specified in __init__ + input_table = read_table( + self.h5file, + self.table_location + f"/tel_{tel_id:03d}", + ) + print(self.table_location + f"/tel_{tel_id:03d}") + self.add_table(tel_id, input_table) class PointingInterpolator(Interpolator): + table_location = "/dl0/monitoring/telescope/pointing" + def __call__(self, tel_id, time): """ Interpolate alt/az for given time and tel_id. @@ -137,11 +143,11 @@ def __call__(self, tel_id, time): interpolated azimuth angle """ - super().__call__(tel_id, time) + self._check_interpolators(tel_id) mjd = time.tai.mjd - az = u.Quantity(self._interpolators[tel_id](mjd), u.rad, copy=False) - alt = u.Quantity(self._secondary_interpolators[tel_id](mjd), u.rad, copy=False) + az = u.Quantity(self._interpolators[tel_id]["az"](mjd), u.rad, copy=False) + alt = u.Quantity(self._interpolators[tel_id]["alt"](mjd), u.rad, copy=False) return alt, az def add_table(self, tel_id, input_table): @@ -181,20 +187,14 @@ def add_table(self, tel_id, input_table): az = np.unwrap(az) alt = input_table["altitude"].quantity.to_value(u.rad) mjd = input_table["time"].tai.mjd - self._interpolators[tel_id] = interp1d(mjd, az, **self.interp_options) - self._secondary_interpolators[tel_id] = interp1d( - mjd, alt, **self.interp_options - ) - - def _read_parameter_table(self, tel_id): - input_table = read_table( - self.h5file, - f"/dl0/monitoring/telescope/pointing/tel_{tel_id:03d}", - ) - self.add_table(tel_id, input_table) + self._interpolators[tel_id] = {} + self._interpolators[tel_id]["az"] = interp1d(mjd, az, **self.interp_options) + self._interpolators[tel_id]["alt"] = interp1d(mjd, alt, **self.interp_options) class CalibrationInterpolator(Interpolator): + table_location = "dl1/calibration" # TBD + def __call__(self, tel_id, time): """ Interpolate pedestal or gain for a given time and tel_id. @@ -212,9 +212,9 @@ def __call__(self, tel_id, time): interpolated calibration quantity """ - super().__call__(tel_id, time) + self._check_interpolators(tel_id) - cal = self._interpolators[tel_id](time) + cal = self._interpolators[tel_id][self.par_name](time) return cal def add_table(self, tel_id, input_table, par_name="pedestal"): @@ -240,14 +240,12 @@ def add_table(self, tel_id, input_table, par_name="pedestal"): if len(missing) > 0: raise ValueError(f"Table is missing required column(s): {missing}") + self.par_name = par_name + input_table.sort("time") time = input_table["time"] cal = input_table[par_name] - self._interpolators[tel_id] = StepFunction(time, cal, **self.interp_options) - - def _read_parameter_table(self, tel_id): - input_table = read_table( - self.h5file, - f"/dl0/calibration/tel_{tel_id:03d}", # needs to be updated once we determine where calibration data is put, might remove method from base class + self._interpolators[tel_id] = {} + self._interpolators[tel_id][par_name] = StepFunction( + time, cal, **self.interp_options ) - self.add_table(tel_id, input_table) diff --git a/src/ctapipe/io/tableloader.py b/src/ctapipe/io/tableloader.py index 771058ffad7..6da74b6f081 100644 --- a/src/ctapipe/io/tableloader.py +++ b/src/ctapipe/io/tableloader.py @@ -14,7 +14,7 @@ from ..core import Component, Provenance, traits from ..instrument import FocalLengthKind, SubarrayDescription from .astropy_helpers import join_allow_empty, read_table -from .interpolating import PointingInterpolator +from .interpolation import PointingInterpolator __all__ = ["TableLoader"] diff --git a/src/ctapipe/io/tests/test_interpolator.py b/src/ctapipe/io/tests/test_interpolator.py index 1fa9ca52972..9883069893d 100644 --- a/src/ctapipe/io/tests/test_interpolator.py +++ b/src/ctapipe/io/tests/test_interpolator.py @@ -5,7 +5,7 @@ from astropy.table import Table from astropy.time import Time -from ctapipe.io.interpolating import CalibrationInterpolator, PointingInterpolator +from ctapipe.io.interpolation import CalibrationInterpolator, PointingInterpolator t0 = Time("2022-01-01T00:00:00") @@ -88,7 +88,6 @@ def test_hdf5(tmp_path): def test_bounds(): """Test invalid pointing tables raise nice errors""" - from ctapipe.io.interpolating import PointingInterpolator table_pointing = Table( { From 681ecca3b9070b858d7967dea5c4f669dcc4aae3 Mon Sep 17 00:00:00 2001 From: Christoph Toennis Date: Thu, 25 Jul 2024 10:45:21 +0200 Subject: [PATCH 12/23] Adding some docustrings --- src/ctapipe/io/interpolation.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/ctapipe/io/interpolation.py b/src/ctapipe/io/interpolation.py index cc4a091fc79..9183bfb3c89 100644 --- a/src/ctapipe/io/interpolation.py +++ b/src/ctapipe/io/interpolation.py @@ -122,6 +122,10 @@ def _read_parameter_table(self, tel_id): class PointingInterpolator(Interpolator): + """ + Interpolator for pointing and pointing correction data + """ + table_location = "/dl0/monitoring/telescope/pointing" def __call__(self, tel_id, time): @@ -193,6 +197,10 @@ def add_table(self, tel_id, input_table): class CalibrationInterpolator(Interpolator): + """ + Interpolator for calibration data + """ + table_location = "dl1/calibration" # TBD def __call__(self, tel_id, time): From 29c361cc6ca15319fbf673b059851d9f18b8f0d0 Mon Sep 17 00:00:00 2001 From: Christoph Toennis Date: Thu, 25 Jul 2024 11:11:24 +0200 Subject: [PATCH 13/23] Added some import in __init__ --- src/ctapipe/io/__init__.py | 3 ++- src/ctapipe/io/interpolation.py | 18 +++++++++--------- 2 files changed, 11 insertions(+), 10 deletions(-) diff --git a/src/ctapipe/io/__init__.py b/src/ctapipe/io/__init__.py index f6e7a0b12ba..baa2ec2e339 100644 --- a/src/ctapipe/io/__init__.py +++ b/src/ctapipe/io/__init__.py @@ -18,7 +18,7 @@ from .datawriter import DATA_MODEL_VERSION, DataWriter -from .interpolation import PointingInterpolator, CalibrationInterpolator +from .interpolation import Interpolator, PointingInterpolator, CalibrationInterpolator __all__ = [ "HDF5TableWriter", @@ -37,6 +37,7 @@ "DataWriter", "DATA_MODEL_VERSION", "get_hdf5_datalevels", + "Interpolator", "PointingInterpolator", "CalibrationInterpolator", ] diff --git a/src/ctapipe/io/interpolation.py b/src/ctapipe/io/interpolation.py index 9183bfb3c89..e8ecff92f5b 100644 --- a/src/ctapipe/io/interpolation.py +++ b/src/ctapipe/io/interpolation.py @@ -60,15 +60,6 @@ def __call__(self, point): class Interpolator(Component): - """ - Interpolator parent class. - - Parameters - ---------- - h5file : None | tables.File - A open hdf5 file with read access. - """ - bounds_error = traits.Bool( default_value=True, help="If true, raises an exception when trying to extrapolate out of the given table", @@ -83,6 +74,15 @@ class Interpolator(Component): table_location = None def __init__(self, h5file=None, **kwargs): + """ + Interpolator parent class. + + Parameters + ---------- + h5file : None | tables.File + A open hdf5 file with read access. + """ + super().__init__(**kwargs) if h5file is not None and not isinstance(h5file, tables.File): From 5a05359c386af86c68ef9d0c44e5f7deb7865cd1 Mon Sep 17 00:00:00 2001 From: Christoph Toennis Date: Thu, 25 Jul 2024 11:31:27 +0200 Subject: [PATCH 14/23] moving docustring --- src/ctapipe/io/interpolation.py | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/ctapipe/io/interpolation.py b/src/ctapipe/io/interpolation.py index e8ecff92f5b..9183bfb3c89 100644 --- a/src/ctapipe/io/interpolation.py +++ b/src/ctapipe/io/interpolation.py @@ -60,6 +60,15 @@ def __call__(self, point): class Interpolator(Component): + """ + Interpolator parent class. + + Parameters + ---------- + h5file : None | tables.File + A open hdf5 file with read access. + """ + bounds_error = traits.Bool( default_value=True, help="If true, raises an exception when trying to extrapolate out of the given table", @@ -74,15 +83,6 @@ class Interpolator(Component): table_location = None def __init__(self, h5file=None, **kwargs): - """ - Interpolator parent class. - - Parameters - ---------- - h5file : None | tables.File - A open hdf5 file with read access. - """ - super().__init__(**kwargs) if h5file is not None and not isinstance(h5file, tables.File): From 52e71604edf9c0bcc7fc24bf247e85d33d8396a1 Mon Sep 17 00:00:00 2001 From: Christoph Toennis Date: Thu, 25 Jul 2024 14:23:51 +0200 Subject: [PATCH 15/23] Implementing latest comments --- src/ctapipe/io/interpolation.py | 54 ++++++++++++++++++++++++--------- 1 file changed, 39 insertions(+), 15 deletions(-) diff --git a/src/ctapipe/io/interpolation.py b/src/ctapipe/io/interpolation.py index 9183bfb3c89..ba41d7f24ea 100644 --- a/src/ctapipe/io/interpolation.py +++ b/src/ctapipe/io/interpolation.py @@ -1,4 +1,4 @@ -from abc import abstractmethod +from abc import ABCMeta, abstractmethod from typing import Any import astropy.units as u @@ -59,7 +59,7 @@ def __call__(self, point): return self.values[i - 1] -class Interpolator(Component): +class Interpolator(Component, metaclass=ABCMeta): """ Interpolator parent class. @@ -81,6 +81,8 @@ class Interpolator(Component): ).tag(config=True) table_location = None + required_columns = set() + expected_units = {} def __init__(self, h5file=None, **kwargs): super().__init__(**kwargs) @@ -103,8 +105,35 @@ def __init__(self, h5file=None, **kwargs): @abstractmethod def add_table(self, tel_id, input_table): + """ + Add a table to this interpolator + This method reads input tables and creates instances of the needed interpolators + to be added to _interpolators. The first index of _interpolators needs to be + tel_id, the second needs to be the name of the parameter that is to be interpolated + + Parameters + ---------- + tel_id : int + Telescope id + input_table : astropy.table.Table + Table of pointing values, expected columns + are always ``time`` as ``Time`` column and other columns for the data + that is to be interpolated + """ + pass + def _check_tables(self, input_table): + missing = self.required_columns - set(input_table.colnames) + if len(missing) > 0: + raise ValueError(f"Table is missing required column(s): {missing}") + for col in self.expected_units: + unit = input_table[col].unit + if unit is None or not self.expected_units[col].is_equivalent(unit): + raise ValueError( + f"{col} must have units compatible with '{self.expected_units[col].name}'" + ) + def _check_interpolators(self, tel_id): if tel_id not in self._interpolators: if self.h5file is not None: @@ -115,9 +144,8 @@ def _check_interpolators(self, tel_id): def _read_parameter_table(self, tel_id): input_table = read_table( self.h5file, - self.table_location + f"/tel_{tel_id:03d}", + f"{self.table_location}/tel_{tel_id:03d}", ) - print(self.table_location + f"/tel_{tel_id:03d}") self.add_table(tel_id, input_table) @@ -127,6 +155,8 @@ class PointingInterpolator(Interpolator): """ table_location = "/dl0/monitoring/telescope/pointing" + required_columns = {"time", "azimuth", "altitude"} + expected_units = {"azimuth": u.rad, "altitude": u.rad} def __call__(self, tel_id, time): """ @@ -168,13 +198,7 @@ def add_table(self, tel_id, input_table): as quantity columns for pointing and pointing correction data. """ - missing = {"time", "azimuth", "altitude"} - set(input_table.colnames) - if len(missing) > 0: - raise ValueError(f"Table is missing required column(s): {missing}") - for col in ("azimuth", "altitude"): - unit = input_table[col].unit - if unit is None or not u.rad.is_equivalent(unit): - raise ValueError(f"{col} must have units compatible with 'rad'") + self._check_tables(input_table) if not isinstance(input_table["time"], Time): raise TypeError("'time' column of pointing table must be astropy.time.Time") @@ -244,12 +268,12 @@ def add_table(self, tel_id, input_table, par_name="pedestal"): interpolate the content of StatisticsContainers """ - missing = {"time", par_name} - set(input_table.colnames) - if len(missing) > 0: - raise ValueError(f"Table is missing required column(s): {missing}") - self.par_name = par_name + self.required_columns = {"time", par_name} + + self._check_tables(input_table) + input_table.sort("time") time = input_table["time"] cal = input_table[par_name] From 8680178d1bde8f02af4f38623097e1597e26fe8d Mon Sep 17 00:00:00 2001 From: Christoph Toennis Date: Wed, 31 Jul 2024 09:53:59 +0200 Subject: [PATCH 16/23] I added the option to require interpolated dat not to have units --- src/ctapipe/io/interpolation.py | 18 +++++++++++++----- 1 file changed, 13 insertions(+), 5 deletions(-) diff --git a/src/ctapipe/io/interpolation.py b/src/ctapipe/io/interpolation.py index ba41d7f24ea..e92e7278460 100644 --- a/src/ctapipe/io/interpolation.py +++ b/src/ctapipe/io/interpolation.py @@ -112,7 +112,7 @@ def add_table(self, tel_id, input_table): tel_id, the second needs to be the name of the parameter that is to be interpolated Parameters - ---------- + -.--------- tel_id : int Telescope id input_table : astropy.table.Table @@ -129,10 +129,18 @@ def _check_tables(self, input_table): raise ValueError(f"Table is missing required column(s): {missing}") for col in self.expected_units: unit = input_table[col].unit - if unit is None or not self.expected_units[col].is_equivalent(unit): - raise ValueError( - f"{col} must have units compatible with '{self.expected_units[col].name}'" - ) + if unit is None: + if self.expected_units[col] is not None: + raise ValueError( + f"{col} must have units compatible with '{self.expected_units[col].name}'" + ) + elif not self.expected_units[col].is_equivalent(unit): + if self.expected_units[col] is None: + raise ValueError(f"{col} must have units compatible with 'None'") + else: + raise ValueError( + f"{col} must have units compatible with '{self.expected_units[col].name}'" + ) def _check_interpolators(self, tel_id): if tel_id not in self._interpolators: From db1c983ddbeb629f09917db49ab5e3cf39fe4e0d Mon Sep 17 00:00:00 2001 From: Christoph Toennis Date: Wed, 31 Jul 2024 14:55:26 +0200 Subject: [PATCH 17/23] I split the CalibrationInterpolator --- src/ctapipe/io/__init__.py | 10 ++- src/ctapipe/io/interpolation.py | 97 +++++++++++++++++------ src/ctapipe/io/tests/test_interpolator.py | 45 ++++++++--- 3 files changed, 116 insertions(+), 36 deletions(-) diff --git a/src/ctapipe/io/__init__.py b/src/ctapipe/io/__init__.py index baa2ec2e339..dbf555a7bd6 100644 --- a/src/ctapipe/io/__init__.py +++ b/src/ctapipe/io/__init__.py @@ -18,7 +18,12 @@ from .datawriter import DATA_MODEL_VERSION, DataWriter -from .interpolation import Interpolator, PointingInterpolator, CalibrationInterpolator +from .interpolation import ( + Interpolator, + PointingInterpolator, + GainInterpolator, + PedestalInterpolator, +) __all__ = [ "HDF5TableWriter", @@ -39,5 +44,6 @@ "get_hdf5_datalevels", "Interpolator", "PointingInterpolator", - "CalibrationInterpolator", + "PedestalInterpolator", + "GainInterpolator", ] diff --git a/src/ctapipe/io/interpolation.py b/src/ctapipe/io/interpolation.py index e92e7278460..dc0953e3892 100644 --- a/src/ctapipe/io/interpolation.py +++ b/src/ctapipe/io/interpolation.py @@ -16,6 +16,8 @@ class StepFunction: """ Step function Interpolator for the gain and pedestals + Interpolates data so that for each time the value from the closest previous + point given. Parameters ---------- @@ -80,7 +82,7 @@ class Interpolator(Component, metaclass=ABCMeta): default_value=False, ).tag(config=True) - table_location = None + telescope_data_group = None required_columns = set() expected_units = {} @@ -152,7 +154,7 @@ def _check_interpolators(self, tel_id): def _read_parameter_table(self, tel_id): input_table = read_table( self.h5file, - f"{self.table_location}/tel_{tel_id:03d}", + f"{self.telescope_data_group}/tel_{tel_id:03d}", ) self.add_table(tel_id, input_table) @@ -162,7 +164,7 @@ class PointingInterpolator(Interpolator): Interpolator for pointing and pointing correction data """ - table_location = "/dl0/monitoring/telescope/pointing" + telescope_data_group = "/dl0/monitoring/telescope/pointing" required_columns = {"time", "azimuth", "altitude"} expected_units = {"azimuth": u.rad, "altitude": u.rad} @@ -228,12 +230,12 @@ def add_table(self, tel_id, input_table): self._interpolators[tel_id]["alt"] = interp1d(mjd, alt, **self.interp_options) -class CalibrationInterpolator(Interpolator): +class GainInterpolator(Interpolator): """ - Interpolator for calibration data + Interpolator for relative gain data """ - table_location = "dl1/calibration" # TBD + telescope_data_group = "dl1/calibration/gain" # TBD def __call__(self, tel_id, time): """ @@ -248,16 +250,16 @@ def __call__(self, tel_id, time): Returns ------- - cal : array [float] - interpolated calibration quantity + gain : array [float] + interpolated relative gain """ self._check_interpolators(tel_id) - cal = self._interpolators[tel_id][self.par_name](time) - return cal + gain = self._interpolators[tel_id]["gain"](time) + return gain - def add_table(self, tel_id, input_table, par_name="pedestal"): + def add_table(self, tel_id, input_table): """ Add a table to this interpolator @@ -267,25 +269,74 @@ def add_table(self, tel_id, input_table, par_name="pedestal"): Telescope id input_table : astropy.table.Table Table of pointing values, expected columns - are ``time`` as ``Time`` column. The calibration - parameter column is given through the variable ``par_name`` - par_name : str - Name of the parameter that is to be interpolated - ``pedestal`` is used for pedestals, ``gain`` for gain - can also be the name of statistical parameters to - interpolate the content of StatisticsContainers + are ``time`` as ``Time`` column and "gain" + for the relative gain data + """ + + self.required_columns = {"time", "gain"} + + self._check_tables(input_table) + + input_table.sort("time") + time = input_table["time"] + gain = input_table["gain"] + self._interpolators[tel_id] = {} + self._interpolators[tel_id]["gain"] = StepFunction( + time, gain, **self.interp_options + ) + + +class PedestalInterpolator(Interpolator): + """ + Interpolator for Pedestal data + """ + + telescope_data_group = "dl1/calibration/pedestal" # TBD + + def __call__(self, tel_id, time): """ + Interpolate pedestal or gain for a given time and tel_id. - self.par_name = par_name + Parameters + ---------- + tel_id : int + telescope id + time : astropy.time.Time + time for which to interpolate the calibration data + + Returns + ------- + pedestal : array [float] + interpolated pedestal values + """ + + self._check_interpolators(tel_id) + + pedestal = self._interpolators[tel_id]["pedestal"](time) + return pedestal + + def add_table(self, tel_id, input_table): + """ + Add a table to this interpolator + + Parameters + ---------- + tel_id : int + Telescope id + input_table : astropy.table.Table + Table of pointing values, expected columns + are ``time`` as ``Time`` column and "pedestal" + for the pedestal data + """ - self.required_columns = {"time", par_name} + self.required_columns = {"time", "pedestal"} self._check_tables(input_table) input_table.sort("time") time = input_table["time"] - cal = input_table[par_name] + pedestal = input_table["pedestal"] self._interpolators[tel_id] = {} - self._interpolators[tel_id][par_name] = StepFunction( - time, cal, **self.interp_options + self._interpolators[tel_id]["pedestal"] = StepFunction( + time, pedestal, **self.interp_options ) diff --git a/src/ctapipe/io/tests/test_interpolator.py b/src/ctapipe/io/tests/test_interpolator.py index 9883069893d..1279dc5c784 100644 --- a/src/ctapipe/io/tests/test_interpolator.py +++ b/src/ctapipe/io/tests/test_interpolator.py @@ -5,7 +5,11 @@ from astropy.table import Table from astropy.time import Time -from ctapipe.io.interpolation import CalibrationInterpolator, PointingInterpolator +from ctapipe.io.interpolation import ( + GainInterpolator, + PedestalInterpolator, + PointingInterpolator, +) t0 = Time("2022-01-01T00:00:00") @@ -97,17 +101,26 @@ def test_bounds(): }, ) - table_cal = Table( + table_pedestal = Table( { "time": np.arange(0.0, 10.1, 2.0), "pedestal": np.reshape(np.random.normal(4.0, 1.0, 1850 * 6), (6, 1850)), }, ) + table_gain = Table( + { + "time": np.arange(0.0, 10.1, 2.0), + "gain": np.reshape(np.random.normal(1.0, 1.0, 1850 * 6), (6, 1850)), + }, + ) + interpolator_pointing = PointingInterpolator() - interpolator_cal = CalibrationInterpolator() + interpolator_pedestal = PedestalInterpolator() + interpolator_gain = GainInterpolator() interpolator_pointing.add_table(1, table_pointing) - interpolator_cal.add_table(1, table_cal, "pedestal") + interpolator_pedestal.add_table(1, table_pedestal) + interpolator_gain.add_table(1, table_gain) error_message = "below the interpolation range" @@ -115,7 +128,10 @@ def test_bounds(): interpolator_pointing(tel_id=1, time=t0 - 0.1 * u.s) with pytest.raises(ValueError, match=error_message): - interpolator_cal(tel_id=1, time=-0.1) + interpolator_pedestal(tel_id=1, time=-0.1) + + with pytest.raises(ValueError, match=error_message): + interpolator_gain(tel_id=1, time=-0.1) with pytest.raises(ValueError, match="above the interpolation range"): interpolator_pointing(tel_id=1, time=t0 + 10.2 * u.s) @@ -124,24 +140,31 @@ def test_bounds(): assert u.isclose(alt, 69 * u.deg) assert u.isclose(az, 1 * u.deg) - pedestal = interpolator_cal(tel_id=1, time=1.0) - assert all(pedestal == table_cal["pedestal"][0]) + pedestal = interpolator_pedestal(tel_id=1, time=1.0) + assert all(pedestal == table_pedestal["pedestal"][0]) + gain = interpolator_gain(tel_id=1, time=1.0) + assert all(gain == table_gain["gain"][0]) with pytest.raises(KeyError): interpolator_pointing(tel_id=2, time=t0 + 1 * u.s) with pytest.raises(KeyError): - interpolator_cal(tel_id=2, time=1.0) + interpolator_pedestal(tel_id=2, time=1.0) + with pytest.raises(KeyError): + interpolator_gain(tel_id=2, time=1.0) interpolator_pointing = PointingInterpolator(bounds_error=False) - interpolator_cal = CalibrationInterpolator(bounds_error=False) + interpolator_pedestal = PedestalInterpolator(bounds_error=False) + interpolator_gain = GainInterpolator(bounds_error=False) interpolator_pointing.add_table(1, table_pointing) - interpolator_cal.add_table(1, table_cal, "pedestal") + interpolator_pedestal.add_table(1, table_pedestal) + interpolator_gain.add_table(1, table_gain) for dt in (-0.1, 10.1) * u.s: alt, az = interpolator_pointing(tel_id=1, time=t0 + dt) assert np.isnan(alt.value) assert np.isnan(az.value) - assert all(np.isnan(interpolator_cal(tel_id=1, time=-0.1))) + assert all(np.isnan(interpolator_pedestal(tel_id=1, time=-0.1))) + assert all(np.isnan(interpolator_gain(tel_id=1, time=-0.1))) interpolator_pointing = PointingInterpolator(bounds_error=False, extrapolate=True) interpolator_pointing.add_table(1, table_pointing) From 86849d1d40430b9bd5df67d12523bb9f0a2fe9db Mon Sep 17 00:00:00 2001 From: Christoph Toennis Date: Mon, 5 Aug 2024 12:17:00 +0200 Subject: [PATCH 18/23] changing the name of GainInterpolator to FlatFieldInterpolator --- src/ctapipe/io/__init__.py | 4 ++-- src/ctapipe/io/interpolation.py | 22 +++++++++++----------- src/ctapipe/io/tests/test_interpolator.py | 22 +++++++++++----------- 3 files changed, 24 insertions(+), 24 deletions(-) diff --git a/src/ctapipe/io/__init__.py b/src/ctapipe/io/__init__.py index dbf555a7bd6..7974f1ffaaa 100644 --- a/src/ctapipe/io/__init__.py +++ b/src/ctapipe/io/__init__.py @@ -21,7 +21,7 @@ from .interpolation import ( Interpolator, PointingInterpolator, - GainInterpolator, + FlatFieldInterpolator, PedestalInterpolator, ) @@ -45,5 +45,5 @@ "Interpolator", "PointingInterpolator", "PedestalInterpolator", - "GainInterpolator", + "FlatFieldInterpolator", ] diff --git a/src/ctapipe/io/interpolation.py b/src/ctapipe/io/interpolation.py index dc0953e3892..e82d96b2e9b 100644 --- a/src/ctapipe/io/interpolation.py +++ b/src/ctapipe/io/interpolation.py @@ -165,7 +165,7 @@ class PointingInterpolator(Interpolator): """ telescope_data_group = "/dl0/monitoring/telescope/pointing" - required_columns = {"time", "azimuth", "altitude"} + required_columns = frozenset(["time", "azimuth", "altitude"]) expected_units = {"azimuth": u.rad, "altitude": u.rad} def __call__(self, tel_id, time): @@ -230,16 +230,16 @@ def add_table(self, tel_id, input_table): self._interpolators[tel_id]["alt"] = interp1d(mjd, alt, **self.interp_options) -class GainInterpolator(Interpolator): +class FlatFieldInterpolator(Interpolator): """ - Interpolator for relative gain data + Interpolator for flatfield data """ telescope_data_group = "dl1/calibration/gain" # TBD def __call__(self, tel_id, time): """ - Interpolate pedestal or gain for a given time and tel_id. + Interpolate flatfield data for a given time and tel_id. Parameters ---------- @@ -250,14 +250,14 @@ def __call__(self, tel_id, time): Returns ------- - gain : array [float] - interpolated relative gain + ffield : array [float] + interpolated flatfield data """ self._check_interpolators(tel_id) - gain = self._interpolators[tel_id]["gain"](time) - return gain + ffield = self._interpolators[tel_id]["gain"](time) + return ffield def add_table(self, tel_id, input_table): """ @@ -270,10 +270,10 @@ def add_table(self, tel_id, input_table): input_table : astropy.table.Table Table of pointing values, expected columns are ``time`` as ``Time`` column and "gain" - for the relative gain data + for the flatfield data """ - self.required_columns = {"time", "gain"} + self.required_columns = frozenset(["time", "gain"]) # TBD self._check_tables(input_table) @@ -329,7 +329,7 @@ def add_table(self, tel_id, input_table): for the pedestal data """ - self.required_columns = {"time", "pedestal"} + self.required_columns = frozenset(["time", "pedestal"]) # TBD self._check_tables(input_table) diff --git a/src/ctapipe/io/tests/test_interpolator.py b/src/ctapipe/io/tests/test_interpolator.py index 1279dc5c784..4d41ceeff0b 100644 --- a/src/ctapipe/io/tests/test_interpolator.py +++ b/src/ctapipe/io/tests/test_interpolator.py @@ -6,7 +6,7 @@ from astropy.time import Time from ctapipe.io.interpolation import ( - GainInterpolator, + FlatFieldInterpolator, PedestalInterpolator, PointingInterpolator, ) @@ -108,7 +108,7 @@ def test_bounds(): }, ) - table_gain = Table( + table_flatfield = Table( { "time": np.arange(0.0, 10.1, 2.0), "gain": np.reshape(np.random.normal(1.0, 1.0, 1850 * 6), (6, 1850)), @@ -117,10 +117,10 @@ def test_bounds(): interpolator_pointing = PointingInterpolator() interpolator_pedestal = PedestalInterpolator() - interpolator_gain = GainInterpolator() + interpolator_flatfield = FlatFieldInterpolator() interpolator_pointing.add_table(1, table_pointing) interpolator_pedestal.add_table(1, table_pedestal) - interpolator_gain.add_table(1, table_gain) + interpolator_flatfield.add_table(1, table_flatfield) error_message = "below the interpolation range" @@ -131,7 +131,7 @@ def test_bounds(): interpolator_pedestal(tel_id=1, time=-0.1) with pytest.raises(ValueError, match=error_message): - interpolator_gain(tel_id=1, time=-0.1) + interpolator_flatfield(tel_id=1, time=-0.1) with pytest.raises(ValueError, match="above the interpolation range"): interpolator_pointing(tel_id=1, time=t0 + 10.2 * u.s) @@ -142,21 +142,21 @@ def test_bounds(): pedestal = interpolator_pedestal(tel_id=1, time=1.0) assert all(pedestal == table_pedestal["pedestal"][0]) - gain = interpolator_gain(tel_id=1, time=1.0) - assert all(gain == table_gain["gain"][0]) + flatfield = interpolator_flatfield(tel_id=1, time=1.0) + assert all(flatfield == table_flatfield["gain"][0]) with pytest.raises(KeyError): interpolator_pointing(tel_id=2, time=t0 + 1 * u.s) with pytest.raises(KeyError): interpolator_pedestal(tel_id=2, time=1.0) with pytest.raises(KeyError): - interpolator_gain(tel_id=2, time=1.0) + interpolator_flatfield(tel_id=2, time=1.0) interpolator_pointing = PointingInterpolator(bounds_error=False) interpolator_pedestal = PedestalInterpolator(bounds_error=False) - interpolator_gain = GainInterpolator(bounds_error=False) + interpolator_flatfield = FlatFieldInterpolator(bounds_error=False) interpolator_pointing.add_table(1, table_pointing) interpolator_pedestal.add_table(1, table_pedestal) - interpolator_gain.add_table(1, table_gain) + interpolator_flatfield.add_table(1, table_flatfield) for dt in (-0.1, 10.1) * u.s: alt, az = interpolator_pointing(tel_id=1, time=t0 + dt) @@ -164,7 +164,7 @@ def test_bounds(): assert np.isnan(az.value) assert all(np.isnan(interpolator_pedestal(tel_id=1, time=-0.1))) - assert all(np.isnan(interpolator_gain(tel_id=1, time=-0.1))) + assert all(np.isnan(interpolator_flatfield(tel_id=1, time=-0.1))) interpolator_pointing = PointingInterpolator(bounds_error=False, extrapolate=True) interpolator_pointing.add_table(1, table_pointing) From 908f381ea344ccfe5f78b5b98d34ef31ff2a8efe Mon Sep 17 00:00:00 2001 From: Christoph Toennis Date: Mon, 5 Aug 2024 13:24:30 +0200 Subject: [PATCH 19/23] Fixing docstrings --- src/ctapipe/io/interpolation.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/ctapipe/io/interpolation.py b/src/ctapipe/io/interpolation.py index e82d96b2e9b..b68adbf6bf4 100644 --- a/src/ctapipe/io/interpolation.py +++ b/src/ctapipe/io/interpolation.py @@ -114,13 +114,13 @@ def add_table(self, tel_id, input_table): tel_id, the second needs to be the name of the parameter that is to be interpolated Parameters - -.--------- + ---------- tel_id : int Telescope id input_table : astropy.table.Table Table of pointing values, expected columns - are always ``time`` as ``Time`` column and other columns for the data - that is to be interpolated + are always ``time`` as ``Time`` column and + other columns for the data that is to be interpolated """ pass From 04f794a4ec60979d110b37c180adb70146f0f98e Mon Sep 17 00:00:00 2001 From: Christoph Toennis Date: Tue, 6 Aug 2024 10:37:01 +0200 Subject: [PATCH 20/23] Moved some things to avoid shadowing variables --- src/ctapipe/io/interpolation.py | 10 ++++++---- src/ctapipe/io/tests/test_interpolator.py | 6 ++++-- 2 files changed, 10 insertions(+), 6 deletions(-) diff --git a/src/ctapipe/io/interpolation.py b/src/ctapipe/io/interpolation.py index b68adbf6bf4..5870e71a91b 100644 --- a/src/ctapipe/io/interpolation.py +++ b/src/ctapipe/io/interpolation.py @@ -11,6 +11,8 @@ from .astropy_helpers import read_table +dimensionless_unit = u.Unit() + class StepFunction: @@ -236,6 +238,8 @@ class FlatFieldInterpolator(Interpolator): """ telescope_data_group = "dl1/calibration/gain" # TBD + required_columns = frozenset(["time", "gain"]) # TBD + expected_units = {"gain": dimensionless_unit} def __call__(self, tel_id, time): """ @@ -273,8 +277,6 @@ def add_table(self, tel_id, input_table): for the flatfield data """ - self.required_columns = frozenset(["time", "gain"]) # TBD - self._check_tables(input_table) input_table.sort("time") @@ -292,6 +294,8 @@ class PedestalInterpolator(Interpolator): """ telescope_data_group = "dl1/calibration/pedestal" # TBD + required_columns = frozenset(["time", "pedestal"]) # TBD + expected_units = {"pedestal": dimensionless_unit} def __call__(self, tel_id, time): """ @@ -329,8 +333,6 @@ def add_table(self, tel_id, input_table): for the pedestal data """ - self.required_columns = frozenset(["time", "pedestal"]) # TBD - self._check_tables(input_table) input_table.sort("time") diff --git a/src/ctapipe/io/tests/test_interpolator.py b/src/ctapipe/io/tests/test_interpolator.py index 4d41ceeff0b..930e1e7d73c 100644 --- a/src/ctapipe/io/tests/test_interpolator.py +++ b/src/ctapipe/io/tests/test_interpolator.py @@ -104,14 +104,16 @@ def test_bounds(): table_pedestal = Table( { "time": np.arange(0.0, 10.1, 2.0), - "pedestal": np.reshape(np.random.normal(4.0, 1.0, 1850 * 6), (6, 1850)), + "pedestal": np.reshape(np.random.normal(4.0, 1.0, 1850 * 6), (6, 1850)) + * u.Unit(), }, ) table_flatfield = Table( { "time": np.arange(0.0, 10.1, 2.0), - "gain": np.reshape(np.random.normal(1.0, 1.0, 1850 * 6), (6, 1850)), + "gain": np.reshape(np.random.normal(1.0, 1.0, 1850 * 6), (6, 1850)) + * u.Unit(), }, ) From e9a21123b9abd597a301aba9e24bd08fe668a50a Mon Sep 17 00:00:00 2001 From: Christoph Toennis Date: Tue, 6 Aug 2024 11:24:16 +0200 Subject: [PATCH 21/23] interpolator now copies input tables so they are not being altered when interpolating --- src/ctapipe/io/interpolation.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/ctapipe/io/interpolation.py b/src/ctapipe/io/interpolation.py index 5870e71a91b..2f82a708710 100644 --- a/src/ctapipe/io/interpolation.py +++ b/src/ctapipe/io/interpolation.py @@ -215,6 +215,7 @@ def add_table(self, tel_id, input_table): if not isinstance(input_table["time"], Time): raise TypeError("'time' column of pointing table must be astropy.time.Time") + input_table = input_table.copy() input_table.sort("time") az = input_table["azimuth"].quantity.to_value(u.rad) @@ -279,6 +280,7 @@ def add_table(self, tel_id, input_table): self._check_tables(input_table) + input_table = input_table.copy() input_table.sort("time") time = input_table["time"] gain = input_table["gain"] @@ -335,6 +337,7 @@ def add_table(self, tel_id, input_table): self._check_tables(input_table) + input_table = input_table.copy() input_table.sort("time") time = input_table["time"] pedestal = input_table["pedestal"] From a5569c644e535e0b9344efb0f004b298af02b20e Mon Sep 17 00:00:00 2001 From: Christoph Toennis Date: Thu, 8 Aug 2024 10:52:17 +0200 Subject: [PATCH 22/23] I removed the Flatfield and pedestal interpolators --- src/ctapipe/io/__init__.py | 4 - src/ctapipe/io/interpolation.py | 116 ---------------------- src/ctapipe/io/tests/test_interpolator.py | 45 --------- 3 files changed, 165 deletions(-) diff --git a/src/ctapipe/io/__init__.py b/src/ctapipe/io/__init__.py index 7974f1ffaaa..28603643160 100644 --- a/src/ctapipe/io/__init__.py +++ b/src/ctapipe/io/__init__.py @@ -21,8 +21,6 @@ from .interpolation import ( Interpolator, PointingInterpolator, - FlatFieldInterpolator, - PedestalInterpolator, ) __all__ = [ @@ -44,6 +42,4 @@ "get_hdf5_datalevels", "Interpolator", "PointingInterpolator", - "PedestalInterpolator", - "FlatFieldInterpolator", ] diff --git a/src/ctapipe/io/interpolation.py b/src/ctapipe/io/interpolation.py index 2f82a708710..2ebbbbfd2fe 100644 --- a/src/ctapipe/io/interpolation.py +++ b/src/ctapipe/io/interpolation.py @@ -11,8 +11,6 @@ from .astropy_helpers import read_table -dimensionless_unit = u.Unit() - class StepFunction: @@ -231,117 +229,3 @@ def add_table(self, tel_id, input_table): self._interpolators[tel_id] = {} self._interpolators[tel_id]["az"] = interp1d(mjd, az, **self.interp_options) self._interpolators[tel_id]["alt"] = interp1d(mjd, alt, **self.interp_options) - - -class FlatFieldInterpolator(Interpolator): - """ - Interpolator for flatfield data - """ - - telescope_data_group = "dl1/calibration/gain" # TBD - required_columns = frozenset(["time", "gain"]) # TBD - expected_units = {"gain": dimensionless_unit} - - def __call__(self, tel_id, time): - """ - Interpolate flatfield data for a given time and tel_id. - - Parameters - ---------- - tel_id : int - telescope id - time : astropy.time.Time - time for which to interpolate the calibration data - - Returns - ------- - ffield : array [float] - interpolated flatfield data - """ - - self._check_interpolators(tel_id) - - ffield = self._interpolators[tel_id]["gain"](time) - return ffield - - def add_table(self, tel_id, input_table): - """ - Add a table to this interpolator - - Parameters - ---------- - tel_id : int - Telescope id - input_table : astropy.table.Table - Table of pointing values, expected columns - are ``time`` as ``Time`` column and "gain" - for the flatfield data - """ - - self._check_tables(input_table) - - input_table = input_table.copy() - input_table.sort("time") - time = input_table["time"] - gain = input_table["gain"] - self._interpolators[tel_id] = {} - self._interpolators[tel_id]["gain"] = StepFunction( - time, gain, **self.interp_options - ) - - -class PedestalInterpolator(Interpolator): - """ - Interpolator for Pedestal data - """ - - telescope_data_group = "dl1/calibration/pedestal" # TBD - required_columns = frozenset(["time", "pedestal"]) # TBD - expected_units = {"pedestal": dimensionless_unit} - - def __call__(self, tel_id, time): - """ - Interpolate pedestal or gain for a given time and tel_id. - - Parameters - ---------- - tel_id : int - telescope id - time : astropy.time.Time - time for which to interpolate the calibration data - - Returns - ------- - pedestal : array [float] - interpolated pedestal values - """ - - self._check_interpolators(tel_id) - - pedestal = self._interpolators[tel_id]["pedestal"](time) - return pedestal - - def add_table(self, tel_id, input_table): - """ - Add a table to this interpolator - - Parameters - ---------- - tel_id : int - Telescope id - input_table : astropy.table.Table - Table of pointing values, expected columns - are ``time`` as ``Time`` column and "pedestal" - for the pedestal data - """ - - self._check_tables(input_table) - - input_table = input_table.copy() - input_table.sort("time") - time = input_table["time"] - pedestal = input_table["pedestal"] - self._interpolators[tel_id] = {} - self._interpolators[tel_id]["pedestal"] = StepFunction( - time, pedestal, **self.interp_options - ) diff --git a/src/ctapipe/io/tests/test_interpolator.py b/src/ctapipe/io/tests/test_interpolator.py index 930e1e7d73c..02f4c4ce306 100644 --- a/src/ctapipe/io/tests/test_interpolator.py +++ b/src/ctapipe/io/tests/test_interpolator.py @@ -6,8 +6,6 @@ from astropy.time import Time from ctapipe.io.interpolation import ( - FlatFieldInterpolator, - PedestalInterpolator, PointingInterpolator, ) @@ -101,40 +99,13 @@ def test_bounds(): }, ) - table_pedestal = Table( - { - "time": np.arange(0.0, 10.1, 2.0), - "pedestal": np.reshape(np.random.normal(4.0, 1.0, 1850 * 6), (6, 1850)) - * u.Unit(), - }, - ) - - table_flatfield = Table( - { - "time": np.arange(0.0, 10.1, 2.0), - "gain": np.reshape(np.random.normal(1.0, 1.0, 1850 * 6), (6, 1850)) - * u.Unit(), - }, - ) - interpolator_pointing = PointingInterpolator() - interpolator_pedestal = PedestalInterpolator() - interpolator_flatfield = FlatFieldInterpolator() interpolator_pointing.add_table(1, table_pointing) - interpolator_pedestal.add_table(1, table_pedestal) - interpolator_flatfield.add_table(1, table_flatfield) - error_message = "below the interpolation range" with pytest.raises(ValueError, match=error_message): interpolator_pointing(tel_id=1, time=t0 - 0.1 * u.s) - with pytest.raises(ValueError, match=error_message): - interpolator_pedestal(tel_id=1, time=-0.1) - - with pytest.raises(ValueError, match=error_message): - interpolator_flatfield(tel_id=1, time=-0.1) - with pytest.raises(ValueError, match="above the interpolation range"): interpolator_pointing(tel_id=1, time=t0 + 10.2 * u.s) @@ -142,32 +113,16 @@ def test_bounds(): assert u.isclose(alt, 69 * u.deg) assert u.isclose(az, 1 * u.deg) - pedestal = interpolator_pedestal(tel_id=1, time=1.0) - assert all(pedestal == table_pedestal["pedestal"][0]) - flatfield = interpolator_flatfield(tel_id=1, time=1.0) - assert all(flatfield == table_flatfield["gain"][0]) with pytest.raises(KeyError): interpolator_pointing(tel_id=2, time=t0 + 1 * u.s) - with pytest.raises(KeyError): - interpolator_pedestal(tel_id=2, time=1.0) - with pytest.raises(KeyError): - interpolator_flatfield(tel_id=2, time=1.0) interpolator_pointing = PointingInterpolator(bounds_error=False) - interpolator_pedestal = PedestalInterpolator(bounds_error=False) - interpolator_flatfield = FlatFieldInterpolator(bounds_error=False) interpolator_pointing.add_table(1, table_pointing) - interpolator_pedestal.add_table(1, table_pedestal) - interpolator_flatfield.add_table(1, table_flatfield) - for dt in (-0.1, 10.1) * u.s: alt, az = interpolator_pointing(tel_id=1, time=t0 + dt) assert np.isnan(alt.value) assert np.isnan(az.value) - assert all(np.isnan(interpolator_pedestal(tel_id=1, time=-0.1))) - assert all(np.isnan(interpolator_flatfield(tel_id=1, time=-0.1))) - interpolator_pointing = PointingInterpolator(bounds_error=False, extrapolate=True) interpolator_pointing.add_table(1, table_pointing) alt, az = interpolator_pointing(tel_id=1, time=t0 - 1 * u.s) From 77153101d3048830f7587baaef7e69a97072af79 Mon Sep 17 00:00:00 2001 From: Christoph Toennis Date: Mon, 12 Aug 2024 14:05:25 +0200 Subject: [PATCH 23/23] removing refernce to StepInterpolator --- docs/changes/2600.feature.rst | 1 - src/ctapipe/io/interpolation.py | 49 --------------------------------- 2 files changed, 50 deletions(-) diff --git a/docs/changes/2600.feature.rst b/docs/changes/2600.feature.rst index aead92bc84b..1c996397c6d 100644 --- a/docs/changes/2600.feature.rst +++ b/docs/changes/2600.feature.rst @@ -1,2 +1 @@ Add Interpolator class to generalize the PointingInterpolator in the io collection. -Add new interpolator class called StepInterpolator to be used for calibration data. diff --git a/src/ctapipe/io/interpolation.py b/src/ctapipe/io/interpolation.py index 2ebbbbfd2fe..82de9f0bd19 100644 --- a/src/ctapipe/io/interpolation.py +++ b/src/ctapipe/io/interpolation.py @@ -12,55 +12,6 @@ from .astropy_helpers import read_table -class StepFunction: - - """ - Step function Interpolator for the gain and pedestals - Interpolates data so that for each time the value from the closest previous - point given. - - Parameters - ---------- - values : None | np.array - Numpy array of the data that is to be interpolated. - The first dimension needs to be an index over time - times : None | np.array - Time values over which data are to be interpolated - need to be sorted and have same length as first dimension of values - """ - - def __init__( - self, - times, - values, - bounds_error=True, - fill_value="extrapolate", - assume_sorted=True, - copy=False, - ): - self.values = values - self.times = times - self.bounds_error = bounds_error - self.fill_value = fill_value - - def __call__(self, point): - if point < self.times[0]: - if self.bounds_error: - raise ValueError("below the interpolation range") - - if self.fill_value == "extrapolate": - return self.values[0] - - else: - a = np.empty(self.values[0].shape) - a[:] = np.nan - return a - - else: - i = np.searchsorted(self.times, point, side="left") - return self.values[i - 1] - - class Interpolator(Component, metaclass=ABCMeta): """ Interpolator parent class.