From e4d131f1bc88bcf56a779bfbcf4f46eca70f0e55 Mon Sep 17 00:00:00 2001 From: Roberto Rossini <71787608+robomics@users.noreply.github.com> Date: Mon, 11 Nov 2024 15:11:09 +0100 Subject: [PATCH] Fix GCC builds. Move a few templates from .hpp to .cpp files --- src/libnchg/include/nchg/nchg.hpp | 9 +- src/libnchg/nchg.cpp | 312 +++++++++++++++++++++++++++++ src/libnchg/nchg_impl.hpp | 322 ------------------------------ 3 files changed, 316 insertions(+), 327 deletions(-) diff --git a/src/libnchg/include/nchg/nchg.hpp b/src/libnchg/include/nchg/nchg.hpp index fe0e4a7..c9f782b 100644 --- a/src/libnchg/include/nchg/nchg.hpp +++ b/src/libnchg/include/nchg/nchg.hpp @@ -152,11 +152,10 @@ class NCHG { [[nodiscard]] auto aggregate_pixels(const hictk::GenomicInterval& range1, const hictk::GenomicInterval& range2) const; - template - requires arithmetic - [[nodiscard]] static NCHGResult compute_stats(hictk::Pixel pixel, double exp, N obs_sum, N N1, - N N2, double exp_sum, double L1, double L2, - std::vector& buffer); + [[nodiscard]] static NCHGResult compute_stats(hictk::Pixel pixel, double exp, + std::uint64_t obs_sum, std::uint64_t N1, + std::uint64_t N2, double exp_sum, double L1, + double L2, std::vector& buffer); public: template diff --git a/src/libnchg/nchg.cpp b/src/libnchg/nchg.cpp index 0762b94..2369ded 100644 --- a/src/libnchg/nchg.cpp +++ b/src/libnchg/nchg.cpp @@ -57,6 +57,138 @@ bool NCHGResult::operator==(const NCHGResult &other) const noexcept { bool NCHGResult::operator!=(const NCHGResult &other) const noexcept { return !(*this == other); } +[[nodiscard]] static std::pair aggregate_pixels_cis( + const hictk::File &f, const ExpectedMatrixStats &expected_matrix, + const hictk::GenomicInterval &range1, const hictk::GenomicInterval &range2, + std::uint64_t min_delta, std::uint64_t max_delta, const std::vector &bin_mask) { + assert(range1.chrom() == range2.chrom()); + assert(min_delta <= max_delta); + + std::uint64_t obs{}; + double exp{}; + + const auto query_offset1 = + static_cast(f.bins().at(range1.chrom(), range1.start()).rel_id()); + const auto query_offset2 = + static_cast(f.bins().at(range2.chrom(), range2.start()).rel_id()); + + const auto query_height = (range1.size() + f.resolution() - 1) / f.resolution(); + const auto query_width = (range2.size() + f.resolution() - 1) / f.resolution(); + + std::visit( + [&](const auto &f_) { + const auto sel = f_.fetch(range1.chrom().name(), range1.start(), range1.end(), + range2.chrom().name(), range2.start(), range2.end()); + + const hictk::transformers::JoinGenomicCoords jsel( + sel.template begin(), sel.template end(), f.bins_ptr()); + + for (const hictk::Pixel &p : jsel) { + const auto delta = p.coords.bin2.start() - p.coords.bin1.start(); + + const auto i1 = p.coords.bin1.rel_id(); + const auto i2 = p.coords.bin2.rel_id(); + + if (delta < min_delta || delta >= max_delta || bin_mask[i1] || bin_mask[i2]) + [[unlikely]] { + continue; + } + + const auto observed_count = p.count; + const auto expected_count = expected_matrix.at(i1, i2); + + const auto j1 = static_cast(i1) - query_offset1; + const auto j2 = static_cast(i2) - query_offset2; + const auto j3 = static_cast(i2) - query_offset1; + const auto j4 = static_cast(i1) - query_offset2; + + bool added = false; + if (j1 >= 0 && j1 < query_height && j2 >= 0 && j2 < query_width) { + obs += observed_count; + exp += expected_count; + added = true; + } + + if (added && j1 == j3 && j2 == j4) [[unlikely]] { + continue; + } + + if (j3 >= 0 && j3 < query_height && j4 >= 0 && j4 < query_width) { + obs += observed_count; + exp += expected_count; + } + } + }, + f.get()); + + return {obs, exp}; +} + +[[nodiscard]] static std::pair aggregate_pixels_trans( + const hictk::File &f, const ExpectedMatrixStats &expected_matrix, + const hictk::GenomicInterval &range1, const hictk::GenomicInterval &range2, + const std::vector &bin_mask1, const std::vector &bin_mask2) { + assert(range1.chrom() != range2.chrom()); + + std::uint64_t obs{}; + double exp{}; + + const auto chrom_offset1 = f.bins().at(range1.chrom(), 0).id(); + const auto chrom_offset2 = f.bins().at(range2.chrom(), 0).id(); + + std::visit( + [&](const auto &f_) { + const auto sel = f_.fetch(range1.chrom().name(), range1.start(), range1.end(), + range2.chrom().name(), range2.start(), range2.end()); + + std::for_each(sel.template begin(), sel.template end(), + [&](const hictk::ThinPixel &p) { + assert(p.bin1_id >= chrom_offset1); + assert(p.bin2_id >= chrom_offset2); + const auto i1 = p.bin1_id - chrom_offset1; + const auto i2 = p.bin2_id - chrom_offset2; + + if (!bin_mask1[i1] && !bin_mask2[i2]) [[likely]] { + obs += p.count; + exp += expected_matrix.at(i1, i2); + } + }); + }, + f.get()); + + return {obs, exp}; +} + +auto NCHG::aggregate_pixels(const hictk::GenomicInterval &range1, + const hictk::GenomicInterval &range2) const { + assert(_fp); + assert(_exp_matrix); + + struct Result { + std::uint64_t obs{}; + double exp{}; + }; + + if (range1.chrom() == range2.chrom()) { + const auto &mask = _expected_values.bin_mask(range1.chrom()); + assert(mask); + const auto [obs, exp] = aggregate_pixels_cis(*_fp, *_exp_matrix, range1, range2, + params().min_delta, params().max_delta, *mask); + + return Result{obs, exp}; + } + + const auto &[mask1, mask2] = _expected_values.bin_mask(range1.chrom(), range2.chrom()); + + assert(mask1); + assert(mask2); + + const auto [obs, exp] = + aggregate_pixels_trans(*_fp, *_exp_matrix, range1, range2, *mask1, *mask2); + + return Result{obs, exp}; +} + NCHG::NCHG(std::shared_ptr f, const hictk::Chromosome &chrom1, const hictk::Chromosome &chrom2, const Params ¶ms) : NCHG(f, chrom1, chrom2, ExpectedValues::chromosome_pair(f, chrom1, chrom2, params)) {} @@ -266,4 +398,184 @@ double NCHG::compute_pvalue_nchg(std::vector &buffer, std::uint64_t obs, return pvalue < 0 ? 1.0 : pvalue; } +template + requires arithmetic +[[nodiscard]] constexpr double compute_odds_ratio(N n, N total_sum, N sum1, N sum2) { + if constexpr (std::is_floating_point_v) { + if (std::isnan(n)) [[unlikely]] { + return std::numeric_limits::quiet_NaN(); + } + } + + if (sum1 == 0 || sum2 == 0) [[unlikely]] { + return std::numeric_limits::quiet_NaN(); + } + + assert(n <= sum1); + assert(n <= sum2); + assert(sum1 <= total_sum); + assert(sum2 <= total_sum); + assert(sum1 + sum2 <= total_sum + n); + + const auto num = n * (total_sum - sum1 - sum2 + n); + const auto denom = (sum1 - n) * (sum2 - n); + + assert(num >= 0); + assert(denom >= 0); + + return conditional_static_cast(num) / conditional_static_cast(denom); +} + +template , double, std::uint64_t>> + requires arithmetic +constexpr std::pair aggregate_marginals(const hictk::GenomicInterval &range, + std::uint32_t resolution, + const std::vector &marginals, + const std::vector &bin_mask) noexcept { + assert(resolution > 0); + + const auto i1 = range.start() / resolution; + const auto i2 = (range.end() + resolution - 1) / resolution; + + if (i1 == i2) [[unlikely]] { + return {N2{}, N2{}}; + } + + assert(i2 <= marginals.size()); + assert(i2 <= bin_mask.size()); + + std::size_t bin_masked = 0; + N2 sum{}; + for (auto i = i1; i < i2; ++i) { + sum += conditional_static_cast(marginals[i]); + bin_masked += bin_mask[i]; + } + + const auto bin_masked_frac = static_cast(bin_masked) / static_cast(i2 - i1); + + return {bin_masked_frac, sum}; +} + +std::uint64_t NCHG::compute_N1(const hictk::GenomicInterval &range1, + const hictk::GenomicInterval &range2, + double max_bad_bin_threshold) const noexcept { + assert(_obs_matrix); + assert(_fp); + const auto &[bad_bin_frac, sum] = + aggregate_marginals(range1, _fp->resolution(), _obs_matrix->marginals1(), + *_expected_values.bin_mask(range1.chrom(), range2.chrom()).first); + + if (bad_bin_frac >= max_bad_bin_threshold) { + return std::numeric_limits::max(); + } + + return sum; +} + +std::uint64_t NCHG::compute_N2(const hictk::GenomicInterval &range1, + const hictk::GenomicInterval &range2, + double max_bad_bin_threshold) const noexcept { + assert(_obs_matrix); + assert(_fp); + const auto &[bad_bin_frac, sum] = + aggregate_marginals(range2, _fp->resolution(), _obs_matrix->marginals2(), + *_expected_values.bin_mask(range1.chrom(), range2.chrom()).second); + + if (bad_bin_frac >= max_bad_bin_threshold) { + return std::numeric_limits::max(); + } + + return sum; +} + +double NCHG::compute_L1(const hictk::GenomicInterval &range1, const hictk::GenomicInterval &range2, + double max_bad_bin_threshold) const noexcept { + assert(_exp_matrix); + assert(_fp); + const auto &[bad_bin_frac, sum] = + aggregate_marginals(range1, _fp->resolution(), _exp_matrix->marginals1(), + *_expected_values.bin_mask(range1.chrom(), range2.chrom()).first); + + if (bad_bin_frac >= max_bad_bin_threshold) { + return std::numeric_limits::quiet_NaN(); + } + + return sum; +} + +double NCHG::compute_L2(const hictk::GenomicInterval &range1, const hictk::GenomicInterval &range2, + double max_bad_bin_threshold) const noexcept { + assert(_exp_matrix); + assert(_fp); + const auto &[bad_bin_frac, sum] = + aggregate_marginals(range2, _fp->resolution(), _exp_matrix->marginals2(), + *_expected_values.bin_mask(range1.chrom(), range2.chrom()).second); + + if (bad_bin_frac >= max_bad_bin_threshold) { + return std::numeric_limits::quiet_NaN(); + } + + return sum; +} + +NCHGResult NCHG::compute_stats(hictk::Pixel pixel, double exp, std::uint64_t obs_sum, + std::uint64_t N1, std::uint64_t N2, double exp_sum, double L1, + double L2, std::vector &buffer) { + constexpr auto bad_sum = std::numeric_limits::max(); + if (pixel.count == 0 || N1 == bad_sum || N2 == bad_sum) [[unlikely]] { + return {.pixel = std::move(pixel), + .expected = exp, + .pval = 1.0, + .log_ratio = 0.0, + .odds_ratio = 0.0, + .omega = 0.0}; + } + + assert(std::isfinite(exp)); + assert(std::isfinite(L1)); + assert(std::isfinite(L2)); + assert(std::isfinite(exp_sum)); + + const auto intra_matrix = pixel.coords.bin1.chrom() == pixel.coords.bin2.chrom(); + + if (!intra_matrix) { + obs_sum *= 2; + } + + const auto odds_ratio = compute_odds_ratio(pixel.count, obs_sum, N1, N2); + const auto omega = intra_matrix ? compute_odds_ratio(exp, exp_sum, L1, L2) : 1.0; + const auto log_ratio = std::log2(static_cast(pixel.count)) - std::log2(exp); + + if (!std::isfinite(odds_ratio) || !std::isfinite(omega)) [[unlikely]] { + return {.pixel = std::move(pixel), + .expected = exp, + .pval = 1.0, + .log_ratio = log_ratio, + .odds_ratio = odds_ratio, + .omega = omega}; + } + + if (!std::isfinite(omega) || omega > odds_ratio) [[unlikely]] { + return {.pixel = std::move(pixel), + .expected = exp, + .pval = 1.0, + .log_ratio = log_ratio, + .odds_ratio = odds_ratio, + .omega = omega}; + } + + if (!intra_matrix) { + obs_sum /= 2; + } + + const auto pvalue = compute_pvalue_nchg(buffer, pixel.count, N1, N2, obs_sum, omega); + return {.pixel = std::move(pixel), + .expected = exp, + .pval = pvalue, + .log_ratio = log_ratio, + .odds_ratio = odds_ratio, + .omega = omega}; +} + } // namespace nchg diff --git a/src/libnchg/nchg_impl.hpp b/src/libnchg/nchg_impl.hpp index d0b4cad..fba0d88 100644 --- a/src/libnchg/nchg_impl.hpp +++ b/src/libnchg/nchg_impl.hpp @@ -35,328 +35,6 @@ #include "nchg/observed_matrix.hpp" namespace nchg { -namespace internal { - -template - requires arithmetic -[[nodiscard]] constexpr double compute_odds_ratio(N n, N total_sum, N sum1, N sum2) { - if constexpr (std::is_floating_point_v) { - if (std::isnan(n)) [[unlikely]] { - return std::numeric_limits::quiet_NaN(); - } - } - - if (sum1 == 0 || sum2 == 0) [[unlikely]] { - return std::numeric_limits::quiet_NaN(); - } - - assert(n <= sum1); - assert(n <= sum2); - assert(sum1 <= total_sum); - assert(sum2 <= total_sum); - assert(sum1 + sum2 <= total_sum + n); - - const auto num = n * (total_sum - sum1 - sum2 + n); - const auto denom = (sum1 - n) * (sum2 - n); - - assert(num >= 0); - assert(denom >= 0); - - return conditional_static_cast(num) / conditional_static_cast(denom); -} - -// GCC gets confused if we declare this as a member function -// (static, not static, inline, constexpr, noexcept(true), noexcept(false)... It does not matter) -template , double, std::uint64_t>> - requires arithmetic -constexpr std::pair aggregate_marginals(const hictk::GenomicInterval &range, - std::uint32_t resolution, - const std::vector &marginals, - const std::vector &bin_mask) noexcept { - assert(resolution > 0); - - const auto i1 = range.start() / resolution; - const auto i2 = (range.end() + resolution - 1) / resolution; - - if (i1 == i2) [[unlikely]] { - return {N2{}, N2{}}; - } - - assert(i2 <= marginals.size()); - assert(i2 <= bin_mask.size()); - - std::size_t bin_masked = 0; - N2 sum{}; - for (auto i = i1; i < i2; ++i) { - sum += conditional_static_cast(marginals[i]); - bin_masked += bin_mask[i]; - } - - const auto bin_masked_frac = static_cast(bin_masked) / static_cast(i2 - i1); - - return {bin_masked_frac, sum}; -} - -[[nodiscard]] inline std::pair aggregate_pixels_cis( - const hictk::File &f, const ExpectedMatrixStats &expected_matrix, - const hictk::GenomicInterval &range1, const hictk::GenomicInterval &range2, - std::uint64_t min_delta, std::uint64_t max_delta, const std::vector &bin_mask) { - assert(range1.chrom() == range2.chrom()); - assert(min_delta <= max_delta); - - std::uint64_t obs{}; - double exp{}; - - const auto query_offset1 = - static_cast(f.bins().at(range1.chrom(), range1.start()).rel_id()); - const auto query_offset2 = - static_cast(f.bins().at(range2.chrom(), range2.start()).rel_id()); - - const auto query_height = (range1.size() + f.resolution() - 1) / f.resolution(); - const auto query_width = (range2.size() + f.resolution() - 1) / f.resolution(); - - std::visit( - [&](const auto &f_) { - const auto sel = f_.fetch(range1.chrom().name(), range1.start(), range1.end(), - range2.chrom().name(), range2.start(), range2.end()); - - const hictk::transformers::JoinGenomicCoords jsel( - sel.template begin(), sel.template end(), f.bins_ptr()); - - for (const hictk::Pixel &p : jsel) { - const auto delta = p.coords.bin2.start() - p.coords.bin1.start(); - - const auto i1 = p.coords.bin1.rel_id(); - const auto i2 = p.coords.bin2.rel_id(); - - if (delta < min_delta || delta >= max_delta || bin_mask[i1] || bin_mask[i2]) - [[unlikely]] { - continue; - } - - const auto observed_count = p.count; - const auto expected_count = expected_matrix.at(i1, i2); - - const auto j1 = static_cast(i1) - query_offset1; - const auto j2 = static_cast(i2) - query_offset2; - const auto j3 = static_cast(i2) - query_offset1; - const auto j4 = static_cast(i1) - query_offset2; - - bool added = false; - if (j1 >= 0 && j1 < query_height && j2 >= 0 && j2 < query_width) { - obs += observed_count; - exp += expected_count; - added = true; - } - - if (added && j1 == j3 && j2 == j4) [[unlikely]] { - continue; - } - - if (j3 >= 0 && j3 < query_height && j4 >= 0 && j4 < query_width) { - obs += observed_count; - exp += expected_count; - } - } - }, - f.get()); - - return {obs, exp}; -} - -[[nodiscard]] inline std::pair aggregate_pixels_trans( - const hictk::File &f, const ExpectedMatrixStats &expected_matrix, - const hictk::GenomicInterval &range1, const hictk::GenomicInterval &range2, - const std::vector &bin_mask1, const std::vector &bin_mask2) { - assert(range1.chrom() != range2.chrom()); - - std::uint64_t obs{}; - double exp{}; - - const auto chrom_offset1 = f.bins().at(range1.chrom(), 0).id(); - const auto chrom_offset2 = f.bins().at(range2.chrom(), 0).id(); - - std::visit( - [&](const auto &f_) { - const auto sel = f_.fetch(range1.chrom().name(), range1.start(), range1.end(), - range2.chrom().name(), range2.start(), range2.end()); - - std::for_each(sel.template begin(), sel.template end(), - [&](const hictk::ThinPixel &p) { - assert(p.bin1_id >= chrom_offset1); - assert(p.bin2_id >= chrom_offset2); - const auto i1 = p.bin1_id - chrom_offset1; - const auto i2 = p.bin2_id - chrom_offset2; - - if (!bin_mask1[i1] && !bin_mask2[i2]) [[likely]] { - obs += p.count; - exp += expected_matrix.at(i1, i2); - } - }); - }, - f.get()); - - return {obs, exp}; -} - -} // namespace internal - -inline std::uint64_t NCHG::compute_N1(const hictk::GenomicInterval &range1, - const hictk::GenomicInterval &range2, - double max_bad_bin_threshold) const noexcept { - assert(_obs_matrix); - assert(_fp); - const auto &[bad_bin_frac, sum] = internal::aggregate_marginals( - range1, _fp->resolution(), _obs_matrix->marginals1(), - *_expected_values.bin_mask(range1.chrom(), range2.chrom()).first); - - if (bad_bin_frac >= max_bad_bin_threshold) { - return std::numeric_limits::max(); - } - - return sum; -} - -inline std::uint64_t NCHG::compute_N2(const hictk::GenomicInterval &range1, - const hictk::GenomicInterval &range2, - double max_bad_bin_threshold) const noexcept { - assert(_obs_matrix); - assert(_fp); - const auto &[bad_bin_frac, sum] = internal::aggregate_marginals( - range2, _fp->resolution(), _obs_matrix->marginals2(), - *_expected_values.bin_mask(range1.chrom(), range2.chrom()).second); - - if (bad_bin_frac >= max_bad_bin_threshold) { - return std::numeric_limits::max(); - } - - return sum; -} - -inline double NCHG::compute_L1(const hictk::GenomicInterval &range1, - const hictk::GenomicInterval &range2, - double max_bad_bin_threshold) const noexcept { - assert(_exp_matrix); - assert(_fp); - const auto &[bad_bin_frac, sum] = internal::aggregate_marginals( - range1, _fp->resolution(), _exp_matrix->marginals1(), - *_expected_values.bin_mask(range1.chrom(), range2.chrom()).first); - - if (bad_bin_frac >= max_bad_bin_threshold) { - return std::numeric_limits::quiet_NaN(); - } - - return sum; -} - -inline double NCHG::compute_L2(const hictk::GenomicInterval &range1, - const hictk::GenomicInterval &range2, - double max_bad_bin_threshold) const noexcept { - assert(_exp_matrix); - assert(_fp); - const auto &[bad_bin_frac, sum] = internal::aggregate_marginals( - range2, _fp->resolution(), _exp_matrix->marginals2(), - *_expected_values.bin_mask(range1.chrom(), range2.chrom()).second); - - if (bad_bin_frac >= max_bad_bin_threshold) { - return std::numeric_limits::quiet_NaN(); - } - - return sum; -} - -inline auto NCHG::aggregate_pixels(const hictk::GenomicInterval &range1, - const hictk::GenomicInterval &range2) const { - assert(_fp); - assert(_exp_matrix); - - struct Result { - std::uint64_t obs{}; - double exp{}; - }; - - if (range1.chrom() == range2.chrom()) { - const auto &mask = _expected_values.bin_mask(range1.chrom()); - assert(mask); - const auto [obs, exp] = internal::aggregate_pixels_cis( - *_fp, *_exp_matrix, range1, range2, params().min_delta, params().max_delta, *mask); - - return Result{obs, exp}; - } - - const auto &[mask1, mask2] = _expected_values.bin_mask(range1.chrom(), range2.chrom()); - - assert(mask1); - assert(mask2); - - const auto [obs, exp] = - internal::aggregate_pixels_trans(*_fp, *_exp_matrix, range1, range2, *mask1, *mask2); - - return Result{obs, exp}; -} - -template - requires arithmetic -inline NCHGResult NCHG::compute_stats(hictk::Pixel pixel, double exp, N obs_sum, N N1, N N2, - double exp_sum, double L1, double L2, - std::vector &buffer) { - constexpr auto bad_sum = std::numeric_limits::max(); - if (pixel.count == 0 || N1 == bad_sum || N2 == bad_sum) [[unlikely]] { - return {.pixel = std::move(pixel), - .expected = exp, - .pval = 1.0, - .log_ratio = 0.0, - .odds_ratio = 0.0, - .omega = 0.0}; - } - - assert(std::isfinite(exp)); - assert(std::isfinite(L1)); - assert(std::isfinite(L2)); - assert(std::isfinite(exp_sum)); - - const auto intra_matrix = pixel.coords.bin1.chrom() == pixel.coords.bin2.chrom(); - - if (!intra_matrix) { - obs_sum *= 2; - } - - const auto odds_ratio = internal::compute_odds_ratio(pixel.count, obs_sum, N1, N2); - const auto omega = intra_matrix ? internal::compute_odds_ratio(exp, exp_sum, L1, L2) : 1.0; - const auto log_ratio = std::log2(static_cast(pixel.count)) - std::log2(exp); - - if (!std::isfinite(odds_ratio) || !std::isfinite(omega)) [[unlikely]] { - return {.pixel = std::move(pixel), - .expected = exp, - .pval = 1.0, - .log_ratio = log_ratio, - .odds_ratio = odds_ratio, - .omega = omega}; - } - - if (!std::isfinite(omega) || omega > odds_ratio) [[unlikely]] { - return {.pixel = std::move(pixel), - .expected = exp, - .pval = 1.0, - .log_ratio = log_ratio, - .odds_ratio = odds_ratio, - .omega = omega}; - } - - if (!intra_matrix) { - obs_sum /= 2; - } - - const auto pvalue = compute_pvalue_nchg(buffer, pixel.count, static_cast(N1), - static_cast(N2), obs_sum, omega); - return {.pixel = std::move(pixel), - .expected = exp, - .pval = pvalue, - .log_ratio = log_ratio, - .odds_ratio = odds_ratio, - .omega = omega}; -} template inline NCHG::iterator::iterator(PixelSelector selector,