diff --git a/.github/workflows/build_and_deploy_conda.yml b/.github/workflows/build_and_deploy_conda.yml index 65483c3..8917419 100644 --- a/.github/workflows/build_and_deploy_conda.yml +++ b/.github/workflows/build_and_deploy_conda.yml @@ -1,9 +1,9 @@ name: Build pkgs and deploy if on main on: - push: - branches: - - main + release: + types: + - published jobs: build: diff --git a/.github/workflows/build_docs.yml b/.github/workflows/build_docs.yml index 24050a3..4fd23e7 100644 --- a/.github/workflows/build_docs.yml +++ b/.github/workflows/build_docs.yml @@ -2,7 +2,10 @@ name: Build the package using cmake then documentation on: workflow_dispatch: - push: + pull_request: + release: + types: + - published permissions: @@ -40,7 +43,7 @@ jobs: run: | mkdir build cd build - cmake .. -DAARE_SYSTEM_LIBRARIES=ON -DAARE_DOCS=ON + cmake .. -DAARE_SYSTEM_LIBRARIES=ON -DAARE_PYTHON_BINDINGS=ON -DAARE_DOCS=ON make -j 2 make docs @@ -55,7 +58,7 @@ jobs: url: ${{ steps.deployment.outputs.page_url }} runs-on: ubuntu-latest needs: build - if: github.ref == 'refs/heads/main' + if: (github.event_name == 'release' && github.event.action == 'published') || (github.event_name == 'workflow_dispatch' ) steps: - name: Deploy to GitHub Pages id: deployment diff --git a/CMakeLists.txt b/CMakeLists.txt index 64c08ad..191448c 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -79,6 +79,9 @@ endif() if(AARE_VERBOSE) add_compile_definitions(AARE_VERBOSE) + add_compile_definitions(AARE_LOG_LEVEL=aare::logDEBUG5) +else() + add_compile_definitions(AARE_LOG_LEVEL=aare::logERROR) endif() if(AARE_CUSTOM_ASSERT) @@ -90,6 +93,7 @@ if(AARE_BENCHMARKS) endif() + set(CMAKE_EXPORT_COMPILE_COMMANDS ON) if(AARE_FETCH_LMFIT) @@ -446,6 +450,7 @@ if(AARE_TESTS) ${CMAKE_CURRENT_SOURCE_DIR}/src/NumpyFile.test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/NumpyHelpers.test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/RawFile.test.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/src/RawSubFile.test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/utils/task.test.cpp ) diff --git a/RELEASE.md b/RELEASE.md new file mode 100644 index 0000000..afda148 --- /dev/null +++ b/RELEASE.md @@ -0,0 +1,22 @@ +# Release notes + + +### head + +Features: + +- Cluster finder now works with 5x5, 7x7 and 9x9 clusters + + +### 2025.05.22 + +Features: + +- Added scurve fitting + +Bugfixes: + +- Fixed crash when opening raw files with large number of data files + + + diff --git a/VERSION b/VERSION index bd52db8..ae365e4 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -0.0.0 \ No newline at end of file +2025.5.22 \ No newline at end of file diff --git a/benchmarks/calculateeta_benchmark.cpp b/benchmarks/calculateeta_benchmark.cpp index a320188..6c40c5c 100644 --- a/benchmarks/calculateeta_benchmark.cpp +++ b/benchmarks/calculateeta_benchmark.cpp @@ -41,8 +41,8 @@ BENCHMARK_F(ClusterFixture, Calculate2x2Eta)(benchmark::State &st) { } // almost takes double the time -BENCHMARK_F(ClusterFixture, - CalculateGeneralEtaFor2x2Cluster)(benchmark::State &st) { +BENCHMARK_F(ClusterFixture, CalculateGeneralEtaFor2x2Cluster) +(benchmark::State &st) { for (auto _ : st) { // This code gets timed Eta2 eta = calculate_eta2(cluster_2x2); @@ -59,8 +59,8 @@ BENCHMARK_F(ClusterFixture, Calculate3x3Eta)(benchmark::State &st) { } // almost takes double the time -BENCHMARK_F(ClusterFixture, - CalculateGeneralEtaFor3x3Cluster)(benchmark::State &st) { +BENCHMARK_F(ClusterFixture, CalculateGeneralEtaFor3x3Cluster) +(benchmark::State &st) { for (auto _ : st) { // This code gets timed Eta2 eta = calculate_eta2(cluster_3x3); diff --git a/benchmarks/ndarray_benchmark.cpp b/benchmarks/ndarray_benchmark.cpp index 55fa263..91a2d9b 100644 --- a/benchmarks/ndarray_benchmark.cpp +++ b/benchmarks/ndarray_benchmark.cpp @@ -1,136 +1,132 @@ -#include #include "aare/NDArray.hpp" - +#include using aare::NDArray; constexpr ssize_t size = 1024; class TwoArrays : public benchmark::Fixture { -public: - NDArray a{{size,size},0}; - NDArray b{{size,size},0}; - void SetUp(::benchmark::State& state) { - for(uint32_t i = 0; i < size; i++){ - for(uint32_t j = 0; j < size; j++){ - a(i, j)= i*j+1; - b(i, j)= i*j+1; - } + public: + NDArray a{{size, size}, 0}; + NDArray b{{size, size}, 0}; + void SetUp(::benchmark::State &state) { + for (uint32_t i = 0; i < size; i++) { + for (uint32_t j = 0; j < size; j++) { + a(i, j) = i * j + 1; + b(i, j) = i * j + 1; + } + } } - } - // void TearDown(::benchmark::State& state) { - // } + // void TearDown(::benchmark::State& state) { + // } }; - - - -BENCHMARK_F(TwoArrays, AddWithOperator)(benchmark::State& st) { - for (auto _ : st) { - // This code gets timed - NDArray res = a+b; - benchmark::DoNotOptimize(res); - } -} -BENCHMARK_F(TwoArrays, AddWithIndex)(benchmark::State& st) { - for (auto _ : st) { - // This code gets timed - NDArray res(a.shape()); - for (uint32_t i = 0; i < a.size(); i++) { - res(i) = a(i) + b(i); +BENCHMARK_F(TwoArrays, AddWithOperator)(benchmark::State &st) { + for (auto _ : st) { + // This code gets timed + NDArray res = a + b; + benchmark::DoNotOptimize(res); + } +} +BENCHMARK_F(TwoArrays, AddWithIndex)(benchmark::State &st) { + for (auto _ : st) { + // This code gets timed + NDArray res(a.shape()); + for (uint32_t i = 0; i < a.size(); i++) { + res(i) = a(i) + b(i); + } + benchmark::DoNotOptimize(res); } - benchmark::DoNotOptimize(res); - } } -BENCHMARK_F(TwoArrays, SubtractWithOperator)(benchmark::State& st) { - for (auto _ : st) { - // This code gets timed - NDArray res = a-b; - benchmark::DoNotOptimize(res); - } -} -BENCHMARK_F(TwoArrays, SubtractWithIndex)(benchmark::State& st) { - for (auto _ : st) { - // This code gets timed - NDArray res(a.shape()); - for (uint32_t i = 0; i < a.size(); i++) { - res(i) = a(i) - b(i); +BENCHMARK_F(TwoArrays, SubtractWithOperator)(benchmark::State &st) { + for (auto _ : st) { + // This code gets timed + NDArray res = a - b; + benchmark::DoNotOptimize(res); + } +} +BENCHMARK_F(TwoArrays, SubtractWithIndex)(benchmark::State &st) { + for (auto _ : st) { + // This code gets timed + NDArray res(a.shape()); + for (uint32_t i = 0; i < a.size(); i++) { + res(i) = a(i) - b(i); + } + benchmark::DoNotOptimize(res); } - benchmark::DoNotOptimize(res); - } } -BENCHMARK_F(TwoArrays, MultiplyWithOperator)(benchmark::State& st) { - for (auto _ : st) { - // This code gets timed - NDArray res = a*b; - benchmark::DoNotOptimize(res); - } -} -BENCHMARK_F(TwoArrays, MultiplyWithIndex)(benchmark::State& st) { - for (auto _ : st) { - // This code gets timed - NDArray res(a.shape()); - for (uint32_t i = 0; i < a.size(); i++) { - res(i) = a(i) * b(i); +BENCHMARK_F(TwoArrays, MultiplyWithOperator)(benchmark::State &st) { + for (auto _ : st) { + // This code gets timed + NDArray res = a * b; + benchmark::DoNotOptimize(res); + } +} +BENCHMARK_F(TwoArrays, MultiplyWithIndex)(benchmark::State &st) { + for (auto _ : st) { + // This code gets timed + NDArray res(a.shape()); + for (uint32_t i = 0; i < a.size(); i++) { + res(i) = a(i) * b(i); + } + benchmark::DoNotOptimize(res); } - benchmark::DoNotOptimize(res); - } } -BENCHMARK_F(TwoArrays, DivideWithOperator)(benchmark::State& st) { - for (auto _ : st) { - // This code gets timed - NDArray res = a/b; - benchmark::DoNotOptimize(res); - } -} -BENCHMARK_F(TwoArrays, DivideWithIndex)(benchmark::State& st) { - for (auto _ : st) { - // This code gets timed - NDArray res(a.shape()); - for (uint32_t i = 0; i < a.size(); i++) { - res(i) = a(i) / b(i); +BENCHMARK_F(TwoArrays, DivideWithOperator)(benchmark::State &st) { + for (auto _ : st) { + // This code gets timed + NDArray res = a / b; + benchmark::DoNotOptimize(res); + } +} +BENCHMARK_F(TwoArrays, DivideWithIndex)(benchmark::State &st) { + for (auto _ : st) { + // This code gets timed + NDArray res(a.shape()); + for (uint32_t i = 0; i < a.size(); i++) { + res(i) = a(i) / b(i); + } + benchmark::DoNotOptimize(res); } - benchmark::DoNotOptimize(res); - } } -BENCHMARK_F(TwoArrays, FourAddWithOperator)(benchmark::State& st) { - for (auto _ : st) { - // This code gets timed - NDArray res = a+b+a+b; - benchmark::DoNotOptimize(res); - } -} -BENCHMARK_F(TwoArrays, FourAddWithIndex)(benchmark::State& st) { - for (auto _ : st) { - // This code gets timed - NDArray res(a.shape()); - for (uint32_t i = 0; i < a.size(); i++) { - res(i) = a(i) + b(i) + a(i) + b(i); +BENCHMARK_F(TwoArrays, FourAddWithOperator)(benchmark::State &st) { + for (auto _ : st) { + // This code gets timed + NDArray res = a + b + a + b; + benchmark::DoNotOptimize(res); + } +} +BENCHMARK_F(TwoArrays, FourAddWithIndex)(benchmark::State &st) { + for (auto _ : st) { + // This code gets timed + NDArray res(a.shape()); + for (uint32_t i = 0; i < a.size(); i++) { + res(i) = a(i) + b(i) + a(i) + b(i); + } + benchmark::DoNotOptimize(res); } - benchmark::DoNotOptimize(res); - } } -BENCHMARK_F(TwoArrays, MultiplyAddDivideWithOperator)(benchmark::State& st) { - for (auto _ : st) { - // This code gets timed - NDArray res = a*a+b/a; - benchmark::DoNotOptimize(res); - } -} -BENCHMARK_F(TwoArrays, MultiplyAddDivideWithIndex)(benchmark::State& st) { - for (auto _ : st) { - // This code gets timed - NDArray res(a.shape()); - for (uint32_t i = 0; i < a.size(); i++) { - res(i) = a(i) * a(i) + b(i) / a(i); +BENCHMARK_F(TwoArrays, MultiplyAddDivideWithOperator)(benchmark::State &st) { + for (auto _ : st) { + // This code gets timed + NDArray res = a * a + b / a; + benchmark::DoNotOptimize(res); + } +} +BENCHMARK_F(TwoArrays, MultiplyAddDivideWithIndex)(benchmark::State &st) { + for (auto _ : st) { + // This code gets timed + NDArray res(a.shape()); + for (uint32_t i = 0; i < a.size(); i++) { + res(i) = a(i) * a(i) + b(i) / a(i); + } + benchmark::DoNotOptimize(res); } - benchmark::DoNotOptimize(res); - } } BENCHMARK_MAIN(); \ No newline at end of file diff --git a/docs/src/ClusterFile.rst b/docs/src/ClusterFile.rst index 79de086..a2ee162 100644 --- a/docs/src/ClusterFile.rst +++ b/docs/src/ClusterFile.rst @@ -4,4 +4,5 @@ ClusterFile .. doxygenclass:: aare::ClusterFile :members: :undoc-members: - :private-members: \ No newline at end of file + :private-members: + diff --git a/docs/src/Philosophy.rst b/docs/src/Philosophy.rst new file mode 100644 index 0000000..f187bad --- /dev/null +++ b/docs/src/Philosophy.rst @@ -0,0 +1,47 @@ +**************** +Philosophy +**************** + + +Fast code with a simple interface +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +Aare should be fast and efficient, but also easy to use. We strive to keep a simple interface that feels intuitive. +Internally we use C++ for performance and the ability to integrate the library in other programs, but we see most +users using the Python interface. + +Live at head +~~~~~~~~~~~~~~~~~~ + +As a user of the library you should be able to, and is expected to, use the latest version. Bug fixes will rarely be backported +to older releases. By upgrading frequently you will benefit from the latest features and minimize the effort to maintain your scripts/code +by doing several small upgrades instead of one big upgrade. + +API +~~~~~~~~~~~~~~~~~~ + +We aim to keep the API stable and only break it for good reasons. But specially now in the early stages of development +the API will change. On those occasions it will be clearly stated in the release notes. However, the norm should be a +backward compatible API. + +Documentation +~~~~~~~~~~~~~~~~~~ + +Being a library it is important to have a well documented API. We use Doxygen to generate the C++ documentation +and Sphinx for the Python part. Breathe is used to integrate the two into one Sphinx html site. The documentation is built +automatically on release by the CI and published to GitHub pages. In addition to the generated API documentation, +certain classes might need more descriptions of the usage. This is then placed in the .rst files in the docs/src directory. + +.. attention:: + + The code should be well documented, but using descriptive names is more important. In the same spirit + if a function is called `getNumberOfFrames()` you don't need to write a comment saying that it gets the + number of frames. + + +Dependencies +~~~~~~~~~~~~~~~~~~ + +Deployment in the scientific community is often tricky. Either due to old OS versions or the lack of package managers. +We strive to keep the dependencies to a minimum and will vendor some libraries to simplify deployment even though it comes +at a cost of build time. \ No newline at end of file diff --git a/docs/src/Requirements.rst b/docs/src/Requirements.rst index c962f73..b0c370f 100644 --- a/docs/src/Requirements.rst +++ b/docs/src/Requirements.rst @@ -2,18 +2,21 @@ Requirements ============================================== - C++17 compiler (gcc 8/clang 7) -- CMake 3.14+ +- CMake 3.15+ **Internally used libraries** .. note :: - These can also be picked up from the system/conda environment by specifying: + To save compile time some of the dependencies can also be picked up from the system/conda environment by specifying: -DAARE_SYSTEM_LIBRARIES=ON during the cmake configuration. -- pybind11 +To simplify deployment we build and statically link a few libraries. + - fmt +- lmfit - https://jugit.fz-juelich.de/mlz/lmfit - nlohmann_json +- pybind11 - ZeroMQ **Extra dependencies for building documentation** diff --git a/docs/src/Workflow.rst b/docs/src/Workflow.rst new file mode 100644 index 0000000..9b41e74 --- /dev/null +++ b/docs/src/Workflow.rst @@ -0,0 +1,86 @@ +**************** +Workflow +**************** + +This page describes how we develop aare. + +GitHub centric +~~~~~~~~~~~~~~~~~~ + +We use GitHub for all development. Issues and pull requests provide a platform for collaboration as well +as a record of the development process. Even if we discuss things in person, we record the outcome in an issue. +If a particular implementation is chosen over another, the reason should be recorded in the pull request. + + +Branches +~~~~~~~~~~~~~~~~~~ + +We aim for an as lightweight branching strategy as possible. Short-lived feature branches are merged back into main. +The main branch is expected to always be in a releasable state. A release is simply a tag on main which provides a +reference and triggers the CI to build the release artifacts (conda, pypi etc.). For large features consider merging +smaller chunks into main as they are completed, rather than waiting for the entire feature to be finished. Worst case +make sure your feature branch merges with main regularly to avoid large merge conflicts later on. + +.. note:: + + The main branch is expected to always work. Feel free to pull from main instead of sticking to a + release + + +Releases +~~~~~~~~~~~~~~~~~~ + +Release early, release often. As soon as "enough" new features have been implemented, a release is created. +A release should not be a big thing, rather a routine part of development that does not require any special person or +unfamiliar steps. + + + +Checklists for deployment +~~~~~~~~~~~~~~~~~~ + +**Feature:** + +#. Create a new issue for the feature (label feature) +#. Create a new branch from main. +#. Implement the feature including test and documentation +#. Add the feature to RELEASE.md under head +#. Create a pull request linked to the issue +#. Code is reviewed by at least one other person +#. Once approved, the branch is merged into main + + +**BugFix:** + +Essentially the same as for a feature, if possible start with +a failing test that demonstrates the bug. + +#. Create a new issue for the bug (label bug) +#. Create a new branch from main. +#. **Write a test that fails for the bug** +#. Implement the fix +#. **Run the test to ensure it passes** +#. Add the bugfix to RELEASE.md under head +#. Create a pull request linked to the issue. +#. Code is reviewed by at least one other person +#. Once approved, the branch is merged into main + +**Release:** + +#. Once "enough" new features have been implemented, a release is created +#. Update RELEASE.md with the tag of the release and verify that it is complete +#. Create the release in GitHub describing the new features and bug fixes +#. CI makes magic + + +**Update documentation only:** + +.. attention:: + + It's possible to update the documentation without changing the code, but take + care since the docs will reflect the code in main and not the latest release. + +#. Create a PR to main with the documentation changes +#. Create a pull request linked to the issue. +#. Code is reviewed by at least one other person +#. Once merged you can manually trigger the CI workflow for documentation \ No newline at end of file diff --git a/docs/src/index.rst b/docs/src/index.rst index af5e99a..a33b4df 100644 --- a/docs/src/index.rst +++ b/docs/src/index.rst @@ -63,4 +63,6 @@ AARE :caption: Developer :maxdepth: 3 + Philosophy + Workflow Tests \ No newline at end of file diff --git a/docs/src/pyClusterFile.rst b/docs/src/pyClusterFile.rst index bdf898c..cd391e0 100644 --- a/docs/src/pyClusterFile.rst +++ b/docs/src/pyClusterFile.rst @@ -2,9 +2,24 @@ ClusterFile ============ + +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. + .. py:currentmodule:: aare .. autoclass:: ClusterFile + :members: + :undoc-members: + :inherited-members: + + +Below is the API of the ClusterFile_Cluster3x3i but all variants share the same API. + +.. autoclass:: aare._aare.ClusterFile_Cluster3x3i + :special-members: __init__ :members: :undoc-members: :show-inheritance: diff --git a/docs/src/pyClusterVector.rst b/docs/src/pyClusterVector.rst index 4277920..ff115c9 100644 --- a/docs/src/pyClusterVector.rst +++ b/docs/src/pyClusterVector.rst @@ -2,8 +2,10 @@ ClusterVector ================ The ClusterVector, holds clusters from the ClusterFinder. Since it is templated -in C++ we use a suffix indicating the data type in python. The suffix is -``_i`` for integer, ``_f`` for float, and ``_d`` for double. +in C++ we use a suffix indicating the type of cluster it holds. The suffix follows +the same pattern as for ClusterFile i.e. ``ClusterVector_Cluster3x3i`` +for a vector holding 3x3 integer clusters. + At the moment the functionality from python is limited and it is not supported to push_back clusters to the vector. The intended use case is to pass it to @@ -26,7 +28,8 @@ C++ functions that support the ClusterVector or to view it as a numpy array. .. py:currentmodule:: aare -.. autoclass:: ClusterVector_i +.. autoclass:: aare._aare.ClusterVector_Cluster3x3i + :special-members: __init__ :members: :undoc-members: :show-inheritance: diff --git a/etc/dev-env.yml b/etc/dev-env.yml index e580c81..2edfc46 100644 --- a/etc/dev-env.yml +++ b/etc/dev-env.yml @@ -5,9 +5,12 @@ dependencies: - anaconda-client - conda-build - doxygen - - sphinx=7.1.2 + - sphinx - breathe - sphinx_rtd_theme - furo - zeromq + - pybind11 + - numpy + - matplotlib diff --git a/include/aare/ArrayExpr.hpp b/include/aare/ArrayExpr.hpp index d326601..e5fb5d7 100644 --- a/include/aare/ArrayExpr.hpp +++ b/include/aare/ArrayExpr.hpp @@ -1,10 +1,9 @@ #pragma once -#include -#include -#include -#include #include "aare/defs.hpp" - +#include +#include +#include +#include namespace aare { @@ -15,7 +14,9 @@ template class ArrayExpr { auto operator[](size_t i) const { return static_cast(*this)[i]; } auto operator()(size_t i) const { return static_cast(*this)[i]; } auto size() const { return static_cast(*this).size(); } - std::array shape() const { return static_cast(*this).shape(); } + std::array shape() const { + return static_cast(*this).shape(); + } }; template @@ -47,7 +48,7 @@ class ArraySub : public ArrayExpr, Ndim> { }; template -class ArrayMul : public ArrayExpr,Ndim> { +class ArrayMul : public ArrayExpr, Ndim> { const A &arr1_; const B &arr2_; @@ -74,15 +75,13 @@ class ArrayDiv : public ArrayExpr, Ndim> { std::array shape() const { return arr1_.shape(); } }; - - template auto operator+(const ArrayExpr &arr1, const ArrayExpr &arr2) { return ArrayAdd, ArrayExpr, Ndim>(arr1, arr2); } template -auto operator-(const ArrayExpr &arr1, const ArrayExpr &arr2) { +auto operator-(const ArrayExpr &arr1, const ArrayExpr &arr2) { return ArraySub, ArrayExpr, Ndim>(arr1, arr2); } @@ -96,6 +95,4 @@ auto operator/(const ArrayExpr &arr1, const ArrayExpr &arr2) { return ArrayDiv, ArrayExpr, Ndim>(arr1, arr2); } - - } // namespace aare \ No newline at end of file diff --git a/include/aare/CircularFifo.hpp b/include/aare/CircularFifo.hpp index 8098082..853a89f 100644 --- a/include/aare/CircularFifo.hpp +++ b/include/aare/CircularFifo.hpp @@ -17,7 +17,8 @@ template class CircularFifo { public: CircularFifo() : CircularFifo(100){}; - CircularFifo(uint32_t size) : fifo_size(size), free_slots(size + 1), filled_slots(size + 1) { + CircularFifo(uint32_t size) + : fifo_size(size), free_slots(size + 1), filled_slots(size + 1) { // TODO! how do we deal with alignment for writing? alignas??? // Do we give the user a chance to provide memory locations? @@ -55,7 +56,8 @@ template class CircularFifo { bool try_pop_free(ItemType &v) { return free_slots.read(v); } - ItemType pop_value(std::chrono::nanoseconds wait, std::atomic &stopped) { + ItemType pop_value(std::chrono::nanoseconds wait, + std::atomic &stopped) { ItemType v; while (!filled_slots.read(v) && !stopped) { std::this_thread::sleep_for(wait); diff --git a/include/aare/ClusterFile.hpp b/include/aare/ClusterFile.hpp index ef78874..c0eca33 100644 --- a/include/aare/ClusterFile.hpp +++ b/include/aare/ClusterFile.hpp @@ -5,6 +5,8 @@ #include "aare/GainMap.hpp" #include "aare/NDArray.hpp" #include "aare/defs.hpp" +#include "aare/logger.hpp" + #include #include #include @@ -369,11 +371,15 @@ ClusterFile::read_frame_without_cut() { "Could not read number of clusters"); } + 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) != static_cast(n_clusters)) { throw std::runtime_error(LOCATION + "Could not read clusters"); diff --git a/include/aare/ClusterFileSink.hpp b/include/aare/ClusterFileSink.hpp index 810e63c..1900774 100644 --- a/include/aare/ClusterFileSink.hpp +++ b/include/aare/ClusterFileSink.hpp @@ -21,7 +21,7 @@ class ClusterFileSink { void process() { m_stopped = false; - fmt::print("ClusterFileSink started\n"); + LOG(logDEBUG) << "ClusterFileSink started"; while (!m_stop_requested || !m_source->isEmpty()) { if (ClusterVector *clusters = m_source->frontPtr(); clusters != nullptr) { @@ -41,13 +41,16 @@ class ClusterFileSink { std::this_thread::sleep_for(m_default_wait); } } - fmt::print("ClusterFileSink stopped\n"); + LOG(logDEBUG) << "ClusterFileSink stopped"; m_stopped = true; } public: ClusterFileSink(ClusterFinderMT *source, const std::filesystem::path &fname) { + LOG(logDEBUG) << "ClusterFileSink: " + << "source: " << source->sink() + << ", file: " << fname.string(); m_source = source->sink(); m_thread = std::thread(&ClusterFileSink::process, this); m_file.open(fname, std::ios::binary); diff --git a/include/aare/ClusterFinder.hpp b/include/aare/ClusterFinder.hpp index ea11162..069d887 100644 --- a/include/aare/ClusterFinder.hpp +++ b/include/aare/ClusterFinder.hpp @@ -38,7 +38,11 @@ class ClusterFinder { : m_image_size(image_size), m_nSigma(nSigma), c2(sqrt((ClusterSizeY + 1) / 2 * (ClusterSizeX + 1) / 2)), c3(sqrt(ClusterSizeX * ClusterSizeY)), - m_pedestal(image_size[0], image_size[1]), m_clusters(capacity) {}; + m_pedestal(image_size[0], image_size[1]), m_clusters(capacity) { + LOG(logDEBUG) << "ClusterFinder: " + << "image_size: " << image_size[0] << "x" << image_size[1] + << ", nSigma: " << nSigma << ", capacity: " << capacity; + } void push_pedestal_frame(NDView frame) { m_pedestal.push(frame); diff --git a/include/aare/ClusterFinderMT.hpp b/include/aare/ClusterFinderMT.hpp index 2dfb279..0340973 100644 --- a/include/aare/ClusterFinderMT.hpp +++ b/include/aare/ClusterFinderMT.hpp @@ -8,6 +8,7 @@ #include "aare/ClusterFinder.hpp" #include "aare/NDArray.hpp" #include "aare/ProducerConsumerQueue.hpp" +#include "aare/logger.hpp" namespace aare { @@ -123,6 +124,12 @@ class ClusterFinderMT { size_t capacity = 2000, size_t n_threads = 3) : m_n_threads(n_threads) { + LOG(logDEBUG1) << "ClusterFinderMT: " + << "image_size: " << image_size[0] << "x" + << image_size[1] << ", nSigma: " << nSigma + << ", capacity: " << capacity + << ", n_threads: " << n_threads; + for (size_t i = 0; i < n_threads; i++) { m_cluster_finders.push_back( std::make_unique< diff --git a/include/aare/ClusterVector.hpp b/include/aare/ClusterVector.hpp index c8b1ea1..9d575d9 100644 --- a/include/aare/ClusterVector.hpp +++ b/include/aare/ClusterVector.hpp @@ -133,9 +133,9 @@ class ClusterVector> { */ size_t capacity() const { return m_data.capacity(); } - const auto begin() const { return m_data.begin(); } + auto begin() const { return m_data.begin(); } - const auto end() const { return m_data.end(); } + auto end() const { return m_data.end(); } /** * @brief Return the size in bytes of a single cluster diff --git a/include/aare/CtbRawFile.hpp b/include/aare/CtbRawFile.hpp index 68dab23..afae0a2 100644 --- a/include/aare/CtbRawFile.hpp +++ b/include/aare/CtbRawFile.hpp @@ -1,27 +1,27 @@ #pragma once #include "aare/FileInterface.hpp" -#include "aare/RawMasterFile.hpp" #include "aare/Frame.hpp" +#include "aare/RawMasterFile.hpp" #include #include -namespace aare{ +namespace aare { - -class CtbRawFile{ +class CtbRawFile { RawMasterFile m_master; std::ifstream m_file; size_t m_current_frame{0}; size_t m_current_subfile{0}; size_t m_num_subfiles{0}; -public: + + public: CtbRawFile(const std::filesystem::path &fname); - void read_into(std::byte *image_buf, DetectorHeader* header = nullptr); - void seek(size_t frame_index); //!< seek to the given frame index - size_t tell() const; //!< get the frame index of the file pointer + void read_into(std::byte *image_buf, DetectorHeader *header = nullptr); + void seek(size_t frame_index); //!< seek to the given frame index + size_t tell() const; //!< get the frame index of the file pointer // in the specific class we can expose more functionality @@ -29,13 +29,13 @@ public: size_t frames_in_file() const; RawMasterFile master() const; -private: + + private: void find_subfiles(); size_t sub_file_index(size_t frame_index) const { return frame_index / m_master.max_frames_per_file(); } void open_data_file(size_t subfile_index); - }; -} \ No newline at end of file +} // namespace aare \ No newline at end of file diff --git a/include/aare/Dtype.hpp b/include/aare/Dtype.hpp index 7e1e62a..7047264 100644 --- a/include/aare/Dtype.hpp +++ b/include/aare/Dtype.hpp @@ -6,31 +6,37 @@ namespace aare { -// The format descriptor is a single character that specifies the type of the data +// The format descriptor is a single character that specifies the type of the +// data // - python documentation: https://docs.python.org/3/c-api/arg.html#numbers -// - py::format_descriptor::format() (in pybind11) does not return the same format as +// - py::format_descriptor::format() (in pybind11) does not return the same +// format as // written in python.org documentation. -// - numpy also doesn't use the same format. and also numpy associates the format -// with variable bitdepth types. (e.g. long is int64 on linux64 and int32 on win64) -// https://numpy.org/doc/stable/reference/arrays.scalars.html +// - numpy also doesn't use the same format. and also numpy associates the +// format +// with variable bitdepth types. (e.g. long is int64 on linux64 and int32 on +// win64) https://numpy.org/doc/stable/reference/arrays.scalars.html // // github issue discussing this: // https://github.com/pybind/pybind11/issues/1908#issuecomment-658358767 // -// [IN LINUX] the difference is for int64 (long) and uint64 (unsigned long). The format -// descriptor is 'q' and 'Q' respectively and in the documentation it is 'l' and 'k'. +// [IN LINUX] the difference is for int64 (long) and uint64 (unsigned long). The +// format descriptor is 'q' and 'Q' respectively and in the documentation it is +// 'l' and 'k'. // in practice numpy doesn't seem to care when reading buffer info: the library // interprets 'q' or 'l' as int64 and 'Q' or 'L' as uint64. -// for this reason we decided to use the same format descriptor as pybind to avoid -// any further discrepancies. +// for this reason we decided to use the same format descriptor as pybind to +// avoid any further discrepancies. // in the following order: // int8, uint8, int16, uint16, int32, uint32, int64, uint64, float, double -const char DTYPE_FORMAT_DSC[] = {'b', 'B', 'h', 'H', 'i', 'I', 'q', 'Q', 'f', 'd'}; +const char DTYPE_FORMAT_DSC[] = {'b', 'B', 'h', 'H', 'i', + 'I', 'q', 'Q', 'f', 'd'}; // on linux64 & apple -const char NUMPY_FORMAT_DSC[] = {'b', 'B', 'h', 'H', 'i', 'I', 'l', 'L', 'f', 'd'}; +const char NUMPY_FORMAT_DSC[] = {'b', 'B', 'h', 'H', 'i', + 'I', 'l', 'L', 'f', 'd'}; /** * @brief enum class to define the endianess of the system */ @@ -52,12 +58,29 @@ enum class endian { */ class Dtype { public: - enum TypeIndex { INT8, UINT8, INT16, UINT16, INT32, UINT32, INT64, UINT64, FLOAT, DOUBLE, ERROR, NONE }; + enum TypeIndex { + INT8, + UINT8, + INT16, + UINT16, + INT32, + UINT32, + INT64, + UINT64, + FLOAT, + DOUBLE, + ERROR, + NONE + }; uint8_t bitdepth() const; size_t bytes() const; - std::string format_descr() const { return std::string(1, DTYPE_FORMAT_DSC[static_cast(m_type)]); } - std::string numpy_descr() const { return std::string(1, NUMPY_FORMAT_DSC[static_cast(m_type)]); } + std::string format_descr() const { + return std::string(1, DTYPE_FORMAT_DSC[static_cast(m_type)]); + } + std::string numpy_descr() const { + return std::string(1, NUMPY_FORMAT_DSC[static_cast(m_type)]); + } explicit Dtype(const std::type_info &t); explicit Dtype(std::string_view sv); diff --git a/include/aare/File.hpp b/include/aare/File.hpp index 1cef898..e8f1589 100644 --- a/include/aare/File.hpp +++ b/include/aare/File.hpp @@ -5,12 +5,12 @@ namespace aare { /** - * @brief RAII File class for reading, and in the future potentially writing - * image files in various formats. Minimal generic interface. For specail fuctions - * plase use the RawFile or NumpyFile classes directly. - * Wraps FileInterface to abstract the underlying file format - * @note **frame_number** refers the the frame number sent by the detector while **frame_index** - * is the position of the frame in the file + * @brief RAII File class for reading, and in the future potentially writing + * image files in various formats. Minimal generic interface. For specail + * fuctions plase use the RawFile or NumpyFile classes directly. Wraps + * FileInterface to abstract the underlying file format + * @note **frame_number** refers the the frame number sent by the detector while + * **frame_index** is the position of the frame in the file */ class File { std::unique_ptr file_impl; @@ -25,42 +25,46 @@ class File { * @throws std::invalid_argument if the file mode is not supported * */ - File(const std::filesystem::path &fname, const std::string &mode="r", const FileConfig &cfg = {}); - - /**Since the object is responsible for managing the file we disable copy construction */ - File(File const &other) = delete; + File(const std::filesystem::path &fname, const std::string &mode = "r", + const FileConfig &cfg = {}); + + /**Since the object is responsible for managing the file we disable copy + * construction */ + File(File const &other) = delete; /**The same goes for copy assignment */ - File& operator=(File const &other) = delete; + File &operator=(File const &other) = delete; File(File &&other) noexcept; - File& operator=(File &&other) noexcept; + File &operator=(File &&other) noexcept; ~File() = default; // void close(); //!< close the file - - Frame read_frame(); //!< read one frame from the file at the current position - Frame read_frame(size_t frame_index); //!< read one frame at the position given by frame number - std::vector read_n(size_t n_frames); //!< read n_frames from the file at the current position + + Frame + read_frame(); //!< read one frame from the file at the current position + Frame read_frame(size_t frame_index); //!< read one frame at the position + //!< given by frame number + std::vector read_n(size_t n_frames); //!< read n_frames from the file + //!< at the current position void read_into(std::byte *image_buf); void read_into(std::byte *image_buf, size_t n_frames); - - size_t frame_number(); //!< get the frame number at the current position - size_t frame_number(size_t frame_index); //!< get the frame number at the given frame index - size_t bytes_per_frame() const; - size_t pixels_per_frame() const; - size_t bytes_per_pixel() const; + + size_t frame_number(); //!< get the frame number at the current position + size_t frame_number( + size_t frame_index); //!< get the frame number at the given frame index + size_t bytes_per_frame() const; + size_t pixels_per_frame() const; + size_t bytes_per_pixel() const; size_t bitdepth() const; - void seek(size_t frame_index); //!< seek to the given frame index - size_t tell() const; //!< get the frame index of the file pointer + void seek(size_t frame_index); //!< seek to the given frame index + size_t tell() const; //!< get the frame index of the file pointer size_t total_frames() const; size_t rows() const; size_t cols() const; DetectorType detector_type() const; - - }; } // namespace aare \ No newline at end of file diff --git a/include/aare/FileInterface.hpp b/include/aare/FileInterface.hpp index 3736c46..6ca4755 100644 --- a/include/aare/FileInterface.hpp +++ b/include/aare/FileInterface.hpp @@ -20,8 +20,10 @@ struct FileConfig { uint64_t rows{}; uint64_t cols{}; bool operator==(const FileConfig &other) const { - return dtype == other.dtype && rows == other.rows && cols == other.cols && geometry == other.geometry && - detector_type == other.detector_type && max_frames_per_file == other.max_frames_per_file; + return dtype == other.dtype && rows == other.rows && + cols == other.cols && geometry == other.geometry && + detector_type == other.detector_type && + max_frames_per_file == other.max_frames_per_file; } bool operator!=(const FileConfig &other) const { return !(*this == other); } @@ -32,8 +34,11 @@ struct FileConfig { int max_frames_per_file{}; size_t total_frames{}; std::string to_string() const { - return "{ dtype: " + dtype.to_string() + ", rows: " + std::to_string(rows) + ", cols: " + std::to_string(cols) + - ", geometry: " + geometry.to_string() + ", detector_type: " + ToString(detector_type) + + return "{ dtype: " + dtype.to_string() + + ", rows: " + std::to_string(rows) + + ", cols: " + std::to_string(cols) + + ", geometry: " + geometry.to_string() + + ", detector_type: " + ToString(detector_type) + ", max_frames_per_file: " + std::to_string(max_frames_per_file) + ", total_frames: " + std::to_string(total_frames) + " }"; } @@ -42,7 +47,8 @@ struct FileConfig { /** * @brief FileInterface class to define the interface for file operations * @note parent class for NumpyFile and RawFile - * @note all functions are pure virtual and must be implemented by the derived classes + * @note all functions are pure virtual and must be implemented by the derived + * classes */ class FileInterface { public: @@ -64,17 +70,20 @@ class FileInterface { * @param n_frames number of frames to read * @return vector of frames */ - virtual std::vector read_n(size_t n_frames) = 0; // Is this the right interface? + virtual std::vector + read_n(size_t n_frames) = 0; // Is this the right interface? /** - * @brief read one frame from the file at the current position and store it in the provided buffer + * @brief read one frame from the file at the current position and store it + * in the provided buffer * @param image_buf buffer to store the frame * @return void */ virtual void read_into(std::byte *image_buf) = 0; /** - * @brief read n_frames from the file at the current position and store them in the provided buffer + * @brief read n_frames from the file at the current position and store them + * in the provided buffer * @param image_buf buffer to store the frames * @param n_frames number of frames to read * @return void @@ -134,7 +143,6 @@ class FileInterface { */ virtual size_t bitdepth() const = 0; - virtual DetectorType detector_type() const = 0; // function to query the data type of the file diff --git a/include/aare/FilePtr.hpp b/include/aare/FilePtr.hpp index 4ddc76e..2ffc293 100644 --- a/include/aare/FilePtr.hpp +++ b/include/aare/FilePtr.hpp @@ -12,7 +12,7 @@ class FilePtr { public: FilePtr() = default; - FilePtr(const std::filesystem::path& fname, const std::string& mode); + FilePtr(const std::filesystem::path &fname, const std::string &mode); FilePtr(const FilePtr &) = delete; // we don't want a copy FilePtr &operator=(const FilePtr &) = delete; // since we handle a resource FilePtr(FilePtr &&other); diff --git a/include/aare/Fit.hpp b/include/aare/Fit.hpp index eb9ac22..1beec0a 100644 --- a/include/aare/Fit.hpp +++ b/include/aare/Fit.hpp @@ -23,16 +23,19 @@ NDArray scurve2(NDView x, NDView par); } // namespace func - /** * @brief Estimate the initial parameters for a Gaussian fit */ -std::array gaus_init_par(const NDView x, const NDView y); +std::array gaus_init_par(const NDView x, + const NDView y); -std::array pol1_init_par(const NDView x, const NDView y); +std::array pol1_init_par(const NDView x, + const NDView y); -std::array scurve_init_par(const NDView x, const NDView y); -std::array scurve2_init_par(const NDView x, const NDView y); +std::array scurve_init_par(const NDView x, + const NDView y); +std::array scurve2_init_par(const NDView x, + const NDView y); static constexpr int DEFAULT_NUM_THREADS = 4; @@ -43,7 +46,6 @@ static constexpr int DEFAULT_NUM_THREADS = 4; */ NDArray fit_gaus(NDView x, NDView y); - /** * @brief Fit a 1D Gaussian to each pixel. Data layout [row, col, values] * @param x x values @@ -54,9 +56,6 @@ NDArray fit_gaus(NDView x, NDView y); NDArray fit_gaus(NDView x, NDView y, int n_threads = DEFAULT_NUM_THREADS); - - - /** * @brief Fit a 1D Gaussian with error estimates * @param x x values @@ -67,7 +66,7 @@ NDArray fit_gaus(NDView x, NDView y, */ void fit_gaus(NDView x, NDView y, NDView y_err, NDView par_out, NDView par_err_out, - double& chi2); + double &chi2); /** * @brief Fit a 1D Gaussian to each pixel with error estimates. Data layout @@ -80,9 +79,8 @@ void fit_gaus(NDView x, NDView y, NDView y_err, * @param n_threads number of threads to use */ void fit_gaus(NDView x, NDView y, NDView y_err, - NDView par_out, NDView par_err_out, NDView chi2_out, - int n_threads = DEFAULT_NUM_THREADS - ); + NDView par_out, NDView par_err_out, + NDView chi2_out, int n_threads = DEFAULT_NUM_THREADS); NDArray fit_pol1(NDView x, NDView y); @@ -90,26 +88,33 @@ NDArray fit_pol1(NDView x, NDView y, int n_threads = DEFAULT_NUM_THREADS); void fit_pol1(NDView x, NDView y, NDView y_err, - NDView par_out, NDView par_err_out, double& chi2); + NDView par_out, NDView par_err_out, + double &chi2); // TODO! not sure we need to offer the different version in C++ void fit_pol1(NDView x, NDView y, NDView y_err, - NDView par_out, NDView par_err_out,NDView chi2_out, - int n_threads = DEFAULT_NUM_THREADS); + NDView par_out, NDView par_err_out, + NDView chi2_out, int n_threads = DEFAULT_NUM_THREADS); -NDArray fit_scurve(NDView x, NDView y); -NDArray fit_scurve(NDView x, NDView y, int n_threads); -void fit_scurve(NDView x, NDView y, NDView y_err, - NDView par_out, NDView par_err_out, double& chi2); -void fit_scurve(NDView x, NDView y, NDView y_err, - NDView par_out, NDView par_err_out, NDView chi2_out, - int n_threads); +NDArray fit_scurve(NDView x, NDView y); +NDArray fit_scurve(NDView x, NDView y, + int n_threads); +void fit_scurve(NDView x, NDView y, + NDView y_err, NDView par_out, + NDView par_err_out, double &chi2); +void fit_scurve(NDView x, NDView y, + NDView y_err, NDView par_out, + NDView par_err_out, NDView chi2_out, + int n_threads); -NDArray fit_scurve2(NDView x, NDView y); -NDArray fit_scurve2(NDView x, NDView y, int n_threads); -void fit_scurve2(NDView x, NDView y, NDView y_err, - NDView par_out, NDView par_err_out, double& chi2); -void fit_scurve2(NDView x, NDView y, NDView y_err, - NDView par_out, NDView par_err_out, NDView chi2_out, - int n_threads); +NDArray fit_scurve2(NDView x, NDView y); +NDArray fit_scurve2(NDView x, NDView y, + int n_threads); +void fit_scurve2(NDView x, NDView y, + NDView y_err, NDView par_out, + NDView par_err_out, double &chi2); +void fit_scurve2(NDView x, NDView y, + NDView y_err, NDView par_out, + NDView par_err_out, NDView chi2_out, + int n_threads); } // namespace aare \ No newline at end of file diff --git a/include/aare/Frame.hpp b/include/aare/Frame.hpp index 02ea82f..27a2a4a 100644 --- a/include/aare/Frame.hpp +++ b/include/aare/Frame.hpp @@ -19,7 +19,7 @@ class Frame { uint32_t m_cols; Dtype m_dtype; std::byte *m_data; - //TODO! Add frame number? + // TODO! Add frame number? public: /** @@ -39,7 +39,7 @@ class Frame { * @param dtype data type of the pixels */ Frame(const std::byte *bytes, uint32_t rows, uint32_t cols, Dtype dtype); - ~Frame(){ delete[] m_data; }; + ~Frame() { delete[] m_data; }; /** @warning Copy is disabled to ensure performance when passing * frames around. Can discuss enabling it. @@ -52,7 +52,6 @@ class Frame { Frame &operator=(Frame &&other) noexcept; Frame(Frame &&other) noexcept; - Frame clone() const; //<- Explicit copy uint32_t rows() const; @@ -93,7 +92,7 @@ class Frame { if (row >= m_rows || col >= m_cols) { throw std::out_of_range("Invalid row or column index"); } - //TODO! add tests then reimplement using pixel_ptr + // TODO! add tests then reimplement using pixel_ptr T data; std::memcpy(&data, m_data + (row * m_cols + col) * m_dtype.bytes(), m_dtype.bytes()); @@ -102,9 +101,9 @@ class Frame { /** * @brief Return an NDView of the frame. This is the preferred way to access * data in the frame. - * + * * @tparam T type of the pixels - * @return NDView + * @return NDView */ template NDView view() { std::array shape = {static_cast(m_rows), @@ -113,7 +112,7 @@ class Frame { return NDView(data, shape); } - /** + /** * @brief Copy the frame data into a new NDArray. This is a deep copy. */ template NDArray image() { diff --git a/include/aare/Interpolator.hpp b/include/aare/Interpolator.hpp index d2b2322..8e65f38 100644 --- a/include/aare/Interpolator.hpp +++ b/include/aare/Interpolator.hpp @@ -51,7 +51,7 @@ Interpolator::interpolate(const ClusterVector &clusters) { Photon photon; photon.x = cluster.x; photon.y = cluster.y; - photon.energy = eta.sum; + photon.energy = static_cast(eta.sum); // auto ie = nearest_index(m_energy_bins, photon.energy)-1; // auto ix = nearest_index(m_etabinsx, eta.x)-1; @@ -99,7 +99,7 @@ Interpolator::interpolate(const ClusterVector &clusters) { Photon photon; photon.x = cluster.x; photon.y = cluster.y; - photon.energy = eta.sum; + photon.energy = static_cast(eta.sum); // Now do some actual interpolation. // Find which energy bin the cluster is in diff --git a/include/aare/JungfrauDataFile.hpp b/include/aare/JungfrauDataFile.hpp index 9b1bc48..f871b86 100644 --- a/include/aare/JungfrauDataFile.hpp +++ b/include/aare/JungfrauDataFile.hpp @@ -3,104 +3,113 @@ #include #include -#include "aare/FilePtr.hpp" -#include "aare/defs.hpp" -#include "aare/NDArray.hpp" #include "aare/FileInterface.hpp" +#include "aare/FilePtr.hpp" +#include "aare/NDArray.hpp" +#include "aare/defs.hpp" namespace aare { - -struct JungfrauDataHeader{ +struct JungfrauDataHeader { uint64_t framenum; uint64_t bunchid; }; class JungfrauDataFile : public FileInterface { - size_t m_rows{}; //!< number of rows in the image, from find_frame_size(); - size_t m_cols{}; //!< number of columns in the image, from find_frame_size(); + size_t m_rows{}; //!< number of rows in the image, from find_frame_size(); + size_t + m_cols{}; //!< number of columns in the image, from find_frame_size(); size_t m_bytes_per_frame{}; //!< number of bytes per frame excluding header - size_t m_total_frames{}; //!< total number of frames in the series of files - size_t m_offset{}; //!< file index of the first file, allow starting at non zero file - size_t m_current_file_index{}; //!< The index of the open file - size_t m_current_frame_index{}; //!< The index of the current frame (with reference to all files) + size_t m_total_frames{}; //!< total number of frames in the series of files + size_t m_offset{}; //!< file index of the first file, allow starting at non + //!< zero file + size_t m_current_file_index{}; //!< The index of the open file + size_t m_current_frame_index{}; //!< The index of the current frame (with + //!< reference to all files) - std::vector m_last_frame_in_file{}; //!< Used for seeking to the correct file + std::vector + m_last_frame_in_file{}; //!< Used for seeking to the correct file std::filesystem::path m_path; //!< path to the files std::string m_base_name; //!< base name used for formatting file names - FilePtr m_fp; //!< RAII wrapper for a FILE* + FilePtr m_fp; //!< RAII wrapper for a FILE* - using pixel_type = uint16_t; - static constexpr size_t header_size = sizeof(JungfrauDataHeader); - static constexpr size_t n_digits_in_file_index = 6; //!< to format file names + static constexpr size_t header_size = sizeof(JungfrauDataHeader); + static constexpr size_t n_digits_in_file_index = + 6; //!< to format file names public: JungfrauDataFile(const std::filesystem::path &fname); - std::string base_name() const; //!< get the base name of the file (without path and extension) - size_t bytes_per_frame() override; - size_t pixels_per_frame() override; - size_t bytes_per_pixel() const; + std::string base_name() + const; //!< get the base name of the file (without path and extension) + size_t bytes_per_frame() override; + size_t pixels_per_frame() override; + size_t bytes_per_pixel() const; size_t bitdepth() const override; - 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 + 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 size_t total_frames() const override; size_t rows() const override; size_t cols() const override; - std::array shape() const; - size_t n_files() const; //!< get the number of files in the series. + std::array shape() const; + size_t n_files() const; //!< get the number of files in the series. // Extra functions needed for FileInterface Frame read_frame() override; Frame read_frame(size_t frame_number) override; - std::vector read_n(size_t n_frames=0) override; + std::vector read_n(size_t n_frames = 0) override; void read_into(std::byte *image_buf) override; void read_into(std::byte *image_buf, size_t n_frames) override; size_t frame_number(size_t frame_index) override; DetectorType detector_type() const override; /** - * @brief Read a single frame from the file into the given buffer. - * @param image_buf buffer to read the frame into. (Note the caller is responsible for allocating the buffer) + * @brief Read a single frame from the file into the given buffer. + * @param image_buf buffer to read the frame into. (Note the caller is + * responsible for allocating the buffer) * @param header pointer to a JungfrauDataHeader or nullptr to skip header) */ void read_into(std::byte *image_buf, JungfrauDataHeader *header = nullptr); /** - * @brief Read a multiple frames from the file into the given buffer. - * @param image_buf buffer to read the frame into. (Note the caller is responsible for allocating the buffer) + * @brief Read a multiple frames from the file into the given buffer. + * @param image_buf buffer to read the frame into. (Note the caller is + * responsible for allocating the buffer) * @param n_frames number of frames to read * @param header pointer to a JungfrauDataHeader or nullptr to skip header) */ - void read_into(std::byte *image_buf, size_t n_frames, JungfrauDataHeader *header = nullptr); - - /** + void read_into(std::byte *image_buf, size_t n_frames, + JungfrauDataHeader *header = nullptr); + + /** * @brief Read a single frame from the file into the given NDArray * @param image NDArray to read the frame into. */ - void read_into(NDArray* image, JungfrauDataHeader* header = nullptr); + void read_into(NDArray *image, + JungfrauDataHeader *header = nullptr); JungfrauDataHeader read_header(); - std::filesystem::path current_file() const { return fpath(m_current_file_index+m_offset); } + std::filesystem::path current_file() const { + return fpath(m_current_file_index + m_offset); + } - - private: + private: /** - * @brief Find the size of the frame in the file. (256x256, 256x1024, 512x1024) + * @brief Find the size of the frame in the file. (256x256, 256x1024, + * 512x1024) * @param fname path to the file - * @throws std::runtime_error if the file is empty or the size cannot be determined + * @throws std::runtime_error if the file is empty or the size cannot be + * determined */ - void find_frame_size(const std::filesystem::path &fname); + void find_frame_size(const std::filesystem::path &fname); - - void parse_fname(const std::filesystem::path &fname); - void scan_files(); - void open_file(size_t file_index); - std::filesystem::path fpath(size_t frame_index) const; - - - }; + void parse_fname(const std::filesystem::path &fname); + void scan_files(); + void open_file(size_t file_index); + std::filesystem::path fpath(size_t frame_index) const; +}; } // namespace aare \ No newline at end of file diff --git a/include/aare/NumpyFile.hpp b/include/aare/NumpyFile.hpp index 7381a76..481a1a0 100644 --- a/include/aare/NumpyFile.hpp +++ b/include/aare/NumpyFile.hpp @@ -1,9 +1,8 @@ #pragma once #include "aare/Dtype.hpp" -#include "aare/defs.hpp" #include "aare/FileInterface.hpp" #include "aare/NumpyHelpers.hpp" - +#include "aare/defs.hpp" #include #include @@ -11,13 +10,12 @@ namespace aare { - - /** * @brief NumpyFile class to read and write numpy files * @note derived from FileInterface * @note implements all the pure virtual functions from FileInterface - * @note documentation for the functions can also be found in the FileInterface class + * @note documentation for the functions can also be found in the FileInterface + * class */ class NumpyFile : public FileInterface { @@ -28,26 +26,35 @@ class NumpyFile : public FileInterface { * @param mode file mode (r, w) * @param cfg file configuration */ - explicit NumpyFile(const std::filesystem::path &fname, const std::string &mode = "r", FileConfig cfg = {}); + explicit NumpyFile(const std::filesystem::path &fname, + const std::string &mode = "r", FileConfig cfg = {}); void write(Frame &frame); Frame read_frame() override { return get_frame(this->current_frame++); } - Frame read_frame(size_t frame_number) override { return get_frame(frame_number); } + Frame read_frame(size_t frame_number) override { + return get_frame(frame_number); + } std::vector read_n(size_t n_frames) override; - void read_into(std::byte *image_buf) override { return get_frame_into(this->current_frame++, image_buf); } + void read_into(std::byte *image_buf) override { + return get_frame_into(this->current_frame++, image_buf); + } void read_into(std::byte *image_buf, size_t n_frames) override; size_t frame_number(size_t frame_index) override { return frame_index; }; size_t bytes_per_frame() override; size_t pixels_per_frame() override; - void seek(size_t frame_number) override { this->current_frame = frame_number; } + void seek(size_t frame_number) override { + this->current_frame = frame_number; + } size_t tell() override { return this->current_frame; } size_t total_frames() const override { return m_header.shape[0]; } size_t rows() const override { return m_header.shape[1]; } size_t cols() const override { return m_header.shape[2]; } size_t bitdepth() const override { return m_header.dtype.bitdepth(); } - DetectorType detector_type() const override { return DetectorType::Unknown; } + DetectorType detector_type() const override { + return DetectorType::Unknown; + } /** * @brief get the data type of the numpy file @@ -70,7 +77,8 @@ class NumpyFile : public FileInterface { template NDArray load() { NDArray arr(make_shape(m_header.shape)); if (fseek(fp, static_cast(header_size), SEEK_SET)) { - throw std::runtime_error(LOCATION + "Error seeking to the start of the data"); + throw std::runtime_error(LOCATION + + "Error seeking to the start of the data"); } size_t rc = fread(arr.data(), sizeof(T), arr.size(), fp); if (rc != static_cast(arr.size())) { @@ -78,16 +86,20 @@ class NumpyFile : public FileInterface { } return arr; } - template void write(NDView &frame) { + template + void write(NDView &frame) { write_impl(frame.data(), frame.total_bytes()); } - template void write(NDArray &frame) { + template + void write(NDArray &frame) { write_impl(frame.data(), frame.total_bytes()); } - template void write(NDView &&frame) { + template + void write(NDView &&frame) { write_impl(frame.data(), frame.total_bytes()); } - template void write(NDArray &&frame) { + template + void write(NDArray &&frame) { write_impl(frame.data(), frame.total_bytes()); } diff --git a/include/aare/NumpyHelpers.hpp b/include/aare/NumpyHelpers.hpp index 8ed0ec7..2facc4c 100644 --- a/include/aare/NumpyHelpers.hpp +++ b/include/aare/NumpyHelpers.hpp @@ -40,15 +40,18 @@ bool parse_bool(const std::string &in); std::string get_value_from_map(const std::string &mapstr); -std::unordered_map parse_dict(std::string in, const std::vector &keys); +std::unordered_map +parse_dict(std::string in, const std::vector &keys); -template bool in_array(T val, const std::array &arr) { +template +bool in_array(T val, const std::array &arr) { return std::find(std::begin(arr), std::end(arr), val) != std::end(arr); } bool is_digits(const std::string &str); aare::Dtype parse_descr(std::string typestring); -size_t write_header(const std::filesystem::path &fname, const NumpyHeader &header); +size_t write_header(const std::filesystem::path &fname, + const NumpyHeader &header); size_t write_header(std::ostream &out, const NumpyHeader &header); } // namespace NumpyHelpers diff --git a/include/aare/Pedestal.hpp b/include/aare/Pedestal.hpp index d6223c1..0efc4b7 100644 --- a/include/aare/Pedestal.hpp +++ b/include/aare/Pedestal.hpp @@ -18,15 +18,15 @@ template class Pedestal { uint32_t m_samples; NDArray m_cur_samples; - - //TODO! in case of int needs to be changed to uint64_t + + // TODO! in case of int needs to be changed to uint64_t NDArray m_sum; NDArray m_sum2; - //Cache mean since it is used over and over in the ClusterFinder - //This optimization is related to the access pattern of the ClusterFinder - //Relies on having more reads than pushes to the pedestal - NDArray m_mean; + // Cache mean since it is used over and over in the ClusterFinder + // This optimization is related to the access pattern of the ClusterFinder + // Relies on having more reads than pushes to the pedestal + NDArray m_mean; public: Pedestal(uint32_t rows, uint32_t cols, uint32_t n_samples = 1000) @@ -42,9 +42,7 @@ template class Pedestal { } ~Pedestal() = default; - NDArray mean() { - return m_mean; - } + NDArray mean() { return m_mean; } SUM_TYPE mean(const uint32_t row, const uint32_t col) const { return m_mean(row, col); @@ -71,8 +69,6 @@ template class Pedestal { return variance_array; } - - NDArray std() { NDArray standard_deviation_array({m_rows, m_cols}); for (uint32_t i = 0; i < m_rows * m_cols; i++) { @@ -83,8 +79,6 @@ template class Pedestal { return standard_deviation_array; } - - void clear() { m_sum = 0; m_sum2 = 0; @@ -92,16 +86,12 @@ template class Pedestal { m_mean = 0; } - - void clear(const uint32_t row, const uint32_t col) { m_sum(row, col) = 0; m_sum2(row, col) = 0; m_cur_samples(row, col) = 0; m_mean(row, col) = 0; } - - template void push(NDView frame) { assert(frame.size() == m_rows * m_cols); @@ -122,7 +112,7 @@ template class Pedestal { /** * Push but don't update the cached mean. Speeds up the process * when initializing the pedestal. - * + * */ template void push_no_update(NDView frame) { assert(frame.size() == m_rows * m_cols); @@ -140,9 +130,6 @@ template class Pedestal { } } - - - template void push(Frame &frame) { assert(frame.rows() == static_cast(m_rows) && frame.cols() == static_cast(m_cols)); @@ -170,7 +157,8 @@ template class Pedestal { m_sum(row, col) += val - m_sum(row, col) / m_samples; m_sum2(row, col) += val * val - m_sum2(row, col) / m_samples; } - //Since we just did a push we know that m_cur_samples(row, col) is at least 1 + // Since we just did a push we know that m_cur_samples(row, col) is at + // least 1 m_mean(row, col) = m_sum(row, col) / m_cur_samples(row, col); } @@ -183,7 +171,8 @@ template class Pedestal { m_cur_samples(row, col)++; } else { m_sum(row, col) += val - m_sum(row, col) / m_cur_samples(row, col); - m_sum2(row, col) += val * val - m_sum2(row, col) / m_cur_samples(row, col); + m_sum2(row, col) += + val * val - m_sum2(row, col) / m_cur_samples(row, col); } } @@ -191,19 +180,16 @@ template class Pedestal { * @brief Update the mean of the pedestal. This is used after having done * push_no_update. It is not necessary to call this function after push. */ - void update_mean(){ - m_mean = m_sum / m_cur_samples; - } + void update_mean() { m_mean = m_sum / m_cur_samples; } - template - void push_fast(const uint32_t row, const uint32_t col, const T val_){ - //Assume we reached the steady state where all pixels have - //m_samples samples + template + void push_fast(const uint32_t row, const uint32_t col, const T val_) { + // Assume we reached the steady state where all pixels have + // m_samples samples SUM_TYPE val = static_cast(val_); m_sum(row, col) += val - m_sum(row, col) / m_samples; m_sum2(row, col) += val * val - m_sum2(row, col) / m_samples; m_mean(row, col) = m_sum(row, col) / m_samples; } - }; } // namespace aare \ No newline at end of file diff --git a/include/aare/PixelMap.hpp b/include/aare/PixelMap.hpp index 1b7a890..9c30680 100644 --- a/include/aare/PixelMap.hpp +++ b/include/aare/PixelMap.hpp @@ -1,7 +1,7 @@ #pragma once -#include "aare/defs.hpp" #include "aare/NDArray.hpp" +#include "aare/defs.hpp" namespace aare { @@ -10,11 +10,11 @@ NDArray GenerateMoench05PixelMap(); NDArray GenerateMoench05PixelMap1g(); NDArray GenerateMoench05PixelMapOld(); -//Matterhorn02 -NDArrayGenerateMH02SingleCounterPixelMap(); +// Matterhorn02 +NDArray GenerateMH02SingleCounterPixelMap(); NDArray GenerateMH02FourCounterPixelMap(); -//Eiger -NDArrayGenerateEigerFlipRowsPixelMap(); +// Eiger +NDArray GenerateEigerFlipRowsPixelMap(); } // namespace aare \ No newline at end of file diff --git a/include/aare/ProducerConsumerQueue.hpp b/include/aare/ProducerConsumerQueue.hpp index 426b9e2..f189cec 100644 --- a/include/aare/ProducerConsumerQueue.hpp +++ b/include/aare/ProducerConsumerQueue.hpp @@ -18,9 +18,9 @@ // @author Jordan DeLong (delong.j@fb.com) // Changes made by PSD Detector Group: -// Copied: Line 34 constexpr std::size_t hardware_destructive_interference_size = 128; from folly/lang/Align.h -// Changed extension to .hpp -// Changed namespace to aare +// Copied: Line 34 constexpr std::size_t hardware_destructive_interference_size +// = 128; from folly/lang/Align.h Changed extension to .hpp Changed namespace to +// aare #pragma once @@ -45,15 +45,14 @@ template struct ProducerConsumerQueue { ProducerConsumerQueue(const ProducerConsumerQueue &) = delete; ProducerConsumerQueue &operator=(const ProducerConsumerQueue &) = delete; - - ProducerConsumerQueue(ProducerConsumerQueue &&other){ + ProducerConsumerQueue(ProducerConsumerQueue &&other) { size_ = other.size_; records_ = other.records_; other.records_ = nullptr; readIndex_ = other.readIndex_.load(std::memory_order_acquire); writeIndex_ = other.writeIndex_.load(std::memory_order_acquire); } - ProducerConsumerQueue &operator=(ProducerConsumerQueue &&other){ + ProducerConsumerQueue &operator=(ProducerConsumerQueue &&other) { size_ = other.size_; records_ = other.records_; other.records_ = nullptr; @@ -61,16 +60,17 @@ template struct ProducerConsumerQueue { writeIndex_ = other.writeIndex_.load(std::memory_order_acquire); return *this; } - - - ProducerConsumerQueue():ProducerConsumerQueue(2){}; + + ProducerConsumerQueue() : ProducerConsumerQueue(2){}; // size must be >= 2. // // Also, note that the number of usable slots in the queue at any // given time is actually (size-1), so if you start with an empty queue, // isFull() will return true after size-1 insertions. explicit ProducerConsumerQueue(uint32_t size) - : size_(size), records_(static_cast(std::malloc(sizeof(T) * size))), readIndex_(0), writeIndex_(0) { + : size_(size), + records_(static_cast(std::malloc(sizeof(T) * size))), + readIndex_(0), writeIndex_(0) { assert(size >= 2); if (!records_) { throw std::bad_alloc(); @@ -154,7 +154,8 @@ template struct ProducerConsumerQueue { } bool isEmpty() const { - return readIndex_.load(std::memory_order_acquire) == writeIndex_.load(std::memory_order_acquire); + return readIndex_.load(std::memory_order_acquire) == + writeIndex_.load(std::memory_order_acquire); } bool isFull() const { @@ -175,7 +176,8 @@ template struct ProducerConsumerQueue { // be removing items concurrently). // * It is undefined to call this from any other thread. size_t sizeGuess() const { - int ret = writeIndex_.load(std::memory_order_acquire) - readIndex_.load(std::memory_order_acquire); + int ret = writeIndex_.load(std::memory_order_acquire) - + readIndex_.load(std::memory_order_acquire); if (ret < 0) { ret += size_; } @@ -192,7 +194,7 @@ template struct ProducerConsumerQueue { // const uint32_t size_; uint32_t size_; // T *const records_; - T* records_; + T *records_; alignas(hardware_destructive_interference_size) AtomicIndex readIndex_; alignas(hardware_destructive_interference_size) AtomicIndex writeIndex_; diff --git a/include/aare/RawFile.hpp b/include/aare/RawFile.hpp index f744ac2..9ffdb7c 100644 --- a/include/aare/RawFile.hpp +++ b/include/aare/RawFile.hpp @@ -1,11 +1,10 @@ #pragma once #include "aare/FileInterface.hpp" -#include "aare/RawMasterFile.hpp" #include "aare/Frame.hpp" #include "aare/NDArray.hpp" //for pixel map +#include "aare/RawMasterFile.hpp" #include "aare/RawSubFile.hpp" - #include namespace aare { @@ -30,22 +29,11 @@ struct ModuleConfig { * Consider using that unless you need raw file specific functionality. */ class RawFile : public FileInterface { - size_t n_subfiles{}; //f0,f1...fn - size_t n_subfile_parts{}; // d0,d1...dn - //TODO! move to vector of SubFile instead of pointers - std::vector> subfiles; //subfiles[f0,f1...fn][d0,d1...dn] - // std::vector positions; - + std::vector> m_subfiles; ModuleConfig cfg{0, 0}; - RawMasterFile m_master; - size_t m_current_frame{}; - - // std::vector m_module_pixel_0; - // size_t m_rows{}; - // size_t m_cols{}; - + size_t m_current_subfile{}; DetectorGeometry m_geometry; public: @@ -56,7 +44,7 @@ class RawFile : public FileInterface { */ RawFile(const std::filesystem::path &fname, const std::string &mode = "r"); - virtual ~RawFile() override; + virtual ~RawFile() override = default; Frame read_frame() override; Frame read_frame(size_t frame_number) override; @@ -64,10 +52,10 @@ class RawFile : public FileInterface { void read_into(std::byte *image_buf) override; void read_into(std::byte *image_buf, size_t n_frames) override; - //TODO! do we need to adapt the API? + // TODO! do we need to adapt the API? void read_into(std::byte *image_buf, DetectorHeader *header); - void read_into(std::byte *image_buf, size_t n_frames, DetectorHeader *header); - + void read_into(std::byte *image_buf, size_t n_frames, + DetectorHeader *header); size_t frame_number(size_t frame_index) override; size_t bytes_per_frame() override; @@ -80,24 +68,21 @@ class RawFile : public FileInterface { size_t cols() const override; size_t bitdepth() const override; xy geometry(); - size_t n_mod() const; - + size_t n_modules() const; + RawMasterFile master() const; - - DetectorType detector_type() const override; private: - /** * @brief read the frame at the given frame index into the image buffer * @param frame_number frame number to read * @param image_buf buffer to store the frame */ - void get_frame_into(size_t frame_index, std::byte *frame_buffer, DetectorHeader *header = nullptr); - + void get_frame_into(size_t frame_index, std::byte *frame_buffer, + DetectorHeader *header = nullptr); /** * @brief get the frame at the given frame index @@ -106,8 +91,6 @@ class RawFile : public FileInterface { */ Frame get_frame(size_t frame_index); - - /** * @brief read the header of the file * @param fname path to the data subfile @@ -115,12 +98,8 @@ class RawFile : public FileInterface { */ static DetectorHeader read_header(const std::filesystem::path &fname); - // void update_geometry_with_roi(); - int find_number_of_subfiles(); - void open_subfiles(); void find_geometry(); }; - } // namespace aare \ No newline at end of file diff --git a/include/aare/RawMasterFile.hpp b/include/aare/RawMasterFile.hpp index beaeb29..2c64a90 100644 --- a/include/aare/RawMasterFile.hpp +++ b/include/aare/RawMasterFile.hpp @@ -45,7 +45,7 @@ class ScanParameters { int m_start = 0; int m_stop = 0; int m_step = 0; - //TODO! add settleTime, requires string to time conversion + // TODO! add settleTime, requires string to time conversion public: ScanParameters(const std::string &par); @@ -61,7 +61,6 @@ class ScanParameters { void increment_stop(); }; - /** * @brief Class for parsing a master file either in our .json format or the old * .raw format @@ -101,7 +100,6 @@ class RawMasterFile { std::optional m_roi; - public: RawMasterFile(const std::filesystem::path &fpath); @@ -121,6 +119,7 @@ class RawMasterFile { size_t total_frames_expected() const; xy geometry() const; + size_t n_modules() const; std::optional analog_samples() const; std::optional digital_samples() const; @@ -128,10 +127,8 @@ class RawMasterFile { std::optional number_of_rows() const; std::optional quad() const; - std::optional roi() const; - ScanParameters scan_parameters() const; private: diff --git a/include/aare/RawSubFile.hpp b/include/aare/RawSubFile.hpp index 350a475..c38f540 100644 --- a/include/aare/RawSubFile.hpp +++ b/include/aare/RawSubFile.hpp @@ -10,23 +10,34 @@ namespace aare { /** - * @brief Class to read a singe subfile written in .raw format. Used from RawFile to read - * the entire detector. Can be used directly to read part of the image. + * @brief Class to read a singe subfile written in .raw format. Used from + * RawFile to read the entire detector. Can be used directly to read part of the + * image. */ class RawSubFile { protected: std::ifstream m_file; DetectorType m_detector_type; size_t m_bitdepth; - std::filesystem::path m_fname; + std::filesystem::path m_path; //!< path to the subfile + std::string m_base_name; //!< base name used for formatting file names + size_t m_offset{}; //!< file index of the first file, allow starting at non + //!< zero file + size_t m_total_frames{}; //!< total number of frames in the series of files size_t m_rows{}; size_t m_cols{}; size_t m_bytes_per_frame{}; - size_t m_num_frames{}; + + int m_module_index{}; + size_t m_current_file_index{}; //!< The index of the open file + size_t m_current_frame_index{}; //!< The index of the current frame (with + //!< reference to all files) + std::vector + m_last_frame_in_file{}; //!< Used for seeking to the correct file + uint32_t m_pos_row{}; uint32_t m_pos_col{}; - - + std::optional> m_pixel_map; public: @@ -40,12 +51,14 @@ class RawSubFile { * @throws std::invalid_argument if the detector,type pair is not supported */ RawSubFile(const std::filesystem::path &fname, DetectorType detector, - size_t rows, size_t cols, size_t bitdepth, uint32_t pos_row = 0, uint32_t pos_col = 0); + size_t rows, size_t cols, size_t bitdepth, uint32_t pos_row = 0, + uint32_t pos_col = 0); ~RawSubFile() = default; /** * @brief Seek to the given frame number - * @note Puts the file pointer at the start of the header, not the start of the data + * @note Puts the file pointer at the start of the header, not the start of + * the data * @param frame_index frame position in file to seek to * @throws std::runtime_error if the frame number is out of range */ @@ -53,26 +66,30 @@ class RawSubFile { size_t tell(); void read_into(std::byte *image_buf, DetectorHeader *header = nullptr); - void read_into(std::byte *image_buf, size_t n_frames, DetectorHeader *header= nullptr); + void read_into(std::byte *image_buf, size_t n_frames, + DetectorHeader *header = nullptr); void get_part(std::byte *buffer, size_t frame_index); - + void read_header(DetectorHeader *header); - + size_t rows() const; size_t cols() const; - + size_t frame_number(size_t frame_index); size_t bytes_per_frame() const { return m_bytes_per_frame; } size_t pixels_per_frame() const { return m_rows * m_cols; } size_t bytes_per_pixel() const { return m_bitdepth / bits_per_byte; } - size_t frames_in_file() const { return m_num_frames; } + size_t frames_in_file() const { return m_total_frames; } -private: - template - void read_with_map(std::byte *image_buf); + private: + template void read_with_map(std::byte *image_buf); + void parse_fname(const std::filesystem::path &fname); + void scan_files(); + void open_file(size_t file_index); + std::filesystem::path fpath(size_t file_index) const; }; } // namespace aare \ No newline at end of file diff --git a/include/aare/VarClusterFinder.hpp b/include/aare/VarClusterFinder.hpp index 596bf06..3679d39 100644 --- a/include/aare/VarClusterFinder.hpp +++ b/include/aare/VarClusterFinder.hpp @@ -38,11 +38,13 @@ template class VarClusterFinder { bool use_noise_map = false; int peripheralThresholdFactor_ = 5; int current_label; - const std::array di{{0, -1, -1, -1}}; // row ### 8-neighbour by scaning from left to right - const std::array dj{{-1, -1, 0, 1}}; // col ### 8-neighbour by scaning from top to bottom + const std::array di{ + {0, -1, -1, -1}}; // row ### 8-neighbour by scaning from left to right + const std::array dj{ + {-1, -1, 0, 1}}; // col ### 8-neighbour by scaning from top to bottom const std::array di_{{0, 0, -1, 1, -1, 1, -1, 1}}; // row const std::array dj_{{-1, 1, 0, 0, 1, -1, -1, 1}}; // col - std::map child; // heirachy: key: child; val: parent + std::map child; // heirachy: key: child; val: parent std::unordered_map h_size; std::vector hits; // std::vector> row @@ -50,7 +52,8 @@ template class VarClusterFinder { public: VarClusterFinder(Shape<2> shape, T threshold) - : shape_(shape), labeled_(shape, 0), peripheral_labeled_(shape, 0), binary_(shape), threshold_(threshold) { + : shape_(shape), labeled_(shape, 0), peripheral_labeled_(shape, 0), + binary_(shape), threshold_(threshold) { hits.reserve(2000); } @@ -60,7 +63,9 @@ template class VarClusterFinder { noiseMap = noise_map; use_noise_map = true; } - void set_peripheralThresholdFactor(int factor) { peripheralThresholdFactor_ = factor; } + void set_peripheralThresholdFactor(int factor) { + peripheralThresholdFactor_ = factor; + } void find_clusters(NDView img); void find_clusters_X(NDView img); void rec_FillHit(int clusterIndex, int i, int j); @@ -144,7 +149,8 @@ template int VarClusterFinder::check_neighbours(int i, int j) { } } -template void VarClusterFinder::find_clusters(NDView img) { +template +void VarClusterFinder::find_clusters(NDView img) { original_ = img; labeled_ = 0; peripheral_labeled_ = 0; @@ -156,7 +162,8 @@ template void VarClusterFinder::find_clusters(NDView img) store_clusters(); } -template void VarClusterFinder::find_clusters_X(NDView img) { +template +void VarClusterFinder::find_clusters_X(NDView img) { original_ = img; int clusterIndex = 0; for (int i = 0; i < shape_[0]; ++i) { @@ -175,7 +182,8 @@ template void VarClusterFinder::find_clusters_X(NDView img h_size.clear(); } -template void VarClusterFinder::rec_FillHit(int clusterIndex, int i, int j) { +template +void VarClusterFinder::rec_FillHit(int clusterIndex, int i, int j) { // printf("original_(%d, %d)=%f\n", i, j, original_(i,j)); // printf("h_size[%d].size=%d\n", clusterIndex, h_size[clusterIndex].size); if (h_size[clusterIndex].size < MAX_CLUSTER_SIZE) { @@ -203,11 +211,15 @@ template void VarClusterFinder::rec_FillHit(int clusterIndex, in } else { // if (h_size[clusterIndex].size < MAX_CLUSTER_SIZE){ // h_size[clusterIndex].size += 1; - // h_size[clusterIndex].rows[h_size[clusterIndex].size] = row; - // h_size[clusterIndex].cols[h_size[clusterIndex].size] = col; - // h_size[clusterIndex].enes[h_size[clusterIndex].size] = original_(row, col); + // h_size[clusterIndex].rows[h_size[clusterIndex].size] = + // row; h_size[clusterIndex].cols[h_size[clusterIndex].size] + // = col; + // h_size[clusterIndex].enes[h_size[clusterIndex].size] = + // original_(row, col); // }// ? weather to include peripheral pixels - original_(row, col) = 0; // remove peripheral pixels, to avoid potential influence for pedestal updating + original_(row, col) = + 0; // remove peripheral pixels, to avoid potential influence + // for pedestal updating } } } @@ -275,8 +287,8 @@ template void VarClusterFinder::store_clusters() { for (int i = 0; i < shape_[0]; ++i) { for (int j = 0; j < shape_[1]; ++j) { if (labeled_(i, j) != 0 || false - // (i-1 >= 0 and labeled_(i-1, j) != 0) or // another circle of peripheral pixels - // (j-1 >= 0 and labeled_(i, j-1) != 0) or + // (i-1 >= 0 and labeled_(i-1, j) != 0) or // another circle of + // peripheral pixels (j-1 >= 0 and labeled_(i, j-1) != 0) or // (i+1 < shape_[0] and labeled_(i+1, j) != 0) or // (j+1 < shape_[1] and labeled_(i, j+1) != 0) ) { diff --git a/include/aare/algorithm.hpp b/include/aare/algorithm.hpp index fc7d51f..c590e91 100644 --- a/include/aare/algorithm.hpp +++ b/include/aare/algorithm.hpp @@ -1,9 +1,9 @@ #pragma once +#include #include #include #include -#include namespace aare { /** @@ -15,26 +15,24 @@ namespace aare { * @param last iterator to the last element * @param val value to compare * @return index of the last element that is smaller than val - * + * */ template -size_t last_smaller(const T* first, const T* last, T val) { - for (auto iter = first+1; iter != last; ++iter) { +size_t last_smaller(const T *first, const T *last, T val) { + for (auto iter = first + 1; iter != last; ++iter) { if (*iter >= val) { - return std::distance(first, iter-1); + return std::distance(first, iter - 1); } } - return std::distance(first, last-1); + return std::distance(first, last - 1); } -template -size_t last_smaller(const NDArray& arr, T val) { +template size_t last_smaller(const NDArray &arr, T val) { return last_smaller(arr.begin(), arr.end(), val); } -template -size_t last_smaller(const std::vector& vec, T val) { - return last_smaller(vec.data(), vec.data()+vec.size(), val); +template size_t last_smaller(const std::vector &vec, T val) { + return last_smaller(vec.data(), vec.data() + vec.size(), val); } /** @@ -48,64 +46,67 @@ size_t last_smaller(const std::vector& vec, T val) { * @return index of the first element that is larger than val */ template -size_t first_larger(const T* first, const T* last, T val) { +size_t first_larger(const T *first, const T *last, T val) { for (auto iter = first; iter != last; ++iter) { if (*iter > val) { return std::distance(first, iter); } } - return std::distance(first, last-1); + return std::distance(first, last - 1); } -template -size_t first_larger(const NDArray& arr, T val) { +template size_t first_larger(const NDArray &arr, T val) { return first_larger(arr.begin(), arr.end(), val); } -template -size_t first_larger(const std::vector& vec, T val) { - return first_larger(vec.data(), vec.data()+vec.size(), val); +template size_t first_larger(const std::vector &vec, T val) { + return first_larger(vec.data(), vec.data() + vec.size(), val); } /** * @brief Index of the nearest element to val. - * Requires a sorted array. If there is no difference it takes the first element. + * Requires a sorted array. If there is no difference it takes the first + * element. * @param first iterator to the first element * @param last iterator to the last element * @param val value to compare * @return index of the nearest element */ template -size_t nearest_index(const T* first, const T* last, T val) { - auto iter = std::min_element(first, last, - [val](T a, T b) { +size_t nearest_index(const T *first, const T *last, T val) { + auto iter = std::min_element(first, last, [val](T a, T b) { return std::abs(a - val) < std::abs(b - val); }); return std::distance(first, iter); } -template -size_t nearest_index(const NDArray& arr, T val) { +template size_t nearest_index(const NDArray &arr, T val) { return nearest_index(arr.begin(), arr.end(), val); } -template -size_t nearest_index(const std::vector& vec, T val) { - return nearest_index(vec.data(), vec.data()+vec.size(), val); +template size_t nearest_index(const std::vector &vec, T val) { + return nearest_index(vec.data(), vec.data() + vec.size(), val); } template -size_t nearest_index(const std::array& arr, T val) { - return nearest_index(arr.data(), arr.data()+arr.size(), val); +size_t nearest_index(const std::array &arr, T val) { + return nearest_index(arr.data(), arr.data() + arr.size(), val); } -template -std::vector cumsum(const std::vector& vec) { +template std::vector cumsum(const std::vector &vec) { std::vector result(vec.size()); std::partial_sum(vec.begin(), vec.end(), result.begin()); return result; } - +template bool all_equal(const Container &c) { + if (!c.empty() && + std::all_of(begin(c), end(c), + [c](const typename Container::value_type &element) { + return element == c.front(); + })) + return true; + return false; +} } // namespace aare \ No newline at end of file diff --git a/include/aare/decode.hpp b/include/aare/decode.hpp index e784c4a..ec24447 100644 --- a/include/aare/decode.hpp +++ b/include/aare/decode.hpp @@ -1,26 +1,27 @@ #pragma once +#include #include #include -#include namespace aare { - uint16_t adc_sar_05_decode64to16(uint64_t input); uint16_t adc_sar_04_decode64to16(uint64_t input); -void adc_sar_05_decode64to16(NDView input, NDView output); -void adc_sar_04_decode64to16(NDView input, NDView output); - +void adc_sar_05_decode64to16(NDView input, + NDView output); +void adc_sar_04_decode64to16(NDView input, + NDView output); /** - * @brief Apply custom weights to a 16-bit input value. Will sum up weights[i]**i - * for each bit i that is set in the input value. + * @brief Apply custom weights to a 16-bit input value. Will sum up + * weights[i]**i for each bit i that is set in the input value. * @throws std::out_of_range if weights.size() < 16 * @param input 16-bit input value * @param weights vector of weights, size must be less than or equal to 16 */ double apply_custom_weights(uint16_t input, const NDView weights); -void apply_custom_weights(NDView input, NDView output, const NDView weights); +void apply_custom_weights(NDView input, NDView output, + const NDView weights); } // namespace aare diff --git a/include/aare/defs.hpp b/include/aare/defs.hpp index ea02059..6f28965 100644 --- a/include/aare/defs.hpp +++ b/include/aare/defs.hpp @@ -200,6 +200,8 @@ struct DetectorGeometry { int module_gap_row{}; int module_gap_col{}; std::vector module_pixel_0; + + auto size() const { return module_pixel_0.size(); } }; struct ROI { diff --git a/include/aare/geo_helpers.hpp b/include/aare/geo_helpers.hpp index d0d5d1a..c6454a5 100644 --- a/include/aare/geo_helpers.hpp +++ b/include/aare/geo_helpers.hpp @@ -1,16 +1,15 @@ #pragma once -#include "aare/defs.hpp" #include "aare/RawMasterFile.hpp" //ROI refactor away -namespace aare{ +#include "aare/defs.hpp" +namespace aare { /** * @brief Update the detector geometry given a region of interest - * - * @param geo - * @param roi - * @return DetectorGeometry + * + * @param geo + * @param roi + * @return DetectorGeometry */ DetectorGeometry update_geometry_with_roi(DetectorGeometry geo, ROI roi); - } // namespace aare \ No newline at end of file diff --git a/include/aare/logger.hpp b/include/aare/logger.hpp new file mode 100644 index 0000000..0bedd7a --- /dev/null +++ b/include/aare/logger.hpp @@ -0,0 +1,141 @@ +#pragma once +/*Utility to log to console*/ + +#include +#include +#include + +namespace aare { + +#define RED "\x1b[31m" +#define GREEN "\x1b[32m" +#define YELLOW "\x1b[33m" +#define BLUE "\x1b[34m" +#define MAGENTA "\x1b[35m" +#define CYAN "\x1b[36m" +#define GRAY "\x1b[37m" +#define DARKGRAY "\x1b[30m" + +#define BG_BLACK "\x1b[48;5;232m" +#define BG_RED "\x1b[41m" +#define BG_GREEN "\x1b[42m" +#define BG_YELLOW "\x1b[43m" +#define BG_BLUE "\x1b[44m" +#define BG_MAGENTA "\x1b[45m" +#define BG_CYAN "\x1b[46m" +#define RESET "\x1b[0m" +#define BOLD "\x1b[1m" + +enum TLogLevel { + logERROR, + logWARNING, + logINFOBLUE, + logINFOGREEN, + logINFORED, + logINFOCYAN, + logINFOMAGENTA, + logINFO, + logDEBUG, // constructors, destructors etc. should still give too much + // output + logDEBUG1, + logDEBUG2, + logDEBUG3, + logDEBUG4, + logDEBUG5 +}; + +// Compiler should optimize away anything below this value +#ifndef AARE_LOG_LEVEL +#define AARE_LOG_LEVEL \ + "LOG LEVEL NOT SET IN CMAKE" // This is configured in the main + // CMakeLists.txt +#endif + +#define __AT__ \ + std::string(__FILE__) + std::string("::") + std::string(__func__) + \ + std::string("(): ") +#define __SHORT_FORM_OF_FILE__ \ + (strrchr(__FILE__, '/') ? strrchr(__FILE__, '/') + 1 : __FILE__) +#define __SHORT_AT__ \ + std::string(__SHORT_FORM_OF_FILE__) + std::string("::") + \ + std::string(__func__) + std::string("(): ") + +class Logger { + std::ostringstream os; + TLogLevel m_level = AARE_LOG_LEVEL; + + public: + Logger() = default; + explicit Logger(TLogLevel level) : m_level(level){}; + ~Logger() { + // output in the destructor to allow for << syntax + os << RESET << '\n'; + std::clog << os.str() << std::flush; // Single write + } + + static TLogLevel & + ReportingLevel() { // singelton eeh TODO! Do we need a runtime option? + static TLogLevel reportingLevel = logDEBUG5; + return reportingLevel; + } + + // Danger this buffer need as many elements as TLogLevel + static const char *Color(TLogLevel level) noexcept { + static const char *const colors[] = { + RED BOLD, YELLOW BOLD, BLUE, GREEN, RED, CYAN, MAGENTA, + RESET, RESET, RESET, RESET, RESET, RESET, RESET}; + // out of bounds + if (level < 0 || level >= sizeof(colors) / sizeof(colors[0])) { + return RESET; + } + return colors[level]; + } + + // Danger this buffer need as many elements as TLogLevel + static std::string ToString(TLogLevel level) { + static const char *const buffer[] = { + "ERROR", "WARNING", "INFO", "INFO", "INFO", + "INFO", "INFO", "INFO", "DEBUG", "DEBUG1", + "DEBUG2", "DEBUG3", "DEBUG4", "DEBUG5"}; + // out of bounds + if (level < 0 || level >= sizeof(buffer) / sizeof(buffer[0])) { + return "UNKNOWN"; + } + return buffer[level]; + } + + std::ostringstream &Get() { + os << Color(m_level) << "- " << Timestamp() << " " << ToString(m_level) + << ": "; + return os; + } + + static std::string Timestamp() { + constexpr size_t buffer_len = 12; + char buffer[buffer_len]; + time_t t; + ::time(&t); + tm r; + strftime(buffer, buffer_len, "%X", localtime_r(&t, &r)); + buffer[buffer_len - 1] = '\0'; + struct timeval tv; + gettimeofday(&tv, nullptr); + constexpr size_t result_len = 100; + char result[result_len]; + snprintf(result, result_len, "%s.%03ld", buffer, + static_cast(tv.tv_usec) / 1000); + result[result_len - 1] = '\0'; + return result; + } +}; + +// TODO! Do we need to keep the runtime option? +#define LOG(level) \ + if (level > AARE_LOG_LEVEL) \ + ; \ + else if (level > aare::Logger::ReportingLevel()) \ + ; \ + else \ + aare::Logger(level).Get() + +} // namespace aare diff --git a/include/aare/utils/ifstream_helpers.hpp b/include/aare/utils/ifstream_helpers.hpp index 0a842ed..a8d0d21 100644 --- a/include/aare/utils/ifstream_helpers.hpp +++ b/include/aare/utils/ifstream_helpers.hpp @@ -4,9 +4,9 @@ #include namespace aare { -/** +/** * @brief Get the error message from an ifstream object -*/ + */ std::string ifstream_error_msg(std::ifstream &ifs); } // namespace aare \ No newline at end of file diff --git a/include/aare/utils/par.hpp b/include/aare/utils/par.hpp index efb1c77..e52c897 100644 --- a/include/aare/utils/par.hpp +++ b/include/aare/utils/par.hpp @@ -1,18 +1,18 @@ #include -#include #include +#include namespace aare { - template - void RunInParallel(F func, const std::vector>& tasks) { - // auto tasks = split_task(0, y.shape(0), n_threads); - std::vector threads; - for (auto &task : tasks) { - threads.push_back(std::thread(func, task.first, task.second)); - } - for (auto &thread : threads) { - thread.join(); - } +template +void RunInParallel(F func, const std::vector> &tasks) { + // auto tasks = split_task(0, y.shape(0), n_threads); + std::vector threads; + for (auto &task : tasks) { + threads.push_back(std::thread(func, task.first, task.second)); } + for (auto &thread : threads) { + thread.join(); + } +} } // namespace aare \ No newline at end of file diff --git a/python/aare/ClusterFinder.py b/python/aare/ClusterFinder.py index 6e7c352..251d938 100644 --- a/python/aare/ClusterFinder.py +++ b/python/aare/ClusterFinder.py @@ -1,22 +1,39 @@ - -from ._aare import ClusterFinder_Cluster3x3i, ClusterFinder_Cluster2x2i, ClusterFinderMT_Cluster3x3i, ClusterFinderMT_Cluster2x2i, ClusterCollector_Cluster3x3i, ClusterCollector_Cluster2x2i - - -from ._aare import ClusterFileSink_Cluster3x3i, ClusterFileSink_Cluster2x2i +from . import _aare import numpy as np +_supported_cluster_sizes = [(2,2), (3,3), (5,5), (7,7), (9,9),] + +def _type_to_char(dtype): + if dtype == np.int32: + return 'i' + elif dtype == np.float32: + return 'f' + elif dtype == np.float64: + return 'd' + else: + raise ValueError(f"Unsupported dtype: {dtype}. Only np.int32, np.float32, and np.float64 are supported.") + +def _get_class(name, cluster_size, dtype): + """ + Helper function to get the class based on the name, cluster size, and dtype. + """ + try: + class_name = f"{name}_Cluster{cluster_size[0]}x{cluster_size[1]}{_type_to_char(dtype)}" + cls = getattr(_aare, class_name) + except AttributeError: + raise ValueError(f"Unsupported combination of type and cluster size: {dtype}/{cluster_size} when requesting {class_name}") + return cls + + + def ClusterFinder(image_size, cluster_size, n_sigma=5, dtype = np.int32, capacity = 1024): """ Factory function to create a ClusterFinder object. Provides a cleaner syntax for the templated ClusterFinder in C++. """ - if dtype == np.int32 and cluster_size == (3,3): - return ClusterFinder_Cluster3x3i(image_size, n_sigma = n_sigma, capacity=capacity) - elif dtype == np.int32 and cluster_size == (2,2): - return ClusterFinder_Cluster2x2i(image_size, n_sigma = n_sigma, capacity=capacity) - else: - #TODO! add the other formats - raise ValueError(f"Unsupported dtype: {dtype}. Only np.int32 is supported.") + cls = _get_class("ClusterFinder", cluster_size, dtype) + return cls(image_size, n_sigma=n_sigma, capacity=capacity) + def ClusterFinderMT(image_size, cluster_size = (3,3), dtype=np.int32, n_sigma=5, capacity = 1024, n_threads = 3): @@ -25,15 +42,9 @@ def ClusterFinderMT(image_size, cluster_size = (3,3), dtype=np.int32, n_sigma=5, the templated ClusterFinderMT in C++. """ - if dtype == np.int32 and cluster_size == (3,3): - return ClusterFinderMT_Cluster3x3i(image_size, n_sigma = n_sigma, - capacity = capacity, n_threads = n_threads) - elif dtype == np.int32 and cluster_size == (2,2): - return ClusterFinderMT_Cluster2x2i(image_size, n_sigma = n_sigma, - capacity = capacity, n_threads = n_threads) - else: - #TODO! add the other formats - raise ValueError(f"Unsupported dtype: {dtype}. Only np.int32 is supported.") + cls = _get_class("ClusterFinderMT", cluster_size, dtype) + return cls(image_size, n_sigma=n_sigma, capacity=capacity, n_threads=n_threads) + def ClusterCollector(clusterfindermt, cluster_size = (3,3), dtype=np.int32): @@ -42,14 +53,8 @@ def ClusterCollector(clusterfindermt, cluster_size = (3,3), dtype=np.int32): the templated ClusterCollector in C++. """ - if dtype == np.int32 and cluster_size == (3,3): - return ClusterCollector_Cluster3x3i(clusterfindermt) - elif dtype == np.int32 and cluster_size == (2,2): - return ClusterCollector_Cluster2x2i(clusterfindermt) - - else: - #TODO! add the other formats - raise ValueError(f"Unsupported dtype: {dtype}. Only np.int32 is supported.") + cls = _get_class("ClusterCollector", cluster_size, dtype) + return cls(clusterfindermt) def ClusterFileSink(clusterfindermt, cluster_file, dtype=np.int32): """ @@ -57,11 +62,26 @@ def ClusterFileSink(clusterfindermt, cluster_file, dtype=np.int32): the templated ClusterCollector in C++. """ - if dtype == np.int32 and clusterfindermt.cluster_size == (3,3): - return ClusterFileSink_Cluster3x3i(clusterfindermt, cluster_file) - elif dtype == np.int32 and clusterfindermt.cluster_size == (2,2): - return ClusterFileSink_Cluster2x2i(clusterfindermt, cluster_file) - - else: - #TODO! add the other formats - raise ValueError(f"Unsupported dtype: {dtype}. Only np.int32 is supported.") \ No newline at end of file + cls = _get_class("ClusterFileSink", clusterfindermt.cluster_size, dtype) + return cls(clusterfindermt, cluster_file) + + +def ClusterFile(fname, cluster_size=(3,3), dtype=np.int32, chunk_size = 1000): + """ + Factory function to create a ClusterFile object. Provides a cleaner syntax for + the templated ClusterFile in C++. + + .. 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 + + """ + + cls = _get_class("ClusterFile", cluster_size, dtype) + return cls(fname, chunk_size=chunk_size) diff --git a/python/aare/__init__.py b/python/aare/__init__.py index d2bbe0a..0b95702 100644 --- a/python/aare/__init__.py +++ b/python/aare/__init__.py @@ -5,13 +5,12 @@ from . import _aare from ._aare import File, RawMasterFile, RawSubFile, JungfrauDataFile from ._aare import Pedestal_d, Pedestal_f, ClusterFinder_Cluster3x3i, VarClusterFinder from ._aare import DetectorType -from ._aare import ClusterFile_Cluster3x3i as ClusterFile from ._aare import hitmap from ._aare import ROI # from ._aare import ClusterFinderMT, ClusterCollector, ClusterFileSink, ClusterVector_i -from .ClusterFinder import ClusterFinder, ClusterCollector, ClusterFinderMT, ClusterFileSink +from .ClusterFinder import ClusterFinder, ClusterCollector, ClusterFinderMT, ClusterFileSink, ClusterFile from .ClusterVector import ClusterVector diff --git a/python/examples/play.py b/python/examples/play.py index da469dc..0f4feca 100644 --- a/python/examples/play.py +++ b/python/examples/play.py @@ -1,79 +1,89 @@ import sys sys.path.append('/home/l_msdetect/erik/aare/build') -from aare._aare import ClusterVector_i, Interpolator -import pickle -import numpy as np -import matplotlib.pyplot as plt -import boost_histogram as bh -import torch -import math -import time +from aare import RawSubFile, DetectorType, RawFile + +from pathlib import Path +path = Path("/home/l_msdetect/erik/data/aare-test-data/raw/jungfrau/") +f = RawSubFile(path/"jungfrau_single_d0_f0_0.raw", DetectorType.Jungfrau, 512, 1024, 16) + +# f = RawFile(path/"jungfrau_single_master_0.json") + + +# from aare._aare import ClusterVector_i, Interpolator + +# import pickle +# import numpy as np +# import matplotlib.pyplot as plt +# import boost_histogram as bh +# import torch +# import math +# import time -def gaussian_2d(mx, my, sigma = 1, res=100, grid_size = 2): - """ - Generate a 2D gaussian as position mx, my, with sigma=sigma. - The gaussian is placed on a 2x2 pixel matrix with resolution - res in one dimesion. - """ - x = torch.linspace(0, pixel_size*grid_size, res) - x,y = torch.meshgrid(x,x, indexing="ij") - return 1 / (2*math.pi*sigma**2) * \ - torch.exp(-((x - my)**2 / (2*sigma**2) + (y - mx)**2 / (2*sigma**2))) +# def gaussian_2d(mx, my, sigma = 1, res=100, grid_size = 2): +# """ +# Generate a 2D gaussian as position mx, my, with sigma=sigma. +# The gaussian is placed on a 2x2 pixel matrix with resolution +# res in one dimesion. +# """ +# x = torch.linspace(0, pixel_size*grid_size, res) +# x,y = torch.meshgrid(x,x, indexing="ij") +# return 1 / (2*math.pi*sigma**2) * \ +# torch.exp(-((x - my)**2 / (2*sigma**2) + (y - mx)**2 / (2*sigma**2))) -scale = 1000 #Scale factor when converting to integer -pixel_size = 25 #um -grid = 2 -resolution = 100 -sigma_um = 10 -xa = np.linspace(0,grid*pixel_size,resolution) -ticks = [0, 25, 50] +# scale = 1000 #Scale factor when converting to integer +# pixel_size = 25 #um +# grid = 2 +# resolution = 100 +# sigma_um = 10 +# xa = np.linspace(0,grid*pixel_size,resolution) +# ticks = [0, 25, 50] -hit = np.array((20,20)) -etahist_fname = "/home/l_msdetect/erik/tmp/test_hist.pkl" +# hit = np.array((20,20)) +# etahist_fname = "/home/l_msdetect/erik/tmp/test_hist.pkl" -local_resolution = 99 -grid_size = 3 -xaxis = np.linspace(0,grid_size*pixel_size, local_resolution) -t = gaussian_2d(hit[0],hit[1], grid_size = grid_size, sigma = 10, res = local_resolution) -pixels = t.reshape(grid_size, t.shape[0] // grid_size, grid_size, t.shape[1] // grid_size).sum(axis = 3).sum(axis = 1) -pixels = pixels.numpy() -pixels = (pixels*scale).astype(np.int32) -v = ClusterVector_i(3,3) -v.push_back(1,1, pixels) +# local_resolution = 99 +# grid_size = 3 +# xaxis = np.linspace(0,grid_size*pixel_size, local_resolution) +# t = gaussian_2d(hit[0],hit[1], grid_size = grid_size, sigma = 10, res = local_resolution) +# pixels = t.reshape(grid_size, t.shape[0] // grid_size, grid_size, t.shape[1] // grid_size).sum(axis = 3).sum(axis = 1) +# pixels = pixels.numpy() +# pixels = (pixels*scale).astype(np.int32) +# v = ClusterVector_i(3,3) +# v.push_back(1,1, pixels) -with open(etahist_fname, "rb") as f: - hist = pickle.load(f) -eta = hist.view().copy() -etabinsx = np.array(hist.axes.edges.T[0].flat) -etabinsy = np.array(hist.axes.edges.T[1].flat) -ebins = np.array(hist.axes.edges.T[2].flat) -p = Interpolator(eta, etabinsx[0:-1], etabinsy[0:-1], ebins[0:-1]) +# with open(etahist_fname, "rb") as f: +# hist = pickle.load(f) +# eta = hist.view().copy() +# etabinsx = np.array(hist.axes.edges.T[0].flat) +# etabinsy = np.array(hist.axes.edges.T[1].flat) +# ebins = np.array(hist.axes.edges.T[2].flat) +# p = Interpolator(eta, etabinsx[0:-1], etabinsy[0:-1], ebins[0:-1]) -#Generate the hit +# #Generate the hit -tmp = p.interpolate(v) -print(f'tmp:{tmp}') -pos = np.array((tmp['x'], tmp['y']))*25 +# tmp = p.interpolate(v) +# print(f'tmp:{tmp}') +# pos = np.array((tmp['x'], tmp['y']))*25 -print(pixels) -fig, ax = plt.subplots(figsize = (7,7)) -ax.pcolormesh(xaxis, xaxis, t) -ax.plot(*pos, 'o') -ax.set_xticks([0,25,50,75]) -ax.set_yticks([0,25,50,75]) -ax.set_xlim(0,75) -ax.set_ylim(0,75) -ax.grid() -print(f'{hit=}') -print(f'{pos=}') \ No newline at end of file +# print(pixels) +# fig, ax = plt.subplots(figsize = (7,7)) +# ax.pcolormesh(xaxis, xaxis, t) +# ax.plot(*pos, 'o') +# ax.set_xticks([0,25,50,75]) +# ax.set_yticks([0,25,50,75]) +# ax.set_xlim(0,75) +# ax.set_ylim(0,75) +# ax.grid() +# print(f'{hit=}') +# print(f'{pos=}') \ No newline at end of file diff --git a/python/src/bind_Cluster.hpp b/python/src/bind_Cluster.hpp new file mode 100644 index 0000000..690d0e8 --- /dev/null +++ b/python/src/bind_Cluster.hpp @@ -0,0 +1,64 @@ +#include "aare/Cluster.hpp" + +#include +#include +#include +#include +#include +#include +#include + +namespace py = pybind11; +using pd_type = double; + +using namespace aare; + +#pragma GCC diagnostic push +#pragma GCC diagnostic ignored "-Wunused-parameter" + +template +void define_Cluster(py::module &m, const std::string &typestr) { + auto class_name = fmt::format("Cluster{}", typestr); + + py::class_>( + m, class_name.c_str(), py::buffer_protocol()) + + .def(py::init([](uint8_t x, uint8_t y, py::array_t data) { + py::buffer_info buf_info = data.request(); + Cluster cluster; + cluster.x = x; + cluster.y = y; + auto r = data.template unchecked<1>(); // no bounds checks + for (py::ssize_t i = 0; i < data.size(); ++i) { + cluster.data[i] = r(i); + } + return cluster; + })); + + /* + //TODO! Review if to keep or not + .def_property( + "data", + [](ClusterType &c) -> py::array { + return py::array(py::buffer_info( + c.data, sizeof(Type), + py::format_descriptor::format(), // Type + // format + 1, // Number of dimensions + {static_cast(ClusterSizeX * + ClusterSizeY)}, // Shape (flattened) + {sizeof(Type)} // Stride (step size between elements) + )); + }, + [](ClusterType &c, py::array_t arr) { + py::buffer_info buf_info = arr.request(); + Type *ptr = static_cast(buf_info.ptr); + std::copy(ptr, ptr + ClusterSizeX * ClusterSizeY, + c.data); // TODO dont iterate over centers!!! + + }); + */ +} + +#pragma GCC diagnostic pop \ No newline at end of file diff --git a/python/src/bind_ClusterCollector.hpp b/python/src/bind_ClusterCollector.hpp new file mode 100644 index 0000000..84172cb --- /dev/null +++ b/python/src/bind_ClusterCollector.hpp @@ -0,0 +1,44 @@ +#include "aare/ClusterCollector.hpp" +#include "aare/ClusterFileSink.hpp" +#include "aare/ClusterFinder.hpp" +#include "aare/ClusterFinderMT.hpp" +#include "aare/ClusterVector.hpp" +#include "aare/NDView.hpp" +#include "aare/Pedestal.hpp" +#include "np_helper.hpp" + +#include +#include +#include +#include +#include + +namespace py = pybind11; +using pd_type = double; + +using namespace aare; + +#pragma GCC diagnostic push +#pragma GCC diagnostic ignored "-Wunused-parameter" + +template +void define_ClusterCollector(py::module &m, const std::string &typestr) { + auto class_name = fmt::format("ClusterCollector_{}", typestr); + + using ClusterType = Cluster; + + py::class_>(m, class_name.c_str()) + .def(py::init *>()) + .def("stop", &ClusterCollector::stop) + .def( + "steal_clusters", + [](ClusterCollector &self) { + auto v = new std::vector>( + self.steal_clusters()); + return v; // TODO change!!! + }, + py::return_value_policy::take_ownership); +} + +#pragma GCC diagnostic pop \ No newline at end of file diff --git a/python/src/cluster_file.hpp b/python/src/bind_ClusterFile.hpp similarity index 92% rename from python/src/cluster_file.hpp rename to python/src/bind_ClusterFile.hpp index ac384b2..5d8aa88 100644 --- a/python/src/cluster_file.hpp +++ b/python/src/bind_ClusterFile.hpp @@ -21,8 +21,7 @@ using namespace ::aare; template -void define_cluster_file_io_bindings(py::module &m, - const std::string &typestr) { +void define_ClusterFile(py::module &m, const std::string &typestr) { using ClusterType = Cluster; @@ -39,19 +38,20 @@ void define_cluster_file_io_bindings(py::module &m, self.read_clusters(n_clusters)); return v; }, - py::return_value_policy::take_ownership) + 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) + .def("set_roi", &ClusterFile::set_roi, + py::arg("roi")) .def( "set_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")) .def("set_gain_map", [](ClusterFile &self, py::array_t gain_map) { diff --git a/python/src/bind_ClusterFileSink.hpp b/python/src/bind_ClusterFileSink.hpp new file mode 100644 index 0000000..f717de6 --- /dev/null +++ b/python/src/bind_ClusterFileSink.hpp @@ -0,0 +1,37 @@ +#include "aare/ClusterCollector.hpp" +#include "aare/ClusterFileSink.hpp" +#include "aare/ClusterFinder.hpp" +#include "aare/ClusterFinderMT.hpp" +#include "aare/ClusterVector.hpp" +#include "aare/NDView.hpp" +#include "aare/Pedestal.hpp" +#include "np_helper.hpp" + +#include +#include +#include +#include +#include + +namespace py = pybind11; +using pd_type = double; + +using namespace aare; + +#pragma GCC diagnostic push +#pragma GCC diagnostic ignored "-Wunused-parameter" + +template +void define_ClusterFileSink(py::module &m, const std::string &typestr) { + auto class_name = fmt::format("ClusterFileSink_{}", typestr); + + using ClusterType = Cluster; + + py::class_>(m, class_name.c_str()) + .def(py::init *, + const std::filesystem::path &>()) + .def("stop", &ClusterFileSink::stop); +} + +#pragma GCC diagnostic pop diff --git a/python/src/bind_ClusterFinder.hpp b/python/src/bind_ClusterFinder.hpp new file mode 100644 index 0000000..5f0fe8d --- /dev/null +++ b/python/src/bind_ClusterFinder.hpp @@ -0,0 +1,77 @@ +#include "aare/ClusterCollector.hpp" +#include "aare/ClusterFileSink.hpp" +#include "aare/ClusterFinder.hpp" +#include "aare/ClusterFinderMT.hpp" +#include "aare/ClusterVector.hpp" +#include "aare/NDView.hpp" +#include "aare/Pedestal.hpp" +#include "np_helper.hpp" + +#include +#include +#include +#include +#include + +namespace py = pybind11; +using pd_type = double; + +using namespace aare; + +#pragma GCC diagnostic push +#pragma GCC diagnostic ignored "-Wunused-parameter" + +template +void define_ClusterFinder(py::module &m, const std::string &typestr) { + auto class_name = fmt::format("ClusterFinder_{}", typestr); + + using ClusterType = Cluster; + + py::class_>( + m, class_name.c_str()) + .def(py::init, pd_type, size_t>(), py::arg("image_size"), + py::arg("n_sigma") = 5.0, py::arg("capacity") = 1'000'000) + .def("push_pedestal_frame", + [](ClusterFinder &self, + py::array_t frame) { + auto view = make_view_2d(frame); + self.push_pedestal_frame(view); + }) + .def("clear_pedestal", + &ClusterFinder::clear_pedestal) + .def_property_readonly( + "pedestal", + [](ClusterFinder &self) { + auto pd = new NDArray{}; + *pd = self.pedestal(); + return return_image_data(pd); + }) + .def_property_readonly( + "noise", + [](ClusterFinder &self) { + auto arr = new NDArray{}; + *arr = self.noise(); + return return_image_data(arr); + }) + .def( + "steal_clusters", + [](ClusterFinder &self, + bool realloc_same_capacity) { + ClusterVector clusters = + self.steal_clusters(realloc_same_capacity); + return clusters; + }, + py::arg("realloc_same_capacity") = false) + .def( + "find_clusters", + [](ClusterFinder &self, + py::array_t frame, uint64_t frame_number) { + auto view = make_view_2d(frame); + self.find_clusters(view, frame_number); + return; + }, + py::arg(), py::arg("frame_number") = 0); +} + +#pragma GCC diagnostic pop \ No newline at end of file diff --git a/python/src/bind_ClusterFinderMT.hpp b/python/src/bind_ClusterFinderMT.hpp new file mode 100644 index 0000000..0ecbbd1 --- /dev/null +++ b/python/src/bind_ClusterFinderMT.hpp @@ -0,0 +1,81 @@ +#include "aare/ClusterCollector.hpp" +#include "aare/ClusterFileSink.hpp" +#include "aare/ClusterFinder.hpp" +#include "aare/ClusterFinderMT.hpp" +#include "aare/ClusterVector.hpp" +#include "aare/NDView.hpp" +#include "aare/Pedestal.hpp" +#include "np_helper.hpp" + +#include +#include +#include +#include +#include + +namespace py = pybind11; +using pd_type = double; + +using namespace aare; + +#pragma GCC diagnostic push +#pragma GCC diagnostic ignored "-Wunused-parameter" + +template +void define_ClusterFinderMT(py::module &m, const std::string &typestr) { + auto class_name = fmt::format("ClusterFinderMT_{}", typestr); + + using ClusterType = Cluster; + + py::class_>( + m, class_name.c_str()) + .def(py::init, pd_type, size_t, size_t>(), + py::arg("image_size"), py::arg("n_sigma") = 5.0, + py::arg("capacity") = 2048, py::arg("n_threads") = 3) + .def("push_pedestal_frame", + [](ClusterFinderMT &self, + py::array_t frame) { + auto view = make_view_2d(frame); + self.push_pedestal_frame(view); + }) + .def( + "find_clusters", + [](ClusterFinderMT &self, + py::array_t frame, uint64_t frame_number) { + auto view = make_view_2d(frame); + self.find_clusters(view, frame_number); + return; + }, + py::arg(), py::arg("frame_number") = 0) + .def_property_readonly( + "cluster_size", + [](ClusterFinderMT &self) { + return py::make_tuple(ClusterSizeX, ClusterSizeY); + }) + .def("clear_pedestal", + &ClusterFinderMT::clear_pedestal) + .def("sync", &ClusterFinderMT::sync) + .def("stop", &ClusterFinderMT::stop) + .def("start", &ClusterFinderMT::start) + .def( + "pedestal", + [](ClusterFinderMT &self, + size_t thread_index) { + auto pd = new NDArray{}; + *pd = self.pedestal(thread_index); + return return_image_data(pd); + }, + py::arg("thread_index") = 0) + .def( + "noise", + [](ClusterFinderMT &self, + size_t thread_index) { + auto arr = new NDArray{}; + *arr = self.noise(thread_index); + return return_image_data(arr); + }, + py::arg("thread_index") = 0); +} + +#pragma GCC diagnostic pop \ No newline at end of file diff --git a/python/src/bind_ClusterVector.hpp b/python/src/bind_ClusterVector.hpp index db8c8a3..9e9c4ab 100644 --- a/python/src/bind_ClusterVector.hpp +++ b/python/src/bind_ClusterVector.hpp @@ -44,10 +44,11 @@ void define_ClusterVector(py::module &m, const std::string &typestr) { auto *vec = new std::vector(self.sum()); return return_vector(vec); }) - .def("sum_2x2", [](ClusterVector &self){ - auto *vec = new std::vector(self.sum_2x2()); - return return_vector(vec); - }) + .def("sum_2x2", + [](ClusterVector &self) { + auto *vec = new std::vector(self.sum_2x2()); + return return_vector(vec); + }) .def_property_readonly("size", &ClusterVector::size) .def("item_size", &ClusterVector::item_size) .def_property_readonly("fmt", @@ -101,4 +102,6 @@ void define_ClusterVector(py::module &m, const std::string &typestr) { return hitmap; }); -} \ No newline at end of file +} + +#pragma GCC diagnostic pop \ No newline at end of file diff --git a/python/src/cluster.hpp b/python/src/cluster.hpp deleted file mode 100644 index 58f137c..0000000 --- a/python/src/cluster.hpp +++ /dev/null @@ -1,211 +0,0 @@ -#include "aare/ClusterCollector.hpp" -#include "aare/ClusterFileSink.hpp" -#include "aare/ClusterFinder.hpp" -#include "aare/ClusterFinderMT.hpp" -#include "aare/ClusterVector.hpp" -#include "aare/NDView.hpp" -#include "aare/Pedestal.hpp" -#include "np_helper.hpp" - -#include -#include -#include -#include -#include - -namespace py = pybind11; -using pd_type = double; - -using namespace aare; - -#pragma GCC diagnostic push -#pragma GCC diagnostic ignored "-Wunused-parameter" - -template -void define_cluster(py::module &m, const std::string &typestr) { - auto class_name = fmt::format("Cluster{}", typestr); - - py::class_>( - m, class_name.c_str(), py::buffer_protocol()) - - .def(py::init([](uint8_t x, uint8_t y, py::array_t data) { - py::buffer_info buf_info = data.request(); - Cluster cluster; - cluster.x = x; - cluster.y = y; - auto r = data.template unchecked<1>(); // no bounds checks - for (py::ssize_t i = 0; i < data.size(); ++i) { - cluster.data[i] = r(i); - } - return cluster; - })); - - /* - .def_property( - "data", - [](ClusterType &c) -> py::array { - return py::array(py::buffer_info( - c.data, sizeof(Type), - py::format_descriptor::format(), // Type - // format - 1, // Number of dimensions - {static_cast(ClusterSizeX * - ClusterSizeY)}, // Shape (flattened) - {sizeof(Type)} // Stride (step size between elements) - )); - }, - [](ClusterType &c, py::array_t arr) { - py::buffer_info buf_info = arr.request(); - Type *ptr = static_cast(buf_info.ptr); - std::copy(ptr, ptr + ClusterSizeX * ClusterSizeY, - c.data); // TODO dont iterate over centers!!! - - }); - */ -} - -template -void define_cluster_finder_mt_bindings(py::module &m, - const std::string &typestr) { - auto class_name = fmt::format("ClusterFinderMT_{}", typestr); - - using ClusterType = Cluster; - - py::class_>( - m, class_name.c_str()) - .def(py::init, pd_type, size_t, size_t>(), - py::arg("image_size"), py::arg("n_sigma") = 5.0, - py::arg("capacity") = 2048, py::arg("n_threads") = 3) - .def("push_pedestal_frame", - [](ClusterFinderMT &self, - py::array_t frame) { - auto view = make_view_2d(frame); - self.push_pedestal_frame(view); - }) - .def( - "find_clusters", - [](ClusterFinderMT &self, - py::array_t frame, uint64_t frame_number) { - auto view = make_view_2d(frame); - self.find_clusters(view, frame_number); - return; - }, - py::arg(), py::arg("frame_number") = 0) - .def_property_readonly("cluster_size", [](ClusterFinderMT &self){ - return py::make_tuple(ClusterSizeX, ClusterSizeY); - }) - .def("clear_pedestal", - &ClusterFinderMT::clear_pedestal) - .def("sync", &ClusterFinderMT::sync) - .def("stop", &ClusterFinderMT::stop) - .def("start", &ClusterFinderMT::start) - .def( - "pedestal", - [](ClusterFinderMT &self, - size_t thread_index) { - auto pd = new NDArray{}; - *pd = self.pedestal(thread_index); - return return_image_data(pd); - }, - py::arg("thread_index") = 0) - .def( - "noise", - [](ClusterFinderMT &self, - size_t thread_index) { - auto arr = new NDArray{}; - *arr = self.noise(thread_index); - return return_image_data(arr); - }, - py::arg("thread_index") = 0); -} - -template -void define_cluster_collector_bindings(py::module &m, - const std::string &typestr) { - auto class_name = fmt::format("ClusterCollector_{}", typestr); - - using ClusterType = Cluster; - - py::class_>(m, class_name.c_str()) - .def(py::init *>()) - .def("stop", &ClusterCollector::stop) - .def( - "steal_clusters", - [](ClusterCollector &self) { - auto v = new std::vector>( - self.steal_clusters()); - return v; // TODO change!!! - }, - py::return_value_policy::take_ownership); -} - -template -void define_cluster_file_sink_bindings(py::module &m, - const std::string &typestr) { - auto class_name = fmt::format("ClusterFileSink_{}", typestr); - - using ClusterType = Cluster; - - py::class_>(m, class_name.c_str()) - .def(py::init *, - const std::filesystem::path &>()) - .def("stop", &ClusterFileSink::stop); -} - -template -void define_cluster_finder_bindings(py::module &m, const std::string &typestr) { - auto class_name = fmt::format("ClusterFinder_{}", typestr); - - using ClusterType = Cluster; - - py::class_>( - m, class_name.c_str()) - .def(py::init, pd_type, size_t>(), py::arg("image_size"), - py::arg("n_sigma") = 5.0, py::arg("capacity") = 1'000'000) - .def("push_pedestal_frame", - [](ClusterFinder &self, - py::array_t frame) { - auto view = make_view_2d(frame); - self.push_pedestal_frame(view); - }) - .def("clear_pedestal", - &ClusterFinder::clear_pedestal) - .def_property_readonly( - "pedestal", - [](ClusterFinder &self) { - auto pd = new NDArray{}; - *pd = self.pedestal(); - return return_image_data(pd); - }) - .def_property_readonly( - "noise", - [](ClusterFinder &self) { - auto arr = new NDArray{}; - *arr = self.noise(); - return return_image_data(arr); - }) - .def( - "steal_clusters", - [](ClusterFinder &self, - bool realloc_same_capacity) { - ClusterVector clusters = - self.steal_clusters(realloc_same_capacity); - return clusters; - }, - py::arg("realloc_same_capacity") = false) - .def( - "find_clusters", - [](ClusterFinder &self, - py::array_t frame, uint64_t frame_number) { - auto view = make_view_2d(frame); - self.find_clusters(view, frame_number); - return; - }, - py::arg(), py::arg("frame_number") = 0); -} -#pragma GCC diagnostic pop diff --git a/python/src/ctb_raw_file.hpp b/python/src/ctb_raw_file.hpp index c9b5310..5eb9652 100644 --- a/python/src/ctb_raw_file.hpp +++ b/python/src/ctb_raw_file.hpp @@ -6,8 +6,8 @@ #include "aare/RawMasterFile.hpp" #include "aare/RawSubFile.hpp" -#include "aare/defs.hpp" #include "aare/decode.hpp" +#include "aare/defs.hpp" // #include "aare/fClusterFileV2.hpp" #include "np_helper.hpp" @@ -26,95 +26,103 @@ using namespace ::aare; void define_ctb_raw_file_io_bindings(py::module &m) { -m.def("adc_sar_05_decode64to16", [](py::array_t input) { + m.def("adc_sar_05_decode64to16", [](py::array_t input) { + if (input.ndim() != 2) { + throw std::runtime_error( + "Only 2D arrays are supported at this moment"); + } - - if(input.ndim() != 2){ - throw std::runtime_error("Only 2D arrays are supported at this moment"); - } + // Create a 2D output array with the same shape as the input + std::vector shape{input.shape(0), + input.shape(1) / + static_cast(bits_per_byte)}; + py::array_t output(shape); - //Create a 2D output array with the same shape as the input - std::vector shape{input.shape(0), input.shape(1)/static_cast(bits_per_byte)}; - py::array_t output(shape); + // Create a view of the input and output arrays + NDView input_view( + reinterpret_cast(input.mutable_data()), + {output.shape(0), output.shape(1)}); + NDView output_view(output.mutable_data(), + {output.shape(0), output.shape(1)}); - //Create a view of the input and output arrays - NDView input_view(reinterpret_cast(input.mutable_data()), {output.shape(0), output.shape(1)}); - NDView output_view(output.mutable_data(), {output.shape(0), output.shape(1)}); + adc_sar_05_decode64to16(input_view, output_view); - adc_sar_05_decode64to16(input_view, output_view); - - return output; -}); - - -m.def("adc_sar_04_decode64to16", [](py::array_t input) { - - - if(input.ndim() != 2){ - throw std::runtime_error("Only 2D arrays are supported at this moment"); - } - - //Create a 2D output array with the same shape as the input - std::vector shape{input.shape(0), input.shape(1)/static_cast(bits_per_byte)}; - py::array_t output(shape); - - //Create a view of the input and output arrays - NDView input_view(reinterpret_cast(input.mutable_data()), {output.shape(0), output.shape(1)}); - NDView output_view(output.mutable_data(), {output.shape(0), output.shape(1)}); - - adc_sar_04_decode64to16(input_view, output_view); - - return output; -}); - -m.def( - "apply_custom_weights", - [](py::array_t &input, - py::array_t - &weights) { - - - // Create new array with same shape as the input array (uninitialized values) - py::buffer_info buf = input.request(); - py::array_t output(buf.shape); - - // Use NDViews to call into the C++ library - auto weights_view = make_view_1d(weights); - NDView input_view(input.mutable_data(), {input.size()}); - NDView output_view(output.mutable_data(), {output.size()}); - - apply_custom_weights(input_view, output_view, weights_view); return output; }); -py::class_(m, "CtbRawFile") - .def(py::init()) - .def("read_frame", - [](CtbRawFile &self) { - size_t image_size = self.image_size_in_bytes(); - py::array image; - std::vector shape; - shape.reserve(2); - shape.push_back(1); - shape.push_back(image_size); + m.def("adc_sar_04_decode64to16", [](py::array_t input) { + if (input.ndim() != 2) { + throw std::runtime_error( + "Only 2D arrays are supported at this moment"); + } - py::array_t header(1); + // Create a 2D output array with the same shape as the input + std::vector shape{input.shape(0), + input.shape(1) / + static_cast(bits_per_byte)}; + py::array_t output(shape); - // always read bytes - image = py::array_t(shape); + // Create a view of the input and output arrays + NDView input_view( + reinterpret_cast(input.mutable_data()), + {output.shape(0), output.shape(1)}); + NDView output_view(output.mutable_data(), + {output.shape(0), output.shape(1)}); - self.read_into(reinterpret_cast(image.mutable_data()), - header.mutable_data()); + adc_sar_04_decode64to16(input_view, output_view); - return py::make_tuple(header, image); - }) - .def("seek", &CtbRawFile::seek) - .def("tell", &CtbRawFile::tell) - .def("master", &CtbRawFile::master) + return output; + }); - .def_property_readonly("image_size_in_bytes", - &CtbRawFile::image_size_in_bytes) + m.def("apply_custom_weights", + [](py::array_t + &input, + py::array_t + &weights) { + // Create new array with same shape as the input array + // (uninitialized values) + py::buffer_info buf = input.request(); + py::array_t output(buf.shape); - .def_property_readonly("frames_in_file", &CtbRawFile::frames_in_file); + // Use NDViews to call into the C++ library + auto weights_view = make_view_1d(weights); + NDView input_view(input.mutable_data(), + {input.size()}); + NDView output_view(output.mutable_data(), + {output.size()}); + apply_custom_weights(input_view, output_view, weights_view); + return output; + }); + + py::class_(m, "CtbRawFile") + .def(py::init()) + .def("read_frame", + [](CtbRawFile &self) { + size_t image_size = self.image_size_in_bytes(); + py::array image; + std::vector shape; + shape.reserve(2); + shape.push_back(1); + shape.push_back(image_size); + + py::array_t header(1); + + // always read bytes + image = py::array_t(shape); + + self.read_into( + reinterpret_cast(image.mutable_data()), + header.mutable_data()); + + return py::make_tuple(header, image); + }) + .def("seek", &CtbRawFile::seek) + .def("tell", &CtbRawFile::tell) + .def("master", &CtbRawFile::master) + + .def_property_readonly("image_size_in_bytes", + &CtbRawFile::image_size_in_bytes) + + .def_property_readonly("frames_in_file", &CtbRawFile::frames_in_file); } diff --git a/python/src/file.hpp b/python/src/file.hpp index f97db96..262b4f8 100644 --- a/python/src/file.hpp +++ b/python/src/file.hpp @@ -20,17 +20,13 @@ namespace py = pybind11; using namespace ::aare; - - - -//Disable warnings for unused parameters, as we ignore some -//in the __exit__ method +// Disable warnings for unused parameters, as we ignore some +// in the __exit__ method #pragma GCC diagnostic push #pragma GCC diagnostic ignored "-Wunused-parameter" void define_file_io_bindings(py::module &m) { - py::enum_(m, "DetectorType") .value("Jungfrau", DetectorType::Jungfrau) .value("Eiger", DetectorType::Eiger) @@ -41,13 +37,10 @@ void define_file_io_bindings(py::module &m) { .value("ChipTestBoard", DetectorType::ChipTestBoard) .value("Unknown", DetectorType::Unknown); - PYBIND11_NUMPY_DTYPE(DetectorHeader, frameNumber, expLength, packetNumber, bunchId, timestamp, modId, row, column, reserved, debug, roundRNumber, detType, version, packetMask); - - py::class_(m, "File") .def(py::init([](const std::filesystem::path &fname) { return File(fname, "r", {}); @@ -112,45 +105,18 @@ void define_file_io_bindings(py::module &m) { reinterpret_cast(image.mutable_data())); return image; }) - .def("read_n", [](File &self, size_t n_frames) { - //adjust for actual frames left in the file - n_frames = std::min(n_frames, self.total_frames()-self.tell()); - if(n_frames == 0){ - throw std::runtime_error("No frames left in file"); - } - std::vector shape{n_frames, self.rows(), self.cols()}; - - py::array image; - const uint8_t item_size = self.bytes_per_pixel(); - if (item_size == 1) { - image = py::array_t(shape); - } else if (item_size == 2) { - image = py::array_t(shape); - } else if (item_size == 4) { - image = py::array_t(shape); - } - self.read_into(reinterpret_cast(image.mutable_data()), - n_frames); - return image; - }) - .def("__enter__", [](File &self) { return &self; }) - .def("__exit__", - [](File &self, - const std::optional &exc_type, - const std::optional &exc_value, - const std::optional &traceback) { - // self.close(); - }) - .def("__iter__", [](File &self) { return &self; }) - .def("__next__", [](File &self) { + .def("read_n", + [](File &self, size_t n_frames) { + // adjust for actual frames left in the file + n_frames = + std::min(n_frames, self.total_frames() - self.tell()); + if (n_frames == 0) { + throw std::runtime_error("No frames left in file"); + } + std::vector shape{n_frames, self.rows(), self.cols()}; - try{ - const uint8_t item_size = self.bytes_per_pixel(); py::array image; - std::vector shape; - shape.reserve(2); - shape.push_back(self.rows()); - shape.push_back(self.cols()); + const uint8_t item_size = self.bytes_per_pixel(); if (item_size == 1) { image = py::array_t(shape); } else if (item_size == 2) { @@ -159,14 +125,41 @@ void define_file_io_bindings(py::module &m) { image = py::array_t(shape); } self.read_into( - reinterpret_cast(image.mutable_data())); + reinterpret_cast(image.mutable_data()), + n_frames); return image; - }catch(std::runtime_error &e){ + }) + .def("__enter__", [](File &self) { return &self; }) + .def("__exit__", + [](File &self, const std::optional &exc_type, + const std::optional &exc_value, + const std::optional &traceback) { + // self.close(); + }) + .def("__iter__", [](File &self) { return &self; }) + .def("__next__", [](File &self) { + try { + const uint8_t item_size = self.bytes_per_pixel(); + py::array image; + std::vector shape; + shape.reserve(2); + shape.push_back(self.rows()); + shape.push_back(self.cols()); + if (item_size == 1) { + image = py::array_t(shape); + } else if (item_size == 2) { + image = py::array_t(shape); + } else if (item_size == 4) { + image = py::array_t(shape); + } + self.read_into( + reinterpret_cast(image.mutable_data())); + return image; + } catch (std::runtime_error &e) { throw py::stop_iteration(); } }); - py::class_(m, "FileConfig") .def(py::init<>()) .def_readwrite("rows", &FileConfig::rows) @@ -183,8 +176,6 @@ void define_file_io_bindings(py::module &m) { return ""; }); - - py::class_(m, "ScanParameters") .def(py::init()) .def(py::init()) @@ -195,7 +186,6 @@ void define_file_io_bindings(py::module &m) { .def_property_readonly("stop", &ScanParameters::stop) .def_property_readonly("step", &ScanParameters::step); - py::class_(m, "ROI") .def(py::init<>()) .def(py::init(), py::arg("xmin"), @@ -204,23 +194,21 @@ void define_file_io_bindings(py::module &m) { .def_readwrite("xmax", &ROI::xmax) .def_readwrite("ymin", &ROI::ymin) .def_readwrite("ymax", &ROI::ymax) - .def("__str__", [](const ROI& self){ - return fmt::format("ROI: xmin: {} xmax: {} ymin: {} ymax: {}", self.xmin, self.xmax, self.ymin, self.ymax); - }) - .def("__repr__", [](const ROI& self){ - return fmt::format("", self.xmin, self.xmax, self.ymin, self.ymax); - }) + .def("__str__", + [](const ROI &self) { + return fmt::format("ROI: xmin: {} xmax: {} ymin: {} ymax: {}", + self.xmin, self.xmax, self.ymin, self.ymax); + }) + .def("__repr__", + [](const ROI &self) { + return fmt::format( + "", self.xmin, + self.xmax, self.ymin, self.ymax); + }) .def("__iter__", [](const ROI &self) { - return py::make_iterator(&self.xmin, &self.ymax+1); //NOLINT + return py::make_iterator(&self.xmin, &self.ymax + 1); // NOLINT }); - - - - - - - #pragma GCC diagnostic pop // py::class_(m, "ClusterHeader") // .def(py::init<>()) diff --git a/python/src/fit.hpp b/python/src/fit.hpp index 97dafb5..47568d6 100644 --- a/python/src/fit.hpp +++ b/python/src/fit.hpp @@ -9,7 +9,6 @@ namespace py = pybind11; using namespace pybind11::literals; - void define_fit_bindings(py::module &m) { // TODO! Evaluate without converting to double @@ -61,7 +60,8 @@ void define_fit_bindings(py::module &m) { py::array_t par) { auto x_view = make_view_1d(x); auto par_view = make_view_1d(par); - auto y = new NDArray{aare::func::scurve(x_view, par_view)}; + auto y = + new NDArray{aare::func::scurve(x_view, par_view)}; return return_image_data(y); }, R"( @@ -82,7 +82,8 @@ void define_fit_bindings(py::module &m) { py::array_t par) { auto x_view = make_view_1d(x); auto par_view = make_view_1d(par); - auto y = new NDArray{aare::func::scurve2(x_view, par_view)}; + auto y = + new NDArray{aare::func::scurve2(x_view, par_view)}; return return_image_data(y); }, R"( @@ -139,7 +140,6 @@ n_threads : int, optional py::array_t y, py::array_t y_err, int n_threads) { - if (y.ndim() == 3) { // Allocate memory for the output // Need to have pointers to allow python to manage @@ -173,7 +173,6 @@ n_threads : int, optional 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); @@ -248,11 +247,10 @@ n_threads : int, optional 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), + "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({2}); auto par_err = new NDArray({2}); @@ -289,7 +287,7 @@ n_threads : int, optional )", py::arg("x"), py::arg("y"), py::arg("y_err"), py::arg("n_threads") = 4); -//========= + //========= m.def( "fit_scurve", [](py::array_t x, @@ -333,13 +331,12 @@ n_threads : int, optional auto chi2 = new NDArray({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); + par_err->view(), chi2->view(), n_threads); return py::dict("par"_a = return_image_data(par), - "par_err"_a = return_image_data(par_err), + "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({2}); auto par_err = new NDArray({2}); @@ -351,7 +348,7 @@ n_threads : int, optional double chi2 = 0; aare::fit_scurve(x_view, y_view, y_view_err, par->view(), - par_err->view(), chi2); + 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); @@ -375,7 +372,6 @@ 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", @@ -420,13 +416,12 @@ n_threads : int, optional auto chi2 = new NDArray({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); + par_err->view(), chi2->view(), n_threads); return py::dict("par"_a = return_image_data(par), - "par_err"_a = return_image_data(par_err), + "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({6}); auto par_err = new NDArray({6}); @@ -438,7 +433,7 @@ n_threads : int, optional double chi2 = 0; aare::fit_scurve2(x_view, y_view, y_view_err, par->view(), - par_err->view(), chi2); + 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); diff --git a/python/src/jungfrau_data_file.hpp b/python/src/jungfrau_data_file.hpp index 942f6a6..62a95c9 100644 --- a/python/src/jungfrau_data_file.hpp +++ b/python/src/jungfrau_data_file.hpp @@ -21,10 +21,7 @@ using namespace ::aare; auto read_dat_frame(JungfrauDataFile &self) { py::array_t header(1); - py::array_t image({ - self.rows(), - self.cols() - }); + py::array_t image({self.rows(), self.cols()}); self.read_into(reinterpret_cast(image.mutable_data()), header.mutable_data()); @@ -40,9 +37,7 @@ auto read_n_dat_frames(JungfrauDataFile &self, size_t n_frames) { } py::array_t header(n_frames); - py::array_t image({ - n_frames, self.rows(), - self.cols()}); + py::array_t image({n_frames, self.rows(), self.cols()}); self.read_into(reinterpret_cast(image.mutable_data()), n_frames, header.mutable_data()); diff --git a/python/src/module.cpp b/python/src/module.cpp index 946a41b..fc04a9f 100644 --- a/python/src/module.cpp +++ b/python/src/module.cpp @@ -1,22 +1,26 @@ // Files with bindings to the different classes -//New style file naming +// New style file naming +#include "bind_Cluster.hpp" +#include "bind_ClusterCollector.hpp" +#include "bind_ClusterFile.hpp" +#include "bind_ClusterFileSink.hpp" +#include "bind_ClusterFinder.hpp" +#include "bind_ClusterFinderMT.hpp" #include "bind_ClusterVector.hpp" -//TODO! migrate the other names -#include "cluster.hpp" -#include "cluster_file.hpp" +// TODO! migrate the other names #include "ctb_raw_file.hpp" #include "file.hpp" #include "fit.hpp" #include "interpolation.hpp" -#include "raw_sub_file.hpp" -#include "raw_master_file.hpp" -#include "raw_file.hpp" -#include "pixel_map.hpp" -#include "var_cluster.hpp" -#include "pedestal.hpp" #include "jungfrau_data_file.hpp" +#include "pedestal.hpp" +#include "pixel_map.hpp" +#include "raw_file.hpp" +#include "raw_master_file.hpp" +#include "raw_sub_file.hpp" +#include "var_cluster.hpp" // Pybind stuff #include @@ -24,6 +28,26 @@ namespace py = pybind11; +/* MACRO that defines Cluster bindings for a specific size and type + +T - Storage type of the cluster data (int, float, double) +N - Number of rows in the cluster +M - Number of columns in the cluster +U - Type of the pixel data (e.g., uint16_t) +TYPE_CODE - A character representing the type code (e.g., 'i' for int, 'd' for +double, 'f' for float) + +*/ +#define DEFINE_CLUSTER_BINDINGS(T, N, M, U, TYPE_CODE) \ + define_ClusterFile(m, "Cluster" #N "x" #M #TYPE_CODE); \ + define_ClusterVector(m, "Cluster" #N "x" #M #TYPE_CODE); \ + define_ClusterFinder(m, "Cluster" #N "x" #M #TYPE_CODE); \ + define_ClusterFinderMT(m, "Cluster" #N "x" #M #TYPE_CODE); \ + define_ClusterFileSink(m, "Cluster" #N "x" #M #TYPE_CODE); \ + define_ClusterCollector(m, "Cluster" #N "x" #M #TYPE_CODE); \ + define_Cluster(m, #N "x" #M #TYPE_CODE); \ + register_calculate_eta(m); + PYBIND11_MODULE(_aare, m) { define_file_io_bindings(m); define_raw_file_io_bindings(m); @@ -38,59 +62,23 @@ PYBIND11_MODULE(_aare, m) { define_interpolation_bindings(m); define_jungfrau_data_file_io_bindings(m); - define_cluster_file_io_bindings(m, "Cluster3x3i"); - define_cluster_file_io_bindings(m, "Cluster3x3d"); - define_cluster_file_io_bindings(m, "Cluster3x3f"); - define_cluster_file_io_bindings(m, "Cluster2x2i"); - define_cluster_file_io_bindings(m, "Cluster2x2f"); - define_cluster_file_io_bindings(m, "Cluster2x2d"); + DEFINE_CLUSTER_BINDINGS(int, 3, 3, uint16_t, i); + DEFINE_CLUSTER_BINDINGS(double, 3, 3, uint16_t, d); + DEFINE_CLUSTER_BINDINGS(float, 3, 3, uint16_t, f); - define_ClusterVector(m, "Cluster3x3i"); - define_ClusterVector(m, "Cluster3x3d"); - define_ClusterVector(m, "Cluster3x3f"); - define_ClusterVector(m, "Cluster2x2i"); - define_ClusterVector(m, "Cluster2x2d"); - define_ClusterVector(m, "Cluster2x2f"); + DEFINE_CLUSTER_BINDINGS(int, 2, 2, uint16_t, i); + DEFINE_CLUSTER_BINDINGS(double, 2, 2, uint16_t, d); + DEFINE_CLUSTER_BINDINGS(float, 2, 2, uint16_t, f); - define_cluster_finder_bindings(m, "Cluster3x3i"); - define_cluster_finder_bindings(m, "Cluster3x3d"); - define_cluster_finder_bindings(m, "Cluster3x3f"); - define_cluster_finder_bindings(m, "Cluster2x2i"); - define_cluster_finder_bindings(m, "Cluster2x2d"); - define_cluster_finder_bindings(m, "Cluster2x2f"); + DEFINE_CLUSTER_BINDINGS(int, 5, 5, uint16_t, i); + DEFINE_CLUSTER_BINDINGS(double, 5, 5, uint16_t, d); + DEFINE_CLUSTER_BINDINGS(float, 5, 5, uint16_t, f); - define_cluster_finder_mt_bindings(m, "Cluster3x3i"); - define_cluster_finder_mt_bindings(m, "Cluster3x3d"); - define_cluster_finder_mt_bindings(m, "Cluster3x3f"); - define_cluster_finder_mt_bindings(m, "Cluster2x2i"); - define_cluster_finder_mt_bindings(m, "Cluster2x2d"); - define_cluster_finder_mt_bindings(m, "Cluster2x2f"); + DEFINE_CLUSTER_BINDINGS(int, 7, 7, uint16_t, i); + DEFINE_CLUSTER_BINDINGS(double, 7, 7, uint16_t, d); + DEFINE_CLUSTER_BINDINGS(float, 7, 7, uint16_t, f); - define_cluster_file_sink_bindings(m, "Cluster3x3i"); - define_cluster_file_sink_bindings(m, "Cluster3x3d"); - define_cluster_file_sink_bindings(m, "Cluster3x3f"); - define_cluster_file_sink_bindings(m, "Cluster2x2i"); - define_cluster_file_sink_bindings(m, "Cluster2x2d"); - define_cluster_file_sink_bindings(m, "Cluster2x2f"); - - define_cluster_collector_bindings(m, "Cluster3x3i"); - define_cluster_collector_bindings(m, "Cluster3x3f"); - define_cluster_collector_bindings(m, "Cluster3x3d"); - define_cluster_collector_bindings(m, "Cluster2x2i"); - define_cluster_collector_bindings(m, "Cluster2x2f"); - define_cluster_collector_bindings(m, "Cluster2x2d"); - - define_cluster(m, "3x3i"); - define_cluster(m, "3x3f"); - define_cluster(m, "3x3d"); - define_cluster(m, "2x2i"); - define_cluster(m, "2x2f"); - define_cluster(m, "2x2d"); - - register_calculate_eta(m); - register_calculate_eta(m); - register_calculate_eta(m); - register_calculate_eta(m); - register_calculate_eta(m); - register_calculate_eta(m); + DEFINE_CLUSTER_BINDINGS(int, 9, 9, uint16_t, i); + DEFINE_CLUSTER_BINDINGS(double, 9, 9, uint16_t, d); + DEFINE_CLUSTER_BINDINGS(float, 9, 9, uint16_t, f); } diff --git a/python/src/pedestal.hpp b/python/src/pedestal.hpp index 77148dc..23d8247 100644 --- a/python/src/pedestal.hpp +++ b/python/src/pedestal.hpp @@ -9,7 +9,8 @@ namespace py = pybind11; -template void define_pedestal_bindings(py::module &m, const std::string &name) { +template +void define_pedestal_bindings(py::module &m, const std::string &name) { py::class_>(m, name.c_str()) .def(py::init()) .def(py::init()) @@ -19,16 +20,18 @@ template void define_pedestal_bindings(py::module &m, const *mea = self.mean(); return return_image_data(mea); }) - .def("variance", [](Pedestal &self) { - auto var = new NDArray{}; - *var = self.variance(); - return return_image_data(var); - }) - .def("std", [](Pedestal &self) { - auto std = new NDArray{}; - *std = self.std(); - return return_image_data(std); - }) + .def("variance", + [](Pedestal &self) { + auto var = new NDArray{}; + *var = self.variance(); + return return_image_data(var); + }) + .def("std", + [](Pedestal &self) { + auto std = new NDArray{}; + *std = self.std(); + return return_image_data(std); + }) .def("clear", py::overload_cast<>(&Pedestal::clear)) .def_property_readonly("rows", &Pedestal::rows) .def_property_readonly("cols", &Pedestal::cols) @@ -39,14 +42,19 @@ template void define_pedestal_bindings(py::module &m, const [&](Pedestal &pedestal) { return Pedestal(pedestal); }) - //TODO! add push for other data types - .def("push", [](Pedestal &pedestal, py::array_t &f) { - auto v = make_view_2d(f); - pedestal.push(v); - }) - .def("push_no_update", [](Pedestal &pedestal, py::array_t &f) { - auto v = make_view_2d(f); - pedestal.push_no_update(v); - }, py::arg().noconvert()) + // TODO! add push for other data types + .def("push", + [](Pedestal &pedestal, py::array_t &f) { + auto v = make_view_2d(f); + pedestal.push(v); + }) + .def( + "push_no_update", + [](Pedestal &pedestal, + py::array_t &f) { + auto v = make_view_2d(f); + pedestal.push_no_update(v); + }, + py::arg().noconvert()) .def("update_mean", &Pedestal::update_mean); } \ No newline at end of file diff --git a/python/src/pixel_map.hpp b/python/src/pixel_map.hpp index 46b1bc4..986728b 100644 --- a/python/src/pixel_map.hpp +++ b/python/src/pixel_map.hpp @@ -1,41 +1,46 @@ #include "aare/PixelMap.hpp" #include "np_helper.hpp" - #include #include #include #include - namespace py = pybind11; -using namespace::aare; - +using namespace ::aare; void define_pixel_map_bindings(py::module &m) { - m.def("GenerateMoench03PixelMap", []() { - auto ptr = new NDArray(GenerateMoench03PixelMap()); - return return_image_data(ptr); - }) - .def("GenerateMoench05PixelMap", []() { - auto ptr = new NDArray(GenerateMoench05PixelMap()); - return return_image_data(ptr); - }) - .def("GenerateMoench05PixelMap1g", []() { - auto ptr = new NDArray(GenerateMoench05PixelMap1g()); - return return_image_data(ptr); - }) - .def("GenerateMoench05PixelMapOld", []() { - auto ptr = new NDArray(GenerateMoench05PixelMapOld()); - return return_image_data(ptr); - }) - .def("GenerateMH02SingleCounterPixelMap", []() { - auto ptr = new NDArray(GenerateMH02SingleCounterPixelMap()); - return return_image_data(ptr); - }) - .def("GenerateMH02FourCounterPixelMap", []() { - auto ptr = new NDArray(GenerateMH02FourCounterPixelMap()); - return return_image_data(ptr); - }); - + m.def("GenerateMoench03PixelMap", + []() { + auto ptr = new NDArray(GenerateMoench03PixelMap()); + return return_image_data(ptr); + }) + .def("GenerateMoench05PixelMap", + []() { + auto ptr = new NDArray(GenerateMoench05PixelMap()); + return return_image_data(ptr); + }) + .def("GenerateMoench05PixelMap1g", + []() { + auto ptr = + new NDArray(GenerateMoench05PixelMap1g()); + return return_image_data(ptr); + }) + .def("GenerateMoench05PixelMapOld", + []() { + auto ptr = + new NDArray(GenerateMoench05PixelMapOld()); + return return_image_data(ptr); + }) + .def("GenerateMH02SingleCounterPixelMap", + []() { + auto ptr = new NDArray( + GenerateMH02SingleCounterPixelMap()); + return return_image_data(ptr); + }) + .def("GenerateMH02FourCounterPixelMap", []() { + auto ptr = + new NDArray(GenerateMH02FourCounterPixelMap()); + return return_image_data(ptr); + }); } \ No newline at end of file diff --git a/python/src/raw_file.hpp b/python/src/raw_file.hpp index 38b4896..689b84e 100644 --- a/python/src/raw_file.hpp +++ b/python/src/raw_file.hpp @@ -32,7 +32,7 @@ void define_raw_file_io_bindings(py::module &m) { shape.push_back(self.cols()); // return headers from all subfiles - py::array_t header(self.n_mod()); + py::array_t header(self.n_modules()); const uint8_t item_size = self.bytes_per_pixel(); if (item_size == 1) { @@ -58,13 +58,14 @@ void define_raw_file_io_bindings(py::module &m) { throw std::runtime_error("No frames left in file"); } std::vector shape{n_frames, self.rows(), self.cols()}; - + // return headers from all subfiles py::array_t header; - if (self.n_mod() == 1) { + if (self.n_modules() == 1) { header = py::array_t(n_frames); } else { - header = py::array_t({self.n_mod(), n_frames}); + header = py::array_t( + {self.n_modules(), n_frames}); } // py::array_t header({self.n_mod(), n_frames}); @@ -100,7 +101,7 @@ void define_raw_file_io_bindings(py::module &m) { .def_property_readonly("cols", &RawFile::cols) .def_property_readonly("bitdepth", &RawFile::bitdepth) .def_property_readonly("geometry", &RawFile::geometry) - .def_property_readonly("n_mod", &RawFile::n_mod) + .def_property_readonly("n_modules", &RawFile::n_modules) .def_property_readonly("detector_type", &RawFile::detector_type) .def_property_readonly("master", &RawFile::master); } \ No newline at end of file diff --git a/python/src/raw_master_file.hpp b/python/src/raw_master_file.hpp index 943437f..9c2bd17 100644 --- a/python/src/raw_master_file.hpp +++ b/python/src/raw_master_file.hpp @@ -57,7 +57,8 @@ void define_raw_master_file_bindings(py::module &m) { .def_property_readonly("total_frames_expected", &RawMasterFile::total_frames_expected) .def_property_readonly("geometry", &RawMasterFile::geometry) - .def_property_readonly("analog_samples", &RawMasterFile::analog_samples, R"( + .def_property_readonly("analog_samples", &RawMasterFile::analog_samples, + R"( Number of analog samples Returns @@ -66,7 +67,7 @@ void define_raw_master_file_bindings(py::module &m) { The number of analog samples in the file (or None if not enabled) )") .def_property_readonly("digital_samples", - &RawMasterFile::digital_samples, R"( + &RawMasterFile::digital_samples, R"( Number of digital samples Returns diff --git a/python/src/raw_sub_file.hpp b/python/src/raw_sub_file.hpp index 2cb83fc..cff511b 100644 --- a/python/src/raw_sub_file.hpp +++ b/python/src/raw_sub_file.hpp @@ -24,8 +24,8 @@ auto read_frame_from_RawSubFile(RawSubFile &self) { py::array_t header(1); const uint8_t item_size = self.bytes_per_pixel(); std::vector shape{static_cast(self.rows()), - static_cast(self.cols())}; - + static_cast(self.cols())}; + py::array image; if (item_size == 1) { image = py::array_t(shape); @@ -43,12 +43,10 @@ auto read_frame_from_RawSubFile(RawSubFile &self) { auto read_n_frames_from_RawSubFile(RawSubFile &self, size_t n_frames) { py::array_t header(n_frames); const uint8_t item_size = self.bytes_per_pixel(); - std::vector shape{ - static_cast(n_frames), - static_cast(self.rows()), - static_cast(self.cols()) - }; - + std::vector shape{static_cast(n_frames), + static_cast(self.rows()), + static_cast(self.cols())}; + py::array image; if (item_size == 1) { image = py::array_t(shape); @@ -57,15 +55,14 @@ auto read_n_frames_from_RawSubFile(RawSubFile &self, size_t n_frames) { } else if (item_size == 4) { image = py::array_t(shape); } - self.read_into(reinterpret_cast(image.mutable_data()), n_frames, - header.mutable_data()); + self.read_into(reinterpret_cast(image.mutable_data()), + n_frames, header.mutable_data()); return py::make_tuple(header, image); } - -//Disable warnings for unused parameters, as we ignore some -//in the __exit__ method +// Disable warnings for unused parameters, as we ignore some +// in the __exit__ method #pragma GCC diagnostic push #pragma GCC diagnostic ignored "-Wunused-parameter" @@ -76,7 +73,7 @@ void define_raw_sub_file_io_bindings(py::module &m) { .def_property_readonly("bytes_per_frame", &RawSubFile::bytes_per_frame) .def_property_readonly("pixels_per_frame", &RawSubFile::pixels_per_frame) - .def_property_readonly("bytes_per_pixel", &RawSubFile::bytes_per_pixel) + .def_property_readonly("bytes_per_pixel", &RawSubFile::bytes_per_pixel) .def("seek", &RawSubFile::seek) .def("tell", &RawSubFile::tell) .def_property_readonly("rows", &RawSubFile::rows) @@ -84,18 +81,17 @@ void define_raw_sub_file_io_bindings(py::module &m) { .def_property_readonly("frames_in_file", &RawSubFile::frames_in_file) .def("read_frame", &read_frame_from_RawSubFile) .def("read_n", &read_n_frames_from_RawSubFile) - .def("read", [](RawSubFile &self){ - self.seek(0); - auto n_frames = self.frames_in_file(); - return read_n_frames_from_RawSubFile(self, n_frames); - }) + .def("read", + [](RawSubFile &self) { + self.seek(0); + auto n_frames = self.frames_in_file(); + return read_n_frames_from_RawSubFile(self, n_frames); + }) .def("__enter__", [](RawSubFile &self) { return &self; }) .def("__exit__", - [](RawSubFile &self, - const std::optional &exc_type, + [](RawSubFile &self, const std::optional &exc_type, const std::optional &exc_value, - const std::optional &traceback) { - }) + const std::optional &traceback) {}) .def("__iter__", [](RawSubFile &self) { return &self; }) .def("__next__", [](RawSubFile &self) { try { @@ -104,7 +100,6 @@ void define_raw_sub_file_io_bindings(py::module &m) { throw py::stop_iteration(); } }); - } #pragma GCC diagnostic pop \ No newline at end of file diff --git a/python/src/var_cluster.hpp b/python/src/var_cluster.hpp index f7b373f..4e7302d 100644 --- a/python/src/var_cluster.hpp +++ b/python/src/var_cluster.hpp @@ -12,10 +12,8 @@ // #include // #include - namespace py = pybind11; -using namespace::aare; - +using namespace ::aare; void define_var_cluster_finder_bindings(py::module &m) { PYBIND11_NUMPY_DTYPE(VarClusterFinder::Hit, size, row, col, @@ -29,12 +27,12 @@ void define_var_cluster_finder_bindings(py::module &m) { return return_image_data(ptr); }) .def("set_noiseMap", - [](VarClusterFinder &self, + [](VarClusterFinder &self, py::array_t noise_map) { - auto noise_map_span = make_view_2d(noise_map); - self.set_noiseMap(noise_map_span); - }) + auto noise_map_span = make_view_2d(noise_map); + self.set_noiseMap(noise_map_span); + }) .def("set_peripheralThresholdFactor", &VarClusterFinder::set_peripheralThresholdFactor) .def("find_clusters", @@ -65,9 +63,7 @@ void define_var_cluster_finder_bindings(py::module &m) { return return_vector(ptr); }) .def("clear_hits", - [](VarClusterFinder &self) { - self.clear_hits(); - }) + [](VarClusterFinder &self) { self.clear_hits(); }) .def("steal_hits", [](VarClusterFinder &self) { auto ptr = new std::vector::Hit>( @@ -75,5 +71,4 @@ void define_var_cluster_finder_bindings(py::module &m) { return return_vector(ptr); }) .def("total_clusters", &VarClusterFinder::total_clusters); - } \ No newline at end of file diff --git a/python/tests/test_RawSubFile.py b/python/tests/test_RawSubFile.py index a5eea91..aa4721a 100644 --- a/python/tests/test_RawSubFile.py +++ b/python/tests/test_RawSubFile.py @@ -5,32 +5,35 @@ from aare import RawSubFile, DetectorType @pytest.mark.files def test_read_a_jungfrau_RawSubFile(test_data_path): + + # Starting with f1 there is now 7 frames left in the series of files with RawSubFile(test_data_path / "raw/jungfrau/jungfrau_single_d0_f1_0.raw", DetectorType.Jungfrau, 512, 1024, 16) as f: - assert f.frames_in_file == 3 + assert f.frames_in_file == 7 headers, frames = f.read() - assert headers.size == 3 - assert frames.shape == (3, 512, 1024) + assert headers.size == 7 + assert frames.shape == (7, 512, 1024) - # Frame numbers in this file should be 4, 5, 6 - for i,h in zip(range(4,7,1), headers): + + for i,h in zip(range(4,11,1), headers): assert h["frameNumber"] == i # Compare to canned data using numpy data = np.load(test_data_path / "raw/jungfrau/jungfrau_single_0.npy") - assert np.all(data[3:6] == frames) + assert np.all(data[3:] == frames) @pytest.mark.files def test_iterate_over_a_jungfrau_RawSubFile(test_data_path): data = np.load(test_data_path / "raw/jungfrau/jungfrau_single_0.npy") + # Given the first subfile in a series we can read all frames from f0, f1, f2...fN with RawSubFile(test_data_path / "raw/jungfrau/jungfrau_single_d0_f0_0.raw", DetectorType.Jungfrau, 512, 1024, 16) as f: i = 0 for header, frame in f: assert header["frameNumber"] == i+1 assert np.all(frame == data[i]) i += 1 - assert i == 3 - assert header["frameNumber"] == 3 \ No newline at end of file + assert i == 10 + assert header["frameNumber"] == 10 diff --git a/src/ClusterFile.cpp b/src/ClusterFile.cpp new file mode 100644 index 0000000..13f7364 --- /dev/null +++ b/src/ClusterFile.cpp @@ -0,0 +1,395 @@ +#include "aare/ClusterFile.hpp" + +#include + +namespace aare { + +ClusterFile::ClusterFile(const std::filesystem::path &fname, size_t chunk_size, + const std::string &mode) + : m_chunk_size(chunk_size), m_mode(mode) { + + if (mode == "r") { + fp = fopen(fname.c_str(), "rb"); + if (!fp) { + throw std::runtime_error("Could not open file for reading: " + + fname.string()); + } + } else if (mode == "w") { + fp = fopen(fname.c_str(), "wb"); + if (!fp) { + throw std::runtime_error("Could not open file for writing: " + + fname.string()); + } + } else if (mode == "a") { + fp = fopen(fname.c_str(), "ab"); + if (!fp) { + throw std::runtime_error("Could not open file for appending: " + + fname.string()); + } + } else { + throw std::runtime_error("Unsupported mode: " + mode); + } +} + +void ClusterFile::set_roi(ROI roi) { m_roi = roi; } + +void ClusterFile::set_noise_map(const NDView noise_map) { + m_noise_map = NDArray(noise_map); +} + +void ClusterFile::set_gain_map(const NDView gain_map) { + m_gain_map = NDArray(gain_map); + + // Gain map is passed as ADU/keV to avoid dividing in when applying the gain + // map we invert it here + for (auto &item : m_gain_map->view()) { + item = 1.0 / item; + } +} + +ClusterFile::~ClusterFile() { close(); } + +void ClusterFile::close() { + if (fp) { + fclose(fp); + fp = nullptr; + } +} + +void ClusterFile::write_frame(const ClusterVector &clusters) { + if (m_mode != "w" && m_mode != "a") { + throw std::runtime_error("File not opened for writing"); + } + if (!(clusters.cluster_size_x() == 3) && + !(clusters.cluster_size_y() == 3)) { + throw std::runtime_error("Only 3x3 clusters are supported"); + } + // First write the frame number - 4 bytes + int32_t frame_number = clusters.frame_number(); + if (fwrite(&frame_number, sizeof(frame_number), 1, fp) != 1) { + throw std::runtime_error(LOCATION + "Could not write frame number"); + } + + // Then write the number of clusters - 4 bytes + uint32_t n_clusters = clusters.size(); + if (fwrite(&n_clusters, sizeof(n_clusters), 1, fp) != 1) { + throw std::runtime_error(LOCATION + + "Could not write number of clusters"); + } + + // Now write the clusters in the frame + if (fwrite(clusters.data(), clusters.item_size(), clusters.size(), fp) != + clusters.size()) { + throw std::runtime_error(LOCATION + "Could not write clusters"); + } +} + +ClusterVector ClusterFile::read_clusters(size_t n_clusters) { + if (m_mode != "r") { + throw std::runtime_error("File not opened for reading"); + } + if (m_noise_map || m_roi) { + return read_clusters_with_cut(n_clusters); + } else { + return read_clusters_without_cut(n_clusters); + } +} + +ClusterVector +ClusterFile::read_clusters_without_cut(size_t n_clusters) { + if (m_mode != "r") { + throw std::runtime_error("File not opened for reading"); + } + + ClusterVector clusters(3, 3, n_clusters); + + int32_t iframe = 0; // frame number needs to be 4 bytes! + size_t nph_read = 0; + uint32_t nn = m_num_left; + uint32_t nph = m_num_left; // number of clusters in frame needs to be 4 + + // auto buf = reinterpret_cast(clusters.data()); + auto buf = clusters.data(); + // if there are photons left from previous frame read them first + if (nph) { + if (nph > n_clusters) { + // if we have more photons left in the frame then photons to read we + // read directly the requested number + nn = n_clusters; + } else { + nn = nph; + } + nph_read += fread((buf + nph_read * clusters.item_size()), + clusters.item_size(), nn, fp); + m_num_left = nph - nn; // write back the number of photons left + } + + if (nph_read < n_clusters) { + // keep on reading frames and photons until reaching n_clusters + while (fread(&iframe, sizeof(iframe), 1, fp)) { + 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; + + nph_read += fread((buf + nph_read * clusters.item_size()), + clusters.item_size(), nn, fp); + m_num_left = nph - nn; + } + if (nph_read >= n_clusters) + break; + } + } + + // Resize the vector to the number of clusters. + // No new allocation, only change bounds. + clusters.resize(nph_read); + if (m_gain_map) + clusters.apply_gain_map(m_gain_map->view()); + return clusters; +} + +ClusterVector ClusterFile::read_clusters_with_cut(size_t n_clusters) { + ClusterVector clusters(3, 3); + clusters.reserve(n_clusters); + + // if there are photons left from previous frame read them first + if (m_num_left) { + while (m_num_left && clusters.size() < n_clusters) { + Cluster3x3 c = read_one_cluster(); + if (is_selected(c)) { + clusters.push_back(c.x, c.y, + reinterpret_cast(c.data)); + } + } + } + + // we did not have enough clusters left in the previous frame + // keep on reading frames until reaching n_clusters + if (clusters.size() < n_clusters) { + // sanity check + if (m_num_left) { + throw std::runtime_error( + LOCATION + "Entered second loop with clusters left\n"); + } + + 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) { + Cluster3x3 c = read_one_cluster(); + if (is_selected(c)) { + clusters.push_back( + c.x, c.y, reinterpret_cast(c.data)); + } + } + } + + // we have enough clusters, break out of the outer while loop + if (clusters.size() >= n_clusters) + break; + } + } + if (m_gain_map) + clusters.apply_gain_map(m_gain_map->view()); + + return clusters; +} + +Cluster3x3 ClusterFile::read_one_cluster() { + Cluster3x3 c; + auto rc = fread(&c, sizeof(c), 1, fp); + if (rc != 1) { + throw std::runtime_error(LOCATION + "Could not read cluster"); + } + --m_num_left; + return c; +} + +ClusterVector ClusterFile::read_frame() { + 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(); + } +} + +ClusterVector ClusterFile::read_frame_without_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"); + } + int32_t frame_number; + if (fread(&frame_number, sizeof(frame_number), 1, fp) != 1) { + throw std::runtime_error(LOCATION + "Could not read frame number"); + } + + int32_t n_clusters; // Saved as 32bit integer in the cluster file + if (fread(&n_clusters, sizeof(n_clusters), 1, fp) != 1) { + throw std::runtime_error(LOCATION + + "Could not read number of clusters"); + } + + ClusterVector clusters(3, 3, n_clusters); + clusters.set_frame_number(frame_number); + + if (fread(clusters.data(), clusters.item_size(), n_clusters, fp) != + static_cast(n_clusters)) { + throw std::runtime_error(LOCATION + "Could not read clusters"); + } + clusters.resize(n_clusters); + if (m_gain_map) + clusters.apply_gain_map(m_gain_map->view()); + return clusters; +} + +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"); + } + int32_t frame_number; + if (fread(&frame_number, sizeof(frame_number), 1, fp) != 1) { + throw std::runtime_error("Could not read frame number"); + } + + if (fread(&m_num_left, sizeof(m_num_left), 1, fp) != 1) { + throw std::runtime_error("Could not read number of clusters"); + } + + ClusterVector clusters(3, 3); + clusters.reserve(m_num_left); + clusters.set_frame_number(frame_number); + while (m_num_left) { + Cluster3x3 c = read_one_cluster(); + if (is_selected(c)) { + clusters.push_back(c.x, c.y, reinterpret_cast(c.data)); + } + } + if (m_gain_map) + clusters.apply_gain_map(m_gain_map->view()); + return clusters; +} + +bool ClusterFile::is_selected(Cluster3x3 &cl) { + // Should fail fast + if (m_roi) { + if (!(m_roi->contains(cl.x, cl.y))) { + return false; + } + } + if (m_noise_map) { + int32_t sum_1x1 = cl.data[4]; // central pixel + int32_t sum_2x2 = cl.sum_2x2(); // highest sum of 2x2 subclusters + int32_t sum_3x3 = cl.sum(); // sum of all pixels + + auto noise = + (*m_noise_map)(cl.y, cl.x); // TODO! check if this is correct + if (sum_1x1 <= noise || sum_2x2 <= 2 * noise || sum_3x3 <= 3 * noise) { + return false; + } + } + // we passed all checks + return true; +} + +NDArray calculate_eta2(ClusterVector &clusters) { + // TOTO! make work with 2x2 clusters + NDArray eta2({static_cast(clusters.size()), 2}); + + if (clusters.cluster_size_x() == 3 || clusters.cluster_size_y() == 3) { + for (size_t i = 0; i < clusters.size(); i++) { + auto e = calculate_eta2(clusters.at(i)); + eta2(i, 0) = e.x; + eta2(i, 1) = e.y; + } + } else if (clusters.cluster_size_x() == 2 || + clusters.cluster_size_y() == 2) { + for (size_t i = 0; i < clusters.size(); i++) { + auto e = calculate_eta2(clusters.at(i)); + eta2(i, 0) = e.x; + eta2(i, 1) = e.y; + } + } else { + throw std::runtime_error("Only 3x3 and 2x2 clusters are supported"); + } + + return eta2; +} + +/** + * @brief Calculate the eta2 values for a 3x3 cluster and return them in a Eta2 + * struct containing etay, etax and the corner of the cluster. + */ +Eta2 calculate_eta2(Cluster3x3 &cl) { + Eta2 eta{}; + + std::array tot2; + tot2[0] = cl.data[0] + cl.data[1] + cl.data[3] + cl.data[4]; + tot2[1] = cl.data[1] + cl.data[2] + cl.data[4] + cl.data[5]; + tot2[2] = cl.data[3] + cl.data[4] + cl.data[6] + cl.data[7]; + tot2[3] = cl.data[4] + cl.data[5] + cl.data[7] + cl.data[8]; + + auto c = std::max_element(tot2.begin(), tot2.end()) - tot2.begin(); + eta.sum = tot2[c]; + switch (c) { + case cBottomLeft: + if ((cl.data[3] + cl.data[4]) != 0) + eta.x = static_cast(cl.data[4]) / (cl.data[3] + cl.data[4]); + if ((cl.data[1] + cl.data[4]) != 0) + eta.y = static_cast(cl.data[4]) / (cl.data[1] + cl.data[4]); + eta.c = cBottomLeft; + break; + case cBottomRight: + if ((cl.data[2] + cl.data[5]) != 0) + eta.x = static_cast(cl.data[5]) / (cl.data[4] + cl.data[5]); + if ((cl.data[1] + cl.data[4]) != 0) + eta.y = static_cast(cl.data[4]) / (cl.data[1] + cl.data[4]); + eta.c = cBottomRight; + break; + case cTopLeft: + if ((cl.data[7] + cl.data[4]) != 0) + eta.x = static_cast(cl.data[4]) / (cl.data[3] + cl.data[4]); + if ((cl.data[7] + cl.data[4]) != 0) + eta.y = static_cast(cl.data[7]) / (cl.data[7] + cl.data[4]); + eta.c = cTopLeft; + break; + case cTopRight: + if ((cl.data[5] + cl.data[4]) != 0) + eta.x = static_cast(cl.data[5]) / (cl.data[5] + cl.data[4]); + if ((cl.data[7] + cl.data[4]) != 0) + eta.y = static_cast(cl.data[7]) / (cl.data[7] + cl.data[4]); + eta.c = cTopRight; + break; + // no default to allow compiler to warn about missing cases + } + return eta; +} + +Eta2 calculate_eta2(Cluster2x2 &cl) { + Eta2 eta{}; + if ((cl.data[0] + cl.data[1]) != 0) + eta.x = static_cast(cl.data[1]) / (cl.data[0] + cl.data[1]); + if ((cl.data[0] + cl.data[2]) != 0) + eta.y = static_cast(cl.data[2]) / (cl.data[0] + cl.data[2]); + eta.sum = cl.data[0] + cl.data[1] + cl.data[2] + cl.data[3]; + eta.c = cBottomLeft; // TODO! This is not correct, but need to put something + return eta; +} + +} // namespace aare \ No newline at end of file diff --git a/src/ClusterFile.test.cpp b/src/ClusterFile.test.cpp index 6254b5d..68e734d 100644 --- a/src/ClusterFile.test.cpp +++ b/src/ClusterFile.test.cpp @@ -10,9 +10,8 @@ using aare::Cluster; using aare::ClusterFile; using aare::ClusterVector; - TEST_CASE("Read one frame from a cluster file", "[.files]") { - //We know that the frame has 97 clusters + // We know that the frame has 97 clusters auto fpath = test_data_path() / "clust" / "single_frame_97_clustrers.clust"; REQUIRE(std::filesystem::exists(fpath)); @@ -27,7 +26,6 @@ TEST_CASE("Read one frame from a cluster file", "[.files]") { std::begin(expected_cluster_data))); } - TEST_CASE("Read one frame using ROI", "[.files]") { // We know that the frame has 97 clusters auto fpath = test_data_path() / "clust" / "single_frame_97_clustrers.clust"; @@ -60,8 +58,6 @@ TEST_CASE("Read one frame using ROI", "[.files]") { std::begin(expected_cluster_data))); } - - TEST_CASE("Read clusters from single frame file", "[.files]") { // frame_number, num_clusters [135] 97 diff --git a/src/CtbRawFile.cpp b/src/CtbRawFile.cpp index 4d9d895..a6a1d92 100644 --- a/src/CtbRawFile.cpp +++ b/src/CtbRawFile.cpp @@ -14,22 +14,24 @@ CtbRawFile::CtbRawFile(const std::filesystem::path &fname) : m_master(fname) { m_file.open(m_master.data_fname(0, 0), std::ios::binary); } -void CtbRawFile::read_into(std::byte *image_buf, DetectorHeader* header) { - if(m_current_frame >= m_master.frames_in_file()){ +void CtbRawFile::read_into(std::byte *image_buf, DetectorHeader *header) { + if (m_current_frame >= m_master.frames_in_file()) { throw std::runtime_error(LOCATION + " End of file reached"); } - if(m_current_frame != 0 && m_current_frame % m_master.max_frames_per_file() == 0){ - open_data_file(m_current_subfile+1); + if (m_current_frame != 0 && + m_current_frame % m_master.max_frames_per_file() == 0) { + open_data_file(m_current_subfile + 1); } - - if(header){ + + if (header) { m_file.read(reinterpret_cast(header), sizeof(DetectorHeader)); - }else{ + } else { m_file.seekg(sizeof(DetectorHeader), std::ios::cur); } - m_file.read(reinterpret_cast(image_buf), m_master.image_size_in_bytes()); + m_file.read(reinterpret_cast(image_buf), + m_master.image_size_in_bytes()); m_current_frame++; } @@ -38,13 +40,16 @@ void CtbRawFile::seek(size_t frame_number) { open_data_file(index); } size_t frame_number_in_file = frame_number % m_master.max_frames_per_file(); - m_file.seekg((sizeof(DetectorHeader)+m_master.image_size_in_bytes()) * frame_number_in_file); + m_file.seekg((sizeof(DetectorHeader) + m_master.image_size_in_bytes()) * + frame_number_in_file); m_current_frame = frame_number; } size_t CtbRawFile::tell() const { return m_current_frame; } -size_t CtbRawFile::image_size_in_bytes() const { return m_master.image_size_in_bytes(); } +size_t CtbRawFile::image_size_in_bytes() const { + return m_master.image_size_in_bytes(); +} size_t CtbRawFile::frames_in_file() const { return m_master.frames_in_file(); } @@ -63,12 +68,11 @@ void CtbRawFile::open_data_file(size_t subfile_index) { throw std::runtime_error(LOCATION + "Subfile index out of range"); } m_current_subfile = subfile_index; - m_file = std::ifstream(m_master.data_fname(0, subfile_index), std::ios::binary); // only one module for CTB + m_file = std::ifstream(m_master.data_fname(0, subfile_index), + std::ios::binary); // only one module for CTB if (!m_file.is_open()) { throw std::runtime_error(LOCATION + "Could not open data file"); } } - - } // namespace aare \ No newline at end of file diff --git a/src/Dtype.cpp b/src/Dtype.cpp index b818ea3..0fdffec 100644 --- a/src/Dtype.cpp +++ b/src/Dtype.cpp @@ -10,7 +10,8 @@ namespace aare { * @brief Construct a DType object from a type_info object * @param t type_info object * @throw runtime_error if the type is not supported - * @note supported types are: int8_t, uint8_t, int16_t, uint16_t, int32_t, uint32_t, int64_t, uint64_t, float, double + * @note supported types are: int8_t, uint8_t, int16_t, uint16_t, int32_t, + * uint32_t, int64_t, uint64_t, float, double * @note the type_info object is obtained using typeid (e.g. typeid(int)) */ Dtype::Dtype(const std::type_info &t) { @@ -35,7 +36,8 @@ Dtype::Dtype(const std::type_info &t) { else if (t == typeid(double)) m_type = TypeIndex::DOUBLE; else - throw std::runtime_error("Could not construct data type. Type not supported."); + throw std::runtime_error( + "Could not construct data type. Type not supported."); } /** @@ -63,7 +65,8 @@ uint8_t Dtype::bitdepth() const { case TypeIndex::NONE: return 0; default: - throw std::runtime_error(LOCATION + "Could not get bitdepth. Type not supported."); + throw std::runtime_error(LOCATION + + "Could not get bitdepth. Type not supported."); } } @@ -138,7 +141,8 @@ Dtype Dtype::from_bitdepth(uint8_t bitdepth) { case 64: return Dtype(TypeIndex::UINT64); default: - throw std::runtime_error("Could not construct data type from bitdepth."); + throw std::runtime_error( + "Could not construct data type from bitdepth."); } } /** @@ -175,17 +179,27 @@ std::string Dtype::to_string() const { case TypeIndex::DOUBLE: return "f8"; case TypeIndex::ERROR: - throw std::runtime_error("Could not get string representation. Type not supported."); + throw std::runtime_error( + "Could not get string representation. Type not supported."); case TypeIndex::NONE: - throw std::runtime_error("Could not get string representation. Type not supported."); + throw std::runtime_error( + "Could not get string representation. Type not supported."); } return {}; } -bool Dtype::operator==(const Dtype &other) const noexcept { return m_type == other.m_type; } -bool Dtype::operator!=(const Dtype &other) const noexcept { return !(*this == other); } +bool Dtype::operator==(const Dtype &other) const noexcept { + return m_type == other.m_type; +} +bool Dtype::operator!=(const Dtype &other) const noexcept { + return !(*this == other); +} -bool Dtype::operator==(const std::type_info &t) const { return Dtype(t) == *this; } -bool Dtype::operator!=(const std::type_info &t) const { return Dtype(t) != *this; } +bool Dtype::operator==(const std::type_info &t) const { + return Dtype(t) == *this; +} +bool Dtype::operator!=(const std::type_info &t) const { + return Dtype(t) != *this; +} } // namespace aare diff --git a/src/Dtype.test.cpp b/src/Dtype.test.cpp index b252267..256be64 100644 --- a/src/Dtype.test.cpp +++ b/src/Dtype.test.cpp @@ -51,4 +51,6 @@ TEST_CASE("Construct from string with endianess") { REQUIRE_THROWS(Dtype(">i4") == typeid(int32_t)); } -TEST_CASE("Convert to string") { REQUIRE(Dtype(typeid(int)).to_string() == "(fname, mode); - } - else if (fname.extension() == ".npy") { + } else if (fname.extension() == ".npy") { // file_impl = new NumpyFile(fname, mode, cfg); file_impl = std::make_unique(fname, mode, cfg); - }else if(fname.extension() == ".dat"){ + } else if (fname.extension() == ".dat") { file_impl = std::make_unique(fname); } else { throw std::runtime_error("Unsupported file type"); } } +File::File(File &&other) noexcept { std::swap(file_impl, other.file_impl); } -File::File(File &&other) noexcept{ - std::swap(file_impl, other.file_impl); -} - -File& File::operator=(File &&other) noexcept { +File &File::operator=(File &&other) noexcept { if (this != &other) { File tmp(std::move(other)); std::swap(file_impl, tmp.file_impl); @@ -70,15 +66,16 @@ size_t File::frame_number(size_t frame_index) { } size_t File::bytes_per_frame() const { return file_impl->bytes_per_frame(); } -size_t File::pixels_per_frame() const{ return file_impl->pixels_per_frame(); } +size_t File::pixels_per_frame() const { return file_impl->pixels_per_frame(); } void File::seek(size_t frame_index) { file_impl->seek(frame_index); } 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(); } -size_t File::bytes_per_pixel() const { return file_impl->bitdepth() / bits_per_byte; } +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 \ No newline at end of file diff --git a/src/FilePtr.cpp b/src/FilePtr.cpp index e3cdb4b..f850080 100644 --- a/src/FilePtr.cpp +++ b/src/FilePtr.cpp @@ -6,10 +6,12 @@ namespace aare { -FilePtr::FilePtr(const std::filesystem::path& fname, const std::string& mode = "rb") { +FilePtr::FilePtr(const std::filesystem::path &fname, + const std::string &mode = "rb") { fp_ = fopen(fname.c_str(), mode.c_str()); if (!fp_) - throw std::runtime_error(fmt::format("Could not open: {}", fname.c_str())); + throw std::runtime_error( + fmt::format("Could not open: {}", fname.c_str())); } FilePtr::FilePtr(FilePtr &&other) { std::swap(fp_, other.fp_); } @@ -24,15 +26,16 @@ FILE *FilePtr::get() { return fp_; } ssize_t FilePtr::tell() { auto pos = ftell(fp_); if (pos == -1) - throw std::runtime_error(fmt::format("Error getting file position: {}", error_msg())); + throw std::runtime_error( + fmt::format("Error getting file position: {}", error_msg())); return pos; -} +} FilePtr::~FilePtr() { if (fp_) fclose(fp_); // check? } -std::string FilePtr::error_msg(){ +std::string FilePtr::error_msg() { if (feof(fp_)) { return "End of file reached"; } diff --git a/src/Fit.cpp b/src/Fit.cpp index d104675..9d4b70b 100644 --- a/src/Fit.cpp +++ b/src/Fit.cpp @@ -1,13 +1,12 @@ #include "aare/Fit.hpp" -#include "aare/utils/task.hpp" #include "aare/utils/par.hpp" +#include "aare/utils/task.hpp" #include #include #include #include - namespace aare { namespace func { @@ -34,8 +33,10 @@ NDArray pol1(NDView x, NDView par) { 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])); +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 scurve(NDView x, NDView par) { @@ -46,8 +47,10 @@ NDArray scurve(NDView x, NDView par) { 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])); +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 scurve2(NDView x, NDView par) { @@ -91,7 +94,8 @@ NDArray fit_gaus(NDView x, NDView y, return result; } -std::array gaus_init_par(const NDView x, const NDView y) { +std::array gaus_init_par(const NDView x, + const NDView y) { std::array start_par{0, 0, 0}; auto e = std::max_element(y.begin(), y.end()); auto idx = std::distance(y.begin(), e); @@ -103,31 +107,29 @@ std::array gaus_init_par(const NDView x, const NDView *e / 2; }) * - delta / 2.35; + start_par[2] = std::count_if(y.begin(), y.end(), + [e](double val) { return val > *e / 2; }) * + delta / 2.35; return start_par; } +std::array pol1_init_par(const NDView x, + const NDView y) { + // Estimate the initial parameters for the fit + std::array start_par{0, 0}; -std::array pol1_init_par(const NDView x, const NDView y){ - // Estimate the initial parameters for the fit - std::array start_par{0, 0}; + auto y2 = std::max_element(y.begin(), y.end()); + auto x2 = x[std::distance(y.begin(), y2)]; + auto y1 = std::min_element(y.begin(), y.end()); + auto x1 = x[std::distance(y.begin(), y1)]; - - auto y2 = std::max_element(y.begin(), y.end()); - auto x2 = x[std::distance(y.begin(), y2)]; - auto y1 = std::min_element(y.begin(), y.end()); - auto x1 = x[std::distance(y.begin(), y1)]; - - start_par[0] = - (*y2 - *y1) / (x2 - x1); // For amplitude we use the maximum value - start_par[1] = - *y1 - ((*y2 - *y1) / (x2 - x1)) * - x1; // For the mean we use the x value of the maximum value - return start_par; + start_par[0] = + (*y2 - *y1) / (x2 - x1); // For amplitude we use the maximum value + start_par[1] = + *y1 - ((*y2 - *y1) / (x2 - x1)) * + x1; // For the mean we use the x value of the maximum value + return start_par; } void fit_gaus(NDView x, NDView y, NDView y_err, @@ -141,7 +143,6 @@ void fit_gaus(NDView x, NDView y, NDView y_err, "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. */ @@ -153,23 +154,32 @@ void fit_gaus(NDView x, NDView y, NDView y_err, // */ // } lm_status_struct; - lm_status_struct status; par_out = gaus_init_par(x, y); - std::array cov{0, 0, 0, 0, 0, 0, 0 , 0 , 0}; + std::array 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). + // 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, @@ -178,12 +188,14 @@ void fit_gaus(NDView x, NDView y, NDView y_err, // 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); + chi2 += + std::pow((y(i) - func::gaus(x(i), par_out.data())) / y_err(i), 2); } } void fit_gaus(NDView x, NDView y, NDView y_err, - NDView par_out, NDView par_err_out, NDView chi2_out, + NDView par_out, NDView par_err_out, + NDView chi2_out, int n_threads) { @@ -197,10 +209,9 @@ void fit_gaus(NDView x, NDView y, NDView y_err, {par_out.shape(2)}); NDView 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)); - } } }; @@ -210,7 +221,8 @@ void fit_gaus(NDView x, NDView y, NDView y_err, } void fit_pol1(NDView x, NDView y, NDView y_err, - NDView par_out, NDView par_err_out, double& chi2) { + NDView par_out, NDView par_err_out, + double &chi2) { // Check that we have the correct sizes if (y.size() != x.size() || y.size() != y_err.size() || @@ -230,13 +242,14 @@ void fit_pol1(NDView x, NDView y, NDView y_err, // 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); + chi2 += + std::pow((y(i) - func::pol1(x(i), par_out.data())) / y_err(i), 2); } } void fit_pol1(NDView x, NDView y, NDView y_err, - NDView par_out, NDView par_err_out, NDView chi2_out, - int n_threads) { + NDView par_out, NDView par_err_out, + NDView 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++) { @@ -249,15 +262,14 @@ void fit_pol1(NDView x, NDView y, NDView y_err, NDView 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)); - + 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 fit_pol1(NDView x, NDView y) { @@ -300,27 +312,29 @@ NDArray fit_pol1(NDView x, NDView y, // ~~ S-CURVES ~~ // SCURVE -- -std::array scurve_init_par(const NDView x, const NDView y){ - // Estimate the initial parameters for the fit - std::array start_par{0, 0, 0, 0, 0, 0}; +std::array scurve_init_par(const NDView x, + const NDView y) { + // Estimate the initial parameters for the fit + std::array start_par{0, 0, 0, 0, 0, 0}; - auto ymax = std::max_element(y.begin(), y.end()); - auto ymin = std::min_element(y.begin(), y.end()); - start_par[4] = *ymin + (*ymax - *ymin) / 2; - - // Find the first x where the corresponding y value is above the threshold (start_par[4]) - for (ssize_t i = 0; i < y.size(); ++i) { - if (y[i] >= start_par[4]) { - start_par[2] = x[i]; - break; // Exit the loop after finding the first valid x - } + auto ymax = std::max_element(y.begin(), y.end()); + auto ymin = std::min_element(y.begin(), y.end()); + start_par[4] = *ymin + (*ymax - *ymin) / 2; + + // Find the first x where the corresponding y value is above the threshold + // (start_par[4]) + for (ssize_t i = 0; i < y.size(); ++i) { + if (y[i] >= start_par[4]) { + start_par[2] = x[i]; + break; // Exit the loop after finding the first valid x } + } - start_par[3] = 2 * sqrt(start_par[2]); - start_par[0] = 100; - start_par[1] = 0.25; - start_par[5] = 1; - return start_par; + start_par[3] = 2 * sqrt(start_par[2]); + start_par[0] = 100; + start_par[1] = 0.25; + start_par[5] = 1; + return start_par; } // - No error @@ -334,7 +348,8 @@ NDArray fit_scurve(NDView x, NDView y) { return result; } -NDArray fit_scurve(NDView x, NDView y, int n_threads) { +NDArray fit_scurve(NDView x, NDView y, + int n_threads) { NDArray result({y.shape(0), y.shape(1), 6}, 0); auto process = [&x, &y, &result](ssize_t first_row, ssize_t last_row) { @@ -358,8 +373,9 @@ NDArray fit_scurve(NDView x, NDView y, int n_th } // - Error -void fit_scurve(NDView x, NDView y, NDView y_err, - NDView par_out, NDView par_err_out, double& chi2) { +void fit_scurve(NDView x, NDView y, + NDView y_err, NDView par_out, + NDView par_err_out, double &chi2) { // Check that we have the correct sizes if (y.size() != x.size() || y.size() != y_err.size() || @@ -380,13 +396,15 @@ void fit_scurve(NDView x, NDView y, NDView y_er // 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); + chi2 += + std::pow((y(i) - func::pol1(x(i), par_out.data())) / y_err(i), 2); } } -void fit_scurve(NDView x, NDView y, NDView y_err, - NDView par_out, NDView par_err_out, NDView chi2_out, - int n_threads) { +void fit_scurve(NDView x, NDView y, + NDView y_err, NDView par_out, + NDView par_err_out, NDView 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++) { @@ -399,40 +417,41 @@ void fit_scurve(NDView x, NDView y, NDView y_er NDView 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)); - + 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 --- -std::array scurve2_init_par(const NDView x, const NDView y){ - // Estimate the initial parameters for the fit - std::array start_par{0, 0, 0, 0, 0, 0}; +std::array scurve2_init_par(const NDView x, + const NDView y) { + // Estimate the initial parameters for the fit + std::array start_par{0, 0, 0, 0, 0, 0}; - auto ymax = std::max_element(y.begin(), y.end()); - auto ymin = std::min_element(y.begin(), y.end()); - start_par[4] = *ymin + (*ymax - *ymin) / 2; - - // Find the first x where the corresponding y value is above the threshold (start_par[4]) - for (ssize_t i = 0; i < y.size(); ++i) { - if (y[i] <= start_par[4]) { - start_par[2] = x[i]; - break; // Exit the loop after finding the first valid x - } + auto ymax = std::max_element(y.begin(), y.end()); + auto ymin = std::min_element(y.begin(), y.end()); + start_par[4] = *ymin + (*ymax - *ymin) / 2; + + // Find the first x where the corresponding y value is above the threshold + // (start_par[4]) + for (ssize_t i = 0; i < y.size(); ++i) { + if (y[i] <= start_par[4]) { + start_par[2] = x[i]; + break; // Exit the loop after finding the first valid x } + } - start_par[3] = 2 * sqrt(start_par[2]); - start_par[0] = 100; - start_par[1] = 0.25; - start_par[5] = -1; - return start_par; + start_par[3] = 2 * sqrt(start_par[2]); + start_par[0] = 100; + start_par[1] = 0.25; + start_par[5] = -1; + return start_par; } // - No error @@ -446,7 +465,8 @@ NDArray fit_scurve2(NDView x, NDView y) { return result; } -NDArray fit_scurve2(NDView x, NDView y, int n_threads) { +NDArray fit_scurve2(NDView x, NDView y, + int n_threads) { NDArray result({y.shape(0), y.shape(1), 6}, 0); auto process = [&x, &y, &result](ssize_t first_row, ssize_t last_row) { @@ -470,8 +490,9 @@ NDArray fit_scurve2(NDView x, NDView y, int n_t } // - Error -void fit_scurve2(NDView x, NDView y, NDView y_err, - NDView par_out, NDView par_err_out, double& chi2) { +void fit_scurve2(NDView x, NDView y, + NDView y_err, NDView par_out, + NDView par_err_out, double &chi2) { // Check that we have the correct sizes if (y.size() != x.size() || y.size() != y_err.size() || @@ -492,13 +513,15 @@ void fit_scurve2(NDView x, NDView y, NDView y_e // 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); + chi2 += + std::pow((y(i) - func::pol1(x(i), par_out.data())) / y_err(i), 2); } } -void fit_scurve2(NDView x, NDView y, NDView y_err, - NDView par_out, NDView par_err_out, NDView chi2_out, - int n_threads) { +void fit_scurve2(NDView x, NDView y, + NDView y_err, NDView par_out, + NDView par_err_out, NDView 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++) { @@ -511,15 +534,14 @@ void fit_scurve2(NDView x, NDView y, NDView y_e NDView 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)); - + 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); - } } // namespace aare \ No newline at end of file diff --git a/src/Frame.cpp b/src/Frame.cpp index d44bed5..ef7675f 100644 --- a/src/Frame.cpp +++ b/src/Frame.cpp @@ -29,8 +29,7 @@ uint64_t Frame::size() const { return m_rows * m_cols; } size_t Frame::bytes() const { return m_rows * m_cols * m_dtype.bytes(); } std::byte *Frame::data() const { return m_data; } - -std::byte *Frame::pixel_ptr(uint32_t row, uint32_t col) const{ +std::byte *Frame::pixel_ptr(uint32_t row, uint32_t col) const { if ((row >= m_rows) || (col >= m_cols)) { std::cerr << "Invalid row or column index" << '\n'; return nullptr; @@ -38,7 +37,6 @@ std::byte *Frame::pixel_ptr(uint32_t row, uint32_t col) const{ return m_data + (row * m_cols + col) * (m_dtype.bytes()); } - Frame &Frame::operator=(Frame &&other) noexcept { if (this == &other) { return *this; @@ -70,5 +68,4 @@ Frame Frame::clone() const { return frame; } - } // namespace aare diff --git a/src/Frame.test.cpp b/src/Frame.test.cpp index 4063701..bafb39c 100644 --- a/src/Frame.test.cpp +++ b/src/Frame.test.cpp @@ -65,7 +65,8 @@ TEST_CASE("Set a value in a 64 bit frame") { // only the value we did set should be non-zero for (size_t i = 0; i < rows; i++) { for (size_t j = 0; j < cols; j++) { - uint64_t *data = reinterpret_cast(frame.pixel_ptr(i, j)); + uint64_t *data = + reinterpret_cast(frame.pixel_ptr(i, j)); REQUIRE(data != nullptr); if (i == 5 && j == 7) { REQUIRE(*data == value); @@ -150,4 +151,3 @@ TEST_CASE("test explicit copy constructor") { REQUIRE(frame2.bytes() == rows * cols * bitdepth / 8); REQUIRE(frame2.data() != data); } - diff --git a/src/JungfrauDataFile.cpp b/src/JungfrauDataFile.cpp index 59a1a0a..5fb99a6 100644 --- a/src/JungfrauDataFile.cpp +++ b/src/JungfrauDataFile.cpp @@ -19,16 +19,15 @@ JungfrauDataFile::JungfrauDataFile(const std::filesystem::path &fname) { open_file(m_current_file_index); } - // FileInterface -Frame JungfrauDataFile::read_frame(){ +Frame JungfrauDataFile::read_frame() { Frame f(rows(), cols(), Dtype::UINT16); read_into(reinterpret_cast(f.data()), nullptr); return f; } -Frame JungfrauDataFile::read_frame(size_t frame_number){ +Frame JungfrauDataFile::read_frame(size_t frame_number) { seek(frame_number); Frame f(rows(), cols(), Dtype::UINT16); read_into(reinterpret_cast(f.data()), nullptr); @@ -37,7 +36,7 @@ Frame JungfrauDataFile::read_frame(size_t frame_number){ std::vector JungfrauDataFile::read_n(size_t n_frames) { std::vector frames; - for(size_t i = 0; i < n_frames; ++i){ + for (size_t i = 0; i < n_frames; ++i) { frames.push_back(read_frame()); } return frames; @@ -48,7 +47,7 @@ void JungfrauDataFile::read_into(std::byte *image_buf) { } void JungfrauDataFile::read_into(std::byte *image_buf, size_t n_frames) { read_into(image_buf, n_frames, nullptr); -} +} size_t JungfrauDataFile::frame_number(size_t frame_index) { seek(frame_index); @@ -59,7 +58,9 @@ std::array JungfrauDataFile::shape() const { return {static_cast(rows()), static_cast(cols())}; } -DetectorType JungfrauDataFile::detector_type() const { return DetectorType::Jungfrau; } +DetectorType JungfrauDataFile::detector_type() const { + return DetectorType::Jungfrau; +} std::string JungfrauDataFile::base_name() const { return m_base_name; } @@ -195,22 +196,23 @@ void JungfrauDataFile::read_into(std::byte *image_buf, size_t n_frames, JungfrauDataHeader *header) { if (header) { for (size_t i = 0; i < n_frames; ++i) - read_into(image_buf + i * m_bytes_per_frame, header + i); - }else{ + read_into(image_buf + i * m_bytes_per_frame, header + i); + } else { for (size_t i = 0; i < n_frames; ++i) read_into(image_buf + i * m_bytes_per_frame, nullptr); } } -void JungfrauDataFile::read_into(NDArray* image, JungfrauDataHeader* header) { - if(image->shape()!=shape()){ - throw std::runtime_error(LOCATION + - "Image shape does not match file size: " + std::to_string(rows()) + "x" + std::to_string(cols())); +void JungfrauDataFile::read_into(NDArray *image, + JungfrauDataHeader *header) { + if (image->shape() != shape()) { + throw std::runtime_error( + LOCATION + "Image shape does not match file size: " + + std::to_string(rows()) + "x" + std::to_string(cols())); } read_into(reinterpret_cast(image->data()), header); } - JungfrauDataHeader JungfrauDataFile::read_header() { JungfrauDataHeader header; if (auto rc = fread(&header, 1, sizeof(header), m_fp.get()); diff --git a/src/JungfrauDataFile.test.cpp b/src/JungfrauDataFile.test.cpp index ce51168..21a4e32 100644 --- a/src/JungfrauDataFile.test.cpp +++ b/src/JungfrauDataFile.test.cpp @@ -1,21 +1,21 @@ #include "aare/JungfrauDataFile.hpp" -#include #include "test_config.hpp" +#include using aare::JungfrauDataFile; using aare::JungfrauDataHeader; TEST_CASE("Open a Jungfrau data file", "[.files]") { - //we know we have 4 files with 7, 7, 7, and 3 frames - //firs frame number if 1 and the bunch id is frame_number**2 - //so we can check the header + // we know we have 4 files with 7, 7, 7, and 3 frames + // firs frame number if 1 and the bunch id is frame_number**2 + // so we can check the header auto fpath = test_data_path() / "dat" / "AldoJF500k_000000.dat"; REQUIRE(std::filesystem::exists(fpath)); JungfrauDataFile f(fpath); REQUIRE(f.rows() == 512); REQUIRE(f.cols() == 1024); - REQUIRE(f.bytes_per_frame() == 1048576); + REQUIRE(f.bytes_per_frame() == 1048576); REQUIRE(f.pixels_per_frame() == 524288); REQUIRE(f.bytes_per_pixel() == 2); REQUIRE(f.bitdepth() == 16); @@ -25,7 +25,7 @@ TEST_CASE("Open a Jungfrau data file", "[.files]") { REQUIRE(f.total_frames() == 24); REQUIRE(f.current_file() == fpath); - //Check that the frame number and buch id is read correctly + // Check that the frame number and buch id is read correctly for (size_t i = 0; i < 24; ++i) { JungfrauDataHeader header; aare::NDArray image(f.shape()); @@ -37,65 +37,64 @@ TEST_CASE("Open a Jungfrau data file", "[.files]") { } } -TEST_CASE("Seek in a JungfrauDataFile", "[.files]"){ +TEST_CASE("Seek in a JungfrauDataFile", "[.files]") { auto fpath = test_data_path() / "dat" / "AldoJF65k_000000.dat"; REQUIRE(std::filesystem::exists(fpath)); JungfrauDataFile f(fpath); - //The file should have 113 frames + // The file should have 113 frames f.seek(19); REQUIRE(f.tell() == 19); auto h = f.read_header(); - REQUIRE(h.framenum == 19+1); + REQUIRE(h.framenum == 19 + 1); - //Reading again does not change the file pointer + // Reading again does not change the file pointer auto h2 = f.read_header(); - REQUIRE(h2.framenum == 19+1); + REQUIRE(h2.framenum == 19 + 1); f.seek(59); REQUIRE(f.tell() == 59); auto h3 = f.read_header(); - REQUIRE(h3.framenum == 59+1); + REQUIRE(h3.framenum == 59 + 1); JungfrauDataHeader h4; aare::NDArray image(f.shape()); f.read_into(&image, &h4); - REQUIRE(h4.framenum == 59+1); + REQUIRE(h4.framenum == 59 + 1); - //now we should be on the next frame + // now we should be on the next frame REQUIRE(f.tell() == 60); - REQUIRE(f.read_header().framenum == 60+1); + REQUIRE(f.read_header().framenum == 60 + 1); - REQUIRE_THROWS(f.seek(86356)); //out of range + REQUIRE_THROWS(f.seek(86356)); // out of range } -TEST_CASE("Open a Jungfrau data file with non zero file index", "[.files]"){ +TEST_CASE("Open a Jungfrau data file with non zero file index", "[.files]") { auto fpath = test_data_path() / "dat" / "AldoJF65k_000003.dat"; REQUIRE(std::filesystem::exists(fpath)); JungfrauDataFile f(fpath); - //18 files per data file, opening the 3rd file we ignore the first 3 - REQUIRE(f.total_frames() == 113-18*3); + // 18 files per data file, opening the 3rd file we ignore the first 3 + REQUIRE(f.total_frames() == 113 - 18 * 3); REQUIRE(f.tell() == 0); - //Frame numbers start at 1 in the first file - REQUIRE(f.read_header().framenum == 18*3+1); + // Frame numbers start at 1 in the first file + REQUIRE(f.read_header().framenum == 18 * 3 + 1); // moving relative to the third file f.seek(5); - REQUIRE(f.read_header().framenum == 18*3+1+5); + REQUIRE(f.read_header().framenum == 18 * 3 + 1 + 5); // ignoring the first 3 files REQUIRE(f.n_files() == 4); REQUIRE(f.current_file().stem() == "AldoJF65k_000003"); - } -TEST_CASE("Read into throws if size doesn't match", "[.files]"){ +TEST_CASE("Read into throws if size doesn't match", "[.files]") { auto fpath = test_data_path() / "dat" / "AldoJF65k_000000.dat"; REQUIRE(std::filesystem::exists(fpath)); @@ -109,6 +108,4 @@ TEST_CASE("Read into throws if size doesn't match", "[.files]"){ REQUIRE_THROWS(f.read_into(&image)); REQUIRE(f.tell() == 0); - - } \ No newline at end of file diff --git a/src/NDArray.test.cpp b/src/NDArray.test.cpp index 819a1a9..91b5933 100644 --- a/src/NDArray.test.cpp +++ b/src/NDArray.test.cpp @@ -35,7 +35,7 @@ TEST_CASE("Construct from an NDView") { } } -TEST_CASE("3D NDArray from NDView"){ +TEST_CASE("3D NDArray from NDView") { std::vector data(27); std::iota(data.begin(), data.end(), 0); NDView view(data.data(), Shape<3>{3, 3, 3}); @@ -44,9 +44,9 @@ TEST_CASE("3D NDArray from NDView"){ REQUIRE(image.size() == view.size()); REQUIRE(image.data() != view.data()); - for(ssize_t i=0; i MultiplyNDArrayUsingOperator(NDArray &a, NDArray &b) { // return a * a * b * b; - NDArrayc = a*b; + NDArray c = a * b; return c; } @@ -162,7 +162,6 @@ NDArray AddNDArrayUsingIndex(NDArray &a, NDArray &b) { return res; } - TEST_CASE("Compare two images") { NDArray a; NDArray b; @@ -222,7 +221,6 @@ TEST_CASE("Bitwise and on data") { REQUIRE(a(2) == 384); } - TEST_CASE("Elementwise operations on images") { std::array shape{5, 5}; double a_val = 3.0; @@ -258,7 +256,8 @@ TEST_CASE("Elementwise operations on images") { NDArray A(shape, a_val); NDArray B(shape, b_val); NDArray C = A - B; - // auto C = A - B; // This works but the result is a lazy ArraySub object + // auto C = A - B; // This works but the result is a lazy ArraySub + // object // Value of C matches for (uint32_t i = 0; i < C.size(); ++i) { @@ -282,7 +281,8 @@ TEST_CASE("Elementwise operations on images") { SECTION("Multiply two images") { NDArray A(shape, a_val); NDArray B(shape, b_val); - // auto C = A * B; // This works but the result is a lazy ArrayMul object + // auto C = A * B; // This works but the result is a lazy ArrayMul + // object NDArray C = A * B; // Value of C matches @@ -307,7 +307,8 @@ TEST_CASE("Elementwise operations on images") { SECTION("Divide two images") { NDArray A(shape, a_val); NDArray B(shape, b_val); - // auto C = A / B; // This works but the result is a lazy ArrayDiv object + // auto C = A / B; // This works but the result is a lazy ArrayDiv + // object NDArray C = A / B; // Value of C matches diff --git a/src/NDView.test.cpp b/src/NDView.test.cpp index 89e76e9..a65758c 100644 --- a/src/NDView.test.cpp +++ b/src/NDView.test.cpp @@ -2,8 +2,8 @@ #include #include -#include #include +#include using aare::NDView; using aare::Shape; @@ -151,8 +151,10 @@ TEST_CASE("divide with another span") { std::vector vec1{3, 2, 1}; std::vector result{3, 6, 3}; - NDView data0(vec0.data(), Shape<1>{static_cast(vec0.size())}); - NDView data1(vec1.data(), Shape<1>{static_cast(vec1.size())}); + NDView data0(vec0.data(), + Shape<1>{static_cast(vec0.size())}); + NDView data1(vec1.data(), + Shape<1>{static_cast(vec1.size())}); data0 /= data1; @@ -181,8 +183,7 @@ TEST_CASE("compare two views") { REQUIRE((view1 == view2)); } - -TEST_CASE("Create a view over a vector"){ +TEST_CASE("Create a view over a vector") { std::vector vec(12); std::iota(vec.begin(), vec.end(), 0); auto v = aare::make_view(vec); diff --git a/src/NumpyFile.cpp b/src/NumpyFile.cpp index e375ce3..4e0c215 100644 --- a/src/NumpyFile.cpp +++ b/src/NumpyFile.cpp @@ -4,16 +4,16 @@ namespace aare { - - -NumpyFile::NumpyFile(const std::filesystem::path &fname, const std::string &mode, FileConfig cfg) { +NumpyFile::NumpyFile(const std::filesystem::path &fname, + const std::string &mode, FileConfig cfg) { // TODO! add opts to constructor m_mode = mode; if (mode == "r") { fp = fopen(fname.string().c_str(), "rb"); if (!fp) { - throw std::runtime_error(fmt::format("Could not open: {} for reading", fname.string())); + throw std::runtime_error( + fmt::format("Could not open: {} for reading", fname.string())); } load_metadata(); } else if (mode == "w") { @@ -24,11 +24,15 @@ NumpyFile::NumpyFile(const std::filesystem::path &fname, const std::string &mode m_header.shape = {0, cfg.rows, cfg.cols}; fp = fopen(fname.string().c_str(), "wb"); if (!fp) { - throw std::runtime_error(fmt::format("Could not open: {} for reading", fname.string())); + throw std::runtime_error( + fmt::format("Could not open: {} for reading", fname.string())); } - initial_header_len = aare::NumpyHelpers::write_header(std::filesystem::path(fname.c_str()), m_header); + initial_header_len = aare::NumpyHelpers::write_header( + std::filesystem::path(fname.c_str()), m_header); } - m_pixels_per_frame = std::accumulate(m_header.shape.begin() + 1, m_header.shape.end(), 1, std::multiplies<>()); + m_pixels_per_frame = + std::accumulate(m_header.shape.begin() + 1, m_header.shape.end(), 1, + std::multiplies<>()); m_bytes_per_frame = m_header.dtype.bitdepth() / 8 * m_pixels_per_frame; } @@ -63,7 +67,8 @@ void NumpyFile::get_frame_into(size_t frame_number, std::byte *image_buf) { if (frame_number > m_header.shape[0]) { throw std::invalid_argument("Frame number out of range"); } - if (fseek(fp, header_size + frame_number * m_bytes_per_frame, SEEK_SET)) // NOLINT + if (fseek(fp, header_size + frame_number * m_bytes_per_frame, + SEEK_SET)) // NOLINT throw std::runtime_error("Could not seek to frame"); size_t const rc = fread(image_buf, m_bytes_per_frame, 1, fp); @@ -113,7 +118,8 @@ NumpyFile::~NumpyFile() noexcept { // write header size_t const rc = fwrite(header_str.c_str(), header_str.size(), 1, fp); if (rc != 1) { - std::cout << "Error writing header to numpy file in destructor" << std::endl; + std::cout << "Error writing header to numpy file in destructor" + << std::endl; } } @@ -140,8 +146,10 @@ void NumpyFile::load_metadata() { } // read version - rc = fread(reinterpret_cast(&major_ver_), sizeof(major_ver_), 1, fp); - rc += fread(reinterpret_cast(&minor_ver_), sizeof(minor_ver_), 1, fp); + rc = + fread(reinterpret_cast(&major_ver_), sizeof(major_ver_), 1, fp); + rc += + fread(reinterpret_cast(&minor_ver_), sizeof(minor_ver_), 1, fp); if (rc != 2) { throw std::runtime_error("Error reading numpy version"); } @@ -159,7 +167,8 @@ void NumpyFile::load_metadata() { if (rc != 1) { throw std::runtime_error("Error reading header length"); } - header_size = aare::NumpyHelpers::magic_string_length + 2 + header_len_size + header_len; + header_size = aare::NumpyHelpers::magic_string_length + 2 + + header_len_size + header_len; if (header_size % 16 != 0) { fmt::print("Warning: header length is not a multiple of 16\n"); } diff --git a/src/NumpyFile.test.cpp b/src/NumpyFile.test.cpp index bdd6451..e687bea 100644 --- a/src/NumpyFile.test.cpp +++ b/src/NumpyFile.test.cpp @@ -1,8 +1,8 @@ #include "aare/NumpyFile.hpp" #include "aare/NDArray.hpp" -#include #include "test_config.hpp" +#include using aare::Dtype; using aare::NumpyFile; @@ -23,7 +23,7 @@ TEST_CASE("Read a 1D numpy file with int32 data type", "[.integration]") { REQUIRE(data(i) == i); } } - + TEST_CASE("Read a 3D numpy file with np.double data type", "[.integration]") { auto fpath = test_data_path() / "numpy" / "test_3d_double.npy"; diff --git a/src/NumpyHelpers.cpp b/src/NumpyHelpers.cpp index b8414d9..106327e 100644 --- a/src/NumpyHelpers.cpp +++ b/src/NumpyHelpers.cpp @@ -29,7 +29,8 @@ namespace aare { std::string NumpyHeader::to_string() const { std::stringstream sstm; - sstm << "dtype: " << dtype.to_string() << ", fortran_order: " << fortran_order << ' '; + sstm << "dtype: " << dtype.to_string() + << ", fortran_order: " << fortran_order << ' '; sstm << "shape: ("; for (auto item : shape) sstm << item << ','; @@ -37,10 +38,10 @@ std::string NumpyHeader::to_string() const { return sstm.str(); } - namespace NumpyHelpers { -std::unordered_map parse_dict(std::string in, const std::vector &keys) { +std::unordered_map +parse_dict(std::string in, const std::vector &keys) { std::unordered_map map; if (keys.empty()) return map; @@ -100,7 +101,8 @@ aare::Dtype parse_descr(std::string typestring) { constexpr char little_endian_char = '<'; constexpr char big_endian_char = '>'; constexpr char no_endian_char = '|'; - constexpr std::array endian_chars = {little_endian_char, big_endian_char, no_endian_char}; + constexpr std::array endian_chars = { + little_endian_char, big_endian_char, no_endian_char}; constexpr std::array numtype_chars = {'f', 'i', 'u', 'c'}; const char byteorder_c = typestring[0]; @@ -139,7 +141,9 @@ std::string get_value_from_map(const std::string &mapstr) { return trim(tmp); } -bool is_digits(const std::string &str) { return std::all_of(str.begin(), str.end(), ::isdigit); } +bool is_digits(const std::string &str) { + return std::all_of(str.begin(), str.end(), ::isdigit); +} std::vector parse_tuple(std::string in) { std::vector v; @@ -215,20 +219,25 @@ inline std::string write_boolean(bool b) { return "False"; } -inline std::string write_header_dict(const std::string &descr, bool fortran_order, const std::vector &shape) { +inline std::string write_header_dict(const std::string &descr, + bool fortran_order, + const std::vector &shape) { std::string const s_fortran_order = write_boolean(fortran_order); std::string const shape_s = write_tuple(shape); - return "{'descr': '" + descr + "', 'fortran_order': " + s_fortran_order + ", 'shape': " + shape_s + ", }"; + return "{'descr': '" + descr + "', 'fortran_order': " + s_fortran_order + + ", 'shape': " + shape_s + ", }"; } -size_t write_header(const std::filesystem::path &fname, const NumpyHeader &header) { +size_t write_header(const std::filesystem::path &fname, + const NumpyHeader &header) { std::ofstream out(fname, std::ios::binary | std::ios::out); return write_header(out, header); } size_t write_header(std::ostream &out, const NumpyHeader &header) { - std::string const header_dict = write_header_dict(header.dtype.to_string(), header.fortran_order, header.shape); + std::string const header_dict = write_header_dict( + header.dtype.to_string(), header.fortran_order, header.shape); size_t length = magic_string_length + 2 + 2 + header_dict.length() + 1; @@ -247,17 +256,22 @@ size_t write_header(std::ostream &out, const NumpyHeader &header) { // write header length if (version_major == 1 && version_minor == 0) { - auto header_len = static_cast(header_dict.length() + padding.length() + 1); + auto header_len = + static_cast(header_dict.length() + padding.length() + 1); - std::array header_len_le16{static_cast((header_len >> 0) & 0xff), - static_cast((header_len >> 8) & 0xff)}; + std::array header_len_le16{ + static_cast((header_len >> 0) & 0xff), + static_cast((header_len >> 8) & 0xff)}; out.write(reinterpret_cast(header_len_le16.data()), 2); } else { - auto header_len = static_cast(header_dict.length() + padding.length() + 1); + auto header_len = + static_cast(header_dict.length() + padding.length() + 1); std::array header_len_le32{ - static_cast((header_len >> 0) & 0xff), static_cast((header_len >> 8) & 0xff), - static_cast((header_len >> 16) & 0xff), static_cast((header_len >> 24) & 0xff)}; + static_cast((header_len >> 0) & 0xff), + static_cast((header_len >> 8) & 0xff), + static_cast((header_len >> 16) & 0xff), + static_cast((header_len >> 24) & 0xff)}; out.write(reinterpret_cast(header_len_le32.data()), 4); } diff --git a/src/NumpyHelpers.test.cpp b/src/NumpyHelpers.test.cpp index 36fcfe6..ad55b17 100644 --- a/src/NumpyHelpers.test.cpp +++ b/src/NumpyHelpers.test.cpp @@ -19,7 +19,9 @@ TEST_CASE("Check for quotes and return stripped string") { REQUIRE(parse_str("''") == ""); } -TEST_CASE("parsing a string without quotes throws") { REQUIRE_THROWS(parse_str("hej")); } +TEST_CASE("parsing a string without quotes throws") { + REQUIRE_THROWS(parse_str("hej")); +} TEST_CASE("trim whitespace") { REQUIRE(trim(" hej ") == "hej"); @@ -53,7 +55,8 @@ TEST_CASE("is element in array") { } TEST_CASE("Parse numpy dict") { - std::string in = "{'descr': ' keys{"descr", "fortran_order", "shape"}; auto map = parse_dict(in, keys); REQUIRE(map["descr"] == "' #include +#include #include #include @@ -58,7 +57,8 @@ TEST_CASE("test pedestal push") { if (k < 5) { REQUIRE(pedestal.cur_samples()(i, j) == k + 1); REQUIRE(pedestal.get_sum()(i, j) == (k + 1) * (i + j)); - REQUIRE(pedestal.get_sum2()(i, j) == (k + 1) * (i + j) * (i + j)); + REQUIRE(pedestal.get_sum2()(i, j) == + (k + 1) * (i + j) * (i + j)); } else { REQUIRE(pedestal.cur_samples()(i, j) == 5); REQUIRE(pedestal.get_sum()(i, j) == 5 * (i + j)); @@ -95,9 +95,12 @@ TEST_CASE("test pedestal with normal distribution") { for (int i = 0; i < 3; i++) { for (int j = 0; j < 5; j++) { - REQUIRE_THAT(mean(i, j), Catch::Matchers::WithinAbs(MEAN, MEAN * TOLERANCE)); - REQUIRE_THAT(variance(i, j), Catch::Matchers::WithinAbs(VAR, VAR * TOLERANCE)); - REQUIRE_THAT(standard_deviation(i, j), Catch::Matchers::WithinAbs(STD, STD * TOLERANCE)); + REQUIRE_THAT(mean(i, j), + Catch::Matchers::WithinAbs(MEAN, MEAN * TOLERANCE)); + REQUIRE_THAT(variance(i, j), + Catch::Matchers::WithinAbs(VAR, VAR * TOLERANCE)); + REQUIRE_THAT(standard_deviation(i, j), + Catch::Matchers::WithinAbs(STD, STD * TOLERANCE)); } } } \ No newline at end of file diff --git a/src/PixelMap.cpp b/src/PixelMap.cpp index d62759a..5be5ed4 100644 --- a/src/PixelMap.cpp +++ b/src/PixelMap.cpp @@ -31,7 +31,7 @@ NDArray GenerateMoench03PixelMap() { } NDArray GenerateMoench05PixelMap() { - std::array adc_numbers = {5, 9, 1}; + std::array adc_numbers = {5, 9, 1}; NDArray order_map({160, 150}); int n_pixel = 0; for (int row = 0; row < 160; row++) { @@ -40,11 +40,11 @@ NDArray GenerateMoench05PixelMap() { for (int i_sc = 0; i_sc < 3; i_sc++) { int col = 50 * i_sc + i_col; int adc_nr = adc_numbers[i_sc]; - int i_analog = n_pixel * 12 + adc_nr; + int i_analog = n_pixel * 12 + adc_nr; - // analog_frame[row * 150 + col] = analog_data[i_analog] & 0x3FFF; + // analog_frame[row * 150 + col] = analog_data[i_analog] & + // 0x3FFF; order_map(row, col) = i_analog; - } } } @@ -52,7 +52,7 @@ NDArray GenerateMoench05PixelMap() { } NDArray GenerateMoench05PixelMap1g() { - std::array adc_numbers = {1, 2, 0}; + std::array adc_numbers = {1, 2, 0}; NDArray order_map({160, 150}); int n_pixel = 0; for (int row = 0; row < 160; row++) { @@ -61,12 +61,11 @@ NDArray GenerateMoench05PixelMap1g() { for (int i_sc = 0; i_sc < 3; i_sc++) { int col = 50 * i_sc + i_col; int adc_nr = adc_numbers[i_sc]; - int i_analog = n_pixel * 3 + adc_nr; + int i_analog = n_pixel * 3 + adc_nr; - - // analog_frame[row * 150 + col] = analog_data[i_analog] & 0x3FFF; + // analog_frame[row * 150 + col] = analog_data[i_analog] & + // 0x3FFF; order_map(row, col) = i_analog; - } } } @@ -85,42 +84,42 @@ NDArray GenerateMoench05PixelMapOld() { int adc_nr = adc_numbers[i_sc]; int i_analog = n_pixel * 32 + adc_nr; - - // analog_frame[row * 150 + col] = analog_data[i_analog] & 0x3FFF; + // analog_frame[row * 150 + col] = analog_data[i_analog] & + // 0x3FFF; order_map(row, col) = i_analog; - } } } return order_map; } -NDArrayGenerateEigerFlipRowsPixelMap(){ +NDArray GenerateEigerFlipRowsPixelMap() { NDArray order_map({256, 512}); - for(int row = 0; row < 256; row++){ - for(int col = 0; col < 512; col++){ - order_map(row, col) = 255*512-row*512 + col; + for (int row = 0; row < 256; row++) { + for (int col = 0; col < 512; col++) { + order_map(row, col) = 255 * 512 - row * 512 + col; } } return order_map; } -NDArrayGenerateMH02SingleCounterPixelMap(){ +NDArray GenerateMH02SingleCounterPixelMap() { NDArray order_map({48, 48}); - for(int row = 0; row < 48; row++){ - for(int col = 0; col < 48; col++){ - order_map(row, col) = row*48 + col; + for (int row = 0; row < 48; row++) { + for (int col = 0; col < 48; col++) { + order_map(row, col) = row * 48 + col; } } return order_map; } -NDArray GenerateMH02FourCounterPixelMap(){ +NDArray GenerateMH02FourCounterPixelMap() { NDArray order_map({4, 48, 48}); - for (int counter=0; counter<4; counter++){ - for(int row = 0; row < 48; row++){ - for(int col = 0; col < 48; col++){ - order_map(counter, row, col) = counter*48*48 + row*48 + col; + for (int counter = 0; counter < 4; counter++) { + for (int row = 0; row < 48; row++) { + for (int col = 0; col < 48; col++) { + order_map(counter, row, col) = + counter * 48 * 48 + row * 48 + col; } } } diff --git a/src/RawFile.cpp b/src/RawFile.cpp index c576453..22bf8dd 100644 --- a/src/RawFile.cpp +++ b/src/RawFile.cpp @@ -1,7 +1,9 @@ #include "aare/RawFile.hpp" #include "aare/PixelMap.hpp" +#include "aare/algorithm.hpp" #include "aare/defs.hpp" #include "aare/geo_helpers.hpp" +#include "aare/logger.hpp" #include #include @@ -14,23 +16,15 @@ RawFile::RawFile(const std::filesystem::path &fname, const std::string &mode) : m_master(fname) { m_mode = mode; if (mode == "r") { - - n_subfiles = find_number_of_subfiles(); // f0,f1...fn - n_subfile_parts = - m_master.geometry().col * m_master.geometry().row; // d0,d1...dn - - - find_geometry(); - - if (m_master.roi()){ - m_geometry = update_geometry_with_roi(m_geometry, m_master.roi().value()); + if (m_master.roi()) { + m_geometry = + update_geometry_with_roi(m_geometry, m_master.roi().value()); } - open_subfiles(); } else { throw std::runtime_error(LOCATION + - "Unsupported mode. Can only read RawFiles."); + " Unsupported mode. Can only read RawFiles."); } } @@ -54,32 +48,31 @@ void RawFile::read_into(std::byte *image_buf) { return get_frame_into(m_current_frame++, image_buf); } - void RawFile::read_into(std::byte *image_buf, DetectorHeader *header) { return get_frame_into(m_current_frame++, image_buf, header); } -void RawFile::read_into(std::byte *image_buf, size_t n_frames, DetectorHeader *header) { +void RawFile::read_into(std::byte *image_buf, size_t n_frames, + DetectorHeader *header) { // return get_frame_into(m_current_frame++, image_buf, header); for (size_t i = 0; i < n_frames; i++) { this->get_frame_into(m_current_frame++, image_buf, header); image_buf += bytes_per_frame(); - if(header) - header+=n_mod(); + if (header) + header += n_modules(); } - } -size_t RawFile::n_mod() const { return n_subfile_parts; } - +size_t RawFile::n_modules() const { return m_master.n_modules(); } size_t RawFile::bytes_per_frame() { - return m_geometry.pixels_x * m_geometry.pixels_y * m_master.bitdepth() / bits_per_byte; + return m_geometry.pixels_x * m_geometry.pixels_y * m_master.bitdepth() / + bits_per_byte; } -size_t RawFile::pixels_per_frame() { - // return m_rows * m_cols; +size_t RawFile::pixels_per_frame() { + // return m_rows * m_cols; return m_geometry.pixels_x * m_geometry.pixels_y; } @@ -106,17 +99,11 @@ xy RawFile::geometry() { return m_master.geometry(); } void RawFile::open_subfiles() { if (m_mode == "r") - for (size_t i = 0; i != n_subfiles; ++i) { - auto v = std::vector(n_subfile_parts); - for (size_t j = 0; j != n_subfile_parts; ++j) { - auto pos = m_geometry.module_pixel_0[j]; - v[j] = new RawSubFile(m_master.data_fname(j, i), - m_master.detector_type(), pos.height, - pos.width, m_master.bitdepth(), - pos.row_index, pos.col_index); - - } - subfiles.push_back(v); + for (size_t i = 0; i != n_modules(); ++i) { + auto pos = m_geometry.module_pixel_0[i]; + m_subfiles.emplace_back(std::make_unique( + m_master.data_fname(i, 0), m_master.detector_type(), pos.height, + pos.width, m_master.bitdepth(), pos.row_index, pos.col_index)); } else { throw std::runtime_error(LOCATION + @@ -141,39 +128,25 @@ DetectorHeader RawFile::read_header(const std::filesystem::path &fname) { return h; } -int RawFile::find_number_of_subfiles() { - int n_files = 0; - // f0,f1...fn How many files is the data split into? - while (std::filesystem::exists(m_master.data_fname(0, n_files))) - n_files++; // increment after test - -#ifdef AARE_VERBOSE - fmt::print("Found: {} subfiles\n", n_files); -#endif - return n_files; - -} - RawMasterFile RawFile::master() const { return m_master; } /** * @brief Find the geometry of the detector by opening all the subfiles and - * reading the headers. + * reading the headers. */ void RawFile::find_geometry() { - - //Hold the maximal row and column number found - //Later used for calculating the total number of rows and columns + + // Hold the maximal row and column number found + // Later used for calculating the total number of rows and columns uint16_t r{}; uint16_t c{}; - - for (size_t i = 0; i < n_subfile_parts; i++) { + for (size_t i = 0; i < n_modules(); i++) { auto h = read_header(m_master.data_fname(i, 0)); r = std::max(r, h.row); c = std::max(c, h.column); // positions.push_back({h.row, h.column}); - + ModuleGeometry g; g.origin_x = h.column * m_master.pixels_x(); g.origin_y = h.row * m_master.pixels_y(); @@ -182,127 +155,110 @@ void RawFile::find_geometry() { g.width = m_master.pixels_x(); g.height = m_master.pixels_y(); m_geometry.module_pixel_0.push_back(g); - } r++; c++; m_geometry.pixels_y = (r * m_master.pixels_y()); - m_geometry.pixels_x = (c * m_master.pixels_x()); + m_geometry.pixels_x = (c * m_master.pixels_x()); m_geometry.modules_x = c; m_geometry.modules_y = r; m_geometry.pixels_y += static_cast((r - 1) * cfg.module_gap_row); - } - Frame RawFile::get_frame(size_t frame_index) { - auto f = Frame(m_geometry.pixels_y, m_geometry.pixels_x, Dtype::from_bitdepth(m_master.bitdepth())); + auto f = Frame(m_geometry.pixels_y, m_geometry.pixels_x, + Dtype::from_bitdepth(m_master.bitdepth())); std::byte *frame_buffer = f.data(); get_frame_into(frame_index, frame_buffer); return f; } +size_t RawFile::bytes_per_pixel() const { return m_master.bitdepth() / 8; } -size_t RawFile::bytes_per_pixel() const { - return m_master.bitdepth() / 8; -} - -void RawFile::get_frame_into(size_t frame_index, std::byte *frame_buffer, DetectorHeader *header) { +void RawFile::get_frame_into(size_t frame_index, std::byte *frame_buffer, + DetectorHeader *header) { + LOG(logDEBUG) << "RawFile::get_frame_into(" << frame_index << ")"; if (frame_index >= total_frames()) { throw std::runtime_error(LOCATION + "Frame number out of range"); } - std::vector frame_numbers(n_subfile_parts); - std::vector frame_indices(n_subfile_parts, frame_index); - + std::vector frame_numbers(n_modules()); + std::vector frame_indices(n_modules(), frame_index); // sync the frame numbers - if (n_subfile_parts != 1) { - for (size_t part_idx = 0; part_idx != n_subfile_parts; ++part_idx) { - auto subfile_id = frame_index / m_master.max_frames_per_file(); - if (subfile_id >= subfiles.size()) { - throw std::runtime_error(LOCATION + - " Subfile out of range. Possible missing data."); - } + if (n_modules() != 1) { // if we have more than one module + for (size_t part_idx = 0; part_idx != n_modules(); ++part_idx) { frame_numbers[part_idx] = - subfiles[subfile_id][part_idx]->frame_number( - frame_index % m_master.max_frames_per_file()); + m_subfiles[part_idx]->frame_number(frame_index); } + // 1. if frame number vector is the same break - while (std::adjacent_find(frame_numbers.begin(), frame_numbers.end(), - std::not_equal_to<>()) != - frame_numbers.end()) { + while (!all_equal(frame_numbers)) { + // 2. find the index of the minimum frame number, auto min_frame_idx = std::distance( frame_numbers.begin(), std::min_element(frame_numbers.begin(), frame_numbers.end())); + // 3. increase its index and update its respective frame number frame_indices[min_frame_idx]++; + // 4. if we can't increase its index => throw error if (frame_indices[min_frame_idx] >= total_frames()) { throw std::runtime_error(LOCATION + "Frame number out of range"); } - auto subfile_id = - frame_indices[min_frame_idx] / m_master.max_frames_per_file(); + frame_numbers[min_frame_idx] = - subfiles[subfile_id][min_frame_idx]->frame_number( - frame_indices[min_frame_idx] % - m_master.max_frames_per_file()); + m_subfiles[min_frame_idx]->frame_number( + frame_indices[min_frame_idx]); } } if (m_master.geometry().col == 1) { // get the part from each subfile and copy it to the frame - for (size_t part_idx = 0; part_idx != n_subfile_parts; ++part_idx) { + for (size_t part_idx = 0; part_idx != n_modules(); ++part_idx) { auto corrected_idx = frame_indices[part_idx]; - auto subfile_id = corrected_idx / m_master.max_frames_per_file(); - if (subfile_id >= subfiles.size()) { - throw std::runtime_error(LOCATION + - " Subfile out of range. Possible missing data."); - } // This is where we start writing - auto offset = (m_geometry.module_pixel_0[part_idx].origin_y * m_geometry.pixels_x + - m_geometry.module_pixel_0[part_idx].origin_x)*m_master.bitdepth()/8; + auto offset = (m_geometry.module_pixel_0[part_idx].origin_y * + m_geometry.pixels_x + + m_geometry.module_pixel_0[part_idx].origin_x) * + m_master.bitdepth() / 8; - if (m_geometry.module_pixel_0[part_idx].origin_x!=0) - throw std::runtime_error(LOCATION + "Implementation error. x pos not 0."); - - //TODO! Risk for out of range access - subfiles[subfile_id][part_idx]->seek(corrected_idx % m_master.max_frames_per_file()); - subfiles[subfile_id][part_idx]->read_into(frame_buffer + offset, header); + if (m_geometry.module_pixel_0[part_idx].origin_x != 0) + throw std::runtime_error(LOCATION + + " Implementation error. x pos not 0."); + + // TODO! What if the files don't match? + m_subfiles[part_idx]->seek(corrected_idx); + m_subfiles[part_idx]->read_into(frame_buffer + offset, header); if (header) ++header; } } else { - //TODO! should we read row by row? + // TODO! should we read row by row? // create a buffer large enough to hold a full module - auto bytes_per_part = m_master.pixels_y() * m_master.pixels_x() * m_master.bitdepth() / 8; // TODO! replace with image_size_in_bytes + auto *part_buffer = new std::byte[bytes_per_part]; // TODO! if we have many submodules we should reorder them on the module // level - for (size_t part_idx = 0; part_idx != n_subfile_parts; ++part_idx) { + for (size_t part_idx = 0; part_idx != n_modules(); ++part_idx) { auto pos = m_geometry.module_pixel_0[part_idx]; auto corrected_idx = frame_indices[part_idx]; - auto subfile_id = corrected_idx / m_master.max_frames_per_file(); - if (subfile_id >= subfiles.size()) { - throw std::runtime_error(LOCATION + - " Subfile out of range. Possible missing data."); - } - subfiles[subfile_id][part_idx]->seek(corrected_idx % m_master.max_frames_per_file()); - subfiles[subfile_id][part_idx]->read_into(part_buffer, header); - if(header) + m_subfiles[part_idx]->seek(corrected_idx); + m_subfiles[part_idx]->read_into(part_buffer, header); + if (header) ++header; for (size_t cur_row = 0; cur_row < static_cast(pos.height); @@ -313,10 +269,9 @@ void RawFile::get_frame_into(size_t frame_index, std::byte *frame_buffer, Detect auto dest = (irow * this->m_geometry.pixels_x + icol); dest = dest * m_master.bitdepth() / 8; memcpy(frame_buffer + dest, - part_buffer + cur_row * pos.width * - m_master.bitdepth() / 8, + part_buffer + + cur_row * pos.width * m_master.bitdepth() / 8, pos.width * m_master.bitdepth() / 8); - } } delete[] part_buffer; @@ -337,27 +292,7 @@ size_t RawFile::frame_number(size_t frame_index) { if (frame_index >= m_master.frames_in_file()) { throw std::runtime_error(LOCATION + " Frame number out of range"); } - size_t subfile_id = frame_index / m_master.max_frames_per_file(); - if (subfile_id >= subfiles.size()) { - throw std::runtime_error( - LOCATION + " Subfile out of range. Possible missing data."); - } - return subfiles[subfile_id][0]->frame_number( - frame_index % m_master.max_frames_per_file()); + return m_subfiles[0]->frame_number(frame_index); } -RawFile::~RawFile() { - - // TODO! Fix this, for file closing - for (auto &vec : subfiles) { - for (auto *subfile : vec) { - delete subfile; - } - } -} - - - - - } // namespace aare diff --git a/src/RawFile.test.cpp b/src/RawFile.test.cpp index 5f9b2e1..1fb441f 100644 --- a/src/RawFile.test.cpp +++ b/src/RawFile.test.cpp @@ -1,18 +1,18 @@ +#include "aare/RawFile.hpp" #include "aare/File.hpp" #include "aare/RawMasterFile.hpp" //needed for ROI -#include "aare/RawFile.hpp" #include #include #include "test_config.hpp" - using aare::File; TEST_CASE("Read number of frames from a jungfrau raw file", "[.integration]") { - auto fpath = test_data_path() / "jungfrau" / "jungfrau_single_master_0.json"; + auto fpath = + test_data_path() / "jungfrau" / "jungfrau_single_master_0.json"; REQUIRE(std::filesystem::exists(fpath)); File f(fpath, "r"); @@ -20,7 +20,8 @@ TEST_CASE("Read number of frames from a jungfrau raw file", "[.integration]") { } TEST_CASE("Read frame numbers from a jungfrau raw file", "[.integration]") { - auto fpath = test_data_path() / "jungfrau" / "jungfrau_single_master_0.json"; + auto fpath = + test_data_path() / "jungfrau" / "jungfrau_single_master_0.json"; REQUIRE(std::filesystem::exists(fpath)); File f(fpath, "r"); @@ -36,7 +37,8 @@ TEST_CASE("Read frame numbers from a jungfrau raw file", "[.integration]") { } TEST_CASE("Read a frame number too high throws", "[.integration]") { - auto fpath = test_data_path() / "jungfrau" / "jungfrau_single_master_0.json"; + auto fpath = + test_data_path() / "jungfrau" / "jungfrau_single_master_0.json"; REQUIRE(std::filesystem::exists(fpath)); File f(fpath, "r"); @@ -49,8 +51,10 @@ TEST_CASE("Read a frame number too high throws", "[.integration]") { REQUIRE_THROWS(f.frame_number(10)); } -TEST_CASE("Read a frame numbers where the subfile is missing throws", "[.integration]") { - auto fpath = test_data_path() / "jungfrau" / "jungfrau_missing_subfile_master_0.json"; +TEST_CASE("Read a frame numbers where the subfile is missing throws", + "[.integration]") { + auto fpath = test_data_path() / "jungfrau" / + "jungfrau_missing_subfile_master_0.json"; REQUIRE(std::filesystem::exists(fpath)); File f(fpath, "r"); @@ -58,7 +62,7 @@ TEST_CASE("Read a frame numbers where the subfile is missing throws", "[.integra // we know this file has 10 frames with frame numbers 1 to 10 // f0 1,2,3 // f1 4,5,6 - but files f1-f3 are missing - // f2 7,8,9 - gone + // f2 7,8,9 - gone // f3 10 - gone REQUIRE(f.frame_number(0) == 1); REQUIRE(f.frame_number(1) == 2); @@ -69,15 +73,18 @@ TEST_CASE("Read a frame numbers where the subfile is missing throws", "[.integra REQUIRE_THROWS(f.frame_number(10)); } - -TEST_CASE("Read data from a jungfrau 500k single port raw file", "[.integration]") { - auto fpath = test_data_path() / "jungfrau" / "jungfrau_single_master_0.json"; +TEST_CASE("Read data from a jungfrau 500k single port raw file", + "[.integration]") { + auto fpath = + test_data_path() / "jungfrau" / "jungfrau_single_master_0.json"; REQUIRE(std::filesystem::exists(fpath)); File f(fpath, "r"); - // we know this file has 10 frames with pixel 0,0 being: 2123, 2051, 2109, 2117, 2089, 2095, 2072, 2126, 2097, 2102 - std::vector pixel_0_0 = {2123, 2051, 2109, 2117, 2089, 2095, 2072, 2126, 2097, 2102}; + // we know this file has 10 frames with pixel 0,0 being: 2123, 2051, 2109, + // 2117, 2089, 2095, 2072, 2126, 2097, 2102 + std::vector pixel_0_0 = {2123, 2051, 2109, 2117, 2089, + 2095, 2072, 2126, 2097, 2102}; for (size_t i = 0; i < 10; i++) { auto frame = f.read_frame(); CHECK(frame.rows() == 512); @@ -99,11 +106,13 @@ TEST_CASE("Read frame numbers from a raw file", "[.integration]") { } } -TEST_CASE("Compare reading from a numpy file with a raw file", "[.integration]") { - auto fpath_raw = test_data_path() / "jungfrau" / "jungfrau_single_master_0.json"; +TEST_CASE("Compare reading from a numpy file with a raw file", "[.files]") { + auto fpath_raw = + test_data_path() / "raw/jungfrau" / "jungfrau_single_master_0.json"; REQUIRE(std::filesystem::exists(fpath_raw)); - auto fpath_npy = test_data_path() / "jungfrau" / "jungfrau_single_0.npy"; + auto fpath_npy = + test_data_path() / "raw/jungfrau" / "jungfrau_single_0.npy"; REQUIRE(std::filesystem::exists(fpath_npy)); File raw(fpath_raw, "r"); @@ -113,6 +122,7 @@ TEST_CASE("Compare reading from a numpy file with a raw file", "[.integration]") CHECK(npy.total_frames() == 10); for (size_t i = 0; i < 10; ++i) { + CHECK(raw.tell() == i); auto raw_frame = raw.read_frame(); auto npy_frame = npy.read_frame(); CHECK((raw_frame.view() == npy_frame.view())); @@ -120,17 +130,23 @@ TEST_CASE("Compare reading from a numpy file with a raw file", "[.integration]") } TEST_CASE("Read multipart files", "[.integration]") { - auto fpath = test_data_path() / "jungfrau" / "jungfrau_double_master_0.json"; + auto fpath = + test_data_path() / "jungfrau" / "jungfrau_double_master_0.json"; REQUIRE(std::filesystem::exists(fpath)); File f(fpath, "r"); // we know this file has 10 frames check read_multiport.py for the values - std::vector pixel_0_0 = {2099, 2121, 2108, 2084, 2084, 2118, 2066, 2108, 2112, 2116}; - std::vector pixel_0_1 = {2842, 2796, 2865, 2798, 2805, 2817, 2852, 2789, 2792, 2833}; - std::vector pixel_255_1023 = {2149, 2037, 2115, 2102, 2118, 2090, 2036, 2071, 2073, 2142}; - std::vector pixel_511_1023 = {3231, 3169, 3167, 3162, 3168, 3160, 3171, 3171, 3169, 3171}; - std::vector pixel_1_0 = {2748, 2614, 2665, 2629, 2618, 2630, 2631, 2634, 2577, 2598}; + std::vector pixel_0_0 = {2099, 2121, 2108, 2084, 2084, + 2118, 2066, 2108, 2112, 2116}; + std::vector pixel_0_1 = {2842, 2796, 2865, 2798, 2805, + 2817, 2852, 2789, 2792, 2833}; + std::vector pixel_255_1023 = {2149, 2037, 2115, 2102, 2118, + 2090, 2036, 2071, 2073, 2142}; + std::vector pixel_511_1023 = {3231, 3169, 3167, 3162, 3168, + 3160, 3171, 3171, 3169, 3171}; + std::vector pixel_1_0 = {2748, 2614, 2665, 2629, 2618, + 2630, 2631, 2634, 2577, 2598}; for (size_t i = 0; i < 10; i++) { auto frame = f.read_frame(); @@ -145,11 +161,9 @@ TEST_CASE("Read multipart files", "[.integration]") { } TEST_CASE("Read file with unordered frames", "[.integration]") { - //TODO! Better explanation and error message + // TODO! Better explanation and error message auto fpath = test_data_path() / "mythen" / "scan242_master_3.raw"; REQUIRE(std::filesystem::exists(fpath)); File f(fpath); REQUIRE_THROWS((f.read_frame())); } - - diff --git a/src/RawMasterFile.cpp b/src/RawMasterFile.cpp index 33807d4..508b396 100644 --- a/src/RawMasterFile.cpp +++ b/src/RawMasterFile.cpp @@ -1,5 +1,5 @@ #include "aare/RawMasterFile.hpp" -#include +#include namespace aare { RawFileNameComponents::RawFileNameComponents( @@ -37,18 +37,15 @@ std::filesystem::path RawFileNameComponents::master_fname() const { } std::filesystem::path RawFileNameComponents::data_fname(size_t mod_id, - size_t file_id - ) const { - - - + size_t file_id) const { + std::string fmt = "{}_d{}_f{}_{}.raw"; - //Before version X we used to name the data files f000000000000 + // Before version X we used to name the data files f000000000000 if (m_old_scheme) { fmt = "{}_d{}_f{:012}_{}.raw"; } - return m_base_path / fmt::format(fmt, m_base_name, mod_id, - file_id, m_file_index); + return m_base_path / + fmt::format(fmt, m_base_name, mod_id, file_id, m_file_index); } void RawFileNameComponents::set_old_scheme(bool old_scheme) { @@ -65,19 +62,19 @@ const std::string &RawFileNameComponents::ext() const { return m_ext; } int RawFileNameComponents::file_index() const { return m_file_index; } // "[enabled\ndac dac 4\nstart 500\nstop 2200\nstep 5\nsettleTime 100us\n]" -ScanParameters::ScanParameters(const std::string& par){ - std::istringstream iss(par.substr(1, par.size()-2)); +ScanParameters::ScanParameters(const std::string &par) { + std::istringstream iss(par.substr(1, par.size() - 2)); std::string line; - while(std::getline(iss, line)){ - if(line == "enabled"){ + while (std::getline(iss, line)) { + if (line == "enabled") { m_enabled = true; - }else if(line.find("dac") != std::string::npos){ + } else if (line.find("dac") != std::string::npos) { m_dac = line.substr(4); - }else if(line.find("start") != std::string::npos){ + } else if (line.find("start") != std::string::npos) { m_start = std::stoi(line.substr(6)); - }else if(line.find("stop") != std::string::npos){ + } else if (line.find("stop") != std::string::npos) { m_stop = std::stoi(line.substr(5)); - }else if(line.find("step") != std::string::npos){ + } else if (line.find("step") != std::string::npos) { m_step = std::stoi(line.substr(5)); } } @@ -85,14 +82,11 @@ ScanParameters::ScanParameters(const std::string& par){ int ScanParameters::start() const { return m_start; } int ScanParameters::stop() const { return m_stop; } -void ScanParameters::increment_stop(){ - m_stop += 1; -} +void ScanParameters::increment_stop() { m_stop += 1; } int ScanParameters::step() const { return m_step; } const std::string &ScanParameters::dac() const { return m_dac; } bool ScanParameters::enabled() const { return m_enabled; } - RawMasterFile::RawMasterFile(const std::filesystem::path &fpath) : m_fnc(fpath) { if (!std::filesystem::exists(fpath)) { @@ -140,6 +134,10 @@ std::optional RawMasterFile::number_of_rows() const { xy RawMasterFile::geometry() const { return m_geometry; } +size_t RawMasterFile::n_modules() const { + return m_geometry.row * m_geometry.col; +} + std::optional RawMasterFile::quad() const { return m_quad; } // optional values, these may or may not be present in the master file @@ -159,10 +157,8 @@ ScanParameters RawMasterFile::scan_parameters() const { return m_scan_parameters; } - std::optional RawMasterFile::roi() const { return m_roi; } - void RawMasterFile::parse_json(const std::filesystem::path &fpath) { std::ifstream ifs(fpath); json j; @@ -201,17 +197,16 @@ void RawMasterFile::parse_json(const std::filesystem::path &fpath) { // keep the optional empty } - // ---------------------------------------------------------------- // Special treatment of analog flag because of Moench03 - try{ + try { m_analog_flag = j.at("Analog Flag"); - }catch (const json::out_of_range &e) { + } catch (const json::out_of_range &e) { // if it doesn't work still set it to one // to try to decode analog samples (Old Moench03) m_analog_flag = 1; } - try { + try { if (m_analog_flag) { m_analog_samples = j.at("Analog Samples"); } @@ -244,27 +239,27 @@ void RawMasterFile::parse_json(const std::filesystem::path &fpath) { // keep the optional empty } - try{ + try { m_transceiver_flag = j.at("Transceiver Flag"); - if(m_transceiver_flag){ + if (m_transceiver_flag) { m_transceiver_samples = j.at("Transceiver Samples"); } - }catch (const json::out_of_range &e) { + } catch (const json::out_of_range &e) { // keep the optional empty } - try{ + try { std::string scan_parameters = j.at("Scan Parameters"); m_scan_parameters = ScanParameters(scan_parameters); - if(v<7.21){ - m_scan_parameters.increment_stop(); //adjust for endpoint being included - } - }catch (const json::out_of_range &e) { + if (v < 7.21) { + m_scan_parameters + .increment_stop(); // adjust for endpoint being included + } + } catch (const json::out_of_range &e) { // not a scan } - - try{ + try { ROI tmp_roi; auto obj = j.at("Receiver Roi"); tmp_roi.xmin = obj.at("xmin"); @@ -272,37 +267,32 @@ void RawMasterFile::parse_json(const std::filesystem::path &fpath) { tmp_roi.ymin = obj.at("ymin"); tmp_roi.ymax = obj.at("ymax"); - //if any of the values are set update the roi + // if any of the values are set update the roi if (tmp_roi.xmin != 4294967295 || tmp_roi.xmax != 4294967295 || tmp_roi.ymin != 4294967295 || tmp_roi.ymax != 4294967295) { - - if(v<7.21){ + + if (v < 7.21) { tmp_roi.xmax++; tmp_roi.ymax++; } - + m_roi = tmp_roi; } - - }catch (const json::out_of_range &e) { + } catch (const json::out_of_range &e) { // leave the optional empty } - //if we have an roi we need to update the geometry for the subfiles - if (m_roi){ - + // if we have an roi we need to update the geometry for the subfiles + if (m_roi) { } - - - - // Update detector type for Moench - // TODO! How does this work with old .raw master files? - #ifdef AARE_VERBOSE +// Update detector type for Moench +// TODO! How does this work with old .raw master files? +#ifdef AARE_VERBOSE fmt::print("Detecting Moench03: m_pixels_y: {}, m_analog_samples: {}\n", m_pixels_y, m_analog_samples.value_or(0)); - #endif +#endif if (m_type == DetectorType::Moench && !m_analog_samples && m_pixels_y == 400) { m_type = DetectorType::Moench03; @@ -328,19 +318,19 @@ void RawMasterFile::parse_raw(const std::filesystem::path &fpath) { if (key == "Version") { m_version = value; - //TODO!: How old versions can we handle? + // TODO!: How old versions can we handle? auto v = std::stod(value); - //TODO! figure out exactly when we did the change - //This enables padding of f to 12 digits - if (v<4.0) + // TODO! figure out exactly when we did the change + // This enables padding of f to 12 digits + if (v < 4.0) m_fnc.set_old_scheme(true); } else if (key == "TimeStamp") { } else if (key == "Detector Type") { m_type = StringTo(value); - if(m_type==DetectorType::Moench){ + if (m_type == DetectorType::Moench) { m_type = DetectorType::Moench03_old; } } else if (key == "Timing Mode") { @@ -377,10 +367,10 @@ void RawMasterFile::parse_raw(const std::filesystem::path &fpath) { pos = value.find(','); m_pixels_x = std::stoi(value.substr(1, pos)); m_pixels_y = std::stoi(value.substr(pos + 1)); - }else if(key == "row"){ + } else if (key == "row") { pos = value.find('p'); m_pixels_y = std::stoi(value.substr(0, pos)); - }else if(key == "col"){ + } else if (key == "col") { pos = value.find('p'); m_pixels_x = std::stoi(value.substr(0, pos)); } else if (key == "Total Frames") { @@ -391,8 +381,8 @@ void RawMasterFile::parse_raw(const std::filesystem::path &fpath) { m_quad = std::stoi(value); } else if (key == "Max Frames Per File") { m_max_frames_per_file = std::stoi(value); - }else if(key == "Max. Frames Per File"){ - //Version 3.0 way of writing it + } else if (key == "Max. Frames Per File") { + // Version 3.0 way of writing it m_max_frames_per_file = std::stoi(value); } else if (key == "Geometry") { pos = value.find(','); @@ -406,15 +396,14 @@ void RawMasterFile::parse_raw(const std::filesystem::path &fpath) { m_type = DetectorType::Moench03_old; } - - //TODO! Look for d0, d1...dn and update geometry - if(m_geometry.col == 0 && m_geometry.row == 0){ - m_geometry = {1,1}; + // TODO! Look for d0, d1...dn and update geometry + if (m_geometry.col == 0 && m_geometry.row == 0) { + m_geometry = {1, 1}; fmt::print("Warning: No geometry found in master file. Assuming 1x1\n"); } - //TODO! Read files and find actual frames - if(m_frames_in_file==0) + // TODO! Read files and find actual frames + if (m_frames_in_file == 0) m_frames_in_file = m_total_frames_expected; } } // namespace aare diff --git a/src/RawMasterFile.test.cpp b/src/RawMasterFile.test.cpp index 560217f..0e75ec4 100644 --- a/src/RawMasterFile.test.cpp +++ b/src/RawMasterFile.test.cpp @@ -1,12 +1,11 @@ #include "aare/RawMasterFile.hpp" -#include #include "test_config.hpp" +#include using namespace aare; - -TEST_CASE("Parse a master file fname"){ +TEST_CASE("Parse a master file fname") { RawFileNameComponents m("test_master_1.json"); REQUIRE(m.base_name() == "test"); REQUIRE(m.ext() == ".json"); @@ -14,7 +13,7 @@ TEST_CASE("Parse a master file fname"){ REQUIRE(m.base_path() == ""); } -TEST_CASE("Extraction of base path works"){ +TEST_CASE("Extraction of base path works") { RawFileNameComponents m("some/path/test_master_73.json"); REQUIRE(m.base_name() == "test"); REQUIRE(m.ext() == ".json"); @@ -22,7 +21,7 @@ TEST_CASE("Extraction of base path works"){ REQUIRE(m.base_path() == "some/path"); } -TEST_CASE("Construction of master file name and data files"){ +TEST_CASE("Construction of master file name and data files") { RawFileNameComponents m("test_master_1.json"); REQUIRE(m.master_fname() == "test_master_1.json"); REQUIRE(m.data_fname(0, 0) == "test_d0_f0_1.raw"); @@ -31,7 +30,7 @@ TEST_CASE("Construction of master file name and data files"){ REQUIRE(m.data_fname(1, 1) == "test_d1_f1_1.raw"); } -TEST_CASE("Construction of master file name and data files using old scheme"){ +TEST_CASE("Construction of master file name and data files using old scheme") { RawFileNameComponents m("test_master_1.raw"); m.set_old_scheme(true); REQUIRE(m.master_fname() == "test_master_1.raw"); @@ -41,16 +40,15 @@ TEST_CASE("Construction of master file name and data files using old scheme"){ REQUIRE(m.data_fname(1, 1) == "test_d1_f000000000001_1.raw"); } -TEST_CASE("Master file name does not fit pattern"){ +TEST_CASE("Master file name does not fit pattern") { REQUIRE_THROWS(RawFileNameComponents("somefile.json")); REQUIRE_THROWS(RawFileNameComponents("another_test_d0_f0_1.raw")); REQUIRE_THROWS(RawFileNameComponents("test_master_1.txt")); } - - -TEST_CASE("Parse scan parameters"){ - ScanParameters s("[enabled\ndac dac 4\nstart 500\nstop 2200\nstep 5\nsettleTime 100us\n]"); +TEST_CASE("Parse scan parameters") { + ScanParameters s("[enabled\ndac dac 4\nstart 500\nstop 2200\nstep " + "5\nsettleTime 100us\n]"); REQUIRE(s.enabled()); REQUIRE(s.dac() == "dac 4"); REQUIRE(s.start() == 500); @@ -58,7 +56,7 @@ TEST_CASE("Parse scan parameters"){ REQUIRE(s.step() == 5); } -TEST_CASE("A disabled scan"){ +TEST_CASE("A disabled scan") { ScanParameters s("[disabled]"); REQUIRE_FALSE(s.enabled()); REQUIRE(s.dac() == ""); @@ -67,9 +65,9 @@ TEST_CASE("A disabled scan"){ REQUIRE(s.step() == 0); } - -TEST_CASE("Parse a master file in .json format", "[.integration]"){ - auto fpath = test_data_path() / "jungfrau" / "jungfrau_single_master_0.json"; +TEST_CASE("Parse a master file in .json format", "[.integration]") { + auto fpath = + test_data_path() / "jungfrau" / "jungfrau_single_master_0.json"; REQUIRE(std::filesystem::exists(fpath)); RawMasterFile f(fpath); @@ -80,7 +78,7 @@ TEST_CASE("Parse a master file in .json format", "[.integration]"){ REQUIRE(f.detector_type() == DetectorType::Jungfrau); // "Timing Mode": "auto", REQUIRE(f.timing_mode() == TimingMode::Auto); - + // "Geometry": { // "x": 1, // "y": 1 @@ -100,10 +98,9 @@ TEST_CASE("Parse a master file in .json format", "[.integration]"){ // "Max Frames Per File": 3, REQUIRE(f.max_frames_per_file() == 3); - //Jungfrau doesn't write but it is 16 + // Jungfrau doesn't write but it is 16 REQUIRE(f.bitdepth() == 16); - // "Frame Discard Policy": "nodiscard", // "Frame Padding": 1, @@ -125,33 +122,35 @@ TEST_CASE("Parse a master file in .json format", "[.integration]"){ // "Frames in File": 10, REQUIRE(f.frames_in_file() == 10); - //TODO! Should we parse this? - // "Frame Header Format": { - // "Frame Number": "8 bytes", - // "SubFrame Number/ExpLength": "4 bytes", - // "Packet Number": "4 bytes", - // "Bunch ID": "8 bytes", - // "Timestamp": "8 bytes", - // "Module Id": "2 bytes", - // "Row": "2 bytes", - // "Column": "2 bytes", - // "Reserved": "2 bytes", - // "Debug": "4 bytes", - // "Round Robin Number": "2 bytes", - // "Detector Type": "1 byte", - // "Header Version": "1 byte", - // "Packets Caught Mask": "64 bytes" - // } - // } + // TODO! Should we parse this? + // "Frame Header Format": { + // "Frame Number": "8 bytes", + // "SubFrame Number/ExpLength": "4 bytes", + // "Packet Number": "4 bytes", + // "Bunch ID": "8 bytes", + // "Timestamp": "8 bytes", + // "Module Id": "2 bytes", + // "Row": "2 bytes", + // "Column": "2 bytes", + // "Reserved": "2 bytes", + // "Debug": "4 bytes", + // "Round Robin Number": "2 bytes", + // "Detector Type": "1 byte", + // "Header Version": "1 byte", + // "Packets Caught Mask": "64 bytes" + // } + // } REQUIRE_FALSE(f.analog_samples()); - REQUIRE_FALSE(f.digital_samples()); - + REQUIRE_FALSE(f.digital_samples()); } -TEST_CASE("Parse a master file in .raw format", "[.integration]"){ - - auto fpath = test_data_path() / "moench/moench04_noise_200V_sto_both_100us_no_light_thresh_900_master_0.raw"; +TEST_CASE("Parse a master file in .raw format", "[.integration]") { + + auto fpath = + test_data_path() / + "moench/" + "moench04_noise_200V_sto_both_100us_no_light_thresh_900_master_0.raw"; REQUIRE(std::filesystem::exists(fpath)); RawMasterFile f(fpath); @@ -209,80 +208,74 @@ TEST_CASE("Parse a master file in .raw format", "[.integration]"){ // Detector Type : 1 byte // Header Version : 1 byte // Packets Caught Mask : 64 bytes - - } - -TEST_CASE("Read eiger master file", "[.integration]"){ -auto fpath = test_data_path() / "eiger" / "eiger_500k_32bit_master_0.json"; +TEST_CASE("Read eiger master file", "[.integration]") { + auto fpath = test_data_path() / "eiger" / "eiger_500k_32bit_master_0.json"; REQUIRE(std::filesystem::exists(fpath)); RawMasterFile f(fpath); - -// { -// "Version": 7.2, -REQUIRE(f.version() == "7.2"); -// "Timestamp": "Tue Mar 26 17:24:34 2024", -// "Detector Type": "Eiger", -REQUIRE(f.detector_type() == DetectorType::Eiger); -// "Timing Mode": "auto", -REQUIRE(f.timing_mode() == TimingMode::Auto); -// "Geometry": { -// "x": 2, -// "y": 2 -// }, -// "Image Size in bytes": 524288, -REQUIRE(f.image_size_in_bytes() == 524288); -// "Pixels": { -// "x": 512, -REQUIRE(f.pixels_x() == 512); -// "y": 256 -REQUIRE(f.pixels_y() == 256); -// }, -// "Max Frames Per File": 10000, -REQUIRE(f.max_frames_per_file() == 10000); -// "Frame Discard Policy": "nodiscard", -REQUIRE(f.frame_discard_policy() == FrameDiscardPolicy::NoDiscard); -// "Frame Padding": 1, -REQUIRE(f.frame_padding() == 1); - -// "Scan Parameters": "[disabled]", -// "Total Frames": 3, -// "Receiver Roi": { -// "xmin": 4294967295, -// "xmax": 4294967295, -// "ymin": 4294967295, -// "ymax": 4294967295 -// }, -// "Dynamic Range": 32, -// "Ten Giga": 0, -// "Exptime": "5s", -// "Period": "1s", -// "Threshold Energy": -1, -// "Sub Exptime": "2.62144ms", -// "Sub Period": "2.62144ms", -// "Quad": 0, -// "Number of rows": 256, -// "Rate Corrections": "[0, 0]", -// "Frames in File": 3, -// "Frame Header Format": { -// "Frame Number": "8 bytes", -// "SubFrame Number/ExpLength": "4 bytes", -// "Packet Number": "4 bytes", -// "Bunch ID": "8 bytes", -// "Timestamp": "8 bytes", -// "Module Id": "2 bytes", -// "Row": "2 bytes", -// "Column": "2 bytes", -// "Reserved": "2 bytes", -// "Debug": "4 bytes", -// "Round Robin Number": "2 bytes", -// "Detector Type": "1 byte", -// "Header Version": "1 byte", -// "Packets Caught Mask": "64 bytes" -// } -// } - + // { + // "Version": 7.2, + REQUIRE(f.version() == "7.2"); + // "Timestamp": "Tue Mar 26 17:24:34 2024", + // "Detector Type": "Eiger", + REQUIRE(f.detector_type() == DetectorType::Eiger); + // "Timing Mode": "auto", + REQUIRE(f.timing_mode() == TimingMode::Auto); + // "Geometry": { + // "x": 2, + // "y": 2 + // }, + // "Image Size in bytes": 524288, + REQUIRE(f.image_size_in_bytes() == 524288); + // "Pixels": { + // "x": 512, + REQUIRE(f.pixels_x() == 512); + // "y": 256 + REQUIRE(f.pixels_y() == 256); + // }, + // "Max Frames Per File": 10000, + REQUIRE(f.max_frames_per_file() == 10000); + // "Frame Discard Policy": "nodiscard", + REQUIRE(f.frame_discard_policy() == FrameDiscardPolicy::NoDiscard); + // "Frame Padding": 1, + REQUIRE(f.frame_padding() == 1); + // "Scan Parameters": "[disabled]", + // "Total Frames": 3, + // "Receiver Roi": { + // "xmin": 4294967295, + // "xmax": 4294967295, + // "ymin": 4294967295, + // "ymax": 4294967295 + // }, + // "Dynamic Range": 32, + // "Ten Giga": 0, + // "Exptime": "5s", + // "Period": "1s", + // "Threshold Energy": -1, + // "Sub Exptime": "2.62144ms", + // "Sub Period": "2.62144ms", + // "Quad": 0, + // "Number of rows": 256, + // "Rate Corrections": "[0, 0]", + // "Frames in File": 3, + // "Frame Header Format": { + // "Frame Number": "8 bytes", + // "SubFrame Number/ExpLength": "4 bytes", + // "Packet Number": "4 bytes", + // "Bunch ID": "8 bytes", + // "Timestamp": "8 bytes", + // "Module Id": "2 bytes", + // "Row": "2 bytes", + // "Column": "2 bytes", + // "Reserved": "2 bytes", + // "Debug": "4 bytes", + // "Round Robin Number": "2 bytes", + // "Detector Type": "1 byte", + // "Header Version": "1 byte", + // "Packets Caught Mask": "64 bytes" + // } + // } } \ No newline at end of file diff --git a/src/RawSubFile.cpp b/src/RawSubFile.cpp index 9e7a421..3ed2c6f 100644 --- a/src/RawSubFile.cpp +++ b/src/RawSubFile.cpp @@ -1,69 +1,70 @@ #include "aare/RawSubFile.hpp" #include "aare/PixelMap.hpp" +#include "aare/algorithm.hpp" +#include "aare/logger.hpp" #include "aare/utils/ifstream_helpers.hpp" + #include // memcpy #include #include - - +#include namespace aare { RawSubFile::RawSubFile(const std::filesystem::path &fname, DetectorType detector, size_t rows, size_t cols, size_t bitdepth, uint32_t pos_row, uint32_t pos_col) - : m_detector_type(detector), m_bitdepth(bitdepth), m_fname(fname), - m_rows(rows), m_cols(cols), - m_bytes_per_frame((m_bitdepth / 8) * m_rows * m_cols), m_pos_row(pos_row), - m_pos_col(pos_col) { + : m_detector_type(detector), m_bitdepth(bitdepth), m_rows(rows), + m_cols(cols), m_bytes_per_frame((m_bitdepth / 8) * m_rows * m_cols), + m_pos_row(pos_row), m_pos_col(pos_col) { + + LOG(logDEBUG) << "RawSubFile::RawSubFile()"; if (m_detector_type == DetectorType::Moench03_old) { m_pixel_map = GenerateMoench03PixelMap(); } else if (m_detector_type == DetectorType::Eiger && m_pos_row % 2 == 0) { m_pixel_map = GenerateEigerFlipRowsPixelMap(); } - if (std::filesystem::exists(fname)) { - m_num_frames = std::filesystem::file_size(fname) / - (sizeof(DetectorHeader) + rows * cols * bitdepth / 8); - } else { - throw std::runtime_error( - LOCATION + fmt::format("File {} does not exist", m_fname.string())); - } - - // fp = fopen(m_fname.string().c_str(), "rb"); - m_file.open(m_fname, std::ios::binary); - if (!m_file.is_open()) { - throw std::runtime_error( - LOCATION + fmt::format("Could not open file {}", m_fname.string())); - } - -#ifdef AARE_VERBOSE - fmt::print("Opened file: {} with {} frames\n", m_fname.string(), m_num_frames); - fmt::print("m_rows: {}, m_cols: {}, m_bitdepth: {}\n", m_rows, m_cols, - m_bitdepth); - fmt::print("file size: {}\n", std::filesystem::file_size(fname)); -#endif + parse_fname(fname); + scan_files(); + open_file(m_current_file_index); // open the first file } void RawSubFile::seek(size_t frame_index) { - if (frame_index >= m_num_frames) { - throw std::runtime_error(LOCATION + fmt::format("Frame index {} out of range in a file with {} frames", frame_index, m_num_frames)); + LOG(logDEBUG) << "RawSubFile::seek(" << frame_index << ")"; + if (frame_index >= m_total_frames) { + throw std::runtime_error(LOCATION + " Frame index out of range: " + + std::to_string(frame_index)); } - m_file.seekg((sizeof(DetectorHeader) + bytes_per_frame()) * frame_index); + m_current_frame_index = frame_index; + auto file_index = first_larger(m_last_frame_in_file, frame_index); + + if (file_index != m_current_file_index) + open_file(file_index); + + auto frame_offset = (file_index) + ? frame_index - m_last_frame_in_file[file_index - 1] + : frame_index; + auto byte_offset = + frame_offset * (m_bytes_per_frame + sizeof(DetectorHeader)); + m_file.seekg(byte_offset); } size_t RawSubFile::tell() { - return m_file.tellg() / (sizeof(DetectorHeader) + bytes_per_frame()); + LOG(logDEBUG) << "RawSubFile::tell():" << m_current_frame_index; + return m_current_frame_index; } void RawSubFile::read_into(std::byte *image_buf, DetectorHeader *header) { + LOG(logDEBUG) << "RawSubFile::read_into()"; + if (header) { m_file.read(reinterpret_cast(header), sizeof(DetectorHeader)); } else { m_file.seekg(sizeof(DetectorHeader), std::ios::cur); } - if (m_file.fail()){ + if (m_file.fail()) { throw std::runtime_error(LOCATION + ifstream_error_msg(m_file)); } @@ -72,14 +73,15 @@ void RawSubFile::read_into(std::byte *image_buf, DetectorHeader *header) { // read into a temporary buffer and then copy the data to the buffer // in the correct order // TODO! add 4 bit support - if(m_bitdepth == 8){ + if (m_bitdepth == 8) { read_with_map(image_buf); - }else if (m_bitdepth == 16) { + } else if (m_bitdepth == 16) { read_with_map(image_buf); } else if (m_bitdepth == 32) { read_with_map(image_buf); - }else{ - throw std::runtime_error("Unsupported bitdepth for read with pixel map"); + } else { + throw std::runtime_error( + "Unsupported bitdepth for read with pixel map"); } } else { @@ -87,12 +89,20 @@ void RawSubFile::read_into(std::byte *image_buf, DetectorHeader *header) { m_file.read(reinterpret_cast(image_buf), bytes_per_frame()); } - if (m_file.fail()){ + if (m_file.fail()) { throw std::runtime_error(LOCATION + ifstream_error_msg(m_file)); } + + ++m_current_frame_index; + if (m_current_frame_index >= m_last_frame_in_file[m_current_file_index] && + (m_current_frame_index < m_total_frames)) { + ++m_current_file_index; + open_file(m_current_file_index); + } } -void RawSubFile::read_into(std::byte *image_buf, size_t n_frames, DetectorHeader *header) { +void RawSubFile::read_into(std::byte *image_buf, size_t n_frames, + DetectorHeader *header) { for (size_t i = 0; i < n_frames; i++) { read_into(image_buf, header); image_buf += bytes_per_frame(); @@ -102,10 +112,7 @@ void RawSubFile::read_into(std::byte *image_buf, size_t n_frames, DetectorHeader } } - - -template -void RawSubFile::read_with_map(std::byte *image_buf) { +template void RawSubFile::read_with_map(std::byte *image_buf) { auto part_buffer = new std::byte[bytes_per_frame()]; m_file.read(reinterpret_cast(part_buffer), bytes_per_frame()); auto *data = reinterpret_cast(image_buf); @@ -130,4 +137,74 @@ size_t RawSubFile::frame_number(size_t frame_index) { return h.frameNumber; } +void RawSubFile::parse_fname(const std::filesystem::path &fname) { + LOG(logDEBUG) << "RawSubFile::parse_fname()"; + // data has the format: /path/too/data/jungfrau_single_d0_f1_0.raw + // d0 is the module index, will not change for this file + // f1 is the file index - thi is the one we need + // 0 is the measurement index, will not change + m_path = fname.parent_path(); + m_base_name = fname.filename(); + + // Regex to extract numbers after 'd' and 'f' + std::regex pattern(R"(^(.*_d)(\d+)(_f)(\d+)(_\d+\.raw)$)"); + std::smatch match; + + if (std::regex_match(m_base_name, match, pattern)) { + m_offset = std::stoi(match[4].str()); // find the first file index in + // case of a truncated series + m_base_name = match[1].str() + match[2].str() + match[3].str() + "{}" + + match[5].str(); + LOG(logDEBUG) << "Base name: " << m_base_name; + LOG(logDEBUG) << "Offset: " << m_offset; + LOG(logDEBUG) << "Path: " << m_path.string(); + } else { + throw std::runtime_error( + LOCATION + + fmt::format("Could not parse file name {}", fname.string())); + } +} + +std::filesystem::path RawSubFile::fpath(size_t file_index) const { + auto fname = fmt::format(m_base_name, file_index); + return m_path / fname; +} + +void RawSubFile::open_file(size_t file_index) { + m_file.close(); + auto fname = fpath(file_index + m_offset); + LOG(logDEBUG) << "RawSubFile::open_file(): " << fname.string(); + m_file.open(fname, std::ios::binary); + if (!m_file.is_open()) { + throw std::runtime_error( + LOCATION + + fmt::format("Could not open file {}", fpath(file_index).string())); + } + m_current_file_index = file_index; +} + +void RawSubFile::scan_files() { + LOG(logDEBUG) << "RawSubFile::scan_files()"; + // find how many files we have and the number of frames in each file + m_last_frame_in_file.clear(); + size_t file_index = m_offset; + + while (std::filesystem::exists(fpath(file_index))) { + auto n_frames = std::filesystem::file_size(fpath(file_index)) / + (m_bytes_per_frame + sizeof(DetectorHeader)); + m_last_frame_in_file.push_back(n_frames); + LOG(logDEBUG) << "Found: " << n_frames + << " frames in file: " << fpath(file_index).string(); + ++file_index; + } + + // find where we need to open the next file and total number of frames + m_last_frame_in_file = cumsum(m_last_frame_in_file); + if (m_last_frame_in_file.empty()) { + m_total_frames = 0; + } else { + m_total_frames = m_last_frame_in_file.back(); + } +} + } // namespace aare \ No newline at end of file diff --git a/src/RawSubFile.test.cpp b/src/RawSubFile.test.cpp new file mode 100644 index 0000000..37d071b --- /dev/null +++ b/src/RawSubFile.test.cpp @@ -0,0 +1,79 @@ +#include "aare/RawSubFile.hpp" +#include "aare/File.hpp" +#include "aare/NDArray.hpp" +#include "test_config.hpp" +#include + +using namespace aare; + +TEST_CASE("Read frames directly from a RawSubFile", "[.files]") { + auto fpath_raw = + test_data_path() / "raw/jungfrau" / "jungfrau_single_d0_f0_0.raw"; + REQUIRE(std::filesystem::exists(fpath_raw)); + + RawSubFile f(fpath_raw, DetectorType::Jungfrau, 512, 1024, 16); + REQUIRE(f.rows() == 512); + REQUIRE(f.cols() == 1024); + REQUIRE(f.pixels_per_frame() == 512 * 1024); + REQUIRE(f.bytes_per_frame() == 512 * 1024 * 2); + REQUIRE(f.bytes_per_pixel() == 2); + + auto fpath_npy = + test_data_path() / "raw/jungfrau" / "jungfrau_single_0.npy"; + REQUIRE(std::filesystem::exists(fpath_npy)); + + // Numpy file with the same data to use as reference + File npy(fpath_npy, "r"); + + CHECK(f.frames_in_file() == 10); + CHECK(npy.total_frames() == 10); + + DetectorHeader header{}; + NDArray image( + {static_cast(f.rows()), static_cast(f.cols())}); + for (size_t i = 0; i < 10; ++i) { + CHECK(f.tell() == i); + f.read_into(image.buffer(), &header); + auto npy_frame = npy.read_frame(); + CHECK((image.view() == npy_frame.view())); + } +} + +TEST_CASE("Read frames directly from a RawSubFile starting at the second file", + "[.files]") { + // we know this file has 10 frames with frame numbers 1 to 10 + // f0 1,2,3 + // f1 4,5,6 <-- starting here + // f2 7,8,9 + // f3 10 + + auto fpath_raw = + test_data_path() / "raw/jungfrau" / "jungfrau_single_d0_f1_0.raw"; + REQUIRE(std::filesystem::exists(fpath_raw)); + + RawSubFile f(fpath_raw, DetectorType::Jungfrau, 512, 1024, 16); + + auto fpath_npy = + test_data_path() / "raw/jungfrau" / "jungfrau_single_0.npy"; + REQUIRE(std::filesystem::exists(fpath_npy)); + + // Numpy file with the same data to use as reference + File npy(fpath_npy, "r"); + npy.seek(3); + + CHECK(f.frames_in_file() == 7); + CHECK(npy.total_frames() == 10); + + DetectorHeader header{}; + NDArray image( + {static_cast(f.rows()), static_cast(f.cols())}); + for (size_t i = 0; i < 7; ++i) { + CHECK(f.tell() == i); + f.read_into(image.buffer(), &header); + // frame numbers start at 1 frame index at 0 + // adding 3 + 1 to verify the frame number + CHECK(header.frameNumber == i + 4); + auto npy_frame = npy.read_frame(); + CHECK((image.view() == npy_frame.view())); + } +} \ No newline at end of file diff --git a/src/algorithm.test.cpp b/src/algorithm.test.cpp index 5452fcf..3706aba 100644 --- a/src/algorithm.test.cpp +++ b/src/algorithm.test.cpp @@ -160,3 +160,35 @@ TEST_CASE("cumsum works with negative numbers", "[algorithm]") { REQUIRE(result[3] == -6); REQUIRE(result[4] == -10); } + +TEST_CASE("cumsum on an empty vector", "[algorithm]") { + std::vector vec = {}; + auto result = aare::cumsum(vec); + REQUIRE(result.size() == 0); +} + +TEST_CASE("All equal on an empty vector is false", "[algorithm]") { + std::vector vec = {}; + REQUIRE(aare::all_equal(vec) == false); +} + +TEST_CASE("All equal on a vector with 1 element is true", "[algorithm]") { + std::vector vec = {1}; + REQUIRE(aare::all_equal(vec) == true); +} + +TEST_CASE("All equal on a vector with 2 elements is true", "[algorithm]") { + std::vector vec = {1, 1}; + REQUIRE(aare::all_equal(vec) == true); +} + +TEST_CASE("All equal on a vector with two different elements is false", + "[algorithm]") { + std::vector vec = {1, 2}; + REQUIRE(aare::all_equal(vec) == false); +} + +TEST_CASE("Last element is different", "[algorithm]") { + std::vector vec = {1, 1, 1, 1, 2}; + REQUIRE(aare::all_equal(vec) == false); +} diff --git a/src/decode.cpp b/src/decode.cpp index 436ad7b..0c78b70 100644 --- a/src/decode.cpp +++ b/src/decode.cpp @@ -2,9 +2,9 @@ #include namespace aare { -uint16_t adc_sar_05_decode64to16(uint64_t input){ +uint16_t adc_sar_05_decode64to16(uint64_t input) { - //we want bits 29,19,28,18,31,21,27,20,24,23,25,22 and then pad to 16 + // we want bits 29,19,28,18,31,21,27,20,24,23,25,22 and then pad to 16 uint16_t output = 0; output |= ((input >> 22) & 1) << 11; output |= ((input >> 25) & 1) << 10; @@ -21,19 +21,21 @@ uint16_t adc_sar_05_decode64to16(uint64_t input){ return output; } -void adc_sar_05_decode64to16(NDView input, NDView output){ - if(input.shape() != output.shape()){ - throw std::invalid_argument(LOCATION + " input and output shapes must match"); +void adc_sar_05_decode64to16(NDView input, + NDView output) { + if (input.shape() != output.shape()) { + throw std::invalid_argument(LOCATION + + " input and output shapes must match"); } - for(ssize_t i = 0; i < input.shape(0); i++){ - for(ssize_t j = 0; j < input.shape(1); j++){ - output(i,j) = adc_sar_05_decode64to16(input(i,j)); + for (ssize_t i = 0; i < input.shape(0); i++) { + for (ssize_t j = 0; j < input.shape(1); j++) { + output(i, j) = adc_sar_05_decode64to16(input(i, j)); } } } -uint16_t adc_sar_04_decode64to16(uint64_t input){ +uint16_t adc_sar_04_decode64to16(uint64_t input) { // bit_map = array([15,17,19,21,23,4,6,8,10,12,14,16] LSB->MSB uint16_t output = 0; @@ -52,20 +54,23 @@ uint16_t adc_sar_04_decode64to16(uint64_t input){ return output; } -void adc_sar_04_decode64to16(NDView input, NDView output){ - if(input.shape() != output.shape()){ - throw std::invalid_argument(LOCATION + " input and output shapes must match"); +void adc_sar_04_decode64to16(NDView input, + NDView output) { + if (input.shape() != output.shape()) { + throw std::invalid_argument(LOCATION + + " input and output shapes must match"); } - for(ssize_t i = 0; i < input.shape(0); i++){ - for(ssize_t j = 0; j < input.shape(1); j++){ - output(i,j) = adc_sar_04_decode64to16(input(i,j)); + for (ssize_t i = 0; i < input.shape(0); i++) { + for (ssize_t j = 0; j < input.shape(1); j++) { + output(i, j) = adc_sar_04_decode64to16(input(i, j)); } } } double apply_custom_weights(uint16_t input, const NDView weights) { - if(weights.size() > 16){ - throw std::invalid_argument("weights size must be less than or equal to 16"); + if (weights.size() > 16) { + throw std::invalid_argument( + "weights size must be less than or equal to 16"); } double result = 0.0; @@ -73,30 +78,30 @@ double apply_custom_weights(uint16_t input, const NDView weights) { result += ((input >> i) & 1) * std::pow(weights[i], i); } return result; - } -void apply_custom_weights(NDView input, NDView output, const NDView weights) { - if(input.shape() != output.shape()){ - throw std::invalid_argument(LOCATION + " input and output shapes must match"); +void apply_custom_weights(NDView input, NDView output, + const NDView weights) { + if (input.shape() != output.shape()) { + throw std::invalid_argument(LOCATION + + " input and output shapes must match"); } - //Calculate weights to avoid repeatedly calling std::pow + // Calculate weights to avoid repeatedly calling std::pow std::vector weights_powers(weights.size()); for (ssize_t i = 0; i < weights.size(); ++i) { weights_powers[i] = std::pow(weights[i], i); } // Apply custom weights to each element in the input array - for (ssize_t i = 0; i < input.shape(0); i++) { + for (ssize_t i = 0; i < input.shape(0); i++) { double result = 0.0; - for (size_t bit_index = 0; bit_index < weights_powers.size(); ++bit_index) { + for (size_t bit_index = 0; bit_index < weights_powers.size(); + ++bit_index) { result += ((input(i) >> bit_index) & 1) * weights_powers[bit_index]; } output(i) = result; } } - - } // namespace aare diff --git a/src/decode.test.cpp b/src/decode.test.cpp index 1e4b2fc..b4310ca 100644 --- a/src/decode.test.cpp +++ b/src/decode.test.cpp @@ -1,17 +1,16 @@ #include "aare/decode.hpp" -#include -#include #include "aare/NDArray.hpp" +#include +#include using Catch::Matchers::WithinAbs; #include -TEST_CASE("test_adc_sar_05_decode64to16"){ +TEST_CASE("test_adc_sar_05_decode64to16") { uint64_t input = 0; uint16_t output = aare::adc_sar_05_decode64to16(input); CHECK(output == 0); - // bit 29 on th input is bit 0 on the output input = 1UL << 29; output = aare::adc_sar_05_decode64to16(input); @@ -25,7 +24,6 @@ TEST_CASE("test_adc_sar_05_decode64to16"){ CHECK(output == (1 << i)); } - // test a few "random" values input = 0; input |= (1UL << 29); @@ -34,7 +32,6 @@ TEST_CASE("test_adc_sar_05_decode64to16"){ output = aare::adc_sar_05_decode64to16(input); CHECK(output == 7UL); - input = 0; input |= (1UL << 18); input |= (1UL << 27); @@ -47,10 +44,9 @@ TEST_CASE("test_adc_sar_05_decode64to16"){ input |= (1UL << 22); output = aare::adc_sar_05_decode64to16(input); CHECK(output == 3072UL); - } +} - - TEST_CASE("test_apply_custom_weights") { +TEST_CASE("test_apply_custom_weights") { uint16_t input = 1; aare::NDArray weights_data({3}, 0.0); @@ -60,7 +56,6 @@ TEST_CASE("test_adc_sar_05_decode64to16"){ auto weights = weights_data.view(); - double output = aare::apply_custom_weights(input, weights); CHECK_THAT(output, WithinAbs(1.0, 0.001)); @@ -68,7 +63,6 @@ TEST_CASE("test_adc_sar_05_decode64to16"){ output = aare::apply_custom_weights(input, weights); CHECK_THAT(output, WithinAbs(2.1, 0.001)); - input = 1 << 2; output = aare::apply_custom_weights(input, weights); CHECK_THAT(output, WithinAbs(3.24, 0.001)); @@ -76,5 +70,4 @@ TEST_CASE("test_adc_sar_05_decode64to16"){ input = 0b111; output = aare::apply_custom_weights(input, weights); CHECK_THAT(output, WithinAbs(6.34, 0.001)); - - } \ No newline at end of file +} \ No newline at end of file diff --git a/src/defs.cpp b/src/defs.cpp index 7c7cc45..b93ff06 100644 --- a/src/defs.cpp +++ b/src/defs.cpp @@ -5,15 +5,11 @@ #include namespace aare { - -void assert_failed(const std::string &msg) - { +void assert_failed(const std::string &msg) { fmt::print(msg); exit(1); } - - /** * @brief Convert a DetectorType to a string * @param type DetectorType @@ -40,7 +36,7 @@ template <> std::string ToString(DetectorType arg) { case DetectorType::Xilinx_ChipTestBoard: return "Xilinx_ChipTestBoard"; - //Custom ones + // Custom ones case DetectorType::Moench03: return "Moench03"; case DetectorType::Moench03_old: @@ -48,8 +44,8 @@ template <> std::string ToString(DetectorType arg) { case DetectorType::Unknown: return "Unknown"; - //no default case to trigger compiler warning if not all - //enum values are handled + // no default case to trigger compiler warning if not all + // enum values are handled } throw std::runtime_error("Could not decode detector to string"); } @@ -80,14 +76,14 @@ template <> DetectorType StringTo(const std::string &arg) { if (arg == "Xilinx_ChipTestBoard") return DetectorType::Xilinx_ChipTestBoard; - //Custom ones + // Custom ones if (arg == "Moench03") return DetectorType::Moench03; if (arg == "Moench03_old") return DetectorType::Moench03_old; if (arg == "Unknown") return DetectorType::Unknown; - + throw std::runtime_error("Could not decode detector from: \"" + arg + "\""); } @@ -102,7 +98,8 @@ template <> TimingMode StringTo(const std::string &arg) { return TimingMode::Auto; if (arg == "trigger") return TimingMode::Trigger; - throw std::runtime_error("Could not decode timing mode from: \"" + arg + "\""); + throw std::runtime_error("Could not decode timing mode from: \"" + arg + + "\""); } template <> FrameDiscardPolicy StringTo(const std::string &arg) { @@ -112,7 +109,8 @@ template <> FrameDiscardPolicy StringTo(const std::string &arg) { return FrameDiscardPolicy::Discard; if (arg == "discardpartial") return FrameDiscardPolicy::DiscardPartial; - throw std::runtime_error("Could not decode frame discard policy from: \"" + arg + "\""); + throw std::runtime_error("Could not decode frame discard policy from: \"" + + arg + "\""); } // template <> TimingMode StringTo(std::string mode); diff --git a/src/defs.test.cpp b/src/defs.test.cpp index 6ab9394..2106d86 100644 --- a/src/defs.test.cpp +++ b/src/defs.test.cpp @@ -3,12 +3,12 @@ #include #include -using aare::ToString; using aare::StringTo; +using aare::ToString; TEST_CASE("Enum to string conversion") { - // TODO! By the way I don't think the enum string conversions should be in the defs.hpp file - // but let's use this to show a test + // TODO! By the way I don't think the enum string conversions should be in + // the defs.hpp file but let's use this to show a test REQUIRE(ToString(aare::DetectorType::Generic) == "Generic"); REQUIRE(ToString(aare::DetectorType::Eiger) == "Eiger"); REQUIRE(ToString(aare::DetectorType::Gotthard) == "Gotthard"); @@ -17,30 +17,42 @@ TEST_CASE("Enum to string conversion") { REQUIRE(ToString(aare::DetectorType::Moench) == "Moench"); REQUIRE(ToString(aare::DetectorType::Mythen3) == "Mythen3"); REQUIRE(ToString(aare::DetectorType::Gotthard2) == "Gotthard2"); - REQUIRE(ToString(aare::DetectorType::Xilinx_ChipTestBoard) == "Xilinx_ChipTestBoard"); + REQUIRE(ToString(aare::DetectorType::Xilinx_ChipTestBoard) == + "Xilinx_ChipTestBoard"); REQUIRE(ToString(aare::DetectorType::Moench03) == "Moench03"); REQUIRE(ToString(aare::DetectorType::Moench03_old) == "Moench03_old"); REQUIRE(ToString(aare::DetectorType::Unknown) == "Unknown"); } -TEST_CASE("String to enum"){ - REQUIRE(StringTo("Generic") == aare::DetectorType::Generic); +TEST_CASE("String to enum") { + REQUIRE(StringTo("Generic") == + aare::DetectorType::Generic); REQUIRE(StringTo("Eiger") == aare::DetectorType::Eiger); - REQUIRE(StringTo("Gotthard") == aare::DetectorType::Gotthard); - REQUIRE(StringTo("Jungfrau") == aare::DetectorType::Jungfrau); - REQUIRE(StringTo("ChipTestBoard") == aare::DetectorType::ChipTestBoard); - REQUIRE(StringTo("Moench") == aare::DetectorType::Moench); - REQUIRE(StringTo("Mythen3") == aare::DetectorType::Mythen3); - REQUIRE(StringTo("Gotthard2") == aare::DetectorType::Gotthard2); - REQUIRE(StringTo("Xilinx_ChipTestBoard") == aare::DetectorType::Xilinx_ChipTestBoard); - REQUIRE(StringTo("Moench03") == aare::DetectorType::Moench03); - REQUIRE(StringTo("Moench03_old") == aare::DetectorType::Moench03_old); - REQUIRE(StringTo("Unknown") == aare::DetectorType::Unknown); + REQUIRE(StringTo("Gotthard") == + aare::DetectorType::Gotthard); + REQUIRE(StringTo("Jungfrau") == + aare::DetectorType::Jungfrau); + REQUIRE(StringTo("ChipTestBoard") == + aare::DetectorType::ChipTestBoard); + REQUIRE(StringTo("Moench") == + aare::DetectorType::Moench); + REQUIRE(StringTo("Mythen3") == + aare::DetectorType::Mythen3); + REQUIRE(StringTo("Gotthard2") == + aare::DetectorType::Gotthard2); + REQUIRE(StringTo("Xilinx_ChipTestBoard") == + aare::DetectorType::Xilinx_ChipTestBoard); + REQUIRE(StringTo("Moench03") == + aare::DetectorType::Moench03); + REQUIRE(StringTo("Moench03_old") == + aare::DetectorType::Moench03_old); + REQUIRE(StringTo("Unknown") == + aare::DetectorType::Unknown); } -TEST_CASE("Enum values"){ - //Since some of the enums are written to file we need to make sure - //they match the value in the slsDetectorPackage +TEST_CASE("Enum values") { + // Since some of the enums are written to file we need to make sure + // they match the value in the slsDetectorPackage REQUIRE(static_cast(aare::DetectorType::Generic) == 0); REQUIRE(static_cast(aare::DetectorType::Eiger) == 1); @@ -52,7 +64,7 @@ TEST_CASE("Enum values"){ REQUIRE(static_cast(aare::DetectorType::Gotthard2) == 7); REQUIRE(static_cast(aare::DetectorType::Xilinx_ChipTestBoard) == 8); - //Not included + // Not included REQUIRE(static_cast(aare::DetectorType::Moench03) == 100); } @@ -85,5 +97,6 @@ TEST_CASE("DynamicCluster creation") { // double v3 = c2.get(33 * 44 - 1); // REQUIRE(aare::compare_floats(123.11, v3)); -// REQUIRE_THROWS_AS(c2.set(0, 1), std::invalid_argument); // set int to double +// REQUIRE_THROWS_AS(c2.set(0, 1), std::invalid_argument); // set int to +// double // } \ No newline at end of file diff --git a/src/geo_helpers.cpp b/src/geo_helpers.cpp index 39086ec..96a9056 100644 --- a/src/geo_helpers.cpp +++ b/src/geo_helpers.cpp @@ -2,15 +2,16 @@ #include "aare/geo_helpers.hpp" #include "fmt/core.h" -namespace aare{ +namespace aare { DetectorGeometry update_geometry_with_roi(DetectorGeometry geo, aare::ROI roi) { - #ifdef AARE_VERBOSE +#ifdef AARE_VERBOSE fmt::println("update_geometry_with_roi() called with ROI: {} {} {} {}", roi.xmin, roi.xmax, roi.ymin, roi.ymax); - fmt::println("Geometry: {} {} {} {} {} {}", - geo.modules_x, geo.modules_y, geo.pixels_x, geo.pixels_y, geo.module_gap_row, geo.module_gap_col); - #endif + fmt::println("Geometry: {} {} {} {} {} {}", geo.modules_x, geo.modules_y, + geo.pixels_x, geo.pixels_y, geo.module_gap_row, + geo.module_gap_col); +#endif int pos_y = 0; int pos_y_increment = 0; for (int row = 0; row < geo.modules_y; row++) { @@ -41,9 +42,9 @@ DetectorGeometry update_geometry_with_roi(DetectorGeometry geo, aare::ROI roi) { if (m.origin_y + m.height < roi.ymin) { m.height = 0; } else { - if ((roi.ymin > m.origin_y) && (roi.ymin < m.origin_y + m.height)) { + if ((roi.ymin > m.origin_y) && + (roi.ymin < m.origin_y + m.height)) { m.height -= roi.ymin - m.origin_y; - } if (roi.ymax < m.origin_y + original_height) { m.height -= m.origin_y + original_height - roi.ymax; @@ -51,9 +52,10 @@ DetectorGeometry update_geometry_with_roi(DetectorGeometry geo, aare::ROI roi) { m.origin_y = pos_y; pos_y_increment = m.height; } - #ifdef AARE_VERBOSE - fmt::println("Module {} {} {} {}", m.origin_x, m.origin_y, m.width, m.height); - #endif +#ifdef AARE_VERBOSE + fmt::println("Module {} {} {} {}", m.origin_x, m.origin_y, m.width, + m.height); +#endif } // increment pos_y pos_y += pos_y_increment; @@ -65,7 +67,6 @@ DetectorGeometry update_geometry_with_roi(DetectorGeometry geo, aare::ROI roi) { geo.pixels_y = roi.height(); return geo; - } } // namespace aare \ No newline at end of file diff --git a/src/geo_helpers.test.cpp b/src/geo_helpers.test.cpp index 08ee96c..48ae9cf 100644 --- a/src/geo_helpers.test.cpp +++ b/src/geo_helpers.test.cpp @@ -1,6 +1,6 @@ #include "aare/File.hpp" -#include "aare/RawMasterFile.hpp" //needed for ROI #include "aare/RawFile.hpp" +#include "aare/RawMasterFile.hpp" //needed for ROI #include #include @@ -8,26 +8,24 @@ #include "aare/geo_helpers.hpp" #include "test_config.hpp" -TEST_CASE("Simple ROIs on one module"){ - // DetectorGeometry update_geometry_with_roi(DetectorGeometry geo, aare::ROI roi) +TEST_CASE("Simple ROIs on one module") { + // DetectorGeometry update_geometry_with_roi(DetectorGeometry geo, aare::ROI + // roi) aare::DetectorGeometry geo; - - aare::ModuleGeometry mod; + aare::ModuleGeometry mod; mod.origin_x = 0; mod.origin_y = 0; mod.width = 1024; mod.height = 512; - - geo.pixels_x = 1024; geo.pixels_y = 512; geo.modules_x = 1; geo.modules_y = 1; geo.module_pixel_0.push_back(mod); - SECTION("ROI is the whole module"){ + SECTION("ROI is the whole module") { aare::ROI roi; roi.xmin = 0; roi.xmax = 1024; @@ -42,7 +40,7 @@ TEST_CASE("Simple ROIs on one module"){ REQUIRE(updated_geo.module_pixel_0[0].height == 512); REQUIRE(updated_geo.module_pixel_0[0].width == 1024); } - SECTION("ROI is the top left corner of the module"){ + SECTION("ROI is the top left corner of the module") { aare::ROI roi; roi.xmin = 100; roi.xmax = 200; @@ -58,7 +56,7 @@ TEST_CASE("Simple ROIs on one module"){ REQUIRE(updated_geo.module_pixel_0[0].width == 100); } - SECTION("ROI is a small square"){ + SECTION("ROI is a small square") { aare::ROI roi; roi.xmin = 1000; roi.xmax = 1010; @@ -73,7 +71,7 @@ TEST_CASE("Simple ROIs on one module"){ REQUIRE(updated_geo.module_pixel_0[0].height == 10); REQUIRE(updated_geo.module_pixel_0[0].width == 10); } - SECTION("ROI is a few columns"){ + SECTION("ROI is a few columns") { aare::ROI roi; roi.xmin = 750; roi.xmax = 800; @@ -90,14 +88,12 @@ TEST_CASE("Simple ROIs on one module"){ } } - - -TEST_CASE("Two modules side by side"){ - // DetectorGeometry update_geometry_with_roi(DetectorGeometry geo, aare::ROI roi) +TEST_CASE("Two modules side by side") { + // DetectorGeometry update_geometry_with_roi(DetectorGeometry geo, aare::ROI + // roi) aare::DetectorGeometry geo; - - aare::ModuleGeometry mod; + aare::ModuleGeometry mod; mod.origin_x = 0; mod.origin_y = 0; mod.width = 1024; @@ -112,7 +108,7 @@ TEST_CASE("Two modules side by side"){ mod.origin_x = 1024; geo.module_pixel_0.push_back(mod); - SECTION("ROI is the whole image"){ + SECTION("ROI is the whole image") { aare::ROI roi; roi.xmin = 0; roi.xmax = 2048; @@ -125,7 +121,7 @@ TEST_CASE("Two modules side by side"){ REQUIRE(updated_geo.modules_x == 2); REQUIRE(updated_geo.modules_y == 1); } - SECTION("rectangle on both modules"){ + SECTION("rectangle on both modules") { aare::ROI roi; roi.xmin = 800; roi.xmax = 1300; @@ -141,11 +137,12 @@ TEST_CASE("Two modules side by side"){ REQUIRE(updated_geo.module_pixel_0[0].width == 224); REQUIRE(updated_geo.module_pixel_0[1].height == 299); REQUIRE(updated_geo.module_pixel_0[1].width == 276); - } + } } -TEST_CASE("Three modules side by side"){ - // DetectorGeometry update_geometry_with_roi(DetectorGeometry geo, aare::ROI roi) +TEST_CASE("Three modules side by side") { + // DetectorGeometry update_geometry_with_roi(DetectorGeometry geo, aare::ROI + // roi) aare::DetectorGeometry geo; aare::ROI roi; roi.xmin = 700; @@ -153,7 +150,7 @@ TEST_CASE("Three modules side by side"){ roi.ymin = 0; roi.ymax = 123; - aare::ModuleGeometry mod; + aare::ModuleGeometry mod; mod.origin_x = 0; mod.origin_y = 0; mod.width = 1024; @@ -184,8 +181,9 @@ TEST_CASE("Three modules side by side"){ REQUIRE(updated_geo.module_pixel_0[2].width == 452); } -TEST_CASE("Four modules as a square"){ - // DetectorGeometry update_geometry_with_roi(DetectorGeometry geo, aare::ROI roi) +TEST_CASE("Four modules as a square") { + // DetectorGeometry update_geometry_with_roi(DetectorGeometry geo, aare::ROI + // roi) aare::DetectorGeometry geo; aare::ROI roi; roi.xmin = 500; @@ -193,7 +191,7 @@ TEST_CASE("Four modules as a square"){ roi.ymin = 500; roi.ymax = 600; - aare::ModuleGeometry mod; + aare::ModuleGeometry mod; mod.origin_x = 0; mod.origin_y = 0; mod.width = 1024; diff --git a/src/utils/ifstream_helpers.cpp b/src/utils/ifstream_helpers.cpp index 74c56f3..79f006d 100644 --- a/src/utils/ifstream_helpers.cpp +++ b/src/utils/ifstream_helpers.cpp @@ -10,7 +10,7 @@ std::string ifstream_error_msg(std::ifstream &ifs) { return " Bad file stream"; } else if (state & std::ios_base::failbit) { return " File read failed"; - }else{ + } else { return " Unknown/no error"; } } diff --git a/src/utils/task.test.cpp b/src/utils/task.test.cpp index e19994a..3ca71c7 100644 --- a/src/utils/task.test.cpp +++ b/src/utils/task.test.cpp @@ -1,10 +1,9 @@ #include "aare/utils/task.hpp" -#include #include +#include - -TEST_CASE("Split a range into multiple tasks"){ +TEST_CASE("Split a range into multiple tasks") { auto tasks = aare::split_task(0, 10, 3); REQUIRE(tasks.size() == 3); @@ -22,11 +21,8 @@ TEST_CASE("Split a range into multiple tasks"){ tasks = aare::split_task(0, 10, 10); REQUIRE(tasks.size() == 10); - for (int i = 0; i < 10; i++){ + for (int i = 0; i < 10; i++) { REQUIRE(tasks[i].first == i); - REQUIRE(tasks[i].second == i+1); + REQUIRE(tasks[i].second == i + 1); } - - - } \ No newline at end of file diff --git a/tests/test.cpp b/tests/test.cpp index 513f690..1d4456f 100644 --- a/tests/test.cpp +++ b/tests/test.cpp @@ -2,8 +2,8 @@ #include #include #include -#include #include +#include TEST_CASE("Test suite can find data assets", "[.integration]") { auto fpath = test_data_path() / "numpy" / "test_numpy_file.npy";