diff --git a/CMakeLists.txt b/CMakeLists.txt index d07e60cd..7def8473 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -458,6 +458,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/FilePtr.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 diff --git a/RELEASE.md b/RELEASE.md index 9a8c2060..886f5090 100644 --- a/RELEASE.md +++ b/RELEASE.md @@ -4,6 +4,10 @@ ### New Features: +- Added ``ClusterFile.frames()`` and ``ClusterFile.chunks()`` in C++ and + Python for iteration from the current file position. Frames preserve empty + frames and frame numbers; chunks support an optional positive size override. + Python's default file iteration continues to yield chunks. - Added ``FastPedestal`` in C++ and Python for per-pixel running mean, population variance, and standard deviation. It supports exponentially weighted updates, initialization from files, direct subtraction from NumPy @@ -22,6 +26,14 @@ ### API Changes: +- Added ``ClusterFile::read_frame(ClusterVector&)`` for allocation-reusing C++ + reads. It returns ``false`` at a clean end of file. The value-returning C++ + overload now returns ``std::optional>`` and + Python ``read_frame()`` returns ``None`` at end of file. Incomplete frames + still raise an error. +- Added an explicit boolean conversion to the C++ ``FilePtr`` type. +- Removed the public C++ ``ClusterFile::open`` method. Construct a new + ``ClusterFile`` to reopen a file or change its mode. - ``ClusterFinder`` now uses ``FastPedestal``. It must receive 1000 pedestal frames before cluster finding; ``find_clusters()`` raises an error until initialization is complete. Added ``update_threshold()`` to @@ -46,6 +58,16 @@ - ``TimingMode::Auto`` changed to ``TimingMode::AUTO_TIMING``, ``TimingMode::Trigger`` changed to ``TimingMode::TRIGGER_EXPOSURE`` ### Bugfixes: + +- Fixed a leaked empty ``ClusterVector`` at the end of Python ``ClusterFile`` + iteration. Chunk iteration now rejects a zero chunk size. +- ``ClusterFile::read_clusters`` and Python iteration now report incomplete + frame headers and cluster records instead of treating truncated files as a + clean end of file, with or without ROI or noise filtering. +- ``ClusterFile::write_frame`` now reports incomplete writes instead of + silently continuing with a truncated file. +- Gain-map application now checks the complete cluster footprint, preventing + out-of-bounds access for cluster sizes larger than 3x3. - ``RawFile`` now derives its frame count from the shortest selected raw subfile series across all ROIs. Frame-number reads use the same bounds, and Python ``len(reader)`` returns the adjusted count. A warning is printed when diff --git a/benchmarks/CMakeLists.txt b/benchmarks/CMakeLists.txt index 78031277..e874a003 100644 --- a/benchmarks/CMakeLists.txt +++ b/benchmarks/CMakeLists.txt @@ -18,8 +18,9 @@ FetchContent_MakeAvailable(benchmark) add_executable(benchmarks) target_sources( - benchmarks PRIVATE ndarray_benchmark.cpp ndview_benchmark.cpp - calculateeta_benchmark.cpp reduce_benchmark.cpp) + benchmarks + PRIVATE ndarray_benchmark.cpp ndview_benchmark.cpp calculateeta_benchmark.cpp + reduce_benchmark.cpp clusterfile_benchmark.cpp) # Link Google Benchmark and other necessary libraries target_link_libraries(benchmarks PRIVATE benchmark::benchmark aare_core diff --git a/benchmarks/clusterfile_benchmark.cpp b/benchmarks/clusterfile_benchmark.cpp new file mode 100644 index 00000000..7d757532 --- /dev/null +++ b/benchmarks/clusterfile_benchmark.cpp @@ -0,0 +1,94 @@ +// SPDX-License-Identifier: MPL-2.0 +#include "aare/ClusterFile.hpp" + +#include +#include +#include +#include +#include + +class ClusterFileFixture : public benchmark::Fixture { + using Cluster = aare::Cluster; + using File = aare::ClusterFile; + using Vector = aare::ClusterVector; + + std::filesystem::path m_path; + static constexpr int32_t frame_count = 256; + static constexpr size_t chunk_size = 1000; + + static void observe(Vector &clusters) { + benchmark::DoNotOptimize(clusters.data()); + benchmark::DoNotOptimize(clusters.size()); + } + + public: + void SetUp(const benchmark::State &state) override { + const auto stamp = + std::chrono::steady_clock::now().time_since_epoch().count(); + m_path = std::filesystem::temp_directory_path() / + ("aare-cluster-benchmark-" + std::to_string(stamp) + ".clust"); + File writer(m_path, chunk_size, "w"); + Vector frame(0); + frame.resize(state.range(0)); + for (int32_t number = 0; number < frame_count; ++number) { + frame.set_frame_number(number); + writer.write_frame(frame); + } + } + + void TearDown(const benchmark::State &) override { + std::error_code error; + std::filesystem::remove(m_path, error); + } + + template void read(benchmark::State &state) { + for (auto _ : state) { + File reader(m_path, chunk_size); + if constexpr (Iterate) { + auto consume = [](auto range) { + for (auto &clusters : range) { + observe(clusters); + } + }; + if constexpr (Frames) { + consume(reader.frames()); + } else { + consume(reader.chunks()); + } + } else if constexpr (Frames) { + Vector frame(0); + while (reader.read_frame(frame)) { + observe(frame); + } + } else { + while (true) { + auto chunk = reader.read_clusters(chunk_size); + if (chunk.empty()) { + break; + } + observe(chunk); + } + } + } + state.SetBytesProcessed(state.iterations() * frame_count * + state.range(0) * sizeof(Cluster)); + } +}; + +BENCHMARK_DEFINE_F(ClusterFileFixture, ReadFrames)(benchmark::State &state) { + read(state); +} +BENCHMARK_DEFINE_F(ClusterFileFixture, IterateFrames)(benchmark::State &state) { + read(state); +} +BENCHMARK_DEFINE_F(ClusterFileFixture, ReadChunks)(benchmark::State &state) { + read(state); +} +BENCHMARK_DEFINE_F(ClusterFileFixture, IterateChunks)(benchmark::State &state) { + read(state); +} + +BENCHMARK_REGISTER_F(ClusterFileFixture, ReadFrames)->Arg(64)->Arg(1024); +BENCHMARK_REGISTER_F(ClusterFileFixture, IterateFrames)->Arg(64)->Arg(1024); +BENCHMARK_REGISTER_F(ClusterFileFixture, ReadChunks)->Arg(64)->Arg(1024); +BENCHMARK_REGISTER_F(ClusterFileFixture, IterateChunks)->Arg(64)->Arg(1024); diff --git a/docs/src/ClusterFile.rst b/docs/src/ClusterFile.rst index a2ee162b..d764b687 100644 --- a/docs/src/ClusterFile.rst +++ b/docs/src/ClusterFile.rst @@ -1,8 +1,40 @@ ClusterFile ============= +Use ``frames()`` to iterate over complete frames, including empty frames and +frames whose clusters are all rejected by ROI or noise filtering. Each result +keeps its stored frame number; missing frame numbers are not synthesized. + +.. code-block:: cpp + + using ClusterType = aare::Cluster; + aare::ClusterFile file("clusters.clust"); + for (auto &frame : file.frames()) { + process_frame(frame.frame_number(), frame); + } + +Use ``chunks()`` for the constructor's chunk size, or ``chunks(10000)`` to +request a different size for that traversal. The size must be positive and +counts selected clusters after filtering. Chunks may split or combine frames, +so their frame numbers are not reliable per-cluster metadata. Empty chunks are +not yielded; only the final chunk can contain fewer clusters than requested. +Both ranges apply the configured gain map and report incomplete files as +errors, just like the explicit read methods. + +Ranges borrow the file and consume its current position without rewinding. +Constructing a range does not read; ``begin()`` reads the first result and +increment reads the next. Breaking a loop does not read ahead. Use one +traversal at a time and separate files for independent cursors. Switching to +frame reads after a chunk stops partway through a frame raises an error until +the remaining clusters in that frame have been read. + +The file must outlive its ranges and iterators and must not be moved while +they are in use. Iterators are move only and support C++17 range-based loops, +not algorithms requiring copyable STL input iterators. References to a result +last until advancement or iterator destruction. Move a result out with +``auto retained = std::move(frame)`` to retain it. Frame iteration reuses the +current vector's storage when it has not been moved out. + .. doxygenclass:: aare::ClusterFile :members: :undoc-members: - :private-members: - diff --git a/docs/src/python/file/pyClusterFile.rst b/docs/src/python/file/pyClusterFile.rst index cd391e0e..d6f3c8b9 100644 --- a/docs/src/python/file/pyClusterFile.rst +++ b/docs/src/python/file/pyClusterFile.rst @@ -1,26 +1,78 @@ ClusterFile -============ +=========== +The :func:`ClusterFile` factory is the main interface for reading and writing + cluster files. Use ``mode="r"`` to read, ``mode="w"`` to truncate and +write, or ``mode="a"`` to append. -The :class:`ClusterFile` class is the main interface to read and write clusters in aare. Unfortunately the -old file format does not include metadata like the cluster size and the data type. This means that the -user has to know this information from other sources. Specifying the wrong cluster size or data type -will lead to garbage data being read. +The format does not store the cluster dimensions, value type, coordinate type, +padding, or byte order. The ``cluster_size`` and ``dtype`` arguments must match +the writer; otherwise the bytes are interpreted incorrectly. + +Use ``frames()`` when frame boundaries and the frame number matter: + +.. code-block:: python + + with ClusterFile("clusters.clust") as file: + for frame in file.frames(): + print(frame.frame_number, frame.size) + +Each stored frame is yielded, including empty frames and frames whose clusters +are all rejected by filtering. Missing frame numbers are not synthesized. +The explicit ``read_frame()`` method returns one frame or ``None`` at EOF. + +Use ``chunks()`` to iterate over chunks of ``chunk_size`` clusters. Per default +the chunk size defined in the constructor is used. Alternatively one can pass +an explicit chunk size to the function for that specific traversal: + +.. code-block:: python + + with ClusterFile("clusters.clust", chunk_size=1000) as file: + for clusters in file.chunks(10000): + process(clusters) + +Chunk sizes must be positive and count selected clusters after filtering. +``file.chunks(10000)`` does not overwrite the constructor argument ``chunk_size`` +. Iterating directly +over ``file`` or calling ``next(file)`` still reads chunks using that default. +Chunks and ``read_clusters(n_clusters)`` may split or combine frames, so their +frame numbers are not reliable per-cluster metadata. Empty chunks are not +yielded. +Chunk reads and iteration raise an error if they encounter an incomplete frame +header or cluster record, including when ROI or noise filtering is enabled. +They return fewer clusters than requested only at a clean end of file. + +Both iterators consume the current file position without rewinding. Creating +an iterator does not read; each ``next()`` reads one result without prefetching +the following result. Use one traversal at a time and separate readers for +independent cursors. Switching to frames after a chunk stops partway through a +frame raises an error until the remaining clusters in that frame have been +read. Exhausted iterators continue to raise ``StopIteration``. + +Iterators keep their file object alive, but explicitly closing the file or +leaving its ``with`` block prevents further reads. Every yielded +``ClusterVector`` owns its storage and can be retained after iteration advances +or the file closes. NumPy views also remain valid across those operations; +as usual, resizing or modifying the allocation of their backing vector can +invalidate them. + +When a gain map is configured, it is applied to every cluster whose complete +footprint lies inside the map. Clusters whose footprint extends beyond the +gain-map boundaries remain in the returned cluster vector, but all their data +values are set to zero. .. py:currentmodule:: aare -.. autoclass:: ClusterFile - :members: - :undoc-members: - :inherited-members: +.. autofunction:: ClusterFile -Below is the API of the ClusterFile_Cluster3x3i but all variants share the same API. +Below is the API of ``ClusterFile_Cluster3x3i``; all compiled variants share the +same API. .. autoclass:: aare._aare.ClusterFile_Cluster3x3i :special-members: __init__ :members: :undoc-members: :show-inheritance: - :inherited-members: \ No newline at end of file + :inherited-members: diff --git a/include/aare/ClusterFile.hpp b/include/aare/ClusterFile.hpp index e464a254..507e2c9c 100644 --- a/include/aare/ClusterFile.hpp +++ b/include/aare/ClusterFile.hpp @@ -3,6 +3,7 @@ #include "aare/Cluster.hpp" #include "aare/ClusterVector.hpp" +#include "aare/FilePtr.hpp" #include "aare/GainMap.hpp" #include "aare/NDArray.hpp" #include "aare/ROI.hpp" @@ -11,38 +12,30 @@ #include #include #include +#include +#include namespace aare { -/* -Binary cluster file. Expects data to be laid out as: -int32_t frame_number -uint32_t number_of_clusters -int16_t x, int16_t y, int32_t data[9] x number_of_clusters -int32_t frame_number -uint32_t number_of_clusters -.... -*/ - -// TODO: change to support any type of clusters, e.g. header line with -// clsuter_size_x, cluster_size_y, /** - * @brief Class to read and write cluster files - * Expects data to be laid out as: + * @brief Read and write legacy binary cluster files. * + * Each frame is stored as: * * int32_t frame_number * uint32_t number_of_clusters - * int16_t x, int16_t y, int32_t data[9] * number_of_clusters - * int32_t frame_number - * uint32_t number_of_clusters - * etc. + * ClusterType clusters[number_of_clusters] + * + * The format stores clusters as their native in-memory representation and has + * no metadata describing the cluster dimensions, value type, coordinate type, + * padding, or byte order. Readers must therefore use the same ClusterType and + * a compatible platform ABI as the writer. */ template >> class ClusterFile { - FILE *fp{}; - const std::string m_filename{}; + FilePtr m_fp; + std::string m_filename{}; uint32_t m_num_left{}; /*Number of photons left in frame*/ size_t m_chunk_size{}; /*Number of clusters to read at a time*/ std::string m_mode; /*Mode to open the file in*/ @@ -54,48 +47,150 @@ class ClusterFile { public: /** - * @brief Construct a new Cluster File object - * @param fname path to the file - * @param chunk_size number of clusters to read at a time when iterating - * over the file - * @param mode mode to open the file in. "r" for reading, "w" for writing, - * "a" for appending - * @throws std::runtime_error if the file could not be opened + * @brief Single-pass range over frames or chunks of selected clusters. + * + * The range borrows the file and shares its current position with all + * other reads. Use only one traversal at a time. The file must outlive the + * range and its iterators and must not be moved while they are in use. + * Constructing a range does not read; begin() reads the first result. + * Calling begin() again starts at the file's then-current position. + */ + template class ReadRange { + ClusterFile *m_file; + size_t m_chunk_size; + + friend class ClusterFile; + ReadRange(ClusterFile &file, size_t chunk_size) + : m_file(&file), m_chunk_size(chunk_size) {} + + public: + struct Sentinel {}; + + /** + * @brief Move-only iterator for C++17 range-based loops. + * + * References to the current vector are valid until the iterator is + * advanced or destroyed. Move the vector out to retain its storage. + * Frame iteration reuses storage when the vector is not moved out. + * This is not a copyable STL input iterator. + */ + class Iterator { + ClusterFile *m_file; + size_t m_chunk_size; + ClusterVector m_clusters{0}; + + friend class ReadRange; + Iterator(ClusterFile &file, size_t chunk_size) + : m_file(&file), m_chunk_size(chunk_size) { + ++*this; + } + + public: + using value_type = ClusterVector; + + Iterator(const Iterator &) = delete; + Iterator &operator=(const Iterator &) = delete; + Iterator(Iterator &&) noexcept = default; + Iterator &operator=(Iterator &&) noexcept = default; + ~Iterator() = default; + + value_type &operator*() { return m_clusters; } + value_type *operator->() { return &m_clusters; } + + Iterator &operator++() { + if (m_file) { + if constexpr (ByFrame) { + if (!m_file->read_frame(m_clusters)) { + m_file = nullptr; + } + } else { + m_clusters = m_file->read_clusters(m_chunk_size); + if (m_clusters.empty()) { + m_file = nullptr; + } + } + } + return *this; + } + + void operator++(int) { ++*this; } + + friend bool operator==(const Iterator &it, Sentinel /*end*/) { + return it.m_file == nullptr; + } + friend bool operator!=(const Iterator &it, Sentinel end) { + return !(it == end); + } + friend bool operator==(Sentinel end, const Iterator &it) { + return it == end; + } + friend bool operator!=(Sentinel end, const Iterator &it) { + return !(it == end); + } + }; + + Iterator begin() const { return Iterator(*m_file, m_chunk_size); } + Sentinel end() const { return {}; } + }; + + using FrameRange = ReadRange; + using ChunkRange = ReadRange; + + /** + * @brief Open a cluster file. + * @param fname Path to the file. + * @param chunk_size Maximum number of selected clusters per chunk step. + * @param mode File mode: "r" to read, "w" to truncate and write, or "a" + * to append. + * @throws std::runtime_error If the mode is unsupported or the file cannot + * be opened. */ ClusterFile(const std::filesystem::path &fname, size_t chunk_size = 1000, const std::string &mode = "r") - : m_filename(fname.string()), m_chunk_size(chunk_size), m_mode(mode) { - - if (mode == "r") { - fp = fopen(m_filename.c_str(), "rb"); - if (!fp) { - throw std::runtime_error("Could not open file for reading: " + - m_filename); - } - } else if (mode == "w") { - fp = fopen(m_filename.c_str(), "wb"); - if (!fp) { - throw std::runtime_error("Could not open file for writing: " + - m_filename); - } - } else if (mode == "a") { - fp = fopen(m_filename.c_str(), "ab"); - if (!fp) { - throw std::runtime_error("Could not open file for appending: " + - m_filename); - } - } else { - throw std::runtime_error("Unsupported mode: " + mode); - } + : m_filename(fname.string()), m_chunk_size(chunk_size) { + open(mode); } - ~ClusterFile() { close(); } + /** + * @brief Iterate over complete frames from the current file position. + * + * Empty and fully filtered frames are yielded with their stored frame + * numbers. Reading stops only at a clean end of file. Advancing uses + * read_frame(), including its error for a prior partial-frame read. + */ + FrameRange frames() & { return FrameRange(*this, 0); } + FrameRange frames() && = delete; + + /** @brief Iterate using the chunk size supplied to the constructor. */ + ChunkRange chunks() & { return chunks(m_chunk_size); } + ChunkRange chunks() && = delete; /** - * @brief Read n_clusters clusters from the file discarding - * frame numbers. If EOF is reached the returned vector will - * have less than n_clusters clusters + * @brief Iterate over chunks from the current file position. + * @param chunk_size Maximum number of selected clusters per step. + * @throws std::invalid_argument If chunk_size is zero. + * @note Chunks can span frames; their frame numbers are not per-cluster + * metadata. Advancing uses read_clusters() and propagates its errors. + * Only the final chunk can be short. Empty chunks are not yielded. + */ + ChunkRange chunks(size_t chunk_size) & { + if (chunk_size == 0) { + throw std::invalid_argument("Chunk size must be greater than zero"); + } + return ChunkRange(*this, chunk_size); + } + ChunkRange chunks(size_t) && = delete; + + /** + * @brief Read up to n_clusters without preserving frame boundaries. + * @param n_clusters Maximum number of selected clusters to return. + * @return A cluster vector that may contain fewer clusters at a clean end + * of file. + * @note The returned vector may combine data from several frames, so its + * frame number must not be used as per-cluster metadata. + * @throws std::runtime_error If the file is not open for reading, an I/O + * error occurs, or an incomplete frame header or cluster record is read. */ ClusterVector read_clusters(size_t n_clusters) { if (m_mode != "r") { @@ -109,76 +204,116 @@ class ClusterFile { } /** - * @brief Read a single frame from the file and return the - * clusters. The cluster vector will have the frame number - * set. - * @throws std::runtime_error if the file is not opened for - * reading or the file pointer not at the beginning of a - * frame + * @brief Read the next complete frame. + * @return The selected clusters with the stored frame number set, or + * std::nullopt at a clean end of file before the next frame. + * @throws std::runtime_error If the file is not open for reading, a prior + * partial-frame read left clusters unread, or the frame is incomplete. + * @note A complete frame produces an engaged optional even when it contains + * no clusters or all of its clusters are removed by the configured filters. */ - ClusterVector read_frame() { + std::optional> read_frame() { + ClusterVector clusters(0); + if (!read_frame(clusters)) { + return std::nullopt; + } + return std::optional>{std::move(clusters)}; + } + + /** + * @brief Read the next complete frame into an existing cluster vector. + * @param clusters Destination whose storage is reused when large enough. + * Existing clusters are replaced after a frame header is read and remain + * unchanged at a clean end of file. + * @return true when a complete frame was read, or false at a clean end of + * file before the next frame. + * @throws std::runtime_error If the file is not open for reading, a prior + * partial-frame read left clusters unread, or a frame is incomplete. + * @note A complete frame is a successful read even when it contains no + * clusters or all of its clusters are removed by the configured filters. + */ + bool read_frame(ClusterVector &clusters) { if (m_mode != "r") { throw std::runtime_error(LOCATION + "File not opened for reading"); } - if (m_noise_map || m_roi) { - return read_frame_with_cut(); - } else { - return read_frame_without_cut(); + if (m_num_left) { + throw std::runtime_error( + LOCATION + "There are still clusters left in the last frame"); } + + if (m_noise_map || m_roi) { + return read_frame_with_cut(clusters); + } + return read_frame_without_cut(clusters); } + /** + * @brief Write one frame to the file. + * @param clusters Clusters to write, including their frame number. + * @throws std::runtime_error If the file is not open for writing or any + * part of the frame cannot be written completely. + */ void write_frame(const ClusterVector &clusters) { if (m_mode != "w" && m_mode != "a") { throw std::runtime_error("File not opened for writing"); } int32_t frame_number = clusters.frame_number(); - fwrite(&frame_number, sizeof(frame_number), 1, fp); + if (fwrite(&frame_number, sizeof(frame_number), 1, m_fp.get()) != 1) { + throw std::runtime_error(LOCATION + "Could not write frame number"); + } + uint32_t n_clusters = clusters.size(); - fwrite(&n_clusters, sizeof(n_clusters), 1, fp); - fwrite(clusters.data(), clusters.item_size(), clusters.size(), fp); + if (fwrite(&n_clusters, sizeof(n_clusters), 1, m_fp.get()) != 1) { + throw std::runtime_error(LOCATION + + "Could not write number of clusters"); + } + + if (fwrite(clusters.data(), clusters.item_size(), clusters.size(), + m_fp.get()) != clusters.size()) { + throw std::runtime_error(LOCATION + "Could not write clusters"); + } } /** - * @brief Return the chunk size + * @brief Return the default number of selected clusters per chunk step. */ size_t chunk_size() const { return m_chunk_size; } /** - * @brief Estimate the number of clusters in the file from its size + * @brief Estimate the number of clusters in the file from its size. * - * Frame headers are included in the estimate, so the result may be slightly - * larger than the actual number of clusters. The file position is not - * changed. - * - * @throws std::runtime_error if the file is not opened for reading + * Frame-header bytes are included in the estimate, so it may exceed the + * actual number of clusters. The file position is not changed. */ size_t estimate_n_clusters() const { return std::filesystem::file_size(m_filename) / sizeof(ClusterType); } /** - * @brief Set the region of interest to use when reading - * clusters. If set only clusters within the ROI will be - * read. + * @brief Select clusters by their center coordinate when reading. + * @param roi Half-open region of interest: [xmin, xmax) x [ymin, ymax). */ void set_roi(ROI roi) { m_roi = roi; } /** - * @brief Set the noise map to use when reading clusters. If - * set clusters below the noise level will be discarded. - * Selection criteria one of: Central pixel above noise, - * highest 2x2 sum above 2 * noise, total sum above 3 * - * noise. + * @brief Discard clusters that do not pass the noise thresholds. + * @param noise_map Per-pixel noise indexed as [y, x]. The map is copied. + * A cluster is retained only when its central pixel exceeds the local + * noise, its highest 2x2 sum exceeds twice the noise, and its total sum + * exceeds three times the noise. + * @warning The map must cover every cluster center coordinate in the file. */ void set_noise_map(const NDView noise_map) { m_noise_map = NDArray(noise_map); } /** - * @brief Set the gain map to use when reading clusters. If set the gain map - * will be applied to the clusters that pass ROI and noise_map selection. - * The gain map is expected to be in ADU/energy. + * @brief Apply a gain map to clusters selected while reading. + * @param gain_map Per-pixel gain in ADU/energy, indexed as [y, x]. The map + * is copied and inverted internally. + * @note Clusters whose complete footprint extends beyond the gain map are + * retained with all cluster data values set to zero. */ void set_gain_map(const NDView gain_map) { m_gain_map = InvertedGainMap(gain_map); @@ -193,65 +328,53 @@ class ClusterFile { } /** - * @brief Close the file. If not closed the file will be - * closed in the destructor + * @brief Close the file. + * + * Calling close more than once is safe. The destructor closes an open file + * automatically. */ void close() { - if (fp) { - fclose(fp); - fp = nullptr; - } + m_fp = FilePtr{}; + m_mode = ""; } /** - * @brief Return the current position in the file (bytes) + * @brief Return the current byte position in the file. + * @throws std::runtime_error If the file is closed or its position cannot + * be determined. */ int64_t tell() { - if (!fp) { + if (!m_fp.get()) { throw std::runtime_error(LOCATION + "File not opened"); } - return ftell(fp); + return m_fp.tell(); } + private: /** @brief Open the file in specific mode * */ void open(const std::string &mode) { - if (fp) { - close(); - } + close(); if (mode == "r") { - fp = fopen(m_filename.c_str(), "rb"); - if (!fp) { - throw std::runtime_error("Could not open file for reading: " + - m_filename); - } + m_fp = FilePtr(m_filename, "rb"); m_mode = "r"; } else if (mode == "w") { - fp = fopen(m_filename.c_str(), "wb"); - if (!fp) { - throw std::runtime_error("Could not open file for writing: " + - m_filename); - } + m_fp = FilePtr(m_filename, "wb"); m_mode = "w"; } else if (mode == "a") { - fp = fopen(m_filename.c_str(), "ab"); - if (!fp) { - throw std::runtime_error("Could not open file for appending: " + - m_filename); - } + m_fp = FilePtr(m_filename, "ab"); m_mode = "a"; } else { throw std::runtime_error("Unsupported mode: " + mode); } } - - private: ClusterVector read_clusters_with_cut(size_t n_clusters); ClusterVector read_clusters_without_cut(size_t n_clusters); - ClusterVector read_frame_with_cut(); - ClusterVector read_frame_without_cut(); + bool read_frame_header(int32_t &frame_number, uint32_t &n_clusters); + bool read_frame_with_cut(ClusterVector &clusters); + bool read_frame_without_cut(ClusterVector &clusters); bool is_selected(ClusterType &cl); ClusterType read_one_cluster(); }; @@ -281,25 +404,31 @@ ClusterFile::read_clusters_without_cut(size_t n_clusters) { } else { nn = nph; } - nph_read += fread((buf + nph_read), clusters.item_size(), nn, fp); + const auto count = + fread((buf + nph_read), clusters.item_size(), nn, m_fp.get()); + if (count != nn) { + throw std::runtime_error(LOCATION + "Could not read clusters"); + } + nph_read += count; m_num_left = nph - nn; // write back the number of photons left } if (nph_read < n_clusters) { // keep on reading frames and photons until reaching n_clusters - while (fread(&iframe, sizeof(iframe), 1, fp)) { + while (read_frame_header(iframe, nph)) { clusters.set_frame_number(iframe); - // read number of clusters in frame - if (fread(&nph, sizeof(nph), 1, fp)) { - if (nph > (n_clusters - nph_read)) - nn = n_clusters - nph_read; - else - nn = nph; + if (nph > (n_clusters - nph_read)) + nn = n_clusters - nph_read; + else + nn = nph; - nph_read += - fread((buf + nph_read), clusters.item_size(), nn, fp); - m_num_left = nph - nn; + const auto count = + fread((buf + nph_read), clusters.item_size(), nn, m_fp.get()); + if (count != nn) { + throw std::runtime_error(LOCATION + "Could not read clusters"); } + nph_read += count; + m_num_left = nph - nn; if (nph_read >= n_clusters) break; } @@ -339,16 +468,15 @@ ClusterFile::read_clusters_with_cut(size_t n_clusters) { } int32_t frame_number = 0; // frame number needs to be 4 bytes! - while (fread(&frame_number, sizeof(frame_number), 1, fp)) { - if (fread(&m_num_left, sizeof(m_num_left), 1, fp)) { - clusters.set_frame_number( - frame_number); // cluster vector will hold the last - // frame number - while (m_num_left && clusters.size() < n_clusters) { - ClusterType c = read_one_cluster(); - if (is_selected(c)) { - clusters.push_back(c); - } + uint32_t n_clusters_in_frame = 0; + while (read_frame_header(frame_number, n_clusters_in_frame)) { + m_num_left = n_clusters_in_frame; + clusters.set_frame_number( + frame_number); // cluster vector will hold the last frame number + while (m_num_left && clusters.size() < n_clusters) { + ClusterType c = read_one_cluster(); + if (is_selected(c)) { + clusters.push_back(c); } } @@ -366,7 +494,7 @@ ClusterFile::read_clusters_with_cut(size_t n_clusters) { template ClusterType ClusterFile::read_one_cluster() { ClusterType c; - auto rc = fread(&c, sizeof(c), 1, fp); + auto rc = fread(&c, sizeof(c), 1, m_fp.get()); if (rc != 1) { throw std::runtime_error(LOCATION + "Could not read cluster"); } @@ -375,73 +503,70 @@ ClusterType ClusterFile::read_one_cluster() { } template -ClusterVector -ClusterFile::read_frame_without_cut() { - if (m_mode != "r") { - throw std::runtime_error(LOCATION + "File not opened for reading"); +bool ClusterFile::read_frame_header(int32_t &frame_number, + uint32_t &n_clusters) { + const auto frame_number_bytes = + fread(&frame_number, 1, sizeof(frame_number), m_fp.get()); + if (frame_number_bytes == 0 && feof(m_fp.get()) && !ferror(m_fp.get())) { + return false; } - if (m_num_left) { - throw std::runtime_error( - LOCATION + "There are still photons left in the last frame"); - } - int32_t frame_number; - if (fread(&frame_number, sizeof(frame_number), 1, fp) != 1) { - if (feof(fp)) - throw std::runtime_error(LOCATION + "Unexpected end of file"); - else if (ferror(fp)) + if (frame_number_bytes != sizeof(frame_number)) { + if (ferror(m_fp.get())) { throw std::runtime_error(LOCATION + "Error reading from file"); - - throw std::runtime_error( - LOCATION + - "Unexpected error (not feof or ferror) when reading frame number"); + } + throw std::runtime_error(LOCATION + "Incomplete frame number"); } + const auto cluster_count_bytes = + fread(&n_clusters, 1, sizeof(n_clusters), m_fp.get()); + if (cluster_count_bytes != sizeof(n_clusters)) { + if (ferror(m_fp.get())) { + throw std::runtime_error(LOCATION + "Error reading from file"); + } + throw std::runtime_error(LOCATION + "Incomplete number of clusters"); + } + return true; +} + +template +bool ClusterFile::read_frame_without_cut( + ClusterVector &clusters) { + int32_t frame_number; uint32_t n_clusters; - if (fread(&n_clusters, sizeof(n_clusters), 1, fp) != 1) { - throw std::runtime_error(LOCATION + - "Could not read number of clusters"); + if (!read_frame_header(frame_number, n_clusters)) { + return false; } LOG(logDEBUG1) << "Reading " << n_clusters << " clusters from frame " << frame_number; - ClusterVector clusters(n_clusters); clusters.set_frame_number(frame_number); clusters.resize(n_clusters); LOG(logDEBUG1) << "clusters.item_size(): " << clusters.item_size(); - if (fread(clusters.data(), clusters.item_size(), n_clusters, fp) != + if (fread(clusters.data(), clusters.item_size(), n_clusters, m_fp.get()) != static_cast(n_clusters)) { throw std::runtime_error(LOCATION + "Could not read clusters"); } if (m_gain_map) m_gain_map->apply_gain_map(clusters); - return clusters; + return true; } template -ClusterVector -ClusterFile::read_frame_with_cut() { - if (m_mode != "r") { - throw std::runtime_error("File not opened for reading"); - } - if (m_num_left) { - throw std::runtime_error( - "There are still photons left in the last frame"); - } +bool ClusterFile::read_frame_with_cut( + ClusterVector &clusters) { int32_t frame_number; - if (fread(&frame_number, sizeof(frame_number), 1, fp) != 1) { - throw std::runtime_error("Could not read frame number"); + uint32_t n_clusters; + if (!read_frame_header(frame_number, n_clusters)) { + return false; } - if (fread(&m_num_left, sizeof(m_num_left), 1, fp) != 1) { - throw std::runtime_error("Could not read number of clusters"); - } - - ClusterVector clusters; - clusters.reserve(m_num_left); + m_num_left = n_clusters; + clusters.resize(0); + clusters.reserve(n_clusters); clusters.set_frame_number(frame_number); while (m_num_left) { ClusterType c = read_one_cluster(); @@ -451,7 +576,7 @@ ClusterFile::read_frame_with_cut() { } if (m_gain_map) m_gain_map->apply_gain_map(clusters); - return clusters; + return true; } template diff --git a/include/aare/FilePtr.hpp b/include/aare/FilePtr.hpp index 2cd08490..2e84b87c 100644 --- a/include/aare/FilePtr.hpp +++ b/include/aare/FilePtr.hpp @@ -18,6 +18,7 @@ class FilePtr { FilePtr &operator=(const FilePtr &) = delete; // since we handle a resource FilePtr(FilePtr &&other); FilePtr &operator=(FilePtr &&other); + explicit operator bool() const noexcept; FILE *get(); ssize_t tell(); void seek(ssize_t offset, int whence = SEEK_SET) { @@ -28,4 +29,4 @@ class FilePtr { ~FilePtr(); }; -} // namespace aare \ No newline at end of file +} // namespace aare diff --git a/include/aare/GainMap.hpp b/include/aare/GainMap.hpp index ff888915..3cc63077 100644 --- a/include/aare/GainMap.hpp +++ b/include/aare/GainMap.hpp @@ -34,29 +34,40 @@ class InvertedGainMap { template >> void apply_gain_map(ClusterVector &clustervec) { - // in principle we need to know the size of the image for this lookup - size_t ClusterSizeX = clustervec.cluster_size_x(); - size_t ClusterSizeY = clustervec.cluster_size_y(); - using T = typename ClusterVector::value_type; - int64_t index_cluster_center_x = ClusterSizeX / 2; - int64_t index_cluster_center_y = ClusterSizeY / 2; + constexpr size_t cluster_size_x = ClusterType::cluster_size_x; + constexpr size_t cluster_size_y = ClusterType::cluster_size_y; + constexpr ssize_t pixel_left_from_center = cluster_size_x / 2; + constexpr ssize_t pixels_right_from_center = + cluster_size_x - pixel_left_from_center - 1; + constexpr ssize_t pixels_top_from_center = cluster_size_y / 2; + constexpr ssize_t pixels_bottom_from_center = + cluster_size_y - pixels_top_from_center - 1; + for (size_t i = 0; i < clustervec.size(); i++) { auto &cl = clustervec[i]; + const auto center_x = static_cast(cl.x); + const auto center_y = static_cast(cl.y); - if (cl.x > 0 && cl.y > 0 && cl.x < m_gain_map.shape(1) - 1 && - cl.y < m_gain_map.shape(0) - 1) { - for (size_t j = 0; j < ClusterSizeX * ClusterSizeY; j++) { - size_t x = cl.x + j % ClusterSizeX - index_cluster_center_x; - size_t y = cl.y + j / ClusterSizeX - index_cluster_center_y; + if (center_x >= pixel_left_from_center && + center_y >= pixels_top_from_center && + center_x < m_gain_map.shape(1) - pixels_right_from_center && + center_y < m_gain_map.shape(0) - pixels_bottom_from_center) { + for (size_t j = 0; j < cluster_size_x * cluster_size_y; j++) { + const auto x = center_x + + static_cast(j % cluster_size_x) - + pixel_left_from_center; + const auto y = center_y + + static_cast(j / cluster_size_x) - + pixels_top_from_center; cl.data[j] = static_cast( static_cast(cl.data[j]) * m_gain_map( y, x)); // cast after conversion to keep precision } } else { - // clear edge clusters + // Clear clusters whose footprint extends beyond the gain map. cl.data.fill(0); } } @@ -66,4 +77,4 @@ class InvertedGainMap { NDArray m_gain_map{}; }; -} // end of namespace aare \ No newline at end of file +} // end of namespace aare diff --git a/python/aare/ClusterFinder.py b/python/aare/ClusterFinder.py index 9c0151b8..d612fcb9 100644 --- a/python/aare/ClusterFinder.py +++ b/python/aare/ClusterFinder.py @@ -89,19 +89,55 @@ def ClusterFileSink(clusterfindermt, cluster_file, dtype=np.int32): def ClusterFile(fname, cluster_size=(3,3), dtype=np.int32, chunk_size = 1000, mode = "r"): - """ - Factory function to create a ClusterFile object. Provides a cleaner syntax for - the templated ClusterFile in C++. + """Create a reader or writer for a binary cluster file. + + Parameters + ---------- + fname : path-like + Cluster file to open. + cluster_size : tuple[int, int], default=(3, 3) + Cluster dimensions stored in the file. + dtype : numpy dtype, default=numpy.int32 + Data type of the cluster values stored in the file. + chunk_size : int, default=1000 + Maximum number of selected clusters returned by ``chunks()`` and + default iteration. Must be positive when iterating over chunks. + mode : {"r", "w", "a"}, default="r" + Open for reading, truncate and write, or append, respectively. + + Returns + ------- + ClusterFile + The compiled ClusterFile specialization matching ``cluster_size`` and + ``dtype``. + + Notes + ----- + The file format contains no cluster shape or data-type metadata. Supplying + values that do not match the file causes its bytes to be interpreted + incorrectly. Use ``frames()`` to iterate over complete frames, including + empty or fully filtered frames with their stored frame numbers. Use + ``chunks()`` or ``chunks(chunk_size)`` to iterate over selected clusters + in batches. Chunks may combine frames, so their frame number is not + reliable per-cluster metadata. + + Iterators consume the current file position without rewinding; use one + traversal at a time. Each result owns its storage and remains valid after + advancing the iterator or closing the file. + + Examples + -------- .. code-block:: python from aare import ClusterFile - - with ClusterFile("clusters.clust", cluster_size=(3,3), dtype=np.int32) as cf: - # cf is now a ClusterFile_Cluster3x3i object but you don't need to know that. - for clusters in cf: - # Loop over clusters in chunks of 1000 - # The type of clusters will be a ClusterVector_Cluster3x3i in this case + + with ClusterFile( + "clusters.clust", cluster_size=(3, 3), dtype=np.int32 + ) as cf: + for clusters in cf.chunks(): + # Process clusters in chunks of at most 1000. + ... """ diff --git a/python/src/bind_ClusterFile.hpp b/python/src/bind_ClusterFile.hpp index 0831d64f..09fc649b 100644 --- a/python/src/bind_ClusterFile.hpp +++ b/python/src/bind_ClusterFile.hpp @@ -11,6 +11,7 @@ #include #include #include +#include // Disable warnings for unused parameters, as we ignore some // in the __exit__ method @@ -20,18 +21,90 @@ namespace py = pybind11; using namespace ::aare; +template class ClusterFileIterator { + Range m_range; + std::optional m_iterator; + bool m_done{false}; + + public: + explicit ClusterFileIterator(Range range) : m_range(std::move(range)) {} + + typename Range::Iterator::value_type next() { + if (m_done) { + throw py::stop_iteration(); + } + // Advance before returning the next result, never after yielding it. + if (m_iterator) { + ++*m_iterator; + } else { + m_iterator.emplace(m_range.begin()); + } + if (*m_iterator == m_range.end()) { + m_done = true; + throw py::stop_iteration(); + } + return std::move(**m_iterator); + } +}; + +template +void define_cluster_file_iterator(py::module &m, const std::string &name) { + py::class_(m, name.c_str()) + .def( + "__iter__", [](Iterator &self) -> Iterator & { return self; }, + py::return_value_policy::reference_internal) + .def("__next__", &Iterator::next); +} + template void define_ClusterFile(py::module &m, const std::string &typestr) { using ClusterType = Cluster; + using File = ClusterFile; + using FrameIterator = ClusterFileIterator; + using ChunkIterator = ClusterFileIterator; auto class_name = fmt::format("ClusterFile_{}", typestr); + define_cluster_file_iterator(m, + class_name + "_FrameIterator"); + define_cluster_file_iterator(m, + class_name + "_ChunkIterator"); - py::class_>(m, class_name.c_str()) + py::class_>( + m, class_name.c_str(), + "Read and write legacy binary cluster files. The format contains no " + "cluster type or shape metadata, so this class must match the file.") .def(py::init(), - py::arg(), py::arg("chunk_size") = 1000, py::arg("mode") = "r") + py::arg("fname"), py::arg("chunk_size") = 1000, + py::arg("mode") = "r", + "Open a cluster file. Mode must be 'r' to read, 'w' to truncate " + "and write, or 'a' to append.") + .def( + "frames", [](File &self) { return FrameIterator(self.frames()); }, + py::keep_alive<0, 1>(), + "Iterate over complete frames from the current position, including " + "empty and fully filtered frames with their stored frame numbers. " + "Results own their storage. The iterator keeps the file alive but " + "shares its cursor; use only one traversal at a time.") + .def( + "chunks", [](File &self) { return ChunkIterator(self.chunks()); }, + py::keep_alive<0, 1>(), + "Iterate from the current position using the constructor's chunk " + "size, which must be positive. Results own their storage. The " + "iterator keeps the file alive but shares its cursor; use only " + "one traversal at a time.") + .def( + "chunks", + [](File &self, size_t chunk_size) { + return ChunkIterator(self.chunks(chunk_size)); + }, + py::arg("chunk_size"), py::keep_alive<0, 1>(), + "Iterate over chunks of up to chunk_size selected clusters. The " + "size must be positive and does not change the constructor's " + "default. Chunks may span frames; their frame numbers are not " + "per-cluster metadata. Only the final chunk can be short.") .def( "read_clusters", [](ClusterFile &self, size_t n_clusters) { @@ -39,32 +112,60 @@ void define_ClusterFile(py::module &m, const std::string &typestr) { self.read_clusters(n_clusters)); return v; }, - py::return_value_policy::take_ownership, py::arg("n_clusters")) - .def("read_frame", - [](ClusterFile &self) { - auto v = new ClusterVector(self.read_frame()); - return v; - }) - .def("set_roi", &ClusterFile::set_roi, py::arg("roi")) - .def("tell", &ClusterFile::tell) + py::return_value_policy::take_ownership, py::arg("n_clusters"), + "Read up to n_clusters without preserving frame boundaries. The " + "result may combine frames, so its frame number is not reliable " + "per-cluster metadata.") + .def( + "read_frame", + [](ClusterFile &self) -> py::object { + auto clusters = self.read_frame(); + if (!clusters) { + return py::none(); + } + return py::cast( + new ClusterVector(std::move(*clusters)), + py::return_value_policy::take_ownership); + }, + "Read and return the next complete frame with its frame number, " + "or None at end of file.") + .def("set_roi", &ClusterFile::set_roi, py::arg("roi"), + "Select clusters whose centers lie within the half-open ROI.") + .def("tell", &ClusterFile::tell, + "Return the current byte position in the file.") .def("estimate_n_clusters", - &ClusterFile::estimate_n_clusters) + &ClusterFile::estimate_n_clusters, + "Estimate the number of clusters from the file size. Frame " + "headers can make this larger than the actual count.") .def( "set_noise_map", - [](ClusterFile &self, py::array_t noise_map) { + [](ClusterFile &self, + py::array_t noise_map) { auto view = make_view_2d(noise_map); self.set_noise_map(view); }, - py::arg("noise_map")) + py::arg("noise_map"), + "Set a two-dimensional, C-contiguous int32 noise map indexed as " + "[y, x]. The map must cover every cluster center coordinate.") - .def("set_gain_map", - [](ClusterFile &self, py::array_t gain_map) { - auto view = make_view_2d(gain_map); - self.set_gain_map(view); - }) + .def( + "set_gain_map", + [](ClusterFile &self, + py::array_t gain_map) { + auto view = make_view_2d(gain_map); + self.set_gain_map(view); + }, + py::arg("gain_map"), + "Set a two-dimensional, C-contiguous float64 gain map in " + "ADU/energy, indexed as [y, x]. Clusters whose complete footprint " + "extends beyond the map are retained with all data values set to " + "zero.") - .def("close", &ClusterFile::close) - .def("write_frame", &ClusterFile::write_frame) + .def("close", &ClusterFile::close, + "Close the file. Calling close more than once is safe.") + .def("write_frame", &ClusterFile::write_frame, + py::arg("clusters"), + "Write one ClusterVector, including its frame number.") .def("__enter__", [](ClusterFile &self) { return &self; }) .def("__exit__", [](ClusterFile &self, @@ -74,14 +175,8 @@ void define_ClusterFile(py::module &m, const std::string &typestr) { self.close(); }) .def("__iter__", [](ClusterFile &self) { return &self; }) - .def("__next__", [](ClusterFile &self) { - auto v = new ClusterVector( - self.read_clusters(self.chunk_size())); - if (v->size() == 0) { - throw py::stop_iteration(); - } - return v; - }); + .def("__next__", + [](File &self) { return ChunkIterator(self.chunks()).next(); }); } #pragma GCC diagnostic pop diff --git a/python/tests/test_ClusterFile.py b/python/tests/test_ClusterFile.py index 3e2fb81d..6ce8a11b 100644 --- a/python/tests/test_ClusterFile.py +++ b/python/tests/test_ClusterFile.py @@ -6,10 +6,279 @@ import boost_histogram as bh import time from pathlib import Path import pickle +import struct +import gc +import weakref -from aare import ClusterFile +from aare import ClusterFile, ROI from conftest import test_data_path + +@pytest.fixture +def iterator_file(tmp_path): + fname = tmp_path / "iteration.clust" + records = [ + (-11, []), + (42, [(5, 6, 10), (7, 8, 20)]), + (103, []), + (104, [(5, 6, 30), (7, 8, 40)]), + (105, []), + ] + fname.write_bytes(b"".join( + struct.pack("@iI", number, len(clusters)) + + b"".join(struct.pack("@HH9i", x, y, *([value] * 9)) + for x, y, value in clusters) + for number, clusters in records + )) + return fname + + +def test_frames_preserve_empty_frames_and_frame_numbers(iterator_file): + with ClusterFile(iterator_file) as reader: + frames = reader.frames() + assert iter(frames) is frames + assert reader.tell() == 0 + results = list(frames) + assert [frame.frame_number for frame in results] == [-11, 42, 103, 104, 105] + assert [frame.size for frame in results] == [0, 2, 0, 2, 0] + for _ in range(2): + with pytest.raises(StopIteration): + next(frames) + assert list(reader.frames()) == [] + + +@pytest.mark.parametrize("chunk_size, sizes", [(1, [1, 1, 1, 1]), (2, [2, 2]), + (3, [3, 1]), (10, [4])]) +@pytest.mark.parametrize("explicit", [False, True]) +def test_chunks_cross_frames(iterator_file, chunk_size, sizes, explicit): + with ClusterFile(iterator_file, chunk_size=chunk_size) as reader: + chunks = reader.chunks(chunk_size=chunk_size) if explicit else reader.chunks() + assert iter(chunks) is chunks + assert reader.tell() == 0 + results = list(chunks) + assert [chunk.size for chunk in results] == sizes + data = np.concatenate([np.asarray(chunk) for chunk in results]) + np.testing.assert_array_equal(data["x"], [5, 7, 5, 7]) + np.testing.assert_array_equal(data["data"][:, 0, 0], [10, 20, 30, 40]) + for _ in range(2): + with pytest.raises(StopIteration): + next(chunks) + + +def test_chunk_override_preserves_default_iteration(iterator_file): + with ClusterFile(iterator_file, chunk_size=2) as reader: + assert iter(reader) is reader + assert next(reader.chunks(1)).size == 1 + assert [chunk.size for chunk in reader] == [2, 1] + + +@pytest.mark.parametrize("use_roi", [False, True]) +@pytest.mark.parametrize("use_noise", [False, True]) +@pytest.mark.parametrize("method", ["frames", "chunks"]) +def test_iterators_apply_filters_and_gain(iterator_file, use_roi, use_noise, method): + with ClusterFile(iterator_file, chunk_size=3) as reader: + if use_roi: + reader.set_roi(ROI(0, 6, 0, 12)) + if use_noise: + reader.set_noise_map(np.full((12, 12), 15, dtype=np.int32)) + reader.set_gain_map(np.full((12, 12), 2.0)) + results = list(getattr(reader, method)()) + if method == "frames": + assert [frame.frame_number for frame in results] == [-11, 42, 103, 104, 105] + expected = [value / 2 for x, value in [(5, 10), (7, 20), (5, 30), (7, 40)] + if (not use_roi or x == 5) and (not use_noise or value > 15)] + data = np.concatenate([np.asarray(result) for result in results]) + np.testing.assert_array_equal(data["data"][:, 0, 0], expected) + + +@pytest.mark.parametrize("method", ["frames", "chunks"]) +def test_iterators_handle_all_rejected_clusters(iterator_file, method): + with ClusterFile(iterator_file) as reader: + reader.set_roi(ROI(0, 1, 0, 1)) + results = list(getattr(reader, method)()) + assert len(results) == (5 if method == "frames" else 0) + assert all(result.size == 0 for result in results) + + +@pytest.mark.parametrize("method", ["frames", "chunks"]) +def test_iterator_keeps_file_alive(iterator_file, method): + reader = ClusterFile(iterator_file) + owner = weakref.ref(reader) + iterator = getattr(reader, method)() + del reader + gc.collect() + assert owner() is not None + results = list(iterator) + assert sum(result.size for result in results) == 4 + del iterator + gc.collect() + assert owner() is None + + +@pytest.mark.parametrize("method", ["frames", "chunks"]) +def test_retained_results_and_numpy_views_own_storage(iterator_file, method): + with ClusterFile(iterator_file, chunk_size=2) as reader: + reader.read_frame() # Skip the initial empty frame. + iterator = getattr(reader, method)() + first = next(iterator) + array = np.asarray(first) + expected = array.copy() + results = list(iterator) + np.testing.assert_array_equal(np.asarray(first), expected) + del first, results, iterator, reader + gc.collect() + np.testing.assert_array_equal(array, expected) + + +def test_frame_iteration_resumes_without_reading_ahead(iterator_file): + with ClusterFile(iterator_file) as reader: + reader.read_frame() + for frame in reader.frames(): + assert frame.frame_number == 42 + break + assert reader.read_frame().frame_number == 103 + assert [frame.frame_number for frame in reader.frames()] == [104, 105] + + +def test_partial_chunk_must_be_completed_before_frame_iteration(iterator_file): + with ClusterFile(iterator_file) as reader: + for chunk in reader.chunks(1): + assert chunk.size == 1 + break + with pytest.raises(RuntimeError, match="clusters left"): + next(reader.frames()) + assert reader.read_clusters(1).size == 1 + assert [frame.frame_number for frame in reader.frames()] == [103, 104, 105] + + +@pytest.mark.parametrize("method", ["frames", "chunks"]) +def test_empty_iterators(tmp_path, method): + fname = tmp_path / "empty.clust" + fname.touch() + with ClusterFile(fname) as reader: + iterator = getattr(reader, method)() + for _ in range(2): + with pytest.raises(StopIteration): + next(iterator) + + +@pytest.mark.parametrize("method", ["frames", "chunks"]) +@pytest.mark.parametrize("started", [False, True]) +def test_iterators_reject_closed_files(iterator_file, method, started): + reader = ClusterFile(iterator_file, chunk_size=1) + iterator = getattr(reader, method)() + if started: + next(iterator) + reader.close() + with pytest.raises(RuntimeError, match="not opened for reading"): + next(iterator) + + +@pytest.mark.parametrize("method", ["frames", "chunks"]) +@pytest.mark.parametrize("mode", ["w", "a"]) +def test_iterators_reject_writing_modes(tmp_path, method, mode): + with ClusterFile(tmp_path / "output.clust", mode=mode) as reader: + with pytest.raises(RuntimeError, match="not opened for reading"): + next(getattr(reader, method)()) + + +def test_chunk_sizes_must_be_positive(iterator_file): + with ClusterFile(iterator_file, chunk_size=0) as reader: + with pytest.raises(ValueError, match="greater than zero"): + reader.chunks() + with pytest.raises(ValueError, match="greater than zero"): + reader.chunks(0) + with pytest.raises(ValueError, match="greater than zero"): + next(reader) + with pytest.raises(TypeError): + reader.chunks(-1) + assert reader.tell() == 0 + assert reader.read_clusters(0).size == 0 + assert len(list(reader.frames())) == 5 + + +@pytest.mark.parametrize("method", ["frames", "chunks"]) +def test_iterators_defer_read_errors_until_next_result(tmp_path, method): + fname = tmp_path / "truncated.clust" + cluster = struct.pack("@HH9i", 5, 6, *range(9)) + frame = struct.pack("@iI", 42, 1) + cluster + fname.write_bytes(frame + frame[:-1]) + with ClusterFile(fname, chunk_size=1) as reader: + iterator = getattr(reader, method)() + assert next(iterator).size == 1 + assert reader.tell() == len(frame) + with pytest.raises(RuntimeError): + next(iterator) + + +@pytest.mark.parametrize( + "shape, dtype", + [(shape, dtype) + for shape in [(2, 2), (3, 3), (5, 5), (7, 7), (9, 9)] + for dtype in [np.int32, np.float32, np.float64]] + + [((3, 3), np.int16)], +) +def test_iterators_are_bound_for_each_cluster_type(tmp_path, shape, dtype): + fname = tmp_path / "typed_frames.clust" + values = list(range(shape[0] * shape[1])) + record = struct.pack("@HH" + np.dtype(dtype).char * len(values), 5, 6, *values) + fname.write_bytes((struct.pack("@iI", -42, 1) + record) * 3) + with ClusterFile(fname, cluster_size=shape, dtype=dtype) as reader: + frame = next(reader.frames()) + assert frame.frame_number == -42 + np.testing.assert_array_equal(np.asarray(frame)["data"].reshape(-1), values) + chunks = list(reader.chunks()) + assert len(chunks) == 1 + assert chunks[0].size == 2 + np.testing.assert_array_equal(np.asarray(chunks[0])["data"].reshape(-1), values * 2) + + +def test_read_frame_returns_none_at_eof(tmp_path): + fname = tmp_path / "empty.clust" + fname.touch() + + with ClusterFile(fname) as f: + assert f.read_frame() is None + + +def test_read_frame_raises_for_malformed_file(tmp_path): + fname = tmp_path / "malformed.clust" + fname.write_bytes(b"\x00") + + with ClusterFile(fname) as f, pytest.raises(RuntimeError): + f.read_frame() + + +@pytest.mark.parametrize("use_roi", [False, True]) +@pytest.mark.parametrize("method", ["read_clusters", "default", "frames", "chunks"]) +@pytest.mark.parametrize( + "size", + [ + pytest.param(3, id="partial-frame-number"), + pytest.param(4, id="missing-cluster-count"), + pytest.param(7, id="partial-cluster-count"), + pytest.param(48, id="missing-cluster-record"), + pytest.param(87, id="partial-cluster-record"), + ], +) +def test_chunk_reads_reject_incomplete_frames(tmp_path, use_roi, method, size): + fname = tmp_path / "incomplete.clust" + cluster = struct.pack("@HH9i", 5, 6, *range(9)) + frame = struct.pack("@iI", 42, 2) + cluster * 2 + fname.write_bytes(frame[:size]) + + with ClusterFile(fname) as reader: + if use_roi: + reader.set_roi(ROI(0, 10, 0, 10)) + with pytest.raises(RuntimeError): + if method == "default": + next(reader) + elif method == "read_clusters": + reader.read_clusters(10) + else: + next(getattr(reader, method)()) + + @pytest.mark.withdata def test_cluster_file(test_data_path): """Test ClusterFile""" diff --git a/src/ClusterFile.test.cpp b/src/ClusterFile.test.cpp index 15caf6bb..6c009ba1 100644 --- a/src/ClusterFile.test.cpp +++ b/src/ClusterFile.test.cpp @@ -5,12 +5,74 @@ #include "aare/defs.hpp" #include #include +#include +#include +#include #include +#include +#include +#include using aare::Cluster; using aare::ClusterFile; using aare::ClusterVector; +namespace { + +class TemporaryClusterFile { + public: + TemporaryClusterFile() { + const auto unique = + std::chrono::steady_clock::now().time_since_epoch().count(); + m_path = std::filesystem::temp_directory_path() / + ("aare-cluster-file-" + std::to_string(unique) + ".clust"); + } + + TemporaryClusterFile(const TemporaryClusterFile &) = delete; + TemporaryClusterFile &operator=(const TemporaryClusterFile &) = delete; + + ~TemporaryClusterFile() { + std::error_code error; + std::filesystem::remove(m_path, error); + } + + const std::filesystem::path &path() const { return m_path; } + + private: + std::filesystem::path m_path; +}; + +using TestCluster = Cluster; + +ClusterVector make_test_frame(int32_t frame_number, + double offset) { + ClusterVector clusters(2, frame_number); + clusters.push_back( + TestCluster{5, + 6, + {offset, offset + 1, offset + 2, offset + 3, offset + 4, + offset + 5, offset + 6, offset + 7, offset + 8}}); + clusters.push_back(TestCluster{7, + 8, + {offset + 9, offset + 10, offset + 11, + offset + 12, offset + 13, offset + 14, + offset + 15, offset + 16, offset + 17}}); + return clusters; +} + +void check_frame(const ClusterVector &actual, + const ClusterVector &expected) { + REQUIRE(actual.frame_number() == expected.frame_number()); + REQUIRE(actual.size() == expected.size()); + for (size_t i = 0; i < expected.size(); ++i) { + CHECK(actual[i].x == expected[i].x); + CHECK(actual[i].y == expected[i].y); + CHECK(actual[i].data == expected[i].data); + } +} + +} // namespace + TEST_CASE("Read one frame from a cluster file", "[.with-data]") { // We know that the frame has 97 clusters auto fpath = test_data_path() / "clust" / "single_frame_97_clustrers.clust"; @@ -20,12 +82,14 @@ TEST_CASE("Read one frame from a cluster file", "[.with-data]") { CHECK(f.estimate_n_clusters() == 97); CHECK(f.tell() == 0); auto clusters = f.read_frame(); - CHECK(clusters.size() == 97); - CHECK(clusters.frame_number() == 135); - CHECK(clusters[0].x == 1); - CHECK(clusters[0].y == 200); + REQUIRE(clusters); + CHECK(clusters->size() == 97); + CHECK(clusters->frame_number() == 135); + CHECK((*clusters)[0].x == 1); + CHECK((*clusters)[0].y == 200); int32_t expected_cluster_data[] = {0, 1, 2, 3, 4, 5, 6, 7, 8}; - CHECK(std::equal(std::begin(clusters[0].data), std::end(clusters[0].data), + CHECK(std::equal(std::begin((*clusters)[0].data), + std::end((*clusters)[0].data), std::begin(expected_cluster_data))); } @@ -42,22 +106,24 @@ TEST_CASE("Read one frame using ROI", "[.with-data]") { roi.ymax = 249; f.set_roi(roi); auto clusters = f.read_frame(); - REQUIRE(clusters.size() == 49); - REQUIRE(clusters.frame_number() == 135); + REQUIRE(clusters); + REQUIRE(clusters->size() == 49); + REQUIRE(clusters->frame_number() == 135); // Check that all clusters are within the ROI - for (size_t i = 0; i < clusters.size(); i++) { - auto c = clusters[i]; + for (size_t i = 0; i < clusters->size(); i++) { + auto c = (*clusters)[i]; REQUIRE(c.x >= roi.xmin); REQUIRE(c.x <= roi.xmax); REQUIRE(c.y >= roi.ymin); REQUIRE(c.y <= roi.ymax); } - CHECK(clusters[0].x == 1); - CHECK(clusters[0].y == 200); + CHECK((*clusters)[0].x == 1); + CHECK((*clusters)[0].y == 200); int32_t expected_cluster_data[] = {0, 1, 2, 3, 4, 5, 6, 7, 8}; - CHECK(std::equal(std::begin(clusters[0].data), std::end(clusters[0].data), + CHECK(std::equal(std::begin((*clusters)[0].data), + std::end((*clusters)[0].data), std::begin(expected_cluster_data))); } @@ -284,70 +350,551 @@ TEST_CASE("Read cluster from multiple frame file", "[.with-data]") { } } -TEST_CASE("Write cluster with potential padding", - "[.with-data][.ClusterFile]") { +TEST_CASE("ClusterFile flushes a frame when the writer is destroyed", + "[ClusterFile]") { + TemporaryClusterFile file; + auto expected = make_test_frame(42, 0.0); - using ClusterType = Cluster; + { + ClusterFile writer(file.path(), 1000, "w"); + writer.write_frame(expected); + } - REQUIRE(std::filesystem::exists(test_data_path() / "clust")); + const auto expected_file_size = sizeof(int32_t) + sizeof(uint32_t) + + expected.size() * sizeof(TestCluster); + CHECK(std::filesystem::file_size(file.path()) == expected_file_size); - auto fpath = test_data_path() / "clust" / "single_frame_2_clusters.clust"; - - ClusterFile file(fpath, 1000, "w"); - - ClusterVector clustervec(2); - uint16_t coordinate = 5; - clustervec.push_back(ClusterType{ - coordinate, coordinate, {0., 0., 0., 0., 0., 0., 0., 0., 0.}}); - clustervec.push_back(ClusterType{ - coordinate, coordinate, {0., 0., 0., 0., 0., 0., 0., 0., 0.}}); - - file.write_frame(clustervec); - - file.close(); - - file.open("r"); - - auto read_cluster_vector = file.read_frame(); - - CHECK(read_cluster_vector.size() == 2); - CHECK(read_cluster_vector.frame_number() == 0); - - CHECK(read_cluster_vector[0].x == clustervec[0].x); - CHECK(read_cluster_vector[0].y == clustervec[0].y); - CHECK(std::equal( - clustervec[0].data.begin(), clustervec[0].data.end(), - read_cluster_vector[0].data.begin(), [](double a, double b) { - return std::abs(a - b) < std::numeric_limits::epsilon(); - })); - - CHECK(read_cluster_vector[1].x == clustervec[1].x); - CHECK(read_cluster_vector[1].y == clustervec[1].y); - CHECK(std::equal( - clustervec[1].data.begin(), clustervec[1].data.end(), - read_cluster_vector[1].data.begin(), [](double a, double b) { - return std::abs(a - b) < std::numeric_limits::epsilon(); - })); + ClusterFile reader(file.path()); + auto actual = reader.read_frame(); + REQUIRE(actual); + check_frame(*actual, expected); } -TEST_CASE("Read frame and modify cluster data", "[.with-data][.ClusterFile]") { +TEST_CASE("ClusterFile appends frames", "[ClusterFile]") { + TemporaryClusterFile file; + auto first_expected = make_test_frame(42, 0.0); + auto second_expected = make_test_frame(43, 100.0); + + { + ClusterFile writer(file.path(), 1000, "w"); + writer.write_frame(first_expected); + } + { + ClusterFile writer(file.path(), 1000, "a"); + writer.write_frame(second_expected); + } + + ClusterFile reader(file.path()); + auto first_actual = reader.read_frame(); + auto second_actual = reader.read_frame(); + REQUIRE(first_actual); + REQUIRE(second_actual); + check_frame(*first_actual, first_expected); + check_frame(*second_actual, second_expected); +} + +TEST_CASE("ClusterFile reads frames into an existing cluster vector", + "[ClusterFile]") { + TemporaryClusterFile file; + auto first_expected = make_test_frame(42, 0.0); + auto second_expected = make_test_frame(43, 100.0); + + { + ClusterFile writer(file.path(), 1000, "w"); + writer.write_frame(first_expected); + writer.write_frame(second_expected); + } + + ClusterFile reader(file.path()); + ClusterVector actual(2); + const auto initial_data = actual.data(); + + REQUIRE(reader.read_frame(actual)); + CHECK(actual.data() == initial_data); + check_frame(actual, first_expected); + + REQUIRE(reader.read_frame(actual)); + CHECK(actual.data() == initial_data); + check_frame(actual, second_expected); + + CHECK_FALSE(reader.read_frame(actual)); + check_frame(actual, second_expected); +} + +TEST_CASE("ClusterFile reports a clean end of file", "[ClusterFile]") { + TemporaryClusterFile file; + { + ClusterFile writer(file.path(), 1000, "w"); + } + + ClusterFile reader(file.path()); + ClusterVector clusters; + CHECK_FALSE(reader.read_frame(clusters)); + CHECK(clusters.empty()); + CHECK_FALSE(reader.read_frame().has_value()); +} + +TEST_CASE("ClusterFile reports complete frames with no output clusters", + "[ClusterFile]") { + TemporaryClusterFile file; + ClusterVector empty_frame(0, 41); + auto expected = make_test_frame(42, 0.0); + { + ClusterFile writer(file.path(), 1000, "w"); + writer.write_frame(empty_frame); + writer.write_frame(expected); + } + + ClusterFile reader(file.path()); + ClusterVector clusters; + auto empty_actual = reader.read_frame(); + REQUIRE(empty_actual); + CHECK(empty_actual->empty()); + CHECK(empty_actual->frame_number() == 41); + + aare::ROI roi; + roi.xmin = 100; + roi.xmax = 200; + roi.ymin = 100; + roi.ymax = 200; + reader.set_roi(roi); + + REQUIRE(reader.read_frame(clusters)); + CHECK(clusters.empty()); + CHECK(clusters.frame_number() == expected.frame_number()); + CHECK_FALSE(reader.read_frame(clusters)); +} + +TEST_CASE("ClusterFile rejects incomplete frames", "[ClusterFile]") { + TemporaryClusterFile file; + auto expected = make_test_frame(42, 0.0); + { + ClusterFile writer(file.path(), 1000, "w"); + writer.write_frame(expected); + } + + SECTION("incomplete frame number") { + std::filesystem::resize_file(file.path(), sizeof(int32_t) - 1); + } + SECTION("incomplete cluster count") { + std::filesystem::resize_file(file.path(), + sizeof(int32_t) + sizeof(uint32_t) - 1); + } + SECTION("incomplete cluster data") { + std::filesystem::resize_file( + file.path(), std::filesystem::file_size(file.path()) - 1); + } + + { + ClusterFile reader(file.path()); + ClusterVector clusters; + CHECK_THROWS_AS(reader.read_frame(clusters), std::runtime_error); + } + { + ClusterFile reader(file.path()); + CHECK_THROWS_AS(reader.read_frame(), std::runtime_error); + } +} + +TEST_CASE("ClusterFile close is idempotent", "[ClusterFile]") { + TemporaryClusterFile file; + auto expected = make_test_frame(42, 0.0); + + ClusterFile writer(file.path(), 1000, "w"); + writer.write_frame(expected); + writer.close(); + + CHECK_NOTHROW(writer.close()); + CHECK_THROWS_AS(writer.tell(), std::runtime_error); + CHECK_THROWS_AS(writer.write_frame(expected), std::runtime_error); + + ClusterFile reader(file.path()); + auto actual = reader.read_frame(); + REQUIRE(actual); + check_frame(*actual, expected); +} + +TEST_CASE("ClusterFile chunk reads reject incomplete frames", "[ClusterFile]") { + const auto use_roi = GENERATE(false, true); + CAPTURE(use_roi); + TemporaryClusterFile file; + auto expected = make_test_frame(42, 0.0); + { + ClusterFile writer(file.path(), 1000, "w"); + writer.write_frame(expected); + } + + SECTION("incomplete frame number") { + std::filesystem::resize_file(file.path(), sizeof(int32_t) - 1); + } + SECTION("missing cluster count") { + std::filesystem::resize_file(file.path(), sizeof(int32_t)); + } + SECTION("incomplete cluster count") { + std::filesystem::resize_file(file.path(), + sizeof(int32_t) + sizeof(uint32_t) - 1); + } + SECTION("missing cluster record") { + std::filesystem::resize_file(file.path(), sizeof(int32_t) + + sizeof(uint32_t) + + sizeof(TestCluster)); + } + SECTION("incomplete cluster record") { + std::filesystem::resize_file( + file.path(), std::filesystem::file_size(file.path()) - 1); + } + + ClusterFile reader(file.path()); + if (use_roi) { + reader.set_roi(aare::ROI{0, 10, 0, 10}); + } + CHECK_THROWS_AS(reader.read_clusters(10), std::runtime_error); +} + +TEST_CASE("ClusterFile chunk reads reject truncated leftover clusters", + "[ClusterFile]") { + const auto use_roi = GENERATE(false, true); + const auto requested = GENERATE(1, 10); + CAPTURE(use_roi, requested); + TemporaryClusterFile file; + auto expected = make_test_frame(42, 0.0); + { + ClusterFile writer(file.path(), 1000, "w"); + writer.write_frame(expected); + } + std::filesystem::resize_file(file.path(), + std::filesystem::file_size(file.path()) - 1); + + ClusterFile reader(file.path()); + if (use_roi) { + reader.set_roi(aare::ROI{0, 10, 0, 10}); + } + auto first = reader.read_clusters(1); + REQUIRE(first.size() == 1); + CHECK(first[0].data == expected[0].data); + CHECK_THROWS_AS(reader.read_clusters(requested), std::runtime_error); +} + +TEST_CASE("ClusterFile chunk reads preserve valid frame traversal", + "[ClusterFile]") { + const auto use_roi = GENERATE(false, true); + CAPTURE(use_roi); + TemporaryClusterFile file; + auto first = make_test_frame(42, 0.0); + auto second = make_test_frame(43, 100.0); + ClusterVector empty(0); + { + ClusterFile writer(file.path(), 1000, "w"); + writer.write_frame(empty); + writer.write_frame(first); + writer.write_frame(empty); + writer.write_frame(second); + writer.write_frame(empty); + } + + ClusterFile reader(file.path()); + if (use_roi) { + reader.set_roi(aare::ROI{0, 10, 0, 10}); + } + CHECK(reader.read_clusters(0).empty()); + CHECK(reader.tell() == 0); + + auto initial = reader.read_clusters(1); + REQUIRE(initial.size() == 1); + CHECK(initial[0].data == first[0].data); + const auto position = reader.tell(); + CHECK(reader.read_clusters(0).empty()); + CHECK(reader.tell() == position); + CHECK_THROWS_AS(reader.read_frame(), std::runtime_error); + + auto crossing = reader.read_clusters(2); + REQUIRE(crossing.size() == 2); + CHECK(crossing[0].data == first[1].data); + CHECK(crossing[1].data == second[0].data); + + auto final = reader.read_clusters(10); + REQUIRE(final.size() == 1); + CHECK(final[0].data == second[1].data); + CHECK(static_cast(reader.tell()) == + std::filesystem::file_size(file.path())); + CHECK(reader.read_clusters(10).empty()); + CHECK_FALSE(reader.read_frame()); +} + +TEST_CASE("ClusterFile frames preserve empty frames and reuse storage", + "[ClusterFile]") { + TemporaryClusterFile file; + auto first = make_test_frame(-42, 0.0); + auto second = make_test_frame(105, 100.0); + { + ClusterFile writer(file.path(), 1000, "w"); + writer.write_frame(ClusterVector(0, -43)); + writer.write_frame(first); + writer.write_frame(ClusterVector(0, 0)); + writer.write_frame(second); + } + + ClusterFile reader(file.path()); + auto frames = reader.frames(); + CHECK(reader.tell() == 0); + const TestCluster *storage = nullptr; + std::vector frame_numbers; + for (auto &frame : frames) { + frame_numbers.push_back(frame.frame_number()); + if (frame.frame_number() == -42) { + check_frame(frame, first); + storage = frame.data(); + } else if (frame.frame_number() == 105) { + check_frame(frame, second); + CHECK(frame.data() == storage); + } else { + CHECK(frame.empty()); + } + } + CHECK(frame_numbers == std::vector{-43, -42, 0, 105}); + CHECK(frames.begin() == frames.end()); +} + +TEST_CASE("ClusterFile chunks preserve cluster order across frame boundaries", + "[ClusterFile]") { + const auto chunk_size = GENERATE(1, 2, 3, 10); + TemporaryClusterFile file; + auto first = make_test_frame(42, 0.0); + auto second = make_test_frame(43, 100.0); + { + ClusterFile writer(file.path(), 1000, "w"); + writer.write_frame(ClusterVector(0)); + writer.write_frame(first); + writer.write_frame(ClusterVector(0)); + writer.write_frame(second); + writer.write_frame(ClusterVector(0)); + } + + ClusterFile reader(file.path(), chunk_size); + auto chunks = reader.chunks(); + auto requested = static_cast(chunk_size); + SECTION("constructor size") {} + SECTION("explicit size overrides the constructor") { + chunks = reader.chunks(3); + requested = 3; + CHECK(reader.chunk_size() == static_cast(chunk_size)); + } + CHECK(reader.tell() == 0); + size_t count = 0; + for (auto &chunk : chunks) { + REQUIRE(chunk.size() == std::min(requested, 4 - count)); + for (const auto &cluster : chunk) { + REQUIRE(count < 4); + const auto &expected = count < 2 ? first[count] : second[count - 2]; + CHECK(cluster.x == expected.x); + CHECK(cluster.y == expected.y); + CHECK(cluster.data == expected.data); + ++count; + } + } + CHECK(count == 4); + CHECK(chunks.begin() == chunks.end()); +} + +TEST_CASE("ClusterFile iterator results can be moved out without copying", + "[ClusterFile]") { + TemporaryClusterFile file; + auto first = make_test_frame(42, 0.0); + auto second = make_test_frame(43, 100.0); + { + ClusterFile writer(file.path(), 1000, "w"); + writer.write_frame(first); + writer.write_frame(second); + } + ClusterFile reader(file.path(), 2); + auto check_moves = [&](auto range) { + auto it = range.begin(); + REQUIRE(it != range.end()); + CHECK(range.end() != it); + const auto storage = it->data(); + auto retained = std::move(*it); + CHECK(retained.data() == storage); + ++it; + REQUIRE(it != range.end()); + check_frame(*it, second); + check_frame(retained, first); + it++; + CHECK(it == range.end()); + CHECK(range.end() == it); + ++it; + CHECK(it == range.end()); + reader.close(); + check_frame(retained, first); + }; + SECTION("frames") { check_moves(reader.frames()); } + SECTION("chunks") { check_moves(reader.chunks()); } +} + +TEST_CASE("ClusterFile iteration preserves filtering and gain correction", + "[ClusterFile]") { + const auto use_roi = GENERATE(false, true); + const auto use_noise = GENERATE(false, true); + TemporaryClusterFile file; + { + ClusterFile writer(file.path(), 1000, "w"); + writer.write_frame(make_test_frame(42, 0.0)); + writer.write_frame(make_test_frame(43, 100.0)); + } + ClusterFile reader(file.path(), 3); + ClusterFile reference(file.path()); + aare::NDArray gain({12, 12}, 2.0); + aare::NDArray noise({12, 12}, 5); + for (auto *file_reader : {&reader, &reference}) { + file_reader->set_gain_map(gain.view()); + if (use_roi) { + file_reader->set_roi(aare::ROI{0, 6, 0, 12}); + } + if (use_noise) { + file_reader->set_noise_map(noise.view()); + } + } + SECTION("frames include fully rejected frames") { + size_t count = 0; + for (auto &frame : reader.frames()) { + auto expected = reference.read_frame(); + REQUIRE(expected); + check_frame(frame, *expected); + ++count; + } + CHECK(count == 2); + CHECK_FALSE(reference.read_frame()); + } + SECTION("chunks count selected clusters") { + for (auto &chunk : reader.chunks()) { + auto expected = reference.read_clusters(3); + check_frame(chunk, expected); + } + CHECK(reference.read_clusters(3).empty()); + } +} + +TEST_CASE("ClusterFile iteration resumes after an early exit", + "[ClusterFile]") { + TemporaryClusterFile file; + auto first = make_test_frame(42, 0.0); + auto second = make_test_frame(43, 100.0); + { + ClusterFile writer(file.path(), 1000, "w"); + writer.write_frame(first); + writer.write_frame(second); + } + ClusterFile reader(file.path()); + SECTION("breaking a frame loop does not consume the following frame") { + for (auto &frame : reader.frames()) { + check_frame(frame, first); + break; + } + auto next = reader.read_frame(); + REQUIRE(next); + check_frame(*next, second); + } + SECTION("partial chunks must be completed before frame iteration") { + for (auto &chunk : reader.chunks(1)) { + REQUIRE(chunk.size() == 1); + CHECK(chunk[0].data == first[0].data); + break; + } + CHECK_THROWS_AS(reader.frames().begin(), std::runtime_error); + auto remainder = reader.read_clusters(1); + REQUIRE(remainder.size() == 1); + CHECK(remainder[0].data == first[1].data); + auto it = reader.frames().begin(); + REQUIRE(it != reader.frames().end()); + check_frame(*it, second); + } +} + +TEST_CASE("ClusterFile iterators distinguish EOF from read errors", + "[ClusterFile]") { + TemporaryClusterFile file; + { + ClusterFile writer(file.path(), 1000, "w"); + writer.write_frame(make_test_frame(42, 0.0)); + } + auto check_error = [&] { + ClusterFile frame_reader(file.path()); + ClusterFile chunk_reader(file.path()); + CHECK_THROWS_AS(frame_reader.frames().begin(), std::runtime_error); + CHECK_THROWS_AS(chunk_reader.chunks().begin(), std::runtime_error); + }; + SECTION("empty file") { + std::filesystem::resize_file(file.path(), 0); + ClusterFile reader(file.path()); + CHECK(reader.frames().begin() == reader.frames().end()); + CHECK(reader.chunks().begin() == reader.chunks().end()); + } + SECTION("partial header") { + std::filesystem::resize_file(file.path(), 7); + check_error(); + } + SECTION("partial record") { + std::filesystem::resize_file( + file.path(), std::filesystem::file_size(file.path()) - 1); + check_error(); + } + SECTION("partial record encountered on increment") { + std::filesystem::resize_file( + file.path(), std::filesystem::file_size(file.path()) - 1); + ClusterFile reader(file.path()); + auto it = reader.chunks(1).begin(); + REQUIRE(it != reader.chunks(1).end()); + CHECK_THROWS_AS(++it, std::runtime_error); + } + SECTION("closed before first read") { + ClusterFile reader(file.path()); + auto frames = reader.frames(); + auto chunks = reader.chunks(); + reader.close(); + CHECK_THROWS_AS(frames.begin(), std::runtime_error); + CHECK_THROWS_AS(chunks.begin(), std::runtime_error); + } + SECTION("closed during traversal") { + ClusterFile reader(file.path()); + auto it = reader.frames().begin(); + reader.close(); + CHECK_THROWS_AS(++it, std::runtime_error); + } + SECTION("write and append modes") { + const auto mode = GENERATE("w", "a"); + ClusterFile writer(file.path(), 1000, mode); + CHECK_THROWS_AS(writer.frames().begin(), std::runtime_error); + CHECK_THROWS_AS(writer.chunks().begin(), std::runtime_error); + } + SECTION("zero chunk sizes are rejected without reading") { + ClusterFile reader(file.path(), 0); + CHECK_THROWS_AS(reader.chunks(), std::invalid_argument); + CHECK_THROWS_AS(reader.chunks(0), std::invalid_argument); + CHECK(reader.tell() == 0); + CHECK(reader.read_clusters(0).empty()); + CHECK(reader.frames().begin() != reader.frames().end()); + } +} + +TEST_CASE("Read frame and modify cluster data", "[.with-data]") { auto fpath = test_data_path() / "clust" / "single_frame_97_clustrers.clust"; REQUIRE(std::filesystem::exists(fpath)); ClusterFile> f(fpath); auto clusters = f.read_frame(); - CHECK(clusters.size() == 97); - CHECK(clusters.frame_number() == 135); + REQUIRE(clusters); + CHECK(clusters->size() == 97); + CHECK(clusters->frame_number() == 135); int32_t expected_cluster_data[] = {0, 1, 2, 3, 4, 5, 6, 7, 8}; - clusters.push_back( + clusters->push_back( Cluster{0, 0, {0, 1, 2, 3, 4, 5, 6, 7, 8}}); - CHECK(clusters.size() == 98); - CHECK(clusters[0].x == 1); - CHECK(clusters[0].y == 200); + CHECK(clusters->size() == 98); + CHECK((*clusters)[0].x == 1); + CHECK((*clusters)[0].y == 200); - CHECK(std::equal(std::begin(clusters[0].data), std::end(clusters[0].data), + CHECK(std::equal(std::begin((*clusters)[0].data), + std::end((*clusters)[0].data), std::begin(expected_cluster_data))); } diff --git a/src/ClusterVector.test.cpp b/src/ClusterVector.test.cpp index b4e3abf6..135719c1 100644 --- a/src/ClusterVector.test.cpp +++ b/src/ClusterVector.test.cpp @@ -1,5 +1,8 @@ // SPDX-License-Identifier: MPL-2.0 #include "aare/ClusterVector.hpp" +#include "aare/GainMap.hpp" +#include "aare/NDArray.hpp" +#include #include #include @@ -374,3 +377,73 @@ TEST_CASE("Gain Map Calculation Index Map") { CHECK(index_map_x == clustertestdata.index_map_x); CHECK(index_map_y == clustertestdata.index_map_y); } + +namespace { + +template +void check_gain_map_cluster_bounds() { + using ClusterType = Cluster; + + constexpr ssize_t rows = 16; + constexpr ssize_t cols = 16; + constexpr uint16_t left = ClusterSizeX / 2; + constexpr uint16_t right = ClusterSizeX - left - 1; + constexpr uint16_t top = ClusterSizeY / 2; + constexpr uint16_t bottom = ClusterSizeY - top - 1; + + aare::NDArray gain_map({rows, cols}, 2.0); + aare::InvertedGainMap inverted_gain_map(gain_map); + ClusterVector clusters(6); + + const auto add_cluster = [&clusters](uint16_t x, uint16_t y) { + ClusterType cluster{}; + cluster.x = x; + cluster.y = y; + cluster.data.fill(2.0); + clusters.push_back(cluster); + }; + + add_cluster(left, top); + add_cluster(cols - right - 1, rows - bottom - 1); + add_cluster(left - 1, top); + add_cluster(left, top - 1); + add_cluster(cols - right, top); + add_cluster(left, rows - bottom); + + inverted_gain_map.apply_gain_map(clusters); + + for (size_t i = 0; i < 2; ++i) { + CHECK(std::all_of(clusters[i].data.begin(), clusters[i].data.end(), + [](double value) { return value == 1.0; })); + } + for (size_t i = 2; i < clusters.size(); ++i) { + CHECK(std::all_of(clusters[i].data.begin(), clusters[i].data.end(), + [](double value) { return value == 0.0; })); + } + + aare::NDArray small_gain_map( + {ClusterSizeY - 1, ClusterSizeX - 1}, 2.0); + aare::InvertedGainMap small_inverted_gain_map(small_gain_map); + ClusterVector oversized_cluster(1); + ClusterType cluster{}; + cluster.x = left; + cluster.y = top; + cluster.data.fill(2.0); + oversized_cluster.push_back(cluster); + + small_inverted_gain_map.apply_gain_map(oversized_cluster); + + CHECK(std::all_of(oversized_cluster[0].data.begin(), + oversized_cluster[0].data.end(), + [](double value) { return value == 0.0; })); +} + +} // namespace + +TEST_CASE("Gain map bounds cover the full cluster footprint", "[GainMap]") { + SECTION("3x3") { check_gain_map_cluster_bounds<3, 3>(); } + SECTION("5x5") { check_gain_map_cluster_bounds<5, 5>(); } + SECTION("7x7") { check_gain_map_cluster_bounds<7, 7>(); } + SECTION("9x9") { check_gain_map_cluster_bounds<9, 9>(); } + SECTION("5x7") { check_gain_map_cluster_bounds<5, 7>(); } +} diff --git a/src/FilePtr.cpp b/src/FilePtr.cpp index dd907bda..45e65d46 100644 --- a/src/FilePtr.cpp +++ b/src/FilePtr.cpp @@ -22,6 +22,8 @@ FilePtr &FilePtr::operator=(FilePtr &&other) { return *this; } +FilePtr::operator bool() const noexcept { return fp_ != nullptr; } + FILE *FilePtr::get() { return fp_; } ssize_t FilePtr::tell() { diff --git a/src/FilePtr.test.cpp b/src/FilePtr.test.cpp new file mode 100644 index 00000000..6b23e62a --- /dev/null +++ b/src/FilePtr.test.cpp @@ -0,0 +1,18 @@ +// SPDX-License-Identifier: MPL-2.0 + +#include "aare/FilePtr.hpp" + +#include +#include + +TEST_CASE("FilePtr converts to its open state", "[FilePtr]") { + aare::FilePtr empty; + CHECK_FALSE(empty); + + aare::FilePtr open(__FILE__, "rb"); + CHECK(open); + + aare::FilePtr moved(std::move(open)); + CHECK(moved); + CHECK_FALSE(open); +}