detectors/aare
25
0
Fork 0
mirror of https://github.com/slsdetectorgroup/aare.git synced 2025-06-05 12:30:39 +02:00
Code Issues Actions Packages Releases Activity

merge from latest developer
All checks were successful
Build on RHEL9 / build (push) Successful in 2m21s
Details
Build on RHEL8 / build (push) Successful in 2m40s
Details

Browse Source
This commit is contained in:
maliakal_d 2025-06-03 11:24:37 +02:00
commit 5681e18403
135 changed files with 9254 additions and 1500 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

42
.clang-tidy Normal file
View File

@ -0,0 +1,42 @@
---
Checks: '*,
-altera-*,
-android-cloexec-fopen,
-cppcoreguidelines-pro-bounds-array-to-pointer-decay,
-cppcoreguidelines-pro-bounds-pointer-arithmetic,
-fuchsia*,
-readability-else-after-return,
-readability-avoid-const-params-in-decls,
-readability-identifier-length,
-cppcoreguidelines-pro-bounds-constant-array-index,
-cppcoreguidelines-pro-type-reinterpret-cast,
-llvm-header-guard,
-modernize-use-nodiscard,
-misc-non-private-member-variables-in-classes,
-readability-static-accessed-through-instance,
-readability-braces-around-statements,
-readability-isolate-declaration,
-readability-implicit-bool-conversion,
-readability-identifier-length,
-readability-identifier-naming,
-hicpp-signed-bitwise,
-hicpp-no-array-decay,
-hicpp-braces-around-statements,
-google-runtime-references,
-google-readability-todo,
-google-readability-braces-around-statements,
-modernize-use-trailing-return-type,
-llvmlibc-*'
HeaderFilterRegex: \.hpp
FormatStyle: none
CheckOptions:
- { key: readability-identifier-naming.NamespaceCase, value: lower_case }
# - { key: readability-identifier-naming.FunctionCase, value: lower_case }
- { key: readability-identifier-naming.ClassCase, value: CamelCase }
# - { key: readability-identifier-naming.MethodCase, value: CamelCase }
# - { key: readability-identifier-naming.StructCase, value: CamelCase }
# - { key: readability-identifier-naming.VariableCase, value: lower_case }
- { key: readability-identifier-naming.GlobalConstantCase, value: UPPER_CASE }
...

58
.gitea/workflows/cmake_build.yml Normal file
View File

@ -0,0 +1,58 @@
name: Build the package using cmake then documentation
on:
workflow_dispatch:
permissions:
contents: read
pages: write
id-token: write
jobs:
build:
strategy:
fail-fast: false
matrix:
platform: [ubuntu-latest, ]
python-version: ["3.12", ]
runs-on: ${{ matrix.platform }}
defaults:
run:
shell: "bash -l {0}"
steps:
- uses: actions/checkout@v4
- name: Setup dev env
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"
- name: Build library
run: |
mkdir build
cd build
cmake .. -DAARE_SYSTEM_LIBRARIES=ON -DAARE_DOCS=ON
make -j 2
make docs

36
.gitea/workflows/rh8-native.yml Normal file
View File

@ -0,0 +1,36 @@
name: Build on RHEL8
on:
push:
workflow_dispatch:
permissions:
contents: read
jobs:
build:
runs-on: "ubuntu-latest"
container:
image: gitea.psi.ch/images/rhel8-developer-gitea-actions
steps:
# workaround until actions/checkout@v4 is available for RH8
# - uses: actions/checkout@v4
- name: Clone repository
run: |
echo Cloning ${{ github.ref_name }}
git clone https://${{secrets.GITHUB_TOKEN}}@gitea.psi.ch/${{ github.repository }}.git --branch=${{ github.ref_name }} .
- name: Install dependencies
run: |
dnf install -y cmake python3.12 python3.12-devel python3.12-pip
- name: Build library
run: |
mkdir build && cd build
cmake .. -DAARE_PYTHON_BINDINGS=ON -DAARE_TESTS=ON -DPython_FIND_VIRTUALENV=FIRST
make -j 2
- name: C++ unit tests
working-directory: ${{gitea.workspace}}/build
run: ctest

31
.gitea/workflows/rh9-native.yml Normal file
View File

@ -0,0 +1,31 @@
name: Build on RHEL9
on:
push:
workflow_dispatch:
permissions:
contents: read
jobs:
build:
runs-on: "ubuntu-latest"
container:
image: gitea.psi.ch/images/rhel9-developer-gitea-actions
steps:
- uses: actions/checkout@v4
- name: Install dependencies
run: |
dnf install -y cmake python3.12 python3.12-devel python3.12-pip
- name: Build library
run: |
mkdir build && cd build
cmake .. -DAARE_PYTHON_BINDINGS=ON -DAARE_TESTS=ON
make -j 2
- name: C++ unit tests
working-directory: ${{gitea.workspace}}/build
run: ctest

16
.github/workflows/deploy.yml → .github/workflows/build_and_deploy_conda.yml vendored
View File

@ -1,10 +1,9 @@
name: Build pkgs and deploy if on main
on:
push:
branches:
- main
- developer
release:
types:
- published
jobs:
build:
@ -25,16 +24,15 @@ jobs:
- uses: actions/checkout@v4
- name: Get conda
uses: conda-incubator/setup-miniconda@v3.0.4
uses: conda-incubator/setup-miniconda@v3
with:
python-version: ${{ matrix.python-version }}
environment-file: etc/dev-env.yml
miniforge-version: latest
channels: conda-forge
- name: Prepare
run: conda install conda-build=24.9 conda-verify pytest anaconda-client
conda-remove-defaults: "true"
- name: Enable upload
if: github.ref == 'refs/heads/main'
run: conda config --set anaconda_upload yes
- name: Build

41
.github/workflows/build_conda.yml vendored Normal file
View File

@ -0,0 +1,41 @@
name: Build pkgs and deploy if on main
on:
push:
branches:
- developer
jobs:
build:
strategy:
fail-fast: false
matrix:
platform: [ubuntu-latest, ] # macos-12, windows-2019]
python-version: ["3.12",]
runs-on: ${{ matrix.platform }}
# The setup-miniconda action needs this to activate miniconda
defaults:
run:
shell: "bash -l {0}"
steps:
- uses: actions/checkout@v4
- 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"
- name: Disable upload
run: conda config --set anaconda_upload no
- name: Build
run: conda build conda-recipe

12
.github/workflows/build_docs.yml vendored
View File

@ -5,7 +5,6 @@ on:
push:
permissions:
contents: read
pages: write
@ -16,12 +15,11 @@ jobs:
strategy:
fail-fast: false
matrix:
platform: [ubuntu-latest, ] # macos-12, windows-2019]
platform: [ubuntu-latest, ]
python-version: ["3.12",]
runs-on: ${{ matrix.platform }}
# The setup-miniconda action needs this to activate miniconda
defaults:
run:
shell: "bash -l {0}"
@ -30,13 +28,13 @@ jobs:
- uses: actions/checkout@v4
- name: Get conda
uses: conda-incubator/setup-miniconda@v3.0.4
uses: conda-incubator/setup-miniconda@v3
with:
python-version: ${{ matrix.python-version }}
environment-file: etc/dev-env.yml
miniforge-version: latest
channels: conda-forge
- name: Prepare
run: conda install doxygen sphinx=7.1.2 breathe pybind11 sphinx_rtd_theme furo nlohmann_json zeromq fmt numpy
conda-remove-defaults: "true"
- name: Build library
run: |

64
.github/workflows/build_wheel.yml vendored Normal file
View File

@ -0,0 +1,64 @@
name: Build wheel
on:
workflow_dispatch:
pull_request:
push:
branches:
- main
release:
types:
- published
jobs:
build_wheels:
name: Build wheels on ${{ matrix.os }}
runs-on: ${{ matrix.os }}
strategy:
matrix:
os: [ubuntu-latest,]
steps:
- uses: actions/checkout@v4
- name: Build wheels
run: pipx run cibuildwheel==2.23.0
- uses: actions/upload-artifact@v4
with:
name: cibw-wheels-${{ matrix.os }}-${{ strategy.job-index }}
path: ./wheelhouse/*.whl
build_sdist:
name: Build source distribution
runs-on: ubuntu-latest
steps:
- uses: actions/checkout@v4
- name: Build sdist
run: pipx run build --sdist
- uses: actions/upload-artifact@v4
with:
name: cibw-sdist
path: dist/*.tar.gz
upload_pypi:
needs: [build_wheels, build_sdist]
runs-on: ubuntu-latest
environment: pypi
permissions:
id-token: write
if: github.event_name == 'release' && github.event.action == 'published'
# or, alternatively, upload to PyPI on every tag starting with 'v' (remove on: release above to use this)
# if: github.event_name == 'push' && startsWith(github.ref, 'refs/tags/v')
steps:
- uses: actions/download-artifact@v4
with:
# unpacks all CIBW artifacts into dist/
pattern: cibw-*
path: dist
merge-multiple: true
- uses: pypa/gh-action-pypi-publish@release/v1

3
.gitignore vendored
View File

@ -17,7 +17,8 @@ Testing/
ctbDict.cpp
ctbDict.h
wheelhouse/
dist/
*.pyc
*/__pycache__/*

149
CMakeLists.txt
View File

@ -1,16 +1,29 @@
cmake_minimum_required(VERSION 3.14)
cmake_minimum_required(VERSION 3.15)
project(aare
VERSION 1.0.0
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")
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_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
)
message(STATUS "Building from git hash: ${GIT_HASH}")
if (${CMAKE_VERSION} VERSION_GREATER "3.24")
cmake_policy(SET CMP0135 NEW) #Fetch content download timestamp
endif()
@ -31,7 +44,7 @@ set(CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/cmake" ${CMAKE_MODULE_PATH})
# General options
option(AARE_PYTHON_BINDINGS "Build python bindings" ON)
option(AARE_PYTHON_BINDINGS "Build python bindings" OFF)
option(AARE_TESTS "Build tests" OFF)
option(AARE_BENCHMARKS "Build benchmarks" OFF)
option(AARE_EXAMPLES "Build examples" OFF)
@ -41,7 +54,7 @@ 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_HDF5 "Hdf5 File Format" OFF)
option(AARE_ASAN "Enable AddressSanitizer" OFF)
# Configure which of the dependencies to use FetchContent for
option(AARE_FETCH_FMT "Use FetchContent to download fmt" ON)
@ -49,6 +62,7 @@ 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)
#Convenience option to use system libraries only (no FetchContent)
@ -60,10 +74,15 @@ if(AARE_SYSTEM_LIBRARIES)
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
endif()
if(AARE_VERBOSE)
add_compile_definitions(AARE_VERBOSE)
add_compile_definitions(AARE_LOG_LEVEL=aare::logDEBUG5)
else()
add_compile_definitions(AARE_LOG_LEVEL=aare::logERROR)
endif()
if(AARE_CUSTOM_ASSERT)
@ -75,18 +94,70 @@ if(AARE_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)
# For cmake < 3.28 we can't supply EXCLUDE_FROM_ALL to FetchContent_Declare
# so we need this workaround
if (${CMAKE_VERSION} VERSION_LESS "3.28")
FetchContent_Declare(
lmfit
GIT_REPOSITORY https://jugit.fz-juelich.de/mlz/lmfit.git
GIT_TAG main
PATCH_COMMAND ${LMFIT_PATCH_COMMAND}
UPDATE_DISCONNECTED 1
)
else()
FetchContent_Declare(
lmfit
GIT_REPOSITORY https://jugit.fz-juelich.de/mlz/lmfit.git
GIT_TAG main
PATCH_COMMAND ${LMFIT_PATCH_COMMAND}
UPDATE_DISCONNECTED 1
EXCLUDE_FROM_ALL 1
)
endif()
#Disable what we don't need from lmfit
set(BUILD_TESTING OFF CACHE BOOL "")
set(LMFIT_CPPTEST OFF CACHE BOOL "")
set(LIB_MAN OFF CACHE BOOL "")
set(LMFIT_CPPTEST OFF CACHE BOOL "")
set(BUILD_SHARED_LIBS OFF CACHE BOOL "")
if (${CMAKE_VERSION} VERSION_LESS "3.28")
if(NOT lmfit_POPULATED)
FetchContent_Populate(lmfit)
add_subdirectory(${lmfit_SOURCE_DIR} ${lmfit_BINARY_DIR} EXCLUDE_FROM_ALL)
endif()
else()
FetchContent_MakeAvailable(lmfit)
endif()
set_property(TARGET lmfit PROPERTY POSITION_INDEPENDENT_CODE ON)
else()
find_package(lmfit REQUIRED)
endif()
if(AARE_FETCH_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")
@ -128,8 +199,8 @@ if (AARE_FETCH_FMT)
LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR}
ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR}
RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR}
INCLUDES DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}
)
INCLUDES DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}
)
else()
find_package(fmt 6 REQUIRED)
endif()
@ -147,7 +218,6 @@ if (AARE_FETCH_JSON)
install(
TARGETS nlohmann_json
EXPORT "${TARGETS_EXPORT_NAME}"
)
message(STATUS "target: ${NLOHMANN_JSON_TARGET_NAME}")
else()
@ -226,13 +296,6 @@ if(CMAKE_BUILD_TYPE STREQUAL "Release")
target_compile_options(aare_compiler_flags INTERFACE -O3)
else()
message(STATUS "Debug build")
target_compile_options(
aare_compiler_flags
INTERFACE
-Og
-ggdb3
)
endif()
# Common flags for GCC and Clang
@ -257,7 +320,21 @@ target_compile_options(
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
)
endif()
@ -273,14 +350,23 @@ endif()
set(PUBLICHEADERS
include/aare/ArrayExpr.hpp
include/aare/CalculateEta.hpp
include/aare/Cluster.hpp
include/aare/ClusterFinder.hpp
include/aare/ClusterFile.hpp
include/aare/CtbRawFile.hpp
include/aare/ClusterVector.hpp
include/aare/decode.hpp
include/aare/defs.hpp
include/aare/Dtype.hpp
include/aare/File.hpp
include/aare/Fit.hpp
include/aare/FileInterface.hpp
include/aare/FilePtr.hpp
include/aare/Frame.hpp
include/aare/GainMap.hpp
include/aare/geo_helpers.hpp
include/aare/JungfrauDataFile.hpp
include/aare/NDArray.hpp
include/aare/NDView.hpp
include/aare/NumpyFile.hpp
@ -291,23 +377,30 @@ set(PUBLICHEADERS
include/aare/RawMasterFile.hpp
include/aare/RawSubFile.hpp
include/aare/VarClusterFinder.hpp
include/aare/utils/task.hpp
)
set(SourceFiles
${CMAKE_CURRENT_SOURCE_DIR}/src/CtbRawFile.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/ClusterFile.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/defs.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/Dtype.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/decode.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/Frame.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/File.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/FilePtr.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/Fit.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/geo_helpers.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/JungfrauDataFile.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/NumpyFile.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/NumpyHelpers.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/Interpolator.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/PixelMap.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/RawFile.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/RawSubFile.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/RawMasterFile.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/utils/task.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/utils/ifstream_helpers.cpp
)
# HDF5
@ -329,9 +422,12 @@ endif (AARE_HDF5)
add_library(aare_core STATIC ${SourceFiles})
target_include_directories(aare_core PUBLIC
"$<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/include>"
"$<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>"
"$<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>"
)
set(THREADS_PREFER_PTHREAD_FLAG ON)
find_package(Threads REQUIRED)
target_link_libraries(
aare_core
PUBLIC
@ -340,6 +436,9 @@ target_link_libraries(
${STD_FS_LIB} # from helpers.cmake
PRIVATE
aare_compiler_flags
Threads::Threads
$<BUILD_INTERFACE:lmfit>
)
if (AARE_HDF5 AND HDF5_FOUND)
@ -363,17 +462,29 @@ endif()
if(AARE_TESTS)
set(TestSources
${CMAKE_CURRENT_SOURCE_DIR}/src/algorithm.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/geo_helpers.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
)
if(HDF5_FOUND)
list (APPEND TestSources
@ -457,4 +568,4 @@ 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)
endif()
endif()

1
VERSION Normal file
View File

@ -0,0 +1 @@
2025.5.22

30
benchmarks/CMakeLists.txt
View File

@ -1,11 +1,27 @@
find_package(benchmark REQUIRED)
add_executable(ndarray_benchmark ndarray_benchmark.cpp)
include(FetchContent)
target_link_libraries(ndarray_benchmark benchmark::benchmark aare_core aare_compiler_flags)
# target_link_libraries(tests PRIVATE aare_core aare_compiler_flags)
set_target_properties(ndarray_benchmark PROPERTIES
RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}
# OUTPUT_NAME run_tests
FetchContent_Declare(
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)
FetchContent_MakeAvailable(benchmark)
add_executable(benchmarks)
target_sources(benchmarks PRIVATE ndarray_benchmark.cpp calculateeta_benchmark.cpp)
# Link Google Benchmark and other necessary libraries
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
)

70
benchmarks/calculateeta_benchmark.cpp Normal file
View File

@ -0,0 +1,70 @@
#include "aare/CalculateEta.hpp"
#include "aare/ClusterFile.hpp"
#include <benchmark/benchmark.h>
using namespace aare;
class ClusterFixture : public benchmark::Fixture {
public:
Cluster<int, 2, 2> cluster_2x2{};
Cluster<int, 3, 3> cluster_3x3{};
private:
using benchmark::Fixture::SetUp;
void SetUp([[maybe_unused]] const benchmark::State &state) override {
int temp_data[4] = {1, 2, 3, 1};
std::copy(std::begin(temp_data), std::end(temp_data),
std::begin(cluster_2x2.data));
cluster_2x2.x = 0;
cluster_2x2.y = 0;
int temp_data2[9] = {1, 2, 3, 1, 3, 4, 5, 1, 20};
std::copy(std::begin(temp_data2), std::end(temp_data2),
std::begin(cluster_3x3.data));
cluster_3x3.x = 0;
cluster_3x3.y = 0;
}
// void TearDown(::benchmark::State& state) {
// }
};
BENCHMARK_F(ClusterFixture, Calculate2x2Eta)(benchmark::State &st) {
for (auto _ : st) {
// This code gets timed
Eta2 eta = calculate_eta2(cluster_2x2);
benchmark::DoNotOptimize(eta);
}
}
// almost takes double the time
BENCHMARK_F(ClusterFixture,
CalculateGeneralEtaFor2x2Cluster)(benchmark::State &st) {
for (auto _ : st) {
// This code gets timed
Eta2 eta = calculate_eta2<int, 2, 2>(cluster_2x2);
benchmark::DoNotOptimize(eta);
}
}
BENCHMARK_F(ClusterFixture, Calculate3x3Eta)(benchmark::State &st) {
for (auto _ : st) {
// This code gets timed
Eta2 eta = calculate_eta2(cluster_3x3);
benchmark::DoNotOptimize(eta);
}
}
// almost takes double the time
BENCHMARK_F(ClusterFixture,
CalculateGeneralEtaFor3x3Cluster)(benchmark::State &st) {
for (auto _ : st) {
// This code gets timed
Eta2 eta = calculate_eta2<int, 3, 3>(cluster_3x3);
benchmark::DoNotOptimize(eta);
}
}
// BENCHMARK_MAIN();

23
conda-recipe/conda_build_config.yaml
View File

@ -1,28 +1,5 @@
python:
- 3.11
- 3.11
- 3.11
- 3.12
- 3.12
- 3.12
- 3.13
numpy:
- 1.26
- 2.0
- 2.1
- 1.26
- 2.0
- 2.1
- 2.1
zip_keys:
- python
- numpy
pin_run_as_build:
numpy: x.x
python: x.x

48
conda-recipe/meta.yaml
View File

@ -1,7 +1,10 @@
source:
path: ../
{% set version = load_file_regex(load_file = 'VERSION', regex_pattern = '(\d+(?:\.\d+)*(?:[\+\w\.]+))').group(1) %}
package:
name: aare
version: 2024.11.28.dev0 #TODO! how to not duplicate this?
version: {{version}}
source:
path: ..
@ -9,44 +12,39 @@ source:
build:
number: 0
script:
- unset CMAKE_GENERATOR && {{ PYTHON }} -m pip install . -vv # [not win]
- {{ PYTHON }} -m pip install . -vv # [win]
- unset CMAKE_GENERATOR && {{ PYTHON }} -m pip install . -vv
requirements:
build:
- python {{python}}
- numpy {{ numpy }}
- {{ compiler('cxx') }}
host:
- cmake
- ninja
- python {{python}}
- numpy {{ numpy }}
host:
- python
- pip
- numpy=2.1
- scikit-build-core
- pybind11 >=2.13.0
- fmt
- zeromq
- nlohmann_json
- catch2
- matplotlib # needed in host to solve the environment for run
run:
- python {{python}}
- numpy {{ numpy }}
- python
- {{ pin_compatible('numpy') }}
- matplotlib
test:
imports:
- aare
# requires:
# - pytest
# source_files:
# - tests
# commands:
# - pytest tests
requires:
- pytest
- boost-histogram
source_files:
- python/tests
commands:
- python -m pytest python/tests
about:
summary: An example project built with pybind11 and scikit-build.
# license_file: LICENSE
summary: Data analysis library for hybrid pixel detectors from PSI

27
docs/CMakeLists.txt
View File

@ -12,33 +12,6 @@ set(SPHINX_BUILD ${CMAKE_CURRENT_BINARY_DIR})
file(GLOB SPHINX_SOURCE_FILES CONFIGURE_DEPENDS "src/*.rst")
# set(SPHINX_SOURCE_FILES
# src/index.rst
# src/Installation.rst
# src/Requirements.rst
# src/NDArray.rst
# src/NDView.rst
# src/File.rst
# src/Frame.rst
# src/Dtype.rst
# src/ClusterFinder.rst
# src/ClusterFile.rst
# src/Pedestal.rst
# src/RawFile.rst
# src/RawSubFile.rst
# src/RawMasterFile.rst
# src/Hdf5File.rst
# src/Hdf5SubFile.rst
# src/Hdf5MasterFile.rst
# src/VarClusterFinder.rst
# src/pyVarClusterFinder.rst
# src/pyFile.rst
# src/pyCtbRawFile.rst
# src/pyRawFile.rst
# src/pyRawMasterFile.rst
# src/pyHdf5File.rst
# src/pyHdf5MasterFile.rst
# )
foreach(filename ${SPHINX_SOURCE_FILES})

1
docs/conf.py.in
View File

@ -29,7 +29,6 @@ version = '@PROJECT_VERSION@'
# extensions coming with Sphinx (named 'sphinx.ext.*') or your custom
# ones.
extensions = ['breathe',
'sphinx_rtd_theme',
'sphinx.ext.autodoc',
'sphinx.ext.napoleon',
]

7
docs/src/ClusterFinderMT.rst Normal file
View File

@ -0,0 +1,7 @@
ClusterFinderMT
==================
.. doxygenclass:: aare::ClusterFinderMT
:members:
:undoc-members:

6
docs/src/ClusterVector.rst Normal file
View File

@ -0,0 +1,6 @@
ClusterVector
=============
.. doxygenclass:: aare::ClusterVector
:members:
:undoc-members:

25
docs/src/JungfrauDataFile.rst Normal file
View File

@ -0,0 +1,25 @@
JungfrauDataFile
==================
JungfrauDataFile is a class to read the .dat files that are produced by Aldo's receiver.
It is mostly used for calibration.
The structure of the file is:
* JungfrauDataHeader
* Binary data (256x256, 256x1024 or 512x1024)
* JungfrauDataHeader
* ...
There is no metadata indicating number of frames or the size of the image, but this
will be infered by this reader.
.. doxygenstruct:: aare::JungfrauDataHeader
:members:
:undoc-members:
:private-members:
.. doxygenclass:: aare::JungfrauDataFile
:members:
:undoc-members:
:private-members:

47
docs/src/Tests.rst Normal file
View File

@ -0,0 +1,47 @@
****************
Tests
****************
We test the code both from the C++ and Python API. By default only tests that does not require image data is run.
C++
~~~~~~~~~~~~~~~~~~
.. code-block:: bash
mkdir build
cd build
cmake .. -DAARE_TESTS=ON
make -j 4
export AARE_TEST_DATA=/path/to/test/data
./run_test [.files] #or using ctest, [.files] is the option to include tests needing data
Python
~~~~~~~~~~~~~~~~~~
.. code-block:: bash
#From the root dir of the library
python -m pytest python/tests --files # passing --files will run the tests needing data
Getting the test data
~~~~~~~~~~~~~~~~~~~~~~~~
.. attention ::
The tests needing the test data are not run by default. To make the data available, you need to set the environment variable
AARE_TEST_DATA to the path of the test data directory. Then pass either [.files] for the C++ tests or --files for Python
The image files needed for the test are large and are not included in the repository. They are stored
using GIT LFS in a separate repository. To get the test data, you need to clone the repository.
To do this, you need to have GIT LFS installed. You can find instructions on how to install it here: https://git-lfs.github.com/
Once you have GIT LFS installed, you can clone the repository like any normal repo using:
.. code-block:: bash
git clone https://gitea.psi.ch/detectors/aare-test-data.git

5
docs/src/algorithm.rst Normal file
View File

@ -0,0 +1,5 @@
algorithm
=============
.. doxygenfile:: algorithm.hpp

17
docs/src/index.rst
View File

@ -20,9 +20,6 @@ AARE
Requirements
Consume
.. toctree::
:caption: Python API
:maxdepth: 1
@ -30,24 +27,32 @@ AARE
pyFile
pyCtbRawFile
pyClusterFile
pyClusterVector
pyJungfrauDataFile
pyRawFile
pyRawMasterFile
pyHdf5File
pyHdf5MasterFile
pyVarClusterFinder
pyFit
.. toctree::
:caption: C++ API
:maxdepth: 1
algorithm
NDArray
NDView
Frame
File
Dtype
ClusterFinder
ClusterFinderMT
ClusterFile
ClusterVector
JungfrauDataFile
Pedestal
RawFile
RawSubFile
@ -58,4 +63,8 @@ AARE
.. toctree::
:caption: Developer
:maxdepth: 3
Tests

33
docs/src/pyClusterVector.rst Normal file
View File

@ -0,0 +1,33 @@
ClusterVector
================
The ClusterVector, holds clusters from the ClusterFinder. Since it is templated
in C++ we use a suffix indicating the data type in python. The suffix is
``_i`` for integer, ``_f`` for float, and ``_d`` for double.
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:: ClusterVector_i
:members:
:undoc-members:
:show-inheritance:
:inherited-members:

19
docs/src/pyFit.rst Normal file
View File

@ -0,0 +1,19 @@
Fit
========
.. py:currentmodule:: aare
**Functions**
.. autofunction:: gaus
.. autofunction:: pol1
**Fitting**
.. autofunction:: fit_gaus
.. autofunction:: fit_pol1

10
docs/src/pyJungfrauDataFile.rst Normal file
View File

@ -0,0 +1,10 @@
JungfrauDataFile
===================
.. py:currentmodule:: aare
.. autoclass:: JungfrauDataFile
:members:
:undoc-members:
:show-inheritance:
:inherited-members:

13
etc/dev-env.yml Normal file
View File

@ -0,0 +1,13 @@
name: dev-environment
channels:
- conda-forge
dependencies:
- anaconda-client
- conda-build
- doxygen
- sphinx=7.1.2
- breathe
- sphinx_rtd_theme
- furo
- zeromq

36
include/aare/ArrayExpr.hpp
View File

@ -1,22 +1,24 @@
#pragma once
#include <cstdint> //int64_t
#include <cstddef> //size_t
#include <cstdint>
#include <cstddef>
#include <array>
#include <cassert>
#include "aare/defs.hpp"
namespace aare {
template <typename E, int64_t Ndim> class ArrayExpr {
template <typename E, ssize_t Ndim> class ArrayExpr {
public:
static constexpr bool is_leaf = false;
auto operator[](size_t i) const { return static_cast<E const &>(*this)[i]; }
auto operator()(size_t i) const { return static_cast<E const &>(*this)[i]; }
auto size() const { return static_cast<E const &>(*this).size(); }
std::array<int64_t, Ndim> shape() const { return static_cast<E const &>(*this).shape(); }
std::array<ssize_t, Ndim> shape() const { return static_cast<E const &>(*this).shape(); }
};
template <typename A, typename B, int64_t Ndim>
template <typename A, typename B, ssize_t Ndim>
class ArrayAdd : public ArrayExpr<ArrayAdd<A, B, Ndim>, Ndim> {
const A &arr1_;
const B &arr2_;
@ -27,10 +29,10 @@ class ArrayAdd : public ArrayExpr<ArrayAdd<A, B, Ndim>, Ndim> {
}
auto operator[](int i) const { return arr1_[i] + arr2_[i]; }
size_t size() const { return arr1_.size(); }
std::array<int64_t, Ndim> shape() const { return arr1_.shape(); }
std::array<ssize_t, Ndim> shape() const { return arr1_.shape(); }
};
template <typename A, typename B, int64_t Ndim>
template <typename A, typename B, ssize_t Ndim>
class ArraySub : public ArrayExpr<ArraySub<A, B, Ndim>, Ndim> {
const A &arr1_;
const B &arr2_;
@ -41,10 +43,10 @@ class ArraySub : public ArrayExpr<ArraySub<A, B, Ndim>, Ndim> {
}
auto operator[](int i) const { return arr1_[i] - arr2_[i]; }
size_t size() const { return arr1_.size(); }
std::array<int64_t, Ndim> shape() const { return arr1_.shape(); }
std::array<ssize_t, Ndim> shape() const { return arr1_.shape(); }
};
template <typename A, typename B, int64_t Ndim>
template <typename A, typename B, ssize_t Ndim>
class ArrayMul : public ArrayExpr<ArrayMul<A, B, Ndim>,Ndim> {
const A &arr1_;
const B &arr2_;
@ -55,10 +57,10 @@ class ArrayMul : public ArrayExpr<ArrayMul<A, B, Ndim>,Ndim> {
}
auto operator[](int i) const { return arr1_[i] * arr2_[i]; }
size_t size() const { return arr1_.size(); }
std::array<int64_t, Ndim> shape() const { return arr1_.shape(); }
std::array<ssize_t, Ndim> shape() const { return arr1_.shape(); }
};
template <typename A, typename B, int64_t Ndim>
template <typename A, typename B, ssize_t Ndim>
class ArrayDiv : public ArrayExpr<ArrayDiv<A, B, Ndim>, Ndim> {
const A &arr1_;
const B &arr2_;
@ -69,27 +71,27 @@ class ArrayDiv : public ArrayExpr<ArrayDiv<A, B, Ndim>, Ndim> {
}
auto operator[](int i) const { return arr1_[i] / arr2_[i]; }
size_t size() const { return arr1_.size(); }
std::array<int64_t, Ndim> shape() const { return arr1_.shape(); }
std::array<ssize_t, Ndim> shape() const { return arr1_.shape(); }
};
template <typename A, typename B, int64_t Ndim>
template <typename A, typename B, ssize_t Ndim>
auto operator+(const ArrayExpr<A, Ndim> &arr1, const ArrayExpr<B, Ndim> &arr2) {
return ArrayAdd<ArrayExpr<A, Ndim>, ArrayExpr<B, Ndim>, Ndim>(arr1, arr2);
}
template <typename A, typename B, int64_t Ndim>
template <typename A, typename B, ssize_t Ndim>
auto operator-(const ArrayExpr<A,Ndim> &arr1, const ArrayExpr<B, Ndim> &arr2) {
return ArraySub<ArrayExpr<A, Ndim>, ArrayExpr<B, Ndim>, Ndim>(arr1, arr2);
}
template <typename A, typename B, int64_t Ndim>
template <typename A, typename B, ssize_t Ndim>
auto operator*(const ArrayExpr<A, Ndim> &arr1, const ArrayExpr<B, Ndim> &arr2) {
return ArrayMul<ArrayExpr<A, Ndim>, ArrayExpr<B, Ndim>, Ndim>(arr1, arr2);
}
template <typename A, typename B, int64_t Ndim>
template <typename A, typename B, ssize_t Ndim>
auto operator/(const ArrayExpr<A, Ndim> &arr1, const ArrayExpr<B, Ndim> &arr2) {
return ArrayDiv<ArrayExpr<A, Ndim>, ArrayExpr<B, Ndim>, Ndim>(arr1, arr2);
}

170
include/aare/CalculateEta.hpp Normal file
View File

@ -0,0 +1,170 @@
#pragma once
#include "aare/Cluster.hpp"
#include "aare/ClusterVector.hpp"
#include "aare/NDArray.hpp"
namespace aare {
enum class corner : int {
cBottomLeft = 0,
cBottomRight = 1,
cTopLeft = 2,
cTopRight = 3
};
enum class pixel : int {
pBottomLeft = 0,
pBottom = 1,
pBottomRight = 2,
pLeft = 3,
pCenter = 4,
pRight = 5,
pTopLeft = 6,
pTop = 7,
pTopRight = 8
};
template <typename T> struct Eta2 {
double x;
double y;
int c;
T sum;
};
/**
* @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});
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;
}
return eta2;
}
/**
* @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.
*/
template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
typename CoordType>
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;
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);
// 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]));
}
eta.c = c; // TODO only supported for 2x2 and 3x3 clusters -> at least no
// underyling enum class
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) {
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]);
if ((cl.data[0] + cl.data[2]) != 0)
eta.y = static_cast<double>(cl.data[2]) / (cl.data[0] + cl.data[2]);
eta.sum = cl.sum();
eta.c = static_cast<int>(corner::cBottomLeft); // TODO! This is not correct,
// but need to put something
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) {
Eta2<T> eta{};
T sum = 0;
std::for_each(std::begin(cl.data), std::end(cl.data),
[&sum](T x) { sum += x; });
eta.sum = sum;
eta.c = corner::cBottomLeft;
if ((cl.data[3] + cl.data[4] + cl.data[5]) != 0)
eta.x = static_cast<double>(-cl.data[3] + cl.data[3 + 2]) /
(cl.data[3] + cl.data[4] + cl.data[5]);
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]);
return eta;
}
} // namespace aare

97
include/aare/CircularFifo.hpp Normal file
View File

@ -0,0 +1,97 @@
#pragma once
#include <chrono>
#include <fmt/color.h>
#include <fmt/format.h>
#include <memory>
#include <thread>
#include "aare/ProducerConsumerQueue.hpp"
namespace aare {
template <class ItemType> class CircularFifo {
uint32_t fifo_size;
aare::ProducerConsumerQueue<ItemType> free_slots;
aare::ProducerConsumerQueue<ItemType> filled_slots;
public:
CircularFifo() : CircularFifo(100){};
CircularFifo(uint32_t size) : fifo_size(size), free_slots(size + 1), filled_slots(size + 1) {
// TODO! how do we deal with alignment for writing? alignas???
// Do we give the user a chance to provide memory locations?
// Templated allocator?
for (size_t i = 0; i < fifo_size; ++i) {
free_slots.write(ItemType{});
}
}
bool next() {
// TODO! avoid default constructing ItemType
ItemType it;
if (!filled_slots.read(it))
return false;
if (!free_slots.write(std::move(it)))
return false;
return true;
}
~CircularFifo() {}
using value_type = ItemType;
auto numFilledSlots() const noexcept { return filled_slots.sizeGuess(); }
auto numFreeSlots() const noexcept { return free_slots.sizeGuess(); }
auto isFull() const noexcept { return filled_slots.isFull(); }
ItemType pop_free() {
ItemType v;
while (!free_slots.read(v))
;
return std::move(v);
// return v;
}
bool try_pop_free(ItemType &v) { return free_slots.read(v); }
ItemType pop_value(std::chrono::nanoseconds wait, std::atomic<bool> &stopped) {
ItemType v;
while (!filled_slots.read(v) && !stopped) {
std::this_thread::sleep_for(wait);
}
return std::move(v);
}
ItemType pop_value() {
ItemType v;
while (!filled_slots.read(v))
;
return std::move(v);
}
ItemType *frontPtr() { return filled_slots.frontPtr(); }
// TODO! Add function to move item from filled to free to be used
// with the frontPtr function
template <class... Args> void push_value(Args &&...recordArgs) {
while (!filled_slots.write(std::forward<Args>(recordArgs)...))
;
}
template <class... Args> bool try_push_value(Args &&...recordArgs) {
return filled_slots.write(std::forward<Args>(recordArgs)...);
}
template <class... Args> void push_free(Args &&...recordArgs) {
while (!free_slots.write(std::forward<Args>(recordArgs)...))
;
}
template <class... Args> bool try_push_free(Args &&...recordArgs) {
return free_slots.write(std::forward<Args>(recordArgs)...);
}
};
} // namespace aare

86
include/aare/Cluster.hpp Normal file
View File

@ -0,0 +1,86 @@
/************************************************
* @file Cluster.hpp
* @short definition of cluster, where CoordType (x,y) give
* the cluster center coordinates and data the actual cluster data
* cluster size is given as template parameters
***********************************************/
#pragma once
#include <algorithm>
#include <array>
#include <cstdint>
#include <numeric>
#include <type_traits>
namespace aare {
// requires clause c++20 maybe update
template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
typename CoordType = int16_t>
struct Cluster {
static_assert(std::is_arithmetic_v<T>, "T needs to be an arithmetic type");
static_assert(std::is_integral_v<CoordType>,
"CoordType needs to be an integral type");
static_assert(ClusterSizeX > 0 && ClusterSizeY > 0,
"Cluster sizes must be bigger than zero");
CoordType x;
CoordType y;
std::array<T, ClusterSizeX * ClusterSizeY> data;
static constexpr uint8_t cluster_size_x = ClusterSizeX;
static constexpr uint8_t cluster_size_y = ClusterSizeY;
using value_type = T;
using coord_type = CoordType;
T sum() const { return std::accumulate(data.begin(), data.end(), T{}); }
std::pair<T, int> max_sum_2x2() const {
if constexpr (cluster_size_x == 3 && cluster_size_y == 3) {
std::array<T, 4> sum_2x2_subclusters;
sum_2x2_subclusters[0] = data[0] + data[1] + data[3] + data[4];
sum_2x2_subclusters[1] = data[1] + data[2] + data[4] + data[5];
sum_2x2_subclusters[2] = data[3] + data[4] + data[6] + data[7];
sum_2x2_subclusters[3] = data[4] + data[5] + data[7] + data[8];
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);
} else if constexpr (cluster_size_x == 2 && cluster_size_y == 2) {
return std::make_pair(data[0] + data[1] + data[2] + data[3], 0);
} else {
constexpr size_t num_2x2_subclusters =
(ClusterSizeX - 1) * (ClusterSizeY - 1);
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];
}
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);
}
}
};
// Type Traits for is_cluster_type
template <typename T>
struct is_cluster : std::false_type {}; // Default case: Not a Cluster
template <typename T, uint8_t X, uint8_t Y, typename CoordType>
struct is_cluster<Cluster<T, X, Y, CoordType>> : std::true_type {}; // Cluster
template <typename T> constexpr bool is_cluster_v = is_cluster<T>::value;
} // namespace aare

58
include/aare/ClusterCollector.hpp Normal file
View File

@ -0,0 +1,58 @@
#pragma once
#include <atomic>
#include <thread>
#include "aare/ClusterFinderMT.hpp"
#include "aare/ClusterVector.hpp"
#include "aare/ProducerConsumerQueue.hpp"
namespace aare {
template <typename ClusterType,
typename = std::enable_if_t<is_cluster_v<ClusterType>>>
class ClusterCollector {
ProducerConsumerQueue<ClusterVector<ClusterType>> *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;
fmt::print("ClusterCollector started\n");
while (!m_stop_requested || !m_source->isEmpty()) {
if (ClusterVector<ClusterType> *clusters = m_source->frontPtr();
clusters != nullptr) {
m_clusters.push_back(std::move(*clusters));
m_source->popFront();
} else {
std::this_thread::sleep_for(m_default_wait);
}
}
fmt::print("ClusterCollector stopped\n");
m_stopped = true;
}
public:
ClusterCollector(ClusterFinderMT<ClusterType, uint16_t, double> *source) {
m_source = source->sink();
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?
}
void stop() {
m_stop_requested = true;
m_thread.join();
}
std::vector<ClusterVector<ClusterType>> steal_clusters() {
if (!m_stopped) {
throw std::runtime_error("ClusterCollector is still running");
}
return std::move(m_clusters);
}
};
} // namespace aare

478
include/aare/ClusterFile.hpp
View File

@ -1,67 +1,449 @@
#pragma once
#include "aare/Cluster.hpp"
#include "aare/ClusterVector.hpp"
#include "aare/GainMap.hpp"
#include "aare/NDArray.hpp"
#include "aare/defs.hpp"
#include <filesystem>
#include <fstream>
#include <optional>
namespace aare {
struct Cluster {
int16_t x;
int16_t y;
int32_t data[9];
};
typedef enum {
cBottomLeft = 0,
cBottomRight = 1,
cTopLeft = 2,
cTopRight = 3
} corner;
typedef enum {
pBottomLeft = 0,
pBottom = 1,
pBottomRight = 2,
pLeft = 3,
pCenter = 4,
pRight = 5,
pTopLeft = 6,
pTop = 7,
pTopRight = 8
} pixel;
struct ClusterAnalysis {
uint32_t c;
int32_t tot;
double etax;
double etay;
};
/*
Binary cluster file. Expects data to be layed 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:
*
*
* 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.
*/
template <typename ClusterType,
typename Enable = std::enable_if_t<is_cluster_v<ClusterType>>>
class ClusterFile {
FILE *fp{};
uint32_t m_num_left{};
size_t m_chunk_size{};
const 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*/
std::optional<ROI> m_roi; /*Region of interest, will be applied if set*/
std::optional<NDArray<int32_t, 2>>
m_noise_map; /*Noise map to cut photons, will be applied if set*/
std::optional<InvertedGainMap> m_gain_map; /*Gain map to apply to the
clusters, will be applied if set*/
public:
ClusterFile(const std::filesystem::path &fname, size_t chunk_size = 1000);
~ClusterFile();
std::vector<Cluster> read_clusters(size_t n_clusters);
std::vector<Cluster> read_frame(int32_t &out_fnum);
std::vector<Cluster>
read_cluster_with_cut(size_t n_clusters, double *noise_map, int nx, int ny);
/**
* @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
*/
ClusterFile(const std::filesystem::path &fname, size_t chunk_size = 1000,
const std::string &mode = "r")
int analyze_data(int32_t *data, int32_t *t2, int32_t *t3, char *quad,
double *eta2x, double *eta2y, double *eta3x, double *eta3y);
int analyze_cluster(Cluster cl, int32_t *t2, int32_t *t3, char *quad,
double *eta2x, double *eta2y, double *eta3x,
double *eta3y);
: 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);
}
}
~ClusterFile() { close(); }
/**
* @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
*/
ClusterVector<ClusterType> read_clusters(size_t n_clusters) {
if (m_mode != "r") {
throw std::runtime_error("File not opened for reading");
}
if (m_noise_map || m_roi) {
return read_clusters_with_cut(n_clusters);
} else {
return read_clusters_without_cut(n_clusters);
}
}
/**
* @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
*/
ClusterVector<ClusterType> read_frame() {
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();
}
}
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);
uint32_t n_clusters = clusters.size();
fwrite(&n_clusters, sizeof(n_clusters), 1, fp);
fwrite(clusters.data(), clusters.item_size(), clusters.size(), fp);
}
/**
* @brief Return the chunk size
*/
size_t chunk_size() const { return m_chunk_size; }
void close();
/**
* @brief Set the region of interest to use when reading
* clusters. If set only clusters within the ROI will be
* read.
*/
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.
*/
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.
*/
void set_gain_map(const NDView<double, 2> gain_map) {
m_gain_map = InvertedGainMap(gain_map);
}
void set_gain_map(const InvertedGainMap &gain_map) {
m_gain_map = gain_map;
}
void set_gain_map(const InvertedGainMap &&gain_map) {
m_gain_map = gain_map;
}
/**
* @brief Close the file. If not closed the file will be
* closed in the destructor
*/
void close() {
if (fp) {
fclose(fp);
fp = nullptr;
}
}
/** @brief Open the file in specific mode
*
*/
void open(const std::string &mode) {
if (fp) {
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_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_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_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 is_selected(ClusterType &cl);
ClusterType read_one_cluster();
};
template <typename ClusterType, typename Enable>
ClusterVector<ClusterType>
ClusterFile<ClusterType, Enable>::read_clusters_without_cut(size_t n_clusters) {
if (m_mode != "r") {
throw std::runtime_error("File not opened for reading");
}
ClusterVector<ClusterType> clusters(n_clusters);
clusters.resize(n_clusters);
int32_t iframe = 0; // frame number needs to be 4 bytes!
size_t nph_read = 0;
uint32_t nn = m_num_left;
uint32_t nph = m_num_left; // number of clusters in frame needs to be 4
auto buf = clusters.data();
// if there are photons left from previous frame read them first
if (nph) {
if (nph > n_clusters) {
// if we have more photons left in the frame then photons to
// read we read directly the requested number
nn = n_clusters;
} else {
nn = nph;
}
nph_read += fread((buf + nph_read), clusters.item_size(), nn, fp);
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)) {
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;
nph_read +=
fread((buf + nph_read), clusters.item_size(), nn, fp);
m_num_left = nph - nn;
}
if (nph_read >= n_clusters)
break;
}
}
// Resize the vector to the number o f clusters.
// No new allocation, only change bounds.
clusters.resize(nph_read);
if (m_gain_map)
m_gain_map->apply_gain_map(clusters);
return clusters;
}
template <typename ClusterType, typename Enable>
ClusterVector<ClusterType>
ClusterFile<ClusterType, Enable>::read_clusters_with_cut(size_t n_clusters) {
ClusterVector<ClusterType> clusters;
clusters.reserve(n_clusters);
// if there are photons left from previous frame read them first
if (m_num_left) {
while (m_num_left && clusters.size() < n_clusters) {
ClusterType c = read_one_cluster();
if (is_selected(c)) {
clusters.push_back(c);
}
}
}
// we did not have enough clusters left in the previous frame
// keep on reading frames until reaching n_clusters
if (clusters.size() < n_clusters) {
// sanity check
if (m_num_left) {
throw std::runtime_error(
LOCATION + "Entered second loop with clusters left\n");
}
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);
}
}
}
// we have enough clusters, break out of the outer while loop
if (clusters.size() >= n_clusters)
break;
}
}
if (m_gain_map)
m_gain_map->apply_gain_map(clusters);
return clusters;
}
template <typename ClusterType, typename Enable>
ClusterType ClusterFile<ClusterType, Enable>::read_one_cluster() {
ClusterType c;
auto rc = fread(&c, sizeof(c), 1, fp);
if (rc != 1) {
throw std::runtime_error(LOCATION + "Could not read cluster");
}
--m_num_left;
return c;
}
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");
}
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");
}
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");
}
ClusterVector<ClusterType> clusters(n_clusters);
clusters.set_frame_number(frame_number);
clusters.resize(n_clusters);
if (fread(clusters.data(), clusters.item_size(), n_clusters, fp) !=
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;
}
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");
}
int32_t frame_number;
if (fread(&frame_number, sizeof(frame_number), 1, fp) != 1) {
throw std::runtime_error("Could not read frame number");
}
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);
clusters.set_frame_number(frame_number);
while (m_num_left) {
ClusterType c = read_one_cluster();
if (is_selected(c)) {
clusters.push_back(c);
}
}
if (m_gain_map)
m_gain_map->apply_gain_map(clusters);
return clusters;
}
template <typename ClusterType, typename Enable>
bool ClusterFile<ClusterType, Enable>::is_selected(ClusterType &cl) {
// Should fail fast
if (m_roi) {
if (!(m_roi->contains(cl.x, cl.y))) {
return false;
}
}
size_t cluster_center_index =
(ClusterType::cluster_size_x / 2) +
(ClusterType::cluster_size_y / 2) * ClusterType::cluster_size_x;
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 noise =
(*m_noise_map)(cl.y, cl.x); // TODO! check if this is correct
if (sum_1x1 <= noise || sum_2x2 <= 2 * noise ||
total_sum <= 3 * noise) {
return false;
}
}
// we passed all checks
return true;
}
} // namespace aare

62
include/aare/ClusterFileSink.hpp Normal file
View File

@ -0,0 +1,62 @@
#pragma once
#include <atomic>
#include <filesystem>
#include <thread>
#include "aare/ClusterFinderMT.hpp"
#include "aare/ClusterVector.hpp"
#include "aare/ProducerConsumerQueue.hpp"
namespace aare {
template <typename ClusterType,
typename = std::enable_if_t<is_cluster_v<ClusterType>>>
class ClusterFileSink {
ProducerConsumerQueue<ClusterVector<ClusterType>> *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::ofstream m_file;
void process() {
m_stopped = false;
fmt::print("ClusterFileSink started\n");
while (!m_stop_requested || !m_source->isEmpty()) {
if (ClusterVector<ClusterType> *clusters = m_source->frontPtr();
clusters != nullptr) {
// Write clusters to file
int32_t frame_number =
clusters->frame_number(); // TODO! Should we store frame
// number already as int?
uint32_t num_clusters = clusters->size();
m_file.write(reinterpret_cast<const char *>(&frame_number),
sizeof(frame_number));
m_file.write(reinterpret_cast<const char *>(&num_clusters),
sizeof(num_clusters));
m_file.write(reinterpret_cast<const char *>(clusters->data()),
clusters->size() * clusters->item_size());
m_source->popFront();
} else {
std::this_thread::sleep_for(m_default_wait);
}
}
fmt::print("ClusterFileSink stopped\n");
m_stopped = true;
}
public:
ClusterFileSink(ClusterFinderMT<ClusterType, uint16_t, double> *source,
const std::filesystem::path &fname) {
m_source = source->sink();
m_thread = std::thread(&ClusterFileSink::process, this);
m_file.open(fname, std::ios::binary);
}
void stop() {
m_stop_requested = true;
m_thread.join();
m_file.close();
}
};
} // namespace aare

148
include/aare/ClusterFileV2.hpp
View File

@ -1,148 +0,0 @@
#pragma once
#include "aare/core/defs.hpp"
#include <filesystem>
#include <string>
#include <fmt/format.h>
namespace aare {
struct ClusterHeader {
int32_t frame_number;
int32_t n_clusters;
std::string to_string() const {
return "frame_number: " + std::to_string(frame_number) + ", n_clusters: " + std::to_string(n_clusters);
}
};
struct ClusterV2_ {
int16_t x;
int16_t y;
std::array<int32_t, 9> data;
std::string to_string(bool detailed = false) const {
if (detailed) {
std::string data_str = "[";
for (auto &d : data) {
data_str += std::to_string(d) + ", ";
}
data_str += "]";
return "x: " + std::to_string(x) + ", y: " + std::to_string(y) + ", data: " + data_str;
}
return "x: " + std::to_string(x) + ", y: " + std::to_string(y);
}
};
struct ClusterV2 {
ClusterV2_ cluster;
int32_t frame_number;
std::string to_string() const {
return "frame_number: " + std::to_string(frame_number) + ", " + cluster.to_string();
}
};
/**
* @brief
* important not: fp always points to the clusters header and does not point to individual clusters
*
*/
class ClusterFileV2 {
std::filesystem::path m_fpath;
std::string m_mode;
FILE *fp{nullptr};
void check_open(){
if (!fp)
throw std::runtime_error(fmt::format("File: {} not open", m_fpath.string()));
}
public:
ClusterFileV2(std::filesystem::path const &fpath, std::string const &mode): m_fpath(fpath), m_mode(mode) {
if (m_mode != "r" && m_mode != "w")
throw std::invalid_argument("mode must be 'r' or 'w'");
if (m_mode == "r" && !std::filesystem::exists(m_fpath))
throw std::invalid_argument("File does not exist");
if (mode == "r") {
fp = fopen(fpath.string().c_str(), "rb");
} else if (mode == "w") {
if (std::filesystem::exists(fpath)) {
fp = fopen(fpath.string().c_str(), "r+b");
} else {
fp = fopen(fpath.string().c_str(), "wb");
}
}
if (fp == nullptr) {
throw std::runtime_error("Failed to open file");
}
}
~ClusterFileV2() { close(); }
std::vector<ClusterV2> read() {
check_open();
ClusterHeader header;
fread(&header, sizeof(ClusterHeader), 1, fp);
std::vector<ClusterV2_> clusters_(header.n_clusters);
fread(clusters_.data(), sizeof(ClusterV2_), header.n_clusters, fp);
std::vector<ClusterV2> clusters;
for (auto &c : clusters_) {
ClusterV2 cluster;
cluster.cluster = std::move(c);
cluster.frame_number = header.frame_number;
clusters.push_back(cluster);
}
return clusters;
}
std::vector<std::vector<ClusterV2>> read(int n_frames) {
std::vector<std::vector<ClusterV2>> clusters;
for (int i = 0; i < n_frames; i++) {
clusters.push_back(read());
}
return clusters;
}
size_t write(std::vector<ClusterV2> const &clusters) {
check_open();
if (m_mode != "w")
throw std::runtime_error("File not opened in write mode");
if (clusters.empty())
return 0;
ClusterHeader header;
header.frame_number = clusters[0].frame_number;
header.n_clusters = clusters.size();
fwrite(&header, sizeof(ClusterHeader), 1, fp);
for (auto &c : clusters) {
fwrite(&c.cluster, sizeof(ClusterV2_), 1, fp);
}
return clusters.size();
}
size_t write(std::vector<std::vector<ClusterV2>> const &clusters) {
check_open();
if (m_mode != "w")
throw std::runtime_error("File not opened in write mode");
size_t n_clusters = 0;
for (auto &c : clusters) {
n_clusters += write(c);
}
return n_clusters;
}
int seek_to_begin() { return fseek(fp, 0, SEEK_SET); }
int seek_to_end() { return fseek(fp, 0, SEEK_END); }
int32_t frame_number() {
auto pos = ftell(fp);
ClusterHeader header;
fread(&header, sizeof(ClusterHeader), 1, fp);
fseek(fp, pos, SEEK_SET);
return header.frame_number;
}
void close() {
if (fp) {
fclose(fp);
fp = nullptr;
}
}
};
} // namespace aare

314
include/aare/ClusterFinder.hpp
View File

@ -1,4 +1,6 @@
#pragma once
#include "aare/ClusterFile.hpp"
#include "aare/ClusterVector.hpp"
#include "aare/Dtype.hpp"
#include "aare/NDArray.hpp"
#include "aare/NDView.hpp"
@ -8,251 +10,147 @@
namespace aare {
/** enum to define the event types */
enum eventType {
PEDESTAL, /** pedestal */
NEIGHBOUR, /** neighbour i.e. below threshold, but in the cluster of a
photon */
PHOTON, /** photon i.e. above threshold */
PHOTON_MAX, /** maximum of a cluster satisfying the photon conditions */
NEGATIVE_PEDESTAL, /** negative value, will not be accounted for as pedestal
in order to avoid drift of the pedestal towards
negative values */
UNDEFINED_EVENT = -1 /** undefined */
};
template <typename FRAME_TYPE = uint16_t, typename PEDESTAL_TYPE = double>
template <typename ClusterType = Cluster<int32_t, 3, 3>,
typename FRAME_TYPE = uint16_t, typename PEDESTAL_TYPE = double>
class ClusterFinder {
Shape<2> m_image_size;
const int m_cluster_sizeX;
const int m_cluster_sizeY;
const double m_threshold;
const double m_nSigma;
const double c2;
const double c3;
const PEDESTAL_TYPE m_nSigma;
const PEDESTAL_TYPE c2;
const PEDESTAL_TYPE c3;
Pedestal<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;
public:
ClusterFinder(Shape<2> image_size, Shape<2>cluster_size, double nSigma = 5.0,
double threshold = 0.0)
: m_image_size(image_size), m_cluster_sizeX(cluster_size[0]), m_cluster_sizeY(cluster_size[1]),
m_threshold(threshold), m_nSigma(nSigma),
c2(sqrt((m_cluster_sizeY + 1) / 2 * (m_cluster_sizeX + 1) / 2)),
c3(sqrt(m_cluster_sizeX * m_cluster_sizeY)),
m_pedestal(image_size[0], image_size[1]) {
// c2 = sqrt((cluster_sizeY + 1) / 2 * (cluster_sizeX + 1) / 2);
// c3 = sqrt(cluster_sizeX * cluster_sizeY);
};
/**
* @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
*
*/
ClusterFinder(Shape<2> image_size, PEDESTAL_TYPE nSigma = 5.0,
size_t capacity = 1000000)
: m_image_size(image_size), m_nSigma(nSigma),
c2(sqrt((ClusterSizeY + 1) / 2 * (ClusterSizeX + 1) / 2)),
c3(sqrt(ClusterSizeX * ClusterSizeY)),
m_pedestal(image_size[0], image_size[1]), m_clusters(capacity) {};
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> 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
*
*/
ClusterVector<ClusterType>
steal_clusters(bool realloc_same_capacity = false) {
ClusterVector<ClusterType> tmp = std::move(m_clusters);
if (realloc_same_capacity)
m_clusters = ClusterVector<ClusterType>(tmp.capacity());
else
m_clusters = ClusterVector<ClusterType>{};
return tmp;
}
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;
std::vector<DynamicCluster>
find_clusters_without_threshold(NDView<FRAME_TYPE, 2> frame,
// Pedestal<PEDESTAL_TYPE> &pedestal,
bool late_update = false) {
struct pedestal_update {
int x;
int y;
FRAME_TYPE value;
};
std::vector<pedestal_update> pedestal_updates;
std::vector<DynamicCluster> clusters;
std::vector<std::vector<eventType>> eventMask;
for (int i = 0; i < frame.shape(0); i++) {
eventMask.push_back(std::vector<eventType>(frame.shape(1)));
}
long double val;
long double max;
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++) {
// initialize max and total
max = std::numeric_limits<FRAME_TYPE>::min();
long double total = 0;
eventMask[iy][ix] = PEDESTAL;
for (short ir = -(m_cluster_sizeY / 2);
ir < (m_cluster_sizeY / 2) + 1; ir++) {
for (short ic = -(m_cluster_sizeX / 2);
ic < (m_cluster_sizeX / 2) + 1; ic++) {
PEDESTAL_TYPE max = std::numeric_limits<FRAME_TYPE>::min();
PEDESTAL_TYPE total = 0;
// 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));
if (value < -m_nSigma * rms)
continue; // NEGATIVE_PEDESTAL go to next pixel
// TODO! No pedestal update???
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)) {
val = frame(iy + ir, ix + ic) -
m_pedestal.mean(iy + ir, ix + ic);
PEDESTAL_TYPE val =
frame(iy + ir, ix + ic) -
m_pedestal.mean(iy + ir, ix + ic);
total += val;
if (val > max) {
max = val;
}
max = std::max(max, val);
}
}
}
auto rms = m_pedestal.std(iy, ix);
if (frame(iy, ix) - m_pedestal.mean(iy, ix) < -m_nSigma * rms) {
eventMask[iy][ix] = NEGATIVE_PEDESTAL;
continue;
} else if (max > m_nSigma * rms) {
eventMask[iy][ix] = PHOTON;
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) {
eventMask[iy][ix] = PHOTON;
// pass
} else {
if (late_update) {
pedestal_updates.push_back({ix, iy, frame(iy, ix)});
} else {
m_pedestal.push(iy, ix, frame(iy, ix));
}
continue;
// 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
}
if (eventMask[iy][ix] == PHOTON &&
(frame(iy, ix) - m_pedestal.mean(iy, ix)) >= max) {
eventMask[iy][ix] = PHOTON_MAX;
DynamicCluster cluster(m_cluster_sizeX, m_cluster_sizeY,
Dtype(typeid(PEDESTAL_TYPE)));
// Store cluster
if (value == max) {
ClusterType cluster{};
cluster.x = ix;
cluster.y = iy;
short i = 0;
for (short ir = -(m_cluster_sizeY / 2);
ir < (m_cluster_sizeY / 2) + 1; ir++) {
for (short ic = -(m_cluster_sizeX / 2);
ic < (m_cluster_sizeX / 2) + 1; ic++) {
// 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)) {
PEDESTAL_TYPE tmp =
static_cast<PEDESTAL_TYPE>(
frame(iy + ir, ix + ic)) -
m_pedestal.mean(iy + ir, ix + ic);
cluster.set<PEDESTAL_TYPE>(i, tmp);
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++;
}
}
}
clusters.push_back(cluster);
// Add the cluster to the output ClusterVector
m_clusters.push_back(cluster);
}
}
}
if (late_update) {
for (auto &update : pedestal_updates) {
m_pedestal.push(update.y, update.x, update.value);
}
}
return clusters;
}
// template <typename FRAME_TYPE, typename PEDESTAL_TYPE>
std::vector<DynamicCluster>
find_clusters_with_threshold(NDView<FRAME_TYPE, 2> frame,
Pedestal<PEDESTAL_TYPE> &pedestal) {
assert(m_threshold > 0);
std::vector<DynamicCluster> clusters;
std::vector<std::vector<eventType>> eventMask;
for (int i = 0; i < frame.shape(0); i++) {
eventMask.push_back(std::vector<eventType>(frame.shape(1)));
}
double tthr, tthr1, tthr2;
NDArray<FRAME_TYPE, 2> rest({frame.shape(0), frame.shape(1)});
NDArray<int, 2> nph({frame.shape(0), frame.shape(1)});
// convert to n photons
// nph = (frame-pedestal.mean()+0.5*m_threshold)/m_threshold; // can be
// optimized with expression templates?
for (int iy = 0; iy < frame.shape(0); iy++) {
for (int ix = 0; ix < frame.shape(1); ix++) {
auto val = frame(iy, ix) - pedestal.mean(iy, ix);
nph(iy, ix) = (val + 0.5 * m_threshold) / m_threshold;
nph(iy, ix) = nph(iy, ix) < 0 ? 0 : nph(iy, ix);
rest(iy, ix) = val - nph(iy, ix) * m_threshold;
}
}
// iterate over frame pixels
for (int iy = 0; iy < frame.shape(0); iy++) {
for (int ix = 0; ix < frame.shape(1); ix++) {
eventMask[iy][ix] = PEDESTAL;
// initialize max and total
FRAME_TYPE max = std::numeric_limits<FRAME_TYPE>::min();
long double total = 0;
if (rest(iy, ix) <= 0.25 * m_threshold) {
pedestal.push(iy, ix, frame(iy, ix));
continue;
}
eventMask[iy][ix] = NEIGHBOUR;
// iterate over cluster pixels around the current pixel (ix,iy)
for (short ir = -(m_cluster_sizeY / 2);
ir < (m_cluster_sizeY / 2) + 1; ir++) {
for (short ic = -(m_cluster_sizeX / 2);
ic < (m_cluster_sizeX / 2) + 1; ic++) {
if (ix + ic >= 0 && ix + ic < frame.shape(1) &&
iy + ir >= 0 && iy + ir < frame.shape(0)) {
auto val = frame(iy + ir, ix + ic) -
pedestal.mean(iy + ir, ix + ic);
total += val;
if (val > max) {
max = val;
}
}
}
}
auto rms = pedestal.std(iy, ix);
if (m_nSigma == 0) {
tthr = m_threshold;
tthr1 = m_threshold;
tthr2 = m_threshold;
} else {
tthr = m_nSigma * rms;
tthr1 = m_nSigma * rms * c3;
tthr2 = m_nSigma * rms * c2;
if (m_threshold > 2 * tthr)
tthr = m_threshold - tthr;
if (m_threshold > 2 * tthr1)
tthr1 = tthr - tthr1;
if (m_threshold > 2 * tthr2)
tthr2 = tthr - tthr2;
}
if (total > tthr1 || max > tthr) {
eventMask[iy][ix] = PHOTON;
nph(iy, ix) += 1;
rest(iy, ix) -= m_threshold;
} else {
pedestal.push(iy, ix, frame(iy, ix));
continue;
}
if (eventMask[iy][ix] == PHOTON &&
frame(iy, ix) - pedestal.mean(iy, ix) >= max) {
eventMask[iy][ix] = PHOTON_MAX;
DynamicCluster cluster(m_cluster_sizeX, m_cluster_sizeY,
Dtype(typeid(FRAME_TYPE)));
cluster.x = ix;
cluster.y = iy;
short i = 0;
for (short ir = -(m_cluster_sizeY / 2);
ir < (m_cluster_sizeY / 2) + 1; ir++) {
for (short ic = -(m_cluster_sizeX / 2);
ic < (m_cluster_sizeX / 2) + 1; ic++) {
if (ix + ic >= 0 && ix + ic < frame.shape(1) &&
iy + ir >= 0 && iy + ir < frame.shape(0)) {
auto tmp = frame(iy + ir, ix + ic) -
pedestal.mean(iy + ir, ix + ic);
cluster.set<FRAME_TYPE>(i, tmp);
i++;
}
}
}
clusters.push_back(cluster);
}
}
}
return clusters;
}
};

277
include/aare/ClusterFinderMT.hpp Normal file
View File

@ -0,0 +1,277 @@
#pragma once
#include <atomic>
#include <cstdint>
#include <memory>
#include <thread>
#include <vector>
#include "aare/ClusterFinder.hpp"
#include "aare/NDArray.hpp"
#include "aare/ProducerConsumerQueue.hpp"
namespace aare {
enum class FrameType {
DATA,
PEDESTAL,
};
struct FrameWrapper {
FrameType type;
uint64_t frame_number;
NDArray<uint16_t, 2> data;
};
/**
* @brief ClusterFinderMT is a multi-threaded version of ClusterFinder. It uses
* a producer-consumer queue to distribute the frames to the threads. The
* clusters are collected in a single output queue.
* @tparam FRAME_TYPE type of the frame data
* @tparam PEDESTAL_TYPE type of the pedestal data
* @tparam CT type of the cluster data
*/
template <typename ClusterType = Cluster<int32_t, 3, 3>,
typename FRAME_TYPE = uint16_t, typename PEDESTAL_TYPE = double>
class ClusterFinderMT {
protected:
using CT = typename ClusterType::value_type;
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 OutputQueue = ProducerConsumerQueue<ClusterVector<ClusterType>>;
std::vector<std::unique_ptr<InputQueue>> m_input_queues;
std::vector<std::unique_ptr<OutputQueue>> m_output_queues;
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};
private:
std::atomic<bool> m_stop_requested{false};
std::atomic<bool> m_processing_threads_stopped{true};
/**
* @brief Function called by the processing threads. It reads the frames
* from the input queue and processes them.
*/
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;
while (!m_stop_requested || !q->isEmpty()) {
if (FrameWrapper *frame = q->frontPtr(); frame != nullptr) {
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));
break;
case FrameType::PEDESTAL:
m_cluster_finders[thread_id]->push_pedestal_frame(
frame->data.view());
break;
}
// frame is processed now discard it
m_input_queues[thread_id]->popFront();
} else {
std::this_thread::sleep_for(m_default_wait);
}
}
}
/**
* @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;
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);
}
queue->popFront();
empty = false;
}
}
}
}
public:
/**
* @brief Construct a new ClusterFinderMT object
* @param image_size size of the image
* @param cluster_size size of the cluster
* @param nSigma number of sigma above the pedestal to consider a photon
* @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
*/
ClusterFinderMT(Shape<2> image_size, PEDESTAL_TYPE nSigma = 5.0,
size_t capacity = 2000, size_t n_threads = 3)
: m_n_threads(n_threads) {
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));
}
for (size_t i = 0; i < n_threads; i++) {
m_input_queues.emplace_back(std::make_unique<InputQueue>(200));
m_output_queues.emplace_back(std::make_unique<OutputQueue>(200));
}
// TODO! Should we start automatically?
start();
}
/**
* @brief Return the sink queue where all the clusters are collected
* @warning You need to empty this queue otherwise the cluster finder will
* wait forever
*/
ProducerConsumerQueue<ClusterVector<ClusterType>> *sink() {
return &m_sink;
}
/**
* @brief Start all processing threads
*/
void start() {
m_processing_threads_stopped = false;
m_stop_requested = false;
for (size_t i = 0; i < m_n_threads; i++) {
m_threads.push_back(
std::thread(&ClusterFinderMT::process, this, i));
}
m_collect_thread = std::thread(&ClusterFinderMT::collect, this);
}
/**
* @brief Stop all processing threads
*/
void stop() {
m_stop_requested = true;
for (auto &thread : m_threads) {
thread.join();
}
m_threads.clear();
m_processing_threads_stopped = true;
m_collect_thread.join();
}
/**
* @brief Wait for all the queues to be empty. Mostly used for timing tests.
*/
void sync() {
for (auto &q : m_input_queues) {
while (!q->isEmpty()) {
std::this_thread::sleep_for(m_default_wait);
}
}
for (auto &q : m_output_queues) {
while (!q->isEmpty()) {
std::this_thread::sleep_for(m_default_wait);
}
}
while (!m_sink.isEmpty()) {
std::this_thread::sleep_for(m_default_wait);
}
}
/**
* @brief Push a pedestal frame to all the cluster finders. The frames is
* 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 (auto &queue : m_input_queues) {
while (!queue->write(fw)) {
std::this_thread::sleep_for(m_default_wait);
}
}
}
/**
* @brief Push the frame to the queue of the next available thread. Function
* returns once the frame is in a queue.
* @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);
}
m_current_thread++;
}
void clear_pedestal() {
if (!m_processing_threads_stopped) {
throw std::runtime_error("ClusterFinderMT is still running");
}
for (auto &cf : m_cluster_finders) {
cf->clear_pedestal();
}
}
/**
* @brief Return the pedestal currently used by the cluster finder
* @param thread_index index of the thread
*/
auto pedestal(size_t thread_index = 0) {
if (m_cluster_finders.empty()) {
throw std::runtime_error("No cluster finders available");
}
if (!m_processing_threads_stopped) {
throw std::runtime_error("ClusterFinderMT is still running");
}
if (thread_index >= m_cluster_finders.size()) {
throw std::runtime_error("Thread index out of range");
}
return m_cluster_finders[thread_index]->pedestal();
}
/**
* @brief Return the noise currently used by the cluster finder
* @param thread_index index of the thread
*/
auto noise(size_t thread_index = 0) {
if (m_cluster_finders.empty()) {
throw std::runtime_error("No cluster finders available");
}
if (!m_processing_threads_stopped) {
throw std::runtime_error("ClusterFinderMT is still running");
}
if (thread_index >= m_cluster_finders.size()) {
throw std::runtime_error("Thread index out of range");
}
return m_cluster_finders[thread_index]->noise();
}
// void push(FrameWrapper&& frame) {
// //TODO! need to loop until we are successful
// auto rc = m_input_queue.write(std::move(frame));
// fmt::print("pushed frame {}\n", rc);
// }
};
} // namespace aare

170
include/aare/ClusterVector.hpp Normal file
View File

@ -0,0 +1,170 @@
#pragma once
#include "aare/Cluster.hpp" //TODO maybe store in seperate file !!!
#include <algorithm>
#include <array>
#include <cstddef>
#include <cstdint>
#include <numeric>
#include <vector>
#include <fmt/core.h>
#include "aare/Cluster.hpp"
#include "aare/NDView.hpp"
namespace aare {
template <typename ClusterType,
typename = std::enable_if_t<is_cluster_v<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
* @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 CoordType data type of the x and y coordinates of the cluster
* (normally int16_t)
*/
template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
typename CoordType>
class ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> {
std::vector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> m_data{};
int32_t m_frame_number{0}; // TODO! Check frame number size and type
public:
using value_type = T;
using ClusterType = Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>;
/**
* @brief Construct a new ClusterVector object
* @param capacity initial capacity of the buffer 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)
: 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();
}
// Move assignment operator
ClusterVector &operator=(ClusterVector &&other) noexcept {
if (this != &other) {
m_data = other.m_data;
m_frame_number = other.m_frame_number;
other.m_data.clear();
other.m_frame_number = 0;
}
return *this;
}
/**
* @brief Sum the pixels in each cluster
* @return std::vector<T> vector of sums for each cluster
*/
std::vector<T> sum() {
std::vector<T> sums(m_data.size());
std::transform(
m_data.begin(), m_data.end(), sums.begin(),
[](const ClusterType &cluster) { return cluster.sum(); });
return sums;
}
/**
* @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
*/
std::vector<T> sum_2x2() {
std::vector<T> 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;
});
return sums_2x2;
}
/**
* @brief Reserve space for at least capacity clusters
* @param capacity number of clusters to reserve space for
* @note If capacity is less than the current capacity, the function does
* nothing.
*/
void reserve(size_t capacity) { m_data.reserve(capacity); }
void resize(size_t size) { m_data.resize(size); }
void push_back(const ClusterType &cluster) { m_data.push_back(cluster); }
ClusterVector &operator+=(const ClusterVector &other) {
m_data.insert(m_data.end(), other.begin(), other.end());
return *this;
}
/**
* @brief Return the number of clusters in the vector
*/
size_t size() const { return m_data.size(); }
uint8_t cluster_size_x() const { return ClusterSizeX; }
uint8_t cluster_size_y() const { return ClusterSizeY; }
/**
* @brief Return the capacity of the buffer in number of clusters. This is
* the number of clusters that can be stored in the current buffer without
* reallocation.
*/
size_t capacity() const { return m_data.capacity(); }
auto begin() const { return m_data.begin(); }
auto end() const { return m_data.end(); }
/**
* @brief Return the size in bytes of a single cluster
*/
size_t item_size() const {
return sizeof(ClusterType); // 2 * sizeof(CoordType) + ClusterSizeX *
// ClusterSizeY * sizeof(T);
}
ClusterType *data() { return m_data.data(); }
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
*/
ClusterType &operator[](size_t i) { return m_data[i]; }
const ClusterType &operator[](size_t i) const { return m_data[i]; }
/**
* @brief Return the frame number of the clusters. 0 is used to indicate
* that the clusters come from many frames
*/
int32_t frame_number() const { return m_frame_number; }
void set_frame_number(int32_t frame_number) {
m_frame_number = frame_number;
}
};
} // namespace aare

3
include/aare/File.hpp
View File

@ -36,6 +36,8 @@ class File {
File(File &&other) noexcept;
File& operator=(File &&other) noexcept;
~File() = default;
// void close(); //!< close the file
Frame read_frame(); //!< read one frame from the file at the current position
Frame read_frame(size_t frame_index); //!< read one frame at the position given by frame number
@ -44,6 +46,7 @@ class File {
void read_into(std::byte *image_buf);
void read_into(std::byte *image_buf, size_t n_frames);
size_t frame_number(); //!< get the frame number at the current position
size_t frame_number(size_t frame_index); //!< get the frame number at the given frame index
size_t bytes_per_frame() const;
size_t pixels_per_frame() const;

30
include/aare/FilePtr.hpp Normal file
View File

@ -0,0 +1,30 @@
#pragma once
#include <cstdio>
#include <filesystem>
namespace aare {
/**
* \brief RAII wrapper for FILE pointer
*/
class FilePtr {
FILE *fp_{nullptr};
public:
FilePtr() = default;
FilePtr(const std::filesystem::path& fname, const std::string& mode);
FilePtr(const FilePtr &) = delete; // we don't want a copy
FilePtr &operator=(const FilePtr &) = delete; // since we handle a resource
FilePtr(FilePtr &&other);
FilePtr &operator=(FilePtr &&other);
FILE *get();
ssize_t tell();
void seek(ssize_t offset, int whence = SEEK_SET) {
if (fseek(fp_, offset, whence) != 0)
throw std::runtime_error("Error seeking in file");
}
std::string error_msg();
~FilePtr();
};
} // namespace aare

115
include/aare/Fit.hpp Normal file
View File

@ -0,0 +1,115 @@
#pragma once
#include <cmath>
#include <fmt/core.h>
#include <vector>
#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
/**
* @brief Estimate the initial parameters for a Gaussian fit
*/
std::array<double, 3> gaus_init_par(const NDView<double, 1> x, const NDView<double, 1> y);
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;
/**
* @brief Fit a 1D Gaussian to data.
* @param data data to fit
* @param x x values
*/
NDArray<double, 1> fit_gaus(NDView<double, 1> x, NDView<double, 1> y);
/**
* @brief Fit a 1D Gaussian to each pixel. Data layout [row, col, values]
* @param x x values
* @param y y values, layout [row, col, values]
* @param n_threads number of threads to use
*/
NDArray<double, 3> fit_gaus(NDView<double, 1> x, NDView<double, 3> y,
int n_threads = DEFAULT_NUM_THREADS);
/**
* @brief Fit a 1D Gaussian with error estimates
* @param x x values
* @param y y values, layout [row, col, values]
* @param y_err error in y, layout [row, col, values]
* @param par_out output parameters
* @param par_err_out output error parameters
*/
void fit_gaus(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
NDView<double, 1> par_out, NDView<double, 1> par_err_out,
double& chi2);
/**
* @brief Fit a 1D Gaussian to each pixel with error estimates. Data layout
* [row, col, values]
* @param x x values
* @param y y values, layout [row, col, values]
* @param y_err error in y, layout [row, col, values]
* @param par_out output parameters, layout [row, col, values]
* @param par_err_out output parameter errors, layout [row, col, values]
* @param n_threads number of threads to use
*/
void fit_gaus(NDView<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_err,
NDView<double, 3> par_out, NDView<double, 3> par_err_out, NDView<double, 2> chi2_out,
int n_threads = DEFAULT_NUM_THREADS
);
NDArray<double, 1> fit_pol1(NDView<double, 1> x, NDView<double, 1> y);
NDArray<double, 3> fit_pol1(NDView<double, 1> x, NDView<double, 3> y,
int n_threads = DEFAULT_NUM_THREADS);
void fit_pol1(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
NDView<double, 1> par_out, NDView<double, 1> par_err_out, double& chi2);
// TODO! not sure we need to offer the different version in C++
void fit_pol1(NDView<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_err,
NDView<double, 3> par_out, NDView<double, 3> par_err_out,NDView<double, 2> chi2_out,
int n_threads = DEFAULT_NUM_THREADS);
NDArray<double, 1> fit_scurve(NDView<double, 1> x, NDView<double, 1> y);
NDArray<double, 3> fit_scurve(NDView<double, 1> x, NDView<double, 3> y, int n_threads);
void fit_scurve(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
NDView<double, 1> par_out, NDView<double, 1> par_err_out, double& chi2);
void fit_scurve(NDView<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_err,
NDView<double, 3> par_out, NDView<double, 3> par_err_out, NDView<double, 2> chi2_out,
int n_threads);
NDArray<double, 1> fit_scurve2(NDView<double, 1> x, NDView<double, 1> y);
NDArray<double, 3> fit_scurve2(NDView<double, 1> x, NDView<double, 3> y, int n_threads);
void fit_scurve2(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
NDView<double, 1> par_out, NDView<double, 1> par_err_out, double& chi2);
void fit_scurve2(NDView<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_err,
NDView<double, 3> par_out, NDView<double, 3> par_err_out, NDView<double, 2> chi2_out,
int n_threads);
} // namespace aare

4
include/aare/Frame.hpp
View File

@ -107,8 +107,8 @@ class Frame {
* @return NDView<T, 2>
*/
template <typename T> NDView<T, 2> view() {
std::array<int64_t, 2> shape = {static_cast<int64_t>(m_rows),
static_cast<int64_t>(m_cols)};
std::array<ssize_t, 2> shape = {static_cast<ssize_t>(m_rows),
static_cast<ssize_t>(m_cols)};
T *data = reinterpret_cast<T *>(m_data);
return NDView<T, 2>(data, shape);
}

68
include/aare/GainMap.hpp Normal file
View File

@ -0,0 +1,68 @@
/************************************************
* @file GainMap.hpp
* @short function to apply gain map of image size to a vector of clusters -
*note stored gainmap is inverted for efficient aaplication to images
***********************************************/
#pragma once
#include "aare/Cluster.hpp"
#include "aare/ClusterVector.hpp"
#include "aare/NDArray.hpp"
#include "aare/NDView.hpp"
#include <memory>
namespace aare {
class InvertedGainMap {
public:
explicit InvertedGainMap(const NDArray<double, 2> &gain_map)
: m_gain_map(gain_map) {
for (auto &item : m_gain_map) {
item = 1.0 / item;
}
};
explicit InvertedGainMap(const NDView<double, 2> gain_map) {
m_gain_map = NDArray<double, 2>(gain_map);
for (auto &item : m_gain_map) {
item = 1.0 / item;
}
}
template <typename ClusterType,
typename = std::enable_if_t<is_cluster_v<ClusterType>>>
void apply_gain_map(ClusterVector<ClusterType> &clustervec) {
// in principle we need to know the size of the image for this lookup
size_t ClusterSizeX = clustervec.cluster_size_x();
size_t ClusterSizeY = clustervec.cluster_size_y();
using T = typename ClusterVector<ClusterType>::value_type;
int64_t index_cluster_center_x = ClusterSizeX / 2;
int64_t index_cluster_center_y = ClusterSizeY / 2;
for (size_t i = 0; i < clustervec.size(); i++) {
auto &cl = clustervec[i];
if (cl.x > 0 && cl.y > 0 && cl.x < m_gain_map.shape(1) - 1 &&
cl.y < m_gain_map.shape(0) - 1) {
for (size_t j = 0; j < ClusterSizeX * ClusterSizeY; j++) {
size_t x = cl.x + j % ClusterSizeX - index_cluster_center_x;
size_t y = cl.y + j / ClusterSizeX - index_cluster_center_y;
cl.data[j] = static_cast<T>(
static_cast<double>(cl.data[j]) *
m_gain_map(
y, x)); // cast after conversion to keep precision
}
} else {
// clear edge clusters
cl.data.fill(0);
}
}
}
private:
NDArray<double, 2> m_gain_map{};
};
} // end of namespace aare

130
include/aare/Interpolator.hpp Normal file
View File

@ -0,0 +1,130 @@
#pragma once
#include "aare/CalculateEta.hpp"
#include "aare/Cluster.hpp"
#include "aare/ClusterFile.hpp" //Cluster_3x3
#include "aare/ClusterVector.hpp"
#include "aare/NDArray.hpp"
#include "aare/NDView.hpp"
#include "aare/algorithm.hpp"
namespace aare {
struct Photon {
double x;
double y;
double energy;
};
class Interpolator {
NDArray<double, 3> m_ietax;
NDArray<double, 3> m_ietay;
NDArray<double, 1> m_etabinsx;
NDArray<double, 1> m_etabinsy;
NDArray<double, 1> m_energy_bins;
public:
Interpolator(NDView<double, 3> etacube, NDView<double, 1> xbins,
NDView<double, 1> ybins, NDView<double, 1> ebins);
NDArray<double, 3> get_ietax() { return m_ietax; }
NDArray<double, 3> get_ietay() { return m_ietay; }
template <typename ClusterType,
typename Eanble = std::enable_if_t<is_cluster_v<ClusterType>>>
std::vector<Photon> interpolate(const ClusterVector<ClusterType> &clusters);
};
// TODO: generalize to support any clustertype!!! otherwise add std::enable_if_t
// to only take Cluster2x2 and Cluster3x3
template <typename ClusterType, typename Enable>
std::vector<Photon>
Interpolator::interpolate(const ClusterVector<ClusterType> &clusters) {
std::vector<Photon> photons;
photons.reserve(clusters.size());
if (clusters.cluster_size_x() == 3 || clusters.cluster_size_y() == 3) {
for (const ClusterType &cluster : clusters) {
auto eta = calculate_eta2(cluster);
Photon photon;
photon.x = cluster.x;
photon.y = cluster.y;
photon.energy = static_cast<decltype(photon.energy)>(eta.sum);
// auto ie = nearest_index(m_energy_bins, photon.energy)-1;
// auto ix = nearest_index(m_etabinsx, eta.x)-1;
// auto iy = nearest_index(m_etabinsy, eta.y)-1;
// Finding the index of the last element that is smaller
// should work fine as long as we have many bins
auto ie = last_smaller(m_energy_bins, photon.energy);
auto ix = last_smaller(m_etabinsx, eta.x);
auto iy = last_smaller(m_etabinsy, eta.y);
// fmt::print("ex: {}, ix: {}, iy: {}\n", ie, ix, iy);
double dX, dY;
// cBottomLeft = 0,
// cBottomRight = 1,
// cTopLeft = 2,
// cTopRight = 3
switch (static_cast<corner>(eta.c)) {
case corner::cTopLeft:
dX = -1.;
dY = 0;
break;
case corner::cTopRight:;
dX = 0;
dY = 0;
break;
case corner::cBottomLeft:
dX = -1.;
dY = -1.;
break;
case corner::cBottomRight:
dX = 0.;
dY = -1.;
break;
}
photon.x += m_ietax(ix, iy, ie) * 2 + dX;
photon.y += m_ietay(ix, iy, ie) * 2 + dY;
photons.push_back(photon);
}
} else if (clusters.cluster_size_x() == 2 ||
clusters.cluster_size_y() == 2) {
for (const ClusterType &cluster : clusters) {
auto eta = calculate_eta2(cluster);
Photon photon;
photon.x = cluster.x;
photon.y = cluster.y;
photon.energy = static_cast<decltype(photon.energy)>(eta.sum);
// Now do some actual interpolation.
// Find which energy bin the cluster is in
// auto ie = nearest_index(m_energy_bins, photon.energy)-1;
// auto ix = nearest_index(m_etabinsx, eta.x)-1;
// auto iy = nearest_index(m_etabinsy, eta.y)-1;
// Finding the index of the last element that is smaller
// should work fine as long as we have many bins
auto ie = last_smaller(m_energy_bins, photon.energy);
auto ix = last_smaller(m_etabinsx, eta.x);
auto iy = last_smaller(m_etabinsy, eta.y);
photon.x += m_ietax(ix, iy, ie) *
2; // eta goes between 0 and 1 but we could move the hit
// anywhere in the 2x2
photon.y += m_ietay(ix, iy, ie) * 2;
photons.push_back(photon);
}
} else {
throw std::runtime_error(
"Only 3x3 and 2x2 clusters are supported for interpolation");
}
return photons;
}
} // namespace aare

106
include/aare/JungfrauDataFile.hpp Normal file
View File

@ -0,0 +1,106 @@
#pragma once
#include <cstdint>
#include <filesystem>
#include <vector>
#include "aare/FilePtr.hpp"
#include "aare/defs.hpp"
#include "aare/NDArray.hpp"
#include "aare/FileInterface.hpp"
namespace aare {
struct JungfrauDataHeader{
uint64_t framenum;
uint64_t bunchid;
};
class JungfrauDataFile : public FileInterface {
size_t m_rows{}; //!< number of rows in the image, from find_frame_size();
size_t m_cols{}; //!< number of columns in the image, from find_frame_size();
size_t m_bytes_per_frame{}; //!< number of bytes per frame excluding header
size_t m_total_frames{}; //!< total number of frames in the series of files
size_t m_offset{}; //!< file index of the first file, allow starting at non zero file
size_t m_current_file_index{}; //!< The index of the open file
size_t m_current_frame_index{}; //!< The index of the current frame (with reference to all files)
std::vector<size_t> m_last_frame_in_file{}; //!< Used for seeking to the correct file
std::filesystem::path m_path; //!< path to the files
std::string m_base_name; //!< base name used for formatting file names
FilePtr m_fp; //!< RAII wrapper for a FILE*
using pixel_type = uint16_t;
static constexpr size_t header_size = sizeof(JungfrauDataHeader);
static constexpr size_t n_digits_in_file_index = 6; //!< to format file names
public:
JungfrauDataFile(const std::filesystem::path &fname);
std::string base_name() const; //!< get the base name of the file (without path and extension)
size_t bytes_per_frame() override;
size_t pixels_per_frame() override;
size_t bytes_per_pixel() const;
size_t bitdepth() const override;
void seek(size_t frame_index) override; //!< seek to the given frame index (note not byte offset)
size_t tell() override; //!< get the frame index of the file pointer
size_t total_frames() const override;
size_t rows() const override;
size_t cols() const override;
std::array<ssize_t,2> shape() const;
size_t n_files() const; //!< get the number of files in the series.
// Extra functions needed for FileInterface
Frame read_frame() override;
Frame read_frame(size_t frame_number) override;
std::vector<Frame> read_n(size_t n_frames=0) override;
void read_into(std::byte *image_buf) override;
void read_into(std::byte *image_buf, size_t n_frames) override;
size_t frame_number(size_t frame_index) override;
DetectorType detector_type() const override;
/**
* @brief Read a single frame from the file into the given buffer.
* @param image_buf buffer to read the frame into. (Note the caller is responsible for allocating the buffer)
* @param header pointer to a JungfrauDataHeader or nullptr to skip header)
*/
void read_into(std::byte *image_buf, JungfrauDataHeader *header = nullptr);
/**
* @brief Read a multiple frames from the file into the given buffer.
* @param image_buf buffer to read the frame into. (Note the caller is responsible for allocating the buffer)
* @param n_frames number of frames to read
* @param header pointer to a JungfrauDataHeader or nullptr to skip header)
*/
void read_into(std::byte *image_buf, size_t n_frames, JungfrauDataHeader *header = nullptr);
/**
* @brief Read a single frame from the file into the given NDArray
* @param image NDArray to read the frame into.
*/
void read_into(NDArray<uint16_t>* image, JungfrauDataHeader* header = nullptr);
JungfrauDataHeader read_header();
std::filesystem::path current_file() const { return fpath(m_current_file_index+m_offset); }
private:
/**
* @brief Find the size of the frame in the file. (256x256, 256x1024, 512x1024)
* @param fname path to the file
* @throws std::runtime_error if the file is empty or the size cannot be determined
*/
void find_frame_size(const std::filesystem::path &fname);
void parse_fname(const std::filesystem::path &fname);
void scan_files();
void open_file(size_t file_index);
std::filesystem::path fpath(size_t frame_index) const;
};
} // namespace aare

124
include/aare/NDArray.hpp
View File

@ -22,10 +22,10 @@ TODO! Add expression templates for operators
namespace aare {
template <typename T, int64_t Ndim = 2>
template <typename T, ssize_t Ndim = 2>
class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
std::array<int64_t, Ndim> shape_;
std::array<int64_t, Ndim> strides_;
std::array<ssize_t, Ndim> shape_;
std::array<ssize_t, Ndim> strides_;
size_t size_{};
T *data_;
@ -42,7 +42,7 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
*
* @param shape shape of the new NDArray
*/
explicit NDArray(std::array<int64_t, Ndim> shape)
explicit NDArray(std::array<ssize_t, Ndim> shape)
: shape_(shape), strides_(c_strides<Ndim>(shape_)),
size_(std::accumulate(shape_.begin(), shape_.end(), 1,
std::multiplies<>())),
@ -55,7 +55,7 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
* @param shape shape of the new array
* @param value value to initialize the array with
*/
NDArray(std::array<int64_t, Ndim> shape, T value) : NDArray(shape) {
NDArray(std::array<ssize_t, Ndim> shape, T value) : NDArray(shape) {
this->operator=(value);
}
@ -69,6 +69,11 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
std::copy(v.begin(), v.end(), begin());
}
template<size_t Size>
NDArray(const std::array<T, Size>& arr) : NDArray<T,1>({Size}) {
std::copy(arr.begin(), arr.end(), begin());
}
// Move constructor
NDArray(NDArray &&other) noexcept
: shape_(other.shape_), strides_(c_strides<Ndim>(shape_)),
@ -87,7 +92,7 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
// Conversion operator from array expression to array
template <typename E>
NDArray(ArrayExpr<E, Ndim> &&expr) : NDArray(expr.shape()) {
for (int i = 0; i < size_; ++i) {
for (size_t i = 0; i < size_; ++i) {
data_[i] = expr[i];
}
}
@ -97,6 +102,9 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
auto begin() { return data_; }
auto end() { return data_ + size_; }
auto begin() const { return data_; }
auto end() const { return data_ + size_; }
using value_type = T;
NDArray &operator=(NDArray &&other) noexcept; // Move assign
@ -105,6 +113,20 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
NDArray &operator-=(const NDArray &other);
NDArray &operator*=(const NDArray &other);
//Write directly to the data array, or create a new one
template<size_t Size>
NDArray<T,1>& operator=(const std::array<T,Size> &other){
if(Size != size_){
delete[] data_;
size_ = Size;
data_ = new T[size_];
}
for (size_t i = 0; i < Size; ++i) {
data_[i] = other[i];
}
return *this;
}
// NDArray& operator/=(const NDArray& other);
template <typename V> NDArray &operator/=(const NDArray<V, Ndim> &other) {
@ -135,6 +157,11 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
NDArray &operator&=(const T & /*mask*/);
void sqrt() {
for (int i = 0; i < size_; ++i) {
data_[i] = std::sqrt(data_[i]);
@ -159,22 +186,22 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
}
// TODO! is int the right type for index?
T &operator()(int i) { return data_[i]; }
const T &operator()(int i) const { return data_[i]; }
T &operator()(ssize_t i) { return data_[i]; }
const T &operator()(ssize_t i) const { return data_[i]; }
T &operator[](int i) { return data_[i]; }
const T &operator[](int i) const { return data_[i]; }
T &operator[](ssize_t i) { return data_[i]; }
const T &operator[](ssize_t i) const { return data_[i]; }
T *data() { return data_; }
std::byte *buffer() { return reinterpret_cast<std::byte *>(data_); }
size_t size() const { return size_; }
ssize_t size() const { return static_cast<ssize_t>(size_); }
size_t total_bytes() const { return size_ * sizeof(T); }
std::array<int64_t, Ndim> shape() const noexcept { return shape_; }
int64_t shape(int64_t i) const noexcept { return shape_[i]; }
std::array<int64_t, Ndim> strides() const noexcept { return strides_; }
std::array<ssize_t, Ndim> shape() const noexcept { return shape_; }
ssize_t shape(ssize_t i) const noexcept { return shape_[i]; }
std::array<ssize_t, Ndim> strides() const noexcept { return strides_; }
size_t bitdepth() const noexcept { return sizeof(T) * 8; }
std::array<int64_t, Ndim> byte_strides() const noexcept {
std::array<ssize_t, Ndim> byte_strides() const noexcept {
auto byte_strides = strides_;
for (auto &val : byte_strides)
val *= sizeof(T);
@ -201,7 +228,7 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
};
// Move assign
template <typename T, int64_t Ndim>
template <typename T, ssize_t Ndim>
NDArray<T, Ndim> &
NDArray<T, Ndim>::operator=(NDArray<T, Ndim> &&other) noexcept {
if (this != &other) {
@ -215,7 +242,7 @@ NDArray<T, Ndim>::operator=(NDArray<T, Ndim> &&other) noexcept {
return *this;
}
template <typename T, int64_t Ndim>
template <typename T, ssize_t Ndim>
NDArray<T, Ndim> &NDArray<T, Ndim>::operator+=(const NDArray<T, Ndim> &other) {
// check shape
if (shape_ == other.shape_) {
@ -227,7 +254,7 @@ NDArray<T, Ndim> &NDArray<T, Ndim>::operator+=(const NDArray<T, Ndim> &other) {
throw(std::runtime_error("Shape of ImageDatas must match"));
}
template <typename T, int64_t Ndim>
template <typename T, ssize_t Ndim>
NDArray<T, Ndim> &NDArray<T, Ndim>::operator-=(const NDArray<T, Ndim> &other) {
// check shape
if (shape_ == other.shape_) {
@ -239,7 +266,7 @@ NDArray<T, Ndim> &NDArray<T, Ndim>::operator-=(const NDArray<T, Ndim> &other) {
throw(std::runtime_error("Shape of ImageDatas must match"));
}
template <typename T, int64_t Ndim>
template <typename T, ssize_t Ndim>
NDArray<T, Ndim> &NDArray<T, Ndim>::operator*=(const NDArray<T, Ndim> &other) {
// check shape
if (shape_ == other.shape_) {
@ -251,14 +278,14 @@ NDArray<T, Ndim> &NDArray<T, Ndim>::operator*=(const NDArray<T, Ndim> &other) {
throw(std::runtime_error("Shape of ImageDatas must match"));
}
template <typename T, int64_t Ndim>
template <typename T, ssize_t Ndim>
NDArray<T, Ndim> &NDArray<T, Ndim>::operator&=(const T &mask) {
for (auto it = begin(); it != end(); ++it)
*it &= mask;
return *this;
}
template <typename T, int64_t Ndim>
template <typename T, ssize_t Ndim>
NDArray<bool, Ndim> NDArray<T, Ndim>::operator>(const NDArray &other) {
if (shape_ == other.shape_) {
NDArray<bool, Ndim> result{shape_};
@ -270,7 +297,7 @@ NDArray<bool, Ndim> NDArray<T, Ndim>::operator>(const NDArray &other) {
throw(std::runtime_error("Shape of ImageDatas must match"));
}
template <typename T, int64_t Ndim>
template <typename T, ssize_t Ndim>
NDArray<T, Ndim> &NDArray<T, Ndim>::operator=(const NDArray<T, Ndim> &other) {
if (this != &other) {
delete[] data_;
@ -283,7 +310,7 @@ NDArray<T, Ndim> &NDArray<T, Ndim>::operator=(const NDArray<T, Ndim> &other) {
return *this;
}
template <typename T, int64_t Ndim>
template <typename T, ssize_t Ndim>
bool NDArray<T, Ndim>::operator==(const NDArray<T, Ndim> &other) const {
if (shape_ != other.shape_)
return false;
@ -295,80 +322,83 @@ bool NDArray<T, Ndim>::operator==(const NDArray<T, Ndim> &other) const {
return true;
}
template <typename T, int64_t Ndim>
template <typename T, ssize_t Ndim>
bool NDArray<T, Ndim>::operator!=(const NDArray<T, Ndim> &other) const {
return !((*this) == other);
}
template <typename T, int64_t Ndim>
template <typename T, ssize_t Ndim>
NDArray<T, Ndim> &NDArray<T, Ndim>::operator++() {
for (uint32_t i = 0; i < size_; ++i)
data_[i] += 1;
return *this;
}
template <typename T, int64_t Ndim>
template <typename T, ssize_t Ndim>
NDArray<T, Ndim> &NDArray<T, Ndim>::operator=(const T &value) {
std::fill_n(data_, size_, value);
return *this;
}
template <typename T, int64_t Ndim>
template <typename T, ssize_t Ndim>
NDArray<T, Ndim> &NDArray<T, Ndim>::operator+=(const T &value) {
for (uint32_t i = 0; i < size_; ++i)
data_[i] += value;
return *this;
}
template <typename T, int64_t Ndim>
template <typename T, ssize_t Ndim>
NDArray<T, Ndim> NDArray<T, Ndim>::operator+(const T &value) {
NDArray result = *this;
result += value;
return result;
}
template <typename T, int64_t Ndim>
template <typename T, ssize_t Ndim>
NDArray<T, Ndim> &NDArray<T, Ndim>::operator-=(const T &value) {
for (uint32_t i = 0; i < size_; ++i)
data_[i] -= value;
return *this;
}
template <typename T, int64_t Ndim>
template <typename T, ssize_t Ndim>
NDArray<T, Ndim> NDArray<T, Ndim>::operator-(const T &value) {
NDArray result = *this;
result -= value;
return result;
}
template <typename T, int64_t Ndim>
template <typename T, ssize_t Ndim>
NDArray<T, Ndim> &NDArray<T, Ndim>::operator/=(const T &value) {
for (uint32_t i = 0; i < size_; ++i)
data_[i] /= value;
return *this;
}
template <typename T, int64_t Ndim>
template <typename T, ssize_t Ndim>
NDArray<T, Ndim> NDArray<T, Ndim>::operator/(const T &value) {
NDArray result = *this;
result /= value;
return result;
}
template <typename T, int64_t Ndim>
template <typename T, ssize_t Ndim>
NDArray<T, Ndim> &NDArray<T, Ndim>::operator*=(const T &value) {
for (uint32_t i = 0; i < size_; ++i)
data_[i] *= value;
return *this;
}
template <typename T, int64_t Ndim>
template <typename T, ssize_t Ndim>
NDArray<T, Ndim> NDArray<T, Ndim>::operator*(const T &value) {
NDArray result = *this;
result *= value;
return result;
}
template <typename T, int64_t Ndim> void NDArray<T, Ndim>::Print() {
if (shape_[0] < 20 && shape_[1] < 20)
Print_all();
else
Print_some();
}
// template <typename T, ssize_t Ndim> void NDArray<T, Ndim>::Print() {
// if (shape_[0] < 20 && shape_[1] < 20)
// Print_all();
// else
// Print_some();
// }
template <typename T, int64_t Ndim>
template <typename T, ssize_t Ndim>
std::ostream &operator<<(std::ostream &os, const NDArray<T, Ndim> &arr) {
for (auto row = 0; row < arr.shape(0); ++row) {
for (auto col = 0; col < arr.shape(1); ++col) {
@ -380,7 +410,7 @@ std::ostream &operator<<(std::ostream &os, const NDArray<T, Ndim> &arr) {
return os;
}
template <typename T, int64_t Ndim> void NDArray<T, Ndim>::Print_all() {
template <typename T, ssize_t Ndim> void NDArray<T, Ndim>::Print_all() {
for (auto row = 0; row < shape_[0]; ++row) {
for (auto col = 0; col < shape_[1]; ++col) {
std::cout << std::setw(3);
@ -389,7 +419,7 @@ template <typename T, int64_t Ndim> void NDArray<T, Ndim>::Print_all() {
std::cout << "\n";
}
}
template <typename T, int64_t Ndim> void NDArray<T, Ndim>::Print_some() {
template <typename T, ssize_t Ndim> void NDArray<T, Ndim>::Print_some() {
for (auto row = 0; row < 5; ++row) {
for (auto col = 0; col < 5; ++col) {
std::cout << std::setw(7);
@ -399,7 +429,7 @@ template <typename T, int64_t Ndim> void NDArray<T, Ndim>::Print_some() {
}
}
template <typename T, int64_t Ndim>
template <typename T, ssize_t Ndim>
void save(NDArray<T, Ndim> &img, std::string &pathname) {
std::ofstream f;
f.open(pathname, std::ios::binary);
@ -407,9 +437,9 @@ void save(NDArray<T, Ndim> &img, std::string &pathname) {
f.close();
}
template <typename T, int64_t Ndim>
template <typename T, ssize_t Ndim>
NDArray<T, Ndim> load(const std::string &pathname,
std::array<int64_t, Ndim> shape) {
std::array<ssize_t, Ndim> shape) {
NDArray<T, Ndim> img{shape};
std::ifstream f;
f.open(pathname, std::ios::binary);
@ -418,4 +448,6 @@ NDArray<T, Ndim> load(const std::string &pathname,
return img;
}
} // namespace aare

60
include/aare/NDView.hpp
View File

@ -1,5 +1,5 @@
#pragma once
#include "aare/defs.hpp"
#include "aare/ArrayExpr.hpp"
#include <algorithm>
@ -14,10 +14,10 @@
#include <vector>
namespace aare {
template <int64_t Ndim> using Shape = std::array<int64_t, Ndim>;
template <ssize_t Ndim> using Shape = std::array<ssize_t, Ndim>;
// TODO! fix mismatch between signed and unsigned
template <int64_t Ndim> Shape<Ndim> make_shape(const std::vector<size_t> &shape) {
template <ssize_t Ndim> Shape<Ndim> make_shape(const std::vector<size_t> &shape) {
if (shape.size() != Ndim)
throw std::runtime_error("Shape size mismatch");
Shape<Ndim> arr;
@ -25,41 +25,41 @@ template <int64_t Ndim> Shape<Ndim> make_shape(const std::vector<size_t> &shape)
return arr;
}
template <int64_t Dim = 0, typename Strides> int64_t element_offset(const Strides & /*unused*/) { return 0; }
template <ssize_t Dim = 0, typename Strides> ssize_t element_offset(const Strides & /*unused*/) { return 0; }
template <int64_t Dim = 0, typename Strides, typename... Ix>
int64_t element_offset(const Strides &strides, int64_t i, Ix... index) {
template <ssize_t Dim = 0, typename Strides, typename... Ix>
ssize_t element_offset(const Strides &strides, ssize_t i, Ix... index) {
return i * strides[Dim] + element_offset<Dim + 1>(strides, index...);
}
template <int64_t Ndim> std::array<int64_t, Ndim> c_strides(const std::array<int64_t, Ndim> &shape) {
std::array<int64_t, Ndim> strides{};
template <ssize_t Ndim> std::array<ssize_t, Ndim> c_strides(const std::array<ssize_t, Ndim> &shape) {
std::array<ssize_t, Ndim> strides{};
std::fill(strides.begin(), strides.end(), 1);
for (int64_t i = Ndim - 1; i > 0; --i) {
for (ssize_t i = Ndim - 1; i > 0; --i) {
strides[i - 1] = strides[i] * shape[i];
}
return strides;
}
template <int64_t Ndim> std::array<int64_t, Ndim> make_array(const std::vector<int64_t> &vec) {
template <ssize_t Ndim> std::array<ssize_t, Ndim> make_array(const std::vector<ssize_t> &vec) {
assert(vec.size() == Ndim);
std::array<int64_t, Ndim> arr{};
std::array<ssize_t, Ndim> arr{};
std::copy_n(vec.begin(), Ndim, arr.begin());
return arr;
}
template <typename T, int64_t Ndim = 2> class NDView : public ArrayExpr<NDView<T, Ndim>, Ndim> {
template <typename T, ssize_t Ndim = 2> class NDView : public ArrayExpr<NDView<T, Ndim>, Ndim> {
public:
NDView() = default;
~NDView() = default;
NDView(const NDView &) = default;
NDView(NDView &&) = default;
NDView(T *buffer, std::array<int64_t, Ndim> shape)
NDView(T *buffer, std::array<ssize_t, Ndim> shape)
: buffer_(buffer), strides_(c_strides<Ndim>(shape)), shape_(shape),
size_(std::accumulate(std::begin(shape), std::end(shape), 1, std::multiplies<>())) {}
// NDView(T *buffer, const std::vector<int64_t> &shape)
// NDView(T *buffer, const std::vector<ssize_t> &shape)
// : buffer_(buffer), strides_(c_strides<Ndim>(make_array<Ndim>(shape))), shape_(make_array<Ndim>(shape)),
// size_(std::accumulate(std::begin(shape), std::end(shape), 1, std::multiplies<>())) {}
@ -71,16 +71,16 @@ template <typename T, int64_t Ndim = 2> class NDView : public ArrayExpr<NDView<T
return buffer_[element_offset(strides_, index...)];
}
size_t size() const { return size_; }
ssize_t size() const { return static_cast<ssize_t>(size_); }
size_t total_bytes() const { return size_ * sizeof(T); }
std::array<int64_t, Ndim> strides() const noexcept { return strides_; }
std::array<ssize_t, Ndim> strides() const noexcept { return strides_; }
T *begin() { return buffer_; }
T *end() { return buffer_ + size_; }
T const *begin() const { return buffer_; }
T const *end() const { return buffer_ + size_; }
T &operator()(int64_t i) const { return buffer_[i]; }
T &operator[](int64_t i) const { return buffer_[i]; }
T &operator()(ssize_t i) const { return buffer_[i]; }
T &operator[](ssize_t i) const { return buffer_[i]; }
bool operator==(const NDView &other) const {
if (size_ != other.size_)
@ -99,6 +99,15 @@ template <typename T, int64_t Ndim = 2> class NDView : public ArrayExpr<NDView<T
NDView &operator/=(const NDView &other) { return elemenwise(other, std::divides<T>()); }
template<size_t Size>
NDView& operator=(const std::array<T, Size> &arr) {
if(size() != static_cast<ssize_t>(arr.size()))
throw std::runtime_error(LOCATION + "Array and NDView size mismatch");
std::copy(arr.begin(), arr.end(), begin());
return *this;
}
NDView &operator=(const T val) {
for (auto it = begin(); it != end(); ++it)
*it = val;
@ -127,15 +136,15 @@ template <typename T, int64_t Ndim = 2> class NDView : public ArrayExpr<NDView<T
}
auto &shape() const { return shape_; }
auto shape(int64_t i) const { return shape_[i]; }
auto shape(ssize_t i) const { return shape_[i]; }
T *data() { return buffer_; }
void print_all() const;
private:
T *buffer_{nullptr};
std::array<int64_t, Ndim> strides_{};
std::array<int64_t, Ndim> shape_{};
std::array<ssize_t, Ndim> strides_{};
std::array<ssize_t, Ndim> shape_{};
uint64_t size_{};
template <class BinaryOperation> NDView &elemenwise(T val, BinaryOperation op) {
@ -151,7 +160,7 @@ template <typename T, int64_t Ndim = 2> class NDView : public ArrayExpr<NDView<T
return *this;
}
};
template <typename T, int64_t Ndim> void NDView<T, Ndim>::print_all() const {
template <typename T, ssize_t Ndim> void NDView<T, Ndim>::print_all() const {
for (auto row = 0; row < shape_[0]; ++row) {
for (auto col = 0; col < shape_[1]; ++col) {
std::cout << std::setw(3);
@ -162,7 +171,7 @@ template <typename T, int64_t Ndim> void NDView<T, Ndim>::print_all() const {
}
template <typename T, int64_t Ndim>
template <typename T, ssize_t Ndim>
std::ostream& operator <<(std::ostream& os, const NDView<T, Ndim>& arr){
for (auto row = 0; row < arr.shape(0); ++row) {
for (auto col = 0; col < arr.shape(1); ++col) {
@ -175,4 +184,9 @@ std::ostream& operator <<(std::ostream& os, const NDView<T, Ndim>& arr){
}
template <typename T>
NDView<T,1> make_view(std::vector<T>& vec){
return NDView<T,1>(vec.data(), {static_cast<ssize_t>(vec.size())});
}
} // namespace aare

2
include/aare/NumpyFile.hpp
View File

@ -69,7 +69,7 @@ class NumpyFile : public FileInterface {
*/
template <typename T, size_t NDim> NDArray<T, NDim> load() {
NDArray<T, NDim> arr(make_shape<NDim>(m_header.shape));
if (fseek(fp, static_cast<int64_t>(header_size), SEEK_SET)) {
if (fseek(fp, static_cast<long>(header_size), SEEK_SET)) {
throw std::runtime_error(LOCATION + "Error seeking to the start of the data");
}
size_t rc = fread(arr.data(), sizeof(T), arr.size(), fp);

138
include/aare/Pedestal.hpp
View File

@ -18,34 +18,48 @@ template <typename SUM_TYPE = double> class Pedestal {
uint32_t m_samples;
NDArray<uint32_t, 2> m_cur_samples;
//TODO! in case of int needs to be changed to uint64_t
NDArray<SUM_TYPE, 2> m_sum;
NDArray<SUM_TYPE, 2> m_sum2;
//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
NDArray<SUM_TYPE, 2> m_mean;
public:
Pedestal(uint32_t rows, uint32_t cols, uint32_t n_samples = 1000)
: m_rows(rows), m_cols(cols), m_samples(n_samples),
m_cur_samples(NDArray<uint32_t, 2>({rows, cols}, 0)),
m_sum(NDArray<SUM_TYPE, 2>({rows, cols})),
m_sum2(NDArray<SUM_TYPE, 2>({rows, cols})) {
m_sum2(NDArray<SUM_TYPE, 2>({rows, cols})),
m_mean(NDArray<SUM_TYPE, 2>({rows, cols})) {
assert(rows > 0 && cols > 0 && n_samples > 0);
m_sum = 0;
m_sum2 = 0;
m_mean = 0;
}
~Pedestal() = default;
NDArray<SUM_TYPE, 2> mean() {
NDArray<SUM_TYPE, 2> mean_array({m_rows, m_cols});
for (uint32_t i = 0; i < m_rows * m_cols; i++) {
mean_array(i / m_cols, i % m_cols) = mean(i / m_cols, i % m_cols);
}
return mean_array;
return m_mean;
}
SUM_TYPE mean(const uint32_t row, const uint32_t col) const {
return m_mean(row, col);
}
SUM_TYPE std(const uint32_t row, const uint32_t col) const {
return std::sqrt(variance(row, col));
}
SUM_TYPE variance(const uint32_t row, const uint32_t col) const {
if (m_cur_samples(row, col) == 0) {
return 0.0;
}
return m_sum(row, col) / m_cur_samples(row, col);
return m_sum2(row, col) / m_cur_samples(row, col) -
mean(row, col) * mean(row, col);
}
NDArray<SUM_TYPE, 2> variance() {
@ -57,13 +71,7 @@ template <typename SUM_TYPE = double> class Pedestal {
return variance_array;
}
SUM_TYPE variance(const uint32_t row, const uint32_t col) const {
if (m_cur_samples(row, col) == 0) {
return 0.0;
}
return m_sum2(row, col) / m_cur_samples(row, col) -
mean(row, col) * mean(row, col);
}
NDArray<SUM_TYPE, 2> std() {
NDArray<SUM_TYPE, 2> standard_deviation_array({m_rows, m_cols});
@ -75,14 +83,13 @@ template <typename SUM_TYPE = double> class Pedestal {
return standard_deviation_array;
}
SUM_TYPE std(const uint32_t row, const uint32_t col) const {
return std::sqrt(variance(row, col));
}
void clear() {
for (uint32_t i = 0; i < m_rows * m_cols; i++) {
clear(i / m_cols, i % m_cols);
}
m_sum = 0;
m_sum2 = 0;
m_cur_samples = 0;
m_mean = 0;
}
@ -91,28 +98,51 @@ template <typename SUM_TYPE = double> class Pedestal {
m_sum(row, col) = 0;
m_sum2(row, col) = 0;
m_cur_samples(row, col) = 0;
m_mean(row, col) = 0;
}
// frame level operations
template <typename T> void push(NDView<T, 2> frame) {
assert(frame.size() == m_rows * m_cols);
// TODO! move away from m_rows, m_cols
if (frame.shape() != std::array<int64_t, 2>{m_rows, m_cols}) {
if (frame.shape() != std::array<ssize_t, 2>{m_rows, m_cols}) {
throw std::runtime_error(
"Frame shape does not match pedestal shape");
}
for (uint32_t row = 0; row < m_rows; row++) {
for (uint32_t col = 0; col < m_cols; col++) {
for (size_t row = 0; row < m_rows; row++) {
for (size_t col = 0; col < m_cols; col++) {
push<T>(row, col, frame(row, col));
}
}
// // TODO: test the effect of #pragma omp parallel for
// for (uint32_t index = 0; index < m_rows * m_cols; index++) {
// push<T>(index / m_cols, index % m_cols, frame(index));
// }
}
/**
* Push but don't update the cached mean. Speeds up the process
* when initializing the pedestal.
*
*/
template <typename T> void push_no_update(NDView<T, 2> frame) {
assert(frame.size() == m_rows * m_cols);
// TODO! move away from m_rows, m_cols
if (frame.shape() != std::array<ssize_t, 2>{m_rows, m_cols}) {
throw std::runtime_error(
"Frame shape does not match pedestal shape");
}
for (size_t row = 0; row < m_rows; row++) {
for (size_t col = 0; col < m_cols; col++) {
push_no_update<T>(row, col, frame(row, col));
}
}
}
template <typename T> void push(Frame &frame) {
assert(frame.rows() == static_cast<size_t>(m_rows) &&
frame.cols() == static_cast<size_t>(m_cols));
@ -132,18 +162,48 @@ template <typename SUM_TYPE = double> class Pedestal {
template <typename T>
void push(const uint32_t row, const uint32_t col, const T val_) {
SUM_TYPE val = static_cast<SUM_TYPE>(val_);
const uint32_t idx = index(row, col);
if (m_cur_samples(idx) < m_samples) {
m_sum(idx) += val;
m_sum2(idx) += val * val;
m_cur_samples(idx)++;
if (m_cur_samples(row, col) < m_samples) {
m_sum(row, col) += val;
m_sum2(row, col) += val * val;
m_cur_samples(row, col)++;
} else {
m_sum(idx) += val - m_sum(idx) / m_cur_samples(idx);
m_sum2(idx) += val * val - m_sum2(idx) / m_cur_samples(idx);
m_sum(row, col) += val - m_sum(row, col) / m_samples;
m_sum2(row, col) += val * val - m_sum2(row, col) / m_samples;
}
//Since we just did a push we know that m_cur_samples(row, col) is at least 1
m_mean(row, col) = m_sum(row, col) / m_cur_samples(row, col);
}
template <typename T>
void push_no_update(const uint32_t row, const uint32_t col, const T val_) {
SUM_TYPE val = static_cast<SUM_TYPE>(val_);
if (m_cur_samples(row, col) < m_samples) {
m_sum(row, col) += val;
m_sum2(row, col) += val * val;
m_cur_samples(row, col)++;
} else {
m_sum(row, col) += val - m_sum(row, col) / m_cur_samples(row, col);
m_sum2(row, col) += val * val - m_sum2(row, col) / m_cur_samples(row, col);
}
}
uint32_t index(const uint32_t row, const uint32_t col) const {
return row * m_cols + col;
};
/**
* @brief Update the mean of the pedestal. This is used after having done
* push_no_update. It is not necessary to call this function after push.
*/
void update_mean(){
m_mean = m_sum / m_cur_samples;
}
template<typename T>
void push_fast(const uint32_t row, const uint32_t col, const T val_){
//Assume we reached the steady state where all pixels have
//m_samples samples
SUM_TYPE val = static_cast<SUM_TYPE>(val_);
m_sum(row, col) += val - m_sum(row, col) / m_samples;
m_sum2(row, col) += val * val - m_sum2(row, col) / m_samples;
m_mean(row, col) = m_sum(row, col) / m_samples;
}
};
} // namespace aare

203
include/aare/ProducerConsumerQueue.hpp Normal file
View File

@ -0,0 +1,203 @@
/*
* Copyright (c) Meta Platforms, Inc. and affiliates.
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
// @author Bo Hu (bhu@fb.com)
// @author Jordan DeLong (delong.j@fb.com)
// Changes made by PSD Detector Group:
// Copied: Line 34 constexpr std::size_t hardware_destructive_interference_size = 128; from folly/lang/Align.h
// Changed extension to .hpp
// Changed namespace to aare
#pragma once
#include <atomic>
#include <cassert>
#include <cstdlib>
#include <memory>
#include <stdexcept>
#include <type_traits>
#include <utility>
constexpr std::size_t hardware_destructive_interference_size = 128;
namespace aare {
/*
* ProducerConsumerQueue is a one producer and one consumer queue
* without locks.
*/
template <class T> struct ProducerConsumerQueue {
typedef T value_type;
ProducerConsumerQueue(const ProducerConsumerQueue &) = delete;
ProducerConsumerQueue &operator=(const ProducerConsumerQueue &) = delete;
ProducerConsumerQueue(ProducerConsumerQueue &&other){
size_ = other.size_;
records_ = other.records_;
other.records_ = nullptr;
readIndex_ = other.readIndex_.load(std::memory_order_acquire);
writeIndex_ = other.writeIndex_.load(std::memory_order_acquire);
}
ProducerConsumerQueue &operator=(ProducerConsumerQueue &&other){
size_ = other.size_;
records_ = other.records_;
other.records_ = nullptr;
readIndex_ = other.readIndex_.load(std::memory_order_acquire);
writeIndex_ = other.writeIndex_.load(std::memory_order_acquire);
return *this;
}
ProducerConsumerQueue():ProducerConsumerQueue(2){};
// size must be >= 2.
//
// Also, note that the number of usable slots in the queue at any
// given time is actually (size-1), so if you start with an empty queue,
// isFull() will return true after size-1 insertions.
explicit ProducerConsumerQueue(uint32_t size)
: size_(size), records_(static_cast<T *>(std::malloc(sizeof(T) * size))), readIndex_(0), writeIndex_(0) {
assert(size >= 2);
if (!records_) {
throw std::bad_alloc();
}
}
~ProducerConsumerQueue() {
// We need to destruct anything that may still exist in our queue.
// (No real synchronization needed at destructor time: only one
// thread can be doing this.)
if (!std::is_trivially_destructible<T>::value) {
size_t readIndex = readIndex_;
size_t endIndex = writeIndex_;
while (readIndex != endIndex) {
records_[readIndex].~T();
if (++readIndex == size_) {
readIndex = 0;
}
}
}
std::free(records_);
}
template <class... Args> bool write(Args &&...recordArgs) {
auto const currentWrite = writeIndex_.load(std::memory_order_relaxed);
auto nextRecord = currentWrite + 1;
if (nextRecord == size_) {
nextRecord = 0;
}
if (nextRecord != readIndex_.load(std::memory_order_acquire)) {
new (&records_[currentWrite]) T(std::forward<Args>(recordArgs)...);
writeIndex_.store(nextRecord, std::memory_order_release);
return true;
}
// queue is full
return false;
}
// move (or copy) the value at the front of the queue to given variable
bool read(T &record) {
auto const currentRead = readIndex_.load(std::memory_order_relaxed);
if (currentRead == writeIndex_.load(std::memory_order_acquire)) {
// queue is empty
return false;
}
auto nextRecord = currentRead + 1;
if (nextRecord == size_) {
nextRecord = 0;
}
record = std::move(records_[currentRead]);
records_[currentRead].~T();
readIndex_.store(nextRecord, std::memory_order_release);
return true;
}
// pointer to the value at the front of the queue (for use in-place) or
// nullptr if empty.
T *frontPtr() {
auto const currentRead = readIndex_.load(std::memory_order_relaxed);
if (currentRead == writeIndex_.load(std::memory_order_acquire)) {
// queue is empty
return nullptr;
}
return &records_[currentRead];
}
// queue must not be empty
void popFront() {
auto const currentRead = readIndex_.load(std::memory_order_relaxed);
assert(currentRead != writeIndex_.load(std::memory_order_acquire));
auto nextRecord = currentRead + 1;
if (nextRecord == size_) {
nextRecord = 0;
}
records_[currentRead].~T();
readIndex_.store(nextRecord, std::memory_order_release);
}
bool isEmpty() const {
return readIndex_.load(std::memory_order_acquire) == writeIndex_.load(std::memory_order_acquire);
}
bool isFull() const {
auto nextRecord = writeIndex_.load(std::memory_order_acquire) + 1;
if (nextRecord == size_) {
nextRecord = 0;
}
if (nextRecord != readIndex_.load(std::memory_order_acquire)) {
return false;
}
// queue is full
return true;
}
// * If called by consumer, then true size may be more (because producer may
// be adding items concurrently).
// * If called by producer, then true size may be less (because consumer may
// be removing items concurrently).
// * It is undefined to call this from any other thread.
size_t sizeGuess() const {
int ret = writeIndex_.load(std::memory_order_acquire) - readIndex_.load(std::memory_order_acquire);
if (ret < 0) {
ret += size_;
}
return ret;
}
// maximum number of items in the queue.
size_t capacity() const { return size_ - 1; }
private:
using AtomicIndex = std::atomic<unsigned int>;
char pad0_[hardware_destructive_interference_size];
// const uint32_t size_;
uint32_t size_;
// T *const records_;
T* records_;
alignas(hardware_destructive_interference_size) AtomicIndex readIndex_;
alignas(hardware_destructive_interference_size) AtomicIndex writeIndex_;
char pad1_[hardware_destructive_interference_size - sizeof(AtomicIndex)];
};
} // namespace aare

21
include/aare/RawFile.hpp
View File

@ -30,19 +30,12 @@ struct ModuleConfig {
* Consider using that unless you need raw file specific functionality.
*/
class RawFile : public FileInterface {
size_t n_subfiles{}; //f0,f1...fn
size_t n_subfile_parts{}; // d0,d1...dn
//TODO! move to vector of SubFile instead of pointers
std::vector<std::vector<RawSubFile *>> subfiles; //subfiles[f0,f1...fn][d0,d1...dn]
std::vector<xy> positions;
std::vector<ModuleGeometry> m_module_pixel_0;
std::vector<std::unique_ptr<RawSubFile>> m_subfiles;
ModuleConfig cfg{0, 0};
RawMasterFile m_master;
size_t m_current_frame{};
size_t m_rows{};
size_t m_cols{};
size_t m_current_subfile{};
DetectorGeometry m_geometry;
public:
/**
@ -52,7 +45,7 @@ class RawFile : public FileInterface {
*/
RawFile(const std::filesystem::path &fname, const std::string &mode = "r");
virtual ~RawFile() override;
virtual ~RawFile() override = default;
Frame read_frame() override;
Frame read_frame(size_t frame_number) override;
@ -76,7 +69,7 @@ class RawFile : public FileInterface {
size_t cols() const override;
size_t bitdepth() const override;
xy geometry();
size_t n_mod() const;
size_t n_modules() const;
RawMasterFile master() const;
@ -111,11 +104,9 @@ class RawFile : public FileInterface {
*/
static DetectorHeader read_header(const std::filesystem::path &fname);
void update_geometry_with_roi();
int find_number_of_subfiles();
void open_subfiles();
void find_geometry();
};
} // namespace aare

12
include/aare/RawMasterFile.hpp
View File

@ -62,17 +62,6 @@ class ScanParameters {
};
struct ROI{
int64_t xmin{};
int64_t xmax{};
int64_t ymin{};
int64_t ymax{};
int64_t height() const { return ymax - ymin; }
int64_t width() const { return xmax - xmin; }
};
/**
* @brief Class for parsing a master file either in our .json format or the old
* .raw format
@ -132,6 +121,7 @@ class RawMasterFile {
size_t total_frames_expected() const;
xy geometry() const;
size_t n_modules() const;
std::optional<size_t> analog_samples() const;
std::optional<size_t> digital_samples() const;

26
include/aare/RawSubFile.hpp
View File

@ -18,11 +18,20 @@ class RawSubFile {
std::ifstream m_file;
DetectorType m_detector_type;
size_t m_bitdepth;
std::filesystem::path m_fname;
std::filesystem::path m_path; //!< path to the subfile
std::string m_base_name; //!< base name used for formatting file names
size_t m_offset{}; //!< file index of the first file, allow starting at non zero file
size_t m_total_frames{}; //!< total number of frames in the series of files
size_t m_rows{};
size_t m_cols{};
size_t m_bytes_per_frame{};
size_t n_frames{};
int m_module_index{};
size_t m_current_file_index{}; //!< The index of the open file
size_t m_current_frame_index{}; //!< The index of the current frame (with reference to all files)
std::vector<size_t> m_last_frame_in_file{}; //!< Used for seeking to the correct file
uint32_t m_pos_row{};
uint32_t m_pos_col{};
@ -53,6 +62,7 @@ class RawSubFile {
size_t tell();
void read_into(std::byte *image_buf, DetectorHeader *header = nullptr);
void read_into(std::byte *image_buf, size_t n_frames, DetectorHeader *header= nullptr);
void get_part(std::byte *buffer, size_t frame_index);
void read_header(DetectorHeader *header);
@ -64,8 +74,18 @@ class RawSubFile {
size_t bytes_per_frame() const { return m_bytes_per_frame; }
size_t pixels_per_frame() const { return m_rows * m_cols; }
size_t bytes_per_pixel() const { return m_bitdepth / 8; }
size_t bytes_per_pixel() const { return m_bitdepth / bits_per_byte; }
size_t frames_in_file() const { return m_total_frames; }
private:
template <typename T>
void read_with_map(std::byte *image_buf);
void parse_fname(const std::filesystem::path &fname);
void scan_files();
void open_file(size_t file_index);
std::filesystem::path fpath(size_t file_index) const;
};

8
include/aare/VarClusterFinder.hpp
View File

@ -7,7 +7,7 @@
#include "aare/NDArray.hpp"
const int MAX_CLUSTER_SIZE = 200;
const int MAX_CLUSTER_SIZE = 50;
namespace aare {
template <typename T> class VarClusterFinder {
@ -28,7 +28,7 @@ template <typename T> class VarClusterFinder {
};
private:
const std::array<int64_t, 2> shape_;
const std::array<ssize_t, 2> shape_;
NDView<T, 2> original_;
NDArray<int, 2> labeled_;
NDArray<int, 2> peripheral_labeled_;
@ -226,7 +226,7 @@ template <typename T> void VarClusterFinder<T>::single_pass(NDView<T, 2> img) {
template <typename T> void VarClusterFinder<T>::first_pass() {
for (size_t i = 0; i < original_.size(); ++i) {
for (ssize_t i = 0; i < original_.size(); ++i) {
if (use_noise_map)
threshold_ = 5 * noiseMap(i);
binary_(i) = (original_(i) > threshold_);
@ -250,7 +250,7 @@ template <typename T> void VarClusterFinder<T>::first_pass() {
template <typename T> void VarClusterFinder<T>::second_pass() {
for (size_t i = 0; i != labeled_.size(); ++i) {
for (ssize_t i = 0; i != labeled_.size(); ++i) {
auto cl = labeled_(i);
if (cl != 0) {
auto it = child.find(cl);

122
include/aare/algorithm.hpp Normal file
View File

@ -0,0 +1,122 @@
#pragma once
#include <algorithm>
#include <array>
#include <vector>
#include <aare/NDArray.hpp>
namespace aare {
/**
* @brief Index of the last element that is smaller than val.
* Requires a sorted array. Uses >= for ordering. If all elements
* are smaller it returns the last element and if all elements are
* larger it returns the first element.
* @param first iterator to the first element
* @param last iterator to the last element
* @param val value to compare
* @return index of the last element that is smaller than val
*
*/
template <typename T>
size_t last_smaller(const T* first, const T* last, T val) {
for (auto iter = first+1; iter != last; ++iter) {
if (*iter >= val) {
return std::distance(first, iter-1);
}
}
return std::distance(first, last-1);
}
template <typename T>
size_t last_smaller(const NDArray<T, 1>& arr, T val) {
return last_smaller(arr.begin(), arr.end(), val);
}
template <typename T>
size_t last_smaller(const std::vector<T>& vec, T val) {
return last_smaller(vec.data(), vec.data()+vec.size(), val);
}
/**
* @brief Index of the first element that is larger than val.
* Requires a sorted array. Uses > for ordering. If all elements
* are larger it returns the first element and if all elements are
* smaller it returns the last element.
* @param first iterator to the first element
* @param last iterator to the last element
* @param val value to compare
* @return index of the first element that is larger than val
*/
template <typename T>
size_t first_larger(const T* first, const T* last, T val) {
for (auto iter = first; iter != last; ++iter) {
if (*iter > val) {
return std::distance(first, iter);
}
}
return std::distance(first, last-1);
}
template <typename T>
size_t first_larger(const NDArray<T, 1>& arr, T val) {
return first_larger(arr.begin(), arr.end(), val);
}
template <typename T>
size_t first_larger(const std::vector<T>& vec, T val) {
return first_larger(vec.data(), vec.data()+vec.size(), val);
}
/**
* @brief Index of the nearest element to val.
* Requires a sorted array. If there is no difference it takes the first element.
* @param first iterator to the first element
* @param last iterator to the last element
* @param val value to compare
* @return index of the nearest element
*/
template <typename T>
size_t nearest_index(const T* first, const T* last, T val) {
auto iter = std::min_element(first, last,
[val](T a, T b) {
return std::abs(a - val) < std::abs(b - val);
});
return std::distance(first, iter);
}
template <typename T>
size_t nearest_index(const NDArray<T, 1>& arr, T val) {
return nearest_index(arr.begin(), arr.end(), val);
}
template <typename T>
size_t nearest_index(const std::vector<T>& vec, T val) {
return nearest_index(vec.data(), vec.data()+vec.size(), val);
}
template <typename T, size_t N>
size_t nearest_index(const std::array<T,N>& arr, T val) {
return nearest_index(arr.data(), arr.data()+arr.size(), val);
}
template <typename T>
std::vector<T> cumsum(const std::vector<T>& vec) {
std::vector<T> result(vec.size());
std::partial_sum(vec.begin(), vec.end(), result.begin());
return result;
}
template <typename Container> bool all_equal(const Container &c) {
if (!c.empty() &&
std::all_of(begin(c), end(c),
[c](const typename Container::value_type &element) {
return element == c.front();
}))
return true;
return false;
}
} // namespace aare

26
include/aare/decode.hpp Normal file
View File

@ -0,0 +1,26 @@
#pragma once
#include <cstdint>
#include <vector>
#include <aare/NDView.hpp>
namespace aare {
uint16_t adc_sar_05_decode64to16(uint64_t input);
uint16_t adc_sar_04_decode64to16(uint64_t input);
void adc_sar_05_decode64to16(NDView<uint64_t, 2> input, NDView<uint16_t,2> output);
void adc_sar_04_decode64to16(NDView<uint64_t, 2> input, NDView<uint16_t,2> output);
/**
* @brief Apply custom weights to a 16-bit input value. Will sum up weights[i]**i
* for each bit i that is set in the input value.
* @throws std::out_of_range if weights.size() < 16
* @param input 16-bit input value
* @param weights vector of weights, size must be less than or equal to 16
*/
double apply_custom_weights(uint16_t input, const NDView<double, 1> weights);
void apply_custom_weights(NDView<uint16_t, 1> input, NDView<double, 1> output, const NDView<double, 1> weights);
} // namespace aare

47
include/aare/defs.hpp
View File

@ -1,11 +1,9 @@
#pragma once
#include "aare/Dtype.hpp"
// #include "aare/utils/logger.hpp"
#include <array>
#include <stdexcept>
#include <cassert>
#include <cstdint>
#include <cstring>
@ -38,16 +36,19 @@
namespace aare {
inline constexpr size_t bits_per_byte = 8;
void assert_failed(const std::string &msg);
class DynamicCluster {
public:
int cluster_sizeX;
int cluster_sizeY;
int16_t x;
int16_t y;
Dtype dt;
Dtype dt; // 4 bytes
private:
std::byte *m_data;
@ -179,15 +180,49 @@ template <typename T> struct t_xy {
using xy = t_xy<uint32_t>;
/**
* @brief Class to hold the geometry of a module. Where pixel 0 is located and the size of the module
*/
struct ModuleGeometry{
int x{};
int y{};
int origin_x{};
int origin_y{};
int height{};
int width{};
int row_index{};
int col_index{};
};
/**
* @brief Class to hold the geometry of a detector. Number of modules, their size and where pixel 0
* for each module is located
*/
struct DetectorGeometry{
int modules_x{};
int modules_y{};
int pixels_x{};
int pixels_y{};
int module_gap_row{};
int module_gap_col{};
std::vector<ModuleGeometry> module_pixel_0;
auto size() const { return module_pixel_0.size(); }
};
using dynamic_shape = std::vector<int64_t>;
struct ROI{
ssize_t xmin{};
ssize_t xmax{};
ssize_t ymin{};
ssize_t ymax{};
ssize_t height() const { return ymax - ymin; }
ssize_t width() const { return xmax - xmin; }
bool contains(ssize_t x, ssize_t y) const {
return x >= xmin && x < xmax && y >= ymin && y < ymax;
}
};
using dynamic_shape = std::vector<ssize_t>;
//TODO! Can we uniform enums between the libraries?

16
include/aare/geo_helpers.hpp Normal file
View File

@ -0,0 +1,16 @@
#pragma once
#include "aare/defs.hpp"
#include "aare/RawMasterFile.hpp" //ROI refactor away
namespace aare{
/**
* @brief Update the detector geometry given a region of interest
*
* @param geo
* @param roi
* @return DetectorGeometry
*/
DetectorGeometry update_geometry_with_roi(DetectorGeometry geo, ROI roi);
} // namespace aare

139
include/aare/logger.hpp Normal file
View File

@ -0,0 +1,139 @@
#pragma once
/*Utility to log to console*/
#include <iostream>
#include <sstream>
#include <sys/time.h>
namespace aare {
#define RED "\x1b[31m"
#define GREEN "\x1b[32m"
#define YELLOW "\x1b[33m"
#define BLUE "\x1b[34m"
#define MAGENTA "\x1b[35m"
#define CYAN "\x1b[36m"
#define GRAY "\x1b[37m"
#define DARKGRAY "\x1b[30m"
#define BG_BLACK "\x1b[48;5;232m"
#define BG_RED "\x1b[41m"
#define BG_GREEN "\x1b[42m"
#define BG_YELLOW "\x1b[43m"
#define BG_BLUE "\x1b[44m"
#define BG_MAGENTA "\x1b[45m"
#define BG_CYAN "\x1b[46m"
#define RESET "\x1b[0m"
#define BOLD "\x1b[1m"
enum TLogLevel {
logERROR,
logWARNING,
logINFOBLUE,
logINFOGREEN,
logINFORED,
logINFOCYAN,
logINFOMAGENTA,
logINFO,
logDEBUG,
logDEBUG1,
logDEBUG2,
logDEBUG3,
logDEBUG4,
logDEBUG5
};
// Compiler should optimize away anything below this value
#ifndef AARE_LOG_LEVEL
#define AARE_LOG_LEVEL "LOG LEVEL NOT SET IN CMAKE" //This is configured in the main CMakeLists.txt
#endif
#define __AT__ \
std::string(__FILE__) + std::string("::") + std::string(__func__) + \
std::string("(): ")
#define __SHORT_FORM_OF_FILE__ \
(strrchr(__FILE__, '/') ? strrchr(__FILE__, '/') + 1 : __FILE__)
#define __SHORT_AT__ \
std::string(__SHORT_FORM_OF_FILE__) + std::string("::") + \
std::string(__func__) + std::string("(): ")
class Logger {
std::ostringstream os;
TLogLevel m_level = AARE_LOG_LEVEL;
public:
Logger() = default;
explicit Logger(TLogLevel level) : m_level(level){};
~Logger() {
// output in the destructor to allow for << syntax
os << RESET << '\n';
std::clog << os.str() << std::flush; // Single write
}
static TLogLevel &ReportingLevel() { // singelton eeh TODO! Do we need a runtime option?
static TLogLevel reportingLevel = logDEBUG5;
return reportingLevel;
}
// Danger this buffer need as many elements as TLogLevel
static const char *Color(TLogLevel level) noexcept {
static const char *const colors[] = {
RED BOLD, YELLOW BOLD, BLUE, GREEN, RED, CYAN, MAGENTA,
RESET, RESET, RESET, RESET, RESET, RESET, RESET};
// out of bounds
if (level < 0 || level >= sizeof(colors) / sizeof(colors[0])) {
return RESET;
}
return colors[level];
}
// Danger this buffer need as many elements as TLogLevel
static std::string ToString(TLogLevel level) {
static const char *const buffer[] = {
"ERROR", "WARNING", "INFO", "INFO", "INFO",
"INFO", "INFO", "INFO", "DEBUG", "DEBUG1",
"DEBUG2", "DEBUG3", "DEBUG4", "DEBUG5"};
// out of bounds
if (level < 0 || level >= sizeof(buffer) / sizeof(buffer[0])) {
return "UNKNOWN";
}
return buffer[level];
}
std::ostringstream &Get() {
os << Color(m_level) << "- " << Timestamp() << " " << ToString(m_level)
<< ": ";
return os;
}
static std::string Timestamp() {
constexpr size_t buffer_len = 12;
char buffer[buffer_len];
time_t t;
::time(&t);
tm r;
strftime(buffer, buffer_len, "%X", localtime_r(&t, &r));
buffer[buffer_len - 1] = '\0';
struct timeval tv;
gettimeofday(&tv, nullptr);
constexpr size_t result_len = 100;
char result[result_len];
snprintf(result, result_len, "%s.%03ld", buffer,
static_cast<long>(tv.tv_usec) / 1000);
result[result_len - 1] = '\0';
return result;
}
};
// TODO! Do we need to keep the runtime option?
#define LOG(level) \
if (level > AARE_LOG_LEVEL) \
; \
else if (level > aare::Logger::ReportingLevel()) \
; \
else \
aare::Logger(level).Get()
} // namespace aare

12
include/aare/utils/ifstream_helpers.hpp Normal file
View File

@ -0,0 +1,12 @@
#pragma once
#include <fstream>
#include <string>
namespace aare {
/**
* @brief Get the error message from an ifstream object
*/
std::string ifstream_error_msg(std::ifstream &ifs);
} // namespace aare

18
include/aare/utils/par.hpp Normal file
View File

@ -0,0 +1,18 @@
#include <thread>
#include <vector>
#include <utility>
namespace aare {
template<typename F>
void RunInParallel(F func, const std::vector<std::pair<int, int>>& tasks) {
// auto tasks = split_task(0, y.shape(0), n_threads);
std::vector<std::thread> threads;
for (auto &task : tasks) {
threads.push_back(std::thread(func, task.first, task.second));
}
for (auto &thread : threads) {
thread.join();
}
}
} // namespace aare

8
include/aare/utils/task.hpp Normal file
View File

@ -0,0 +1,8 @@
#include <utility>
#include <vector>
namespace aare {
std::vector<std::pair<int, int>> split_task(int first, int last, int n_threads);
} // namespace aare

18
patches/libzmq_cmake_version.patch Normal file
View File

@ -0,0 +1,18 @@
diff --git a/CMakeLists.txt b/CMakeLists.txt
index dd3d8eb9..c0187747 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,11 +1,8 @@
# CMake build script for ZeroMQ
project(ZeroMQ)
-if(${CMAKE_SYSTEM_NAME} STREQUAL Darwin)
- cmake_minimum_required(VERSION 3.0.2)
-else()
- cmake_minimum_required(VERSION 2.8.12)
-endif()
+cmake_minimum_required(VERSION 3.15)
+message(STATUS "Patched cmake version")
include(CheckIncludeFiles)
include(CheckCCompilerFlag)

13
patches/lmfit.patch Normal file
View File

@ -0,0 +1,13 @@
diff --git a/lib/CMakeLists.txt b/lib/CMakeLists.txt
index 4efb7ed..6533660 100644
--- a/lib/CMakeLists.txt
+++ b/lib/CMakeLists.txt
@@ -11,7 +11,7 @@ target_compile_definitions(${lib} PRIVATE "LMFIT_EXPORT") # for Windows DLL expo
target_include_directories(${lib}
PUBLIC
- $<BUILD_INTERFACE:${CMAKE_SOURCE_DIR}/>
+ $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/>
$<INSTALL_INTERFACE:include/>
)

34
pyproject.toml
View File

@ -1,15 +1,41 @@
[tool.scikit-build.metadata.version]
provider = "scikit_build_core.metadata.regex"
input = "VERSION"
regex = '^(?P<version>\d+(?:\.\d+)*(?:[\.\+\w]+)?)$'
result = "{version}"
[build-system]
requires = ["scikit-build-core>=0.10", "pybind11", "numpy"]
build-backend = "scikit_build_core.build"
[project]
name = "aare"
version = "2024.11.28.dev0"
dynamic = ["version"]
requires-python = ">=3.11"
dependencies = [
"numpy",
"matplotlib",
]
[tool.cibuildwheel]
build = "cp{311,312,313}-manylinux_x86_64"
[tool.scikit-build]
cmake.verbose = true
build.verbose = true
cmake.build-type = "Release"
install.components = ["python"]
[tool.scikit-build.cmake.define]
AARE_PYTHON_BINDINGS = "ON"
AARE_SYSTEM_LIBRARIES = "ON"
AARE_INSTALL_PYTHONEXT = "ON"
AARE_INSTALL_PYTHONEXT = "ON"
[tool.pytest.ini_options]
markers = [
"files: marks tests that need additional data (deselect with '-m \"not files\"')",
]

46
python/CMakeLists.txt
View File

@ -1,12 +1,13 @@
find_package (Python 3.10 COMPONENTS Interpreter Development REQUIRED)
find_package (Python 3.10 COMPONENTS Interpreter Development.Module REQUIRED)
set(PYBIND11_FINDPYTHON ON) # Needed for RH8
# Download or find pybind11 depending on configuration
if(AARE_FETCH_PYBIND11)
FetchContent_Declare(
pybind11
GIT_REPOSITORY https://github.com/pybind/pybind11
GIT_TAG v2.13.0
GIT_TAG v2.13.6
)
FetchContent_MakeAvailable(pybind11)
else()
@ -28,6 +29,10 @@ target_link_libraries(_aare PRIVATE aare_core aare_compiler_flags)
set( PYTHON_FILES
aare/__init__.py
aare/CtbRawFile.py
aare/ClusterFinder.py
aare/ClusterVector.py
aare/func.py
aare/RawFile.py
aare/transform.py
aare/ScanParameters.py
@ -35,22 +40,6 @@ set( PYTHON_FILES
)
#HDF5
if (AARE_HDF5)
find_package(HDF5 1.10 COMPONENTS CXX REQUIRED)
add_definitions(
${HDF5_DEFINITIONS}
)
list(APPEND PYTHON_FILES
aare/Hdf5File.py
)
if(HDF5_FOUND)
add_definitions(-DHDF5_FOUND)
target_link_libraries(_aare PUBLIC ${HDF5_LIBRARIES})
target_include_directories(_aare PUBLIC ${HDF5_INCLUDE_DIRS})
endif()
endif()
# Copy the python files to the build directory
foreach(FILE ${PYTHON_FILES})
configure_file(${FILE} ${CMAKE_BINARY_DIR}/${FILE} )
@ -60,17 +49,30 @@ set_target_properties(_aare PROPERTIES
LIBRARY_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/aare
)
set(PYTHON_EXAMPLES
examples/play.py
examples/fits.py
)
# Copy the examples/scripts to the build directory
configure_file(examples/play.py ${CMAKE_BINARY_DIR}/play.py)
# Copy the python examples to the build directory
foreach(FILE ${PYTHON_EXAMPLES})
configure_file(${FILE} ${CMAKE_BINARY_DIR}/${FILE} )
message(STATUS "Copying ${FILE} to ${CMAKE_BINARY_DIR}/${FILE}")
endforeach(FILE ${PYTHON_EXAMPLES})
if(AARE_INSTALL_PYTHONEXT)
install(TARGETS _aare
install(
TARGETS _aare
EXPORT "${TARGETS_EXPORT_NAME}"
LIBRARY DESTINATION aare
COMPONENT python
)
install(FILES ${PYTHON_FILES} DESTINATION aare)
install(
FILES ${PYTHON_FILES}
DESTINATION aare
COMPONENT python
)
endif()

67
python/aare/ClusterFinder.py Normal file
View File

@ -0,0 +1,67 @@
from ._aare import ClusterFinder_Cluster3x3i, ClusterFinder_Cluster2x2i, ClusterFinderMT_Cluster3x3i, ClusterFinderMT_Cluster2x2i, ClusterCollector_Cluster3x3i, ClusterCollector_Cluster2x2i
from ._aare import ClusterFileSink_Cluster3x3i, ClusterFileSink_Cluster2x2i
import numpy as np
def ClusterFinder(image_size, cluster_size, n_sigma=5, dtype = np.int32, capacity = 1024):
"""
Factory function to create a ClusterFinder object. Provides a cleaner syntax for
the templated ClusterFinder in C++.
"""
if dtype == np.int32 and cluster_size == (3,3):
return ClusterFinder_Cluster3x3i(image_size, n_sigma = n_sigma, capacity=capacity)
elif dtype == np.int32 and cluster_size == (2,2):
return ClusterFinder_Cluster2x2i(image_size, n_sigma = n_sigma, capacity=capacity)
else:
#TODO! add the other formats
raise ValueError(f"Unsupported dtype: {dtype}. Only np.int32 is supported.")
def ClusterFinderMT(image_size, cluster_size = (3,3), dtype=np.int32, n_sigma=5, capacity = 1024, n_threads = 3):
"""
Factory function to create a ClusterFinderMT object. Provides a cleaner syntax for
the templated ClusterFinderMT in C++.
"""
if dtype == np.int32 and cluster_size == (3,3):
return ClusterFinderMT_Cluster3x3i(image_size, n_sigma = n_sigma,
capacity = capacity, n_threads = n_threads)
elif dtype == np.int32 and cluster_size == (2,2):
return ClusterFinderMT_Cluster2x2i(image_size, n_sigma = n_sigma,
capacity = capacity, n_threads = n_threads)
else:
#TODO! add the other formats
raise ValueError(f"Unsupported dtype: {dtype}. Only np.int32 is supported.")
def ClusterCollector(clusterfindermt, cluster_size = (3,3), dtype=np.int32):
"""
Factory function to create a ClusterCollector object. Provides a cleaner syntax for
the templated ClusterCollector in C++.
"""
if dtype == np.int32 and cluster_size == (3,3):
return ClusterCollector_Cluster3x3i(clusterfindermt)
elif dtype == np.int32 and cluster_size == (2,2):
return ClusterCollector_Cluster2x2i(clusterfindermt)
else:
#TODO! add the other formats
raise ValueError(f"Unsupported dtype: {dtype}. Only np.int32 is supported.")
def ClusterFileSink(clusterfindermt, cluster_file, dtype=np.int32):
"""
Factory function to create a ClusterCollector object. Provides a cleaner syntax for
the templated ClusterCollector in C++.
"""
if dtype == np.int32 and clusterfindermt.cluster_size == (3,3):
return ClusterFileSink_Cluster3x3i(clusterfindermt, cluster_file)
elif dtype == np.int32 and clusterfindermt.cluster_size == (2,2):
return ClusterFileSink_Cluster2x2i(clusterfindermt, cluster_file)
else:
#TODO! add the other formats
raise ValueError(f"Unsupported dtype: {dtype}. Only np.int32 is supported.")

11
python/aare/ClusterVector.py Normal file
View File

@ -0,0 +1,11 @@
from ._aare import ClusterVector_Cluster3x3i
import numpy as np
def ClusterVector(cluster_size, dtype = np.int32):
if dtype == np.int32 and cluster_size == (3,3):
return ClusterVector_Cluster3x3i()
else:
raise ValueError(f"Unsupported dtype: {dtype}. Only np.int32 is supported.")

26
python/aare/__init__.py
View File

@ -2,15 +2,33 @@
from . import _aare
from ._aare import File, RawMasterFile, RawSubFile, Hdf5MasterFile
from ._aare import Pedestal, ClusterFinder, VarClusterFinder
from ._aare import File, RawMasterFile, RawSubFile, Hdf5MasterFile, JungfrauDataFile
from ._aare import Pedestal_d, Pedestal_f, ClusterFinder_Cluster3x3i, VarClusterFinder
from ._aare import DetectorType
from ._aare import ClusterFile
from ._aare import ClusterFile_Cluster3x3i as ClusterFile
from ._aare import hitmap
from ._aare import ROI
# from ._aare import ClusterFinderMT, ClusterCollector, ClusterFileSink, ClusterVector_i
from .ClusterFinder import ClusterFinder, ClusterCollector, ClusterFinderMT, ClusterFileSink
from .ClusterVector import ClusterVector
from ._aare import fit_gaus, fit_pol1, fit_scurve, fit_scurve2
from ._aare import Interpolator
from ._aare import calculate_eta2
from ._aare import apply_custom_weights
from .CtbRawFile import CtbRawFile
from .RawFile import RawFile
from .Hdf5File import Hdf5File
from .ScanParameters import ScanParameters
from .utils import random_pixels, random_pixel
from .utils import random_pixels, random_pixel, flat_list, add_colorbar
#make functions available in the top level API
from .func import *

1
python/aare/func.py Normal file
View File

@ -0,0 +1 @@
from ._aare import gaus, pol1, scurve, scurve2

12
python/aare/transform.py
View File

@ -2,6 +2,14 @@ import numpy as np
from . import _aare
class AdcSar04Transform64to16:
def __call__(self, data):
return _aare.adc_sar_04_decode64to16(data)
class AdcSar05Transform64to16:
def __call__(self, data):
return _aare.adc_sar_05_decode64to16(data)
class Moench05Transform:
#Could be moved to C++ without changing the interface
def __init__(self):
@ -45,4 +53,6 @@ class Matterhorn02Transform:
moench05 = Moench05Transform()
moench05_1g = Moench05Transform1g()
moench05_old = Moench05TransformOld()
matterhorn02 = Matterhorn02Transform()
matterhorn02 = Matterhorn02Transform()
adc_sar_04_64to16 = AdcSar04Transform64to16()
adc_sar_05_64to16 = AdcSar05Transform64to16()

15
python/aare/utils.py
View File

@ -1,4 +1,6 @@
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable
def random_pixels(n_pixels, xmin=0, xmax=512, ymin=0, ymax=1024):
"""Return a list of random pixels.
@ -20,4 +22,15 @@ def random_pixel(xmin=0, xmax=512, ymin=0, ymax=1024):
Returns:
tuple: (row, col)
"""
return random_pixels(1, xmin, xmax, ymin, ymax)[0]
return random_pixels(1, xmin, xmax, ymin, ymax)[0]
def flat_list(xss):
"""Flatten a list of lists."""
return [x for xs in xss for x in xs]
def add_colorbar(ax, im, size="5%", pad=0.05):
"""Add a colorbar with the same height as the image."""
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size=size, pad=pad)
plt.colorbar(im, cax=cax)
return ax, im, cax

79
python/examples/fits.py Normal file
View File

@ -0,0 +1,79 @@
import matplotlib.pyplot as plt
import numpy as np
from aare import fit_gaus, fit_pol1
from aare import gaus, pol1
textpm = f"±" #
textmu = f"μ" #
textsigma = f"σ" #
# ================================= Gauss fit =================================
# Parameters
mu = np.random.uniform(1, 100) # Mean of Gaussian
sigma = np.random.uniform(4, 20) # Standard deviation
num_points = 10000 # Number of points for smooth distribution
noise_sigma = 100
# Generate Gaussian distribution
data = np.random.normal(mu, sigma, num_points)
# Generate errors for each point
errors = np.abs(np.random.normal(0, sigma, num_points)) # Errors with mean 0, std 0.5
# Create subplot
fig0, ax0 = plt.subplots(1, 1, num=0, figsize=(12, 8))
x = np.histogram(data, bins=30)[1][:-1] + 0.05
y = np.histogram(data, bins=30)[0]
yerr = errors[:30]
# Add the errors as error bars in the step plot
ax0.errorbar(x, y, yerr=yerr, fmt=". ", capsize=5)
ax0.grid()
par, err = fit_gaus(x, y, yerr)
print(par, err)
x = np.linspace(x[0], x[-1], 1000)
ax0.plot(x, gaus(x, par), marker="")
ax0.set(xlabel="x", ylabel="Counts", title=f"A0 = {par[0]:0.2f}{textpm}{err[0]:0.2f}\n"
f"{textmu} = {par[1]:0.2f}{textpm}{err[1]:0.2f}\n"
f"{textsigma} = {par[2]:0.2f}{textpm}{err[2]:0.2f}\n"
f"(init: {textmu}: {mu:0.2f}, {textsigma}: {sigma:0.2f})")
fig0.tight_layout()
# ================================= pol1 fit =================================
# Parameters
n_points = 40
# Generate random slope and intercept (origin)
slope = np.random.uniform(-10, 10) # Random slope between 0.5 and 2.0
intercept = np.random.uniform(-10, 10) # Random intercept between -10 and 10
# Generate random x values
x_values = np.random.uniform(-10, 10, n_points)
# Calculate y values based on the linear function y = mx + b + error
errors = np.abs(np.random.normal(0, np.random.uniform(1, 5), n_points))
var_points = np.random.normal(0, np.random.uniform(0.1, 2), n_points)
y_values = slope * x_values + intercept + var_points
fig1, ax1 = plt.subplots(1, 1, num=1, figsize=(12, 8))
ax1.errorbar(x_values, y_values, yerr=errors, fmt=". ", capsize=5)
par, err = fit_pol1(x_values, y_values, errors)
x = np.linspace(np.min(x_values), np.max(x_values), 1000)
ax1.plot(x, pol1(x, par), marker="")
ax1.set(xlabel="x", ylabel="y", title=f"a = {par[0]:0.2f}{textpm}{err[0]:0.2f}\n"
f"b = {par[1]:0.2f}{textpm}{err[1]:0.2f}\n"
f"(init: {slope:0.2f}, {intercept:0.2f})")
fig1.tight_layout()
plt.show()

96
python/examples/play.py
View File

@ -1,15 +1,89 @@
import matplotlib.pyplot as plt
import numpy as np
plt.ion()
import sys
sys.path.append('/home/l_msdetect/erik/aare/build')
from aare import RawSubFile, DetectorType, RawFile
from pathlib import Path
from aare import ClusterFile
path = Path("/home/l_msdetect/erik/data/aare-test-data/raw/jungfrau/")
f = RawSubFile(path/"jungfrau_single_d0_f0_0.raw", DetectorType.Jungfrau, 512, 1024, 16)
base = Path('~/data/aare_test_data/clusters').expanduser()
f = ClusterFile(base / 'beam_En700eV_-40deg_300V_10us_d0_f0_100.clust')
# f = ClusterFile(base / 'single_frame_97_clustrers.clust')
# f = RawFile(path/"jungfrau_single_master_0.json")
for i in range(10):
fn, cl = f.read_frame()
print(fn, cl.size)
# from aare._aare import ClusterVector_i, Interpolator
# import pickle
# import numpy as np
# import matplotlib.pyplot as plt
# import boost_histogram as bh
# import torch
# import math
# import time
# def gaussian_2d(mx, my, sigma = 1, res=100, grid_size = 2):
# """
# Generate a 2D gaussian as position mx, my, with sigma=sigma.
# The gaussian is placed on a 2x2 pixel matrix with resolution
# res in one dimesion.
# """
# x = torch.linspace(0, pixel_size*grid_size, res)
# x,y = torch.meshgrid(x,x, indexing="ij")
# return 1 / (2*math.pi*sigma**2) * \
# torch.exp(-((x - my)**2 / (2*sigma**2) + (y - mx)**2 / (2*sigma**2)))
# scale = 1000 #Scale factor when converting to integer
# pixel_size = 25 #um
# grid = 2
# resolution = 100
# sigma_um = 10
# xa = np.linspace(0,grid*pixel_size,resolution)
# ticks = [0, 25, 50]
# hit = np.array((20,20))
# etahist_fname = "/home/l_msdetect/erik/tmp/test_hist.pkl"
# local_resolution = 99
# grid_size = 3
# xaxis = np.linspace(0,grid_size*pixel_size, local_resolution)
# t = gaussian_2d(hit[0],hit[1], grid_size = grid_size, sigma = 10, res = local_resolution)
# pixels = t.reshape(grid_size, t.shape[0] // grid_size, grid_size, t.shape[1] // grid_size).sum(axis = 3).sum(axis = 1)
# pixels = pixels.numpy()
# pixels = (pixels*scale).astype(np.int32)
# v = ClusterVector_i(3,3)
# v.push_back(1,1, pixels)
# with open(etahist_fname, "rb") as f:
# hist = pickle.load(f)
# eta = hist.view().copy()
# etabinsx = np.array(hist.axes.edges.T[0].flat)
# etabinsy = np.array(hist.axes.edges.T[1].flat)
# ebins = np.array(hist.axes.edges.T[2].flat)
# p = Interpolator(eta, etabinsx[0:-1], etabinsy[0:-1], ebins[0:-1])
# #Generate the hit
# tmp = p.interpolate(v)
# print(f'tmp:{tmp}')
# pos = np.array((tmp['x'], tmp['y']))*25
# print(pixels)
# fig, ax = plt.subplots(figsize = (7,7))
# ax.pcolormesh(xaxis, xaxis, t)
# ax.plot(*pos, 'o')
# ax.set_xticks([0,25,50,75])
# ax.set_yticks([0,25,50,75])
# ax.set_xlim(0,75)
# ax.set_ylim(0,75)
# ax.grid()
# print(f'{hit=}')
# print(f'{pos=}')

64
python/src/bind_Cluster.hpp Normal file
View File

@ -0,0 +1,64 @@
#include "aare/Cluster.hpp"
#include <cstdint>
#include <fmt/format.h>
#include <filesystem>
#include <pybind11/pybind11.h>
#include <pybind11/numpy.h>
#include <pybind11/stl.h>
#include <pybind11/stl_bind.h>
namespace py = pybind11;
using pd_type = double;
using namespace aare;
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wunused-parameter"
template <typename Type, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
typename CoordType>
void define_Cluster(py::module &m, const std::string &typestr) {
auto class_name = fmt::format("Cluster{}", typestr);
py::class_<Cluster<Type, ClusterSizeX, ClusterSizeY, CoordType>>(
m, class_name.c_str(), py::buffer_protocol())
.def(py::init([](uint8_t x, uint8_t y, py::array_t<Type> data) {
py::buffer_info buf_info = data.request();
Cluster<Type, ClusterSizeX, ClusterSizeY, CoordType> cluster;
cluster.x = x;
cluster.y = y;
auto r = data.template unchecked<1>(); // no bounds checks
for (py::ssize_t i = 0; i < data.size(); ++i) {
cluster.data[i] = r(i);
}
return cluster;
}));
/*
//TODO! Review if to keep or not
.def_property(
"data",
[](ClusterType &c) -> py::array {
return py::array(py::buffer_info(
c.data, sizeof(Type),
py::format_descriptor<Type>::format(), // Type
// format
1, // Number of dimensions
{static_cast<ssize_t>(ClusterSizeX *
ClusterSizeY)}, // Shape (flattened)
{sizeof(Type)} // Stride (step size between elements)
));
},
[](ClusterType &c, py::array_t<Type> arr) {
py::buffer_info buf_info = arr.request();
Type *ptr = static_cast<Type *>(buf_info.ptr);
std::copy(ptr, ptr + ClusterSizeX * ClusterSizeY,
c.data); // TODO dont iterate over centers!!!
});
*/
}
#pragma GCC diagnostic pop

46
python/src/bind_ClusterCollector.hpp Normal file
View File

@ -0,0 +1,46 @@
#include "aare/ClusterCollector.hpp"
#include "aare/ClusterFileSink.hpp"
#include "aare/ClusterFinder.hpp"
#include "aare/ClusterFinderMT.hpp"
#include "aare/ClusterVector.hpp"
#include "aare/NDView.hpp"
#include "aare/Pedestal.hpp"
#include "np_helper.hpp"
#include <cstdint>
#include <filesystem>
#include <pybind11/pybind11.h>
#include <pybind11/stl.h>
#include <pybind11/stl_bind.h>
namespace py = pybind11;
using pd_type = double;
using namespace aare;
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wunused-parameter"
template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
typename CoordType = uint16_t>
void define_ClusterCollector(py::module &m,
const std::string &typestr) {
auto class_name = fmt::format("ClusterCollector_{}", typestr);
using ClusterType = Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>;
py::class_<ClusterCollector<ClusterType>>(m, class_name.c_str())
.def(py::init<ClusterFinderMT<ClusterType, uint16_t, double> *>())
.def("stop", &ClusterCollector<ClusterType>::stop)
.def(
"steal_clusters",
[](ClusterCollector<ClusterType> &self) {
auto v = new std::vector<ClusterVector<ClusterType>>(
self.steal_clusters());
return v; // TODO change!!!
},
py::return_value_policy::take_ownership);
}
#pragma GCC diagnostic pop

94
python/src/bind_ClusterFile.hpp Normal file
View File

@ -0,0 +1,94 @@
#include "aare/CalculateEta.hpp"
#include "aare/ClusterFile.hpp"
#include "aare/defs.hpp"
#include <cstdint>
#include <filesystem>
#include <pybind11/iostream.h>
#include <pybind11/numpy.h>
#include <pybind11/pybind11.h>
#include <pybind11/stl.h>
#include <pybind11/stl/filesystem.h>
#include <string>
// Disable warnings for unused parameters, as we ignore some
// in the __exit__ method
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wunused-parameter"
namespace py = pybind11;
using namespace ::aare;
template <typename Type, uint8_t CoordSizeX, uint8_t CoordSizeY,
typename CoordType = uint16_t>
void define_ClusterFile(py::module &m,
const std::string &typestr) {
using ClusterType = Cluster<Type, CoordSizeX, CoordSizeY, CoordType>;
auto class_name = fmt::format("ClusterFile_{}", typestr);
py::class_<ClusterFile<ClusterType>>(m, class_name.c_str())
.def(py::init<const std::filesystem::path &, size_t,
const std::string &>(),
py::arg(), py::arg("chunk_size") = 1000, py::arg("mode") = "r")
.def(
"read_clusters",
[](ClusterFile<ClusterType> &self, size_t n_clusters) {
auto v = new ClusterVector<ClusterType>(
self.read_clusters(n_clusters));
return v;
},
py::return_value_policy::take_ownership)
.def("read_frame",
[](ClusterFile<ClusterType> &self) {
auto v = new ClusterVector<ClusterType>(self.read_frame());
return v;
})
.def("set_roi", &ClusterFile<ClusterType>::set_roi)
.def(
"set_noise_map",
[](ClusterFile<ClusterType> &self, py::array_t<int32_t> noise_map) {
auto view = make_view_2d(noise_map);
self.set_noise_map(view);
})
.def("set_gain_map",
[](ClusterFile<ClusterType> &self, py::array_t<double> gain_map) {
auto view = make_view_2d(gain_map);
self.set_gain_map(view);
})
.def("close", &ClusterFile<ClusterType>::close)
.def("write_frame", &ClusterFile<ClusterType>::write_frame)
.def("__enter__", [](ClusterFile<ClusterType> &self) { return &self; })
.def("__exit__",
[](ClusterFile<ClusterType> &self,
const std::optional<pybind11::type> &exc_type,
const std::optional<pybind11::object> &exc_value,
const std::optional<pybind11::object> &traceback) {
self.close();
})
.def("__iter__", [](ClusterFile<ClusterType> &self) { return &self; })
.def("__next__", [](ClusterFile<ClusterType> &self) {
auto v = new ClusterVector<ClusterType>(
self.read_clusters(self.chunk_size()));
if (v->size() == 0) {
throw py::stop_iteration();
}
return v;
});
}
template <typename Type, uint8_t CoordSizeX, uint8_t CoordSizeY,
typename CoordType = uint16_t>
void register_calculate_eta(py::module &m) {
using ClusterType = Cluster<Type, CoordSizeX, CoordSizeY, CoordType>;
m.def("calculate_eta2",
[](const aare::ClusterVector<ClusterType> &clusters) {
auto eta2 = new NDArray<double, 2>(calculate_eta2(clusters));
return return_image_data(eta2);
});
}
#pragma GCC diagnostic pop

44
python/src/bind_ClusterFileSink.hpp Normal file
View File

@ -0,0 +1,44 @@
#include "aare/ClusterCollector.hpp"
#include "aare/ClusterFileSink.hpp"
#include "aare/ClusterFinder.hpp"
#include "aare/ClusterFinderMT.hpp"
#include "aare/ClusterVector.hpp"
#include "aare/NDView.hpp"
#include "aare/Pedestal.hpp"
#include "np_helper.hpp"
#include <cstdint>
#include <filesystem>
#include <pybind11/pybind11.h>
#include <pybind11/stl.h>
#include <pybind11/stl_bind.h>
namespace py = pybind11;
using pd_type = double;
using namespace aare;
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wunused-parameter"
template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
typename CoordType = uint16_t>
void define_ClusterFileSink(py::module &m,
const std::string &typestr) {
auto class_name = fmt::format("ClusterFileSink_{}", typestr);
using ClusterType = Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>;
py::class_<ClusterFileSink<ClusterType>>(m, class_name.c_str())
.def(py::init<ClusterFinderMT<ClusterType, uint16_t, double> *,
const std::filesystem::path &>())
.def("stop", &ClusterFileSink<ClusterType>::stop);
}
#pragma GCC diagnostic pop

77
python/src/bind_ClusterFinder.hpp Normal file
View File

@ -0,0 +1,77 @@
#include "aare/ClusterCollector.hpp"
#include "aare/ClusterFileSink.hpp"
#include "aare/ClusterFinder.hpp"
#include "aare/ClusterFinderMT.hpp"
#include "aare/ClusterVector.hpp"
#include "aare/NDView.hpp"
#include "aare/Pedestal.hpp"
#include "np_helper.hpp"
#include <cstdint>
#include <filesystem>
#include <pybind11/pybind11.h>
#include <pybind11/stl.h>
#include <pybind11/stl_bind.h>
namespace py = pybind11;
using pd_type = double;
using namespace aare;
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wunused-parameter"
template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
typename CoordType = uint16_t>
void define_ClusterFinder(py::module &m, const std::string &typestr) {
auto class_name = fmt::format("ClusterFinder_{}", typestr);
using ClusterType = Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>;
py::class_<ClusterFinder<ClusterType, uint16_t, pd_type>>(
m, class_name.c_str())
.def(py::init<Shape<2>, pd_type, size_t>(), py::arg("image_size"),
py::arg("n_sigma") = 5.0, py::arg("capacity") = 1'000'000)
.def("push_pedestal_frame",
[](ClusterFinder<ClusterType, uint16_t, pd_type> &self,
py::array_t<uint16_t> frame) {
auto view = make_view_2d(frame);
self.push_pedestal_frame(view);
})
.def("clear_pedestal",
&ClusterFinder<ClusterType, uint16_t, pd_type>::clear_pedestal)
.def_property_readonly(
"pedestal",
[](ClusterFinder<ClusterType, uint16_t, pd_type> &self) {
auto pd = new NDArray<pd_type, 2>{};
*pd = self.pedestal();
return return_image_data(pd);
})
.def_property_readonly(
"noise",
[](ClusterFinder<ClusterType, uint16_t, pd_type> &self) {
auto arr = new NDArray<pd_type, 2>{};
*arr = self.noise();
return return_image_data(arr);
})
.def(
"steal_clusters",
[](ClusterFinder<ClusterType, uint16_t, pd_type> &self,
bool realloc_same_capacity) {
ClusterVector<ClusterType> clusters =
self.steal_clusters(realloc_same_capacity);
return clusters;
},
py::arg("realloc_same_capacity") = false)
.def(
"find_clusters",
[](ClusterFinder<ClusterType, uint16_t, pd_type> &self,
py::array_t<uint16_t> frame, uint64_t frame_number) {
auto view = make_view_2d(frame);
self.find_clusters(view, frame_number);
return;
},
py::arg(), py::arg("frame_number") = 0);
}
#pragma GCC diagnostic pop

81
python/src/bind_ClusterFinderMT.hpp Normal file
View File

@ -0,0 +1,81 @@
#include "aare/ClusterCollector.hpp"
#include "aare/ClusterFileSink.hpp"
#include "aare/ClusterFinder.hpp"
#include "aare/ClusterFinderMT.hpp"
#include "aare/ClusterVector.hpp"
#include "aare/NDView.hpp"
#include "aare/Pedestal.hpp"
#include "np_helper.hpp"
#include <cstdint>
#include <filesystem>
#include <pybind11/pybind11.h>
#include <pybind11/stl.h>
#include <pybind11/stl_bind.h>
namespace py = pybind11;
using pd_type = double;
using namespace aare;
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wunused-parameter"
template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
typename CoordType = uint16_t>
void define_ClusterFinderMT(py::module &m,
const std::string &typestr) {
auto class_name = fmt::format("ClusterFinderMT_{}", typestr);
using ClusterType = Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>;
py::class_<ClusterFinderMT<ClusterType, uint16_t, pd_type>>(
m, class_name.c_str())
.def(py::init<Shape<2>, pd_type, size_t, size_t>(),
py::arg("image_size"), py::arg("n_sigma") = 5.0,
py::arg("capacity") = 2048, py::arg("n_threads") = 3)
.def("push_pedestal_frame",
[](ClusterFinderMT<ClusterType, uint16_t, pd_type> &self,
py::array_t<uint16_t> frame) {
auto view = make_view_2d(frame);
self.push_pedestal_frame(view);
})
.def(
"find_clusters",
[](ClusterFinderMT<ClusterType, uint16_t, pd_type> &self,
py::array_t<uint16_t> frame, uint64_t frame_number) {
auto view = make_view_2d(frame);
self.find_clusters(view, frame_number);
return;
},
py::arg(), py::arg("frame_number") = 0)
.def_property_readonly("cluster_size", [](ClusterFinderMT<ClusterType, uint16_t, pd_type> &self){
return py::make_tuple(ClusterSizeX, ClusterSizeY);
})
.def("clear_pedestal",
&ClusterFinderMT<ClusterType, uint16_t, pd_type>::clear_pedestal)
.def("sync", &ClusterFinderMT<ClusterType, uint16_t, pd_type>::sync)
.def("stop", &ClusterFinderMT<ClusterType, uint16_t, pd_type>::stop)
.def("start", &ClusterFinderMT<ClusterType, uint16_t, pd_type>::start)
.def(
"pedestal",
[](ClusterFinderMT<ClusterType, uint16_t, pd_type> &self,
size_t thread_index) {
auto pd = new NDArray<pd_type, 2>{};
*pd = self.pedestal(thread_index);
return return_image_data(pd);
},
py::arg("thread_index") = 0)
.def(
"noise",
[](ClusterFinderMT<ClusterType, uint16_t, pd_type> &self,
size_t thread_index) {
auto arr = new NDArray<pd_type, 2>{};
*arr = self.noise(thread_index);
return return_image_data(arr);
},
py::arg("thread_index") = 0);
}
#pragma GCC diagnostic pop

106
python/src/bind_ClusterVector.hpp Normal file
View File

@ -0,0 +1,106 @@
#include "aare/ClusterCollector.hpp"
#include "aare/ClusterFileSink.hpp"
#include "aare/ClusterFinder.hpp"
#include "aare/ClusterFinderMT.hpp"
#include "aare/ClusterVector.hpp"
#include "aare/NDView.hpp"
#include "aare/Pedestal.hpp"
#include "np_helper.hpp"
#include <cstdint>
#include <filesystem>
#include <pybind11/pybind11.h>
#include <pybind11/stl.h>
#include <pybind11/stl_bind.h>
namespace py = pybind11;
using pd_type = double;
using namespace aare;
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wunused-parameter"
template <typename Type, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
typename CoordType = uint16_t>
void define_ClusterVector(py::module &m, const std::string &typestr) {
using ClusterType = Cluster<Type, ClusterSizeX, ClusterSizeY, CoordType>;
auto class_name = fmt::format("ClusterVector_{}", typestr);
py::class_<ClusterVector<
Cluster<Type, ClusterSizeX, ClusterSizeY, CoordType>, void>>(
m, class_name.c_str(),
py::buffer_protocol())
.def(py::init()) // TODO change!!!
.def("push_back",
[](ClusterVector<ClusterType> &self, const ClusterType &cluster) {
self.push_back(cluster);
})
.def("sum",
[](ClusterVector<ClusterType> &self) {
auto *vec = new std::vector<Type>(self.sum());
return return_vector(vec);
})
.def("sum_2x2", [](ClusterVector<ClusterType> &self){
auto *vec = new std::vector<Type>(self.sum_2x2());
return return_vector(vec);
})
.def_property_readonly("size", &ClusterVector<ClusterType>::size)
.def("item_size", &ClusterVector<ClusterType>::item_size)
.def_property_readonly("fmt",
[typestr](ClusterVector<ClusterType> &self) {
return fmt_format<ClusterType>;
})
.def_property_readonly("cluster_size_x",
&ClusterVector<ClusterType>::cluster_size_x)
.def_property_readonly("cluster_size_y",
&ClusterVector<ClusterType>::cluster_size_y)
.def_property_readonly("capacity",
&ClusterVector<ClusterType>::capacity)
.def_property("frame_number", &ClusterVector<ClusterType>::frame_number,
&ClusterVector<ClusterType>::set_frame_number)
.def_buffer(
[typestr](ClusterVector<ClusterType> &self) -> py::buffer_info {
return py::buffer_info(
self.data(), /* Pointer to buffer */
self.item_size(), /* Size of one scalar */
fmt_format<ClusterType>, /* Format descriptor */
1, /* Number of dimensions */
{self.size()}, /* Buffer dimensions */
{self.item_size()} /* Strides (in bytes) for each index */
);
});
// Free functions using ClusterVector
m.def("hitmap",
[](std::array<size_t, 2> image_size, ClusterVector<ClusterType> &cv) {
// Create a numpy array to hold the hitmap
// The shape of the array is (image_size[0], image_size[1])
// note that the python array is passed as [row, col] which
// is the opposite of the clusters [x,y]
py::array_t<int32_t> hitmap(image_size);
auto r = hitmap.mutable_unchecked<2>();
// Initialize hitmap to 0
for (py::ssize_t i = 0; i < r.shape(0); i++)
for (py::ssize_t j = 0; j < r.shape(1); j++)
r(i, j) = 0;
// Loop over the clusters and increment the hitmap
// Skip out of bound clusters
for (const auto &cluster : cv) {
auto x = cluster.x;
auto y = cluster.y;
if (x < image_size[1] && y < image_size[0])
r(cluster.y, cluster.x) += 1;
}
return hitmap;
});
}
#pragma GCC diagnostic pop

52
python/src/cluster.hpp
View File

@ -1,52 +0,0 @@
#include "aare/ClusterFinder.hpp"
#include "aare/NDView.hpp"
#include "aare/Pedestal.hpp"
#include "np_helper.hpp"
#include <cstdint>
#include <filesystem>
#include <pybind11/pybind11.h>
#include <pybind11/stl.h>
namespace py = pybind11;
void define_cluster_finder_bindings(py::module &m) {
py::class_<ClusterFinder<uint16_t, double>>(m, "ClusterFinder")
.def(py::init<Shape<2>, Shape<2>>())
.def("push_pedestal_frame",
[](ClusterFinder<uint16_t, double> &self,
py::array_t<uint16_t> frame) {
auto view = make_view_2d(frame);
self.push_pedestal_frame(view);
})
.def("pedestal",
[](ClusterFinder<uint16_t, double> &self) {
auto pd = new NDArray<double, 2>{};
*pd = self.pedestal();
return return_image_data(pd);
})
.def("find_clusters_without_threshold",
[](ClusterFinder<uint16_t, double> &self,
py::array_t<uint16_t> frame) {
auto view = make_view_2d(frame);
auto clusters = self.find_clusters_without_threshold(view);
return clusters;
});
py::class_<DynamicCluster>(m, "DynamicCluster", py::buffer_protocol())
.def(py::init<int, int, Dtype>())
.def("size", &DynamicCluster::size)
.def("begin", &DynamicCluster::begin)
.def("end", &DynamicCluster::end)
.def_readwrite("x", &DynamicCluster::x)
.def_readwrite("y", &DynamicCluster::y)
.def_buffer([](DynamicCluster &c) -> py::buffer_info {
return py::buffer_info(c.data(), c.dt.bytes(), c.dt.format_descr(),
1, {c.size()}, {c.dt.bytes()});
})
.def("__repr__", [](const DynamicCluster &a) {
return "<DynamicCluster: x: " + std::to_string(a.x) +
", y: " + std::to_string(a.y) + ">";
});
}

50
python/src/cluster_file.hpp
View File

@ -1,50 +0,0 @@
#include "aare/ClusterFile.hpp"
#include "aare/defs.hpp"
#include <cstdint>
#include <filesystem>
#include <pybind11/iostream.h>
#include <pybind11/numpy.h>
#include <pybind11/pybind11.h>
#include <pybind11/stl.h>
#include <pybind11/stl/filesystem.h>
#include <string>
namespace py = pybind11;
using namespace ::aare;
void define_cluster_file_io_bindings(py::module &m) {
PYBIND11_NUMPY_DTYPE(Cluster, x, y, data);
py::class_<ClusterFile>(m, "ClusterFile")
.def(py::init<const std::filesystem::path &, size_t>(), py::arg(), py::arg("chunk_size") = 1000)
.def("read_clusters",
[](ClusterFile &self, size_t n_clusters) {
auto* vec = new std::vector<Cluster>(self.read_clusters(n_clusters));
return return_vector(vec);
})
.def("read_frame",
[](ClusterFile &self) {
int32_t frame_number;
auto* vec = new std::vector<Cluster>(self.read_frame(frame_number));
return py::make_tuple(frame_number, return_vector(vec));
})
.def("read_cluster_with_cut",
[](ClusterFile &self, size_t n_clusters, py::array_t<double> noise_map, int nx, int ny) {
auto view = make_view_2d(noise_map);
auto* vec = new std::vector<Cluster>(self.read_cluster_with_cut(n_clusters, view.data(), nx, ny));
return return_vector(vec);
})
.def("__enter__", [](ClusterFile &self) { return &self; })
.def("__exit__", [](ClusterFile &self) { self.close();})
.def("__iter__", [](ClusterFile &self) { return &self; })
.def("__next__", [](ClusterFile &self) {
auto vec = new std::vector<Cluster>(self.read_clusters(self.chunk_size()));
if(vec->size() == 0) {
throw py::stop_iteration();
}
return return_vector(vec);
});
}

113
python/src/ctb_raw_file.hpp
View File

@ -7,8 +7,11 @@
#include "aare/RawSubFile.hpp"
#include "aare/defs.hpp"
#include "aare/decode.hpp"
// #include "aare/fClusterFileV2.hpp"
#include "np_helper.hpp"
#include <cstdint>
#include <filesystem>
#include <pybind11/iostream.h>
@ -23,35 +26,95 @@ using namespace ::aare;
void define_ctb_raw_file_io_bindings(py::module &m) {
py::class_<CtbRawFile>(m, "CtbRawFile")
.def(py::init<const std::filesystem::path &>())
.def("read_frame",
[](CtbRawFile &self) {
size_t image_size = self.image_size_in_bytes();
py::array image;
std::vector<ssize_t> shape;
shape.reserve(2);
shape.push_back(1);
shape.push_back(image_size);
m.def("adc_sar_05_decode64to16", [](py::array_t<uint8_t> input) {
py::array_t<DetectorHeader> header(1);
if(input.ndim() != 2){
throw std::runtime_error("Only 2D arrays are supported at this moment");
}
// always read bytes
image = py::array_t<uint8_t>(shape);
//Create a 2D output array with the same shape as the input
std::vector<ssize_t> shape{input.shape(0), input.shape(1)/static_cast<ssize_t>(bits_per_byte)};
py::array_t<uint16_t> output(shape);
self.read_into(
reinterpret_cast<std::byte *>(image.mutable_data()),
header.mutable_data());
//Create a view of the input and output arrays
NDView<uint64_t, 2> input_view(reinterpret_cast<uint64_t*>(input.mutable_data()), {output.shape(0), output.shape(1)});
NDView<uint16_t, 2> output_view(output.mutable_data(), {output.shape(0), output.shape(1)});
return py::make_tuple(header, image);
})
.def("seek", &CtbRawFile::seek)
.def("tell", &CtbRawFile::tell)
.def("master", &CtbRawFile::master)
adc_sar_05_decode64to16(input_view, output_view);
.def_property_readonly("image_size_in_bytes",
&CtbRawFile::image_size_in_bytes)
return output;
});
.def_property_readonly("frames_in_file", &CtbRawFile::frames_in_file);
}
m.def("adc_sar_04_decode64to16", [](py::array_t<uint8_t> input) {
if(input.ndim() != 2){
throw std::runtime_error("Only 2D arrays are supported at this moment");
}
//Create a 2D output array with the same shape as the input
std::vector<ssize_t> shape{input.shape(0), input.shape(1)/static_cast<ssize_t>(bits_per_byte)};
py::array_t<uint16_t> output(shape);
//Create a view of the input and output arrays
NDView<uint64_t, 2> input_view(reinterpret_cast<uint64_t*>(input.mutable_data()), {output.shape(0), output.shape(1)});
NDView<uint16_t, 2> output_view(output.mutable_data(), {output.shape(0), output.shape(1)});
adc_sar_04_decode64to16(input_view, output_view);
return output;
});
m.def(
"apply_custom_weights",
[](py::array_t<uint16_t, py::array::c_style | py::array::forcecast> &input,
py::array_t<double, py::array::c_style | py::array::forcecast>
&weights) {
// Create new array with same shape as the input array (uninitialized values)
py::buffer_info buf = input.request();
py::array_t<double> output(buf.shape);
// Use NDViews to call into the C++ library
auto weights_view = make_view_1d(weights);
NDView<uint16_t, 1> input_view(input.mutable_data(), {input.size()});
NDView<double, 1> output_view(output.mutable_data(), {output.size()});
apply_custom_weights(input_view, output_view, weights_view);
return output;
});
py::class_<CtbRawFile>(m, "CtbRawFile")
.def(py::init<const std::filesystem::path &>())
.def("read_frame",
[](CtbRawFile &self) {
size_t image_size = self.image_size_in_bytes();
py::array image;
std::vector<ssize_t> shape;
shape.reserve(2);
shape.push_back(1);
shape.push_back(image_size);
py::array_t<DetectorHeader> header(1);
// always read bytes
image = py::array_t<uint8_t>(shape);
self.read_into(reinterpret_cast<std::byte *>(image.mutable_data()),
header.mutable_data());
return py::make_tuple(header, image);
})
.def("seek", &CtbRawFile::seek)
.def("tell", &CtbRawFile::tell)
.def("master", &CtbRawFile::master)
.def_property_readonly("image_size_in_bytes",
&CtbRawFile::image_size_in_bytes)
.def_property_readonly("frames_in_file", &CtbRawFile::frames_in_file);
}

79
python/src/file.hpp
View File

@ -25,6 +25,14 @@
namespace py = pybind11;
using namespace ::aare;
//Disable warnings for unused parameters, as we ignore some
//in the __exit__ method
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wunused-parameter"
void define_file_io_bindings(py::module &m) {
@ -56,7 +64,8 @@ void define_file_io_bindings(py::module &m) {
.def(py::init<const std::filesystem::path &, const std::string &,
const FileConfig &>())
.def("frame_number", &File::frame_number)
.def("frame_number", py::overload_cast<>(&File::frame_number))
.def("frame_number", py::overload_cast<size_t>(&File::frame_number))
.def_property_readonly("bytes_per_frame", &File::bytes_per_frame)
.def_property_readonly("pixels_per_frame", &File::pixels_per_frame)
.def("seek", &File::seek)
@ -128,8 +137,41 @@ void define_file_io_bindings(py::module &m) {
self.read_into(reinterpret_cast<std::byte *>(image.mutable_data()),
n_frames);
return image;
})
.def("__enter__", [](File &self) { return &self; })
.def("__exit__",
[](File &self,
const std::optional<pybind11::type> &exc_type,
const std::optional<pybind11::object> &exc_value,
const std::optional<pybind11::object> &traceback) {
// self.close();
})
.def("__iter__", [](File &self) { return &self; })
.def("__next__", [](File &self) {
try{
const uint8_t item_size = self.bytes_per_pixel();
py::array image;
std::vector<ssize_t> shape;
shape.reserve(2);
shape.push_back(self.rows());
shape.push_back(self.cols());
if (item_size == 1) {
image = py::array_t<uint8_t>(shape);
} else if (item_size == 2) {
image = py::array_t<uint16_t>(shape);
} else if (item_size == 4) {
image = py::array_t<uint32_t>(shape);
}
self.read_into(
reinterpret_cast<std::byte *>(image.mutable_data()));
return image;
}catch(std::runtime_error &e){
throw py::stop_iteration();
}
});
py::class_<FileConfig>(m, "FileConfig")
.def(py::init<>())
.def_readwrite("rows", &FileConfig::rows)
@ -161,6 +203,8 @@ void define_file_io_bindings(py::module &m) {
py::class_<ROI>(m, "ROI")
.def(py::init<>())
.def(py::init<ssize_t, ssize_t, ssize_t, ssize_t>(), py::arg("xmin"),
py::arg("xmax"), py::arg("ymin"), py::arg("ymax"))
.def_readwrite("xmin", &ROI::xmin)
.def_readwrite("xmax", &ROI::xmax)
.def_readwrite("ymin", &ROI::ymin)
@ -178,38 +222,11 @@ void define_file_io_bindings(py::module &m) {
py::class_<RawSubFile>(m, "RawSubFile")
.def(py::init<const std::filesystem::path &, DetectorType, size_t,
size_t, size_t>())
.def_property_readonly("bytes_per_frame", &RawSubFile::bytes_per_frame)
.def_property_readonly("pixels_per_frame",
&RawSubFile::pixels_per_frame)
.def("seek", &RawSubFile::seek)
.def("tell", &RawSubFile::tell)
.def_property_readonly("rows", &RawSubFile::rows)
.def_property_readonly("cols", &RawSubFile::cols)
.def("read_frame",
[](RawSubFile &self) {
const uint8_t item_size = self.bytes_per_pixel();
py::array image;
std::vector<ssize_t> shape;
shape.reserve(2);
shape.push_back(self.rows());
shape.push_back(self.cols());
if (item_size == 1) {
image = py::array_t<uint8_t>(shape);
} else if (item_size == 2) {
image = py::array_t<uint16_t>(shape);
} else if (item_size == 4) {
image = py::array_t<uint32_t>(shape);
}
fmt::print("item_size: {} rows: {} cols: {}\n", item_size, self.rows(), self.cols());
self.read_into(
reinterpret_cast<std::byte *>(image.mutable_data()));
return image;
});
#pragma GCC diagnostic pop
// py::class_<ClusterHeader>(m, "ClusterHeader")
// .def(py::init<>())
// .def_readwrite("frame_number", &ClusterHeader::frame_number)

465
python/src/fit.hpp Normal file
View File

@ -0,0 +1,465 @@
#include <cstdint>
#include <filesystem>
#include <pybind11/pybind11.h>
#include <pybind11/stl.h>
#include <pybind11/stl_bind.h>
#include "aare/Fit.hpp"
namespace py = pybind11;
using namespace pybind11::literals;
void define_fit_bindings(py::module &m) {
// TODO! Evaluate without converting to double
m.def(
"gaus",
[](py::array_t<double, py::array::c_style | py::array::forcecast> x,
py::array_t<double, py::array::c_style | py::array::forcecast> par) {
auto x_view = make_view_1d(x);
auto par_view = make_view_1d(par);
auto y = new NDArray<double, 1>{aare::func::gaus(x_view, par_view)};
return return_image_data(y);
},
R"(
Evaluate a 1D Gaussian function for all points in x using parameters par.
Parameters
----------
x : array_like
The points at which to evaluate the Gaussian function.
par : array_like
The parameters of the Gaussian function. The first element is the amplitude, the second element is the mean, and the third element is the standard deviation.
)",
py::arg("x"), py::arg("par"));
m.def(
"pol1",
[](py::array_t<double, py::array::c_style | py::array::forcecast> x,
py::array_t<double, py::array::c_style | py::array::forcecast> par) {
auto x_view = make_view_1d(x);
auto par_view = make_view_1d(par);
auto y = new NDArray<double, 1>{aare::func::pol1(x_view, par_view)};
return return_image_data(y);
},
R"(
Evaluate a 1D polynomial function for all points in x using parameters par. (p0+p1*x)
Parameters
----------
x : array_like
The points at which to evaluate the polynomial function.
par : array_like
The parameters of the polynomial function. The first element is the intercept, and the second element is the slope.
)",
py::arg("x"), py::arg("par"));
m.def(
"scurve",
[](py::array_t<double, py::array::c_style | py::array::forcecast> x,
py::array_t<double, py::array::c_style | py::array::forcecast> par) {
auto x_view = make_view_1d(x);
auto par_view = make_view_1d(par);
auto y = new NDArray<double, 1>{aare::func::scurve(x_view, par_view)};
return return_image_data(y);
},
R"(
Evaluate a 1D scurve function for all points in x using parameters par.
Parameters
----------
x : array_like
The points at which to evaluate the scurve function.
par : array_like
The parameters of the scurve function. The first element is the background slope, the second element is the background intercept, the third element is the mean, the fourth element is the standard deviation, the fifth element is inflexion point count number, and the sixth element is C.
)",
py::arg("x"), py::arg("par"));
m.def(
"scurve2",
[](py::array_t<double, py::array::c_style | py::array::forcecast> x,
py::array_t<double, py::array::c_style | py::array::forcecast> par) {
auto x_view = make_view_1d(x);
auto par_view = make_view_1d(par);
auto y = new NDArray<double, 1>{aare::func::scurve2(x_view, par_view)};
return return_image_data(y);
},
R"(
Evaluate a 1D scurve2 function for all points in x using parameters par.
Parameters
----------
x : array_like
The points at which to evaluate the scurve function.
par : array_like
The parameters of the scurve2 function. The first element is the background slope, the second element is the background intercept, the third element is the mean, the fourth element is the standard deviation, the fifth element is inflexion point count number, and the sixth element is C.
)",
py::arg("x"), py::arg("par"));
m.def(
"fit_gaus",
[](py::array_t<double, py::array::c_style | py::array::forcecast> x,
py::array_t<double, py::array::c_style | py::array::forcecast> y,
int n_threads) {
if (y.ndim() == 3) {
auto par = new NDArray<double, 3>{};
auto y_view = make_view_3d(y);
auto x_view = make_view_1d(x);
*par = aare::fit_gaus(x_view, y_view, n_threads);
return return_image_data(par);
} else if (y.ndim() == 1) {
auto par = new NDArray<double, 1>{};
auto y_view = make_view_1d(y);
auto x_view = make_view_1d(x);
*par = aare::fit_gaus(x_view, y_view);
return return_image_data(par);
} else {
throw std::runtime_error("Data must be 1D or 3D");
}
},
R"(
Fit a 1D Gaussian to data.
Parameters
----------
x : array_like
The x values.
y : array_like
The y values.
n_threads : int, optional
The number of threads to use. Default is 4.
)",
py::arg("x"), py::arg("y"), py::arg("n_threads") = 4);
m.def(
"fit_gaus",
[](py::array_t<double, py::array::c_style | py::array::forcecast> x,
py::array_t<double, py::array::c_style | py::array::forcecast> y,
py::array_t<double, py::array::c_style | py::array::forcecast> y_err,
int n_threads) {
if (y.ndim() == 3) {
// Allocate memory for the output
// Need to have pointers to allow python to manage
// the memory
auto par = new NDArray<double, 3>({y.shape(0), y.shape(1), 3});
auto par_err =
new NDArray<double, 3>({y.shape(0), y.shape(1), 3});
auto chi2 = new NDArray<double, 2>({y.shape(0), y.shape(1)});
// Make views of the numpy arrays
auto y_view = make_view_3d(y);
auto y_view_err = make_view_3d(y_err);
auto x_view = make_view_1d(x);
aare::fit_gaus(x_view, y_view, y_view_err, par->view(),
par_err->view(), chi2->view(), n_threads);
return py::dict("par"_a = return_image_data(par),
"par_err"_a = return_image_data(par_err),
"chi2"_a = return_image_data(chi2),
"Ndf"_a = y.shape(2) - 3);
} else if (y.ndim() == 1) {
// Allocate memory for the output
// Need to have pointers to allow python to manage
// the memory
auto par = new NDArray<double, 1>({3});
auto par_err = new NDArray<double, 1>({3});
// Decode the numpy arrays
auto y_view = make_view_1d(y);
auto y_view_err = make_view_1d(y_err);
auto x_view = make_view_1d(x);
double chi2 = 0;
aare::fit_gaus(x_view, y_view, y_view_err, par->view(),
par_err->view(), chi2);
return py::dict("par"_a = return_image_data(par),
"par_err"_a = return_image_data(par_err),
"chi2"_a = chi2, "Ndf"_a = y.size() - 3);
} else {
throw std::runtime_error("Data must be 1D or 3D");
}
},
R"(
Fit a 1D Gaussian to data with error estimates.
Parameters
----------
x : array_like
The x values.
y : array_like
The y values.
y_err : array_like
The error in the y values.
n_threads : int, optional
The number of threads to use. Default is 4.
)",
py::arg("x"), py::arg("y"), py::arg("y_err"), py::arg("n_threads") = 4);
m.def(
"fit_pol1",
[](py::array_t<double, py::array::c_style | py::array::forcecast> x,
py::array_t<double, py::array::c_style | py::array::forcecast> y,
int n_threads) {
if (y.ndim() == 3) {
auto par = new NDArray<double, 3>{};
auto x_view = make_view_1d(x);
auto y_view = make_view_3d(y);
*par = aare::fit_pol1(x_view, y_view, n_threads);
return return_image_data(par);
} else if (y.ndim() == 1) {
auto par = new NDArray<double, 1>{};
auto x_view = make_view_1d(x);
auto y_view = make_view_1d(y);
*par = aare::fit_pol1(x_view, y_view);
return return_image_data(par);
} else {
throw std::runtime_error("Data must be 1D or 3D");
}
},
py::arg("x"), py::arg("y"), py::arg("n_threads") = 4);
m.def(
"fit_pol1",
[](py::array_t<double, py::array::c_style | py::array::forcecast> x,
py::array_t<double, py::array::c_style | py::array::forcecast> y,
py::array_t<double, py::array::c_style | py::array::forcecast> y_err,
int n_threads) {
if (y.ndim() == 3) {
auto par = new NDArray<double, 3>({y.shape(0), y.shape(1), 2});
auto par_err =
new NDArray<double, 3>({y.shape(0), y.shape(1), 2});
auto y_view = make_view_3d(y);
auto y_view_err = make_view_3d(y_err);
auto x_view = make_view_1d(x);
auto chi2 = new NDArray<double, 2>({y.shape(0), y.shape(1)});
aare::fit_pol1(x_view, y_view, y_view_err, par->view(),
par_err->view(), chi2->view(), n_threads);
return py::dict("par"_a = return_image_data(par),
"par_err"_a = return_image_data(par_err),
"chi2"_a = return_image_data(chi2),
"Ndf"_a = y.shape(2) - 2);
} else if (y.ndim() == 1) {
auto par = new NDArray<double, 1>({2});
auto par_err = new NDArray<double, 1>({2});
auto y_view = make_view_1d(y);
auto y_view_err = make_view_1d(y_err);
auto x_view = make_view_1d(x);
double chi2 = 0;
aare::fit_pol1(x_view, y_view, y_view_err, par->view(),
par_err->view(), chi2);
return py::dict("par"_a = return_image_data(par),
"par_err"_a = return_image_data(par_err),
"chi2"_a = chi2, "Ndf"_a = y.size() - 2);
} else {
throw std::runtime_error("Data must be 1D or 3D");
}
},
R"(
Fit a 1D polynomial to data with error estimates.
Parameters
----------
x : array_like
The x values.
y : array_like
The y values.
y_err : array_like
The error in the y values.
n_threads : int, optional
The number of threads to use. Default is 4.
)",
py::arg("x"), py::arg("y"), py::arg("y_err"), py::arg("n_threads") = 4);
//=========
m.def(
"fit_scurve",
[](py::array_t<double, py::array::c_style | py::array::forcecast> x,
py::array_t<double, py::array::c_style | py::array::forcecast> y,
int n_threads) {
if (y.ndim() == 3) {
auto par = new NDArray<double, 3>{};
auto x_view = make_view_1d(x);
auto y_view = make_view_3d(y);
*par = aare::fit_scurve(x_view, y_view, n_threads);
return return_image_data(par);
} else if (y.ndim() == 1) {
auto par = new NDArray<double, 1>{};
auto x_view = make_view_1d(x);
auto y_view = make_view_1d(y);
*par = aare::fit_scurve(x_view, y_view);
return return_image_data(par);
} else {
throw std::runtime_error("Data must be 1D or 3D");
}
},
py::arg("x"), py::arg("y"), py::arg("n_threads") = 4);
m.def(
"fit_scurve",
[](py::array_t<double, py::array::c_style | py::array::forcecast> x,
py::array_t<double, py::array::c_style | py::array::forcecast> y,
py::array_t<double, py::array::c_style | py::array::forcecast> y_err,
int n_threads) {
if (y.ndim() == 3) {
auto par = new NDArray<double, 3>({y.shape(0), y.shape(1), 6});
auto par_err =
new NDArray<double, 3>({y.shape(0), y.shape(1), 6});
auto y_view = make_view_3d(y);
auto y_view_err = make_view_3d(y_err);
auto x_view = make_view_1d(x);
auto chi2 = new NDArray<double, 2>({y.shape(0), y.shape(1)});
aare::fit_scurve(x_view, y_view, y_view_err, par->view(),
par_err->view(), chi2->view(), n_threads);
return py::dict("par"_a = return_image_data(par),
"par_err"_a = return_image_data(par_err),
"chi2"_a = return_image_data(chi2),
"Ndf"_a = y.shape(2) - 2);
} else if (y.ndim() == 1) {
auto par = new NDArray<double, 1>({2});
auto par_err = new NDArray<double, 1>({2});
auto y_view = make_view_1d(y);
auto y_view_err = make_view_1d(y_err);
auto x_view = make_view_1d(x);
double chi2 = 0;
aare::fit_scurve(x_view, y_view, y_view_err, par->view(),
par_err->view(), chi2);
return py::dict("par"_a = return_image_data(par),
"par_err"_a = return_image_data(par_err),
"chi2"_a = chi2, "Ndf"_a = y.size() - 2);
} else {
throw std::runtime_error("Data must be 1D or 3D");
}
},
R"(
Fit a 1D polynomial to data with error estimates.
Parameters
----------
x : array_like
The x values.
y : array_like
The y values.
y_err : array_like
The error in the y values.
n_threads : int, optional
The number of threads to use. Default is 4.
)",
py::arg("x"), py::arg("y"), py::arg("y_err"), py::arg("n_threads") = 4);
m.def(
"fit_scurve2",
[](py::array_t<double, py::array::c_style | py::array::forcecast> x,
py::array_t<double, py::array::c_style | py::array::forcecast> y,
int n_threads) {
if (y.ndim() == 3) {
auto par = new NDArray<double, 3>{};
auto x_view = make_view_1d(x);
auto y_view = make_view_3d(y);
*par = aare::fit_scurve2(x_view, y_view, n_threads);
return return_image_data(par);
} else if (y.ndim() == 1) {
auto par = new NDArray<double, 1>{};
auto x_view = make_view_1d(x);
auto y_view = make_view_1d(y);
*par = aare::fit_scurve2(x_view, y_view);
return return_image_data(par);
} else {
throw std::runtime_error("Data must be 1D or 3D");
}
},
py::arg("x"), py::arg("y"), py::arg("n_threads") = 4);
m.def(
"fit_scurve2",
[](py::array_t<double, py::array::c_style | py::array::forcecast> x,
py::array_t<double, py::array::c_style | py::array::forcecast> y,
py::array_t<double, py::array::c_style | py::array::forcecast> y_err,
int n_threads) {
if (y.ndim() == 3) {
auto par = new NDArray<double, 3>({y.shape(0), y.shape(1), 6});
auto par_err =
new NDArray<double, 3>({y.shape(0), y.shape(1), 6});
auto y_view = make_view_3d(y);
auto y_view_err = make_view_3d(y_err);
auto x_view = make_view_1d(x);
auto chi2 = new NDArray<double, 2>({y.shape(0), y.shape(1)});
aare::fit_scurve2(x_view, y_view, y_view_err, par->view(),
par_err->view(), chi2->view(), n_threads);
return py::dict("par"_a = return_image_data(par),
"par_err"_a = return_image_data(par_err),
"chi2"_a = return_image_data(chi2),
"Ndf"_a = y.shape(2) - 2);
} else if (y.ndim() == 1) {
auto par = new NDArray<double, 1>({6});
auto par_err = new NDArray<double, 1>({6});
auto y_view = make_view_1d(y);
auto y_view_err = make_view_1d(y_err);
auto x_view = make_view_1d(x);
double chi2 = 0;
aare::fit_scurve2(x_view, y_view, y_view_err, par->view(),
par_err->view(), chi2);
return py::dict("par"_a = return_image_data(par),
"par_err"_a = return_image_data(par_err),
"chi2"_a = chi2, "Ndf"_a = y.size() - 2);
} else {
throw std::runtime_error("Data must be 1D or 3D");
}
},
R"(
Fit a 1D polynomial to data with error estimates.
Parameters
----------
x : array_like
The x values.
y : array_like
The y values.
y_err : array_like
The error in the y values.
n_threads : int, optional
The number of threads to use. Default is 4.
)",
py::arg("x"), py::arg("y"), py::arg("y_err"), py::arg("n_threads") = 4);
}

82
python/src/interpolation.hpp Normal file
View File

@ -0,0 +1,82 @@
#include "aare/Interpolator.hpp"
#include "aare/NDArray.hpp"
#include "aare/NDView.hpp"
#include "np_helper.hpp"
#include <cstdint>
#include <filesystem>
#include <pybind11/pybind11.h>
#include <pybind11/stl.h>
namespace py = pybind11;
template <typename Type, uint8_t CoordSizeX, uint8_t CoordSizeY,
typename CoordType = uint16_t>
void register_interpolate(py::class_<aare::Interpolator> &interpolator) {
using ClusterType = Cluster<Type, CoordSizeX, CoordSizeY, CoordType>;
interpolator.def("interpolate",
[](aare::Interpolator &self,
const ClusterVector<ClusterType> &clusters) {
auto photons = self.interpolate<ClusterType>(clusters);
auto *ptr = new std::vector<Photon>{photons};
return return_vector(ptr);
});
}
void define_interpolation_bindings(py::module &m) {
PYBIND11_NUMPY_DTYPE(aare::Photon, x, y, energy);
auto interpolator =
py::class_<aare::Interpolator>(m, "Interpolator")
.def(py::init([](py::array_t<double, py::array::c_style |
py::array::forcecast>
etacube,
py::array_t<double> xbins,
py::array_t<double> ybins,
py::array_t<double> ebins) {
return Interpolator(make_view_3d(etacube), make_view_1d(xbins),
make_view_1d(ybins), make_view_1d(ebins));
}))
.def("get_ietax",
[](Interpolator &self) {
auto *ptr = new NDArray<double, 3>{};
*ptr = self.get_ietax();
return return_image_data(ptr);
})
.def("get_ietay", [](Interpolator &self) {
auto *ptr = new NDArray<double, 3>{};
*ptr = self.get_ietay();
return return_image_data(ptr);
});
register_interpolate<int, 3, 3, uint16_t>(interpolator);
register_interpolate<float, 3, 3, uint16_t>(interpolator);
register_interpolate<double, 3, 3, uint16_t>(interpolator);
register_interpolate<int, 2, 2, uint16_t>(interpolator);
register_interpolate<float, 2, 2, uint16_t>(interpolator);
register_interpolate<double, 2, 2, uint16_t>(interpolator);
// TODO! Evaluate without converting to double
m.def(
"hej",
[]() {
// auto boost_histogram = py::module_::import("boost_histogram");
// py::object axis =
// boost_histogram.attr("axis").attr("Regular")(10, 0.0, 10.0);
// py::object histogram = boost_histogram.attr("Histogram")(axis);
// return histogram;
// return h;
},
R"(
Evaluate a 1D Gaussian function for all points in x using parameters par.
Parameters
----------
x : array_like
The points at which to evaluate the Gaussian function.
par : array_like
The parameters of the Gaussian function. The first element is the amplitude, the second element is the mean, and the third element is the standard deviation.
)");
}

116
python/src/jungfrau_data_file.hpp Normal file
View File

@ -0,0 +1,116 @@
#include "aare/JungfrauDataFile.hpp"
#include "aare/defs.hpp"
#include <cstdint>
#include <filesystem>
#include <pybind11/iostream.h>
#include <pybind11/numpy.h>
#include <pybind11/pybind11.h>
#include <pybind11/stl.h>
#include <pybind11/stl/filesystem.h>
#include <string>
namespace py = pybind11;
using namespace ::aare;
// Disable warnings for unused parameters, as we ignore some
// in the __exit__ method
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wunused-parameter"
auto read_dat_frame(JungfrauDataFile &self) {
py::array_t<JungfrauDataHeader> header(1);
py::array_t<uint16_t> image({
self.rows(),
self.cols()
});
self.read_into(reinterpret_cast<std::byte *>(image.mutable_data()),
header.mutable_data());
return py::make_tuple(header, image);
}
auto read_n_dat_frames(JungfrauDataFile &self, size_t n_frames) {
// adjust for actual frames left in the file
n_frames = std::min(n_frames, self.total_frames() - self.tell());
if (n_frames == 0) {
throw std::runtime_error("No frames left in file");
}
py::array_t<JungfrauDataHeader> header(n_frames);
py::array_t<uint16_t> image({
n_frames, self.rows(),
self.cols()});
self.read_into(reinterpret_cast<std::byte *>(image.mutable_data()),
n_frames, header.mutable_data());
return py::make_tuple(header, image);
}
void define_jungfrau_data_file_io_bindings(py::module &m) {
// Make the JungfrauDataHeader usable from numpy
PYBIND11_NUMPY_DTYPE(JungfrauDataHeader, framenum, bunchid);
py::class_<JungfrauDataFile>(m, "JungfrauDataFile")
.def(py::init<const std::filesystem::path &>())
.def("seek", &JungfrauDataFile::seek,
R"(
Seek to the given frame index.
)")
.def("tell", &JungfrauDataFile::tell,
R"(
Get the current frame index.
)")
.def_property_readonly("rows", &JungfrauDataFile::rows)
.def_property_readonly("cols", &JungfrauDataFile::cols)
.def_property_readonly("base_name", &JungfrauDataFile::base_name)
.def_property_readonly("bytes_per_frame",
&JungfrauDataFile::bytes_per_frame)
.def_property_readonly("pixels_per_frame",
&JungfrauDataFile::pixels_per_frame)
.def_property_readonly("bytes_per_pixel",
&JungfrauDataFile::bytes_per_pixel)
.def_property_readonly("bitdepth", &JungfrauDataFile::bitdepth)
.def_property_readonly("current_file", &JungfrauDataFile::current_file)
.def_property_readonly("total_frames", &JungfrauDataFile::total_frames)
.def_property_readonly("n_files", &JungfrauDataFile::n_files)
.def("read_frame", &read_dat_frame,
R"(
Read a single frame from the file.
)")
.def("read_n", &read_n_dat_frames,
R"(
Read maximum n_frames frames from the file.
)")
.def(
"read",
[](JungfrauDataFile &self) {
self.seek(0);
auto n_frames = self.total_frames();
return read_n_dat_frames(self, n_frames);
},
R"(
Read all frames from the file. Seeks to the beginning before reading.
)")
.def("__enter__", [](JungfrauDataFile &self) { return &self; })
.def("__exit__",
[](JungfrauDataFile &self,
const std::optional<pybind11::type> &exc_type,
const std::optional<pybind11::object> &exc_value,
const std::optional<pybind11::object> &traceback) {
// self.close();
})
.def("__iter__", [](JungfrauDataFile &self) { return &self; })
.def("__next__", [](JungfrauDataFile &self) {
try {
return read_dat_frame(self);
} catch (std::runtime_error &e) {
throw py::stop_iteration();
}
});
}
#pragma GCC diagnostic pop

95
python/src/module.cpp
View File

@ -1,19 +1,32 @@
//Files with bindings to the different classes
#include "file.hpp"
#include "raw_file.hpp"
// Files with bindings to the different classes
//New style file naming
#include "bind_Cluster.hpp"
#include "bind_ClusterCollector.hpp"
#include "bind_ClusterFinder.hpp"
#include "bind_ClusterFinderMT.hpp"
#include "bind_ClusterFile.hpp"
#include "bind_ClusterFileSink.hpp"
#include "bind_ClusterVector.hpp"
//TODO! migrate the other names
#include "ctb_raw_file.hpp"
#include "file.hpp"
#include "fit.hpp"
#include "interpolation.hpp"
#include "raw_sub_file.hpp"
#include "raw_master_file.hpp"
#include "raw_file.hpp"
#ifdef HDF5_FOUND
#include "hdf5_file.hpp"
#include "hdf5_master_file.hpp"
#endif
#include "var_cluster.hpp"
#include "pixel_map.hpp"
#include "var_cluster.hpp"
#include "pedestal.hpp"
#include "cluster.hpp"
#include "cluster_file.hpp"
#include "jungfrau_data_file.hpp"
//Pybind stuff
// Pybind stuff
#include <pybind11/pybind11.h>
#include <pybind11/stl.h>
@ -22,6 +35,7 @@ namespace py = pybind11;
PYBIND11_MODULE(_aare, m) {
define_file_io_bindings(m);
define_raw_file_io_bindings(m);
define_raw_sub_file_io_bindings(m);
define_ctb_raw_file_io_bindings(m);
define_raw_master_file_bindings(m);
#ifdef HDF5_FOUND
@ -30,8 +44,65 @@ PYBIND11_MODULE(_aare, m) {
#endif
define_var_cluster_finder_bindings(m);
define_pixel_map_bindings(m);
define_pedestal_bindings<double>(m, "Pedestal");
define_pedestal_bindings<float>(m, "Pedestal_float32");
define_cluster_finder_bindings(m);
define_cluster_file_io_bindings(m);
}
define_pedestal_bindings<double>(m, "Pedestal_d");
define_pedestal_bindings<float>(m, "Pedestal_f");
define_fit_bindings(m);
define_interpolation_bindings(m);
define_jungfrau_data_file_io_bindings(m);
define_ClusterFile<int, 3, 3, uint16_t>(m, "Cluster3x3i");
define_ClusterFile<double, 3, 3, uint16_t>(m, "Cluster3x3d");
define_ClusterFile<float, 3, 3, uint16_t>(m, "Cluster3x3f");
define_ClusterFile<int, 2, 2, uint16_t>(m, "Cluster2x2i");
define_ClusterFile<float, 2, 2, uint16_t>(m, "Cluster2x2f");
define_ClusterFile<double, 2, 2, uint16_t>(m, "Cluster2x2d");
define_ClusterVector<int, 3, 3, uint16_t>(m, "Cluster3x3i");
define_ClusterVector<double, 3, 3, uint16_t>(m, "Cluster3x3d");
define_ClusterVector<float, 3, 3, uint16_t>(m, "Cluster3x3f");
define_ClusterVector<int, 2, 2, uint16_t>(m, "Cluster2x2i");
define_ClusterVector<double, 2, 2, uint16_t>(m, "Cluster2x2d");
define_ClusterVector<float, 2, 2, uint16_t>(m, "Cluster2x2f");
define_ClusterFinder<int, 3, 3, uint16_t>(m, "Cluster3x3i");
define_ClusterFinder<double, 3, 3, uint16_t>(m, "Cluster3x3d");
define_ClusterFinder<float, 3, 3, uint16_t>(m, "Cluster3x3f");
define_ClusterFinder<int, 2, 2, uint16_t>(m, "Cluster2x2i");
define_ClusterFinder<double, 2, 2, uint16_t>(m, "Cluster2x2d");
define_ClusterFinder<float, 2, 2, uint16_t>(m, "Cluster2x2f");
define_ClusterFinderMT<int, 3, 3, uint16_t>(m, "Cluster3x3i");
define_ClusterFinderMT<double, 3, 3, uint16_t>(m, "Cluster3x3d");
define_ClusterFinderMT<float, 3, 3, uint16_t>(m, "Cluster3x3f");
define_ClusterFinderMT<int, 2, 2, uint16_t>(m, "Cluster2x2i");
define_ClusterFinderMT<double, 2, 2, uint16_t>(m, "Cluster2x2d");
define_ClusterFinderMT<float, 2, 2, uint16_t>(m, "Cluster2x2f");
define_ClusterFileSink<int, 3, 3, uint16_t>(m, "Cluster3x3i");
define_ClusterFileSink<double, 3, 3, uint16_t>(m, "Cluster3x3d");
define_ClusterFileSink<float, 3, 3, uint16_t>(m, "Cluster3x3f");
define_ClusterFileSink<int, 2, 2, uint16_t>(m, "Cluster2x2i");
define_ClusterFileSink<double, 2, 2, uint16_t>(m, "Cluster2x2d");
define_ClusterFileSink<float, 2, 2, uint16_t>(m, "Cluster2x2f");
define_ClusterCollector<int, 3, 3, uint16_t>(m, "Cluster3x3i");
define_ClusterCollector<double, 3, 3, uint16_t>(m, "Cluster3x3d");
define_ClusterCollector<float, 3, 3, uint16_t>(m, "Cluster3x3f");
define_ClusterCollector<int, 2, 2, uint16_t>(m, "Cluster2x2i");
define_ClusterCollector<double, 2, 2, uint16_t>(m, "Cluster2x2d");
define_ClusterCollector<float, 2, 2, uint16_t>(m, "Cluster2x2f");
define_Cluster<int, 3, 3, uint16_t>(m, "3x3i");
define_Cluster<float, 3, 3, uint16_t>(m, "3x3f");
define_Cluster<double, 3, 3, uint16_t>(m, "3x3d");
define_Cluster<int, 2, 2, uint16_t>(m, "2x2i");
define_Cluster<float, 2, 2, uint16_t>(m, "2x2f");
define_Cluster<double, 2, 2, uint16_t>(m, "2x2d");
register_calculate_eta<int, 3, 3, uint16_t>(m);
register_calculate_eta<float, 3, 3, uint16_t>(m);
register_calculate_eta<double, 3, 3, uint16_t>(m);
register_calculate_eta<int, 2, 2, uint16_t>(m);
register_calculate_eta<float, 2, 2, uint16_t>(m);
register_calculate_eta<double, 2, 2, uint16_t>(m);
}

100
python/src/np_helper.hpp
View File

@ -10,9 +10,10 @@
#include "aare/NDView.hpp"
namespace py = pybind11;
using namespace aare;
// Pass image data back to python as a numpy array
template <typename T, int64_t Ndim>
template <typename T, ssize_t Ndim>
py::array return_image_data(aare::NDArray<T, Ndim> *image) {
py::capsule free_when_done(image, [](void *f) {
@ -39,78 +40,47 @@ template <typename T> py::array return_vector(std::vector<T> *vec) {
free_when_done); // numpy array references this parent
}
// template <typename Reader> py::array do_read(Reader &r, size_t n_frames) {
// py::array image;
// if (n_frames == 0)
// n_frames = r.total_frames();
// std::array<ssize_t, 3> shape{static_cast<ssize_t>(n_frames), r.rows(),
// r.cols()};
// const uint8_t item_size = r.bytes_per_pixel();
// if (item_size == 1) {
// image = py::array_t<uint8_t, py::array::c_style | py::array::forcecast>(
// shape);
// } else if (item_size == 2) {
// image =
// py::array_t<uint16_t, py::array::c_style | py::array::forcecast>(
// shape);
// } else if (item_size == 4) {
// image =
// py::array_t<uint32_t, py::array::c_style | py::array::forcecast>(
// shape);
// }
// r.read_into(reinterpret_cast<std::byte *>(image.mutable_data()), n_frames);
// return image;
// }
// py::array return_frame(pl::Frame *ptr) {
// py::capsule free_when_done(ptr, [](void *f) {
// pl::Frame *foo = reinterpret_cast<pl::Frame *>(f);
// delete foo;
// });
// const uint8_t item_size = ptr->bytes_per_pixel();
// std::vector<ssize_t> shape;
// for (auto val : ptr->shape())
// if (val > 1)
// shape.push_back(val);
// std::vector<ssize_t> strides;
// if (shape.size() == 1)
// strides.push_back(item_size);
// else if (shape.size() == 2) {
// strides.push_back(item_size * shape[1]);
// strides.push_back(item_size);
// }
// if (item_size == 1)
// return py::array_t<uint8_t>(
// shape, strides,
// reinterpret_cast<uint8_t *>(ptr->data()), free_when_done);
// else if (item_size == 2)
// return py::array_t<uint16_t>(shape, strides,
// reinterpret_cast<uint16_t *>(ptr->data()),
// free_when_done);
// else if (item_size == 4)
// return py::array_t<uint32_t>(shape, strides,
// reinterpret_cast<uint32_t *>(ptr->data()),
// free_when_done);
// return {};
// }
// todo rewrite generic
template <class T, int Flags> auto get_shape_3d(py::array_t<T, Flags> arr) {
template <class T, int Flags>
auto get_shape_3d(const py::array_t<T, Flags> &arr) {
return aare::Shape<3>{arr.shape(0), arr.shape(1), arr.shape(2)};
}
template <class T, int Flags> auto make_view_3d(py::array_t<T, Flags> arr) {
template <class T, int Flags> auto make_view_3d(py::array_t<T, Flags> &arr) {
return aare::NDView<T, 3>(arr.mutable_data(), get_shape_3d<T, Flags>(arr));
}
template <class T, int Flags> auto get_shape_2d(py::array_t<T, Flags> arr) {
template <class T, int Flags>
auto get_shape_2d(const py::array_t<T, Flags> &arr) {
return aare::Shape<2>{arr.shape(0), arr.shape(1)};
}
template <class T, int Flags> auto make_view_2d(py::array_t<T, Flags> arr) {
template <class T, int Flags>
auto get_shape_1d(const py::array_t<T, Flags> &arr) {
return aare::Shape<1>{arr.shape(0)};
}
template <class T, int Flags> auto make_view_2d(py::array_t<T, Flags> &arr) {
return aare::NDView<T, 2>(arr.mutable_data(), get_shape_2d<T, Flags>(arr));
}
}
template <class T, int Flags> auto make_view_1d(py::array_t<T, Flags> &arr) {
return aare::NDView<T, 1>(arr.mutable_data(), get_shape_1d<T, Flags>(arr));
}
template <typename ClusterType> struct fmt_format_trait; // forward declaration
template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
typename CoordType>
struct fmt_format_trait<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> {
static std::string value() {
return fmt::format("T{{{}:x:{}:y:{}:data:}}",
py::format_descriptor<CoordType>::format(),
py::format_descriptor<CoordType>::format(),
fmt::format("({},{}){}", ClusterSizeX, ClusterSizeY,
py::format_descriptor<T>::format()));
}
};
template <typename ClusterType>
auto fmt_format = fmt_format_trait<ClusterType>::value();

7
python/src/pedestal.hpp
View File

@ -43,5 +43,10 @@ template <typename SUM_TYPE> void define_pedestal_bindings(py::module &m, const
.def("push", [](Pedestal<SUM_TYPE> &pedestal, py::array_t<uint16_t> &f) {
auto v = make_view_2d(f);
pedestal.push(v);
});
})
.def("push_no_update", [](Pedestal<SUM_TYPE> &pedestal, py::array_t<uint16_t, py::array::c_style> &f) {
auto v = make_view_2d(f);
pedestal.push_no_update(v);
}, py::arg().noconvert())
.def("update_mean", &Pedestal<SUM_TYPE>::update_mean);
}

8
python/src/raw_file.hpp
View File

@ -32,7 +32,7 @@ void define_raw_file_io_bindings(py::module &m) {
shape.push_back(self.cols());
// return headers from all subfiles
py::array_t<DetectorHeader> header(self.n_mod());
py::array_t<DetectorHeader> header(self.n_modules());
const uint8_t item_size = self.bytes_per_pixel();
if (item_size == 1) {
@ -61,10 +61,10 @@ void define_raw_file_io_bindings(py::module &m) {
// return headers from all subfiles
py::array_t<DetectorHeader> header;
if (self.n_mod() == 1) {
if (self.n_modules() == 1) {
header = py::array_t<DetectorHeader>(n_frames);
} else {
header = py::array_t<DetectorHeader>({self.n_mod(), n_frames});
header = py::array_t<DetectorHeader>({self.n_modules(), n_frames});
}
// py::array_t<DetectorHeader> header({self.n_mod(), n_frames});
@ -100,7 +100,7 @@ void define_raw_file_io_bindings(py::module &m) {
.def_property_readonly("cols", &RawFile::cols)
.def_property_readonly("bitdepth", &RawFile::bitdepth)
.def_property_readonly("geometry", &RawFile::geometry)
.def_property_readonly("n_mod", &RawFile::n_mod)
.def_property_readonly("n_modules", &RawFile::n_modules)
.def_property_readonly("detector_type", &RawFile::detector_type)
.def_property_readonly("master", &RawFile::master);
}

110
python/src/raw_sub_file.hpp Normal file
View File

@ -0,0 +1,110 @@
#include "aare/CtbRawFile.hpp"
#include "aare/File.hpp"
#include "aare/Frame.hpp"
#include "aare/RawFile.hpp"
#include "aare/RawMasterFile.hpp"
#include "aare/RawSubFile.hpp"
#include "aare/defs.hpp"
// #include "aare/fClusterFileV2.hpp"
#include <cstdint>
#include <filesystem>
#include <pybind11/iostream.h>
#include <pybind11/numpy.h>
#include <pybind11/pybind11.h>
#include <pybind11/stl.h>
#include <pybind11/stl/filesystem.h>
#include <string>
namespace py = pybind11;
using namespace ::aare;
auto read_frame_from_RawSubFile(RawSubFile &self) {
py::array_t<DetectorHeader> header(1);
const uint8_t item_size = self.bytes_per_pixel();
std::vector<ssize_t> shape{static_cast<ssize_t>(self.rows()),
static_cast<ssize_t>(self.cols())};
py::array image;
if (item_size == 1) {
image = py::array_t<uint8_t>(shape);
} else if (item_size == 2) {
image = py::array_t<uint16_t>(shape);
} else if (item_size == 4) {
image = py::array_t<uint32_t>(shape);
}
self.read_into(reinterpret_cast<std::byte *>(image.mutable_data()),
header.mutable_data());
return py::make_tuple(header, image);
}
auto read_n_frames_from_RawSubFile(RawSubFile &self, size_t n_frames) {
py::array_t<DetectorHeader> header(n_frames);
const uint8_t item_size = self.bytes_per_pixel();
std::vector<ssize_t> shape{
static_cast<ssize_t>(n_frames),
static_cast<ssize_t>(self.rows()),
static_cast<ssize_t>(self.cols())
};
py::array image;
if (item_size == 1) {
image = py::array_t<uint8_t>(shape);
} else if (item_size == 2) {
image = py::array_t<uint16_t>(shape);
} else if (item_size == 4) {
image = py::array_t<uint32_t>(shape);
}
self.read_into(reinterpret_cast<std::byte *>(image.mutable_data()), n_frames,
header.mutable_data());
return py::make_tuple(header, image);
}
//Disable warnings for unused parameters, as we ignore some
//in the __exit__ method
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wunused-parameter"
void define_raw_sub_file_io_bindings(py::module &m) {
py::class_<RawSubFile>(m, "RawSubFile")
.def(py::init<const std::filesystem::path &, DetectorType, size_t,
size_t, size_t>())
.def_property_readonly("bytes_per_frame", &RawSubFile::bytes_per_frame)
.def_property_readonly("pixels_per_frame",
&RawSubFile::pixels_per_frame)
.def_property_readonly("bytes_per_pixel", &RawSubFile::bytes_per_pixel)
.def("seek", &RawSubFile::seek)
.def("tell", &RawSubFile::tell)
.def_property_readonly("rows", &RawSubFile::rows)
.def_property_readonly("cols", &RawSubFile::cols)
.def_property_readonly("frames_in_file", &RawSubFile::frames_in_file)
.def("read_frame", &read_frame_from_RawSubFile)
.def("read_n", &read_n_frames_from_RawSubFile)
.def("read", [](RawSubFile &self){
self.seek(0);
auto n_frames = self.frames_in_file();
return read_n_frames_from_RawSubFile(self, n_frames);
})
.def("__enter__", [](RawSubFile &self) { return &self; })
.def("__exit__",
[](RawSubFile &self,
const std::optional<pybind11::type> &exc_type,
const std::optional<pybind11::object> &exc_value,
const std::optional<pybind11::object> &traceback) {
})
.def("__iter__", [](RawSubFile &self) { return &self; })
.def("__next__", [](RawSubFile &self) {
try {
return read_frame_from_RawSubFile(self);
} catch (std::runtime_error &e) {
throw py::stop_iteration();
}
});
}
#pragma GCC diagnostic pop

37
python/src/var_cluster.hpp
View File

@ -19,15 +19,24 @@ using namespace::aare;
void define_var_cluster_finder_bindings(py::module &m) {
PYBIND11_NUMPY_DTYPE(VarClusterFinder<double>::Hit, size, row, col,
reserved, energy, max);
reserved, energy, max, rows, cols, enes);
py::class_<VarClusterFinder<double>>(m, "VarClusterFinder")
.def(py::init<Shape<2>, double>())
.def("labeled",
[](VarClusterFinder<double> &self) {
auto ptr = new NDArray<int, 2>(self.labeled());
auto *ptr = new NDArray<int, 2>(self.labeled());
return return_image_data(ptr);
})
.def("set_noiseMap",
[](VarClusterFinder<double> &self,
py::array_t<double, py::array::c_style | py::array::forcecast>
noise_map) {
auto noise_map_span = make_view_2d(noise_map);
self.set_noiseMap(noise_map_span);
})
.def("set_peripheralThresholdFactor",
&VarClusterFinder<double>::set_peripheralThresholdFactor)
.def("find_clusters",
[](VarClusterFinder<double> &self,
py::array_t<double, py::array::c_style | py::array::forcecast>
@ -35,6 +44,30 @@ void define_var_cluster_finder_bindings(py::module &m) {
auto view = make_view_2d(img);
self.find_clusters(view);
})
.def("find_clusters_X",
[](VarClusterFinder<double> &self,
py::array_t<double, py::array::c_style | py::array::forcecast>
img) {
auto img_span = make_view_2d(img);
self.find_clusters_X(img_span);
})
.def("single_pass",
[](VarClusterFinder<double> &self,
py::array_t<double, py::array::c_style | py::array::forcecast>
img) {
auto img_span = make_view_2d(img);
self.single_pass(img_span);
})
.def("hits",
[](VarClusterFinder<double> &self) {
auto ptr = new std::vector<VarClusterFinder<double>::Hit>(
self.steal_hits());
return return_vector(ptr);
})
.def("clear_hits",
[](VarClusterFinder<double> &self) {
self.clear_hits();
})
.def("steal_hits",
[](VarClusterFinder<double> &self) {
auto ptr = new std::vector<VarClusterFinder<double>::Hit>(

34
python/tests/conftest.py Normal file
View File

@ -0,0 +1,34 @@
import os
from pathlib import Path
import pytest
def pytest_addoption(parser):
parser.addoption(
"--files", action="store_true", default=False, help="run slow tests"
)
def pytest_configure(config):
config.addinivalue_line("markers", "files: mark test as needing image files to run")
def pytest_collection_modifyitems(config, items):
if config.getoption("--files"):
return
skip = pytest.mark.skip(reason="need --files option to run")
for item in items:
if "files" in item.keywords:
item.add_marker(skip)
@pytest.fixture
def test_data_path():
env_value = os.environ.get("AARE_TEST_DATA")
if not env_value:
raise RuntimeError("Environment variable AARE_TEST_DATA is not set or is empty")
return Path(env_value)

110
python/tests/test_Cluster.py Normal file
View File

@ -0,0 +1,110 @@
import pytest
import numpy as np
from aare import _aare #import the C++ module
from conftest import test_data_path
def test_cluster_vector_can_be_converted_to_numpy():
cv = _aare.ClusterVector_Cluster3x3i()
arr = np.array(cv, copy=False)
assert arr.shape == (0,) # 4 for x, y, size, energy and 9 for the cluster data
def test_ClusterVector():
"""Test ClusterVector"""
clustervector = _aare.ClusterVector_Cluster3x3i()
assert clustervector.cluster_size_x == 3
assert clustervector.cluster_size_y == 3
assert clustervector.item_size() == 4+9*4
assert clustervector.frame_number == 0
assert clustervector.size == 0
cluster = _aare.Cluster3x3i(0,0,np.ones(9, dtype=np.int32))
clustervector.push_back(cluster)
assert clustervector.size == 1
with pytest.raises(TypeError): # Or use the appropriate exception type
clustervector.push_back(_aare.Cluster2x2i(0,0,np.ones(4, dtype=np.int32)))
with pytest.raises(TypeError):
clustervector.push_back(_aare.Cluster3x3f(0,0,np.ones(9, dtype=np.float32)))
def test_Interpolator():
"""Test Interpolator"""
ebins = np.linspace(0,10, 20, dtype=np.float64)
xbins = np.linspace(0, 5, 30, dtype=np.float64)
ybins = np.linspace(0, 5, 30, dtype=np.float64)
etacube = np.zeros(shape=[30, 30, 20], dtype=np.float64)
interpolator = _aare.Interpolator(etacube, xbins, ybins, ebins)
assert interpolator.get_ietax().shape == (30,30,20)
assert interpolator.get_ietay().shape == (30,30,20)
clustervector = _aare.ClusterVector_Cluster3x3i()
cluster = _aare.Cluster3x3i(0,0, np.ones(9, dtype=np.int32))
clustervector.push_back(cluster)
interpolated_photons = interpolator.interpolate(clustervector)
assert interpolated_photons.size == 1
assert interpolated_photons[0]["x"] == -1
assert interpolated_photons[0]["y"] == -1
assert interpolated_photons[0]["energy"] == 4 #eta_sum = 4, dx, dy = -1,-1 m_ietax = 0, m_ietay = 0
clustervector = _aare.ClusterVector_Cluster2x2i()
cluster = _aare.Cluster2x2i(0,0, np.ones(4, dtype=np.int32))
clustervector.push_back(cluster)
interpolated_photons = interpolator.interpolate(clustervector)
assert interpolated_photons.size == 1
assert interpolated_photons[0]["x"] == 0
assert interpolated_photons[0]["y"] == 0
assert interpolated_photons[0]["energy"] == 4
def test_calculate_eta():
"""Calculate Eta"""
clusters = _aare.ClusterVector_Cluster3x3i()
clusters.push_back(_aare.Cluster3x3i(0,0, np.ones(9, dtype=np.int32)))
clusters.push_back(_aare.Cluster3x3i(0,0, np.array([1,1,1,2,2,2,3,3,3])))
eta2 = _aare.calculate_eta2(clusters)
assert eta2.shape == (2,2)
assert eta2[0,0] == 0.5
assert eta2[0,1] == 0.5
assert eta2[1,0] == 0.5
assert eta2[1,1] == 0.6 #1/5
def test_cluster_finder():
"""Test ClusterFinder"""
clusterfinder = _aare.ClusterFinder_Cluster3x3i([100,100])
#frame = np.random.rand(100,100)
frame = np.zeros(shape=[100,100])
clusterfinder.find_clusters(frame)
clusters = clusterfinder.steal_clusters(False) #conversion does not work
assert clusters.size == 0

64
python/tests/test_ClusterFile.py Normal file
View File

@ -0,0 +1,64 @@
import pytest
import numpy as np
import boost_histogram as bh
import time
from pathlib import Path
import pickle
from aare import ClusterFile
from conftest import test_data_path
@pytest.mark.files
def test_cluster_file(test_data_path):
"""Test ClusterFile"""
f = ClusterFile(test_data_path / "clust/single_frame_97_clustrers.clust")
cv = f.read_clusters(10) #conversion does not work
assert cv.frame_number == 135
assert cv.size == 10
#Known data
#frame_number, num_clusters [135] 97
#[ 1 200] [0 1 2 3 4 5 6 7 8]
#[ 2 201] [ 9 10 11 12 13 14 15 16 17]
#[ 3 202] [18 19 20 21 22 23 24 25 26]
#[ 4 203] [27 28 29 30 31 32 33 34 35]
#[ 5 204] [36 37 38 39 40 41 42 43 44]
#[ 6 205] [45 46 47 48 49 50 51 52 53]
#[ 7 206] [54 55 56 57 58 59 60 61 62]
#[ 8 207] [63 64 65 66 67 68 69 70 71]
#[ 9 208] [72 73 74 75 76 77 78 79 80]
#[ 10 209] [81 82 83 84 85 86 87 88 89]
#conversion to numpy array
arr = np.array(cv, copy = False)
assert arr.size == 10
for i in range(10):
assert arr[i]['x'] == i+1
@pytest.mark.files
def test_read_clusters_and_fill_histogram(test_data_path):
# Create the histogram
n_bins = 100
xmin = -100
xmax = 1e4
hist_aare = bh.Histogram(bh.axis.Regular(n_bins, xmin, xmax))
fname = test_data_path / "clust/beam_En700eV_-40deg_300V_10us_d0_f0_100.clust"
#Read clusters and fill the histogram with pixel values
with ClusterFile(fname, chunk_size = 10000) as f:
for clusters in f:
arr = np.array(clusters, copy = False)
hist_aare.fill(arr['data'].flat)
#Load the histogram from the pickle file
with open(fname.with_suffix('.pkl'), 'rb') as f:
hist_py = pickle.load(f)
#Compare the two histograms
assert hist_aare == hist_py

54
python/tests/test_ClusterVector.py Normal file
View File

@ -0,0 +1,54 @@
import pytest
import numpy as np
import boost_histogram as bh
import time
from pathlib import Path
import pickle
from aare import ClusterFile
from aare import _aare
from conftest import test_data_path
def test_create_cluster_vector():
cv = _aare.ClusterVector_Cluster3x3i()
assert cv.cluster_size_x == 3
assert cv.cluster_size_y == 3
assert cv.size == 0
def test_push_back_on_cluster_vector():
cv = _aare.ClusterVector_Cluster2x2i()
assert cv.cluster_size_x == 2
assert cv.cluster_size_y == 2
assert cv.size == 0
cluster = _aare.Cluster2x2i(19, 22, np.ones(4, dtype=np.int32))
cv.push_back(cluster)
assert cv.size == 1
arr = np.array(cv, copy=False)
assert arr[0]['x'] == 19
assert arr[0]['y'] == 22
def test_make_a_hitmap_from_cluster_vector():
cv = _aare.ClusterVector_Cluster3x3i()
# Push back 4 clusters with different positions
cv.push_back(_aare.Cluster3x3i(0, 0, np.ones(9, dtype=np.int32)))
cv.push_back(_aare.Cluster3x3i(1, 1, np.ones(9, dtype=np.int32)))
cv.push_back(_aare.Cluster3x3i(1, 1, np.ones(9, dtype=np.int32)))
cv.push_back(_aare.Cluster3x3i(2, 2, np.ones(9, dtype=np.int32)))
ref = np.zeros((5, 5), dtype=np.int32)
ref[0,0] = 1
ref[1,1] = 2
ref[2,2] = 1
img = _aare.hitmap((5,5), cv)
# print(img)
# print(ref)
assert (img == ref).all()

39
python/tests/test_RawSubFile.py Normal file
View File

@ -0,0 +1,39 @@
import pytest
import numpy as np
from aare import RawSubFile, DetectorType
@pytest.mark.files
def test_read_a_jungfrau_RawSubFile(test_data_path):
# Starting with f1 there is now 7 frames left in the series of files
with RawSubFile(test_data_path / "raw/jungfrau/jungfrau_single_d0_f1_0.raw", DetectorType.Jungfrau, 512, 1024, 16) as f:
assert f.frames_in_file == 7
headers, frames = f.read()
assert headers.size == 7
assert frames.shape == (7, 512, 1024)
for i,h in zip(range(4,11,1), headers):
assert h["frameNumber"] == i
# Compare to canned data using numpy
data = np.load(test_data_path / "raw/jungfrau/jungfrau_single_0.npy")
assert np.all(data[3:] == frames)
@pytest.mark.files
def test_iterate_over_a_jungfrau_RawSubFile(test_data_path):
data = np.load(test_data_path / "raw/jungfrau/jungfrau_single_0.npy")
# Given the first subfile in a series we can read all frames from f0, f1, f2...fN
with RawSubFile(test_data_path / "raw/jungfrau/jungfrau_single_d0_f0_0.raw", DetectorType.Jungfrau, 512, 1024, 16) as f:
i = 0
for header, frame in f:
assert header["frameNumber"] == i+1
assert np.all(frame == data[i])
i += 1
assert i == 10
assert header["frameNumber"] == 10

92
python/tests/test_jungfrau_dat_files.py Normal file
View File

@ -0,0 +1,92 @@
import pytest
import numpy as np
from aare import JungfrauDataFile
@pytest.mark.files
def test_jfungfrau_dat_read_number_of_frames(test_data_path):
with JungfrauDataFile(test_data_path / "dat/AldoJF500k_000000.dat") as dat_file:
assert dat_file.total_frames == 24
with JungfrauDataFile(test_data_path / "dat/AldoJF250k_000000.dat") as dat_file:
assert dat_file.total_frames == 53
with JungfrauDataFile(test_data_path / "dat/AldoJF65k_000000.dat") as dat_file:
assert dat_file.total_frames == 113
@pytest.mark.files
def test_jfungfrau_dat_read_number_of_file(test_data_path):
with JungfrauDataFile(test_data_path / "dat/AldoJF500k_000000.dat") as dat_file:
assert dat_file.n_files == 4
with JungfrauDataFile(test_data_path / "dat/AldoJF250k_000000.dat") as dat_file:
assert dat_file.n_files == 7
with JungfrauDataFile(test_data_path / "dat/AldoJF65k_000000.dat") as dat_file:
assert dat_file.n_files == 7
@pytest.mark.files
def test_read_module(test_data_path):
"""
Read all frames from the series of .dat files. Compare to canned data in npz format.
"""
# Read all frames from the .dat file
with JungfrauDataFile(test_data_path / "dat/AldoJF500k_000000.dat") as f:
header, data = f.read()
#Sanity check
n_frames = 24
assert header.size == n_frames
assert data.shape == (n_frames, 512, 1024)
# Read reference data using numpy
with np.load(test_data_path / "dat/AldoJF500k.npz") as f:
ref_header = f["headers"]
ref_data = f["frames"]
# Check that the data is the same
assert np.all(ref_header == header)
assert np.all(ref_data == data)
@pytest.mark.files
def test_read_half_module(test_data_path):
# Read all frames from the .dat file
with JungfrauDataFile(test_data_path / "dat/AldoJF250k_000000.dat") as f:
header, data = f.read()
n_frames = 53
assert header.size == n_frames
assert data.shape == (n_frames, 256, 1024)
# Read reference data using numpy
with np.load(test_data_path / "dat/AldoJF250k.npz") as f:
ref_header = f["headers"]
ref_data = f["frames"]
# Check that the data is the same
assert np.all(ref_header == header)
assert np.all(ref_data == data)
@pytest.mark.files
def test_read_single_chip(test_data_path):
# Read all frames from the .dat file
with JungfrauDataFile(test_data_path / "dat/AldoJF65k_000000.dat") as f:
header, data = f.read()
n_frames = 113
assert header.size == n_frames
assert data.shape == (n_frames, 256, 256)
# Read reference data using numpy
with np.load(test_data_path / "dat/AldoJF65k.npz") as f:
ref_header = f["headers"]
ref_data = f["frames"]
# Check that the data is the same
assert np.all(ref_header == header)
assert np.all(ref_data == data)

Some files were not shown because too many files have changed in this diff Show More