mod pedestal

CMakeLists.txt

@@ -40,7 +40,7 @@ option(AARE_DOCS "Build documentation" OFF)
 option(AARE_VERBOSE "Verbose output" OFF)
 option(AARE_CUSTOM_ASSERT "Use custom assert" OFF)
 option(AARE_INSTALL_PYTHONEXT "Install the python extension in the install tree under CMAKE_INSTALL_PREFIX/aare/" OFF)
-
+option(AARE_ASAN "Enable AddressSanitizer" OFF)
 
 # Configure which of the dependencies to use FetchContent for
 option(AARE_FETCH_FMT "Use FetchContent to download fmt" ON)
@@ -225,13 +225,6 @@ if(CMAKE_BUILD_TYPE STREQUAL "Release")
     target_compile_options(aare_compiler_flags INTERFACE -O3)
 else()
     message(STATUS "Debug build")
-    target_compile_options(
-        aare_compiler_flags 
-        INTERFACE 
-            -Og 
-            -ggdb3 
-    )
-
 endif()
 
 # Common flags for GCC and Clang
@@ -256,7 +249,21 @@ target_compile_options(
 endif() #GCC/Clang specific
 
 
-
+if(AARE_ASAN)
+    message(STATUS "AddressSanitizer enabled")
+    target_compile_options(
+        aare_compiler_flags 
+        INTERFACE 
+            -fsanitize=address,undefined,pointer-compare 
+            -fno-omit-frame-pointer
+    )
+    target_link_libraries(
+        aare_compiler_flags 
+        INTERFACE 
+            -fsanitize=address,undefined,pointer-compare 
+            -fno-omit-frame-pointer
+    )
+endif()
 
 
 
@@ -316,6 +323,8 @@ target_include_directories(aare_core PUBLIC
     "$<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>"
 )
 
+
+
 target_link_libraries(
     aare_core 
     PUBLIC 

include/aare/ClusterFinder.hpp

@@ -82,8 +82,8 @@ class ClusterFinder {
         // // TODO! deal with even size clusters
         // // currently 3,3 -> +/- 1
         // //  4,4 -> +/- 2
-        short dy = m_cluster_sizeY / 2;
+        int dy = m_cluster_sizeY / 2;
-        short dx = m_cluster_sizeX / 2;
+        int dx = m_cluster_sizeX / 2;
 
         std::vector<CT> cluster_data(m_cluster_sizeX * m_cluster_sizeY);
         for (int iy = 0; iy < frame.shape(0); iy++) {
@@ -100,8 +100,8 @@ class ClusterFinder {
                     continue; // NEGATIVE_PEDESTAL go to next pixel
                               // TODO! No pedestal update???
 
-                for (short ir = -dy; ir < dy + 1; ir++) {
+                for (int ir = -dy; ir < dy + 1; ir++) {
-                    for (short ic = -dx; ic < dx + 1; ic++) {
+                    for (int ic = -dx; ic < dx + 1; ic++) {
                         if (ix + ic >= 0 && ix + ic < frame.shape(1) &&
                             iy + ir >= 0 && iy + ir < frame.shape(0)) {
                             PEDESTAL_TYPE val =

include/aare/ClusterVector.hpp

@@ -13,12 +13,12 @@ namespace aare {
  * contiguous memory buffer to store the clusters. 
  * @note push_back can invalidate pointers to elements in the container
  * @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> class ClusterVector {
+template <typename T, typename CoordType=int16_t> class ClusterVector {
     using value_type = T;
-    using coord_t = int16_t;
+    size_t m_cluster_size_x;
-    coord_t m_cluster_size_x;
+    size_t m_cluster_size_y;
-    coord_t m_cluster_size_y;
     std::byte *m_data{};
     size_t m_size{0};
     size_t m_capacity;
@@ -39,7 +39,7 @@ template <typename T> class ClusterVector {
      * @param cluster_size_y size of the cluster in y direction
      * @param capacity initial capacity of the buffer in number of clusters
      */
-    ClusterVector(coord_t cluster_size_x, coord_t cluster_size_y,
+    ClusterVector(size_t cluster_size_x, size_t cluster_size_y,
                   size_t capacity = 1024)
         : m_cluster_size_x(cluster_size_x), m_cluster_size_y(cluster_size_y),
           m_capacity(capacity) {
@@ -94,21 +94,22 @@ template <typename T> class ClusterVector {
      * @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(coord_t x, coord_t y, const std::byte *data) {
+    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<coord_t *>(ptr) = x;
+        *reinterpret_cast<CoordType *>(ptr) = x;
-        ptr += sizeof(coord_t);
+        ptr += sizeof(CoordType);
-        *reinterpret_cast<coord_t *>(ptr) = y;
+        *reinterpret_cast<CoordType *>(ptr) = y;
-        ptr += sizeof(coord_t);
+        ptr += sizeof(CoordType);
 
         std::copy(data, data + m_cluster_size_x * m_cluster_size_y * sizeof(T),
                   ptr);
         m_size++;
     }
 
+
     /**
      * @brief Sum the pixels in each cluster
      * @return std::vector<T> vector of sums for each cluster
@@ -117,7 +118,7 @@ template <typename T> class ClusterVector {
         std::vector<T> sums(m_size);
         const size_t stride = element_offset();
         const size_t n_pixels = m_cluster_size_x * m_cluster_size_y;
-        std::byte *ptr = m_data + 2 * sizeof(coord_t); // skip x and y
+        std::byte *ptr = m_data + 2 * sizeof(CoordType); // skip x and y
 
         for (size_t i = 0; i < m_size; i++) {
             sums[i] =
@@ -135,7 +136,7 @@ template <typename T> class ClusterVector {
      * @brief Return the offset in bytes for a single cluster
      */
     size_t element_offset() const {
-        return sizeof(m_cluster_size_x) + sizeof(m_cluster_size_y) +
+        return 2*sizeof(CoordType)  +
                m_cluster_size_x * m_cluster_size_y * sizeof(T);
     }
     /**
@@ -148,8 +149,8 @@ template <typename T> class ClusterVector {
      */
     std::byte *element_ptr(size_t i) { return m_data + element_offset(i); }
 
-    int16_t cluster_size_x() const { return m_cluster_size_x; }
+    size_t cluster_size_x() const { return m_cluster_size_x; }
-    int16_t cluster_size_y() const { return m_cluster_size_y; }
+    size_t cluster_size_y() const { return m_cluster_size_y; }
 
     std::byte *data() { return m_data; }
     const std::byte *data() const { return m_data; }

include/aare/NDArray.hpp

@@ -87,7 +87,7 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
     // Conversion operator from array expression to array
     template <typename E>
     NDArray(ArrayExpr<E, Ndim> &&expr) : NDArray(expr.shape()) {
-        for (int i = 0; i < size_; ++i) {
+        for (size_t i = 0; i < size_; ++i) {
             data_[i] = expr[i];
         }
     }
@@ -159,11 +159,11 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
     }
 
     // TODO! is int the right type for index?
-    T &operator()(int i) { return data_[i]; }
+    T &operator()(int64_t i) { return data_[i]; }
-    const T &operator()(int i) const { return data_[i]; }
+    const T &operator()(int64_t i) const { return data_[i]; }
 
-    T &operator[](int i) { return data_[i]; }
+    T &operator[](int64_t i) { return data_[i]; }
-    const T &operator[](int i) const { return data_[i]; }
+    const T &operator[](int64_t i) const { return data_[i]; }
 
     T *data() { return data_; }
     std::byte *buffer() { return reinterpret_cast<std::byte *>(data_); }

include/aare/Pedestal.hpp

@@ -23,31 +23,43 @@ template <typename SUM_TYPE = double> class Pedestal {
     NDArray<SUM_TYPE, 2> m_sum;
     NDArray<SUM_TYPE, 2> m_sum2;
 
+    //Cache mean since it is used over and over in the ClusterFinder
+    //This optimization is related to the access pattern of the ClusterFinder
+    //Relies on having more reads than pushes to the pedestal
+    NDArray<SUM_TYPE, 2> m_mean; 
+
   public:
     Pedestal(uint32_t rows, uint32_t cols, uint32_t n_samples = 1000)
         : m_rows(rows), m_cols(cols), m_samples(n_samples),
           m_cur_samples(NDArray<uint32_t, 2>({rows, cols}, 0)),
           m_sum(NDArray<SUM_TYPE, 2>({rows, cols})),
-          m_sum2(NDArray<SUM_TYPE, 2>({rows, cols})) {
+          m_sum2(NDArray<SUM_TYPE, 2>({rows, cols})),
+          m_mean(NDArray<SUM_TYPE, 2>({rows, cols})) {
         assert(rows > 0 && cols > 0 && n_samples > 0);
         m_sum = 0;
         m_sum2 = 0;
+        m_mean = 0;
     }
     ~Pedestal() = default;
 
     NDArray<SUM_TYPE, 2> mean() {
-        NDArray<SUM_TYPE, 2> mean_array({m_rows, m_cols});
+        return m_mean;
-        for (uint32_t i = 0; i < m_rows * m_cols; i++) {
-            mean_array(i / m_cols, i % m_cols) = mean(i / m_cols, i % m_cols);
-        }
-        return mean_array;
     }
 
     SUM_TYPE mean(const uint32_t row, const uint32_t col) const {
+        return m_mean(row, col);
+    }
+
+    SUM_TYPE std(const uint32_t row, const uint32_t col) const {
+        return std::sqrt(variance(row, col));
+    }
+
+    SUM_TYPE variance(const uint32_t row, const uint32_t col) const {
         if (m_cur_samples(row, col) == 0) {
             return 0.0;
         }
-        return m_sum(row, col) / m_cur_samples(row, col);
+        return m_sum2(row, col) / m_cur_samples(row, col) -
+               mean(row, col) * mean(row, col);
     }
 
     NDArray<SUM_TYPE, 2> variance() {
@@ -59,13 +71,7 @@ template <typename SUM_TYPE = double> class Pedestal {
         return variance_array;
     }
 
-    SUM_TYPE variance(const uint32_t row, const uint32_t col) const {
+    
-        if (m_cur_samples(row, col) == 0) {
-            return 0.0;
-        }
-        return m_sum2(row, col) / m_cur_samples(row, col) -
-               mean(row, col) * mean(row, col);
-    }
 
     NDArray<SUM_TYPE, 2> std() {
         NDArray<SUM_TYPE, 2> standard_deviation_array({m_rows, m_cols});
@@ -77,14 +83,12 @@ template <typename SUM_TYPE = double> class Pedestal {
         return standard_deviation_array;
     }
 
-    SUM_TYPE std(const uint32_t row, const uint32_t col) const {
+    
-        return std::sqrt(variance(row, col));
-    }
 
     void clear() {
-        for (uint32_t i = 0; i < m_rows * m_cols; i++) {
+        m_sum = 0;
-            clear(i / m_cols, i % m_cols);
+        m_sum2 = 0;
-        }
+        m_cur_samples = 0;
     }
 
     
@@ -104,8 +108,8 @@ template <typename SUM_TYPE = double> class Pedestal {
                 "Frame shape does not match pedestal shape");
         }
 
-        for (uint32_t row = 0; row < m_rows; row++) {
+        for (size_t row = 0; row < m_rows; row++) {
-            for (uint32_t col = 0; col < m_cols; col++) {
+            for (size_t col = 0; col < m_cols; col++) {
                 push<T>(row, col, frame(row, col));
             }
         }
@@ -134,18 +138,17 @@ template <typename SUM_TYPE = double> class Pedestal {
     template <typename T>
     void push(const uint32_t row, const uint32_t col, const T val_) {
         SUM_TYPE val = static_cast<SUM_TYPE>(val_);
-        const uint32_t idx = index(row, col);
+        if (m_cur_samples(row, col) < m_samples) {
-        if (m_cur_samples(idx) < m_samples) {
+            m_sum(row, col) += val;
-            m_sum(idx) += val;
+            m_sum2(row, col) += val * val;
-            m_sum2(idx) += val * val;
+            m_cur_samples(row, col)++;
-            m_cur_samples(idx)++;
         } else {
-            m_sum(idx) += val - m_sum(idx) / m_cur_samples(idx);
+            m_sum(row, col) += val - m_sum(row, col) / m_cur_samples(row, col);
-            m_sum2(idx) += val * val - m_sum2(idx) / m_cur_samples(idx);
+            m_sum2(row, col) += val * val - m_sum2(row, col) / m_cur_samples(row, col);
         }
+        //Since we just did a push we know that m_cur_samples(row, col) is at least 1
+        m_mean(row, col) = m_sum(row, col) / m_cur_samples(row, col);
     }
-    uint32_t index(const uint32_t row, const uint32_t col) const {
+
-        return row * m_cols + col;
-    };
 };
 } // namespace aare

python/aare/__init__.py

@@ -3,9 +3,10 @@ from . import _aare
 
 
 from ._aare import File, RawMasterFile, RawSubFile
-from ._aare import Pedestal, ClusterFinder, VarClusterFinder
+from ._aare import Pedestal_d, Pedestal_f, ClusterFinder, VarClusterFinder
 from ._aare import DetectorType
 from ._aare import ClusterFile
+from ._aare import hitmap
 
 from .CtbRawFile import CtbRawFile
 from .RawFile import RawFile

python/src/cluster.hpp

@@ -80,7 +80,26 @@ void define_cluster_finder_bindings(py::module &m) {
                  return;
              });
 
+    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.element_offset();
+        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");

python/src/cluster_file.hpp

@@ -10,6 +10,12 @@
 #include <pybind11/stl/filesystem.h>
 #include <string>
 
+//Disable warnings for unused parameters, as we ignore some
+//in the __exit__ method
+#pragma GCC diagnostic push
+#pragma GCC diagnostic ignored "-Wunused-parameter"
+
+
 namespace py = pybind11;
 using namespace ::aare;
 
@@ -61,3 +67,5 @@ void define_cluster_file_io_bindings(py::module &m) {
             return return_vector(vec);
         });
 }
+
+#pragma GCC diagnostic pop

python/src/module.cpp

@@ -22,8 +22,8 @@ PYBIND11_MODULE(_aare, m) {
     define_raw_master_file_bindings(m);
     define_var_cluster_finder_bindings(m);
     define_pixel_map_bindings(m);
-    define_pedestal_bindings<double>(m, "Pedestal");
+    define_pedestal_bindings<double>(m, "Pedestal_d");
-    define_pedestal_bindings<float>(m, "Pedestal_float32");
+    define_pedestal_bindings<float>(m, "Pedestal_f");
     define_cluster_finder_bindings(m);
     define_cluster_file_io_bindings(m);
 }

src/ClusterVector.test.cpp

@@ -1,7 +1,7 @@
 #include <cstdint>
 #include "aare/ClusterVector.hpp"
 
-// #include <catch2/matchers/catch_matchers_floating_point.hpp>
+#include <catch2/matchers/catch_matchers_floating_point.hpp>
 #include <catch2/catch_test_macros.hpp>
 
 using aare::ClusterVector;
@@ -44,10 +44,10 @@ TEST_CASE("Summing 3x1 clusters of int64"){
     struct Cluster_l3x1{
         int16_t x;
         int16_t y;
-        int64_t data[3];
+        int32_t data[3];
     };
 
-    ClusterVector<int64_t> cv(3, 1, 2);
+    ClusterVector<int32_t> cv(3, 1, 2);
     REQUIRE(cv.capacity() == 2);
     REQUIRE(cv.size() == 0);
     REQUIRE(cv.cluster_size_x() == 3);
@@ -75,3 +75,34 @@ TEST_CASE("Summing 3x1 clusters of int64"){
     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, 2);
+    REQUIRE(cv.capacity() == 2);
+    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]));
+    REQUIRE(cv.capacity() == 2);
+    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]));
+    REQUIRE(cv.capacity() == 2);
+    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));
+}

src/Frame.test.cpp

@@ -19,7 +19,7 @@ TEST_CASE("Construct a frame") {
     // data should be initialized to 0
     for (size_t i = 0; i < rows; i++) {
         for (size_t j = 0; j < cols; j++) {
-            uint8_t *data = (uint8_t *)frame.pixel_ptr(i, j);
+            uint8_t *data = reinterpret_cast<uint8_t *>(frame.pixel_ptr(i, j));
             REQUIRE(data != nullptr);
             REQUIRE(*data == 0);
         }
@@ -40,7 +40,7 @@ TEST_CASE("Set a value in a 8 bit frame") {
     // only the value we did set should be non-zero
     for (size_t i = 0; i < rows; i++) {
         for (size_t j = 0; j < cols; j++) {
-            uint8_t *data = (uint8_t *)frame.pixel_ptr(i, j);
+            uint8_t *data = reinterpret_cast<uint8_t *>(frame.pixel_ptr(i, j));
             REQUIRE(data != nullptr);
             if (i == 5 && j == 7) {
                 REQUIRE(*data == value);
@@ -65,7 +65,7 @@ TEST_CASE("Set a value in a 64 bit frame") {
     // only the value we did set should be non-zero
     for (size_t i = 0; i < rows; i++) {
         for (size_t j = 0; j < cols; j++) {
-            uint64_t *data = (uint64_t *)frame.pixel_ptr(i, j);
+            uint64_t *data = reinterpret_cast<uint64_t *>(frame.pixel_ptr(i, j));
             REQUIRE(data != nullptr);
             if (i == 5 && j == 7) {
                 REQUIRE(*data == value);
@@ -150,3 +150,4 @@ TEST_CASE("test explicit copy constructor") {
     REQUIRE(frame2.bytes() == rows * cols * bitdepth / 8);
     REQUIRE(frame2.data() != data);
 }
+

tests/CMakeLists.txt

@@ -17,8 +17,8 @@ endif()
 list(APPEND CMAKE_MODULE_PATH ${Catch2_SOURCE_DIR}/extras)
 
 add_executable(tests test.cpp)
-target_link_libraries(tests PRIVATE Catch2::Catch2WithMain)
+target_link_libraries(tests PRIVATE Catch2::Catch2WithMain aare_core aare_compiler_flags)
-
+# target_compile_options(tests PRIVATE -fno-omit-frame-pointer -fsanitize=address)
 set_target_properties(tests PROPERTIES
         RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}
         OUTPUT_NAME run_tests
@@ -34,7 +34,7 @@ set(TestSources
 target_sources(tests PRIVATE ${TestSources} )
 
 #Work around to remove, this is not the way to do it =) 
-target_link_libraries(tests PRIVATE aare_core aare_compiler_flags)
+# target_link_libraries(tests PRIVATE aare_core aare_compiler_flags)
 
 
 #configure a header to pass test file paths

tests/test.cpp

@@ -3,6 +3,7 @@
 #include <climits>
 #include <filesystem>
 #include <fstream>
+#include <fmt/core.h>
 
 TEST_CASE("Test suite can find data assets", "[.integration]") {
     auto fpath = test_data_path() / "numpy" / "test_numpy_file.npy";
@@ -19,3 +20,19 @@ TEST_CASE("Test float32 and char8") {
     REQUIRE(sizeof(float) == 4);
     REQUIRE(CHAR_BIT == 8);
 }
+
+/**
+ * Uncomment the following tests to verify that asan is working
+ */
+
+// TEST_CASE("trigger asan stack"){
+//     int arr[5] =  {1,2,3,4,5};
+//     int val = arr[7];
+//     fmt::print("val: {}\n", val);
+// }
+
+// TEST_CASE("trigger asan heap"){
+//     auto *ptr = new int[5];
+//     ptr[70] = 5;
+//     fmt::print("ptr: {}\n", ptr[70]);
+// }