Merge branch 'main' into dev/strixels/remap_simple

.gitignore

@@ -127,6 +127,8 @@ cover/
 *.mo
 *.pot
 
+.*
+
 # Django stuff:
 *.log
 local_settings.py
@@ -192,4 +194,4 @@ cython_debug/
 .ruff_cache/
 
 # user defined
-wheelhouse/
+wheelhouse/

CMakeLists.txt

@@ -62,7 +62,6 @@ option(AARE_FETCH_PYBIND11 "Use FetchContent to download pybind11" ON)
 option(AARE_FETCH_CATCH "Use FetchContent to download catch2" ON)
 option(AARE_FETCH_JSON "Use FetchContent to download nlohmann::json" ON)
 option(AARE_FETCH_ZMQ "Use FetchContent to download libzmq" ON)
-option(AARE_FETCH_LMFIT "Use FetchContent to download lmfit" ON)
 option(AARE_FETCH_MINUIT2 "Use FetchContent to download Minuit2" ON)
 
 option(AARE_WARNINGS_AS_ERRORS "Treat warnings as errors during compilation"
@@ -87,8 +86,8 @@ if(AARE_SYSTEM_LIBRARIES)
   set(AARE_FETCH_ZMQ
       OFF
       CACHE BOOL "Disabled FetchContent for libzmq" FORCE)
-  # Still fetch lmfit and Minuit2 when setting AARE_SYSTEM_LIBRARIES since these
-  # are not available on conda-forge
+  # Still fetch Minuit2 when setting AARE_SYSTEM_LIBRARIES since it is not
+  # available on conda-forge
 endif()
 
 if(AARE_BENCHMARKS)
@@ -97,61 +96,6 @@ endif()
 
 set(CMAKE_EXPORT_COMPILE_COMMANDS ON)
 
-if(AARE_FETCH_LMFIT)
-  # TODO! Should we fetch lmfit from the web or inlcude a tar.gz in the repo?
-  set(LMFIT_PATCH_COMMAND git apply
-                          ${CMAKE_CURRENT_SOURCE_DIR}/patches/lmfit.patch)
-
-  # For cmake < 3.28 we can't supply EXCLUDE_FROM_ALL to FetchContent_Declare so
-  # we need this workaround
-  if(${CMAKE_VERSION} VERSION_LESS "3.28")
-    FetchContent_Declare(
-      lmfit
-      GIT_REPOSITORY https://jugit.fz-juelich.de/mlz/lmfit.git
-      GIT_TAG main
-      PATCH_COMMAND ${LMFIT_PATCH_COMMAND}
-      UPDATE_DISCONNECTED 1)
-  else()
-    FetchContent_Declare(
-      lmfit
-      GIT_REPOSITORY https://jugit.fz-juelich.de/mlz/lmfit.git
-      GIT_TAG main
-      PATCH_COMMAND ${LMFIT_PATCH_COMMAND}
-      UPDATE_DISCONNECTED 1
-      EXCLUDE_FROM_ALL 1)
-  endif()
-
-  # Disable what we don't need from lmfit
-  set(BUILD_TESTING
-      OFF
-      CACHE BOOL "")
-  set(LMFIT_CPPTEST
-      OFF
-      CACHE BOOL "")
-  set(LIB_MAN
-      OFF
-      CACHE BOOL "")
-  set(LMFIT_CPPTEST
-      OFF
-      CACHE BOOL "")
-  set(BUILD_SHARED_LIBS
-      OFF
-      CACHE BOOL "")
-
-  if(${CMAKE_VERSION} VERSION_LESS "3.28")
-    if(NOT lmfit_POPULATED)
-      FetchContent_Populate(lmfit)
-      add_subdirectory(${lmfit_SOURCE_DIR} ${lmfit_BINARY_DIR} EXCLUDE_FROM_ALL)
-    endif()
-  else()
-    FetchContent_MakeAvailable(lmfit)
-  endif()
-
-  set_property(TARGET lmfit PROPERTY POSITION_INDEPENDENT_CODE ON)
-else()
-  find_package(lmfit REQUIRED)
-endif()
-
 if(AARE_FETCH_MINUIT2)
   # We are building Minuit2 from sources.
 
@@ -396,6 +340,7 @@ set(PUBLICHEADERS
     include/aare/FileInterface.hpp
     include/aare/FilePtr.hpp
     include/aare/Frame.hpp
+    include/aare/MultiThreadedFileReader.hpp
     include/aare/hist/PixelHistogram.hpp
     include/aare/hist/PixelHistogramImpl.hpp
     include/aare/hist/PedestalTrackingPixelHistogram.hpp
@@ -420,7 +365,9 @@ set(PUBLICHEADERS
     include/aare/RawMasterFile.hpp
     include/aare/RawSubFile.hpp
     include/aare/VarClusterFinder.hpp
-    include/aare/utils/task.hpp)
+    include/aare/utils/task.hpp
+    include/aare/utils/ifstream_helpers.hpp
+    include/aare/utils/math_helpers.hpp)
 
 set(SourceFiles
     ${CMAKE_CURRENT_SOURCE_DIR}/src/calibration.cpp
@@ -437,6 +384,7 @@ set(SourceFiles
     ${CMAKE_CURRENT_SOURCE_DIR}/src/Frame.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/src/Interpolator.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/src/JungfrauDataFile.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/src/MultiThreadedFileReader.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/src/NumpyFile.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/src/NumpyHelpers.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/src/PixelMap.cpp
@@ -461,8 +409,7 @@ target_link_libraries(
   aare_core
   PUBLIC fmt::fmt nlohmann_json::nlohmann_json ${STD_FS_LIB} # from
                                                              # helpers.cmake
-  PRIVATE aare_compiler_flags Threads::Threads $<BUILD_INTERFACE:lmfit>
-          $<BUILD_INTERFACE:aare::Minuit2>)
+  PRIVATE aare_compiler_flags Threads::Threads $<BUILD_INTERFACE:aare::Minuit2>)
 
 target_include_directories(
   aare_core SYSTEM
@@ -494,6 +441,7 @@ if(AARE_TESTS)
       ${CMAKE_CURRENT_SOURCE_DIR}/src/defs.test.cpp
       ${CMAKE_CURRENT_SOURCE_DIR}/src/decode.test.cpp
       ${CMAKE_CURRENT_SOURCE_DIR}/src/Dtype.test.cpp
+      ${CMAKE_CURRENT_SOURCE_DIR}/src/Fit.test.cpp
       ${CMAKE_CURRENT_SOURCE_DIR}/src/Frame.test.cpp
       ${CMAKE_CURRENT_SOURCE_DIR}/src/DetectorGeometry.test.cpp
       ${CMAKE_CURRENT_SOURCE_DIR}/src/Interpolation.test.cpp
@@ -510,6 +458,7 @@ if(AARE_TESTS)
       ${CMAKE_CURRENT_SOURCE_DIR}/src/hist/PixelHistogramImpl.test.cpp
       ${CMAKE_CURRENT_SOURCE_DIR}/src/hist/PixelHistogram.test.cpp
       ${CMAKE_CURRENT_SOURCE_DIR}/src/JungfrauDataFile.test.cpp
+      ${CMAKE_CURRENT_SOURCE_DIR}/src/MultiThreadedFileReader.test.cpp
       ${CMAKE_CURRENT_SOURCE_DIR}/src/NumpyFile.test.cpp
       ${CMAKE_CURRENT_SOURCE_DIR}/src/NumpyHelpers.test.cpp
       ${CMAKE_CURRENT_SOURCE_DIR}/src/RawFile.test.cpp

RELEASE.md

@@ -1,5 +1,17 @@
 # Release notes
 
+## Next
+
+### API Changes:
+
+- Removed the lmfit dependency and the legacy ``fit_gaus``, ``fit_pol1``,
+  ``fit_scurve``, and ``fit_scurve2`` APIs. Use ``Gaussian``, ``Pol1``,
+  ``RisingScurve``, or ``FallingScurve`` and call ``model.fit(...)`` (or
+  ``fit(model, ...)``) instead.
+- Removed the legacy ``gaus``, ``pol1``, ``scurve``, and ``scurve2`` function
+  evaluators. Model objects are callable and provide the replacement, for
+  example ``Gaussian()(x, par)``.
+
 ## 2026.7.2
 
 
@@ -151,4 +163,3 @@ dhanya.thattil@psi.ch
 
 
 
-

benchmarks/fit_benchmark.cpp

@@ -75,24 +75,7 @@ static void report_accuracy(benchmark::State &state, const TestCase &tc,
 // Benchmarks
 // ----------
 
-// 1. lmcurve
-static void BM_FitGausLm(benchmark::State &state) {
-    const auto &tc = get_test_cases()[state.range(0)];
-    auto data = generate_gaussian_data(tc);
-    auto xv = data.x.view();
-    auto yv = data.y.view();
-
-    aare::NDArray<double, 1> result;
-    for (auto _ : state) {
-        result = aare::fit_gaus(xv, yv);
-        benchmark::DoNotOptimize(result.data());
-    }
-
-    report_accuracy(state, tc, result);
-    state.SetLabel(tc.name);
-}
-
-// 2. Minuit2, analytic gradient (no Hesse)
+// Minuit2, analytic gradient (no Hesse)
 static void BM_FitGausMinuitGrad(benchmark::State &state) {
     const auto &tc = get_test_cases()[state.range(0)];
     auto data = generate_gaussian_data(tc);
@@ -115,7 +98,7 @@ static void BM_FitGausMinuitGrad(benchmark::State &state) {
     state.SetLabel(tc.name);
 }
 
-// 3. Minuit2, analytic gradient + Hesse
+// Minuit2, analytic gradient + Hesse
 static void BM_FitGausMinuitGradHesse(benchmark::State &state) {
     const auto &tc = get_test_cases()[state.range(0)];
     auto data = generate_gaussian_data(tc);
@@ -145,8 +128,6 @@ static void BM_FitGausMinuitGradHesse(benchmark::State &state) {
     state.SetLabel(tc.name);
 }
 
-BENCHMARK(BM_FitGausLm)->DenseRange(0, 5)->Unit(benchmark::kMicrosecond);
-
 BENCHMARK(BM_FitGausMinuitGrad)
     ->DenseRange(0, 5)
     ->Unit(benchmark::kMicrosecond);
@@ -155,4 +136,4 @@ BENCHMARK(BM_FitGausMinuitGradHesse)
     ->DenseRange(0, 5)
     ->Unit(benchmark::kMicrosecond);
 
-BENCHMARK_MAIN();
+BENCHMARK_MAIN();

docs/src/MultiThreadedFileReader.rst

@@ -0,0 +1,37 @@
+MultiThreadedFileReader
+=======================
+
+``MultiThreadedFileReader`` reads one chunk per worker on each call. Each
+worker owns an independent :cpp:class:`aare::File`, while all workers write
+into non-overlapping regions of the destination buffer. The output is ordered
+by frame index and successive calls advance through the file.
+
+.. code-block:: cpp
+
+   #include "aare/MultiThreadedFileReader.hpp"
+
+   // Four workers, chunks of 128 frames, and at most 10,000 frames.
+   aare::experimental::MultiThreadedFileReader reader(path, 4, 128, 10'000);
+   while (reader.remaining_frames() != 0) {
+       // Contains at most 4 * 128 frames.
+       auto batch = reader.read();
+       process(batch);
+   }
+
+Omit the final argument to read every frame in the source. An explicit value
+of zero requests an empty result. The low-level ``read_into`` overload avoids
+an allocation when the caller already owns a buffer of at least
+``reader.next_read_bytes()`` bytes. Use ``read_all()`` to read every frame
+remaining from the current position, and ``seek()`` to reposition the reader.
+Call ``close()`` to release all worker file handles early.
+
+.. note::
+
+   Multiple workers do not guarantee faster reads. Performance depends on the
+   storage device and file format, so the thread count and chunk size should be
+   benchmark-driven.
+
+.. doxygenclass:: aare::experimental::MultiThreadedFileReader
+   :members:
+   :undoc-members:
+   :private-members:

docs/src/Requirements.rst

@@ -14,7 +14,7 @@ Requirements
 To simplify deployment we build and statically link a few libraries.
 
 - fmt
-- lmfit - https://jugit.fz-juelich.de/mlz/lmfit
+- Minuit2
 - nlohmann_json
 - pybind11
 - ZeroMQ
@@ -23,4 +23,4 @@ To simplify deployment we build and statically link a few libraries.
 
 - Sphinx
 - Breathe
-- Doxygen
+- Doxygen

docs/src/index.rst

@@ -28,6 +28,7 @@ AARE
     pycalibration
     python/cluster/index
     python/file/index
+    python/experimental/index
     python/histogram/index
     pyFit
     
@@ -42,6 +43,7 @@ AARE
     NDView
     Frame
     File
+    MultiThreadedFileReader
     Dtype
     Cluster
     ClusterFinder

docs/src/pyFit.rst

@@ -1,19 +1,30 @@
-
-Fit
-========
+Fitting
+-------
 
 .. py:currentmodule:: aare
 
+Aare fits one-dimensional scans and three-dimensional pixel data with
+Minuit2. Create a model object and call its :meth:`fit` method::
 
-**Functions**
+    model = Gaussian(compute_errors=True)
+    result = model.fit(x, y, y_err)
 
-.. autofunction:: gaus
+The model object is also callable, which evaluates it at the supplied points::
 
-.. autofunction:: pol1
+    fitted_y = model(x, result["par"])
 
+The available models are ``Gaussian``, ``GaussianErfcPlateau``,
+``GaussianChargeSharing``, ``GaussianChargeSharingKb``, ``Pol1``, ``Pol2``,
+``RisingScurve``, and ``FallingScurve``. The module-level :func:`fit` function
+accepts the same model objects when a functional interface is preferred.
 
-**Fitting**
+For three-dimensional data, pass an array with shape
+``(rows, columns, scan_points)`` and select the worker count with
+``n_threads``::
 
-.. autofunction:: fit_gaus
+    result = model.fit(x, image_data, image_errors, n_threads=8)
 
-.. autofunction:: fit_pol1
+The result dictionary contains ``par`` and ``chi2``. It also contains
+``par_err`` when ``compute_errors`` is enabled.
+
+.. autofunction:: fit

docs/src/python/experimental/MultiThreadedFileReader.rst

@@ -0,0 +1,30 @@
+MultiThreadedFileReader
+=======================
+
+.. py:currentmodule:: aare.experimental
+
+The reader returns a NumPy array with shape ``(frames, rows, cols)`` and
+preserves the source pixel dtype. Each iteration reads at most
+``n_threads * chunk_size`` frames—one chunk per worker. File I/O runs with the
+Python GIL released.
+
+.. code-block:: python
+
+    from aare.experimental import MultiThreadedFileReader
+
+    with MultiThreadedFileReader(
+        "frames.npy", n_threads=4, chunk_size=128, total_frames=10_000
+    ) as reader:
+        for frames in reader:
+            process(frames)
+
+Call ``read()`` directly for the next batch, or ``read_all()`` for all frames
+remaining from the current position. ``tell()`` and ``seek()`` expose the
+iteration position. The context manager closes all worker files on exit.
+``close()`` is also available for explicit cleanup and may be called
+repeatedly.
+
+.. autoclass:: MultiThreadedFileReader
+    :members:
+    :undoc-members:
+    :show-inheritance:

docs/src/python/experimental/index.rst

@@ -0,0 +1,10 @@
+Experimental
+============
+
+APIs in this module may change without notice.
+
+.. toctree::
+   :caption: Experimental
+   :maxdepth: 1
+
+   MultiThreadedFileReader

include/aare/File.hpp

@@ -64,8 +64,9 @@ class File {
     size_t total_frames() const;
     size_t rows() const;
     size_t cols() const;
+    Dtype dtype() const;
 
     DetectorType detector_type() const;
 };
 
-} // namespace aare
+} // namespace aare

include/aare/FileInterface.hpp

@@ -145,10 +145,15 @@ class FileInterface {
      */
     virtual size_t bitdepth() const = 0;
 
-    virtual DetectorType detector_type() const = 0;
+    /**
+     * @brief get the data type of the pixels
+     * @return pixel data type
+     */
+    virtual Dtype dtype() const {
+        return Dtype::from_bitdepth(static_cast<uint8_t>(bitdepth()));
+    }
 
-    // function to query the data type of the file
-    /*virtual DataType dtype = 0; */
+    virtual DetectorType detector_type() const = 0;
 
     virtual ~FileInterface() = default;
 
@@ -168,4 +173,4 @@ class FileInterface {
     // size_t current_frame{};
 };
 
-} // namespace aare
+} // namespace aare

include/aare/Fit.hpp

@@ -1,111 +1,11 @@
 // SPDX-License-Identifier: MPL-2.0
 #pragma once
 
-#include <cmath>
-#include <vector>
-
 #include "aare/FitModel.hpp"
 #include "aare/NDArray.hpp"
-#include "aare/utils/par.hpp"
-#include "aare/utils/task.hpp"
 
 namespace aare {
 
-namespace func {
-double gaus(const double x, const double *par);
-NDArray<double, 1> gaus(NDView<double, 1> x, NDView<double, 1> par);
-
-double pol1(const double x, const double *par);
-NDArray<double, 1> pol1(NDView<double, 1> x, NDView<double, 1> par);
-
-double scurve(const double x, const double *par);
-NDArray<double, 1> scurve(NDView<double, 1> x, NDView<double, 1> par);
-
-double scurve2(const double x, const double *par);
-NDArray<double, 1> scurve2(NDView<double, 1> x, NDView<double, 1> par);
-
-} // namespace func
-
-static constexpr int DEFAULT_NUM_THREADS = 4;
-
-/**
- * @brief Fit a 1D Gaussian to data.
- * @param data data to fit
- * @param x x values
- */
-NDArray<double, 1> fit_gaus(NDView<double, 1> x, NDView<double, 1> y);
-
-/**
- * @brief Fit a 1D Gaussian to each pixel. Data layout [row, col, values]
- * @param x x values
- * @param y y values, layout [row, col, values]
- * @param n_threads number of threads to use
- */
-NDArray<double, 3> fit_gaus(NDView<double, 1> x, NDView<double, 3> y,
-                            int n_threads = DEFAULT_NUM_THREADS);
-
-/**
- * @brief Fit a 1D Gaussian with error estimates
- * @param x x values
- * @param y y values, layout [row, col, values]
- * @param y_err error in y, layout [row, col, values]
- * @param par_out output parameters
- * @param par_err_out output error parameters
- */
-void fit_gaus(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
-              NDView<double, 1> par_out, NDView<double, 1> par_err_out,
-              double &chi2);
-
-/**
- * @brief Fit a 1D Gaussian to each pixel with error estimates. Data layout
- * [row, col, values]
- * @param x x values
- * @param y y values, layout [row, col, values]
- * @param y_err error in y, layout [row, col, values]
- * @param par_out output parameters, layout [row, col, values]
- * @param par_err_out output parameter errors, layout [row, col, values]
- * @param n_threads number of threads to use
- */
-void fit_gaus(NDView<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_err,
-              NDView<double, 3> par_out, NDView<double, 3> par_err_out,
-              NDView<double, 2> chi2_out, int n_threads = DEFAULT_NUM_THREADS);
-
-NDArray<double, 1> fit_pol1(NDView<double, 1> x, NDView<double, 1> y);
-
-NDArray<double, 3> fit_pol1(NDView<double, 1> x, NDView<double, 3> y,
-                            int n_threads = DEFAULT_NUM_THREADS);
-
-void fit_pol1(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
-              NDView<double, 1> par_out, NDView<double, 1> par_err_out,
-              double &chi2);
-
-// TODO! not sure we need to offer the different version in C++
-void fit_pol1(NDView<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_err,
-              NDView<double, 3> par_out, NDView<double, 3> par_err_out,
-              NDView<double, 2> chi2_out, int n_threads = DEFAULT_NUM_THREADS);
-
-NDArray<double, 1> fit_scurve(NDView<double, 1> x, NDView<double, 1> y);
-NDArray<double, 3> fit_scurve(NDView<double, 1> x, NDView<double, 3> y,
-                              int n_threads);
-void fit_scurve(NDView<double, 1> x, NDView<double, 1> y,
-                NDView<double, 1> y_err, NDView<double, 1> par_out,
-                NDView<double, 1> par_err_out, double &chi2);
-void fit_scurve(NDView<double, 1> x, NDView<double, 3> y,
-                NDView<double, 3> y_err, NDView<double, 3> par_out,
-                NDView<double, 3> par_err_out, NDView<double, 2> chi2_out,
-                int n_threads);
-
-NDArray<double, 1> fit_scurve2(NDView<double, 1> x, NDView<double, 1> y);
-NDArray<double, 3> fit_scurve2(NDView<double, 1> x, NDView<double, 3> y,
-                               int n_threads);
-void fit_scurve2(NDView<double, 1> x, NDView<double, 1> y,
-                 NDView<double, 1> y_err, NDView<double, 1> par_out,
-                 NDView<double, 1> par_err_out, double &chi2);
-void fit_scurve2(NDView<double, 1> x, NDView<double, 3> y,
-                 NDView<double, 3> y_err, NDView<double, 3> par_out,
-                 NDView<double, 3> par_err_out, NDView<double, 2> chi2_out,
-                 int n_threads);
-
 // ---------------------------------------------------------------------------
 // Minuit2-based pixel fitting.
 // Template bodies and explicit instantiations live in src/Fit.cpp.

include/aare/JungfrauDataFile.hpp

@@ -49,6 +49,7 @@ class JungfrauDataFile : public FileInterface {
     size_t pixels_per_frame() override;
     size_t bytes_per_pixel() const;
     size_t bitdepth() const override;
+    Dtype dtype() const override { return Dtype::UINT16; }
     void seek(size_t frame_index)
         override; //!< seek to the given frame index (note not byte offset)
     size_t tell() override; //!< get the frame index of the file pointer
@@ -113,4 +114,4 @@ class JungfrauDataFile : public FileInterface {
     std::filesystem::path fpath(size_t frame_index) const;
 };
 
-} // namespace aare
+} // namespace aare

include/aare/MultiThreadedFileReader.hpp

@@ -0,0 +1,104 @@
+// SPDX-License-Identifier: MPL-2.0
+#pragma once
+
+#include "aare/Dtype.hpp"
+#include "aare/File.hpp"
+
+#include <cstddef>
+#include <filesystem>
+#include <optional>
+#include <vector>
+
+namespace aare::experimental {
+
+/**
+ * @brief Read independent chunks of a file in parallel.
+ *
+ * Each worker opens its own File instance, so seeking and reading do not share
+ * mutable file state. Chunks are written directly to their position in the
+ * destination buffer and the resulting frame order is the same as in the file.
+ */
+class MultiThreadedFileReader {
+  public:
+    /**
+     * @param fname path accepted by File
+     * @param n_threads maximum number of worker threads
+     * @param chunk_size number of frames claimed by a worker at a time
+     * @param total_frames number of frames to read, or all frames when omitted
+     */
+    MultiThreadedFileReader(std::filesystem::path fname, size_t n_threads,
+                            size_t chunk_size,
+                            std::optional<size_t> total_frames = std::nullopt);
+
+    MultiThreadedFileReader(const MultiThreadedFileReader &) = delete;
+    MultiThreadedFileReader &
+    operator=(const MultiThreadedFileReader &) = delete;
+    MultiThreadedFileReader(MultiThreadedFileReader &&) noexcept = default;
+    MultiThreadedFileReader &
+    operator=(MultiThreadedFileReader &&) noexcept = default;
+
+    /**
+     * @brief Read one chunk per active worker into a caller-owned buffer.
+     *
+     * The buffer must hold at least next_read_bytes() bytes. The reader's
+     * position advances by the returned number of frames. At the end of the
+     * configured range this function returns zero and does not access the
+     * destination.
+     */
+    size_t read_into(std::byte *destination);
+
+    /** @brief Read the next wave of chunks into an owned byte buffer. */
+    std::vector<std::byte> read();
+
+    /** @brief Read every frame remaining from the current position. */
+    std::vector<std::byte> read_all();
+
+    /** @brief Set the next frame index to read. The end position is valid. */
+    void seek(size_t frame_index);
+
+    /** @brief Return the next frame index to read. */
+    size_t tell() const noexcept { return m_current_frame; }
+
+    /** @brief Close all worker files. Safe to call more than once. */
+    void close() noexcept { m_files.clear(); }
+
+    /** @brief Return whether the worker files are open. */
+    bool is_open() const noexcept { return !m_files.empty(); }
+
+    size_t n_threads() const noexcept { return m_n_threads; }
+    size_t chunk_size() const noexcept { return m_chunk_size; }
+    size_t total_frames() const noexcept { return m_total_frames; }
+    size_t source_total_frames() const noexcept {
+        return m_source_total_frames;
+    }
+    size_t rows() const noexcept { return m_rows; }
+    size_t cols() const noexcept { return m_cols; }
+    size_t bitdepth() const noexcept { return m_bitdepth; }
+    Dtype dtype() const noexcept { return m_dtype; }
+    size_t bytes_per_frame() const noexcept { return m_bytes_per_frame; }
+    size_t total_bytes() const noexcept { return m_total_bytes; }
+    size_t remaining_frames() const noexcept;
+    size_t next_read_frames() const noexcept;
+    size_t next_read_bytes() const noexcept {
+        return next_read_frames() * m_bytes_per_frame;
+    }
+
+  private:
+    std::filesystem::path m_fname;
+    size_t m_n_threads;
+    size_t m_chunk_size;
+    size_t m_total_frames;
+    size_t m_source_total_frames;
+    size_t m_rows;
+    size_t m_cols;
+    size_t m_bitdepth;
+    Dtype m_dtype;
+    size_t m_bytes_per_frame;
+    size_t m_total_bytes;
+    size_t m_current_frame;
+    std::vector<File> m_files;
+
+    void ensure_open() const;
+};
+
+} // namespace aare::experimental

include/aare/NumpyFile.hpp

@@ -61,7 +61,7 @@ class NumpyFile : public FileInterface {
      * @brief get the data type of the numpy file
      * @return DType
      */
-    Dtype dtype() const { return m_header.dtype; }
+    Dtype dtype() const override { return m_header.dtype; }
 
     /**
      * @brief get the shape of the numpy file
@@ -129,4 +129,4 @@ class NumpyFile : public FileInterface {
     void write_impl(void *data, uint64_t size);
 };
 
-} // namespace aare
+} // namespace aare

include/aare/RawFile.hpp

@@ -118,6 +118,7 @@ class RawFile : public FileInterface {
      */
     size_t cols(const size_t roi_index) const;
     size_t bitdepth() const override;
+    Dtype dtype() const override { return Dtype::from_bitdepth(bitdepth()); }
     size_t n_modules() const;
 
     /**
@@ -170,4 +171,4 @@ class RawFile : public FileInterface {
     void open_subfiles(const size_t roi_index);
 };
 
-} // namespace aare
+} // namespace aare

include/aare/defs.hpp

@@ -38,87 +38,6 @@ inline constexpr size_t bits_per_byte = 8;
 
 void assert_failed(const std::string &msg);
 
-class DynamicCluster {
-  public:
-    int cluster_sizeX;
-    int cluster_sizeY;
-    int16_t x;
-    int16_t y;
-    Dtype dt; // 4 bytes
-
-  private:
-    std::byte *m_data;
-
-  public:
-    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()]{};
-    }
-    DynamicCluster() : DynamicCluster(3, 3) {}
-    DynamicCluster(const DynamicCluster &other)
-        : 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());
-    }
-    DynamicCluster &operator=(const DynamicCluster &other) {
-        if (this == &other)
-            return *this;
-        this->~DynamicCluster();
-        new (this) DynamicCluster(other);
-        return *this;
-    }
-    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);
-    }
-    ~DynamicCluster() { delete[] m_data; }
-    template <typename T> T get(int idx) {
-        (sizeof(T) == dt.bytes())
-            ? 0
-            : throw std::invalid_argument("[ERROR] Type size mismatch");
-        return *reinterpret_cast<T *>(m_data + idx * dt.bytes());
-    }
-    template <typename T> auto set(int idx, T val) {
-        (sizeof(T) == dt.bytes())
-            ? 0
-            : throw std::invalid_argument("[ERROR] Type size mismatch");
-        return memcpy(m_data + idx * dt.bytes(), &val, dt.bytes());
-    }
-
-    template <typename T> std::string to_string() const {
-        (sizeof(T) == dt.bytes())
-            ? 0
-            : throw std::invalid_argument("[ERROR] Type size mismatch");
-        std::string s = "x: " + std::to_string(x) + " y: " + std::to_string(y) +
-                        "\nm_data: [";
-        for (int i = 0; i < cluster_sizeX * cluster_sizeY; i++) {
-            s += std::to_string(
-                     *reinterpret_cast<T *>(m_data + i * dt.bytes())) +
-                 " ";
-        }
-        s += "]";
-        return s;
-    }
-    /**
-     * @brief size of the cluster in bytes when saved to a file
-     */
-    size_t size() const { return cluster_sizeX * cluster_sizeY; }
-    size_t bytes() const { return cluster_sizeX * cluster_sizeY * dt.bytes(); }
-    auto begin() const { return m_data; }
-    auto end() const {
-        return m_data + cluster_sizeX * cluster_sizeY * dt.bytes();
-    }
-    std::byte *data() { return m_data; }
-};
-
 /**
  * @brief header contained in parts of frames
  */

include/aare/utils/math_helpers.hpp

@@ -0,0 +1,14 @@
+#pragma once
+#include <cstddef>
+
+namespace aare {
+
+/**
+ * @brief Compute the ceiling of the integer division of n by d.
+ * @param n The numerator.
+ * @param d The denominator.
+ * @return The ceiling of the integer division.
+ */
+constexpr size_t ceil_div(size_t n, size_t d) { return n / d + (n % d != 0); }
+
+} // namespace aare

patches/lmfit.patch

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

python/CMakeLists.txt

@@ -40,7 +40,7 @@ set(PYTHON_FILES
     aare/ClusterVector.py
     aare/Cluster.py
     aare/calibration.py
-    aare/func.py
+    aare/experimental.py
     aare/RawFile.py
     aare/transform.py
     aare/ScanParameters.py

python/aare/__init__.py

@@ -3,8 +3,14 @@
 from . import _aare
 
 from . import transform
+from . import experimental
 
-from ._aare import File, RawMasterFile, RawSubFile, JungfrauDataFile
+from ._aare import (
+    File,
+    JungfrauDataFile,
+    RawMasterFile,
+    RawSubFile,
+)
 from ._aare import Pedestal_d, Pedestal_f, ClusterFinder_Cluster3x3i, VarClusterFinder
 from ._aare import DetectorType, ReadoutMode 
 from ._aare import hitmap
@@ -20,7 +26,6 @@ from .Cluster import Cluster
 
 from ._aare import Gaussian, RisingScurve, FallingScurve, Pol1, Pol2, GaussianErfcPlateau, GaussianChargeSharing, GaussianChargeSharingKb
 from ._aare import fit
-from ._aare import fit_gaus, fit_pol1, fit_scurve, fit_scurve2
 from ._aare import Interpolator
 from ._aare import calculate_eta2, calculate_eta3, calculate_cross_eta3, calculate_full_eta2
 from ._aare import reduce_to_2x2, reduce_to_3x3
@@ -36,9 +41,6 @@ from .ScanParameters import ScanParameters
 from .utils import random_pixels, random_pixel, flat_list, add_colorbar, Timer
 
 
-#make functions available in the top level API
-from .func import *
-
 from .calibration import *
 from ._aare import apply_calibration, count_switching_pixels
 from ._aare import calculate_pedestal, calculate_pedestal_float, calculate_pedestal_g0, calculate_pedestal_g0_float

python/aare/experimental.py

@@ -0,0 +1,6 @@
+# SPDX-License-Identifier: MPL-2.0
+"""Experimental APIs that may change without notice."""
+
+from ._aare.experimental import MultiThreadedFileReader
+
+__all__ = ["MultiThreadedFileReader"]

python/aare/func.py

@@ -1,2 +0,0 @@
-# SPDX-License-Identifier: MPL-2.0
-from ._aare import gaus, pol1, scurve, scurve2

python/examples/fits.py

@@ -1,106 +1,78 @@
 # SPDX-License-Identifier: MPL-2.0
 import matplotlib.pyplot as plt
 import numpy as np
-import sys
-sys.path.insert(0, '/home/kferjaoui/sw/aare/build')
-from aare import fit_gaus, fit_pol1
-from aare import Gaussian, fit
-from aare import pol1
 
-textpm = f"±"  #
-textmu = f"μ"  #
-textsigma = f"σ"  #
+from aare import Gaussian, Pol1
 
+textpm = "±"
+textmu = "μ"
+textsigma = "σ"
 
 
 # ================================= Gauss fit =================================
-# Parameters
-mu = np.random.uniform(1, 100)  # Mean of Gaussian
-sigma = np.random.uniform(4, 20)  # Standard deviation
-num_points = 10000  # Number of points for smooth distribution
-# noise_sigma = 10
-
-# Generate Gaussian distribution
-data = np.random.normal(mu, sigma, num_points)
+mu = np.random.uniform(1, 100)
+sigma = np.random.uniform(4, 20)
+data = np.random.normal(mu, sigma, 10000)
 counts, edges = np.histogram(data, bins=100)
 
-x = 0.5 * (edges[:-1] + edges[1:])       # proper bin centers
+x = 0.5 * (edges[:-1] + edges[1:])
 y = counts.astype(np.float64)
-
-# Poisson noise
 yerr = np.sqrt(np.maximum(y, 1))
-# yerr = np.abs(np.random.normal(0, noise_sigma, len(x)))
 
-# Create subplot
 fig0, ax0 = plt.subplots(1, 1, num=0, figsize=(12, 8))
-
-# Add the errors as error bars in the step plot
 ax0.errorbar(x, y, yerr=yerr, fmt=". ", capsize=5)
 ax0.grid()
 
-# Fit with lmfit
-result_lm = fit_gaus(x, y, yerr)
-par_lm = result_lm["par"]
-err_lm = result_lm["par_err"]
-chi2_lm = result_lm["chi2"] 
-print("[lmfit] fit_gaus:            ", par_lm, err_lm, chi2_lm)
+gaussian = Gaussian(compute_errors=True)
+result = gaussian.fit(x, y, yerr)
+par = result["par"]
+err = result["par_err"]
+chi2 = result["chi2"]
+print(f"Gaussian.fit: par={par}, err={err}, chi2={chi2}")
 
-# Fit with Minuit2 + analytic gradient + Hesse errors
-gaussian = Gaussian()
-gaussian.compute_errors = True
-result_m2 = gaussian.fit(x, y, yerr)
-par_m2 = result_m2['par']
-err_m2 = result_m2['par_err']
-chi2_m2 = result_m2['chi2']
-print(f"[minuit2] gaussian.fit: par={par_m2}, err={err_m2}, chi2={chi2_m2}")
-
-x = np.linspace(x[0], x[-1], 1000)
-ax0.plot(x, gaussian(x, par_lm), marker="", label="fit_gaus")
-ax0.plot(x, gaussian(x, par_m2), marker="", linestyle=":", label="fit_gaus_minuit_grad")
+x_plot = np.linspace(x[0], x[-1], 1000)
+ax0.plot(x_plot, gaussian(x_plot, par), marker="", label="Gaussian.fit")
 ax0.legend()
-ax0.set(xlabel="x", ylabel="Counts", 
+ax0.set(
+    xlabel="x",
+    ylabel="Counts",
     title=(
-        f"fit_gaus:       A={par_lm[0]:0.2f}{textpm}{err_lm[0]:0.2f}  "
-        f"{textmu}={par_lm[1]:0.2f}{textpm}{err_lm[1]:0.2f}  "
-        f"{textsigma}={par_lm[2]:0.2f}{textpm}{err_lm[2]:0.2f}\n"
-        f"minuit_grad: A={par_m2[0]:0.2f}{textpm}{err_m2[0]:0.2f}  "
-        f"{textmu}={par_m2[1]:0.2f}{textpm}{err_m2[1]:0.2f}  "
-        f"{textsigma}={par_m2[2]:0.2f}{textpm}{err_m2[2]:0.2f}\n"
+        f"A={par[0]:0.2f}{textpm}{err[0]:0.2f}  "
+        f"{textmu}={par[1]:0.2f}{textpm}{err[1]:0.2f}  "
+        f"{textsigma}={par[2]:0.2f}{textpm}{err[2]:0.2f}\n"
         f"(truth: {textmu}={mu:0.2f}, {textsigma}={sigma:0.2f})"
     ),
 )
 fig0.tight_layout()
 
 
-
-# ================================= pol1 fit =================================
-# Parameters
+# ================================= Pol1 fit =================================
 n_points = 40
-
-# Generate random slope and intercept (origin)
-slope = np.random.uniform(-10, 10)  # Random slope between 0.5 and 2.0
-intercept = np.random.uniform(-10, 10)  # Random intercept between -10 and 10
-
-# Generate random x values
+slope = np.random.uniform(-10, 10)
+intercept = np.random.uniform(-10, 10)
 x_values = np.random.uniform(-10, 10, n_points)
-
-# Calculate y values based on the linear function y = mx + b + error
 errors = np.abs(np.random.normal(0, np.random.uniform(1, 5), n_points))
-var_points = np.random.normal(0, np.random.uniform(0.1, 2), n_points)
-y_values = slope * x_values + intercept + var_points
+y_values = slope * x_values + intercept + np.random.normal(0, 1, n_points)
 
 fig1, ax1 = plt.subplots(1, 1, num=1, figsize=(12, 8))
 ax1.errorbar(x_values, y_values, yerr=errors, fmt=". ", capsize=5)
-result_pol = fit_pol1(x_values, y_values, errors)
-par = result_pol["par"]
-err = result_pol["par_err"]
 
-x = np.linspace(np.min(x_values), np.max(x_values), 1000)
-ax1.plot(x, pol1(x, par), marker="")
-ax1.set(xlabel="x", ylabel="y", title=f"a = {par[0]:0.2f}{textpm}{err[0]:0.2f}\n"
-                                      f"b = {par[1]:0.2f}{textpm}{err[1]:0.2f}\n"
-                                      f"(init: {slope:0.2f}, {intercept:0.2f})")
+pol1 = Pol1(compute_errors=True)
+result = pol1.fit(x_values, y_values, errors)
+par = result["par"]
+err = result["par_err"]
+
+x_plot = np.linspace(np.min(x_values), np.max(x_values), 1000)
+ax1.plot(x_plot, pol1(x_plot, par), marker="")
+ax1.set(
+    xlabel="x",
+    ylabel="y",
+    title=(
+        f"intercept = {par[0]:0.2f}{textpm}{err[0]:0.2f}\n"
+        f"slope = {par[1]:0.2f}{textpm}{err[1]:0.2f}\n"
+        f"(truth: {intercept:0.2f}, {slope:0.2f})"
+    ),
+)
 fig1.tight_layout()
 
 plt.show()
-

python/src/bind_MultiThreadedFileReader.hpp

@@ -0,0 +1,175 @@
+// SPDX-License-Identifier: MPL-2.0
+#pragma once
+
+#include "aare/MultiThreadedFileReader.hpp"
+
+#include <cstddef>
+#include <cstdint>
+#include <filesystem>
+#include <optional>
+#include <stdexcept>
+#include <vector>
+
+#include <pybind11/numpy.h>
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+#include <pybind11/stl/filesystem.h>
+
+namespace py = pybind11;
+
+inline py::dtype multi_threaded_reader_numpy_dtype(const aare::Dtype &dtype) {
+    using aare::Dtype;
+    if (dtype == Dtype::INT8)
+        return py::dtype::of<int8_t>();
+    if (dtype == Dtype::UINT8)
+        return py::dtype::of<uint8_t>();
+    if (dtype == Dtype::INT16)
+        return py::dtype::of<int16_t>();
+    if (dtype == Dtype::UINT16)
+        return py::dtype::of<uint16_t>();
+    if (dtype == Dtype::INT32)
+        return py::dtype::of<int32_t>();
+    if (dtype == Dtype::UINT32)
+        return py::dtype::of<uint32_t>();
+    if (dtype == Dtype::INT64)
+        return py::dtype::of<int64_t>();
+    if (dtype == Dtype::UINT64)
+        return py::dtype::of<uint64_t>();
+    if (dtype == Dtype::FLOAT)
+        return py::dtype::of<float>();
+    if (dtype == Dtype::DOUBLE)
+        return py::dtype::of<double>();
+    throw std::runtime_error("Unsupported pixel data type");
+}
+
+inline py::array
+multi_threaded_reader_read(aare::experimental::MultiThreadedFileReader &reader,
+                           bool read_all) {
+    const size_t n_frames =
+        read_all ? reader.remaining_frames() : reader.next_read_frames();
+    const std::vector<py::ssize_t> shape{
+        static_cast<py::ssize_t>(n_frames),
+        static_cast<py::ssize_t>(reader.rows()),
+        static_cast<py::ssize_t>(reader.cols())};
+
+    py::array image(multi_threaded_reader_numpy_dtype(reader.dtype()), shape);
+    auto *destination = reinterpret_cast<std::byte *>(image.mutable_data());
+    {
+        py::gil_scoped_release release;
+        if (read_all) {
+            size_t offset = 0;
+            while (reader.remaining_frames() != 0) {
+                const size_t frames_read =
+                    reader.read_into(destination + offset);
+                offset += frames_read * reader.bytes_per_frame();
+            }
+        } else {
+            reader.read_into(destination);
+        }
+    }
+    return image;
+}
+
+inline void define_multi_threaded_file_reader_bindings(py::module_ &m) {
+    using aare::experimental::MultiThreadedFileReader;
+
+    auto reader =
+        py::class_<MultiThreadedFileReader>(m, "MultiThreadedFileReader");
+    reader.attr("__module__") = "aare.experimental";
+    reader
+        .def(py::init<std::filesystem::path, size_t, size_t,
+                      std::optional<size_t>>(),
+             py::arg("fname"), py::arg("n_threads"), py::arg("chunk_size"),
+             py::arg("total_frames") = py::none(),
+             R"doc(
+                 Read chunks of detector frames concurrently.
+
+                 Each worker opens an independent File. The returned array is
+                 ordered by frame index even though chunks are read in parallel.
+
+                 Args:
+                     fname: Path accepted by File.
+                     n_threads: Maximum number of worker threads.
+                     chunk_size: Number of frames read per claimed chunk.
+                     total_frames: Optional frame limit. None reads all frames.
+                 )doc")
+        .def(
+            "read",
+            [](MultiThreadedFileReader &self) {
+                return multi_threaded_reader_read(self, false);
+            },
+            R"doc(
+                 Read one chunk per active worker into a NumPy array.
+
+                 Returns:
+                     An array containing at most n_threads * chunk_size frames.
+                     An empty array is returned at the end of the configured
+                     frame range. The GIL is released while file data is read.
+                 )doc")
+        .def(
+            "read_all",
+            [](MultiThreadedFileReader &self) {
+                return multi_threaded_reader_read(self, true);
+            },
+            R"doc(
+                 Read all frames remaining from the current position.
+                 )doc")
+        .def_property_readonly("n_threads", &MultiThreadedFileReader::n_threads)
+        .def_property_readonly("chunk_size",
+                               &MultiThreadedFileReader::chunk_size)
+        .def_property_readonly("total_frames",
+                               &MultiThreadedFileReader::total_frames)
+        .def_property_readonly("source_total_frames",
+                               &MultiThreadedFileReader::source_total_frames)
+        .def_property_readonly("rows", &MultiThreadedFileReader::rows)
+        .def_property_readonly("cols", &MultiThreadedFileReader::cols)
+        .def_property_readonly("bitdepth", &MultiThreadedFileReader::bitdepth)
+        .def_property_readonly("dtype",
+                               [](const MultiThreadedFileReader &self) {
+                                   return multi_threaded_reader_numpy_dtype(
+                                       self.dtype());
+                               })
+        .def_property_readonly("bytes_per_frame",
+                               &MultiThreadedFileReader::bytes_per_frame)
+        .def_property_readonly("total_bytes",
+                               &MultiThreadedFileReader::total_bytes)
+        .def_property_readonly("remaining_frames",
+                               &MultiThreadedFileReader::remaining_frames)
+        .def_property_readonly("next_read_frames",
+                               &MultiThreadedFileReader::next_read_frames)
+        .def_property_readonly("next_read_bytes",
+                               &MultiThreadedFileReader::next_read_bytes)
+        .def("seek", &MultiThreadedFileReader::seek, py::arg("frame_index"))
+        .def("tell", &MultiThreadedFileReader::tell)
+        .def("close", &MultiThreadedFileReader::close,
+             "Close all worker files. Safe to call more than once.")
+        .def_property_readonly(
+            "closed",
+            [](const MultiThreadedFileReader &self) { return !self.is_open(); })
+        .def("__len__", &MultiThreadedFileReader::total_frames)
+        .def(
+            "__enter__",
+            [](MultiThreadedFileReader &self) -> MultiThreadedFileReader * {
+                if (!self.is_open()) {
+                    throw std::runtime_error(
+                        "Cannot enter a closed MultiThreadedFileReader");
+                }
+                return &self;
+            },
+            py::return_value_policy::reference_internal)
+        .def("__exit__",
+             [](MultiThreadedFileReader &self, const py::object &,
+                const py::object &, const py::object &) {
+                 self.close();
+                 return false;
+             })
+        .def(
+            "__iter__", [](MultiThreadedFileReader &self) { return &self; },
+            py::return_value_policy::reference_internal)
+        .def("__next__", [](MultiThreadedFileReader &self) {
+            if (self.remaining_frames() == 0) {
+                throw py::stop_iteration();
+            }
+            return multi_threaded_reader_read(self, false);
+        });
+}

python/src/ctb_raw_file.hpp

@@ -9,7 +9,6 @@
 
 #include "aare/decode.hpp"
 #include "aare/defs.hpp"
-// #include "aare/fClusterFileV2.hpp"
 
 #include "np_helper.hpp"
 

python/src/file.hpp

@@ -7,7 +7,6 @@
 #include "aare/RawSubFile.hpp"
 
 #include "aare/defs.hpp"
-// #include "aare/fClusterFileV2.hpp"
 
 #include <cstdint>
 #include <filesystem>
@@ -195,52 +194,4 @@ void define_file_io_bindings(py::module &m) {
         });
 
 #pragma GCC diagnostic pop
-    // py::class_<ClusterHeader>(m, "ClusterHeader")
-    //     .def(py::init<>())
-    //     .def_readwrite("frame_number", &ClusterHeader::frame_number)
-    //     .def_readwrite("n_clusters", &ClusterHeader::n_clusters)
-    //     .def("__repr__", [](const ClusterHeader &a) { return "<ClusterHeader:
-    //     " + a.to_string() + ">"; });
-
-    // py::class_<ClusterV2_>(m, "ClusterV2_")
-    //     .def(py::init<>())
-    //     .def_readwrite("x", &ClusterV2_::x)
-    //     .def_readwrite("y", &ClusterV2_::y)
-    //     .def_readwrite("data", &ClusterV2_::data)
-    //     .def("__repr__", [](const ClusterV2_ &a) { return "<ClusterV2_: " +
-    //     a.to_string(false) + ">"; });
-
-    // py::class_<ClusterV2>(m, "ClusterV2")
-    //     .def(py::init<>())
-    //     .def_readwrite("cluster", &ClusterV2::cluster)
-    //     .def_readwrite("frame_number", &ClusterV2::frame_number)
-    //     .def("__repr__", [](const ClusterV2 &a) { return "<ClusterV2: " +
-    //     a.to_string() + ">"; });
-
-    // py::class_<ClusterFileV2>(m, "ClusterFileV2")
-    //     .def(py::init<const std::filesystem::path &, const std::string &>())
-    //     .def("read", py::overload_cast<>(&ClusterFileV2::read))
-    //     .def("read", py::overload_cast<int>(&ClusterFileV2::read))
-    //     .def("frame_number", &ClusterFileV2::frame_number)
-    //     .def("write", py::overload_cast<std::vector<ClusterV2> const
-    //     &>(&ClusterFileV2::write))
-
-    //     .def("close", &ClusterFileV2::close);
-
-    // m.def("to_clustV2", [](std::vector<DynamicCluster> &clusters, const int
-    // frame_number) {
-    //     std::vector<ClusterV2> clusters_;
-    //     for (auto &c : clusters) {
-    //         ClusterV2 cluster;
-    //         cluster.cluster.x = c.x;
-    //         cluster.cluster.y = c.y;
-    //         int i=0;
-    //         for(auto &d : cluster.cluster.data) {
-    //             d=c.get<double>(i++);
-    //         }
-    //         cluster.frame_number = frame_number;
-    //         clusters_.push_back(cluster);
-    //     }
-    //     return clusters_;
-    // });
-}
+}

python/src/fit.hpp

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

python/src/module.cpp

@@ -12,6 +12,7 @@
 #include "bind_Defs.hpp"
 #include "bind_Eta.hpp"
 #include "bind_Interpolator.hpp"
+#include "bind_MultiThreadedFileReader.hpp"
 #include "bind_PedestalTrackingPixelHistogram.hpp"
 #include "bind_PixelHistogram.hpp"
 #include "bind_PixelMap.hpp"
@@ -59,7 +60,11 @@ double, 'f' for float)
     define_ClusterCollector<T, N, M, U>(m, "Cluster" #N "x" #M #TYPE_CODE);
 
 PYBIND11_MODULE(_aare, m) {
+    auto experimental = m.def_submodule(
+        "experimental", "Experimental APIs that may change without notice");
+
     define_file_io_bindings(m);
+    define_multi_threaded_file_reader_bindings(experimental);
     define_raw_file_io_bindings(m);
     define_raw_sub_file_io_bindings(m);
     define_ctb_raw_file_io_bindings(m);

python/src/raw_master_file.hpp

@@ -8,7 +8,6 @@
 #include "aare/RawSubFile.hpp"
 
 #include "aare/defs.hpp"
-// #include "aare/fClusterFileV2.hpp"
 
 #include <cstdint>
 #include <filesystem>

python/src/raw_sub_file.hpp

@@ -7,7 +7,6 @@
 #include "aare/RawSubFile.hpp"
 
 #include "aare/defs.hpp"
-// #include "aare/fClusterFileV2.hpp"
 
 #include <cstdint>
 #include <filesystem>

python/src/var_cluster.hpp

@@ -1,17 +1,11 @@
 // SPDX-License-Identifier: MPL-2.0
 #include "aare/VarClusterFinder.hpp"
 #include "np_helper.hpp"
-// #include "aare/defs.hpp"
-// #include "aare/fClusterFileV2.hpp"
 
 #include <cstdint>
-// #include <filesystem>
 #include <pybind11/numpy.h>
-// #include <pybind11/iostream.h>
 #include <pybind11/pybind11.h>
 #include <pybind11/stl.h>
-// #include <pybind11/stl/filesystem.h>
-// #include <string>
 
 namespace py = pybind11;
 using namespace ::aare;

python/tests/test_Fit.py

@@ -0,0 +1,15 @@
+# SPDX-License-Identifier: MPL-2.0
+import numpy as np
+import aare
+
+
+
+def test_gaussian_model_evaluates_and_fits_data():
+    x = np.linspace(-5.0, 5.0, 51)
+    expected = np.array([20.0, 0.5, 1.2])
+    model = aare.Gaussian()
+
+    y = model(x, expected)
+    result = model.fit(x, y)
+
+    np.testing.assert_allclose(result["par"], expected, atol=2e-3)

python/tests/test_MultiThreadedFileReader.py

@@ -0,0 +1,147 @@
+# SPDX-License-Identifier: MPL-2.0
+import numpy as np
+import pytest
+
+import aare
+from aare.experimental import MultiThreadedFileReader
+
+
+@pytest.fixture
+def frame_file(tmp_path):
+    data = np.arange(10 * 2 * 3, dtype=np.uint16).reshape(10, 2, 3)
+    path = tmp_path / "frames.npy"
+    np.save(path, data)
+    return path, data
+
+
+def test_experimental_import_path():
+    assert aare.experimental.MultiThreadedFileReader is MultiThreadedFileReader
+    assert MultiThreadedFileReader.__module__ == "aare.experimental"
+    assert not hasattr(aare, "MultiThreadedFileReader")
+
+
+def test_reads_all_frames_in_order(frame_file):
+    path, expected = frame_file
+    reader = MultiThreadedFileReader(path, n_threads=2, chunk_size=3)
+
+    first = reader.read()
+    second = reader.read()
+    exhausted = reader.read()
+
+    assert np.array_equal(first, expected[:6])
+    assert np.array_equal(second, expected[6:])
+    assert exhausted.shape == (0, 2, 3)
+    assert first.dtype == np.uint16
+    assert reader.n_threads == 2
+    assert reader.chunk_size == 3
+    assert reader.total_frames == 10
+    assert reader.source_total_frames == 10
+    assert reader.rows == 2
+    assert reader.cols == 3
+    assert reader.bitdepth == 16
+    assert reader.dtype == np.dtype(np.uint16)
+    assert reader.bytes_per_frame == 12
+    assert reader.total_bytes == expected.nbytes
+    assert len(reader) == 10
+    assert reader.tell() == 10
+    assert reader.remaining_frames == 0
+    assert reader.next_read_frames == 0
+    assert reader.next_read_bytes == 0
+
+
+def test_total_frame_limit(frame_file):
+    path, expected = frame_file
+    reader = MultiThreadedFileReader(
+        path, n_threads=2, chunk_size=2, total_frames=7
+    )
+
+    assert np.array_equal(reader.read(), expected[:4])
+    assert np.array_equal(reader.read_all(), expected[4:7])
+
+
+def test_iteration_yields_one_chunk_per_thread(frame_file):
+    path, expected = frame_file
+    reader = MultiThreadedFileReader(path, n_threads=2, chunk_size=2)
+
+    batches = list(reader)
+
+    assert [len(batch) for batch in batches] == [4, 4, 2]
+    assert np.array_equal(np.concatenate(batches), expected)
+
+
+def test_seek_resets_iteration_position(frame_file):
+    path, expected = frame_file
+    reader = MultiThreadedFileReader(path, n_threads=2, chunk_size=2)
+
+    reader.read()
+    assert reader.tell() == 4
+
+    reader.seek(1)
+    assert reader.tell() == 1
+    assert np.array_equal(reader.read(), expected[1:5])
+
+    with pytest.raises(IndexError):
+        reader.seek(11)
+
+
+def test_context_manager_closes_worker_files(frame_file):
+    path, expected = frame_file
+
+    with MultiThreadedFileReader(path, n_threads=2, chunk_size=2) as reader:
+        assert not reader.closed
+        assert np.array_equal(reader.read(), expected[:4])
+
+    assert reader.closed
+    reader.close()
+    with pytest.raises(RuntimeError):
+        reader.read()
+    with pytest.raises(RuntimeError):
+        reader.seek(0)
+    with pytest.raises(RuntimeError):
+        with reader:
+            pass
+
+
+def test_explicit_zero_frame_limit(frame_file):
+    path, expected = frame_file
+    reader = MultiThreadedFileReader(
+        path, n_threads=8, chunk_size=3, total_frames=0
+    )
+
+    actual = reader.read()
+    assert actual.shape == (0, *expected.shape[1:])
+    assert actual.dtype == expected.dtype
+
+
+@pytest.mark.parametrize(
+    "dtype", [np.int8, np.int32, np.uint64, np.float32, np.float64]
+)
+def test_preserves_numpy_dtype(tmp_path, dtype):
+    expected = np.arange(4 * 2 * 3, dtype=dtype).reshape(4, 2, 3)
+    path = tmp_path / "typed-frames.npy"
+    np.save(path, expected)
+
+    reader = MultiThreadedFileReader(path, n_threads=2, chunk_size=3)
+    actual = reader.read()
+
+    assert actual.dtype == expected.dtype
+    assert reader.dtype == expected.dtype
+    assert np.array_equal(actual, expected)
+
+
+@pytest.mark.parametrize(
+    ("n_threads", "chunk_size", "total_frames"),
+    [(0, 1, None), (1, 0, None), (1, 1, 11)],
+)
+def test_invalid_configuration(
+    frame_file, n_threads, chunk_size, total_frames
+):
+    path, _ = frame_file
+
+    with pytest.raises(ValueError):
+        MultiThreadedFileReader(
+            path,
+            n_threads=n_threads,
+            chunk_size=chunk_size,
+            total_frames=total_frames,
+        )

python/tests/test_PixelHistogram.py

@@ -0,0 +1,101 @@
+# SPDX-License-Identifier: MPL-2.0
+import numpy as np
+import pytest
+
+from aare import PixelHistogram
+
+
+def _random_frames(rows, cols, n, xmin, xmax, seed=0):
+    rng = np.random.default_rng(seed)
+    return [rng.uniform(xmin - 0.25, xmax + 0.25, size=(rows, cols)).astype(np.float64)
+            for _ in range(n)]
+
+
+def _reference_hdata(frames, rows, cols, n_bins, xmin, xmax):
+    expected = np.zeros((rows, cols, n_bins), dtype=np.uint16)
+    inv_range = n_bins / (xmax - xmin)
+    for img in frames:
+        for r in range(rows):
+            for c in range(cols):
+                v = float(img[r, c])
+                if not (xmin <= v < xmax):
+                    continue
+                b = int((v - xmin) * inv_range)
+                if b >= n_bins:
+                    b = n_bins - 1
+                expected[r, c, b] += 1
+    return expected
+
+
+def test_async_fill_matches_reference():
+    rows, cols, n_bins = 5, 7, 8
+    xmin, xmax = 0.0, 2.0
+    frames = _random_frames(rows, cols, n=3, xmin=xmin, xmax=xmax, seed=1)
+
+    hist = PixelHistogram(rows=rows, cols=cols, n_bins=n_bins, xmin=xmin, xmax=xmax)
+    for img in frames:
+        hist.fill_async(img)
+
+    np.testing.assert_array_equal(
+        hist.values(),
+        _reference_hdata(frames, rows, cols, n_bins, xmin, xmax),
+    )
+
+
+
+
+def test_fill_async_copies_buffer():
+    # After fill_async returns, the caller should be free to mutate the
+    # numpy array without affecting the pending fill.
+    rows, cols, n_bins = 4, 4, 4
+    xmin, xmax = 0.0, 1.0
+    hist = PixelHistogram(rows=rows, cols=cols, n_bins=n_bins, xmin=xmin, xmax=xmax, n_threads=1, max_pending=8)
+
+    img = np.full((rows, cols), 0.1, dtype=np.float64)  # falls in bin 0
+    hist.fill_async(img)
+    # Mutate the original array immediately; this must not affect the
+    # value that was already enqueued.
+    img[:] = 0.9  # would be bin 3
+    hist.flush()
+
+    h = hist.values()
+    assert h.shape == (rows, cols, n_bins)
+    # Every pixel saw one value in bin 0, none elsewhere.
+    assert (h[:, :, 0] == 1).all()
+    assert (h[:, :, 1:] == 0).all()
+
+
+    def test_fill_async_rejects_wrong_shape():
+        hist = PixelHistogram(8, 8, 4, 0.0, 1.0)
+        bad = np.zeros((4, 4), dtype=np.float32)
+        with pytest.raises(ValueError):
+            hist.fill_async(bad)
+
+
+def test_hdata_flushes_pending():
+    # Submit several frames with a tiny queue and read hdata() without an
+    # explicit flush(); hdata() must drain everything first.
+    rows, cols, n_bins = 3, 3, 4
+    xmin, xmax = 0.0, 1.0
+    hist = PixelHistogram(rows=rows, cols=cols, n_bins=n_bins, xmin=xmin, xmax=xmax,
+                          n_threads=1, max_pending=1)
+    frames = _random_frames(rows, cols, n=8, xmin=xmin, xmax=xmax, seed=3)
+    for img in frames:
+        hist.fill_async(img)
+
+    h = hist.values()  # no explicit flush()
+    np.testing.assert_array_equal(
+        h, _reference_hdata(frames, rows, cols, n_bins, xmin, xmax)
+    )
+
+
+def test_bin_centers_and_edges():
+    n_bins = 5
+    xmin, xmax = 0.0, 1.0
+    hist = PixelHistogram(rows=2, cols=2, n_bins=n_bins, xmin=xmin, xmax=xmax)
+    edges = hist.bin_edges()
+    centers = hist.bin_centers()
+    assert edges.shape == (n_bins + 1,)
+    assert centers.shape == (n_bins,)
+    np.testing.assert_allclose(edges, np.linspace(xmin, xmax, n_bins + 1), atol=1e-6)
+    np.testing.assert_allclose(centers, 0.5 * (edges[:-1] + edges[1:]), atol=1e-6)

src/Dtype.cpp

@@ -92,6 +92,12 @@ Dtype::Dtype(Dtype::TypeIndex ti) : m_type(ti) {}
  */
 Dtype::Dtype(std::string_view sv) {
 
+    // NumPy uses '|' for data types whose byte order is not applicable,
+    // notably one-byte integer types.
+    if (!sv.empty() && sv.front() == '|') {
+        sv.remove_prefix(1);
+    }
+
     // Check if the file is using our native endianess
     if (auto pos = sv.find_first_of("<>"); pos != std::string_view::npos) {
         const auto endianess = [](const char c) {

src/Dtype.test.cpp

@@ -12,6 +12,9 @@ TEST_CASE("Construct from typeid") {
 }
 
 TEST_CASE("Construct from string") {
+    REQUIRE(Dtype("|i1") == typeid(int8_t));
+    REQUIRE(Dtype("|u1") == typeid(uint8_t));
+
     if (endian::native == endian::little) {
         REQUIRE(Dtype("<i1") == typeid(int8_t));
         REQUIRE(Dtype("<u1") == typeid(uint8_t));
@@ -53,4 +56,4 @@ TEST_CASE("Construct from string with endianess") {
 
 TEST_CASE("Convert to string") {
     REQUIRE(Dtype(typeid(int)).to_string() == "<i4");
-}
+}

src/File.cpp

@@ -73,10 +73,11 @@ size_t File::tell() const { return file_impl->tell(); }
 size_t File::rows() const { return file_impl->rows(); }
 size_t File::cols() const { return file_impl->cols(); }
 size_t File::bitdepth() const { return file_impl->bitdepth(); }
+Dtype File::dtype() const { return file_impl->dtype(); }
 size_t File::bytes_per_pixel() const {
     return file_impl->bitdepth() / bits_per_byte;
 }
 
 DetectorType File::detector_type() const { return file_impl->detector_type(); }
 
-} // namespace aare
+} // namespace aare

src/Fit.cpp

@@ -4,7 +4,6 @@
 #include "Minuit2/FunctionMinimum.h"
 #include "Minuit2/MnHesse.h"
 #include "Minuit2/MnMigrad.h"
-#include "Minuit2/MnPrint.h"
 #include "Minuit2/MnStrategy.h"
 #include "Minuit2/MnUserParameters.h"
 #include "aare/Models.hpp"
@@ -12,461 +11,11 @@
 #include "aare/utils/task.hpp"
 #include <array>
 #include <cmath>
-#include <lmcurve2.h>
-#include <lmfit.hpp>
 #include <memory>
 #include <stdexcept>
-#include <thread>
 #include <type_traits>
 
 namespace aare {
-
-namespace func {
-
-double gaus(const double x, const double *par) {
-    return par[0] * exp(-pow(x - par[1], 2) / (2 * pow(par[2], 2)));
-}
-
-NDArray<double, 1> gaus(NDView<double, 1> x, NDView<double, 1> par) {
-    NDArray<double, 1> y({x.shape(0)}, 0);
-    for (ssize_t i = 0; i < x.size(); i++) {
-        y(i) = gaus(x(i), par.data());
-    }
-    return y;
-}
-
-double pol1(const double x, const double *par) { return par[0] * x + par[1]; }
-
-NDArray<double, 1> pol1(NDView<double, 1> x, NDView<double, 1> par) {
-    NDArray<double, 1> y({x.shape()}, 0);
-    for (ssize_t i = 0; i < x.size(); i++) {
-        y(i) = pol1(x(i), par.data());
-    }
-    return y;
-}
-
-double scurve(const double x, const double *par) {
-    return (par[0] + par[1] * x) +
-           0.5 * (1 + erf((x - par[2]) / (sqrt(2) * par[3]))) *
-               (par[4] + par[5] * (x - par[2]));
-}
-
-NDArray<double, 1> scurve(NDView<double, 1> x, NDView<double, 1> par) {
-    NDArray<double, 1> y({x.shape()}, 0);
-    for (ssize_t i = 0; i < x.size(); i++) {
-        y(i) = scurve(x(i), par.data());
-    }
-    return y;
-}
-
-double scurve2(const double x, const double *par) {
-    return (par[0] + par[1] * x) +
-           0.5 * (1 - erf((x - par[2]) / (sqrt(2) * par[3]))) *
-               (par[4] + par[5] * (x - par[2]));
-}
-
-NDArray<double, 1> scurve2(NDView<double, 1> x, NDView<double, 1> par) {
-    NDArray<double, 1> y({x.shape()}, 0);
-    for (ssize_t i = 0; i < x.size(); i++) {
-        y(i) = scurve2(x(i), par.data());
-    }
-    return y;
-}
-
-} // namespace func
-
-NDArray<double, 1> fit_gaus(NDView<double, 1> x, NDView<double, 1> y) {
-    NDArray<double, 1> result = model::Gaussian::estimate_par(x, y);
-    lm_status_struct status;
-
-    lmcurve(result.size(), result.data(), x.size(), x.data(), y.data(),
-            aare::func::gaus, &lm_control_double, &status);
-
-    return result;
-}
-
-NDArray<double, 3> fit_gaus(NDView<double, 1> x, NDView<double, 3> y,
-                            int n_threads) {
-    NDArray<double, 3> result({y.shape(0), y.shape(1), 3}, 0);
-
-    auto process = [&x, &y, &result](ssize_t first_row, ssize_t last_row) {
-        for (ssize_t row = first_row; row < last_row; row++) {
-            for (ssize_t col = 0; col < y.shape(1); col++) {
-                NDView<double, 1> values(&y(row, col, 0), {y.shape(2)});
-                auto res = fit_gaus(x, values);
-                result(row, col, 0) = res(0);
-                result(row, col, 1) = res(1);
-                result(row, col, 2) = res(2);
-            }
-        }
-    };
-
-    auto tasks = split_task(0, y.shape(0), n_threads);
-    RunInParallel(process, tasks);
-    return result;
-}
-
-void fit_gaus(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
-              NDView<double, 1> par_out, NDView<double, 1> par_err_out,
-              double &chi2) {
-
-    // Check that we have the correct sizes
-    if (y.size() != x.size() || y.size() != y_err.size() ||
-        par_out.size() != 3 || par_err_out.size() != 3) {
-        throw std::runtime_error("Data, x, data_err must have the same size "
-                                 "and par_out, par_err_out must have size 3");
-    }
-
-    // /* Collection of output parameters for status info. */
-    // typedef struct {
-    //     double fnorm;  /* norm of the residue vector fvec. */
-    //     int nfev;      /* actual number of iterations. */
-    //     int outcome;   /* Status indicator. Nonnegative values are used as
-    //     index
-    //                     for the message text lm_infmsg, set in lmmin.c. */
-    //     int userbreak; /* Set when function evaluation requests termination.
-    //     */
-    // } lm_status_struct;
-
-    lm_status_struct status;
-    par_out = model::Gaussian::estimate_par(x, y);
-    std::array<double, 9> cov{0, 0, 0, 0, 0, 0, 0, 0, 0};
-
-    // void lmcurve2( const int n_par, double *par, double *parerr, double
-    // *covar, const int m_dat, const double *t, const double *y, const double
-    // *dy, double (*f)( const double ti, const double *par ), const
-    // lm_control_struct *control, lm_status_struct *status); n_par - Number of
-    // free variables. Length of parameter vector par. par - Parameter vector.
-    // On input, it must contain a reasonable guess. On output, it contains the
-    // solution found to minimize ||r||. parerr - Parameter uncertainties
-    // vector. Array of length n_par or NULL. On output, unless it or covar is
-    // NULL, it contains the weighted parameter uncertainties for the found
-    // parameters. covar - Covariance matrix. Array of length n_par * n_par or
-    // NULL. On output, unless it is NULL, it contains the covariance matrix.
-    // m_dat - Number of data points. Length of vectors t, y, dy. Must statisfy
-    // n_par <= m_dat. t - Array of length m_dat. Contains the abcissae (time,
-    // or "x") for which function f will be evaluated. y - Array of length
-    // m_dat. Contains the ordinate values that shall be fitted. dy - Array of
-    // length m_dat. Contains the standard deviations of the values y. f - A
-    // user-supplied parametric function f(ti;par). control - Parameter
-    // collection for tuning the fit procedure. In most cases, the default
-    // &lm_control_double is adequate. If f is only computed with
-    // single-precision accuracy, &lm_control_float should be used. Parameters
-    // are explained in lmmin2(3). status - A record used to return information
-    // about the minimization process: For details, see lmmin2(3).
-
-    lmcurve2(par_out.size(), par_out.data(), par_err_out.data(), cov.data(),
-             x.size(), x.data(), y.data(), y_err.data(), aare::func::gaus,
-             &lm_control_double, &status);
-
-    // Calculate chi2
-    chi2 = 0;
-    for (ssize_t i = 0; i < y.size(); i++) {
-        chi2 +=
-            std::pow((y(i) - func::gaus(x(i), par_out.data())) / y_err(i), 2);
-    }
-}
-
-void fit_gaus(NDView<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_err,
-              NDView<double, 3> par_out, NDView<double, 3> par_err_out,
-              NDView<double, 2> chi2_out,
-
-              int n_threads) {
-
-    auto process = [&](ssize_t first_row, ssize_t last_row) {
-        for (ssize_t row = first_row; row < last_row; row++) {
-            for (ssize_t col = 0; col < y.shape(1); col++) {
-                NDView<double, 1> y_view(&y(row, col, 0), {y.shape(2)});
-                NDView<double, 1> y_err_view(&y_err(row, col, 0),
-                                             {y_err.shape(2)});
-                NDView<double, 1> par_out_view(&par_out(row, col, 0),
-                                               {par_out.shape(2)});
-                NDView<double, 1> par_err_out_view(&par_err_out(row, col, 0),
-                                                   {par_err_out.shape(2)});
-
-                fit_gaus(x, y_view, y_err_view, par_out_view, par_err_out_view,
-                         chi2_out(row, col));
-            }
-        }
-    };
-
-    auto tasks = split_task(0, y.shape(0), n_threads);
-    RunInParallel(process, tasks);
-}
-
-void fit_pol1(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
-              NDView<double, 1> par_out, NDView<double, 1> par_err_out,
-              double &chi2) {
-
-    // Check that we have the correct sizes
-    if (y.size() != x.size() || y.size() != y_err.size() ||
-        par_out.size() != 2 || par_err_out.size() != 2) {
-        throw std::runtime_error("Data, x, data_err must have the same size "
-                                 "and par_out, par_err_out must have size 2");
-    }
-
-    lm_status_struct status;
-    par_out = model::Pol1::estimate_par(x, y);
-    std::array<double, 4> cov{0, 0, 0, 0};
-
-    lmcurve2(par_out.size(), par_out.data(), par_err_out.data(), cov.data(),
-             x.size(), x.data(), y.data(), y_err.data(), aare::func::pol1,
-             &lm_control_double, &status);
-
-    // Calculate chi2
-    chi2 = 0;
-    for (ssize_t i = 0; i < y.size(); i++) {
-        chi2 +=
-            std::pow((y(i) - func::pol1(x(i), par_out.data())) / y_err(i), 2);
-    }
-}
-
-void fit_pol1(NDView<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_err,
-              NDView<double, 3> par_out, NDView<double, 3> par_err_out,
-              NDView<double, 2> chi2_out, int n_threads) {
-
-    auto process = [&](ssize_t first_row, ssize_t last_row) {
-        for (ssize_t row = first_row; row < last_row; row++) {
-            for (ssize_t col = 0; col < y.shape(1); col++) {
-                NDView<double, 1> y_view(&y(row, col, 0), {y.shape(2)});
-                NDView<double, 1> y_err_view(&y_err(row, col, 0),
-                                             {y_err.shape(2)});
-                NDView<double, 1> par_out_view(&par_out(row, col, 0),
-                                               {par_out.shape(2)});
-                NDView<double, 1> par_err_out_view(&par_err_out(row, col, 0),
-                                                   {par_err_out.shape(2)});
-
-                fit_pol1(x, y_view, y_err_view, par_out_view, par_err_out_view,
-                         chi2_out(row, col));
-            }
-        }
-    };
-
-    auto tasks = split_task(0, y.shape(0), n_threads);
-    RunInParallel(process, tasks);
-}
-
-NDArray<double, 1> fit_pol1(NDView<double, 1> x, NDView<double, 1> y) {
-    // // Check that we have the correct sizes
-    // if (y.size() != x.size() || y.size() != y_err.size() ||
-    // par_out.size() != 2 || par_err_out.size() != 2) {
-    // throw std::runtime_error("Data, x, data_err must have the same size "
-    //                      "and par_out, par_err_out must have size 2");
-    // }
-    NDArray<double, 1> par = model::Pol1::estimate_par(x, y);
-
-    lm_status_struct status;
-    lmcurve(par.size(), par.data(), x.size(), x.data(), y.data(),
-            aare::func::pol1, &lm_control_double, &status);
-
-    return par;
-}
-
-NDArray<double, 3> fit_pol1(NDView<double, 1> x, NDView<double, 3> y,
-                            int n_threads) {
-    NDArray<double, 3> result({y.shape(0), y.shape(1), 2}, 0);
-
-    auto process = [&](ssize_t first_row, ssize_t last_row) {
-        for (ssize_t row = first_row; row < last_row; row++) {
-            for (ssize_t col = 0; col < y.shape(1); col++) {
-                NDView<double, 1> values(&y(row, col, 0), {y.shape(2)});
-                auto res = fit_pol1(x, values);
-                result(row, col, 0) = res(0);
-                result(row, col, 1) = res(1);
-            }
-        }
-    };
-
-    auto tasks = split_task(0, y.shape(0), n_threads);
-
-    RunInParallel(process, tasks);
-    return result;
-}
-
-// ~~ S-CURVES ~~
-
-// - No error
-NDArray<double, 1> fit_scurve(NDView<double, 1> x, NDView<double, 1> y) {
-    NDArray<double, 1> result = model::RisingScurve::estimate_par(x, y);
-    lm_status_struct status;
-
-    lmcurve(result.size(), result.data(), x.size(), x.data(), y.data(),
-            aare::func::scurve, &lm_control_double, &status);
-
-    return result;
-}
-
-NDArray<double, 3> fit_scurve(NDView<double, 1> x, NDView<double, 3> y,
-                              int n_threads) {
-    NDArray<double, 3> result({y.shape(0), y.shape(1), 6}, 0);
-
-    auto process = [&x, &y, &result](ssize_t first_row, ssize_t last_row) {
-        for (ssize_t row = first_row; row < last_row; row++) {
-            for (ssize_t col = 0; col < y.shape(1); col++) {
-                NDView<double, 1> values(&y(row, col, 0), {y.shape(2)});
-                auto res = fit_scurve(x, values);
-                result(row, col, 0) = res(0);
-                result(row, col, 1) = res(1);
-                result(row, col, 2) = res(2);
-                result(row, col, 3) = res(3);
-                result(row, col, 4) = res(4);
-                result(row, col, 5) = res(5);
-            }
-        }
-    };
-
-    auto tasks = split_task(0, y.shape(0), n_threads);
-    RunInParallel(process, tasks);
-    return result;
-}
-
-// - Error
-void fit_scurve(NDView<double, 1> x, NDView<double, 1> y,
-                NDView<double, 1> y_err, NDView<double, 1> par_out,
-                NDView<double, 1> par_err_out, double &chi2) {
-
-    // Check that we have the correct sizes
-    if (y.size() != x.size() || y.size() != y_err.size() ||
-        par_out.size() != 6 || par_err_out.size() != 6) {
-        throw std::runtime_error("Data, x, data_err must have the same size "
-                                 "and par_out, par_err_out must have size 6");
-    }
-
-    lm_status_struct status;
-    par_out = model::RisingScurve::estimate_par(x, y);
-    std::array<double, 36> cov = {0}; // size 6x6
-    // std::array<double, 4> cov{0, 0, 0, 0};
-
-    lmcurve2(par_out.size(), par_out.data(), par_err_out.data(), cov.data(),
-             x.size(), x.data(), y.data(), y_err.data(), aare::func::scurve,
-             &lm_control_double, &status);
-
-    // Calculate chi2
-    chi2 = 0;
-    for (ssize_t i = 0; i < y.size(); i++) {
-        chi2 +=
-            std::pow((y(i) - func::pol1(x(i), par_out.data())) / y_err(i), 2);
-    }
-}
-
-void fit_scurve(NDView<double, 1> x, NDView<double, 3> y,
-                NDView<double, 3> y_err, NDView<double, 3> par_out,
-                NDView<double, 3> par_err_out, NDView<double, 2> chi2_out,
-                int n_threads) {
-
-    auto process = [&](ssize_t first_row, ssize_t last_row) {
-        for (ssize_t row = first_row; row < last_row; row++) {
-            for (ssize_t col = 0; col < y.shape(1); col++) {
-                NDView<double, 1> y_view(&y(row, col, 0), {y.shape(2)});
-                NDView<double, 1> y_err_view(&y_err(row, col, 0),
-                                             {y_err.shape(2)});
-                NDView<double, 1> par_out_view(&par_out(row, col, 0),
-                                               {par_out.shape(2)});
-                NDView<double, 1> par_err_out_view(&par_err_out(row, col, 0),
-                                                   {par_err_out.shape(2)});
-
-                fit_scurve(x, y_view, y_err_view, par_out_view,
-                           par_err_out_view, chi2_out(row, col));
-            }
-        }
-    };
-
-    auto tasks = split_task(0, y.shape(0), n_threads);
-    RunInParallel(process, tasks);
-}
-
-// SCURVE2 ---
-
-// - No error
-NDArray<double, 1> fit_scurve2(NDView<double, 1> x, NDView<double, 1> y) {
-    NDArray<double, 1> result = model::FallingScurve::estimate_par(x, y);
-    lm_status_struct status;
-
-    lmcurve(result.size(), result.data(), x.size(), x.data(), y.data(),
-            aare::func::scurve2, &lm_control_double, &status);
-
-    return result;
-}
-
-NDArray<double, 3> fit_scurve2(NDView<double, 1> x, NDView<double, 3> y,
-                               int n_threads) {
-    NDArray<double, 3> result({y.shape(0), y.shape(1), 6}, 0);
-
-    auto process = [&x, &y, &result](ssize_t first_row, ssize_t last_row) {
-        for (ssize_t row = first_row; row < last_row; row++) {
-            for (ssize_t col = 0; col < y.shape(1); col++) {
-                NDView<double, 1> values(&y(row, col, 0), {y.shape(2)});
-                auto res = fit_scurve2(x, values);
-                result(row, col, 0) = res(0);
-                result(row, col, 1) = res(1);
-                result(row, col, 2) = res(2);
-                result(row, col, 3) = res(3);
-                result(row, col, 4) = res(4);
-                result(row, col, 5) = res(5);
-            }
-        }
-    };
-
-    auto tasks = split_task(0, y.shape(0), n_threads);
-    RunInParallel(process, tasks);
-    return result;
-}
-
-// - Error
-void fit_scurve2(NDView<double, 1> x, NDView<double, 1> y,
-                 NDView<double, 1> y_err, NDView<double, 1> par_out,
-                 NDView<double, 1> par_err_out, double &chi2) {
-
-    // Check that we have the correct sizes
-    if (y.size() != x.size() || y.size() != y_err.size() ||
-        par_out.size() != 6 || par_err_out.size() != 6) {
-        throw std::runtime_error("Data, x, data_err must have the same size "
-                                 "and par_out, par_err_out must have size 6");
-    }
-
-    lm_status_struct status;
-    par_out = model::FallingScurve::estimate_par(x, y);
-    std::array<double, 36> cov = {0}; // size 6x6
-    // std::array<double, 4> cov{0, 0, 0, 0};
-
-    lmcurve2(par_out.size(), par_out.data(), par_err_out.data(), cov.data(),
-             x.size(), x.data(), y.data(), y_err.data(), aare::func::scurve2,
-             &lm_control_double, &status);
-
-    // Calculate chi2
-    chi2 = 0;
-    for (ssize_t i = 0; i < y.size(); i++) {
-        chi2 +=
-            std::pow((y(i) - func::pol1(x(i), par_out.data())) / y_err(i), 2);
-    }
-}
-
-void fit_scurve2(NDView<double, 1> x, NDView<double, 3> y,
-                 NDView<double, 3> y_err, NDView<double, 3> par_out,
-                 NDView<double, 3> par_err_out, NDView<double, 2> chi2_out,
-                 int n_threads) {
-
-    auto process = [&](ssize_t first_row, ssize_t last_row) {
-        for (ssize_t row = first_row; row < last_row; row++) {
-            for (ssize_t col = 0; col < y.shape(1); col++) {
-                NDView<double, 1> y_view(&y(row, col, 0), {y.shape(2)});
-                NDView<double, 1> y_err_view(&y_err(row, col, 0),
-                                             {y_err.shape(2)});
-                NDView<double, 1> par_out_view(&par_out(row, col, 0),
-                                               {par_out.shape(2)});
-                NDView<double, 1> par_err_out_view(&par_err_out(row, col, 0),
-                                                   {par_err_out.shape(2)});
-
-                fit_scurve2(x, y_view, y_err_view, par_out_view,
-                            par_err_out_view, chi2_out(row, col));
-            }
-        }
-    };
-
-    auto tasks = split_task(0, y.shape(0), n_threads);
-    RunInParallel(process, tasks);
-}
-
 // ============================================================================
 // FitModel<Model> — method definitions
 // (constructor, destructor, copy, and all methods that touch Minuit2 state)
@@ -756,4 +305,4 @@ AARE_INSTANTIATE_FIT(model::FallingScurve)
 #undef AARE_INSTANTIATE_FIT
 // NOLINTEND
 
-} // namespace aare
+} // namespace aare

src/Fit.test.cpp

@@ -0,0 +1,81 @@
+// SPDX-License-Identifier: MPL-2.0
+#include "aare/Fit.hpp"
+#include "aare/FitModel.hpp"
+#include "aare/Models.hpp"
+
+#include <catch2/catch_approx.hpp>
+#include <catch2/catch_test_macros.hpp>
+#include <cmath>
+
+namespace {
+
+constexpr ssize_t n_points = 61;
+
+void fill_gaussian(aare::NDArray<double, 1> &x, aare::NDArray<double, 1> &y,
+                   double amplitude, double mean, double sigma) {
+    for (ssize_t i = 0; i < n_points; ++i) {
+        x(i) = -6.0 + 0.2 * static_cast<double>(i);
+        const double z = (x(i) - mean) / sigma;
+        y(i) = amplitude * std::exp(-0.5 * z * z);
+    }
+}
+
+} // namespace
+
+TEST_CASE("Minuit2 fits weighted and unweighted Gaussian data", "[fit]") {
+    aare::NDArray<double, 1> x({n_points});
+    aare::NDArray<double, 1> y({n_points});
+    aare::NDArray<double, 1> y_err({n_points}, 1.0);
+    fill_gaussian(x, y, 120.0, 0.8, 1.3);
+
+    const aare::FitModel<aare::model::Gaussian> unweighted_model;
+    const auto unweighted =
+        aare::fit_pixel(unweighted_model, x.view(), y.view());
+
+    REQUIRE(unweighted.size() == 4);
+    CHECK(unweighted(0) == Catch::Approx(120.0).epsilon(1e-5));
+    CHECK(unweighted(1) == Catch::Approx(0.8).epsilon(1e-5));
+    CHECK(unweighted(2) == Catch::Approx(1.3).epsilon(1e-5));
+    CHECK(unweighted(3) == Catch::Approx(0.0).margin(1e-4));
+
+    const aare::FitModel<aare::model::Gaussian> weighted_model(0, 100, 0.5,
+                                                               true);
+    const auto weighted =
+        aare::fit_pixel(weighted_model, x.view(), y.view(), y_err.view());
+
+    REQUIRE(weighted.size() == 7);
+    CHECK(weighted(0) == Catch::Approx(120.0).epsilon(1e-5));
+    CHECK(weighted(1) == Catch::Approx(0.8).epsilon(1e-5));
+    CHECK(weighted(2) == Catch::Approx(1.3).epsilon(1e-5));
+    CHECK(weighted(6) == Catch::Approx(0.0).margin(1e-4));
+}
+
+TEST_CASE("Minuit2 fits a Gaussian data cube in parallel", "[fit]") {
+    aare::NDArray<double, 1> x({n_points});
+    aare::NDArray<double, 1> values({n_points});
+    fill_gaussian(x, values, 80.0, -0.6, 0.9);
+
+    aare::NDArray<double, 3> y({2, 2, n_points});
+    for (ssize_t row = 0; row < 2; ++row) {
+        for (ssize_t col = 0; col < 2; ++col) {
+            for (ssize_t i = 0; i < n_points; ++i)
+                y(row, col, i) = values(i);
+        }
+    }
+
+    aare::NDArray<double, 3> par({2, 2, 3});
+    aare::NDArray<double, 2> chi2({2, 2});
+    const aare::FitModel<aare::model::Gaussian> model;
+
+    aare::fit_3d(model, x.view(), y.view(), aare::NDView<double, 3>{},
+                 par.view(), aare::NDView<double, 3>{}, chi2.view(), 2);
+
+    for (ssize_t row = 0; row < 2; ++row) {
+        for (ssize_t col = 0; col < 2; ++col) {
+            CHECK(par(row, col, 0) == Catch::Approx(80.0).epsilon(1e-5));
+            CHECK(par(row, col, 1) == Catch::Approx(-0.6).epsilon(1e-5));
+            CHECK(par(row, col, 2) == Catch::Approx(0.9).epsilon(1e-5));
+            CHECK(chi2(row, col) == Catch::Approx(0.0).margin(1e-8));
+        }
+    }
+}

src/MultiThreadedFileReader.cpp

@@ -0,0 +1,151 @@
+// SPDX-License-Identifier: MPL-2.0
+#include "aare/MultiThreadedFileReader.hpp"
+#include "aare/File.hpp"
+#include "aare/utils/math_helpers.hpp"
+
+#include <algorithm>
+#include <future>
+#include <limits>
+#include <stdexcept>
+#include <utility>
+
+namespace aare::experimental {
+namespace {
+
+size_t checked_product(size_t lhs, size_t rhs) {
+    if (lhs != 0 && rhs > std::numeric_limits<size_t>::max() / lhs) {
+        throw std::overflow_error(
+            "MultiThreadedFileReader buffer size overflow");
+    }
+    return lhs * rhs;
+}
+
+} // namespace
+
+MultiThreadedFileReader::MultiThreadedFileReader(
+    std::filesystem::path fname, size_t n_threads, size_t chunk_size,
+    std::optional<size_t> total_frames)
+    : m_fname(std::move(fname)), m_n_threads(n_threads),
+      m_chunk_size(chunk_size), m_total_frames(0), m_source_total_frames(0),
+      m_rows(0), m_cols(0), m_bitdepth(0), m_dtype(Dtype::NONE),
+      m_bytes_per_frame(0), m_total_bytes(0), m_current_frame(0) {
+    if (m_n_threads == 0) {
+        throw std::invalid_argument(
+            "MultiThreadedFileReader requires at least one thread");
+    }
+    if (m_chunk_size == 0) {
+        throw std::invalid_argument(
+            "MultiThreadedFileReader chunk size must be greater than zero");
+    }
+
+    File file(m_fname);
+    m_source_total_frames = file.total_frames();
+    m_total_frames = total_frames.value_or(m_source_total_frames);
+    if (m_total_frames > m_source_total_frames) {
+        throw std::invalid_argument(
+            "Requested frame count exceeds the number of frames in the file");
+    }
+
+    m_rows = file.rows();
+    m_cols = file.cols();
+    m_bitdepth = file.bitdepth();
+    m_dtype = file.dtype();
+    m_bytes_per_frame = file.bytes_per_frame();
+    m_total_bytes = checked_product(m_total_frames, m_bytes_per_frame);
+
+    m_files.reserve(m_n_threads);
+    m_files.push_back(std::move(file));
+    for (size_t i = 1; i < m_n_threads; ++i) {
+        m_files.emplace_back(m_fname);
+    }
+}
+
+size_t MultiThreadedFileReader::remaining_frames() const noexcept {
+    return m_total_frames - m_current_frame;
+}
+
+size_t MultiThreadedFileReader::next_read_frames() const noexcept {
+    const size_t remaining = remaining_frames();
+    if (remaining == 0) {
+        return 0;
+    }
+
+    const size_t chunks_remaining = ceil_div(remaining, m_chunk_size);
+    if (m_n_threads >= chunks_remaining) {
+        return remaining;
+    }
+    // This multiplication is safe: in this branch n_threads * chunk_size is
+    // strictly smaller than remaining.
+    return m_n_threads * m_chunk_size;
+}
+
+void MultiThreadedFileReader::ensure_open() const {
+    if (!is_open()) {
+        throw std::runtime_error("MultiThreadedFileReader is closed");
+    }
+}
+
+size_t MultiThreadedFileReader::read_into(std::byte *destination) {
+    ensure_open();
+    const size_t frames_to_read = next_read_frames();
+    if (frames_to_read == 0) {
+        return 0;
+    }
+    if (destination == nullptr) {
+        throw std::invalid_argument(
+            "MultiThreadedFileReader destination must not be null");
+    }
+
+    const size_t first_frame = m_current_frame;
+    const size_t active_threads = ceil_div(frames_to_read, m_chunk_size);
+    auto worker = [&](size_t worker_index) {
+        File &file = m_files[worker_index];
+        const size_t batch_offset = worker_index * m_chunk_size;
+        const size_t begin = first_frame + batch_offset;
+        const size_t count =
+            std::min(m_chunk_size, frames_to_read - batch_offset);
+        file.seek(begin);
+        file.read_into(destination + batch_offset * m_bytes_per_frame, count);
+    };
+
+    std::vector<std::future<void>> workers;
+    workers.reserve(active_threads);
+    for (size_t i = 0; i < active_threads; ++i) {
+        workers.emplace_back(std::async(std::launch::async, worker, i));
+    }
+
+    for (auto &future : workers) {
+        future.get();
+    }
+    m_current_frame += frames_to_read;
+    return frames_to_read;
+}
+
+std::vector<std::byte> MultiThreadedFileReader::read() {
+    std::vector<std::byte> data(next_read_bytes());
+    read_into(data.data());
+    return data;
+}
+
+std::vector<std::byte> MultiThreadedFileReader::read_all() {
+    ensure_open();
+    std::vector<std::byte> data(
+        checked_product(remaining_frames(), m_bytes_per_frame));
+    size_t offset = 0;
+    while (remaining_frames() != 0) {
+        const size_t frames_read = read_into(data.data() + offset);
+        offset += frames_read * m_bytes_per_frame;
+    }
+    return data;
+}
+
+void MultiThreadedFileReader::seek(size_t frame_index) {
+    ensure_open();
+    if (frame_index > m_total_frames) {
+        throw std::out_of_range(
+            "MultiThreadedFileReader frame index is out of range");
+    }
+    m_current_frame = frame_index;
+}
+
+} // namespace aare::experimental

src/MultiThreadedFileReader.test.cpp

@@ -0,0 +1,189 @@
+// SPDX-License-Identifier: MPL-2.0
+#include "aare/MultiThreadedFileReader.hpp"
+
+#include "aare/Dtype.hpp"
+#include "aare/File.hpp"
+#include "aare/FileInterface.hpp"
+#include "aare/Frame.hpp"
+#include "aare/NumpyFile.hpp"
+
+#include "test_config.hpp"
+#include <catch2/catch_test_macros.hpp>
+
+#include <chrono>
+#include <cstddef>
+#include <cstdint>
+#include <filesystem>
+#include <stdexcept>
+#include <string>
+#include <vector>
+
+using aare::File;
+using aare::FileConfig;
+using aare::Frame;
+using aare::NumpyFile;
+using aare::experimental::MultiThreadedFileReader;
+
+namespace {
+
+class TemporaryNumpyFile {
+  public:
+    TemporaryNumpyFile() {
+        const auto unique =
+            std::chrono::steady_clock::now().time_since_epoch().count();
+        m_path = std::filesystem::temp_directory_path() /
+                 ("aare-mt-reader-" + std::to_string(unique) + ".npy");
+
+        FileConfig cfg;
+        cfg.dtype = aare::Dtype::UINT16;
+        cfg.rows = 2;
+        cfg.cols = 3;
+        NumpyFile file(m_path, "w", cfg);
+        for (uint16_t frame_index = 0; frame_index < 10; ++frame_index) {
+            Frame frame(cfg.rows, cfg.cols, cfg.dtype);
+            auto image = frame.view<uint16_t>();
+            for (ssize_t row = 0; row < image.shape(0); ++row) {
+                for (ssize_t col = 0; col < image.shape(1); ++col) {
+                    image(row, col) = static_cast<uint16_t>(frame_index * 100 +
+                                                            row * 10 + col);
+                }
+            }
+            file.write(frame);
+        }
+    }
+
+    TemporaryNumpyFile(const TemporaryNumpyFile &) = delete;
+    TemporaryNumpyFile &operator=(const TemporaryNumpyFile &) = delete;
+
+    ~TemporaryNumpyFile() { std::filesystem::remove(m_path); }
+
+    const std::filesystem::path &path() const { return m_path; }
+    void truncate() { std::filesystem::resize_file(m_path, 0); }
+
+  private:
+    std::filesystem::path m_path;
+};
+
+std::vector<std::byte> read_reference(const std::filesystem::path &fpath,
+                                      size_t n_frames) {
+    File file(fpath);
+    std::vector<std::byte> data(n_frames * file.bytes_per_frame());
+    if (n_frames != 0) {
+        file.read_into(data.data(), n_frames);
+    }
+    return data;
+}
+
+} // namespace
+
+TEST_CASE("Multi-threaded reader preserves numpy frame order",
+          "[MultiThreadedFileReader]") {
+    TemporaryNumpyFile file;
+    const auto &fpath = file.path();
+
+    MultiThreadedFileReader reader(fpath, 2, 3);
+
+    CHECK(reader.n_threads() == 2);
+    CHECK(reader.chunk_size() == 3);
+    CHECK(reader.total_frames() == 10);
+    CHECK(reader.source_total_frames() == 10);
+    CHECK(reader.rows() == 2);
+    CHECK(reader.cols() == 3);
+    CHECK(reader.bitdepth() == 16);
+    CHECK(reader.total_bytes() ==
+          reader.total_frames() * reader.bytes_per_frame());
+
+    const auto reference = read_reference(fpath, reader.total_frames());
+    auto first = reader.read();
+    auto second = reader.read();
+
+    CHECK(first == std::vector<std::byte>(reference.begin(),
+                                          reference.begin() + 6 * 12));
+    CHECK(second ==
+          std::vector<std::byte>(reference.begin() + 6 * 12, reference.end()));
+    CHECK(reader.read().empty());
+    CHECK(reader.tell() == 10);
+    CHECK(reader.remaining_frames() == 0);
+}
+
+TEST_CASE("Multi-threaded reader handles uneven raw chunks and frame limits",
+          "[.with-data][MultiThreadedFileReader]") {
+    const auto fpath =
+        test_data_path() / "raw/jungfrau/jungfrau_single_master_0.json";
+    REQUIRE(std::filesystem::exists(fpath));
+
+    MultiThreadedFileReader reader(fpath, 2, 3, 9);
+
+    CHECK(reader.total_frames() == 9);
+    CHECK(reader.source_total_frames() == 10);
+    CHECK(reader.read_all() == read_reference(fpath, 9));
+}
+
+TEST_CASE("Multi-threaded reader can seek and reread",
+          "[MultiThreadedFileReader]") {
+    TemporaryNumpyFile file;
+
+    MultiThreadedFileReader reader(file.path(), 2, 1);
+    const auto first = reader.read();
+    CHECK(reader.tell() == 2);
+    CHECK(reader.next_read_frames() == 2);
+
+    reader.seek(0);
+    CHECK(reader.tell() == 0);
+    CHECK(reader.read() == first);
+    CHECK_THROWS_AS(reader.seek(11), std::out_of_range);
+}
+
+TEST_CASE("Multi-threaded reader validates its configuration",
+          "[MultiThreadedFileReader]") {
+    TemporaryNumpyFile file;
+    const auto &fpath = file.path();
+
+    CHECK_THROWS_AS(MultiThreadedFileReader(fpath, 0, 1),
+                    std::invalid_argument);
+    CHECK_THROWS_AS(MultiThreadedFileReader(fpath, 1, 0),
+                    std::invalid_argument);
+    CHECK_THROWS_AS(MultiThreadedFileReader(fpath, 1, 1, 11),
+                    std::invalid_argument);
+
+    MultiThreadedFileReader reader(fpath, 2, 2);
+    CHECK_THROWS_AS(reader.read_into(nullptr), std::invalid_argument);
+}
+
+TEST_CASE("An explicit zero frame limit produces an empty read",
+          "[MultiThreadedFileReader]") {
+    TemporaryNumpyFile file;
+
+    MultiThreadedFileReader reader(file.path(), 8, 3, 0);
+    CHECK(reader.total_frames() == 0);
+    CHECK(reader.total_bytes() == 0);
+    CHECK(reader.read().empty());
+    CHECK_NOTHROW(reader.read_into(nullptr));
+}
+
+TEST_CASE("read_into reads at most one chunk per worker",
+          "[MultiThreadedFileReader]") {
+    TemporaryNumpyFile file;
+    MultiThreadedFileReader reader(file.path(), 3, 2);
+    std::vector<std::byte> data(reader.next_read_bytes());
+
+    CHECK(reader.next_read_frames() == 6);
+    CHECK(reader.read_into(data.data()) == 6);
+    CHECK(reader.tell() == 6);
+    CHECK(reader.next_read_frames() == 4);
+    CHECK(reader.next_read_bytes() == 4 * reader.bytes_per_frame());
+}
+
+TEST_CASE("Multi-threaded reader can close its worker files",
+          "[MultiThreadedFileReader]") {
+    TemporaryNumpyFile file;
+    MultiThreadedFileReader reader(file.path(), 2, 2);
+
+    CHECK(reader.is_open());
+    reader.close();
+    CHECK_FALSE(reader.is_open());
+    CHECK_NOTHROW(reader.close());
+    CHECK_THROWS_AS(reader.read(), std::runtime_error);
+    CHECK_THROWS_AS(reader.read_all(), std::runtime_error);
+    CHECK_THROWS_AS(reader.seek(0), std::runtime_error);
+}

src/NumpyHelpers.test.cpp

@@ -34,6 +34,9 @@ TEST_CASE("trim whitespace") {
 }
 
 TEST_CASE("parse data type descriptions") {
+    REQUIRE(parse_descr("|i1") == aare::Dtype::INT8);
+    REQUIRE(parse_descr("|u1") == aare::Dtype::UINT8);
+
     REQUIRE(parse_descr("<i1") == aare::Dtype::INT8);
     REQUIRE(parse_descr("<i2") == aare::Dtype::INT16);
     REQUIRE(parse_descr("<i4") == aare::Dtype::INT32);
@@ -63,4 +66,4 @@ TEST_CASE("Parse numpy dict") {
     REQUIRE(map["descr"] == "'<f4'");
     REQUIRE(map["fortran_order"] == "False");
     REQUIRE(map["shape"] == "(3, 4)");
-}
+}

src/defs.test.cpp

@@ -22,20 +22,6 @@ TEST_CASE("Enum values") {
     REQUIRE(static_cast<int>(aare::DetectorType::Moench03) == 100);
 }
 
-TEST_CASE("DynamicCluster creation") {
-    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::DynamicCluster c2(c);
-    REQUIRE(c2.cluster_sizeX == 13);
-    REQUIRE(c2.cluster_sizeY == 15);
-    REQUIRE(c2.dt == aare::Dtype(typeid(int32_t)));
-    REQUIRE(c2.data() != nullptr);
-}
-
 TEST_CASE("Basic ops on BitOffset") {
     REQUIRE_THROWS(aare::BitOffset(10));