adt_quality_control.hpp

Go to the documentation of this file.

#ifndef SCRAN_QC_ADT_QUALITY_CONTROL_HPP
#define SCRAN_QC_ADT_QUALITY_CONTROL_HPP
 
#include <vector>
#include <limits>
#include <algorithm>
#include <type_traits>
 
#include "tatami/tatami.hpp"
 
#include "find_median_mad.hpp"
#include "per_cell_qc_metrics.hpp"
#include "choose_filter_thresholds.hpp"
 
namespace scran_qc {
 
struct ComputeAdtQcMetricsOptions {
    int num_threads = 1;
};
 
template<typename Sum_ = double, typename Detected_ = int>
struct ComputeAdtQcMetricsBuffers {
    Sum_* sum;
 
    Detected_* detected;
 
    std::vector<Sum_*> subset_sum;
};
 
template<typename Value_, typename Index_, typename Subset_, typename Sum_, typename Detected_>
void compute_adt_qc_metrics(
    const tatami::Matrix<Value_, Index_>& mat,
    const std::vector<Subset_>& subsets,
    const ComputeAdtQcMetricsBuffers<Sum_, Detected_>& output,
    const ComputeAdtQcMetricsOptions& options)
{
    PerCellQcMetricsBuffers<Sum_, Detected_, Value_, Index_> tmp;
    tmp.sum = output.sum;
    tmp.detected = output.detected;
    tmp.subset_sum = output.subset_sum;
 
    PerCellQcMetricsOptions opt;
    opt.num_threads = options.num_threads;
    per_cell_qc_metrics(mat, subsets, tmp, opt);
}
 
template<typename Sum_ = double, typename Detected_ = int>
struct ComputeAdtQcMetricsResults {
    std::vector<Sum_> sum;
 
    std::vector<Detected_> detected;
 
    std::vector<std::vector<Sum_> > subset_sum;
};
 
template<typename Sum_ = double, typename Detected_ = int, typename Value_ = double, typename Index_ = int, typename Subset_ = const uint8_t*>
ComputeAdtQcMetricsResults<Sum_, Detected_> compute_adt_qc_metrics(
    const tatami::Matrix<Value_, Index_>& mat,
    const std::vector<Subset_>& subsets,
    const ComputeAdtQcMetricsOptions& options)
{
    auto NC = mat.ncol();
    ComputeAdtQcMetricsBuffers<Sum_, Detected_> x;
    ComputeAdtQcMetricsResults<Sum_, Detected_> output;
 
    output.sum.resize(NC
#ifdef SCRAN_QC_TEST_INIT
        , SCRAN_QC_TEST_INIT
#endif
    );
    x.sum = output.sum.data();
 
    output.detected.resize(NC
#ifdef SCRAN_QC_TEST_INIT
        , SCRAN_QC_TEST_INIT
#endif
    );
    x.detected = output.detected.data();
 
    size_t nsubsets = subsets.size();
    x.subset_sum.resize(nsubsets);
    output.subset_sum.resize(nsubsets);
    for (size_t s = 0; s < nsubsets; ++s) {
        output.subset_sum[s].resize(NC
#ifdef SCRAN_QC_TEST_INIT
            , SCRAN_QC_TEST_INIT
#endif
        );
        x.subset_sum[s] = output.subset_sum[s].data();
    }
 
    compute_adt_qc_metrics(mat, subsets, x, options);
    return output;
}
 
struct ComputeAdtQcFiltersOptions {
    double detected_num_mads = 3;
 
    double detected_min_drop = 0.1;
 
    double subset_sum_num_mads = 3;
};
 
namespace internal {
 
template<typename Float_, class Host_, typename Sum_, typename Detected_, typename BlockSource_>
void adt_populate(Host_& host, size_t n, const ComputeAdtQcMetricsBuffers<Sum_, Detected_>& res, BlockSource_ block, const ComputeAdtQcFiltersOptions& options) {
    constexpr bool unblocked = std::is_same<BlockSource_, bool>::value;
    auto buffer = [&]() {
        if constexpr(unblocked) {
            return std::vector<Float_>(n);
        } else {
            return FindMedianMadWorkspace<Float_, size_t>(n, block);
        }
    }();
 
    {
        ChooseFilterThresholdsOptions opts;
        opts.num_mads = options.detected_num_mads;
        opts.log = true;
        opts.upper = false;
        opts.min_diff = -std::log(1 - options.detected_min_drop);
        host.get_detected() = [&]() {
            if constexpr(unblocked) {
                return choose_filter_thresholds(n, res.detected, buffer.data(), opts).lower;
            } else {
                return internal::strip_threshold<true>(choose_filter_thresholds_blocked(n, res.detected, block, &buffer, opts));
            }
        }();
    }
 
    {
        size_t nsubsets = res.subset_sum.size();
        host.get_subset_sum().resize(nsubsets);
 
        ChooseFilterThresholdsOptions opts;
        opts.num_mads = options.subset_sum_num_mads;
        opts.log = true;
        opts.lower = false;
 
        for (size_t s = 0; s < nsubsets; ++s) {
            auto sub = res.subset_sum[s];
            host.get_subset_sum()[s] = [&]() {
                if constexpr(unblocked) {
                    return choose_filter_thresholds(n, sub, buffer.data(), opts).upper;
                } else {
                    return internal::strip_threshold<false>(choose_filter_thresholds_blocked(n, sub, block, &buffer, opts));
                }
            }();
        }
    }
}
 
template<class Host_, typename Sum_, typename Detected_, typename BlockSource_, typename Output_>
void adt_filter(const Host_& host, size_t n, const ComputeAdtQcMetricsBuffers<Sum_, Detected_>& metrics, BlockSource_ block, Output_* output) {
    constexpr bool unblocked = std::is_same<BlockSource_, bool>::value;
    std::fill_n(output, n, 1);
 
    const auto& detected = host.get_detected();
    for (size_t i = 0; i < n; ++i) {
        auto thresh = [&]() {
            if constexpr(unblocked) {
                return detected;
            } else {
                return detected[block[i]];
            }
        }();
        output[i] = output[i] && (metrics.detected[i] >= thresh);
    }
 
    size_t nsubsets = metrics.subset_sum.size();
    for (size_t s = 0; s < nsubsets; ++s) {
        auto sub = metrics.subset_sum[s];
        const auto& sthresh = host.get_subset_sum()[s];
        for (size_t i = 0; i < n; ++i) {
            auto thresh = [&]() {
                if constexpr(unblocked) {
                    return sthresh;
                } else {
                    return sthresh[block[i]];
                }
            }();
            output[i] = output[i] && (sub[i] <= thresh);
        }
    }
}
 
template<typename Sum_, typename Detected_>
ComputeAdtQcMetricsBuffers<const Sum_, const Detected_> to_buffer(const ComputeAdtQcMetricsResults<Sum_, Detected_>& metrics) {
    ComputeAdtQcMetricsBuffers<const Sum_, const Detected_> buffer;
    buffer.sum = metrics.sum.data();
    buffer.detected = metrics.detected.data();
    buffer.subset_sum.reserve(metrics.subset_sum.size());
    for (const auto& s : metrics.subset_sum) {
        buffer.subset_sum.push_back(s.data());
    }
    return buffer;
}
 
}
template<typename Float_ = double>
class AdtQcFilters {
public:
    Float_ get_detected() const {
        return my_detected;
    }
 
    const std::vector<Float_>& get_subset_sum() const {
        return my_subset_sum;
    }
 
    Float_& get_detected() {
        return my_detected;
    }
 
    std::vector<Float_>& get_subset_sum() {
        return my_subset_sum;
    }
 
private:
    Float_ my_detected = 0;
    std::vector<Float_> my_subset_sum;
 
public:
    template<typename Sum_, typename Detected_, typename Output_>
    void filter(size_t num, const ComputeAdtQcMetricsBuffers<Sum_, Detected_>& metrics, Output_* output) const {
        internal::adt_filter(*this, num, metrics, false, output);
    }
 
    template<typename Sum_, typename Detected_, typename Output_>
    void filter(const ComputeAdtQcMetricsResults<Sum_, Detected_>& metrics, Output_* output) const {
        return filter(metrics.detected.size(), internal::to_buffer(metrics), output);
    }
 
    template<typename Output_ = uint8_t, typename Sum_ = double, typename Detected_ = int>
    std::vector<Output_> filter(const ComputeAdtQcMetricsResults<Sum_, Detected_>& metrics) const {
        std::vector<Output_> output(metrics.detected.size()
#ifdef SCRAN_QC_TEST_INIT
            , SCRAN_QC_TEST_INIT
#endif
        );
        filter(metrics, output.data());
        return output;
    }
};
 
template<typename Float_ = double, typename Sum_, typename Detected_>
AdtQcFilters<Float_> compute_adt_qc_filters(size_t num, const ComputeAdtQcMetricsBuffers<Sum_, Detected_>& metrics, const ComputeAdtQcFiltersOptions& options) {
    AdtQcFilters<Float_> output;
    internal::adt_populate<Float_>(output, num, metrics, false, options);
    return output;
}
 
template<typename Float_ = double, typename Sum_, typename Detected_>
AdtQcFilters<Float_> compute_adt_qc_filters(const ComputeAdtQcMetricsResults<Sum_, Detected_>& metrics, const ComputeAdtQcFiltersOptions& options) {
    return compute_adt_qc_filters(metrics.detected.size(), internal::to_buffer(metrics), options);
}
 
template<typename Float_ = double>
class AdtQcBlockedFilters {
public:
    const std::vector<Float_>& get_detected() const {
        return my_detected;
    }
 
    const std::vector<std::vector<Float_> >& get_subset_sum() const {
        return my_subset_sum;
    }
 
    std::vector<Float_>& get_detected() {
        return my_detected;
    }
 
    std::vector<std::vector<Float_> >& get_subset_sum() {
        return my_subset_sum;
    }
 
private:
    std::vector<Float_> my_sum;
    std::vector<Float_> my_detected;
    std::vector<std::vector<Float_> > my_subset_sum;
 
public:
    template<typename Index_, typename Sum_, typename Detected_, typename Block_, typename Output_>
    void filter(Index_ num, const ComputeAdtQcMetricsBuffers<Sum_, Detected_>& metrics, const Block_* block, Output_* output) const {
        internal::adt_filter(*this, num, metrics, block, output);
    }
 
    template<typename Sum_, typename Detected_, typename Block_, typename Output_>
    void filter(const ComputeAdtQcMetricsResults<Sum_, Detected_>& metrics, const Block_* block, Output_* output) const {
        return filter(metrics.detected.size(), internal::to_buffer(metrics), block, output);
    }
 
    template<typename Output_ = uint8_t, typename Sum_ = double, typename Detected_ = int, typename Block_ = int>
    std::vector<Output_> filter(const ComputeAdtQcMetricsResults<Sum_, Detected_>& metrics, const Block_* block) const {
        std::vector<Output_> output(metrics.detected.size()
#ifdef SCRAN_QC_TEST_INIT
            , SCRAN_QC_TEST_INIT
#endif
        );
        filter(metrics, block, output.data());
        return output;
    }
};
 
template<typename Float_ = double, typename Sum_, typename Detected_, typename Block_>
AdtQcBlockedFilters<Float_> compute_adt_qc_filters_blocked(
    size_t num,
    const ComputeAdtQcMetricsBuffers<Sum_, Detected_>& metrics,
    const Block_* block,
    const ComputeAdtQcFiltersOptions& options)
{
    AdtQcBlockedFilters<Float_> output;
    internal::adt_populate<Float_>(output, num, metrics, block, options);
    return output;
}
 
template<typename Float_ = double, typename Sum_, typename Detected_, typename Block_>
AdtQcBlockedFilters<Float_> compute_adt_qc_filters_blocked(
    const ComputeAdtQcMetricsResults<Sum_, Detected_>& metrics,
    const Block_* block,
    const ComputeAdtQcFiltersOptions& options)
{
    return compute_adt_qc_filters_blocked(metrics.detected.size(), internal::to_buffer(metrics), block, options);
}
 
}
 
#endif