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

For 2025.5.22 release (#181)
All checks were successful
Build on RHEL9 / build (push) Successful in 2m22s
Details
Build on RHEL8 / build (push) Successful in 2m29s
Details

Browse Source 
Co-authored-by: Patrick <patrick.sieberer@psi.ch>
Co-authored-by: JulianHeymes <julian.heymes@psi.ch>
Co-authored-by: Dhanya Thattil <dhanya.thattil@psi.ch>
Co-authored-by: Xiangyu Xie <45243914+xiangyuxie@users.noreply.github.com>
Co-authored-by: xiangyu.xie <xiangyu.xie@psi.ch>
Co-authored-by: AliceMazzoleni99 <alice.mazzoleni@psi.ch>
Co-authored-by: Mazzoleni Alice Francesca <mazzol_a@pc17378.psi.ch>
Co-authored-by: siebsi <sieb.patr@gmail.com>
This commit is contained in:
Erik Fröjdh 2025-05-22 11:40:39 +02:00 committed by GitHub
parent fd0196f2fd
commit 94ac58b09e
No known key found for this signature in database
GPG Key ID: B5690EEEBB952194
78 changed files with 3865 additions and 1461 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

14
.github/workflows/build_and_deploy_conda.yml vendored
View File

@ -1,9 +1,9 @@
name: Build pkgs and deploy if on main
on:
push:
branches:
- main
release:
types:
- published
jobs:
build:
@ -24,13 +24,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 conda-build=24.9 conda-verify pytest anaconda-client
conda-remove-defaults: "true"
- name: Enable upload
run: conda config --set anaconda_upload yes

9
.github/workflows/build_conda.yml vendored
View File

@ -24,14 +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
conda-remove-defaults: "true"
- name: Prepare
run: conda install conda-build=24.9 conda-verify pytest anaconda-client
- name: Disable upload
run: conda config --set anaconda_upload no

32
CMakeLists.txt
View File

@ -1,12 +1,17 @@
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)
@ -39,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)
@ -74,6 +79,9 @@ 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)
@ -85,6 +93,7 @@ if(AARE_BENCHMARKS)
endif()
set(CMAKE_EXPORT_COMPILE_COMMANDS ON)
if(AARE_FETCH_LMFIT)
@ -340,6 +349,8 @@ 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
@ -352,6 +363,7 @@ set(PUBLICHEADERS
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
@ -365,13 +377,11 @@ set(PUBLICHEADERS
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
@ -388,19 +398,18 @@ set(SourceFiles
${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
)
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
@ -410,7 +419,8 @@ target_link_libraries(
${STD_FS_LIB} # from helpers.cmake
PRIVATE
aare_compiler_flags
"$<BUILD_INTERFACE:lmfit>"
Threads::Threads
$<BUILD_INTERFACE:lmfit>
)
@ -436,12 +446,16 @@ if(AARE_TESTS)
${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
)

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

51
conda-recipe/meta.yaml
View File

@ -1,11 +1,10 @@
source:
path: ../
{% set version = load_file_regex(load_file = 'VERSION', regex_pattern = '(\d+(?:\.\d+)*(?:[\+\w\.]+))').group(1) %}
package:
name: aare
version: 2025.4.22 #TODO! how to not duplicate this?
version: {{version}}
source:
path: ..
@ -13,45 +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

6
etc/dev-env.yml
View File

@ -3,13 +3,11 @@ channels:
- conda-forge
dependencies:
- anaconda-client
- conda-build
- doxygen
- sphinx=7.1.2
- breathe
- pybind11
- sphinx_rtd_theme
- furo
- nlohmann_json
- zeromq
- fmt
- numpy

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

94
include/aare/Cluster.hpp
View File

@ -1,36 +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 <cstddef>
#include <cstdint>
#include <numeric>
#include <type_traits>
namespace aare {
//TODO! Template this?
struct Cluster3x3 {
int16_t x;
int16_t y;
int32_t data[9];
// requires clause c++20 maybe update
template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
typename CoordType = int16_t>
struct Cluster {
int32_t sum_2x2() const{
std::array<int32_t, 4> total;
total[0] = data[0] + data[1] + data[3] + data[4];
total[1] = data[1] + data[2] + data[4] + data[5];
total[2] = data[3] + data[4] + data[6] + data[7];
total[3] = data[4] + data[5] + data[7] + data[8];
return *std::max_element(total.begin(), total.end());
}
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");
int32_t sum() const{
return std::accumulate(data, data + 9, 0);
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);
}
}
};
struct Cluster2x2 {
int16_t x;
int16_t y;
int32_t data[4];
};
} // namespace aare
// 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

60
include/aare/ClusterCollector.hpp
View File

@ -2,29 +2,31 @@
#include <atomic>
#include <thread>
#include "aare/ProducerConsumerQueue.hpp"
#include "aare/ClusterVector.hpp"
#include "aare/ClusterFinderMT.hpp"
#include "aare/ClusterVector.hpp"
#include "aare/ProducerConsumerQueue.hpp"
namespace aare {
class ClusterCollector{
ProducerConsumerQueue<ClusterVector<int>>* 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<int>> m_clusters;
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(){
void process() {
m_stopped = false;
fmt::print("ClusterCollector started\n");
while (!m_stop_requested || !m_source->isEmpty()) {
if (ClusterVector<int> *clusters = m_source->frontPtr();
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{
} else {
std::this_thread::sleep_for(m_default_wait);
}
}
@ -32,21 +34,25 @@ class ClusterCollector{
m_stopped = true;
}
public:
ClusterCollector(ClusterFinderMT<uint16_t, double, int32_t>* source){
m_source = source->sink();
m_thread = std::thread(&ClusterCollector::process, this);
}
void stop(){
m_stop_requested = true;
m_thread.join();
}
std::vector<ClusterVector<int>> steal_clusters(){
if(!m_stopped){
throw std::runtime_error("ClusterCollector is still running");
}
return std::move(m_clusters);
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

469
include/aare/ClusterFile.hpp
View File

@ -2,6 +2,7 @@
#include "aare/Cluster.hpp"
#include "aare/ClusterVector.hpp"
#include "aare/GainMap.hpp"
#include "aare/NDArray.hpp"
#include "aare/defs.hpp"
#include <filesystem>
@ -10,43 +11,18 @@
namespace aare {
/*
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! Legacy enums, migrate to enum class
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 Eta2 {
double x;
double y;
corner c;
int32_t sum;
};
struct ClusterAnalysis {
uint32_t c;
int32_t tot;
double etax;
double etay;
};
// 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:
@ -59,14 +35,19 @@ struct ClusterAnalysis {
* 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{}; /*Number of photons left in frame*/
size_t m_chunk_size{}; /*Number of clusters to read at a time*/
const 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<NDArray<double, 2>> m_gain_map; /*Gain map to apply to the clusters, will be applied if set*/
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:
/**
@ -79,74 +60,390 @@ class ClusterFile {
* @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");
~ClusterFile();
const std::string &mode = "r")
: m_filename(fname.string()), m_chunk_size(chunk_size), m_mode(mode) {
if (mode == "r") {
fp = fopen(m_filename.c_str(), "rb");
if (!fp) {
throw std::runtime_error("Could not open file for reading: " +
m_filename);
}
} else if (mode == "w") {
fp = fopen(m_filename.c_str(), "wb");
if (!fp) {
throw std::runtime_error("Could not open file for writing: " +
m_filename);
}
} else if (mode == "a") {
fp = fopen(m_filename.c_str(), "ab");
if (!fp) {
throw std::runtime_error("Could not open file for appending: " +
m_filename);
}
} else {
throw std::runtime_error("Unsupported mode: " + mode);
}
}
~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
* @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<int32_t> read_clusters(size_t n_clusters);
ClusterVector<int32_t> read_clusters(size_t n_clusters, ROI roi);
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
* @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<int32_t> read_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);
}
void write_frame(const ClusterVector<int32_t> &clusters);
/**
* @brief Return the chunk size
*/
size_t chunk_size() const { return m_chunk_size; }
/**
* @brief Set the region of interest to use when reading clusters. If set only clusters within
* the ROI will be read.
* @brief 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);
void set_roi(ROI roi) { m_roi = roi; }
/**
* @brief Set the noise map to use when reading clusters. If set clusters below the noise
* level will be discarded. Selection criteria one of: Central pixel above noise, highest
* 2x2 sum above 2 * noise, total sum above 3 * noise.
* @brief 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);
void set_noise_map(const NDView<int32_t, 2> noise_map) {
m_noise_map = NDArray<int32_t, 2>(noise_map);
}
/**
* @brief Set the gain map to use when reading clusters. If set the gain map will be applied
* to the clusters that pass ROI and noise_map selection. The gain map is expected to be in ADU/energy.
* @brief 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);
/**
* @brief Close the file. If not closed the file will be closed in the destructor
*/
void close();
void set_gain_map(const NDView<double, 2> gain_map) {
m_gain_map = InvertedGainMap(gain_map);
}
private:
ClusterVector<int32_t> read_clusters_with_cut(size_t n_clusters);
ClusterVector<int32_t> read_clusters_without_cut(size_t n_clusters);
ClusterVector<int32_t> read_frame_with_cut();
ClusterVector<int32_t> read_frame_without_cut();
bool is_selected(Cluster3x3 &cl);
Cluster3x3 read_one_cluster();
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();
};
//TODO! helper functions that doesn't really belong here
NDArray<double, 2> calculate_eta2(ClusterVector<int> &clusters);
Eta2 calculate_eta2(Cluster3x3 &cl);
Eta2 calculate_eta2(Cluster2x2 &cl);
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

58
include/aare/ClusterFileSink.hpp
View File

@ -3,35 +3,41 @@
#include <filesystem>
#include <thread>
#include "aare/ProducerConsumerQueue.hpp"
#include "aare/ClusterVector.hpp"
#include "aare/ClusterFinderMT.hpp"
#include "aare/ClusterVector.hpp"
#include "aare/ProducerConsumerQueue.hpp"
namespace aare{
namespace aare {
class ClusterFileSink{
ProducerConsumerQueue<ClusterVector<int>>* m_source;
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(){
void process() {
m_stopped = false;
fmt::print("ClusterFileSink started\n");
while (!m_stop_requested || !m_source->isEmpty()) {
if (ClusterVector<int> *clusters = m_source->frontPtr();
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?
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_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{
} else {
std::this_thread::sleep_for(m_default_wait);
}
}
@ -39,18 +45,18 @@ class ClusterFileSink{
m_stopped = true;
}
public:
ClusterFileSink(ClusterFinderMT<uint16_t, double, int32_t>* 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();
}
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

78
include/aare/ClusterFinder.hpp
View File

@ -10,17 +10,19 @@
namespace aare {
template <typename FRAME_TYPE = uint16_t, typename PEDESTAL_TYPE = double,
typename CT = int32_t>
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 PEDESTAL_TYPE m_nSigma;
const PEDESTAL_TYPE c2;
const PEDESTAL_TYPE c3;
Pedestal<PEDESTAL_TYPE> m_pedestal;
ClusterVector<CT> m_clusters;
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:
/**
@ -31,15 +33,12 @@ class ClusterFinder {
* @param capacity initial capacity of the cluster vector
*
*/
ClusterFinder(Shape<2> image_size, Shape<2> cluster_size,
PEDESTAL_TYPE nSigma = 5.0, size_t capacity = 1000000)
: m_image_size(image_size), m_cluster_sizeX(cluster_size[0]),
m_cluster_sizeY(cluster_size[1]),
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]),
m_clusters(m_cluster_sizeX, m_cluster_sizeY, capacity) {};
ClusterFinder(Shape<2> image_size, PEDESTAL_TYPE nSigma = 5.0,
size_t capacity = 1000000)
: m_image_size(image_size), m_nSigma(nSigma),
c2(sqrt((ClusterSizeY + 1) / 2 * (ClusterSizeX + 1) / 2)),
c3(sqrt(ClusterSizeX * ClusterSizeY)),
m_pedestal(image_size[0], image_size[1]), m_clusters(capacity) {};
void push_pedestal_frame(NDView<FRAME_TYPE, 2> frame) {
m_pedestal.push(frame);
@ -56,23 +55,28 @@ class ClusterFinder {
* same capacity as the old one
*
*/
ClusterVector<CT> steal_clusters(bool realloc_same_capacity = false) {
ClusterVector<CT> tmp = std::move(m_clusters);
ClusterVector<ClusterType>
steal_clusters(bool realloc_same_capacity = false) {
ClusterVector<ClusterType> tmp = std::move(m_clusters);
if (realloc_same_capacity)
m_clusters = ClusterVector<CT>(m_cluster_sizeX, m_cluster_sizeY,
tmp.capacity());
m_clusters = ClusterVector<ClusterType>(tmp.capacity());
else
m_clusters = ClusterVector<CT>(m_cluster_sizeX, m_cluster_sizeY);
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 = m_cluster_sizeY / 2;
int dx = m_cluster_sizeX / 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;
m_clusters.set_frame_number(frame_number);
std::vector<CT> cluster_data(m_cluster_sizeX * m_cluster_sizeY);
for (int iy = 0; iy < frame.shape(0); iy++) {
for (int ix = 0; ix < frame.shape(1); ix++) {
@ -87,8 +91,8 @@ class ClusterFinder {
continue; // NEGATIVE_PEDESTAL go to next pixel
// TODO! No pedestal update???
for (int ir = -dy; ir < dy + 1; ir++) {
for (int ic = -dx; ic < dx + 1; ic++) {
for (int ir = -dy; ir < dy + has_center_pixel_y; ir++) {
for (int ic = -dx; ic < dx + has_center_pixel_x; ic++) {
if (ix + ic >= 0 && ix + ic < frame.shape(1) &&
iy + ir >= 0 && iy + ir < frame.shape(0)) {
PEDESTAL_TYPE val =
@ -109,27 +113,33 @@ class ClusterFinder {
// pass
} else {
// m_pedestal.push(iy, ix, frame(iy, ix)); // Safe option
m_pedestal.push_fast(iy, ix, frame(iy, ix)); // Assume we have reached n_samples in the pedestal, slight performance improvement
continue; // It was a pedestal value nothing to store
m_pedestal.push_fast(
iy, ix,
frame(iy,
ix)); // Assume we have reached n_samples in the
// pedestal, slight performance improvement
continue; // It was a pedestal value nothing to store
}
// Store cluster
if (value == max) {
// Zero out the cluster data
std::fill(cluster_data.begin(), cluster_data.end(), 0);
ClusterType cluster{};
cluster.x = ix;
cluster.y = iy;
// Fill the cluster data since we have a photon to store
// It's worth redoing the look since most of the time we
// don't have a photon
int i = 0;
for (int ir = -dy; ir < dy + 1; ir++) {
for (int ic = -dx; ic < dx + 1; ic++) {
for (int ir = -dy; ir < dy + has_center_pixel_y; ir++) {
for (int ic = -dx; ic < dx + has_center_pixel_y; ic++) {
if (ix + ic >= 0 && ix + ic < frame.shape(1) &&
iy + ir >= 0 && iy + ir < frame.shape(0)) {
CT tmp =
static_cast<CT>(frame(iy + ir, ix + ic)) -
m_pedestal.mean(iy + ir, ix + ic);
cluster_data[i] =
static_cast<CT>(
m_pedestal.mean(iy + ir, ix + ic));
cluster.data[i] =
tmp; // Watch for out of bounds access
i++;
}
@ -137,9 +147,7 @@ class ClusterFinder {
}
// Add the cluster to the output ClusterVector
m_clusters.push_back(
ix, iy,
reinterpret_cast<std::byte *>(cluster_data.data()));
m_clusters.push_back(cluster);
}
}
}

35
include/aare/ClusterFinderMT.hpp
View File

@ -30,14 +30,17 @@ struct FrameWrapper {
* @tparam PEDESTAL_TYPE type of the pedestal data
* @tparam CT type of the cluster data
*/
template <typename FRAME_TYPE = uint16_t, typename PEDESTAL_TYPE = double,
typename CT = int32_t>
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<FRAME_TYPE, PEDESTAL_TYPE, CT>;
using Finder = ClusterFinder<ClusterType, FRAME_TYPE, PEDESTAL_TYPE>;
using InputQueue = ProducerConsumerQueue<FrameWrapper>;
using OutputQueue = ProducerConsumerQueue<ClusterVector<int>>;
using OutputQueue = ProducerConsumerQueue<ClusterVector<ClusterType>>;
std::vector<std::unique_ptr<InputQueue>> m_input_queues;
std::vector<std::unique_ptr<OutputQueue>> m_output_queues;
@ -48,6 +51,7 @@ class ClusterFinderMT {
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};
@ -66,7 +70,8 @@ class ClusterFinderMT {
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));
m_output_queues[thread_id]->write(
cf->steal_clusters(realloc_same_capacity));
break;
case FrameType::PEDESTAL:
@ -114,28 +119,32 @@ class ClusterFinderMT {
* expected number of clusters in a frame per frame.
* @param n_threads number of threads to use
*/
ClusterFinderMT(Shape<2> image_size, Shape<2> cluster_size,
PEDESTAL_TYPE nSigma = 5.0, size_t capacity = 2000,
size_t n_threads = 3)
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<FRAME_TYPE, PEDESTAL_TYPE, CT>>(
image_size, cluster_size, nSigma, capacity));
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?
// 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
* @warning You need to empty this queue otherwise the cluster finder will
* wait forever
*/
ProducerConsumerQueue<ClusterVector<int>> *sink() { return &m_sink; }
ProducerConsumerQueue<ClusterVector<ClusterType>> *sink() {
return &m_sink;
}
/**
* @brief Start all processing threads

280
include/aare/ClusterVector.hpp
View File

@ -1,4 +1,5 @@
#pragma once
#include "aare/Cluster.hpp" //TODO maybe store in seperate file !!!
#include <algorithm>
#include <array>
#include <cstddef>
@ -13,292 +14,157 @@
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.
* @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.
* @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, typename CoordType = int16_t> class ClusterVector {
using value_type = T;
size_t m_cluster_size_x;
size_t m_cluster_size_y;
std::byte *m_data{};
size_t m_size{0};
size_t m_capacity;
uint64_t m_frame_number{0}; // TODO! Check frame number size and type
/*
Format string used in the python bindings to create a numpy
array from the buffer
= - native byte order
h - short
d - double
i - int
*/
constexpr static char m_fmt_base[] = "=h:x:\nh:y:\n({},{}){}:data:";
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 cluster_size_x size of the cluster in x direction
* @param cluster_size_y size of the cluster in y direction
* @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 cluster_size_x = 3, size_t cluster_size_y = 3,
size_t capacity = 1024, uint64_t frame_number = 0)
: m_cluster_size_x(cluster_size_x), m_cluster_size_y(cluster_size_y),
m_capacity(capacity), m_frame_number(frame_number) {
allocate_buffer(capacity);
ClusterVector(size_t capacity = 1024, uint64_t frame_number = 0)
: m_frame_number(frame_number) {
m_data.reserve(capacity);
}
~ClusterVector() { delete[] m_data; }
// Move constructor
ClusterVector(ClusterVector &&other) noexcept
: m_cluster_size_x(other.m_cluster_size_x),
m_cluster_size_y(other.m_cluster_size_y), m_data(other.m_data),
m_size(other.m_size), m_capacity(other.m_capacity),
m_frame_number(other.m_frame_number) {
other.m_data = nullptr;
other.m_size = 0;
other.m_capacity = 0;
: 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) {
delete[] m_data;
m_cluster_size_x = other.m_cluster_size_x;
m_cluster_size_y = other.m_cluster_size_y;
m_data = other.m_data;
m_size = other.m_size;
m_capacity = other.m_capacity;
m_frame_number = other.m_frame_number;
other.m_data = nullptr;
other.m_size = 0;
other.m_capacity = 0;
other.m_data.clear();
other.m_frame_number = 0;
}
return *this;
}
/**
* @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) {
if (capacity > m_capacity) {
allocate_buffer(capacity);
}
}
/**
* @brief Add a cluster to the vector
* @param x x-coordinate of the cluster
* @param y y-coordinate of the cluster
* @param data pointer to the data of the cluster
* @warning The data pointer must point to a buffer of size cluster_size_x *
* cluster_size_y * sizeof(T)
*/
void push_back(CoordType x, CoordType y, const std::byte *data) {
if (m_size == m_capacity) {
allocate_buffer(m_capacity * 2);
}
std::byte *ptr = element_ptr(m_size);
*reinterpret_cast<CoordType *>(ptr) = x;
ptr += sizeof(CoordType);
*reinterpret_cast<CoordType *>(ptr) = y;
ptr += sizeof(CoordType);
std::copy(data, data + m_cluster_size_x * m_cluster_size_y * sizeof(T),
ptr);
m_size++;
}
ClusterVector &operator+=(const ClusterVector &other) {
if (m_size + other.m_size > m_capacity) {
allocate_buffer(m_capacity + other.m_size);
}
std::copy(other.m_data, other.m_data + other.m_size * item_size(),
m_data + m_size * item_size());
m_size += other.m_size;
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_size);
const size_t stride = item_size();
const size_t n_pixels = m_cluster_size_x * m_cluster_size_y;
std::byte *ptr = m_data + 2 * sizeof(CoordType); // skip x and y
std::vector<T> sums(m_data.size());
std::transform(
m_data.begin(), m_data.end(), sums.begin(),
[](const ClusterType &cluster) { return cluster.sum(); });
for (size_t i = 0; i < m_size; i++) {
sums[i] =
std::accumulate(reinterpret_cast<T *>(ptr),
reinterpret_cast<T *>(ptr) + n_pixels, T{});
ptr += stride;
}
return sums;
}
/**
* @brief Return the maximum sum of the 2x2 subclusters in each cluster
* @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
* @throws std::runtime_error if the cluster size is not 3x3
* @warning Only 3x3 clusters are supported for the 2x2 sum.
*/
std::vector<T> sum_2x2() {
std::vector<T> sums(m_size);
const size_t stride = item_size();
std::vector<T> sums_2x2(m_data.size());
if (m_cluster_size_x != 3 || m_cluster_size_y != 3) {
throw std::runtime_error(
"Only 3x3 clusters are supported for the 2x2 sum.");
}
std::byte *ptr = m_data + 2 * sizeof(CoordType); // skip x and y
std::transform(m_data.begin(), m_data.end(), sums_2x2.begin(),
[](const ClusterType &cluster) {
return cluster.max_sum_2x2().first;
});
for (size_t i = 0; i < m_size; i++) {
std::array<T, 4> total;
auto T_ptr = reinterpret_cast<T *>(ptr);
total[0] = T_ptr[0] + T_ptr[1] + T_ptr[3] + T_ptr[4];
total[1] = T_ptr[1] + T_ptr[2] + T_ptr[4] + T_ptr[5];
total[2] = T_ptr[3] + T_ptr[4] + T_ptr[6] + T_ptr[7];
total[3] = T_ptr[4] + T_ptr[5] + T_ptr[7] + T_ptr[8];
return sums_2x2;
}
sums[i] = *std::max_element(total.begin(), total.end());
ptr += stride;
}
/**
* @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); }
return sums;
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_size; }
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_capacity; }
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 2 * sizeof(CoordType) +
m_cluster_size_x * m_cluster_size_y * sizeof(T);
return sizeof(ClusterType); // 2 * sizeof(CoordType) + ClusterSizeX *
// ClusterSizeY * sizeof(T);
}
/**
* @brief Return the offset in bytes for the i-th cluster
*/
size_t element_offset(size_t i) const { return item_size() * i; }
/**
* @brief Return a pointer to the i-th cluster
*/
std::byte *element_ptr(size_t i) { return m_data + element_offset(i); }
/**
* @brief Return a pointer to the i-th cluster
*/
const std::byte *element_ptr(size_t i) const {
return m_data + element_offset(i);
}
size_t cluster_size_x() const { return m_cluster_size_x; }
size_t cluster_size_y() const { return m_cluster_size_y; }
std::byte *data() { return m_data; }
std::byte const *data() const { return m_data; }
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
*/
template <typename V> V &at(size_t i) {
return *reinterpret_cast<V *>(element_ptr(i));
}
ClusterType &operator[](size_t i) { return m_data[i]; }
template <typename V> const V &at(size_t i) const {
return *reinterpret_cast<const V *>(element_ptr(i));
}
const std::string_view fmt_base() const {
// TODO! how do we match on coord_t?
return m_fmt_base;
}
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
*/
uint64_t frame_number() const { return m_frame_number; }
int32_t frame_number() const { return m_frame_number; }
void set_frame_number(uint64_t frame_number) {
void set_frame_number(int32_t frame_number) {
m_frame_number = frame_number;
}
/**
* @brief Resize the vector to contain new_size clusters. If new_size is
* greater than the current capacity, a new buffer is allocated. If the size
* is smaller no memory is freed, size is just updated.
* @param new_size new size of the vector
* @warning The additional clusters are not initialized
*/
void resize(size_t new_size) {
// TODO! Should we initialize the new clusters?
if (new_size > m_capacity) {
allocate_buffer(new_size);
}
m_size = new_size;
}
void apply_gain_map(const NDView<double> gain_map){
//in principle we need to know the size of the image for this lookup
//TODO! check orientations
std::array<int64_t, 9> xcorr = {-1, 0, 1, -1, 0, 1, -1, 0, 1};
std::array<int64_t, 9> ycorr = {-1, -1, -1, 0, 0, 0, 1, 1, 1};
for (size_t i=0; i<m_size; i++){
auto& cl = at<Cluster3x3>(i);
if (cl.x > 0 && cl.y > 0 && cl.x < gain_map.shape(1)-1 && cl.y < gain_map.shape(0)-1){
for (size_t j=0; j<9; j++){
size_t x = cl.x + xcorr[j];
size_t y = cl.y + ycorr[j];
cl.data[j] = static_cast<T>(cl.data[j] * gain_map(y, x));
}
}else{
memset(cl.data, 0, 9*sizeof(T)); //clear edge clusters
}
}
}
private:
void allocate_buffer(size_t new_capacity) {
size_t num_bytes = item_size() * new_capacity;
std::byte *new_data = new std::byte[num_bytes]{};
std::copy(m_data, m_data + item_size() * m_size, new_data);
delete[] m_data;
m_data = new_data;
m_capacity = new_capacity;
}
};
} // namespace aare

4
include/aare/FilePtr.hpp
View File

@ -18,8 +18,8 @@ class FilePtr {
FilePtr(FilePtr &&other);
FilePtr &operator=(FilePtr &&other);
FILE *get();
int64_t tell();
void seek(int64_t offset, int whence = SEEK_SET) {
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");
}

29
include/aare/Fit.hpp
View File

@ -15,6 +15,12 @@ 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
@ -25,6 +31,9 @@ std::array<double, 3> gaus_init_par(const NDView<double, 1> x, const NDView<doub
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;
/**
@ -38,7 +47,7 @@ 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 vales, layout [row, col, values]
* @param y y values, layout [row, col, values]
* @param n_threads number of threads to use
*/
@ -51,7 +60,7 @@ NDArray<double, 3> fit_gaus(NDView<double, 1> x, NDView<double, 3> y,
/**
* @brief Fit a 1D Gaussian with error estimates
* @param x x values
* @param y y vales, layout [row, col, 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
@ -64,7 +73,7 @@ void fit_gaus(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
* @brief Fit a 1D Gaussian to each pixel with error estimates. Data layout
* [row, col, values]
* @param x x values
* @param y y vales, layout [row, col, 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]
@ -88,5 +97,19 @@ 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

121
include/aare/Interpolator.hpp
View File

@ -1,29 +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/ClusterVector.hpp"
#include "aare/ClusterFile.hpp" //Cluster_3x3
namespace aare{
#include "aare/algorithm.hpp"
struct Photon{
namespace aare {
struct Photon {
double x;
double y;
double energy;
};
class Interpolator{
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;}
std::vector<Photon> interpolate(const ClusterVector<int32_t>& clusters);
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

78
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);
}
@ -186,22 +186,22 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
}
// TODO! is int the right type for index?
T &operator()(int64_t i) { return data_[i]; }
const T &operator()(int64_t 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[](int64_t i) { return data_[i]; }
const T &operator[](int64_t 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_); }
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);
@ -228,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) {
@ -242,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_) {
@ -254,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_) {
@ -266,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_) {
@ -278,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_};
@ -297,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_;
@ -310,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;
@ -322,23 +322,23 @@ 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;
@ -348,57 +348,57 @@ NDArray<T, Ndim> &NDArray<T, Ndim>::operator+=(const T &value) {
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() {
// 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) {
@ -410,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);
@ -419,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);
@ -429,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);
@ -437,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);

44
include/aare/NDView.hpp
View File

@ -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<>())) {}
@ -73,14 +73,14 @@ template <typename T, int64_t Ndim = 2> class NDView : public ArrayExpr<NDView<T
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_)
@ -136,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) {
@ -160,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);
@ -171,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) {
@ -186,7 +186,7 @@ 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<int64_t>(vec.size())});
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);

4
include/aare/Pedestal.hpp
View File

@ -107,7 +107,7 @@ template <typename SUM_TYPE = double> class Pedestal {
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");
}
@ -128,7 +128,7 @@ template <typename SUM_TYPE = double> class Pedestal {
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");
}

22
include/aare/RawFile.hpp
View File

@ -30,22 +30,11 @@ 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<std::unique_ptr<RawSubFile>> m_subfiles;
ModuleConfig cfg{0, 0};
RawMasterFile m_master;
size_t m_current_frame{};
// std::vector<ModuleGeometry> m_module_pixel_0;
// size_t m_rows{};
// size_t m_cols{};
size_t m_current_subfile{};
DetectorGeometry m_geometry;
public:
@ -56,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;
@ -80,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;
@ -115,9 +104,6 @@ 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();
};

1
include/aare/RawMasterFile.hpp
View File

@ -121,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;

20
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 m_num_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{};
@ -67,12 +76,17 @@ class RawSubFile {
size_t pixels_per_frame() const { return m_rows * m_cols; }
size_t bytes_per_pixel() const { return m_bitdepth / bits_per_byte; }
size_t frames_in_file() const { return m_num_frames; }
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;
};
} // namespace aare

2
include/aare/VarClusterFinder.hpp
View File

@ -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_;

11
include/aare/algorithm.hpp
View File

@ -107,5 +107,16 @@ std::vector<T> cumsum(const std::vector<T>& vec) {
}
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

18
include/aare/defs.hpp
View File

@ -204,23 +204,25 @@ struct DetectorGeometry{
int module_gap_row{};
int module_gap_col{};
std::vector<ModuleGeometry> module_pixel_0;
auto size() const { return module_pixel_0.size(); }
};
struct ROI{
int64_t xmin{};
int64_t xmax{};
int64_t ymin{};
int64_t ymax{};
ssize_t xmin{};
ssize_t xmax{};
ssize_t ymin{};
ssize_t ymax{};
int64_t height() const { return ymax - ymin; }
int64_t width() const { return xmax - xmin; }
bool contains(int64_t x, int64_t y) const {
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<int64_t>;
using dynamic_shape = std::vector<ssize_t>;
//TODO! Can we uniform enums between the libraries?

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

8
pyproject.toml
View File

@ -1,10 +1,16 @@
[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 = "2025.4.22"
dynamic = ["version"]
requires-python = ">=3.11"
dependencies = [
"numpy",

4
python/CMakeLists.txt
View File

@ -29,6 +29,9 @@ 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
@ -36,6 +39,7 @@ set( PYTHON_FILES
aare/utils.py
)
# Copy the python files to the build directory
foreach(FILE ${PYTHON_FILES})
configure_file(${FILE} ${CMAKE_BINARY_DIR}/${FILE} )

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.")

13
python/aare/__init__.py
View File

@ -3,16 +3,21 @@ from . import _aare
from ._aare import File, RawMasterFile, RawSubFile, JungfrauDataFile
from ._aare import Pedestal_d, Pedestal_f, ClusterFinder, VarClusterFinder
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 ._aare import ClusterFinderMT, ClusterCollector, ClusterFileSink, ClusterVector_i
from ._aare import fit_gaus, fit_pol1
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

2
python/aare/func.py
View File

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

126
python/examples/play.py
View File

@ -1,79 +1,89 @@
import sys
sys.path.append('/home/l_msdetect/erik/aare/build')
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
from aare import RawSubFile, DetectorType, RawFile
from pathlib import Path
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)
# f = RawFile(path/"jungfrau_single_master_0.json")
# 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)))
# 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]
# 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"
# 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)
# 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])
# 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
# #Generate the hit
tmp = p.interpolate(v)
print(f'tmp:{tmp}')
pos = np.array((tmp['x'], tmp['y']))*25
# 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=}')
# 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=}')

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

@ -0,0 +1,104 @@
#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;
});
}

292
python/src/cluster.hpp
View File

@ -16,194 +16,196 @@
namespace py = pybind11;
using pd_type = double;
template <typename T>
void define_cluster_vector(py::module &m, const std::string &typestr) {
auto class_name = fmt::format("ClusterVector_{}", typestr);
py::class_<ClusterVector<T>>(m, class_name.c_str(), py::buffer_protocol())
.def(py::init<int, int>(),
py::arg("cluster_size_x") = 3, py::arg("cluster_size_y") = 3)
.def("push_back",
[](ClusterVector<T> &self, int x, int y, py::array_t<T> data) {
// auto view = make_view_2d(data);
self.push_back(x, y, reinterpret_cast<const std::byte*>(data.data()));
})
.def_property_readonly("size", &ClusterVector<T>::size)
.def("item_size", &ClusterVector<T>::item_size)
.def_property_readonly("fmt",
[typestr](ClusterVector<T> &self) {
return fmt::format(
self.fmt_base(), self.cluster_size_x(),
self.cluster_size_y(), typestr);
})
.def("sum",
[](ClusterVector<T> &self) {
auto *vec = new std::vector<T>(self.sum());
return return_vector(vec);
})
.def("sum_2x2", [](ClusterVector<T> &self) {
auto *vec = new std::vector<T>(self.sum_2x2());
return return_vector(vec);
})
.def_property_readonly("cluster_size_x", &ClusterVector<T>::cluster_size_x)
.def_property_readonly("cluster_size_y", &ClusterVector<T>::cluster_size_y)
.def_property_readonly("capacity", &ClusterVector<T>::capacity)
.def_property("frame_number", &ClusterVector<T>::frame_number,
&ClusterVector<T>::set_frame_number)
.def_buffer([typestr](ClusterVector<T> &self) -> py::buffer_info {
return py::buffer_info(
self.data(), /* Pointer to buffer */
self.item_size(), /* Size of one scalar */
fmt::format(self.fmt_base(), self.cluster_size_x(),
self.cluster_size_y(),
typestr), /* Format descriptor */
1, /* Number of dimensions */
{self.size()}, /* Buffer dimensions */
{self.item_size()} /* Strides (in bytes) for each index */
);
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;
}));
/*
.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!!!
});
*/
}
void define_cluster_finder_mt_bindings(py::module &m) {
py::class_<ClusterFinderMT<uint16_t, pd_type>>(m, "ClusterFinderMT")
.def(py::init<Shape<2>, Shape<2>, pd_type, size_t, size_t>(),
py::arg("image_size"), py::arg("cluster_size"),
py::arg("n_sigma") = 5.0, py::arg("capacity") = 2048,
py::arg("n_threads") = 3)
template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
typename CoordType = uint16_t>
void define_cluster_finder_mt_bindings(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<uint16_t, pd_type> &self,
[](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<uint16_t, pd_type> &self,
[](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("clear_pedestal", &ClusterFinderMT<uint16_t, pd_type>::clear_pedestal)
.def("sync", &ClusterFinderMT<uint16_t, pd_type>::sync)
.def("stop", &ClusterFinderMT<uint16_t, pd_type>::stop)
.def("start", &ClusterFinderMT<uint16_t, pd_type>::start)
.def("pedestal",
[](ClusterFinderMT<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<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);
.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);
}
void define_cluster_collector_bindings(py::module &m) {
py::class_<ClusterCollector>(m, "ClusterCollector")
.def(py::init<ClusterFinderMT<uint16_t, double, int32_t> *>())
.def("stop", &ClusterCollector::stop)
template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
typename CoordType = uint16_t>
void define_cluster_collector_bindings(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 &self) {
auto v =
new std::vector<ClusterVector<int>>(self.steal_clusters());
return v;
[](ClusterCollector<ClusterType> &self) {
auto v = new std::vector<ClusterVector<ClusterType>>(
self.steal_clusters());
return v; // TODO change!!!
},
py::return_value_policy::take_ownership);
}
void define_cluster_file_sink_bindings(py::module &m) {
py::class_<ClusterFileSink>(m, "ClusterFileSink")
.def(py::init<ClusterFinderMT<uint16_t, double, int32_t> *,
template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
typename CoordType = uint16_t>
void define_cluster_file_sink_bindings(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::stop);
.def("stop", &ClusterFileSink<ClusterType>::stop);
}
void define_cluster_finder_bindings(py::module &m) {
py::class_<ClusterFinder<uint16_t, pd_type>>(m, "ClusterFinder")
.def(py::init<Shape<2>, Shape<2>, pd_type, size_t>(),
py::arg("image_size"), py::arg("cluster_size"),
template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
typename CoordType = uint16_t>
void define_cluster_finder_bindings(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<uint16_t, pd_type> &self,
[](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<uint16_t, pd_type>::clear_pedestal)
.def_property_readonly("pedestal",
[](ClusterFinder<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<uint16_t, pd_type> &self) {
auto arr = new NDArray<pd_type, 2>{};
*arr = self.noise();
return return_image_data(arr);
})
.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<uint16_t, pd_type> &self,
[](ClusterFinder<ClusterType, uint16_t, pd_type> &self,
bool realloc_same_capacity) {
auto v = new ClusterVector<int>(
self.steal_clusters(realloc_same_capacity));
return v;
ClusterVector<ClusterType> clusters =
self.steal_clusters(realloc_same_capacity);
return clusters;
},
py::arg("realloc_same_capacity") = false)
.def(
"find_clusters",
[](ClusterFinder<uint16_t, pd_type> &self,
[](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);
m.def("hitmap",
[](std::array<size_t, 2> image_size, ClusterVector<int32_t> &cv) {
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;
size_t stride = cv.item_size();
auto ptr = cv.data();
for (size_t i = 0; i < cv.size(); i++) {
auto x = *reinterpret_cast<int16_t *>(ptr);
auto y = *reinterpret_cast<int16_t *>(ptr + sizeof(int16_t));
r(y, x) += 1;
ptr += stride;
}
return hitmap;
});
define_cluster_vector<int>(m, "i");
define_cluster_vector<double>(m, "d");
define_cluster_vector<float>(m, "f");
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) + ">";
});
}
}
#pragma GCC diagnostic pop

87
python/src/cluster_file.hpp
View File

@ -1,3 +1,4 @@
#include "aare/CalculateEta.hpp"
#include "aare/ClusterFile.hpp"
#include "aare/defs.hpp"
@ -10,64 +11,84 @@
#include <pybind11/stl/filesystem.h>
#include <string>
//Disable warnings for unused parameters, as we ignore some
//in the __exit__ method
// 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;
void define_cluster_file_io_bindings(py::module &m) {
PYBIND11_NUMPY_DTYPE(Cluster3x3, x, y, data);
template <typename Type, uint8_t CoordSizeX, uint8_t CoordSizeY,
typename CoordType = uint16_t>
void define_cluster_file_io_bindings(py::module &m,
const std::string &typestr) {
py::class_<ClusterFile>(m, "ClusterFile")
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 &self, size_t n_clusters) {
auto v = new ClusterVector<int32_t>(self.read_clusters(n_clusters));
.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)
},
py::return_value_policy::take_ownership)
.def("read_frame",
[](ClusterFile &self) {
auto v = new ClusterVector<int32_t>(self.read_frame());
return v;
[](ClusterFile<ClusterType> &self) {
auto v = new ClusterVector<ClusterType>(self.read_frame());
return v;
})
.def("set_roi", &ClusterFile::set_roi)
.def("set_noise_map", [](ClusterFile &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 &self, py::array_t<double> gain_map) {
auto view = make_view_2d(gain_map);
self.set_gain_map(view);
})
.def("close", &ClusterFile::close)
.def("write_frame", &ClusterFile::write_frame)
.def("__enter__", [](ClusterFile &self) { return &self; })
.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 &self,
[](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 &self) { return &self; })
.def("__next__", [](ClusterFile &self) {
auto v = new ClusterVector<int32_t>(self.read_clusters(self.chunk_size()));
.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;
});
}
m.def("calculate_eta2", []( aare::ClusterVector<int32_t> &clusters) {
auto eta2 = new NDArray<double, 2>(calculate_eta2(clusters));
return return_image_data(eta2);
});
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

4
python/src/ctb_raw_file.hpp
View File

@ -34,7 +34,7 @@ m.def("adc_sar_05_decode64to16", [](py::array_t<uint8_t> input) {
}
//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<int64_t>(bits_per_byte)};
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
@ -55,7 +55,7 @@ m.def("adc_sar_04_decode64to16", [](py::array_t<uint8_t> input) {
}
//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<int64_t>(bits_per_byte)};
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

2
python/src/file.hpp
View File

@ -198,7 +198,7 @@ void define_file_io_bindings(py::module &m) {
py::class_<ROI>(m, "ROI")
.def(py::init<>())
.def(py::init<int64_t, int64_t, int64_t, int64_t>(), py::arg("xmin"),
.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)

215
python/src/fit.hpp
View File

@ -55,6 +55,47 @@ void define_fit_bindings(py::module &m) {
)",
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",
@ -235,6 +276,180 @@ n_threads : int, optional
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

68
python/src/interpolation.hpp
View File

@ -8,31 +8,55 @@
#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);
PYBIND11_NUMPY_DTYPE(aare::Photon, x, y, energy);
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);
})
.def("interpolate", [](Interpolator& self, const ClusterVector<int32_t>& clusters){
auto photons = self.interpolate(clusters);
auto* ptr = new std::vector<Photon>{photons};
return return_vector(ptr);
});
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(

86
python/src/module.cpp
View File

@ -1,20 +1,24 @@
//Files with bindings to the different classes
#include "file.hpp"
#include "raw_file.hpp"
#include "ctb_raw_file.hpp"
#include "raw_master_file.hpp"
#include "var_cluster.hpp"
#include "pixel_map.hpp"
#include "pedestal.hpp"
// Files with bindings to the different classes
//New style file naming
#include "bind_ClusterVector.hpp"
//TODO! migrate the other names
#include "cluster.hpp"
#include "cluster_file.hpp"
#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"
#include "pixel_map.hpp"
#include "var_cluster.hpp"
#include "pedestal.hpp"
#include "jungfrau_data_file.hpp"
//Pybind stuff
// Pybind stuff
#include <pybind11/pybind11.h>
#include <pybind11/stl.h>
@ -30,13 +34,63 @@ PYBIND11_MODULE(_aare, m) {
define_pixel_map_bindings(m);
define_pedestal_bindings<double>(m, "Pedestal_d");
define_pedestal_bindings<float>(m, "Pedestal_f");
define_cluster_finder_bindings(m);
define_cluster_finder_mt_bindings(m);
define_cluster_file_io_bindings(m);
define_cluster_collector_bindings(m);
define_cluster_file_sink_bindings(m);
define_fit_bindings(m);
define_interpolation_bindings(m);
define_jungfrau_data_file_io_bindings(m);
}
define_cluster_file_io_bindings<int, 3, 3, uint16_t>(m, "Cluster3x3i");
define_cluster_file_io_bindings<double, 3, 3, uint16_t>(m, "Cluster3x3d");
define_cluster_file_io_bindings<float, 3, 3, uint16_t>(m, "Cluster3x3f");
define_cluster_file_io_bindings<int, 2, 2, uint16_t>(m, "Cluster2x2i");
define_cluster_file_io_bindings<float, 2, 2, uint16_t>(m, "Cluster2x2f");
define_cluster_file_io_bindings<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_cluster_finder_bindings<int, 3, 3, uint16_t>(m, "Cluster3x3i");
define_cluster_finder_bindings<double, 3, 3, uint16_t>(m, "Cluster3x3d");
define_cluster_finder_bindings<float, 3, 3, uint16_t>(m, "Cluster3x3f");
define_cluster_finder_bindings<int, 2, 2, uint16_t>(m, "Cluster2x2i");
define_cluster_finder_bindings<double, 2, 2, uint16_t>(m, "Cluster2x2d");
define_cluster_finder_bindings<float, 2, 2, uint16_t>(m, "Cluster2x2f");
define_cluster_finder_mt_bindings<int, 3, 3, uint16_t>(m, "Cluster3x3i");
define_cluster_finder_mt_bindings<double, 3, 3, uint16_t>(m, "Cluster3x3d");
define_cluster_finder_mt_bindings<float, 3, 3, uint16_t>(m, "Cluster3x3f");
define_cluster_finder_mt_bindings<int, 2, 2, uint16_t>(m, "Cluster2x2i");
define_cluster_finder_mt_bindings<double, 2, 2, uint16_t>(m, "Cluster2x2d");
define_cluster_finder_mt_bindings<float, 2, 2, uint16_t>(m, "Cluster2x2f");
define_cluster_file_sink_bindings<int, 3, 3, uint16_t>(m, "Cluster3x3i");
define_cluster_file_sink_bindings<double, 3, 3, uint16_t>(m, "Cluster3x3d");
define_cluster_file_sink_bindings<float, 3, 3, uint16_t>(m, "Cluster3x3f");
define_cluster_file_sink_bindings<int, 2, 2, uint16_t>(m, "Cluster2x2i");
define_cluster_file_sink_bindings<double, 2, 2, uint16_t>(m, "Cluster2x2d");
define_cluster_file_sink_bindings<float, 2, 2, uint16_t>(m, "Cluster2x2f");
define_cluster_collector_bindings<int, 3, 3, uint16_t>(m, "Cluster3x3i");
define_cluster_collector_bindings<double, 3, 3, uint16_t>(m, "Cluster3x3f");
define_cluster_collector_bindings<float, 3, 3, uint16_t>(m, "Cluster3x3d");
define_cluster_collector_bindings<int, 2, 2, uint16_t>(m, "Cluster2x2i");
define_cluster_collector_bindings<double, 2, 2, uint16_t>(m, "Cluster2x2f");
define_cluster_collector_bindings<float, 2, 2, uint16_t>(m, "Cluster2x2d");
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);
}

38
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) {
@ -40,25 +41,46 @@ template <typename T> py::array return_vector(std::vector<T> *vec) {
}
// todo rewrite generic
template <class T, int Flags> auto get_shape_3d(const 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(const 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 get_shape_1d(const 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) {
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) {
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();

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);
}

7
python/tests/conftest.py
View File

@ -25,5 +25,10 @@ def pytest_collection_modifyitems(config, items):
@pytest.fixture
def test_data_path():
return Path(os.environ["AARE_TEST_DATA"])
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()

19
python/tests/test_RawSubFile.py
View File

@ -5,32 +5,35 @@ 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 == 3
assert f.frames_in_file == 7
headers, frames = f.read()
assert headers.size == 3
assert frames.shape == (3, 512, 1024)
assert headers.size == 7
assert frames.shape == (7, 512, 1024)
# Frame numbers in this file should be 4, 5, 6
for i,h in zip(range(4,7,1), headers):
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:6] == frames)
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 == 3
assert header["frameNumber"] == 3
assert i == 10
assert header["frameNumber"] == 10

127
src/CalculateEta.test.cpp Normal file
View File

@ -0,0 +1,127 @@
/************************************************
* @file CalculateEta.test.cpp
* @short test case to calculate_eta2
***********************************************/
#include "aare/CalculateEta.hpp"
#include "aare/Cluster.hpp"
#include "aare/ClusterFile.hpp"
// #include "catch.hpp"
#include <array>
#include <catch2/catch_all.hpp>
#include <catch2/catch_test_macros.hpp>
using namespace aare;
using ClusterTypes =
std::variant<Cluster<int, 2, 2>, Cluster<int, 3, 3>, Cluster<int, 5, 5>,
Cluster<int, 4, 2>, Cluster<int, 2, 3>>;
auto get_test_parameters() {
return GENERATE(
std::make_tuple(ClusterTypes{Cluster<int, 2, 2>{0, 0, {1, 2, 3, 1}}},
Eta2<int>{2. / 3, 3. / 4,
static_cast<int>(corner::cBottomLeft), 7}),
std::make_tuple(
ClusterTypes{Cluster<int, 3, 3>{0, 0, {1, 2, 3, 4, 5, 6, 1, 2, 7}}},
Eta2<int>{6. / 11, 2. / 7, static_cast<int>(corner::cTopRight),
20}),
std::make_tuple(ClusterTypes{Cluster<int, 5, 5>{
0, 0, {1, 6, 7, 6, 5, 4, 3, 2, 1, 2, 8, 9, 8,
1, 4, 5, 6, 7, 8, 4, 1, 1, 1, 1, 1}}},
Eta2<int>{8. / 17, 7. / 15, 9, 30}),
std::make_tuple(
ClusterTypes{Cluster<int, 4, 2>{0, 0, {1, 4, 7, 2, 5, 6, 4, 3}}},
Eta2<int>{4. / 10, 4. / 11, 1, 21}),
std::make_tuple(
ClusterTypes{Cluster<int, 2, 3>{0, 0, {1, 3, 2, 3, 4, 2}}},
Eta2<int>{3. / 5, 2. / 5, 1, 11}));
}
TEST_CASE("compute_largest_2x2_subcluster", "[eta_calculation]") {
auto [cluster, expected_eta] = get_test_parameters();
auto [sum, index] = std::visit(
[](const auto &clustertype) { return clustertype.max_sum_2x2(); },
cluster);
CHECK(expected_eta.c == index);
CHECK(expected_eta.sum == sum);
}
TEST_CASE("calculate_eta2", "[eta_calculation]") {
auto [cluster, expected_eta] = get_test_parameters();
auto eta = std::visit(
[](const auto &clustertype) { return calculate_eta2(clustertype); },
cluster);
CHECK(eta.x == expected_eta.x);
CHECK(eta.y == expected_eta.y);
CHECK(eta.c == expected_eta.c);
CHECK(eta.sum == expected_eta.sum);
}
// 3x3 cluster layout (rotated to match the cBottomLeft enum):
// 6, 7, 8
// 3, 4, 5
// 0, 1, 2
TEST_CASE("Calculate eta2 for a 3x3 int32 cluster with the largest 2x2 sum in "
"the bottom left",
"[eta_calculation]") {
// Create a 3x3 cluster
Cluster<int32_t, 3, 3> cl;
cl.x = 0;
cl.y = 0;
cl.data[0] = 30;
cl.data[1] = 23;
cl.data[2] = 5;
cl.data[3] = 20;
cl.data[4] = 50;
cl.data[5] = 3;
cl.data[6] = 8;
cl.data[7] = 2;
cl.data[8] = 3;
// 8, 2, 3
// 20, 50, 3
// 30, 23, 5
auto eta = calculate_eta2(cl);
CHECK(eta.c == static_cast<int>(corner::cBottomLeft));
CHECK(eta.x == 50.0 / (20 + 50)); // 4/(3+4)
CHECK(eta.y == 50.0 / (23 + 50)); // 4/(1+4)
CHECK(eta.sum == 30 + 23 + 20 + 50);
}
TEST_CASE("Calculate eta2 for a 3x3 int32 cluster with the largest 2x2 sum in "
"the top left",
"[eta_calculation]") {
// Create a 3x3 cluster
Cluster<int32_t, 3, 3> cl;
cl.x = 0;
cl.y = 0;
cl.data[0] = 8;
cl.data[1] = 12;
cl.data[2] = 5;
cl.data[3] = 77;
cl.data[4] = 80;
cl.data[5] = 3;
cl.data[6] = 82;
cl.data[7] = 91;
cl.data[8] = 3;
// 82, 91, 3
// 77, 80, 3
// 8, 12, 5
auto eta = calculate_eta2(cl);
CHECK(eta.c == static_cast<int>(corner::cTopLeft));
CHECK(eta.x == 80. / (77 + 80)); // 4/(3+4)
CHECK(eta.y == 91.0 / (91 + 80)); // 7/(7+4)
CHECK(eta.sum == 77 + 80 + 82 + 91);
}

21
src/Cluster.test.cpp Normal file
View File

@ -0,0 +1,21 @@
/************************************************
* @file test-Cluster.cpp
* @short test case for generic Cluster, ClusterVector, and calculate_eta2
***********************************************/
#include "aare/Cluster.hpp"
#include "aare/CalculateEta.hpp"
#include "aare/ClusterFile.hpp"
// #include "catch.hpp"
#include <array>
#include <catch2/catch_all.hpp>
#include <catch2/catch_test_macros.hpp>
using namespace aare;
TEST_CASE("Test sum of Cluster", "[.cluster]") {
Cluster<int, 2, 2> cluster{0, 0, {1, 2, 3, 4}};
CHECK(cluster.sum() == 10);
}

301
src/ClusterFile.test.cpp
View File

@ -1,35 +1,39 @@
#include "aare/ClusterFile.hpp"
#include "test_config.hpp"
#include "aare/defs.hpp"
#include <algorithm>
#include <catch2/catch_test_macros.hpp>
#include <filesystem>
using aare::Cluster;
using aare::ClusterFile;
using aare::ClusterVector;
TEST_CASE("Read one frame from a a cluster file", "[.files]") {
TEST_CASE("Read one frame from a cluster file", "[.files]") {
//We know that the frame has 97 clusters
auto fpath = test_data_path() / "clust" / "single_frame_97_clustrers.clust";
REQUIRE(std::filesystem::exists(fpath));
ClusterFile f(fpath);
ClusterFile<Cluster<int32_t, 3, 3>> f(fpath);
auto clusters = f.read_frame();
REQUIRE(clusters.size() == 97);
REQUIRE(clusters.frame_number() == 135);
CHECK(clusters.size() == 97);
CHECK(clusters.frame_number() == 135);
CHECK(clusters[0].x == 1);
CHECK(clusters[0].y == 200);
int32_t expected_cluster_data[] = {0, 1, 2, 3, 4, 5, 6, 7, 8};
CHECK(std::equal(std::begin(clusters[0].data), std::end(clusters[0].data),
std::begin(expected_cluster_data)));
}
TEST_CASE("Read one frame using ROI", "[.files]") {
//We know that the frame has 97 clusters
// We know that the frame has 97 clusters
auto fpath = test_data_path() / "clust" / "single_frame_97_clustrers.clust";
REQUIRE(std::filesystem::exists(fpath));
ClusterFile f(fpath);
ClusterFile<Cluster<int32_t, 3, 3>> f(fpath);
aare::ROI roi;
roi.xmin = 0;
roi.xmax = 50;
@ -40,45 +44,308 @@ TEST_CASE("Read one frame using ROI", "[.files]") {
REQUIRE(clusters.size() == 49);
REQUIRE(clusters.frame_number() == 135);
//Check that all clusters are within the ROI
// Check that all clusters are within the ROI
for (size_t i = 0; i < clusters.size(); i++) {
auto c = clusters.at<aare::Cluster3x3>(i);
auto c = clusters[i];
REQUIRE(c.x >= roi.xmin);
REQUIRE(c.x <= roi.xmax);
REQUIRE(c.y >= roi.ymin);
REQUIRE(c.y <= roi.ymax);
}
CHECK(clusters[0].x == 1);
CHECK(clusters[0].y == 200);
int32_t expected_cluster_data[] = {0, 1, 2, 3, 4, 5, 6, 7, 8};
CHECK(std::equal(std::begin(clusters[0].data), std::end(clusters[0].data),
std::begin(expected_cluster_data)));
}
TEST_CASE("Read clusters from single frame file", "[.files]") {
// 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]
// [ 11 210] [90 91 92 93 94 95 96 97 98]
// [ 12 211] [ 99 100 101 102 103 104 105 106 107]
// [ 13 212] [108 109 110 111 112 113 114 115 116]
// [ 14 213] [117 118 119 120 121 122 123 124 125]
// [ 15 214] [126 127 128 129 130 131 132 133 134]
// [ 16 215] [135 136 137 138 139 140 141 142 143]
// [ 17 216] [144 145 146 147 148 149 150 151 152]
// [ 18 217] [153 154 155 156 157 158 159 160 161]
// [ 19 218] [162 163 164 165 166 167 168 169 170]
// [ 20 219] [171 172 173 174 175 176 177 178 179]
// [ 21 220] [180 181 182 183 184 185 186 187 188]
// [ 22 221] [189 190 191 192 193 194 195 196 197]
// [ 23 222] [198 199 200 201 202 203 204 205 206]
// [ 24 223] [207 208 209 210 211 212 213 214 215]
// [ 25 224] [216 217 218 219 220 221 222 223 224]
// [ 26 225] [225 226 227 228 229 230 231 232 233]
// [ 27 226] [234 235 236 237 238 239 240 241 242]
// [ 28 227] [243 244 245 246 247 248 249 250 251]
// [ 29 228] [252 253 254 255 256 257 258 259 260]
// [ 30 229] [261 262 263 264 265 266 267 268 269]
// [ 31 230] [270 271 272 273 274 275 276 277 278]
// [ 32 231] [279 280 281 282 283 284 285 286 287]
// [ 33 232] [288 289 290 291 292 293 294 295 296]
// [ 34 233] [297 298 299 300 301 302 303 304 305]
// [ 35 234] [306 307 308 309 310 311 312 313 314]
// [ 36 235] [315 316 317 318 319 320 321 322 323]
// [ 37 236] [324 325 326 327 328 329 330 331 332]
// [ 38 237] [333 334 335 336 337 338 339 340 341]
// [ 39 238] [342 343 344 345 346 347 348 349 350]
// [ 40 239] [351 352 353 354 355 356 357 358 359]
// [ 41 240] [360 361 362 363 364 365 366 367 368]
// [ 42 241] [369 370 371 372 373 374 375 376 377]
// [ 43 242] [378 379 380 381 382 383 384 385 386]
// [ 44 243] [387 388 389 390 391 392 393 394 395]
// [ 45 244] [396 397 398 399 400 401 402 403 404]
// [ 46 245] [405 406 407 408 409 410 411 412 413]
// [ 47 246] [414 415 416 417 418 419 420 421 422]
// [ 48 247] [423 424 425 426 427 428 429 430 431]
// [ 49 248] [432 433 434 435 436 437 438 439 440]
// [ 50 249] [441 442 443 444 445 446 447 448 449]
// [ 51 250] [450 451 452 453 454 455 456 457 458]
// [ 52 251] [459 460 461 462 463 464 465 466 467]
// [ 53 252] [468 469 470 471 472 473 474 475 476]
// [ 54 253] [477 478 479 480 481 482 483 484 485]
// [ 55 254] [486 487 488 489 490 491 492 493 494]
// [ 56 255] [495 496 497 498 499 500 501 502 503]
// [ 57 256] [504 505 506 507 508 509 510 511 512]
// [ 58 257] [513 514 515 516 517 518 519 520 521]
// [ 59 258] [522 523 524 525 526 527 528 529 530]
// [ 60 259] [531 532 533 534 535 536 537 538 539]
// [ 61 260] [540 541 542 543 544 545 546 547 548]
// [ 62 261] [549 550 551 552 553 554 555 556 557]
// [ 63 262] [558 559 560 561 562 563 564 565 566]
// [ 64 263] [567 568 569 570 571 572 573 574 575]
// [ 65 264] [576 577 578 579 580 581 582 583 584]
// [ 66 265] [585 586 587 588 589 590 591 592 593]
// [ 67 266] [594 595 596 597 598 599 600 601 602]
// [ 68 267] [603 604 605 606 607 608 609 610 611]
// [ 69 268] [612 613 614 615 616 617 618 619 620]
// [ 70 269] [621 622 623 624 625 626 627 628 629]
// [ 71 270] [630 631 632 633 634 635 636 637 638]
// [ 72 271] [639 640 641 642 643 644 645 646 647]
// [ 73 272] [648 649 650 651 652 653 654 655 656]
// [ 74 273] [657 658 659 660 661 662 663 664 665]
// [ 75 274] [666 667 668 669 670 671 672 673 674]
// [ 76 275] [675 676 677 678 679 680 681 682 683]
// [ 77 276] [684 685 686 687 688 689 690 691 692]
// [ 78 277] [693 694 695 696 697 698 699 700 701]
// [ 79 278] [702 703 704 705 706 707 708 709 710]
// [ 80 279] [711 712 713 714 715 716 717 718 719]
// [ 81 280] [720 721 722 723 724 725 726 727 728]
// [ 82 281] [729 730 731 732 733 734 735 736 737]
// [ 83 282] [738 739 740 741 742 743 744 745 746]
// [ 84 283] [747 748 749 750 751 752 753 754 755]
// [ 85 284] [756 757 758 759 760 761 762 763 764]
// [ 86 285] [765 766 767 768 769 770 771 772 773]
// [ 87 286] [774 775 776 777 778 779 780 781 782]
// [ 88 287] [783 784 785 786 787 788 789 790 791]
// [ 89 288] [792 793 794 795 796 797 798 799 800]
// [ 90 289] [801 802 803 804 805 806 807 808 809]
// [ 91 290] [810 811 812 813 814 815 816 817 818]
// [ 92 291] [819 820 821 822 823 824 825 826 827]
// [ 93 292] [828 829 830 831 832 833 834 835 836]
// [ 94 293] [837 838 839 840 841 842 843 844 845]
// [ 95 294] [846 847 848 849 850 851 852 853 854]
// [ 96 295] [855 856 857 858 859 860 861 862 863]
// [ 97 296] [864 865 866 867 868 869 870 871 872]
auto fpath = test_data_path() / "clust" / "single_frame_97_clustrers.clust";
REQUIRE(std::filesystem::exists(fpath));
SECTION("Read fewer clusters than available") {
ClusterFile f(fpath);
ClusterFile<Cluster<int32_t, 3, 3>> f(fpath);
auto clusters = f.read_clusters(50);
REQUIRE(clusters.size() == 50);
REQUIRE(clusters.frame_number() == 135);
REQUIRE(clusters.frame_number() == 135);
int32_t expected_cluster_data[] = {0, 1, 2, 3, 4, 5, 6, 7, 8};
REQUIRE(clusters[0].x == 1);
REQUIRE(clusters[0].y == 200);
CHECK(std::equal(std::begin(clusters[0].data),
std::end(clusters[0].data),
std::begin(expected_cluster_data)));
}
SECTION("Read more clusters than available") {
ClusterFile f(fpath);
ClusterFile<Cluster<int32_t, 3, 3>> f(fpath);
// 100 is the maximum number of clusters read
auto clusters = f.read_clusters(100);
REQUIRE(clusters.size() == 97);
REQUIRE(clusters.frame_number() == 135);
int32_t expected_cluster_data[] = {0, 1, 2, 3, 4, 5, 6, 7, 8};
REQUIRE(clusters[0].x == 1);
REQUIRE(clusters[0].y == 200);
CHECK(std::equal(std::begin(clusters[0].data),
std::end(clusters[0].data),
std::begin(expected_cluster_data)));
}
SECTION("Read all clusters") {
ClusterFile f(fpath);
ClusterFile<Cluster<int32_t, 3, 3>> f(fpath);
auto clusters = f.read_clusters(97);
REQUIRE(clusters.size() == 97);
REQUIRE(clusters.frame_number() == 135);
REQUIRE(clusters[0].x == 1);
REQUIRE(clusters[0].y == 200);
int32_t expected_cluster_data[] = {0, 1, 2, 3, 4, 5, 6, 7, 8};
CHECK(std::equal(std::begin(clusters[0].data),
std::end(clusters[0].data),
std::begin(expected_cluster_data)));
}
}
TEST_CASE("Read clusters from single frame file with ROI", "[.files]") {
auto fpath = test_data_path() / "clust" / "single_frame_97_clustrers.clust";
REQUIRE(std::filesystem::exists(fpath));
ClusterFile<Cluster<int32_t, 3, 3>> f(fpath);
aare::ROI roi;
roi.xmin = 0;
roi.xmax = 50;
roi.ymin = 200;
roi.ymax = 249;
f.set_roi(roi);
auto clusters = f.read_clusters(10);
CHECK(clusters.size() == 10);
CHECK(clusters.frame_number() == 135);
CHECK(clusters[0].x == 1);
CHECK(clusters[0].y == 200);
int32_t expected_cluster_data[] = {0, 1, 2, 3, 4, 5, 6, 7, 8};
CHECK(std::equal(std::begin(clusters[0].data), std::end(clusters[0].data),
std::begin(expected_cluster_data)));
}
TEST_CASE("Read cluster from multiple frame file", "[.files]") {
using ClusterType = Cluster<double, 2, 2>;
auto fpath =
test_data_path() / "clust" / "Two_frames_2x2double_test_clusters.clust";
REQUIRE(std::filesystem::exists(fpath));
// Two_frames_2x2double_test_clusters.clust
// frame number, num_clusters 0, 4
//[10, 20], {0. ,0., 0., 0.}
//[11, 30], {1., 1., 1., 1.}
//[12, 40], {2., 2., 2., 2.}
//[13, 50], {3., 3., 3., 3.}
// 1,4
//[10, 20], {4., 4., 4., 4.}
//[11, 30], {5., 5., 5., 5.}
//[12, 40], {6., 6., 6., 6.}
//[13, 50], {7., 7., 7., 7.}
SECTION("Read clusters from both frames") {
ClusterFile<ClusterType> f(fpath);
auto clusters = f.read_clusters(2);
REQUIRE(clusters.size() == 2);
REQUIRE(clusters.frame_number() == 0);
auto clusters1 = f.read_clusters(3);
REQUIRE(clusters1.size() == 3);
REQUIRE(clusters1.frame_number() == 1);
}
SECTION("Read all clusters") {
ClusterFile<ClusterType> f(fpath);
auto clusters = f.read_clusters(8);
REQUIRE(clusters.size() == 8);
REQUIRE(clusters.frame_number() == 1);
}
SECTION("Read clusters from one frame") {
ClusterFile<ClusterType> f(fpath);
auto clusters = f.read_clusters(2);
REQUIRE(clusters.size() == 2);
REQUIRE(clusters.frame_number() == 0);
auto clusters1 = f.read_clusters(1);
REQUIRE(clusters1.size() == 1);
REQUIRE(clusters1.frame_number() == 0);
}
}
TEST_CASE("Write cluster with potential padding", "[.files][.ClusterFile]") {
using ClusterType = Cluster<double, 3, 3>;
REQUIRE(std::filesystem::exists(test_data_path() / "clust"));
auto fpath = test_data_path() / "clust" / "single_frame_2_clusters.clust";
ClusterFile<ClusterType> file(fpath, 1000, "w");
ClusterVector<ClusterType> clustervec(2);
int16_t coordinate = 5;
clustervec.push_back(ClusterType{
coordinate, coordinate, {0., 0., 0., 0., 0., 0., 0., 0., 0.}});
clustervec.push_back(ClusterType{
coordinate, coordinate, {0., 0., 0., 0., 0., 0., 0., 0., 0.}});
file.write_frame(clustervec);
file.close();
file.open("r");
auto read_cluster_vector = file.read_frame();
CHECK(read_cluster_vector.size() == 2);
CHECK(read_cluster_vector.frame_number() == 0);
CHECK(read_cluster_vector[0].x == clustervec[0].x);
CHECK(read_cluster_vector[0].y == clustervec[0].y);
CHECK(std::equal(
clustervec[0].data.begin(), clustervec[0].data.end(),
read_cluster_vector[0].data.begin(), [](double a, double b) {
return std::abs(a - b) < std::numeric_limits<double>::epsilon();
}));
CHECK(read_cluster_vector[1].x == clustervec[1].x);
CHECK(read_cluster_vector[1].y == clustervec[1].y);
CHECK(std::equal(
clustervec[1].data.begin(), clustervec[1].data.end(),
read_cluster_vector[1].data.begin(), [](double a, double b) {
return std::abs(a - b) < std::numeric_limits<double>::epsilon();
}));
}
TEST_CASE("Read frame and modify cluster data", "[.files][.ClusterFile]") {
auto fpath = test_data_path() / "clust" / "single_frame_97_clustrers.clust";
REQUIRE(std::filesystem::exists(fpath));
ClusterFile<Cluster<int32_t, 3, 3>> f(fpath);
auto clusters = f.read_frame();
CHECK(clusters.size() == 97);
CHECK(clusters.frame_number() == 135);
int32_t expected_cluster_data[] = {0, 1, 2, 3, 4, 5, 6, 7, 8};
clusters.push_back(
Cluster<int32_t, 3, 3>{0, 0, {0, 1, 2, 3, 4, 5, 6, 7, 8}});
CHECK(clusters.size() == 98);
CHECK(clusters[0].x == 1);
CHECK(clusters[0].y == 200);
CHECK(std::equal(std::begin(clusters[0].data), std::end(clusters[0].data),
std::begin(expected_cluster_data)));
}

24
src/ClusterFinder.test.cpp
View File

@ -1,19 +1,18 @@
#include "aare/ClusterFinder.hpp"
#include "aare/Pedestal.hpp"
#include <catch2/matchers/catch_matchers_floating_point.hpp>
#include <catch2/catch_test_macros.hpp>
#include <catch2/matchers/catch_matchers_floating_point.hpp>
#include <chrono>
#include <random>
using namespace aare;
//TODO! Find a way to test the cluster finder
// TODO! Find a way to test the cluster finder
// class ClusterFinderUnitTest : public ClusterFinder {
// public:
// ClusterFinderUnitTest(int cluster_sizeX, int cluster_sizeY, double nSigma = 5.0, double threshold = 0.0)
// ClusterFinderUnitTest(int cluster_sizeX, int cluster_sizeY, double nSigma
// = 5.0, double threshold = 0.0)
// : ClusterFinder(cluster_sizeX, cluster_sizeY, nSigma, threshold) {}
// double get_c2() { return c2; }
// double get_c3() { return c3; }
@ -37,8 +36,8 @@ using namespace aare;
// REQUIRE_THAT(cf.get_c3(), Catch::Matchers::WithinRel(c3, 1e-9));
// }
TEST_CASE("Construct a cluster finder"){
ClusterFinder clusterFinder({400,400}, {3,3});
TEST_CASE("Construct a cluster finder") {
ClusterFinder clusterFinder({400, 400});
// REQUIRE(clusterFinder.get_cluster_sizeX() == 3);
// REQUIRE(clusterFinder.get_cluster_sizeY() == 3);
// REQUIRE(clusterFinder.get_threshold() == 1);
@ -49,16 +48,17 @@ TEST_CASE("Construct a cluster finder"){
// aare::Pedestal pedestal(10, 10, 5);
// NDArray<double, 2> frame({10, 10});
// frame = 0;
// ClusterFinder clusterFinder(3, 3, 1, 1); // 3x3 cluster, 1 nSigma, 1 threshold
// ClusterFinder clusterFinder(3, 3, 1, 1); // 3x3 cluster, 1 nSigma, 1
// threshold
// auto clusters = clusterFinder.find_clusters_without_threshold(frame.span(), pedestal);
// auto clusters =
// clusterFinder.find_clusters_without_threshold(frame.span(), pedestal);
// REQUIRE(clusters.size() == 0);
// frame(5, 5) = 10;
// clusters = clusterFinder.find_clusters_without_threshold(frame.span(), pedestal);
// REQUIRE(clusters.size() == 1);
// REQUIRE(clusters[0].x == 5);
// clusters = clusterFinder.find_clusters_without_threshold(frame.span(),
// pedestal); REQUIRE(clusters.size() == 1); REQUIRE(clusters[0].x == 5);
// REQUIRE(clusters[0].y == 5);
// for (int i = 0; i < 3; i++) {
// for (int j = 0; j < 3; j++) {

99
src/ClusterFinderMT.test.cpp Normal file
View File

@ -0,0 +1,99 @@
#include "aare/ClusterFinderMT.hpp"
#include "aare/Cluster.hpp"
#include "aare/ClusterCollector.hpp"
#include "aare/File.hpp"
#include "test_config.hpp"
#include <catch2/catch_test_macros.hpp>
#include <filesystem>
#include <memory>
using namespace aare;
// wrapper function to access private member variables for testing
template <typename ClusterType, typename FRAME_TYPE = uint16_t,
typename PEDESTAL_TYPE = double>
class ClusterFinderMTWrapper
: public ClusterFinderMT<ClusterType, FRAME_TYPE, PEDESTAL_TYPE> {
public:
ClusterFinderMTWrapper(Shape<2> image_size, PEDESTAL_TYPE nSigma = 5.0,
size_t capacity = 2000, size_t n_threads = 3)
: ClusterFinderMT<ClusterType, FRAME_TYPE, PEDESTAL_TYPE>(
image_size, nSigma, capacity, n_threads) {}
size_t get_m_input_queues_size() const {
return this->m_input_queues.size();
}
size_t get_m_output_queues_size() const {
return this->m_output_queues.size();
}
size_t get_m_cluster_finders_size() const {
return this->m_cluster_finders.size();
}
bool m_output_queues_are_empty() const {
for (auto &queue : this->m_output_queues) {
if (!queue->isEmpty())
return false;
}
return true;
}
bool m_input_queues_are_empty() const {
for (auto &queue : this->m_input_queues) {
if (!queue->isEmpty())
return false;
}
return true;
}
bool m_sink_is_empty() const { return this->m_sink.isEmpty(); }
size_t m_sink_size() const { return this->m_sink.sizeGuess(); }
};
TEST_CASE("multithreaded cluster finder", "[.files][.ClusterFinder]") {
auto fpath = "/mnt/sls_det_storage/matterhorn_data/aare_test_data/"
"Moench03new/cu_half_speed_master_4.json";
File file(fpath);
size_t n_threads = 2;
size_t n_frames_pd = 10;
using ClusterType = Cluster<int32_t, 3, 3>;
ClusterFinderMTWrapper<ClusterType> cf(
{static_cast<int64_t>(file.rows()), static_cast<int64_t>(file.cols())},
5, 2000, n_threads); // no idea what frame type is!!! default uint16_t
CHECK(cf.get_m_input_queues_size() == n_threads);
CHECK(cf.get_m_output_queues_size() == n_threads);
CHECK(cf.get_m_cluster_finders_size() == n_threads);
CHECK(cf.m_output_queues_are_empty() == true);
CHECK(cf.m_input_queues_are_empty() == true);
for (size_t i = 0; i < n_frames_pd; ++i) {
cf.find_clusters(file.read_frame().view<uint16_t>());
}
cf.stop();
CHECK(cf.m_output_queues_are_empty() == true);
CHECK(cf.m_input_queues_are_empty() == true);
CHECK(cf.m_sink_size() == n_frames_pd);
ClusterCollector<ClusterType> clustercollector(&cf);
clustercollector.stop();
CHECK(cf.m_sink_size() == 0);
auto clustervec = clustercollector.steal_clusters();
// CHECK(clustervec.size() == ) //dont know how many clusters to expect
}

240
src/ClusterVector.test.cpp
View File

@ -1,21 +1,52 @@
#include <cstdint>
#include "aare/ClusterVector.hpp"
#include <cstdint>
#include <catch2/matchers/catch_matchers_floating_point.hpp>
#include <catch2/catch_all.hpp>
#include <catch2/catch_test_macros.hpp>
#include <catch2/matchers/catch_matchers_floating_point.hpp>
using aare::Cluster;
using aare::ClusterVector;
struct Cluster_i2x2 {
int16_t x;
int16_t y;
int32_t data[4];
};
TEST_CASE("item_size return the size of the cluster stored") {
using C1 = Cluster<int32_t, 2, 2>;
ClusterVector<C1> cv(4);
CHECK(cv.item_size() == sizeof(C1));
TEST_CASE("ClusterVector 2x2 int32_t capacity 4, push back then read") {
// Sanity check
// 2*2*4 = 16 bytes of data for the cluster
// 2*2 = 4 bytes for the x and y coordinates
REQUIRE(cv.item_size() == 20);
ClusterVector<int32_t> cv(2, 2, 4);
using C2 = Cluster<int32_t, 3, 3>;
ClusterVector<C2> cv2(4);
CHECK(cv2.item_size() == sizeof(C2));
using C3 = Cluster<double, 2, 3>;
ClusterVector<C3> cv3(4);
CHECK(cv3.item_size() == sizeof(C3));
using C4 = Cluster<char, 10, 5>;
ClusterVector<C4> cv4(4);
CHECK(cv4.item_size() == sizeof(C4));
using C5 = Cluster<int32_t, 2, 3>;
ClusterVector<C5> cv5(4);
CHECK(cv5.item_size() == sizeof(C5));
using C6 = Cluster<double, 5, 5>;
ClusterVector<C6> cv6(4);
CHECK(cv6.item_size() == sizeof(C6));
using C7 = Cluster<double, 3, 3>;
ClusterVector<C7> cv7(4);
CHECK(cv7.item_size() == sizeof(C7));
}
TEST_CASE("ClusterVector 2x2 int32_t capacity 4, push back then read",
"[.ClusterVector]") {
ClusterVector<Cluster<int32_t, 2, 2>> cv(4);
REQUIRE(cv.capacity() == 4);
REQUIRE(cv.size() == 0);
REQUIRE(cv.cluster_size_x() == 2);
@ -23,112 +54,102 @@ TEST_CASE("ClusterVector 2x2 int32_t capacity 4, push back then read") {
// int16_t, int16_t, 2x2 int32_t = 20 bytes
REQUIRE(cv.item_size() == 20);
//Create a cluster and push back into the vector
Cluster_i2x2 c1 = {1, 2, {3, 4, 5, 6}};
cv.push_back(c1.x, c1.y, reinterpret_cast<std::byte*>(&c1.data[0]));
// Create a cluster and push back into the vector
Cluster<int32_t, 2, 2> c1 = {1, 2, {3, 4, 5, 6}};
cv.push_back(c1);
REQUIRE(cv.size() == 1);
REQUIRE(cv.capacity() == 4);
//Read the cluster back out using copy. TODO! Can we improve the API?
Cluster_i2x2 c2;
std::byte *ptr = cv.element_ptr(0);
std::copy(ptr, ptr + cv.item_size(), reinterpret_cast<std::byte*>(&c2));
auto c2 = cv[0];
//Check that the data is the same
// Check that the data is the same
REQUIRE(c1.x == c2.x);
REQUIRE(c1.y == c2.y);
for(size_t i = 0; i < 4; i++) {
for (size_t i = 0; i < 4; i++) {
REQUIRE(c1.data[i] == c2.data[i]);
}
}
TEST_CASE("Summing 3x1 clusters of int64"){
struct Cluster_l3x1{
int16_t x;
int16_t y;
int32_t data[3];
};
ClusterVector<int32_t> cv(3, 1, 2);
TEST_CASE("Summing 3x1 clusters of int64", "[.ClusterVector]") {
ClusterVector<Cluster<int32_t, 3, 1>> cv(2);
REQUIRE(cv.capacity() == 2);
REQUIRE(cv.size() == 0);
REQUIRE(cv.cluster_size_x() == 3);
REQUIRE(cv.cluster_size_y() == 1);
//Create a cluster and push back into the vector
Cluster_l3x1 c1 = {1, 2, {3, 4, 5}};
cv.push_back(c1.x, c1.y, reinterpret_cast<std::byte*>(&c1.data[0]));
// Create a cluster and push back into the vector
Cluster<int32_t, 3, 1> c1 = {1, 2, {3, 4, 5}};
cv.push_back(c1);
REQUIRE(cv.capacity() == 2);
REQUIRE(cv.size() == 1);
Cluster_l3x1 c2 = {6, 7, {8, 9, 10}};
cv.push_back(c2.x, c2.y, reinterpret_cast<std::byte*>(&c2.data[0]));
Cluster<int32_t, 3, 1> c2 = {6, 7, {8, 9, 10}};
cv.push_back(c2);
REQUIRE(cv.capacity() == 2);
REQUIRE(cv.size() == 2);
Cluster_l3x1 c3 = {11, 12, {13, 14, 15}};
cv.push_back(c3.x, c3.y, reinterpret_cast<std::byte*>(&c3.data[0]));
Cluster<int32_t, 3, 1> c3 = {11, 12, {13, 14, 15}};
cv.push_back(c3);
REQUIRE(cv.capacity() == 4);
REQUIRE(cv.size() == 3);
/*
auto sums = cv.sum();
REQUIRE(sums.size() == 3);
REQUIRE(sums[0] == 12);
REQUIRE(sums[1] == 27);
REQUIRE(sums[2] == 42);
*/
}
TEST_CASE("Storing floats"){
struct Cluster_f4x2{
int16_t x;
int16_t y;
float data[8];
};
ClusterVector<float> cv(2, 4, 10);
TEST_CASE("Storing floats", "[.ClusterVector]") {
ClusterVector<Cluster<float, 2, 4>> cv(10);
REQUIRE(cv.capacity() == 10);
REQUIRE(cv.size() == 0);
REQUIRE(cv.cluster_size_x() == 2);
REQUIRE(cv.cluster_size_y() == 4);
//Create a cluster and push back into the vector
Cluster_f4x2 c1 = {1, 2, {3.0, 4.0, 5.0, 6.0,3.0, 4.0, 5.0, 6.0}};
cv.push_back(c1.x, c1.y, reinterpret_cast<std::byte*>(&c1.data[0]));
// Create a cluster and push back into the vector
Cluster<float, 2, 4> c1 = {1, 2, {3.0, 4.0, 5.0, 6.0, 3.0, 4.0, 5.0, 6.0}};
cv.push_back(c1);
REQUIRE(cv.capacity() == 10);
REQUIRE(cv.size() == 1);
Cluster_f4x2 c2 = {6, 7, {8.0, 9.0, 10.0, 11.0,8.0, 9.0, 10.0, 11.0}};
cv.push_back(c2.x, c2.y, reinterpret_cast<std::byte*>(&c2.data[0]));
Cluster<float, 2, 4> c2 = {
6, 7, {8.0, 9.0, 10.0, 11.0, 8.0, 9.0, 10.0, 11.0}};
cv.push_back(c2);
REQUIRE(cv.capacity() == 10);
REQUIRE(cv.size() == 2);
/*
auto sums = cv.sum();
REQUIRE(sums.size() == 2);
REQUIRE_THAT(sums[0], Catch::Matchers::WithinAbs(36.0, 1e-6));
REQUIRE_THAT(sums[1], Catch::Matchers::WithinAbs(76.0, 1e-6));
*/
}
TEST_CASE("Push back more than initial capacity"){
ClusterVector<int32_t> cv(2, 2, 2);
TEST_CASE("Push back more than initial capacity", "[.ClusterVector]") {
ClusterVector<Cluster<int32_t, 2, 2>> cv(2);
auto initial_data = cv.data();
Cluster_i2x2 c1 = {1, 2, {3, 4, 5, 6}};
cv.push_back(c1.x, c1.y, reinterpret_cast<std::byte*>(&c1.data[0]));
Cluster<int32_t, 2, 2> c1 = {1, 2, {3, 4, 5, 6}};
cv.push_back(c1);
REQUIRE(cv.size() == 1);
REQUIRE(cv.capacity() == 2);
Cluster_i2x2 c2 = {6, 7, {8, 9, 10, 11}};
cv.push_back(c2.x, c2.y, reinterpret_cast<std::byte*>(&c2.data[0]));
Cluster<int32_t, 2, 2> c2 = {6, 7, {8, 9, 10, 11}};
cv.push_back(c2);
REQUIRE(cv.size() == 2);
REQUIRE(cv.capacity() == 2);
Cluster_i2x2 c3 = {11, 12, {13, 14, 15, 16}};
cv.push_back(c3.x, c3.y, reinterpret_cast<std::byte*>(&c3.data[0]));
REQUIRE(cv.size() == 3);
Cluster<int32_t, 2, 2> c3 = {11, 12, {13, 14, 15, 16}};
cv.push_back(c3);
REQUIRE(cv.size() == 3);
REQUIRE(cv.capacity() == 4);
Cluster_i2x2* ptr = reinterpret_cast<Cluster_i2x2*>(cv.data());
Cluster<int32_t, 2, 2> *ptr =
reinterpret_cast<Cluster<int32_t, 2, 2> *>(cv.data());
REQUIRE(ptr[0].x == 1);
REQUIRE(ptr[0].y == 2);
REQUIRE(ptr[1].x == 6);
@ -136,29 +157,31 @@ TEST_CASE("Push back more than initial capacity"){
REQUIRE(ptr[2].x == 11);
REQUIRE(ptr[2].y == 12);
//We should have allocated a new buffer, since we outgrew the initial capacity
// We should have allocated a new buffer, since we outgrew the initial
// capacity
REQUIRE(initial_data != cv.data());
}
TEST_CASE("Concatenate two cluster vectors where the first has enough capacity"){
ClusterVector<int32_t> cv1(2, 2, 12);
Cluster_i2x2 c1 = {1, 2, {3, 4, 5, 6}};
cv1.push_back(c1.x, c1.y, reinterpret_cast<std::byte*>(&c1.data[0]));
Cluster_i2x2 c2 = {6, 7, {8, 9, 10, 11}};
cv1.push_back(c2.x, c2.y, reinterpret_cast<std::byte*>(&c2.data[0]));
TEST_CASE("Concatenate two cluster vectors where the first has enough capacity",
"[.ClusterVector]") {
ClusterVector<Cluster<int32_t, 2, 2>> cv1(12);
Cluster<int32_t, 2, 2> c1 = {1, 2, {3, 4, 5, 6}};
cv1.push_back(c1);
Cluster<int32_t, 2, 2> c2 = {6, 7, {8, 9, 10, 11}};
cv1.push_back(c2);
ClusterVector<int32_t> cv2(2, 2, 2);
Cluster_i2x2 c3 = {11, 12, {13, 14, 15, 16}};
cv2.push_back(c3.x, c3.y, reinterpret_cast<std::byte*>(&c3.data[0]));
Cluster_i2x2 c4 = {16, 17, {18, 19, 20, 21}};
cv2.push_back(c4.x, c4.y, reinterpret_cast<std::byte*>(&c4.data[0]));
ClusterVector<Cluster<int32_t, 2, 2>> cv2(2);
Cluster<int32_t, 2, 2> c3 = {11, 12, {13, 14, 15, 16}};
cv2.push_back(c3);
Cluster<int32_t, 2, 2> c4 = {16, 17, {18, 19, 20, 21}};
cv2.push_back(c4);
cv1 += cv2;
REQUIRE(cv1.size() == 4);
REQUIRE(cv1.capacity() == 12);
Cluster_i2x2* ptr = reinterpret_cast<Cluster_i2x2*>(cv1.data());
Cluster<int32_t, 2, 2> *ptr =
reinterpret_cast<Cluster<int32_t, 2, 2> *>(cv1.data());
REQUIRE(ptr[0].x == 1);
REQUIRE(ptr[0].y == 2);
REQUIRE(ptr[1].x == 6);
@ -169,24 +192,26 @@ TEST_CASE("Concatenate two cluster vectors where the first has enough capacity")
REQUIRE(ptr[3].y == 17);
}
TEST_CASE("Concatenate two cluster vectors where we need to allocate"){
ClusterVector<int32_t> cv1(2, 2, 2);
Cluster_i2x2 c1 = {1, 2, {3, 4, 5, 6}};
cv1.push_back(c1.x, c1.y, reinterpret_cast<std::byte*>(&c1.data[0]));
Cluster_i2x2 c2 = {6, 7, {8, 9, 10, 11}};
cv1.push_back(c2.x, c2.y, reinterpret_cast<std::byte*>(&c2.data[0]));
TEST_CASE("Concatenate two cluster vectors where we need to allocate",
"[.ClusterVector]") {
ClusterVector<Cluster<int32_t, 2, 2>> cv1(2);
Cluster<int32_t, 2, 2> c1 = {1, 2, {3, 4, 5, 6}};
cv1.push_back(c1);
Cluster<int32_t, 2, 2> c2 = {6, 7, {8, 9, 10, 11}};
cv1.push_back(c2);
ClusterVector<int32_t> cv2(2, 2, 2);
Cluster_i2x2 c3 = {11, 12, {13, 14, 15, 16}};
cv2.push_back(c3.x, c3.y, reinterpret_cast<std::byte*>(&c3.data[0]));
Cluster_i2x2 c4 = {16, 17, {18, 19, 20, 21}};
cv2.push_back(c4.x, c4.y, reinterpret_cast<std::byte*>(&c4.data[0]));
ClusterVector<Cluster<int32_t, 2, 2>> cv2(2);
Cluster<int32_t, 2, 2> c3 = {11, 12, {13, 14, 15, 16}};
cv2.push_back(c3);
Cluster<int32_t, 2, 2> c4 = {16, 17, {18, 19, 20, 21}};
cv2.push_back(c4);
cv1 += cv2;
REQUIRE(cv1.size() == 4);
REQUIRE(cv1.capacity() == 4);
Cluster_i2x2* ptr = reinterpret_cast<Cluster_i2x2*>(cv1.data());
Cluster<int32_t, 2, 2> *ptr =
reinterpret_cast<Cluster<int32_t, 2, 2> *>(cv1.data());
REQUIRE(ptr[0].x == 1);
REQUIRE(ptr[0].y == 2);
REQUIRE(ptr[1].x == 6);
@ -195,4 +220,49 @@ TEST_CASE("Concatenate two cluster vectors where we need to allocate"){
REQUIRE(ptr[2].y == 12);
REQUIRE(ptr[3].x == 16);
REQUIRE(ptr[3].y == 17);
}
struct ClusterTestData {
uint8_t ClusterSizeX;
uint8_t ClusterSizeY;
std::vector<int64_t> index_map_x;
std::vector<int64_t> index_map_y;
};
TEST_CASE("Gain Map Calculation Index Map", "[.ClusterVector][.gain_map]") {
auto clustertestdata = GENERATE(
ClusterTestData{3,
3,
{-1, 0, 1, -1, 0, 1, -1, 0, 1},
{-1, -1, -1, 0, 0, 0, 1, 1, 1}},
ClusterTestData{
4,
4,
{-2, -1, 0, 1, -2, -1, 0, 1, -2, -1, 0, 1, -2, -1, 0, 1},
{-2, -2, -2, -2, -1, -1, -1, -1, 0, 0, 0, 0, 1, 1, 1, 1}},
ClusterTestData{2, 2, {-1, 0, -1, 0}, {-1, -1, 0, 0}},
ClusterTestData{5,
5,
{-2, -1, 0, 1, 2, -2, -1, 0, 1, 2, -2, -1, 0,
1, 2, -2, -1, 0, 1, 2, -2, -1, 0, 1, 2},
{-2, -2, -2, -2, -2, -1, -1, -1, -1, -1, 0, 0, 0,
0, 0, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2}});
uint8_t ClusterSizeX = clustertestdata.ClusterSizeX;
uint8_t ClusterSizeY = clustertestdata.ClusterSizeY;
std::vector<int64_t> index_map_x(ClusterSizeX * ClusterSizeY);
std::vector<int64_t> index_map_y(ClusterSizeX * ClusterSizeY);
int64_t index_cluster_center_x = ClusterSizeX / 2;
int64_t index_cluster_center_y = ClusterSizeY / 2;
for (size_t j = 0; j < ClusterSizeX * ClusterSizeY; j++) {
index_map_x[j] = j % ClusterSizeX - index_cluster_center_x;
index_map_y[j] = j / ClusterSizeX - index_cluster_center_y;
}
CHECK(index_map_x == clustertestdata.index_map_x);
CHECK(index_map_y == clustertestdata.index_map_y);
}

2
src/FilePtr.cpp
View File

@ -21,7 +21,7 @@ FilePtr &FilePtr::operator=(FilePtr &&other) {
FILE *FilePtr::get() { return fp_; }
int64_t FilePtr::tell() {
ssize_t FilePtr::tell() {
auto pos = ftell(fp_);
if (pos == -1)
throw std::runtime_error(fmt::format("Error getting file position: {}", error_msg()));

251
src/Fit.cpp
View File

@ -34,6 +34,30 @@ NDArray<double, 1> pol1(NDView<double, 1> x, NDView<double, 1> par) {
return y;
}
double scurve(const double x, const double * par) {
return (par[0] + par[1] * x) + 0.5 * (1 + erf((x - par[2]) / (sqrt(2) * par[3]))) * (par[4] + par[5] * (x - par[2]));
}
NDArray<double, 1> scurve(NDView<double, 1> x, NDView<double, 1> par) {
NDArray<double, 1> y({x.shape()}, 0);
for (ssize_t i = 0; i < x.size(); i++) {
y(i) = scurve(x(i), par.data());
}
return y;
}
double scurve2(const double x, const double * par) {
return (par[0] + par[1] * x) + 0.5 * (1 - erf((x - par[2]) / (sqrt(2) * par[3]))) * (par[4] + par[5] * (x - par[2]));
}
NDArray<double, 1> scurve2(NDView<double, 1> x, NDView<double, 1> par) {
NDArray<double, 1> y({x.shape()}, 0);
for (ssize_t i = 0; i < x.size(); i++) {
y(i) = scurve2(x(i), par.data());
}
return y;
}
} // namespace func
NDArray<double, 1> fit_gaus(NDView<double, 1> x, NDView<double, 1> y) {
@ -81,7 +105,7 @@ std::array<double, 3> gaus_init_par(const NDView<double, 1> x, const NDView<doub
auto delta = x[1] - x[0];
start_par[2] =
std::count_if(y.begin(), y.end(),
[e, delta](double val) { return val > *e / 2; }) *
[e](double val) { return val > *e / 2; }) *
delta / 2.35;
return start_par;
@ -273,4 +297,229 @@ NDArray<double, 3> fit_pol1(NDView<double, 1> x, NDView<double, 3> y,
return result;
}
// ~~ S-CURVES ~~
// SCURVE --
std::array<double, 6> scurve_init_par(const NDView<double, 1> x, const NDView<double, 1> y){
// Estimate the initial parameters for the fit
std::array<double, 6> start_par{0, 0, 0, 0, 0, 0};
auto ymax = std::max_element(y.begin(), y.end());
auto ymin = std::min_element(y.begin(), y.end());
start_par[4] = *ymin + (*ymax - *ymin) / 2;
// Find the first x where the corresponding y value is above the threshold (start_par[4])
for (ssize_t i = 0; i < y.size(); ++i) {
if (y[i] >= start_par[4]) {
start_par[2] = x[i];
break; // Exit the loop after finding the first valid x
}
}
start_par[3] = 2 * sqrt(start_par[2]);
start_par[0] = 100;
start_par[1] = 0.25;
start_par[5] = 1;
return start_par;
}
// - No error
NDArray<double, 1> fit_scurve(NDView<double, 1> x, NDView<double, 1> y) {
NDArray<double, 1> result = scurve_init_par(x, y);
lm_status_struct status;
lmcurve(result.size(), result.data(), x.size(), x.data(), y.data(),
aare::func::scurve, &lm_control_double, &status);
return result;
}
NDArray<double, 3> fit_scurve(NDView<double, 1> x, NDView<double, 3> y, int n_threads) {
NDArray<double, 3> result({y.shape(0), y.shape(1), 6}, 0);
auto process = [&x, &y, &result](ssize_t first_row, ssize_t last_row) {
for (ssize_t row = first_row; row < last_row; row++) {
for (ssize_t col = 0; col < y.shape(1); col++) {
NDView<double, 1> values(&y(row, col, 0), {y.shape(2)});
auto res = fit_scurve(x, values);
result(row, col, 0) = res(0);
result(row, col, 1) = res(1);
result(row, col, 2) = res(2);
result(row, col, 3) = res(3);
result(row, col, 4) = res(4);
result(row, col, 5) = res(5);
}
}
};
auto tasks = split_task(0, y.shape(0), n_threads);
RunInParallel(process, tasks);
return result;
}
// - Error
void fit_scurve(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
NDView<double, 1> par_out, NDView<double, 1> par_err_out, double& chi2) {
// Check that we have the correct sizes
if (y.size() != x.size() || y.size() != y_err.size() ||
par_out.size() != 6 || par_err_out.size() != 6) {
throw std::runtime_error("Data, x, data_err must have the same size "
"and par_out, par_err_out must have size 6");
}
lm_status_struct status;
par_out = scurve_init_par(x, y);
std::array<double, 36> cov = {0}; // size 6x6
// std::array<double, 4> cov{0, 0, 0, 0};
lmcurve2(par_out.size(), par_out.data(), par_err_out.data(), cov.data(),
x.size(), x.data(), y.data(), y_err.data(), aare::func::scurve,
&lm_control_double, &status);
// Calculate chi2
chi2 = 0;
for (ssize_t i = 0; i < y.size(); i++) {
chi2 += std::pow((y(i) - func::pol1(x(i), par_out.data())) / y_err(i), 2);
}
}
void fit_scurve(NDView<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_err,
NDView<double, 3> par_out, NDView<double, 3> par_err_out, NDView<double, 2> chi2_out,
int n_threads) {
auto process = [&](ssize_t first_row, ssize_t last_row) {
for (ssize_t row = first_row; row < last_row; row++) {
for (ssize_t col = 0; col < y.shape(1); col++) {
NDView<double, 1> y_view(&y(row, col, 0), {y.shape(2)});
NDView<double, 1> y_err_view(&y_err(row, col, 0),
{y_err.shape(2)});
NDView<double, 1> par_out_view(&par_out(row, col, 0),
{par_out.shape(2)});
NDView<double, 1> par_err_out_view(&par_err_out(row, col, 0),
{par_err_out.shape(2)});
fit_scurve(x, y_view, y_err_view, par_out_view, par_err_out_view, chi2_out(row, col));
}
}
};
auto tasks = split_task(0, y.shape(0), n_threads);
RunInParallel(process, tasks);
}
// SCURVE2 ---
std::array<double, 6> scurve2_init_par(const NDView<double, 1> x, const NDView<double, 1> y){
// Estimate the initial parameters for the fit
std::array<double, 6> start_par{0, 0, 0, 0, 0, 0};
auto ymax = std::max_element(y.begin(), y.end());
auto ymin = std::min_element(y.begin(), y.end());
start_par[4] = *ymin + (*ymax - *ymin) / 2;
// Find the first x where the corresponding y value is above the threshold (start_par[4])
for (ssize_t i = 0; i < y.size(); ++i) {
if (y[i] <= start_par[4]) {
start_par[2] = x[i];
break; // Exit the loop after finding the first valid x
}
}
start_par[3] = 2 * sqrt(start_par[2]);
start_par[0] = 100;
start_par[1] = 0.25;
start_par[5] = -1;
return start_par;
}
// - No error
NDArray<double, 1> fit_scurve2(NDView<double, 1> x, NDView<double, 1> y) {
NDArray<double, 1> result = scurve2_init_par(x, y);
lm_status_struct status;
lmcurve(result.size(), result.data(), x.size(), x.data(), y.data(),
aare::func::scurve2, &lm_control_double, &status);
return result;
}
NDArray<double, 3> fit_scurve2(NDView<double, 1> x, NDView<double, 3> y, int n_threads) {
NDArray<double, 3> result({y.shape(0), y.shape(1), 6}, 0);
auto process = [&x, &y, &result](ssize_t first_row, ssize_t last_row) {
for (ssize_t row = first_row; row < last_row; row++) {
for (ssize_t col = 0; col < y.shape(1); col++) {
NDView<double, 1> values(&y(row, col, 0), {y.shape(2)});
auto res = fit_scurve2(x, values);
result(row, col, 0) = res(0);
result(row, col, 1) = res(1);
result(row, col, 2) = res(2);
result(row, col, 3) = res(3);
result(row, col, 4) = res(4);
result(row, col, 5) = res(5);
}
}
};
auto tasks = split_task(0, y.shape(0), n_threads);
RunInParallel(process, tasks);
return result;
}
// - Error
void fit_scurve2(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
NDView<double, 1> par_out, NDView<double, 1> par_err_out, double& chi2) {
// Check that we have the correct sizes
if (y.size() != x.size() || y.size() != y_err.size() ||
par_out.size() != 6 || par_err_out.size() != 6) {
throw std::runtime_error("Data, x, data_err must have the same size "
"and par_out, par_err_out must have size 6");
}
lm_status_struct status;
par_out = scurve2_init_par(x, y);
std::array<double, 36> cov = {0}; // size 6x6
// std::array<double, 4> cov{0, 0, 0, 0};
lmcurve2(par_out.size(), par_out.data(), par_err_out.data(), cov.data(),
x.size(), x.data(), y.data(), y_err.data(), aare::func::scurve2,
&lm_control_double, &status);
// Calculate chi2
chi2 = 0;
for (ssize_t i = 0; i < y.size(); i++) {
chi2 += std::pow((y(i) - func::pol1(x(i), par_out.data())) / y_err(i), 2);
}
}
void fit_scurve2(NDView<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_err,
NDView<double, 3> par_out, NDView<double, 3> par_err_out, NDView<double, 2> chi2_out,
int n_threads) {
auto process = [&](ssize_t first_row, ssize_t last_row) {
for (ssize_t row = first_row; row < last_row; row++) {
for (ssize_t col = 0; col < y.shape(1); col++) {
NDView<double, 1> y_view(&y(row, col, 0), {y.shape(2)});
NDView<double, 1> y_err_view(&y_err(row, col, 0),
{y_err.shape(2)});
NDView<double, 1> par_out_view(&par_out(row, col, 0),
{par_out.shape(2)});
NDView<double, 1> par_err_out_view(&par_err_out(row, col, 0),
{par_err_out.shape(2)});
fit_scurve2(x, y_view, y_err_view, par_out_view, par_err_out_view, chi2_out(row, col));
}
}
};
auto tasks = split_task(0, y.shape(0), n_threads);
RunInParallel(process, tasks);
}
} // namespace aare

87
src/Interpolator.cpp
View File

@ -1,11 +1,11 @@
#include "aare/Interpolator.hpp"
#include "aare/algorithm.hpp"
namespace aare {
Interpolator::Interpolator(NDView<double, 3> etacube, NDView<double, 1> xbins,
NDView<double, 1> ybins, NDView<double, 1> ebins)
: m_ietax(etacube), m_ietay(etacube), m_etabinsx(xbins), m_etabinsy(ybins), m_energy_bins(ebins) {
: m_ietax(etacube), m_ietay(etacube), m_etabinsx(xbins), m_etabinsy(ybins),
m_energy_bins(ebins) {
if (etacube.shape(0) != xbins.size() || etacube.shape(1) != ybins.size() ||
etacube.shape(2) != ebins.size()) {
throw std::invalid_argument(
@ -53,87 +53,4 @@ Interpolator::Interpolator(NDView<double, 3> etacube, NDView<double, 1> xbins,
}
}
std::vector<Photon> Interpolator::interpolate(const ClusterVector<int32_t>& clusters) {
std::vector<Photon> photons;
photons.reserve(clusters.size());
if (clusters.cluster_size_x() == 3 || clusters.cluster_size_y() == 3) {
for (size_t i = 0; i<clusters.size(); i++){
auto cluster = clusters.at<Cluster3x3>(i);
Eta2 eta= calculate_eta2(cluster);
Photon photon;
photon.x = cluster.x;
photon.y = cluster.y;
photon.energy = eta.sum;
//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);
double dX{}, dY{};
// cBottomLeft = 0,
// cBottomRight = 1,
// cTopLeft = 2,
// cTopRight = 3
switch (eta.c) {
case cTopLeft:
dX = -1.;
dY = 0.;
break;
case cTopRight:;
dX = 0.;
dY = 0.;
break;
case cBottomLeft:
dX = -1.;
dY = -1.;
break;
case 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 (size_t i = 0; i<clusters.size(); i++){
auto cluster = clusters.at<Cluster2x2>(i);
Eta2 eta= calculate_eta2(cluster);
Photon photon;
photon.x = cluster.x;
photon.y = cluster.y;
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

4
src/JungfrauDataFile.cpp
View File

@ -89,7 +89,7 @@ void JungfrauDataFile::seek(size_t frame_index) {
: frame_index;
auto byte_offset = frame_offset * (m_bytes_per_frame + header_size);
m_fp.seek(byte_offset);
};
}
size_t JungfrauDataFile::tell() { return m_current_frame_index; }
size_t JungfrauDataFile::total_frames() const { return m_total_frames; }
@ -235,4 +235,4 @@ std::filesystem::path JungfrauDataFile::fpath(size_t file_index) const {
return m_path / fname;
}
} // namespace aare
} // namespace aare

20
src/NDArray.test.cpp
View File

@ -44,9 +44,9 @@ TEST_CASE("3D NDArray from NDView"){
REQUIRE(image.size() == view.size());
REQUIRE(image.data() != view.data());
for(int64_t i=0; i<image.shape(0); i++){
for(int64_t j=0; j<image.shape(1); j++){
for(int64_t k=0; k<image.shape(2); k++){
for(ssize_t i=0; i<image.shape(0); i++){
for(ssize_t j=0; j<image.shape(1); j++){
for(ssize_t k=0; k<image.shape(2); k++){
REQUIRE(image(i, j, k) == view(i, j, k));
}
}
@ -54,7 +54,7 @@ TEST_CASE("3D NDArray from NDView"){
}
TEST_CASE("1D image") {
std::array<int64_t, 1> shape{{20}};
std::array<ssize_t, 1> shape{{20}};
NDArray<short, 1> img(shape, 3);
REQUIRE(img.size() == 20);
REQUIRE(img(5) == 3);
@ -71,7 +71,7 @@ TEST_CASE("Accessing a const object") {
}
TEST_CASE("Indexing of a 2D image") {
std::array<int64_t, 2> shape{{3, 7}};
std::array<ssize_t, 2> shape{{3, 7}};
NDArray<long> img(shape, 5);
for (uint32_t i = 0; i != img.size(); ++i) {
REQUIRE(img(i) == 5);
@ -114,7 +114,7 @@ TEST_CASE("Divide double by int") {
}
TEST_CASE("Elementwise multiplication of 3D image") {
std::array<int64_t, 3> shape{3, 4, 2};
std::array<ssize_t, 3> shape{3, 4, 2};
NDArray<double, 3> a{shape};
NDArray<double, 3> b{shape};
for (uint32_t i = 0; i != a.size(); ++i) {
@ -179,9 +179,9 @@ TEST_CASE("Compare two images") {
}
TEST_CASE("Size and shape matches") {
int64_t w = 15;
int64_t h = 75;
std::array<int64_t, 2> shape{w, h};
ssize_t w = 15;
ssize_t h = 75;
std::array<ssize_t, 2> shape{w, h};
NDArray<double> a{shape};
REQUIRE(a.size() == w * h);
REQUIRE(a.shape() == shape);
@ -224,7 +224,7 @@ TEST_CASE("Bitwise and on data") {
TEST_CASE("Elementwise operations on images") {
std::array<int64_t, 2> shape{5, 5};
std::array<ssize_t, 2> shape{5, 5};
double a_val = 3.0;
double b_val = 8.0;

6
src/NDView.test.cpp
View File

@ -142,7 +142,7 @@ TEST_CASE("iterators") {
// for (int i = 0; i != 12; ++i) {
// vec.push_back(i);
// }
// std::vector<int64_t> shape{3, 4};
// std::vector<ssize_t> shape{3, 4};
// NDView<int, 2> data(vec.data(), shape);
// }
@ -151,8 +151,8 @@ TEST_CASE("divide with another span") {
std::vector<int> vec1{3, 2, 1};
std::vector<int> result{3, 6, 3};
NDView<int, 1> data0(vec0.data(), Shape<1>{static_cast<int64_t>(vec0.size())});
NDView<int, 1> data1(vec1.data(), Shape<1>{static_cast<int64_t>(vec1.size())});
NDView<int, 1> data0(vec0.data(), Shape<1>{static_cast<ssize_t>(vec0.size())});
NDView<int, 1> data1(vec1.data(), Shape<1>{static_cast<ssize_t>(vec1.size())});
data0 /= data1;

6
src/NumpyFile.cpp
View File

@ -72,8 +72,8 @@ void NumpyFile::get_frame_into(size_t frame_number, std::byte *image_buf) {
}
}
size_t NumpyFile::pixels_per_frame() { return m_pixels_per_frame; };
size_t NumpyFile::bytes_per_frame() { return m_bytes_per_frame; };
size_t NumpyFile::pixels_per_frame() { return m_pixels_per_frame; }
size_t NumpyFile::bytes_per_frame() { return m_bytes_per_frame; }
std::vector<Frame> NumpyFile::read_n(size_t n_frames) {
// TODO: implement this in a more efficient way
@ -197,4 +197,4 @@ void NumpyFile::load_metadata() {
m_header = {dtype, fortran_order, shape};
}
} // namespace aare
} // namespace aare

146
src/RawFile.cpp
View File

@ -1,6 +1,8 @@
#include "aare/RawFile.hpp"
#include "aare/algorithm.hpp"
#include "aare/PixelMap.hpp"
#include "aare/defs.hpp"
#include "aare/logger.hpp"
#include "aare/geo_helpers.hpp"
#include <fmt/format.h>
@ -14,27 +16,18 @@ RawFile::RawFile(const std::filesystem::path &fname, const std::string &mode)
: m_master(fname) {
m_mode = mode;
if (mode == "r") {
n_subfiles = find_number_of_subfiles(); // f0,f1...fn
n_subfile_parts =
m_master.geometry().col * m_master.geometry().row; // d0,d1...dn
find_geometry();
if (m_master.roi()){
m_geometry = update_geometry_with_roi(m_geometry, m_master.roi().value());
}
open_subfiles();
} else {
throw std::runtime_error(LOCATION +
"Unsupported mode. Can only read RawFiles.");
" Unsupported mode. Can only read RawFiles.");
}
}
Frame RawFile::read_frame() { return get_frame(m_current_frame++); };
Frame RawFile::read_frame() { return get_frame(m_current_frame++); }
Frame RawFile::read_frame(size_t frame_number) {
seek(frame_number);
@ -52,13 +45,13 @@ void RawFile::read_into(std::byte *image_buf, size_t n_frames) {
void RawFile::read_into(std::byte *image_buf) {
return get_frame_into(m_current_frame++, image_buf);
};
}
void RawFile::read_into(std::byte *image_buf, DetectorHeader *header) {
return get_frame_into(m_current_frame++, image_buf, header);
};
}
void RawFile::read_into(std::byte *image_buf, size_t n_frames, DetectorHeader *header) {
// return get_frame_into(m_current_frame++, image_buf, header);
@ -67,12 +60,12 @@ void RawFile::read_into(std::byte *image_buf, size_t n_frames, DetectorHeader *h
this->get_frame_into(m_current_frame++, image_buf, header);
image_buf += bytes_per_frame();
if(header)
header+=n_mod();
header+=n_modules();
}
};
}
size_t RawFile::n_mod() const { return n_subfile_parts; }
size_t RawFile::n_modules() const { return m_master.n_modules(); }
size_t RawFile::bytes_per_frame() {
@ -94,9 +87,9 @@ void RawFile::seek(size_t frame_index) {
frame_index, total_frames()));
}
m_current_frame = frame_index;
};
}
size_t RawFile::tell() { return m_current_frame; };
size_t RawFile::tell() { return m_current_frame; }
size_t RawFile::total_frames() const { return m_master.frames_in_file(); }
size_t RawFile::rows() const { return m_geometry.pixels_y; }
@ -106,17 +99,11 @@ xy RawFile::geometry() { return m_master.geometry(); }
void RawFile::open_subfiles() {
if (m_mode == "r")
for (size_t i = 0; i != n_subfiles; ++i) {
auto v = std::vector<RawSubFile *>(n_subfile_parts);
for (size_t j = 0; j != n_subfile_parts; ++j) {
auto pos = m_geometry.module_pixel_0[j];
v[j] = new RawSubFile(m_master.data_fname(j, i),
m_master.detector_type(), pos.height,
pos.width, m_master.bitdepth(),
pos.row_index, pos.col_index);
}
subfiles.push_back(v);
for (size_t i = 0; i != n_modules(); ++i) {
auto pos = m_geometry.module_pixel_0[i];
m_subfiles.emplace_back(std::make_unique<RawSubFile>(
m_master.data_fname(i, 0), m_master.detector_type(), pos.height,
pos.width, m_master.bitdepth(), pos.row_index, pos.col_index));
}
else {
throw std::runtime_error(LOCATION +
@ -141,18 +128,6 @@ DetectorHeader RawFile::read_header(const std::filesystem::path &fname) {
return h;
}
int RawFile::find_number_of_subfiles() {
int n_files = 0;
// f0,f1...fn How many files is the data split into?
while (std::filesystem::exists(m_master.data_fname(0, n_files)))
n_files++; // increment after test
#ifdef AARE_VERBOSE
fmt::print("Found: {} subfiles\n", n_files);
#endif
return n_files;
}
RawMasterFile RawFile::master() const { return m_master; }
@ -168,7 +143,7 @@ void RawFile::find_geometry() {
uint16_t c{};
for (size_t i = 0; i < n_subfile_parts; i++) {
for (size_t i = 0; i < n_modules(); i++) {
auto h = read_header(m_master.data_fname(i, 0));
r = std::max(r, h.row);
c = std::max(c, h.column);
@ -210,70 +185,58 @@ size_t RawFile::bytes_per_pixel() const {
}
void RawFile::get_frame_into(size_t frame_index, std::byte *frame_buffer, DetectorHeader *header) {
LOG(logDEBUG) << "RawFile::get_frame_into(" << frame_index << ")";
if (frame_index >= total_frames()) {
throw std::runtime_error(LOCATION + "Frame number out of range");
}
std::vector<size_t> frame_numbers(n_subfile_parts);
std::vector<size_t> frame_indices(n_subfile_parts, frame_index);
std::vector<size_t> frame_numbers(n_modules());
std::vector<size_t> frame_indices(n_modules(), frame_index);
// sync the frame numbers
if (n_subfile_parts != 1) {
for (size_t part_idx = 0; part_idx != n_subfile_parts; ++part_idx) {
auto subfile_id = frame_index / m_master.max_frames_per_file();
if (subfile_id >= subfiles.size()) {
throw std::runtime_error(LOCATION +
" Subfile out of range. Possible missing data.");
}
frame_numbers[part_idx] =
subfiles[subfile_id][part_idx]->frame_number(
frame_index % m_master.max_frames_per_file());
if (n_modules() != 1) { //if we have more than one module
for (size_t part_idx = 0; part_idx != n_modules(); ++part_idx) {
frame_numbers[part_idx] = m_subfiles[part_idx]->frame_number(frame_index);
}
// 1. if frame number vector is the same break
while (std::adjacent_find(frame_numbers.begin(), frame_numbers.end(),
std::not_equal_to<>()) !=
frame_numbers.end()) {
while (!all_equal(frame_numbers)) {
// 2. find the index of the minimum frame number,
auto min_frame_idx = std::distance(
frame_numbers.begin(),
std::min_element(frame_numbers.begin(), frame_numbers.end()));
// 3. increase its index and update its respective frame number
frame_indices[min_frame_idx]++;
// 4. if we can't increase its index => throw error
if (frame_indices[min_frame_idx] >= total_frames()) {
throw std::runtime_error(LOCATION +
"Frame number out of range");
}
auto subfile_id =
frame_indices[min_frame_idx] / m_master.max_frames_per_file();
frame_numbers[min_frame_idx] =
subfiles[subfile_id][min_frame_idx]->frame_number(
frame_indices[min_frame_idx] %
m_master.max_frames_per_file());
m_subfiles[min_frame_idx]->frame_number(frame_indices[min_frame_idx]);
}
}
if (m_master.geometry().col == 1) {
// get the part from each subfile and copy it to the frame
for (size_t part_idx = 0; part_idx != n_subfile_parts; ++part_idx) {
for (size_t part_idx = 0; part_idx != n_modules(); ++part_idx) {
auto corrected_idx = frame_indices[part_idx];
auto subfile_id = corrected_idx / m_master.max_frames_per_file();
if (subfile_id >= subfiles.size()) {
throw std::runtime_error(LOCATION +
" Subfile out of range. Possible missing data.");
}
// This is where we start writing
auto offset = (m_geometry.module_pixel_0[part_idx].origin_y * m_geometry.pixels_x +
m_geometry.module_pixel_0[part_idx].origin_x)*m_master.bitdepth()/8;
if (m_geometry.module_pixel_0[part_idx].origin_x!=0)
throw std::runtime_error(LOCATION + "Implementation error. x pos not 0.");
throw std::runtime_error(LOCATION + " Implementation error. x pos not 0.");
//TODO! Risk for out of range access
subfiles[subfile_id][part_idx]->seek(corrected_idx % m_master.max_frames_per_file());
subfiles[subfile_id][part_idx]->read_into(frame_buffer + offset, header);
//TODO! What if the files don't match?
m_subfiles[part_idx]->seek(corrected_idx);
m_subfiles[part_idx]->read_into(frame_buffer + offset, header);
if (header)
++header;
}
@ -282,26 +245,21 @@ void RawFile::get_frame_into(size_t frame_index, std::byte *frame_buffer, Detect
//TODO! should we read row by row?
// create a buffer large enough to hold a full module
auto bytes_per_part = m_master.pixels_y() * m_master.pixels_x() *
m_master.bitdepth() /
8; // TODO! replace with image_size_in_bytes
auto *part_buffer = new std::byte[bytes_per_part];
// TODO! if we have many submodules we should reorder them on the module
// level
for (size_t part_idx = 0; part_idx != n_subfile_parts; ++part_idx) {
for (size_t part_idx = 0; part_idx != n_modules(); ++part_idx) {
auto pos = m_geometry.module_pixel_0[part_idx];
auto corrected_idx = frame_indices[part_idx];
auto subfile_id = corrected_idx / m_master.max_frames_per_file();
if (subfile_id >= subfiles.size()) {
throw std::runtime_error(LOCATION +
" Subfile out of range. Possible missing data.");
}
subfiles[subfile_id][part_idx]->seek(corrected_idx % m_master.max_frames_per_file());
subfiles[subfile_id][part_idx]->read_into(part_buffer, header);
m_subfiles[part_idx]->seek(corrected_idx);
m_subfiles[part_idx]->read_into(part_buffer, header);
if(header)
++header;
@ -321,6 +279,7 @@ void RawFile::get_frame_into(size_t frame_index, std::byte *frame_buffer, Detect
}
delete[] part_buffer;
}
}
std::vector<Frame> RawFile::read_n(size_t n_frames) {
@ -337,27 +296,8 @@ size_t RawFile::frame_number(size_t frame_index) {
if (frame_index >= m_master.frames_in_file()) {
throw std::runtime_error(LOCATION + " Frame number out of range");
}
size_t subfile_id = frame_index / m_master.max_frames_per_file();
if (subfile_id >= subfiles.size()) {
throw std::runtime_error(
LOCATION + " Subfile out of range. Possible missing data.");
}
return subfiles[subfile_id][0]->frame_number(
frame_index % m_master.max_frames_per_file());
}
RawFile::~RawFile() {
// TODO! Fix this, for file closing
for (auto &vec : subfiles) {
for (auto *subfile : vec) {
delete subfile;
}
}
return m_subfiles[0]->frame_number(frame_index);
}
} // namespace aare
} // namespace aare

7
src/RawFile.test.cpp
View File

@ -99,11 +99,11 @@ TEST_CASE("Read frame numbers from a raw file", "[.integration]") {
}
}
TEST_CASE("Compare reading from a numpy file with a raw file", "[.integration]") {
auto fpath_raw = test_data_path() / "jungfrau" / "jungfrau_single_master_0.json";
TEST_CASE("Compare reading from a numpy file with a raw file", "[.files]") {
auto fpath_raw = test_data_path() / "raw/jungfrau" / "jungfrau_single_master_0.json";
REQUIRE(std::filesystem::exists(fpath_raw));
auto fpath_npy = test_data_path() / "jungfrau" / "jungfrau_single_0.npy";
auto fpath_npy = test_data_path() / "raw/jungfrau" / "jungfrau_single_0.npy";
REQUIRE(std::filesystem::exists(fpath_npy));
File raw(fpath_raw, "r");
@ -113,6 +113,7 @@ TEST_CASE("Compare reading from a numpy file with a raw file", "[.integration]")
CHECK(npy.total_frames() == 10);
for (size_t i = 0; i < 10; ++i) {
CHECK(raw.tell() == i);
auto raw_frame = raw.read_frame();
auto npy_frame = npy.read_frame();
CHECK((raw_frame.view<uint16_t>() == npy_frame.view<uint16_t>()));

8
src/RawMasterFile.cpp
View File

@ -87,7 +87,7 @@ int ScanParameters::start() const { return m_start; }
int ScanParameters::stop() const { return m_stop; }
void ScanParameters::increment_stop(){
m_stop += 1;
};
}
int ScanParameters::step() const { return m_step; }
const std::string &ScanParameters::dac() const { return m_dac; }
bool ScanParameters::enabled() const { return m_enabled; }
@ -140,6 +140,10 @@ std::optional<size_t> RawMasterFile::number_of_rows() const {
xy RawMasterFile::geometry() const { return m_geometry; }
size_t RawMasterFile::n_modules() const {
return m_geometry.row * m_geometry.col;
}
std::optional<uint8_t> RawMasterFile::quad() const { return m_quad; }
// optional values, these may or may not be present in the master file
@ -417,4 +421,4 @@ void RawMasterFile::parse_raw(const std::filesystem::path &fpath) {
if(m_frames_in_file==0)
m_frames_in_file = m_total_frames_expected;
}
} // namespace aare
} // namespace aare

128
src/RawSubFile.cpp
View File

@ -1,9 +1,15 @@
#include "aare/RawSubFile.hpp"
#include "aare/PixelMap.hpp"
#include "aare/algorithm.hpp"
#include "aare/utils/ifstream_helpers.hpp"
#include "aare/logger.hpp"
#include <cstring> // memcpy
#include <fmt/core.h>
#include <iostream>
#include <regex>
@ -12,51 +18,51 @@ namespace aare {
RawSubFile::RawSubFile(const std::filesystem::path &fname,
DetectorType detector, size_t rows, size_t cols,
size_t bitdepth, uint32_t pos_row, uint32_t pos_col)
: m_detector_type(detector), m_bitdepth(bitdepth), m_fname(fname),
: m_detector_type(detector), m_bitdepth(bitdepth),
m_rows(rows), m_cols(cols),
m_bytes_per_frame((m_bitdepth / 8) * m_rows * m_cols), m_pos_row(pos_row),
m_pos_col(pos_col) {
LOG(logDEBUG) << "RawSubFile::RawSubFile()";
if (m_detector_type == DetectorType::Moench03_old) {
m_pixel_map = GenerateMoench03PixelMap();
} else if (m_detector_type == DetectorType::Eiger && m_pos_row % 2 == 0) {
m_pixel_map = GenerateEigerFlipRowsPixelMap();
}
if (std::filesystem::exists(fname)) {
m_num_frames = std::filesystem::file_size(fname) /
(sizeof(DetectorHeader) + rows * cols * bitdepth / 8);
} else {
throw std::runtime_error(
LOCATION + fmt::format("File {} does not exist", m_fname.string()));
}
// fp = fopen(m_fname.string().c_str(), "rb");
m_file.open(m_fname, std::ios::binary);
if (!m_file.is_open()) {
throw std::runtime_error(
LOCATION + fmt::format("Could not open file {}", m_fname.string()));
}
#ifdef AARE_VERBOSE
fmt::print("Opened file: {} with {} frames\n", m_fname.string(), m_num_frames);
fmt::print("m_rows: {}, m_cols: {}, m_bitdepth: {}\n", m_rows, m_cols,
m_bitdepth);
fmt::print("file size: {}\n", std::filesystem::file_size(fname));
#endif
parse_fname(fname);
scan_files();
open_file(m_current_file_index); // open the first file
}
void RawSubFile::seek(size_t frame_index) {
if (frame_index >= m_num_frames) {
throw std::runtime_error(LOCATION + fmt::format("Frame index {} out of range in a file with {} frames", frame_index, m_num_frames));
LOG(logDEBUG) << "RawSubFile::seek(" << frame_index << ")";
if (frame_index >= m_total_frames) {
throw std::runtime_error(LOCATION + " Frame index out of range: " +
std::to_string(frame_index));
}
m_file.seekg((sizeof(DetectorHeader) + bytes_per_frame()) * frame_index);
m_current_frame_index = frame_index;
auto file_index = first_larger(m_last_frame_in_file, frame_index);
if (file_index != m_current_file_index)
open_file(file_index);
auto frame_offset = (file_index)
? frame_index - m_last_frame_in_file[file_index - 1]
: frame_index;
auto byte_offset = frame_offset * (m_bytes_per_frame + sizeof(DetectorHeader));
m_file.seekg(byte_offset);
}
size_t RawSubFile::tell() {
return m_file.tellg() / (sizeof(DetectorHeader) + bytes_per_frame());
LOG(logDEBUG) << "RawSubFile::tell():" << m_current_frame_index;
return m_current_frame_index;
}
void RawSubFile::read_into(std::byte *image_buf, DetectorHeader *header) {
LOG(logDEBUG) << "RawSubFile::read_into()";
if (header) {
m_file.read(reinterpret_cast<char *>(header), sizeof(DetectorHeader));
} else {
@ -90,6 +96,13 @@ void RawSubFile::read_into(std::byte *image_buf, DetectorHeader *header) {
if (m_file.fail()){
throw std::runtime_error(LOCATION + ifstream_error_msg(m_file));
}
++ m_current_frame_index;
if (m_current_frame_index >= m_last_frame_in_file[m_current_file_index] &&
(m_current_frame_index < m_total_frames)) {
++m_current_file_index;
open_file(m_current_file_index);
}
}
void RawSubFile::read_into(std::byte *image_buf, size_t n_frames, DetectorHeader *header) {
@ -130,4 +143,69 @@ size_t RawSubFile::frame_number(size_t frame_index) {
return h.frameNumber;
}
void RawSubFile::parse_fname(const std::filesystem::path &fname) {
LOG(logDEBUG) << "RawSubFile::parse_fname()";
// data has the format: /path/too/data/jungfrau_single_d0_f1_0.raw
// d0 is the module index, will not change for this file
// f1 is the file index - thi is the one we need
// 0 is the measurement index, will not change
m_path = fname.parent_path();
m_base_name = fname.filename();
// Regex to extract numbers after 'd' and 'f'
std::regex pattern(R"(^(.*_d)(\d+)(_f)(\d+)(_\d+\.raw)$)");
std::smatch match;
if (std::regex_match(m_base_name, match, pattern)) {
m_offset = std::stoi(match[4].str()); // find the first file index in case of a truncated series
m_base_name = match[1].str() + match[2].str() + match[3].str() + "{}" + match[5].str();
LOG(logDEBUG) << "Base name: " << m_base_name;
LOG(logDEBUG) << "Offset: " << m_offset;
LOG(logDEBUG) << "Path: " << m_path.string();
} else {
throw std::runtime_error(
LOCATION + fmt::format("Could not parse file name {}", fname.string()));
}
}
std::filesystem::path RawSubFile::fpath(size_t file_index) const {
auto fname = fmt::format(m_base_name, file_index);
return m_path / fname;
}
void RawSubFile::open_file(size_t file_index) {
m_file.close();
auto fname = fpath(file_index+m_offset);
LOG(logDEBUG) << "RawSubFile::open_file(): " << fname.string();
m_file.open(fname, std::ios::binary);
if (!m_file.is_open()) {
throw std::runtime_error(
LOCATION + fmt::format("Could not open file {}", fpath(file_index).string()));
}
m_current_file_index = file_index;
}
void RawSubFile::scan_files() {
LOG(logDEBUG) << "RawSubFile::scan_files()";
// find how many files we have and the number of frames in each file
m_last_frame_in_file.clear();
size_t file_index = m_offset;
while (std::filesystem::exists(fpath(file_index))) {
auto n_frames = std::filesystem::file_size(fpath(file_index)) /
(m_bytes_per_frame + sizeof(DetectorHeader));
m_last_frame_in_file.push_back(n_frames);
LOG(logDEBUG) << "Found: " << n_frames << " frames in file: " << fpath(file_index).string();
++file_index;
}
// find where we need to open the next file and total number of frames
m_last_frame_in_file = cumsum(m_last_frame_in_file);
if(m_last_frame_in_file.empty()){
m_total_frames = 0;
}else{
m_total_frames = m_last_frame_in_file.back();
}
}
} // namespace aare

76
src/RawSubFile.test.cpp Normal file
View File

@ -0,0 +1,76 @@
#include "aare/RawSubFile.hpp"
#include "aare/File.hpp"
#include "aare/NDArray.hpp"
#include <catch2/catch_test_macros.hpp>
#include "test_config.hpp"
using namespace aare;
TEST_CASE("Read frames directly from a RawSubFile", "[.files]"){
auto fpath_raw = test_data_path() / "raw/jungfrau" / "jungfrau_single_d0_f0_0.raw";
REQUIRE(std::filesystem::exists(fpath_raw));
RawSubFile f(fpath_raw, DetectorType::Jungfrau, 512, 1024, 16);
REQUIRE(f.rows() == 512);
REQUIRE(f.cols() == 1024);
REQUIRE(f.pixels_per_frame() == 512 * 1024);
REQUIRE(f.bytes_per_frame() == 512 * 1024 * 2);
REQUIRE(f.bytes_per_pixel() == 2);
auto fpath_npy = test_data_path() / "raw/jungfrau" / "jungfrau_single_0.npy";
REQUIRE(std::filesystem::exists(fpath_npy));
//Numpy file with the same data to use as reference
File npy(fpath_npy, "r");
CHECK(f.frames_in_file() == 10);
CHECK(npy.total_frames() == 10);
DetectorHeader header{};
NDArray<uint16_t, 2> image({static_cast<ssize_t>(f.rows()), static_cast<ssize_t>(f.cols())});
for (size_t i = 0; i < 10; ++i) {
CHECK(f.tell() == i);
f.read_into(image.buffer(), &header);
auto npy_frame = npy.read_frame();
CHECK((image.view() == npy_frame.view<uint16_t>()));
}
}
TEST_CASE("Read frames directly from a RawSubFile starting at the second file", "[.files]"){
// we know this file has 10 frames with frame numbers 1 to 10
// f0 1,2,3
// f1 4,5,6 <-- starting here
// f2 7,8,9
// f3 10
auto fpath_raw = test_data_path() / "raw/jungfrau" / "jungfrau_single_d0_f1_0.raw";
REQUIRE(std::filesystem::exists(fpath_raw));
RawSubFile f(fpath_raw, DetectorType::Jungfrau, 512, 1024, 16);
auto fpath_npy = test_data_path() / "raw/jungfrau" / "jungfrau_single_0.npy";
REQUIRE(std::filesystem::exists(fpath_npy));
//Numpy file with the same data to use as reference
File npy(fpath_npy, "r");
npy.seek(3);
CHECK(f.frames_in_file() == 7);
CHECK(npy.total_frames() == 10);
DetectorHeader header{};
NDArray<uint16_t, 2> image({static_cast<ssize_t>(f.rows()), static_cast<ssize_t>(f.cols())});
for (size_t i = 0; i < 7; ++i) {
CHECK(f.tell() == i);
f.read_into(image.buffer(), &header);
// frame numbers start at 1 frame index at 0
// adding 3 + 1 to verify the frame number
CHECK(header.frameNumber == i + 4);
auto npy_frame = npy.read_frame();
CHECK((image.view() == npy_frame.view<uint16_t>()));
}
}

88
src/algorithm.test.cpp
View File

@ -1,8 +1,7 @@
#include <catch2/catch_test_macros.hpp>
#include <aare/algorithm.hpp>
#include <catch2/catch_test_macros.hpp>
TEST_CASE("Find the closed index in a 1D array", "[algorithm]") {
aare::NDArray<double, 1> arr({5});
@ -17,7 +16,7 @@ TEST_CASE("Find the closed index in a 1D array", "[algorithm]") {
REQUIRE(aare::nearest_index(arr, -1.0) == 0);
}
TEST_CASE("Passing integers to nearest_index works", "[algorithm]"){
TEST_CASE("Passing integers to nearest_index works", "[algorithm]") {
aare::NDArray<int, 1> arr({5});
for (ssize_t i = 0; i < arr.size(); i++) {
arr[i] = i;
@ -30,8 +29,7 @@ TEST_CASE("Passing integers to nearest_index works", "[algorithm]"){
REQUIRE(aare::nearest_index(arr, -1) == 0);
}
TEST_CASE("nearest_index works with std::vector", "[algorithm]"){
TEST_CASE("nearest_index works with std::vector", "[algorithm]") {
std::vector<double> vec = {0, 1, 2, 3, 4};
REQUIRE(aare::nearest_index(vec, 2.123) == 2);
REQUIRE(aare::nearest_index(vec, 2.66) == 3);
@ -40,7 +38,7 @@ TEST_CASE("nearest_index works with std::vector", "[algorithm]"){
REQUIRE(aare::nearest_index(vec, -10.0) == 0);
}
TEST_CASE("nearest index works with std::array", "[algorithm]"){
TEST_CASE("nearest index works with std::array", "[algorithm]") {
std::array<double, 5> arr = {0, 1, 2, 3, 4};
REQUIRE(aare::nearest_index(arr, 2.123) == 2);
REQUIRE(aare::nearest_index(arr, 2.501) == 3);
@ -49,18 +47,20 @@ TEST_CASE("nearest index works with std::array", "[algorithm]"){
REQUIRE(aare::nearest_index(arr, -10.0) == 0);
}
TEST_CASE("nearest index when there is no different uses the first element", "[algorithm]"){
TEST_CASE("nearest index when there is no different uses the first element",
"[algorithm]") {
std::vector<int> vec = {5, 5, 5, 5, 5};
REQUIRE(aare::nearest_index(vec, 5) == 0);
}
TEST_CASE("nearest index when there is no different uses the first element also when all smaller", "[algorithm]"){
TEST_CASE("nearest index when there is no different uses the first element "
"also when all smaller",
"[algorithm]") {
std::vector<int> vec = {5, 5, 5, 5, 5};
REQUIRE(aare::nearest_index(vec, 10) == 0);
}
TEST_CASE("last smaller", "[algorithm]"){
TEST_CASE("last smaller", "[algorithm]") {
aare::NDArray<double, 1> arr({5});
for (ssize_t i = 0; i < arr.size(); i++) {
arr[i] = i;
@ -72,17 +72,17 @@ TEST_CASE("last smaller", "[algorithm]"){
REQUIRE(aare::last_smaller(arr, 253.) == 4);
}
TEST_CASE("returns last bin strictly smaller", "[algorithm]"){
TEST_CASE("returns last bin strictly smaller", "[algorithm]") {
aare::NDArray<double, 1> arr({5});
for (ssize_t i = 0; i < arr.size(); i++) {
arr[i] = i;
}
// arr 0, 1, 2, 3, 4
REQUIRE(aare::last_smaller(arr, 2.0) == 1);
}
TEST_CASE("last_smaller with all elements smaller returns last element", "[algorithm]"){
TEST_CASE("last_smaller with all elements smaller returns last element",
"[algorithm]") {
aare::NDArray<double, 1> arr({5});
for (ssize_t i = 0; i < arr.size(); i++) {
arr[i] = i;
@ -91,7 +91,8 @@ TEST_CASE("last_smaller with all elements smaller returns last element", "[algor
REQUIRE(aare::last_smaller(arr, 50.) == 4);
}
TEST_CASE("last_smaller with all elements bigger returns first element", "[algorithm]"){
TEST_CASE("last_smaller with all elements bigger returns first element",
"[algorithm]") {
aare::NDArray<double, 1> arr({5});
for (ssize_t i = 0; i < arr.size(); i++) {
arr[i] = i;
@ -100,38 +101,41 @@ TEST_CASE("last_smaller with all elements bigger returns first element", "[algor
REQUIRE(aare::last_smaller(arr, -50.) == 0);
}
TEST_CASE("last smaller with all elements equal returns the first element", "[algorithm]"){
std::vector<int> vec = {5,5,5,5,5,5,5};
TEST_CASE("last smaller with all elements equal returns the first element",
"[algorithm]") {
std::vector<int> vec = {5, 5, 5, 5, 5, 5, 5};
REQUIRE(aare::last_smaller(vec, 5) == 0);
}
TEST_CASE("first_lager with vector", "[algorithm]"){
TEST_CASE("first_lager with vector", "[algorithm]") {
std::vector<double> vec = {0, 1, 2, 3, 4};
REQUIRE(aare::first_larger(vec, 2.5) == 3);
}
TEST_CASE("first_lager with all elements smaller returns last element", "[algorithm]"){
TEST_CASE("first_lager with all elements smaller returns last element",
"[algorithm]") {
std::vector<double> vec = {0, 1, 2, 3, 4};
REQUIRE(aare::first_larger(vec, 50.) == 4);
}
TEST_CASE("first_lager with all elements bigger returns first element", "[algorithm]"){
TEST_CASE("first_lager with all elements bigger returns first element",
"[algorithm]") {
std::vector<double> vec = {0, 1, 2, 3, 4};
REQUIRE(aare::first_larger(vec, -50.) == 0);
}
TEST_CASE("first_lager with all elements the same as the check returns last", "[algorithm]"){
TEST_CASE("first_lager with all elements the same as the check returns last",
"[algorithm]") {
std::vector<int> vec = {14, 14, 14, 14, 14};
REQUIRE(aare::first_larger(vec, 14) == 4);
}
TEST_CASE("first larger with the same element", "[algorithm]"){
std::vector<int> vec = {7,8,9,10,11};
TEST_CASE("first larger with the same element", "[algorithm]") {
std::vector<int> vec = {7, 8, 9, 10, 11};
REQUIRE(aare::first_larger(vec, 9) == 3);
}
TEST_CASE("cumsum works", "[algorithm]"){
TEST_CASE("cumsum works", "[algorithm]") {
std::vector<double> vec = {0, 1, 2, 3, 4};
auto result = aare::cumsum(vec);
REQUIRE(result.size() == vec.size());
@ -141,12 +145,12 @@ TEST_CASE("cumsum works", "[algorithm]"){
REQUIRE(result[3] == 6);
REQUIRE(result[4] == 10);
}
TEST_CASE("cumsum works with empty vector", "[algorithm]"){
TEST_CASE("cumsum works with empty vector", "[algorithm]") {
std::vector<double> vec = {};
auto result = aare::cumsum(vec);
REQUIRE(result.size() == 0);
}
TEST_CASE("cumsum works with negative numbers", "[algorithm]"){
TEST_CASE("cumsum works with negative numbers", "[algorithm]") {
std::vector<double> vec = {0, -1, -2, -3, -4};
auto result = aare::cumsum(vec);
REQUIRE(result.size() == vec.size());
@ -157,3 +161,35 @@ TEST_CASE("cumsum works with negative numbers", "[algorithm]"){
REQUIRE(result[4] == -10);
}
TEST_CASE("cumsum on an empty vector", "[algorithm]") {
std::vector<double> vec = {};
auto result = aare::cumsum(vec);
REQUIRE(result.size() == 0);
}
TEST_CASE("All equal on an empty vector is false", "[algorithm]") {
std::vector<int> vec = {};
REQUIRE(aare::all_equal(vec) == false);
}
TEST_CASE("All equal on a vector with 1 element is true", "[algorithm]") {
std::vector<int> vec = {1};
REQUIRE(aare::all_equal(vec) == true);
}
TEST_CASE("All equal on a vector with 2 elements is true", "[algorithm]") {
std::vector<int> vec = {1, 1};
REQUIRE(aare::all_equal(vec) == true);
}
TEST_CASE("All equal on a vector with two different elements is false", "[algorithm]") {
std::vector<int> vec = {1, 2};
REQUIRE(aare::all_equal(vec) == false);
}
TEST_CASE("Last element is different", "[algorithm]") {
std::vector<int> vec = {1, 1, 1, 1, 2};
REQUIRE(aare::all_equal(vec) == false);
}

8
src/decode.cpp
View File

@ -26,8 +26,8 @@ void adc_sar_05_decode64to16(NDView<uint64_t, 2> input, NDView<uint16_t,2> outpu
throw std::invalid_argument(LOCATION + " input and output shapes must match");
}
for(int64_t i = 0; i < input.shape(0); i++){
for(int64_t j = 0; j < input.shape(1); j++){
for(ssize_t i = 0; i < input.shape(0); i++){
for(ssize_t j = 0; j < input.shape(1); j++){
output(i,j) = adc_sar_05_decode64to16(input(i,j));
}
}
@ -56,8 +56,8 @@ void adc_sar_04_decode64to16(NDView<uint64_t, 2> input, NDView<uint16_t,2> outpu
if(input.shape() != output.shape()){
throw std::invalid_argument(LOCATION + " input and output shapes must match");
}
for(int64_t i = 0; i < input.shape(0); i++){
for(int64_t j = 0; j < input.shape(1); j++){
for(ssize_t i = 0; i < input.shape(0); i++){
for(ssize_t j = 0; j < input.shape(1); j++){
output(i,j) = adc_sar_04_decode64to16(input(i,j));
}
}

57
update_version.py Normal file
View File

@ -0,0 +1,57 @@
# SPDX-License-Identifier: LGPL-3.0-or-other
# Copyright (C) 2021 Contributors to the Aare Package
"""
Script to update VERSION file with semantic versioning if provided as an argument, or with 0.0.0 if no argument is provided.
"""
import sys
import os
import re
from packaging.version import Version, InvalidVersion
SCRIPT_DIR = os.path.dirname(os.path.abspath(__file__))
def is_integer(value):
try:
int(value)
except ValueError:
return False
else:
return True
def get_version():
# Check at least one argument is passed
if len(sys.argv) < 2:
return "0.0.0"
version = sys.argv[1]
try:
v = Version(version) # normalize check if version follows PEP 440 specification
version_normalized = version.replace("-", ".")
version_normalized = re.sub(r'0*(\d+)', lambda m : str(int(m.group(0))), version_normalized) #remove leading zeros
return version_normalized
except InvalidVersion as e:
print(f"Invalid version {version}. Version format must follow semantic versioning format of python PEP 440 version identification specification.")
sys.exit(1)
def write_version_to_file(version):
version_file_path = os.path.join(SCRIPT_DIR, "VERSION")
with open(version_file_path, "w") as version_file:
version_file.write(version)
print(f"Version {version} written to VERSION file.")
# Main script
if __name__ == "__main__":
version = get_version()
write_version_to_file(version)