diff --git a/src/libhictk/hic/include/hictk/hic/impl/pixel_selector_impl.hpp b/src/libhictk/hic/include/hictk/hic/impl/pixel_selector_impl.hpp index aac7b08f..202764e7 100644 --- a/src/libhictk/hic/include/hictk/hic/impl/pixel_selector_impl.hpp +++ b/src/libhictk/hic/include/hictk/hic/impl/pixel_selector_impl.hpp @@ -211,12 +211,6 @@ inline bool PixelSelector::is_intra() const noexcept { return chrom1() == chrom2 inline bool PixelSelector::is_inter() const noexcept { return !is_intra(); } -template -inline N PixelSelector::sum() const noexcept { - return _reader.sum(); -} -inline double PixelSelector::avg() const noexcept { return _reader.avg(); } - inline std::size_t PixelSelector::estimate_optimal_cache_size( [[maybe_unused]] std::size_t num_samples) const { if (_reader.index().empty()) { diff --git a/src/libhictk/hic/include/hictk/hic/pixel_selector.hpp b/src/libhictk/hic/include/hictk/hic/pixel_selector.hpp index cfe0a298..f96f3070 100644 --- a/src/libhictk/hic/include/hictk/hic/pixel_selector.hpp +++ b/src/libhictk/hic/include/hictk/hic/pixel_selector.hpp @@ -97,9 +97,6 @@ class PixelSelector { [[nodiscard]] bool is_inter() const noexcept; [[nodiscard]] bool is_intra() const noexcept; - template - [[nodiscard]] N sum() const noexcept; - [[nodiscard]] double avg() const noexcept; [[nodiscard]] std::size_t estimate_optimal_cache_size(std::size_t num_samples = 500) const; void clear_cache() const; diff --git a/src/libhictk/transformers/include/hictk/transformers.hpp b/src/libhictk/transformers/include/hictk/transformers.hpp index 31459c42..2778cfa5 100644 --- a/src/libhictk/transformers/include/hictk/transformers.hpp +++ b/src/libhictk/transformers/include/hictk/transformers.hpp @@ -6,3 +6,4 @@ #include "hictk/transformers/coarsen.hpp" #include "hictk/transformers/join_genomic_coords.hpp" +#include "hictk/transformers/stats.hpp" diff --git a/src/libhictk/transformers/include/hictk/transformers/impl/stats_impl.hpp b/src/libhictk/transformers/include/hictk/transformers/impl/stats_impl.hpp new file mode 100644 index 00000000..3390118d --- /dev/null +++ b/src/libhictk/transformers/include/hictk/transformers/impl/stats_impl.hpp @@ -0,0 +1,46 @@ +// Copyright (C) 2023 Roberto Rossini +// +// SPDX-License-Identifier: MIT + +#pragma once + +#include +#include + +namespace hictk::transformers { + +template +inline N sum(PixelIt first, PixelIt last) { + return std::accumulate(first, last, N(0), + [&](const N accumulator, const auto& p) { return accumulator + p.count; }); +} + +template +inline N max(PixelIt first, PixelIt last) { + // I prefer using for_each instead of max_element to avoid keeping around a copy + // of the iterator for the current max_element + N max_ = 0; + std::for_each(first, last, [&](const auto& p) { max_ = std::max(max_, p.count); }); + return max_; +} + +template +inline std::size_t nnz(PixelIt first, PixelIt last) { + return static_cast(std::distance(first, last)); +} + +template +inline double avg(PixelIt first, PixelIt last) { + std::size_t nnz = 0; + const auto sum = + std::accumulate(first, last, double(0), [&](const double accumulator, const auto& p) { + ++nnz; + return accumulator + p.count; + }); + if (nnz == 0) { + return 0.0; + } + return sum / static_cast(nnz); +} + +} // namespace hictk::transformers diff --git a/src/libhictk/transformers/include/hictk/transformers/stats.hpp b/src/libhictk/transformers/include/hictk/transformers/stats.hpp new file mode 100644 index 00000000..c628e73a --- /dev/null +++ b/src/libhictk/transformers/include/hictk/transformers/stats.hpp @@ -0,0 +1,27 @@ +// Copyright (C) 2023 Roberto Rossini +// +// SPDX-License-Identifier: MIT + +#pragma once + +#include +#include +#include + +namespace hictk::transformers { + +template +[[nodiscard]] double avg(PixelIt first, PixelIt last); + +template ()->count)>> +[[nodiscard]] N max(PixelIt first, PixelIt last); + +template +[[nodiscard]] std::size_t nnz(PixelIt first, PixelIt last); + +template ()->count)>> +[[nodiscard]] N sum(PixelIt first, PixelIt last); + +} // namespace hictk::transformers + +#include "./impl/stats_impl.hpp" diff --git a/test/units/transformers/transformers_test.cpp b/test/units/transformers/transformers_test.cpp index 89882a0e..52e82baf 100644 --- a/test/units/transformers/transformers_test.cpp +++ b/test/units/transformers/transformers_test.cpp @@ -7,6 +7,7 @@ #include #include +#include #include "hictk/cooler/cooler.hpp" #include "hictk/hic.hpp" @@ -94,6 +95,19 @@ TEST_CASE("Transformers (cooler)", "[transformers][short]") { CHECK(v1[i] == v2[i].to_thin()); } } + + SECTION("stats") { + const auto path = datadir / "cooler/ENCFF993FGR.2500000.cool"; + auto clr = cooler::File::open(path.string()); + auto sel = clr.fetch("chr1"); + auto first = sel.begin(); + auto last = sel.end(); + + CHECK_THAT(avg(first, last), Catch::Matchers::WithinRel(25231.981858902574)); + CHECK(nnz(first, last) == 4'465); + CHECK(max(first, last) == 1'357'124); + CHECK(sum(first, last) == 112'660'799); + } } // NOLINTNEXTLINE(readability-function-cognitive-complexity)