choose_pseudo_count.hpp

Go to the documentation of this file.

#ifndef SCRAN_NORM_CHOOSE_PSEUDO_COUNT_HPP
#define SCRAN_NORM_CHOOSE_PSEUDO_COUNT_HPP
 
#include <algorithm>
#include <vector>
#include <cstddef>
 
namespace scran_norm {
 
struct ChoosePseudoCountOptions {
    double quantile = 0.05;
 
    double max_bias = 0.1;
 
    double min_value = 1;
};
 
namespace internal {
 
template<typename Float_>
Float_ find_quantile(Float_ quantile, std::size_t n, Float_* ptr) {
    double raw = static_cast<double>(n - 1) * quantile;
    std::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));
}
 
}
template<typename Float_>
Float_ choose_pseudo_count_raw(std::size_t num, Float_* size_factors, const ChoosePseudoCountOptions& options) {
    if (num <= 1) {
        return options.min_value;
    }
 
    // Avoid problems with zeros.
    decltype(num) counter = 0;
    for (decltype(num) i = 0; i < num; ++i) {
        auto val = size_factors[i];
        if (std::isfinite(val) && val > 0) {
            if (i != counter) {
                size_factors[counter] = val;
            }
            ++counter;
        }
    }
    num = counter;
 
    if (num <= 1) {
        return options.min_value;
    }
 
    double lower_sf, upper_sf;
    if (options.quantile == 0) {
        lower_sf = *std::min_element(size_factors, size_factors + num);
        upper_sf = *std::max_element(size_factors, size_factors + num);
    } else {
        lower_sf = internal::find_quantile(options.quantile, num, size_factors);
        upper_sf = internal::find_quantile(1 - options.quantile, num, size_factors);
    }
 
    // Very confusing formulation in Equation 3, but whatever.
    Float_ pseudo_count = (1.0 / lower_sf - 1.0 / upper_sf) / (8 * options.max_bias);
 
    return std::max(options.min_value, pseudo_count);
}
 
template<typename Float_>
Float_ choose_pseudo_count(std::size_t num, const Float_* size_factors, const ChoosePseudoCountOptions& options) {
    std::vector<Float_> buffer(size_factors, size_factors + num);
    return choose_pseudo_count_raw(num, buffer.data(), options);
}
 
}
 
#endif