SuggestCrisprQcFilters.hpp

Go to the documentation of this file.

#ifndef SCRAN_SUGGEST_CRISPR_QC_FILTERS_H
#define SCRAN_SUGGEST_CRISPR_QC_FILTERS_H
 
#include "../utils/macros.hpp"
 
#include <vector>
#include <cmath>
 
#include "utils.hpp"
#include "PerCellCrisprQcMetrics.hpp"
#include "ComputeMedianMad.hpp"
#include "ChooseOutlierFilters.hpp"
 
namespace scran {
 
class SuggestCrisprQcFilters {
public:
    struct Defaults {
        static constexpr double num_mads = 3;
    };
 
private:
    double num_mads = Defaults::num_mads;
 
public:
    SuggestCrisprQcFilters& set_num_mads(double n = Defaults::num_mads) {
        num_mads = n;
        return *this;
    }
 
public:
    struct Thresholds {
        std::vector<double> max_count;
 
    public:
        template<bool overwrite = true, typename Float, typename Integer, typename Output>
        void filter(size_t n, const PerCellCrisprQcMetrics::Buffers<Float, Integer>& buffers, Output* output) const {
            if (max_count.size() != 1) {
                throw std::runtime_error("should use filter_blocked() for multiple batches");
            }
            filter_<overwrite>(n, buffers, output, [](size_t i) -> size_t { return 0; });
        }
 
        template<typename Output = uint8_t>
        std::vector<Output> filter(const PerCellCrisprQcMetrics::Results& metrics) const {
            std::vector<Output> output(metrics.sums.size());
            filter(output.size(), metrics.buffers(), output.data());
            return output;
        }
 
    public:
        template<bool overwrite = true, typename Block, typename Float, typename Integer, typename Output>
        void filter_blocked(size_t n, const Block* block, const PerCellCrisprQcMetrics::Buffers<Float, Integer>& buffers, Output* output) const {
            if (block) {
                filter_<overwrite>(n, buffers, output, [&](size_t i) -> Block { return block[i]; });
            } else {
                filter<overwrite>(n, buffers, output);
            }
        }
 
        template<typename Output = uint8_t, typename Block>
        std::vector<Output> filter_blocked(const PerCellCrisprQcMetrics::Results& metrics, const Block* block) const {
            std::vector<Output> output(metrics.sums.size());
            filter_blocked(output.size(), block, metrics.buffers(), output.data());
            return output;
        }
 
    private:
        template<bool overwrite, typename Float, typename Integer, typename Output, typename Function>
        void filter_(size_t n, const PerCellCrisprQcMetrics::Buffers<Float, Integer>& buffers, Output* output, Function indexer) const {
            for (size_t i = 0; i < n; ++i) {
                auto b = indexer(i);
                auto candidate = buffers.max_proportion[i] * buffers.sums[i];
 
                if (quality_control::is_less_than<double>(candidate, max_count[b])) {
                    output[i] = true;
                    continue;
                }
 
                if constexpr(overwrite) {
                    output[i] = false;
                }
            }
 
            return;
        }
    };
 
public:
    template<typename Float, typename Integer>
    Thresholds run(size_t n, const PerCellCrisprQcMetrics::Buffers<Float, Integer>& buffers) const {
        return run_blocked(n, static_cast<int*>(NULL), buffers);
    }
 
    Thresholds run(const PerCellCrisprQcMetrics::Results& metrics) const {
        return run(metrics.max_proportion.size(), metrics.buffers());
    }
 
public:
    template<typename Block, typename Float, typename Integer>
    Thresholds run_blocked(size_t n, const Block* block, const PerCellCrisprQcMetrics::Buffers<Float, Integer>& buffers) const {
        Thresholds output;
        std::vector<int> starts;
        std::vector<Block> subblock;
        if (block) {
            starts = ComputeMedianMad::compute_block_starts(n, block);
            subblock.reserve(n);
        }
 
        // Subsetting to the observations in the top 50% of proportions.
        std::vector<double> workspace(n);
        std::vector<double> subvalues;
        subvalues.reserve(n);
 
        {
            ComputeMedianMad meddler;
            meddler.set_median_only(true);
            auto prop_res = meddler.run_blocked(n, block, starts, buffers.max_proportion, workspace.data());
 
            if (block) {
                for (size_t i = 0; i < n; ++i) {
                    auto p = buffers.max_proportion[i];
                    if (quality_control::is_greater_than_or_equal_to(p, prop_res.medians[block[i]])) {
                        subvalues.push_back(p * buffers.sums[i]);
                        subblock.push_back(block[i]);
                    }
                }
            } else {
                for (size_t i = 0; i < n; ++i) {
                    auto p = buffers.max_proportion[i];
                    if (quality_control::is_greater_than_or_equal_to(p, prop_res.medians[0])) {
                        subvalues.push_back(p * buffers.sums[i]);
                    }
                }
            }
        }
 
        // Filtering on the max counts.
        {
            ComputeMedianMad meddler;
            meddler.set_log(true);
 
            const Block* bptr = NULL;
            starts.clear();
            if (block) {
                bptr = subblock.data();
                starts = ComputeMedianMad::compute_block_starts(subblock.size(), bptr);
            }
 
            auto sums_res = meddler.run_blocked(subvalues.size(), bptr, starts, subvalues.data(), workspace.data());
 
            ChooseOutlierFilters filter;
            filter.set_num_mads(num_mads);
            filter.set_upper(false);
            auto sums_filt = filter.run(std::move(sums_res));
 
            output.max_count = std::move(sums_filt.lower);
        }
 
        return output;
    }
 
    template<typename Block>
    Thresholds run_blocked(const PerCellCrisprQcMetrics::Results& metrics, const Block* block) const {
        return run_blocked(metrics.max_proportion.size(), block, metrics.buffers());
    }
};
 
}
 
#endif