ChoosePseudoCount.hpp

Go to the documentation of this file.

#ifndef SCRAN_CHOOSE_PSEUDO_COUNT_HPP
#define SCRAN_CHOOSE_PSEUDO_COUNT_HPP
 
#include <algorithm>
#include <vector>
#include <cmath>
 
namespace scran {
 
class ChoosePseudoCount {
public:
    struct Defaults {
        static constexpr double quantile = 0.05;
 
        static constexpr double max_bias = 0.1;
 
        static constexpr double min_value = 1;
    };
 
private:
    double quantile = Defaults::quantile;
    double max_bias = Defaults::max_bias;
    double min_value = Defaults::min_value;
 
public:
    ChoosePseudoCount& set_quantile(double q = Defaults::quantile) {
        quantile = q;
        return *this;
    }
 
    ChoosePseudoCount& set_max_bias(double b = Defaults::max_bias) {
        max_bias = b;
        return *this;
    }
 
    ChoosePseudoCount& set_min_value(double v = Defaults::min_value) {
        min_value = v;
        return *this;
    }
 
public:
    static double find_quantile(double quantile, size_t n, double* ptr) {
        double raw = static_cast<double>(n - 1) * quantile;
        size_t index = std::ceil(raw);
        std::nth_element(ptr, ptr + index, ptr + n);
        double upper = *(ptr + index);
        std::nth_element(ptr, ptr + index - 1, ptr + index);
        double lower = *(ptr + index - 1);
        return lower * (index - raw) + upper * (raw - (index - 1));
    }
public:
    double run(size_t n, const double* size_factors, double* buffer) const {
        if (n <= 1) {
            return min_value;
        }
 
        // Avoid problems with zeros.
        size_t counter = 0;
        for (size_t i = 0; i < n; ++i) {
            if (size_factors[i] > 0) {
                buffer[counter] = size_factors[i];
                ++counter;
            }
        }
        n = counter;
 
        if (n <= 1) {
            return min_value;
        }
 
        double lower_sf, upper_sf;
        if (quantile == 0) {
            lower_sf = *std::min_element(buffer, buffer + n);
            upper_sf = *std::max_element(buffer, buffer + n);
        } else {
            lower_sf = find_quantile(quantile, n, buffer);
            upper_sf = find_quantile(1 - quantile, n, buffer);
        }
 
        // Very confusing formulation in Equation 3, but whatever.
        double pseudo_count = (1 / lower_sf - 1 / upper_sf) / (8 * max_bias);
 
        return std::max(min_value, pseudo_count);
    }
 
    double run(size_t n, const double* size_factors) const {
        std::vector<double> buffer(n);
        return run(n, size_factors, buffer.data());
    }
};
 
}
 
#endif