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 <algorithm>
#include <cstddef>
#include <cassert>
 
#include "sanisizer/sanisizer.hpp"
#include "quickstats/quickstats.hpp"
 
#include "utils.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;
};
 
template<typename Float_>
Float_ unlog_threshold(const Float_ val, const bool was_logged) {
    if (was_logged) {
        if (std::isinf(val)) {
            if (val < 0) {
                return 0;
            }
        } else {
            return std::exp(val);
        }
    }
    return val;
}
 
template<typename Float_>
ChooseFilterThresholdsResults<Float_> choose_filter_thresholds_internal(
    std::size_t num_cells,
    Float_* metrics,
    const ChooseFilterThresholdsOptions& options
) {
    static_assert(std::is_floating_point<Float_>::value);
 
    // Rotate all the NaNs to the front of the buffer and ignore them.
    I<decltype(num_cells)> lost = 0;
    for (I<decltype(num_cells)> i = 0; i < num_cells; ++i) {
        if (std::isnan(metrics[i])) {
            std::swap(metrics[i], metrics[lost]);
            ++lost;
        }
    }
    metrics += lost;
    num_cells -= lost;
 
    // Maybe we do some log-transformation, if that's requested.
    if (options.log) {
        for (I<decltype(num_cells)> i = 0; i < num_cells; ++i) {
            auto& val = metrics[i];
            if (val > 0) {
                val = std::log(val);
            } else if (val == 0) {
                val = -std::numeric_limits<double>::infinity();
            } else {
                throw std::runtime_error("cannot log-transform negative values");
            }
        }
    }
 
    ChooseFilterThresholdsResults<Float_> output;
    Float_& lthresh = output.lower;
    Float_& uthresh = output.upper;
    lthresh = -std::numeric_limits<Float_>::infinity();
    uthresh = std::numeric_limits<Float_>::infinity();
 
    quickstats::MedianOptions<Float_> medopt;
    medopt.placeholder = std::numeric_limits<Float_>::quiet_NaN();
    const auto median = quickstats::median<Float_>(num_cells, metrics, medopt);
 
    quickstats::MadOptions<Float_> madopt;
    madopt.placeholder = std::numeric_limits<Float_>::quiet_NaN();
    madopt.difference_between_infinities_is_zero = true; // for sane handling of log(0).
    const auto mad = quickstats::scale_mad_to_sd(quickstats::mad<Float_>(num_cells, metrics, median, madopt));
 
    if (!std::isnan(mad)) {
        const auto delta = std::max(static_cast<Float_>(options.min_diff), static_cast<Float_>(options.num_mads * mad));
        if (options.lower) {
            const auto threshold = median - delta;
            if (!std::isnan(threshold)) {
                lthresh = unlog_threshold(threshold, options.log);
            }
        }
        if (options.upper) {
            const auto threshold = median + delta;
            if (!std::isnan(threshold)) {
                uthresh = unlog_threshold(threshold, options.log);
            }
        }
    }
 
    return output;
}
template<typename Value_, typename Float_>
ChooseFilterThresholdsResults<Float_> choose_filter_thresholds(
    const std::size_t num_cells,
    const Value_* const metrics,
    Float_* const buffer,
    const ChooseFilterThresholdsOptions& options
) {
    // Only copying it if it's not exactly the same.
    [&](){
        if constexpr(std::is_same<Value_, Float_>::value) {
            if (metrics == buffer) {
                return;
            }
        }
        std::copy_n(metrics, num_cells, buffer);
    }();
    return choose_filter_thresholds_internal(num_cells, buffer, options);
}
 
template<typename Float_>
struct ChooseFilterThresholdsBlockedWorkspace {
    template<typename Block_>
    ChooseFilterThresholdsBlockedWorkspace(const std::size_t num_cells, const Block_* const block, const std::size_t num_blocks) {
        reset_choose_filter_thresholds_blocked_workspace(*this, num_cells, block, num_blocks);
    }
 
    ChooseFilterThresholdsBlockedWorkspace() = default;
 
public:
    std::vector<Float_> buffer;
    std::vector<std::size_t> block_starts;
    std::vector<std::size_t> block_offsets;
};
 
template<typename Float_, typename Block_>
void reset_choose_filter_thresholds_blocked_workspace(
    ChooseFilterThresholdsBlockedWorkspace<Float_>& work,
    const std::size_t num_cells,
    const Block_* const block,
    const std::size_t num_blocks
) {
    work.block_starts.clear();
 
    sanisizer::resize(work.block_starts, num_blocks);
    for (I<decltype(num_cells)> i = 0; i < num_cells; ++i) {
        ++work.block_starts[block[i]];
    }
 
    std::size_t sofar = 0;
    for (auto& s : work.block_starts) {
        const auto last = sofar;
        sofar += s;
        s = last;
    }
 
    sanisizer::resize(work.buffer, num_cells
#ifdef SCRAN_QC_TEST_INIT
        , SCRAN_QC_TEST_INIT
#endif
    );
 
    sanisizer::resize(work.block_offsets, num_blocks
#ifdef SCRAN_QC_TEST_INIT
        , SCRAN_QC_TEST_INIT
#endif
    );
}
 
template<typename Value_, typename Block_, typename Float_>
std::vector<ChooseFilterThresholdsResults<Float_> > choose_filter_thresholds_blocked(
    const std::size_t num_cells,
    const Value_* const metrics,
    const Block_* const block,
    const std::size_t num_blocks,
    ChooseFilterThresholdsBlockedWorkspace<Float_>& workspace,
    const ChooseFilterThresholdsOptions& options
) {
    std::vector<ChooseFilterThresholdsResults<Float_> > output;
    output.reserve(num_blocks);
    process_blocks_for_choose_filter_thresholds(
        num_cells,
        metrics,
        block,
        num_blocks,
        workspace,
        [&](const std::size_t len, Float_* const ptr) -> void {
            output.push_back(choose_filter_thresholds_internal<Float_>(len, ptr, options));
        }
    );
    return output;
}
 
template<typename Value_, typename Block_, typename Float_, class Function_>
void process_blocks_for_choose_filter_thresholds(
    const std::size_t num_cells,
    const Value_* const metrics,
    const Block_* const block,
    const std::size_t num_blocks,
    ChooseFilterThresholdsBlockedWorkspace<Float_>& workspace,
    Function_ fun 
) {
    assert(num_cells == workspace.buffer.size());
    assert(num_blocks == workspace.block_starts.size());
 
    auto& buffer = workspace.buffer;
    const auto& starts = workspace.block_starts;
    auto& offsets = workspace.block_offsets;
    std::copy(starts.begin(), starts.end(), offsets.begin());
    for (I<decltype(num_cells)> i = 0; i < num_cells; ++i) {
        auto& pos = offsets[block[i]];
        buffer[pos] = metrics[i];
        ++pos;
    }
 
    for (I<decltype(num_blocks)> g = 0; g < num_blocks; ++g) {
        fun(offsets[g] - starts[g], buffer.data() + starts[g]);
    }
}
 
 
 
template<bool lower_, typename Float_>
std::vector<Float_> extract_filter_thresholds(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;
}
}
 
#endif