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>
#include <cassert>
#include <array>
 
#include "quickstats/quickstats.hpp"
#include "sanisizer/sanisizer.hpp"
 
#include "utils.hpp"
 
namespace scran_norm {
 
struct ChoosePseudoCountOptions {
    double quantile = 0.05;
 
    double max_bias = 0.1;
 
    double min_value = 1;
};
 
template<typename Float_>
Float_ choose_pseudo_count_raw(const std::size_t num, Float_* const size_factors, const ChoosePseudoCountOptions& options) {
    // Avoid problems with zeros.
    I<decltype(num)> actual_num = 0;
    for (I<decltype(num)> i = 0; i < num; ++i) {
        const auto val = size_factors[i];
        if (std::isfinite(val) && val > 0) {
            if (i != actual_num) {
                size_factors[actual_num] = val;
            }
            ++actual_num;
        }
    }
 
    if (actual_num == 0) {
        return options.min_value;
    }
 
    if (options.quantile < 0 || options.quantile >= 0.5) {
        throw std::runtime_error("'quantile' should be in [0, 0.5)");
    }
    quickstats::MultipleQuantilesFixedNumber<Float_> comp(actual_num, std::vector<double>{ options.quantile, 1.0 - options.quantile });
    std::array<Float_, 2> sf;
    comp(size_factors, [&](std::size_t i, Float_ val) -> void {
        sf[i] = val;
    });
    const Float_ lower_sf = sf[0];
    const Float_ upper_sf = sf[1];
 
    // Very confusing formulation in Equation 3, but whatever.
    const 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(const std::size_t num, const Float_* const size_factors, const ChoosePseudoCountOptions& options) {
    std::vector<Float_> buffer(size_factors, size_factors + num);
    return choose_pseudo_count_raw(num, buffer.data(), options);
}
 
}
 
#endif