Merge pull request #176 from paulsengroup/feature/transformers

Refactor code to return interactions as sparse/dense Eigen matrices
paulsengroup · May 22, 2024 · 9f426e7 · 9f426e7
2 parents a019f1f + 4bad2f5
commit 9f426e7
Show file tree

Hide file tree

Showing 17 changed files with 406 additions and 333 deletions.
diff --git a/src/libhictk/cooler/include/hictk/cooler/impl/pixel_selector_impl.hpp b/src/libhictk/cooler/include/hictk/cooler/impl/pixel_selector_impl.hpp
@@ -124,65 +124,6 @@ inline std::vector<Pixel<N>> PixelSelector::read_all() const {
   return buff;
 }
 
-#ifdef HICTK_WITH_EIGEN
-template <typename N>
-inline Eigen::SparseMatrix<N> PixelSelector::read_sparse() const {
-  const auto bin_size = _bins->resolution();
-  const auto span1 = coord1().bin2.end() - coord1().bin1.start();
-  const auto span2 = coord2().bin2.end() - coord2().bin1.start();
-  const auto num_rows = static_cast<std::int64_t>((span1 + bin_size - 1) / bin_size);
-  const auto num_cols = static_cast<std::int64_t>((span2 + bin_size - 1) / bin_size);
-
-  const auto offset1 = coord1().bin1.id();
-  const auto offset2 = coord2().bin1.id();
-
-  Eigen::SparseMatrix<N> matrix(num_rows, num_cols);
-  std::for_each(begin<N>(), end<N>(), [&](const ThinPixel<N> &p) {
-    matrix.insert(static_cast<std::int64_t>(p.bin1_id - offset1),
-                  static_cast<std::int64_t>(p.bin2_id - offset2)) = p.count;
-  });
-  matrix.makeCompressed();
-  return matrix;
-}
-
-template <typename N>
-[[nodiscard]] Eigen::Matrix<N, Eigen::Dynamic, Eigen::Dynamic> PixelSelector::read_dense() const {
-  const auto bin_size = _bins->resolution();
-  const auto span1 = coord1().bin2.end() - coord1().bin1.start();
-  const auto span2 = coord2().bin2.end() - coord2().bin1.start();
-  const auto num_rows = static_cast<std::int64_t>((span1 + bin_size - 1) / bin_size);
-  const auto num_cols = static_cast<std::int64_t>((span2 + bin_size - 1) / bin_size);
-
-  const auto offset1 = coord1().bin1.id();
-  const auto offset2 = coord2().bin1.id();
-
-  const auto mirror_matrix = coord1().bin1.chrom() == coord2().bin1.chrom();
-
-  using MatrixT = Eigen::Matrix<N, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
-  MatrixT matrix = MatrixT::Zero(num_rows, num_cols);
-  std::for_each(begin<N>(), end<N>(), [&](const ThinPixel<N> &p) {
-    const auto i1 = static_cast<std::int64_t>(p.bin1_id - offset1);
-    const auto i2 = static_cast<std::int64_t>(p.bin2_id - offset2);
-    matrix(i1, i2) = p.count;
-
-    if (mirror_matrix) {
-      const auto delta = i2 - i1;
-      if (delta >= 0 && delta < num_rows && i1 < num_cols && i2 < num_rows) {
-        matrix(i2, i1) = p.count;
-      } else if ((delta < 0 || delta > num_cols) && i1 < num_cols && i2 < num_rows) {
-        const auto i3 = static_cast<std::int64_t>(p.bin2_id - offset1);
-        const auto i4 = static_cast<std::int64_t>(p.bin1_id - offset2);
-
-        if (i3 >= 0 && i3 < num_rows && i4 >= 0 && i4 < num_cols) {
-          matrix(i3, i4) = p.count;
-        }
-      }
-    }
-  });
-  return matrix;
-}
-#endif
-
 inline const PixelCoordinates &PixelSelector::coord1() const noexcept { return _coord1; }
 
 inline const PixelCoordinates &PixelSelector::coord2() const noexcept { return _coord2; }

diff --git a/src/libhictk/cooler/include/hictk/cooler/pixel_selector.hpp b/src/libhictk/cooler/include/hictk/cooler/pixel_selector.hpp
@@ -6,10 +6,6 @@
 
 // IWYU pragma: private, include "hictk/cooler.hpp"
 
-#ifdef HICTK_WITH_EIGEN
-#include <Eigen/Dense>
-#include <Eigen/SparseCore>
-#endif
 #include <cstddef>
 #include <cstdint>
 #include <iterator>
@@ -70,12 +66,6 @@ class PixelSelector {
 
   template <typename N>
   [[nodiscard]] std::vector<Pixel<N>> read_all() const;
-#ifdef HICTK_WITH_EIGEN
-  template <typename N>
-  [[nodiscard]] Eigen::SparseMatrix<N> read_sparse() const;
-  template <typename N>
-  [[nodiscard]] Eigen::Matrix<N, Eigen::Dynamic, Eigen::Dynamic> read_dense() const;
-#endif
 
   [[nodiscard]] const PixelCoordinates &coord1() const noexcept;
   [[nodiscard]] const PixelCoordinates &coord2() const noexcept;

diff --git a/src/libhictk/file/include/hictk/file.hpp b/src/libhictk/file/include/hictk/file.hpp
@@ -4,11 +4,6 @@
 
 #pragma once
 
-#ifdef HICTK_WITH_EIGEN
-#include <Eigen/Dense>
-#include <Eigen/SparseCore>
-#endif
-
 #include <cstddef>
 #include <cstdint>
 #include <iterator>
@@ -55,13 +50,6 @@ class PixelSelector {
   template <typename N>
   [[nodiscard]] std::vector<Pixel<N>> read_all() const;
 
-#ifdef HICTK_WITH_EIGEN
-  template <typename N>
-  [[nodiscard]] Eigen::SparseMatrix<N> read_sparse() const;
-  template <typename N>
-  [[nodiscard]] Eigen::Matrix<N, Eigen::Dynamic, Eigen::Dynamic> read_dense() const;
-#endif
-
   [[nodiscard]] const PixelCoordinates &coord1() const;
   [[nodiscard]] const PixelCoordinates &coord2() const;
 

diff --git a/src/libhictk/file/include/hictk/impl/file_impl.hpp b/src/libhictk/file/include/hictk/impl/file_impl.hpp
@@ -68,18 +68,6 @@ inline std::vector<Pixel<N>> PixelSelector::read_all() const {
   return std::visit([&](const auto& sel) { return sel.template read_all<N>(); }, _sel);
 }
 
-#ifdef HICTK_WITH_EIGEN
-template <typename N>
-inline Eigen::SparseMatrix<N> PixelSelector::read_sparse() const {
-  return std::visit([&](const auto& sel) { return sel.template read_sparse<N>(); }, _sel);
-}
-
-template <typename N>
-inline Eigen::Matrix<N, Eigen::Dynamic, Eigen::Dynamic> PixelSelector::read_dense() const {
-  return std::visit([&](const auto& sel) { return sel.template read_dense<N>(); }, _sel);
-}
-#endif
-
 inline const PixelCoordinates& PixelSelector::coord1() const {
   return std::visit(
       [&](const auto& sel) -> const PixelCoordinates& {

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
@@ -153,65 +153,6 @@ inline std::vector<Pixel<N>> PixelSelector::read_all() const {
   return buff;
 }
 
-#ifdef HICTK_WITH_EIGEN
-template <typename N>
-inline Eigen::SparseMatrix<N> PixelSelector::read_sparse() const {
-  const auto bin_size = bins().resolution();
-  const auto span1 = coord1().bin2.end() - coord1().bin1.start();
-  const auto span2 = coord2().bin2.end() - coord2().bin1.start();
-  const auto num_rows = static_cast<std::int64_t>((span1 + bin_size - 1) / bin_size);
-  const auto num_cols = static_cast<std::int64_t>((span2 + bin_size - 1) / bin_size);
-
-  const auto offset1 = coord1().bin1.id();
-  const auto offset2 = coord2().bin1.id();
-
-  Eigen::SparseMatrix<N> matrix(num_rows, num_cols);
-  std::for_each(begin<N>(), end<N>(), [&](const ThinPixel<N> &p) {
-    matrix.insert(static_cast<std::int64_t>(p.bin1_id - offset1),
-                  static_cast<std::int64_t>(p.bin2_id - offset2)) = p.count;
-  });
-  matrix.makeCompressed();
-  return matrix;
-}
-
-template <typename N>
-[[nodiscard]] Eigen::Matrix<N, Eigen::Dynamic, Eigen::Dynamic> PixelSelector::read_dense() const {
-  const auto bin_size = bins().resolution();
-  const auto span1 = coord1().bin2.end() - coord1().bin1.start();
-  const auto span2 = coord2().bin2.end() - coord2().bin1.start();
-  const auto num_rows = static_cast<std::int64_t>((span1 + bin_size - 1) / bin_size);
-  const auto num_cols = static_cast<std::int64_t>((span2 + bin_size - 1) / bin_size);
-
-  const auto offset1 = coord1().bin1.id();
-  const auto offset2 = coord2().bin1.id();
-
-  const auto mirror_matrix = coord1().bin1.chrom() == coord2().bin1.chrom();
-
-  using MatrixT = Eigen::Matrix<N, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
-  MatrixT matrix = MatrixT::Zero(num_rows, num_cols);
-  std::for_each(begin<N>(), end<N>(), [&](const ThinPixel<N> &p) {
-    const auto i1 = static_cast<std::int64_t>(p.bin1_id - offset1);
-    const auto i2 = static_cast<std::int64_t>(p.bin2_id - offset2);
-    matrix(i1, i2) = p.count;
-
-    if (mirror_matrix) {
-      const auto delta = i2 - i1;
-      if (delta >= 0 && delta < num_rows && i1 < num_cols && i2 < num_rows) {
-        matrix(i2, i1) = p.count;
-      } else if ((delta < 0 || delta > num_cols) && i1 < num_cols && i2 < num_rows) {
-        const auto i3 = static_cast<std::int64_t>(p.bin2_id - offset1);
-        const auto i4 = static_cast<std::int64_t>(p.bin1_id - offset2);
-
-        if (i3 >= 0 && i3 < num_rows && i4 >= 0 && i4 < num_cols) {
-          matrix(i3, i4) = p.count;
-        }
-      }
-    }
-  });
-  return matrix;
-}
-#endif
-
 inline const PixelCoordinates &PixelSelector::coord1() const noexcept { return *_coord1; }
 inline const PixelCoordinates &PixelSelector::coord2() const noexcept { return *_coord2; }
 inline MatrixType PixelSelector::matrix_type() const noexcept { return metadata().matrix_type; }
@@ -749,35 +690,6 @@ inline std::vector<Pixel<N>> PixelSelectorAll::read_all() const {
   return buff;
 }
 
-#ifdef HICTK_WITH_EIGEN
-template <typename N>
-inline Eigen::SparseMatrix<N> PixelSelectorAll::read_sparse() const {
-  const auto num_bins = static_cast<std::int64_t>(bins().size());
-  Eigen::SparseMatrix<N> matrix(num_bins, num_bins);
-  std::for_each(begin<N>(), end<N>(), [&](const ThinPixel<N> &p) {
-    matrix.insert(static_cast<std::int64_t>(p.bin1_id), static_cast<std::int64_t>(p.bin2_id)) =
-        p.count;
-  });
-  matrix.makeCompressed();
-  return matrix;
-}
-
-template <typename N>
-[[nodiscard]] Eigen::Matrix<N, Eigen::Dynamic, Eigen::Dynamic> PixelSelectorAll::read_dense()
-    const {
-  const auto num_bins = static_cast<std::int64_t>(bins().size());
-  using MatrixT = Eigen::Matrix<N, Eigen::Dynamic, Eigen::Dynamic>;
-  MatrixT matrix = MatrixT::Zero(num_bins, num_bins);
-  std::for_each(begin<N>(false), end<N>(), [&](const ThinPixel<N> &p) {
-    const auto i1 = static_cast<std::int64_t>(p.bin1_id);
-    const auto i2 = static_cast<std::int64_t>(p.bin2_id);
-    matrix(i1, i2) = p.count;
-    matrix(i2, i1) = p.count;
-  });
-  return matrix;
-}
-#endif
-
 inline MatrixType PixelSelectorAll::matrix_type() const noexcept {
   return _selectors.front().matrix_type();
 }

diff --git a/src/libhictk/hic/include/hictk/hic/pixel_selector.hpp b/src/libhictk/hic/include/hictk/hic/pixel_selector.hpp
@@ -10,10 +10,6 @@
 
 #include <cstddef>
 #include <cstdint>
-#ifdef HICTK_WITH_EIGEN
-#include <Eigen/Dense>
-#include <Eigen/SparseCore>
-#endif
 #include <functional>
 #include <iterator>
 #include <memory>
@@ -79,13 +75,6 @@ class PixelSelector {
   template <typename N>
   [[nodiscard]] std::vector<Pixel<N>> read_all() const;
 
-#ifdef HICTK_WITH_EIGEN
-  template <typename N>
-  [[nodiscard]] Eigen::SparseMatrix<N> read_sparse() const;
-  template <typename N>
-  [[nodiscard]] Eigen::Matrix<N, Eigen::Dynamic, Eigen::Dynamic> read_dense() const;
-#endif
-
   [[nodiscard]] const PixelCoordinates &coord1() const noexcept;
   [[nodiscard]] const PixelCoordinates &coord2() const noexcept;
 
@@ -209,13 +198,6 @@ class PixelSelectorAll {
   template <typename N>
   [[nodiscard]] std::vector<Pixel<N>> read_all() const;
 
-#ifdef HICTK_WITH_EIGEN
-  template <typename N>
-  [[nodiscard]] Eigen::SparseMatrix<N> read_sparse() const;
-  template <typename N>
-  [[nodiscard]] Eigen::Matrix<N, Eigen::Dynamic, Eigen::Dynamic> read_dense() const;
-#endif
-
   [[nodiscard]] MatrixType matrix_type() const noexcept;
   [[nodiscard]] balancing::Method normalization() const noexcept;
   [[nodiscard]] MatrixUnit unit() const noexcept;

diff --git a/src/libhictk/transformers/include/hictk/transformers.hpp b/src/libhictk/transformers/include/hictk/transformers.hpp
@@ -8,3 +8,5 @@
 #include "hictk/transformers/join_genomic_coords.hpp"
 #include "hictk/transformers/pixel_merger.hpp"
 #include "hictk/transformers/stats.hpp"
+#include "hictk/transformers/to_dense_matrix.hpp"
+#include "hictk/transformers/to_sparse_matrix.hpp"
diff --git a/src/libhictk/transformers/include/hictk/transformers/impl/common.hpp b/src/libhictk/transformers/include/hictk/transformers/impl/common.hpp
@@ -0,0 +1,16 @@
+// Copyright (C) 2024 Roberto Rossini <roberros@uio.no>
+//
+// SPDX-License-Identifier: MIT
+
+#pragma once
+
+#include <type_traits>
+
+namespace hictk::transformers::internal {
+template <typename T, typename = std::void_t<>>
+inline constexpr bool has_coord1_member_fx = false;
+
+template <typename T>
+inline constexpr bool has_coord1_member_fx<T, std::void_t<decltype(std::declval<T>().coord1())>> =
+    true;
+}