Implement transformers to compute basic stats

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 <typename N>
-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()) {

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 <typename N = double>
-  [[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;

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"

src/libhictk/transformers/include/hictk/transformers/impl/stats_impl.hpp

@@ -0,0 +1,46 @@
+// Copyright (C) 2023 Roberto Rossini <[email protected]>
+//
+// SPDX-License-Identifier: MIT
+
+#pragma once
+
+#include <algorithm>
+#include <numeric>
+
+namespace hictk::transformers {
+
+template <typename PixelIt, typename N>
+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 <typename PixelIt, typename N>
+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 <typename PixelIt>
+inline std::size_t nnz(PixelIt first, PixelIt last) {
+  return static_cast<std::size_t>(std::distance(first, last));
+}
+
+template <typename PixelIt>
+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<double>(nnz);
+}
+
+}  // namespace hictk::transformers

src/libhictk/transformers/include/hictk/transformers/stats.hpp

@@ -0,0 +1,27 @@
+// Copyright (C) 2023 Roberto Rossini <[email protected]>
+//
+// SPDX-License-Identifier: MIT
+
+#pragma once
+
+#include <cstddef>
+#include <iterator>
+#include <type_traits>
+
+namespace hictk::transformers {
+
+template <typename PixelIt>
+[[nodiscard]] double avg(PixelIt first, PixelIt last);
+
+template <typename PixelIt, typename N = std::decay_t<decltype(std::declval<PixelIt>()->count)>>
+[[nodiscard]] N max(PixelIt first, PixelIt last);
+
+template <typename PixelIt>
+[[nodiscard]] std::size_t nnz(PixelIt first, PixelIt last);
+
+template <typename PixelIt, typename N = std::decay_t<decltype(std::declval<PixelIt>()->count)>>
+[[nodiscard]] N sum(PixelIt first, PixelIt last);
+
+}  // namespace hictk::transformers
+
+#include "./impl/stats_impl.hpp"

test/units/transformers/transformers_test.cpp

@@ -7,6 +7,7 @@
 #include <fmt/format.h>
 
 #include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_floating_point.hpp>
 
 #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<std::int32_t>();
+    auto last = sel.end<std::int32_t>();
+
+    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)