diff --git a/docs/src/index.rst b/docs/src/index.rst index 8b9af28..d021764 100644 --- a/docs/src/index.rst +++ b/docs/src/index.rst @@ -29,6 +29,7 @@ AARE python/cluster/index python/file/index python/histogram/index + python/pedestal/index pyFit diff --git a/docs/src/python/pedestal/index.rst b/docs/src/python/pedestal/index.rst new file mode 100644 index 0000000..4e0a05e --- /dev/null +++ b/docs/src/python/pedestal/index.rst @@ -0,0 +1,9 @@ +Pedestal +======== + +.. toctree:: + :caption: Pedestal + :maxdepth: 1 + + pyFastPedestal + pyPedestal diff --git a/docs/src/python/pedestal/pyFastPedestal.rst b/docs/src/python/pedestal/pyFastPedestal.rst new file mode 100644 index 0000000..00ada7f --- /dev/null +++ b/docs/src/python/pedestal/pyFastPedestal.rst @@ -0,0 +1,83 @@ +FastPedestal +============ + +``FastPedestal`` calculates a running mean, variance and standard deviation for each pixel in a +series of frames. The python binding only exposes ``uint16`` input but the underlying +C++ class is templated. Initialize it with ``n_samples`` frames using +``push_init()``. Once ``ready`` is true, use ``push()`` for steady-state +updates. + +.. warning:: + + FastPedestal is not usable until you have pushed ``n_samples`` initial frames with ``push_init(raw)``. + You can check the state with ``ready``. + +The public factory selects the bound C++ specialization from ``dtype``: + +* ``numpy.float64`` creates ``FastPedestal_d`` +* ``numpy.float32`` creates ``FastPedestal_f`` +* ``numpy.int16`` creates ``FastPedestal_i16`` + +The internal calculations are done with double, but the cached mean and on demand var and std are returned in the specified type. + +Factory +------- + +.. py:currentmodule:: aare + +.. autofunction:: FastPedestal + +Loading from a file +------------------- + +``FastPedestal.from_file()`` initializes the pedestal from ``n_samples`` +frames after ``skip_first``, then applies steady-state updates for any frames +remaining in the file. The input frames must contain ``uint16`` data; ``dtype`` +selects the output type of the pedestal statistics. + +.. autofunction:: aare.FastPedestal.from_file + +.. code-block:: python + + pedestal = FastPedestal.from_file( + "frames.npy", n_samples=100, skip_first=10, dtype=np.float32 + ) + +Example +------- + +.. code-block:: python + + import numpy as np + from aare import FastPedestal + + pedestal = FastPedestal(512, 1024, n_samples=100, dtype=np.float32) + + # Initialize with n_samples frames + for frame in initialization_frames: + pedestal.push_init(frame) + + # Now we can push a frame for pedestal update + if pedestal.ready: + pedestal.push(next_frame) + + # Mean and std are also ready + mean = pedestal.mean() + noise = pedestal.std() + + # Direct pedestal subtraction is also supported + for frame in raw_data: + image = frame - pedestal + +Complete API +------------ + +The API below is for the ``float64`` specialization. All dtype variants share +the same API. + +.. autoclass:: aare._aare.FastPedestal_d + :special-members: __init__ + :members: + :undoc-members: + :show-inheritance: + :inherited-members: diff --git a/docs/src/python/pedestal/pyPedestal.rst b/docs/src/python/pedestal/pyPedestal.rst new file mode 100644 index 0000000..7d6aecb --- /dev/null +++ b/docs/src/python/pedestal/pyPedestal.rst @@ -0,0 +1,42 @@ +Pedestal +======== + +``Pedestal`` calculates a running mean and variance for each pixel in a series +of ``uint16`` frames. ``push()`` updates the cached mean immediately. For +faster batch initialization, use ``push_no_update()`` for each frame and call +``update_mean()`` after the batch. + +Three specializations are available from :mod:`aare`: + +* ``Pedestal_d`` uses ``float64`` storage +* ``Pedestal_f`` uses ``float32`` storage +* ``Pedestal_i16`` uses ``int16`` storage + +Example +------- + +.. code-block:: python + + from aare import Pedestal_d + + pedestal = Pedestal_d(512, 1024, 100) + + for frame in initialization_frames: + pedestal.push_no_update(frame) + + pedestal.update_mean() + mean = pedestal.mean() + noise = pedestal.std() + +Complete API +------------ + +The API below is for the ``float64`` specialization. All dtype variants share +the same API. + +.. autoclass:: aare._aare.Pedestal_d + :special-members: __init__ + :members: + :undoc-members: + :show-inheritance: + :inherited-members: diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt index 0377548..20932b4 100644 --- a/python/CMakeLists.txt +++ b/python/CMakeLists.txt @@ -39,7 +39,9 @@ set(PYTHON_FILES aare/ClusterFinder.py aare/ClusterVector.py aare/Cluster.py + aare/FastPedestal.py aare/calibration.py + aare/factory.py aare/func.py aare/RawFile.py aare/transform.py diff --git a/python/aare/Cluster.py b/python/aare/Cluster.py index f2ea07e..7946567 100644 --- a/python/aare/Cluster.py +++ b/python/aare/Cluster.py @@ -1,6 +1,6 @@ from . import _aare import numpy as np -from .ClusterFinder import _type_to_char +from .factory import _type_to_char def Cluster(x : int, y : int, data, cluster_size=(3,3), dtype = np.int32): @@ -21,4 +21,4 @@ def Cluster(x : int, y : int, data, cluster_size=(3,3), dtype = np.int32): except AttributeError: raise ValueError(f"Unsupported combination of type and cluster size: {dtype}/{cluster_size} when requesting {class_name}") - return cls(x, y, data) \ No newline at end of file + return cls(x, y, data) diff --git a/python/aare/ClusterFinder.py b/python/aare/ClusterFinder.py index eae7ad7..93af877 100644 --- a/python/aare/ClusterFinder.py +++ b/python/aare/ClusterFinder.py @@ -1,21 +1,10 @@ # SPDX-License-Identifier: MPL-2.0 from . import _aare import numpy as np +from .factory import _type_to_char _supported_cluster_sizes = [(2,2), (3,3), (5,5), (7,7), (9,9),] -def _type_to_char(dtype): - if dtype == np.int32: - return 'i' - elif dtype == np.float32: - return 'f' - elif dtype == np.float64: - return 'd' - elif dtype == np.int16: - return 'i16' - else: - raise ValueError(f"Unsupported dtype: {dtype}. Only np.int32, np.float32, and np.float64 are supported.") - def _get_class(name, cluster_size, dtype): """ Helper function to get the class based on the name, cluster size, and dtype. diff --git a/python/aare/FastPedestal.py b/python/aare/FastPedestal.py new file mode 100644 index 0000000..710e343 --- /dev/null +++ b/python/aare/FastPedestal.py @@ -0,0 +1,34 @@ +# SPDX-License-Identifier: MPL-2.0 + +import numpy as np + +from .factory import _get_typed_class + + +def _get_fast_pedestal_class(dtype): + return _get_typed_class("FastPedestal", dtype) + + +def FastPedestal(rows, cols, n_samples=1000, dtype=np.float64): + """Create a FastPedestal with the requested output dtype. + + This factory hides the dtype suffix used by the templated C++ bindings. + Supported dtypes are ``np.float64``, ``np.float32``, and ``np.int16``. + """ + cls = _get_fast_pedestal_class(dtype) + return cls(rows, cols, n_samples) + + +def from_file(filename, n_samples=1000, skip_first=0, dtype=np.float64): + """Create a FastPedestal from frames in a file. + + The input frames must contain uint16 data. ``dtype`` selects the output + type of the pedestal mean, variance, and standard deviation. + """ + cls = _get_fast_pedestal_class(dtype) + return cls.from_file( + filename, n_samples=n_samples, skip_first=skip_first + ) + + +FastPedestal.from_file = from_file diff --git a/python/aare/__init__.py b/python/aare/__init__.py index 92c3652..005af54 100644 --- a/python/aare/__init__.py +++ b/python/aare/__init__.py @@ -23,6 +23,7 @@ from ._aare import corner # from ._aare import ClusterFinderMT, ClusterCollector, ClusterFileSink, ClusterVector_i from ._version import __version__ +from .FastPedestal import FastPedestal from .ClusterFinder import ClusterFinder, ClusterCollector, ClusterFinderMT, ClusterFileSink, ClusterFile from .ClusterVector import ClusterVector from .Cluster import Cluster diff --git a/python/aare/factory.py b/python/aare/factory.py new file mode 100644 index 0000000..4da3fdb --- /dev/null +++ b/python/aare/factory.py @@ -0,0 +1,36 @@ +# SPDX-License-Identifier: MPL-2.0 + +import numpy as np + +from . import _aare + + +_TYPE_TO_CHAR = { + np.dtype(np.int32): "i", + np.dtype(np.float32): "f", + np.dtype(np.float64): "d", + np.dtype(np.int16): "i16", +} + + +def _type_to_char(dtype): + """Return the suffix used by bindings instantiated for ``dtype``.""" + try: + return _TYPE_TO_CHAR[np.dtype(dtype)] + except (KeyError, TypeError): + supported = ", ".join(str(dtype) for dtype in _TYPE_TO_CHAR) + raise ValueError( + f"Unsupported dtype: {dtype}. Supported dtypes are {supported}." + ) from None + + +def _get_typed_class(name, dtype): + """Return a bound class named ``_``.""" + class_name = f"{name}_{_type_to_char(dtype)}" + try: + return getattr(_aare, class_name) + except AttributeError: + raise ValueError( + f"Unsupported dtype for {name}: {dtype} " + f"(binding {class_name} is not available)." + ) from None diff --git a/python/src/bind_FastPedestal.hpp b/python/src/bind_FastPedestal.hpp index a701d8d..0ca238e 100644 --- a/python/src/bind_FastPedestal.hpp +++ b/python/src/bind_FastPedestal.hpp @@ -14,8 +14,9 @@ template void define_fast_pedestal_bindings(py::module &m, const std::string &name) { py::class_>(m, name.c_str(), py::buffer_protocol()) - .def(py::init()) - .def(py::init()) + .def(py::init(), py::arg("rows"), + py::arg("cols"), py::arg("n_samples")) + .def(py::init(), py::arg("rows"), py::arg("cols")) .def( "mean", @@ -24,7 +25,7 @@ void define_fast_pedestal_bindings(py::module &m, const std::string &name) { *mean = self.mean(); return return_image_data(mean); }, - "Return a copy of the mean of the pedestal as a NumPy array") + "Return a copy of the cached mean as a NumPy array") .def( "var", [](FastPedestal &self) { @@ -32,7 +33,7 @@ void define_fast_pedestal_bindings(py::module &m, const std::string &name) { *variance = self.variance(); return return_image_data(variance); }, - "Return a copy of the variance of the pedestal as a NumPy array") + "Calculate the variance and return it as a NumPy array") .def( "std", [](FastPedestal &self) { @@ -40,14 +41,15 @@ void define_fast_pedestal_bindings(py::module &m, const std::string &name) { *standard_deviation = self.std(); return return_image_data(standard_deviation); }, - "Return a copy of the standard deviation of the pedestal as a " + "Calculate the standard deviation and return it as a " "NumPy array") .def( "view", [](py::object self_py) { return py::module_::import("numpy").attr("asarray")(self_py); }, - "Return non-owning, non-writable view of the pedestal as a NumPy " + "Return non-owning, non-writable view of the pedestal mean as a " + "NumPy " "array") // We need to buffer protocol to allow for numpy operations using the @@ -95,7 +97,8 @@ void define_fast_pedestal_bindings(py::module &m, const std::string &name) { "pushed at least n_samples samples)") .def_property_readonly( "n_samples", &FastPedestal::n_samples, - "Return the number of samples to push to the pedestal to be ready") + "Return the number of samples to push to the pedestal to be ready. " + "This value also affect how fast it updates in the steady state.") .def("clone", [](FastPedestal &pedestal) { return FastPedestal(pedestal); diff --git a/python/tests/test_FastPedestal.py b/python/tests/test_FastPedestal.py index b5d774f..b3e40dc 100644 --- a/python/tests/test_FastPedestal.py +++ b/python/tests/test_FastPedestal.py @@ -1,7 +1,88 @@ import numpy as np import pytest -from aare import FastPedestal_d, FastPedestal_f +from aare import ( + FastPedestal, + FastPedestal_d, + FastPedestal_f, + FastPedestal_i16, +) + + +@pytest.mark.parametrize( + ("dtype", "pedestal_type"), + [ + (np.float64, FastPedestal_d), + (np.float32, FastPedestal_f), + (np.int16, FastPedestal_i16), + ], +) +def test_fast_pedestal_factory(dtype, pedestal_type): + pedestal = FastPedestal(2, 3, n_samples=4, dtype=dtype) + + assert isinstance(pedestal, pedestal_type) + assert pedestal.rows == 2 + assert pedestal.cols == 3 + assert pedestal.n_samples == 4 + + +def test_fast_pedestal_factory_defaults_to_double(): + assert isinstance(FastPedestal(2, 3), FastPedestal_d) + + +def test_fast_pedestal_factory_rejects_unbound_dtype(): + with pytest.raises(ValueError, match="Unsupported dtype for FastPedestal"): + FastPedestal(2, 3, dtype=np.int32) + + +@pytest.mark.parametrize( + ("kwargs", "expected_n_samples"), + [ + ({"rows": 2, "cols": 3}, 1000), + ({"rows": 2, "cols": 3, "n_samples": 4}, 4), + ], +) +def test_fast_pedestal_binding_accepts_constructor_keywords( + kwargs, expected_n_samples +): + pedestal = FastPedestal_d(**kwargs) + + assert pedestal.rows == 2 + assert pedestal.cols == 3 + assert pedestal.n_samples == expected_n_samples + + +@pytest.mark.parametrize( + ("dtype", "pedestal_type", "expected_dtype"), + [ + (np.float64, FastPedestal_d, np.float64), + (np.float32, FastPedestal_f, np.float32), + (np.int16, FastPedestal_i16, np.int16), + ], +) +def test_fast_pedestal_factory_from_file( + tmp_path, dtype, pedestal_type, expected_dtype +): + frames = np.array( + [[[100, 100]], [[2, 4]], [[4, 6]], [[5, 7]]], dtype=np.uint16 + ) + filename = tmp_path / "frames.npy" + np.save(filename, frames) + + pedestal = FastPedestal.from_file( + filename, n_samples=2, skip_first=1, dtype=dtype + ) + + assert isinstance(pedestal, pedestal_type) + assert pedestal.ready + assert pedestal.cur_samples == 2 + assert pedestal.mean().dtype == expected_dtype + np.testing.assert_array_equal(pedestal.mean(), [[4, 6]]) + + +def test_fast_pedestal_factory_from_file_rejects_unbound_dtype(): + with pytest.raises(ValueError, match="Unsupported dtype for FastPedestal"): + FastPedestal.from_file("unused.npy", dtype=np.int32) @pytest.mark.parametrize(