Cluster reading, expression templates (#99)

Co-authored-by: froejdh_e <erik.frojdh@psi.ch>

CMakeLists.txt

@@ -29,6 +29,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_TESTS "Build tests" OFF)
+option(AARE_BENCHMARKS "Build benchmarks" OFF)
 option(AARE_EXAMPLES "Build examples" OFF)
 option(AARE_IN_GITHUB_ACTIONS "Running in Github Actions" OFF)
 option(AARE_DOCS "Build documentation" OFF)
@@ -63,6 +64,10 @@ if(AARE_CUSTOM_ASSERT)
     add_compile_definitions(AARE_CUSTOM_ASSERT)
 endif()
 
+if(AARE_BENCHMARKS)
+    add_subdirectory(benchmarks)
+endif()
+
 
 set(CMAKE_EXPORT_COMPILE_COMMANDS ON)
 
@@ -242,7 +247,9 @@ endif()
 ###------------------------------------------------------------------------------------------
 
 set(PUBLICHEADERS
+    include/aare/ArrayExpr.hpp
     include/aare/ClusterFinder.hpp
+    include/aare/ClusterFile.hpp
     include/aare/CtbRawFile.hpp
     include/aare/defs.hpp
     include/aare/Dtype.hpp
@@ -265,6 +272,7 @@ set(PUBLICHEADERS
 
 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/Frame.cpp

benchmarks/CMakeLists.txt

@@ -0,0 +1,11 @@
+find_package(benchmark REQUIRED)
+
+add_executable(ndarray_benchmark ndarray_benchmark.cpp)
+
+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
+)

benchmarks/ndarray_benchmark.cpp

@@ -0,0 +1,136 @@
+#include <benchmark/benchmark.h>
+#include "aare/NDArray.hpp"
+
+
+using aare::NDArray;
+
+constexpr ssize_t size = 1024;
+class TwoArrays : public benchmark::Fixture {
+public:
+    NDArray<int,2> a{{size,size},0};
+    NDArray<int,2> b{{size,size},0};
+  void SetUp(::benchmark::State& state) {
+    for(uint32_t i = 0; i < size; i++){
+      for(uint32_t j = 0; j < size; j++){
+        a(i, j)= i*j+1;
+        b(i, j)= i*j+1;
+      }
+    }
+  }
+
+  // void TearDown(::benchmark::State& state) {
+  // }
+};
+
+
+
+
+BENCHMARK_F(TwoArrays, AddWithOperator)(benchmark::State& st) {
+   for (auto _ : st) {
+    // This code gets timed
+    NDArray<int,2> res = a+b;
+    benchmark::DoNotOptimize(res);
+  }
+}
+BENCHMARK_F(TwoArrays, AddWithIndex)(benchmark::State& st) {
+   for (auto _ : st) {
+    // This code gets timed
+    NDArray<int,2> res(a.shape());
+    for (uint32_t i = 0; i < a.size(); i++) {
+        res(i) = a(i) + b(i);
+    }
+    benchmark::DoNotOptimize(res);
+  }
+}
+
+BENCHMARK_F(TwoArrays, SubtractWithOperator)(benchmark::State& st) {
+   for (auto _ : st) {
+    // This code gets timed
+    NDArray<int,2> res = a-b;
+    benchmark::DoNotOptimize(res);
+  }
+}
+BENCHMARK_F(TwoArrays, SubtractWithIndex)(benchmark::State& st) {
+   for (auto _ : st) {
+    // This code gets timed
+    NDArray<int,2> res(a.shape());
+    for (uint32_t i = 0; i < a.size(); i++) {
+        res(i) = a(i) - b(i);
+    }
+    benchmark::DoNotOptimize(res);
+  }
+}
+
+BENCHMARK_F(TwoArrays, MultiplyWithOperator)(benchmark::State& st) {
+   for (auto _ : st) {
+    // This code gets timed
+    NDArray<int,2> res = a*b;
+    benchmark::DoNotOptimize(res);
+  }
+}
+BENCHMARK_F(TwoArrays, MultiplyWithIndex)(benchmark::State& st) {
+   for (auto _ : st) {
+    // This code gets timed
+    NDArray<int,2> res(a.shape());
+    for (uint32_t i = 0; i < a.size(); i++) {
+        res(i) = a(i) * b(i);
+    }
+    benchmark::DoNotOptimize(res);
+  }
+}
+
+BENCHMARK_F(TwoArrays, DivideWithOperator)(benchmark::State& st) {
+   for (auto _ : st) {
+    // This code gets timed
+    NDArray<int,2> res = a/b;
+    benchmark::DoNotOptimize(res);
+  }
+}
+BENCHMARK_F(TwoArrays, DivideWithIndex)(benchmark::State& st) {
+   for (auto _ : st) {
+    // This code gets timed
+    NDArray<int,2> res(a.shape());
+    for (uint32_t i = 0; i < a.size(); i++) {
+        res(i) = a(i) / b(i);
+    }
+    benchmark::DoNotOptimize(res);
+  }
+}
+
+BENCHMARK_F(TwoArrays, FourAddWithOperator)(benchmark::State& st) {
+   for (auto _ : st) {
+    // This code gets timed
+    NDArray<int,2> res = a+b+a+b;
+    benchmark::DoNotOptimize(res);
+  }
+}
+BENCHMARK_F(TwoArrays, FourAddWithIndex)(benchmark::State& st) {
+   for (auto _ : st) {
+    // This code gets timed
+    NDArray<int,2> res(a.shape());
+    for (uint32_t i = 0; i < a.size(); i++) {
+        res(i) = a(i) + b(i) + a(i) + b(i);
+    }
+    benchmark::DoNotOptimize(res);
+  }
+}
+
+BENCHMARK_F(TwoArrays, MultiplyAddDivideWithOperator)(benchmark::State& st) {
+   for (auto _ : st) {
+    // This code gets timed
+    NDArray<int,2> res = a*a+b/a;
+    benchmark::DoNotOptimize(res);
+  }
+}
+BENCHMARK_F(TwoArrays, MultiplyAddDivideWithIndex)(benchmark::State& st) {
+   for (auto _ : st) {
+    // This code gets timed
+    NDArray<int,2> res(a.shape());
+    for (uint32_t i = 0; i < a.size(); i++) {
+        res(i) = a(i) * a(i) + b(i) / a(i);
+    }
+    benchmark::DoNotOptimize(res);
+  }
+}
+
+BENCHMARK_MAIN();

conda-recipe/meta.yaml

@@ -1,6 +1,6 @@
 package:
   name: aare
-  version: 2024.11.11.dev0 #TODO! how to not duplicate this?
+  version: 2024.11.15.dev0 #TODO! how to not duplicate this?
 
 
 source:

docs/CMakeLists.txt

@@ -10,31 +10,36 @@ configure_file(${DOXYGEN_IN} ${DOXYGEN_OUT} @ONLY)
 set(SPHINX_SOURCE ${CMAKE_CURRENT_SOURCE_DIR}/src)
 set(SPHINX_BUILD ${CMAKE_CURRENT_BINARY_DIR})
 
-set(SPHINX_SOURCE_FILES
+
-  src/index.rst
+file(GLOB SPHINX_SOURCE_FILES CONFIGURE_DEPENDS "src/*.rst")
-  src/Installation.rst
+# set(SPHINX_SOURCE_FILES
-  src/Requirements.rst
+#   src/index.rst
-  src/NDArray.rst
+#   src/Installation.rst
-  src/NDView.rst
+#   src/Requirements.rst
-  src/File.rst
+#   src/NDArray.rst
-  src/Frame.rst
+#   src/NDView.rst
-  src/Dtype.rst
+#   src/File.rst
-  src/ClusterFinder.rst
+#   src/Frame.rst
-  src/Pedestal.rst
+#   src/Dtype.rst
-  src/RawFile.rst
+#   src/ClusterFinder.rst
-  src/RawSubFile.rst
+#   src/ClusterFile.rst
-  src/RawMasterFile.rst
+#   src/Pedestal.rst
-  src/VarClusterFinder.rst
+#   src/RawFile.rst
-  src/pyVarClusterFinder.rst
+#   src/RawSubFile.rst
-  src/pyFile.rst
+#   src/RawMasterFile.rst
-  src/pyCtbRawFile.rst
+#   src/VarClusterFinder.rst
-  src/pyRawFile.rst
+#   src/pyVarClusterFinder.rst
-  src/pyRawMasterFile.rst
+#   src/pyFile.rst
-)
+#   src/pyCtbRawFile.rst
+#   src/pyRawFile.rst
+#   src/pyRawMasterFile.rst
+# )
+
 
 foreach(filename ${SPHINX_SOURCE_FILES})
-configure_file(${CMAKE_CURRENT_SOURCE_DIR}/${filename} 
+    get_filename_component(fname ${filename} NAME)
-              "${SPHINX_BUILD}/${filename}")
+    message(STATUS "Copying ${filename} to ${SPHINX_BUILD}/src/${fname}")
+    configure_file(${filename} "${SPHINX_BUILD}/src/${fname}")  
 endforeach(filename ${SPHINX_SOURCE_FILES})
 
 configure_file(

docs/src/ClusterFile.rst

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

docs/src/index.rst

@@ -14,27 +14,19 @@ AARE
     Requirements
 
 
+
 .. toctree::
     :caption: Python API
     :maxdepth: 1
     
     pyFile
     pyCtbRawFile
+    pyClusterFile
     pyRawFile
     pyRawMasterFile
     pyVarClusterFinder
 
 
-.. toctree::
-    :caption: Python API
-    :maxdepth: 1
-    
-    pyFile
-    pyCtbRawFile
-    pyRawMasterFile
-    pyVarClusterFinder
-
-
 .. toctree::
     :caption: C++ API
     :maxdepth: 1
@@ -45,6 +37,7 @@ AARE
     File
     Dtype
     ClusterFinder
+    ClusterFile
     Pedestal
     RawFile
     RawSubFile

docs/src/pyClusterFile.rst

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

include/aare/ArrayExpr.hpp

@@ -0,0 +1,99 @@
+#pragma once
+#include <cstdint> //int64_t
+#include <cstddef> //size_t
+#include <array> 
+
+#include <cassert>
+namespace aare {
+
+template <typename E, int64_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(); }
+};
+
+template <typename A, typename B, int64_t Ndim>
+class ArrayAdd : public ArrayExpr<ArrayAdd<A, B, Ndim>, Ndim> {
+    const A &arr1_;
+    const B &arr2_;
+
+  public:
+    ArrayAdd(const A &arr1, const B &arr2) : arr1_(arr1), arr2_(arr2) {
+        assert(arr1.size() == arr2.size());
+    }
+    auto operator[](int i) const { return arr1_[i] + arr2_[i]; }
+    size_t size() const { return arr1_.size(); }
+    std::array<int64_t, Ndim> shape() const { return arr1_.shape(); }
+};
+
+template <typename A, typename B, int64_t Ndim>
+class ArraySub : public ArrayExpr<ArraySub<A, B, Ndim>, Ndim> {
+    const A &arr1_;
+    const B &arr2_;
+
+  public:
+    ArraySub(const A &arr1, const B &arr2) : arr1_(arr1), arr2_(arr2) {
+        assert(arr1.size() == arr2.size());
+    }
+    auto operator[](int i) const { return arr1_[i] - arr2_[i]; }
+    size_t size() const { return arr1_.size(); }
+    std::array<int64_t, Ndim> shape() const { return arr1_.shape(); }
+};
+
+template <typename A, typename B, int64_t Ndim>
+class ArrayMul : public ArrayExpr<ArrayMul<A, B, Ndim>,Ndim> {
+    const A &arr1_;
+    const B &arr2_;
+
+  public:
+    ArrayMul(const A &arr1, const B &arr2) : arr1_(arr1), arr2_(arr2) {
+        assert(arr1.size() == arr2.size());
+    }
+    auto operator[](int i) const { return arr1_[i] * arr2_[i]; }
+    size_t size() const { return arr1_.size(); }
+    std::array<int64_t, Ndim> shape() const { return arr1_.shape(); }
+};
+
+template <typename A, typename B, int64_t Ndim>
+class ArrayDiv : public ArrayExpr<ArrayDiv<A, B, Ndim>, Ndim> {
+    const A &arr1_;
+    const B &arr2_;
+
+  public:
+    ArrayDiv(const A &arr1, const B &arr2) : arr1_(arr1), arr2_(arr2) {
+        assert(arr1.size() == arr2.size());
+    }
+    auto operator[](int i) const { return arr1_[i] / arr2_[i]; }
+    size_t size() const { return arr1_.size(); }
+    std::array<int64_t, Ndim> shape() const { return arr1_.shape(); }
+};
+
+
+
+template <typename A, typename B, int64_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>
+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>
+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>
+auto operator/(const ArrayExpr<A, Ndim> &arr1, const ArrayExpr<B, Ndim> &arr2) {
+    return ArrayDiv<ArrayExpr<A, Ndim>, ArrayExpr<B, Ndim>, Ndim>(arr1, arr2);
+}
+
+
+
+} // namespace aare

include/aare/ClusterFile.hpp

@@ -0,0 +1,65 @@
+
+
+#include "aare/defs.hpp"
+#include <filesystem>
+#include <fstream>
+
+namespace aare {
+
+struct Cluster {
+    int16_t x;
+    int16_t y;
+    int32_t data[9];
+};
+
+typedef enum {
+    cBottomLeft = 0,
+    cBottomRight = 1,
+    cTopLeft = 2,
+    cTopRight = 3
+} corner;
+
+typedef enum {
+    pBottomLeft = 0,
+    pBottom = 1,
+    pBottomRight = 2,
+    pLeft = 3,
+    pCenter = 4,
+    pRight = 5,
+    pTopLeft = 6,
+    pTop = 7,
+    pTopRight = 8
+} pixel;
+
+struct ClusterAnalysis {
+    uint32_t c;
+    int32_t tot;
+    double etax;
+    double etay;
+};
+
+
+
+class ClusterFile {
+    FILE *fp{};
+    uint32_t m_num_left{};
+    size_t m_chunk_size{};
+
+  public:
+    ClusterFile(const std::filesystem::path &fname, size_t chunk_size = 1000);
+    std::vector<Cluster> read_clusters(size_t n_clusters);
+    std::vector<Cluster> read_frame(int32_t &out_fnum);
+    std::vector<Cluster>
+    read_cluster_with_cut(size_t n_clusters, double *noise_map, int nx, int ny);
+
+    int analyze_data(int32_t *data, int32_t *t2, int32_t *t3, char *quad,
+                 double *eta2x, double *eta2y, double *eta3x, double *eta3y);
+    int analyze_cluster(Cluster cl, int32_t *t2, int32_t *t3, char *quad,
+                    double *eta2x, double *eta2y, double *eta3x,
+                    double *eta3y);
+
+    size_t chunk_size() const { return m_chunk_size; }
+
+};
+
+} // namespace aare

include/aare/ClusterFinder.hpp

@@ -53,7 +53,7 @@ class ClusterFinder {
         return m_pedestal.mean();
     }
 
-    std::vector<Cluster>
+    std::vector<DynamicCluster>
     find_clusters_without_threshold(NDView<FRAME_TYPE, 2> frame,
                                     // Pedestal<PEDESTAL_TYPE> &pedestal,
                                     bool late_update = false) {
@@ -64,7 +64,7 @@ class ClusterFinder {
         };
         std::vector<pedestal_update> pedestal_updates;
 
-        std::vector<Cluster> clusters;
+        std::vector<DynamicCluster> clusters;
         std::vector<std::vector<eventType>> eventMask;
         for (int i = 0; i < frame.shape(0); i++) {
             eventMask.push_back(std::vector<eventType>(frame.shape(1)));
@@ -115,7 +115,7 @@ class ClusterFinder {
                 if (eventMask[iy][ix] == PHOTON &&
                     (frame(iy, ix) - m_pedestal.mean(iy, ix)) >= max) {
                     eventMask[iy][ix] = PHOTON_MAX;
-                    Cluster cluster(m_cluster_sizeX, m_cluster_sizeY,
+                    DynamicCluster cluster(m_cluster_sizeX, m_cluster_sizeY,
                                     Dtype(typeid(PEDESTAL_TYPE)));
                     cluster.x = ix;
                     cluster.y = iy;
@@ -149,11 +149,11 @@ class ClusterFinder {
     }
 
     // template <typename FRAME_TYPE, typename PEDESTAL_TYPE>
-    std::vector<Cluster>
+    std::vector<DynamicCluster>
     find_clusters_with_threshold(NDView<FRAME_TYPE, 2> frame,
                                  Pedestal<PEDESTAL_TYPE> &pedestal) {
         assert(m_threshold > 0);
-        std::vector<Cluster> clusters;
+        std::vector<DynamicCluster> clusters;
         std::vector<std::vector<eventType>> eventMask;
         for (int i = 0; i < frame.shape(0); i++) {
             eventMask.push_back(std::vector<eventType>(frame.shape(1)));
@@ -230,7 +230,7 @@ class ClusterFinder {
                 if (eventMask[iy][ix] == PHOTON &&
                     frame(iy, ix) - pedestal.mean(iy, ix) >= max) {
                     eventMask[iy][ix] = PHOTON_MAX;
-                    Cluster cluster(m_cluster_sizeX, m_cluster_sizeY,
+                    DynamicCluster cluster(m_cluster_sizeX, m_cluster_sizeY,
                                     Dtype(typeid(FRAME_TYPE)));
                     cluster.x = ix;
                     cluster.y = iy;

include/aare/NDArray.hpp

@@ -7,6 +7,7 @@ memory.
 TODO! Add expression templates for operators
 
 */
+#include "aare/ArrayExpr.hpp"
 #include "aare/NDView.hpp"
 
 #include <algorithm>
@@ -20,7 +21,9 @@ TODO! Add expression templates for operators
 
 namespace aare {
 
-template <typename T, int64_t Ndim = 2> class NDArray {
+
+template <typename T, int64_t Ndim = 2>
+class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
     std::array<int64_t, Ndim> shape_;
     std::array<int64_t, Ndim> strides_;
     size_t size_{};
@@ -43,7 +46,8 @@ template <typename T, int64_t Ndim = 2> class NDArray {
         : shape_(shape), strides_(c_strides<Ndim>(shape_)),
           size_(std::accumulate(shape_.begin(), shape_.end(), 1,
                                 std::multiplies<>())),
-          data_(new T[size_]) {};
+          data_(new T[size_]) {}
+
 
     /**
      * @brief Construct a new NDArray object with a shape and value.
@@ -69,8 +73,8 @@ template <typename T, int64_t Ndim = 2> class NDArray {
     NDArray(NDArray &&other) noexcept
         : shape_(other.shape_), strides_(c_strides<Ndim>(shape_)),
           size_(other.size_), data_(other.data_) {
+        other.reset(); // TODO! is this necessary?
 
-        other.reset();
     }
 
     // Copy constructor
@@ -78,7 +82,14 @@ template <typename T, int64_t Ndim = 2> class NDArray {
         : shape_(other.shape_), strides_(c_strides<Ndim>(shape_)),
           size_(other.size_), data_(new T[size_]) {
         std::copy(other.data_, other.data_ + size_, data_);
-        // fmt::print("NDArray(const NDArray &other)\n");
+    }
+
+    // Conversion operator from array expression to array
+    template <typename E>
+    NDArray(ArrayExpr<E, Ndim> &&expr) : NDArray(expr.shape()) {
+        for (int i = 0; i < size_; ++i) {
+            data_[i] = expr[i];
+        }
     }
 
     ~NDArray() { delete[] data_; }
@@ -90,14 +101,10 @@ template <typename T, int64_t Ndim = 2> class NDArray {
 
     NDArray &operator=(NDArray &&other) noexcept; // Move assign
     NDArray &operator=(const NDArray &other);     // Copy assign
-
-    NDArray operator+(const NDArray &other);
     NDArray &operator+=(const NDArray &other);
-    NDArray operator-(const NDArray &other);
     NDArray &operator-=(const NDArray &other);
-    NDArray operator*(const NDArray &other);
     NDArray &operator*=(const NDArray &other);
-    NDArray operator/(const NDArray &other);
+
     // NDArray& operator/=(const NDArray& other);
 
     template <typename V> NDArray &operator/=(const NDArray<V, Ndim> &other) {
@@ -151,12 +158,16 @@ template <typename T, int64_t Ndim = 2> class NDArray {
         return data_[element_offset(strides_, index...)];
     }
 
+    // TODO! is int the right type for index?
     T &operator()(int i) { return data_[i]; }
     const T &operator()(int i) const { return data_[i]; }
 
+    T &operator[](int i) { return data_[i]; }
+    const T &operator[](int i) const { return data_[i]; }
+
     T *data() { return data_; }
     std::byte *buffer() { return reinterpret_cast<std::byte *>(data_); }
-    uint64_t size() const { return size_; }
+    size_t size() const { return 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]; }
@@ -204,17 +215,11 @@ NDArray<T, Ndim>::operator=(NDArray<T, Ndim> &&other) noexcept {
     return *this;
 }
 
-template <typename T, int64_t Ndim>
-NDArray<T, Ndim> NDArray<T, Ndim>::operator+(const NDArray &other) {
-    NDArray result(*this);
-    result += other;
-    return result;
-}
 template <typename T, int64_t Ndim>
 NDArray<T, Ndim> &NDArray<T, Ndim>::operator+=(const NDArray<T, Ndim> &other) {
     // check shape
     if (shape_ == other.shape_) {
-        for (uint32_t i = 0; i < size_; ++i) {
+        for (size_t i = 0; i < size_; ++i) {
             data_[i] += other.data_[i];
         }
         return *this;
@@ -222,13 +227,6 @@ 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>
-NDArray<T, Ndim> NDArray<T, Ndim>::operator-(const NDArray &other) {
-    NDArray result{*this};
-    result -= other;
-    return result;
-}
-
 template <typename T, int64_t Ndim>
 NDArray<T, Ndim> &NDArray<T, Ndim>::operator-=(const NDArray<T, Ndim> &other) {
     // check shape
@@ -240,12 +238,6 @@ 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>
-NDArray<T, Ndim> NDArray<T, Ndim>::operator*(const NDArray &other) {
-    NDArray result = *this;
-    result *= other;
-    return result;
-}
 
 template <typename T, int64_t Ndim>
 NDArray<T, Ndim> &NDArray<T, Ndim>::operator*=(const NDArray<T, Ndim> &other) {
@@ -259,13 +251,6 @@ 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>
-NDArray<T, Ndim> NDArray<T, Ndim>::operator/(const NDArray &other) {
-    NDArray result = *this;
-    result /= other;
-    return result;
-}
-
 template <typename T, int64_t Ndim>
 NDArray<T, Ndim> &NDArray<T, Ndim>::operator&=(const T &mask) {
     for (auto it = begin(); it != end(); ++it)
@@ -276,7 +261,7 @@ NDArray<T, Ndim> &NDArray<T, Ndim>::operator&=(const T &mask) {
 template <typename T, int64_t Ndim>
 NDArray<bool, Ndim> NDArray<T, Ndim>::operator>(const NDArray &other) {
     if (shape_ == other.shape_) {
-        NDArray<bool> result{shape_};
+        NDArray<bool, Ndim> result{shape_};
         for (int i = 0; i < size_; ++i) {
             result(i) = (data_[i] > other.data_[i]);
         }

include/aare/NDView.hpp

@@ -1,4 +1,7 @@
 #pragma once
+
+#include "aare/ArrayExpr.hpp"
+
 #include <algorithm>
 #include <array>
 #include <cassert>
@@ -45,7 +48,7 @@ template <int64_t Ndim> std::array<int64_t, Ndim> make_array(const std::vector<i
     return arr;
 }
 
-template <typename T, int64_t Ndim = 2> class NDView {
+template <typename T, int64_t Ndim = 2> class NDView : public ArrayExpr<NDView<T, Ndim>, Ndim> {
   public:
     NDView() = default;
     ~NDView() = default;
@@ -68,7 +71,7 @@ template <typename T, int64_t Ndim = 2> class NDView {
         return buffer_[element_offset(strides_, index...)];
     }
 
-    uint64_t size() const { return size_; }
+    size_t size() const { return size_; }
     size_t total_bytes() const { return size_ * sizeof(T); }
     std::array<int64_t, Ndim> strides() const noexcept { return strides_; }
 

include/aare/defs.hpp

@@ -41,7 +41,7 @@ namespace aare {
 void assert_failed(const std::string &msg);
 
 
-class Cluster {
+class DynamicCluster {
   public:
     int cluster_sizeX;
     int cluster_sizeY;
@@ -53,36 +53,36 @@ class Cluster {
     std::byte *m_data;
 
   public:
-    Cluster(int cluster_sizeX_, int cluster_sizeY_,
+    DynamicCluster(int cluster_sizeX_, int cluster_sizeY_,
             Dtype dt_ = Dtype(typeid(int32_t)))
         : cluster_sizeX(cluster_sizeX_), cluster_sizeY(cluster_sizeY_),
           dt(dt_) {
         m_data = new std::byte[cluster_sizeX * cluster_sizeY * dt.bytes()]{};
     }
-    Cluster() : Cluster(3, 3) {}
+    DynamicCluster() : DynamicCluster(3, 3) {}
-    Cluster(const Cluster &other)
+    DynamicCluster(const DynamicCluster &other)
-        : Cluster(other.cluster_sizeX, other.cluster_sizeY, other.dt) {
+        : DynamicCluster(other.cluster_sizeX, other.cluster_sizeY, other.dt) {
         if (this == &other)
             return;
         x = other.x;
         y = other.y;
         memcpy(m_data, other.m_data, other.bytes());
     }
-    Cluster &operator=(const Cluster &other) {
+    DynamicCluster &operator=(const DynamicCluster &other) {
         if (this == &other)
             return *this;
-        this->~Cluster();
+        this->~DynamicCluster();
-        new (this) Cluster(other);
+        new (this) DynamicCluster(other);
         return *this;
     }
-    Cluster(Cluster &&other) noexcept
+    DynamicCluster(DynamicCluster &&other) noexcept
         : cluster_sizeX(other.cluster_sizeX),
           cluster_sizeY(other.cluster_sizeY), x(other.x), y(other.y),
           dt(other.dt), m_data(other.m_data) {
         other.m_data = nullptr;
         other.dt = Dtype(Dtype::TypeIndex::ERROR);
     }
-    ~Cluster() { delete[] m_data; }
+    ~DynamicCluster() { delete[] m_data; }
     template <typename T> T get(int idx) {
         (sizeof(T) == dt.bytes())
             ? 0
@@ -95,10 +95,6 @@ class Cluster {
             : throw std::invalid_argument("[ERROR] Type size mismatch");
         return memcpy(m_data + idx * dt.bytes(), &val, (size_t)dt.bytes());
     }
-    // auto x() const { return x; }
-    // auto y() const { return y; }
-    // auto x(int16_t x_) { return x = x_; }
-    // auto y(int16_t y_) { return y = y_; }
 
     template <typename T> std::string to_string() const {
         (sizeof(T) == dt.bytes())

pyproject.toml

@@ -4,7 +4,7 @@ build-backend = "scikit_build_core.build"
 
 [project]
 name = "aare"
-version = "2024.11.11.dev0"
+version = "2024.11.15.dev0"
 
 
 [tool.scikit-build]

python/aare/__init__.py

@@ -5,6 +5,7 @@ from . import _aare
 from ._aare import File, RawFile, RawMasterFile, RawSubFile
 from ._aare import Pedestal, ClusterFinder, VarClusterFinder
 from ._aare import DetectorType
+from ._aare import ClusterFile
 
 from .CtbRawFile import CtbRawFile
 from .ScanParameters import ScanParameters

python/examples/play.py

@@ -1,95 +1,15 @@
 import matplotlib.pyplot as plt
 import numpy as np
 plt.ion()
-
-import aare
-from aare import CtbRawFile
-print('aare imported')
-from aare import transform
-print('transform imported')
 from pathlib import Path
+from aare import ClusterFile
 
-import json
+base = Path('~/data/aare_test_data/clusters').expanduser()
 
-def decode(frames, rawdata):
+f = ClusterFile(base / 'beam_En700eV_-40deg_300V_10us_d0_f0_100.clust')
-    # rawdata = np.fromfile(f, dtype = np.uint16)
+# f = ClusterFile(base / 'single_frame_97_clustrers.clust')
-    counters = int((np.shape(rawdata)[0]/frames-56)/(48*48))
-    print('Counters:', counters)
-    rawdata = rawdata.reshape(frames,-1)[:,56:]
-    rawdata = rawdata.reshape(frames,576*counters,4) #Data come in "blocks" of 4 pixels/receiver
-    tr1 = rawdata[:,0:576*counters:2] #Transceiver1
-    tr1=tr1.reshape((frames,48*counters,24)) 
-
-    tr2 = rawdata[:,1:576*counters:2] #Transceiver2
-    tr2=tr2.reshape((frames,48*counters,24))     
-    
-    data = np.append(tr1,tr2,axis=2)
-    return data
-
-def get_Mh02_frames(fname):
-    # this function gives you the data from a file that is not a scan
-    # it returns a (frames,48*counters,48)
-
-    jsonf = open(fname)
-    jsonpar = json.load(jsonf)
-    jsonf.close()
-
-    frames=jsonpar["Frames in File"]
-    print('Frames:', frames)
-
-    rawf = fname.replace('master','d0_f0')	
-    rawf = rawf.replace('.json','.raw')	
-
-    with open(rawf, 'rb') as f:
-        rawdata = np.fromfile(f, dtype = np.uint16)
-        data = decode(frames, rawdata) 
-        print('Data:', np.shape(data))
-        
-        return data
 
 
-#target format
+for i in range(10):
-# [frame, counter, row, col]
+    fn, cl = f.read_frame()
-# plt.imshow(data[0,0])
+    print(fn, cl.size)
-
-
-base = Path('/mnt/sls_det_storage/matterhorn_data/aare_test_data/ci/aare_test_data')
-# p = Path(base / 'jungfrau/jungfrau_single_master_0.json')
-
-# f = aare.File(p)
-# for i in range(10):
-#     frame = f.read_frame()
-
-
-
-# # f2 = aare.CtbRawFile(fpath, transform=transform.matterhorn02)
-# # header, data = f2.read()
-# # plt.plot(data[:,0,20,20])
-# from aare import RawMasterFile, File, RawSubFile, DetectorType, RawFile
-# base = Path('/mnt/sls_det_storage/matterhorn_data/aare_test_data/Jungfrau10/Jungfrau_DoubleModule_1UDP_ROI/SideBySide/')
-# fpath = Path('241019_JF_12keV_Si_FF_GaAs_FF_7p88mmFilter_PedestalStart_ZPos_5.5_master_0.json')
-# raw = Path('241019_JF_12keV_Si_FF_GaAs_FF_7p88mmFilter_PedestalStart_ZPos_5.5_d0_f0_0.raw')
-
-
-
-# m = RawMasterFile(base / fpath)
-# # roi = m.roi
-# # rows = roi.ymax-roi.ymin+1
-# # cols = 1024-roi.xmin
-# # sf = RawSubFile(base / raw, DetectorType.Jungfrau, rows, cols, 16)
-
-from aare import RawFile
-
-
-
-from aare import RawFile, File
-
-base = Path('/mnt/sls_det_storage/matterhorn_data/aare_test_data/Jungfrau10/Jungfrau_DoubleModule_1UDP_ROI/')
-fname = base / Path('SideBySide/241019_JF_12keV_Si_FF_GaAs_FF_7p88mmFilter_PedestalStart_ZPos_5.5_master_0.json')
-# fname = Path(base / 'jungfrau/jungfrau_single_master_0.json')
-# fname = base / 'Stacked/241024_JF10_m450_m367_KnifeEdge_TestBesom_9keV_750umFilter_PedestalStart_ZPos_-6_master_0.json'
-
-
-f = RawFile(fname)
-h,img = f.read_frame()
-print(f'{h["frameNumber"]}')

python/src/cluster.hpp

@@ -33,20 +33,20 @@ void define_cluster_finder_bindings(py::module &m) {
                  return clusters;
              });
 
-    py::class_<Cluster>(m, "Cluster", py::buffer_protocol())
+    py::class_<DynamicCluster>(m, "DynamicCluster", py::buffer_protocol())
         .def(py::init<int, int, Dtype>())
-        .def("size", &Cluster::size)
+        .def("size", &DynamicCluster::size)
-        .def("begin", &Cluster::begin)
+        .def("begin", &DynamicCluster::begin)
-        .def("end", &Cluster::end)
+        .def("end", &DynamicCluster::end)
-        .def_readwrite("x", &Cluster::x)
+        .def_readwrite("x", &DynamicCluster::x)
-        .def_readwrite("y", &Cluster::y)
+        .def_readwrite("y", &DynamicCluster::y)
-        .def_buffer([](Cluster &c) -> py::buffer_info {
+        .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 Cluster &a) {
+        .def("__repr__", [](const DynamicCluster &a) {
-            return "<Cluster: x: " + std::to_string(a.x) +
+            return "<DynamicCluster: x: " + std::to_string(a.x) +
                    ", y: " + std::to_string(a.y) + ">";
         });
 }

python/src/cluster_file.hpp

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

python/src/file.hpp

@@ -264,7 +264,7 @@ void define_file_io_bindings(py::module &m) {
 
     //     .def("close", &ClusterFileV2::close);
 
-    // m.def("to_clustV2", [](std::vector<Cluster> &clusters, const int
+    // m.def("to_clustV2", [](std::vector<DynamicCluster> &clusters, const int
     // frame_number) {
     //     std::vector<ClusterV2> clusters_;
     //     for (auto &c : clusters) {

python/src/module.cpp

@@ -6,6 +6,7 @@
 #include "pixel_map.hpp"
 #include "pedestal.hpp"
 #include "cluster.hpp"
+#include "cluster_file.hpp"
 
 //Pybind stuff
 #include <pybind11/pybind11.h>
@@ -22,4 +23,5 @@ PYBIND11_MODULE(_aare, m) {
     define_pedestal_bindings<double>(m, "Pedestal");
     define_pedestal_bindings<float>(m, "Pedestal_float32");
     define_cluster_finder_bindings(m);
+    define_cluster_file_io_bindings(m);
 }

src/ClusterFile.cpp

@@ -0,0 +1,312 @@
+#include "aare/ClusterFile.hpp"
+
+namespace aare {
+
+ClusterFile::ClusterFile(const std::filesystem::path &fname, size_t chunk_size): m_chunk_size(chunk_size) {
+    fp = fopen(fname.c_str(), "rb");
+    if (!fp) {
+        throw std::runtime_error("Could not open file: " + fname.string());
+    }
+}
+
+std::vector<Cluster> ClusterFile::read_clusters(size_t n_clusters) {
+    std::vector<Cluster> clusters(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 = reinterpret_cast<Cluster *>(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((void *)(buf + nph_read), sizeof(Cluster), 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)) {
+            // 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((void *)(buf + nph_read), sizeof(Cluster), nn, fp);
+                m_num_left = nph - nn;
+            }
+            if (nph_read >= n_clusters)
+                break;
+        }
+    }
+
+    // Resize the vector to the number of clusters.
+    // No new allocation, only change bounds.
+    clusters.resize(nph_read);
+    return clusters;
+}
+
+std::vector<Cluster> ClusterFile::read_frame(int32_t &out_fnum) {
+    if (m_num_left) {
+        throw std::runtime_error("There are still photons left in the last frame");
+    }
+
+    if (fread(&out_fnum, sizeof(out_fnum), 1, fp) != 1) {
+        throw std::runtime_error("Could not read frame number");
+    }
+
+    int n_clusters;
+    if (fread(&n_clusters, sizeof(n_clusters), 1, fp) != 1) {
+        throw std::runtime_error("Could not read number of clusters");
+    }
+    std::vector<Cluster> clusters(n_clusters);
+
+    if (fread(clusters.data(), sizeof(Cluster), n_clusters, fp) != n_clusters) {
+        throw std::runtime_error("Could not read clusters");
+    }
+    return clusters;
+
+}
+
+std::vector<Cluster> ClusterFile::read_cluster_with_cut(size_t n_clusters,
+                                                        double *noise_map,
+                                                        int nx, int ny) {
+
+    std::vector<Cluster> clusters(n_clusters);
+    // size_t read_clusters_with_cut(FILE *fp, size_t n_clusters, Cluster *buf,
+    //                               uint32_t *n_left, double *noise_map, int
+    //                               nx, int ny) {
+    int iframe = 0;
+    //     uint32_t nph = *n_left;
+    uint32_t nph = m_num_left;
+    //     uint32_t nn = *n_left;
+    uint32_t nn = m_num_left;
+    size_t nph_read = 0;
+
+    int32_t t2max, tot1;
+    int32_t tot3;
+    // Cluster *ptr = buf;
+    Cluster *ptr = clusters.data();
+    int good = 1;
+    double noise;
+    // read photons left from previous frame
+    if (noise_map)
+        printf("Using noise map\n");
+
+    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;
+        }
+        for (size_t iph = 0; iph < nn; iph++) {
+            // read photons 1 by 1
+            size_t n_read = fread((void *)(ptr), sizeof(Cluster), 1, fp);
+            if (n_read != 1) {
+                clusters.resize(nph_read);
+                return clusters;
+            }
+            // TODO! error handling on read
+            good = 1;
+            if (noise_map) {
+                if (ptr->x >= 0 && ptr->x < nx && ptr->y >= 0 && ptr->y < ny) {
+                    tot1 = ptr->data[4];
+                    analyze_cluster(*ptr, &t2max, &tot3, NULL, NULL, NULL, NULL,
+                                    NULL);
+                    noise = noise_map[ptr->y * nx + ptr->x];
+                    if (tot1 > noise || t2max > 2 * noise || tot3 > 3 * noise) {
+                        ;
+                    } else {
+                        good = 0;
+                        printf("%d %d %f %d %d %d\n", ptr->x, ptr->y, noise,
+                               tot1, t2max, tot3);
+                    }
+                } else {
+                    printf("Bad pixel number %d %d\n", ptr->x, ptr->y);
+                    good = 0;
+                }
+            }
+            if (good) {
+                ptr++;
+                nph_read++;
+            }
+            (m_num_left)--;
+            if (nph_read >= n_clusters)
+                break;
+        }
+    }
+        if (nph_read < n_clusters) {
+    //         // keep on reading frames and photons until reaching n_clusters
+            while (fread(&iframe, sizeof(iframe), 1, fp)) {
+    //             // printf("%d\n",nph_read);
+
+                if (fread(&nph, sizeof(nph), 1, fp)) {
+    //                 // printf("** %d\n",nph);
+                    m_num_left = nph;
+                    for (size_t iph = 0; iph < nph; iph++) {
+    //                     // read photons 1 by 1
+                        size_t n_read =
+                            fread((void *)(ptr), sizeof(Cluster), 1, fp);
+                        if (n_read != 1) {
+                            clusters.resize(nph_read);
+                            return clusters;
+                            // return nph_read;
+                        }
+                        good = 1;
+                        if (noise_map) {
+                            if (ptr->x >= 0 && ptr->x < nx && ptr->y >= 0 &&
+                                ptr->y < ny) {
+                                tot1 = ptr->data[4];
+                                analyze_cluster(*ptr, &t2max, &tot3, NULL,
+                                NULL,
+                                                NULL, NULL, NULL);
+                                // noise = noise_map[ptr->y * nx + ptr->x];
+                                noise = noise_map[ptr->y + ny * ptr->x];
+    			    if (tot1 > noise || t2max > 2 * noise ||
+                                    tot3 > 3 * noise) {
+                                    ;
+                                } else
+                                    good = 0;
+                            } else {
+                                printf("Bad pixel number %d %d\n", ptr->x,
+                                ptr->y); good = 0;
+                            }
+                        }
+                        if (good) {
+                            ptr++;
+    			nph_read++;
+                        }
+    		    (m_num_left)--;
+                        if (nph_read >= n_clusters)
+                            break;
+                    }
+                }
+                if (nph_read >= n_clusters)
+                    break;
+            }
+        }
+        // printf("%d\n",nph_read);
+        clusters.resize(nph_read);
+        return clusters;
+
+}
+
+int ClusterFile::analyze_cluster(Cluster cl, int32_t *t2, int32_t *t3, char *quad,
+                    double *eta2x, double *eta2y, double *eta3x,
+                    double *eta3y) {
+
+    return analyze_data(cl.data, t2, t3, quad, eta2x, eta2y, eta3x, eta3y);
+}
+
+int ClusterFile::analyze_data(int32_t *data, int32_t *t2, int32_t *t3, char *quad,
+                 double *eta2x, double *eta2y, double *eta3x, double *eta3y) {
+
+    int ok = 1;
+
+    int32_t tot2[4];
+    int32_t t2max = 0;
+    char c = 0;
+    int32_t val, tot3;
+
+    tot3 = 0;
+    for (int i = 0; i < 4; i++)
+        tot2[i] = 0;
+
+    for (int ix = 0; ix < 3; ix++) {
+        for (int iy = 0; iy < 3; iy++) {
+            val = data[iy * 3 + ix];
+            //	printf ("%d ",data[iy * 3 + ix]);
+            tot3 += val;
+            if (ix <= 1 && iy <= 1)
+                tot2[cBottomLeft] += val;
+            if (ix >= 1 && iy <= 1)
+                tot2[cBottomRight] += val;
+            if (ix <= 1 && iy >= 1)
+                tot2[cTopLeft] += val;
+            if (ix >= 1 && iy >= 1)
+                tot2[cTopRight] += val;
+        }
+        //	printf ("\n");
+    }
+    // printf ("\n");
+
+    if (t2 || quad) {
+
+        t2max = tot2[0];
+        c = cBottomLeft;
+        for (int i = 1; i < 4; i++) {
+            if (tot2[i] > t2max) {
+                t2max = tot2[i];
+                c = i;
+            }
+        }
+	//printf("*** %d %d %d %d -- %d\n",tot2[0],tot2[1],tot2[2],tot2[3],t2max);
+        if (quad)
+            *quad = c;
+        if (t2)
+            *t2 = t2max;
+    }
+    if (t3)
+        *t3 = tot3;
+
+    if (eta2x || eta2y) {
+        if (eta2x)
+            *eta2x = 0;
+        if (eta2y)
+            *eta2y = 0;
+        switch (c) {
+        case cBottomLeft:
+            if (eta2x && (data[3] + data[4]) != 0)
+                *eta2x = (double)(data[4]) / (data[3] + data[4]);
+            if (eta2y && (data[1] + data[4]) != 0)
+                *eta2y = (double)(data[4]) / (data[1] + data[4]);
+            break;
+        case cBottomRight:
+            if (eta2x && (data[2] + data[5]) != 0)
+                *eta2x = (double)(data[5]) / (data[4] + data[5]);
+            if (eta2y && (data[1] + data[4]) != 0)
+                *eta2y = (double)(data[4]) / (data[1] + data[4]);
+            break;
+        case cTopLeft:
+            if (eta2x && (data[7] + data[4]) != 0)
+                *eta2x = (double)(data[4]) / (data[3] + data[4]);
+            if (eta2y && (data[7] + data[4]) != 0)
+                *eta2y = (double)(data[7]) / (data[7] + data[4]);
+            break;
+        case cTopRight:
+            if (eta2x && t2max != 0)
+                *eta2x = (double)(data[5]) / (data[5] + data[4]);
+            if (eta2y && t2max != 0)
+                *eta2y = (double)(data[7]) / (data[7] + data[4]);
+            break;
+        default:;
+        }
+    }
+
+    if (eta3x || eta3y) {
+        if (eta3x && (data[3] + data[4] + data[5]) != 0)
+            *eta3x = (double)(-data[3] + data[3 + 2]) /
+                     (data[3] + data[4] + data[5]);
+        if (eta3y && (data[1] + data[4] + data[7]) != 0)
+            *eta3y = (double)(-data[1] + data[2 * 3 + 1]) /
+                     (data[1] + data[4] + data[7]);
+    }
+
+    return ok;
+}
+
+
+
+} // namespace aare

src/NDArray.test.cpp

@@ -1,5 +1,6 @@
 #include "aare/NDArray.hpp"
 #include <array>
+#include <catch2/benchmark/catch_benchmark.hpp>
 #include <catch2/catch_test_macros.hpp>
 
 using aare::NDArray;
@@ -101,7 +102,8 @@ TEST_CASE("Elementwise multiplication of 3D image") {
         a(i) = i;
         b(i) = i;
     }
-    auto c = a * b;
+    // auto c = a * b; // This works but the result is a lazy ArrayMul object
+    NDArray<double, 3> c = a * b;
     REQUIRE(c(0, 0, 0) == 0 * 0);
     REQUIRE(c(0, 0, 1) == 1 * 1);
     REQUIRE(c(0, 1, 1) == 3 * 3);
@@ -109,6 +111,39 @@ TEST_CASE("Elementwise multiplication of 3D image") {
     REQUIRE(c(2, 3, 1) == 23 * 23);
 }
 
+NDArray<int> MultiplyNDArrayUsingOperator(NDArray<int> &a, NDArray<int> &b) {
+    // return a * a * b * b;
+    NDArray<int>c =  a*b;
+    return c;
+}
+
+NDArray<int> MultiplyNDArrayUsingIndex(NDArray<int> &a, NDArray<int> &b) {
+    NDArray<int> res(a.shape());
+    for (uint32_t i = 0; i < a.size(); i++) {
+        // res(i) = a(i) * a(i) * b(i) * b(i);
+        res(i) = a(i) * b(i);
+    }
+    return res;
+}
+
+NDArray<int> AddNDArrayUsingOperator(NDArray<int> &a, NDArray<int> &b) {
+    // return a * a * b * b;
+    // NDArray<int>c = a+b;
+    NDArray<int> c(a.shape());
+    c = a + b;
+    return c;
+}
+
+NDArray<int> AddNDArrayUsingIndex(NDArray<int> &a, NDArray<int> &b) {
+    NDArray<int> res(a.shape());
+    for (uint32_t i = 0; i < a.size(); i++) {
+        // res(i) = a(i) * a(i) * b(i) * b(i);
+        res(i) = a(i) + b(i);
+    }
+    return res;
+}
+
+
 TEST_CASE("Compare two images") {
     NDArray<int> a;
     NDArray<int> b;
@@ -178,7 +213,8 @@ TEST_CASE("Elementwise operations on images") {
         NDArray<double> A(shape, a_val);
         NDArray<double> B(shape, b_val);
 
-        auto C = A + B;
+        NDArray<double> C = A + B;
+        // auto C = A+B; // This works but the result is a lazy ArraySum object
 
         // Value of C matches
         for (uint32_t i = 0; i < C.size(); ++i) {
@@ -202,7 +238,8 @@ TEST_CASE("Elementwise operations on images") {
     SECTION("Subtract two images") {
         NDArray<double> A(shape, a_val);
         NDArray<double> B(shape, b_val);
-        auto C = A - B;
+        NDArray<double> C = A - B;
+        // auto C = A - B; // This works but the result is a lazy ArraySub object
 
         // Value of C matches
         for (uint32_t i = 0; i < C.size(); ++i) {
@@ -226,7 +263,8 @@ TEST_CASE("Elementwise operations on images") {
     SECTION("Multiply two images") {
         NDArray<double> A(shape, a_val);
         NDArray<double> B(shape, b_val);
-        auto C = A * B;
+        // auto C = A * B; // This works but the result is a lazy ArrayMul object
+        NDArray<double> C = A * B;
 
         // Value of C matches
         for (uint32_t i = 0; i < C.size(); ++i) {
@@ -250,7 +288,8 @@ TEST_CASE("Elementwise operations on images") {
     SECTION("Divide two images") {
         NDArray<double> A(shape, a_val);
         NDArray<double> B(shape, b_val);
-        auto C = A / B;
+        // auto C = A / B; // This works but the result is a lazy ArrayDiv object
+        NDArray<double> C = A / B;
 
         // Value of C matches
         for (uint32_t i = 0; i < C.size(); ++i) {

src/NumpyFile.test.cpp

@@ -6,7 +6,7 @@
 
 using aare::Dtype;
 using aare::NumpyFile;
-TEST_CASE("Read a 1D numpy file with int32 data type") {
+TEST_CASE("Read a 1D numpy file with int32 data type", "[.integration]") {
 
     auto fpath = test_data_path() / "numpy" / "test_1d_int32.npy";
     REQUIRE(std::filesystem::exists(fpath));
@@ -24,7 +24,7 @@ TEST_CASE("Read a 1D numpy file with int32 data type") {
     }
 }
  
-TEST_CASE("Read a 3D numpy file with np.double data type") {
+TEST_CASE("Read a 3D numpy file with np.double data type", "[.integration]") {
 
     auto fpath = test_data_path() / "numpy" / "test_3d_double.npy";
     REQUIRE(std::filesystem::exists(fpath));

src/RawFile.test.cpp

@@ -7,7 +7,7 @@
 
 using aare::File;
 
-TEST_CASE("Read number of frames from a jungfrau raw file") {
+TEST_CASE("Read number of frames from a jungfrau raw file", "[.integration]") {
 
     auto fpath = test_data_path() / "jungfrau" / "jungfrau_single_master_0.json";
     REQUIRE(std::filesystem::exists(fpath));
@@ -16,7 +16,7 @@ TEST_CASE("Read number of frames from a jungfrau raw file") {
     REQUIRE(f.total_frames() == 10);
 }
 
-TEST_CASE("Read frame numbers from a jungfrau raw file") {
+TEST_CASE("Read frame numbers from a jungfrau raw file", "[.integration]") {
     auto fpath = test_data_path() / "jungfrau" / "jungfrau_single_master_0.json";
     REQUIRE(std::filesystem::exists(fpath));
 
@@ -32,7 +32,7 @@ TEST_CASE("Read frame numbers from a jungfrau raw file") {
     }
 }
 
-TEST_CASE("Read a frame number too high throws") {
+TEST_CASE("Read a frame number too high throws", "[.integration]") {
     auto fpath = test_data_path() / "jungfrau" / "jungfrau_single_master_0.json";
     REQUIRE(std::filesystem::exists(fpath));
 
@@ -46,7 +46,7 @@ TEST_CASE("Read a frame number too high throws") {
     REQUIRE_THROWS(f.frame_number(10));
 }
 
-TEST_CASE("Read a frame numbers where the subfile is missing throws") {
+TEST_CASE("Read a frame numbers where the subfile is missing throws", "[.integration]") {
     auto fpath = test_data_path() / "jungfrau" / "jungfrau_missing_subfile_master_0.json";
     REQUIRE(std::filesystem::exists(fpath));
 
@@ -67,7 +67,7 @@ TEST_CASE("Read a frame numbers where the subfile is missing throws") {
 }
 
 
-TEST_CASE("Read data from a jungfrau 500k single port raw file") {
+TEST_CASE("Read data from a jungfrau 500k single port raw file", "[.integration]") {
     auto fpath = test_data_path() / "jungfrau" / "jungfrau_single_master_0.json";
     REQUIRE(std::filesystem::exists(fpath));
 
@@ -83,7 +83,7 @@ TEST_CASE("Read data from a jungfrau 500k single port raw file") {
     }
 }
 
-TEST_CASE("Read frame numbers from a raw file") {
+TEST_CASE("Read frame numbers from a raw file", "[.integration]") {
     auto fpath = test_data_path() / "eiger" / "eiger_500k_16bit_master_0.json";
     REQUIRE(std::filesystem::exists(fpath));
 
@@ -96,7 +96,7 @@ TEST_CASE("Read frame numbers from a raw file") {
     }
 }
 
-TEST_CASE("Compare reading from a numpy file with a raw file") {
+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";
     REQUIRE(std::filesystem::exists(fpath_raw));
 
@@ -116,7 +116,7 @@ TEST_CASE("Compare reading from a numpy file with a raw file") {
     }
 }
 
-TEST_CASE("Read multipart files") {
+TEST_CASE("Read multipart files", "[.integration]") {
     auto fpath = test_data_path() / "jungfrau" / "jungfrau_double_master_0.json";
     REQUIRE(std::filesystem::exists(fpath));
 
@@ -141,7 +141,7 @@ TEST_CASE("Read multipart files") {
     }
 }
 
-TEST_CASE("Read file with unordered frames") {
+TEST_CASE("Read file with unordered frames", "[.integration]") {
     //TODO! Better explanation and error message
     auto fpath = test_data_path() / "mythen" / "scan242_master_3.raw";
     REQUIRE(std::filesystem::exists(fpath));

src/RawMasterFile.test.cpp

@@ -58,7 +58,7 @@ TEST_CASE("A disabled scan"){
 }
 
 
-TEST_CASE("Parse a master file in .json format"){
+TEST_CASE("Parse a master file in .json format", "[.integration]"){
     auto fpath = test_data_path() / "jungfrau" / "jungfrau_single_master_0.json";
     REQUIRE(std::filesystem::exists(fpath));
     RawMasterFile f(fpath);
@@ -139,7 +139,7 @@ TEST_CASE("Parse a master file in .json format"){
 
 }
 
-TEST_CASE("Parse a master file in .raw format"){
+TEST_CASE("Parse a master file in .raw format", "[.integration]"){
     
     auto fpath = test_data_path() / "moench/moench04_noise_200V_sto_both_100us_no_light_thresh_900_master_0.raw";
     REQUIRE(std::filesystem::exists(fpath));
@@ -204,7 +204,7 @@ TEST_CASE("Parse a master file in .raw format"){
 }
 
 
-TEST_CASE("Read eiger master file"){
+TEST_CASE("Read eiger master file", "[.integration]"){
 auto fpath = test_data_path() / "eiger" / "eiger_500k_32bit_master_0.json";
     REQUIRE(std::filesystem::exists(fpath));
     RawMasterFile f(fpath);

src/defs.test.cpp

@@ -9,14 +9,14 @@ TEST_CASE("Enum to string conversion") {
     REQUIRE(ToString(aare::DetectorType::Jungfrau) == "Jungfrau");
 }
 
-TEST_CASE("Cluster creation") {
+TEST_CASE("DynamicCluster creation") {
-    aare::Cluster c(13, 15);
+    aare::DynamicCluster c(13, 15);
     REQUIRE(c.cluster_sizeX == 13);
     REQUIRE(c.cluster_sizeY == 15);
     REQUIRE(c.dt == aare::Dtype(typeid(int32_t)));
     REQUIRE(c.data() != nullptr);
 
-    aare::Cluster c2(c);
+    aare::DynamicCluster c2(c);
     REQUIRE(c2.cluster_sizeX == 13);
     REQUIRE(c2.cluster_sizeY == 15);
     REQUIRE(c2.dt == aare::Dtype(typeid(int32_t)));
@@ -25,7 +25,7 @@ TEST_CASE("Cluster creation") {
 
 // TEST_CASE("cluster set and get data") {
 
-//     aare::Cluster c2(33, 44, aare::Dtype(typeid(double)));
+//     aare::DynamicCluster c2(33, 44, aare::Dtype(typeid(double)));
 //     REQUIRE(c2.cluster_sizeX == 33);
 //     REQUIRE(c2.cluster_sizeY == 44);
 //     REQUIRE(c2.dt == aare::Dtype::DOUBLE);

tests/test.cpp

@@ -4,12 +4,12 @@
 #include <filesystem>
 #include <fstream>
 
-TEST_CASE("Test suite can find data assets") {
+TEST_CASE("Test suite can find data assets", "[.integration]") {
     auto fpath = test_data_path() / "numpy" / "test_numpy_file.npy";
     REQUIRE(std::filesystem::exists(fpath));
 }
 
-TEST_CASE("Test suite can open data assets") {
+TEST_CASE("Test suite can open data assets", "[.integration]") {
     auto fpath = test_data_path() / "numpy" / "test_numpy_file.npy";
     auto f = std::ifstream(fpath, std::ios::binary);
     REQUIRE(f.is_open());