ChooseOutlierFilters.hpp

Go to the documentation of this file.

#ifndef SCRAN_CHOOSE_OUTLIER_FILTERS_HPP
#define SCRAN_CHOOSE_OUTLIER_FILTERS_HPP
 
#include "../utils/macros.hpp"
 
#include <vector>
#include <limits>
#include <cmath>
 
#include "utils.hpp"
#include "ComputeMedianMad.hpp"
 
namespace scran {
 
struct ChooseOutlierFilters {
public:
    struct Defaults {
        static constexpr bool lower = true;
 
        static constexpr bool upper = true;
 
        static constexpr double num_mads = 3;
 
        static constexpr double min_diff = 0;
    };
 
private:
    bool lower = Defaults::lower;
    bool upper = Defaults::upper;
    double num_mads = Defaults::num_mads;
    double min_diff = Defaults::min_diff;
 
public:
    ChooseOutlierFilters& set_lower(bool l = Defaults::lower) {
        lower = l;
        return *this;
    }
 
    ChooseOutlierFilters& set_upper(bool u = Defaults::upper) {
        upper = u;
        return *this;
    }
 
    ChooseOutlierFilters& set_num_mads(double n = Defaults::num_mads) {
        num_mads = n;
        return *this;
    }
 
    ChooseOutlierFilters& set_min_diff(double m = Defaults::min_diff) {
        min_diff = m;
        return *this;
    }
 
private:
    double sanitize(double proposed, bool log) const {
        if (log) {
            if (std::isinf(proposed)) {
                if (proposed < 0) {
                    proposed = 0;
                }
            } else {
                proposed = std::exp(proposed);
            }
        } 
        return proposed;
    }
 
public:
    struct Thresholds {
        Thresholds(size_t nblocks = 0) : lower(nblocks), upper(nblocks) {}
        std::vector<double> lower;
 
        std::vector<double> upper;
 
    public:
        template<bool overwrite = true, typename Input, typename Output>
        void filter(size_t n, Input* input, Output* output) const {
            if (upper.size() > 1 || lower.size() > 1) {
                throw std::runtime_error("should use filter_blocked() for multiple batches");
            }
            filter_<overwrite>(n, input, output, [](size_t i) -> size_t { return 0; });
        }
 
        template<typename Output = uint8_t, typename Input>
        std::vector<Output> filter(size_t n, const Input* input) const {
            std::vector<Output> output(n);
            filter(n, input, output.data());
            return output;
        }
 
    public:
        template<bool overwrite = true, typename Block, typename Input, typename Output>
        void filter_blocked(size_t n, const Block* block, const Input* input, Output* output) const {
            if (block) {
                filter_<overwrite>(n, input, output, [&](size_t i) -> Block { return block[i]; });
            } else {
                filter<overwrite>(n, input, output);
            }
        }
 
        template<typename Output = uint8_t, typename Block, typename Input>
        std::vector<Output> filter_blocked(size_t n, const Block* block, const Input* input) const {
            std::vector<Output> output(n);
            filter_blocked(n, block, input, output.data());
            return output;
        }
 
    private:
        template<bool overwrite, typename Input, typename Output, typename Function>
        void filter_(size_t n, const Input* input, Output* output, Function indexer) const {
            for (size_t i = 0; i < n; ++i) {
                auto b = indexer(i);
 
                if (!lower.empty()) {
                    if (quality_control::is_less_than(input[i], lower[b])) {
                        output[i] = true;
                        continue;
                    }
                }
 
                if (!upper.empty()) {
                    if (quality_control::is_greater_than(input[i], upper[b])) {
                        output[i] = true;
                        continue;
                    }
                }
 
                if constexpr(overwrite) {
                    output[i] = false;
                }
            }
 
            return;
        }
    };
 
public:
    Thresholds run(ComputeMedianMad::Results x) const {
        size_t nblocks = x.medians.size();
 
        for (size_t b = 0; b < nblocks; ++b) {
            auto med = x.medians[b];
            auto mad = x.mads[b];
 
            double& lthresh = x.medians[b];
            double& uthresh = x.mads[b];
 
            if (std::isnan(med)) {
                if (lower) {
                    lthresh = std::numeric_limits<double>::quiet_NaN();
                }
                if (upper) {
                    uthresh = std::numeric_limits<double>::quiet_NaN();
                }
            } else {
                auto delta = std::max(min_diff, num_mads * mad);
                if (lower) {
                    lthresh = sanitize(med - delta, x.log);
                }
                if (upper) {
                    uthresh = sanitize(med + delta, x.log);
                }
            }
        }
 
        Thresholds output;
        if (lower) {
            output.lower = std::move(x.medians);
        }
        if (upper) {
            output.upper = std::move(x.mads);
        }
        return output;
    }
};
 
}
 
#endif