Merge pull request #98 from paulsengroup/bug/hic-fetch-expected

Fix bug when dumping expected .hic matrices
paulsengroup · Jan 9, 2024 · 40f8743 · 40f8743
2 parents 446d174 + 48ef784
commit 40f8743
Show file tree

Hide file tree

Showing 4 changed files with 12 additions and 9 deletions.
diff --git a/src/libhictk/config/version.hpp.in b/src/libhictk/config/version.hpp.in
@@ -46,7 +46,7 @@ namespace internal {
   }
   buff.append(is_dirty_suffix);
 
-  if (buff.front() == '-') {
+  if (!buff.empty() && buff.front() == '-') {
     buff.erase(0, 1);
   }
 

diff --git a/src/libhictk/hic/include/hictk/hic/impl/block_reader_impl.hpp b/src/libhictk/hic/include/hictk/hic/impl/block_reader_impl.hpp
@@ -69,8 +69,11 @@ inline double HiCBlockReader::avg() const {
   }
 
   const auto bin_size = bins().bin_size();
-  const auto num_bins1 = (chrom1().size() + bin_size - 1) / bin_size;
-  const auto num_bins2 = (chrom2().size() + bin_size - 1) / bin_size;
+  // We round down for two reasons:
+  // - to be consistent with straw
+  // - because the last bin is usually smaller than the bin_size
+  const auto num_bins1 = chrom1().size() / bin_size;
+  const auto num_bins2 = chrom2().size() / bin_size;
 
   return sum() / double(num_bins1 * num_bins2);
 }

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
@@ -124,8 +124,8 @@ inline ThinPixel<N> PixelSelector::transform_pixel(ThinPixel<float> pixel) const
       return float(_reader->avg());
     }
 
-    const auto i = (bin2 - bin1);
-    assert(i < expected.size());
+    const auto i =
+        std::min(bin2 - bin1, conditional_static_cast<std::uint64_t>(expected.size() - 1));
     return float(expected[i]);
   }();
 
@@ -694,8 +694,8 @@ inline ThinPixel<N> PixelSelector::iterator<N>::transform_pixel(ThinPixel<float>
       return float(_reader->avg());
     }
 
-    const auto i = (bin2 - bin1);
-    assert(i < expected.size());
+    const auto i =
+        std::min(bin2 - bin1, conditional_static_cast<std::uint64_t>(expected.size() - 1));
     return float(expected[i]);
   }();
 

diff --git a/test/units/hic/pixel_selector_test.cpp b/test/units/hic/pixel_selector_test.cpp
@@ -359,7 +359,7 @@ TEST_CASE("HiC: pixel selector fetch (expected NONE BP 10000)", "[hic][long]") {
       }
       SECTION("inter-chromosomal") {
         constexpr std::size_t expected_size = 56743;
-        constexpr double expected_sum = 12610.80619812;
+        constexpr double expected_sum = 12710.32078149;
         auto sel = File(path, 10'000, MatrixType::expected, MatrixUnit::BP)
                        .fetch("chr2L", "chr4", hictk::balancing::Method::NONE());
         const auto buffer = sel.read_all<double>();
@@ -387,7 +387,7 @@ TEST_CASE("HiC: pixel selector fetch (oe NONE BP 10000)", "[hic][long]") {
     }
     SECTION("inter-chromosomal") {
       constexpr std::size_t expected_size = 56743;
-      constexpr double expected_sum = 317520.00459671;
+      constexpr double expected_sum = 315034.01705551;
       auto sel = File(path, 10'000, MatrixType::oe, MatrixUnit::BP)
                      .fetch("chr2L", "chr4", hictk::balancing::Method::NONE());
       const auto buffer = sel.read_all<double>();