detectors/aare
25
0
Fork 0
mirror of https://github.com/slsdetectorgroup/aare.git synced 2026-09-07 10:42:37 +02:00
Code Issues Actions Packages Releases Activity

Merge branch 'main' into feature/cuda_clusterfinder
Build on RHEL8 / build (push) Successful in 3m34s
Details
Build on RHEL9 / build (push) Successful in 3m42s
Details
Run tests using data on local RHEL8 / build (push) Failing after 4m46s
Details

Browse Source
This commit is contained in:
Khalil Ferjaoui
2026-09-02 14:18:37 +02:00
committed by GitHub
parent 8567d3f00e eddb919328
commit 1d963a778b
100 changed files with 4748 additions and 2464 deletions
Download Patch File Download Diff File Expand all files Collapse all files
.gitignore
+3 -1
View File
@@ -127,6 +127,8 @@ cover/
*.mo
*.pot
.*
# Django stuff:
*.log
local_settings.py
@@ -192,4 +194,4 @@ cython_debug/
.ruff_cache/
# user defined
wheelhouse/
wheelhouse/
AGENTS.md
+166
View File
@@ -0,0 +1,166 @@
# Repository guidelines
## Scope and project goals
These instructions apply to the entire repository unless a more specific
`AGENTS.md` exists in a subdirectory.
Aare is a data-analysis library for PSI hybrid detectors. The C++17 core is
the canonical implementation, while Python is the main user-facing interface
through pybind11. Changes should preserve the project's priorities: fast and
efficient processing, simple interfaces, API stability, and a small dependency
footprint.
## Repository layout
- `include/aare/` contains the public C++ API and header-defined templates.
- `src/` contains C++ implementations and colocated Catch2 tests named
`*.test.cpp`.
- `tests/` contains the C++ test executable, configuration, and shared helpers.
- `python/src/` contains pybind11 bindings and module registration.
- `python/aare/` contains Python facades, convenience APIs, and public exports.
- `python/tests/` contains the pytest suite.
- `docs/src/` contains Sphinx/reStructuredText documentation. Doxygen and
Breathe integrate the C++ API into the Sphinx site.
- `benchmarks/` contains Google Benchmark programs.
- `cmake/`, `conda-recipe/`, and `pyproject.toml` support builds and packaging.
- `RELEASE.md` contains pending and published release notes.
## Build environment
The project requires CMake 3.15 or newer, C++17 with compiler extensions
disabled, and Python 3.11 or newer. `etc/dev-env.yml` defines the Conda
development environment.
CMake fetches several dependencies by default. Use
`-DAARE_SYSTEM_LIBRARIES=ON` only when all required system or Conda packages
are available.
Use this configuration for normal development:
```bash
cmake -S . -B build \
-DCMAKE_BUILD_TYPE=Debug \
-DAARE_TESTS=ON \
-DAARE_PYTHON_BINDINGS=ON
cmake --build build -j4
```
Useful optional settings include `AARE_DOCS`, `AARE_BENCHMARKS`, `AARE_ASAN`,
and `AARE_WARNINGS_AS_ERRORS`. Reconfigure an existing build directory instead
of creating alternate in-tree build layouts unless isolation is needed.
## Implementation conventions
- Put public declarations in `include/aare/` and private implementation details
in `src/`.
- Add new compiled headers, sources, and C++ tests to the explicit lists in the
root `CMakeLists.txt`.
- Keep template implementations in headers unless the supported types are
explicitly instantiated.
- Follow the naming in adjacent code. Broadly, use the `aare` namespace,
CamelCase types, and snake_case functions.
- Treat ownership, lifetime, const-correctness, array shapes, and buffer
contiguity as part of the API when working with `NDArray`, `NDView`, or NumPy
bindings.
- Avoid unnecessary allocations and copies in detector-data and per-pixel
processing paths. Add or update a benchmark when performance is central to a
change.
- Place APIs that may change without notice under the existing experimental
namespace/module.
- Start new source files with `SPDX-License-Identifier: MPL-2.0`, using the
appropriate comment syntax.
- Prefer descriptive names to comments that only restate the code.
- Do not add a dependency unless the benefit justifies the packaging and
deployment cost.
## Python-facing changes
A Python-facing feature can require coordinated changes in several layers:
1. Update the C++ public API and implementation.
2. Add or update its binding in `python/src/`.
3. Register new bindings in `python/src/module.cpp`.
4. Update the facade or public exports in `python/aare/`.
5. If adding a Python module, add it to `PYTHON_FILES` in
`python/CMakeLists.txt` so it is copied and installed.
6. Add Python tests and update user documentation.
Bindings must validate NumPy dimensions and data types before constructing
views. Do not return a view whose backing C++ or Python storage can expire while
the view remains reachable.
## Tests
Run tests that do not require external detector data with:
```bash
ctest --test-dir build --output-on-failure -j4
PYTHONPATH="$PWD/build" python -m pytest python/tests
```
For focused runs, use a Catch2 tag or an individual pytest file/test:
```bash
build/run_tests "[tag]"
PYTHONPATH="$PWD/build" python -m pytest python/tests/test_example.py
```
Large detector test files live outside this repository. To include data-backed
tests, set `AARE_TEST_DATA` and opt in explicitly:
```bash
export AARE_TEST_DATA=/path/to/aare-test-data
build/run_tests "[.with-data]"
PYTHONPATH="$PWD/build" python -m pytest python/tests --with-data
```
- Start bug fixes with a failing regression test when practical.
- Put C++ tests beside the relevant implementation as `src/Thing.test.cpp`.
- Mark C++ tests requiring external files with `[.with-data]`.
- Mark Python tests requiring external files with `@pytest.mark.withdata`.
- Run both suites for changes that cross the C++/Python boundary.
- Do not silently skip required data-backed coverage. Report when the external
test data is unavailable.
## Formatting and static analysis
Before handing off a broad change, run the relevant checks:
```bash
pre-commit run --all-files
cmake --build build --target check-format
cmake --build build --target clang-tidy
```
C++ formatting follows `.clang-format` (four-space indentation and an 80-column
limit). CMake files are checked by `cmake-format`. Avoid formatting unrelated
code as part of a focused change.
## Documentation and release notes
- Update the relevant `.rst` pages for public behavior or API changes.
- Add new pages to `docs/src/index.rst` or the relevant nested toctree.
- Update `RELEASE.md` under `## Next` for user-visible features, bug fixes, and
API changes.
- Do not change `VERSION` unless performing an explicitly requested release.
- Preserve compatibility with existing detector formats and older recorded
files where practical. Call out intentional API or format incompatibilities.
Build the documentation with:
```bash
cmake -S . -B build \
-DAARE_DOCS=ON \
-DAARE_PYTHON_BINDINGS=ON
cmake --build build --target docs
```
## Working practices and handoff
- Inspect the adjacent implementation, tests, and documentation before editing.
- Keep changes focused and preserve unrelated modifications in the worktree.
- Do not edit generated files or fetched dependency sources under `build/`.
- For large features, prefer independently testable increments.
- At handoff, summarize the behavior and important files changed, checks run,
checks not run and why, and any compatibility or performance considerations.
CMakeLists.txt
+39 -78
View File
@@ -55,6 +55,10 @@ option(
"Install the python extension in the install tree under CMAKE_INSTALL_PREFIX/aare/"
OFF)
option(AARE_ASAN "Enable AddressSanitizer" OFF)
option(
AARE_TUNE_LOCAL
"Optimize for the building machine's CPU (-march=native -mtune=native). Not portable, the resulting binaries might not run on other machines."
OFF)
option(AARE_CUDA "Build CUDA cluster finder backend" OFF)
set(AARE_CUDA_ARCHITECTURES
@@ -67,7 +71,6 @@ option(AARE_FETCH_PYBIND11 "Use FetchContent to download pybind11" ON)
option(AARE_FETCH_CATCH "Use FetchContent to download catch2" ON)
option(AARE_FETCH_JSON "Use FetchContent to download nlohmann::json" ON)
option(AARE_FETCH_ZMQ "Use FetchContent to download libzmq" ON)
option(AARE_FETCH_LMFIT "Use FetchContent to download lmfit" ON)
option(AARE_FETCH_MINUIT2 "Use FetchContent to download Minuit2" ON)
option(AARE_WARNINGS_AS_ERRORS "Treat warnings as errors during compilation"
@@ -92,8 +95,8 @@ if(AARE_SYSTEM_LIBRARIES)
set(AARE_FETCH_ZMQ
OFF
CACHE BOOL "Disabled FetchContent for libzmq" FORCE)
# Still fetch lmfit and Minuit2 when setting AARE_SYSTEM_LIBRARIES since these
# are not available on conda-forge
# Still fetch Minuit2 when setting AARE_SYSTEM_LIBRARIES since it is not
# available on conda-forge
endif()
if(AARE_CUDA)
@@ -122,61 +125,6 @@ endif()
set(CMAKE_EXPORT_COMPILE_COMMANDS ON)
if(AARE_FETCH_LMFIT)
# TODO! Should we fetch lmfit from the web or inlcude a tar.gz in the repo?
set(LMFIT_PATCH_COMMAND git apply
${CMAKE_CURRENT_SOURCE_DIR}/patches/lmfit.patch)
# For cmake < 3.28 we can't supply EXCLUDE_FROM_ALL to FetchContent_Declare so
# we need this workaround
if(${CMAKE_VERSION} VERSION_LESS "3.28")
FetchContent_Declare(
lmfit
GIT_REPOSITORY https://jugit.fz-juelich.de/mlz/lmfit.git
GIT_TAG main
PATCH_COMMAND ${LMFIT_PATCH_COMMAND}
UPDATE_DISCONNECTED 1)
else()
FetchContent_Declare(
lmfit
GIT_REPOSITORY https://jugit.fz-juelich.de/mlz/lmfit.git
GIT_TAG main
PATCH_COMMAND ${LMFIT_PATCH_COMMAND}
UPDATE_DISCONNECTED 1
EXCLUDE_FROM_ALL 1)
endif()
# Disable what we don't need from lmfit
set(BUILD_TESTING
OFF
CACHE BOOL "")
set(LMFIT_CPPTEST
OFF
CACHE BOOL "")
set(LIB_MAN
OFF
CACHE BOOL "")
set(LMFIT_CPPTEST
OFF
CACHE BOOL "")
set(BUILD_SHARED_LIBS
OFF
CACHE BOOL "")
if(${CMAKE_VERSION} VERSION_LESS "3.28")
if(NOT lmfit_POPULATED)
FetchContent_Populate(lmfit)
add_subdirectory(${lmfit_SOURCE_DIR} ${lmfit_BINARY_DIR} EXCLUDE_FROM_ALL)
endif()
else()
FetchContent_MakeAvailable(lmfit)
endif()
set_property(TARGET lmfit PROPERTY POSITION_INDEPENDENT_CODE ON)
else()
find_package(lmfit REQUIRED)
endif()
if(AARE_FETCH_MINUIT2)
# We are building Minuit2 from sources.
@@ -380,6 +328,22 @@ else()
target_compile_options(aare_compiler_flags INTERFACE -Werror)
endif()
if(AARE_TUNE_LOCAL)
include(CheckCXXCompilerFlag)
check_cxx_compiler_flag("-march=native" AARE_HAS_MARCH_NATIVE)
check_cxx_compiler_flag("-mtune=native" AARE_HAS_MTUNE_NATIVE)
if(AARE_HAS_MARCH_NATIVE AND AARE_HAS_MTUNE_NATIVE)
message(STATUS "Tuning for the local CPU: -march=native -mtune=native")
target_compile_options(aare_compiler_flags INTERFACE -march=native
-mtune=native)
else()
message(
WARNING
"AARE_TUNE_LOCAL requested but the compiler does not support -march=native/-mtune=native. Ignoring."
)
endif()
endif()
endif() # GCC/Clang specific
if(AARE_PYTHON_BINDINGS)
@@ -423,7 +387,9 @@ set(PUBLICHEADERS
include/aare/Models.hpp
include/aare/FileInterface.hpp
include/aare/FilePtr.hpp
include/aare/FastPedestal.hpp
include/aare/Frame.hpp
include/aare/MultiThreadedFileReader.hpp
include/aare/hist/PixelHistogram.hpp
include/aare/hist/PixelHistogramImpl.hpp
include/aare/hist/PedestalTrackingPixelHistogram.hpp
@@ -442,7 +408,9 @@ set(PUBLICHEADERS
include/aare/RawMasterFile.hpp
include/aare/RawSubFile.hpp
include/aare/VarClusterFinder.hpp
include/aare/utils/task.hpp)
include/aare/utils/task.hpp
include/aare/utils/ifstream_helpers.hpp
include/aare/utils/math_helpers.hpp)
if(AARE_CUDA)
list(APPEND PUBLICHEADERS include/aare/ClusterFinderCUDA.hpp
@@ -464,6 +432,7 @@ set(SourceFiles
${CMAKE_CURRENT_SOURCE_DIR}/src/Frame.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/Interpolator.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/JungfrauDataFile.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/MultiThreadedFileReader.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/NumpyFile.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/NumpyHelpers.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/PixelMap.cpp
@@ -486,8 +455,7 @@ target_link_libraries(
aare_core
PUBLIC fmt::fmt nlohmann_json::nlohmann_json ${STD_FS_LIB} # from
# helpers.cmake
PRIVATE aare_compiler_flags Threads::Threads $<BUILD_INTERFACE:lmfit>
$<BUILD_INTERFACE:aare::Minuit2>)
PRIVATE aare_compiler_flags Threads::Threads $<BUILD_INTERFACE:aare::Minuit2>)
target_include_directories(
aare_core SYSTEM
@@ -530,6 +498,7 @@ if(AARE_TESTS)
${CMAKE_CURRENT_SOURCE_DIR}/src/defs.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/decode.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/Dtype.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/Fit.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/Frame.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/DetectorGeometry.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/Interpolation.test.cpp
@@ -545,7 +514,9 @@ if(AARE_TESTS)
${CMAKE_CURRENT_SOURCE_DIR}/src/Pedestal.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/hist/PixelHistogramImpl.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/hist/PixelHistogram.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/hist/PedestalTrackingPixelHistogram.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/JungfrauDataFile.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/MultiThreadedFileReader.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/NumpyFile.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/NumpyHelpers.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/RawFile.test.cpp
@@ -591,9 +562,8 @@ endif()
add_custom_target(
check-format
COMMAND
find \\ (-name "*.cpp" -o -name "*.hpp" \\) -not -path "./build/*" | xargs
-I {} -n 1 -P 10 bash -c
"clang-format -Werror -style=\"file:.clang-format\" {} | diff {} -"
bash -c
[=[ find \( -name "*.cpp" -o -name "*.hpp" \) -not -path "./build/*" | xargs -I {} -n 1 -P 10 bash -c "clang-format -Werror -style=\"file:.clang-format\" {} | diff {} -" ]=]
WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}
COMMENT "Checking code formatting with clang-format"
VERBATIM)
@@ -601,8 +571,8 @@ add_custom_target(
add_custom_target(
format-files
COMMAND
find \\ (-name "*.cpp" -o -name "*.hpp" \\) -not -path "./build/*" | xargs
-I {} -n 1 -P 10 bash -c "clang-format -i -style=\"file:.clang-format\" {}"
bash -c
[=[ find \( -name "*.cpp" -o -name "*.hpp" \) -not -path "./build/*" | xargs -I {} -n 1 -P 10 bash -c "clang-format -i -style=\"file:.clang-format\" {}" ]=]
WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}
COMMENT "Formatting with clang-format"
VERBATIM)
@@ -617,18 +587,9 @@ endif()
add_custom_target(
clang-tidy
COMMAND
find \\ (-path
"./src/*"
-a
-not
-path
"./src/python/*"
-a
\\
(-name "*.cpp" -not -name "*.test.cpp" \\)
\\) -not -name "CircularFifo.hpp" -not -name
"ProducerConsumerQueue.hpp" -not -name "VariableSizeClusterFinder.hpp" |
xargs -I {} -n 1 -P 10 bash -c
find -path "./src/*" -not -path "./src/python/*" -not -name "*.test.cpp"
-not -name "CircularFifo.hpp" -not -name "ProducerConsumerQueue.hpp" -not
-name "VariableSizeClusterFinder.hpp" | xargs -I {} -n 1 -P 10 bash -c
"${CLANG_TIDY_COMMAND} --config-file=.clang-tidy -p build {}"
WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}
COMMENT "linting with clang-tidy"
RELEASE.md
+40 -4
View File
@@ -1,5 +1,45 @@
# Release notes
## Next
### New Features:
- Added ``FastPedestal`` in C++ and Python for per-pixel running mean,
population variance, and standard deviation. It supports exponentially
weighted updates, initialization from files, direct subtraction from NumPy
arrays, and ``float64``, ``float32``, and ``int16`` output types.
- Added the ``Pedestal_i16`` Python binding alongside
``FastPedestal_d``, ``FastPedestal_f``, and ``FastPedestal_i16``.
- ``PedestalTrackingPixelHistogram.fill_from_file()`` now uses parallel,
double-buffered file reading and accepts ``reader_threads`` and
``reader_chunk_size`` tuning parameters.
- Added the ``AARE_TUNE_LOCAL`` CMake option to build with ``-march=native``
and ``-mtune=native`` when supported. Binaries built with this option are
specific to the local CPU and may not be portable.
### API Changes:
- ``ClusterFinder`` now uses ``FastPedestal``. It must receive 1000 pedestal
frames before cluster finding; ``find_clusters()`` raises
an error until initialization is complete. Added ``update_threshold()`` to
recompute the per-pixel detection thresholds.
- Added the ``queue_depth`` constructor argument to ``ClusterFinderMT`` to
configure the number of preallocated frame buffers per worker thread.
- Exposed ``ClusterVector.empty()`` in the Python API.
- Added ClusterVector.estimate_n_clusters
- Removed the lmfit dependency and the legacy ``fit_gaus``, ``fit_pol1``,
``fit_scurve``, and ``fit_scurve2`` APIs. Use ``Gaussian``, ``Pol1``,
``RisingScurve``, or ``FallingScurve`` and call ``model.fit(...)`` (or
``fit(model, ...)``) instead.
- Removed the legacy ``gaus``, ``pol1``, ``scurve``, and ``scurve2`` function
evaluators. Model objects are callable and provide the replacement, for
example ``Gaussian()(x, par)``.
- ``NDView<T, Ndim>`` now converts to ``NDView<const T, Ndim>``;
``expand4to8bit`` and ``expand24to32bit`` accept const input views.
### Bugfixes:
- Fixed broken reading of old (pre reordering) Moench03
## 2026.7.2
@@ -148,7 +188,3 @@ dhanya.thattil@psi.ch
benchmarks/CMakeLists.txt
+2 -2
View File
@@ -18,8 +18,8 @@ FetchContent_MakeAvailable(benchmark)
add_executable(benchmarks)
target_sources(
benchmarks PRIVATE ndarray_benchmark.cpp calculateeta_benchmark.cpp
reduce_benchmark.cpp)
benchmarks PRIVATE ndarray_benchmark.cpp ndview_benchmark.cpp
calculateeta_benchmark.cpp reduce_benchmark.cpp)
# Link Google Benchmark and other necessary libraries
target_link_libraries(benchmarks PRIVATE benchmark::benchmark aare_core
benchmarks/fit_benchmark.cpp
+3 -22
View File
@@ -75,24 +75,7 @@ static void report_accuracy(benchmark::State &state, const TestCase &tc,
// Benchmarks
// ----------
// 1. lmcurve
static void BM_FitGausLm(benchmark::State &state) {
const auto &tc = get_test_cases()[state.range(0)];
auto data = generate_gaussian_data(tc);
auto xv = data.x.view();
auto yv = data.y.view();
aare::NDArray<double, 1> result;
for (auto _ : state) {
result = aare::fit_gaus(xv, yv);
benchmark::DoNotOptimize(result.data());
}
report_accuracy(state, tc, result);
state.SetLabel(tc.name);
}
// 2. Minuit2, analytic gradient (no Hesse)
// Minuit2, analytic gradient (no Hesse)
static void BM_FitGausMinuitGrad(benchmark::State &state) {
const auto &tc = get_test_cases()[state.range(0)];
auto data = generate_gaussian_data(tc);
@@ -115,7 +98,7 @@ static void BM_FitGausMinuitGrad(benchmark::State &state) {
state.SetLabel(tc.name);
}
// 3. Minuit2, analytic gradient + Hesse
// Minuit2, analytic gradient + Hesse
static void BM_FitGausMinuitGradHesse(benchmark::State &state) {
const auto &tc = get_test_cases()[state.range(0)];
auto data = generate_gaussian_data(tc);
@@ -145,8 +128,6 @@ static void BM_FitGausMinuitGradHesse(benchmark::State &state) {
state.SetLabel(tc.name);
}
BENCHMARK(BM_FitGausLm)->DenseRange(0, 5)->Unit(benchmark::kMicrosecond);
BENCHMARK(BM_FitGausMinuitGrad)
->DenseRange(0, 5)
->Unit(benchmark::kMicrosecond);
@@ -155,4 +136,4 @@ BENCHMARK(BM_FitGausMinuitGradHesse)
->DenseRange(0, 5)
->Unit(benchmark::kMicrosecond);
BENCHMARK_MAIN();
BENCHMARK_MAIN();
benchmarks/ndview_benchmark.cpp
+17
View File
@@ -0,0 +1,17 @@
// SPDX-License-Identifier: MPL-2.0
#include "aare/NDArray.hpp"
#include "aare/NDView.hpp"
#include <benchmark/benchmark.h>
using aare::NDArray;
using aare::NDView;
static void BM_CreateNDView(benchmark::State &st) {
NDArray<int, 2> arr{{1024, 1024}, 0};
for (auto _ : st) {
// This code gets timed
auto res = arr.view();
benchmark::DoNotOptimize(res);
}
}
BENCHMARK(BM_CreateNDView);
docs/src/MultiThreadedFileReader.rst
+37
View File
@@ -0,0 +1,37 @@
MultiThreadedFileReader
=======================
``MultiThreadedFileReader`` reads one chunk per worker on each call. Each
worker owns an independent :cpp:class:`aare::File`, while all workers write
into non-overlapping regions of the destination buffer. The output is ordered
by frame index and successive calls advance through the file.
.. code-block:: cpp
#include "aare/MultiThreadedFileReader.hpp"
// Four workers, chunks of 128 frames, and at most 10,000 frames.
aare::experimental::MultiThreadedFileReader reader(path, 4, 128, 10'000);
while (reader.remaining_frames() != 0) {
// Contains at most 4 * 128 frames.
auto batch = reader.read();
process(batch);
}
Omit the final argument to read every frame in the source. An explicit value
of zero requests an empty result. The low-level ``read_into`` overload avoids
an allocation when the caller already owns a buffer of at least
``reader.next_read_bytes()`` bytes. Use ``read_all()`` to read every frame
remaining from the current position, and ``seek()`` to reposition the reader.
Call ``close()`` to release all worker file handles early.
.. note::
Multiple workers do not guarantee faster reads. Performance depends on the
storage device and file format, so the thread count and chunk size should be
benchmark-driven.
.. doxygenclass:: aare::experimental::MultiThreadedFileReader
:members:
:undoc-members:
:private-members:
docs/src/Requirements.rst
+2 -2
View File
@@ -14,7 +14,7 @@ Requirements
To simplify deployment we build and statically link a few libraries.
- fmt
- lmfit - https://jugit.fz-juelich.de/mlz/lmfit
- Minuit2
- nlohmann_json
- pybind11
- ZeroMQ
@@ -23,4 +23,4 @@ To simplify deployment we build and statically link a few libraries.
- Sphinx
- Breathe
- Doxygen
- Doxygen
docs/src/index.rst
+3
View File
@@ -28,7 +28,9 @@ AARE
pycalibration
python/cluster/index
python/file/index
python/experimental/index
python/histogram/index
python/pedestal/index
pyFit
@@ -42,6 +44,7 @@ AARE
NDView
Frame
File
MultiThreadedFileReader
Dtype
Cluster
ClusterFinder
docs/src/pyFit.rst
+20 -9
View File
@@ -1,19 +1,30 @@
Fit
========
Fitting
-------
.. py:currentmodule:: aare
Aare fits one-dimensional scans and three-dimensional pixel data with
Minuit2. Create a model object and call its :meth:`fit` method::
**Functions**
model = Gaussian(compute_errors=True)
result = model.fit(x, y, y_err)
.. autofunction:: gaus
The model object is also callable, which evaluates it at the supplied points::
.. autofunction:: pol1
fitted_y = model(x, result["par"])
The available models are ``Gaussian``, ``GaussianErfcPlateau``,
``GaussianChargeSharing``, ``GaussianChargeSharingKb``, ``Pol1``, ``Pol2``,
``RisingScurve``, and ``FallingScurve``. The module-level :func:`fit` function
accepts the same model objects when a functional interface is preferred.
**Fitting**
For three-dimensional data, pass an array with shape
``(rows, columns, scan_points)`` and select the worker count with
``n_threads``::
.. autofunction:: fit_gaus
result = model.fit(x, image_data, image_errors, n_threads=8)
.. autofunction:: fit_pol1
The result dictionary contains ``par`` and ``chi2``. It also contains
``par_err`` when ``compute_errors`` is enabled.
.. autofunction:: fit
docs/src/python/experimental/MultiThreadedFileReader.rst
+30
View File
@@ -0,0 +1,30 @@
MultiThreadedFileReader
=======================
.. py:currentmodule:: aare.experimental
The reader returns a NumPy array with shape ``(frames, rows, cols)`` and
preserves the source pixel dtype. Each iteration reads at most
``n_threads * chunk_size`` frames—one chunk per worker. File I/O runs with the
Python GIL released.
.. code-block:: python
from aare.experimental import MultiThreadedFileReader
with MultiThreadedFileReader(
"frames.npy", n_threads=4, chunk_size=128, total_frames=10_000
) as reader:
for frames in reader:
process(frames)
Call ``read()`` directly for the next batch, or ``read_all()`` for all frames
remaining from the current position. ``tell()`` and ``seek()`` expose the
iteration position. The context manager closes all worker files on exit.
``close()`` is also available for explicit cleanup and may be called
repeatedly.
.. autoclass:: MultiThreadedFileReader
:members:
:undoc-members:
:show-inheritance:
docs/src/python/experimental/index.rst
+10
View File
@@ -0,0 +1,10 @@
Experimental
============
APIs in this module may change without notice.
.. toctree::
:caption: Experimental
:maxdepth: 1
MultiThreadedFileReader
docs/src/python/histogram/pyPedestalTrackingPixelHistogram.rst
+8
View File
@@ -15,6 +15,14 @@ Use ``push_pedestal_no_update()`` to seed the pedestal estimate, then
asynchronous fills are drained by ``flush()``, and snapshot methods such as
``values()`` and ``pedestal_mean()`` return numpy arrays.
``fill_from_file()`` uses parallel file-reader workers and a double-buffered
pipeline. After the initial batch has been read, histogram processing of one
batch overlaps reading of the next. ``reader_threads`` and
``reader_chunk_size`` tune the I/O stage independently of the histogram worker
count. Two fixed-capacity buffers are allocated once and reused by alternating
their read and histogram roles. Their approximate memory use is ``2 *
reader_threads * reader_chunk_size * rows * cols * sizeof(uint16)``.
.. py:currentmodule:: aare
.. autoclass:: PedestalTrackingPixelHistogram
docs/src/python/pedestal/index.rst
+9
View File
@@ -0,0 +1,9 @@
Pedestal
========
.. toctree::
:caption: Pedestal
:maxdepth: 1
pyFastPedestal
pyPedestal
docs/src/python/pedestal/pyFastPedestal.rst
+83
View File
@@ -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
``add_init_frame()``. Once ``ready`` is true, use ``push_ema()`` to update the exponential
moving average initialized by the mean and with smoothing factor 1/n_samples.
.. warning::
FastPedestal is not usable until you have added ``n_samples`` initial frames with ``add_init_frame(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.add_init_frame(frame)
# Now we can push a frame for pedestal update
if pedestal.ready:
pedestal.push_ema(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:
docs/src/python/pedestal/pyPedestal.rst
+42
View File
@@ -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:
include/aare/Backoff.hpp
+48
View File
@@ -0,0 +1,48 @@
// SPDX-License-Identifier: MPL-2.0
#pragma once
#include <chrono>
#include <thread>
namespace aare {
/**
* @brief Hint the CPU that this is a spin-wait. Prefer this over yield() when
* the wait is expected to be short.
*/
inline void cpu_relax() noexcept {
#if defined(__x86_64__) || defined(__i386__)
__builtin_ia32_pause();
#elif defined(__aarch64__)
asm volatile("yield" ::: "memory");
#else
std::this_thread::yield();
#endif
}
/**
* @brief Escalating wait for producer/consumer idle and backpressure loops.
*
* Starts with a pause instruction (sub-microsecond, no syscall), then yield,
* then a short sleep. Call reset() whenever work arrives so a busy pipeline
* never leaves the pause tier.
*/
class Backoff {
int m_count{0};
public:
void reset() noexcept { m_count = 0; }
void pause() noexcept {
if (m_count < 64) {
cpu_relax();
} else if (m_count < 256) {
std::this_thread::yield();
} else {
std::this_thread::sleep_for(std::chrono::microseconds(50));
}
++m_count;
}
};
} // namespace aare
include/aare/CircularFifo.hpp
+42 -2
View File
@@ -1,11 +1,13 @@
// SPDX-License-Identifier: MPL-2.0
#pragma once
#include <cassert>
#include <chrono>
#include <fmt/color.h>
#include <fmt/format.h>
#include <memory>
#include <thread>
#include <utility>
#include "aare/ProducerConsumerQueue.hpp"
@@ -29,6 +31,23 @@ template <class ItemType> class CircularFifo {
}
}
/**
* @brief Construct a fifo and seed the free list using a factory.
* @param size number of items circulating in the fifo
* @param make_item callable invoked as make_item(i) for each slot index i
*
* Use this instead of CircularFifo(size) when the items need to be
* initialized, for example to hold a preallocated buffer or to carry
* their own slot index.
*/
template <class F>
CircularFifo(uint32_t size, F make_item)
: fifo_size(size), free_slots(size + 1), filled_slots(size + 1) {
for (size_t i = 0; i < fifo_size; ++i) {
free_slots.write(make_item(i));
}
}
bool next() {
// TODO! avoid default constructing ItemType
ItemType it;
@@ -47,6 +66,13 @@ template <class ItemType> class CircularFifo {
auto numFreeSlots() const noexcept { return free_slots.sizeGuess(); }
auto isFull() const noexcept { return filled_slots.isFull(); }
/**
* @brief True if there are no filled slots waiting to be consumed.
* @note Prefer this over numFilledSlots() == 0 since sizeGuess() may
* under-report when called from the producing thread.
*/
auto isEmpty() const noexcept { return filled_slots.isEmpty(); }
ItemType pop_free() {
ItemType v;
while (!free_slots.read(v))
@@ -75,8 +101,22 @@ template <class ItemType> class CircularFifo {
ItemType *frontPtr() { return filled_slots.frontPtr(); }
// TODO! Add function to move item from filled to free to be used
// with the frontPtr function
/**
* @brief Return the front filled item to the free list. To be used
* together with frontPtr() once the item has been consumed in place.
* @warning The fifo must not be empty when calling this.
*
* The item is written to the free list before it is popped from the
* filled list, so it can never be dropped. The write cannot fail: both
* queues hold size + 1 slots while only size items circulate.
*/
void recycle_front() {
ItemType *it = filled_slots.frontPtr();
assert(it != nullptr);
[[maybe_unused]] const bool ok = free_slots.write(std::move(*it));
assert(ok);
filled_slots.popFront();
}
template <class... Args> void push_value(Args &&...recordArgs) {
while (!filled_slots.write(std::forward<Args>(recordArgs)...))
include/aare/ClusterCollector.hpp
+33 -13
View File
@@ -2,7 +2,11 @@
#pragma once
#include <atomic>
#include <thread>
#include <type_traits>
#include <utility>
#include <vector>
#include "aare/Backoff.hpp"
#include "aare/ClusterFinderMT.hpp"
#include "aare/ClusterVector.hpp"
#include "aare/ProducerConsumerQueue.hpp"
@@ -13,48 +17,64 @@ namespace aare {
template <typename ClusterType,
typename = std::enable_if_t<is_cluster_v<ClusterType>>>
class ClusterCollector {
ProducerConsumerQueue<ClusterVector<ClusterType>> *m_source;
using SourceQueue = ProducerConsumerQueue<ClusterVector<ClusterType>>;
SourceQueue *m_source;
std::atomic<bool> m_stop_requested{false};
std::atomic<bool> m_stopped{true};
std::chrono::milliseconds m_default_wait{1};
std::thread m_thread;
std::vector<ClusterVector<ClusterType>> m_clusters;
void process() {
m_stopped = false;
m_stopped.store(false, std::memory_order_release);
fmt::print("ClusterCollector started\n");
while (!m_stop_requested || !m_source->isEmpty()) {
Backoff backoff;
while (true) {
if (ClusterVector<ClusterType> *clusters = m_source->frontPtr();
clusters != nullptr) {
backoff.reset();
m_clusters.push_back(std::move(*clusters));
m_source->popFront();
} else {
std::this_thread::sleep_for(m_default_wait);
continue;
}
if (m_stop_requested.load(std::memory_order_acquire))
break;
backoff.pause();
}
fmt::print("ClusterCollector stopped\n");
m_stopped = true;
m_stopped.store(true, std::memory_order_release);
}
public:
ClusterCollector(ClusterFinderMT<ClusterType, uint16_t, double> *source) {
m_source = source->sink();
explicit ClusterCollector(SourceQueue *source) : m_source(source) {
m_thread =
std::thread(&ClusterCollector::process,
this); // only one process does that so why isnt it
// automatically written to m_cluster in collect
// - instead of writing first to m_sink?
}
template <typename Finder,
std::enable_if_t<
std::is_same_v<decltype(std::declval<Finder &>().sink()),
SourceQueue *>,
int> = 0>
explicit ClusterCollector(Finder *source)
: ClusterCollector(source->sink()) {}
void stop() {
m_stop_requested = true;
m_thread.join();
m_stop_requested.store(true, std::memory_order_release);
if (m_thread.joinable())
m_thread.join();
}
std::vector<ClusterVector<ClusterType>> steal_clusters() {
if (!m_stopped) {
if (!m_stopped.load(std::memory_order_acquire)) {
throw std::runtime_error("ClusterCollector is still running");
}
return std::move(m_clusters);
}
};
} // namespace aare
} // namespace aare
include/aare/ClusterFile.hpp
+13
View File
@@ -144,6 +144,19 @@ class ClusterFile {
*/
size_t chunk_size() const { return m_chunk_size; }
/**
* @brief Estimate the number of clusters in the file from its size
*
* Frame headers are included in the estimate, so the result may be slightly
* larger than the actual number of clusters. The file position is not
* changed.
*
* @throws std::runtime_error if the file is not opened for reading
*/
size_t estimate_n_clusters() const {
return std::filesystem::file_size(m_filename) / sizeof(ClusterType);
}
/**
* @brief Set the region of interest to use when reading
* clusters. If set only clusters within the ROI will be
include/aare/ClusterFinder.hpp
+223 -110
View File
@@ -3,6 +3,7 @@
#include "aare/ClusterFile.hpp"
#include "aare/ClusterVector.hpp"
#include "aare/Dtype.hpp"
#include "aare/FastPedestal.hpp"
#include "aare/NDArray.hpp"
#include "aare/NDView.hpp"
#include "aare/Pedestal.hpp"
@@ -18,6 +19,13 @@ struct no_2x2_cluster {
ClusterType::cluster_size_x > 2 && ClusterType::cluster_size_y > 2;
};
/**
* @brief Find fixed-size photon clusters using a per-pixel pedestal and noise
* threshold.
* @tparam ClusterType Output cluster type; both dimensions must exceed 2.
* @tparam FRAME_TYPE Input pixel type.
* @tparam PEDESTAL_TYPE Type used for pedestal and threshold calculations.
*/
template <typename ClusterType = Cluster<int32_t, 3, 3>,
typename FRAME_TYPE = uint16_t, typename PEDESTAL_TYPE = double,
typename = std::enable_if_t<no_2x2_cluster<ClusterType>::value>>
@@ -26,51 +34,86 @@ class ClusterFinder {
PEDESTAL_TYPE m_nSigma;
const PEDESTAL_TYPE c2;
const PEDESTAL_TYPE c3;
Pedestal<PEDESTAL_TYPE> m_pedestal;
FastPedestal<PEDESTAL_TYPE> m_pedestal;
ClusterVector<ClusterType> m_clusters;
static const uint8_t ClusterSizeX = ClusterType::cluster_size_x;
static const uint8_t ClusterSizeY = ClusterType::cluster_size_y;
using CT = typename ClusterType::value_type;
NDArray<PEDESTAL_TYPE, 2> m_threshold;
NDArray<PEDESTAL_TYPE, 2> m_pd_corrected_frame;
public:
/**
* @brief Construct a new ClusterFinder object
* @param image_size size of the image
* @param cluster_size size of the cluster (x, y)
* @param nSigma number of sigma above the pedestal to consider a photon
* @param capacity initial capacity of the cluster vector
*
* @brief Construct a cluster finder with an empty pedestal.
* @param image_size Image shape as (rows, columns).
* @param nSigma Per-pixel noise threshold multiplier.
* @param capacity Initial cluster-vector capacity.
*/
ClusterFinder(Shape<2> image_size, PEDESTAL_TYPE nSigma = 5.0,
size_t capacity = 1000000)
: m_image_size(image_size), m_nSigma(nSigma),
c2(sqrt((ClusterSizeY + 1) / 2 * (ClusterSizeX + 1) / 2)),
c3(sqrt(ClusterSizeX * ClusterSizeY)),
m_pedestal(image_size[0], image_size[1]), m_clusters(capacity) {
m_pedestal(image_size[0], image_size[1]), m_clusters(capacity),
m_threshold({image_size[0], image_size[1]}, 0),
m_pd_corrected_frame({image_size[0], image_size[1]}, 0) {
LOG(logDEBUG) << "ClusterFinder: "
<< "image_size: " << image_size[0] << "x" << image_size[1]
<< ", nSigma: " << nSigma << ", capacity: " << capacity;
}
void set_nSigma(PEDESTAL_TYPE nSigma) { m_nSigma = nSigma; }
PEDESTAL_TYPE get_nSigma() const { return m_nSigma; }
void push_pedestal_frame(NDView<FRAME_TYPE, 2> frame) {
m_pedestal.push(frame);
/**
* @brief Set the noise multiplier used for threshold calculation and
* recompute the threshold.
* @param nSigma New per-pixel noise multiplier.
*/
void set_nSigma(PEDESTAL_TYPE nSigma) {
m_nSigma = nSigma;
update_threshold();
}
NDArray<PEDESTAL_TYPE, 2> pedestal() { return m_pedestal.mean(); }
NDArray<PEDESTAL_TYPE, 2> noise() { return m_pedestal.std(); }
void clear_pedestal() { m_pedestal.clear(); }
/** @brief Return the current noise multiplier used for threshold
* calculation. */
PEDESTAL_TYPE get_nSigma() const { return m_nSigma; }
/**
* @brief Move the clusters from the ClusterVector in the ClusterFinder to a
* new ClusterVector and return it.
* @param realloc_same_capacity if true the new ClusterVector will have the
* same capacity as the old one
* @brief Add a dark frame to the pedestal estimator.
*
* The threshold is initialized automatically when the pedestal first
* becomes ready. Later frames update the ready pedestal.
* @param frame Dark frame matching the configured image shape.
* @throws std::runtime_error if the frame shape does not match.
*/
void push_pedestal_frame(NDView<FRAME_TYPE, 2> frame) {
if (!m_pedestal.ready()) {
m_pedestal.add_init_frame(frame);
// Initialize the threshold when the pedestal becomes ready.
if (m_pedestal.ready()) {
update_threshold();
}
} else {
m_pedestal.push_ema(frame);
}
}
/** @brief Return a copy of the per-pixel pedestal mean. */
NDArray<PEDESTAL_TYPE, 2> pedestal() { return m_pedestal.mean(); }
/** @brief Return the per-pixel pedestal standard deviation (noise). */
NDArray<PEDESTAL_TYPE, 2> noise() { return m_pedestal.std(); }
/** @brief Clear the pedestal and mark it as not ready. */
void clear_pedestal() { m_pedestal.clear(); }
/** @brief Recompute the threshold as noise multiplied by nSigma. */
void update_threshold() { m_threshold = m_pedestal.std() * m_nSigma; }
/**
* @brief Move out all accumulated clusters and reset the internal vector.
* @param realloc_same_capacity Preserve the previous capacity when true.
* @return The accumulated clusters and their frame metadata.
*/
ClusterVector<ClusterType>
steal_clusters(bool realloc_same_capacity = false) {
@@ -81,109 +124,179 @@ class ClusterFinder {
m_clusters = ClusterVector<ClusterType>{};
return tmp;
}
private:
/**
* @brief Process a single pixel: scan its cluster window, decide whether it
* is a photon or a pedestal value, and store the cluster if needed.
* @tparam CheckBounds Skip out-of-image neighbours when true; assume the
* complete window is in bounds when false. Skipped cluster values remain 0.
*/
template <bool CheckBounds>
void process_pixel(const NDView<FRAME_TYPE, 2> &frame, const int iy,
const int ix) {
constexpr int dy = ClusterSizeY / 2;
constexpr int dx = ClusterSizeX / 2;
constexpr int has_center_pixel_x = ClusterSizeX % 2;
constexpr int has_center_pixel_y = ClusterSizeY % 2;
PEDESTAL_TYPE max = std::numeric_limits<PEDESTAL_TYPE>::lowest();
PEDESTAL_TYPE total = 0;
const int cols = static_cast<int>(frame.shape(1));
const int rows = static_cast<int>(frame.shape(0));
const auto center =
(static_cast<std::size_t>(iy) * static_cast<std::size_t>(cols)) +
static_cast<std::size_t>(ix);
const auto *corrected = m_pd_corrected_frame.data();
const PEDESTAL_TYPE threshold = m_threshold.data()[center];
const PEDESTAL_TYPE value = corrected[center];
if (value < -threshold)
return; // NEGATIVE_PEDESTAL, nothing to do for this pixel
// TODO! No pedestal update???
for (int ir = -dy; ir < dy + has_center_pixel_y; ir++) {
const int y = iy + ir;
if constexpr (CheckBounds) {
if (y < 0 || y >= rows)
continue;
}
const auto *row = corrected + static_cast<std::size_t>(y) * cols;
for (int ic = -dx; ic < dx + has_center_pixel_x; ic++) {
const int x = ix + ic;
if constexpr (CheckBounds) {
if (x < 0 || x >= cols)
continue;
}
const PEDESTAL_TYPE val = row[x];
total += val;
max = std::max(max, val);
}
}
if ((max > threshold)) {
if (value < max)
return; // Not max go to the next pixel, no pedestal update
} else if (total > c3 * threshold) {
// pass, store the cluster below
} else {
m_pedestal.push_ema_unchecked(center, frame.data()[center]);
return; // It was a pedestal value nothing to store
}
// Store cluster
if (value == max) {
ClusterType cluster{};
cluster.x = ix;
cluster.y = iy;
int i = 0;
for (int ir = -dy; ir < dy + has_center_pixel_y; ir++) {
const int y = iy + ir;
for (int ic = -dx; ic < dx + has_center_pixel_x; ic++, i++) {
const int x = ix + ic;
if constexpr (CheckBounds) {
if (x < 0 || x >= cols || y < 0 || y >= rows)
continue;
}
const PEDESTAL_TYPE corrected_value =
corrected[(static_cast<std::size_t>(y) * cols) + x];
// If the cluster type is an integral type, and the
// pedestal is a floating point type then we need to
// round the value before storing it
if constexpr (std::is_integral_v<CT> &&
std::is_floating_point_v<PEDESTAL_TYPE>) {
cluster.data[i] =
static_cast<CT>(std::lround(corrected_value));
}
// On the other hand if both are floating point or both
// are integral then we can just static cast directly
else {
cluster.data[i] = static_cast<CT>(corrected_value);
}
}
}
// Add the cluster to the output ClusterVector
m_clusters.push_back(cluster);
}
}
public:
/**
* @brief Find clusters in one frame and update eligible pedestal pixels.
* @param frame Input frame matching the configured image shape.
* @param frame_number Metadata assigned to the accumulated clusters.
* @pre frame has the same shape as image_size passed to the constructor.
* @throws std::runtime_error if the pedestal is not ready.
* @note Clusters accumulate until steal_clusters() is called.
*/
void find_clusters(NDView<FRAME_TYPE, 2> frame, uint64_t frame_number = 0) {
// // TODO! deal with even size clusters
// // currently 3,3 -> +/- 1
// // 4,4 -> +/- 2
int dy = ClusterSizeY / 2;
int dx = ClusterSizeX / 2;
int has_center_pixel_x =
ClusterSizeX %
2; // for even sized clusters there is no proper cluster center and
// even amount of pixels around the center
int has_center_pixel_y = ClusterSizeY % 2;
if (!m_pedestal.ready()) {
throw std::runtime_error(
"Pedestal is not ready, cannot find clusters");
}
constexpr int dy = ClusterSizeY / 2;
constexpr int dx = ClusterSizeX / 2;
constexpr int has_center_pixel_x = ClusterSizeX % 2;
constexpr int has_center_pixel_y = ClusterSizeY % 2;
// Largest neighbour offset below/right of the current pixel. Pixels
// further than this from an edge have their whole window in bounds.
constexpr int down = dy + has_center_pixel_y - 1;
constexpr int right = dx + has_center_pixel_x - 1;
m_clusters.set_frame_number(frame_number);
for (int iy = 0; iy < frame.shape(0); iy++) {
for (int ix = 0; ix < frame.shape(1); ix++) {
PEDESTAL_TYPE max = std::numeric_limits<FRAME_TYPE>::min();
PEDESTAL_TYPE total = 0;
const int rows = static_cast<int>(frame.shape(0));
const int cols = static_cast<int>(frame.shape(1));
// What can we short circuit here?
PEDESTAL_TYPE rms = m_pedestal.std(iy, ix);
PEDESTAL_TYPE value = (frame(iy, ix) - m_pedestal.mean(iy, ix));
// TODO! See if we can get the same performace using the operator-
// m_pd_corrected_frame = frame - m_pedestal.view();
if (value < -m_nSigma * rms)
continue; // NEGATIVE_PEDESTAL go to next pixel
// TODO! No pedestal update???
// here we should be able to safely assume that the frame and corrected
// frame have the same size
auto n_pixels = frame.size();
auto pd = m_pedestal.view().data();
auto corrected = m_pd_corrected_frame.data();
auto frame_data = frame.data();
for (ssize_t i = 0; i < n_pixels; i++) {
corrected[i] = static_cast<PEDESTAL_TYPE>(frame_data[i]) - pd[i];
}
for (int ir = -dy; ir < dy + has_center_pixel_y; ir++) {
for (int ic = -dx; ic < dx + has_center_pixel_x; ic++) {
if (ix + ic >= 0 && ix + ic < frame.shape(1) &&
iy + ir >= 0 && iy + ir < frame.shape(0)) {
PEDESTAL_TYPE val =
frame(iy + ir, ix + ic) -
m_pedestal.mean(iy + ir, ix + ic);
// Interior pixels can skip the per-neighbour bounds checks; pixels
// within dx/dy of an edge take the bounds-checked path. Iteration order
// (row-major, increasing ix) is preserved so results are identical.
const int ix_begin = dx;
const int ix_end = cols - right; // exclusive
total += val;
max = std::max(max, val);
}
}
}
for (int iy = 0; iy < rows; iy++) {
const bool interior_row = iy >= dy && iy < rows - down;
if ((max > m_nSigma * rms)) {
if (value < max)
continue; // Not max go to the next pixel
// but also no pedestal update
} else if (total > c3 * m_nSigma * rms) {
// pass
} else {
// m_pedestal.push(iy, ix, frame(iy, ix)); // Safe option
m_pedestal.push_fast(
iy, ix,
frame(iy,
ix)); // Assume we have reached n_samples in the
// pedestal, slight performance improvement
continue; // It was a pedestal value nothing to store
}
// Store cluster
if (value == max) {
ClusterType cluster{};
cluster.x = ix;
cluster.y = iy;
// Fill the cluster data since we have a photon to store
// It's worth redoing the look since most of the time we
// don't have a photon
int i = 0;
for (int ir = -dy; ir < dy + has_center_pixel_y; ir++) {
for (int ic = -dx; ic < dx + has_center_pixel_x; ic++) {
if (ix + ic >= 0 && ix + ic < frame.shape(1) &&
iy + ir >= 0 && iy + ir < frame.shape(0)) {
// If the cluster type is an integral type, and
// the pedestal is a floating point type then we
// need to round the value before storing it
if constexpr (std::is_integral_v<CT> &&
std::is_floating_point_v<
PEDESTAL_TYPE>) {
auto tmp = std::lround(
frame(iy + ir, ix + ic) -
m_pedestal.mean(iy + ir, ix + ic));
cluster.data[i] = static_cast<CT>(tmp);
}
// On the other hand if both are floating point
// or both are integral then we can just static
// cast directly
else {
auto tmp =
frame(iy + ir, ix + ic) -
m_pedestal.mean(iy + ir, ix + ic);
cluster.data[i] = static_cast<CT>(tmp);
}
}
i++;
}
}
// Add the cluster to the output ClusterVector
m_clusters.push_back(cluster);
}
if (!interior_row || ix_begin >= ix_end) {
for (int ix = 0; ix < cols; ix++)
process_pixel<true>(frame, iy, ix);
continue;
}
for (int ix = 0; ix < ix_begin; ix++)
process_pixel<true>(frame, iy, ix);
for (int ix = ix_begin; ix < ix_end; ix++)
process_pixel<false>(frame, iy, ix);
for (int ix = ix_end; ix < cols; ix++)
process_pixel<true>(frame, iy, ix);
}
}
};
} // namespace aare
} // namespace aare
include/aare/ClusterFinderMT.hpp
+134 -39
View File
@@ -6,11 +6,13 @@
#include <thread>
#include <vector>
#include "aare/Backoff.hpp"
#include "aare/CircularFifo.hpp"
#include "aare/ClusterFinder.hpp"
#include "aare/NDArray.hpp"
#include "aare/ProducerConsumerQueue.hpp"
#include "aare/logger.hpp"
#include <ctime>
namespace aare {
enum class FrameType {
@@ -18,10 +20,40 @@ enum class FrameType {
PEDESTAL,
};
struct FrameWrapper {
FrameType type;
uint64_t frame_number;
NDArray<uint16_t, 2> data;
/**
* @brief Ticket identifying a frame buffer in a FramePool. Trivially
* copyable, the buffer itself never travels through the queues.
*/
struct FrameRef {
FrameType type{};
uint32_t slot{};
uint64_t frame_number{};
};
/**
* @brief Fixed set of frame buffers allocated once at construction.
*
* Buffers are addressed by slot index and are never moved or reallocated, so
* the queues only need to pass a FrameRef around. This keeps the per frame
* allocation out of the hot path entirely.
*/
class FramePool {
std::vector<NDArray<uint16_t, 2>> m_buffers;
public:
FramePool(size_t depth, Shape<2> shape) {
m_buffers.reserve(depth); // no reallocation, slots stay stable
for (size_t i = 0; i < depth; ++i) {
m_buffers.emplace_back(shape);
}
}
NDArray<uint16_t, 2> &operator[](uint32_t slot) { return m_buffers[slot]; }
const NDArray<uint16_t, 2> &operator[](uint32_t slot) const {
return m_buffers[slot];
}
size_t size() const { return m_buffers.size(); }
};
/**
@@ -42,17 +74,18 @@ class ClusterFinderMT {
size_t m_current_thread{0};
size_t m_n_threads{0};
using Finder = ClusterFinder<ClusterType, FRAME_TYPE, PEDESTAL_TYPE>;
using InputQueue = ProducerConsumerQueue<FrameWrapper>;
using InputQueue = CircularFifo<FrameRef>;
using OutputQueue = ProducerConsumerQueue<ClusterVector<ClusterType>>;
std::vector<std::unique_ptr<InputQueue>> m_input_queues;
std::vector<std::unique_ptr<OutputQueue>> m_output_queues;
std::vector<std::unique_ptr<FramePool>> m_frame_pools;
OutputQueue m_sink{1000}; // All clusters go into this queue
std::vector<std::unique_ptr<Finder>> m_cluster_finders;
std::vector<std::thread> m_threads;
std::thread m_collect_thread;
std::chrono::milliseconds m_default_wait{1};
std::chrono::microseconds m_default_wait{50};
private:
std::atomic<bool> m_stop_requested{false};
@@ -65,31 +98,35 @@ class ClusterFinderMT {
void process(int thread_id) {
auto cf = m_cluster_finders[thread_id].get();
auto q = m_input_queues[thread_id].get();
bool realloc_same_capacity = true;
auto *pool = m_frame_pools[thread_id].get();
Backoff backoff;
while (!m_stop_requested || !q->isEmpty()) {
if (FrameWrapper *frame = q->frontPtr(); frame != nullptr) {
if (FrameRef *ref = q->frontPtr(); ref != nullptr) {
backoff.reset();
auto view = (*pool)[ref->slot].view();
switch (frame->type) {
case FrameType::DATA:
cf->find_clusters(frame->data.view(), frame->frame_number);
switch (ref->type) {
case FrameType::DATA: {
cf->find_clusters(view, ref->frame_number);
// Steal before the write so a failed write cannot drop the
// clusters by re-stealing an empty vector on retry.
auto clusters = cf->steal_clusters(true);
while (!m_output_queues[thread_id]->write(
cf->steal_clusters(realloc_same_capacity))) {
std::this_thread::sleep_for(m_default_wait);
std::move(clusters))) {
backoff.pause();
}
break;
}
case FrameType::PEDESTAL:
m_cluster_finders[thread_id]->push_pedestal_frame(
frame->data.view());
m_cluster_finders[thread_id]->push_pedestal_frame(view);
break;
}
// frame is processed now discard it
m_input_queues[thread_id]->popFront();
// frame is processed, hand the buffer back to the free list
q->recycle_front();
} else {
std::this_thread::sleep_for(m_default_wait);
backoff.pause();
}
}
}
@@ -100,18 +137,24 @@ class ClusterFinderMT {
*/
void collect() {
bool empty = true;
Backoff backoff;
while (!m_stop_requested || !empty || !m_processing_threads_stopped) {
empty = true;
bool moved_any = false;
for (auto &queue : m_output_queues) {
if (!queue->isEmpty()) {
while (!m_sink.write(std::move(*queue->frontPtr()))) {
std::this_thread::sleep_for(m_default_wait);
while (auto *front = queue->frontPtr()) {
while (!m_sink.write(std::move(*front))) {
backoff.pause();
}
queue->popFront();
empty = false;
moved_any = true;
}
}
empty = !moved_any;
if (moved_any) {
backoff.reset();
} else {
backoff.pause();
}
}
}
@@ -124,16 +167,22 @@ class ClusterFinderMT {
* @param capacity initial capacity of the cluster vector. Should match
* expected number of clusters in a frame per frame.
* @param n_threads number of threads to use
* @param queue_depth number of frame buffers per thread. These are
* allocated once and recycled, so the total resident frame memory is
* n_threads * queue_depth * frame size. Keeping the in flight data below
* the L3 size keeps the per frame copy cheap.
*/
ClusterFinderMT(Shape<2> image_size, PEDESTAL_TYPE nSigma = 5.0,
size_t capacity = 2000, size_t n_threads = 3)
size_t capacity = 2000, size_t n_threads = 3,
size_t queue_depth = 16)
: m_n_threads(n_threads) {
LOG(logDEBUG1) << "ClusterFinderMT: "
<< "image_size: " << image_size[0] << "x"
<< image_size[1] << ", nSigma: " << nSigma
<< ", capacity: " << capacity
<< ", n_threads: " << n_threads;
<< ", n_threads: " << n_threads
<< ", queue_depth: " << queue_depth;
for (size_t i = 0; i < n_threads; i++) {
m_cluster_finders.push_back(
@@ -142,7 +191,13 @@ class ClusterFinderMT {
image_size, nSigma, capacity));
}
for (size_t i = 0; i < n_threads; i++) {
m_input_queues.emplace_back(std::make_unique<InputQueue>(200));
m_frame_pools.emplace_back(
std::make_unique<FramePool>(queue_depth, image_size));
m_input_queues.emplace_back(std::make_unique<InputQueue>(
static_cast<uint32_t>(queue_depth), [](size_t slot) {
return FrameRef{FrameType::DATA,
static_cast<uint32_t>(slot), 0};
}));
m_output_queues.emplace_back(std::make_unique<OutputQueue>(200));
}
// TODO! Should we start automatically?
@@ -212,13 +267,22 @@ class ClusterFinderMT {
* expected to be dark. No photon finding is done. Just pedestal update.
*/
void push_pedestal_frame(NDView<FRAME_TYPE, 2> frame) {
FrameWrapper fw{FrameType::PEDESTAL, 0,
NDArray(frame)}; // TODO! copies the data!
for (size_t i = 0; i < m_n_threads; ++i) {
auto *q = m_input_queues[i].get();
Backoff backoff;
for (auto &queue : m_input_queues) {
while (!queue->write(fw)) {
std::this_thread::sleep_for(m_default_wait);
FrameRef ref;
while (!q->try_pop_free(ref)) {
backoff.pause();
}
ref.type = FrameType::PEDESTAL;
ref.frame_number = 0;
(*m_frame_pools[i])[ref.slot].copy_from(frame);
// Cannot fail, the free list is what limits how many frames are
// in flight so there is always room in the filled list.
q->try_push_value(ref);
}
}
@@ -228,11 +292,29 @@ class ClusterFinderMT {
* @note Spin locks with a default wait if the queue is full.
*/
void find_clusters(NDView<FRAME_TYPE, 2> frame, uint64_t frame_number = 0) {
FrameWrapper fw{FrameType::DATA, frame_number,
NDArray(frame)}; // TODO! copies the data!
while (!m_input_queues[m_current_thread % m_n_threads]->write(fw)) {
std::this_thread::sleep_for(m_default_wait);
const size_t tid = m_current_thread % m_n_threads;
auto *q = m_input_queues[tid].get();
Backoff backoff;
FrameRef ref;
while (!q->try_pop_free(ref)) {
backoff.pause();
}
ref.type = FrameType::DATA;
ref.frame_number = frame_number;
// DualTimer dt;
(*m_frame_pools[tid])[ref.slot].copy_from(frame);
// auto [wall_ns, cpu_ns] = dt.elapsed();
// std::cerr << "ClusterFinderMT: find_clusters: copied frame "
// << frame_number << " took " << wall_ns/1000.0 << " wall_us
// and "
// << cpu_ns/1000.0 << " cpu_us" << std::endl;
// Cannot fail, the free list is what limits how many frames are in
// flight so there is always room in the filled list.
q->try_push_value(ref);
m_current_thread++;
}
@@ -245,6 +327,19 @@ class ClusterFinderMT {
}
}
/**
* @brief Recompute the threshold (nSigma * pedestal std) on all cluster
* finders. Requires the processing threads to be stopped.
*/
void update_threshold() {
if (!m_processing_threads_stopped) {
throw std::runtime_error("ClusterFinderMT is still running");
}
for (auto &cf : m_cluster_finders) {
cf->update_threshold();
}
}
/**
* @brief Return the pedestal currently used by the cluster finder
* @param thread_index index of the thread
include/aare/ClusterVector.hpp
+12 -17
View File
@@ -53,13 +53,8 @@ class ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> {
m_data.reserve(capacity);
}
// Move constructor
ClusterVector(ClusterVector &&other) noexcept
: m_data(std::move(other.m_data)),
m_frame_number(other.m_frame_number) {
// other.m_data.clear();
other.m_frame_number = 0;
}
ClusterVector(ClusterVector &&other) noexcept = default;
ClusterVector &operator=(ClusterVector &&other) noexcept = default;
/**
* @brief Create a copy of the clustervector by filtering clusters in the
@@ -82,16 +77,16 @@ class ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> {
return result;
}
// Move assignment operator
ClusterVector &operator=(ClusterVector &&other) noexcept {
if (this != &other) {
m_data = std::move(other.m_data);
m_frame_number = other.m_frame_number;
other.m_data.clear();
other.m_frame_number = 0;
}
return *this;
}
// // Move assignment operator
// ClusterVector &operator=(ClusterVector &&other) noexcept {
// if (this != &other) {
// m_data = other.m_data;
// m_frame_number = other.m_frame_number;
// other.m_data.clear();
// other.m_frame_number = 0;
// }
// return *this;
// }
/**
* @brief Sum the pixels in each cluster
include/aare/FastPedestal.hpp
+424
View File
@@ -0,0 +1,424 @@
// SPDX-License-Identifier: MPL-2.0
#pragma once
#include "aare/File.hpp"
#include "aare/Frame.hpp"
#include "aare/NDArray.hpp"
#include "aare/NDView.hpp"
#include <cstddef>
#include <sys/types.h>
namespace aare {
/**
* @brief Maintain per-pixel mean, population variance, and standard deviation.
*
* Initialization accumulates exactly n_samples frames. Subsequent push_ema()
* update the exponential moving average initialized with the mean using
* a smoothing factor of (1/ n_samples). Internal moments are stored
* in double precision.
* @tparam PEDESTAL_TYPE Type returned for the mean, variance, and standard
* deviation.
*/
template <typename PEDESTAL_TYPE> class FastPedestal {
// Did we accumulate enough samples and updated the mean?
bool m_ready = false;
uint32_t m_rows;
uint32_t m_cols;
uint32_t m_samples;
double m_inv_samples; // precompute 1/m_samples for faster division
uint32_t m_cur_samples = 0; // number of samples accumulated so far
// For cache locality we want to keep sum and sum2 close. Should
// improve performance for random access.
struct Entry {
double sum;
double sum2;
};
NDArray<Entry, 2> m_sum;
// Cache mean since it is used over and over in the ClusterFinder
// This optimization is related to the access pattern of the ClusterFinder
// Relies on having more reads than pushes to the pedestal
// But also makes sense when subtracting the pedestal from the frame
NDArray<PEDESTAL_TYPE, 2> m_mean;
// Helper function to convert row and column indices to a flat index
// used to provide both row column and flat index access to the pedestal
size_t rc_to_index(uint32_t row, uint32_t col) const {
return (static_cast<std::size_t>(row) * m_cols) + col;
}
public:
/**
* @brief Construct an empty pedestal that becomes ready after n_samples
* initialization frames.
* @param rows Number of image rows.
* @param cols Number of image columns.
* @param n_samples Number of initialization frames and reciprocal of the
* weight assigned to each subsequent value.
* @throws std::runtime_error if rows, cols, or n_samples is zero.
*/
FastPedestal(uint32_t rows, uint32_t cols, uint32_t n_samples = 1000)
: m_rows(rows), m_cols(cols), m_samples(n_samples),
m_inv_samples(1.0 / n_samples), m_sum({rows, cols}, Entry{0, 0}),
m_mean({rows, cols}, PEDESTAL_TYPE(0)) {
if (!(rows > 0 && cols > 0 && n_samples > 0)) {
throw std::runtime_error(
fmt::format("Invalid parameters for FastPedestal: rows={}, "
"cols={}, n_samples={} need to be positive",
rows, cols, n_samples));
}
}
~FastPedestal() = default;
/**
* @brief Return a non-owning view of the cached mean.
* @throws std::runtime_error if ready() is false.
* @note The caller must treat the data as read-only and must not retain the
* view after this object is destroyed, moved, or assigned.
*/
NDView<const PEDESTAL_TYPE, 2> view() const {
if (!ready()) {
throw std::runtime_error(
"Pedestal is not ready, cannot return view");
}
return m_mean.view();
}
/**
* @brief Return a copy of the cached mean.
* @throws std::runtime_error if ready() is false.
*/
NDArray<PEDESTAL_TYPE, 2> mean() {
if (!ready()) {
throw std::runtime_error(
"Pedestal is not ready, cannot return mean");
}
return m_mean;
}
/**
* @brief Return the cached mean at (row, col).
* @throws std::runtime_error if ready() is false or either index is out of
* range.
*/
PEDESTAL_TYPE mean(uint32_t row, uint32_t col) const {
if (!ready()) {
throw std::runtime_error(
"Pedestal is not ready, cannot return mean");
}
if (row >= m_rows || col >= m_cols) {
throw std::runtime_error(
fmt::format("Invalid indices for FastPedestal mean: row={}, "
"col={} must be in [0, {}), [0, {})",
row, col, m_rows, m_cols));
}
return m_mean(row, col);
}
/**
* @brief Return the cached mean at a flat row-major index.
* @pre ready() is true, the cache is current, and index is valid; the index
* is not checked.
*/
PEDESTAL_TYPE mean_unchecked(ssize_t index) const { return m_mean[index]; }
/**
* @brief Calculate and return the population variance of every pixel.
* @throws std::runtime_error if ready() is false.
* @note The result is normalized by n_samples.
*/
NDArray<PEDESTAL_TYPE, 2> variance() {
if (!ready()) {
throw std::runtime_error(
"Pedestal is not ready, cannot return variance");
}
NDArray<PEDESTAL_TYPE, 2> res({m_rows, m_cols});
for (ssize_t i = 0; i < m_sum.size(); ++i) {
res[i] = variance_unchecked(i);
}
return res;
}
/**
* @brief Calculate the population variance at (row, col).
* @throws std::runtime_error if ready() is false or either index is out of
* range.
*/
PEDESTAL_TYPE variance(const uint32_t row, const uint32_t col) const {
if (!ready()) {
throw std::runtime_error(
"Pedestal is not ready, cannot return variance");
}
if (row >= m_rows || col >= m_cols) {
throw std::runtime_error(fmt::format(
"Invalid indices for FastPedestal variance: row={}, "
"col={} must be in [0, {}), [0, {})",
row, col, m_rows, m_cols));
}
return variance_unchecked(rc_to_index(row, col));
}
/**
* @brief Calculate the population variance at a flat row-major index.
* @pre ready() is true and index is valid; the index is not checked.
*/
PEDESTAL_TYPE variance_unchecked(ssize_t index) const {
const auto &entry = m_sum[index];
const auto m = entry.sum * m_inv_samples;
return std::fma(-m, m, entry.sum2 * m_inv_samples);
}
/**
* @brief Calculate and return the population standard deviation of every
* pixel.
* @throws std::runtime_error if ready() is false.
*/
NDArray<PEDESTAL_TYPE, 2> std() {
if (!ready()) {
throw std::runtime_error(
"Pedestal is not ready, cannot return std");
}
NDArray<PEDESTAL_TYPE, 2> res({m_rows, m_cols});
for (ssize_t i = 0; i < m_sum.size(); ++i) {
res[i] = std_unchecked(i);
}
return res;
}
/**
* @brief Calculate the population standard deviation at (row, col).
* @throws std::runtime_error if ready() is false or either index is out of
* range.
*/
PEDESTAL_TYPE std(const uint32_t row, const uint32_t col) const {
if (!ready()) {
throw std::runtime_error(
"Pedestal is not ready, cannot return std");
}
if (row >= m_rows || col >= m_cols) {
throw std::runtime_error(
fmt::format("Invalid indices for FastPedestal std: row={}, "
"col={} must be in [0, {}), [0, {})",
row, col, m_rows, m_cols));
}
return std::sqrt(variance(row, col));
}
/**
* @brief Calculate the population standard deviation at a flat row-major
* index.
* @pre ready() is true and index is valid; the index is not checked.
*/
PEDESTAL_TYPE std_unchecked(ssize_t index) const {
return std::sqrt(variance_unchecked(index));
}
/**
* @brief Return whether initialization is complete (cur_samples() equals
* n_samples()).
*/
bool ready() const { return m_ready; }
/**
* @brief Return the stored number of accumulated initialization frames.
* @note The value is in [0, n_samples] and does not change during
* steady-state pushes.
*/
uint32_t cur_samples() const { return m_cur_samples; }
/**
* @brief Zero the moments and cached mean, and mark the pedestal not ready.
*/
void clear() {
m_sum = Entry{0., 0.};
m_mean = PEDESTAL_TYPE(0.);
m_cur_samples = 0;
m_ready = false;
}
/**
* @brief Update every pixel using the steady-state exponential estimator,
* giving the new value weight 1 / n_samples.
* @param frame Frame whose shape must exactly match the pedestal.
* @throws std::runtime_error if the shape differs or ready() is false.
*/
template <typename T> void push_ema(NDView<T, 2> frame) {
if (frame.shape() != std::array<ssize_t, 2>{m_rows, m_cols}) {
throw std::runtime_error(
"Frame shape does not match pedestal shape");
}
if (!ready()) {
throw std::runtime_error("Pedestal is not ready, cannot push");
}
const auto size = static_cast<std::size_t>(m_rows) * m_cols;
const auto *data = frame.data();
for (std::size_t index = 0; index < size; ++index) {
push_ema_unchecked(index, data[index]);
}
}
/**
* @brief Update every pixel from a Frame using the steady-state estimator.
* @tparam T Actual pixel type stored in frame; this is not runtime-checked.
* @param frame Frame whose shape must exactly match the pedestal.
* @throws std::runtime_error if the shape differs or ready() is false.
*/
template <typename T> void push_ema(Frame &frame) {
push_ema<T>(frame.view<T>());
}
/**
* @brief Update the exponential moving average with smoothing factor
* 1/n_samples
* @param row Pixel row.
* @param col Pixel column.
* @param val New pixel value.
* @pre row and col are valid; indices are not checked.
* @throws std::runtime_error if ready() is false.
*/
template <typename T>
void push_ema(const uint32_t row, const uint32_t col, const T val) {
if (!ready()) {
throw std::runtime_error("Pedestal is not ready, cannot push");
}
push_ema_unchecked(rc_to_index(row, col), val);
}
/**
* @brief Update one pixel and its cached mean without runtime checks in
* release builds.
* @param index Flat row-major pixel index.
* @param value New pixel value, with weight 1 / n_samples.
* @pre ready() is true and index is a valid flat row-major index. These
* preconditions are asserted only in debug builds.
*/
template <typename T>
void push_ema_unchecked(const std::size_t index, const T value) noexcept {
assert(m_ready);
assert(index < static_cast<std::size_t>(m_sum.size()));
const auto val = static_cast<double>(value);
auto &entry = m_sum[index];
entry.sum += val - entry.sum * m_inv_samples;
entry.sum2 += val * val - entry.sum2 * m_inv_samples;
m_mean[index] = static_cast<PEDESTAL_TYPE>(entry.sum * m_inv_samples);
}
/**
* @brief Accumulate one initialization frame.
* @param frame Frame whose shape must exactly match the pedestal.
* @throws std::runtime_error if the shape differs or n_samples frames have
* already been accumulated.
* @note The statistics can be accessed and ready() becomes true only after
* the n_samples frames have been added.
*/
template <typename T> void add_init_frame(NDView<T, 2> frame) {
if (frame.shape() != std::array<ssize_t, 2>{m_rows, m_cols}) {
throw std::runtime_error(
"Frame shape does not match pedestal shape");
}
// if the pedestal is already initialized we cannot add more frames
if (ready()) {
throw std::runtime_error("Pedestal initialization is already done");
}
for (ssize_t i = 0; i < m_sum.size(); ++i) {
const auto val = static_cast<double>(frame[i]);
auto &entry = m_sum[i];
entry.sum += val;
entry.sum2 += val * val;
}
m_cur_samples += 1;
if (m_cur_samples == m_samples) {
update_mean();
m_ready = true;
}
}
/**
* @brief Initialize from n_samples frames after skip_first, then apply all
* remaining file frames as steady-state updates.
* @tparam T Pixel representation stored in the file; this is not checked.
* @param filename Input image file. Its dimensions define the pedestal
* shape.
* @param n_samples Number of initialization frames.
* @param skip_first Number of leading frames to ignore.
* @throws std::runtime_error if fewer than n_samples frames remain after
* skip_first, or if any constructor argument is invalid.
*/
template <typename T>
static FastPedestal from_file(const std::filesystem::path &filename,
uint32_t n_samples = 1000,
uint32_t skip_first = 0) {
File f(filename);
const auto total_frames = f.total_frames();
const auto first_frame = static_cast<size_t>(skip_first);
const auto initialization_frames = static_cast<size_t>(n_samples);
if (first_frame > total_frames ||
initialization_frames > total_frames - first_frame) {
throw std::runtime_error(
"File has less frames than the number of samples needed to "
"initialize the pedestal");
}
if (skip_first > 0) {
f.seek(static_cast<size_t>(skip_first));
}
const auto rows = static_cast<uint32_t>(f.rows());
const auto cols = static_cast<uint32_t>(f.cols());
FastPedestal pedestal(rows, cols, n_samples);
NDArray<T, 2> frame({rows, cols});
auto frame_index = first_frame;
const auto initialization_end = first_frame + initialization_frames;
while (frame_index < initialization_end) {
f.read_into(frame.buffer());
pedestal.template add_init_frame<T>(frame.view());
frame_index++;
}
// read the rest of the file
while (frame_index < total_frames) {
f.read_into(frame.buffer());
pedestal.template push_ema<T>(frame.view());
frame_index++;
}
return pedestal;
}
/** @brief Return the number of image rows. */
uint32_t rows() const { return m_rows; }
/** @brief Return the number of image columns. */
uint32_t cols() const { return m_cols; }
/**
* @brief Return the initialization frame count and steady-state
* update-weight denominator.
*/
uint32_t n_samples() const { return m_samples; }
private:
/**
* @brief Write the cached mean after the final add_init_frame. All other
* (non initialization) pushes update the cached mean immediately.
*/
void update_mean() {
for (ssize_t i = 0; i < m_sum.size(); i++) {
auto &entry = m_sum[i];
m_mean[i] = static_cast<PEDESTAL_TYPE>(entry.sum * m_inv_samples);
}
}
};
} // namespace aare
include/aare/File.hpp
+2 -1
View File
@@ -64,8 +64,9 @@ class File {
size_t total_frames() const;
size_t rows() const;
size_t cols() const;
Dtype dtype() const;
DetectorType detector_type() const;
};
} // namespace aare
} // namespace aare
include/aare/FileInterface.hpp
+9 -4
View File
@@ -145,10 +145,15 @@ class FileInterface {
*/
virtual size_t bitdepth() const = 0;
virtual DetectorType detector_type() const = 0;
/**
* @brief get the data type of the pixels
* @return pixel data type
*/
virtual Dtype dtype() const {
return Dtype::from_bitdepth(static_cast<uint8_t>(bitdepth()));
}
// function to query the data type of the file
/*virtual DataType dtype = 0; */
virtual DetectorType detector_type() const = 0;
virtual ~FileInterface() = default;
@@ -168,4 +173,4 @@ class FileInterface {
// size_t current_frame{};
};
} // namespace aare
} // namespace aare
include/aare/Fit.hpp
-100
View File
@@ -1,111 +1,11 @@
// SPDX-License-Identifier: MPL-2.0
#pragma once
#include <cmath>
#include <vector>
#include "aare/FitModel.hpp"
#include "aare/NDArray.hpp"
#include "aare/utils/par.hpp"
#include "aare/utils/task.hpp"
namespace aare {
namespace func {
double gaus(const double x, const double *par);
NDArray<double, 1> gaus(NDView<double, 1> x, NDView<double, 1> par);
double pol1(const double x, const double *par);
NDArray<double, 1> pol1(NDView<double, 1> x, NDView<double, 1> par);
double scurve(const double x, const double *par);
NDArray<double, 1> scurve(NDView<double, 1> x, NDView<double, 1> par);
double scurve2(const double x, const double *par);
NDArray<double, 1> scurve2(NDView<double, 1> x, NDView<double, 1> par);
} // namespace func
static constexpr int DEFAULT_NUM_THREADS = 4;
/**
* @brief Fit a 1D Gaussian to data.
* @param data data to fit
* @param x x values
*/
NDArray<double, 1> fit_gaus(NDView<double, 1> x, NDView<double, 1> y);
/**
* @brief Fit a 1D Gaussian to each pixel. Data layout [row, col, values]
* @param x x values
* @param y y values, layout [row, col, values]
* @param n_threads number of threads to use
*/
NDArray<double, 3> fit_gaus(NDView<double, 1> x, NDView<double, 3> y,
int n_threads = DEFAULT_NUM_THREADS);
/**
* @brief Fit a 1D Gaussian with error estimates
* @param x x values
* @param y y values, layout [row, col, values]
* @param y_err error in y, layout [row, col, values]
* @param par_out output parameters
* @param par_err_out output error parameters
*/
void fit_gaus(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
NDView<double, 1> par_out, NDView<double, 1> par_err_out,
double &chi2);
/**
* @brief Fit a 1D Gaussian to each pixel with error estimates. Data layout
* [row, col, values]
* @param x x values
* @param y y values, layout [row, col, values]
* @param y_err error in y, layout [row, col, values]
* @param par_out output parameters, layout [row, col, values]
* @param par_err_out output parameter errors, layout [row, col, values]
* @param n_threads number of threads to use
*/
void fit_gaus(NDView<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_err,
NDView<double, 3> par_out, NDView<double, 3> par_err_out,
NDView<double, 2> chi2_out, int n_threads = DEFAULT_NUM_THREADS);
NDArray<double, 1> fit_pol1(NDView<double, 1> x, NDView<double, 1> y);
NDArray<double, 3> fit_pol1(NDView<double, 1> x, NDView<double, 3> y,
int n_threads = DEFAULT_NUM_THREADS);
void fit_pol1(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
NDView<double, 1> par_out, NDView<double, 1> par_err_out,
double &chi2);
// TODO! not sure we need to offer the different version in C++
void fit_pol1(NDView<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_err,
NDView<double, 3> par_out, NDView<double, 3> par_err_out,
NDView<double, 2> chi2_out, int n_threads = DEFAULT_NUM_THREADS);
NDArray<double, 1> fit_scurve(NDView<double, 1> x, NDView<double, 1> y);
NDArray<double, 3> fit_scurve(NDView<double, 1> x, NDView<double, 3> y,
int n_threads);
void fit_scurve(NDView<double, 1> x, NDView<double, 1> y,
NDView<double, 1> y_err, NDView<double, 1> par_out,
NDView<double, 1> par_err_out, double &chi2);
void fit_scurve(NDView<double, 1> x, NDView<double, 3> y,
NDView<double, 3> y_err, NDView<double, 3> par_out,
NDView<double, 3> par_err_out, NDView<double, 2> chi2_out,
int n_threads);
NDArray<double, 1> fit_scurve2(NDView<double, 1> x, NDView<double, 1> y);
NDArray<double, 3> fit_scurve2(NDView<double, 1> x, NDView<double, 3> y,
int n_threads);
void fit_scurve2(NDView<double, 1> x, NDView<double, 1> y,
NDView<double, 1> y_err, NDView<double, 1> par_out,
NDView<double, 1> par_err_out, double &chi2);
void fit_scurve2(NDView<double, 1> x, NDView<double, 3> y,
NDView<double, 3> y_err, NDView<double, 3> par_out,
NDView<double, 3> par_err_out, NDView<double, 2> chi2_out,
int n_threads);
// ---------------------------------------------------------------------------
// Minuit2-based pixel fitting.
// Template bodies and explicit instantiations live in src/Fit.cpp.
include/aare/JungfrauDataFile.hpp
+2 -1
View File
@@ -49,6 +49,7 @@ class JungfrauDataFile : public FileInterface {
size_t pixels_per_frame() override;
size_t bytes_per_pixel() const;
size_t bitdepth() const override;
Dtype dtype() const override { return Dtype::UINT16; }
void seek(size_t frame_index)
override; //!< seek to the given frame index (note not byte offset)
size_t tell() override; //!< get the frame index of the file pointer
@@ -113,4 +114,4 @@ class JungfrauDataFile : public FileInterface {
std::filesystem::path fpath(size_t frame_index) const;
};
} // namespace aare
} // namespace aare
include/aare/MultiThreadedFileReader.hpp
+104
View File
@@ -0,0 +1,104 @@
// SPDX-License-Identifier: MPL-2.0
#pragma once
#include "aare/Dtype.hpp"
#include "aare/File.hpp"
#include <cstddef>
#include <filesystem>
#include <optional>
#include <vector>
namespace aare::experimental {
/**
* @brief Read independent chunks of a file in parallel.
*
* Each worker opens its own File instance, so seeking and reading do not share
* mutable file state. Chunks are written directly to their position in the
* destination buffer and the resulting frame order is the same as in the file.
*/
class MultiThreadedFileReader {
public:
/**
* @param fname path accepted by File
* @param n_threads maximum number of worker threads
* @param chunk_size number of frames claimed by a worker at a time
* @param total_frames number of frames to read, or all frames when omitted
*/
MultiThreadedFileReader(std::filesystem::path fname, size_t n_threads,
size_t chunk_size,
std::optional<size_t> total_frames = std::nullopt);
MultiThreadedFileReader(const MultiThreadedFileReader &) = delete;
MultiThreadedFileReader &
operator=(const MultiThreadedFileReader &) = delete;
MultiThreadedFileReader(MultiThreadedFileReader &&) noexcept = default;
MultiThreadedFileReader &
operator=(MultiThreadedFileReader &&) noexcept = default;
/**
* @brief Read one chunk per active worker into a caller-owned buffer.
*
* The buffer must hold at least next_read_bytes() bytes. The reader's
* position advances by the returned number of frames. At the end of the
* configured range this function returns zero and does not access the
* destination.
*/
size_t read_into(std::byte *destination);
/** @brief Read the next wave of chunks into an owned byte buffer. */
std::vector<std::byte> read();
/** @brief Read every frame remaining from the current position. */
std::vector<std::byte> read_all();
/** @brief Set the next frame index to read. The end position is valid. */
void seek(size_t frame_index);
/** @brief Return the next frame index to read. */
size_t tell() const noexcept { return m_current_frame; }
/** @brief Close all worker files. Safe to call more than once. */
void close() noexcept { m_files.clear(); }
/** @brief Return whether the worker files are open. */
bool is_open() const noexcept { return !m_files.empty(); }
size_t n_threads() const noexcept { return m_n_threads; }
size_t chunk_size() const noexcept { return m_chunk_size; }
size_t total_frames() const noexcept { return m_total_frames; }
size_t source_total_frames() const noexcept {
return m_source_total_frames;
}
size_t rows() const noexcept { return m_rows; }
size_t cols() const noexcept { return m_cols; }
size_t bitdepth() const noexcept { return m_bitdepth; }
Dtype dtype() const noexcept { return m_dtype; }
size_t bytes_per_frame() const noexcept { return m_bytes_per_frame; }
size_t total_bytes() const noexcept { return m_total_bytes; }
size_t remaining_frames() const noexcept;
size_t next_read_frames() const noexcept;
size_t next_read_bytes() const noexcept {
return next_read_frames() * m_bytes_per_frame;
}
private:
std::filesystem::path m_fname;
size_t m_n_threads;
size_t m_chunk_size;
size_t m_total_frames;
size_t m_source_total_frames;
size_t m_rows;
size_t m_cols;
size_t m_bitdepth;
Dtype m_dtype;
size_t m_bytes_per_frame;
size_t m_total_bytes;
size_t m_current_frame;
std::vector<File> m_files;
void ensure_open() const;
};
} // namespace aare::experimental
include/aare/NDArray.hpp
+24 -6
View File
@@ -146,6 +146,23 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
*/
~NDArray() { delete[] data_; }
/**
* @brief Copy data from a view of matching shape into this array without
* reallocating.
* @param v view to copy from, must have the same shape as this array
* @throws std::runtime_error if the shapes differ
*
* Use this instead of assigning a new NDArray when the buffer needs to be
* kept, for example when the array is part of a preallocated pool.
*/
void copy_from(const NDView<T, Ndim> v) {
if (v.shape() != shape_) {
throw std::runtime_error(LOCATION +
"Shape mismatch in NDArray::copy_from");
}
std::copy(v.begin(), v.end(), begin());
}
///////////////////////////////////////////////////////////////////////////////
// Iterators and indexing
//
@@ -506,12 +523,13 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
return *this;
}
/**
* @brief Create a view of the NDArray.
*
* @return NDView<T, Ndim>
*/
NDView<T, Ndim> view() const { return NDView<T, Ndim>{data_, shape_}; }
/** @brief Create a mutable view of the NDArray. */
NDView<T, Ndim> view() { return NDView<T, Ndim>{data_, shape_}; }
/** @brief Create a read-only view of a const NDArray. */
NDView<const T, Ndim> view() const {
return NDView<const T, Ndim>{data_, shape_};
}
private:
/**
include/aare/NDView.hpp
+10 -1
View File
@@ -12,6 +12,7 @@
#include <iostream>
#include <numeric>
#include <stdexcept>
#include <type_traits>
#include <vector>
namespace aare {
@@ -92,6 +93,12 @@ class NDView : public ArrayExpr<NDView<T, Ndim>, Ndim> {
: buffer_(buffer), strides_(c_strides<Ndim>(shape)), shape_(shape),
size_(num_elements(shape)) {}
template <typename U, std::enable_if_t<
std::is_convertible_v<U (*)[], T (*)[]>, int> = 0>
NDView(const NDView<U, Ndim> &other) noexcept
: buffer_(other.buffer_), strides_(other.strides_),
shape_(other.shape_), size_(other.size_) {}
template <typename... Ix>
std::enable_if_t<sizeof...(Ix) == Ndim, T &> operator()(Ix... index) {
return buffer_[element_offset(strides_, index...)];
@@ -216,6 +223,8 @@ class NDView : public ArrayExpr<NDView<T, Ndim>, Ndim> {
}
private:
template <typename, ssize_t> friend class NDView;
T *buffer_{nullptr};
std::array<ssize_t, Ndim> strides_{};
std::array<ssize_t, Ndim> shape_{};
@@ -263,4 +272,4 @@ template <typename T> NDView<T, 1> make_view(std::vector<T> &vec) {
return NDView<T, 1>(vec.data(), {static_cast<ssize_t>(vec.size())});
}
} // namespace aare
} // namespace aare
include/aare/NumpyFile.hpp
+2 -2
View File
@@ -61,7 +61,7 @@ class NumpyFile : public FileInterface {
* @brief get the data type of the numpy file
* @return DType
*/
Dtype dtype() const { return m_header.dtype; }
Dtype dtype() const override { return m_header.dtype; }
/**
* @brief get the shape of the numpy file
@@ -129,4 +129,4 @@ class NumpyFile : public FileInterface {
void write_impl(void *data, uint64_t size);
};
} // namespace aare
} // namespace aare
include/aare/Pedestal.hpp
+27 -3
View File
@@ -29,28 +29,39 @@ template <typename SUM_TYPE = double> class Pedestal {
// Relies on having more reads than pushes to the pedestal
NDArray<SUM_TYPE, 2> m_mean;
// Cache std. Only refreshed via update_std() to keep push() cheap.
NDArray<SUM_TYPE, 2> m_std;
public:
Pedestal(uint32_t rows, uint32_t cols, uint32_t n_samples = 1000)
: m_rows(rows), m_cols(cols), m_samples(n_samples),
m_cur_samples(NDArray<uint32_t, 2>({rows, cols}, 0)),
m_sum(NDArray<SUM_TYPE, 2>({rows, cols})),
m_sum2(NDArray<SUM_TYPE, 2>({rows, cols})),
m_mean(NDArray<SUM_TYPE, 2>({rows, cols})) {
m_mean(NDArray<SUM_TYPE, 2>({rows, cols})),
m_std(NDArray<SUM_TYPE, 2>({rows, cols})) {
assert(rows > 0 && cols > 0 && n_samples > 0);
m_sum = 0;
m_sum2 = 0;
m_mean = 0;
m_std = 0;
}
~Pedestal() = default;
NDArray<SUM_TYPE, 2> mean() { return m_mean; }
const NDView<SUM_TYPE, 2> view() const { return m_mean.view(); }
NDView<const SUM_TYPE, 2> view() const { return m_mean.view(); }
SUM_TYPE mean(const uint32_t row, const uint32_t col) const {
return m_mean(row, col);
}
NDArray<SUM_TYPE, 2> cached_std() { return m_std; }
SUM_TYPE cached_std(const uint32_t row, const uint32_t col) const {
return m_std(row, col);
}
SUM_TYPE std(const uint32_t row, const uint32_t col) const {
return std::sqrt(variance(row, col));
}
@@ -87,6 +98,7 @@ template <typename SUM_TYPE = double> class Pedestal {
m_sum2 = 0;
m_cur_samples = 0;
m_mean = 0;
m_std = 0;
}
void clear(const uint32_t row, const uint32_t col) {
@@ -94,6 +106,7 @@ template <typename SUM_TYPE = double> class Pedestal {
m_sum2(row, col) = 0;
m_cur_samples(row, col) = 0;
m_mean(row, col) = 0;
m_std(row, col) = 0;
}
template <typename T> void push(NDView<T, 2> frame) {
@@ -211,6 +224,17 @@ template <typename SUM_TYPE = double> class Pedestal {
*/
void update_mean() { m_mean = m_sum / m_cur_samples; }
/**
* @brief Refresh the cached std for all pixels from the current sums.
* Kept separate from push() so pushes stay cheap; call before reading
* cached_std() (analogous to update_mean()).
*/
void update_std() {
for (uint32_t i = 0; i < m_rows * m_cols; i++) {
m_std(i / m_cols, i % m_cols) = std(i / m_cols, i % m_cols);
}
}
template <typename T>
void push_fast(const uint32_t row, const uint32_t col, const T val_) {
// Assume we reached the steady state where all pixels have
@@ -221,4 +245,4 @@ template <typename SUM_TYPE = double> class Pedestal {
m_mean(row, col) = m_sum(row, col) / m_samples;
}
};
} // namespace aare
} // namespace aare
include/aare/RawFile.hpp
+2 -1
View File
@@ -118,6 +118,7 @@ class RawFile : public FileInterface {
*/
size_t cols(const size_t roi_index) const;
size_t bitdepth() const override;
Dtype dtype() const override { return Dtype::from_bitdepth(bitdepth()); }
size_t n_modules() const;
/**
@@ -170,4 +171,4 @@ class RawFile : public FileInterface {
void open_subfiles(const size_t roi_index);
};
} // namespace aare
} // namespace aare
include/aare/decode.hpp
+2 -2
View File
@@ -35,7 +35,7 @@ uint32_t mask32to24bits(uint32_t input, BitOffset offset = {});
* @param offset Offset within the first byte to where the data starts (0-7
* bits)
*/
void expand24to32bit(NDView<uint8_t, 1> input, NDView<uint32_t, 1> output,
void expand24to32bit(NDView<const uint8_t, 1> input, NDView<uint32_t, 1> output,
BitOffset offset = {});
/**
@@ -43,7 +43,7 @@ void expand24to32bit(NDView<uint8_t, 1> input, NDView<uint32_t, 1> output,
* @param input input buffer with 4 bit values packed into 8 bit
* @param output output buffer with 8 bit values
*/
void expand4to8bit(NDView<uint8_t, 1> input, NDView<uint8_t, 1> output);
void expand4to8bit(NDView<const uint8_t, 1> input, NDView<uint8_t, 1> output);
/**
* @brief Apply custom weights to a 16-bit input value. Will sum up
include/aare/defs.hpp
-81
View File
@@ -38,87 +38,6 @@ inline constexpr size_t bits_per_byte = 8;
void assert_failed(const std::string &msg);
class DynamicCluster {
public:
int cluster_sizeX;
int cluster_sizeY;
int16_t x;
int16_t y;
Dtype dt; // 4 bytes
private:
std::byte *m_data;
public:
DynamicCluster(int cluster_sizeX_, int cluster_sizeY_,
Dtype dt_ = Dtype(typeid(int32_t)))
: cluster_sizeX(cluster_sizeX_), cluster_sizeY(cluster_sizeY_),
dt(dt_) {
m_data = new std::byte[cluster_sizeX * cluster_sizeY * dt.bytes()]{};
}
DynamicCluster() : DynamicCluster(3, 3) {}
DynamicCluster(const DynamicCluster &other)
: DynamicCluster(other.cluster_sizeX, other.cluster_sizeY, other.dt) {
if (this == &other)
return;
x = other.x;
y = other.y;
memcpy(m_data, other.m_data, other.bytes());
}
DynamicCluster &operator=(const DynamicCluster &other) {
if (this == &other)
return *this;
this->~DynamicCluster();
new (this) DynamicCluster(other);
return *this;
}
DynamicCluster(DynamicCluster &&other) noexcept
: cluster_sizeX(other.cluster_sizeX),
cluster_sizeY(other.cluster_sizeY), x(other.x), y(other.y),
dt(other.dt), m_data(other.m_data) {
other.m_data = nullptr;
other.dt = Dtype(Dtype::TypeIndex::ERROR);
}
~DynamicCluster() { delete[] m_data; }
template <typename T> T get(int idx) {
(sizeof(T) == dt.bytes())
? 0
: throw std::invalid_argument("[ERROR] Type size mismatch");
return *reinterpret_cast<T *>(m_data + idx * dt.bytes());
}
template <typename T> auto set(int idx, T val) {
(sizeof(T) == dt.bytes())
? 0
: throw std::invalid_argument("[ERROR] Type size mismatch");
return memcpy(m_data + idx * dt.bytes(), &val, dt.bytes());
}
template <typename T> std::string to_string() const {
(sizeof(T) == dt.bytes())
? 0
: throw std::invalid_argument("[ERROR] Type size mismatch");
std::string s = "x: " + std::to_string(x) + " y: " + std::to_string(y) +
"\nm_data: [";
for (int i = 0; i < cluster_sizeX * cluster_sizeY; i++) {
s += std::to_string(
*reinterpret_cast<T *>(m_data + i * dt.bytes())) +
" ";
}
s += "]";
return s;
}
/**
* @brief size of the cluster in bytes when saved to a file
*/
size_t size() const { return cluster_sizeX * cluster_sizeY; }
size_t bytes() const { return cluster_sizeX * cluster_sizeY * dt.bytes(); }
auto begin() const { return m_data; }
auto end() const {
return m_data + cluster_sizeX * cluster_sizeY * dt.bytes();
}
std::byte *data() { return m_data; }
};
/**
* @brief header contained in parts of frames
*/
include/aare/hist/PedestalTrackingPixelHistogram.hpp
+23 -2
View File
@@ -1,4 +1,5 @@
#pragma once
#include "aare/FastPedestal.hpp"
#include "aare/NDArray.hpp"
#include "aare/NDView.hpp"
#include "aare/Pedestal.hpp"
@@ -44,7 +45,7 @@ class PedestalTrackingPixelHistogram {
// worker using the LOCAL row index (i.e. 0..row_count(t)-1), NOT the
// global row index. Owned exclusively by worker `t` during a
// dispatched fan-out.
std::vector<Pedestal<AxisType>> partial_pedestals_;
std::vector<FastPedestal<AxisType>> partial_pedestals_;
std::vector<NDArray<AxisType, 2>> partial_std_; // cached for pedestal
// tracking
@@ -65,6 +66,11 @@ class PedestalTrackingPixelHistogram {
// Always the outermost lock; work_mutex_ is taken briefly inside it.
mutable std::mutex fill_mutex_;
// Serialises producers of async_queue_ (which is SPSC) and gives
// fill_from_file exclusive ownership of the ingestion path while its
// direct, double-buffered dispatch is active.
std::mutex ingestion_mutex_;
// Async producer/consumer pipeline. SPSC queue feeds the coordinator
// thread, which batches queued images before dispatching them.
std::unique_ptr<AsyncQueue> async_queue_;
@@ -107,8 +113,23 @@ class PedestalTrackingPixelHistogram {
void fill_async(NDArray<FrameType, 2> &&image);
/**
* @brief Fill from ordered, parallel-read batches.
*
* After the first batch, reading batch N+1 overlaps histogram processing
* for batch N. Two fixed-capacity read buffers are allocated once and
* reused for the duration of the call.
*
* @param fname input file accepted by File
* @param max_frames maximum frames to process, or all frames for -1
* @param verbose print periodic progress information
* @param reader_threads number of MultiThreadedFileReader workers
* @param reader_chunk_size frames claimed per reader worker and batch
*/
void fill_from_file(const std::filesystem::path &fname,
ssize_t max_frames = -1, bool verbose = false);
ssize_t max_frames = -1, bool verbose = false,
std::size_t reader_threads = 2,
std::size_t reader_chunk_size = 4);
void process_pedestal_file(const std::filesystem::path &fname,
ssize_t max_frames = -1, bool verbose = false);
include/aare/hist/PixelHistogramImpl.hpp
+20 -9
View File
@@ -10,6 +10,7 @@ silently dropped.
#include "aare/NDArray.hpp"
#include "aare/NDView.hpp"
#include <algorithm>
#include <cstddef>
#include <limits>
#include <stdexcept>
@@ -37,12 +38,15 @@ template <typename T, typename StorageType> class PixelHistogramImpl {
void fill(const NDView<T, 2> &frame);
void fill(int row, int col, T value);
void fill_unchecked(int row, int col, T value);
// Fill using a row-major pixel index, avoiding repeated 3D stride
// calculation in tiled worker loops. The index is intentionally unchecked.
void fill_flat_unchecked(std::size_t pixel, T value);
NDArray<StorageType, 3> values() const;
// Zero-copy view of the underlying [rows x cols x n_bins] storage.
// Lifetime is tied to *this. Use for low-level merge/stitching paths;
// prefer values() for the public API where you want an owned copy.
NDView<StorageType, 3> view() const;
NDView<const StorageType, 3> view() const;
NDArray<T, 1> bin_centers() const;
NDArray<T, 1> bin_edges() const;
};
@@ -98,22 +102,29 @@ void PixelHistogramImpl<T, StorageType>::fill(int row, int col, T value) {
template <typename T, typename StorageType>
void PixelHistogramImpl<T, StorageType>::fill_unchecked(int row, int col,
T value) {
const auto pixel =
static_cast<std::size_t>(row) * static_cast<std::size_t>(m_cols) +
static_cast<std::size_t>(col);
fill_flat_unchecked(pixel, value);
}
template <typename T, typename StorageType>
void PixelHistogramImpl<T, StorageType>::fill_flat_unchecked(std::size_t pixel,
T value) {
if (value < m_xmin || value >= m_xmax) {
return;
}
int bin = static_cast<int>((value - m_xmin) * m_scale);
// Guard against floating-point rounding pushing val just below
// xmax to bin == n_bins.
if (bin >= m_n_bins) {
bin = m_n_bins - 1;
}
bin = std::clamp(bin, 0, m_n_bins - 1);
auto &cell = m_values.data()[pixel * static_cast<std::size_t>(m_n_bins) +
static_cast<std::size_t>(bin)];
if constexpr (std::is_integral_v<StorageType>) {
if (m_values(row, col, bin) >=
std::numeric_limits<StorageType>::max()) {
if (cell >= std::numeric_limits<StorageType>::max()) {
return;
}
}
++m_values(row, col, bin);
++cell;
}
template <typename T, typename StorageType>
@@ -122,7 +133,7 @@ NDArray<StorageType, 3> PixelHistogramImpl<T, StorageType>::values() const {
}
template <typename T, typename StorageType>
NDView<StorageType, 3> PixelHistogramImpl<T, StorageType>::view() const {
NDView<const StorageType, 3> PixelHistogramImpl<T, StorageType>::view() const {
return m_values.view();
}
include/aare/utils/math_helpers.hpp
+14
View File
@@ -0,0 +1,14 @@
#pragma once
#include <cstddef>
namespace aare {
/**
* @brief Compute the ceiling of the integer division of n by d.
* @param n The numerator.
* @param d The denominator.
* @return The ceiling of the integer division.
*/
constexpr size_t ceil_div(size_t n, size_t d) { return n / d + (n % d != 0); }
} // namespace aare
patches/lmfit.patch
-13
View File
@@ -1,13 +0,0 @@
diff --git a/lib/CMakeLists.txt b/lib/CMakeLists.txt
index 4efb7ed..6533660 100644
--- a/lib/CMakeLists.txt
+++ b/lib/CMakeLists.txt
@@ -11,7 +11,7 @@ target_compile_definitions(${lib} PRIVATE "LMFIT_EXPORT") # for Windows DLL expo
target_include_directories(${lib}
PUBLIC
- $<BUILD_INTERFACE:${CMAKE_SOURCE_DIR}/>
+ $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/>
$<INSTALL_INTERFACE:include/>
)
python/CMakeLists.txt
+3 -1
View File
@@ -68,8 +68,10 @@ set(PYTHON_FILES
aare/ClusterFinder.py
aare/ClusterVector.py
aare/Cluster.py
aare/FastPedestal.py
aare/calibration.py
aare/func.py
aare/factory.py
aare/experimental.py
aare/RawFile.py
aare/transform.py
aare/ScanParameters.py
python/aare/Cluster.py
+2 -2
View File
@@ -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)
return cls(x, y, data)
python/aare/ClusterFinder.py
+1 -12
View File
@@ -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.
python/aare/FastPedestal.py
+56
View File
@@ -0,0 +1,56 @@
# 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 an empty per-pixel running pedestal.
This factory hides the dtype suffix used by the templated C++ bindings.
Call ``add_init_frame()`` exactly ``n_samples`` times before using the
statistics or calling ``push_ema()``. Subsequent frames have weight
``1 / n_samples`` in the running mean and population variance.
Args:
rows: Number of image rows.
cols: Number of image columns.
n_samples: Initialization frame count and steady-state update-weight
denominator.
dtype: Output dtype for the mean, variance, and standard deviation.
Supported values 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.
After ignoring ``skip_first`` frames, the next ``n_samples`` frames
initialize the pedestal. Every remaining frame is then applied as a
steady-state update. Input frames are read as uint16 data.
Args:
filename: Input image file.
n_samples: Number of frames used for initialization.
skip_first: Number of leading frames to ignore.
dtype: Output dtype for the mean, variance, and standard deviation.
Raises:
RuntimeError: If fewer than ``n_samples`` frames remain after
``skip_first`` or ``n_samples`` is zero.
"""
cls = _get_fast_pedestal_class(dtype)
return cls.from_file(
filename, n_samples=n_samples, skip_first=skip_first
)
FastPedestal.from_file = from_file
python/aare/__init__.py
+20 -7
View File
@@ -20,8 +20,24 @@ except ImportError:
pass
from . import transform
from . import experimental
from ._aare import File, RawMasterFile, RawSubFile, JungfrauDataFile
from ._aare import (
FastPedestal_d,
FastPedestal_f,
FastPedestal_i16,
Pedestal_d,
Pedestal_f,
Pedestal_i16,
ClusterFinder_Cluster3x3i,
VarClusterFinder,
)
from ._aare import (
File,
JungfrauDataFile,
RawMasterFile,
RawSubFile,
)
from ._aare import Pedestal_d, Pedestal_f, ClusterFinder_Cluster3x3i, VarClusterFinder
from ._aare import DetectorType, ReadoutMode
from ._aare import hitmap
@@ -31,14 +47,14 @@ from ._aare import corner
# from ._aare import ClusterFinderMT, ClusterCollector, ClusterFileSink, ClusterVector_i
from ._version import __version__
from .ClusterFinder import ClusterFinder, ClusterFinderFrozen, ClusterCollector, ClusterFinderMT, ClusterFileSink, ClusterFile
from .ClusterFinder import ClusterFinderCUDA, ClusterFinderCUDAGraph, _cuda_available, find_cluster_views_batched_iter
from .FastPedestal import FastPedestal
from .ClusterFinder import (ClusterFinder, ClusterFinderFrozen, ClusterCollector, ClusterFinderMT, ClusterFileSink, ClusterFile,
ClusterFinderCUDA, ClusterFinderCUDAGraph, _cuda_available, find_cluster_views_batched_iter)
from .ClusterVector import ClusterVector
from .Cluster import Cluster
from ._aare import Gaussian, RisingScurve, FallingScurve, Pol1, Pol2, GaussianErfcPlateau, GaussianChargeSharing, GaussianChargeSharingKb
from ._aare import fit
from ._aare import fit_gaus, fit_pol1, fit_scurve, fit_scurve2
from ._aare import Interpolator
from ._aare import calculate_eta2, calculate_eta3, calculate_cross_eta3, calculate_full_eta2
from ._aare import reduce_to_2x2, reduce_to_3x3
@@ -54,9 +70,6 @@ from .ScanParameters import ScanParameters
from .utils import random_pixels, random_pixel, flat_list, add_colorbar, Timer
#make functions available in the top level API
from .func import *
from .calibration import *
from ._aare import apply_calibration, count_switching_pixels
from ._aare import calculate_pedestal, calculate_pedestal_float, calculate_pedestal_g0, calculate_pedestal_g0_float
python/aare/experimental.py
+6
View File
@@ -0,0 +1,6 @@
# SPDX-License-Identifier: MPL-2.0
"""Experimental APIs that may change without notice."""
from ._aare.experimental import MultiThreadedFileReader
__all__ = ["MultiThreadedFileReader"]
python/aare/factory.py
+36
View File
@@ -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 ``<name>_<dtype suffix>``."""
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
python/aare/func.py
-2
View File
@@ -1,2 +0,0 @@
# SPDX-License-Identifier: MPL-2.0
from ._aare import gaus, pol1, scurve, scurve2
python/examples/fits.py
+42 -70
View File
@@ -1,106 +1,78 @@
# SPDX-License-Identifier: MPL-2.0
import matplotlib.pyplot as plt
import numpy as np
import sys
sys.path.insert(0, '/home/kferjaoui/sw/aare/build')
from aare import fit_gaus, fit_pol1
from aare import Gaussian, fit
from aare import pol1
textpm = f"±" #
textmu = f"μ" #
textsigma = f"σ" #
from aare import Gaussian, Pol1
textpm = "±"
textmu = "μ"
textsigma = "σ"
# ================================= Gauss fit =================================
# Parameters
mu = np.random.uniform(1, 100) # Mean of Gaussian
sigma = np.random.uniform(4, 20) # Standard deviation
num_points = 10000 # Number of points for smooth distribution
# noise_sigma = 10
# Generate Gaussian distribution
data = np.random.normal(mu, sigma, num_points)
mu = np.random.uniform(1, 100)
sigma = np.random.uniform(4, 20)
data = np.random.normal(mu, sigma, 10000)
counts, edges = np.histogram(data, bins=100)
x = 0.5 * (edges[:-1] + edges[1:]) # proper bin centers
x = 0.5 * (edges[:-1] + edges[1:])
y = counts.astype(np.float64)
# Poisson noise
yerr = np.sqrt(np.maximum(y, 1))
# yerr = np.abs(np.random.normal(0, noise_sigma, len(x)))
# Create subplot
fig0, ax0 = plt.subplots(1, 1, num=0, figsize=(12, 8))
# Add the errors as error bars in the step plot
ax0.errorbar(x, y, yerr=yerr, fmt=". ", capsize=5)
ax0.grid()
# Fit with lmfit
result_lm = fit_gaus(x, y, yerr)
par_lm = result_lm["par"]
err_lm = result_lm["par_err"]
chi2_lm = result_lm["chi2"]
print("[lmfit] fit_gaus: ", par_lm, err_lm, chi2_lm)
gaussian = Gaussian(compute_errors=True)
result = gaussian.fit(x, y, yerr)
par = result["par"]
err = result["par_err"]
chi2 = result["chi2"]
print(f"Gaussian.fit: par={par}, err={err}, chi2={chi2}")
# Fit with Minuit2 + analytic gradient + Hesse errors
gaussian = Gaussian()
gaussian.compute_errors = True
result_m2 = gaussian.fit(x, y, yerr)
par_m2 = result_m2['par']
err_m2 = result_m2['par_err']
chi2_m2 = result_m2['chi2']
print(f"[minuit2] gaussian.fit: par={par_m2}, err={err_m2}, chi2={chi2_m2}")
x = np.linspace(x[0], x[-1], 1000)
ax0.plot(x, gaussian(x, par_lm), marker="", label="fit_gaus")
ax0.plot(x, gaussian(x, par_m2), marker="", linestyle=":", label="fit_gaus_minuit_grad")
x_plot = np.linspace(x[0], x[-1], 1000)
ax0.plot(x_plot, gaussian(x_plot, par), marker="", label="Gaussian.fit")
ax0.legend()
ax0.set(xlabel="x", ylabel="Counts",
ax0.set(
xlabel="x",
ylabel="Counts",
title=(
f"fit_gaus: A={par_lm[0]:0.2f}{textpm}{err_lm[0]:0.2f} "
f"{textmu}={par_lm[1]:0.2f}{textpm}{err_lm[1]:0.2f} "
f"{textsigma}={par_lm[2]:0.2f}{textpm}{err_lm[2]:0.2f}\n"
f"minuit_grad: A={par_m2[0]:0.2f}{textpm}{err_m2[0]:0.2f} "
f"{textmu}={par_m2[1]:0.2f}{textpm}{err_m2[1]:0.2f} "
f"{textsigma}={par_m2[2]:0.2f}{textpm}{err_m2[2]:0.2f}\n"
f"A={par[0]:0.2f}{textpm}{err[0]:0.2f} "
f"{textmu}={par[1]:0.2f}{textpm}{err[1]:0.2f} "
f"{textsigma}={par[2]:0.2f}{textpm}{err[2]:0.2f}\n"
f"(truth: {textmu}={mu:0.2f}, {textsigma}={sigma:0.2f})"
),
)
fig0.tight_layout()
# ================================= pol1 fit =================================
# Parameters
# ================================= Pol1 fit =================================
n_points = 40
# Generate random slope and intercept (origin)
slope = np.random.uniform(-10, 10) # Random slope between 0.5 and 2.0
intercept = np.random.uniform(-10, 10) # Random intercept between -10 and 10
# Generate random x values
slope = np.random.uniform(-10, 10)
intercept = np.random.uniform(-10, 10)
x_values = np.random.uniform(-10, 10, n_points)
# Calculate y values based on the linear function y = mx + b + error
errors = np.abs(np.random.normal(0, np.random.uniform(1, 5), n_points))
var_points = np.random.normal(0, np.random.uniform(0.1, 2), n_points)
y_values = slope * x_values + intercept + var_points
y_values = slope * x_values + intercept + np.random.normal(0, 1, n_points)
fig1, ax1 = plt.subplots(1, 1, num=1, figsize=(12, 8))
ax1.errorbar(x_values, y_values, yerr=errors, fmt=". ", capsize=5)
result_pol = fit_pol1(x_values, y_values, errors)
par = result_pol["par"]
err = result_pol["par_err"]
x = np.linspace(np.min(x_values), np.max(x_values), 1000)
ax1.plot(x, pol1(x, par), marker="")
ax1.set(xlabel="x", ylabel="y", title=f"a = {par[0]:0.2f}{textpm}{err[0]:0.2f}\n"
f"b = {par[1]:0.2f}{textpm}{err[1]:0.2f}\n"
f"(init: {slope:0.2f}, {intercept:0.2f})")
pol1 = Pol1(compute_errors=True)
result = pol1.fit(x_values, y_values, errors)
par = result["par"]
err = result["par_err"]
x_plot = np.linspace(np.min(x_values), np.max(x_values), 1000)
ax1.plot(x_plot, pol1(x_plot, par), marker="")
ax1.set(
xlabel="x",
ylabel="y",
title=(
f"intercept = {par[0]:0.2f}{textpm}{err[0]:0.2f}\n"
f"slope = {par[1]:0.2f}{textpm}{err[1]:0.2f}\n"
f"(truth: {intercept:0.2f}, {slope:0.2f})"
),
)
fig1.tight_layout()
plt.show()
python/src/bind_Cluster.hpp
-1
View File
@@ -10,7 +10,6 @@
#include <pybind11/stl_bind.h>
namespace py = pybind11;
using pd_type = double;
using namespace aare;
python/src/bind_ClusterCollector.hpp
+3 -2
View File
@@ -6,6 +6,8 @@
#include "aare/ClusterVector.hpp"
#include "aare/NDView.hpp"
#include "aare/Pedestal.hpp"
#include "module_config.hpp"
#include "np_helper.hpp"
#include <cstdint>
@@ -15,7 +17,6 @@
#include <pybind11/stl_bind.h>
namespace py = pybind11;
using pd_type = double;
using namespace aare;
@@ -30,7 +31,7 @@ void define_ClusterCollector(py::module &m, const std::string &typestr) {
using ClusterType = Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>;
py::class_<ClusterCollector<ClusterType>>(m, class_name.c_str())
.def(py::init<ClusterFinderMT<ClusterType, uint16_t, double> *>())
.def(py::init<ClusterFinderMT<ClusterType, uint16_t, pd_type> *>())
.def("stop", &ClusterCollector<ClusterType>::stop)
.def(
"steal_clusters",
python/src/bind_ClusterFile.hpp
+3 -1
View File
@@ -47,6 +47,8 @@ void define_ClusterFile(py::module &m, const std::string &typestr) {
})
.def("set_roi", &ClusterFile<ClusterType>::set_roi, py::arg("roi"))
.def("tell", &ClusterFile<ClusterType>::tell)
.def("estimate_n_clusters",
&ClusterFile<ClusterType>::estimate_n_clusters)
.def(
"set_noise_map",
[](ClusterFile<ClusterType> &self, py::array_t<int32_t> noise_map) {
@@ -82,4 +84,4 @@ void define_ClusterFile(py::module &m, const std::string &typestr) {
});
}
#pragma GCC diagnostic pop
#pragma GCC diagnostic pop
python/src/bind_ClusterFileSink.hpp
+2 -2
View File
@@ -15,8 +15,6 @@
#include <pybind11/stl_bind.h>
namespace py = pybind11;
using pd_type = double;
using namespace aare;
#pragma GCC diagnostic push
@@ -29,6 +27,8 @@ void define_ClusterFileSink(py::module &m, const std::string &typestr) {
using ClusterType = Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>;
// TODO! adapt to set pedestal type (needs templating of ClusterFileSink)
// or maybe access through base class?
py::class_<ClusterFileSink<ClusterType>>(m, class_name.c_str())
.def(py::init<ClusterFinderMT<ClusterType, uint16_t, double> *,
const std::filesystem::path &>())
python/src/bind_ClusterFinder.hpp
+4 -1
View File
@@ -6,6 +6,8 @@
#include "aare/ClusterVector.hpp"
#include "aare/NDView.hpp"
#include "aare/Pedestal.hpp"
#include "module_config.hpp"
#include "np_helper.hpp"
#include <cstdint>
@@ -15,7 +17,6 @@
#include <pybind11/stl_bind.h>
namespace py = pybind11;
using pd_type = double;
using namespace aare;
@@ -48,6 +49,8 @@ void define_ClusterFinder(py::module &m, const std::string &typestr) {
})
.def("clear_pedestal",
&ClusterFinder<ClusterType, uint16_t, pd_type>::clear_pedestal)
.def("update_threshold",
&ClusterFinder<ClusterType, uint16_t, pd_type>::update_threshold)
.def_property_readonly(
"pedestal",
[](ClusterFinder<ClusterType, uint16_t, pd_type> &self) {
python/src/bind_ClusterFinderMT.hpp
+7 -4
View File
@@ -6,6 +6,8 @@
#include "aare/ClusterVector.hpp"
#include "aare/NDView.hpp"
#include "aare/Pedestal.hpp"
#include "module_config.hpp"
#include "np_helper.hpp"
#include <cstdint>
@@ -15,7 +17,6 @@
#include <pybind11/stl_bind.h>
namespace py = pybind11;
using pd_type = double;
using namespace aare;
@@ -31,9 +32,10 @@ void define_ClusterFinderMT(py::module &m, const std::string &typestr) {
py::class_<ClusterFinderMT<ClusterType, uint16_t, pd_type>>(
m, class_name.c_str())
.def(py::init<Shape<2>, pd_type, size_t, size_t>(),
.def(py::init<Shape<2>, pd_type, size_t, size_t, size_t>(),
py::arg("image_size"), py::arg("n_sigma") = 5.0,
py::arg("capacity") = 2048, py::arg("n_threads") = 3)
py::arg("capacity") = 2048, py::arg("n_threads") = 3,
py::arg("queue_depth") = 16)
.def("push_pedestal_frame",
[](ClusterFinderMT<ClusterType, uint16_t, pd_type> &self,
py::array_t<uint16_t> frame) {
@@ -46,7 +48,6 @@ void define_ClusterFinderMT(py::module &m, const std::string &typestr) {
py::array_t<uint16_t> frame, uint64_t frame_number) {
auto view = make_view_2d(frame);
self.find_clusters(view, frame_number);
return;
},
py::arg(), py::arg("frame_number") = 0)
.def_property_readonly(
@@ -56,6 +57,8 @@ void define_ClusterFinderMT(py::module &m, const std::string &typestr) {
})
.def("clear_pedestal",
&ClusterFinderMT<ClusterType, uint16_t, pd_type>::clear_pedestal)
.def("update_threshold",
&ClusterFinderMT<ClusterType, uint16_t, pd_type>::update_threshold)
.def("sync", &ClusterFinderMT<ClusterType, uint16_t, pd_type>::sync)
.def("stop", &ClusterFinderMT<ClusterType, uint16_t, pd_type>::stop)
.def("start", &ClusterFinderMT<ClusterType, uint16_t, pd_type>::start)
python/src/bind_ClusterVector.hpp
+2 -2
View File
@@ -15,7 +15,6 @@
#include <pybind11/stl_bind.h>
namespace py = pybind11;
using pd_type = double;
using namespace aare;
@@ -72,6 +71,7 @@ void define_ClusterVector(py::module &m, const std::string &typestr) {
R"(calculates sum of 2x2 subcluster with highest energy and index relative to cluster center 0: top_left, 1: top_right, 2: bottom_left, 3: bottom_right
)")
.def_property_readonly("size", &ClusterVector<ClusterType>::size)
.def("empty", &ClusterVector<ClusterType>::empty)
.def("item_size", &ClusterVector<ClusterType>::item_size)
.def_property_readonly("fmt",
[typestr](ClusterVector<ClusterType> &self) {
@@ -169,4 +169,4 @@ void define_3x3_reduction(py::module &m) {
py::arg("clustervector"));
}
#pragma GCC diagnostic pop
#pragma GCC diagnostic pop
python/src/bind_FastPedestal.hpp
+147
View File
@@ -0,0 +1,147 @@
// SPDX-License-Identifier: MPL-2.0
#include "aare/FastPedestal.hpp"
#include "np_helper.hpp"
#include <cstdint>
#include <pybind11/numpy.h>
#include <pybind11/pybind11.h>
#include <pybind11/stl.h>
namespace py = pybind11;
template <typename SUM_TYPE>
void define_fast_pedestal_bindings(py::module &m, const std::string &name) {
py::class_<FastPedestal<SUM_TYPE>>(
m, name.c_str(),
"Maintain a per-pixel running mean, population variance, and "
"standard deviation.",
py::buffer_protocol())
.def(py::init<uint32_t, uint32_t, uint32_t>(), py::arg("rows"),
py::arg("cols"), py::arg("n_samples"),
"Construct an empty pedestal. It becomes ready after n_samples "
"calls to add_init_frame().")
.def(py::init<uint32_t, uint32_t>(), py::arg("rows"), py::arg("cols"),
"Construct an empty pedestal with n_samples=1000.")
.def(
"mean",
[](FastPedestal<SUM_TYPE> &self) {
auto mean = new NDArray<SUM_TYPE, 2>{};
*mean = self.mean();
return return_image_data(mean);
},
"Return a copy of the cached mean. The pedestal must be ready.")
.def(
"var",
[](FastPedestal<SUM_TYPE> &self) {
auto variance = new NDArray<SUM_TYPE, 2>{};
*variance = self.variance();
return return_image_data(variance);
},
"Return the population variance, normalized by n_samples, as a "
"NumPy array. The pedestal must be ready.")
.def(
"std",
[](FastPedestal<SUM_TYPE> &self) {
auto standard_deviation = new NDArray<SUM_TYPE, 2>{};
*standard_deviation = self.std();
return return_image_data(standard_deviation);
},
"Return the population standard deviation as a NumPy array. The "
"pedestal must be ready.")
.def(
"view",
[](py::object self_py) {
return py::module_::import("numpy").attr("asarray")(self_py);
},
"Return a non-owning, non-writable NumPy view of the cached mean. "
"The pedestal must be ready.")
// We need to buffer protocol to allow for numpy operations using the
// pedestal mean
.def_buffer([](FastPedestal<SUM_TYPE> &self) {
auto mean = self.view();
return py::buffer_info(
const_cast<SUM_TYPE *>(mean.data()), sizeof(SUM_TYPE),
py::format_descriptor<SUM_TYPE>::format(), 2,
{static_cast<py::ssize_t>(mean.shape(0)),
static_cast<py::ssize_t>(mean.shape(1))},
{static_cast<py::ssize_t>(mean.strides()[0] * sizeof(SUM_TYPE)),
static_cast<py::ssize_t>(mean.strides()[1] *
sizeof(SUM_TYPE))},
true);
})
// Subtracting a FastPedestal from a NumPy array
.def(
"__array_ufunc__",
[](py::object self, py::object ufunc, const std::string &method,
py::args inputs, py::kwargs kwargs) -> py::object {
if (method != "__call__" || inputs.size() != 2 ||
inputs[1].ptr() != self.ptr() ||
py::cast<std::string>(ufunc.attr("__name__")) !=
"subtract") {
return py::reinterpret_borrow<py::object>(
Py_NotImplemented);
}
auto mean =
py::module_::import("builtins").attr("memoryview")(self);
return ufunc(inputs[0], mean, **kwargs);
},
"Support subtracting a FastPedestal from a NumPy array.")
.def("clear", py::overload_cast<>(&FastPedestal<SUM_TYPE>::clear),
"Reset all statistics and initialization state to zero.")
.def_property_readonly("rows", &FastPedestal<SUM_TYPE>::rows,
"Number of image rows.")
.def_property_readonly("cols", &FastPedestal<SUM_TYPE>::cols,
"Number of image columns.")
.def_property_readonly("cur_samples",
&FastPedestal<SUM_TYPE>::cur_samples,
"Number of initialization frames accumulated. "
"Steady-state pushes do not change it.")
.def_property_readonly(
"ready", &FastPedestal<SUM_TYPE>::ready,
"Whether n_samples initialization frames have been accumulated.")
.def_property_readonly(
"n_samples", &FastPedestal<SUM_TYPE>::n_samples,
"Initialization frame count and steady-state update-weight "
"denominator.")
.def(
"clone",
[](FastPedestal<SUM_TYPE> &pedestal) {
return FastPedestal<SUM_TYPE>(pedestal);
},
"Return an independent copy of the pedestal and its state.")
.def(
"push_ema",
[](FastPedestal<SUM_TYPE> &pedestal,
py::array_t<uint16_t, py::array::c_style> &frame) {
pedestal.push_ema(make_view_2d(frame));
},
py::arg("frame").noconvert(),
"Update exponential moving average. The pedstal must "
"already be ready for this update.")
.def(
"add_init_frame",
[](FastPedestal<SUM_TYPE> &pedestal,
py::array_t<uint16_t, py::array::c_style> &frame) {
pedestal.add_init_frame(make_view_2d(frame));
},
py::arg("frame").noconvert(),
"Accumulate one uint16 initialization frame. Call exactly "
"n_samples times to make the pedestal ready.")
.def_static(
"from_file",
[](const std::filesystem::path &filename, uint32_t n_samples,
uint32_t skip_first) {
return FastPedestal<SUM_TYPE>::template from_file<uint16_t>(
filename, n_samples, skip_first);
},
py::arg("filename"), py::arg("n_samples") = 1000,
py::arg("skip_first") = 0,
"Create a pedestal from a uint16 file. Skip skip_first frames, use "
"the next n_samples for initialization, then apply every remaining "
"frame as a steady-state update.");
}
python/src/bind_MultiThreadedFileReader.hpp
+175
View File
@@ -0,0 +1,175 @@
// SPDX-License-Identifier: MPL-2.0
#pragma once
#include "aare/MultiThreadedFileReader.hpp"
#include <cstddef>
#include <cstdint>
#include <filesystem>
#include <optional>
#include <stdexcept>
#include <vector>
#include <pybind11/numpy.h>
#include <pybind11/pybind11.h>
#include <pybind11/stl.h>
#include <pybind11/stl/filesystem.h>
namespace py = pybind11;
inline py::dtype multi_threaded_reader_numpy_dtype(const aare::Dtype &dtype) {
using aare::Dtype;
if (dtype == Dtype::INT8)
return py::dtype::of<int8_t>();
if (dtype == Dtype::UINT8)
return py::dtype::of<uint8_t>();
if (dtype == Dtype::INT16)
return py::dtype::of<int16_t>();
if (dtype == Dtype::UINT16)
return py::dtype::of<uint16_t>();
if (dtype == Dtype::INT32)
return py::dtype::of<int32_t>();
if (dtype == Dtype::UINT32)
return py::dtype::of<uint32_t>();
if (dtype == Dtype::INT64)
return py::dtype::of<int64_t>();
if (dtype == Dtype::UINT64)
return py::dtype::of<uint64_t>();
if (dtype == Dtype::FLOAT)
return py::dtype::of<float>();
if (dtype == Dtype::DOUBLE)
return py::dtype::of<double>();
throw std::runtime_error("Unsupported pixel data type");
}
inline py::array
multi_threaded_reader_read(aare::experimental::MultiThreadedFileReader &reader,
bool read_all) {
const size_t n_frames =
read_all ? reader.remaining_frames() : reader.next_read_frames();
const std::vector<py::ssize_t> shape{
static_cast<py::ssize_t>(n_frames),
static_cast<py::ssize_t>(reader.rows()),
static_cast<py::ssize_t>(reader.cols())};
py::array image(multi_threaded_reader_numpy_dtype(reader.dtype()), shape);
auto *destination = reinterpret_cast<std::byte *>(image.mutable_data());
{
py::gil_scoped_release release;
if (read_all) {
size_t offset = 0;
while (reader.remaining_frames() != 0) {
const size_t frames_read =
reader.read_into(destination + offset);
offset += frames_read * reader.bytes_per_frame();
}
} else {
reader.read_into(destination);
}
}
return image;
}
inline void define_multi_threaded_file_reader_bindings(py::module_ &m) {
using aare::experimental::MultiThreadedFileReader;
auto reader =
py::class_<MultiThreadedFileReader>(m, "MultiThreadedFileReader");
reader.attr("__module__") = "aare.experimental";
reader
.def(py::init<std::filesystem::path, size_t, size_t,
std::optional<size_t>>(),
py::arg("fname"), py::arg("n_threads"), py::arg("chunk_size"),
py::arg("total_frames") = py::none(),
R"doc(
Read chunks of detector frames concurrently.
Each worker opens an independent File. The returned array is
ordered by frame index even though chunks are read in parallel.
Args:
fname: Path accepted by File.
n_threads: Maximum number of worker threads.
chunk_size: Number of frames read per claimed chunk.
total_frames: Optional frame limit. None reads all frames.
)doc")
.def(
"read",
[](MultiThreadedFileReader &self) {
return multi_threaded_reader_read(self, false);
},
R"doc(
Read one chunk per active worker into a NumPy array.
Returns:
An array containing at most n_threads * chunk_size frames.
An empty array is returned at the end of the configured
frame range. The GIL is released while file data is read.
)doc")
.def(
"read_all",
[](MultiThreadedFileReader &self) {
return multi_threaded_reader_read(self, true);
},
R"doc(
Read all frames remaining from the current position.
)doc")
.def_property_readonly("n_threads", &MultiThreadedFileReader::n_threads)
.def_property_readonly("chunk_size",
&MultiThreadedFileReader::chunk_size)
.def_property_readonly("total_frames",
&MultiThreadedFileReader::total_frames)
.def_property_readonly("source_total_frames",
&MultiThreadedFileReader::source_total_frames)
.def_property_readonly("rows", &MultiThreadedFileReader::rows)
.def_property_readonly("cols", &MultiThreadedFileReader::cols)
.def_property_readonly("bitdepth", &MultiThreadedFileReader::bitdepth)
.def_property_readonly("dtype",
[](const MultiThreadedFileReader &self) {
return multi_threaded_reader_numpy_dtype(
self.dtype());
})
.def_property_readonly("bytes_per_frame",
&MultiThreadedFileReader::bytes_per_frame)
.def_property_readonly("total_bytes",
&MultiThreadedFileReader::total_bytes)
.def_property_readonly("remaining_frames",
&MultiThreadedFileReader::remaining_frames)
.def_property_readonly("next_read_frames",
&MultiThreadedFileReader::next_read_frames)
.def_property_readonly("next_read_bytes",
&MultiThreadedFileReader::next_read_bytes)
.def("seek", &MultiThreadedFileReader::seek, py::arg("frame_index"))
.def("tell", &MultiThreadedFileReader::tell)
.def("close", &MultiThreadedFileReader::close,
"Close all worker files. Safe to call more than once.")
.def_property_readonly(
"closed",
[](const MultiThreadedFileReader &self) { return !self.is_open(); })
.def("__len__", &MultiThreadedFileReader::total_frames)
.def(
"__enter__",
[](MultiThreadedFileReader &self) -> MultiThreadedFileReader * {
if (!self.is_open()) {
throw std::runtime_error(
"Cannot enter a closed MultiThreadedFileReader");
}
return &self;
},
py::return_value_policy::reference_internal)
.def("__exit__",
[](MultiThreadedFileReader &self, const py::object &,
const py::object &, const py::object &) {
self.close();
return false;
})
.def(
"__iter__", [](MultiThreadedFileReader &self) { return &self; },
py::return_value_policy::reference_internal)
.def("__next__", [](MultiThreadedFileReader &self) {
if (self.remaining_frames() == 0) {
throw py::stop_iteration();
}
return multi_threaded_reader_read(self, false);
});
}
python/src/bind_PedestalTrackingPixelHistogram.hpp
+5 -1
View File
@@ -153,9 +153,13 @@ void define_pedestal_tracking_pixel_histogram_bindings(py::module &m) {
Args:
file_path: Path to the file to fill from
max_frames: Maximum number of frames to fill from the file (default: -1)
reader_threads: Number of parallel file reader workers (default: 2)
reader_chunk_size: Frames claimed by each reader worker per batch (default: 4)
)",
py::call_guard<py::gil_scoped_release>(), py::arg("fname"),
py::arg("max_frames") = -1, py::arg("verbose") = false)
py::arg("max_frames") = -1, py::arg("verbose") = false,
py::arg("reader_threads") = std::size_t{2},
py::arg("reader_chunk_size") = std::size_t{4})
.def("process_pedestal_file",
&PedestalTrackingPixelHistogram::process_pedestal_file,
R"(
python/src/ctb_raw_file.hpp
-1
View File
@@ -9,7 +9,6 @@
#include "aare/decode.hpp"
#include "aare/defs.hpp"
// #include "aare/fClusterFileV2.hpp"
#include "np_helper.hpp"
python/src/file.hpp
+1 -50
View File
@@ -7,7 +7,6 @@
#include "aare/RawSubFile.hpp"
#include "aare/defs.hpp"
// #include "aare/fClusterFileV2.hpp"
#include <cstdint>
#include <filesystem>
@@ -195,52 +194,4 @@ void define_file_io_bindings(py::module &m) {
});
#pragma GCC diagnostic pop
// py::class_<ClusterHeader>(m, "ClusterHeader")
// .def(py::init<>())
// .def_readwrite("frame_number", &ClusterHeader::frame_number)
// .def_readwrite("n_clusters", &ClusterHeader::n_clusters)
// .def("__repr__", [](const ClusterHeader &a) { return "<ClusterHeader:
// " + a.to_string() + ">"; });
// py::class_<ClusterV2_>(m, "ClusterV2_")
// .def(py::init<>())
// .def_readwrite("x", &ClusterV2_::x)
// .def_readwrite("y", &ClusterV2_::y)
// .def_readwrite("data", &ClusterV2_::data)
// .def("__repr__", [](const ClusterV2_ &a) { return "<ClusterV2_: " +
// a.to_string(false) + ">"; });
// py::class_<ClusterV2>(m, "ClusterV2")
// .def(py::init<>())
// .def_readwrite("cluster", &ClusterV2::cluster)
// .def_readwrite("frame_number", &ClusterV2::frame_number)
// .def("__repr__", [](const ClusterV2 &a) { return "<ClusterV2: " +
// a.to_string() + ">"; });
// py::class_<ClusterFileV2>(m, "ClusterFileV2")
// .def(py::init<const std::filesystem::path &, const std::string &>())
// .def("read", py::overload_cast<>(&ClusterFileV2::read))
// .def("read", py::overload_cast<int>(&ClusterFileV2::read))
// .def("frame_number", &ClusterFileV2::frame_number)
// .def("write", py::overload_cast<std::vector<ClusterV2> const
// &>(&ClusterFileV2::write))
// .def("close", &ClusterFileV2::close);
// m.def("to_clustV2", [](std::vector<DynamicCluster> &clusters, const int
// frame_number) {
// std::vector<ClusterV2> clusters_;
// for (auto &c : clusters) {
// ClusterV2 cluster;
// cluster.cluster.x = c.x;
// cluster.cluster.y = c.y;
// int i=0;
// for(auto &d : cluster.cluster.data) {
// d=c.get<double>(i++);
// }
// cluster.frame_number = frame_number;
// clusters_.push_back(cluster);
// }
// return clusters_;
// });
}
}
python/src/fit.hpp
+1 -447
View File
@@ -227,452 +227,6 @@ fit_dispatch(const aare::FitModel<Model> &model,
}
void define_fit_bindings(py::module &m) {
// TODO! Evaluate without converting to double
m.def(
"gaus",
[](py::array_t<double, py::array::c_style | py::array::forcecast> x,
py::array_t<double, py::array::c_style | py::array::forcecast> par) {
auto x_view = make_view_1d(x);
auto par_view = make_view_1d(par);
auto y = new NDArray<double, 1>{aare::func::gaus(x_view, par_view)};
return return_image_data(y);
},
R"(
Evaluate a 1D Gaussian function for all points in x using parameters par.
Parameters
----------
x : array_like
The points at which to evaluate the Gaussian function.
par : array_like
The parameters of the Gaussian function. The first element is the amplitude, the second element is the mean, and the third element is the standard deviation.
)",
py::arg("x"), py::arg("par"));
m.def(
"pol1",
[](py::array_t<double, py::array::c_style | py::array::forcecast> x,
py::array_t<double, py::array::c_style | py::array::forcecast> par) {
auto x_view = make_view_1d(x);
auto par_view = make_view_1d(par);
auto y = new NDArray<double, 1>{aare::func::pol1(x_view, par_view)};
return return_image_data(y);
},
R"(
Evaluate a 1D polynomial function for all points in x using parameters par. (p0+p1*x)
Parameters
----------
x : array_like
The points at which to evaluate the polynomial function.
par : array_like
The parameters of the polynomial function. The first element is the intercept, and the second element is the slope.
)",
py::arg("x"), py::arg("par"));
m.def(
"scurve",
[](py::array_t<double, py::array::c_style | py::array::forcecast> x,
py::array_t<double, py::array::c_style | py::array::forcecast> par) {
auto x_view = make_view_1d(x);
auto par_view = make_view_1d(par);
auto y =
new NDArray<double, 1>{aare::func::scurve(x_view, par_view)};
return return_image_data(y);
},
R"(
Evaluate a 1D scurve function for all points in x using parameters par.
Parameters
----------
x : array_like
The points at which to evaluate the scurve function.
par : array_like
The parameters of the scurve function. The first element is the background slope, the second element is the background intercept, the third element is the mean, the fourth element is the standard deviation, the fifth element is inflexion point count number, and the sixth element is C.
)",
py::arg("x"), py::arg("par"));
m.def(
"scurve2",
[](py::array_t<double, py::array::c_style | py::array::forcecast> x,
py::array_t<double, py::array::c_style | py::array::forcecast> par) {
auto x_view = make_view_1d(x);
auto par_view = make_view_1d(par);
auto y =
new NDArray<double, 1>{aare::func::scurve2(x_view, par_view)};
return return_image_data(y);
},
R"(
Evaluate a 1D scurve2 function for all points in x using parameters par.
Parameters
----------
x : array_like
The points at which to evaluate the scurve function.
par : array_like
The parameters of the scurve2 function. The first element is the background slope, the second element is the background intercept, the third element is the mean, the fourth element is the standard deviation, the fifth element is inflexion point count number, and the sixth element is C.
)",
py::arg("x"), py::arg("par"));
m.def(
"fit_gaus",
[](py::array_t<double, py::array::c_style | py::array::forcecast> x,
py::array_t<double, py::array::c_style | py::array::forcecast> y,
int n_threads) {
if (y.ndim() == 3) {
auto par = new NDArray<double, 3>{};
auto y_view = make_view_3d(y);
auto x_view = make_view_1d(x);
*par = aare::fit_gaus(x_view, y_view, n_threads);
return return_image_data(par);
} else if (y.ndim() == 1) {
auto par = new NDArray<double, 1>{};
auto y_view = make_view_1d(y);
auto x_view = make_view_1d(x);
*par = aare::fit_gaus(x_view, y_view);
return return_image_data(par);
} else {
throw std::runtime_error("Data must be 1D or 3D");
}
},
R"(
Fit a 1D Gaussian to data.
Parameters
----------
x : array_like
The x values.
y : array_like
The y values.
n_threads : int, optional
The number of threads to use. Default is 4.
)",
py::arg("x"), py::arg("y"), py::arg("n_threads") = 4);
m.def(
"fit_gaus",
[](py::array_t<double, py::array::c_style | py::array::forcecast> x,
py::array_t<double, py::array::c_style | py::array::forcecast> y,
py::array_t<double, py::array::c_style | py::array::forcecast> y_err,
int n_threads) {
if (y.ndim() == 3) {
// Allocate memory for the output
// Need to have pointers to allow python to manage
// the memory
auto par = new NDArray<double, 3>({y.shape(0), y.shape(1), 3});
auto par_err =
new NDArray<double, 3>({y.shape(0), y.shape(1), 3});
auto chi2 = new NDArray<double, 2>({y.shape(0), y.shape(1)});
// Make views of the numpy arrays
auto y_view = make_view_3d(y);
auto y_view_err = make_view_3d(y_err);
auto x_view = make_view_1d(x);
aare::fit_gaus(x_view, y_view, y_view_err, par->view(),
par_err->view(), chi2->view(), n_threads);
return py::dict("par"_a = return_image_data(par),
"par_err"_a = return_image_data(par_err),
"chi2"_a = return_image_data(chi2),
"Ndf"_a = y.shape(2) - 3);
} else if (y.ndim() == 1) {
// Allocate memory for the output
// Need to have pointers to allow python to manage
// the memory
auto par = new NDArray<double, 1>({3});
auto par_err = new NDArray<double, 1>({3});
// Decode the numpy arrays
auto y_view = make_view_1d(y);
auto y_view_err = make_view_1d(y_err);
auto x_view = make_view_1d(x);
double chi2 = 0;
aare::fit_gaus(x_view, y_view, y_view_err, par->view(),
par_err->view(), chi2);
return py::dict("par"_a = return_image_data(par),
"par_err"_a = return_image_data(par_err),
"chi2"_a = chi2, "Ndf"_a = y.size() - 3);
} else {
throw std::runtime_error("Data must be 1D or 3D");
}
},
R"(
Fit a 1D Gaussian to data with error estimates.
Parameters
----------
x : array_like
The x values.
y : array_like
The y values.
y_err : array_like
The error in the y values.
n_threads : int, optional
The number of threads to use. Default is 4.
)",
py::arg("x"), py::arg("y"), py::arg("y_err"), py::arg("n_threads") = 4);
m.def(
"fit_pol1",
[](py::array_t<double, py::array::c_style | py::array::forcecast> x,
py::array_t<double, py::array::c_style | py::array::forcecast> y,
int n_threads) {
if (y.ndim() == 3) {
auto par = new NDArray<double, 3>{};
auto x_view = make_view_1d(x);
auto y_view = make_view_3d(y);
*par = aare::fit_pol1(x_view, y_view, n_threads);
return return_image_data(par);
} else if (y.ndim() == 1) {
auto par = new NDArray<double, 1>{};
auto x_view = make_view_1d(x);
auto y_view = make_view_1d(y);
*par = aare::fit_pol1(x_view, y_view);
return return_image_data(par);
} else {
throw std::runtime_error("Data must be 1D or 3D");
}
},
py::arg("x"), py::arg("y"), py::arg("n_threads") = 4);
m.def(
"fit_pol1",
[](py::array_t<double, py::array::c_style | py::array::forcecast> x,
py::array_t<double, py::array::c_style | py::array::forcecast> y,
py::array_t<double, py::array::c_style | py::array::forcecast> y_err,
int n_threads) {
if (y.ndim() == 3) {
auto par = new NDArray<double, 3>({y.shape(0), y.shape(1), 2});
auto par_err =
new NDArray<double, 3>({y.shape(0), y.shape(1), 2});
auto y_view = make_view_3d(y);
auto y_view_err = make_view_3d(y_err);
auto x_view = make_view_1d(x);
auto chi2 = new NDArray<double, 2>({y.shape(0), y.shape(1)});
aare::fit_pol1(x_view, y_view, y_view_err, par->view(),
par_err->view(), chi2->view(), n_threads);
return py::dict("par"_a = return_image_data(par),
"par_err"_a = return_image_data(par_err),
"chi2"_a = return_image_data(chi2),
"Ndf"_a = y.shape(2) - 2);
} else if (y.ndim() == 1) {
auto par = new NDArray<double, 1>({2});
auto par_err = new NDArray<double, 1>({2});
auto y_view = make_view_1d(y);
auto y_view_err = make_view_1d(y_err);
auto x_view = make_view_1d(x);
double chi2 = 0;
aare::fit_pol1(x_view, y_view, y_view_err, par->view(),
par_err->view(), chi2);
return py::dict("par"_a = return_image_data(par),
"par_err"_a = return_image_data(par_err),
"chi2"_a = chi2, "Ndf"_a = y.size() - 2);
} else {
throw std::runtime_error("Data must be 1D or 3D");
}
},
R"(
Fit a 1D polynomial to data with error estimates.
Parameters
----------
x : array_like
The x values.
y : array_like
The y values.
y_err : array_like
The error in the y values.
n_threads : int, optional
The number of threads to use. Default is 4.
)",
py::arg("x"), py::arg("y"), py::arg("y_err"), py::arg("n_threads") = 4);
//=========
m.def(
"fit_scurve",
[](py::array_t<double, py::array::c_style | py::array::forcecast> x,
py::array_t<double, py::array::c_style | py::array::forcecast> y,
int n_threads) {
if (y.ndim() == 3) {
auto par = new NDArray<double, 3>{};
auto x_view = make_view_1d(x);
auto y_view = make_view_3d(y);
*par = aare::fit_scurve(x_view, y_view, n_threads);
return return_image_data(par);
} else if (y.ndim() == 1) {
auto par = new NDArray<double, 1>{};
auto x_view = make_view_1d(x);
auto y_view = make_view_1d(y);
*par = aare::fit_scurve(x_view, y_view);
return return_image_data(par);
} else {
throw std::runtime_error("Data must be 1D or 3D");
}
},
py::arg("x"), py::arg("y"), py::arg("n_threads") = 4);
m.def(
"fit_scurve",
[](py::array_t<double, py::array::c_style | py::array::forcecast> x,
py::array_t<double, py::array::c_style | py::array::forcecast> y,
py::array_t<double, py::array::c_style | py::array::forcecast> y_err,
int n_threads) {
if (y.ndim() == 3) {
auto par = new NDArray<double, 3>({y.shape(0), y.shape(1), 6});
auto par_err =
new NDArray<double, 3>({y.shape(0), y.shape(1), 6});
auto y_view = make_view_3d(y);
auto y_view_err = make_view_3d(y_err);
auto x_view = make_view_1d(x);
auto chi2 = new NDArray<double, 2>({y.shape(0), y.shape(1)});
aare::fit_scurve(x_view, y_view, y_view_err, par->view(),
par_err->view(), chi2->view(), n_threads);
return py::dict("par"_a = return_image_data(par),
"par_err"_a = return_image_data(par_err),
"chi2"_a = return_image_data(chi2),
"Ndf"_a = y.shape(2) - 2);
} else if (y.ndim() == 1) {
auto par = new NDArray<double, 1>({2});
auto par_err = new NDArray<double, 1>({2});
auto y_view = make_view_1d(y);
auto y_view_err = make_view_1d(y_err);
auto x_view = make_view_1d(x);
double chi2 = 0;
aare::fit_scurve(x_view, y_view, y_view_err, par->view(),
par_err->view(), chi2);
return py::dict("par"_a = return_image_data(par),
"par_err"_a = return_image_data(par_err),
"chi2"_a = chi2, "Ndf"_a = y.size() - 2);
} else {
throw std::runtime_error("Data must be 1D or 3D");
}
},
R"(
Fit a 1D polynomial to data with error estimates.
Parameters
----------
x : array_like
The x values.
y : array_like
The y values.
y_err : array_like
The error in the y values.
n_threads : int, optional
The number of threads to use. Default is 4.
)",
py::arg("x"), py::arg("y"), py::arg("y_err"), py::arg("n_threads") = 4);
m.def(
"fit_scurve2",
[](py::array_t<double, py::array::c_style | py::array::forcecast> x,
py::array_t<double, py::array::c_style | py::array::forcecast> y,
int n_threads) {
if (y.ndim() == 3) {
auto par = new NDArray<double, 3>{};
auto x_view = make_view_1d(x);
auto y_view = make_view_3d(y);
*par = aare::fit_scurve2(x_view, y_view, n_threads);
return return_image_data(par);
} else if (y.ndim() == 1) {
auto par = new NDArray<double, 1>{};
auto x_view = make_view_1d(x);
auto y_view = make_view_1d(y);
*par = aare::fit_scurve2(x_view, y_view);
return return_image_data(par);
} else {
throw std::runtime_error("Data must be 1D or 3D");
}
},
py::arg("x"), py::arg("y"), py::arg("n_threads") = 4);
m.def(
"fit_scurve2",
[](py::array_t<double, py::array::c_style | py::array::forcecast> x,
py::array_t<double, py::array::c_style | py::array::forcecast> y,
py::array_t<double, py::array::c_style | py::array::forcecast> y_err,
int n_threads) {
if (y.ndim() == 3) {
auto par = new NDArray<double, 3>({y.shape(0), y.shape(1), 6});
auto par_err =
new NDArray<double, 3>({y.shape(0), y.shape(1), 6});
auto y_view = make_view_3d(y);
auto y_view_err = make_view_3d(y_err);
auto x_view = make_view_1d(x);
auto chi2 = new NDArray<double, 2>({y.shape(0), y.shape(1)});
aare::fit_scurve2(x_view, y_view, y_view_err, par->view(),
par_err->view(), chi2->view(), n_threads);
return py::dict("par"_a = return_image_data(par),
"par_err"_a = return_image_data(par_err),
"chi2"_a = return_image_data(chi2),
"Ndf"_a = y.shape(2) - 2);
} else if (y.ndim() == 1) {
auto par = new NDArray<double, 1>({6});
auto par_err = new NDArray<double, 1>({6});
auto y_view = make_view_1d(y);
auto y_view_err = make_view_1d(y_err);
auto x_view = make_view_1d(x);
double chi2 = 0;
aare::fit_scurve2(x_view, y_view, y_view_err, par->view(),
par_err->view(), chi2);
return py::dict("par"_a = return_image_data(par),
"par_err"_a = return_image_data(par_err),
"chi2"_a = chi2, "Ndf"_a = y.size() - 2);
} else {
throw std::runtime_error("Data must be 1D or 3D");
}
},
R"(
Fit a 1D polynomial to data with error estimates.
Parameters
----------
x : array_like
The x values.
y : array_like
The y values.
y_err : array_like
The error in the y values.
n_threads : int, optional
The number of threads to use. Default is 4.
)",
py::arg("x"), py::arg("y"), py::arg("y_err"), py::arg("n_threads") = 4);
// ── Bind model classes ──────────────────────────────────────────
bind_fit_model<aare::model::Gaussian>(m, "Gaussian");
bind_fit_model<aare::model::GaussianErfcPlateau>(m, "GaussianErfcPlateau");
@@ -794,4 +348,4 @@ void define_fit_bindings(py::module &m) {
)",
py::arg("model"), py::arg("x"), py::arg("y"),
py::arg("y_err") = py::none(), py::arg("n_threads") = 4);
}
}
python/src/module.cpp
+14
View File
@@ -1,6 +1,8 @@
// SPDX-License-Identifier: MPL-2.0
// Files with bindings to the different classes
#include "module_config.hpp"
// New style file naming
#include "bind_Cluster.hpp"
#include "bind_ClusterCollector.hpp"
@@ -12,7 +14,9 @@
#include "bind_ClusterVector.hpp"
#include "bind_Defs.hpp"
#include "bind_Eta.hpp"
#include "bind_FastPedestal.hpp"
#include "bind_Interpolator.hpp"
#include "bind_MultiThreadedFileReader.hpp"
#include "bind_PedestalTrackingPixelHistogram.hpp"
#include "bind_PixelHistogram.hpp"
#include "bind_PixelMap.hpp"
@@ -30,6 +34,7 @@
#include "var_cluster.hpp"
// Pybind stuff
#include <cstdint>
#include <pybind11/pybind11.h>
#include <pybind11/stl.h>
@@ -61,7 +66,11 @@ double, 'f' for float)
define_ClusterCollector<T, N, M, U>(m, "Cluster" #N "x" #M #TYPE_CODE);
PYBIND11_MODULE(_aare, m) {
auto experimental = m.def_submodule(
"experimental", "Experimental APIs that may change without notice");
define_file_io_bindings(m);
define_multi_threaded_file_reader_bindings(experimental);
define_raw_file_io_bindings(m);
define_raw_sub_file_io_bindings(m);
define_ctb_raw_file_io_bindings(m);
@@ -72,6 +81,10 @@ PYBIND11_MODULE(_aare, m) {
define_pedestal_tracking_pixel_histogram_bindings(m);
define_pedestal_bindings<double>(m, "Pedestal_d");
define_pedestal_bindings<float>(m, "Pedestal_f");
define_pedestal_bindings<int16_t>(m, "Pedestal_i16");
define_fast_pedestal_bindings<double>(m, "FastPedestal_d");
define_fast_pedestal_bindings<float>(m, "FastPedestal_f");
define_fast_pedestal_bindings<int16_t>(m, "FastPedestal_i16");
define_fit_bindings(m);
define_interpolation_bindings(m);
define_jungfrau_data_file_io_bindings(m);
@@ -103,6 +116,7 @@ PYBIND11_MODULE(_aare, m) {
DEFINE_BINDINGS_CLUSTERFINDER(int, 3, 3, uint16_t, i);
DEFINE_BINDINGS_CLUSTERFINDER(double, 3, 3, uint16_t, d);
DEFINE_BINDINGS_CLUSTERFINDER(float, 3, 3, uint16_t, f);
DEFINE_BINDINGS_CLUSTERFINDER(int16_t, 3, 3, uint16_t, i16);
DEFINE_BINDINGS_CLUSTERFINDER(int, 5, 5, uint16_t, i);
DEFINE_BINDINGS_CLUSTERFINDER(double, 5, 5, uint16_t, d);
python/src/module_config.hpp
+5
View File
@@ -0,0 +1,5 @@
#pragma once
#include <cstdint>
// Configure module wide pedestal type for cluster finding
using pd_type = double;
python/src/pedestal.hpp
+6
View File
@@ -52,6 +52,12 @@ void define_pedestal_bindings(py::module &m, const std::string &name) {
*std = self.std();
return return_image_data(std);
})
.def("cached_std",
[](Pedestal<SUM_TYPE> &self) {
auto standard_deviation = new NDArray<SUM_TYPE, 2>{};
*standard_deviation = self.cached_std();
return return_image_data(standard_deviation);
})
.def(
"__array_ufunc__",
[](py::object self, py::object ufunc, const std::string &method,
python/src/raw_master_file.hpp
-1
View File
@@ -8,7 +8,6 @@
#include "aare/RawSubFile.hpp"
#include "aare/defs.hpp"
// #include "aare/fClusterFileV2.hpp"
#include <cstdint>
#include <filesystem>
python/src/raw_sub_file.hpp
-1
View File
@@ -7,7 +7,6 @@
#include "aare/RawSubFile.hpp"
#include "aare/defs.hpp"
// #include "aare/fClusterFileV2.hpp"
#include <cstdint>
#include <filesystem>
python/src/var_cluster.hpp
-6
View File
@@ -1,17 +1,11 @@
// SPDX-License-Identifier: MPL-2.0
#include "aare/VarClusterFinder.hpp"
#include "np_helper.hpp"
// #include "aare/defs.hpp"
// #include "aare/fClusterFileV2.hpp"
#include <cstdint>
// #include <filesystem>
#include <pybind11/numpy.h>
// #include <pybind11/iostream.h>
#include <pybind11/pybind11.h>
#include <pybind11/stl.h>
// #include <pybind11/stl/filesystem.h>
// #include <string>
namespace py = pybind11;
using namespace ::aare;
python/tests/test_Cluster.py
+10 -7
View File
@@ -103,15 +103,18 @@ def test_max_sum():
def test_cluster_finder():
"""Test ClusterFinder"""
shape = [100,100]
cf = _aare.ClusterFinder_Cluster3x3i(shape)
clusterfinder = _aare.ClusterFinder_Cluster3x3i([100,100])
#Push 1000 frames to the pedestal
for i in range(1000):
frame = np.random.normal(loc = 100, scale = 5, size = shape).astype(np.uint16)
cf.push_pedestal_frame(frame)
cf.update_threshold()
frame = np.zeros(shape=shape, dtype=np.uint16)
cf.find_clusters(frame)
#frame = np.random.rand(100,100)
frame = np.zeros(shape=[100,100])
clusterfinder.find_clusters(frame)
clusters = clusterfinder.steal_clusters(False) #conversion does not work
clusters = cf.steal_clusters(False) #conversion does not work
assert clusters.size == 0
python/tests/test_ClusterFile.py
+3 -1
View File
@@ -14,6 +14,8 @@ from conftest import test_data_path
def test_cluster_file(test_data_path):
"""Test ClusterFile"""
f = ClusterFile(test_data_path / "clust/single_frame_97_clustrers.clust")
assert f.estimate_n_clusters() == 97
assert f.tell() == 0
cv = f.read_clusters(10) #conversion does not work
@@ -62,4 +64,4 @@ def test_read_clusters_and_fill_histogram(test_data_path):
hist_py = pickle.load(f)
#Compare the two histograms
assert hist_aare == hist_py
assert hist_aare == hist_py
python/tests/test_ClusterVector.py
+3 -1
View File
@@ -16,6 +16,7 @@ def test_create_cluster_vector():
assert cv.cluster_size_x == 3
assert cv.cluster_size_y == 3
assert cv.size == 0
assert cv.empty()
def test_push_back_on_cluster_vector():
@@ -27,6 +28,7 @@ def test_push_back_on_cluster_vector():
cluster = _aare.Cluster2x2i(19, 22, np.ones(4, dtype=np.int32))
cv.push_back(cluster)
assert cv.size == 1
assert not cv.empty()
arr = np.array(cv, copy=False)
assert arr[0]['x'] == 19
@@ -127,4 +129,4 @@ def test_masking():
assert cv_masked_array[0]["x"] == 1
assert cv_masked_array[0]["y"] == 2
assert (cv_masked_array[0]["data"] == np.ones((3,3),dtype=np.int32)).all()
assert (cv_masked_array[0]["data"] == np.ones((3,3),dtype=np.int32)).all()
python/tests/test_FastPedestal.py
+145
View File
@@ -0,0 +1,145 @@
import numpy as np
import pytest
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)
def test_fast_pedestal_from_file_rejects_skip_beyond_end(tmp_path):
filename = tmp_path / "frames.npy"
np.save(filename, np.zeros((1, 1, 1), dtype=np.uint16))
with pytest.raises(RuntimeError, match="less frames"):
FastPedestal.from_file(filename, n_samples=1, skip_first=2)
@pytest.mark.parametrize(
("pedestal_type", "expected_dtype"),
[(FastPedestal_d, np.float64), (FastPedestal_f, np.float32)],
)
def test_fast_pedestal_initialization(pedestal_type, expected_dtype):
pedestal = pedestal_type(2, 3, 2)
first = np.array([[2, 4, 6], [8, 10, 12]], dtype=np.uint16)
second = np.array([[4, 6, 8], [10, 12, 14]], dtype=np.uint16)
pedestal.add_init_frame(first)
pedestal.add_init_frame(second)
expected_mean = np.array(
[[3, 5, 7], [9, 11, 13]], dtype=expected_dtype
)
np.testing.assert_array_equal(pedestal.mean(), expected_mean)
np.testing.assert_array_equal(pedestal.std(), np.ones((2, 3)))
def test_fast_pedestal_steady_state_push_ema():
pedestal = FastPedestal_d(1, 2, 2)
pedestal.add_init_frame(np.array([[2, 4]], dtype=np.uint16))
pedestal.add_init_frame(np.array([[4, 6]], dtype=np.uint16))
pedestal.push_ema(np.array([[6, 8]], dtype=np.uint16))
np.testing.assert_array_equal(pedestal.mean(), [[4.5, 6.5]])
def test_fast_pedestal_exposes_read_only_buffer_and_subtraction():
pedestal = FastPedestal_d(1, 2, 1)
pedestal.add_init_frame(np.array([[2, 4]], dtype=np.uint16))
view = np.asarray(pedestal)
result = np.array([[12, 14]], dtype=np.uint16) - pedestal
np.testing.assert_array_equal(view, [[2, 4]])
np.testing.assert_array_equal(result, [[10, 10]])
assert np.shares_memory(view, pedestal.view())
assert not view.flags.writeable
def test_fast_pedestal_rejects_wrong_shape():
pedestal = FastPedestal_d(2, 3)
with pytest.raises(RuntimeError, match="shape"):
pedestal.add_init_frame(np.zeros((2, 2), dtype=np.uint16))
python/tests/test_Fit.py
+15
View File
@@ -0,0 +1,15 @@
# SPDX-License-Identifier: MPL-2.0
import numpy as np
import aare
def test_gaussian_model_evaluates_and_fits_data():
x = np.linspace(-5.0, 5.0, 51)
expected = np.array([20.0, 0.5, 1.2])
model = aare.Gaussian()
y = model(x, expected)
result = model.fit(x, y)
np.testing.assert_allclose(result["par"], expected, atol=2e-3)
python/tests/test_Minuit2_gauss.ipynb
+408 -483
View File
File diff suppressed because one or more lines are too long
python/tests/test_Minuit2_scurve_1d.ipynb
+24 -136
View File
File diff suppressed because one or more lines are too long
python/tests/test_Minuit2_scurve_3d.ipynb
+22 -121
View File
File diff suppressed because one or more lines are too long
python/tests/test_MultiThreadedFileReader.py
+147
View File
@@ -0,0 +1,147 @@
# SPDX-License-Identifier: MPL-2.0
import numpy as np
import pytest
import aare
from aare.experimental import MultiThreadedFileReader
@pytest.fixture
def frame_file(tmp_path):
data = np.arange(10 * 2 * 3, dtype=np.uint16).reshape(10, 2, 3)
path = tmp_path / "frames.npy"
np.save(path, data)
return path, data
def test_experimental_import_path():
assert aare.experimental.MultiThreadedFileReader is MultiThreadedFileReader
assert MultiThreadedFileReader.__module__ == "aare.experimental"
assert not hasattr(aare, "MultiThreadedFileReader")
def test_reads_all_frames_in_order(frame_file):
path, expected = frame_file
reader = MultiThreadedFileReader(path, n_threads=2, chunk_size=3)
first = reader.read()
second = reader.read()
exhausted = reader.read()
assert np.array_equal(first, expected[:6])
assert np.array_equal(second, expected[6:])
assert exhausted.shape == (0, 2, 3)
assert first.dtype == np.uint16
assert reader.n_threads == 2
assert reader.chunk_size == 3
assert reader.total_frames == 10
assert reader.source_total_frames == 10
assert reader.rows == 2
assert reader.cols == 3
assert reader.bitdepth == 16
assert reader.dtype == np.dtype(np.uint16)
assert reader.bytes_per_frame == 12
assert reader.total_bytes == expected.nbytes
assert len(reader) == 10
assert reader.tell() == 10
assert reader.remaining_frames == 0
assert reader.next_read_frames == 0
assert reader.next_read_bytes == 0
def test_total_frame_limit(frame_file):
path, expected = frame_file
reader = MultiThreadedFileReader(
path, n_threads=2, chunk_size=2, total_frames=7
)
assert np.array_equal(reader.read(), expected[:4])
assert np.array_equal(reader.read_all(), expected[4:7])
def test_iteration_yields_one_chunk_per_thread(frame_file):
path, expected = frame_file
reader = MultiThreadedFileReader(path, n_threads=2, chunk_size=2)
batches = list(reader)
assert [len(batch) for batch in batches] == [4, 4, 2]
assert np.array_equal(np.concatenate(batches), expected)
def test_seek_resets_iteration_position(frame_file):
path, expected = frame_file
reader = MultiThreadedFileReader(path, n_threads=2, chunk_size=2)
reader.read()
assert reader.tell() == 4
reader.seek(1)
assert reader.tell() == 1
assert np.array_equal(reader.read(), expected[1:5])
with pytest.raises(IndexError):
reader.seek(11)
def test_context_manager_closes_worker_files(frame_file):
path, expected = frame_file
with MultiThreadedFileReader(path, n_threads=2, chunk_size=2) as reader:
assert not reader.closed
assert np.array_equal(reader.read(), expected[:4])
assert reader.closed
reader.close()
with pytest.raises(RuntimeError):
reader.read()
with pytest.raises(RuntimeError):
reader.seek(0)
with pytest.raises(RuntimeError):
with reader:
pass
def test_explicit_zero_frame_limit(frame_file):
path, expected = frame_file
reader = MultiThreadedFileReader(
path, n_threads=8, chunk_size=3, total_frames=0
)
actual = reader.read()
assert actual.shape == (0, *expected.shape[1:])
assert actual.dtype == expected.dtype
@pytest.mark.parametrize(
"dtype", [np.int8, np.int32, np.uint64, np.float32, np.float64]
)
def test_preserves_numpy_dtype(tmp_path, dtype):
expected = np.arange(4 * 2 * 3, dtype=dtype).reshape(4, 2, 3)
path = tmp_path / "typed-frames.npy"
np.save(path, expected)
reader = MultiThreadedFileReader(path, n_threads=2, chunk_size=3)
actual = reader.read()
assert actual.dtype == expected.dtype
assert reader.dtype == expected.dtype
assert np.array_equal(actual, expected)
@pytest.mark.parametrize(
("n_threads", "chunk_size", "total_frames"),
[(0, 1, None), (1, 0, None), (1, 1, 11)],
)
def test_invalid_configuration(
frame_file, n_threads, chunk_size, total_frames
):
path, _ = frame_file
with pytest.raises(ValueError):
MultiThreadedFileReader(
path,
n_threads=n_threads,
chunk_size=chunk_size,
total_frames=total_frames,
)
python/tests/test_PixelHistogram.py
+101
View File
@@ -0,0 +1,101 @@
# SPDX-License-Identifier: MPL-2.0
import numpy as np
import pytest
from aare import PixelHistogram
def _random_frames(rows, cols, n, xmin, xmax, seed=0):
rng = np.random.default_rng(seed)
return [rng.uniform(xmin - 0.25, xmax + 0.25, size=(rows, cols)).astype(np.float64)
for _ in range(n)]
def _reference_hdata(frames, rows, cols, n_bins, xmin, xmax):
expected = np.zeros((rows, cols, n_bins), dtype=np.uint16)
inv_range = n_bins / (xmax - xmin)
for img in frames:
for r in range(rows):
for c in range(cols):
v = float(img[r, c])
if not (xmin <= v < xmax):
continue
b = int((v - xmin) * inv_range)
if b >= n_bins:
b = n_bins - 1
expected[r, c, b] += 1
return expected
def test_async_fill_matches_reference():
rows, cols, n_bins = 5, 7, 8
xmin, xmax = 0.0, 2.0
frames = _random_frames(rows, cols, n=3, xmin=xmin, xmax=xmax, seed=1)
hist = PixelHistogram(rows=rows, cols=cols, n_bins=n_bins, xmin=xmin, xmax=xmax)
for img in frames:
hist.fill_async(img)
np.testing.assert_array_equal(
hist.values(),
_reference_hdata(frames, rows, cols, n_bins, xmin, xmax),
)
def test_fill_async_copies_buffer():
# After fill_async returns, the caller should be free to mutate the
# numpy array without affecting the pending fill.
rows, cols, n_bins = 4, 4, 4
xmin, xmax = 0.0, 1.0
hist = PixelHistogram(rows=rows, cols=cols, n_bins=n_bins, xmin=xmin, xmax=xmax, n_threads=1, max_pending=8)
img = np.full((rows, cols), 0.1, dtype=np.float64) # falls in bin 0
hist.fill_async(img)
# Mutate the original array immediately; this must not affect the
# value that was already enqueued.
img[:] = 0.9 # would be bin 3
hist.flush()
h = hist.values()
assert h.shape == (rows, cols, n_bins)
# Every pixel saw one value in bin 0, none elsewhere.
assert (h[:, :, 0] == 1).all()
assert (h[:, :, 1:] == 0).all()
def test_fill_async_rejects_wrong_shape():
hist = PixelHistogram(8, 8, 4, 0.0, 1.0)
bad = np.zeros((4, 4), dtype=np.float32)
with pytest.raises(ValueError):
hist.fill_async(bad)
def test_hdata_flushes_pending():
# Submit several frames with a tiny queue and read hdata() without an
# explicit flush(); hdata() must drain everything first.
rows, cols, n_bins = 3, 3, 4
xmin, xmax = 0.0, 1.0
hist = PixelHistogram(rows=rows, cols=cols, n_bins=n_bins, xmin=xmin, xmax=xmax,
n_threads=1, max_pending=1)
frames = _random_frames(rows, cols, n=8, xmin=xmin, xmax=xmax, seed=3)
for img in frames:
hist.fill_async(img)
h = hist.values() # no explicit flush()
np.testing.assert_array_equal(
h, _reference_hdata(frames, rows, cols, n_bins, xmin, xmax)
)
def test_bin_centers_and_edges():
n_bins = 5
xmin, xmax = 0.0, 1.0
hist = PixelHistogram(rows=2, cols=2, n_bins=n_bins, xmin=xmin, xmax=xmax)
edges = hist.bin_edges()
centers = hist.bin_centers()
assert edges.shape == (n_bins + 1,)
assert centers.shape == (n_bins,)
np.testing.assert_allclose(edges, np.linspace(xmin, xmax, n_bins + 1), atol=1e-6)
np.testing.assert_allclose(centers, 0.5 * (edges[:-1] + edges[1:]), atol=1e-6)
src/ClusterFile.test.cpp
+4
View File
@@ -17,6 +17,8 @@ TEST_CASE("Read one frame from a cluster file", "[.with-data]") {
REQUIRE(std::filesystem::exists(fpath));
ClusterFile<Cluster<int32_t, 3, 3>> f(fpath);
CHECK(f.estimate_n_clusters() == 97);
CHECK(f.tell() == 0);
auto clusters = f.read_frame();
CHECK(clusters.size() == 97);
CHECK(clusters.frame_number() == 135);
@@ -250,6 +252,8 @@ TEST_CASE("Read cluster from multiple frame file", "[.with-data]") {
SECTION("Read clusters from both frames") {
ClusterFile<ClusterType> f(fpath);
CHECK(f.estimate_n_clusters() == 8);
CHECK(f.tell() == 0);
auto clusters = f.read_clusters(2);
REQUIRE(clusters.size() == 2);
REQUIRE(clusters.frame_number() == 0);
src/ClusterFinderMT.test.cpp
+79 -2
View File
@@ -21,14 +21,21 @@ class ClusterFinderMTWrapper
public:
ClusterFinderMTWrapper(Shape<2> image_size, PEDESTAL_TYPE nSigma = 5.0,
size_t capacity = 2000, size_t n_threads = 3)
size_t capacity = 2000, size_t n_threads = 3,
size_t queue_depth = 16)
: ClusterFinderMT<ClusterType, FRAME_TYPE, PEDESTAL_TYPE>(
image_size, nSigma, capacity, n_threads) {}
image_size, nSigma, capacity, n_threads, queue_depth) {}
size_t get_m_input_queues_size() const {
return this->m_input_queues.size();
}
size_t get_m_frame_pools_size() const { return this->m_frame_pools.size(); }
size_t get_frame_pool_depth(size_t thread_index) const {
return this->m_frame_pools[thread_index]->size();
}
size_t get_m_output_queues_size() const {
return this->m_output_queues.size();
}
@@ -100,3 +107,73 @@ TEST_CASE("multithreaded cluster finder", "[.with-data]") {
auto clustervec = clustercollector.steal_clusters();
// CHECK(clustervec.size() == ) //dont know how many clusters to expect
}
TEST_CASE("frame buffers are recycled when pushing more frames than the pool "
"holds",
"[.files]") {
using ClusterType = Cluster<int32_t, 3, 3>;
const size_t n_threads = 2;
const size_t queue_depth = 4;
const size_t n_frames = 5 * queue_depth;
const Shape<2> image_size{10, 10};
ClusterFinderMTWrapper<ClusterType> cf(image_size, 5, 200, n_threads,
queue_depth);
CHECK(cf.get_m_frame_pools_size() == n_threads);
for (size_t i = 0; i < n_threads; ++i) {
CHECK(cf.get_frame_pool_depth(i) == queue_depth);
}
NDArray<uint16_t, 2> frame(image_size, 0);
// More frames than the pool holds, so this only completes if the workers
// return the buffers to the free list.
for (size_t i = 0; i < n_frames; ++i) {
cf.find_clusters(frame.view(), i);
}
cf.stop();
CHECK(cf.m_input_queues_are_empty() == true);
}
TEST_CASE("cluster collector accepts finders with a matching cluster type") {
using ClusterType = Cluster<int32_t, 3, 3>;
using Finder = ClusterFinderMTWrapper<ClusterType, uint16_t, float>;
using OtherFinder =
ClusterFinderMTWrapper<Cluster<int32_t, 5, 5>, uint16_t, float>;
static_assert(
std::is_constructible_v<ClusterCollector<ClusterType>, Finder *>);
static_assert(
!std::is_constructible_v<ClusterCollector<ClusterType>, OtherFinder *>);
Finder cf({10, 10});
cf.stop();
ClusterCollector<ClusterType> collector(&cf);
collector.stop();
CHECK(collector.steal_clusters().empty());
}
TEST_CASE("cluster collector drains queued clusters when stopped") {
using ClusterType = Cluster<int32_t, 3, 3>;
ProducerConsumerQueue<ClusterVector<ClusterType>> source(4);
for (uint64_t frame_number = 1; frame_number <= 3; ++frame_number) {
REQUIRE(source.write(ClusterVector<ClusterType>(4, frame_number)));
}
ClusterCollector<ClusterType> collector(&source);
collector.stop();
auto clusters = collector.steal_clusters();
REQUIRE(clusters.size() == 3);
CHECK(source.isEmpty());
for (size_t i = 0; i < clusters.size(); ++i) {
CHECK(clusters[i].frame_number() == static_cast<int32_t>(i + 1));
}
}
src/ClusterVector.test.cpp
+25 -1
View File
@@ -1,6 +1,7 @@
// SPDX-License-Identifier: MPL-2.0
#include "aare/ClusterVector.hpp"
#include <cstdint>
#include <utility>
#include <catch2/catch_all.hpp>
#include <catch2/catch_test_macros.hpp>
@@ -21,6 +22,29 @@ TEST_CASE("After pushing back one element the ClusterVector is not empty") {
REQUIRE(!cv.empty());
}
TEST_CASE("ClusterVector move constructor preserves contents") {
ClusterVector<C1> source(4, 42);
source.push_back(C1{1, 2, {3, 4, 5, 6}});
ClusterVector<C1> moved(std::move(source));
REQUIRE(moved.size() == 1);
CHECK(moved.frame_number() == 42);
CHECK(moved[0].data[3] == 6);
}
TEST_CASE("ClusterVector move assignment preserves contents") {
ClusterVector<C1> source(4, 42);
source.push_back(C1{1, 2, {3, 4, 5, 6}});
ClusterVector<C1> moved(1);
moved = std::move(source);
REQUIRE(moved.size() == 1);
CHECK(moved.frame_number() == 42);
CHECK(moved[0].data[3] == 6);
}
TEST_CASE("item_size return the size of the cluster stored") {
ClusterVector<C1> cv(4);
CHECK(cv.item_size() == sizeof(C1));
@@ -275,4 +299,4 @@ TEST_CASE("Gain Map Calculation Index Map") {
CHECK(index_map_x == clustertestdata.index_map_x);
CHECK(index_map_y == clustertestdata.index_map_y);
}
}
src/Dtype.cpp
+6
View File
@@ -92,6 +92,12 @@ Dtype::Dtype(Dtype::TypeIndex ti) : m_type(ti) {}
*/
Dtype::Dtype(std::string_view sv) {
// NumPy uses '|' for data types whose byte order is not applicable,
// notably one-byte integer types.
if (!sv.empty() && sv.front() == '|') {
sv.remove_prefix(1);
}
// Check if the file is using our native endianess
if (auto pos = sv.find_first_of("<>"); pos != std::string_view::npos) {
const auto endianess = [](const char c) {
src/Dtype.test.cpp
+4 -1
View File
@@ -12,6 +12,9 @@ TEST_CASE("Construct from typeid") {
}
TEST_CASE("Construct from string") {
REQUIRE(Dtype("|i1") == typeid(int8_t));
REQUIRE(Dtype("|u1") == typeid(uint8_t));
if (endian::native == endian::little) {
REQUIRE(Dtype("<i1") == typeid(int8_t));
REQUIRE(Dtype("<u1") == typeid(uint8_t));
@@ -53,4 +56,4 @@ TEST_CASE("Construct from string with endianess") {
TEST_CASE("Convert to string") {
REQUIRE(Dtype(typeid(int)).to_string() == "<i4");
}
}
src/File.cpp
+2 -1
View File
@@ -73,10 +73,11 @@ size_t File::tell() const { return file_impl->tell(); }
size_t File::rows() const { return file_impl->rows(); }
size_t File::cols() const { return file_impl->cols(); }
size_t File::bitdepth() const { return file_impl->bitdepth(); }
Dtype File::dtype() const { return file_impl->dtype(); }
size_t File::bytes_per_pixel() const {
return file_impl->bitdepth() / bits_per_byte;
}
DetectorType File::detector_type() const { return file_impl->detector_type(); }
} // namespace aare
} // namespace aare
src/Fit.cpp
+1 -452
View File
@@ -4,7 +4,6 @@
#include "Minuit2/FunctionMinimum.h"
#include "Minuit2/MnHesse.h"
#include "Minuit2/MnMigrad.h"
#include "Minuit2/MnPrint.h"
#include "Minuit2/MnStrategy.h"
#include "Minuit2/MnUserParameters.h"
#include "aare/Models.hpp"
@@ -12,461 +11,11 @@
#include "aare/utils/task.hpp"
#include <array>
#include <cmath>
#include <lmcurve2.h>
#include <lmfit.hpp>
#include <memory>
#include <stdexcept>
#include <thread>
#include <type_traits>
namespace aare {
namespace func {
double gaus(const double x, const double *par) {
return par[0] * exp(-pow(x - par[1], 2) / (2 * pow(par[2], 2)));
}
NDArray<double, 1> gaus(NDView<double, 1> x, NDView<double, 1> par) {
NDArray<double, 1> y({x.shape(0)}, 0);
for (ssize_t i = 0; i < x.size(); i++) {
y(i) = gaus(x(i), par.data());
}
return y;
}
double pol1(const double x, const double *par) { return par[0] * x + par[1]; }
NDArray<double, 1> pol1(NDView<double, 1> x, NDView<double, 1> par) {
NDArray<double, 1> y({x.shape()}, 0);
for (ssize_t i = 0; i < x.size(); i++) {
y(i) = pol1(x(i), par.data());
}
return y;
}
double scurve(const double x, const double *par) {
return (par[0] + par[1] * x) +
0.5 * (1 + erf((x - par[2]) / (sqrt(2) * par[3]))) *
(par[4] + par[5] * (x - par[2]));
}
NDArray<double, 1> scurve(NDView<double, 1> x, NDView<double, 1> par) {
NDArray<double, 1> y({x.shape()}, 0);
for (ssize_t i = 0; i < x.size(); i++) {
y(i) = scurve(x(i), par.data());
}
return y;
}
double scurve2(const double x, const double *par) {
return (par[0] + par[1] * x) +
0.5 * (1 - erf((x - par[2]) / (sqrt(2) * par[3]))) *
(par[4] + par[5] * (x - par[2]));
}
NDArray<double, 1> scurve2(NDView<double, 1> x, NDView<double, 1> par) {
NDArray<double, 1> y({x.shape()}, 0);
for (ssize_t i = 0; i < x.size(); i++) {
y(i) = scurve2(x(i), par.data());
}
return y;
}
} // namespace func
NDArray<double, 1> fit_gaus(NDView<double, 1> x, NDView<double, 1> y) {
NDArray<double, 1> result = model::Gaussian::estimate_par(x, y);
lm_status_struct status;
lmcurve(result.size(), result.data(), x.size(), x.data(), y.data(),
aare::func::gaus, &lm_control_double, &status);
return result;
}
NDArray<double, 3> fit_gaus(NDView<double, 1> x, NDView<double, 3> y,
int n_threads) {
NDArray<double, 3> result({y.shape(0), y.shape(1), 3}, 0);
auto process = [&x, &y, &result](ssize_t first_row, ssize_t last_row) {
for (ssize_t row = first_row; row < last_row; row++) {
for (ssize_t col = 0; col < y.shape(1); col++) {
NDView<double, 1> values(&y(row, col, 0), {y.shape(2)});
auto res = fit_gaus(x, values);
result(row, col, 0) = res(0);
result(row, col, 1) = res(1);
result(row, col, 2) = res(2);
}
}
};
auto tasks = split_task(0, y.shape(0), n_threads);
RunInParallel(process, tasks);
return result;
}
void fit_gaus(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
NDView<double, 1> par_out, NDView<double, 1> par_err_out,
double &chi2) {
// Check that we have the correct sizes
if (y.size() != x.size() || y.size() != y_err.size() ||
par_out.size() != 3 || par_err_out.size() != 3) {
throw std::runtime_error("Data, x, data_err must have the same size "
"and par_out, par_err_out must have size 3");
}
// /* Collection of output parameters for status info. */
// typedef struct {
// double fnorm; /* norm of the residue vector fvec. */
// int nfev; /* actual number of iterations. */
// int outcome; /* Status indicator. Nonnegative values are used as
// index
// for the message text lm_infmsg, set in lmmin.c. */
// int userbreak; /* Set when function evaluation requests termination.
// */
// } lm_status_struct;
lm_status_struct status;
par_out = model::Gaussian::estimate_par(x, y);
std::array<double, 9> cov{0, 0, 0, 0, 0, 0, 0, 0, 0};
// void lmcurve2( const int n_par, double *par, double *parerr, double
// *covar, const int m_dat, const double *t, const double *y, const double
// *dy, double (*f)( const double ti, const double *par ), const
// lm_control_struct *control, lm_status_struct *status); n_par - Number of
// free variables. Length of parameter vector par. par - Parameter vector.
// On input, it must contain a reasonable guess. On output, it contains the
// solution found to minimize ||r||. parerr - Parameter uncertainties
// vector. Array of length n_par or NULL. On output, unless it or covar is
// NULL, it contains the weighted parameter uncertainties for the found
// parameters. covar - Covariance matrix. Array of length n_par * n_par or
// NULL. On output, unless it is NULL, it contains the covariance matrix.
// m_dat - Number of data points. Length of vectors t, y, dy. Must statisfy
// n_par <= m_dat. t - Array of length m_dat. Contains the abcissae (time,
// or "x") for which function f will be evaluated. y - Array of length
// m_dat. Contains the ordinate values that shall be fitted. dy - Array of
// length m_dat. Contains the standard deviations of the values y. f - A
// user-supplied parametric function f(ti;par). control - Parameter
// collection for tuning the fit procedure. In most cases, the default
// &lm_control_double is adequate. If f is only computed with
// single-precision accuracy, &lm_control_float should be used. Parameters
// are explained in lmmin2(3). status - A record used to return information
// about the minimization process: For details, see lmmin2(3).
lmcurve2(par_out.size(), par_out.data(), par_err_out.data(), cov.data(),
x.size(), x.data(), y.data(), y_err.data(), aare::func::gaus,
&lm_control_double, &status);
// Calculate chi2
chi2 = 0;
for (ssize_t i = 0; i < y.size(); i++) {
chi2 +=
std::pow((y(i) - func::gaus(x(i), par_out.data())) / y_err(i), 2);
}
}
void fit_gaus(NDView<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_err,
NDView<double, 3> par_out, NDView<double, 3> par_err_out,
NDView<double, 2> chi2_out,
int n_threads) {
auto process = [&](ssize_t first_row, ssize_t last_row) {
for (ssize_t row = first_row; row < last_row; row++) {
for (ssize_t col = 0; col < y.shape(1); col++) {
NDView<double, 1> y_view(&y(row, col, 0), {y.shape(2)});
NDView<double, 1> y_err_view(&y_err(row, col, 0),
{y_err.shape(2)});
NDView<double, 1> par_out_view(&par_out(row, col, 0),
{par_out.shape(2)});
NDView<double, 1> par_err_out_view(&par_err_out(row, col, 0),
{par_err_out.shape(2)});
fit_gaus(x, y_view, y_err_view, par_out_view, par_err_out_view,
chi2_out(row, col));
}
}
};
auto tasks = split_task(0, y.shape(0), n_threads);
RunInParallel(process, tasks);
}
void fit_pol1(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
NDView<double, 1> par_out, NDView<double, 1> par_err_out,
double &chi2) {
// Check that we have the correct sizes
if (y.size() != x.size() || y.size() != y_err.size() ||
par_out.size() != 2 || par_err_out.size() != 2) {
throw std::runtime_error("Data, x, data_err must have the same size "
"and par_out, par_err_out must have size 2");
}
lm_status_struct status;
par_out = model::Pol1::estimate_par(x, y);
std::array<double, 4> cov{0, 0, 0, 0};
lmcurve2(par_out.size(), par_out.data(), par_err_out.data(), cov.data(),
x.size(), x.data(), y.data(), y_err.data(), aare::func::pol1,
&lm_control_double, &status);
// Calculate chi2
chi2 = 0;
for (ssize_t i = 0; i < y.size(); i++) {
chi2 +=
std::pow((y(i) - func::pol1(x(i), par_out.data())) / y_err(i), 2);
}
}
void fit_pol1(NDView<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_err,
NDView<double, 3> par_out, NDView<double, 3> par_err_out,
NDView<double, 2> chi2_out, int n_threads) {
auto process = [&](ssize_t first_row, ssize_t last_row) {
for (ssize_t row = first_row; row < last_row; row++) {
for (ssize_t col = 0; col < y.shape(1); col++) {
NDView<double, 1> y_view(&y(row, col, 0), {y.shape(2)});
NDView<double, 1> y_err_view(&y_err(row, col, 0),
{y_err.shape(2)});
NDView<double, 1> par_out_view(&par_out(row, col, 0),
{par_out.shape(2)});
NDView<double, 1> par_err_out_view(&par_err_out(row, col, 0),
{par_err_out.shape(2)});
fit_pol1(x, y_view, y_err_view, par_out_view, par_err_out_view,
chi2_out(row, col));
}
}
};
auto tasks = split_task(0, y.shape(0), n_threads);
RunInParallel(process, tasks);
}
NDArray<double, 1> fit_pol1(NDView<double, 1> x, NDView<double, 1> y) {
// // Check that we have the correct sizes
// if (y.size() != x.size() || y.size() != y_err.size() ||
// par_out.size() != 2 || par_err_out.size() != 2) {
// throw std::runtime_error("Data, x, data_err must have the same size "
// "and par_out, par_err_out must have size 2");
// }
NDArray<double, 1> par = model::Pol1::estimate_par(x, y);
lm_status_struct status;
lmcurve(par.size(), par.data(), x.size(), x.data(), y.data(),
aare::func::pol1, &lm_control_double, &status);
return par;
}
NDArray<double, 3> fit_pol1(NDView<double, 1> x, NDView<double, 3> y,
int n_threads) {
NDArray<double, 3> result({y.shape(0), y.shape(1), 2}, 0);
auto process = [&](ssize_t first_row, ssize_t last_row) {
for (ssize_t row = first_row; row < last_row; row++) {
for (ssize_t col = 0; col < y.shape(1); col++) {
NDView<double, 1> values(&y(row, col, 0), {y.shape(2)});
auto res = fit_pol1(x, values);
result(row, col, 0) = res(0);
result(row, col, 1) = res(1);
}
}
};
auto tasks = split_task(0, y.shape(0), n_threads);
RunInParallel(process, tasks);
return result;
}
// ~~ S-CURVES ~~
// - No error
NDArray<double, 1> fit_scurve(NDView<double, 1> x, NDView<double, 1> y) {
NDArray<double, 1> result = model::RisingScurve::estimate_par(x, y);
lm_status_struct status;
lmcurve(result.size(), result.data(), x.size(), x.data(), y.data(),
aare::func::scurve, &lm_control_double, &status);
return result;
}
NDArray<double, 3> fit_scurve(NDView<double, 1> x, NDView<double, 3> y,
int n_threads) {
NDArray<double, 3> result({y.shape(0), y.shape(1), 6}, 0);
auto process = [&x, &y, &result](ssize_t first_row, ssize_t last_row) {
for (ssize_t row = first_row; row < last_row; row++) {
for (ssize_t col = 0; col < y.shape(1); col++) {
NDView<double, 1> values(&y(row, col, 0), {y.shape(2)});
auto res = fit_scurve(x, values);
result(row, col, 0) = res(0);
result(row, col, 1) = res(1);
result(row, col, 2) = res(2);
result(row, col, 3) = res(3);
result(row, col, 4) = res(4);
result(row, col, 5) = res(5);
}
}
};
auto tasks = split_task(0, y.shape(0), n_threads);
RunInParallel(process, tasks);
return result;
}
// - Error
void fit_scurve(NDView<double, 1> x, NDView<double, 1> y,
NDView<double, 1> y_err, NDView<double, 1> par_out,
NDView<double, 1> par_err_out, double &chi2) {
// Check that we have the correct sizes
if (y.size() != x.size() || y.size() != y_err.size() ||
par_out.size() != 6 || par_err_out.size() != 6) {
throw std::runtime_error("Data, x, data_err must have the same size "
"and par_out, par_err_out must have size 6");
}
lm_status_struct status;
par_out = model::RisingScurve::estimate_par(x, y);
std::array<double, 36> cov = {0}; // size 6x6
// std::array<double, 4> cov{0, 0, 0, 0};
lmcurve2(par_out.size(), par_out.data(), par_err_out.data(), cov.data(),
x.size(), x.data(), y.data(), y_err.data(), aare::func::scurve,
&lm_control_double, &status);
// Calculate chi2
chi2 = 0;
for (ssize_t i = 0; i < y.size(); i++) {
chi2 +=
std::pow((y(i) - func::pol1(x(i), par_out.data())) / y_err(i), 2);
}
}
void fit_scurve(NDView<double, 1> x, NDView<double, 3> y,
NDView<double, 3> y_err, NDView<double, 3> par_out,
NDView<double, 3> par_err_out, NDView<double, 2> chi2_out,
int n_threads) {
auto process = [&](ssize_t first_row, ssize_t last_row) {
for (ssize_t row = first_row; row < last_row; row++) {
for (ssize_t col = 0; col < y.shape(1); col++) {
NDView<double, 1> y_view(&y(row, col, 0), {y.shape(2)});
NDView<double, 1> y_err_view(&y_err(row, col, 0),
{y_err.shape(2)});
NDView<double, 1> par_out_view(&par_out(row, col, 0),
{par_out.shape(2)});
NDView<double, 1> par_err_out_view(&par_err_out(row, col, 0),
{par_err_out.shape(2)});
fit_scurve(x, y_view, y_err_view, par_out_view,
par_err_out_view, chi2_out(row, col));
}
}
};
auto tasks = split_task(0, y.shape(0), n_threads);
RunInParallel(process, tasks);
}
// SCURVE2 ---
// - No error
NDArray<double, 1> fit_scurve2(NDView<double, 1> x, NDView<double, 1> y) {
NDArray<double, 1> result = model::FallingScurve::estimate_par(x, y);
lm_status_struct status;
lmcurve(result.size(), result.data(), x.size(), x.data(), y.data(),
aare::func::scurve2, &lm_control_double, &status);
return result;
}
NDArray<double, 3> fit_scurve2(NDView<double, 1> x, NDView<double, 3> y,
int n_threads) {
NDArray<double, 3> result({y.shape(0), y.shape(1), 6}, 0);
auto process = [&x, &y, &result](ssize_t first_row, ssize_t last_row) {
for (ssize_t row = first_row; row < last_row; row++) {
for (ssize_t col = 0; col < y.shape(1); col++) {
NDView<double, 1> values(&y(row, col, 0), {y.shape(2)});
auto res = fit_scurve2(x, values);
result(row, col, 0) = res(0);
result(row, col, 1) = res(1);
result(row, col, 2) = res(2);
result(row, col, 3) = res(3);
result(row, col, 4) = res(4);
result(row, col, 5) = res(5);
}
}
};
auto tasks = split_task(0, y.shape(0), n_threads);
RunInParallel(process, tasks);
return result;
}
// - Error
void fit_scurve2(NDView<double, 1> x, NDView<double, 1> y,
NDView<double, 1> y_err, NDView<double, 1> par_out,
NDView<double, 1> par_err_out, double &chi2) {
// Check that we have the correct sizes
if (y.size() != x.size() || y.size() != y_err.size() ||
par_out.size() != 6 || par_err_out.size() != 6) {
throw std::runtime_error("Data, x, data_err must have the same size "
"and par_out, par_err_out must have size 6");
}
lm_status_struct status;
par_out = model::FallingScurve::estimate_par(x, y);
std::array<double, 36> cov = {0}; // size 6x6
// std::array<double, 4> cov{0, 0, 0, 0};
lmcurve2(par_out.size(), par_out.data(), par_err_out.data(), cov.data(),
x.size(), x.data(), y.data(), y_err.data(), aare::func::scurve2,
&lm_control_double, &status);
// Calculate chi2
chi2 = 0;
for (ssize_t i = 0; i < y.size(); i++) {
chi2 +=
std::pow((y(i) - func::pol1(x(i), par_out.data())) / y_err(i), 2);
}
}
void fit_scurve2(NDView<double, 1> x, NDView<double, 3> y,
NDView<double, 3> y_err, NDView<double, 3> par_out,
NDView<double, 3> par_err_out, NDView<double, 2> chi2_out,
int n_threads) {
auto process = [&](ssize_t first_row, ssize_t last_row) {
for (ssize_t row = first_row; row < last_row; row++) {
for (ssize_t col = 0; col < y.shape(1); col++) {
NDView<double, 1> y_view(&y(row, col, 0), {y.shape(2)});
NDView<double, 1> y_err_view(&y_err(row, col, 0),
{y_err.shape(2)});
NDView<double, 1> par_out_view(&par_out(row, col, 0),
{par_out.shape(2)});
NDView<double, 1> par_err_out_view(&par_err_out(row, col, 0),
{par_err_out.shape(2)});
fit_scurve2(x, y_view, y_err_view, par_out_view,
par_err_out_view, chi2_out(row, col));
}
}
};
auto tasks = split_task(0, y.shape(0), n_threads);
RunInParallel(process, tasks);
}
// ============================================================================
// FitModel<Model> — method definitions
// (constructor, destructor, copy, and all methods that touch Minuit2 state)
@@ -756,4 +305,4 @@ AARE_INSTANTIATE_FIT(model::FallingScurve)
#undef AARE_INSTANTIATE_FIT
// NOLINTEND
} // namespace aare
} // namespace aare
src/Fit.test.cpp
+81
View File
@@ -0,0 +1,81 @@
// SPDX-License-Identifier: MPL-2.0
#include "aare/Fit.hpp"
#include "aare/FitModel.hpp"
#include "aare/Models.hpp"
#include <catch2/catch_approx.hpp>
#include <catch2/catch_test_macros.hpp>
#include <cmath>
namespace {
constexpr ssize_t n_points = 61;
void fill_gaussian(aare::NDArray<double, 1> &x, aare::NDArray<double, 1> &y,
double amplitude, double mean, double sigma) {
for (ssize_t i = 0; i < n_points; ++i) {
x(i) = -6.0 + 0.2 * static_cast<double>(i);
const double z = (x(i) - mean) / sigma;
y(i) = amplitude * std::exp(-0.5 * z * z);
}
}
} // namespace
TEST_CASE("Minuit2 fits weighted and unweighted Gaussian data", "[fit]") {
aare::NDArray<double, 1> x({n_points});
aare::NDArray<double, 1> y({n_points});
aare::NDArray<double, 1> y_err({n_points}, 1.0);
fill_gaussian(x, y, 120.0, 0.8, 1.3);
const aare::FitModel<aare::model::Gaussian> unweighted_model;
const auto unweighted =
aare::fit_pixel(unweighted_model, x.view(), y.view());
REQUIRE(unweighted.size() == 4);
CHECK(unweighted(0) == Catch::Approx(120.0).epsilon(1e-5));
CHECK(unweighted(1) == Catch::Approx(0.8).epsilon(1e-5));
CHECK(unweighted(2) == Catch::Approx(1.3).epsilon(1e-5));
CHECK(unweighted(3) == Catch::Approx(0.0).margin(1e-4));
const aare::FitModel<aare::model::Gaussian> weighted_model(0, 100, 0.5,
true);
const auto weighted =
aare::fit_pixel(weighted_model, x.view(), y.view(), y_err.view());
REQUIRE(weighted.size() == 7);
CHECK(weighted(0) == Catch::Approx(120.0).epsilon(1e-5));
CHECK(weighted(1) == Catch::Approx(0.8).epsilon(1e-5));
CHECK(weighted(2) == Catch::Approx(1.3).epsilon(1e-5));
CHECK(weighted(6) == Catch::Approx(0.0).margin(1e-4));
}
TEST_CASE("Minuit2 fits a Gaussian data cube in parallel", "[fit]") {
aare::NDArray<double, 1> x({n_points});
aare::NDArray<double, 1> values({n_points});
fill_gaussian(x, values, 80.0, -0.6, 0.9);
aare::NDArray<double, 3> y({2, 2, n_points});
for (ssize_t row = 0; row < 2; ++row) {
for (ssize_t col = 0; col < 2; ++col) {
for (ssize_t i = 0; i < n_points; ++i)
y(row, col, i) = values(i);
}
}
aare::NDArray<double, 3> par({2, 2, 3});
aare::NDArray<double, 2> chi2({2, 2});
const aare::FitModel<aare::model::Gaussian> model;
aare::fit_3d(model, x.view(), y.view(), aare::NDView<double, 3>{},
par.view(), aare::NDView<double, 3>{}, chi2.view(), 2);
for (ssize_t row = 0; row < 2; ++row) {
for (ssize_t col = 0; col < 2; ++col) {
CHECK(par(row, col, 0) == Catch::Approx(80.0).epsilon(1e-5));
CHECK(par(row, col, 1) == Catch::Approx(-0.6).epsilon(1e-5));
CHECK(par(row, col, 2) == Catch::Approx(0.9).epsilon(1e-5));
CHECK(chi2(row, col) == Catch::Approx(0.0).margin(1e-8));
}
}
}
src/MultiThreadedFileReader.cpp
+151
View File
@@ -0,0 +1,151 @@
// SPDX-License-Identifier: MPL-2.0
#include "aare/MultiThreadedFileReader.hpp"
#include "aare/File.hpp"
#include "aare/utils/math_helpers.hpp"
#include <algorithm>
#include <future>
#include <limits>
#include <stdexcept>
#include <utility>
namespace aare::experimental {
namespace {
size_t checked_product(size_t lhs, size_t rhs) {
if (lhs != 0 && rhs > std::numeric_limits<size_t>::max() / lhs) {
throw std::overflow_error(
"MultiThreadedFileReader buffer size overflow");
}
return lhs * rhs;
}
} // namespace
MultiThreadedFileReader::MultiThreadedFileReader(
std::filesystem::path fname, size_t n_threads, size_t chunk_size,
std::optional<size_t> total_frames)
: m_fname(std::move(fname)), m_n_threads(n_threads),
m_chunk_size(chunk_size), m_total_frames(0), m_source_total_frames(0),
m_rows(0), m_cols(0), m_bitdepth(0), m_dtype(Dtype::NONE),
m_bytes_per_frame(0), m_total_bytes(0), m_current_frame(0) {
if (m_n_threads == 0) {
throw std::invalid_argument(
"MultiThreadedFileReader requires at least one thread");
}
if (m_chunk_size == 0) {
throw std::invalid_argument(
"MultiThreadedFileReader chunk size must be greater than zero");
}
File file(m_fname);
m_source_total_frames = file.total_frames();
m_total_frames = total_frames.value_or(m_source_total_frames);
if (m_total_frames > m_source_total_frames) {
throw std::invalid_argument(
"Requested frame count exceeds the number of frames in the file");
}
m_rows = file.rows();
m_cols = file.cols();
m_bitdepth = file.bitdepth();
m_dtype = file.dtype();
m_bytes_per_frame = file.bytes_per_frame();
m_total_bytes = checked_product(m_total_frames, m_bytes_per_frame);
m_files.reserve(m_n_threads);
m_files.push_back(std::move(file));
for (size_t i = 1; i < m_n_threads; ++i) {
m_files.emplace_back(m_fname);
}
}
size_t MultiThreadedFileReader::remaining_frames() const noexcept {
return m_total_frames - m_current_frame;
}
size_t MultiThreadedFileReader::next_read_frames() const noexcept {
const size_t remaining = remaining_frames();
if (remaining == 0) {
return 0;
}
const size_t chunks_remaining = ceil_div(remaining, m_chunk_size);
if (m_n_threads >= chunks_remaining) {
return remaining;
}
// This multiplication is safe: in this branch n_threads * chunk_size is
// strictly smaller than remaining.
return m_n_threads * m_chunk_size;
}
void MultiThreadedFileReader::ensure_open() const {
if (!is_open()) {
throw std::runtime_error("MultiThreadedFileReader is closed");
}
}
size_t MultiThreadedFileReader::read_into(std::byte *destination) {
ensure_open();
const size_t frames_to_read = next_read_frames();
if (frames_to_read == 0) {
return 0;
}
if (destination == nullptr) {
throw std::invalid_argument(
"MultiThreadedFileReader destination must not be null");
}
const size_t first_frame = m_current_frame;
const size_t active_threads = ceil_div(frames_to_read, m_chunk_size);
auto worker = [&](size_t worker_index) {
File &file = m_files[worker_index];
const size_t batch_offset = worker_index * m_chunk_size;
const size_t begin = first_frame + batch_offset;
const size_t count =
std::min(m_chunk_size, frames_to_read - batch_offset);
file.seek(begin);
file.read_into(destination + batch_offset * m_bytes_per_frame, count);
};
std::vector<std::future<void>> workers;
workers.reserve(active_threads);
for (size_t i = 0; i < active_threads; ++i) {
workers.emplace_back(std::async(std::launch::async, worker, i));
}
for (auto &future : workers) {
future.get();
}
m_current_frame += frames_to_read;
return frames_to_read;
}
std::vector<std::byte> MultiThreadedFileReader::read() {
std::vector<std::byte> data(next_read_bytes());
read_into(data.data());
return data;
}
std::vector<std::byte> MultiThreadedFileReader::read_all() {
ensure_open();
std::vector<std::byte> data(
checked_product(remaining_frames(), m_bytes_per_frame));
size_t offset = 0;
while (remaining_frames() != 0) {
const size_t frames_read = read_into(data.data() + offset);
offset += frames_read * m_bytes_per_frame;
}
return data;
}
void MultiThreadedFileReader::seek(size_t frame_index) {
ensure_open();
if (frame_index > m_total_frames) {
throw std::out_of_range(
"MultiThreadedFileReader frame index is out of range");
}
m_current_frame = frame_index;
}
} // namespace aare::experimental
src/MultiThreadedFileReader.test.cpp
+189
View File
@@ -0,0 +1,189 @@
// SPDX-License-Identifier: MPL-2.0
#include "aare/MultiThreadedFileReader.hpp"
#include "aare/Dtype.hpp"
#include "aare/File.hpp"
#include "aare/FileInterface.hpp"
#include "aare/Frame.hpp"
#include "aare/NumpyFile.hpp"
#include "test_config.hpp"
#include <catch2/catch_test_macros.hpp>
#include <chrono>
#include <cstddef>
#include <cstdint>
#include <filesystem>
#include <stdexcept>
#include <string>
#include <vector>
using aare::File;
using aare::FileConfig;
using aare::Frame;
using aare::NumpyFile;
using aare::experimental::MultiThreadedFileReader;
namespace {
class TemporaryNumpyFile {
public:
TemporaryNumpyFile() {
const auto unique =
std::chrono::steady_clock::now().time_since_epoch().count();
m_path = std::filesystem::temp_directory_path() /
("aare-mt-reader-" + std::to_string(unique) + ".npy");
FileConfig cfg;
cfg.dtype = aare::Dtype::UINT16;
cfg.rows = 2;
cfg.cols = 3;
NumpyFile file(m_path, "w", cfg);
for (uint16_t frame_index = 0; frame_index < 10; ++frame_index) {
Frame frame(cfg.rows, cfg.cols, cfg.dtype);
auto image = frame.view<uint16_t>();
for (ssize_t row = 0; row < image.shape(0); ++row) {
for (ssize_t col = 0; col < image.shape(1); ++col) {
image(row, col) = static_cast<uint16_t>(frame_index * 100 +
row * 10 + col);
}
}
file.write(frame);
}
}
TemporaryNumpyFile(const TemporaryNumpyFile &) = delete;
TemporaryNumpyFile &operator=(const TemporaryNumpyFile &) = delete;
~TemporaryNumpyFile() { std::filesystem::remove(m_path); }
const std::filesystem::path &path() const { return m_path; }
void truncate() { std::filesystem::resize_file(m_path, 0); }
private:
std::filesystem::path m_path;
};
std::vector<std::byte> read_reference(const std::filesystem::path &fpath,
size_t n_frames) {
File file(fpath);
std::vector<std::byte> data(n_frames * file.bytes_per_frame());
if (n_frames != 0) {
file.read_into(data.data(), n_frames);
}
return data;
}
} // namespace
TEST_CASE("Multi-threaded reader preserves numpy frame order",
"[MultiThreadedFileReader]") {
TemporaryNumpyFile file;
const auto &fpath = file.path();
MultiThreadedFileReader reader(fpath, 2, 3);
CHECK(reader.n_threads() == 2);
CHECK(reader.chunk_size() == 3);
CHECK(reader.total_frames() == 10);
CHECK(reader.source_total_frames() == 10);
CHECK(reader.rows() == 2);
CHECK(reader.cols() == 3);
CHECK(reader.bitdepth() == 16);
CHECK(reader.total_bytes() ==
reader.total_frames() * reader.bytes_per_frame());
const auto reference = read_reference(fpath, reader.total_frames());
auto first = reader.read();
auto second = reader.read();
CHECK(first == std::vector<std::byte>(reference.begin(),
reference.begin() + 6 * 12));
CHECK(second ==
std::vector<std::byte>(reference.begin() + 6 * 12, reference.end()));
CHECK(reader.read().empty());
CHECK(reader.tell() == 10);
CHECK(reader.remaining_frames() == 0);
}
TEST_CASE("Multi-threaded reader handles uneven raw chunks and frame limits",
"[.with-data][MultiThreadedFileReader]") {
const auto fpath =
test_data_path() / "raw/jungfrau/jungfrau_single_master_0.json";
REQUIRE(std::filesystem::exists(fpath));
MultiThreadedFileReader reader(fpath, 2, 3, 9);
CHECK(reader.total_frames() == 9);
CHECK(reader.source_total_frames() == 10);
CHECK(reader.read_all() == read_reference(fpath, 9));
}
TEST_CASE("Multi-threaded reader can seek and reread",
"[MultiThreadedFileReader]") {
TemporaryNumpyFile file;
MultiThreadedFileReader reader(file.path(), 2, 1);
const auto first = reader.read();
CHECK(reader.tell() == 2);
CHECK(reader.next_read_frames() == 2);
reader.seek(0);
CHECK(reader.tell() == 0);
CHECK(reader.read() == first);
CHECK_THROWS_AS(reader.seek(11), std::out_of_range);
}
TEST_CASE("Multi-threaded reader validates its configuration",
"[MultiThreadedFileReader]") {
TemporaryNumpyFile file;
const auto &fpath = file.path();
CHECK_THROWS_AS(MultiThreadedFileReader(fpath, 0, 1),
std::invalid_argument);
CHECK_THROWS_AS(MultiThreadedFileReader(fpath, 1, 0),
std::invalid_argument);
CHECK_THROWS_AS(MultiThreadedFileReader(fpath, 1, 1, 11),
std::invalid_argument);
MultiThreadedFileReader reader(fpath, 2, 2);
CHECK_THROWS_AS(reader.read_into(nullptr), std::invalid_argument);
}
TEST_CASE("An explicit zero frame limit produces an empty read",
"[MultiThreadedFileReader]") {
TemporaryNumpyFile file;
MultiThreadedFileReader reader(file.path(), 8, 3, 0);
CHECK(reader.total_frames() == 0);
CHECK(reader.total_bytes() == 0);
CHECK(reader.read().empty());
CHECK_NOTHROW(reader.read_into(nullptr));
}
TEST_CASE("read_into reads at most one chunk per worker",
"[MultiThreadedFileReader]") {
TemporaryNumpyFile file;
MultiThreadedFileReader reader(file.path(), 3, 2);
std::vector<std::byte> data(reader.next_read_bytes());
CHECK(reader.next_read_frames() == 6);
CHECK(reader.read_into(data.data()) == 6);
CHECK(reader.tell() == 6);
CHECK(reader.next_read_frames() == 4);
CHECK(reader.next_read_bytes() == 4 * reader.bytes_per_frame());
}
TEST_CASE("Multi-threaded reader can close its worker files",
"[MultiThreadedFileReader]") {
TemporaryNumpyFile file;
MultiThreadedFileReader reader(file.path(), 2, 2);
CHECK(reader.is_open());
reader.close();
CHECK_FALSE(reader.is_open());
CHECK_NOTHROW(reader.close());
CHECK_THROWS_AS(reader.read(), std::runtime_error);
CHECK_THROWS_AS(reader.read_all(), std::runtime_error);
CHECK_THROWS_AS(reader.seek(0), std::runtime_error);
}
src/NDArray.test.cpp
+15
View File
@@ -4,6 +4,7 @@
#include <catch2/benchmark/catch_benchmark.hpp>
#include <catch2/catch_test_macros.hpp>
#include <numeric>
#include <type_traits>
using aare::NDArray;
using aare::NDView;
@@ -71,6 +72,20 @@ TEST_CASE("Accessing a const object") {
REQUIRE(img.shape(2) == 5);
}
TEST_CASE("A const NDArray returns a read-only view") {
NDArray<int, 2> mutable_array({2, 3}, 1);
const auto &const_array = mutable_array;
static_assert(
std::is_same_v<decltype(mutable_array.view()), NDView<int, 2>>);
static_assert(
std::is_same_v<decltype(const_array.view()), NDView<const int, 2>>);
auto view = const_array.view();
static_assert(std::is_same_v<decltype(view.data()), const int *>);
REQUIRE(view(1, 2) == 1);
}
TEST_CASE("Indexing of a 2D image") {
std::array<ssize_t, 2> shape{{3, 7}};
NDArray<long> img(shape, 5);
src/NDView.test.cpp
+16 -1
View File
@@ -5,12 +5,16 @@
#include <cstddef>
#include <iostream>
#include <numeric>
#include <type_traits>
#include <vector>
using aare::NDView;
using aare::num_elements;
using aare::Shape;
static_assert(std::is_convertible_v<NDView<int, 2>, NDView<const int, 2>>);
static_assert(!std::is_convertible_v<NDView<const int, 2>, NDView<int, 2>>);
TEST_CASE("Calculate size from a shape") {
Shape<3> shape{2, 3, 4};
REQUIRE(num_elements(shape) == 24);
@@ -88,6 +92,17 @@ TEST_CASE("Element reference 1D with a const NDView") {
}
}
TEST_CASE("Convert a mutable NDView to a const NDView") {
std::vector<int> vec{1, 2, 3, 4};
NDView<int, 2> mutable_view(vec.data(), Shape<2>{2, 2});
NDView<const int, 2> const_view = mutable_view;
REQUIRE(const_view.data() == mutable_view.data());
REQUIRE(const_view.shape() == mutable_view.shape());
REQUIRE(const_view.strides() == mutable_view.strides());
}
TEST_CASE("Element reference 2D") {
std::vector<int> vec(12);
std::iota(vec.begin(), vec.end(), 0);
@@ -279,4 +294,4 @@ TEST_CASE("NDView over byte") {
auto v = aare::make_view(buf);
REQUIRE(v.shape()[0] == 5);
REQUIRE(v[0] == std::byte{0});
}
}
src/NumpyHelpers.test.cpp
+4 -1
View File
@@ -34,6 +34,9 @@ TEST_CASE("trim whitespace") {
}
TEST_CASE("parse data type descriptions") {
REQUIRE(parse_descr("|i1") == aare::Dtype::INT8);
REQUIRE(parse_descr("|u1") == aare::Dtype::UINT8);
REQUIRE(parse_descr("<i1") == aare::Dtype::INT8);
REQUIRE(parse_descr("<i2") == aare::Dtype::INT16);
REQUIRE(parse_descr("<i4") == aare::Dtype::INT32);
@@ -63,4 +66,4 @@ TEST_CASE("Parse numpy dict") {
REQUIRE(map["descr"] == "'<f4'");
REQUIRE(map["fortran_order"] == "False");
REQUIRE(map["shape"] == "(3, 4)");
}
}
src/RawMasterFile.cpp
+21 -8
View File
@@ -303,18 +303,31 @@ void RawMasterFile::parse_json(std::istream &is) {
}
// ----------------------------------------------------------------
// Special treatment of analog flag because of Moench03
m_analog_flag = v < 8.0 && (m_type == DetectorType::Moench);
// Special treatment of analog flag because of Moench03.
// Before SW 8.0.0 (NOT json file version v8!) Moench03 had Analog Samples
// but no Analog Flag. Therefore we need to set analog flag to true and try
// to read analog samples The detector type is later set depending on analog
// samples and number of pixels.
try {
m_analog_flag = static_cast<bool>(j.at("Analog Flag").get<int>());
if (m_analog_flag) {
m_analog_samples = j.at("Analog Samples");
if (m_type == DetectorType::Moench) {
if (auto it = j.find("Analog Samples"); it != j.end()) {
m_analog_flag = true;
m_analog_samples = it->get<size_t>();
} else {
m_analog_flag = false;
}
} else {
if (auto it = j.find("Analog Flag"); it != j.end()) {
m_analog_flag = static_cast<bool>(it->get<int>());
if (m_analog_flag) {
m_analog_samples = j.at("Analog Samples");
}
} else {
m_analog_flag = false;
}
} catch (const json::out_of_range &e) {
// keep the optional empty
}
//-----------------------------------------------------------------
try {
m_quad = j.at("Quad");
} catch (const json::out_of_range &e) {
src/RawMasterFile.test.cpp
+548 -1
View File
@@ -411,6 +411,10 @@ TEST_CASE("Parse EIGER 7.2 master from string stream") {
REQUIRE(f.frame_padding() == 1);
REQUIRE(f.total_frames_expected() == 3);
REQUIRE(f.quad() == 0);
REQUIRE(f.number_of_rows() == 256);
REQUIRE(f.n_modules() == 4);
REQUIRE(f.bitdepth() == 32);
REQUIRE(f.frames_in_file() == 3);
@@ -489,9 +493,473 @@ TEST_CASE("Parse JUNGFRAU 7.2 master from string stream") {
REQUIRE(f.number_of_rows() == 512);
REQUIRE(f.frames_in_file() == 10);
REQUIRE(f.quad() == 0);
REQUIRE(f.n_modules() == 2);
REQUIRE(f.udp_interfaces_per_module() == xy{2, 1});
}
TEST_CASE(
"Parse CTB 7.2 master (SW 7.0.3) with analog samples from string stream") {
std::string master_content = R"({
"Version": 7.2,
"Timestamp": "Wed Aug 19 09:45:29 2026",
"Detector Type": "ChipTestBoard",
"Timing Mode": "auto",
"Geometry": {
"x": 1,
"y": 1
},
"Image Size in bytes": 192000,
"Pixels": {
"x": 32,
"y": 1
},
"Max Frames Per File": 20000,
"Frame Discard Policy": "nodiscard",
"Frame Padding": 1,
"Scan Parameters": "[disabled]",
"Total Frames": 1,
"Receiver Roi": {
"xmin": 4294967295,
"xmax": 4294967295,
"ymin": 4294967295,
"ymax": 4294967295
},
"Exptime": "0ns",
"Period": "1ms",
"Ten Giga": 0,
"ADC Mask": "0xffffffff",
"Analog Flag": 1,
"Analog Samples": 3000,
"Digital Flag": 0,
"Digital Samples": 2000,
"Dbit Offset": 0,
"Dbit Bitset": 0,
"Frames in File": 1,
"Frame Header Format": {
"Frame Number": "8 bytes",
"SubFrame Number/ExpLength": "4 bytes",
"Packet Number": "4 bytes",
"Bunch ID": "8 bytes",
"Timestamp": "8 bytes",
"Module Id": "2 bytes",
"Row": "2 bytes",
"Column": "2 bytes",
"Reserved": "2 bytes",
"Debug": "4 bytes",
"Round Robin Number": "2 bytes",
"Detector Type": "1 byte",
"Header Version": "1 byte",
"Packets Caught Mask": "64 bytes"
}
})";
std::istringstream iss(master_content);
RawMasterFile f(iss, "test_master_0.json");
REQUIRE(f.version() == "7.2");
REQUIRE(f.detector_type() == DetectorType::ChipTestBoard);
REQUIRE(f.timing_mode() == TimingMode::Auto);
REQUIRE(f.geometry() == xy{1, 1});
REQUIRE(f.image_size_in_bytes() == 192000);
REQUIRE(f.pixels_x() == 32);
REQUIRE(f.pixels_y() == 1);
REQUIRE(f.max_frames_per_file() == 20000);
REQUIRE(f.frame_discard_policy() == FrameDiscardPolicy::NoDiscard);
REQUIRE(f.frame_padding() == 1);
REQUIRE(f.total_frames_expected() == 1);
REQUIRE(f.exptime() == std::chrono::nanoseconds(0));
REQUIRE(f.period() == std::chrono::milliseconds(1));
REQUIRE(f.analog_samples() == 3000);
REQUIRE(f.digital_samples() == std::nullopt);
REQUIRE(f.transceiver_samples() == std::nullopt);
REQUIRE(f.frames_in_file() == 1);
}
TEST_CASE(
"Parse CTB 7.2 master (SW 7.0.3) with digital samples from string stream") {
std::string master_content = R"({
"Version": 7.2,
"Timestamp": "Wed Aug 19 09:48:34 2026",
"Detector Type": "ChipTestBoard",
"Timing Mode": "auto",
"Geometry": {
"x": 1,
"y": 1
},
"Image Size in bytes": 16000,
"Pixels": {
"x": 64,
"y": 1
},
"Max Frames Per File": 20000,
"Frame Discard Policy": "nodiscard",
"Frame Padding": 1,
"Scan Parameters": "[disabled]",
"Total Frames": 1,
"Receiver Roi": {
"xmin": 4294967295,
"xmax": 4294967295,
"ymin": 4294967295,
"ymax": 4294967295
},
"Exptime": "0ns",
"Period": "1ms",
"Ten Giga": 0,
"ADC Mask": "0xffffffff",
"Analog Flag": 0,
"Analog Samples": 3000,
"Digital Flag": 1,
"Digital Samples": 2000,
"Dbit Offset": 0,
"Dbit Bitset": 0,
"Frames in File": 1,
"Frame Header Format": {
"Frame Number": "8 bytes",
"SubFrame Number/ExpLength": "4 bytes",
"Packet Number": "4 bytes",
"Bunch ID": "8 bytes",
"Timestamp": "8 bytes",
"Module Id": "2 bytes",
"Row": "2 bytes",
"Column": "2 bytes",
"Reserved": "2 bytes",
"Debug": "4 bytes",
"Round Robin Number": "2 bytes",
"Detector Type": "1 byte",
"Header Version": "1 byte",
"Packets Caught Mask": "64 bytes"
}
})";
std::istringstream iss(master_content);
RawMasterFile f(iss, "test_master_0.json");
REQUIRE(f.version() == "7.2");
REQUIRE(f.detector_type() == DetectorType::ChipTestBoard);
REQUIRE(f.timing_mode() == TimingMode::Auto);
REQUIRE(f.geometry() == xy{1, 1});
REQUIRE(f.image_size_in_bytes() == 16000);
REQUIRE(f.pixels_x() == 64);
REQUIRE(f.pixels_y() == 1);
REQUIRE(f.max_frames_per_file() == 20000);
REQUIRE(f.frame_discard_policy() == FrameDiscardPolicy::NoDiscard);
REQUIRE(f.frame_padding() == 1);
REQUIRE(f.total_frames_expected() == 1);
REQUIRE(f.exptime() == std::chrono::nanoseconds(0));
REQUIRE(f.period() == std::chrono::milliseconds(1));
REQUIRE(f.analog_samples() == std::nullopt);
REQUIRE(f.digital_samples() == 2000);
REQUIRE(f.transceiver_samples() == std::nullopt);
REQUIRE(f.frames_in_file() == 1);
}
TEST_CASE("Parse Moench 7.2 master (SW 7.0.3) from string stream") {
std::string master_content = R"({
"Version": 7.2,
"Timestamp": "Wed Aug 19 10:32:06 2026",
"Detector Type": "Moench",
"Timing Mode": "auto",
"Geometry": {
"x": 1,
"y": 1
},
"Image Size in bytes": 320000,
"Pixels": {
"x": 400,
"y": 400
},
"Max Frames Per File": 100000,
"Frame Discard Policy": "nodiscard",
"Frame Padding": 1,
"Scan Parameters": "[disabled]",
"Total Frames": 1,
"Receiver Roi": {
"xmin": 4294967295,
"xmax": 4294967295,
"ymin": 4294967295,
"ymax": 4294967295
},
"Exptime": "20us",
"Period": "2ms",
"Ten Giga": 0,
"ADC Mask": "0xffffffff",
"Analog Samples": 5000,
"Frames in File": 1,
"Frame Header Format": {
"Frame Number": "8 bytes",
"SubFrame Number/ExpLength": "4 bytes",
"Packet Number": "4 bytes",
"Bunch ID": "8 bytes",
"Timestamp": "8 bytes",
"Module Id": "2 bytes",
"Row": "2 bytes",
"Column": "2 bytes",
"Reserved": "2 bytes",
"Debug": "4 bytes",
"Round Robin Number": "2 bytes",
"Detector Type": "1 byte",
"Header Version": "1 byte",
"Packets Caught Mask": "64 bytes"
}
})";
std::istringstream iss(master_content);
RawMasterFile f(iss, "test_master_0.json");
REQUIRE(f.version() == "7.2");
REQUIRE(f.detector_type() == DetectorType::Moench03_old);
REQUIRE(f.timing_mode() == TimingMode::Auto);
REQUIRE(f.geometry() == xy{1, 1});
REQUIRE(f.image_size_in_bytes() == 320000);
REQUIRE(f.pixels_x() == 400);
REQUIRE(f.pixels_y() == 400);
REQUIRE(f.max_frames_per_file() == 100000);
REQUIRE(f.frame_discard_policy() == FrameDiscardPolicy::NoDiscard);
REQUIRE(f.frame_padding() == 1);
REQUIRE(f.total_frames_expected() == 1);
REQUIRE(f.exptime() == std::chrono::microseconds(20));
REQUIRE(f.period() == std::chrono::milliseconds(2));
REQUIRE(f.analog_samples() == 5000);
REQUIRE(f.digital_samples() == std::nullopt);
REQUIRE(f.transceiver_samples() == std::nullopt);
REQUIRE(f.frames_in_file() == 1);
}
TEST_CASE("Parse Moench 7.2 master (SW 8.0.0) from string stream") {
std::string master_content = R"({
"Version": 7.2,
"Timestamp": "Wed Aug 19 09:54:53 2026",
"Detector Type": "Moench",
"Timing Mode": "auto",
"Geometry": {
"x": 1,
"y": 1
},
"Image Size in bytes": 320000,
"Pixels": {
"x": 400,
"y": 400
},
"Max Frames Per File": 100000,
"Frame Discard Policy": "discardpartial",
"Frame Padding": 1,
"Scan Parameters": "[disabled]",
"Total Frames": 1,
"Receiver Roi": {
"xmin": 4294967295,
"xmax": 4294967295,
"ymin": 4294967295,
"ymax": 4294967295
},
"Exptime": "10us",
"Period": "2ms",
"Number of UDP Interfaces": 1,
"Number of rows": 400,
"Frames in File": 1,
"Frame Header Format": {
"Frame Number": "8 bytes",
"SubFrame Number/ExpLength": "4 bytes",
"Packet Number": "4 bytes",
"Bunch ID": "8 bytes",
"Timestamp": "8 bytes",
"Module Id": "2 bytes",
"Row": "2 bytes",
"Column": "2 bytes",
"Reserved": "2 bytes",
"Debug": "4 bytes",
"Round Robin Number": "2 bytes",
"Detector Type": "1 byte",
"Header Version": "1 byte",
"Packets Caught Mask": "64 bytes"
}
})";
std::istringstream iss(master_content);
RawMasterFile f(iss, "test_master_0.json");
REQUIRE(f.version() == "7.2");
REQUIRE(f.detector_type() == DetectorType::Moench03);
REQUIRE(f.timing_mode() == TimingMode::Auto);
REQUIRE(f.geometry() == xy{1, 1});
REQUIRE(f.image_size_in_bytes() == 320000);
REQUIRE(f.pixels_x() == 400);
REQUIRE(f.pixels_y() == 400);
REQUIRE(f.max_frames_per_file() == 100000);
REQUIRE(f.frame_discard_policy() == FrameDiscardPolicy::DiscardPartial);
REQUIRE(f.frame_padding() == 1);
REQUIRE(f.total_frames_expected() == 1);
REQUIRE(f.exptime() == std::chrono::microseconds(10));
REQUIRE(f.period() == std::chrono::milliseconds(2));
REQUIRE(f.number_of_rows() == 400);
REQUIRE(f.frames_in_file() == 1);
REQUIRE(f.analog_samples() == std::nullopt);
REQUIRE(f.digital_samples() == std::nullopt);
REQUIRE(f.transceiver_samples() == std::nullopt);
REQUIRE(f.udp_interfaces_per_module() == xy{1, 1});
}
TEST_CASE("Parse CTB 7.2 master (SW 8.0.0) from string stream") {
std::string master_content = R"({
"Version": 7.2,
"Timestamp": "Wed Aug 19 09:56:30 2026",
"Detector Type": "ChipTestBoard",
"Timing Mode": "auto",
"Geometry": {
"x": 1,
"y": 1
},
"Image Size in bytes": 192000,
"Pixels": {
"x": 32,
"y": 1
},
"Max Frames Per File": 20000,
"Frame Discard Policy": "nodiscard",
"Frame Padding": 1,
"Scan Parameters": "[disabled]",
"Total Frames": 1,
"Receiver Roi": {
"xmin": 4294967295,
"xmax": 4294967295,
"ymin": 4294967295,
"ymax": 4294967295
},
"Exptime": "0ns",
"Period": "1ms",
"Ten Giga": 0,
"ADC Mask": "0xffffffff",
"Analog Flag": 1,
"Analog Samples": 3000,
"Digital Flag": 0,
"Digital Samples": 2000,
"Dbit Offset": 0,
"Dbit Bitset": 0,
"Transceiver Mask": "0x3",
"Transceiver Flag": 0,
"Transceiver Samples": 1,
"Frames in File": 1,
"Frame Header Format": {
"Frame Number": "8 bytes",
"SubFrame Number/ExpLength": "4 bytes",
"Packet Number": "4 bytes",
"Bunch ID": "8 bytes",
"Timestamp": "8 bytes",
"Module Id": "2 bytes",
"Row": "2 bytes",
"Column": "2 bytes",
"Reserved": "2 bytes",
"Debug": "4 bytes",
"Round Robin Number": "2 bytes",
"Detector Type": "1 byte",
"Header Version": "1 byte",
"Packets Caught Mask": "64 bytes"
}
})";
std::istringstream iss(master_content);
RawMasterFile f(iss, "test_master_0.json");
REQUIRE(f.version() == "7.2");
REQUIRE(f.detector_type() == DetectorType::ChipTestBoard);
REQUIRE(f.timing_mode() == TimingMode::Auto);
REQUIRE(f.geometry() == xy{1, 1});
REQUIRE(f.image_size_in_bytes() == 192000);
REQUIRE(f.pixels_x() == 32);
REQUIRE(f.pixels_y() == 1);
REQUIRE(f.max_frames_per_file() == 20000);
REQUIRE(f.frame_discard_policy() == FrameDiscardPolicy::NoDiscard);
REQUIRE(f.frame_padding() == 1);
REQUIRE(f.total_frames_expected() == 1);
REQUIRE(f.exptime() == std::chrono::nanoseconds(0));
REQUIRE(f.period() == std::chrono::milliseconds(1));
REQUIRE(f.analog_samples() == 3000);
REQUIRE(f.digital_samples() == std::nullopt);
REQUIRE(f.transceiver_samples() == std::nullopt);
REQUIRE(f.frames_in_file() == 1);
}
TEST_CASE(
"Parse CTB 7.2 master (SW 8.0.0) with digital samples from string stream") {
std::string master_content = R"({
"Version": 7.2,
"Timestamp": "Wed Aug 19 09:58:22 2026",
"Detector Type": "ChipTestBoard",
"Timing Mode": "auto",
"Geometry": {
"x": 1,
"y": 1
},
"Image Size in bytes": 16000,
"Pixels": {
"x": 64,
"y": 1
},
"Max Frames Per File": 20000,
"Frame Discard Policy": "nodiscard",
"Frame Padding": 1,
"Scan Parameters": "[disabled]",
"Total Frames": 1,
"Receiver Roi": {
"xmin": 4294967295,
"xmax": 4294967295,
"ymin": 4294967295,
"ymax": 4294967295
},
"Exptime": "0ns",
"Period": "1ms",
"Ten Giga": 0,
"ADC Mask": "0xffffffff",
"Analog Flag": 0,
"Analog Samples": 3000,
"Digital Flag": 1,
"Digital Samples": 2000,
"Dbit Offset": 0,
"Dbit Bitset": 0,
"Transceiver Mask": "0x3",
"Transceiver Flag": 0,
"Transceiver Samples": 1,
"Frames in File": 1,
"Frame Header Format": {
"Frame Number": "8 bytes",
"SubFrame Number/ExpLength": "4 bytes",
"Packet Number": "4 bytes",
"Bunch ID": "8 bytes",
"Timestamp": "8 bytes",
"Module Id": "2 bytes",
"Row": "2 bytes",
"Column": "2 bytes",
"Reserved": "2 bytes",
"Debug": "4 bytes",
"Round Robin Number": "2 bytes",
"Detector Type": "1 byte",
"Header Version": "1 byte",
"Packets Caught Mask": "64 bytes"
}
})";
std::istringstream iss(master_content);
RawMasterFile f(iss, "test_master_0.json");
REQUIRE(f.version() == "7.2");
REQUIRE(f.detector_type() == DetectorType::ChipTestBoard);
REQUIRE(f.timing_mode() == TimingMode::Auto);
REQUIRE(f.geometry() == xy{1, 1});
REQUIRE(f.image_size_in_bytes() == 16000);
REQUIRE(f.pixels_x() == 64);
REQUIRE(f.pixels_y() == 1);
REQUIRE(f.max_frames_per_file() == 20000);
REQUIRE(f.frame_discard_policy() == FrameDiscardPolicy::NoDiscard);
REQUIRE(f.frame_padding() == 1);
REQUIRE(f.total_frames_expected() == 1);
REQUIRE(f.exptime() == std::chrono::nanoseconds(0));
REQUIRE(f.period() == std::chrono::milliseconds(1));
REQUIRE(f.analog_samples() == std::nullopt);
REQUIRE(f.digital_samples() == 2000);
REQUIRE(f.transceiver_samples() == std::nullopt);
REQUIRE(f.frames_in_file() == 1);
}
TEST_CASE("Parse a CTB file from stream") {
std::string master_content = R"({
"Version": 8.0,
@@ -738,4 +1206,83 @@ TEST_CASE("Parse a v7.1 Mythen3 from stream") {
// Period is ok though
REQUIRE(f.period() == std::chrono::milliseconds(2));
}
}
TEST_CASE("Parse old Moench03 from stream") {
std::string master_content = R"(
{
"Version": 7.1,
"Timestamp": "Mon Mar 25 10:07:02 2024",
"Detector Type": "Moench",
"Timing Mode": "auto",
"Geometry": {
"x": 1,
"y": 1
},
"Image Size in bytes": 320000,
"Pixels": {
"x": 400,
"y": 400
},
"Max Frames Per File": 100000,
"Frame Discard Policy": "discardpartial",
"Frame Padding": 1,
"Scan Parameters": "[disabled]",
"Total Frames": 1000000,
"Receiver Roi": {
"xmin": 4294967295,
"xmax": 4294967295,
"ymin": 4294967295,
"ymax": 4294967295
},
"Exptime": "50us",
"Period": "600us",
"Ten Giga": 1,
"ADC Mask": "0xffffffff",
"Analog Samples": 5000,
"Additional Json Header": "{detectorMode: analog, frameMode: newPedestal}",
"Frames in File": 999995,
"Frame Header Format": {
"Frame Number": "8 bytes",
"SubFrame Number/ExpLength": "4 bytes",
"Packet Number": "4 bytes",
"Bunch ID": "8 bytes",
"Timestamp": "8 bytes",
"Module Id": "2 bytes",
"Row": "2 bytes",
"Column": "2 bytes",
"Reserved": "2 bytes",
"Debug": "4 bytes",
"Round Robin Number": "2 bytes",
"Detector Type": "1 byte",
"Header Version": "1 byte",
"Packets Caught Mask": "64 bytes"
}
}
)";
std::istringstream iss(master_content);
RawMasterFile f(iss, "test_master_0.json");
REQUIRE(f.version() == "7.1");
REQUIRE(f.detector_type() == DetectorType::Moench03_old);
REQUIRE(f.timing_mode() == TimingMode::Auto);
REQUIRE(f.geometry().col == 1);
REQUIRE(f.geometry().row == 1);
REQUIRE(f.image_size_in_bytes() == 320000);
REQUIRE(f.pixels_x() == 400);
REQUIRE(f.pixels_y() == 400);
REQUIRE(f.max_frames_per_file() == 100000);
REQUIRE((f.total_frames_expected() ==
1000000)); // This is Total Frames in the master file
REQUIRE(f.frames_in_file() == 999995);
REQUIRE(f.frame_discard_policy() == FrameDiscardPolicy::DiscardPartial);
REQUIRE(f.frame_padding() == 1);
REQUIRE(f.n_modules() == 1);
REQUIRE(f.quad() == 0);
REQUIRE(f.bitdepth() == 16);
REQUIRE(f.exptime() == std::chrono::microseconds(50));
REQUIRE(f.period() == std::chrono::microseconds(600));
REQUIRE(f.analog_samples() == 5000);
}
src/decode.cpp
+2 -2
View File
@@ -144,7 +144,7 @@ uint32_t mask32to24bits(uint32_t input, BitOffset offset) {
return (input >> offset.value()) & mask24bits;
}
void expand4to8bit(NDView<uint8_t, 1> input, NDView<uint8_t, 1> output) {
void expand4to8bit(NDView<const uint8_t, 1> input, NDView<uint8_t, 1> output) {
if (2 * input.size() != output.size())
throw std::runtime_error(
@@ -162,7 +162,7 @@ void expand4to8bit(NDView<uint8_t, 1> input, NDView<uint8_t, 1> output) {
}
}
void expand24to32bit(NDView<uint8_t, 1> input, NDView<uint32_t, 1> output,
void expand24to32bit(NDView<const uint8_t, 1> input, NDView<uint32_t, 1> output,
BitOffset bit_offset) {
ssize_t bytes_per_channel = 3; // 24bit
src/decode.test.cpp
+36 -1
View File
@@ -154,6 +154,20 @@ TEST_CASE("Expand container with 24 bit data to 32") {
}
}
TEST_CASE("Expand 24 bit values to 32 bit values from a const buffer") {
const uint8_t buffer[] = {
0x0F, 0x00, 0x00, 0xFF, 0x00, 0x00, 0xFF, 0xFF, 0xFF,
};
aare::NDView<const uint8_t, 1> input(buffer, {9});
aare::NDArray<uint32_t, 1> out({3});
aare::expand24to32bit(input, out.view());
CHECK(out(0) == 0xF);
CHECK(out(1) == 0xFF);
CHECK(out(2) == 0xFFFFFF);
}
TEST_CASE("Expand 4 bit values packed into 8 bit to 8 bit values") {
{
uint8_t buffer[] = {
@@ -172,4 +186,25 @@ TEST_CASE("Expand 4 bit values packed into 8 bit to 8 bit values") {
CHECK(out(i) == expected_output[i]);
}
}
}
}
TEST_CASE("Expand 4 bit values packed into 8 bit to 8 bit values from a const "
"buffer") {
{
const uint8_t buffer[] = {
0x00, 0xF0, 0xFF, 0x00, 0xF0, 0xFF,
};
aare::NDView<const uint8_t, 1> input(&buffer[0], {6});
aare::NDArray<uint8_t, 1> out({12});
aare::expand4to8bit(input, out.view());
uint8_t expected_output[] = {
0x0, 0x0, 0x0, 0xF, 0xF, 0xF,
0x0, 0x0, 0x0, 0xF, 0xF, 0xF}; // assuming little endian
for (size_t i = 0; i < 12; ++i) {
CHECK(out(i) == expected_output[i]);
}
}
}
src/defs.test.cpp
-14
View File
@@ -22,20 +22,6 @@ TEST_CASE("Enum values") {
REQUIRE(static_cast<int>(aare::DetectorType::Moench03) == 100);
}
TEST_CASE("DynamicCluster creation") {
aare::DynamicCluster c(13, 15);
REQUIRE(c.cluster_sizeX == 13);
REQUIRE(c.cluster_sizeY == 15);
REQUIRE(c.dt == aare::Dtype(typeid(int32_t)));
REQUIRE(c.data() != nullptr);
aare::DynamicCluster c2(c);
REQUIRE(c2.cluster_sizeX == 13);
REQUIRE(c2.cluster_sizeY == 15);
REQUIRE(c2.dt == aare::Dtype(typeid(int32_t)));
REQUIRE(c2.data() != nullptr);
}
TEST_CASE("Basic ops on BitOffset") {
REQUIRE_THROWS(aare::BitOffset(10));
src/hist/PedestalTrackingPixelHistogram.cpp
+161 -91
View File
@@ -1,11 +1,16 @@
#include "aare/hist/PedestalTrackingPixelHistogram.hpp"
#include "aare/File.hpp"
#include "aare/MultiThreadedFileReader.hpp"
#include <algorithm>
#include <array>
#include <chrono>
#include <cmath>
#include <cstring>
#include <functional>
#include <future>
#include <stdexcept>
#include <type_traits>
#include <utility>
#include <vector>
@@ -202,17 +207,8 @@ void PedestalTrackingPixelHistogram::worker_loop(int thread_id) {
switch (kind) {
case WorkKind::PushPedestal: {
// Accumulate raw frame values into this thread's pedestal
// shard. Uses the pixel-level push_no_update which only
// touches m_sum/m_sum2/m_cur_samples (no m_mean writes).
for (int local_row = 0; local_row < local_rows; ++local_row) {
const auto row = static_cast<ssize_t>(first_row + local_row);
for (ssize_t col = 0; col < image->shape(1); ++col) {
my_pedestal.template push_no_update<FrameType>(
static_cast<uint32_t>(local_row),
static_cast<uint32_t>(col), (*image)(row, col));
}
}
auto frame = image->sub_view(first_row, first_row + local_rows);
my_pedestal.add_init_frame(frame);
break;
}
case WorkKind::UpdateMean: {
@@ -220,7 +216,6 @@ void PedestalTrackingPixelHistogram::worker_loop(int thread_id) {
// thread's shard. Also refresh the cached per-pixel std so
// FillWithThreshold can read it without recomputing on the
// hot path.
my_pedestal.update_mean();
auto &my_std = partial_std_[thread_id];
for (int local_row = 0; local_row < local_rows; ++local_row) {
for (int col = 0; col < cols_; ++col) {
@@ -240,59 +235,63 @@ void PedestalTrackingPixelHistogram::worker_loop(int thread_id) {
// tracking gate entirely. The [xmin, xmax) histogram gate
// lives inside PixelHistogramImpl::fill.
const auto n_sigma = n_sigma_.load(std::memory_order_relaxed);
if (n_sigma <= AxisType{0.0}) {
// Fill without pedestal tracking.
for (const auto &frame : *images) {
const auto cols = frame.shape(1);
for (int local_row = 0; local_row < local_rows;
++local_row) {
const auto row =
static_cast<ssize_t>(first_row + local_row);
for (ssize_t col = 0; col < cols; ++col) {
const FrameType raw = frame(row, col);
const AxisType val =
static_cast<AxisType>(raw) -
static_cast<AxisType>(my_pedestal.mean(
static_cast<uint32_t>(local_row),
static_cast<uint32_t>(col)));
my_hist.fill_unchecked(local_row,
static_cast<int>(col), val);
}
}
}
break;
} else {
// Do pedestal tracking. Duplicated code for clean hot path.
constexpr std::size_t pixel_tile_size = 256;
const auto local_pixels = static_cast<std::size_t>(local_rows) *
static_cast<std::size_t>(cols_);
const auto global_pixel_begin =
static_cast<std::size_t>(first_row) *
static_cast<std::size_t>(cols_);
const auto &frames = *images;
auto &my_std = partial_std_[thread_id];
for (const auto &frame : *images) {
const auto cols = frame.shape(1);
for (int local_row = 0; local_row < local_rows;
++local_row) {
const auto row =
static_cast<ssize_t>(first_row + local_row);
for (ssize_t col = 0; col < cols; ++col) {
const FrameType raw = frame(row, col);
// Compile the shared traversal once for each mode. if constexpr
// removes the tracking gate entirely from the histogram-only hot
// path, while keeping the runtime mode decision outside the
// per-pixel loops.
const auto fill_tiles = [&](auto tracking,
const AxisType *std_data) {
constexpr bool track_pedestal = decltype(tracking)::value;
// Each worker retains exclusive ownership of its row shard.
// Pixel tiling keeps a bounded part of that shard's
// [pixel x bin] storage hot while input remains contiguous.
for (std::size_t tile_begin = 0; tile_begin < local_pixels;
tile_begin += pixel_tile_size) {
const auto tile_end =
std::min(tile_begin + pixel_tile_size, local_pixels);
for (const auto &frame : frames) {
const auto *input =
frame.data() + global_pixel_begin + tile_begin;
for (std::size_t local_pixel = tile_begin;
local_pixel < tile_end; ++local_pixel, ++input) {
const FrameType raw = *input;
const AxisType val =
static_cast<AxisType>(raw) -
static_cast<AxisType>(my_pedestal.mean(
static_cast<uint32_t>(local_row),
static_cast<uint32_t>(col)));
my_hist.fill_unchecked(local_row,
static_cast<int>(col), val);
const AxisType sigma = my_std(local_row, col);
if (sigma > AxisType{0.0} &&
std::abs(static_cast<AxisType>(val)) <
n_sigma * sigma) {
my_pedestal.template push<FrameType>(
static_cast<uint32_t>(local_row),
static_cast<uint32_t>(col), raw);
my_pedestal.mean_unchecked(
static_cast<ssize_t>(local_pixel));
my_hist.fill_flat_unchecked(local_pixel, val);
if constexpr (track_pedestal) {
const AxisType sigma = std_data[local_pixel];
if (sigma > AxisType{0.0} &&
std::abs(val) < n_sigma * sigma) {
my_pedestal.push_ema_unchecked<FrameType>(
local_pixel, raw);
}
}
}
}
}
break;
};
if (n_sigma <= AxisType{0.0}) {
fill_tiles(std::false_type{}, nullptr);
} else {
// Frames stay chronological for every pixel. Accepted raw
// values update that pixel's sums and cached mean immediately;
// pixels are otherwise independent.
fill_tiles(std::true_type{}, partial_std_[thread_id].data());
}
break;
}
}
@@ -395,6 +394,11 @@ void PedestalTrackingPixelHistogram::fill_async(NDArray<FrameType, 2> &&image) {
"constructor shape");
}
// ProducerConsumerQueue is SPSC. Serialising this short producer-side
// operation also prevents fill_async from racing fill_from_file's direct
// batch dispatch.
std::lock_guard<std::mutex> ingestion_lock(ingestion_mutex_);
// SPSC backpressure: spin with a short sleep until a slot frees up.
// The std::move only consumes `image` on the iteration that succeeds
// (placement-new inside write() runs only when the slot is free).
@@ -467,35 +471,42 @@ PedestalTrackingPixelHistogram::bin_edges() const {
}
void PedestalTrackingPixelHistogram::fill_from_file(
const std::filesystem::path &fname, ssize_t max_frames, bool verbose) {
const std::filesystem::path &fname, ssize_t max_frames, bool verbose,
std::size_t reader_threads, std::size_t reader_chunk_size) {
constexpr std::size_t progress_interval = 66;
auto last = std::chrono::steady_clock::now();
const auto completed_start =
completed_async_fills_.load(std::memory_order_acquire);
std::size_t last_reported = 0;
const auto completed_for_this_file = [&]() {
return completed_async_fills_.load(std::memory_order_acquire) -
completed_start;
};
const auto wait_for_completed = [&](std::size_t target) {
while (completed_for_this_file() < target) {
std::this_thread::sleep_for(async_wait_);
}
};
if (max_frames < -1) {
throw std::invalid_argument(
"PedestalTrackingPixelHistogram max_frames must be -1 or "
"non-negative");
}
File f(fname);
// check that row col matches constructor
if (f.rows() != static_cast<size_t>(rows_) ||
f.cols() != static_cast<size_t>(cols_)) {
// Preserve the old max_frames behaviour: values beyond EOF are clamped
// rather than rejected by MultiThreadedFileReader.
const auto source_total_frames = File(fname).total_frames();
const auto n_frames = max_frames == -1
? source_total_frames
: std::min(static_cast<std::size_t>(max_frames),
source_total_frames);
experimental::MultiThreadedFileReader reader(fname, reader_threads,
reader_chunk_size, n_frames);
if (reader.rows() != static_cast<size_t>(rows_) ||
reader.cols() != static_cast<size_t>(cols_)) {
throw std::invalid_argument("PedestalTrackingPixelHistogram: Frame in "
"file {} has shape ({}, {}) does not match "
"constructor shape");
}
if (reader.dtype() != Dtype::UINT16 ||
reader.bytes_per_frame() != static_cast<std::size_t>(rows_) *
static_cast<std::size_t>(cols_) *
sizeof(FrameType)) {
throw std::invalid_argument(
"PedestalTrackingPixelHistogram requires uint16 file frames");
}
const ssize_t total_frames = f.total_frames();
const ssize_t n_frames =
max_frames == -1 ? total_frames : std::min(max_frames, total_frames);
const auto print_progress = [&](std::size_t done) {
const auto now = std::chrono::steady_clock::now();
const double dt = std::chrono::duration<double>(now - last).count();
@@ -503,30 +514,89 @@ void PedestalTrackingPixelHistogram::fill_from_file(
const double fps =
dt > 0.0 ? static_cast<double>(done_in_interval) / dt : 0.0;
fmt::print(
"\rProgress: {}/{} ({:.1f}%) {:.1f} FPS ", done, n_frames,
100.0 * static_cast<double>(done) / static_cast<double>(n_frames),
fps);
fmt::print("\rProgress: {}/{} ({:.1f}%) {:.1f} FPS ", done,
n_frames,
n_frames == 0 ? 100.0
: 100.0 * static_cast<double>(done) /
static_cast<double>(n_frames),
fps);
std::fflush(stdout);
last = now;
last_reported = done;
};
for (ssize_t i = 0; i < n_frames; ++i) {
aare::NDArray<uint16_t> frame({rows_, cols_});
f.read_into(reinterpret_cast<std::byte *>(frame.data()));
fill_async(std::move(frame));
using ReadBuffer = NDArray<FrameType, 3>;
const auto read_into_buffer = [&reader](ReadBuffer &buffer) {
const auto frame_count = reader.next_read_frames();
const auto frames_read = reader.read_into(buffer.buffer());
if (frames_read != frame_count) {
throw std::runtime_error(
"MultiThreadedFileReader returned an incomplete batch");
}
return frames_read;
};
if (verbose && (i + 1) % progress_interval == 0) {
wait_for_completed(static_cast<std::size_t>(i + 1));
print_progress(static_cast<std::size_t>(i + 1));
const auto fill_batch = [this](ReadBuffer &buffer,
std::size_t frame_count) {
std::vector<NDView<FrameType, 2>> views;
views.reserve(frame_count);
auto batch_view = buffer.view();
for (std::size_t i = 0; i < frame_count; ++i) {
views.push_back(batch_view(i));
}
std::lock_guard<std::mutex> fill_lock(fill_mutex_);
dispatch_fill_batch_(views);
};
// Exclude other queue producers for the duration. Drain anything already
// submitted before bypassing the coordinator with direct batch dispatch.
std::lock_guard<std::mutex> ingestion_lock(ingestion_mutex_);
flush();
std::size_t completed = 0;
if (n_frames != 0) {
// Allocate the maximum wave size once per buffer. The last read may
// contain fewer frames; its returned count limits the views passed to
// the histogram workers, leaving the unused tail untouched.
const auto buffer_capacity = reader.next_read_frames();
std::array<ReadBuffer, 2> buffers{
ReadBuffer({static_cast<ssize_t>(buffer_capacity), rows_, cols_}),
ReadBuffer({static_cast<ssize_t>(buffer_capacity), rows_, cols_})};
std::size_t current_index = 0;
auto current_frames = read_into_buffer(buffers[current_index]);
while (true) {
const bool has_next = completed + current_frames < n_frames;
const std::size_t next_index = current_index ^ std::size_t{1};
// MultiThreadedFileReader::read_into is blocking. Run the next
// read on a dedicated prefetch task while the current batch is
// processed by the histogram worker pool.
std::future<std::size_t> next;
if (has_next) {
next = std::async(std::launch::async, read_into_buffer,
std::ref(buffers[next_index]));
}
fill_batch(buffers[current_index], current_frames);
completed += current_frames;
if (verbose && (completed - last_reported >= progress_interval ||
completed == n_frames)) {
print_progress(completed);
}
if (!has_next) {
break;
}
current_frames = next.get();
current_index = next_index;
}
}
flush();
if (verbose) {
const auto done = completed_for_this_file();
if (done > last_reported) {
print_progress(done);
if (completed > last_reported || n_frames == 0) {
print_progress(completed);
}
fmt::print("\n\n");
std::fflush(stdout);
src/hist/PedestalTrackingPixelHistogram.test.cpp
+289
View File
@@ -0,0 +1,289 @@
// SPDX-License-Identifier: MPL-2.0
#include "aare/hist/PedestalTrackingPixelHistogram.hpp"
#include "aare/File.hpp"
#include "aare/Frame.hpp"
#include "aare/NumpyFile.hpp"
#include <catch2/catch_test_macros.hpp>
#include <algorithm>
#include <chrono>
#include <cstdint>
#include <filesystem>
#include <functional>
#include <stdexcept>
#include <string>
#include <utility>
using aare::FileConfig;
using aare::Frame;
using aare::NumpyFile;
using aare::PedestalTrackingPixelHistogram;
namespace {
class TemporaryHistogramFile {
public:
using Generator =
std::function<std::uint16_t(std::size_t, std::size_t, std::size_t)>;
TemporaryHistogramFile(std::size_t rows = 2, std::size_t cols = 3,
std::size_t frames = 10, Generator generator = {}) {
const auto unique =
std::chrono::steady_clock::now().time_since_epoch().count();
path_ = std::filesystem::temp_directory_path() /
("aare-pedestal-hist-" + std::to_string(unique) + ".npy");
FileConfig config;
config.dtype = aare::Dtype::UINT16;
config.rows = rows;
config.cols = cols;
NumpyFile file(path_, "w", config);
for (std::size_t frame_index = 0; frame_index < frames; ++frame_index) {
Frame frame(rows, cols, config.dtype);
auto image = frame.view<std::uint16_t>();
for (ssize_t row = 0; row < image.shape(0); ++row) {
for (ssize_t col = 0; col < image.shape(1); ++col) {
image(row, col) =
generator ? generator(frame_index,
static_cast<std::size_t>(row),
static_cast<std::size_t>(col))
: static_cast<std::uint16_t>(frame_index);
}
}
file.write(frame);
}
}
~TemporaryHistogramFile() { std::filesystem::remove(path_); }
TemporaryHistogramFile(const TemporaryHistogramFile &) = delete;
TemporaryHistogramFile &operator=(const TemporaryHistogramFile &) = delete;
const std::filesystem::path &path() const { return path_; }
private:
std::filesystem::path path_;
};
} // namespace
TEST_CASE("Pedestal tracking histogram fills ordered multi-reader batches",
"[PedestalTrackingPixelHistogram]") {
TemporaryHistogramFile file;
PedestalTrackingPixelHistogram histogram(2, 3, 10, 0.0f, 10.0f, 2, 4, 0.0f);
// Two reader workers claiming two frames each produces a full batch of
// four followed by a partial batch of three.
histogram.fill_from_file(file.path(), 7, false, 2, 2);
const auto values = histogram.values();
for (ssize_t row = 0; row < 2; ++row) {
for (ssize_t col = 0; col < 3; ++col) {
for (ssize_t bin = 0; bin < 10; ++bin) {
CHECK(values(row, col, bin) == (bin < 7 ? 1 : 0));
}
}
}
}
TEST_CASE("Pedestal tracking file fill handles limits and reader options",
"[PedestalTrackingPixelHistogram]") {
TemporaryHistogramFile file;
PedestalTrackingPixelHistogram histogram(2, 3, 10, 0.0f, 10.0f, 2, 4, 0.0f);
CHECK_NOTHROW(histogram.fill_from_file(file.path(), 0, false, 2, 2));
CHECK_THROWS_AS(histogram.fill_from_file(file.path(), -2, false, 2, 2),
std::invalid_argument);
CHECK_THROWS_AS(histogram.fill_from_file(file.path(), -1, false, 0, 2),
std::invalid_argument);
CHECK_THROWS_AS(histogram.fill_from_file(file.path(), -1, false, 2, 0),
std::invalid_argument);
// Preserve the previous API's clamp-at-EOF behaviour.
histogram.fill_from_file(file.path(), 100, false, 2, 3);
const auto values = histogram.values();
for (ssize_t row = 0; row < 2; ++row) {
for (ssize_t col = 0; col < 3; ++col) {
for (ssize_t bin = 0; bin < 10; ++bin) {
CHECK(values(row, col, bin) == 1);
}
}
}
}
TEST_CASE("Pedestal tracking file fill validates frame metadata",
"[PedestalTrackingPixelHistogram]") {
TemporaryHistogramFile wrong_shape(3, 3, 1);
PedestalTrackingPixelHistogram histogram(2, 3, 10, 0.0f, 10.0f, 1, 4, 0.0f);
CHECK_THROWS_AS(
histogram.fill_from_file(wrong_shape.path(), -1, false, 2, 1),
std::invalid_argument);
}
TEST_CASE("Histogram-only fill crosses a pixel tile without changing pedestal",
"[PedestalTrackingPixelHistogram]") {
constexpr int rows = 1;
constexpr int cols = 513;
constexpr std::size_t pedestal_samples = 1000;
const auto baseline = [](std::size_t col) {
return static_cast<std::uint16_t>(100 + col % 17);
};
TemporaryHistogramFile file(
rows, cols, 2,
[baseline](std::size_t frame, std::size_t, std::size_t col) {
return static_cast<std::uint16_t>(baseline(col) + frame + 1);
});
// Negative n_sigma exercises the complete histogram-only condition; zero
// is covered by the ordered batch tests above.
PedestalTrackingPixelHistogram histogram(rows, cols, 4, 0.0f, 4.0f, 1, 6,
-1.0f);
for (std::size_t seed = 0; seed < pedestal_samples; ++seed) {
aare::NDArray<std::uint16_t, 2> frame({rows, cols});
const int offset = seed % 2 == 0 ? -1 : 1;
for (ssize_t col = 0; col < cols; ++col) {
frame(0, col) = static_cast<std::uint16_t>(
static_cast<int>(baseline(static_cast<std::size_t>(col))) +
offset);
}
histogram.push_pedestal_no_update(frame.view());
}
histogram.update_mean();
const auto mean_before = histogram.pedestal_mean();
// One two-frame batch traverses the 512-pixel tile and its one-pixel tail.
histogram.fill_from_file(file.path(), -1, false, 1, 2);
const auto values = histogram.values();
for (ssize_t col = 0; col < cols; ++col) {
for (ssize_t bin = 0; bin < 4; ++bin) {
CHECK(values(0, col, bin) == ((bin == 1 || bin == 2) ? 1 : 0));
}
}
const auto mean_after = histogram.pedestal_mean();
CHECK(
std::equal(mean_before.begin(), mean_before.end(), mean_after.begin()));
}
TEST_CASE("Pedestal tracking threshold is strict and rejects zero sigma",
"[PedestalTrackingPixelHistogram]") {
constexpr std::size_t pedestal_samples = 1000;
PedestalTrackingPixelHistogram histogram(1, 3, 8, -4.0f, 4.0f, 1, 4, 2.0f);
// Pixels 0 and 1 have mean 100 and population sigma 1. Pixel 2 has the
// same mean and zero sigma.
for (std::size_t seed = 0; seed < pedestal_samples; ++seed) {
aare::NDArray<std::uint16_t, 2> frame({1, 3});
const auto noisy = static_cast<std::uint16_t>(seed % 2 == 0 ? 99 : 101);
frame(0, 0) = noisy;
frame(0, 1) = noisy;
frame(0, 2) = 100;
histogram.push_pedestal_no_update(frame.view());
}
histogram.update_mean();
const auto mean_before = histogram.pedestal_mean();
REQUIRE(mean_before(0, 0) == 100.0f);
REQUIRE(mean_before(0, 1) == 100.0f);
REQUIRE(mean_before(0, 2) == 100.0f);
aare::NDArray<std::uint16_t, 2> frame({1, 3});
frame(0, 0) = 102; // residual == 2 * sigma: excluded
frame(0, 1) = 101; // residual < 2 * sigma: included
frame(0, 2) = 101; // sigma == 0: excluded
histogram.fill_async(std::move(frame));
histogram.flush();
const auto mean_after = histogram.pedestal_mean();
CHECK(mean_after(0, 0) == 100.0f);
CHECK(mean_after(0, 1) == static_cast<float>(100001.0 * (1.0 / 1000.0)));
CHECK(mean_after(0, 2) == 100.0f);
// Histogramming uses each residual before a possible EMA update.
const auto values = histogram.values();
for (ssize_t bin = 0; bin < 8; ++bin) {
CHECK(values(0, 0, bin) == (bin == 6 ? 1 : 0));
CHECK(values(0, 1, bin) == (bin == 5 ? 1 : 0));
CHECK(values(0, 2, bin) == (bin == 5 ? 1 : 0));
}
}
TEST_CASE("Tiled pedestal tracking matches chronological single-frame fills",
"[PedestalTrackingPixelHistogram]") {
constexpr int rows = 5;
constexpr int cols = 257;
constexpr std::size_t frames = 40;
constexpr int bins = 32;
constexpr float xmin = -16.0f;
constexpr float xmax = 16.0f;
const auto baseline = [](std::size_t row, std::size_t col) {
return static_cast<std::uint16_t>(100 + (row + col) % 5);
};
const auto sample = [baseline](std::size_t frame, std::size_t row,
std::size_t col) {
const int offsets[] = {1, -1, 8, 0};
return static_cast<std::uint16_t>(static_cast<int>(baseline(row, col)) +
offsets[frame % 4]);
};
TemporaryHistogramFile file(rows, cols, frames, sample);
PedestalTrackingPixelHistogram tiled(rows, cols, bins, xmin, xmax, 2, 16,
2.0f);
PedestalTrackingPixelHistogram chronological(rows, cols, bins, xmin, xmax,
2, 16, 2.0f);
// Seed a non-zero cached standard deviation. Five rows split over two
// workers make both row shards cross the 512-pixel tile boundary. Forty
// data frames exercise chronological processing across a sizable batch.
for (std::size_t seed = 0; seed < 1000; ++seed) {
aare::NDArray<std::uint16_t, 2> frame({rows, cols});
const int offset = seed % 2 == 0 ? -2 : 2;
for (ssize_t row = 0; row < rows; ++row) {
for (ssize_t col = 0; col < cols; ++col) {
frame(row, col) = static_cast<std::uint16_t>(
static_cast<int>(baseline(row, col)) + offset);
}
}
tiled.push_pedestal_no_update(frame.view());
chronological.push_pedestal_no_update(frame.view());
}
tiled.update_mean();
chronological.update_mean();
// One forty-frame reader wave exercises tiled batch processing.
tiled.fill_from_file(file.path(), -1, false, 2, 20);
// The reference path establishes the same result one chronological frame
// at a time, without any cross-frame traversal reordering.
for (std::size_t frame_index = 0; frame_index < frames; ++frame_index) {
aare::NDArray<std::uint16_t, 2> frame({rows, cols});
for (ssize_t row = 0; row < rows; ++row) {
for (ssize_t col = 0; col < cols; ++col) {
frame(row, col) =
sample(frame_index, static_cast<std::size_t>(row),
static_cast<std::size_t>(col));
}
}
chronological.fill_async(std::move(frame));
chronological.flush();
}
const auto tiled_values = tiled.values();
const auto chronological_values = chronological.values();
REQUIRE(tiled_values.shape() == chronological_values.shape());
CHECK(std::equal(tiled_values.begin(), tiled_values.end(),
chronological_values.begin()));
const auto tiled_mean = tiled.pedestal_mean();
const auto chronological_mean = chronological.pedestal_mean();
REQUIRE(tiled_mean.shape() == chronological_mean.shape());
CHECK(std::equal(tiled_mean.begin(), tiled_mean.end(),
chronological_mean.begin()));
}
src/hist/PixelHistogramImpl.test.cpp
+11
View File
@@ -72,6 +72,17 @@ TEST_CASE("Fill a small histogram from an NDArray") {
REQUIRE(v(1, 1, 4) == 1);
}
TEST_CASE("Flat pixel filling uses pixel-major histogram storage") {
aare::PixelHistogramImpl<float, uint16_t> hist(2, 3, 4, 0.0f, 4.0f);
hist.fill_flat_unchecked(0, 0.5f);
hist.fill_flat_unchecked(4, 2.5f);
const auto values = hist.view();
CHECK(values(0, 0, 0) == 1);
CHECK(values(1, 1, 2) == 1);
}
TEST_CASE("Check that pixel histogram does not overflow") {
int rows = 1;
int cols = 1;