pedestal (#67)

* add config files for multimodule receiving

* read subfiles with unordered and missing frames

* save work debugging

* Revert "save work debugging"

This reverts commit e791992a05efd754f93a80c980d17397eb4b6045.

* Revert "read subfiles with unordered and missing frames"

This reverts commit 1177fd129d3690db92e9597ccda62598e5a44d41.

* throw when two frames have different frame numbers

* write single part RawFile (working beta)

* correct total number of frames in master file

* add new mythen file with syncd frames

* save work

* save work for receive multimodule
multimodule config results in too much packet loss. needs more debugging.

* setup Task Distributiosn/ parallelization programming model

* read frames with same frame number

* clang-tidy fixes, formatting, add tests

* added second receiver

* Synchronize between zmq streams and merge frames

* improve readability in loop

* fix failing tests

* add simple test for merge frames

* restructure files and use top-level header

* working pedestal + tests

* test_pedestal statistics

* QA test pedestal, fix clang-tidy errors

---------

Co-authored-by: Bechir <bechir.brahem420@gmail.com>
Co-authored-by: Erik Frojdh <erik.frojdh@gmail.com>

.clang-tidy

@@ -43,7 +43,8 @@ Checks: '*,
 HeaderFilterRegex: '*.hpp'
 ImplementationFileExtensions: [cpp]
 AnalyzeTemporaryDtors: false
-WarningsAsErrors: '*'
+WarningsAsErrors: '*,
+-bugprone-implicit-widening-of-multiplication-result'
 FormatStyle:     none
 CheckOptions:
   - { key: readability-identifier-naming.NamespaceCase,       value: lower_case }

.github/workflows/common-workflow.yml

@@ -71,7 +71,14 @@ jobs:
           pwd
           export AARE_ROOT_DIR="$PWD"
           ls build/examples/*_example
-          find build/examples -name "*_example" -not -name "zmq_*" | xargs -I {} -n 1 -t  bash -c {}
+          # examples to run
+          cd build/examples
+          examples=(raw_example cluster_example json_example logger_example multiport_example mythen_example)
+          examples+=(numpy_read_example numpy_write_example)
+          for example in "${examples[@]}"; do
+            echo "Running example: $example"
+            ./$example
+          done
 
 
 

CMakeLists.txt

@@ -173,15 +173,12 @@ if(AARE_TESTS)
     add_subdirectory(tests)
 endif()
 
-add_subdirectory(core)
-add_subdirectory(file_io)
-add_subdirectory(utils)
-add_subdirectory(network_io)
+add_subdirectory(src)
 
 
 #Overall target to link to when using the library
 add_library(aare INTERFACE) 
-target_link_libraries(aare INTERFACE core file_io utils network_io ${Boost_LIBRARIES})
+target_link_libraries(aare INTERFACE core file_io utils network_io processing ${Boost_LIBRARIES})
 target_include_directories(aare INTERFACE 
     $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/include>
     $<INSTALL_INTERFACE:include>
@@ -190,9 +187,6 @@ target_include_directories(aare INTERFACE
 add_subdirectory(examples)
 
 
-if(AARE_PYTHON_BINDINGS)
-    add_subdirectory(python)
-endif()
 
 
 # custom target to run check formatting with clang-format
@@ -222,7 +216,7 @@ endif()
 
 add_custom_target(
     clang-tidy
-    COMMAND find \( -path "./core/*" -o -path "./file_io/*" -o -path "./network_io/*" -o -path "./utils/*" \) -name "*.cpp" -not -name "*.test.cpp" -not -name "CircularFifo.hpp" -not -name "ProducerConsumerQueue.hpp" -not -name "VariableSizeClusterFinder.hpp" | xargs -I {} -n 1 -P 10 bash -c "${CLANG_TIDY_COMMAND} --config-file=.clang-tidy -p build {}"
+    COMMAND find \( -path "./src/*" -a -not -path "./src/python/*" -a \( -name "*.cpp"   -not -name "*.test.cpp"\)  \)    -not -name "CircularFifo.hpp" -not -name "ProducerConsumerQueue.hpp" -not -name "VariableSizeClusterFinder.hpp" | xargs -I {} -n 1 -P 10 bash -c "${CLANG_TIDY_COMMAND} --config-file=.clang-tidy -p build {}"
     WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}
     COMMENT "linting with clang-tidy"
     VERBATIM

examples/CMakeLists.txt

@@ -2,10 +2,12 @@
 set(EXAMPLE_LIST "json_example;logger_example;numpy_read_example;multiport_example;raw_example;zmq_restream_example")
 set(EXAMPLE_LIST "${EXAMPLE_LIST};mythen_example;numpy_write_example;zmq_receiver_example;zmq_sender_example;")    
 set(EXAMPLE_LIST "${EXAMPLE_LIST};cluster_example;zmq_multi_receiver;zmq_task_ventilator;zmq_worker;zmq_sink")
+set(EXAMPLE_LIST "${EXAMPLE_LIST};testing_pedestal;")
+
 foreach(example ${EXAMPLE_LIST})
     add_executable(${example} ${example}.cpp)
     target_include_directories(${example} PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}/include)
-    target_link_libraries(${example} PUBLIC aare PRIVATE aare_compiler_flags libzmq)
+    target_link_libraries(${example} PUBLIC aare PRIVATE aare_compiler_flags libzmq )
 endforeach()
 
 

examples/cluster_example.cpp

@@ -1,5 +1,6 @@
-#include "aare/core/defs.hpp"
-#include "aare/file_io/ClusterFile.hpp"
+#include "aare.hpp"
+#include "aare/examples/defs.hpp"
+
 #include <cassert>
 #include <iostream>
 

examples/json_example.cpp

@@ -1,7 +1,6 @@
-// Your First C++ Program
+#include "aare.hpp"
 #include "aare/examples/defs.hpp"
-#include "aare/file_io/File.hpp"
-#include "aare/utils/logger.hpp"
+
 #include <iostream>
 
 using aare::File;

examples/logger_example.cpp

@@ -1,5 +1,5 @@
+#include "aare.hpp"
 #include "aare/examples/defs.hpp"
-#include "aare/utils/logger.hpp"
 
 #include <fstream>
 #include <iostream>

examples/multiport_example.cpp

@@ -1,7 +1,5 @@
-// Your First C++ Program
+#include "aare.hpp"
 #include "aare/examples/defs.hpp"
-#include "aare/file_io/File.hpp"
-#include "aare/utils/logger.hpp"
 
 #include <iostream>
 

examples/mythen_example.cpp

@@ -1,7 +1,5 @@
-// Your First C++ Program
+#include "aare.hpp"
 #include "aare/examples/defs.hpp"
-#include "aare/file_io/File.hpp"
-#include "aare/utils/logger.hpp"
 
 #include <iostream>
 

examples/numpy_read_example.cpp

@@ -1,6 +1,5 @@
-// Your First C++ Program
+#include "aare.hpp"
 #include "aare/examples/defs.hpp"
-#include "aare/file_io/File.hpp"
 
 #include <iostream>
 

examples/numpy_write_example.cpp

@@ -1,7 +1,5 @@
-// Your First C++ Program
-#include "aare/core/Frame.hpp"
+#include "aare.hpp"
 #include "aare/examples/defs.hpp"
-#include "aare/file_io/File.hpp"
 
 #include <iostream>
 

examples/raw_example.cpp

@@ -1,6 +1,6 @@
-// Your First C++ Program
-#include "aare/aare.hpp"
+#include "aare.hpp"
 #include "aare/examples/defs.hpp"
+
 #include <iostream>
 
 using namespace aare;

examples/testing_pedestal.cpp

@@ -0,0 +1,112 @@
+#include "aare.hpp"
+#include <iostream>
+
+using namespace std;
+using namespace aare;
+#include <algorithm>
+#include <chrono>
+#include <numeric>
+#include <random>
+#include <vector>
+
+void print_vpair(std::vector<std::pair<int, double>> &v) {
+    std::cout << "[ ";
+    for (int i = 0; i < v.size(); i++) {
+        std::cout << "(" << v[i].first << "," << v[i].second << "), ";
+    }
+    std::cout << "]" << std::endl;
+}
+int range(int min, int max, int i, int steps) { return min + (max - min) * i / steps; }
+int main() {
+    const int rows = 1, cols = 1;
+    double MEAN = 5.0, STD = 1.0, VAR = STD * STD, TOLERANCE = 0.1;
+
+    unsigned seed = std::chrono::system_clock::now().time_since_epoch().count();
+    std::default_random_engine generator(seed);
+    std::normal_distribution<double> distribution(MEAN, STD);
+
+    Pedestal pedestal(rows, cols, 1000);
+    std::vector<double> values;
+    std::vector<double> vmean, vvariance, pmean, pvariance, pstandard_deviation;
+    std::vector<int> samples;
+    std::vector<std::pair<int, double>> cur_mean, cur_variance;
+
+    // int steps = 1000;
+    // for (int n_iters = 0; n_iters < steps; n_iters++) {
+    //     int x=range(1000, 1000000, n_iters,steps);
+    //     std::cout<<x<<std::endl;
+
+    //     for (int i = 0; i < x; i++) {
+    //         Frame frame(rows, cols, 64);
+    //         for (int j = 0; j < rows; j++)
+    //             for (int k = 0; k < cols; k++) {
+    //                 double val = distribution(generator);
+    //                 frame.set<double>(j, k, val);
+    //             }
+    //         pedestal.push<double>(frame);
+    //     }
+    //     pmean.push_back({x, pedestal.mean(0, 0)});
+    //     pvariance.push_back({x, pedestal.variance(0, 0)});
+    // }
+
+    long double sum = 0;
+    long double sum2 = 0;
+    // fill 1000 first values of pedestal
+    for (int x = 0; x < 1000; x++) {
+        Frame frame(rows, cols, 64);
+        double val = distribution(generator);
+        frame.set<double>(0, 0, val);
+        pedestal.push<double>(frame);
+        values.push_back(val);
+        sum += val;
+        sum2 += val * val;
+    }
+
+    for (int x = 0, aa = 0; x < 100000; x++, aa++) {
+        Frame frame(rows, cols, 64);
+        double val = distribution(generator);
+        frame.set<double>(0, 0, val);
+        pedestal.push<double>(frame);
+
+        values.push_back(val);
+        auto old_val = values[values.size() - 1000 - 1];
+        sum += val - old_val;
+        sum2 += val * val - old_val * old_val;
+        if (aa % 100 == 1) {
+            aa = 2;
+            samples.push_back(x);
+            vmean.push_back(sum / 1000);
+            vvariance.push_back(sum2 / (1000) - (sum / (1000)) * (sum / (1000)));
+            pmean.push_back(pedestal.mean(0, 0));
+            pvariance.push_back(pedestal.variance(0, 0));
+        }
+        if (x % 1000 == 999) {
+            MEAN *= 1.1;
+            STD *= 1.1;
+
+            distribution.param(std::normal_distribution<double>::param_type(MEAN, STD));
+            cur_mean.push_back({x, MEAN});
+            cur_variance.push_back({x, STD * STD});
+        }
+    }
+
+    logger::info("x6=", samples);
+    logger::info("pmean6=", pmean);
+    logger::info("pvar6=", pvariance);
+    logger::info("vmean6=", vmean);
+    logger::info("vvar6=", vvariance);
+    std::cout << "cur_mean6=";
+    print_vpair(cur_mean);
+    std::cout << "cur_variance6=";
+    print_vpair(cur_variance);
+
+    // std::cout << "PEDESTAL" << std::endl;
+    // std::cout << "Mean: " << pmean(0, 0) << std::endl;
+    // std::cout << "Variance: " << pvariance(0, 0) << std::endl;
+    // std::cout << "Standard Deviation: " << pstandard_deviation(0, 0) << std::endl;
+
+    // std::cout << "VALUES" << std::endl;
+    // std::cout << "Mean: " << std::accumulate(values.begin(), values.end(), 0.0) / values.size() << std::endl;
+    // std::cout << "Variance: " << variance<double>(values) << std::endl;
+    // std::cout << "Standard Deviation: " << std::sqrt(variance<double>(values)) << std::endl;
+}

examples/zmq_multi_receiver.cpp

@@ -1,7 +1,5 @@
+#include "aare.hpp"
 #include "aare/examples/defs.hpp"
-#include "aare/network_io/ZmqMultiReceiver.hpp"
-#include "aare/network_io/ZmqSocketReceiver.hpp"
-#include "aare/network_io/defs.hpp"
 
 #include <boost/program_options.hpp>
 #include <cassert>

examples/zmq_receiver_example.cpp

@@ -1,6 +1,5 @@
+#include "aare.hpp"
 #include "aare/examples/defs.hpp"
-#include "aare/network_io/ZmqSocketReceiver.hpp"
-#include "aare/network_io/defs.hpp"
 
 #include <boost/program_options.hpp>
 #include <cassert>

examples/zmq_restream_example.cpp

@@ -1,11 +1,9 @@
+#include "aare.hpp"
 #include "aare/examples/defs.hpp"
-#include <chrono>
-#include <thread>
-
-#include "aare/file_io/File.hpp"
-#include "aare/network_io/ZmqSocketSender.hpp"
 
 #include <boost/program_options.hpp>
+#include <chrono>
+#include <thread>
 
 using namespace aare;
 using namespace std;

examples/zmq_sender_example.cpp

@@ -1,9 +1,5 @@
-#include "aare/core/Frame.hpp"
+#include "aare.hpp"
 #include "aare/examples/defs.hpp"
-#include "aare/network_io/ZmqHeader.hpp"
-#include "aare/network_io/ZmqSocketSender.hpp"
-#include "aare/network_io/defs.hpp"
-#include "aare/utils/logger.hpp"
 
 #include <ctime> // std::time
 #include <fmt/core.h>

examples/zmq_sink.cpp

@@ -1,8 +1,5 @@
-
-#include "aare/network_io/ZmqSink.hpp"
-#include "aare/network_io/ZmqSocketReceiver.hpp"
-#include "aare/network_io/ZmqVentilator.hpp"
-#include "aare/network_io/ZmqWorker.hpp"
+#include "aare.hpp"
+#include "aare/examples/defs.hpp"
 
 #include "zmq.h"
 #include <boost/program_options.hpp>

examples/zmq_task_ventilator.cpp

@@ -1,5 +1,5 @@
-#include "aare/network_io/ZmqSocketReceiver.hpp"
-#include "aare/network_io/ZmqVentilator.hpp"
+#include "aare.hpp"
+#include "aare/examples/defs.hpp"
 
 #include "zmq.h"
 #include <boost/program_options.hpp>

examples/zmq_worker.cpp

@@ -1,8 +1,6 @@
 
-#include "aare/network_io/ZmqSink.hpp"
-#include "aare/network_io/ZmqSocketReceiver.hpp"
-#include "aare/network_io/ZmqVentilator.hpp"
-#include "aare/network_io/ZmqWorker.hpp"
+#include "aare.hpp"
+#include "aare/examples/defs.hpp"
 
 #include "zmq.h"
 #include <boost/program_options.hpp>

include/aare.hpp

@@ -4,4 +4,5 @@
 #include "aare/core.hpp"
 #include "aare/file_io.hpp"
 #include "aare/network_io.hpp"
+#include "aare/processing.hpp"
 #include "aare/utils.hpp"

include/aare/file_io.hpp

@@ -1,4 +1,3 @@
-// ClusterFile.hpp  File.hpp  FileInterface.hpp  NumpyFile.hpp  NumpyHelpers.hpp  RawFile.hpp  SubFile.hpp
 #include "aare/file_io/ClusterFile.hpp"
 #include "aare/file_io/File.hpp"
 #include "aare/file_io/FileInterface.hpp"

include/aare/network_io.hpp

@@ -1,4 +1,5 @@
 #include "aare/network_io/ZmqHeader.hpp"
+#include "aare/network_io/ZmqMultiReceiver.hpp"
 #include "aare/network_io/ZmqSink.hpp"
 #include "aare/network_io/ZmqSocket.hpp"
 #include "aare/network_io/ZmqSocketReceiver.hpp"

include/aare/processing.hpp

@@ -0,0 +1 @@
+#include "aare/processing/Pedestal.hpp"

src/CMakeLists.txt

@@ -0,0 +1,8 @@
+add_subdirectory(core)
+add_subdirectory(file_io)
+add_subdirectory(utils)
+add_subdirectory(network_io)
+add_subdirectory(processing)
+if(AARE_PYTHON_BINDINGS)
+    add_subdirectory(python)
+endif()

src/core/include/aare/core/Frame.hpp

@@ -79,10 +79,10 @@ class Frame {
         std::memcpy(m_data, other.m_data, m_rows * m_cols * m_bitdepth / 8);
     }
 
-    template <typename T> NDView<T> view() {
-        std::vector<ssize_t> shape = {static_cast<ssize_t>(m_rows), static_cast<ssize_t>(m_cols)};
+    template <typename T> NDView<T, 2> view() {
+        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>(data, shape);
+        return NDView<T, 2>(data, shape);
     }
 
     template <typename T> NDArray<T> image() { return NDArray<T>(this->view<T>()); }

src/core/include/aare/core/NDView.hpp

@@ -54,9 +54,9 @@ template <typename T, ssize_t Ndim = 2> class NDView {
         : 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<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<>())) {}
+    // 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<>())) {}
 
     template <typename... Ix> std::enable_if_t<sizeof...(Ix) == Ndim, T &> operator()(Ix... index) {
         return buffer_[element_offset(strides_, index...)];
@@ -70,8 +70,8 @@ template <typename T, ssize_t Ndim = 2> class NDView {
 
     T *begin() { return buffer_; }
     T *end() { return buffer_ + size_; }
-    T &operator()(ssize_t i) { return buffer_[i]; }
-    T &operator[](ssize_t i) { 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_)
@@ -142,6 +142,4 @@ template <typename T, ssize_t Ndim = 2> class NDView {
     }
 };
 
-template class NDView<uint16_t, 2>;
-
 } // namespace aare

src/core/include/aare/core/defs.hpp

@@ -52,7 +52,7 @@ template <typename T> struct t_xy {
     bool operator!=(const t_xy &other) const { return !(*this == other); }
     std::string to_string() const { return "{ x: " + std::to_string(row) + " y: " + std::to_string(col) + " }"; }
 };
-typedef t_xy<uint32_t> xy;
+using xy = t_xy<uint32_t>;
 
 using dynamic_shape = std::vector<ssize_t>;
 

src/core/test/NDView.test.cpp

@@ -142,14 +142,14 @@ TEST_CASE("iterators") {
     }
 }
 
-TEST_CASE("shape from vector") {
-    std::vector<int> vec;
-    for (int i = 0; i != 12; ++i) {
-        vec.push_back(i);
-    }
-    std::vector<ssize_t> shape{3, 4};
-    NDView<int, 2> data(vec.data(), shape);
-}
+// TEST_CASE("shape from vector") {
+//     std::vector<int> vec;
+//     for (int i = 0; i != 12; ++i) {
+//         vec.push_back(i);
+//     }
+//     std::vector<ssize_t> shape{3, 4};
+//     NDView<int, 2> data(vec.data(), shape);
+// }
 
 TEST_CASE("divide with another span") {
     std::vector<int> vec0{9, 12, 3};

src/core/test/wrappers.test.cpp

@@ -32,7 +32,7 @@ TEST_CASE("NDView") {
         data[i] = i;
     }
     SECTION("constructors") {
-        NDView<uint16_t, 2> ds(data, std::vector<ssize_t>({10, 10}));
+        NDView<uint16_t, 2> ds(data, std::array<ssize_t, 2>({10, 10}));
         for (int i = 0; i < 100; i++) {
             REQUIRE(ds(i / 10, i % 10) == data[i]);
         }

src/network_io/include/aare/network_io/ZmqHeader.hpp

@@ -144,7 +144,7 @@ template <typename T, int N, typename SIMDJSON_VALUE> std::array<T, N> simd_conv
     std::array<T, N> arr{};
     int i = 0;
     for (auto v : simd_array) {
-        int64_t tmp;
+        int64_t tmp{};
         err = v.get(tmp);
         if (err)
             throw std::runtime_error("error converting simdjson::ondemand::value");

src/network_io/include/aare/network_io/defs.hpp

@@ -28,13 +28,13 @@ struct Task {
     size_t data_size{};
 #pragma GCC diagnostic push
 #pragma GCC diagnostic ignored "-Wpedantic"
-    std::byte payload[];
+    std::byte payload[]; // NOLINT
 #pragma GCC diagnostic pop
 
     static const size_t MAX_DATA_SIZE = 1024 * 1024; // 1MB
     size_t size() const { return sizeof(Task) + data_size; }
     static Task *init(std::byte *data, size_t data_size) {
-        Task *task = (Task *)new std::byte[sizeof(Task) + data_size];
+        Task *task = reinterpret_cast<Task *>(new std::byte[sizeof(Task) + data_size]);
         task->data_size = data_size;
         if (data_size > 0)
             memcpy(task->payload, data, data_size);

src/processing/CMakeLists.txt

@@ -0,0 +1,17 @@
+
+add_library(processing STATIC src/Pedestal.cpp)
+target_include_directories(processing PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}/include)
+target_link_libraries(processing PUBLIC core)
+
+
+if(AARE_PYTHON_BINDINGS)
+set_property(TARGET processing PROPERTY POSITION_INDEPENDENT_CODE ON)
+endif()
+
+if(AARE_TESTS)
+    set(TestSources
+            ${CMAKE_CURRENT_SOURCE_DIR}/test/pedestal.test.cpp
+        )
+    target_sources(tests PRIVATE ${TestSources} )
+    target_link_libraries(tests PRIVATE processing)
+endif()

src/processing/include/aare/processing/Pedestal.hpp

@@ -0,0 +1,64 @@
+#pragma once
+#include "aare/core/Frame.hpp"
+#include "aare/core/NDArray.hpp"
+#include "aare/core/NDView.hpp"
+#include <cstddef>
+
+namespace aare {
+class Pedestal {
+  public:
+    Pedestal(int rows, int cols, int n_samples = 1000);
+    ~Pedestal();
+
+    // frame level operations
+    template <typename T> void push(NDView<T, 2> frame) {
+        assert(frame.size() == m_rows * m_cols);
+        // TODO: test the effect of #pragma omp parallel for
+        for (int index = 0; index < m_rows * m_cols; index++) {
+            push<T>(index % m_cols, index / m_cols, frame(index));
+        }
+    }
+    template <typename T> void push(Frame &frame) {
+        assert(frame.rows() == static_cast<size_t>(m_rows) && frame.cols() == static_cast<size_t>(m_cols));
+        push<T>(frame.view<T>());
+    }
+    NDArray<double> mean();
+    NDArray<double> variance();
+    NDArray<double> standard_deviation();
+    void clear();
+
+    // getter functions
+    inline int rows() const { return m_rows; }
+    inline int cols() const { return m_cols; }
+    inline int n_samples() const { return m_samples; }
+    inline uint32_t *cur_samples() const { return m_cur_samples; }
+    inline double *get_sum() const { return m_sum; }
+    inline double *get_sum2() const { return m_sum2; }
+
+    // pixel level operations (should be refactored to allow users to implement their own pixel level operations)
+    template <typename T> inline void push(const int row, const int col, const T val) {
+        const int idx = index(row, col);
+        if (m_cur_samples[idx] < m_samples) {
+            m_sum[idx] += val;
+            m_sum2[idx] += val * val;
+            m_cur_samples[idx]++;
+        } else {
+            m_sum[idx] += val - m_sum[idx] / m_cur_samples[idx];
+            m_sum2[idx] += val * val - m_sum2[idx] / m_cur_samples[idx];
+        }
+    }
+    double mean(const int row, const int col) const;
+    double variance(const int row, const int col) const;
+    double standard_deviation(const int row, const int col) const;
+    inline int index(const int row, const int col) const { return row * m_cols + col; };
+    void clear(const int row, const int col);
+
+  private:
+    int m_rows;
+    int m_cols;
+    int m_samples;
+    uint32_t *m_cur_samples{nullptr};
+    double *m_sum{nullptr};
+    double *m_sum2{nullptr};
+};
+} // namespace aare

src/processing/src/Pedestal.cpp

@@ -0,0 +1,70 @@
+#include "aare/processing/Pedestal.hpp"
+#include <cmath>
+#include <cstddef>
+
+namespace aare {
+Pedestal::Pedestal(int rows, int cols, int n_samples)
+    : m_rows(rows), m_cols(cols), m_samples(n_samples), m_sum(new double[rows * cols]{}),
+      m_sum2(new double[rows * cols]{}), m_cur_samples(new uint32_t[rows * cols]{}) {
+    assert(rows > 0 && cols > 0 && n_samples > 0);
+}
+
+NDArray<double, 2> Pedestal::mean() {
+    NDArray<double, 2> mean_array({m_rows, m_cols});
+    for (int 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;
+}
+
+NDArray<double, 2> Pedestal::variance() {
+    NDArray<double, 2> variance_array({m_rows, m_cols});
+    for (int i = 0; i < m_rows * m_cols; i++) {
+        variance_array(i % m_cols, i / m_cols) = variance(i % m_cols, i / m_cols);
+    }
+    return variance_array;
+}
+
+NDArray<double, 2> Pedestal::standard_deviation() {
+    NDArray<double, 2> standard_deviation_array({m_rows, m_cols});
+    for (int i = 0; i < m_rows * m_cols; i++) {
+        standard_deviation_array(i % m_cols, i / m_cols) = standard_deviation(i % m_cols, i / m_cols);
+    }
+    return standard_deviation_array;
+}
+void Pedestal::clear() {
+    for (int i = 0; i < m_rows * m_cols; i++) {
+        clear(i % m_cols, i / m_cols);
+    }
+}
+
+/*
+ * index level operations
+ */
+double Pedestal::mean(const int row, const int col) const {
+    if (m_cur_samples[index(row, col)] == 0) {
+        return 0.0;
+    }
+    return m_sum[index(row, col)] / m_cur_samples[index(row, col)];
+}
+double Pedestal::variance(const int row, const int col) const {
+    if (m_cur_samples[index(row, col)] == 0) {
+        return 0.0;
+    }
+    return m_sum2[index(row, col)] / m_cur_samples[index(row, col)] - mean(row, col) * mean(row, col);
+}
+double Pedestal::standard_deviation(const int row, const int col) const { return std::sqrt(variance(row, col)); }
+
+void Pedestal::clear(const int row, const int col) {
+    m_sum[index(row, col)] = 0;
+    m_sum2[index(row, col)] = 0;
+    m_cur_samples[index(row, col)] = 0;
+}
+
+Pedestal::~Pedestal() {
+    delete[] m_sum;
+    delete[] m_sum2;
+    delete[] m_cur_samples;
+}
+
+} // namespace aare

src/processing/test/pedestal.test.cpp

@@ -0,0 +1,106 @@
+#include "aare/processing/Pedestal.hpp"
+#include "aare/utils/floats.hpp"
+#include <catch2/catch_test_macros.hpp>
+#include <chrono>
+#include <random>
+
+using namespace aare;
+TEST_CASE("test pedestal constructor") {
+    aare::Pedestal pedestal(10, 10, 5);
+    REQUIRE(pedestal.rows() == 10);
+    REQUIRE(pedestal.cols() == 10);
+    REQUIRE(pedestal.n_samples() == 5);
+    REQUIRE(pedestal.cur_samples() != nullptr);
+    REQUIRE(pedestal.get_sum() != nullptr);
+    REQUIRE(pedestal.get_sum2() != nullptr);
+    for (int i = 0; i < 10; i++) {
+        for (int j = 0; j < 10; j++) {
+            REQUIRE(pedestal.get_sum()[pedestal.index(i, j)] == 0);
+            REQUIRE(pedestal.get_sum2()[pedestal.index(i, j)] == 0);
+            REQUIRE(pedestal.cur_samples()[pedestal.index(i, j)] == 0);
+        }
+    }
+}
+
+TEST_CASE("test pedestal push") {
+    aare::Pedestal pedestal(10, 10, 5);
+    aare::Frame frame(10, 10, 16);
+    for (int i = 0; i < 10; i++) {
+        for (int j = 0; j < 10; j++) {
+            frame.set<uint16_t>(i, j, i + j);
+        }
+    }
+
+    // test single push
+    pedestal.push<uint16_t>(frame);
+    for (int i = 0; i < 10; i++) {
+        for (int j = 0; j < 10; j++) {
+            REQUIRE(pedestal.get_sum()[pedestal.index(i, j)] == i + j);
+            REQUIRE(pedestal.get_sum2()[pedestal.index(i, j)] == (i + j) * (i + j));
+            REQUIRE(pedestal.cur_samples()[pedestal.index(i, j)] == 1);
+        }
+    }
+
+    // test clear
+    pedestal.clear();
+    for (int i = 0; i < 10; i++) {
+        for (int j = 0; j < 10; j++) {
+            REQUIRE(pedestal.get_sum()[pedestal.index(i, j)] == 0);
+            REQUIRE(pedestal.get_sum2()[pedestal.index(i, j)] == 0);
+            REQUIRE(pedestal.cur_samples()[pedestal.index(i, j)] == 0);
+        }
+    }
+
+    // test number of samples after multiple push
+    for (int k = 0; k < 50; k++) {
+        pedestal.push<uint16_t>(frame);
+        for (int i = 0; i < 10; i++) {
+            for (int j = 0; j < 10; j++) {
+                if (k < 5) {
+                    REQUIRE(pedestal.cur_samples()[pedestal.index(i, j)] == k + 1);
+                    REQUIRE(pedestal.get_sum()[pedestal.index(i, j)] == (k + 1) * (i + j));
+                    REQUIRE(pedestal.get_sum2()[pedestal.index(i, j)] == (k + 1) * (i + j) * (i + j));
+                } else {
+                    REQUIRE(pedestal.cur_samples()[pedestal.index(i, j)] == 5);
+                    REQUIRE(pedestal.get_sum()[pedestal.index(i, j)] == 5 * (i + j));
+                    REQUIRE(pedestal.get_sum2()[pedestal.index(i, j)] == 5 * (i + j) * (i + j));
+                }
+                REQUIRE(pedestal.mean(i, j) == (i + j));
+                REQUIRE(pedestal.variance(i, j) == 0);
+                REQUIRE(pedestal.standard_deviation(i, j) == 0);
+            }
+        }
+    }
+}
+
+TEST_CASE("test pedestal with normal distribution") {
+    const double MEAN = 5.0, STD = 2.0, VAR = STD * STD, TOLERANCE = 0.1;
+
+    unsigned seed = std::chrono::system_clock::now().time_since_epoch().count();
+    std::default_random_engine generator(seed);
+    std::normal_distribution<double> distribution(MEAN, STD);
+
+    aare::Pedestal pedestal(4, 4, 10000);
+    for (int i = 0; i < 10000; i++) {
+        aare::Frame frame(4, 4, 64);
+        for (int j = 0; j < 4; j++) {
+            for (int k = 0; k < 4; k++) {
+                frame.set<double>(j, k, distribution(generator));
+            }
+        }
+        pedestal.push<double>(frame);
+    }
+
+    auto mean = pedestal.mean();
+    auto variance = pedestal.variance();
+    auto standard_deviation = pedestal.standard_deviation();
+
+    for (int i = 0; i < 4; i++) {
+        for (int j = 0; j < 4; j++) {
+            // 10% tolerance
+            compare_floats<double>(mean(i, j), MEAN, MEAN * TOLERANCE);
+            compare_floats<double>(variance(i, j), VAR, VAR * TOLERANCE);
+            compare_floats<double>(standard_deviation(i, j), STD, STD * TOLERANCE);
+        }
+    }
+}

src/utils/CMakeLists.txt

@@ -8,7 +8,9 @@ target_link_libraries(utils PUBLIC core)
 if(AARE_TESTS)
     set(TestSources
             ${CMAKE_CURRENT_SOURCE_DIR}/test/merge_frames.test.cpp
+            ${CMAKE_CURRENT_SOURCE_DIR}/test/floats.test.cpp
         )
     target_sources(tests PRIVATE ${TestSources} )
     target_link_libraries(tests PRIVATE utils)
+
 endif()

src/utils/include/aare/utils/compare_files.hpp

@@ -13,7 +13,7 @@ namespace aare {
  * @param p2 path to the second file
  * @return true if the files are the same, false otherwise
  */
-bool compare_files(const std::string &p1, const std::string &p2) {
+bool inline compare_files(const std::string &p1, const std::string &p2) {
     std::ifstream f1(p1, std::ifstream::binary | std::ifstream::ate);
     std::ifstream f2(p2, std::ifstream::binary | std::ifstream::ate);
 
@@ -31,7 +31,7 @@ bool compare_files(const std::string &p1, const std::string &p2) {
     return std::equal(std::istreambuf_iterator<char>(f1.rdbuf()), std::istreambuf_iterator<char>(),
                       std::istreambuf_iterator<char>(f2.rdbuf()));
 }
-bool compare_files(const std::filesystem::path &p1, const std::filesystem::path &p2) {
+bool inline compare_files(const std::filesystem::path &p1, const std::filesystem::path &p2) {
     return compare_files(p1.string(), p2.string());
 }
 } // namespace aare

src/utils/include/aare/utils/floats.hpp

@@ -0,0 +1,25 @@
+#include <cfloat>
+#include <cmath>
+#include <iostream>
+namespace aare {
+template <typename T = float> bool compare_floats(T A, T B, T maxRelDiff = -1) {
+    bool user_defined_diff = true;
+    if (maxRelDiff < 0) {
+        user_defined_diff = false;
+        maxRelDiff = std::numeric_limits<T>::epsilon();
+    }
+    // Calculate the difference.
+    T diff = fabs(A - B);
+    A = fabs(A);
+    B = fabs(B);
+    // Find the largest
+    T largest = (B > A) ? B : A;
+
+    // if user defined maxRelDiff then compare to it without scaling
+    if (user_defined_diff) {
+        return diff <= maxRelDiff;
+    }
+    return diff <= largest * maxRelDiff;
+}
+
+} // namespace aare