choose_filter_thresholds.hpp

Go to the documentation of this file.

#ifndef SCRAN_QC_CHOOSE_FILTER_THRESHOLDS_HPP
#define SCRAN_QC_CHOOSE_FILTER_THRESHOLDS_HPP
 
#include <vector>
#include <limits>
#include <cmath>
 
#include "find_median_mad.hpp"
 
namespace scran_qc {
 
struct ChooseFilterThresholdsOptions {
    bool lower = true;
 
    bool upper = true;
 
    double num_mads = 3;
 
    double min_diff = 0;
 
    bool log = false;
};
 
template<typename Float_>
struct ChooseFilterThresholdsResults {
    Float_ lower = 0;
 
    Float_ upper = 0;
};
 
namespace internal {
 
template<typename Float_>
Float_ unlog_threshold(Float_ val, bool was_logged) {
    if (was_logged) {
        if (std::isinf(val)) {
            if (val < 0) {
                return 0;
            }
        } else {
            return std::exp(val);
        }
    }
    return val;
}
 
template<bool lower_, typename Float_>
std::vector<Float_> strip_threshold(const std::vector<ChooseFilterThresholdsResults<Float_> >& res) {
    std::vector<Float_> output;
    output.reserve(res.size());
    for (const auto& r : res) {
        if constexpr(lower_) {
            output.push_back(r.lower);
        } else {
            output.push_back(r.upper);
        }
    }
    return output;
}
 
}
template<typename Float_>
ChooseFilterThresholdsResults<Float_> choose_filter_thresholds(const FindMedianMadResults<Float_>& mm, const ChooseFilterThresholdsOptions& options) {
    static_assert(std::is_floating_point<Float_>::value);
    ChooseFilterThresholdsResults<Float_> output;
    Float_& lthresh = output.lower;
    Float_& uthresh = output.upper;
    lthresh = -std::numeric_limits<Float_>::infinity();
    uthresh = std::numeric_limits<double>::infinity();
 
    auto median = mm.median;
    auto mad = mm.mad;
    if (!std::isnan(median) && !std::isnan(mad)) {
        auto delta = std::max(static_cast<Float_>(options.min_diff), options.num_mads * mad);
        if (options.lower) {
            lthresh = internal::unlog_threshold(median - delta, options.log);
        }
        if (options.upper) {
            uthresh = internal::unlog_threshold(median + delta, options.log);
        }
    }
 
    return output;
}
 
template<typename Index_, typename Float_>
ChooseFilterThresholdsResults<Float_> choose_filter_thresholds(Index_ num, Float_* metrics, const ChooseFilterThresholdsOptions& options) {
    FindMedianMadOptions fopt;
    fopt.log = options.log;
    auto mm = find_median_mad(num, metrics, fopt);
    return choose_filter_thresholds(mm, options);
}
 
template<typename Index_, typename Value_, typename Float_>
ChooseFilterThresholdsResults<Float_> choose_filter_thresholds(Index_ num, const Value_* metrics, Float_* buffer, const ChooseFilterThresholdsOptions& options) {
    FindMedianMadOptions fopt;
    fopt.log = options.log;
    auto mm = find_median_mad(num, metrics, buffer, fopt);
    return choose_filter_thresholds(mm, options);
}
 
template<typename Float_>
std::vector<ChooseFilterThresholdsResults<Float_> > choose_filter_thresholds_blocked(
    const std::vector<FindMedianMadResults<Float_> > mms,
    const ChooseFilterThresholdsOptions& options)
{
    std::vector<ChooseFilterThresholdsResults<Float_> > output;
    output.reserve(mms.size());
    for (auto& mm : mms) {
        output.emplace_back(choose_filter_thresholds(mm, options));
    }
    return output;
}
 
template<typename Index_, typename Value_, typename Block_, typename Float_>
std::vector<ChooseFilterThresholdsResults<Float_> > choose_filter_thresholds_blocked(
    Index_ num,
    const Value_* metrics,
    const Block_* block,
    FindMedianMadWorkspace<Float_, Index_>* workspace,
    const ChooseFilterThresholdsOptions& options)
{
    FindMedianMadOptions fopt;
    fopt.log = options.log;
    auto mms = find_median_mad_blocked(num, metrics, block, workspace, fopt);
    return choose_filter_thresholds_blocked(mms, options);
}
 
}
 
#endif