detectors/aare
25
0 Fork 0
mirror of https://github.com/slsdetectorgroup/aare.git synced 2026-10-30 12:31:28 +01:00
Code Issues Actions Packages Releases Activity

Merge branch 'main' into fix/spsc-move-assign-leak
Build on RHEL9 / build (push) Successful in 3m12s
Details
Build on RHEL8 / build (push) Successful in 3m44s
Details
Run tests using data on local RHEL8 / build (push) Successful in 4m26s
Details

Browse Source
This commit is contained in:
Erik Fröjdh authored and GitHub committed 2026-09-23 16:17:54 +02:00
2 parents e412e72e5b 1398a74877
commit e5cd15dc28
277 files changed
+25657 -4532
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

.gitea/workflows/cmake_build.yml
+11 -10
View File
@@ -32,21 +32,22 @@ jobs:
run: |
sudo apt-get update
sudo apt-get -y install cmake gcc g++
- name: Get conda
uses: conda-incubator/setup-miniconda@v3
with:
python-version: ${{ matrix.python-version }}
environment-file: etc/dev-env.yml
miniforge-version: latest
channels: conda-forge
conda-remove-defaults: "true"
sudo apt-get -y install python3.12 python3.12-dev python3.12-venv python3-pip
sudo apt-get -y install doxygen
python3.12 -m venv venv
source venv/bin/activate
pip install --upgrade pip
pip install breathe
pip install sphinx_rtd_theme sphinx
pip install numpy
pip install furo
- name: Build library
run: |
source venv/bin/activate
mkdir build
cd build
cmake .. -DAARE_SYSTEM_LIBRARIES=ON -DAARE_DOCS=ON
cmake .. -DAARE_PYTHON_BINDINGS=ON -DAARE_DOCS=ON -DAARE_WARNINGS_AS_ERRORS=ON
make -j 2
make docs
.gitea/workflows/rh8-data-tests-local.yml
+47
View File
@@ -0,0 +1,47 @@
name: Run tests using data on local RHEL8
on:
push: # pull_request only works if pull_request in gitea
workflow_dispatch:
permissions:
contents: read
jobs:
build:
runs-on: "detectors-software-RH8"
steps:
- uses: actions/checkout@v4
- name: Clone aare-test-data repository
run: |
git lfs install
git clone https://gitea.psi.ch/detectors/aare-test-data.git ../aare-test-data
cd ../aare-test-data
git lfs pull
- name: Build library
run: |
source /home/gitea_runner/.bashrc
conda activate aare_test
mkdir build && cd build
cmake .. -DAARE_PYTHON_BINDINGS=ON -DAARE_TESTS=ON
make -j 4
- name: C++ unit tests
working-directory: ${{github.workspace}}/build
env:
AARE_TEST_DATA: ${{github.workspace}}/../aare-test-data
run: ./run_tests [.with-data] # TODO: should we run all tests?
- name: Python unit tests
working-directory: ${{github.workspace}}/build
env:
AARE_TEST_DATA: ${{github.workspace}}/../aare-test-data
run: |
source /home/gitea_runner/.bashrc
conda activate aare_test
python -m pytest ${{github.workspace}}/python/tests/ --with-data # runs all tests
.gitea/workflows/rh8-local.yml
+30
View File
@@ -0,0 +1,30 @@
name: Build on local RHEL8
on:
push:
branches:
- main
workflow_dispatch:
permissions:
contents: read
jobs:
build:
runs-on: "detectors-software-RH8"
steps:
- uses: actions/checkout@v4
- name: Build library
run: |
source /home/gitea_runner/.bashrc
conda activate det
mkdir build && cd build
conda activate det
cmake .. -DAARE_PYTHON_BINDINGS=ON
make -j 2
- name: Deploy to NFS update server
if: gitea.ref == 'refs/heads/main'
run: |
sftp -r gitea_runner@mpc2935:/ <<< $'put build/aare'
.gitea/workflows/rh8-native.yml
+15 -4
View File
@@ -11,7 +11,7 @@ jobs:
build:
runs-on: "ubuntu-latest"
container:
image: gitea.psi.ch/images/rhel8-developer-gitea-actions
image: gitea.psi.ch/detectors/rhel8-detectors-dev
steps:
# workaround until actions/checkout@v4 is available for RH8
# - uses: actions/checkout@v4
@@ -21,9 +21,15 @@ jobs:
git clone https://${{secrets.GITHUB_TOKEN}}@gitea.psi.ch/${{ github.repository }}.git --branch=${{ github.ref_name }} .
- name: Install dependencies
- name: Install Python dependencies
run: |
dnf install -y cmake python3.12 python3.12-devel python3.12-pip
python3.12 -m pip install --upgrade pip
python3.12 -m pip install pytest
python3.12 -m pip install numpy
python3.12 -m pip install pytest-check
python3.12 -m pip install matplotlib
python3.12 -m pip install boost-histogram
- name: Build library
run: |
@@ -33,4 +39,9 @@ jobs:
- name: C++ unit tests
working-directory: ${{gitea.workspace}}/build
run: ctest
run: ctest
- name: Python unit tests
working-directory: ${{gitea.workspace}}/build
run: |
python3.12 -m pytest ${{gitea.workspace}}/python/tests/
.gitea/workflows/rh9-native.yml
+14 -4
View File
@@ -11,14 +11,19 @@ jobs:
build:
runs-on: "ubuntu-latest"
container:
image: gitea.psi.ch/images/rhel9-developer-gitea-actions
image: gitea.psi.ch/detectors/rhel9-detectors-dev
steps:
- uses: actions/checkout@v4
- name: Install dependencies
- name: Install Python dependencies
run: |
dnf install -y cmake python3.12 python3.12-devel python3.12-pip
python3.12 -m pip install --upgrade pip
python3.12 -m pip install pytest
python3.12 -m pip install numpy
python3.12 -m pip install pytest-check
python3.12 -m pip install matplotlib
python3.12 -m pip install boost-histogram
- name: Build library
run: |
@@ -28,4 +33,9 @@ jobs:
- name: C++ unit tests
working-directory: ${{gitea.workspace}}/build
run: ctest
run: ctest
- name: Python unit tests
working-directory: ${{gitea.workspace}}/build
run: |
python3.12 -m pytest ${{gitea.workspace}}/python/tests/
.github/workflows/build_conda.yml
+10 -7
View File
@@ -1,16 +1,20 @@
name: Build pkgs and deploy if on main
name: Build conda pkgs
on:
push:
branches:
- developer
- main
pull_request:
branches:
- main
jobs:
build:
strategy:
fail-fast: false
matrix:
platform: [ubuntu-latest, ] # macos-12, windows-2019]
platform: [ubuntu-latest, macos-latest]
python-version: ["3.12",]
runs-on: ${{ matrix.platform }}
@@ -21,10 +25,10 @@ jobs:
shell: "bash -l {0}"
steps:
- uses: actions/checkout@v4
- uses: actions/checkout@v6
- name: Get conda
uses: conda-incubator/setup-miniconda@v3
uses: conda-incubator/setup-miniconda@v4
with:
python-version: ${{ matrix.python-version }}
environment-file: etc/dev-env.yml
@@ -32,10 +36,9 @@ jobs:
channels: conda-forge
conda-remove-defaults: "true"
- name: Disable upload
run: conda config --set anaconda_upload no
- name: Build
run: conda build conda-recipe
run: conda-build conda-recipe
.github/workflows/build_and_deploy_conda.yml → .github/workflows/build_conda_and_deploy_on_main.yml
+4 -4
View File
@@ -10,7 +10,7 @@ jobs:
strategy:
fail-fast: false
matrix:
platform: [ubuntu-latest, ] # macos-12, windows-2019]
platform: [ubuntu-latest, macos-latest] # macos-12, windows-2019]
python-version: ["3.12",]
runs-on: ${{ matrix.platform }}
@@ -21,10 +21,10 @@ jobs:
shell: "bash -l {0}"
steps:
- uses: actions/checkout@v4
- uses: actions/checkout@v6
- name: Get conda
uses: conda-incubator/setup-miniconda@v3
uses: conda-incubator/setup-miniconda@v4
with:
python-version: ${{ matrix.python-version }}
environment-file: etc/dev-env.yml
@@ -38,5 +38,5 @@ jobs:
- name: Build
env:
CONDA_TOKEN: ${{ secrets.CONDA_TOKEN }}
run: conda build conda-recipe --user slsdetectorgroup --token ${CONDA_TOKEN}
run: conda-build conda-recipe --user slsdetectorgroup --token ${CONDA_TOKEN}
.github/workflows/build_wheel.yml
+2 -2
View File
@@ -17,13 +17,13 @@ jobs:
runs-on: ${{ matrix.os }}
strategy:
matrix:
os: [ubuntu-latest,]
os: [ubuntu-latest,macos-latest]
steps:
- uses: actions/checkout@v4
- name: Build wheels
run: pipx run cibuildwheel==2.23.0
run: pipx run cibuildwheel==3.4.0
- uses: actions/upload-artifact@v4
with:
.github/workflows/build_docs.yml → .github/workflows/build_with_docs.yml
+18 -4
View File
@@ -2,11 +2,15 @@ name: Build the package using cmake then documentation
on:
workflow_dispatch:
push:
pull_request:
release:
types:
- published
env:
# Customize the CMake build type here (Release, Debug, RelWithDebInfo, etc.)
BUILD_TYPE: Debug
permissions:
contents: read
@@ -18,7 +22,7 @@ jobs:
strategy:
fail-fast: false
matrix:
platform: [ubuntu-latest, ]
platform: [ubuntu-latest, macos-latest]
python-version: ["3.12",]
runs-on: ${{ matrix.platform }}
@@ -39,15 +43,25 @@ jobs:
channels: conda-forge
conda-remove-defaults: "true"
- name: Build library
- name: Build library and docs
run: |
mkdir build
cd build
cmake .. -DAARE_SYSTEM_LIBRARIES=ON -DAARE_PYTHON_BINDINGS=ON -DAARE_DOCS=ON
make -j 2
cmake .. -DCMAKE_BUILD_TYPE=${{env.BUILD_TYPE}} -DAARE_SYSTEM_LIBRARIES=ON -DAARE_PYTHON_BINDINGS=ON -DAARE_DOCS=ON -DAARE_TESTS=ON -DAARE_WARNINGS_AS_ERRORS=ON
make -j 4
make docs
- name: C++ unit tests
working-directory: ${{github.workspace}}/build
run: ctest -C ${{env.BUILD_TYPE}} -j4
- name: Python unit tests
working-directory: ${{github.workspace}}/build
run: python -m pytest ${{github.workspace}}/python/tests/
- name: Upload static files as artifact
if: matrix.platform == 'ubuntu-latest'
id: deployment
uses: actions/upload-pages-artifact@v3
with:
.github/workflows/pre-commit.yml
+12
View File
@@ -0,0 +1,12 @@
name: pre-commit
on:
pull_request:
jobs:
pre-commit:
name: pre-commit
runs-on: ubuntu-latest
steps:
- uses: actions/checkout@v6
- uses: pre-commit/action@v3.0.1
.gitignore
+188 -16
View File
@@ -1,25 +1,197 @@
install/
.cproject
.project
bin/
.settings
*.aux
*.log
*.out
*.toc
# C++ .gitignore template (https://github.com/github/gitignore/blob/main/C%2B%2B.gitignore)
# Prerequisites
*.d
# Compiled Object files
*.slo
*.lo
*.o
*.obj
# Precompiled Headers
*.gch
*.pch
# Linker files
*.ilk
# Debugger Files
*.pdb
# Compiled Dynamic libraries
*.so
.*
*.dylib
*.dll
*.so.*
# Fortran module files
*.mod
*.smod
# Compiled Static libraries
*.lai
*.la
*.a
*.lib
# Executables
*.exe
*.out
*.app
# Build directories
build/
RELEASE.txt
Build/
build-*/
# CMake generated files
CMakeFiles/
CMakeCache.txt
cmake_install.cmake
Makefile
install_manifest.txt
compile_commands.json
# Temporary files
*.tmp
*.log
*.bak
*.swp
# vcpkg
vcpkg_installed/
# debug information files
*.dwo
# test output & cache
Testing/
.cache/
ctbDict.cpp
ctbDict.h
# Python .gitignore template (https://github.com/github/gitignore/blob/main/Python.gitignore)
# Byte-compiled / optimized / DLL files
__pycache__/
*.py[codz]
*$py.class
wheelhouse/
# C extensions
*.so
# Distribution / packaging
.Python
build/
develop-eggs/
dist/
downloads/
eggs/
.eggs/
lib/
lib64/
parts/
sdist/
var/
wheels/
share/python-wheels/
*.egg-info/
.installed.cfg
*.egg
MANIFEST
*.pyc
*/__pycache__/*
# PyInstaller
# Usually these files are written by a python script from a template
# before PyInstaller builds the exe, so as to inject date/other infos into it.
*.manifest
*.spec
# Installer logs
pip-log.txt
pip-delete-this-directory.txt
# Unit test / coverage reports
htmlcov/
.tox/
.nox/
.coverage
.coverage.*
.cache
nosetests.xml
coverage.xml
*.cover
*.py.cover
.hypothesis/
.pytest_cache/
cover/
# Translations
*.mo
*.pot
.*
# Django stuff:
*.log
local_settings.py
db.sqlite3
db.sqlite3-journal
# Flask stuff:
instance/
.webassets-cache
# Scrapy stuff:
.scrapy
# Sphinx documentation
docs/_build/
# PyBuilder
.pybuilder/
target/
# Jupyter Notebook
.ipynb_checkpoints
# IPython
profile_default/
ipython_config.py
# PEP 582; used by e.g. github.com/David-OConnor/pyflow and github.com/pdm-project/pdm
__pypackages__/
# Environments
.env
.envrc
.venv
env/
venv/
ENV/
env.bak/
venv.bak/
# mypy
.mypy_cache/
.dmypy.json
dmypy.json
# Pyre type checker
.pyre/
# pytype static type analyzer
.pytype/
# Cython debug symbols
cython_debug/
# Visual Studio Code
# Visual Studio Code specific template is maintained in a separate VisualStudioCode.gitignore
# that can be found at https://github.com/github/gitignore/blob/main/Global/VisualStudioCode.gitignore
# and can be added to the global gitignore or merged into this file. However, if you prefer,
# you could uncomment the following to ignore the entire vscode folder
.vscode/
# Ruff stuff:
.ruff_cache/
# user defined
wheelhouse/
.pre-commit-config.yaml
+29
View File
@@ -0,0 +1,29 @@
fail_fast: true
repos:
- repo: https://github.com/pre-commit/pre-commit-hooks
rev: v6.0.0
hooks:
- id: check-added-large-files
- id: check-yaml
exclude: conda-recipe/
- id: check-toml
- id: check-json
- id: mixed-line-ending
- repo: https://github.com/crate-ci/typos
rev: v1.45.0
hooks:
- id: typos
files: docs
args: [] # empty, to remove write-changes from the default arguments.
- repo: https://github.com/pre-commit/mirrors-clang-format
rev: v22.1.3
hooks:
- id: clang-format
name: clang-format .cpp and .hpp files
- repo: https://github.com/zultron/cmake-format-precommit
rev: v0.6.14
hooks:
- id: cmake-format
AGENTS.md
+307
View File
@@ -0,0 +1,307 @@
# 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`,
`AARE_WARNINGS_AS_ERRORS`, `AARE_VERBOSE` (raises the compile-time log
level), and `AARE_TUNE_LOCAL` (`-march=native`, not portable). Reconfigure an
existing build directory instead of creating alternate in-tree build layouts
unless isolation is needed; `build/CMakeCache.txt` records the options and
interpreter it was configured with.
`-Werror=return-type` is always enabled. New code should also compile cleanly
under the project's `-Wall -Wextra -pedantic -Wshadow -Wold-style-cast
-Wdouble-promotion` flags, which `AARE_WARNINGS_AS_ERRORS` turns into errors.
The build copies `python/aare/*.py` and the compiled `_aare` extension into
`build/aare/`, so `PYTHONPATH=build` makes the package importable. The
extension only imports into the interpreter CMake was configured with, which
is normally the Conda or venv Python active at configure time. Activate that
environment, or resolve the interpreter from the cache:
```bash
PY=$(grep -E '^_?Python_EXECUTABLE' build/CMakeCache.txt | head -1 | cut -d= -f2)
PYTHONPATH="$PWD/build" "$PY" -m pytest python/tests
```
Do not fall back to a different system Python.
With `AARE_BENCHMARKS=ON`, Google Benchmark programs build as
`build/run_benchmarks` and `build/fit_benchmark`.
## Architecture
### Core containers
- `NDArray<T, Ndim>` owns contiguous C-order memory; `NDView<T, Ndim>` is a
non-owning view with the same shape and stride API. Both derive from
`ArrayExpr` (expression templates for element-wise math). Shapes are
`std::array<ssize_t, Ndim>`. `NDView<T>` converts to `NDView<const T>`.
- `Frame` is a byte buffer plus rows, cols, and `Dtype`, the lowest common
denominator returned by file readers. `Dtype` is a runtime type tag built
from a bitdepth, a NumPy descriptor, or `typeid`.
- `defs.hpp` holds detector enums (`DetectorType`, `TimingMode`,
`FrameDiscardPolicy`, `ReadoutMode`), the on-disk `DetectorHeader`, chip
constants (`Moench04`, `Matterhorn02`, ...), the `LOCATION` macro used in
error messages, and `AARE_ASSERT` (active only with `AARE_CUSTOM_ASSERT`).
- `logger.hpp` provides `LOG(logLEVEL) << ...`; levels above `AARE_LOG_LEVEL`
(set by CMake) are compiled out.
### File reading
`File` is the generic RAII facade. It holds a `unique_ptr<FileInterface>` and
dispatches on file extension to one of the concrete readers:
- `RawFile` reads slsDetector `.raw` data with a `_master_N.json` (or legacy
`.raw`) master file. It owns a `RawMasterFile` (master JSON, older raw-text
masters, scan parameters, ROIs, disabled UDP ports) and a two-dimensional
table of `RawSubFile` indexed `[roi][module]`. Each `RawSubFile` handles one
module's sequence of `_dN_fM_` data files and per-frame `DetectorHeader`s.
Module placement into the assembled image goes through `DetectorGeometry`,
`ModuleGeometry`, and `ROIGeometry` (module gaps, ROIs, UDP-port halves).
Frame numbers from the detector are distinct from frame indices in the
file; the API uses both and names them accordingly.
- `NumpyFile` (`.npy`), `JungfrauDataFile` (legacy `.dat`), and `CtbRawFile`
(chip test board; raw ADC samples decoded through `decode.hpp` and the
`PixelMap` generators, exposed in Python as `aare.transform` callables).
- `experimental::MultiThreadedFileReader` opens one `File` per worker and
fills a caller-owned buffer in file order.
Some chip-specific decoding (Moench03/05 pixel maps, ADC SAR transforms)
lives in `python/aare/transform.py` on top of C++ generators.
### Pedestal and cluster finding
- `Pedestal<T>` and `FastPedestal<T>` maintain per-pixel running mean and
standard deviation, always accumulating in `double`. `FastPedestal` must
see `min_pedestal_samples` frames before `ready()`.
- `ClusterFinder<ClusterType, FRAME_TYPE, PEDESTAL_TYPE>` is single-threaded:
push pedestal frames, then `find_clusters(view, frame_number)` appends to an
internal `ClusterVector` that callers take with `steal_clusters()`.
Thresholds are `n_sigma * std` and are recomputed by `update_threshold()`.
- `ClusterFinderMT` wraps N `ClusterFinder`s, one per thread, each with a
`CircularFifo<FrameRef>` input queue backed by a preallocated `FramePool`
and a `ProducerConsumerQueue` output. A collector thread funnels all
outputs into a single sink queue. Consumers of that sink run their own
thread: `ClusterCollector` (accumulate in memory, `steal_clusters()`) or
`ClusterFileSink` (write the binary cluster file). The `stop()` and
`sync()` ordering matters; read the existing tests before changing
lifecycle code.
- `Cluster<T, SizeX, SizeY, CoordType>` is a POD (x, y, data[]) and
`ClusterVector<Cluster<...>>` is a frame-tagged `std::vector` of them.
`ClusterFile` reads and writes the binary format (`int32 frame_number`,
`uint32 n_clusters`, raw cluster bytes). The format carries no shape or
dtype metadata, so the caller must supply matching template arguments.
`CalculateEta.hpp` and `Interpolator` operate on `ClusterVector`s;
`algorithm.hpp` holds small 1-D search helpers.
- `hist/PixelHistogram` and `PedestalTrackingPixelHistogram` build per-pixel
histograms across a frame series, multithreaded by row bands.
### Fitting
`Fit.hpp`, `FitModel<Model>`, and `Models.hpp` (`Gaussian`, `Pol1`, `Pol2`,
`RisingScurve`, `FallingScurve`, `GaussianErfcPlateau`,
`GaussianChargeSharing`, ...) wrap Minuit2. Template bodies and explicit
instantiations live in `src/Fit.cpp`; `FitModelImpl` is a pimpl so Minuit2
headers never leak into the public API. Minuit2 is a private,
`BUILD_INTERFACE`-only dependency of `aare_core`.
### Python layer
- `python/src/module.cpp` registers everything. Templated classes are
instantiated per value type and cluster size through the
`DEFINE_CLUSTER_BINDINGS` and `DEFINE_BINDINGS_CLUSTERFINDER` macros,
producing names such as `ClusterFinder_Cluster3x3i`, `Pedestal_d`, and
`FastPedestal_i16`. Supported cluster sizes are 2x2, 3x3, 5x5, 7x7, and
9x9; value types are `i`, `f`, `d`, and `i16` (3x3 only).
`module_config.hpp` fixes the module-wide pedestal type (`pd_type`).
- `python/aare/factory.py` maps NumPy dtypes to those suffixes. The facades
in `python/aare/*.py` are thin factory functions or subclasses that select
the right binding and add Pythonic conveniences (iteration, context
managers, `ScanParameters`).
- `np_helper.hpp` provides `return_image_data` (moves an `NDArray` into a
NumPy array with a capsule deleter) and `make_view_1d/2d/3d` (views into a
NumPy buffer).
- `aare.experimental` and `_aare.experimental` hold unstable APIs.
## 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 test name, or an individual pytest
file or test:
```bash
build/run_tests "[tag]"
build/run_tests "Exact test case name"
PYTHONPATH="$PWD/build" python -m pytest python/tests/test_example.py
PYTHONPATH="$PWD/build" python -m pytest python/tests/test_example.py -k name
```
Large detector test files live outside this repository in
`https://gitea.psi.ch/detectors/aare-test-data.git` (Git LFS). 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`
and add them to `TestSources` in the root `CMakeLists.txt`.
- Mark C++ tests requiring external files with `[.with-data]`.
- Mark Python tests requiring external files with `@pytest.mark.withdata` and
use the `test_data_path` fixture from `python/tests/conftest.py`.
- Prefer synthesizing small input files over depending on external data.
`tests/raw_file_helpers.hpp` builds throwaway raw files in C++, and many
Python tests write small raw files into `tmp_path`.
- `tests/test_config.hpp.in` provides `test_data_path()`;
`tests/test_macros.hpp` and `tests/friend_test.hpp` (`TEST_CASE_PRIVATE`
with `FRIEND_TEST`) let a test reach private members.
- 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
The project uses pre-commit. Install the git hook once per clone with
`pre-commit install` so it runs on every commit. The hooks in
`.pre-commit-config.yaml` are clang-format on C++, cmake-format on CMake
files, `typos` on `docs/` (project words are whitelisted in `typos.toml`), and
basic YAML, TOML, JSON, and large-file checks. CI runs the same hooks on pull
requests.
Before handing off a broad change, run the relevant checks:
```bash
pre-commit install # once per clone
pre-commit run --all-files # what CI runs
cmake --build build --target check-format # diff-based clang-format check
cmake --build build --target format-files # apply clang-format to all .cpp/.hpp
cmake --build build --target clang-tidy # uses build/compile_commands.json
```
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. C++
pages live at `docs/src/<Class>.rst`; Python pages live under
`docs/src/python/<area>/`.
- 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.
CLAUDE.md
+8
View File
@@ -0,0 +1,8 @@
# CLAUDE.md
This file provides guidance to Claude Code (claude.ai/code) when working with code in this repository.
All project guidance lives in `AGENTS.md`: layout, build and test commands,
architecture, conventions, formatting, documentation, and release notes.
Read it and follow it. Keep new guidance there rather than here so that every
agent and contributor sees the same instructions.
CMakeLists.txt
+345 -318
View File
@@ -1,10 +1,11 @@
# SPDX-License-Identifier: MPL-2.0
cmake_minimum_required(VERSION 3.15)
project(aare
DESCRIPTION "Data processing library for PSI detectors"
HOMEPAGE_URL "https://github.com/slsdetectorgroup/aare"
LANGUAGES C CXX
)
project(
aare
DESCRIPTION "Data processing library for PSI detectors"
HOMEPAGE_URL "https://github.com/slsdetectorgroup/aare"
LANGUAGES C CXX)
# Read VERSION file into project version
set(VERSION_FILE "${CMAKE_CURRENT_SOURCE_DIR}/VERSION")
@@ -12,37 +13,34 @@ file(READ "${VERSION_FILE}" VERSION_CONTENT)
string(STRIP "${VERSION_CONTENT}" PROJECT_VERSION_STRING)
set(PROJECT_VERSION ${PROJECT_VERSION_STRING})
set(CMAKE_CXX_STANDARD 17)
set(CMAKE_CXX_STANDARD 17)
set(CMAKE_CXX_STANDARD_REQUIRED ON)
set(CMAKE_CXX_EXTENSIONS OFF)
execute_process(
COMMAND git log -1 --format=%h
WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}
OUTPUT_VARIABLE GIT_HASH
OUTPUT_STRIP_TRAILING_WHITESPACE
)
COMMAND git log -1 --format=%h
WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}
OUTPUT_VARIABLE GIT_HASH
OUTPUT_STRIP_TRAILING_WHITESPACE)
message(STATUS "Building from git hash: ${GIT_HASH}")
if (${CMAKE_VERSION} VERSION_GREATER "3.24")
cmake_policy(SET CMP0135 NEW) #Fetch content download timestamp
if(${CMAKE_VERSION} VERSION_GREATER "3.24")
cmake_policy(SET CMP0135 NEW) # Fetch content download timestamp
endif()
cmake_policy(SET CMP0079 NEW)
include(GNUInstallDirs)
include(FetchContent)
#Set default build type if none was specified
# Set default build type if none was specified
include(cmake/helpers.cmake)
default_build_type("Release")
set_std_fs_lib()
message(STATUS "Extra linking to fs lib:${STD_FS_LIB}")
set(CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/cmake" ${CMAKE_MODULE_PATH})
# General options
option(AARE_PYTHON_BINDINGS "Build python bindings" OFF)
option(AARE_TESTS "Build tests" OFF)
@@ -52,8 +50,15 @@ option(AARE_IN_GITHUB_ACTIONS "Running in Github Actions" OFF)
option(AARE_DOCS "Build documentation" OFF)
option(AARE_VERBOSE "Verbose output" OFF)
option(AARE_CUSTOM_ASSERT "Use custom assert" OFF)
option(AARE_INSTALL_PYTHONEXT "Install the python extension in the install tree under CMAKE_INSTALL_PREFIX/aare/" OFF)
option(
AARE_INSTALL_PYTHONEXT
"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)
# Configure which of the dependencies to use FetchContent for
option(AARE_FETCH_FMT "Use FetchContent to download fmt" ON)
@@ -61,155 +66,160 @@ 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"
OFF)
#Convenience option to use system libraries only (no FetchContent)
# Convenience option to use system libraries only (no FetchContent)
option(AARE_SYSTEM_LIBRARIES "Use system libraries" OFF)
if(AARE_SYSTEM_LIBRARIES)
message(STATUS "Build using system libraries")
set(AARE_FETCH_FMT OFF CACHE BOOL "Disabled FetchContent for FMT" FORCE)
set(AARE_FETCH_PYBIND11 OFF CACHE BOOL "Disabled FetchContent for pybind11" FORCE)
set(AARE_FETCH_CATCH OFF CACHE BOOL "Disabled FetchContent for catch2" FORCE)
set(AARE_FETCH_JSON OFF CACHE BOOL "Disabled FetchContent for nlohmann::json" FORCE)
set(AARE_FETCH_ZMQ OFF CACHE BOOL "Disabled FetchContent for libzmq" FORCE)
# Still fetch lmfit when setting AARE_SYSTEM_LIBRARIES since this is not available
# on conda-forge
message(STATUS "Build using system libraries")
set(AARE_FETCH_FMT
OFF
CACHE BOOL "Disabled FetchContent for FMT" FORCE)
set(AARE_FETCH_PYBIND11
OFF
CACHE BOOL "Disabled FetchContent for pybind11" FORCE)
set(AARE_FETCH_CATCH
OFF
CACHE BOOL "Disabled FetchContent for catch2" FORCE)
set(AARE_FETCH_JSON
OFF
CACHE BOOL "Disabled FetchContent for nlohmann::json" FORCE)
set(AARE_FETCH_ZMQ
OFF
CACHE BOOL "Disabled FetchContent for libzmq" FORCE)
# Still fetch Minuit2 when setting AARE_SYSTEM_LIBRARIES since it is not
# available on conda-forge
endif()
if(AARE_BENCHMARKS)
add_subdirectory(benchmarks)
add_subdirectory(benchmarks)
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)
if(AARE_FETCH_MINUIT2)
# We are building Minuit2 from sources.
# 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()
FetchContent_Declare(
Minuit2
GIT_REPOSITORY https://github.com/GooFit/Minuit2.git
GIT_TAG master
UPDATE_DISCONNECTED 1)
# Disable Minuit2 extras we don't need
set(minuit2_mpi
OFF
CACHE BOOL "")
set(minuit2_omp
OFF
CACHE BOOL "")
set(BUILD_TESTING
OFF
CACHE BOOL "")
set(MINUIT2_INSTALL
OFF
CACHE BOOL "")
FetchContent_MakeAvailable(Minuit2)
#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 "")
set_property(TARGET Minuit2Math PROPERTY POSITION_INDEPENDENT_CODE ON)
set_property(TARGET Minuit2 PROPERTY POSITION_INDEPENDENT_CODE ON)
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)
add_library(aare::Minuit2 ALIAS Minuit2)
else()
find_package(lmfit REQUIRED)
find_package(Minuit2 CONFIG QUIET)
if(NOT Minuit2_FOUND)
message(
STATUS
"Stand alone Minuit2 not found. Trying to find it as a part of ROOT")
find_package(
ROOT
COMPONENTS Minuit2
REQUIRED)
add_library(aare::Minuit2 ALIAS ROOT::Minuit2)
message(STATUS "Found Minuit2 as a part of ROOT: ${ROOT_DIR}")
endif()
endif()
if(AARE_FETCH_ZMQ)
# Fetchcontent_Declare is deprecated need to find a way to update this
# for now setting the policy to old is enough
if (${CMAKE_VERSION} VERSION_GREATER_EQUAL "3.30")
cmake_policy(SET CMP0169 OLD)
endif()
set(ZMQ_PATCH_COMMAND git apply ${CMAKE_CURRENT_SOURCE_DIR}/patches/libzmq_cmake_version.patch)
FetchContent_Declare(
libzmq
GIT_REPOSITORY https://github.com/zeromq/libzmq.git
GIT_TAG v4.3.4
PATCH_COMMAND ${ZMQ_PATCH_COMMAND}
UPDATE_DISCONNECTED 1
)
# Disable unwanted options from libzmq
set(BUILD_TESTS OFF CACHE BOOL "Switch off libzmq test build")
set(BUILD_SHARED OFF CACHE BOOL "Switch off libzmq shared libs")
set(WITH_PERF_TOOL OFF CACHE BOOL "")
set(ENABLE_CPACK OFF CACHE BOOL "")
set(ENABLE_CLANG OFF CACHE BOOL "")
set(ENABLE_CURVE OFF CACHE BOOL "")
set(ENABLE_DRAFTS OFF CACHE BOOL "")
# Fetchcontent_Declare is deprecated need to find a way to update this for now
# setting the policy to old is enough
if(${CMAKE_VERSION} VERSION_GREATER_EQUAL "3.30")
cmake_policy(SET CMP0169 OLD)
endif()
set(ZMQ_PATCH_COMMAND
git apply ${CMAKE_CURRENT_SOURCE_DIR}/patches/libzmq_cmake_version.patch)
FetchContent_Declare(
libzmq
GIT_REPOSITORY https://github.com/zeromq/libzmq.git
GIT_TAG v4.3.4
PATCH_COMMAND ${ZMQ_PATCH_COMMAND}
UPDATE_DISCONNECTED 1)
# Disable unwanted options from libzmq
set(BUILD_TESTS
OFF
CACHE BOOL "Switch off libzmq test build")
set(BUILD_SHARED
OFF
CACHE BOOL "Switch off libzmq shared libs")
set(WITH_PERF_TOOL
OFF
CACHE BOOL "")
set(ENABLE_CPACK
OFF
CACHE BOOL "")
set(ENABLE_CLANG
OFF
CACHE BOOL "")
set(ENABLE_CURVE
OFF
CACHE BOOL "")
set(ENABLE_DRAFTS
OFF
CACHE BOOL "")
# TODO! Verify that this is what we want to do in aare
# Using GetProperties and Populate to be able to exclude zmq
# from install (not possible with FetchContent_MakeAvailable(libzmq))
FetchContent_GetProperties(libzmq)
if(NOT libzmq_POPULATED)
FetchContent_Populate(libzmq)
add_subdirectory(${libzmq_SOURCE_DIR} ${libzmq_BINARY_DIR} EXCLUDE_FROM_ALL)
endif()
# TODO! Verify that this is what we want to do in aare Using GetProperties and
# Populate to be able to exclude zmq from install (not possible with
# FetchContent_MakeAvailable(libzmq))
FetchContent_GetProperties(libzmq)
if(NOT libzmq_POPULATED)
FetchContent_Populate(libzmq)
add_subdirectory(${libzmq_SOURCE_DIR} ${libzmq_BINARY_DIR} EXCLUDE_FROM_ALL)
endif()
else()
find_package(ZeroMQ 4 REQUIRED)
find_package(ZeroMQ 4 REQUIRED)
endif()
if (AARE_FETCH_FMT)
set(FMT_TEST OFF CACHE INTERNAL "disabling fmt tests")
FetchContent_Declare(
fmt
GIT_REPOSITORY https://github.com/fmtlib/fmt.git
GIT_TAG 10.2.1
GIT_PROGRESS TRUE
USES_TERMINAL_DOWNLOAD TRUE
)
set(FMT_INSTALL ON CACHE BOOL "")
# set(FMT_CMAKE_DIR "")
FetchContent_MakeAvailable(fmt)
set_property(TARGET fmt PROPERTY POSITION_INDEPENDENT_CODE ON)
install(TARGETS fmt
EXPORT ${project}-targets
LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR}
ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR}
RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR}
INCLUDES DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}
)
if(AARE_FETCH_FMT)
set(FMT_TEST
OFF
CACHE INTERNAL "disabling fmt tests")
FetchContent_Declare(
fmt
GIT_REPOSITORY https://github.com/fmtlib/fmt.git
GIT_TAG 12.1.0
GIT_PROGRESS TRUE
USES_TERMINAL_DOWNLOAD TRUE)
set(FMT_INSTALL
ON
CACHE BOOL "")
FetchContent_MakeAvailable(fmt)
set_property(TARGET fmt PROPERTY POSITION_INDEPENDENT_CODE ON)
else()
find_package(fmt 6 REQUIRED)
find_package(fmt 8 REQUIRED)
endif()
if (AARE_FETCH_JSON)
FetchContent_Declare(
json
URL https://github.com/nlohmann/json/releases/download/v3.11.3/json.tar.xz
)
set(JSON_Install ON CACHE BOOL "")
FetchContent_MakeAvailable(json)
set(NLOHMANN_JSON_TARGET_NAME nlohmann_json)
install(
TARGETS nlohmann_json
EXPORT "${TARGETS_EXPORT_NAME}"
)
message(STATUS "target: ${NLOHMANN_JSON_TARGET_NAME}")
if(AARE_FETCH_JSON)
FetchContent_Declare(
json
URL https://github.com/nlohmann/json/releases/download/v3.11.3/json.tar.xz)
set(JSON_Install
ON
CACHE BOOL "")
FetchContent_MakeAvailable(json)
else()
find_package(nlohmann_json 3.11.3 REQUIRED)
find_package(nlohmann_json 3.11.3 REQUIRED)
endif()
include(GNUInstallDirs)
@@ -223,7 +233,6 @@ endif()
string(TOUPPER "${PROJECT_NAME}" PROJECT_NAME_UPPER)
string(TOLOWER "${PROJECT_NAME}" PROJECT_NAME_LOWER)
# Set targets export name (used by slsDetectorPackage and dependencies)
set(TARGETS_EXPORT_NAME "${PROJECT_NAME_LOWER}-targets")
set(namespace "aare::")
@@ -232,106 +241,106 @@ set(CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/cmake" ${CMAKE_MODULE_PATH})
# Check if project is being used directly or via add_subdirectory
set(AARE_MASTER_PROJECT OFF)
if (CMAKE_CURRENT_SOURCE_DIR STREQUAL CMAKE_SOURCE_DIR)
set(AARE_MASTER_PROJECT ON)
if(CMAKE_CURRENT_SOURCE_DIR STREQUAL CMAKE_SOURCE_DIR)
set(AARE_MASTER_PROJECT ON)
endif()
add_library(aare_compiler_flags INTERFACE)
target_compile_features(aare_compiler_flags INTERFACE cxx_std_17)
if(AARE_PYTHON_BINDINGS)
add_subdirectory(python)
endif()
#################
# ##############################################################################
# MSVC specific #
#################
# ##############################################################################
if(MSVC)
add_compile_definitions(AARE_MSVC)
if(CMAKE_BUILD_TYPE STREQUAL "Release")
add_compile_definitions(AARE_MSVC)
if(CMAKE_BUILD_TYPE STREQUAL "Release")
message(STATUS "Release build")
target_compile_options(aare_compiler_flags INTERFACE /O2)
else()
else()
message(STATUS "Debug build")
target_compile_options(
aare_compiler_flags
INTERFACE
/Od
/Zi
/MDd
/D_ITERATOR_DEBUG_LEVEL=2
)
target_link_options(
aare_compiler_flags
INTERFACE
/DEBUG:FULL
)
endif()
target_compile_options(
aare_compiler_flags
INTERFACE
/w # disable warnings
)
target_compile_options(aare_compiler_flags
INTERFACE /Od /Zi /MDd /D_ITERATOR_DEBUG_LEVEL=2)
target_link_options(aare_compiler_flags INTERFACE /DEBUG:FULL)
endif()
target_compile_options(aare_compiler_flags INTERFACE /w # disable warnings
)
else()
######################
# GCC/Clang specific #
######################
# ############################################################################
# GCC/Clang specific #
# ############################################################################
if(CMAKE_BUILD_TYPE STREQUAL "Release")
if(CMAKE_BUILD_TYPE STREQUAL "Release")
message(STATUS "Release build")
target_compile_options(aare_compiler_flags INTERFACE -O3)
else()
else()
message(STATUS "Debug build")
endif()
# Common flags for GCC and Clang
target_compile_options(
aare_compiler_flags
INTERFACE -Wall
-Wextra
-pedantic
-Wshadow
-Wold-style-cast
-Wnon-virtual-dtor
-Woverloaded-virtual
-Wdouble-promotion
-Wformat=2
-Wredundant-decls
-Wvla
-Wdouble-promotion
-Werror=return-type # important can cause segfault in optimzed
# builds
)
if(AARE_WARNINGS_AS_ERRORS)
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)
add_subdirectory(python)
endif()
# Common flags for GCC and Clang
target_compile_options(
aare_compiler_flags
INTERFACE
-Wall
-Wextra
-pedantic
-Wshadow
-Wold-style-cast
-Wnon-virtual-dtor
-Woverloaded-virtual
-Wdouble-promotion
-Wformat=2
-Wredundant-decls
-Wvla
-Wdouble-promotion
-Werror=return-type #important can cause segfault in optimzed builds
)
endif() #GCC/Clang specific
if(AARE_ASAN)
message(STATUS "AddressSanitizer enabled")
target_compile_options(
aare_compiler_flags
INTERFACE
-fsanitize=address,undefined,pointer-compare
-fno-omit-frame-pointer
)
target_link_libraries(
aare_compiler_flags
INTERFACE
-fsanitize=address,undefined,pointer-compare
-fno-omit-frame-pointer
)
message(STATUS "AddressSanitizer enabled")
target_compile_options(
aare_compiler_flags INTERFACE -fsanitize=address,undefined,pointer-compare
-fno-omit-frame-pointer)
target_link_libraries(
aare_compiler_flags INTERFACE -fsanitize=address,undefined,pointer-compare
-fno-omit-frame-pointer)
endif()
if(AARE_TESTS)
enable_testing()
add_subdirectory(tests)
target_compile_definitions(tests PRIVATE AARE_TESTS)
enable_testing()
add_subdirectory(tests)
target_compile_definitions(tests PRIVATE AARE_TESTS)
endif()
###------------------------------------------------------------------------------MAIN LIBRARY
###------------------------------------------------------------------------------------------
# ------------------------------------------------------------------------------MAIN
# LIBRARY
# ------------------------------------------------------------------------------------------
set(PUBLICHEADERS
include/aare/ArrayExpr.hpp
@@ -346,10 +355,19 @@ set(PUBLICHEADERS
include/aare/Dtype.hpp
include/aare/File.hpp
include/aare/Fit.hpp
include/aare/FitModel.hpp
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
include/aare/GainMap.hpp
include/aare/ROI.hpp
include/aare/ROIGeometry.hpp
include/aare/DetectorGeometry.hpp
include/aare/JungfrauDataFile.hpp
include/aare/logger.hpp
@@ -363,16 +381,20 @@ set(PUBLICHEADERS
include/aare/RawMasterFile.hpp
include/aare/RawSubFile.hpp
include/aare/VarClusterFinder.hpp
include/aare/to_string.hpp
include/aare/utils/task.hpp
)
include/aare/utils/ifstream_helpers.hpp
include/aare/utils/math_helpers.hpp)
set(SourceFiles
${CMAKE_CURRENT_SOURCE_DIR}/src/calibration.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/CtbRawFile.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/decode.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/defs.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/ROI.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/ROIGeometry.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/DetectorGeometry.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/hist/PedestalTrackingPixelHistogram.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/Dtype.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/File.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/FilePtr.cpp
@@ -380,155 +402,160 @@ 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
${CMAKE_CURRENT_SOURCE_DIR}/src/RawFile.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/RawMasterFile.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/RawSubFile.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/to_string.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/utils/task.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/utils/ifstream_helpers.cpp
)
${CMAKE_CURRENT_SOURCE_DIR}/src/utils/ifstream_helpers.cpp)
add_library(aare_core STATIC ${SourceFiles})
target_include_directories(aare_core PUBLIC
"$<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/include>"
"$<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>"
)
add_library(aare_core STATIC ${SourceFiles} ${PUBLICHEADERS})
target_include_directories(
aare_core PUBLIC "$<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/include>"
"$<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>")
set(THREADS_PREFER_PTHREAD_FLAG ON)
find_package(Threads REQUIRED)
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>
aare_core
PUBLIC fmt::fmt nlohmann_json::nlohmann_json ${STD_FS_LIB} # from
# helpers.cmake
PRIVATE aare_compiler_flags Threads::Threads $<BUILD_INTERFACE:aare::Minuit2>)
)
target_include_directories(
aare_core SYSTEM
PRIVATE $<TARGET_PROPERTY:aare::Minuit2,INTERFACE_INCLUDE_DIRECTORIES>)
set_property(TARGET aare_core PROPERTY POSITION_INDEPENDENT_CODE ON)
if(AARE_TESTS)
target_compile_definitions(aare_core PRIVATE AARE_TESTS)
target_compile_definitions(aare_core PRIVATE AARE_TESTS)
endif()
if(AARE_VERBOSE)
target_compile_definitions(aare_core PUBLIC AARE_VERBOSE)
target_compile_definitions(aare_core PUBLIC AARE_LOG_LEVEL=aare::logDEBUG5)
target_compile_definitions(aare_core PUBLIC AARE_VERBOSE)
target_compile_definitions(aare_core PUBLIC AARE_LOG_LEVEL=aare::logDEBUG5)
else()
target_compile_definitions(aare_core PUBLIC AARE_LOG_LEVEL=aare::logERROR)
target_compile_definitions(aare_core PUBLIC AARE_LOG_LEVEL=aare::logWARNING)
endif()
if(AARE_CUSTOM_ASSERT)
target_compile_definitions(aare_core PUBLIC AARE_CUSTOM_ASSERT)
target_compile_definitions(aare_core PUBLIC AARE_CUSTOM_ASSERT)
endif()
set_target_properties(aare_core PROPERTIES
ARCHIVE_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}
PUBLIC_HEADER "${PUBLICHEADERS}"
)
if (AARE_PYTHON_BINDINGS)
set_property(TARGET aare_core PROPERTY POSITION_INDEPENDENT_CODE ON)
endif()
set_target_properties(aare_core PROPERTIES ARCHIVE_OUTPUT_DIRECTORY
${CMAKE_BINARY_DIR})
if(AARE_TESTS)
set(TestSources
${CMAKE_CURRENT_SOURCE_DIR}/src/algorithm.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/calibration.test.cpp
${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/Frame.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/DetectorGeometry.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/RawMasterFile.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/NDArray.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/NDView.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/ClusterFinder.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/ClusterVector.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/Cluster.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/CalculateEta.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/ClusterFile.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/ClusterFinderMT.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/Pedestal.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/JungfrauDataFile.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
${CMAKE_CURRENT_SOURCE_DIR}/src/RawSubFile.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/utils/task.test.cpp
)
target_sources(tests PRIVATE ${TestSources} )
set(TestSources
${CMAKE_CURRENT_SOURCE_DIR}/tests/raw_file_helpers.hpp
${CMAKE_CURRENT_SOURCE_DIR}/src/algorithm.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/calibration.test.cpp
${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/FilePtr.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
${CMAKE_CURRENT_SOURCE_DIR}/src/RawMasterFile.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/NDArray.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/NDView.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/ClusterFinder.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/ClusterVector.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/Cluster.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/CalculateEta.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/ClusterFile.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/ClusterFinderMT.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/Pedestal.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/FastPedestal.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
${CMAKE_CURRENT_SOURCE_DIR}/src/RawSubFile.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/ROI.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/utils/task.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/to_string.test.cpp)
target_sources(tests PRIVATE ${TestSources})
endif()
if(AARE_MASTER_PROJECT)
install(TARGETS aare_core aare_compiler_flags
EXPORT "${TARGETS_EXPORT_NAME}"
LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR}
ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR}
PUBLIC_HEADER DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/aare
)
install(
TARGETS aare_core aare_compiler_flags
EXPORT "${TARGETS_EXPORT_NAME}"
LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR}
ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR})
install(
DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/include/aare/
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/aare
FILES_MATCHING
PATTERN "*.hpp")
endif()
set(CMAKE_POSITION_INDEPENDENT_CODE ON)
set(CMAKE_INSTALL_RPATH $ORIGIN)
set(CMAKE_BUILD_WITH_INSTALL_RPATH FALSE)
# #Overall target to link to when using the library
# add_library(aare INTERFACE)
# #Overall target to link to when using the library add_library(aare INTERFACE)
# target_link_libraries(aare INTERFACE aare_core aare_compiler_flags)
# target_include_directories(aare INTERFACE
# $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/include>
# $<INSTALL_INTERFACE:include>
# )
# target_include_directories(aare INTERFACE
# $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/include>
# $<INSTALL_INTERFACE:include> )
# add_subdirectory(examples)
if(AARE_DOCS)
add_subdirectory(docs)
add_subdirectory(docs)
endif()
# custom target to run check formatting with clang-format
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 {} -"
WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}
COMMENT "Checking code formatting with clang-format"
VERBATIM
)
check-format
COMMAND
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)
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\" {}"
WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}
COMMENT "Formatting with clang-format"
VERBATIM
)
format-files
COMMAND
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)
if (AARE_IN_GITHUB_ACTIONS)
message(STATUS "Running in Github Actions")
set(CLANG_TIDY_COMMAND "clang-tidy-17")
if(AARE_IN_GITHUB_ACTIONS)
message(STATUS "Running in Github Actions")
set(CLANG_TIDY_COMMAND "clang-tidy-17")
else()
set(CLANG_TIDY_COMMAND "clang-tidy")
set(CLANG_TIDY_COMMAND "clang-tidy")
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 "${CLANG_TIDY_COMMAND} --config-file=.clang-tidy -p build {}"
WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}
COMMENT "linting with clang-tidy"
VERBATIM
)
clang-tidy
COMMAND
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"
VERBATIM)
if(AARE_MASTER_PROJECT)
set(CMAKE_INSTALL_DIR "share/cmake/${PROJECT_NAME}")
set(PROJECT_LIBRARIES aare-core aare-compiler-flags )
include(cmake/package_config.cmake)
set(CMAKE_INSTALL_DIR "share/cmake/${PROJECT_NAME}")
set(PROJECT_LIBRARIES aare-core aare-compiler-flags)
include(cmake/package_config.cmake)
endif()
LICENSE
+373
View File
@@ -0,0 +1,373 @@
Mozilla Public License Version 2.0
==================================
1. Definitions
--------------
1.1. "Contributor"
means each individual or legal entity that creates, contributes to
the creation of, or owns Covered Software.
1.2. "Contributor Version"
means the combination of the Contributions of others (if any) used
by a Contributor and that particular Contributor's Contribution.
1.3. "Contribution"
means Covered Software of a particular Contributor.
1.4. "Covered Software"
means Source Code Form to which the initial Contributor has attached
the notice in Exhibit A, the Executable Form of such Source Code
Form, and Modifications of such Source Code Form, in each case
including portions thereof.
1.5. "Incompatible With Secondary Licenses"
means
(a) that the initial Contributor has attached the notice described
in Exhibit B to the Covered Software; or
(b) that the Covered Software was made available under the terms of
version 1.1 or earlier of the License, but not also under the
terms of a Secondary License.
1.6. "Executable Form"
means any form of the work other than Source Code Form.
1.7. "Larger Work"
means a work that combines Covered Software with other material, in
a separate file or files, that is not Covered Software.
1.8. "License"
means this document.
1.9. "Licensable"
means having the right to grant, to the maximum extent possible,
whether at the time of the initial grant or subsequently, any and
all of the rights conveyed by this License.
1.10. "Modifications"
means any of the following:
(a) any file in Source Code Form that results from an addition to,
deletion from, or modification of the contents of Covered
Software; or
(b) any new file in Source Code Form that contains any Covered
Software.
1.11. "Patent Claims" of a Contributor
means any patent claim(s), including without limitation, method,
process, and apparatus claims, in any patent Licensable by such
Contributor that would be infringed, but for the grant of the
License, by the making, using, selling, offering for sale, having
made, import, or transfer of either its Contributions or its
Contributor Version.
1.12. "Secondary License"
means either the GNU General Public License, Version 2.0, the GNU
Lesser General Public License, Version 2.1, the GNU Affero General
Public License, Version 3.0, or any later versions of those
licenses.
1.13. "Source Code Form"
means the form of the work preferred for making modifications.
1.14. "You" (or "Your")
means an individual or a legal entity exercising rights under this
License. For legal entities, "You" includes any entity that
controls, is controlled by, or is under common control with You. For
purposes of this definition, "control" means (a) the power, direct
or indirect, to cause the direction or management of such entity,
whether by contract or otherwise, or (b) ownership of more than
fifty percent (50%) of the outstanding shares or beneficial
ownership of such entity.
2. License Grants and Conditions
--------------------------------
2.1. Grants
Each Contributor hereby grants You a world-wide, royalty-free,
non-exclusive license:
(a) under intellectual property rights (other than patent or trademark)
Licensable by such Contributor to use, reproduce, make available,
modify, display, perform, distribute, and otherwise exploit its
Contributions, either on an unmodified basis, with Modifications, or
as part of a Larger Work; and
(b) under Patent Claims of such Contributor to make, use, sell, offer
for sale, have made, import, and otherwise transfer either its
Contributions or its Contributor Version.
2.2. Effective Date
The licenses granted in Section 2.1 with respect to any Contribution
become effective for each Contribution on the date the Contributor first
distributes such Contribution.
2.3. Limitations on Grant Scope
The licenses granted in this Section 2 are the only rights granted under
this License. No additional rights or licenses will be implied from the
distribution or licensing of Covered Software under this License.
Notwithstanding Section 2.1(b) above, no patent license is granted by a
Contributor:
(a) for any code that a Contributor has removed from Covered Software;
or
(b) for infringements caused by: (i) Your and any other third party's
modifications of Covered Software, or (ii) the combination of its
Contributions with other software (except as part of its Contributor
Version); or
(c) under Patent Claims infringed by Covered Software in the absence of
its Contributions.
This License does not grant any rights in the trademarks, service marks,
or logos of any Contributor (except as may be necessary to comply with
the notice requirements in Section 3.4).
2.4. Subsequent Licenses
No Contributor makes additional grants as a result of Your choice to
distribute the Covered Software under a subsequent version of this
License (see Section 10.2) or under the terms of a Secondary License (if
permitted under the terms of Section 3.3).
2.5. Representation
Each Contributor represents that the Contributor believes its
Contributions are its original creation(s) or it has sufficient rights
to grant the rights to its Contributions conveyed by this License.
2.6. Fair Use
This License is not intended to limit any rights You have under
applicable copyright doctrines of fair use, fair dealing, or other
equivalents.
2.7. Conditions
Sections 3.1, 3.2, 3.3, and 3.4 are conditions of the licenses granted
in Section 2.1.
3. Responsibilities
-------------------
3.1. Distribution of Source Form
All distribution of Covered Software in Source Code Form, including any
Modifications that You create or to which You contribute, must be under
the terms of this License. You must inform recipients that the Source
Code Form of the Covered Software is governed by the terms of this
License, and how they can obtain a copy of this License. You may not
attempt to alter or restrict the recipients' rights in the Source Code
Form.
3.2. Distribution of Executable Form
If You distribute Covered Software in Executable Form then:
(a) such Covered Software must also be made available in Source Code
Form, as described in Section 3.1, and You must inform recipients of
the Executable Form how they can obtain a copy of such Source Code
Form by reasonable means in a timely manner, at a charge no more
than the cost of distribution to the recipient; and
(b) You may distribute such Executable Form under the terms of this
License, or sublicense it under different terms, provided that the
license for the Executable Form does not attempt to limit or alter
the recipients' rights in the Source Code Form under this License.
3.3. Distribution of a Larger Work
You may create and distribute a Larger Work under terms of Your choice,
provided that You also comply with the requirements of this License for
the Covered Software. If the Larger Work is a combination of Covered
Software with a work governed by one or more Secondary Licenses, and the
Covered Software is not Incompatible With Secondary Licenses, this
License permits You to additionally distribute such Covered Software
under the terms of such Secondary License(s), so that the recipient of
the Larger Work may, at their option, further distribute the Covered
Software under the terms of either this License or such Secondary
License(s).
3.4. Notices
You may not remove or alter the substance of any license notices
(including copyright notices, patent notices, disclaimers of warranty,
or limitations of liability) contained within the Source Code Form of
the Covered Software, except that You may alter any license notices to
the extent required to remedy known factual inaccuracies.
3.5. Application of Additional Terms
You may choose to offer, and to charge a fee for, warranty, support,
indemnity or liability obligations to one or more recipients of Covered
Software. However, You may do so only on Your own behalf, and not on
behalf of any Contributor. You must make it absolutely clear that any
such warranty, support, indemnity, or liability obligation is offered by
You alone, and You hereby agree to indemnify every Contributor for any
liability incurred by such Contributor as a result of warranty, support,
indemnity or liability terms You offer. You may include additional
disclaimers of warranty and limitations of liability specific to any
jurisdiction.
4. Inability to Comply Due to Statute or Regulation
---------------------------------------------------
If it is impossible for You to comply with any of the terms of this
License with respect to some or all of the Covered Software due to
statute, judicial order, or regulation then You must: (a) comply with
the terms of this License to the maximum extent possible; and (b)
describe the limitations and the code they affect. Such description must
be placed in a text file included with all distributions of the Covered
Software under this License. Except to the extent prohibited by statute
or regulation, such description must be sufficiently detailed for a
recipient of ordinary skill to be able to understand it.
5. Termination
--------------
5.1. The rights granted under this License will terminate automatically
if You fail to comply with any of its terms. However, if You become
compliant, then the rights granted under this License from a particular
Contributor are reinstated (a) provisionally, unless and until such
Contributor explicitly and finally terminates Your grants, and (b) on an
ongoing basis, if such Contributor fails to notify You of the
non-compliance by some reasonable means prior to 60 days after You have
come back into compliance. Moreover, Your grants from a particular
Contributor are reinstated on an ongoing basis if such Contributor
notifies You of the non-compliance by some reasonable means, this is the
first time You have received notice of non-compliance with this License
from such Contributor, and You become compliant prior to 30 days after
Your receipt of the notice.
5.2. If You initiate litigation against any entity by asserting a patent
infringement claim (excluding declaratory judgment actions,
counter-claims, and cross-claims) alleging that a Contributor Version
directly or indirectly infringes any patent, then the rights granted to
You by any and all Contributors for the Covered Software under Section
2.1 of this License shall terminate.
5.3. In the event of termination under Sections 5.1 or 5.2 above, all
end user license agreements (excluding distributors and resellers) which
have been validly granted by You or Your distributors under this License
prior to termination shall survive termination.
************************************************************************
* *
* 6. Disclaimer of Warranty *
* ------------------------- *
* *
* Covered Software is provided under this License on an "as is" *
* basis, without warranty of any kind, either expressed, implied, or *
* statutory, including, without limitation, warranties that the *
* Covered Software is free of defects, merchantable, fit for a *
* particular purpose or non-infringing. The entire risk as to the *
* quality and performance of the Covered Software is with You. *
* Should any Covered Software prove defective in any respect, You *
* (not any Contributor) assume the cost of any necessary servicing, *
* repair, or correction. This disclaimer of warranty constitutes an *
* essential part of this License. No use of any Covered Software is *
* authorized under this License except under this disclaimer. *
* *
************************************************************************
************************************************************************
* *
* 7. Limitation of Liability *
* -------------------------- *
* *
* Under no circumstances and under no legal theory, whether tort *
* (including negligence), contract, or otherwise, shall any *
* Contributor, or anyone who distributes Covered Software as *
* permitted above, be liable to You for any direct, indirect, *
* special, incidental, or consequential damages of any character *
* including, without limitation, damages for lost profits, loss of *
* goodwill, work stoppage, computer failure or malfunction, or any *
* and all other commercial damages or losses, even if such party *
* shall have been informed of the possibility of such damages. This *
* limitation of liability shall not apply to liability for death or *
* personal injury resulting from such party's negligence to the *
* extent applicable law prohibits such limitation. Some *
* jurisdictions do not allow the exclusion or limitation of *
* incidental or consequential damages, so this exclusion and *
* limitation may not apply to You. *
* *
************************************************************************
8. Litigation
-------------
Any litigation relating to this License may be brought only in the
courts of a jurisdiction where the defendant maintains its principal
place of business and such litigation shall be governed by laws of that
jurisdiction, without reference to its conflict-of-law provisions.
Nothing in this Section shall prevent a party's ability to bring
cross-claims or counter-claims.
9. Miscellaneous
----------------
This License represents the complete agreement concerning the subject
matter hereof. If any provision of this License is held to be
unenforceable, such provision shall be reformed only to the extent
necessary to make it enforceable. Any law or regulation which provides
that the language of a contract shall be construed against the drafter
shall not be used to construe this License against a Contributor.
10. Versions of the License
---------------------------
10.1. New Versions
Mozilla Foundation is the license steward. Except as provided in Section
10.3, no one other than the license steward has the right to modify or
publish new versions of this License. Each version will be given a
distinguishing version number.
10.2. Effect of New Versions
You may distribute the Covered Software under the terms of the version
of the License under which You originally received the Covered Software,
or under the terms of any subsequent version published by the license
steward.
10.3. Modified Versions
If you create software not governed by this License, and you want to
create a new license for such software, you may create and use a
modified version of this License if you rename the license and remove
any references to the name of the license steward (except to note that
such modified license differs from this License).
10.4. Distributing Source Code Form that is Incompatible With Secondary
Licenses
If You choose to distribute Source Code Form that is Incompatible With
Secondary Licenses under the terms of this version of the License, the
notice described in Exhibit B of this License must be attached.
Exhibit A - Source Code Form License Notice
-------------------------------------------
This Source Code Form is subject to the terms of the Mozilla Public
License, v. 2.0. If a copy of the MPL was not distributed with this
file, You can obtain one at http://mozilla.org/MPL/2.0/.
If it is not possible or desirable to put the notice in a particular
file, then You may include the notice in a location (such as a LICENSE
file in a relevant directory) where a recipient would be likely to look
for such a notice.
You may add additional accurate notices of copyright ownership.
Exhibit B - "Incompatible With Secondary Licenses" Notice
---------------------------------------------------------
This Source Code Form is "Incompatible With Secondary Licenses", as
defined by the Mozilla Public License, v. 2.0.
README.md
+27 -5
View File
@@ -1,6 +1,14 @@
# aare
Data analysis library for PSI hybrid detectors
## Documentation
Detailed documentation including installation can be found in [Documentation](https://slsdetectorgroup.github.io/aare/)
## License
This project is licensed under the MPL-2.0 license.
See the LICENSE file or https://www.mozilla.org/en-US/MPL/ for details.
## Build and install
@@ -17,27 +25,30 @@ mkdir build
cd build
#configure using cmake
cmake ../aare
cmake ../aare -DAARE_PYTHON_BINDINGS=ON
#build (replace 4 with the number of threads you want to use)
make -j4
```
Now you can use the Python module from your build directory
Now you can use the Python module from your build directory
```python
import aare
f = aare.File('Some/File/I/Want_to_open_master_0.json')
```
To run form other folders either add the path to your conda environment using conda-build or add it to your PYTHONPATH
To run from other folders either add the path to your conda environment using conda-build or add the module to your PYTHONPATH
```bash
export PYTHONPATH=path_to_aare/aare/build:$PYTHONPATH
```
### Install using conda/mamba
```bash
#enable your env first!
conda install aare=2024.10.29.dev0 -c slsdetectorgroup
conda install aare -c slsdetectorgroup # installs latest version
```
### Install to a custom location and use in your project
@@ -68,4 +79,15 @@ make install
```bash
conda build . --variants="{python: [3.11, 3.12, 3.13]}"
```
```
## Developer's guide
We are looking forward to your contributions via pull requests!
If you want to fix an existing bug or propose a new feature:
1. Install `pre-commit` python package and setup it `pre-commit install`
2. Create a new branch with `git branch branch_name`
3. Implement your changes and make a commit (`pre-commit` will check your code automatically)
4. Push your commit and open a pull request if needed
RELEASE.md
+230
View File
@@ -1,5 +1,214 @@
# Release notes
## Next
### New Features:
- Added the Python ``FrameDiscardPolicy`` enum with ``NoDiscard``, ``Discard``,
and ``DiscardPartial``, enabling access to ``RawMasterFile.frame_discard_policy``.
- Added the Python ``Pedestal`` factory with ``dtype`` selection, matching
``FastPedestal`` and defaulting to ``float64`` output.
- Added ``FastPedestal`` in C++ and Python for per-pixel running mean
and population standard deviation. It supports exponentially
- Added ``ClusterFile.frames()`` and ``ClusterFile.chunks()`` in C++ and
Python for iteration from the current file position. Frames preserve empty
frames and frame numbers; chunks support an optional positive size override.
Python's default file iteration continues to yield chunks.
- 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.
- Added string representator in python for Cluster and Eta
- Added roi slice method in python for easy slicing of numpy arrays ``array[roi.slice()]``.
- added context manager for ``aare.RawMasterFile``
### API Changes:
- ``FastPedestal`` variance is now a private ``double`` intermediate. Removed
the C++ ``variance()``/``variance_unchecked()`` APIs and Python ``var()``.
Standard deviation is calculated before conversion to the output type,
avoiding overflow of intermediate variance for ``int16`` output. Negative
variance from floating-point roundoff is clamped to zero.
- ``Pedestal`` now always accumulates sums and sums of squares in ``double``,
like ``FastPedestal``. Mean and standard deviation output types
are unchanged. Removed the C++ ``get_sum()``/``get_sum2()`` getters and
Python ``sum``/``sum2`` properties; internal sums are no longer exposed.
- Removed the public ``Pedestal.variance()`` and ``cached_std()`` APIs and
C++ ``update_std()``. Use ``std()`` to calculate the current population
standard deviation; variance is now an internal implementation detail.
- Added ``ClusterFile::read_frame(ClusterVector&)`` for allocation-reusing C++
reads. It returns ``false`` at a clean end of file. The value-returning C++
overload now returns ``std::optional<ClusterVector<ClusterType>>`` and
Python ``read_frame()`` returns ``None`` at end of file. Incomplete frames
still raise an error.
- Added an explicit boolean conversion to the C++ ``FilePtr`` type.
- Removed the public C++ ``ClusterFile::open`` method. Construct a new
``ClusterFile`` to reopen a file or change its mode.
- ``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.
- ``RawMasterFile::geometry()`` is deprecetad and returns full detector geometry information including module geometry. Use
``RawMasterFile::module_layout()`` to get num_modules in x an y
- ``RawMasterFile::rois()`` always returns a list of rois (no optional). Per default it returns a list of one ROI element spawing the entire detector
- ``TimingMode::Auto`` changed to ``TimingMode::AUTO_TIMING``, ``TimingMode::Trigger`` changed to ``TimingMode::TRIGGER_EXPOSURE``
### Bugfixes:
- Mismatched operators inhibited vectorization in gcc of NDArray math operators
- ``RawFile`` and ``File`` reject raw files with frame padding disabled unless
the frame discard policy is ``discardpartial``. The constructor reports the
master path before opening data subfiles. Legacy ``.raw`` master files now
parse the frame discard policy so unpadded ``discardpartial`` files remain
readable.
- Fixed ``CtbRawFile.read_frame(index)`` and reads after ``seek(index)`` at
subfile boundaries skipping a subfile, returning the wrong frame or raising
``Subfile index out of range``.
- ``Pedestal`` reports mismatched frame shapes with exceptions in all push
overloads, including Debug builds, instead of aborting on assertions.
- Python ``Pedestal`` constructors reject negative dimensions and sample
counts instead of converting them to large unsigned values.
- ``Pedestal`` clamps negative variance from floating-point roundoff to zero,
preventing NaN standard deviations for nearly constant inputs.
- Python ``Pedestal.push()`` now requires C-contiguous ``uint16`` frames
without implicit conversion. Both ``push()`` and ``push_with_threshold()``
validate that frames and thresholds are two-dimensional before constructing
views, preventing incorrect results from unsupported array layouts or ranks.
- Fixed a leaked empty ``ClusterVector`` at the end of Python ``ClusterFile``
iteration. Chunk iteration now rejects a zero chunk size.
- ``ClusterFile::read_clusters`` and Python iteration now report incomplete
frame headers and cluster records instead of treating truncated files as a
clean end of file, with or without ROI or noise filtering.
- ``ClusterFile::write_frame`` now reports incomplete writes instead of
silently continuing with a truncated file.
- Gain-map application now checks the complete cluster footprint, preventing
out-of-bounds access for cluster sizes larger than 3x3.
- ``RawFile`` now derives its frame count from the shortest selected raw
subfile series across all ROIs. Frame-number reads use the same bounds, and
Python ``len(reader)`` returns the adjusted count. A warning is printed when
subfile counts differ or their minimum differs from the recorded master
count. The expected frame count is not used and master metadata is preserved.
Warning-level logging is now enabled in non-verbose builds.
- Raw frame, batch, ROI, and frame-number read errors now include the attempted
frame index and file path in C++ and Python, preserving subfile error details.
Subfile errors omit redundant master-file context and C++ source locations;
out-of-range errors include the available frame count.
Synchronization errors identify the last raw data file for the affected module.
Top-level frame bounds errors state the total frame count and that indices
are zero-based.
- Fixed ``ClusterVector`` move operations to transfer storage instead of
copying every cluster.
- Validate that ``ClusterVector`` masks are one-dimensional, C-contiguous
Boolean arrays.
- Preserve signed ``ClusterVector`` frame numbers when filtering or reducing
cluster dimensions.
- Fixed broken reading of old (pre reordering) Moench03
- Supports reading all timing modes supported in slsDetectorPackage (auto, trigger, gating, burst_trigger, trigger_gating)
## 2026.7.2
### New Features:
- Added a new Minuit2-based fitting framework for ``Gaussian``, ``RisingScurve``, ``FallingScurve``, ``Pol1`` and ``Pol2`` models.
- setter and getter for nSigma for ClusterFinder ``aare.ClusterFinder().nSigma = 2``, ``aare.ClusterFinderMT().set_nSigma(2)``
- mask opeartor for ClusterVector ``masked_clustervector = aare.ClusterVector()(mask)``
- passing pre computed eta values to ``aare.Interpolator.interpolate`` alongside clusters
- Added ``PixelHistogram`` and ``PedestalTrackingPixelHistogram``
- ``aare.transfrom.Matterhorn10Transform`` handles counter artefact in chip. Mind that enabling only one counter or three counters or enabling the wromg two counters e.g. 0,1 will still lead to erreneous data.
- ``aare.transfrom.Matterhorn10Transform`` reshapes data such that first dimension is number of counters
- Added support for len() for files. Returns the number of frames
- Added support for direct subtraction of Pedestal from numpy array
- Added support to read files with disabled udp ports
### Bugfixes:
- Fixed ``split_task(first, last, n_threads)`` so task ranges now correctly respect the ``first`` offset. Previously, non-zero starting indices could generate incorrect subranges.
- Fixed overflow issue causing failed allocations for NDArrays abouve ~2GB
- Fixed libfmt failures due to consteval when building with C++20
## 2026.3.17
### New Features:
- Decoding transceiver data from Matterhorn10 ``transformed_data = aare.transform.Matterhorn10Transform(num_counters=2, dynamic_range=16)(data)``
- Expanding 24 to 32 bit data ``aare._aare.expand24to32bit(data, offset=4)``
- Decoding digital data from Mythen 302 ``transformed_data = aare.transform.Mythen302Transform(offset=4)(data)``
- added ``aare.Interpolator.transform_eta_values``. Function transforms $`\eta`$-values to uniform spatial coordinates. Should only be used for easier debugging.
- New ``to_string``, ``string_to``
- Added exptime and period members to RawMasterFile including decoding
- Removed redundant ``arr.value(ix,iy...)`` on NDArray use ``arr(ix,iy...)``
- Removed Print/Print_some/Print_all form NDArray (operator ``<<`` still works)
- Added const* version of .data()
- reading multiple ROI's supported for aare.
- Use ``aare.RawFile.read_roi(roi_index=0)`` to read a specific ROI for the current frame
- Use ``aare.RawFile.read_rois()`` to read multiple ROIs for the current frame
- Use ``aare.RawFile.read_n_with_roi(num_frames = 2, roi_index = 0)`` to read multiple frames for a specific ROI.
- Note ``read_frame`` and ``read_n`` is not supported for multiple ROI's.
- Building conda/pypi pkgs for python 3.14. Removing 3.11 builds.
### Bugfixes:
- multi threaded cluster finder doesnt drop frames if queues are full
- Round before casting in the cluster finder to avoid biasing clusters by truncating
### 2025.11.21
### New Features:
- Added SPDX-License-Identifier: MPL-2.0 to source files
- Calculate Eta3 supports all cluster types
- interpolation class supports using cross eta3x3 and eta3x3 on full cluster as well as eta2x2 on full cluster
- interpolation class has option to calculate the rosenblatt transform
- reduction operations to reduce Clusters of general size to 2x2 or 3x3 clusters
- `max_sum_2x2` including index of subcluster with highest energy is now available from Python API
- interpolation supports bilinear interpolation of eta values for more fine grained transformed uniform coordinates
- Interpolation is documented
- Added tell to ClusterFile. Returns position in bytes for debugging
### Resolved Features:
- calculate_eta coincides with theoretical definition
### Bugfixes:
- eta calculation assumes correct photon center
- eta transformation to uniform coordinates starts at 0
- Bug in interpolation
- File supports reading new master json file format (multiple ROI's not supported yet)
### API Changes:
- ClusterFinder for 2x2 Cluster disabled
- eta stores corner as enum class cTopLeft, cTopRight, BottomLeft, cBottomRight indicating 2x2 subcluster with largest energy relative to cluster center
- max_sum_2x2 returns corner as index
### 2025.8.22
@@ -40,6 +249,27 @@ Bugfixes:
- Fixed crash when opening raw files with large number of data files
## Download, Documentation & Support
### Download
The Source Code:
https://github.com/slsdetectorgroup/aare
### Documentation
Documentation including installation details:
https://github.com/slsdetectorgroup/aare
### Support
erik.frojdh@psi.ch \
alice.mazzoleni@psi.ch \
dhanya.thattil@psi.ch
VERSION
+1 -1
View File
@@ -1 +1 @@
2025.8.22
2026.7.2
benchmarks/CMakeLists.txt
+25 -11
View File
@@ -1,27 +1,41 @@
# SPDX-License-Identifier: MPL-2.0
include(FetchContent)
FetchContent_Declare(
benchmark
GIT_REPOSITORY https://github.com/google/benchmark.git
GIT_TAG v1.8.3 # Change to the latest version if needed
benchmark
GIT_REPOSITORY https://github.com/google/benchmark.git
GIT_TAG v1.8.3 # Change to the latest version if needed
)
# Ensure Google Benchmark is built correctly
set(BENCHMARK_ENABLE_TESTING OFF CACHE BOOL "" FORCE)
set(BENCHMARK_ENABLE_TESTING
OFF
CACHE BOOL "" FORCE)
FetchContent_MakeAvailable(benchmark)
add_executable(benchmarks)
target_sources(benchmarks PRIVATE ndarray_benchmark.cpp calculateeta_benchmark.cpp reduce_benchmark.cpp)
target_sources(
benchmarks
PRIVATE ndarray_benchmark.cpp ndview_benchmark.cpp calculateeta_benchmark.cpp
reduce_benchmark.cpp clusterfile_benchmark.cpp)
# Link Google Benchmark and other necessary libraries
target_link_libraries(benchmarks PRIVATE benchmark::benchmark aare_core aare_compiler_flags)
target_link_libraries(benchmarks PRIVATE benchmark::benchmark aare_core
aare_compiler_flags)
# Set output properties
set_target_properties(benchmarks PROPERTIES
RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}
OUTPUT_NAME run_benchmarks
)
set_target_properties(
benchmarks PROPERTIES RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}
OUTPUT_NAME run_benchmarks)
add_executable(fit_benchmark fit_benchmark.cpp)
target_link_libraries(fit_benchmark PRIVATE benchmark::benchmark aare_core
aare_compiler_flags)
target_include_directories(
fit_benchmark SYSTEM
PRIVATE $<TARGET_PROPERTY:Minuit2::Minuit2,INTERFACE_INCLUDE_DIRECTORIES>)
set_target_properties(fit_benchmark PROPERTIES RUNTIME_OUTPUT_DIRECTORY
${CMAKE_BINARY_DIR})
benchmarks/calculateeta_benchmark.cpp
+34
View File
@@ -1,3 +1,4 @@
// SPDX-License-Identifier: MPL-2.0
#include "aare/CalculateEta.hpp"
#include "aare/ClusterFile.hpp"
#include <benchmark/benchmark.h>
@@ -8,6 +9,7 @@ class ClusterFixture : public benchmark::Fixture {
public:
Cluster<int, 2, 2> cluster_2x2{};
Cluster<int, 3, 3> cluster_3x3{};
Cluster<int, 4, 4> cluster_4x4{};
private:
using benchmark::Fixture::SetUp;
@@ -26,6 +28,13 @@ class ClusterFixture : public benchmark::Fixture {
cluster_3x3.x = 0;
cluster_3x3.y = 0;
int temp_data3[16] = {1, 2, 3, 4, 5, 6, 7, 8,
9, 10, 11, 12, 13, 14, 15, 16};
std::copy(std::begin(temp_data3), std::end(temp_data3),
std::begin(cluster_4x4.data));
cluster_4x4.x = 0;
cluster_4x4.y = 0;
}
// void TearDown(::benchmark::State& state) {
@@ -67,4 +76,29 @@ BENCHMARK_F(ClusterFixture, CalculateGeneralEtaFor3x3Cluster)
benchmark::DoNotOptimize(eta);
}
}
BENCHMARK_F(ClusterFixture, Calculate2x2Etawithreduction)
(benchmark::State &st) {
for (auto _ : st) {
// This code gets timed
auto reduced_cluster = reduce_to_2x2(cluster_4x4);
Eta2 eta = calculate_eta2(reduced_cluster);
auto reduced_cluster_from_3x3 = reduce_to_2x2(cluster_3x3);
Eta2 eta2 = calculate_eta2(reduced_cluster_from_3x3);
benchmark::DoNotOptimize(eta);
benchmark::DoNotOptimize(eta2);
}
}
BENCHMARK_F(ClusterFixture, Calculate2x2Etawithoutreduction)
(benchmark::State &st) {
for (auto _ : st) {
// This code gets timed
Eta2 eta = calculate_eta2(cluster_4x4);
Eta2 eta2 = calculate_eta2(cluster_3x3);
benchmark::DoNotOptimize(eta);
benchmark::DoNotOptimize(eta2);
}
}
// BENCHMARK_MAIN();
benchmarks/clusterfile_benchmark.cpp
+94
View File
@@ -0,0 +1,94 @@
// SPDX-License-Identifier: MPL-2.0
#include "aare/ClusterFile.hpp"
#include <benchmark/benchmark.h>
#include <chrono>
#include <filesystem>
#include <string>
#include <system_error>
class ClusterFileFixture : public benchmark::Fixture {
using Cluster = aare::Cluster<int32_t, 3, 3>;
using File = aare::ClusterFile<Cluster>;
using Vector = aare::ClusterVector<Cluster>;
std::filesystem::path m_path;
static constexpr int32_t frame_count = 256;
static constexpr size_t chunk_size = 1000;
static void observe(Vector &clusters) {
benchmark::DoNotOptimize(clusters.data());
benchmark::DoNotOptimize(clusters.size());
}
public:
void SetUp(const benchmark::State &state) override {
const auto stamp =
std::chrono::steady_clock::now().time_since_epoch().count();
m_path = std::filesystem::temp_directory_path() /
("aare-cluster-benchmark-" + std::to_string(stamp) + ".clust");
File writer(m_path, chunk_size, "w");
Vector frame(0);
frame.resize(state.range(0));
for (int32_t number = 0; number < frame_count; ++number) {
frame.set_frame_number(number);
writer.write_frame(frame);
}
}
void TearDown(const benchmark::State &) override {
std::error_code error;
std::filesystem::remove(m_path, error);
}
template <bool Frames, bool Iterate> void read(benchmark::State &state) {
for (auto _ : state) {
File reader(m_path, chunk_size);
if constexpr (Iterate) {
auto consume = [](auto range) {
for (auto &clusters : range) {
observe(clusters);
}
};
if constexpr (Frames) {
consume(reader.frames());
} else {
consume(reader.chunks());
}
} else if constexpr (Frames) {
Vector frame(0);
while (reader.read_frame(frame)) {
observe(frame);
}
} else {
while (true) {
auto chunk = reader.read_clusters(chunk_size);
if (chunk.empty()) {
break;
}
observe(chunk);
}
}
}
state.SetBytesProcessed(state.iterations() * frame_count *
state.range(0) * sizeof(Cluster));
}
};
BENCHMARK_DEFINE_F(ClusterFileFixture, ReadFrames)(benchmark::State &state) {
read<true, false>(state);
}
BENCHMARK_DEFINE_F(ClusterFileFixture, IterateFrames)(benchmark::State &state) {
read<true, true>(state);
}
BENCHMARK_DEFINE_F(ClusterFileFixture, ReadChunks)(benchmark::State &state) {
read<false, false>(state);
}
BENCHMARK_DEFINE_F(ClusterFileFixture, IterateChunks)(benchmark::State &state) {
read<false, true>(state);
}
BENCHMARK_REGISTER_F(ClusterFileFixture, ReadFrames)->Arg(64)->Arg(1024);
BENCHMARK_REGISTER_F(ClusterFileFixture, IterateFrames)->Arg(64)->Arg(1024);
BENCHMARK_REGISTER_F(ClusterFileFixture, ReadChunks)->Arg(64)->Arg(1024);
BENCHMARK_REGISTER_F(ClusterFileFixture, IterateChunks)->Arg(64)->Arg(1024);
benchmarks/fit_benchmark.cpp
+139
View File
@@ -0,0 +1,139 @@
// SPDX-License-Identifier: MPL-2.0
#include "aare/Fit.hpp"
#include "aare/FitModel.hpp"
#include "aare/Models.hpp"
#include <benchmark/benchmark.h>
#include <cmath>
#include <random>
#include <string>
#include <vector>
struct TestCase {
std::string name;
double true_A;
double true_mu;
double true_sig;
double noise_frac;
};
static const std::vector<TestCase> &get_test_cases() {
static const std::vector<TestCase> cases = {
{"Clean_signal", 1000.0, 50.0, 5.0, 0.02},
{"Moderate_noise", 1000.0, 50.0, 5.0, 0.10},
{"High_noise", 1000.0, 50.0, 5.0, 0.30},
{"Narrow_peak", 500.0, 25.0, 1.0, 0.05},
{"Wide_peak", 200.0, 100.0, 20.0, 0.05},
{"Off_center_peak", 800.0, -15.0, 3.0, 0.05},
};
return cases;
}
// ----------------------------------------------------------------
// Synthetic data generation (deterministic per test case)
// ----------------------------------------------------------------
static constexpr ssize_t N_POINTS = 100;
static constexpr unsigned SEED = 42;
struct GeneratedData {
aare::NDArray<double, 1> x;
aare::NDArray<double, 1> y;
aare::NDArray<double, 1> y_err;
GeneratedData() : x({N_POINTS}), y({N_POINTS}), y_err({N_POINTS}) {}
};
static GeneratedData generate_gaussian_data(const TestCase &tc) {
GeneratedData d;
double x_min = tc.true_mu - 5.0 * tc.true_sig;
double x_max = tc.true_mu + 5.0 * tc.true_sig;
double dx = (x_max - x_min) / (N_POINTS - 1);
std::mt19937 rng(SEED);
double noise_sigma = tc.noise_frac * tc.true_A;
std::normal_distribution<double> noise(0.0, noise_sigma);
for (ssize_t i = 0; i < N_POINTS; ++i) {
d.x[i] = x_min + i * dx;
double clean = tc.true_A * std::exp(-std::pow(d.x[i] - tc.true_mu, 2) /
(2.0 * std::pow(tc.true_sig, 2)));
d.y[i] = clean + noise(rng);
d.y_err[i] = noise_sigma;
}
return d;
}
static void report_accuracy(benchmark::State &state, const TestCase &tc,
const aare::NDArray<double, 1> &result) {
state.counters["dA"] = result(0) - tc.true_A;
state.counters["dMu"] = result(1) - tc.true_mu;
state.counters["dSig"] = result(2) - tc.true_sig;
}
// ----------
// Benchmarks
// ----------
// 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);
auto xv = data.x.view();
auto yv = data.y.view();
const auto model = aare::FitModel<aare::model::Gaussian>(
/*strategy = */ 0,
/*max_calls = */ 500, // increase for noisy signals
/*tolerance = */ 0.5,
/*compute_errors = */ false);
aare::NDArray<double, 1> result;
for (auto _ : state) {
result = aare::fit_pixel<aare::model::Gaussian>(model, xv, yv);
benchmark::DoNotOptimize(result.data());
}
report_accuracy(state, tc, result);
state.SetLabel(tc.name);
}
// 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);
auto xv = data.x.view();
auto yv = data.y.view();
auto ev = data.y_err.view();
const auto model = aare::FitModel<aare::model::Gaussian>(
0, 500, 0.5, true); // compute_errors = true -> Runs Hesse and provides
// errors on fitted params
aare::NDArray<double, 1> result;
for (auto _ : state) {
result = aare::fit_pixel<aare::model::Gaussian>(model, xv, yv, ev);
benchmark::DoNotOptimize(result.data());
}
// result has 6 elements: [A, mu, sig, err_A, err_mu, err_sig]
report_accuracy(state, tc, result);
// Also report Hesse uncertainties
if (result.size() >= 6) {
state.counters["errA"] = result(3);
state.counters["errMu"] = result(4);
state.counters["errSig"] = result(5);
}
state.SetLabel(tc.name);
}
BENCHMARK(BM_FitGausMinuitGrad)
->DenseRange(0, 5)
->Unit(benchmark::kMicrosecond);
BENCHMARK(BM_FitGausMinuitGradHesse)
->DenseRange(0, 5)
->Unit(benchmark::kMicrosecond);
BENCHMARK_MAIN();
benchmarks/ndarray_benchmark.cpp
+1
View File
@@ -1,3 +1,4 @@
// SPDX-License-Identifier: MPL-2.0
#include "aare/NDArray.hpp"
#include <benchmark/benchmark.h>
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);
benchmarks/reduce_benchmark.cpp
+5 -4
View File
@@ -1,3 +1,4 @@
// SPDX-License-Identifier: MPL-2.0
#include "aare/Cluster.hpp"
#include <benchmark/benchmark.h>
@@ -33,8 +34,8 @@ class ClustersForReduceFixture : public benchmark::Fixture {
};
template <typename T>
Cluster<T, 3, 3, int16_t> reduce_to_3x3(const Cluster<T, 5, 5, int16_t> &c) {
Cluster<T, 3, 3, int16_t> result;
Cluster<T, 3, 3, uint16_t> reduce_to_3x3(const Cluster<T, 5, 5, uint16_t> &c) {
Cluster<T, 3, 3, uint16_t> result;
// Write out the sums in the hope that the compiler can optimize this
std::array<T, 9> sum_3x3_subclusters;
@@ -140,7 +141,7 @@ Cluster<T, 3, 3, int16_t> reduce_to_3x3(const Cluster<T, 5, 5, int16_t> &c) {
BENCHMARK_F(ClustersForReduceFixture, Reduce2x2)(benchmark::State &st) {
for (auto _ : st) {
// This code gets timed
benchmark::DoNotOptimize(reduce_to_2x2<int, 3, 3, int16_t>(
benchmark::DoNotOptimize(reduce_to_2x2<int, 3, 3, uint16_t>(
cluster_3x3)); // make sure compiler evaluates the expression
}
}
@@ -156,7 +157,7 @@ BENCHMARK_F(ClustersForReduceFixture, Reduce3x3)(benchmark::State &st) {
for (auto _ : st) {
// This code gets timed
benchmark::DoNotOptimize(
reduce_to_3x3<int, 5, 5, int16_t>(cluster_5x5));
reduce_to_3x3<int, 5, 5, uint16_t>(cluster_5x5));
}
}
cmake/FindSphinx.cmake
+10 -10
View File
@@ -1,11 +1,11 @@
#Look for an executable called sphinx-build
find_program(SPHINX_EXECUTABLE
NAMES sphinx-build sphinx-build-3.6
DOC "Path to sphinx-build executable")
# Look for an executable called sphinx-build
find_program(
SPHINX_EXECUTABLE
NAMES sphinx-build sphinx-build-3.6
DOC "Path to sphinx-build executable")
include(FindPackageHandleStandardArgs)
#Handle standard arguments to find_package like REQUIRED and QUIET
find_package_handle_standard_args(Sphinx
"Failed to find sphinx-build executable"
SPHINX_EXECUTABLE)
# Handle standard arguments to find_package like REQUIRED and QUIET
find_package_handle_standard_args(
Sphinx "Failed to find sphinx-build executable" SPHINX_EXECUTABLE)
cmake/helpers.cmake
+48 -41
View File
@@ -1,46 +1,53 @@
function(default_build_type val)
if (NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES)
message(STATUS "No build type selected, default to Release")
set(CMAKE_BUILD_TYPE ${val} CACHE STRING "Build type (default ${val})" FORCE)
endif()
if(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES)
message(STATUS "No build type selected, default to Release")
set(CMAKE_BUILD_TYPE
${val}
CACHE STRING "Build type (default ${val})" FORCE)
endif()
endfunction()
function(set_std_fs_lib)
# from pybind11
# Check if we need to add -lstdc++fs or -lc++fs or nothing
if(DEFINED CMAKE_CXX_STANDARD AND CMAKE_CXX_STANDARD LESS 17)
set(STD_FS_NO_LIB_NEEDED TRUE)
elseif(MSVC)
set(STD_FS_NO_LIB_NEEDED TRUE)
else()
file(
WRITE ${CMAKE_CURRENT_BINARY_DIR}/main.cpp
"#include <filesystem>\nint main(int argc, char ** argv) {\n std::filesystem::path p(argv[0]);\n return p.string().length();\n}"
)
try_compile(
STD_FS_NO_LIB_NEEDED ${CMAKE_CURRENT_BINARY_DIR}
SOURCES ${CMAKE_CURRENT_BINARY_DIR}/main.cpp
COMPILE_DEFINITIONS -std=c++17)
try_compile(
STD_FS_NEEDS_STDCXXFS ${CMAKE_CURRENT_BINARY_DIR}
SOURCES ${CMAKE_CURRENT_BINARY_DIR}/main.cpp
COMPILE_DEFINITIONS -std=c++17
LINK_LIBRARIES stdc++fs)
try_compile(
STD_FS_NEEDS_CXXFS ${CMAKE_CURRENT_BINARY_DIR}
SOURCES ${CMAKE_CURRENT_BINARY_DIR}/main.cpp
COMPILE_DEFINITIONS -std=c++17
LINK_LIBRARIES c++fs)
endif()
# from pybind11 Check if we need to add -lstdc++fs or -lc++fs or nothing
if(DEFINED CMAKE_CXX_STANDARD AND CMAKE_CXX_STANDARD LESS 17)
set(STD_FS_NO_LIB_NEEDED TRUE)
elseif(MSVC)
set(STD_FS_NO_LIB_NEEDED TRUE)
else()
file(
WRITE ${CMAKE_CURRENT_BINARY_DIR}/main.cpp
"#include <filesystem>\nint main(int argc, char ** argv) {\n std::filesystem::path p(argv[0]);\n return p.string().length();\n}"
)
try_compile(
STD_FS_NO_LIB_NEEDED ${CMAKE_CURRENT_BINARY_DIR}
SOURCES ${CMAKE_CURRENT_BINARY_DIR}/main.cpp
COMPILE_DEFINITIONS -std=c++17)
try_compile(
STD_FS_NEEDS_STDCXXFS ${CMAKE_CURRENT_BINARY_DIR}
SOURCES ${CMAKE_CURRENT_BINARY_DIR}/main.cpp
COMPILE_DEFINITIONS -std=c++17
LINK_LIBRARIES stdc++fs)
try_compile(
STD_FS_NEEDS_CXXFS ${CMAKE_CURRENT_BINARY_DIR}
SOURCES ${CMAKE_CURRENT_BINARY_DIR}/main.cpp
COMPILE_DEFINITIONS -std=c++17
LINK_LIBRARIES c++fs)
endif()
if(${STD_FS_NEEDS_STDCXXFS})
set(STD_FS_LIB stdc++fs PARENT_SCOPE)
elseif(${STD_FS_NEEDS_CXXFS})
set(STD_FS_LIB c++fs PARENT_SCOPE)
elseif(${STD_FS_NO_LIB_NEEDED})
set(STD_FS_LIB "" PARENT_SCOPE)
else()
message(WARNING "Unknown C++17 compiler - not passing -lstdc++fs")
set(STD_FS_LIB "")
endif()
endfunction()
if(${STD_FS_NEEDS_STDCXXFS})
set(STD_FS_LIB
stdc++fs
PARENT_SCOPE)
elseif(${STD_FS_NEEDS_CXXFS})
set(STD_FS_LIB
c++fs
PARENT_SCOPE)
elseif(${STD_FS_NO_LIB_NEEDED})
set(STD_FS_LIB
""
PARENT_SCOPE)
else()
message(WARNING "Unknown C++17 compiler - not passing -lstdc++fs")
set(STD_FS_LIB "")
endif()
endfunction()
cmake/package_config.cmake
+8 -12
View File
@@ -15,21 +15,17 @@ configure_package_config_file(
write_basic_package_version_file(
"${PROJECT_BINARY_DIR}/${PROJECT_NAME_LOWER}-config-version.cmake"
VERSION ${PROJECT_VERSION}
COMPATIBILITY SameMajorVersion
)
COMPATIBILITY SameMajorVersion)
install(FILES
"${PROJECT_BINARY_DIR}/${PROJECT_NAME_LOWER}-config.cmake"
"${PROJECT_BINARY_DIR}/${PROJECT_NAME_LOWER}-config-version.cmake"
install(
FILES "${PROJECT_BINARY_DIR}/${PROJECT_NAME_LOWER}-config.cmake"
"${PROJECT_BINARY_DIR}/${PROJECT_NAME_LOWER}-config-version.cmake"
COMPONENT devel
DESTINATION ${CMAKE_INSTALL_DIR}
)
DESTINATION ${CMAKE_INSTALL_DIR})
if (PROJECT_LIBRARIES OR PROJECT_STATIC_LIBRARIES)
if(PROJECT_LIBRARIES OR PROJECT_STATIC_LIBRARIES)
install(
EXPORT "${TARGETS_EXPORT_NAME}"
FILE ${PROJECT_NAME_LOWER}-targets.cmake
DESTINATION ${CMAKE_INSTALL_DIR}
)
endif ()
DESTINATION ${CMAKE_INSTALL_DIR})
endif()
conda-recipe/conda_build_config.yaml
+6 -2
View File
@@ -1,16 +1,20 @@
python:
- 3.11
- 3.12
- 3.13
- 3.14
c_compiler:
- gcc # [linux]
- clang # [osx]
c_stdlib:
- sysroot # [linux]
- macosx_deployment_target # [osx]
cxx_compiler:
- gxx # [linux]
- clangxx # [osx]
c_stdlib_version: # [linux]
c_stdlib_version:
- 2.17 # [linux]
- "11.0" # [osx]
conda-recipe/meta.yaml
+4 -2
View File
@@ -25,9 +25,9 @@ requirements:
host:
- python
- pip
- numpy=2.1
- numpy >=2.1
- scikit-build-core
- pybind11 >=2.13.0
- pybind11 >=3.0.3
- matplotlib # needed in host to solve the environment for run
run:
@@ -42,6 +42,7 @@ test:
- aare
requires:
- pytest
- pytest_check
- boost-histogram
source_files:
- python/tests
@@ -49,4 +50,5 @@ test:
- python -m pytest python/tests
about:
license: SPDX-License-Identifier MPL-2.0
summary: Data analysis library for hybrid pixel detectors from PSI
docs/CMakeLists.txt
+25 -30
View File
@@ -1,47 +1,42 @@
# SPDX-License-Identifier: MPL-2.0
find_package(Doxygen REQUIRED)
find_package(Sphinx REQUIRED)
#Doxygen
# Doxygen
set(DOXYGEN_IN ${CMAKE_CURRENT_SOURCE_DIR}/Doxyfile.in)
set(DOXYGEN_OUT ${CMAKE_CURRENT_BINARY_DIR}/Doxyfile)
configure_file(${DOXYGEN_IN} ${DOXYGEN_OUT} @ONLY)
#Sphinx
# Sphinx
set(SPHINX_SOURCE ${CMAKE_CURRENT_SOURCE_DIR}/src)
set(SPHINX_BUILD ${CMAKE_CURRENT_BINARY_DIR})
file(
GLOB_RECURSE SPHINX_SOURCE_FILES CONFIGURE_DEPENDS
RELATIVE "${CMAKE_CURRENT_SOURCE_DIR}/src"
"${CMAKE_CURRENT_SOURCE_DIR}/src/*.rst")
file(GLOB SPHINX_SOURCE_FILES CONFIGURE_DEPENDS "src/*.rst")
foreach(relpath IN LISTS SPHINX_SOURCE_FILES)
set(src "${CMAKE_CURRENT_SOURCE_DIR}/src/${relpath}")
set(dst "${SPHINX_BUILD}/src/${relpath}")
message(STATUS "Copying ${src} to ${dst}")
configure_file("${src}" "${dst}" COPYONLY)
endforeach()
configure_file("${CMAKE_CURRENT_SOURCE_DIR}/conf.py.in"
"${SPHINX_BUILD}/conf.py" @ONLY)
foreach(filename ${SPHINX_SOURCE_FILES})
get_filename_component(fname ${filename} NAME)
message(STATUS "Copying ${filename} to ${SPHINX_BUILD}/src/${fname}")
configure_file(${filename} "${SPHINX_BUILD}/src/${fname}")
endforeach(filename ${SPHINX_SOURCE_FILES})
file(COPY "${CMAKE_CURRENT_SOURCE_DIR}/figures" DESTINATION "${SPHINX_BUILD}")
configure_file(
"${CMAKE_CURRENT_SOURCE_DIR}/conf.py.in"
"${SPHINX_BUILD}/conf.py"
@ONLY
)
configure_file(
"${CMAKE_CURRENT_SOURCE_DIR}/static/extra.css"
"${SPHINX_BUILD}/static/css/extra.css"
@ONLY
)
configure_file("${CMAKE_CURRENT_SOURCE_DIR}/static/extra.css"
"${SPHINX_BUILD}/static/css/extra.css" @ONLY)
add_custom_target(
docs
COMMAND ${DOXYGEN_EXECUTABLE} ${DOXYGEN_OUT}
COMMAND ${SPHINX_EXECUTABLE} -a -b html
-Dbreathe_projects.aare=${CMAKE_CURRENT_BINARY_DIR}/xml
-c "${SPHINX_BUILD}"
${SPHINX_BUILD}/src
${SPHINX_BUILD}/html
COMMENT "Generating documentation with Sphinx"
)
docs
COMMAND ${DOXYGEN_EXECUTABLE} ${DOXYGEN_OUT}
COMMAND
${SPHINX_EXECUTABLE} -W -a -b html
-Dbreathe_projects.aare=${CMAKE_CURRENT_BINARY_DIR}/xml -c "${SPHINX_BUILD}"
${SPHINX_BUILD}/src ${SPHINX_BUILD}/html
COMMENT "Generating documentation with Sphinx")
docs/figures/Eta2x2.pdf
BIN
View File
Binary file not shown.
docs/figures/Eta2x2.png
BIN
View File
Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 9.3 KiB

docs/figures/Eta2x2Full.pdf
BIN
View File
Binary file not shown.
docs/figures/Eta2x2Full.png
BIN
View File
Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 14 KiB

docs/figures/Eta3x3.pdf
BIN
View File
Binary file not shown.
docs/figures/Eta3x3.png
BIN
View File
Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 20 KiB

docs/figures/Eta3x3Cross.pdf
BIN
View File
Binary file not shown.
docs/figures/Eta3x3Cross.png
BIN
View File
Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 9.5 KiB

docs/src/Cluster.rst
+15
View File
@@ -0,0 +1,15 @@
Cluster
========
.. doxygenstruct:: aare::Cluster
:members:
:undoc-members:
:private-members:
**Free Functions:**
.. doxygenfunction:: aare::reduce_to_3x3(const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>&)
.. doxygenfunction:: aare::reduce_to_2x2(const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>&)
docs/src/ClusterFile.rst
+34 -2
View File
@@ -1,8 +1,40 @@
ClusterFile
=============
Use ``frames()`` to iterate over complete frames, including empty frames and
frames whose clusters are all rejected by ROI or noise filtering. Each result
keeps its stored frame number; missing frame numbers are not synthesized.
.. code-block:: cpp
using ClusterType = aare::Cluster<int32_t, 3, 3>;
aare::ClusterFile<ClusterType> file("clusters.clust");
for (auto &frame : file.frames()) {
process_frame(frame.frame_number(), frame);
}
Use ``chunks()`` for the constructor's chunk size, or ``chunks(10000)`` to
request a different size for that traversal. The size must be positive and
counts selected clusters after filtering. Chunks may split or combine frames,
so their frame numbers are not reliable per-cluster metadata. Empty chunks are
not yielded; only the final chunk can contain fewer clusters than requested.
Both ranges apply the configured gain map and report incomplete files as
errors, just like the explicit read methods.
Ranges borrow the file and consume its current position without rewinding.
Constructing a range does not read; ``begin()`` reads the first result and
increment reads the next. Breaking a loop does not read ahead. Use one
traversal at a time and separate files for independent cursors. Switching to
frame reads after a chunk stops partway through a frame raises an error until
the remaining clusters in that frame have been read.
The file must outlive its ranges and iterators and must not be moved while
they are in use. Iterators are move only and support C++17 range-based loops,
not algorithms requiring copyable STL input iterators. References to a result
last until advancement or iterator destruction. Move a result out with
``auto retained = std::move(frame)`` to retain it. Frame iteration reuses the
current vector's storage when it has not been moved out.
.. doxygenclass:: aare::ClusterFile
:members:
:undoc-members:
:private-members:
docs/src/Installation.rst
+25 -4
View File
@@ -15,9 +15,16 @@ dependencies on older distributions.
.. note ::
aare is developing rapidly. Check for the latest release by
aare is developing rapidly. Check for available releases by
using: **conda search aare -c slsdetectorgroup**
To install the latest version of aare using conda, run:
.. code-block:: bash
conda install aare -c slsdetectorgroup
To install a specific version of aare, run:
.. code-block:: bash
@@ -28,7 +35,7 @@ dependencies on older distributions.
cmake build (development install)
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
If you are working on aare or want to test our a version that doesn't yet have
If you are working on aare or want to test out a version that doesn't yet have
a conda package. Build using cmake and then run from the build folder.
.. code-block:: bash
@@ -47,6 +54,20 @@ a conda package. Build using cmake and then run from the build folder.
# add the build folder to your PYTHONPATH and then you should be able to
# import aare in python
To build the python bindings and use aare in python compile with option ``-DAARE_PYTHON_BINDINGS=ON``.
.. code-block:: bash
#in build folder
cmake ../aare -DAARE_PYTHON_BINDINGS=ON
make -j4
Make sure to add the build folder to your PYTHONPATH and then you should be able to import aare in python.
.. code-block:: bash
export PYTHONPATH=/path/to/aare/build:$PYTHONPATH
cmake build + install and use in your C++ project
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -79,9 +100,9 @@ For detailed options see the CMakeLists.txt file in the root directory of the pr
Use system libraries instead of using FetchContent to pull in dependencies. Default option is off.
**AARE_PYTHON_BINDINGS "Build python bindings" ON**
**AARE_PYTHON_BINDINGS "Build python bindings" OFF**
Build the Python bindings. Default option is on.
Build the Python bindings. Default option is off.
.. warning ::
docs/src/Interpolation.rst
+160
View File
@@ -0,0 +1,160 @@
.. _Interpolation_C++API:
Interpolation
==============
The Interpolation class implements the :math:`\eta`-interpolation method.
This interpolation technique is based on charge sharing: for detected photon hits (e.g. clusters), it refines the estimated photon hit using information from neighboring pixels.
The method relies on the so-called :math:`\eta`-functions, which describe the relationship between the energy measured in the central cluster pixel (the initially estimated photon hit) and the energies measured in its neighboring pixels.
Depending on how much energy each neighboring pixel receives relative to the central pixel, the estimated photon hit is shifted toward that neighbor by a certain offset to the actual photon hit position in the pixel :math:`(x, y)`.
The mapping between the :math:`\eta` values and the corresponding spatial photon position :math:`(x,y)` can be viewed as an optimal transport problem.
One can readily compute the probability distribution :math:`P_{\eta}` of the :math:`\eta` values by forming a 2D histogram.
However, the probability distribution :math:`P_{x,y}` of the true photon positions is generally unknown unless the detector is illuminated uniformly (i.e. under flat-field conditions).
In a flat-field, the photon positions are uniformly distributed.
With this assumption, the problem reduces to determining a transport map :math:`T:(\eta_x,\eta_y) \rightarrow (x,y)`, that pushes forward the distribution of :math:`(\eta_x, \eta_y)` to the known uniform distribution of photon positions of a flatfield.
The map :math:`T` is given by:
.. math::
\begin{align*}
T_1: & F_{x}^{-1} F_{\eta_x|\eta_y} \\
T_2: & F_{y}^{-1} F_{\eta_y|\eta_x},
\end{align*}
where :math:`F_{\eta_x|\eta_y}` and :math:`F_{\eta_y|\eta_x}` are the conditional cumulative distribution functions e.g. :math:`F_{\eta_x|\eta_y}(\eta_x', \eta_y') = P_{\eta_x, \eta_y}(\eta_x \leq \eta_x' | \eta_y = \eta_y')`.
And :math:`F_{x}` and :math:`F_{y}` are the cumulative distribution functions of :math:`x` and :math:`y`. Note as :math:`x` and :math:`y` are uniformly distributed :math:`F_{x}` and :math:`F_{y}` are the identity functions. The map :math:`T` thus simplifies to
.. math::
\begin{align*}
T_1: & F_{\eta_x|\eta_y} \\
T_2: & F_{\eta_y|\eta_x}.
\end{align*}
Note that for the implementation :math:`P_{\eta}` is not only a distribution of :math:`\eta_x`, :math:`\eta_y` but also of the estimated photon energy :math:`e`.
The energy level correlates slightly with the z-depth. Higher z-depth leads to more charge sharing and a different :math:`\eta` distribution. Thus we create a mapping :math:`T` for each energy level.
:math:`\eta`-Functions:
---------------------------
.. doxygenstruct:: aare::Eta2
:members:
:undoc-members:
:private-members:
.. note::
The corner value ``c`` is only relevant when one uses ``calculate_eta_2`` or ``calculate_full_eta2``. Otherwise its default value is ``cTopLeft``.
Supported are the following :math:`\eta`-functions:
:math:`\eta`-Function on 2x2 Clusters:
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
.. image:: ../figures/Eta2x2.png
:target: ../figures/Eta2x2.png
:width: 650px
:align: center
:alt: Eta2x2
.. math::
\begin{equation*}
{\color{blue}{\eta_x}} = \frac{Q_{1,1}}{Q_{1,0} + Q_{1,1}} \quad \quad
{\color{green}{\eta_y}} = \frac{Q_{1,1}}{Q_{0,1} + Q_{1,1}}
\end{equation*}
The :math:`\eta` values can range between 0,1. Note they only range between 0,1 because the position of the center pixel (red) can change.
If the center pixel is in the bottom left pixel :math:`\eta_x` will be close to zero. If the center pixel is in the bottom right pixel :math:`\eta_y` will be close to 1.
One can apply this :math:`\eta` not only on 2x2 clusters but on clusters with any size. Then the 2x2 subcluster with maximum energy is chosen and the :math:`\eta` function applied on the subcluster.
.. doxygenfunction:: aare::calculate_eta2(const ClusterVector<ClusterType>&)
.. doxygenfunction:: aare::calculate_eta2(const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>&)
Full :math:`\eta`-Function on 2x2 Clusters:
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
.. image:: ../figures/Eta2x2Full.png
:target: ../figures/Eta2x2Full.png
:width: 650px
:align: center
:alt: Eta2x2 Full
.. math::
\begin{equation*}
{\color{blue}{\eta_x}} = \frac{Q_{0,1} + Q_{1,1}}{\sum_i^{1}\sum_j^{1}Q_{i,j}} \quad \quad
{\textcolor{green}{\eta_y}} = \frac{Q_{1,0} + Q_{1,1}}{\sum_i^{1}\sum_j^{1}Q_{i,j}}
\end{equation*}
The :math:`\eta` values can range between 0,1. Note they only range between 0,1 because the position of the center pixel (red) can change.
If the center pixel is in the bottom left pixel :math:`\eta_x` will be close to zero. If the center pixel is in the bottom right pixel :math:`\eta_y` will be close to 1.
.. doxygenfunction:: aare::calculate_full_eta2(const ClusterVector<ClusterType>&)
.. doxygenfunction:: aare::calculate_full_eta2(const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>&)
Full :math:`\eta`-Function on 3x3 Clusters:
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
.. image:: ../figures/Eta3x3.png
:target: ../figures/Eta3x3.png
:width: 650px
:align: center
:alt: Eta3x3
.. math::
\begin{equation*}
{\color{blue}{\eta_x}} = \frac{\sum_{i=0}^{2} Q_{i,2} - \sum_{i=0}^{2} Q_{i,0}}{\sum_{i=0}^{2}\sum_{j=0}^{2} Q_{i,j}} \quad \quad
{\color{green}{\eta_y}} = \frac{\sum_{j=0}^{2} Q_{2,j} - \sum_{j=0}^{2} Q_{0,j}}{\sum_{i=0}^{2}\sum_{j=0}^{2} Q_{i,j}}
\end{equation*}
The :math:`\eta` values can range between -0.5,0.5.
.. doxygenfunction:: aare::calculate_eta3(const ClusterVector<ClusterType>&)
.. doxygenfunction:: aare::calculate_eta3(const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>&)
Cross :math:`\eta`-Function on 3x3 Clusters:
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
.. image:: ../figures/Eta3x3Cross.png
:target: ../figures/Eta3x3Cross.png
:width: 650px
:align: center
:alt: Cross Eta3x3
.. math::
\begin{equation*}
{\color{blue}{\eta_x}} = \frac{Q_{1,2} - Q_{1,0}}{Q_{1,0} + Q_{1,1} + Q_{1,2}} \quad \quad
{\color{green}{\eta_y}} = \frac{Q_{0,2} - Q_{0,1}}{Q_{0,1} + Q_{1,1} + Q_{2,1}}
\end{equation*}
The :math:`\eta` values can range between -0.5,0.5.
.. doxygenfunction:: aare::calculate_cross_eta3(const ClusterVector<ClusterType>&)
.. doxygenfunction:: aare::calculate_cross_eta3(const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>&)
Interpolation class:
---------------------
.. warning::
The interpolation might lead to erroneous photon positions for clusters at the borders of a frame. Make sure to filter out such cases.
.. Warning::
Make sure to use the same :math:`\eta`-function during interpolation as given by the joint :math:`\eta`-distribution passed to the constructor.
.. Note::
Make sure to use reasonable energy bins, when constructing the joint distribution. If data is too sparse for a given energy the interpolation will lead to erroneous results.
.. doxygenclass:: aare::Interpolator
:members:
:undoc-members:
:private-members:
docs/src/JungfrauDataFile.rst
+1 -1
View File
@@ -12,7 +12,7 @@ The structure of the file is:
* ...
There is no metadata indicating number of frames or the size of the image, but this
will be infered by this reader.
will be inferred by this reader.
.. doxygenstruct:: aare::JungfrauDataHeader
:members:
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/RawFile.rst
+41 -1
View File
@@ -1,8 +1,48 @@
RawFile
===============
Reading requires frame padding to be enabled
or ``Frame Discard Policy`` to be ``discardpartial``.
The constructor checks the
parsed master metadata before opening data subfiles and throws
``std::runtime_error`` with the master path if this requirement is not met.
This applies to JSON and legacy ``.raw`` masters, including when opened through
``aare::File``. ``RawMasterFile`` can still be used to inspect their metadata.
``total_frames()`` is the minimum actual frame count across all selected
``RawSubFile`` objects in all ROIs. Each subfile count includes its complete
series of ``.raw`` files and is determined when the file is opened. Disabled
UDP ports are excluded. An empty subfile makes the count zero; opening a file
with no selected subfiles raises an error.
A warning is logged once on opening if the subfile counts differ or their
minimum differs from ``master().frames_in_file()``. It reports the subfile count
range, the master count, and the count used for reading. The master metadata
is preserved. Its expected frame count does not affect the count used for
reading.
Frame, batch, ROI, and frame-number reads use this minimum as their bound.
Frame-number synchronization may still fail within this bound if detector
frame numbers differ between subfiles.
Read errors include the attempted zero-based frame index and the failing file
path. Errors from a data subfile identify that subfile's path and preserve
the underlying failure reason without adding a master-file wrapper or C++
source location. Out-of-range subfile errors also report the available frame
count across that subfile's series. This applies to frame, batch, ROI, and
frame-number reads, including reads through ``aare::File``.
The frame index is the position in the file, not the frame number recorded by
the detector. Batch errors report the index of the frame that failed.
Synchronization errors identify the last ``.raw`` data file read for the module
that ran out of frames. Requests rejected as out of range before accessing a
subfile report the master path, since there is no single failing data file.
The top-level frame bounds error states the number of frames in the file and
that indices are zero-based.
.. doxygenclass:: aare::RawFile
:members:
:undoc-members:
:private-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/_install.rst
+4 -1
View File
@@ -6,7 +6,10 @@
cd build
#configure using cmake
cmake ../aare -DCMAKE_INSTALL_PREFIX=/where/to/put/aare
cmake ../aare -DCMAKE_INSTALL_PREFIX=/where/to/put/aare
#install python bindings
#cmake .. -DCMAKE_INSTALL_PREFIX=/where/to/put/aare -DAARE_INSTALL_PYTHONEXT=ON -DAARE_PYTHON_BINDINGS=ON
#build (replace 4 with the number of threads you want to use)
make -j4
docs/src/index.rst
+13 -12
View File
@@ -22,19 +22,17 @@ AARE
.. toctree::
:caption: Python API
:maxdepth: 1
pyFile
pycalibration
pyCtbRawFile
pyClusterFile
pyClusterVector
pyJungfrauDataFile
pyRawFile
pyRawMasterFile
pyVarClusterFinder
:maxdepth: 3
:hidden:
pycalibration
python/cluster/index
python/file/index
python/experimental/index
python/histogram/index
python/pedestal/index
pyFit
.. toctree::
@@ -46,11 +44,14 @@ AARE
NDView
Frame
File
MultiThreadedFileReader
Dtype
Cluster
ClusterFinder
ClusterFinderMT
ClusterFile
ClusterVector
Interpolation
JungfrauDataFile
Pedestal
RawFile
@@ -66,4 +67,4 @@ AARE
Philosophy
Workflow
Tests
Tests
docs/src/pyClusterFile.rst
-26
View File
@@ -1,26 +0,0 @@
ClusterFile
============
The :class:`ClusterFile` class is the main interface to read and write clusters in aare. Unfortunately the
old file format does not include metadata like the cluster size and the data type. This means that the
user has to know this information from other sources. Specifying the wrong cluster size or data type
will lead to garbage data being read.
.. py:currentmodule:: aare
.. autoclass:: ClusterFile
:members:
:undoc-members:
:inherited-members:
Below is the API of the ClusterFile_Cluster3x3i but all variants share the same API.
.. autoclass:: aare._aare.ClusterFile_Cluster3x3i
:special-members: __init__
:members:
:undoc-members:
:show-inheritance:
:inherited-members:
docs/src/pyClusterVector.rst
-49
View File
@@ -1,49 +0,0 @@
ClusterVector
================
The ClusterVector, holds clusters from the ClusterFinder. Since it is templated
in C++ we use a suffix indicating the type of cluster it holds. The suffix follows
the same pattern as for ClusterFile i.e. ``ClusterVector_Cluster3x3i``
for a vector holding 3x3 integer clusters.
At the moment the functionality from python is limited and it is not supported
to push_back clusters to the vector. The intended use case is to pass it to
C++ functions that support the ClusterVector or to view it as a numpy array.
**View ClusterVector as numpy array**
.. code:: python
from aare import ClusterFile
with ClusterFile("path/to/file") as f:
cluster_vector = f.read_frame()
# Create a copy of the cluster data in a numpy array
clusters = np.array(cluster_vector)
# Avoid copying the data by passing copy=False
clusters = np.array(cluster_vector, copy = False)
.. py:currentmodule:: aare
.. autoclass:: aare._aare.ClusterVector_Cluster3x3i
:special-members: __init__
:members:
:undoc-members:
:show-inheritance:
:inherited-members:
**Free Functions:**
.. autofunction:: reduce_to_3x3
:noindex:
Reduce a single Cluster to 3x3 by taking the 3x3 subcluster with highest photon energy.
.. autofunction:: reduce_to_2x2
:noindex:
Reduce a single Cluster to 2x2 by taking the 2x2 subcluster with highest photon energy.
docs/src/pyCtbRawFile.rst
-11
View File
@@ -1,11 +0,0 @@
CtbRawFile
============
.. py:currentmodule:: aare
.. autoclass:: CtbRawFile
:members:
:undoc-members:
:show-inheritance:
:inherited-members:
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/pyRawFile.rst
-10
View File
@@ -1,10 +0,0 @@
RawFile
===================
.. py:currentmodule:: aare
.. autoclass:: RawFile
:members:
:undoc-members:
:show-inheritance:
:inherited-members:
docs/src/python/cluster/index.rst
+11
View File
@@ -0,0 +1,11 @@
Cluster & Interpolation
==========================
.. toctree::
:caption: Cluster & Interpolation
:maxdepth: 1
pyCluster
pyClusterVector
pyInterpolation
pyVarClusterFinder
docs/src/python/cluster/pyCluster.rst
+23
View File
@@ -0,0 +1,23 @@
Cluster
========
.. py:currentmodule:: aare
.. autoclass:: Cluster
:members:
:undoc-members:
:inherited-members:
Below is the API of a cluster of size :math:`3\times 3` and type ``int`` but all variants share the same API.
.. autoclass:: aare._aare.Cluster3x3i
:special-members: __init__
:members:
:undoc-members:
:show-inheritance:
:inherited-members:
.. note::
More functions can be found in the :ref:`ClusterVector <py_clustervector>` documentation. Generally apply functions directly on the ``ClusterVector`` instead of looping over individual clusters.
docs/src/python/cluster/pyClusterVector.rst
+92
View File
@@ -0,0 +1,92 @@
.. _py_clustervector:
ClusterVector
================
A ClusterVector stores fixed-size clusters contiguously. Since it is templated
in C++, each bound class has a suffix indicating the cluster type. The suffix
follows the same pattern as ClusterFile; for example,
``ClusterVector_Cluster3x3i`` stores 3x3 clusters with 32-bit integer pixels.
The intended use case is to pass a ClusterVector to C++ functions that support
it or to view it as a NumPy array.
**View ClusterVector as a NumPy array**
.. code:: python
from aare import ClusterFile
with ClusterFile("path/to/file") as f:
cluster_vector = f.read_frame()
# Create a copy of the cluster data in a numpy array
clusters = np.array(cluster_vector)
# Avoid copying the data by passing copy=False
clusters = np.array(cluster_vector, copy=False)
.. warning::
A NumPy array created with ``copy=False`` is a view of the ClusterVector's
current storage. Do not call ``push_back`` or otherwise change the
ClusterVector while using the view. A ``push_back`` that reallocates the
backing buffer leaves existing NumPy views pointing to invalid memory, and
appending without reallocation does not update their shape. Use
``copy=True`` if the ClusterVector may change after creating the array.
.. py:currentmodule:: aare
.. autoclass:: ClusterVector
:members:
:undoc-members:
:inherited-members:
Below is the API of ``ClusterVector_Cluster3x3i``. All variants share the same
API.
.. autoclass:: aare._aare.ClusterVector_Cluster3x3i
:special-members: __init__, __call__
:members:
:undoc-members:
:show-inheritance:
:inherited-members:
**Free Functions:**
.. py:function:: hitmap(image_size, clusters)
:noindex:
Count cluster centers into an ``int32`` image. ``image_size`` is given as
``(rows, columns)``, and output element ``[y, x]`` contains the number of
photon hits at that coordinate. Out-of-bounds hits are ignored. All registered
ClusterVector variants are accepted.
:param tuple[int, int] image_size: Shape of the output image.
:param ClusterVector clusters: Clusters whose centers are counted.
:return: Hit counts with shape ``image_size``.
:rtype: numpy.ndarray
.. py:function:: reduce_to_3x3(clustervector)
:noindex:
Return a new vector containing the central 3x3 block of every input cluster.
Cluster order, coordinates, frame number, and pixel dtype are preserved.
:param ClusterVector clustervector: Input clusters with both dimensions at
least 3 and at least one dimension greater than 3.
:return: Reduced 3x3 clusters.
:rtype: ClusterVector
.. py:function:: reduce_to_2x2(clustervector)
:noindex:
Return a new vector containing the highest-sum center-adjacent 2x2 block of
every input cluster. Cluster order, coordinates, frame number, and pixel
dtype are preserved.
:param ClusterVector clustervector: Input clusters of size 2x2 or larger.
:return: Reduced 2x2 clusters.
:rtype: ClusterVector
docs/src/python/cluster/pyInterpolation.rst
+124
View File
@@ -0,0 +1,124 @@
Interpolation
==============
The Interpolation class implements the :math:`\eta`-interpolation method.
This interpolation technique is based on charge sharing: for detected photon hits (e.g. clusters), it refines the estimated photon hit using information from neighboring pixels.
See :ref:`Interpolation_C++API` for a more elaborate documentation and explanation of the method.
:math:`\eta`-Functions:
--------------------------
Below is an example of the Eta class of type ``double``. Supported are ``Etaf`` of type ``float`` and ``Etai`` of type ``int``.
.. autoclass:: aare._aare.Etad
:members:
:private-members:
.. note::
The corner value ``c`` is only relevant when one uses ``calculate_eta_2`` or ``calculate_full_eta2``. Otherwise its default value is ``cTopLeft``.
Supported are the following :math:`\eta`-functions:
:math:`\eta`-Function on 2x2 Clusters:
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
.. py:currentmodule:: aare
.. image:: ../../../figures/Eta2x2.png
:target: ../../../figures/Eta2x2.png
:width: 650px
:align: center
:alt: Eta2x2
.. math::
\begin{equation*}
{\color{blue}{\eta_x}} = \frac{Q_{1,1}}{Q_{1,0} + Q_{1,1}} \quad \quad
{\color{green}{\eta_y}} = \frac{Q_{1,1}}{Q_{0,1} + Q_{1,1}}
\end{equation*}
The :math:`\eta` values can range between 0,1. Note they only range between 0,1 because the position of the center pixel (red) can change.
If the center pixel is in the bottom left pixel :math:`\eta_x` will be close to zero. If the center pixel is in the bottom right pixel :math:`\eta_y` will be close to 1.
.. autofunction:: calculate_eta2
Full :math:`\eta`-Function on 2x2 Clusters:
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
.. image:: ../../../figures/Eta2x2Full.png
:target: ../../../figures/Eta2x2Full.png
:width: 650px
:align: center
:alt: Eta2x2 Full
.. math::
\begin{equation*}
{\color{blue}{\eta_x}} = \frac{Q_{0,1} + Q_{1,1}}{\sum_{i=0}^{1}\sum_{j=0}^{1}Q_{i,j}} \quad \quad
{\textcolor{green}{\eta_y}} = \frac{Q_{1,0} + Q_{1,1}}{\sum_{i=0}^{1}\sum_{j=0}^{1}Q_{i,j}}
\end{equation*}
The :math:`\eta` values can range between 0,1. Note they only range between 0,1 because the position of the center pixel (red) can change.
If the center pixel is in the bottom left pixel :math:`\eta_x` will be close to zero. If the center pixel is in the bottom right pixel :math:`\eta_y` will be close to 1.
.. autofunction:: calculate_full_eta2
Full :math:`\eta`-Function on 3x3 Clusters:
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
.. image:: ../../../figures/Eta3x3.png
:target: ../../../figures/Eta3x3.png
:width: 650px
:align: center
:alt: Eta3x3
.. math::
\begin{equation*}
{\color{blue}{\eta_x}} = \frac{\sum_{i=0}^{2} Q_{i,2} - \sum_{i=0}^{2} Q_{i,0}}{\sum_{i=0}^{2}\sum_{j}^{3} Q_{i,j}} \quad \quad
{\color{green}{\eta_y}} = \frac{\sum_{j=0}^{2} Q_{2,j} - \sum_{j=0}^{2} Q_{0,j}}{\sum_{i=0}^{2}\sum_{j}^{3} Q_{i,j}}
\end{equation*}
The :math:`\eta` values can range between -0.5,0.5.
.. autofunction:: calculate_eta3
Cross :math:`\eta`-Function on 3x3 Clusters:
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
.. image:: ../../../figures/Eta3x3Cross.png
:target: ../../../figures/Eta3x3Cross.png
:width: 650px
:align: center
:alt: Cross Eta3x3
.. math::
\begin{equation*}
{\color{blue}{\eta_x}} = \frac{Q_{1,2} - Q_{1,0}}{Q_{1,0} + Q_{1,1} + Q_{1,2}} \quad \quad
{\color{green}{\eta_y}} = \frac{Q_{0,2} - Q_{0,1}}{Q_{0,1} + Q_{1,1} + Q_{2,1}}
\end{equation*}
The :math:`\eta` values can range between -0.5,0.5.
.. autofunction:: calculate_cross_eta3
Interpolation class for :math:`\eta`-Interpolation
----------------------------------------------------
.. Warning::
Make sure to use the same :math:`\eta`-function during interpolation as given by the joint :math:`\eta`-distribution passed to the constructor.
.. Warning::
The interpolation might lead to erroneous photon positions for clusters at the boarders of a frame. Make sure to filter out such cases.
.. Note::
Make sure to use reasonable energy bins, when constructing the joint distribution. If data is too sparse for a given energy the interpolation will lead to erreneous results.
.. py:currentmodule:: aare
.. autoclass:: Interpolator
:special-members: __init__
:members:
:undoc-members:
:inherited-members:
docs/src/pyVarClusterFinder.rst → docs/src/python/cluster/pyVarClusterFinder.rst
View File
File renamed without changes.
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/file/index.rst
+14
View File
@@ -0,0 +1,14 @@
File I/O
===================
.. toctree::
:caption: File I/O
:maxdepth: 1
pyClusterFile
pyCtbRawFile
pyFile
pyJungfrauDataFile
pyRawFile
pyRawMasterFile
pyTransform
docs/src/python/file/pyClusterFile.rst
+78
View File
@@ -0,0 +1,78 @@
ClusterFile
===========
The :func:`ClusterFile` factory is the main interface for reading and writing
cluster files. Use ``mode="r"`` to read, ``mode="w"`` to truncate and
write, or ``mode="a"`` to append.
The format does not store the cluster dimensions, value type, coordinate type,
padding, or byte order. The ``cluster_size`` and ``dtype`` arguments must match
the writer; otherwise the bytes are interpreted incorrectly.
Use ``frames()`` when frame boundaries and the frame number matter:
.. code-block:: python
with ClusterFile("clusters.clust") as file:
for frame in file.frames():
print(frame.frame_number, frame.size)
Each stored frame is yielded, including empty frames and frames whose clusters
are all rejected by filtering. Missing frame numbers are not synthesized.
The explicit ``read_frame()`` method returns one frame or ``None`` at EOF.
Use ``chunks()`` to iterate over chunks of ``chunk_size`` clusters. Per default
the chunk size defined in the constructor is used. Alternatively one can pass
an explicit chunk size to the function for that specific traversal:
.. code-block:: python
with ClusterFile("clusters.clust", chunk_size=1000) as file:
for clusters in file.chunks(10000):
process(clusters)
Chunk sizes must be positive and count selected clusters after filtering.
``file.chunks(10000)`` does not overwrite the constructor argument ``chunk_size``
. Iterating directly
over ``file`` or calling ``next(file)`` still reads chunks using that default.
Chunks and ``read_clusters(n_clusters)`` may split or combine frames, so their
frame numbers are not reliable per-cluster metadata. Empty chunks are not
yielded.
Chunk reads and iteration raise an error if they encounter an incomplete frame
header or cluster record, including when ROI or noise filtering is enabled.
They return fewer clusters than requested only at a clean end of file.
Both iterators consume the current file position without rewinding. Creating
an iterator does not read; each ``next()`` reads one result without prefetching
the following result. Use one traversal at a time and separate readers for
independent cursors. Switching to frames after a chunk stops partway through a
frame raises an error until the remaining clusters in that frame have been
read. Exhausted iterators continue to raise ``StopIteration``.
Iterators keep their file object alive, but explicitly closing the file or
leaving its ``with`` block prevents further reads. Every yielded
``ClusterVector`` owns its storage and can be retained after iteration advances
or the file closes. NumPy views also remain valid across those operations;
as usual, resizing or modifying the allocation of their backing vector can
invalidate them.
When a gain map is configured, it is applied to every cluster whose complete
footprint lies inside the map. Clusters whose footprint extends beyond the
gain-map boundaries remain in the returned cluster vector, but all their data
values are set to zero.
.. py:currentmodule:: aare
.. autofunction:: ClusterFile
Below is the API of ``ClusterFile_Cluster3x3i``; all compiled variants share the
same API.
.. autoclass:: aare._aare.ClusterFile_Cluster3x3i
:special-members: __init__
:members:
:undoc-members:
:show-inheritance:
:inherited-members:
docs/src/python/file/pyCtbRawFile.rst
+29
View File
@@ -0,0 +1,29 @@
CtbRawFile
============
Read analog, digital and transceiver samples from a raw file containing
data from the Chip Test Board. Uses :mod:`aare.transform` to decode the
data into a format that the user can work with.
``read_frame(index)`` reads the frame at the zero-based index across the raw
subfiles and advances the position to ``index + 1``. A subsequent
``read_frame()`` reads the following frame, including across subfile boundaries.
.. code:: python
import aare
from aare.transform import Mythen302Transform
my302 = Mythen302Transform(offset = 4)
with aare.CtbRawFile(fname, transform = my302) as f:
for header, data in f:
#do something with the data
.. py:currentmodule:: aare
.. autoclass:: CtbRawFile
:members:
:undoc-members:
:show-inheritance:
:inherited-members:
docs/src/pyFile.rst → docs/src/python/file/pyFile.rst
View File
File renamed without changes.
docs/src/pyJungfrauDataFile.rst → docs/src/python/file/pyJungfrauDataFile.rst
View File
File renamed without changes.
docs/src/python/file/pyRawFile.rst
+47
View File
@@ -0,0 +1,47 @@
RawFile
===================
.. py:currentmodule:: aare
Reading requires frame padding to be enabled
or ``Frame Discard Policy`` to be ``discardpartial``. Construction raises
``RuntimeError`` with the master path before opening data subfiles if this
requirement is not met. This applies to JSON and legacy ``.raw`` masters,
including when opened through ``File``. ``RawMasterFile`` can still be used
to inspect their metadata.
``total_frames`` and ``len(reader)`` report the minimum actual frame count
across all selected subfiles and ROIs, including each subfile's complete
``.raw`` series. ``read()`` uses this count, and ``read_n()`` and
``read_n_with_roi()`` limit batches to the remaining frames. Counts are
determined when the reader is constructed.
An empty subfile makes the count zero, and construction fails if no subfiles
are selected.
The C++ logger prints one warning to stderr on opening if the subfile counts
differ or their minimum differs from ``reader.master.frames_in_file``. The
warning includes the count range, the recorded master count, and the count
used for reading. Master metadata is preserved; its expected frame count
does not affect reading bounds. Detector frame-number mismatches can still
cause synchronization errors within these bounds.
Reading errors raise ``RuntimeError`` with the attempted zero-based frame index
and the failing file path. If a data subfile fails, its path and the underlying
failure reason are reported. Out-of-range subfile errors include the available frame count across
that subfile's series. This applies to ``read_frame()``, ``read_n()``,
and the ROI reading methods. The frame index is the position in the file, not
the frame number recorded by the detector.
Synchronization errors also identify the last ``.raw`` data file read for the
module that ran out of frames. An out-of-range request rejected before accessing
a subfile reports the master path.
For a file containing five frames, the valid indices are 0 through 4. Attempting
to read index 5 raises an error stating: ``Frame index 5 is out of range: file
contains 5 frames (indices are zero-based)``.
.. autoclass:: RawFile
:members:
:undoc-members:
:show-inheritance:
:inherited-members:
docs/src/pyRawMasterFile.rst → docs/src/python/file/pyRawMasterFile.rst
View File
File renamed without changes.
docs/src/python/file/pyTransform.rst
+27
View File
@@ -0,0 +1,27 @@
Transform
===================
The transform module takes data read by :class:`aare.CtbRawFile` and decodes it
to a useful image format. Depending on detector it supports both analog
and digital samples.
For convenience the following transform objects are defined with a short name
.. code:: python
moench05 = Moench05Transform()
moench05_1g = Moench05Transform1g()
moench05_old = Moench05TransformOld()
matterhorn02 = Matterhorn02Transform()
adc_sar_04_64to16 = AdcSar04Transform64to16()
adc_sar_05_64to16 = AdcSar05Transform64to16()
.. py:currentmodule:: aare
.. automodule:: aare.transform
:members:
:undoc-members:
:private-members:
:special-members: __call__
:show-inheritance:
:inherited-members:
docs/src/python/histogram/index.rst
+9
View File
@@ -0,0 +1,9 @@
Histogram
=========
.. toctree::
:caption: Histogram
:maxdepth: 1
pyPixelHistogram
pyPedestalTrackingPixelHistogram
docs/src/python/histogram/pyPedestalTrackingPixelHistogram.rst
+32
View File
@@ -0,0 +1,32 @@
PedestalTrackingPixelHistogram
==============================
.. warning::
``PedestalTrackingPixelHistogram`` is specifically designed for use in the Jungfrau calibration
pipeline. Make sure you understand the behaviour before using it in other contexts.
``PedestalTrackingPixelHistogram`` accumulates a pixel-wise histogram of
``frame - pedestal`` residuals while maintaining a running per-pixel pedestal
estimate.
Use ``push_pedestal_no_update()`` to seed the pedestal estimate, then
``update_mean()`` before submitting frames with ``fill_async()``. Pending
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
:members:
:undoc-members:
:show-inheritance:
:inherited-members:
docs/src/python/histogram/pyPixelHistogram.rst
+41
View File
@@ -0,0 +1,41 @@
PixelHistogram
==============
``PixelHistogram`` accumulates one histogram per detector pixel from 2D
``float64`` images. The public ``PixelHistogram`` name is an alias for
``PixelHistogram_d``, which stores bin counts as ``float64``.
.. note::
``PixelHistogram`` is designed for fast filling from 2D images, utilizing multiple threads,
and contiguous storage of a specific type. No over/underflow bins are provided. For generic use cases consider
boost-histogram. (https://boost-histogram.readthedocs.io/en/latest/)
Storage-specific variants are also available when a smaller or integer count
type is preferred:
* ``PixelHistogram_d``: ``float64`` storage
* ``PixelHistogram_f``: ``float32`` storage
* ``PixelHistogram_u64``: ``uint64`` storage
* ``PixelHistogram_u32``: ``uint32`` storage
* ``PixelHistogram_u16``: ``uint16`` storage
* ``PixelHistogram_u8``: ``uint8`` storage
.. py:currentmodule:: aare
.. autoclass:: PixelHistogram
:members:
:undoc-members:
:show-inheritance:
:inherited-members:
Showing API for PixelHistogram_d, all variant share the same API.
.. autoclass:: aare._aare.PixelHistogram_d
:members:
:undoc-members:
:show-inheritance:
:inherited-members:
:noindex:
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
+86
View File
@@ -0,0 +1,86 @@
FastPedestal
============
``FastPedestal`` calculates a running mean 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``
Internal moments and variance are calculated in double precision. Variance is
private and stays in double precision through the square root; the cached mean
and on-demand standard deviation are returned in the specified type. Negative
variance caused by floating-point roundoff is clamped to zero.
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
+69
View File
@@ -0,0 +1,69 @@
Pedestal
========
``Pedestal`` calculates a running mean and population standard deviation for
each pixel in a series of ``uint16`` frames. ``push()`` updates the cached mean
immediately; ``std()`` calculates the noise from the current statistics.
``push()`` and ``push_with_threshold()`` require C-contiguous, two-dimensional
NumPy frames with dtype ``uint16`` and shape matching the pedestal.
``push_with_threshold()`` also requires a C-contiguous, two-dimensional
threshold array with the pedestal's output dtype and shape. Noncontiguous
inputs raise ``TypeError``; inputs with the wrong number of dimensions raise
``ValueError``. Mismatched frame or threshold shapes raise ``RuntimeError``.
Use ``numpy.ascontiguousarray()`` to copy a sliced or transposed array into the
required layout when needed.
Internal sums and sums of squares always use ``float64``. Three
specializations are available from :mod:`aare` for the mean and standard
deviation output types:
* ``Pedestal_d`` returns ``float64``
* ``Pedestal_f`` returns ``float32``
* ``Pedestal_i16`` returns ``int16``
The public ``Pedestal`` factory selects the specialization from ``dtype``,
defaulting to ``numpy.float64``.
Constructor dimensions and ``n_samples`` must be positive integers. Negative
values raise ``TypeError`` and zero values raise ``RuntimeError``.
Internally, negative variance caused by floating-point roundoff is clamped to
zero before taking its square root, keeping the standard deviation finite for
nearly constant inputs. Only the final standard deviation is converted to the
output dtype.
Factory
-------
.. py:currentmodule:: aare
.. autofunction:: Pedestal
Example
-------
.. code-block:: python
import numpy as np
from aare import Pedestal
pedestal = Pedestal(512, 1024, n_samples=100, dtype=np.float32)
for frame in initialization_frames:
pedestal.push(frame)
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:
etc/add_license.py
+103
View File
@@ -0,0 +1,103 @@
#!/usr/bin/env python3
import argparse
import fnmatch
import os
from pathlib import Path
CPP_PATTERNS = ["*.h", "*.hpp", "*.cpp"]
PY_PATTERNS = ["*.py"]
CMAKE_PATTERNS = ["CMakeLists.txt"]
FILE_PATTERNS = CPP_PATTERNS + PY_PATTERNS + CMAKE_PATTERNS
LICENSE_TEXT = "SPDX-License-Identifier: MPL-2.0"
def get_comment_prefix(filename: str) -> str | None:
if any(fnmatch.fnmatch(filename, p) for p in CPP_PATTERNS):
return "// "
if any(fnmatch.fnmatch(filename, p) for p in (PY_PATTERNS + CMAKE_PATTERNS)):
return "# "
return None
def matches_pattern(filename: str) -> bool:
return any(fnmatch.fnmatch(filename, p) for p in FILE_PATTERNS)
def process_file(filepath: Path) -> bool:
filename = filepath.name
prefix = get_comment_prefix(filename)
if not prefix:
return False
license_line = f"{prefix}{LICENSE_TEXT}\n"
try:
with filepath.open("r", encoding="utf-8") as f:
lines = f.readlines()
except Exception as e:
print(f"⚠️ Error reading {filepath}: {e}")
return False
# Skip if SPDX already present anywhere in the file
if any("SPDX-License-Identifier" in line for line in lines):
return False
insert_index = 0
# For Python, keep shebang on the very first line
if filename.endswith(".py") and lines:
if lines[0].startswith("#!"):
insert_index = 1
lines.insert(insert_index, license_line)
try:
with filepath.open("w", encoding="utf-8") as f:
f.writelines(lines)
except Exception as e:
print(f"⚠️ Error writing {filepath}: {e}")
return False
return True
def main() -> None:
parser = argparse.ArgumentParser(
description="Add SPDX-License-Identifier: MPL-2.0 to source files."
)
parser.add_argument(
"path",
help="Root directory to recursively process "
"(*.h, *.cpp, *.py, and CMakeLists.txt).",
)
args = parser.parse_args()
root_path = Path(args.path).expanduser().resolve()
if not root_path.exists():
print(f"Error: Path does not exist: {root_path}")
raise SystemExit(1)
if not root_path.is_dir():
print(f"Error: Path is not a directory: {root_path}")
raise SystemExit(1)
print(f"Processing directory: {root_path}")
modified = 0
for dirpath, _, files in os.walk(root_path):
dirpath = Path(dirpath)
for name in files:
if matches_pattern(name):
fullpath = dirpath / name
if process_file(fullpath):
print(f"✔ Added SPDX: {fullpath}")
modified += 1
print(f"\nDone. Updated {modified} file(s).")
if __name__ == "__main__":
main()
etc/dev-env.yml
+10 -5
View File
@@ -3,14 +3,19 @@ channels:
- conda-forge
dependencies:
- anaconda-client
- catch2
- conda-build
- doxygen
- doxygen
- sphinx
- breathe
- sphinx_rtd_theme
- furo
- zeromq
- breathe
- sphinx_rtd_theme
- furo
- zeromq
- pybind11
- numpy
- matplotlib
- nlohmann_json
- pytest
- pytest-check
- boost-histogram
update_version.py → etc/update_version.py
+5 -3
View File
@@ -1,4 +1,4 @@
# SPDX-License-Identifier: LGPL-3.0-or-other
# SPDX-License-Identifier: MPL-2.0
# Copyright (C) 2021 Contributors to the Aare Package
"""
Script to update VERSION file with semantic versioning if provided as an argument, or with 0.0.0 if no argument is provided.
@@ -8,11 +8,12 @@ import sys
import os
import re
from datetime import datetime
from pathlib import Path
from packaging.version import Version, InvalidVersion
SCRIPT_DIR = os.path.dirname(os.path.abspath(__file__))
SCRIPT_DIR = Path(__file__).absolute().parent.parent
def is_integer(value):
try:
@@ -46,7 +47,8 @@ def get_version():
def write_version_to_file(version):
version_file_path = os.path.join(SCRIPT_DIR, "VERSION")
version_file_path = SCRIPT_DIR/"VERSION"
print(version_file_path)
with open(version_file_path, "w") as version_file:
version_file.write(version)
print(f"Version {version} written to VERSION file.")
include/aare/ArrayExpr.hpp
+5 -4
View File
@@ -1,3 +1,4 @@
// SPDX-License-Identifier: MPL-2.0
#pragma once
#include "aare/defs.hpp"
#include <array>
@@ -28,7 +29,7 @@ class ArrayAdd : public ArrayExpr<ArrayAdd<A, B, Ndim>, Ndim> {
ArrayAdd(const A &arr1, const B &arr2) : arr1_(arr1), arr2_(arr2) {
assert(arr1.size() == arr2.size());
}
auto operator[](int i) const { return arr1_[i] + arr2_[i]; }
auto operator[](size_t i) const { return arr1_[i] + arr2_[i]; }
size_t size() const { return arr1_.size(); }
std::array<ssize_t, Ndim> shape() const { return arr1_.shape(); }
};
@@ -42,7 +43,7 @@ class ArraySub : public ArrayExpr<ArraySub<A, B, Ndim>, Ndim> {
ArraySub(const A &arr1, const B &arr2) : arr1_(arr1), arr2_(arr2) {
assert(arr1.size() == arr2.size());
}
auto operator[](int i) const { return arr1_[i] - arr2_[i]; }
auto operator[](size_t i) const { return arr1_[i] - arr2_[i]; }
size_t size() const { return arr1_.size(); }
std::array<ssize_t, Ndim> shape() const { return arr1_.shape(); }
};
@@ -56,7 +57,7 @@ class ArrayMul : public ArrayExpr<ArrayMul<A, B, Ndim>, Ndim> {
ArrayMul(const A &arr1, const B &arr2) : arr1_(arr1), arr2_(arr2) {
assert(arr1.size() == arr2.size());
}
auto operator[](int i) const { return arr1_[i] * arr2_[i]; }
auto operator[](size_t i) const { return arr1_[i] * arr2_[i]; }
size_t size() const { return arr1_.size(); }
std::array<ssize_t, Ndim> shape() const { return arr1_.shape(); }
};
@@ -70,7 +71,7 @@ class ArrayDiv : public ArrayExpr<ArrayDiv<A, B, Ndim>, Ndim> {
ArrayDiv(const A &arr1, const B &arr2) : arr1_(arr1), arr2_(arr2) {
assert(arr1.size() == arr2.size());
}
auto operator[](int i) const { return arr1_[i] / arr2_[i]; }
auto operator[](size_t i) const { return arr1_[i] / arr2_[i]; }
size_t size() const { return arr1_.size(); }
std::array<ssize_t, Ndim> shape() const { return arr1_.shape(); }
};
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/CalculateEta.hpp
+369 -95
View File
@@ -1,18 +1,13 @@
// SPDX-License-Identifier: MPL-2.0
#pragma once
#include "aare/Cluster.hpp"
#include "aare/ClusterVector.hpp"
#include "aare/NDArray.hpp"
#include "aare/defs.hpp"
namespace aare {
enum class corner : int {
cBottomLeft = 0,
cBottomRight = 1,
cTopLeft = 2,
cTopRight = 3
};
enum class pixel : int {
pBottomLeft = 0,
pBottom = 1,
@@ -25,149 +20,428 @@ enum class pixel : int {
pTopRight = 8
};
// TODO: better to have sum after x,y
/**
* eta struct
*/
template <typename T> struct Eta2 {
double x;
double y;
int c{0};
T sum;
/// @brief eta in x direction
double x{};
/// @brief eta in y direction
double y{};
/// @brief index of subcluster with highest energy value (given as corner
/// relative to cluster center)
corner c{0};
/// @brief photon energy (cluster sum)
T sum{};
};
/**
* @brief Calculate the eta2 values for all clusters in a Clustervector
* @brief Calculate the eta2 values for all clusters in a ClusterVector
*/
template <typename ClusterType,
typename = std::enable_if_t<is_cluster_v<ClusterType>>>
NDArray<double, 2> calculate_eta2(const ClusterVector<ClusterType> &clusters) {
NDArray<double, 2> eta2({static_cast<int64_t>(clusters.size()), 2});
std::vector<Eta2<typename ClusterType::value_type>>
calculate_eta2(const ClusterVector<ClusterType> &clusters) {
std::vector<Eta2<typename ClusterType::value_type>> eta2{};
eta2.reserve(clusters.size());
for (size_t i = 0; i < clusters.size(); i++) {
auto e = calculate_eta2(clusters[i]);
eta2(i, 0) = e.x;
eta2(i, 1) = e.y;
eta2.push_back(e);
}
return eta2;
}
/**
* @brief Calculate the full eta2 values for all clusters in a ClusterVector
*/
template <typename ClusterType,
typename = std::enable_if_t<is_cluster_v<ClusterType>>>
std::vector<Eta2<typename ClusterType::value_type>>
calculate_full_eta2(const ClusterVector<ClusterType> &clusters) {
std::vector<Eta2<typename ClusterType::value_type>> eta2{};
eta2.reserve(clusters.size());
for (size_t i = 0; i < clusters.size(); i++) {
auto e = calculate_full_eta2(clusters[i]);
eta2.push_back(e);
}
return eta2;
}
/**
* @brief Calculate eta3 for all 3x3 clusters in a ClusterVector
*/
template <typename ClusterType,
typename = std::enable_if_t<is_cluster_v<ClusterType>>>
std::vector<Eta2<typename ClusterType::value_type>>
calculate_eta3(const ClusterVector<ClusterType> &clusters) {
std::vector<Eta2<typename ClusterType::value_type>> eta2{};
eta2.reserve(clusters.size());
for (size_t i = 0; i < clusters.size(); i++) {
auto e = calculate_eta3(clusters[i]);
eta2.push_back(e);
}
return eta2;
}
/**
* @brief Calculate cross eta3 for all 3x3 clusters in a ClusterVector
*/
template <typename ClusterType,
typename = std::enable_if_t<is_cluster_v<ClusterType>>>
std::vector<Eta2<typename ClusterType::value_type>>
calculate_cross_eta3(const ClusterVector<ClusterType> &clusters) {
std::vector<Eta2<typename ClusterType::value_type>> eta2{};
eta2.reserve(clusters.size());
for (size_t i = 0; i < clusters.size(); i++) {
auto e = calculate_cross_eta3(clusters[i]);
eta2.push_back(e);
}
return eta2;
}
/**
* @brief helper function to calculate eta2 x and y values
* @param eta reference to the Eta2 object to update
* @param left_x value of the left pixel
* @param right_x value of the right pixel
* @param bottom_y value of the bottom pixel
* @param top_y value of the top pixel
*/
template <typename T>
inline void calculate_eta2(Eta2<T> &eta, const T left_x, const T right_x,
const T bottom_y, const T top_y) {
if ((right_x + left_x) != 0)
eta.x = static_cast<double>(right_x) /
static_cast<double>(right_x + left_x); // between (0,1) the
// closer to zero left
// value probably larger
if ((top_y + bottom_y) != 0)
eta.y = static_cast<double>(top_y) /
static_cast<double>(top_y + bottom_y); // between (0,1) the
// closer to zero bottom
// value probably larger
}
/**
* @brief Calculate the eta2 values for a generic sized cluster and return them
* in a Eta2 struct containing etay, etax and the index of the respective 2x2
* subcluster.
* in a Eta2 struct containing etay, etax and the index (as corner) of the
* respective 2x2 subcluster relative to the cluster center.
*/
template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
typename CoordType>
typename CoordType = uint16_t>
Eta2<T>
calculate_eta2(const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType> &cl) {
Eta2<T> eta{};
auto max_sum = cl.max_sum_2x2();
eta.sum = max_sum.first;
auto c = max_sum.second;
static_assert(ClusterSizeX > 1 && ClusterSizeY > 1);
Eta2<T> eta{};
size_t cluster_center_index =
(ClusterSizeX / 2) + (ClusterSizeY / 2) * ClusterSizeX;
size_t index_bottom_left_max_2x2_subcluster =
(int(c / (ClusterSizeX - 1))) * ClusterSizeX + c % (ClusterSizeX - 1);
auto max_sum = cl.max_sum_2x2();
eta.sum = max_sum.sum;
corner c = max_sum.index;
// calculate direction of gradient
// check that cluster center is in max subcluster
if (cluster_center_index != index_bottom_left_max_2x2_subcluster &&
cluster_center_index != index_bottom_left_max_2x2_subcluster + 1 &&
cluster_center_index !=
index_bottom_left_max_2x2_subcluster + ClusterSizeX &&
cluster_center_index !=
index_bottom_left_max_2x2_subcluster + ClusterSizeX + 1)
throw std::runtime_error("Photon center is not in max 2x2_subcluster");
if ((cluster_center_index - index_bottom_left_max_2x2_subcluster) %
ClusterSizeX ==
0) {
if ((cl.data[cluster_center_index + 1] +
cl.data[cluster_center_index]) != 0)
eta.x = static_cast<double>(cl.data[cluster_center_index + 1]) /
static_cast<double>((cl.data[cluster_center_index + 1] +
cl.data[cluster_center_index]));
} else {
if ((cl.data[cluster_center_index] +
cl.data[cluster_center_index - 1]) != 0)
eta.x = static_cast<double>(cl.data[cluster_center_index]) /
static_cast<double>((cl.data[cluster_center_index - 1] +
cl.data[cluster_center_index]));
}
if ((cluster_center_index - index_bottom_left_max_2x2_subcluster) /
ClusterSizeX <
1) {
assert(cluster_center_index + ClusterSizeX <
ClusterSizeX * ClusterSizeY); // suppress warning
if ((cl.data[cluster_center_index] +
cl.data[cluster_center_index + ClusterSizeX]) != 0)
eta.y = static_cast<double>(
cl.data[cluster_center_index + ClusterSizeX]) /
static_cast<double>(
(cl.data[cluster_center_index] +
cl.data[cluster_center_index + ClusterSizeX]));
} else {
if ((cl.data[cluster_center_index] +
cl.data[cluster_center_index - ClusterSizeX]) != 0)
eta.y = static_cast<double>(cl.data[cluster_center_index]) /
static_cast<double>(
(cl.data[cluster_center_index] +
cl.data[cluster_center_index - ClusterSizeX]));
// subcluster top right from center
switch (c) {
case (corner::cTopLeft):
calculate_eta2(eta, cl.data[cluster_center_index - 1],
cl.data[cluster_center_index],
cl.data[cluster_center_index - ClusterSizeX],
cl.data[cluster_center_index]);
// dx = -1
// dy = -1
break;
case (corner::cTopRight):
calculate_eta2(eta, cl.data[cluster_center_index],
cl.data[cluster_center_index + 1],
cl.data[cluster_center_index - ClusterSizeX],
cl.data[cluster_center_index]);
// dx = 0
// dy = -1
break;
case (corner::cBottomLeft):
calculate_eta2(eta, cl.data[cluster_center_index - 1],
cl.data[cluster_center_index],
cl.data[cluster_center_index],
cl.data[cluster_center_index + ClusterSizeX]);
// dx = -1
// dy = 0
break;
case (corner::cBottomRight):
calculate_eta2(eta, cl.data[cluster_center_index],
cl.data[cluster_center_index + 1],
cl.data[cluster_center_index],
cl.data[cluster_center_index + ClusterSizeX]);
// dx = 0
// dy = 0
break;
}
eta.c = c; // TODO only supported for 2x2 and 3x3 clusters -> at least no
// underyling enum class
eta.c = c;
return eta;
}
// TODO! Look up eta2 calculation - photon center should be top right corner
template <typename T>
Eta2<T> calculate_eta2(const Cluster<T, 2, 2, int16_t> &cl) {
/**
* @brief Calculate the eta2 values for a generic sized cluster and return them
* in a Eta2 struct containing etay, etax and the index (as corner) of the
* respective 2x2 subcluster relative to the cluster center.
*/
template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
typename CoordType>
Eta2<T> calculate_full_eta2(
const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType> &cl) {
static_assert(ClusterSizeX > 1 && ClusterSizeY > 1);
Eta2<T> eta{};
if ((cl.data[0] + cl.data[1]) != 0)
eta.x = static_cast<double>(cl.data[1]) /
(cl.data[0] + cl.data[1]); // between (0,1) the closer to zero
// left value probably larger
if ((cl.data[0] + cl.data[2]) != 0)
eta.y = static_cast<double>(cl.data[2]) /
(cl.data[0] + cl.data[2]); // between (0,1) the closer to zero
// bottom value probably larger
constexpr size_t cluster_center_index =
(ClusterSizeX / 2) + (ClusterSizeY / 2) * ClusterSizeX;
auto max_sum = cl.max_sum_2x2();
eta.sum = max_sum.sum;
corner c = max_sum.index;
// subcluster top right from center
switch (c) {
case (corner::cTopLeft):
if (eta.sum != 0) {
eta.x = static_cast<double>(
cl.data[cluster_center_index] +
cl.data[cluster_center_index - ClusterSizeX]) /
static_cast<double>(eta.sum);
eta.y = static_cast<double>(cl.data[cluster_center_index - 1] +
cl.data[cluster_center_index]) /
static_cast<double>(eta.sum);
}
// dx = -1
// dy = -1
break;
case (corner::cTopRight):
if (eta.sum != 0) {
eta.x = static_cast<double>(
cl.data[cluster_center_index + 1] +
cl.data[cluster_center_index - ClusterSizeX + 1]) /
static_cast<double>(eta.sum);
eta.y = static_cast<double>(cl.data[cluster_center_index] +
cl.data[cluster_center_index + 1]) /
static_cast<double>(eta.sum);
}
// dx = 0
// dy = -1
break;
case (corner::cBottomLeft):
if (eta.sum != 0) {
eta.x = static_cast<double>(
cl.data[cluster_center_index] +
cl.data[cluster_center_index + ClusterSizeX]) /
static_cast<double>(eta.sum);
eta.y = static_cast<double>(
cl.data[cluster_center_index + ClusterSizeX] +
cl.data[cluster_center_index + ClusterSizeX - 1]) /
static_cast<double>(eta.sum);
}
// dx = -1
// dy = 0
break;
case (corner::cBottomRight):
if (eta.sum != 0) {
eta.x = static_cast<double>(
cl.data[cluster_center_index + 1] +
cl.data[cluster_center_index + ClusterSizeX + 1]) /
static_cast<double>(eta.sum);
eta.y = static_cast<double>(
cl.data[cluster_center_index + ClusterSizeX] +
cl.data[cluster_center_index + ClusterSizeX + 1]) /
static_cast<double>(eta.sum);
}
// dx = 0
// dy = 0
break;
}
eta.c = c;
return eta;
}
template <typename T>
Eta2<T> calculate_eta2(const Cluster<T, 2, 2, uint16_t> &cl) {
Eta2<T> eta{};
// TODO: maybe have as member function of cluster
const uint8_t photon_hit_index =
std::max_element(cl.data.begin(), cl.data.end()) - cl.data.begin();
eta.c = static_cast<corner>(3 - photon_hit_index);
switch (eta.c) {
case corner::cTopLeft:
calculate_eta2(eta, cl.data[2], cl.data[3], cl.data[1], cl.data[3]);
break;
case corner::cTopRight:
calculate_eta2(eta, cl.data[2], cl.data[3], cl.data[0], cl.data[2]);
break;
case corner::cBottomLeft:
calculate_eta2(eta, cl.data[0], cl.data[1], cl.data[1], cl.data[3]);
break;
case corner::cBottomRight:
calculate_eta2(eta, cl.data[0], cl.data[1], cl.data[0], cl.data[2]);
break;
}
eta.sum = cl.sum();
return eta;
}
// calculates Eta3 for 3x3 cluster based on code from analyze_cluster
// TODO only supported for 3x3 Clusters
template <typename T> Eta2<T> calculate_eta3(const Cluster<T, 3, 3> &cl) {
template <typename T>
Eta2<T> calculate_full_eta2(const Cluster<T, 2, 2, uint16_t> &cl) {
Eta2<T> eta{};
T sum = 0;
eta.sum = cl.sum();
std::for_each(std::begin(cl.data), std::end(cl.data),
[&sum](T x) { sum += x; });
const uint8_t photon_hit_index =
std::max_element(cl.data.begin(), cl.data.end()) - cl.data.begin();
eta.sum = sum;
eta.c = static_cast<corner>(3 - photon_hit_index);
if (eta.sum != 0) {
eta.x = static_cast<double>(cl.data[1] + cl.data[3]) /
static_cast<double>(eta.sum);
eta.y = static_cast<double>(cl.data[2] + cl.data[3]) /
static_cast<double>(eta.sum);
}
return eta;
}
// TODO generalize
template <typename T>
Eta2<T> calculate_eta2(const Cluster<T, 1, 2, uint16_t> &cl) {
Eta2<T> eta{};
eta.x = 0;
eta.y = static_cast<double>(cl.data[1]) / cl.data[0];
eta.sum = cl.sum();
}
template <typename T>
Eta2<T> calculate_eta2(const Cluster<T, 2, 1, uint16_t> &cl) {
Eta2<T> eta{};
eta.x = static_cast<double>(cl.data[1]) / cl.data[0];
eta.y = 0;
eta.sum = cl.sum();
}
/**
* @brief calculates cross Eta3 for 3x3 cluster
* cross Eta3 calculates the eta by taking into account only the cross pixels
* {top, bottom, left, right, center}
*/
template <typename T, typename CoordType = uint16_t>
Eta2<T> calculate_cross_eta3(const Cluster<T, 3, 3, CoordType> &cl) {
Eta2<T> eta{};
T photon_energy = cl.sum();
eta.sum = photon_energy;
if ((cl.data[3] + cl.data[4] + cl.data[5]) != 0)
eta.x = static_cast<double>(-cl.data[3] + cl.data[3 + 2]) /
eta.x =
static_cast<double>(-cl.data[3] + cl.data[3 + 2]) /
(cl.data[3] + cl.data[4] + cl.data[5]); // (-1,1)
static_cast<double>(cl.data[3] + cl.data[4] + cl.data[5]); // (-1,1)
if ((cl.data[1] + cl.data[4] + cl.data[7]) != 0)
eta.y = static_cast<double>(-cl.data[1] + cl.data[2 * 3 + 1]) /
(cl.data[1] + cl.data[4] + cl.data[7]);
static_cast<double>(cl.data[1] + cl.data[4] + cl.data[7]);
return eta;
}
template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
typename CoordType = uint16_t>
Eta2<T> calculate_cross_eta3(
const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType> &cl) {
static_assert(ClusterSizeX > 2 && ClusterSizeY > 2,
"calculate_eta3 only defined for clusters larger than 2x2");
if constexpr (ClusterSizeX != 3 || ClusterSizeY != 3) {
auto reduced_cluster = reduce_cluster_to_3x3(cl);
return calculate_cross_eta3(reduced_cluster);
} else {
return calculate_cross_eta3(cl);
}
}
/**
* @brief calculates Eta3 for 3x3 cluster
* It calculates the eta by taking into account all pixels in the 3x3 cluster
*/
template <typename T, typename CoordType = uint16_t>
Eta2<T> calculate_eta3(const Cluster<T, 3, 3, CoordType> &cl) {
Eta2<T> eta{};
T photon_energy = cl.sum();
eta.sum = photon_energy;
// TODO: how do we handle potential arithmetic overflows? - T could be
// uint16
if (photon_energy != 0) {
std::array<T, 2> column_sums{
static_cast<T>(cl.data[0] + cl.data[3] + cl.data[6]),
static_cast<T>(cl.data[2] + cl.data[5] + cl.data[8])};
eta.x = static_cast<double>(-column_sums[0] + column_sums[1]) /
static_cast<double>(photon_energy);
std::array<T, 2> row_sums{
static_cast<T>(cl.data[0] + cl.data[1] + cl.data[2]),
static_cast<T>(cl.data[6] + cl.data[7] + cl.data[8])};
eta.y = static_cast<double>(-row_sums[0] + row_sums[1]) /
static_cast<double>(photon_energy);
}
return eta;
}
template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
typename CoordType = uint16_t>
Eta2<T>
calculate_eta3(const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType> &cl) {
static_assert(ClusterSizeX > 2 && ClusterSizeY > 2,
"calculate_eta3 only defined for clusters larger than 2x2");
if constexpr (ClusterSizeX != 3 || ClusterSizeY != 3) {
auto reduced_cluster = reduce_cluster_to_3x3(cl);
return calculate_eta3(reduced_cluster);
} else {
return calculate_eta3(cl);
}
}
} // namespace aare
include/aare/CircularFifo.hpp
+44 -3
View File
@@ -1,10 +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"
@@ -16,7 +19,7 @@ template <class ItemType> class CircularFifo {
aare::ProducerConsumerQueue<ItemType> filled_slots;
public:
CircularFifo() : CircularFifo(100){};
CircularFifo() : CircularFifo(100) {};
CircularFifo(uint32_t size)
: fifo_size(size), free_slots(size + 1), filled_slots(size + 1) {
@@ -28,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;
@@ -46,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))
@@ -74,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/Cluster.hpp
+113 -117
View File
@@ -1,3 +1,4 @@
// SPDX-License-Identifier: MPL-2.0
/************************************************
* @file Cluster.hpp
@@ -8,7 +9,7 @@
#pragma once
#include "logger.hpp"
#include "defs.hpp"
#include <algorithm>
#include <array>
#include <cstdint>
@@ -18,8 +19,12 @@
namespace aare {
// requires clause c++20 maybe update
/**
* @brief Cluster struct
*/
template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
typename CoordType = int16_t>
typename CoordType = uint16_t>
struct Cluster {
static_assert(std::is_arithmetic_v<T>, "T needs to be an arithmetic type");
@@ -28,8 +33,11 @@ struct Cluster {
static_assert(ClusterSizeX > 0 && ClusterSizeY > 0,
"Cluster sizes must be bigger than zero");
/// @brief Cluster center x coordinate (in pixel coordinates)
CoordType x;
/// @brief Cluster center y coordinate (in pixel coordinates)
CoordType y;
/// @brief Cluster data stored in row-major order starting from top-left
std::array<T, ClusterSizeX * ClusterSizeY> data;
static constexpr uint8_t cluster_size_x = ClusterSizeX;
@@ -37,9 +45,18 @@ struct Cluster {
using value_type = T;
using coord_type = CoordType;
/**
* @brief Sum of all elements in the cluster
*/
T sum() const { return std::accumulate(data.begin(), data.end(), T{}); }
std::pair<T, int> max_sum_2x2() const {
// TODO: handle 1 dimensional clusters
/**
* @brief Find the highest-sum center-adjacent 2x2 subcluster.
* @return Sum and corner index of the selected 2x2 subcluster relative to
* the cluster center
*/
Sum_index_pair<T, corner> max_sum_2x2() const {
if constexpr (cluster_size_x == 3 && cluster_size_y == 3) {
std::array<T, 4> sum_2x2_subclusters;
@@ -50,96 +67,122 @@ struct Cluster {
int index = std::max_element(sum_2x2_subclusters.begin(),
sum_2x2_subclusters.end()) -
sum_2x2_subclusters.begin();
return std::make_pair(sum_2x2_subclusters[index], index);
return Sum_index_pair<T, corner>{sum_2x2_subclusters[index],
corner{index}};
} else if constexpr (cluster_size_x == 2 && cluster_size_y == 2) {
return std::make_pair(data[0] + data[1] + data[2] + data[3], 0);
return Sum_index_pair<T, corner>{
data[0] + data[1] + data[2] + data[3], corner{0}};
} else {
constexpr size_t num_2x2_subclusters =
(ClusterSizeX - 1) * (ClusterSizeY - 1);
constexpr size_t cluster_center_index =
(ClusterSizeX / 2) + (ClusterSizeY / 2) * ClusterSizeX;
std::array<T, num_2x2_subclusters> sum_2x2_subcluster;
for (size_t i = 0; i < ClusterSizeY - 1; ++i) {
for (size_t j = 0; j < ClusterSizeX - 1; ++j)
sum_2x2_subcluster[i * (ClusterSizeX - 1) + j] =
data[i * ClusterSizeX + j] +
data[i * ClusterSizeX + j + 1] +
data[(i + 1) * ClusterSizeX + j] +
data[(i + 1) * ClusterSizeX + j + 1];
std::array<T, 4> sum_2x2_subcluster{0};
// subcluster top left from center
sum_2x2_subcluster[0] =
data[cluster_center_index] + data[cluster_center_index - 1] +
data[cluster_center_index - ClusterSizeX] +
data[cluster_center_index - 1 - ClusterSizeX];
// subcluster top right from center
if (ClusterSizeX > 2) {
sum_2x2_subcluster[1] =
data[cluster_center_index] +
data[cluster_center_index + 1] +
data[cluster_center_index - ClusterSizeX] +
data[cluster_center_index - ClusterSizeX + 1];
}
// subcluster bottom left from center
if (ClusterSizeY > 2) {
sum_2x2_subcluster[2] =
data[cluster_center_index] +
data[cluster_center_index - 1] +
data[cluster_center_index + ClusterSizeX] +
data[cluster_center_index + ClusterSizeX - 1];
}
// subcluster bottom right from center
if (ClusterSizeX > 2 && ClusterSizeY > 2) {
sum_2x2_subcluster[3] =
data[cluster_center_index] +
data[cluster_center_index + 1] +
data[cluster_center_index + ClusterSizeX] +
data[cluster_center_index + ClusterSizeX + 1];
}
int index = std::max_element(sum_2x2_subcluster.begin(),
sum_2x2_subcluster.end()) -
sum_2x2_subcluster.begin();
return std::make_pair(sum_2x2_subcluster[index], index);
return Sum_index_pair<T, corner>{sum_2x2_subcluster[index],
corner{index}};
}
}
};
/**
* @brief Reduce a cluster to a 2x2 cluster by selecting the 2x2 block with the
* highest sum.
* @brief Reduce a cluster to its highest-sum center-adjacent 2x2 block.
* @param c Cluster to reduce
* @return reduced cluster
* @return Reduced cluster with the input coordinates
* @note The cluster is filled using row major ordering starting at the top-left
*/
template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
typename CoordType = int16_t>
typename CoordType = uint16_t>
Cluster<T, 2, 2, CoordType>
reduce_to_2x2(const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType> &c) {
static_assert(ClusterSizeX >= 2 && ClusterSizeY >= 2,
"Cluster sizes must be at least 2x2 for reduction to 2x2");
// TODO maybe add sanity check and check that center is in max subcluster
Cluster<T, 2, 2, CoordType> result;
Cluster<T, 2, 2, CoordType> result{};
auto [sum, index] = c.max_sum_2x2();
int16_t cluster_center_index =
constexpr int16_t cluster_center_index =
(ClusterSizeX / 2) + (ClusterSizeY / 2) * ClusterSizeX;
int16_t index_bottom_left_max_2x2_subcluster =
(int(index / (ClusterSizeX - 1))) * ClusterSizeX +
index % (ClusterSizeX - 1);
int16_t index_top_left_max_2x2_subcluster = cluster_center_index;
switch (index) {
case corner::cTopLeft:
index_top_left_max_2x2_subcluster -= (ClusterSizeX + 1);
break;
case corner::cTopRight:
index_top_left_max_2x2_subcluster -= ClusterSizeX;
break;
case corner::cBottomLeft:
index_top_left_max_2x2_subcluster -= 1;
break;
case corner::cBottomRight:
// no change needed
break;
}
result.x =
c.x + (index_bottom_left_max_2x2_subcluster - cluster_center_index) %
ClusterSizeX;
result.x = c.x;
result.y = c.y;
result.y =
c.y - (index_bottom_left_max_2x2_subcluster - cluster_center_index) /
ClusterSizeX;
result.data = {
c.data[index_bottom_left_max_2x2_subcluster],
c.data[index_bottom_left_max_2x2_subcluster + 1],
c.data[index_bottom_left_max_2x2_subcluster + ClusterSizeX],
c.data[index_bottom_left_max_2x2_subcluster + ClusterSizeX + 1]};
c.data[index_top_left_max_2x2_subcluster],
c.data[index_top_left_max_2x2_subcluster + 1],
c.data[index_top_left_max_2x2_subcluster + ClusterSizeX],
c.data[index_top_left_max_2x2_subcluster + ClusterSizeX + 1]};
return result;
}
template <typename T>
Cluster<T, 2, 2, int16_t> reduce_to_2x2(const Cluster<T, 3, 3, int16_t> &c) {
Cluster<T, 2, 2, int16_t> result;
Cluster<T, 2, 2, uint16_t> reduce_to_2x2(const Cluster<T, 3, 3, uint16_t> &c) {
Cluster<T, 2, 2, uint16_t> result{};
auto [s, i] = c.max_sum_2x2();
result.x = c.x;
result.y = c.y;
switch (i) {
case 0:
result.x = c.x - 1;
result.y = c.y + 1;
case corner::cTopLeft:
result.data = {c.data[0], c.data[1], c.data[3], c.data[4]};
break;
case 1:
result.x = c.x;
result.y = c.y + 1;
case corner::cTopRight:
result.data = {c.data[1], c.data[2], c.data[4], c.data[5]};
break;
case 2:
result.x = c.x - 1;
result.y = c.y;
case corner::cBottomLeft:
result.data = {c.data[3], c.data[4], c.data[6], c.data[7]};
break;
case 3:
result.x = c.x;
result.y = c.y;
case corner::cBottomRight:
result.data = {c.data[4], c.data[5], c.data[7], c.data[8]};
break;
}
@@ -147,88 +190,41 @@ Cluster<T, 2, 2, int16_t> reduce_to_2x2(const Cluster<T, 3, 3, int16_t> &c) {
return result;
}
template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
typename CoordType = int16_t>
inline std::pair<T, uint16_t>
max_3x3_sum(const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType> &cluster) {
if constexpr (ClusterSizeX == 3 && ClusterSizeY == 3) {
return std::make_pair(cluster.sum(), 0);
} else {
size_t index = 0;
T max_3x3_subcluster_sum = 0;
for (size_t i = 0; i < ClusterSizeY - 2; ++i) {
for (size_t j = 0; j < ClusterSizeX - 2; ++j) {
T sum = cluster.data[i * ClusterSizeX + j] +
cluster.data[i * ClusterSizeX + j + 1] +
cluster.data[i * ClusterSizeX + j + 2] +
cluster.data[(i + 1) * ClusterSizeX + j] +
cluster.data[(i + 1) * ClusterSizeX + j + 1] +
cluster.data[(i + 1) * ClusterSizeX + j + 2] +
cluster.data[(i + 2) * ClusterSizeX + j] +
cluster.data[(i + 2) * ClusterSizeX + j + 1] +
cluster.data[(i + 2) * ClusterSizeX + j + 2];
if (sum > max_3x3_subcluster_sum) {
max_3x3_subcluster_sum = sum;
index = i * (ClusterSizeX - 2) + j;
}
}
}
return std::make_pair(max_3x3_subcluster_sum, index);
}
}
/**
* @brief Reduce a cluster to a 3x3 cluster by selecting the 3x3 block with the
* highest sum.
* @brief Reduce a cluster to the 3x3 block around its center index.
* @param c Cluster to reduce
* @return reduced cluster
* @pre ClusterSizeX and ClusterSizeY must both be at least 3, and at least one
* must be greater than 3.
* @return Reduced cluster with the input coordinates are preserved.
*/
template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
typename CoordType = int16_t>
Cluster<T, 3, 3, CoordType>
reduce_to_3x3(const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType> &c) {
static_assert(ClusterSizeX >= 3 && ClusterSizeY >= 3,
"Cluster sizes must be at least 3x3 for reduction to 3x3");
static_assert(ClusterSizeX >= 3 && ClusterSizeY >= 3 &&
(ClusterSizeX > 3 || ClusterSizeY > 3),
"Cluster sizes must both be at least 3, and at least one "
"must be greater than 3 for reduction to 3x3");
Cluster<T, 3, 3, CoordType> result;
// TODO maybe add sanity check and check that center is in max subcluster
auto [sum, index] = max_3x3_sum(c);
Cluster<T, 3, 3, CoordType> result{};
int16_t cluster_center_index =
(ClusterSizeX / 2) + (ClusterSizeY / 2) * ClusterSizeX;
int16_t index_center_max_3x3_subcluster =
(int(index / (ClusterSizeX - 2))) * ClusterSizeX + ClusterSizeX +
index % (ClusterSizeX - 2) + 1;
result.x = c.x;
result.y = c.y;
int16_t index_3x3_subcluster_cluster_center =
int((cluster_center_index - 1 - ClusterSizeX) / ClusterSizeX) *
(ClusterSizeX - 2) +
(cluster_center_index - 1 - ClusterSizeX) % ClusterSizeX;
result.data = {c.data[cluster_center_index - ClusterSizeX - 1],
c.data[cluster_center_index - ClusterSizeX],
c.data[cluster_center_index - ClusterSizeX + 1],
c.data[cluster_center_index - 1],
c.data[cluster_center_index],
c.data[cluster_center_index + 1],
c.data[cluster_center_index + ClusterSizeX - 1],
c.data[cluster_center_index + ClusterSizeX],
c.data[cluster_center_index + ClusterSizeX + 1]};
result.x =
c.x + (index % (ClusterSizeX - 2) -
(index_3x3_subcluster_cluster_center % (ClusterSizeX - 2)));
result.y =
c.y - (index / (ClusterSizeX - 2) -
(index_3x3_subcluster_cluster_center / (ClusterSizeX - 2)));
result.data = {c.data[index_center_max_3x3_subcluster - ClusterSizeX - 1],
c.data[index_center_max_3x3_subcluster - ClusterSizeX],
c.data[index_center_max_3x3_subcluster - ClusterSizeX + 1],
c.data[index_center_max_3x3_subcluster - 1],
c.data[index_center_max_3x3_subcluster],
c.data[index_center_max_3x3_subcluster + 1],
c.data[index_center_max_3x3_subcluster + ClusterSizeX - 1],
c.data[index_center_max_3x3_subcluster + ClusterSizeX],
c.data[index_center_max_3x3_subcluster + ClusterSizeX + 1]};
return result;
}
include/aare/ClusterCollector.hpp
+35 -13
View File
@@ -1,58 +1,80 @@
// SPDX-License-Identifier: MPL-2.0
#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"
#include "aare/defs.hpp"
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
+325 -169
View File
@@ -1,47 +1,41 @@
// SPDX-License-Identifier: MPL-2.0
#pragma once
#include "aare/Cluster.hpp"
#include "aare/ClusterVector.hpp"
#include "aare/FilePtr.hpp"
#include "aare/GainMap.hpp"
#include "aare/NDArray.hpp"
#include "aare/defs.hpp"
#include "aare/ROI.hpp"
#include "aare/logger.hpp"
#include <filesystem>
#include <fstream>
#include <optional>
#include <stdexcept>
#include <utility>
namespace aare {
/*
Binary cluster file. Expects data to be laid out as:
int32_t frame_number
uint32_t number_of_clusters
int16_t x, int16_t y, int32_t data[9] x number_of_clusters
int32_t frame_number
uint32_t number_of_clusters
....
*/
// TODO: change to support any type of clusters, e.g. header line with
// clsuter_size_x, cluster_size_y,
/**
* @brief Class to read and write cluster files
* Expects data to be laid out as:
* @brief Read and write legacy binary cluster files.
*
* Each frame is stored as:
*
* int32_t frame_number
* uint32_t number_of_clusters
* int16_t x, int16_t y, int32_t data[9] * number_of_clusters
* int32_t frame_number
* uint32_t number_of_clusters
* etc.
* ClusterType clusters[number_of_clusters]
*
* The format stores clusters as their native in-memory representation and has
* no metadata describing the cluster dimensions, value type, coordinate type,
* padding, or byte order. Readers must therefore use the same ClusterType and
* a compatible platform ABI as the writer.
*/
template <typename ClusterType,
typename Enable = std::enable_if_t<is_cluster_v<ClusterType>>>
class ClusterFile {
FILE *fp{};
const std::string m_filename{};
FilePtr m_fp;
std::string m_filename{};
uint32_t m_num_left{}; /*Number of photons left in frame*/
size_t m_chunk_size{}; /*Number of clusters to read at a time*/
std::string m_mode; /*Mode to open the file in*/
@@ -53,48 +47,150 @@ class ClusterFile {
public:
/**
* @brief Construct a new Cluster File object
* @param fname path to the file
* @param chunk_size number of clusters to read at a time when iterating
* over the file
* @param mode mode to open the file in. "r" for reading, "w" for writing,
* "a" for appending
* @throws std::runtime_error if the file could not be opened
* @brief Single-pass range over frames or chunks of selected clusters.
*
* The range borrows the file and shares its current position with all
* other reads. Use only one traversal at a time. The file must outlive the
* range and its iterators and must not be moved while they are in use.
* Constructing a range does not read; begin() reads the first result.
* Calling begin() again starts at the file's then-current position.
*/
template <bool ByFrame> class ReadRange {
ClusterFile *m_file;
size_t m_chunk_size;
friend class ClusterFile;
ReadRange(ClusterFile &file, size_t chunk_size)
: m_file(&file), m_chunk_size(chunk_size) {}
public:
struct Sentinel {};
/**
* @brief Move-only iterator for C++17 range-based loops.
*
* References to the current vector are valid until the iterator is
* advanced or destroyed. Move the vector out to retain its storage.
* Frame iteration reuses storage when the vector is not moved out.
* This is not a copyable STL input iterator.
*/
class Iterator {
ClusterFile *m_file;
size_t m_chunk_size;
ClusterVector<ClusterType> m_clusters{0};
friend class ReadRange;
Iterator(ClusterFile &file, size_t chunk_size)
: m_file(&file), m_chunk_size(chunk_size) {
++*this;
}
public:
using value_type = ClusterVector<ClusterType>;
Iterator(const Iterator &) = delete;
Iterator &operator=(const Iterator &) = delete;
Iterator(Iterator &&) noexcept = default;
Iterator &operator=(Iterator &&) noexcept = default;
~Iterator() = default;
value_type &operator*() { return m_clusters; }
value_type *operator->() { return &m_clusters; }
Iterator &operator++() {
if (m_file) {
if constexpr (ByFrame) {
if (!m_file->read_frame(m_clusters)) {
m_file = nullptr;
}
} else {
m_clusters = m_file->read_clusters(m_chunk_size);
if (m_clusters.empty()) {
m_file = nullptr;
}
}
}
return *this;
}
void operator++(int) { ++*this; }
friend bool operator==(const Iterator &it, Sentinel /*end*/) {
return it.m_file == nullptr;
}
friend bool operator!=(const Iterator &it, Sentinel end) {
return !(it == end);
}
friend bool operator==(Sentinel end, const Iterator &it) {
return it == end;
}
friend bool operator!=(Sentinel end, const Iterator &it) {
return !(it == end);
}
};
Iterator begin() const { return Iterator(*m_file, m_chunk_size); }
Sentinel end() const { return {}; }
};
using FrameRange = ReadRange<true>;
using ChunkRange = ReadRange<false>;
/**
* @brief Open a cluster file.
* @param fname Path to the file.
* @param chunk_size Maximum number of selected clusters per chunk step.
* @param mode File mode: "r" to read, "w" to truncate and write, or "a"
* to append.
* @throws std::runtime_error If the mode is unsupported or the file cannot
* be opened.
*/
ClusterFile(const std::filesystem::path &fname, size_t chunk_size = 1000,
const std::string &mode = "r")
: m_filename(fname.string()), m_chunk_size(chunk_size), m_mode(mode) {
if (mode == "r") {
fp = fopen(m_filename.c_str(), "rb");
if (!fp) {
throw std::runtime_error("Could not open file for reading: " +
m_filename);
}
} else if (mode == "w") {
fp = fopen(m_filename.c_str(), "wb");
if (!fp) {
throw std::runtime_error("Could not open file for writing: " +
m_filename);
}
} else if (mode == "a") {
fp = fopen(m_filename.c_str(), "ab");
if (!fp) {
throw std::runtime_error("Could not open file for appending: " +
m_filename);
}
} else {
throw std::runtime_error("Unsupported mode: " + mode);
}
: m_filename(fname.string()), m_chunk_size(chunk_size) {
open(mode);
}
~ClusterFile() { close(); }
/**
* @brief Iterate over complete frames from the current file position.
*
* Empty and fully filtered frames are yielded with their stored frame
* numbers. Reading stops only at a clean end of file. Advancing uses
* read_frame(), including its error for a prior partial-frame read.
*/
FrameRange frames() & { return FrameRange(*this, 0); }
FrameRange frames() && = delete;
/** @brief Iterate using the chunk size supplied to the constructor. */
ChunkRange chunks() & { return chunks(m_chunk_size); }
ChunkRange chunks() && = delete;
/**
* @brief Read n_clusters clusters from the file discarding
* frame numbers. If EOF is reached the returned vector will
* have less than n_clusters clusters
* @brief Iterate over chunks from the current file position.
* @param chunk_size Maximum number of selected clusters per step.
* @throws std::invalid_argument If chunk_size is zero.
* @note Chunks can span frames; their frame numbers are not per-cluster
* metadata. Advancing uses read_clusters() and propagates its errors.
* Only the final chunk can be short. Empty chunks are not yielded.
*/
ChunkRange chunks(size_t chunk_size) & {
if (chunk_size == 0) {
throw std::invalid_argument("Chunk size must be greater than zero");
}
return ChunkRange(*this, chunk_size);
}
ChunkRange chunks(size_t) && = delete;
/**
* @brief Read up to n_clusters without preserving frame boundaries.
* @param n_clusters Maximum number of selected clusters to return.
* @return A cluster vector that may contain fewer clusters at a clean end
* of file.
* @note The returned vector may combine data from several frames, so its
* frame number must not be used as per-cluster metadata.
* @throws std::runtime_error If the file is not open for reading, an I/O
* error occurs, or an incomplete frame header or cluster record is read.
*/
ClusterVector<ClusterType> read_clusters(size_t n_clusters) {
if (m_mode != "r") {
@@ -108,63 +204,116 @@ class ClusterFile {
}
/**
* @brief Read a single frame from the file and return the
* clusters. The cluster vector will have the frame number
* set.
* @throws std::runtime_error if the file is not opened for
* reading or the file pointer not at the beginning of a
* frame
* @brief Read the next complete frame.
* @return The selected clusters with the stored frame number set, or
* std::nullopt at a clean end of file before the next frame.
* @throws std::runtime_error If the file is not open for reading, a prior
* partial-frame read left clusters unread, or the frame is incomplete.
* @note A complete frame produces an engaged optional even when it contains
* no clusters or all of its clusters are removed by the configured filters.
*/
ClusterVector<ClusterType> read_frame() {
std::optional<ClusterVector<ClusterType>> read_frame() {
ClusterVector<ClusterType> clusters(0);
if (!read_frame(clusters)) {
return std::nullopt;
}
return std::optional<ClusterVector<ClusterType>>{std::move(clusters)};
}
/**
* @brief Read the next complete frame into an existing cluster vector.
* @param clusters Destination whose storage is reused when large enough.
* Existing clusters are replaced after a frame header is read and remain
* unchanged at a clean end of file.
* @return true when a complete frame was read, or false at a clean end of
* file before the next frame.
* @throws std::runtime_error If the file is not open for reading, a prior
* partial-frame read left clusters unread, or a frame is incomplete.
* @note A complete frame is a successful read even when it contains no
* clusters or all of its clusters are removed by the configured filters.
*/
bool read_frame(ClusterVector<ClusterType> &clusters) {
if (m_mode != "r") {
throw std::runtime_error(LOCATION + "File not opened for reading");
}
if (m_noise_map || m_roi) {
return read_frame_with_cut();
} else {
return read_frame_without_cut();
if (m_num_left) {
throw std::runtime_error(
LOCATION + "There are still clusters left in the last frame");
}
if (m_noise_map || m_roi) {
return read_frame_with_cut(clusters);
}
return read_frame_without_cut(clusters);
}
/**
* @brief Write one frame to the file.
* @param clusters Clusters to write, including their frame number.
* @throws std::runtime_error If the file is not open for writing or any
* part of the frame cannot be written completely.
*/
void write_frame(const ClusterVector<ClusterType> &clusters) {
if (m_mode != "w" && m_mode != "a") {
throw std::runtime_error("File not opened for writing");
}
int32_t frame_number = clusters.frame_number();
fwrite(&frame_number, sizeof(frame_number), 1, fp);
if (fwrite(&frame_number, sizeof(frame_number), 1, m_fp.get()) != 1) {
throw std::runtime_error(LOCATION + "Could not write frame number");
}
uint32_t n_clusters = clusters.size();
fwrite(&n_clusters, sizeof(n_clusters), 1, fp);
fwrite(clusters.data(), clusters.item_size(), clusters.size(), fp);
if (fwrite(&n_clusters, sizeof(n_clusters), 1, m_fp.get()) != 1) {
throw std::runtime_error(LOCATION +
"Could not write number of clusters");
}
if (fwrite(clusters.data(), clusters.item_size(), clusters.size(),
m_fp.get()) != clusters.size()) {
throw std::runtime_error(LOCATION + "Could not write clusters");
}
}
/**
* @brief Return the chunk size
* @brief Return the default number of selected clusters per chunk step.
*/
size_t chunk_size() const { return m_chunk_size; }
/**
* @brief Set the region of interest to use when reading
* clusters. If set only clusters within the ROI will be
* read.
* @brief Estimate the number of clusters in the file from its size.
*
* Frame-header bytes are included in the estimate, so it may exceed the
* actual number of clusters. The file position is not changed.
*/
size_t estimate_n_clusters() const {
return std::filesystem::file_size(m_filename) / sizeof(ClusterType);
}
/**
* @brief Select clusters by their center coordinate when reading.
* @param roi Half-open region of interest: [xmin, xmax) x [ymin, ymax).
*/
void set_roi(ROI roi) { m_roi = roi; }
/**
* @brief Set the noise map to use when reading clusters. If
* set clusters below the noise level will be discarded.
* Selection criteria one of: Central pixel above noise,
* highest 2x2 sum above 2 * noise, total sum above 3 *
* noise.
* @brief Discard clusters that do not pass the noise thresholds.
* @param noise_map Per-pixel noise indexed as [y, x]. The map is copied.
* A cluster is retained only when its central pixel exceeds the local
* noise, its highest 2x2 sum exceeds twice the noise, and its total sum
* exceeds three times the noise.
* @warning The map must cover every cluster center coordinate in the file.
*/
void set_noise_map(const NDView<int32_t, 2> noise_map) {
m_noise_map = NDArray<int32_t, 2>(noise_map);
}
/**
* @brief Set the gain map to use when reading clusters. If set the gain map
* will be applied to the clusters that pass ROI and noise_map selection.
* The gain map is expected to be in ADU/energy.
* @brief Apply a gain map to clusters selected while reading.
* @param gain_map Per-pixel gain in ADU/energy, indexed as [y, x]. The map
* is copied and inverted internally.
* @note Clusters whose complete footprint extends beyond the gain map are
* retained with all cluster data values set to zero.
*/
void set_gain_map(const NDView<double, 2> gain_map) {
m_gain_map = InvertedGainMap(gain_map);
@@ -179,55 +328,53 @@ class ClusterFile {
}
/**
* @brief Close the file. If not closed the file will be
* closed in the destructor
* @brief Close the file.
*
* Calling close more than once is safe. The destructor closes an open file
* automatically.
*/
void close() {
if (fp) {
fclose(fp);
fp = nullptr;
}
m_fp = FilePtr{};
m_mode = "";
}
/**
* @brief Return the current byte position in the file.
* @throws std::runtime_error If the file is closed or its position cannot
* be determined.
*/
int64_t tell() {
if (!m_fp.get()) {
throw std::runtime_error(LOCATION + "File not opened");
}
return m_fp.tell();
}
private:
/** @brief Open the file in specific mode
*
*/
void open(const std::string &mode) {
if (fp) {
close();
}
close();
if (mode == "r") {
fp = fopen(m_filename.c_str(), "rb");
if (!fp) {
throw std::runtime_error("Could not open file for reading: " +
m_filename);
}
m_fp = FilePtr(m_filename, "rb");
m_mode = "r";
} else if (mode == "w") {
fp = fopen(m_filename.c_str(), "wb");
if (!fp) {
throw std::runtime_error("Could not open file for writing: " +
m_filename);
}
m_fp = FilePtr(m_filename, "wb");
m_mode = "w";
} else if (mode == "a") {
fp = fopen(m_filename.c_str(), "ab");
if (!fp) {
throw std::runtime_error("Could not open file for appending: " +
m_filename);
}
m_fp = FilePtr(m_filename, "ab");
m_mode = "a";
} else {
throw std::runtime_error("Unsupported mode: " + mode);
}
}
private:
ClusterVector<ClusterType> read_clusters_with_cut(size_t n_clusters);
ClusterVector<ClusterType> read_clusters_without_cut(size_t n_clusters);
ClusterVector<ClusterType> read_frame_with_cut();
ClusterVector<ClusterType> read_frame_without_cut();
bool read_frame_header(int32_t &frame_number, uint32_t &n_clusters);
bool read_frame_with_cut(ClusterVector<ClusterType> &clusters);
bool read_frame_without_cut(ClusterVector<ClusterType> &clusters);
bool is_selected(ClusterType &cl);
ClusterType read_one_cluster();
};
@@ -257,25 +404,31 @@ ClusterFile<ClusterType, Enable>::read_clusters_without_cut(size_t n_clusters) {
} else {
nn = nph;
}
nph_read += fread((buf + nph_read), clusters.item_size(), nn, fp);
const auto count =
fread((buf + nph_read), clusters.item_size(), nn, m_fp.get());
if (count != nn) {
throw std::runtime_error(LOCATION + "Could not read clusters");
}
nph_read += count;
m_num_left = nph - nn; // write back the number of photons left
}
if (nph_read < n_clusters) {
// keep on reading frames and photons until reaching n_clusters
while (fread(&iframe, sizeof(iframe), 1, fp)) {
while (read_frame_header(iframe, nph)) {
clusters.set_frame_number(iframe);
// read number of clusters in frame
if (fread(&nph, sizeof(nph), 1, fp)) {
if (nph > (n_clusters - nph_read))
nn = n_clusters - nph_read;
else
nn = nph;
if (nph > (n_clusters - nph_read))
nn = n_clusters - nph_read;
else
nn = nph;
nph_read +=
fread((buf + nph_read), clusters.item_size(), nn, fp);
m_num_left = nph - nn;
const auto count =
fread((buf + nph_read), clusters.item_size(), nn, m_fp.get());
if (count != nn) {
throw std::runtime_error(LOCATION + "Could not read clusters");
}
nph_read += count;
m_num_left = nph - nn;
if (nph_read >= n_clusters)
break;
}
@@ -315,16 +468,15 @@ ClusterFile<ClusterType, Enable>::read_clusters_with_cut(size_t n_clusters) {
}
int32_t frame_number = 0; // frame number needs to be 4 bytes!
while (fread(&frame_number, sizeof(frame_number), 1, fp)) {
if (fread(&m_num_left, sizeof(m_num_left), 1, fp)) {
clusters.set_frame_number(
frame_number); // cluster vector will hold the last
// frame number
while (m_num_left && clusters.size() < n_clusters) {
ClusterType c = read_one_cluster();
if (is_selected(c)) {
clusters.push_back(c);
}
uint32_t n_clusters_in_frame = 0;
while (read_frame_header(frame_number, n_clusters_in_frame)) {
m_num_left = n_clusters_in_frame;
clusters.set_frame_number(
frame_number); // cluster vector will hold the last frame number
while (m_num_left && clusters.size() < n_clusters) {
ClusterType c = read_one_cluster();
if (is_selected(c)) {
clusters.push_back(c);
}
}
@@ -342,7 +494,7 @@ ClusterFile<ClusterType, Enable>::read_clusters_with_cut(size_t n_clusters) {
template <typename ClusterType, typename Enable>
ClusterType ClusterFile<ClusterType, Enable>::read_one_cluster() {
ClusterType c;
auto rc = fread(&c, sizeof(c), 1, fp);
auto rc = fread(&c, sizeof(c), 1, m_fp.get());
if (rc != 1) {
throw std::runtime_error(LOCATION + "Could not read cluster");
}
@@ -351,66 +503,70 @@ ClusterType ClusterFile<ClusterType, Enable>::read_one_cluster() {
}
template <typename ClusterType, typename Enable>
ClusterVector<ClusterType>
ClusterFile<ClusterType, Enable>::read_frame_without_cut() {
if (m_mode != "r") {
throw std::runtime_error("File not opened for reading");
bool ClusterFile<ClusterType, Enable>::read_frame_header(int32_t &frame_number,
uint32_t &n_clusters) {
const auto frame_number_bytes =
fread(&frame_number, 1, sizeof(frame_number), m_fp.get());
if (frame_number_bytes == 0 && feof(m_fp.get()) && !ferror(m_fp.get())) {
return false;
}
if (m_num_left) {
throw std::runtime_error(
"There are still photons left in the last frame");
}
int32_t frame_number;
if (fread(&frame_number, sizeof(frame_number), 1, fp) != 1) {
throw std::runtime_error(LOCATION + "Could not read frame number");
if (frame_number_bytes != sizeof(frame_number)) {
if (ferror(m_fp.get())) {
throw std::runtime_error(LOCATION + "Error reading from file");
}
throw std::runtime_error(LOCATION + "Incomplete frame number");
}
int32_t n_clusters; // Saved as 32bit integer in the cluster file
if (fread(&n_clusters, sizeof(n_clusters), 1, fp) != 1) {
throw std::runtime_error(LOCATION +
"Could not read number of clusters");
const auto cluster_count_bytes =
fread(&n_clusters, 1, sizeof(n_clusters), m_fp.get());
if (cluster_count_bytes != sizeof(n_clusters)) {
if (ferror(m_fp.get())) {
throw std::runtime_error(LOCATION + "Error reading from file");
}
throw std::runtime_error(LOCATION + "Incomplete number of clusters");
}
return true;
}
template <typename ClusterType, typename Enable>
bool ClusterFile<ClusterType, Enable>::read_frame_without_cut(
ClusterVector<ClusterType> &clusters) {
int32_t frame_number;
uint32_t n_clusters;
if (!read_frame_header(frame_number, n_clusters)) {
return false;
}
LOG(logDEBUG1) << "Reading " << n_clusters << " clusters from frame "
<< frame_number;
ClusterVector<ClusterType> clusters(n_clusters);
clusters.set_frame_number(frame_number);
clusters.resize(n_clusters);
LOG(logDEBUG1) << "clusters.item_size(): " << clusters.item_size();
if (fread(clusters.data(), clusters.item_size(), n_clusters, fp) !=
if (fread(clusters.data(), clusters.item_size(), n_clusters, m_fp.get()) !=
static_cast<size_t>(n_clusters)) {
throw std::runtime_error(LOCATION + "Could not read clusters");
}
if (m_gain_map)
m_gain_map->apply_gain_map(clusters);
return clusters;
return true;
}
template <typename ClusterType, typename Enable>
ClusterVector<ClusterType>
ClusterFile<ClusterType, Enable>::read_frame_with_cut() {
if (m_mode != "r") {
throw std::runtime_error("File not opened for reading");
}
if (m_num_left) {
throw std::runtime_error(
"There are still photons left in the last frame");
}
bool ClusterFile<ClusterType, Enable>::read_frame_with_cut(
ClusterVector<ClusterType> &clusters) {
int32_t frame_number;
if (fread(&frame_number, sizeof(frame_number), 1, fp) != 1) {
throw std::runtime_error("Could not read frame number");
uint32_t n_clusters;
if (!read_frame_header(frame_number, n_clusters)) {
return false;
}
if (fread(&m_num_left, sizeof(m_num_left), 1, fp) != 1) {
throw std::runtime_error("Could not read number of clusters");
}
ClusterVector<ClusterType> clusters;
clusters.reserve(m_num_left);
m_num_left = n_clusters;
clusters.resize(0);
clusters.reserve(n_clusters);
clusters.set_frame_number(frame_number);
while (m_num_left) {
ClusterType c = read_one_cluster();
@@ -420,7 +576,7 @@ ClusterFile<ClusterType, Enable>::read_frame_with_cut() {
}
if (m_gain_map)
m_gain_map->apply_gain_map(clusters);
return clusters;
return true;
}
template <typename ClusterType, typename Enable>
@@ -438,8 +594,8 @@ bool ClusterFile<ClusterType, Enable>::is_selected(ClusterType &cl) {
if (m_noise_map) {
auto sum_1x1 = cl.data[cluster_center_index]; // central pixel
auto sum_2x2 = cl.max_sum_2x2().first; // highest sum of 2x2 subclusters
auto total_sum = cl.sum(); // sum of all pixels
auto sum_2x2 = cl.max_sum_2x2().sum; // highest sum of 2x2 subclusters
auto total_sum = cl.sum(); // sum of all pixels
auto noise =
(*m_noise_map)(cl.y, cl.x); // TODO! check if this is correct
include/aare/ClusterFileSink.hpp
+3 -1
View File
@@ -1,3 +1,4 @@
// SPDX-License-Identifier: MPL-2.0
#pragma once
#include <atomic>
#include <filesystem>
@@ -10,7 +11,8 @@
namespace aare {
template <typename ClusterType,
typename = std::enable_if_t<is_cluster_v<ClusterType>>>
typename = std::enable_if_t<is_cluster_v<ClusterType>>,
typename = std::enable_if_t<no_2x2_cluster<ClusterType>::value>>
class ClusterFileSink {
ProducerConsumerQueue<ClusterVector<ClusterType>> *m_source;
std::atomic<bool> m_stop_requested{false};
include/aare/ClusterFinder.hpp
+242 -98
View File
@@ -1,7 +1,9 @@
// SPDX-License-Identifier: MPL-2.0
#pragma once
#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"
@@ -10,54 +12,111 @@
namespace aare {
template <typename ClusterType,
typename = std::enable_if_t<is_cluster_v<ClusterType>>>
struct no_2x2_cluster {
constexpr static bool value =
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 FRAME_TYPE = uint16_t, typename PEDESTAL_TYPE = double,
typename = std::enable_if_t<no_2x2_cluster<ClusterType>::value>>
class ClusterFinder {
Shape<2> m_image_size;
const PEDESTAL_TYPE m_nSigma;
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.
* @param min_pedestal_samples Minimum number of pedestal frames to
* accumulate to get reasonable statistics
*/
ClusterFinder(Shape<2> image_size, PEDESTAL_TYPE nSigma = 5.0,
size_t capacity = 1000000)
size_t capacity = 1000000, size_t min_pedestal_samples = 1000)
: 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], min_pedestal_samples),
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;
<< ", nSigma: " << nSigma << ", capacity: " << capacity
<< ", min_pedestal_samples: " << min_pedestal_samples;
}
void push_pedestal_frame(NDView<FRAME_TYPE, 2> frame) {
m_pedestal.push(frame);
}
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 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 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();
}
/** @brief Return the current noise multiplier used for threshold
* calculation. */
PEDESTAL_TYPE get_nSigma() const { return m_nSigma; }
/**
* @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) {
@@ -68,94 +127,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_y; ic++) {
if (ix + ic >= 0 && ix + ic < frame.shape(1) &&
iy + ir >= 0 && iy + ir < frame.shape(0)) {
CT tmp =
static_cast<CT>(frame(iy + ir, ix + ic)) -
static_cast<CT>(
m_pedestal.mean(iy + ir, ix + ic));
cluster.data[i] =
tmp; // Watch for out of bounds access
}
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
+170 -42
View File
@@ -1,3 +1,4 @@
// SPDX-License-Identifier: MPL-2.0
#pragma once
#include <atomic>
#include <cstdint>
@@ -5,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 {
@@ -17,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(); }
};
/**
@@ -32,7 +65,8 @@ struct FrameWrapper {
* @tparam CT type of the cluster data
*/
template <typename ClusterType = Cluster<int32_t, 3, 3>,
typename FRAME_TYPE = uint16_t, typename PEDESTAL_TYPE = double>
typename FRAME_TYPE = uint16_t, typename PEDESTAL_TYPE = double,
typename = std::enable_if_t<no_2x2_cluster<ClusterType>::value>>
class ClusterFinderMT {
protected:
@@ -40,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};
@@ -63,50 +98,71 @@ 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);
m_output_queues[thread_id]->write(
cf->steal_clusters(realloc_same_capacity));
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(
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();
}
}
}
bool output_queues_are_empty() const {
for (auto &q : m_output_queues) {
if (!q->isEmpty()) {
return false;
}
}
return true;
}
/**
* @brief Collect all the clusters from the output queues and write them to
* the sink
*/
void collect() {
bool empty = true;
while (!m_stop_requested || !empty || !m_processing_threads_stopped) {
empty = true;
Backoff backoff;
while (!m_stop_requested || !output_queues_are_empty() ||
!m_processing_threads_stopped) {
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);
if (auto *front = queue->frontPtr(); front != nullptr) {
while (!m_sink.write(std::move(*front))) {
backoff.pause();
}
queue->popFront();
empty = false;
moved_any = true;
}
}
if (moved_any) {
backoff.reset();
} else {
backoff.pause();
}
}
}
@@ -119,25 +175,40 @@ 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.
* @param min_pedestal_samples minimum number of pedestal samples to
* accumulate before using the pedestal
*/
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, size_t min_pedestal_samples = 1000)
: 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;
<< ", min_pedestal_samples: " << min_pedestal_samples
<< ", n_threads: " << n_threads
<< ", queue_depth: " << queue_depth;
for (size_t i = 0; i < n_threads; i++) {
m_cluster_finders.push_back(
std::make_unique<
ClusterFinder<ClusterType, FRAME_TYPE, PEDESTAL_TYPE>>(
image_size, nSigma, capacity));
image_size, nSigma, capacity, min_pedestal_samples));
}
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?
@@ -207,13 +278,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);
}
}
@@ -223,11 +303,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++;
}
@@ -240,6 +338,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
@@ -279,6 +390,23 @@ class ClusterFinderMT {
// auto rc = m_input_queue.write(std::move(frame));
// fmt::print("pushed frame {}\n", rc);
// }
/**
* @brief Set the nSigma value for all the cluster finders.
* @param nSigma number of sigma above the pedestal to consider a photon
* during cluster finding.
*/
void set_nSigma(const PEDESTAL_TYPE nSigma) {
// Wait for all queues to be empty before changing the sigma
for (auto &q : m_input_queues) {
while (!q->isEmpty()) {
std::this_thread::sleep_for(m_default_wait);
}
}
for (auto &cf : m_cluster_finders) {
cf->set_nSigma(nSigma);
}
}
};
} // namespace aare
include/aare/ClusterVector.hpp
+99 -47
View File
@@ -1,3 +1,4 @@
// SPDX-License-Identifier: MPL-2.0
#pragma once
#include "aare/Cluster.hpp" //TODO maybe store in seperate file !!!
#include <algorithm>
@@ -19,16 +20,21 @@ template <typename ClusterType,
class ClusterVector; // Forward declaration
/**
* @brief ClusterVector is a container for clusters of various sizes. It
* uses a contiguous memory buffer to store the clusters. It is templated on
* the data type and the coordinate type of the clusters.
* @note push_back can invalidate pointers to elements in the container
* @brief Move-only container that stores fixed-size clusters contiguously.
*
* The pixel type, cluster dimensions, and coordinate type are determined by
* the Cluster specialization supplied as the template argument.
*
* @note push_back, reserve, and resize can invalidate pointers, references, and
* iterators to elements in the container.
* @warning ClusterVector is currently move only to catch unintended copies,
* but this might change since there are probably use cases where copying is
* needed.
* @tparam T data type of the pixels in the cluster
* @tparam ClusterSizeX cluster size in the x dimension
* @tparam ClusterSizeY cluster size in the y dimension
* @tparam CoordType data type of the x and y coordinates of the cluster
* (normally int16_t)
* (normally uint16_t)
*/
template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
typename CoordType>
@@ -43,35 +49,48 @@ class ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> {
/**
* @brief Construct a new ClusterVector object
* @param capacity initial capacity of the buffer in number of clusters
* @param capacity minimum initial capacity in number of clusters
* @param frame_number frame number of the clusters. Default is 0, which is
* also used to indicate that the clusters come from many frames
*/
ClusterVector(size_t capacity = 1024, uint64_t frame_number = 0)
ClusterVector(size_t capacity = 1024, int32_t frame_number = 0)
: m_frame_number(frame_number) {
m_data.reserve(capacity);
}
// Move constructor
ClusterVector(ClusterVector &&other) noexcept
: m_data(other.m_data), m_frame_number(other.m_frame_number) {
other.m_data.clear();
}
ClusterVector(ClusterVector &&other) noexcept = default;
ClusterVector &operator=(ClusterVector &&other) noexcept = default;
// 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;
/**
* @brief Return a filtered copy selected by a one-dimensional Boolean mask.
* @param mask boolean 1d mask true if element in ClusterVector will be
* included
* @return ClusterVector containing the selected clusters in their original
* order and with the original frame number
* @throws std::runtime_error if the mask length differs from size()
*/
ClusterVector operator()(NDView<bool, 1> mask) {
if (static_cast<size_t>(mask.size()) != m_data.size()) {
throw std::runtime_error(
LOCATION + "Mask size does not match number of clusters");
}
return *this;
if (m_data.empty()) {
return ClusterVector(0, frame_number());
}
const auto selected =
static_cast<size_t>(std::count(mask.begin(), mask.end(), true));
ClusterVector result(selected, frame_number());
for (size_t i = 0; i < m_data.size(); ++i) {
if (mask(i)) {
result.push_back(m_data[i]);
}
}
return result;
}
/**
* @brief Sum the pixels in each cluster
* @return std::vector<T> vector of sums for each cluster
* @brief Sum the pixels in each cluster.
* @return One sum for every cluster, in container order
*/
std::vector<T> sum() {
std::vector<T> sums(m_data.size());
@@ -84,17 +103,17 @@ class ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> {
}
/**
* @brief Sum the pixels in the 2x2 subcluster with the biggest pixel sum in
* each cluster
* @return std::vector<T> vector of sums for each cluster
* @brief Find the highest-sum center-adjacent 2x2 subcluster in each
* cluster.
* @return One sum and corner-index pair for every cluster, in container
* order
*/
std::vector<T> sum_2x2() {
std::vector<T> sums_2x2(m_data.size());
std::vector<Sum_index_pair<T, corner>> sum_2x2() {
std::vector<Sum_index_pair<T, corner>> sums_2x2(m_data.size());
std::transform(m_data.begin(), m_data.end(), sums_2x2.begin(),
[](const ClusterType &cluster) {
return cluster.max_sum_2x2().first;
});
std::transform(
m_data.begin(), m_data.end(), sums_2x2.begin(),
[](const ClusterType &cluster) { return cluster.max_sum_2x2(); });
return sums_2x2;
}
@@ -107,10 +126,29 @@ class ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> {
*/
void reserve(size_t capacity) { m_data.reserve(capacity); }
/**
* @brief Change the number of stored clusters.
* @param size new number of clusters
* @note Growing the vector value-initializes new clusters and can
* invalidate pointers, references, and iterators.
*/
void resize(size_t size) { m_data.resize(size); }
/**
* @brief Append a cluster to the vector.
* @param cluster cluster to append
* @note Reallocation invalidates pointers, references, iterators, and
* zero-copy NumPy views of the storage.
*/
void push_back(const ClusterType &cluster) { m_data.push_back(cluster); }
/**
* @brief Append all clusters from another vector.
* @param other vector whose clusters are appended
* @return Reference to this vector
* @note The frame number of this vector is unchanged.
* @warning other must not refer to this vector.
*/
ClusterVector &operator+=(const ClusterVector &other) {
m_data.insert(m_data.end(), other.begin(), other.end());
@@ -127,8 +165,10 @@ class ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> {
*/
bool empty() const { return m_data.empty(); }
/** @brief Return the cluster size in the x dimension. */
uint8_t cluster_size_x() const { return ClusterSizeX; }
/** @brief Return the cluster size in the y dimension. */
uint8_t cluster_size_y() const { return ClusterSizeY; }
/**
@@ -143,7 +183,7 @@ class ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> {
auto end() const { return m_data.end(); }
/**
* @brief Return the size in bytes of a single cluster
* @brief Return the size in bytes of one stored cluster, including padding.
*/
size_t item_size() const {
return sizeof(ClusterType); // 2 * sizeof(CoordType) + ClusterSizeX *
@@ -154,8 +194,8 @@ class ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> {
ClusterType const *data() const { return m_data.data(); }
/**
* @brief Return a reference to the i-th cluster casted to type V
* @tparam V type of the cluster
* @brief Return a reference to the i-th cluster without bounds checking.
* @param i zero-based cluster index
*/
ClusterType &operator[](size_t i) { return m_data[i]; }
@@ -167,23 +207,31 @@ class ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> {
*/
int32_t frame_number() const { return m_frame_number; }
/**
* @brief Set the signed 32-bit frame number associated with the clusters.
* @param frame_number frame number, or 0 for clusters from multiple frames
*/
void set_frame_number(int32_t frame_number) {
m_frame_number = frame_number;
}
};
/**
* @brief Reduce a cluster to a 2x2 cluster by selecting the 2x2 block with the
* highest sum.
* @param cv Clustervector containing clusters to reduce
* @return Clustervector with reduced clusters
* @brief Reduce every cluster to its highest-sum center-adjacent 2x2 block.
* @param cv ClusterVector containing clusters to reduce
* @return ClusterVector of 2x2 clusters in the original order and with the
* original frame number
* @note The cluster is filled using row major ordering starting at the top-left
* (thus for a max subcluster in the top left cornern the photon hit is at
* the fourth position)
*/
template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
typename CoordType = uint16_t>
typename CoordType>
ClusterVector<Cluster<T, 2, 2, CoordType>> reduce_to_2x2(
const ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>>
&cv) {
ClusterVector<Cluster<T, 2, 2, CoordType>> result;
ClusterVector<Cluster<T, 2, 2, CoordType>> result(cv.size(),
cv.frame_number());
for (const auto &c : cv) {
result.push_back(reduce_to_2x2(c));
}
@@ -191,21 +239,25 @@ ClusterVector<Cluster<T, 2, 2, CoordType>> reduce_to_2x2(
}
/**
* @brief Reduce a cluster to a 3x3 cluster by selecting the 3x3 block with the
* highest sum.
* @param cv Clustervector containing clusters to reduce
* @return Clustervector with reduced clusters
* @brief Reduce every cluster to the 3x3 block around its center index.
* @param cv ClusterVector containing clusters to reduce
* @pre ClusterSizeX and ClusterSizeY must both be at least 3, and at least one
* must be greater than 3.
* @return ClusterVector of 3x3 clusters in the original order and with the
* original frame number
* @note Coordinates are preserved.
*/
template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
typename CoordType = uint16_t>
typename CoordType>
ClusterVector<Cluster<T, 3, 3, CoordType>> reduce_to_3x3(
const ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>>
&cv) {
ClusterVector<Cluster<T, 3, 3, CoordType>> result;
ClusterVector<Cluster<T, 3, 3, CoordType>> result(cv.size(),
cv.frame_number());
for (const auto &c : cv) {
result.push_back(reduce_to_3x3(c));
}
return result;
}
} // namespace aare
} // namespace aare
include/aare/CtbRawFile.hpp
+2
View File
@@ -1,3 +1,4 @@
// SPDX-License-Identifier: MPL-2.0
#pragma once
#include "aare/FileInterface.hpp"
@@ -27,6 +28,7 @@ class CtbRawFile {
size_t image_size_in_bytes() const;
size_t frames_in_file() const;
size_t total_frames() const;
RawMasterFile master() const;
include/aare/DetectorGeometry.hpp
+72 -12
View File
@@ -1,6 +1,9 @@
// SPDX-License-Identifier: MPL-2.0
#pragma once
#include "aare/RawMasterFile.hpp" //ROI refactor away
#include "aare/defs.hpp"
#include "aare/ROI.hpp"
#include "aare/ROIGeometry.hpp"
#include <iostream>
namespace aare {
struct ModuleConfig {
@@ -17,16 +20,69 @@ struct ModuleConfig {
};
/**
* @brief Class to hold the geometry of a module. Where pixel 0 is located and
* the size of the module
* @brief Class to hold the geometry of a module. Where pixel 0 is located
* relative to full Detector/ROI and the size of the module e.g. size of ROI
* occupied by the module
*/
struct ModuleGeometry {
/// @brief start x coordinate of the module in the full detector/ROI
int origin_x{};
/// @brief start y coordinate of the module in the full detector/ROI
int origin_y{};
/// @brief height of the module in pixels
int height{};
/// @brief width of the module in pixels
int width{};
/// @brief module index of the module in the detector
int row_index{};
/// @brief module index of the module in the detector
int col_index{};
/**
* @brief checks if module is in the region of interest
*/
bool module_in_roi(const ROI &roi) {
// module is to the left of the roi
bool to_the_left = origin_x + width - 1 < roi.xmin;
bool to_the_right = origin_x >= roi.xmax;
bool above = origin_y + height - 1 < roi.ymin;
bool below = origin_y >= roi.ymax;
return !(to_the_left || to_the_right || below || above);
}
/**
* @brief Update the module geometry given a region of interest
*/
void update_geometry_with_roi(const ROI &roi) {
if (!(module_in_roi(roi))) {
return; // do nothing
}
int new_width = std::min(origin_x + width, static_cast<int>(roi.xmax)) -
std::max(origin_x, static_cast<int>(roi.xmin));
// ROI starts inside module
if (roi.xmin >= origin_x && roi.xmin < origin_x + width) {
origin_x = 0;
} else {
origin_x -= roi.xmin;
}
int new_height =
std::min(origin_y + height, static_cast<int>(roi.ymax)) -
std::max(origin_y, static_cast<int>(roi.ymin));
if (roi.ymin >= origin_y && roi.ymin < origin_y + height) {
origin_y = 0;
} else {
origin_y -= roi.ymin;
}
height = new_height;
width = new_width;
}
};
/**
@@ -40,6 +96,8 @@ class DetectorGeometry {
const xy udp_interfaces_per_module = xy{1, 1},
const bool quad = false);
DetectorGeometry() = default;
~DetectorGeometry() = default;
/**
@@ -48,34 +106,36 @@ class DetectorGeometry {
* @param roi
* @return DetectorGeometry
*/
void update_geometry_with_roi(ROI roi);
size_t n_modules() const;
size_t n_modules_in_roi() const;
size_t pixels_x() const;
size_t pixels_y() const;
size_t modules_x() const;
size_t modules_y() const;
const std::vector<ssize_t> &get_modules_in_roi() const;
ssize_t get_modules_in_roi(const size_t index) const;
xy udp_interfaces_per_module() const;
const std::vector<ModuleGeometry> &get_module_geometries() const;
const ModuleGeometry &get_module_geometries(const size_t index) const;
ModuleGeometry &get_module_geometries(const size_t index);
private:
size_t m_modules_x{};
size_t m_modules_y{};
size_t m_pixels_x{};
size_t m_pixels_y{};
xy m_udp_interfaces_per_module{};
static constexpr ModuleConfig cfg{0, 0};
std::vector<ModuleGeometry> module_geometries{};
std::vector<ssize_t> modules_in_roi{};
// TODO: maybe remove - should be a member in ROIGeometry - in particular
// only works as no overlap between ROIs allowed? maybe just have another
// additional vector of ModuleGeometry for each ROIGeometry - some
// duplication but nicer design?
std::vector<ModuleGeometry> module_geometries;
};
} // namespace aare
include/aare/Dtype.hpp
+1
View File
@@ -1,3 +1,4 @@
// SPDX-License-Identifier: MPL-2.0
#pragma once
#include <cstdint>
#include <map>
include/aare/FastPedestal.hpp
+386
View File
@@ -0,0 +1,386 @@
// 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 <algorithm>
#include <cstddef>
#include <sys/types.h>
namespace aare {
/**
* @brief Maintain per-pixel mean and population 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. Variance is a private double-precision intermediate.
* @tparam PEDESTAL_TYPE Type returned for the mean 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() const {
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 standard deviation of every
* pixel.
* @throws std::runtime_error if ready() is false.
*/
NDArray<PEDESTAL_TYPE, 2> std() const {
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_unchecked(rc_to_index(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:
double variance_unchecked(ssize_t index) const {
const auto &entry = m_sum[index];
const auto mean = entry.sum * m_inv_samples;
const auto variance = std::fma(-mean, mean, entry.sum2 * m_inv_samples);
// Roundoff in the moments can make a near-zero variance negative.
return std::max(variance, 0.0);
}
/**
* @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
+3 -1
View File
@@ -1,3 +1,4 @@
// SPDX-License-Identifier: MPL-2.0
#pragma once
#include "aare/FileInterface.hpp"
#include <memory>
@@ -63,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
+20 -13
View File
@@ -1,3 +1,4 @@
// SPDX-License-Identifier: MPL-2.0
#pragma once
#include "aare/Dtype.hpp"
#include "aare/Frame.hpp"
@@ -33,15 +34,16 @@ struct FileConfig {
DetectorType detector_type{DetectorType::Unknown};
int max_frames_per_file{};
size_t total_frames{};
std::string to_string() const {
return "{ dtype: " + dtype.to_string() +
", rows: " + std::to_string(rows) +
", cols: " + std::to_string(cols) +
", geometry: " + geometry.to_string() +
", detector_type: " + ToString(detector_type) +
", max_frames_per_file: " + std::to_string(max_frames_per_file) +
", total_frames: " + std::to_string(total_frames) + " }";
}
// std::string to_string() const {
// return "{ dtype: " + dtype.to_string() +
// ", rows: " + std::to_string(rows) +
// ", cols: " + std::to_string(cols) +
// ", geometry: " + geometry.to_string() +
// ", detector_type: " + ToString(detector_type) +
// ", max_frames_per_file: " +
// std::to_string(max_frames_per_file) +
// ", total_frames: " + std::to_string(total_frames) + " }";
// }
};
/**
@@ -143,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;
@@ -166,4 +173,4 @@ class FileInterface {
// size_t current_frame{};
};
} // namespace aare
} // namespace aare
include/aare/FilePtr.hpp
+3 -1
View File
@@ -1,3 +1,4 @@
// SPDX-License-Identifier: MPL-2.0
#pragma once
#include <cstdio>
#include <filesystem>
@@ -17,6 +18,7 @@ class FilePtr {
FilePtr &operator=(const FilePtr &) = delete; // since we handle a resource
FilePtr(FilePtr &&other);
FilePtr &operator=(FilePtr &&other);
explicit operator bool() const noexcept;
FILE *get();
ssize_t tell();
void seek(ssize_t offset, int whence = SEEK_SET) {
@@ -27,4 +29,4 @@ class FilePtr {
~FilePtr();
};
} // namespace aare
} // namespace aare
include/aare/Fit.hpp
+51 -106
View File
@@ -1,120 +1,65 @@
// SPDX-License-Identifier: MPL-2.0
#pragma once
#include <cmath>
#include <fmt/core.h>
#include <vector>
#include "aare/FitModel.hpp"
#include "aare/NDArray.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
// ---------------------------------------------------------------------------
// Minuit2-based pixel fitting.
// Template bodies and explicit instantiations live in src/Fit.cpp.
// ---------------------------------------------------------------------------
/**
* @brief Estimate the initial parameters for a Gaussian fit
* @brief Fit a single pixel's data using Minuit2.
*
* User-precedence rules:
* - Fixed parameters: untouched (value and fixed flag preserved from model).
* - User-set start: value preserved, step size auto-filled.
* - Neither: both value and step size auto-filled from data.
*
* @tparam Model Model struct (Gaussian, RisingScurve, …).
*
* @param model The FitModel configuration (read-only).
* @param x Scan points.
* @param y Measured values for this pixel.
* @param y_err Per-point uncertainties (empty view -> unweighted fit).
*
* @return NDArray<double,1> of size:
* - compute_errors: [p0..pN, err0..errN, chi2] -> 2*npar + 1
* - otherwise: [p0..pN, chi2] -> npar + 1
*/
std::array<double, 3> gaus_init_par(const NDView<double, 1> x,
const NDView<double, 1> y);
template <typename Model>
NDArray<double, 1> fit_pixel(const FitModel<Model> &model, NDView<double, 1> x,
NDView<double, 1> y, NDView<double, 1> y_err);
std::array<double, 2> pol1_init_par(const NDView<double, 1> x,
const NDView<double, 1> y);
std::array<double, 6> scurve_init_par(const NDView<double, 1> x,
const NDView<double, 1> y);
std::array<double, 6> scurve2_init_par(const NDView<double, 1> x,
const NDView<double, 1> y);
static constexpr int DEFAULT_NUM_THREADS = 4;
// Overload: uncertainties not provided
template <typename Model>
NDArray<double, 1> fit_pixel(const FitModel<Model> &model, NDView<double, 1> x,
NDView<double, 1> y);
/**
* @brief Fit a 1D Gaussian to data.
* @param data data to fit
* @param x x values
* @brief Fit all pixels in a 3D data cube (rows x cols x n_scan).
*
* @tparam Model Model struct.
*
* @param model Fit configuration shared by all pixels.
* @param x Scan points, shape `(n_scan)`.
* @param y Measured values, shape `(rows, cols, n_scan)`.
* @param y_err Uncertainties, same shape as y, or empty for unweighted
* fits.
* @param par_out Output parameters, shape `(rows, cols, npar)`.
* @param err_out Output parameter errors, shape `(rows, cols, npar)`, if
* used.
* @param chi2_out Output chi-squared / objective values, shape `(rows,
* cols)`.
* @param n_threads Number of threads used to split rows.
*/
NDArray<double, 1> fit_gaus(NDView<double, 1> x, NDView<double, 1> y);
template <typename Model>
void fit_3d(const FitModel<Model> &model, NDView<double, 1> x,
NDView<double, 3> y, NDView<double, 3> y_err,
NDView<double, 3> par_out, NDView<double, 3> err_out,
NDView<double, 2> chi2_out, int n_threads);
/**
* @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);
} // namespace aare
} // namespace aare
include/aare/FitModel.hpp
+102
View File
@@ -0,0 +1,102 @@
// SPDX-License-Identifier: MPL-2.0
#pragma once
#include "aare/Models.hpp"
#include <array>
#include <memory>
#include <stdexcept>
#include <string>
#include <vector>
namespace aare {
template <typename Model> class FitModel {
// Forward declaration only — full definition lives in src/FitModelImpl.hpp
// and is never installed. Downstream code cannot dereference this pointer.
struct FitModelImpl; // stores Minuit2 user-set parameters
std::unique_ptr<FitModelImpl> impl_;
unsigned int max_calls_;
double tolerance_;
bool compute_errors_;
std::array<bool, Model::npar> user_fixed_{};
std::array<bool, Model::npar> user_start_{};
/** @brief Safely resolve a parameter name to its index. */
unsigned int checked_index(const std::string &name) const;
public:
static constexpr std::size_t npar = Model::npar;
/**
* @brief Construct a fit model with sensible defaults.
*
* @param strategy Minuit2 strategy level (0 = fast/gradient, 1 =
* default).
* @param max_calls Maximum FCN calls per pixel minimisation.
* @param tolerance Minuit2 EDM tolerance.
* @param compute_errors If true, run MnHesse after minimisation.
*/
explicit FitModel(unsigned int strategy = 0, unsigned int max_calls = 100,
double tolerance = 0.5, bool compute_errors = false);
// Destructor must be defined in Fit.cpp where FitModelImpl is complete.
~FitModel();
FitModel(const FitModel &);
FitModel &operator=(const FitModel &);
FitModel(FitModel &&) noexcept = default;
FitModel &operator=(FitModel &&) noexcept = default;
/** @brief Set lower and upper bounds for parameter idx.*/
void SetParLimits(unsigned int idx, double lo, double hi);
/**
* @brief Fix parameter idx at value val.
*
* Excluded from minimisation. Automatic estimates will not touch it.
*/
void FixParameter(unsigned int idx, double val);
/** @brief Release a previously fixed parameter, re-enabling auto estimates.
*/
void ReleaseParameter(unsigned int idx);
void ReleaseParameter(const std::string &name);
/** @brief Set an explicit starting value for parameter idx.*/
void SetParameter(unsigned int idx, double val);
void SetParameter(const std::string &name, double val);
void FixParameter(const std::string &name, double val);
void SetParLimits(const std::string &name, double lo, double hi);
std::string GetParName(unsigned int idx) const;
std::vector<std::string> GetParNames() const;
static constexpr std::size_t GetNpar() noexcept { return npar; }
void SetMaxCalls(unsigned int n) { max_calls_ = n; }
void SetTolerance(double t) { tolerance_ = t; }
void SetComputeErrors(bool b) { compute_errors_ = b; }
unsigned int max_calls() const { return max_calls_; }
double tolerance() const { return tolerance_; }
bool compute_errors() const { return compute_errors_; }
bool is_user_fixed(unsigned int idx) const { return user_fixed_[idx]; }
bool is_user_start(unsigned int idx) const { return user_start_[idx]; }
// Returns the internal Minuit2 state. FitModelImpl is an incomplete type
// here; only callers that include src/FitModelImpl.hpp can dereference it.
FitModelImpl *impl() const { return impl_.get(); }
};
// Suppress implicit instantiation for all supported model types.
// Definitions live in src/Fit.cpp.
extern template class FitModel<model::Gaussian>;
extern template class FitModel<model::GaussianErfcPlateau>;
extern template class FitModel<model::GaussianChargeSharing>;
extern template class FitModel<model::GaussianChargeSharingKb>;
extern template class FitModel<model::Pol1>;
extern template class FitModel<model::Pol2>;
extern template class FitModel<model::RisingScurve>;
extern template class FitModel<model::FallingScurve>;
} // namespace aare
include/aare/Frame.hpp
+1
View File
@@ -1,3 +1,4 @@
// SPDX-License-Identifier: MPL-2.0
#pragma once
#include "aare/Dtype.hpp"
#include "aare/NDArray.hpp"
Loaded 100 of 277 files, more files were not shown because too many files have changed in this diff. Show more